2.1.

Characterization of Spatially Correlated Noise

First, we consider a two-dimensional noise map $N(\mathbf{x})$, where $\mathbf{x}$ denotes the position of the pixels. Pixel regions with signal from the object of interest are excluded. The statistical properties of the noise are assumed to be uniform in the noise image, which appears to be valid in the interferometric image (We discuss this in Sec. 2.2.) In most of the literature, the noise in the radio interferometric image is quantified and reported with the rms of the noise map $N(\mathbf{x})$:

where the brackets denote the expected value for each pixel, which is practically estimated by averaging over the noise map. The mean of the noise in the image $\mu \equiv ⟨N(x)⟩\approx 0$, since the most of the noise represented by the system temperature $T_{\mathrm{sys}}$ is not correlated in a pair of antennas and the power of the noise does not appear in the correlator output of interferometers, such as ALMA. Extended background emission, such as the cosmic microwave background, is resolved out without total power observation. However, these sources of noise contribute to the random noise associated with visibility measurements that propagate into the noise on the image by the Fourier transform. In the reminder of this paper, we assume that the mean of the noise map to be zero or already have been subtracted in other cases, and thus the rms and the standard deviation of the noise can be used interchangeably. Figure 1 shows the example ALMA band 7 noise map and its histogram from the observation targeting the hyper luminous infrared galaxy BRI 1335-0417 at redshift of 4.4, which will be used in the later sections. The noise map is created by eliminating astronomical sources by 4 sigma clipping as well as removing pixels adjacent to these clipped regions out to 3 times the full-width of half-maximum (FWHM) of the synthesized beam. The histogram of the pixel values in the noise map is well fitted by a Gaussian function.

When noise can be assumed to be Gaussian, the statistical and correlation properties of noise are fully quantified by the noise ACF:²

where the expected value is estimated by averaging all pairs of pixels with the relative distance ${\mathbf{x}}_{i,j}=(i,j)$ in the noise image. The value of the ACF noise at zero lag, ${\mathbf{x}}_{i,j}=\mathbf{0}$, is equal to the variance of the noise as $\xi (0)=⟨N(\mathbf{x}{)}^2⟩={\sigma }_N^2$. When the noise has no interpixel correlations, the noise ACF becomes To evaluate the statistical uncertainty of the derived noise ACF, we first considered the number of independent pixel pairs $N_{\text{pair}}$ in the number of all available pairs $N_{\text{pair}}^{\prime }$ used to evaluate the bracket in Eq. (2), since the pixels within a beam area are expected to be strongly correlated and not independent. We estimated the number of independent pixel pairs $N_{\text{pair}}$ as the ratio of the number of all pixel pairs $N_{\text{pair}}^{\prime }$ and the number of pixels in the beam (beam area in pixels) $A_{\text{beam}}$. [The beam area in pixels is typically estimated by $2\pi b_{\mathrm{maj}}{b}_{\min }/8\ln 2$, where $b_{\mathrm{maj}}$ and $b_{\min }$ are the major and minor FWHMs of the mainlobe of the synthesized beam (the “CLEAN” beam)]: Then the associated statistical uncertainty of the noise ACF $\mathrm{\Delta }\xi ({\mathbf{x}}_{i,j})$ is calculated as the usual standard error of the mean but with an independent sample size $N_{\text{pair}}$, that is, the standard deviation of the multiplication of the values across all pairs of pixels with separation ${\mathbf{x}}_{i,j}$ divided by the root of the number of independent pixel pairs $N_{\text{pair}}$:

Figure 2 shows the results of the noise ACF [Eq. (2)] computed for the noise map shown in Fig. 1, and the synthesized beam of the observation, both of which are normalized so that the central value is one. The noise ACF shows a pattern similar to that of the synthesized beam, with a correlation signal near the center and a correlation signal away from the center corresponding to the main lobe and side lobe of the synthesized beam, respectively. This suggests that most of the correlation of the noise originates from the discrete Fourier transform involved in the interferometric imaging, which is illustrated in the following section.

Note that the noise ACF is measured for the noise maps of the band 7 continuum image (shown in Fig. 2), [Cii] line (velocity integrated over the velocity range of $-400$ to $400\mathrm{km}{\mathrm{s}}^{-1}$, where the velocity is computed with respect to the redshifted [Cii] line frequency with the galaxy’s redshift of 4.4074⁹) moment 0 map, and the [Cii] line cubes (each velocity channel map). These maps are primary beam uncorrected and cleaned images produced by the CLEAN algorithm in CASA (see details in Ref. 8).

2.2.

Origin of the Noise Correlation

In interferometric observations, measurable quantities are visibility (Fourier amplitude and phase) of the astronomical image at the given spatial frequencies $(u,v)=\mathbf{D}/\lambda$, which are related to the antenna baseline vector $\mathbf{D}$ (separation vector of pairs of antennas) projected onto the plane of the sky and the observed wavelength $\lambda$. The image is then computed by the Fourier transform of the measured visibility.

To explore the origin of the noise correlation in the image, seen in (Fig. 2), we start with the ideal case in which the observation measures visibilities at all spatial frequencies $(u,v)$. The visibility of the source of interest is $V(u,v)$, which is the Fourier transform of the true flux distribution of the source in the image, $\widehat{S}(x,y)=\mathrm{FT}[V(u,v)]$, where FT denotes the Fourier transform. A measurement of $V(u,v)$ usually involves uncorrelated random noise, which we describe with the random variable ${\widehat{N}}_{\mathrm{vis}}(u,v)$ with zero mean. We assume that the statistical property of the random variable ${\widehat{N}}_{\mathrm{vis}}(u,v)$ is uniform as a function of $u$ and $v$, that is, the system noise temperature is the same for all antennas. The image obtained $I(x,y)$ is the Fourier transform of the measurement $V(u,v)+{\widehat{N}}_{\mathrm{vis}}(u,v)$

$\widehat{N}(x,y)=\mathrm{FT}({\widehat{N}}_{\mathrm{vis}}(u,v))$ is the noise component of the image, which is a random variable with zero mean. (The Fourier transform of the random variable with zero mean is also random variable with zero mean.) The noise component of the image $\widehat{N}(x,y)$ is due to the random noise associated with visibility measurements ${\widehat{N}}_{\mathrm{vis}}(u,v)$, and the resulting noise map $N(x,y)=\widehat{N}(x,y)$ is not spatially correlated in the ideal case where all spatial frequencies are measured.

In practice, visibilities are measured only at the limited spatial frequencies $\{(u_1,v_1),(u_2,v_2),\dots ,(u_M,v_M)\}$ ($uv$ coverage). The spatial transfer function $W(u,v)$ is used to describe the spatial frequencies $(u,v)$, in which we measure the visibility. This function $W(u,v)$ is nonzero if the visibility at $(u,v)$ is actually measured, which can be expressed as

where $\delta$ is the Dirac delta function. The synthesized beam $b(x,y)$ is the Fourier transform of the spatial transfer function $W(u,v)$, $b(x,y)=\mathrm{FT}[W(u,v)]$. The resulting image $I(x,y)$, decomposed as the signal $S(x,y)$ from the source and noise map $N(x,y)$, is where $*$ represents convolution. As the noise correlation pattern (noise ACF) and the synthesized beam show a similar pattern in Fig. 2, the noise component of the image $N(x,y)$ is the convolution product of the random variable $\widehat{N}(x,y)$ and the synthesized beam $b(x,y)$. Because of this, the noise in the image is well behaved; in particular, its statistical properties are uniform over the image, as assumed in Sec. 2.1, when measuring the noise ACF.

For convenience, by replacing the sky position of $(x,y)$ with the pixel position $\mathbf{x}$, the noise map of the image in Eq. (8) is written as

The autocorrelation of the noise map then becomes² where we used the noise ACF property of uncorrelated noise $\widehat{N}$ [Eq. (3)] for the fourth equality and we have defined $\alpha ({\mathbf{x}}_{i,j})$ as beam autocorrelation:

Equation (10) implies that the noise ACF is related to the ACF of the synthesized beam with a constant multiplicative factor, which is the variance of the noise. In Fig. 3, we compare the ACF of noise and that of the synthesized beam, along with the residuals (noise ACF minus beam ACF). Although the noise ACF and beam ACF show a common characteristic pattern, they do not completely coincide. The difference of the two ACF shows an extended weak positive correlation and a relatively large negative around the main beam in the residual. This disagreement is likely due to not only (1) the remaining contamination by the emission from the sources, but also (2) the process involved in the imaging of the visibility measurements. These are discussed in detail in Appendix A comparing Fig. 3 obtained from the actual data with the one (Fig. 13) obtained from the simulated data with a similar observational setup and realistic noise in the visibilities but without emission in the sky. Note that our interest is on the statistical property of the noise in the image plane, which we characterize by the noise ACF including effects of contamination from the source and the imaging process.

Due to the limited spatial frequency coverage, the synthesized beam $b(x,y)$ has a complex structure with sidelobes that extend from the center to a large radius. The flux from the source is spread out by the side lobes to the distant pixels in the image. The CLEAN algorithm, which is most commonly used in radio imaging, deconvolves the beam pattern $b(x,y)$ for signals with high $S/N$ $(S(x,y)/{\sigma }_N)>3$) and replaces it with a CLEAN beam without side lobes (a Gaussian that approximates the mainlobe of the synthesized beam). The CLEAN algorithm successfully suppresses the influence of the sidelobe and produces a high-fidelity image, but cannot remove the spatial correlations that exist in stochastic noise $N(x,y)$. Therefore, it is important to evaluate their effects on image analysis and signal detection, which we will describe in Sec. 3.

3.1.

Contribution of the Correlated Noise to the Statistical Uncertainty in the Measured Flux

The most fundamental measurement of astronomy is the total flux distributing over some sky region in the images, which are measured by summing the pixel values over the region of interest (i.e., aperture photometry in optical astronomy). In particular, at the submillimeter band of ALMA, the flux of the continuum emission arising primarily from thermal dust, and line emission and absorption by the various atomic and molecular gases are used to estimate the physical properties of the interstellar medium (e.g., dust mass, gas mass, and the energy source of the ionization or excitation). It is important to estimate the uncertainty of the measured quantities. As shown in Fig. 2, the noise in interferometric images correlates significantly between pixels, making the estimation of the noise in the integrated flux difficult. In the previous literature, statistical uncertainties of integrated fluxes were estimated by one of the two methods: (1) randomly placing identical apertures in the noise region of the image, measuring the sum within each aperture and then adopting the rms as the noise in the sum of pixels in the aperture⁵ and (2) assuming that the regions in the image separated with a beam size do not correlate and adopting ${\sigma }_N{A}_{\text{beam}}\sqrt{N_{\text{beam}}}$, where $N_{\text{beam}}$ is the number of beams (independent regions) in the aperture.³ $N_{\text{beam}}$ is estimated as $A_{\text{aperture}}/A_{\text{beam}}$, where $A_{\text{aperture}}$ and $A_{\text{beam}}$ are the aperture area and the CLEAN beam area in pixels, respectively. For convenience, we call methods (1) and (2) “random aperture method” and “independent beam method,” respectively, in this paper.

In the “independent beam method” (${\sigma }_N{A}_{\text{beam}}\sqrt{N_{\text{beam}}}$), the factor ${\sigma }_N{A}_{\text{beam}}$ is the standard deviation of the sum of noise in individual pixels within a beam assuming that the noise perfectly correlates within a beam. Then the standard deviation of the sum of the noise of each independent beam area in the aperture is computed by scaling by the square root of the number of independent beams $N_{\text{beam}}$ within the aperture. The terms $A_{\text{beam}}$ and $N_{\text{beam}}$ in the ${\sigma }_N{A}_{\text{beam}}\sqrt{N_{\text{beam}}}$ denote just the number of data points to be summed. Therefore, we caution readers that ${\sigma }_N{A}_{\text{beam}}\sqrt{N_{\text{beam}}}$ has the same unit as ${\sigma }_N$. Most interferometric maps and measured ${\sigma }_N$ are in brightness units, e.g., $\mathrm{Jy}{\text{beam}}^{-1}\mathrm{km}{\mathrm{s}}^{-1}$ or $\mathrm{Jy}{\text{beam}}^{-1}$. So we need to divide $A_{\text{beam}}$ to compare with the integrated flux or spectral flux density, e.g., $\mathrm{Jy}\mathrm{km}{\mathrm{s}}^{-1}$ or Jy. ${\sigma }_N{A}_{\text{beam}}\sqrt{N_{\text{beam}}}$ is a factor of $A_{\text{beam}}$ different from ${\sigma }_N\sqrt{N_{\text{beam}}}$ described by Alatalo et al.³ due to the unit difference where they assume the quantity in the unit of flux.

This section introduces how to derive the statistical uncertainty associated with the spatially integrated flux directly from the computed noise ACF. We consider adding all the pixel values at pixel positions $\mathbf{x}$ within the sky region of interest $S$. The random noise $N(\mathbf{x})$ in the map is characterized by the noise ACF $\xi ({\mathbf{x}}_{i,j})$. The $1\sigma$ statistical uncertainty associated with the summed value within the pixel region $S$ ${\sigma }_{\mathrm{int}}$ can be estimated as

where $N_{\mathrm{pix}}$ is the number of pixels in the region $S$. Var and Cov indicate variance and covariance, respectively. The second term of the last line is the sum of the noise ACF for all possible pixel separation vectors ${\mathbf{x}}_{i,j}$ between two pixels within the region $S$. If the noise does not have an interpixel correlation, the second term becomes zero, resulting in ${\sigma }_{\mathrm{int}}={\sigma }_{\mathrm{N}}\sqrt{N_{\mathrm{pix}}}$. The method Sun et al.¹ proposed to estimate ${\sigma }_{\mathrm{int}}$ is equivalent to Eq. (12) but approximately substituting the covariance $\xi ({\mathbf{x}}_{i,j})$ by ${\sigma }_N^{2}b({\mathbf{x}}_{i,j})$. However, $\xi ({\mathbf{x}}_{i,j})={\sigma }_N^2\alpha ({\mathbf{x}}_{i,j})$ as shown in Eq. (10).

As an illustration, we used the spatially resolved [Cii] moment 0 map of BRI 1335-0417 taken by ALMA, shown in Fig. 4, where the emission spreads over multiple pixels. We calculated the noise in the integrated flux measured using a variety of apertures $S$ with different sizes. The largest aperture is a dotted line shown in Fig. 4. Figure 5 shows the computed noise from the measured noise ACF compared to the previously used “random aperture” and “independent beam” methods. The noise calculated from the ACF is in excellent agreement with the random aperture method, whereas the independent beam method tends to overestimate ${\sigma }_{\mathrm{int}}$ at smaller apertures (considering the limiting case that the aperture is a single pixel, the obtained value by “independent beam method”, ${\sigma }_N\sqrt{A_{\text{beam}}}$) is clearly different from the fiducial value ${\sigma }_N$) and underestimate ${\sigma }_{\mathrm{int}}$ at large apertures, showing that the assumption of “independent beam” is oversimplified. When the field of view is small (the field of view becomes smaller at higher frequency bands in ALMA) or the aperture area becomes larger, the random aperture method cannot place apertures randomly in the limited area of the emission-free region, and thus the standard error of the estimate increases, as shown in the blue shaded region in Fig. 5. The proposed method of calculating ${\sigma }_{\mathrm{int}}$, however, can provide the best estimate by exploiting all available data to estimate the noise ACF. The total [Cii] flux measured with the aperture shown in Fig. 4 is 29.51 $\pm 1.05\mathrm{Jy}\mathrm{km}{\mathrm{s}}^{-1}$ ($1\sigma$ statistical uncertainty calculated from the noise ACF).

To demonstrate the significant effect of spatially correlated noise, Fig. 6 shows the noise variance ${\sigma }_N^2$ in the integrated flux calculated from the ACF along with the contributions of the noise variance of individual pixels in the aperture [the first term in Eq. (12), this is the value we obtain if we are unaware of the noise correlation] and the interpixel correlation in the aperture [the second term in Eq. (12)]. In the case of an aperture with two pixels, it is not mathematically allowed for the first term to exceed the second term. However, as the number of pixels in the aperture $N_{\mathrm{pix}}$ increases, the contribution of the second term dominates over the first term because the number of the first terms is $N_{\mathrm{pix}}$, whereas the number of the second terms is $N_{\mathrm{pix}}(N_{\mathrm{pix}}-1)$. Ignoring the noise correlation will lead to a significant underestimation of the integrated flux uncertainty.

Figure 7 further divides the variance due to noise correlation into the effects of the mainlobe (correlation due to the mainlobe of the synthesized beam) and the sidelobe (long-range pixel correlation due to the sidelobe of the synthesized beam). To illustrate the significance of the long-range correlation conservatively, we define short-range/long-range correlation components of the noise ACF inside/outside the ellipse with the beam FWHM in radius, respectively. After $N_{\mathrm{pix}}$ exceeds 200, the effect of the sidelobe becomes significant. This explains the deviation of the estimate by the independent beam method from the true value (computed by the noise ACF or the random aperture method; see Fig. 5); the number of independent beams $N_{\text{beam}}$ is estimated by the area of the CLEAN beam, and the long-range correlation due to the sidelobe is not properly taken into account.

Another important measurement in astronomy is the shape of the spectrum integrated over a certain region of interest. Similar to deriving the noise in the integrated flux, we can derive the underlying noise in the spatially integrated spectrum using the noise ACFs computed for every velocity channel of the data cube. Figure 8 shows the spatially integrated [Cii] spectrum enclosed by the aperture shown in Fig. 4 with $1\sigma$ and $3\sigma$ noise levels. The [Cii] spectrum of BRI 1335-0417 is well described by a single Gaussian without deviations from the Gaussian above $3\sigma$. Because emission lines from interesting astronomical phenomena, such as outflows and tidal tails, are faint, it is important to accurately estimate the noise, otherwise, it will lead to false detections.

3.2.

Simulating the Noise Maps

In image analysis, generating random noise based on the statistical properties of the noise is useful to assess the significance of the results. In this section, we describe how to generate random noise based on the noise ACF, which fully characterizes the correlated Gaussian noise. Using an example of simulated noise, we demonstrate how ignoring correlated noise leads to misinterpretation of results.

Once the noise ACF is measured, the noise at the ${\mathbf{x}}_{i,j}$ positions in images with a size of $M\times M$ pixels can be generated randomly by the joint probability distribution, the probability that $N({\mathbf{x}}_{i,j})$ takes the value in small intervals ($N_{i,j}+\mathrm{d}{N}_{i,j}$) given by¹⁰

where $B^{-1}$ is the inverse of the matrix $B$ defined by

Figure 9 shows the comparison of the noise of the observed data, the noise randomly generated from the measured noise ACF using the joint Gaussian probability distribution [Eq. (13); we used the MULTIVARIATE_NORMAL function from the scipy package¹¹], and the spatially uncorrelated Gaussian noise, all of which have the same standard deviation ${\sigma }_N$. The observed noise and the randomly generated noise using the noise ACF are qualitatively similar, whereas the spatially correlated noise and the spatially uncorrelated noise look completely different, illustrating how dangerous it is to assume naively that noise is uncorrelated in image analysis.

In the literature, emission-free regions have been extracted and used as realistic noise maps to estimate the statistical uncertainty associated with the parameters derived by the model fitting.⁶ The field of view of the interferometric observation is generally small, preventing us from obtaining a sufficient number of independent noise maps or cubes to conduct Monte Carlo experiments. (As the $T_{\mathrm{sys}}$ is not strongly variable across a spectral window, we may use the line-free channels with cautions of the channels affected by emission lines in the atmosphere and the usually negligible variation of spatial frequencies $(u,v)$ with spectral frequency.) The proposed method can generate random noise repeatedly from the measured statistical properties of the noise (the noise ACF) in the interferometric image with the best precision limited by the available area of the emission-free region.

There are many potential applications of generating noise. We illustrate the significance of the correlated noise using one example: the image analysis done in Ref. 8, which expands the image with series of logarithmic spirals and identifies the dominant spiral structure in the image. In Fig. 10, we show the Fourier spectrum of the logarithmic spirals in the [Cii] intensity map of the BRI 1335-0417 shown in Fig. 4 (black solid line) and the underlying noise spectra estimated from the simulated noise maps from noise ACF (blue shaded region) and the simulated noise maps with no spatial correlation but the same standard deviation ${\sigma }_N$ (orange shaded region). The spectrum shows the amplitude of the logarithmic spiral with $m$ arms ($m$-fold symmetry) as a function of a variable $p$, which is related to the pitch angle of $\alpha =\arctan (-m/p)$. To calculate the spectrum, the image (Fig. 4) was first deprojected to be viewed face-on with an inclination angle of 37.8 deg and a position angle of 4.5 deg, where the package SCIKIT-IMAGE performed rotation and stretching for the deprojection (which induces an additional noise correlation in the image). Then the amplitude of each Fourier component of the logarithmic spiral $m$ and $\alpha$ is calculated (see Equation S2 in Ref. 8). The noise spectrum was measured by generating 300 noise maps, calculating the amplitude of each Fourier component in the same way as the image, and taking the 84th percentiles. The noise spectrum computed by assuming the uncorrelated Gaussian noise map significantly underestimates the true noise spectrum, which could lead to the false detection of statistically insignificant structures, such as the second peak in $m=2$ and multiple peaks in $m=3,4$. The estimated noise spectra are higher than those reported in the original paper because the image stretching and rotation were not applied to the noise maps in the original paper, resulting in an underestimation of the noise level. However, the result of the original article does not change except for the updated statistical significance of each peak in $m=1,2,3,4$ being $2.1\sigma$, $3.5\sigma$, $1.7\sigma$, and $1.6\sigma$, respectively. The $m=1$ peak with $2.1\sigma$ corresponds to the fact that the northern arm is longer than the southern arm.

3.3.

Fitting a Model to an Interferometric Image Under Spatially Correlated Noise

Fitting an analytic model of the intensity distribution to the 2D observed image is one of the most fundamental tasks in astronomy, allowing us to extract information to characterize astronomical sources from images efficiently. For a simple example, fitting 2D Gaussian to the intensity image provides essential information regarding the source’s position, size, and brightness. This section introduces the construction of a covariance matrix from the measured noise ACF, which is required to calculate ${\chi }^2$, to properly estimate the statistical uncertainty on the parameters derived by fitting the model to the observed interferometric image.

Consider fitting a 2D model $I_{i,j}(\theta )$ to an observed image $I_{i,j}^{\prime }$ with $M\times M$ pixels, where $\theta$ denotes a set of model parameters. We denote the flattened model image and the observed image by $I_m(\theta )$ and $I_m^{\prime }$, respectively, with the flat indices $m=(1,2,3,\dots ,M^2)$ that correspond to the coordinate indices $(i,j)=((\mathrm{1,1}),\dots ,(M,1),(\mathrm{1,2}),\dots ,(M,2),\dots ,(M,M))$ in the images. The formal expression to calculate ${\chi }^2$ is

where $C^{-1}$ is the inverse of the covariance matrix $C$ with the size of $M^2\times M^2$, and $\mathbf{r}(\theta )$ is the residual vector with the elements $r_m(\theta )=I_m^{\prime }-I_m(\theta )$. When noise has no correlation between pixels, the covariance matrix $C={\sigma }_N^{2}E$ ($E$ is the identity matrix) and ${\chi }^2$ leads to When fitting the 2D model to the observed interferometric image, in the literature, the noise correlation has been ignored to estimate the statistical uncertainties on the derived parameters using Eq. (16) instead of Eq. (15); therefore, the statistical uncertainties reported are, in most cases, significantly underestimated.

We can construct the covariance matrix $C$ in Eq. (15) from measured noise ACF by

Here recall that $m$ (and $m^{\prime }$) are flattened indices with a one-to-one relationship with the pixel coordinate $m\to (i,j)$. As $m$ and $m^{\prime }$ are exhangable, the covariance matrix $C$ is a real-valued symmetric matrix, whose inverse can be efficiently calculated using the LINALG.PINVH function in the Scipy package.

Now, the question becomes; what effect does the noise correlation have on the estimated statistical uncertainties of the model parameters derived from the fitting process? To explore this, we fit the 2D Gaussian model to a noiseless 2D Gaussian image and sampled the posterior distribution of the model parameters: total flux of the Gaussian distribution $L$, center of the Gaussian $(x,y)$, major axis and minor axis of the Gaussian $({\sigma }_{\mathrm{maj}},{\sigma }_{\min })$, and position angle (P.A.) from the north to the major axis (counterclockwise). The noise effect perturbing the best-fitting parameters is effectively taken into account by the covariance matrix $C$ in the likelihood [Eq. (15)]. We compared the posterior distributions that resulted in two cases where (1) the noise has no spatial correlation [using Eq. (16)] and (2) the noise has spatial correlation characterized by the noise ACF shown in Fig. 2 [using Eq. (15), and the covariance matrix constructed from the noise ACF]. Note that in both cases, the noise variance in individual pixels is the same, but the latter case has spatial correlation. We used emcee¹² to sample posterior distributions with logarithmic likelihood $L\propto -0.5{\chi }^2$ and a uniform prior. Figure 11 shows the posterior distributions of the parameters of the Gaussian model in two cases. Accounting for the spatial correlation of the noise results in a significantly larger confidence interval than that obtained assuming that there is no spatial correlation of the noise. The model parameters would appear to be constrained too well if we ignored the noise correlation in the fitting process.

Figure 12 compares the posterior distribution sampled with the ${\chi }^2$ covariance matrix (the same as Fig. 11) with the distribution estimated by Monte Carlo resampling. The resampling is performed by repeatedly obtaining bestfitting parameter values after adding the randomly generated correlated noise to the noiseless image (the method described in Sec. 3.2). This confirms that the posterior distributions estimated by the two approaches agree well.

In addition to the 2D fitting problem when investigating emission distributions, many recent studies fit 3D model cubes to observe 3D data cubes (R.A., Dec., and velocity) thanks to the improvement of interferometric imaging in terms of spatial resolution and sensitivity. Radio interferometer has sufficiently high-frequency resolution ($0.01\mathrm{km}{\mathrm{s}}^{-1}$ or $R=3\times {10}^7$ at 110 GHz for ALMA), which is usually further binned to increase the signal to noise. Therefore, the noise can be regarded as independent across the frequency (or velocity) channels, and the above method described for the 2D case can be applied to compute ${\chi }^2$ for the 3D case. One of the main examples of 3D cube fitting is the kinematical modeling of emission lines.

There are several publicly available modeling codes for this purpose, which assume the disk geometry of the emission line, parametrize the rotation curve, velocity, dispersion, and 3D geometry of the disk (position angle, inclination), and produce 3D model cubes that can be directly compared to observations taking into account the resolution and pixel binning (e.g., TIRIFIC,¹³ KINMS,¹⁴ 3D-BAROLO,¹⁵ GALPAK 3D,¹⁶ and QUBEFIT¹⁷). Under dynamical equilibrium, rotation curves are powerful tracers for the mass distribution of galaxies that can be further decomposed into black hole mass, stellar bulge, exponential disk, etc., depending on the data quality. These codes have been widely used to measure galaxy mass distributions and BH masses from distant to nearby galaxies (e.g., Refs. 8 and 14). However, these codes do not consider the correlation of the noise. As the number of data becomes large and the effect of noise correlation becomes significant, estimates of the statistical uncertainty on best-fitting parameters are too small.

Davis et al.⁴ estimate a covariance matrix from a point spread function and use it to calculate the ${\chi }^2$ values. They approximate the point spread function by a single Gaussian function ignoring the side-lobe, which is an oversimplification as seen in Fig. 7. We find that the side-lobe effect can dominate as an image region becomes bigger. As a result, the estimates of the statistical uncertainty of parameters can be underestimated. In addition, as seen in Sec. 2.2, the point spread function does not coincide with the actual noise correlation pattern in the image. Therefore, we recommend constructing the covariance matrix using Eq. (17) from the noise ACF, which fully characterizes the actual noise correlation pattern rather than from the point spread function.

Boizelle et al.⁶ proposed to estimate the formal ${\chi }^2$ by block averaging the data and model cube to form roughly beam-sized cells to avoid the need to calculate the covariance matrix mentioned above. They stated that the block averaging method does not fully mitigate the correlation between neighboring pixels with the presence of the long-range correlation as shown in Fig. 7. They estimated the final statistical uncertainty on the model parameters by Monte Carlo realisation using line-free channels.

The example calculation for the 2D fitting confirms that the two approaches, Monte Carlo resampling (in 2D or 3D),⁶ as well as the use of the covariance matrix during fitting are equivalent and should provide the correct estimation under spatially correlated noise.