Total Power Map to Visibilities (TP2VIS): Joint Deconvolution of ALMA 12m, 7m, and Total Power Array Data

TL;DR

The paper introduces TP2VIS, a new software package that converts single-dish ALMA maps into visibilities for joint deconvolution with other array data, improving imaging of extended and compact structures.

Contribution

TP2VIS is a novel tool that enables the conversion of single-dish maps into visibilities, facilitating joint deconvolution of ALMA data from multiple arrays.

Findings

TP2VIS improves uv coverage and deconvolution quality.

Reproduces model images within 5% error over large flux ranges.

Enhances imaging of both extended and compact astronomical structures.

Abstract

We present a new package for joint deconvolution of ALMA 12m, 7m, and Total Power (TP) data, dubbed ``Total Power Map to Visibilities (TP2VIS)". It converts a TP (single-dish) map into visibilities on the CASA platform, which can be input into deconvolvers (e.g., CLEAN) along with 12m and 7m visibilities. A manual is presented in the Github repository (https://github.com/tp2vis/distribute). Combining data from the different ALMA arrays is a driver for a number of science topics, namely those that probe size scales of extended and compact structures simultaneously. We test TP2VIS using model images, one with a single Gaussian and another that mimics the internal structures of giant molecular clouds. The result shows that the better uv coverage with TP2VIS visibilities helps the deconvolution process and reproduces the model image within errors of only 5% over two orders of magnitude in flux.

Tables

Table 1. Table 1: Beam Terminology

Number	Beam Name	Symbol
1	Primary beam of Total Power (TP) telescopes	${\Omega }_{\mathrm{TP}}$
2	Primary beam of Virtual Interferometer (VI) telescopes	${\Omega }_{\mathrm{vir}}$
3	Dirty beam (output from clean/tclean with niter=0)	${\Omega }_{\mathrm{dirty}}$
	Synthesized beam	${\Omega }_{\mathrm{syn}}$
	Deconvolution beam	${\Omega }_{\mathrm{dec}}$
	Point spread function	${\Omega }_{\mathrm{psf}}$
4	CLEAN beam (typically derived from Gaussian fit to ${\Omega }_{\mathrm{dirty}}$)	${\Omega }_{\mathrm{clean}}$
	Convolution beam	${\Omega }_{\mathrm{conv}}$
	Restoring beam	${\Omega }_{\mathrm{res}}$

Equations

$$B(l,m)=\iint W(u,v)e^{+2\pi i(ul+vm)}dudv.$$

$${\Omega }_{\mathrm{syn}}=\iint B(l,m)dldm=W(0,0).$$

$$\left[\frac{F}{(\mathrm{Jy}/{\Omega }_{\mathrm{conv}})}\right]=\left[\frac{M}{(\mathrm{Jy}/{\Omega }_{\mathrm{conv}})}\right]+\left[\frac{R}{(\mathrm{Jy}/{\Omega }_{\mathrm{syn}})},\right].$$

$${\left(\frac{1}{\Delta S^i}\right)}^2=\sum _k{\left(\frac{1}{\Delta S_k^v}\right)}^2=N_{\mathrm{vis}}{\left(\frac{1}{\Delta S^v}\right)}^2,$$

$$w_k\equiv {\left(\frac{1}{\Delta S_k^v}\right)}^2=\frac{1}{N_{\mathrm{vis}}}{\left(\frac{1}{\Delta S^i}\right)}^2.$$

$$\Delta S_k^v=C_{ij}\sqrt{\frac{T_{\mathrm{sys},i}{T}_{\mathrm{sys},j}}{B\cdot t_{\mathrm{vis}}}}=C_{ij}\sqrt{\frac{T_{\mathrm{sys}}^2}{B\cdot t_{\mathrm{vis}}}},$$

$$C_{ij}=\frac{2k_B}{\sqrt{({\eta }_{a,i}{A}_i)({\eta }_{a,j}{A}_j)}}\frac{1}{\sqrt{2}{\eta }_q},$$

$$\Delta S^i=C_{ij}\sqrt{\frac{T_{\mathrm{sys}}^2}{B\cdot t_{\mathrm{tot}}}},$$

$$\Delta S_k^v=C_{\mathrm{TP}}\sqrt{\frac{T_{\mathrm{sys}}^2}{B\cdot t_{\mathrm{vis}}}}$$

$$\Delta S^i=C_{\mathrm{TP}}\sqrt{\frac{T_{\mathrm{sys}}^2}{B\cdot t_{\mathrm{tot}}}},$$

$$C_{\mathrm{TP}}=\frac{k_B}{{\eta }_{a}A}\frac{1}{{\eta }_{\mathrm{mb}}}\frac{1}{{\eta }_q},$$

$$W(u,v)$$

$$=$$

$$B(l,m)=\frac{1}{1+\beta }\left[B_{\mathrm{INT}}(l,m)+\beta B_{\mathrm{TP}}(l,m)\right].$$

$${\Omega }_{\mathrm{syn}}=W(0,0)=\frac{\beta }{1+\beta }{W}_{\mathrm{TP}}(0,0),$$

$${\Omega }_{\mathrm{conv}}=\frac{\pi b_{\mathrm{maj}}\times b_{\min }}{4\ln 2}.$$

$$\beta =\frac{{\Omega }_{\mathrm{syn}}}{W_{\mathrm{TP}}(0,0)-{\Omega }_{\mathrm{syn}}}=\frac{{\Omega }_{\mathrm{conv}}}{W_{\mathrm{TP}}(0,0)-{\Omega }_{\mathrm{conv}}}.$$

$$W_{\mathrm{TP}}(0,0)=\frac{\pi b_{\mathrm{TP}}^2}{4\ln 2},$$

$$\begin{array}{rl}W& {}_{\mathrm{INT}}(u,v):\beta W_{\mathrm{TP}}(u,v)\\ & =\frac{{\overline{W}}_{\mathrm{INT}}(u,v)}{\iint {\overline{W}}_{\mathrm{INT}}(u,v)dudv}:\beta W_{\mathrm{TP}}(u,v)\\ & ={\overline{W}}_{\mathrm{INT}}(u,v):\beta \left[\iint {\overline{W}}_{\mathrm{INT}}(u,v)dudv\right]W_{\mathrm{TP}}(u,v).\end{array}$$

$$w_k^{\mathrm{TP}}=\frac{\beta }{N_{\mathrm{vis}}^{\mathrm{TP}}}\sum _k{w}_k^{\mathrm{INT}}.$$

$$\left[\frac{\mathrm{D}}{(\mathrm{Jy}/{\Omega }_{\mathrm{syn}})}\right]=\left[\frac{\mathrm{I}}{(\mathrm{Jy}/{\Omega }_{\mathrm{syn}})}\right]+\left[\frac{\mathrm{R}}{(\mathrm{Jy}/{\Omega }_{\mathrm{syn}})}\right],$$

$$\frac{{\Omega }_{\mathrm{conv}}}{{\Omega }_{\mathrm{syn}}}=\frac{F-R}{D-R}.$$

$$\text{Accuracy}=\frac{(\text{Cleaned image})-(\text{Smoothed model})}{(\text{Smoothed model})}.$$

$$\delta (x)\equiv \frac{\rho (x)-{\rho }_0}{{\rho }_0}.$$

$$\mathbb{x}=D\left(\begin{array}[]{c}l/L\\ m/M\\ n/N\end{array}\right),\mathbb{k}=\frac{1}{D}\left(\begin{array}[]{c}u\\ v\\ w\end{array}\right),$$

$$\delta (l,m,n)$$

$$\tilde{\delta }(u,v,w)$$

$$\tilde{\delta }(u,v,w)={\tilde{\delta }}^{\ast }(L-u,M-v,N-w).$$

$$P(k)=\frac{C}{(k^2+k_{\min }^2{)}^{n/2}},$$

$${\sigma }^2=\sum _{k_u,k_v,k_w}P(k)[\Delta k_u\Delta k_v\Delta k_w]=C\sum _{k_u,k_v,k_w}\frac{1}{(k^2+k_{\min }^2{)}^{n/2}}[\Delta k_u\Delta k_v\Delta k_w].$$

$${\sigma }^2$$

Full text

††slugcomment: Accepted to PASP

1 Introduction

The Atacama Large Millimeter/submillimeter Array (ALMA) consists of three sets of telescopes: two interferometers with 12m and 7m dishes and Total Power (TP) 12m single-dish telescopes. The interferometers alone cannot obtain information on an extended emission distribution. Combined, these three sets can recover structures from small to large scales. This emphasizes the importance of single-dish data in the imaging process. Interferometers collect Fourier components of a 2-dimensional brightness distribution in the sky, but cannot obtain full coverage of the $uv$ space (i.e., Fourier space). Hence, the empty part of the $uv$ space has to be filled (or “guessed”) in the deconvolution process, with the filled part of the $uv$ space used as a guideline. Convergence is not guaranteed, but a better sampled $uv$ space helps this process. The combination of interferometer and single-dish data provides better $uv$ coverage, thus leading to better imaging results.

The idea of combining interferometry and single-dish data was introduced early (Bajaja & van Albada, 1979; Ekers & Rots, 1979, who called it “short-spacing corrections”), but in practice the procedure is still under debate. The history and techniques for combination are summarized in Stanimirovic (2002). The techniques can be classified into three general categories: whether the data are combined (1) before, (2) during, or (3) after deconvolution (e.g., Clean). For better $uv$ coverage for deconvolution, (1) and (2) have natural advantages (e.g., Cornwell, 1988). Each of these categories can be split into two subcategories: whether the operations are performed (a) in the Fourier domain ($uv$ domain) or (b) in the image domain.

1.1 Combination After Deconvolution

A common technique of the after-deconvolution combination (category 3a) applies low- and high-pass filters to single-dish and cleaned interferometer maps (or images), respectively. The maps are then added in the Fourier domain and transformed back to a final map. It is now called the Feather technique. Several implementations exist (e.g., the Imerge task in Aips; Immerge task in Miriad; Herbstmeier et al., 1996; Weiß et al., 2001). The Feather task111https://casa.nrao.edu/casadocs/latest/global-task-list/task_feather/about in Casa is one more recent example. It adopts the low-pass filter, $L(u,v)=FT\{B_{\rm TP}(l,m)\}$, and a high-pass filter, $H(u,v)=1-L(u,v)$, where the Fourier transform (FT) of the TP primary beam $B_{\rm TP}(l,m)$ is used 222Feather first deconvolves a TP map with $B_{\rm TP}$, and then multiplies $L(u,v)=FT\{B_{\rm TP}(l,m)\}$ in the Fourier domain. Therefore, it recovers the original TP map. In other words, Feather simply adds the original TP map.. Feather’s choice of $H(u,v)$ extends to $uv\approx 0$ which interferometer data does not cover. To avoid this, Blagrave et al. (2017) used different filters with non-zero values only where the data exist. Since multiplications (or divisions) in the Fourier domain are equivalent to convolutions (deconvolutions) in the image domain, it is possible to implement the same technique both in the Fourier and image domains. Faridani et al. (2018) implemented the Feather technique in the image domain (category 3b).

1.2 Combination During Deconvolution

The TP map can be used as an initial model for deconvolution of interferometer data (category 2a). Practically, such a scheme deconvolves single-dish and interferometer data jointly (Cornwell, 1988). The Clean/Tclean tasks in Casa can take a single-dish map as an input. Dirienzo et al. (2015) employed this technique and produced a combined map. They further applied Feather to ensure that the large-scale components precisely represented the TP data.

1.3 Combination Before Deconvolution

Pseudo visibilities can be generated from a single-dish map. They can be added to interferometer visibilities before deconvolution (category 1a; Vogel et al., 1984). Once TP visibilities are generated from the TP map, they are readily fed into existing deconvolvers (e.g., Clean/Tclean). This technique has also been used, though there are some discussions on how to optimize the distribution and weights of the pseudo visibilities (e.g., Rodríguez-Fernández, 2008; Kurono et al., 2009; Pety & Rodríguez-Fernández, 2010; Koda et al., 2011). For example, Koda et al. (2011) generated a Gaussian visibility distribution and set the visibility weight so that the single-dish primary beam is represented. One of the advantages of pseudo visibilities techniques is that the visibility weight can be easily controlled within single-dish data, as well as with the interferometer data.

The same technique can be implemented in the image domain (category 1b; Stanimirovic et al., 1999). The single-dish and interferometer dirty beams, as well as their maps, are added with a choice of $L(u,v)$ and $H(u,v)$. The combined dirty map can then be deconvolved with the combined dirty beam. This is equivalent to the pseudo visibility techniques, since a FT of a sum of two functions is equal to the sum of the FT of the two functions. The weights can be changed by adjusting $L(u,v)$ and $H(u,v)$. Such implementation is being developed in Casa (U. Rau in preperation).

1.4 This Work

Tools are needed for the combination of ALMA data in the platform of the Common Astronomy Software Applications (Casa; McMullin et al., 2007). Here, we introduce a new package in Casa that converts ALMA’s TP map into pseudo visibilities (dubbed Tp2vis). It generates visibilities as if they were observed by a virtual interferometer with short spacings down to zero spacial frequency. Together with ALMA12m+7m data, the Tp2vis visibilities can be fed into standard deconvolvers for joint deconvolution. The technique follows the one developed by Koda et al. (2011) on the Miriad platform (Sault et al., 1995), but with improvements.

The Tp2vis package provides four functions in the CASA platform. The key function, Tp2vis, generates a CASA measurement set (MS) –visibilities with weights– from a TP map. The weights are calculated based on the root-mean-square (RMS) noise from the TP map. This MS is ready for joint deconvolution, though in some cases one may wish to manipulate the weights using the supplementary function, Tp2viswt. An accessory function, tp2vispl, plots the weights of ALMA 12m, 7m, and TP visibilities for comparison. Another function, tp2vistweak, attempts to fix the problem of beam-size mismatch in the image space after deconvolution (Section 3.2).

2 The TP2VIS procedure

The Tp2vis function takes a sky brightness distribution (TP map) and simulates observations with a virtual interferometer (defined below). Figure 1 shows the flowchart of steps in the tp2vis function. The details of the steps are in the following subsections. From a top view, the procedure can be separated into three parts. First, it converts the TP map into the sky brightness distribution, which will be observed by the virtual interferometer (steps A and B in the figure). Second, it converts the brightness distribution into visibilities (steps C and D). And third, the weights of the TP visibilities are set so that they represent the RMS noise of the original TP map (step E).

Different types of beams are involved in our discussions. For clarity we list the beams in Table 1.

2.1 Step A: Deconvolution with TP Beam $\Omega_{\rm TP}$

This step is a preparation for the virtual interferometer observations. When the TP telescopes observe the sky, the TP telescope beam, $\Omega_{\rm TP}$, is convolved with the brightness distribution in the sky. Therefore, the first step is to deconvolve the TP map with $\Omega_{\rm TP}$ and to obtain the true sky brightness distribution in Jy/beam.

We approximate $\Omega_{\rm TP}$ with a Gaussian with a full width half maximum (FWHM) of $56.6\arcsec$ at 115.2 GHz. We assume that it scales linearly with 1/frequency.

2.2 Step B: Apply Virtual Primary Beam $\Omega_{\rm vir}$

When an interferometer observes the sky, the brightness distribution is attenuated by its primary beam (PB). Hence, the sky brightness distribution from Step A is multiplied by the PB pattern of the virtual interferometer ($\Omega_{\rm vir}$). In case of mosaic observations, this must be done at each pointing position.

In our virtual observations, the coordinates of pointing centers can be set arbitrarily within the TP map. We typically take them from the ALMA 12m data, so that the 12m and TP visibilities have the same sky coverage. The pointing coordinates are listed in an ASCII file and passed to the Tp2vis function with the “ptg=” argument.

2.3 Step C: Generate Gaussian Visibility Distribution

In parallel with steps A and B, a visibility distribution for the virtual observations (hence a MS) is prepared. We generate a two-dimensional Gaussian distribution of visibilities in the $uv$ space, so that its Fourier transformation, i.e., the synthesized beam $\Omega_{\rm syn}$, reproduces the Gaussian beam of the TP telescopes ($\Omega_{\rm TP}$). Adopting a constant weight for all TP visibilities (Section 4), this distribution would represent the sensitivity distribution within the TP beam. The width of the Gaussian visibility distribution is determined from those of the Gaussian beam. The standard deviations of the Gaussian beam in the sky $\sigma$ and the Gaussian distribution in the $uv$ space $\sigma_{\rm F}$ are related as $\sigma_{\rm F}=1/(2\pi\sigma)$.

The full width at half maximum (FWHM) is calculated as FWHM=$2\sqrt{2\ln 2}\sigma$ ($\approx 2.355\sigma$). For a FWHM of $56.6\arcsec$, $\sigma=2.74\times 10^{-4}\,\rm radian$ and $\sigma_{\rm F}=580\lambda$, where $\lambda$ is the wavelength.

2.4 Step D: Fill Visibility Amplitudes and Phases

We already have the sky brightness distribution with the PB of the virtual interferometer applied (from step B), and the Gaussian visibility distribution in the $uv$ space (step C). Now, the sky brightness distribution is Fourier transformed into the $uv$ space. The amplitude and phase of the transformation are then sampled at the visibility positions, which fill the AMPLITUDE and PHASE columns of the visibilities in the MS.

These visibilities are already internally consistent. The Clean/Tclean tasks with natural weighting produce a synthesized beam $\Omega_{\rm syn}$ roughly equal to the TP telescope beam $\Omega_{\rm TP}$. The brightness distribution also reproduces the observed TP brightness distribution.

2.5 Step E: Set Weights of TP Visibilities

The Tp2vis function sets visibility weights according to the RMS noise in the TP map. The RMS value is set manually by the “rms=” argument [in units of Jy/beam of the TP map]. [Note that in the current implementation, the RMS is set by hand, because it is not always trivial to find emission-free channels using an algorithm.] The RMS is converted to the weight of individual visibilities. The equation for this conversion is in Section 4. The sensitivity-based weight is a proper representation of the quality of the TP data and is the default of Tp2vis.

In the history of radio interferometry, the weights are often manipulated. The most common examples are the uniform, robust/briggs, or taper weighting schemes. There is no reason not to manipulate them in some other ways – for example, to match the weights of TP visibilities with those of ALMA 12m and 7m data (their sensitivities may not match when one works on archival data). The supplementary tp2viswt function provides several options to manipulate the weights. An accessory tp2vispl function plots the weight distributions of ALMA 12m, 7m, and TP visibilities for comparison, so that users can decide how the TP visibilities need to be weighted with respect to 12m and 7m data.

3 Beam Disparity and Flux Problem

Imaging of interferometer data involves several different beams. Some of them are often referred with different names by different researchers or in different circumstances. This is very confusing. Furthermore, their differences cause a more practical issue in flux conservation in the process of imaging. In order to discuss the issue, we first explain the beams. For Tp2vis we use four different beams (Table 1).

The first two beams were used in Section 2: the primary beam of the TP telescopes ($\Omega_{\rm TP}$) and that of the virtual interferometer ($\Omega_{\rm vir}$). The third and fourth beams are from visibilities, i.e., a synthesized beam ($\Omega_{\rm syn}$) and convolution beam ($\Omega_{\rm conv}$). The $\Omega_{\rm syn}$ is also called the dirty beam ($\Omega_{\rm dirty}$), deconvolution beam ($\Omega_{\rm dec}$), or point spread function ($\Omega_{\rm psf}$). The $\Omega_{\rm conv}$ is also called the clean beam ($\Omega_{\rm clean}$) or restoring beam ($\Omega_{\rm res}$). For discussions in later sections, we give more precise definitions of $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$ in Section 3.1.

Jörsater & van Moorsel (1995) extensively discussed the problem of flux conservation when there is a disparity between the areas of $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$. The problem becomes more apparent when the TP visibilities are included (Section 3.2). We discuss two techniques to resolve/mitigate this problem: by adjusting the weights of the TP visibilities (Section 4.2) and/or by scaling the residual map (Section 5).

Note that we use the symbols, $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$, to represent the synthesized and convolution beams, respectively, but also to denote their areas.

3.1 The Synthesized and Convolution Beams

The synthesized beam $\Omega_{\rm syn}$ is determined by the distributions of visibility weights in $uv$ space [Note: $\Omega_{\rm syn}=\Omega_{\rm dirty}=\Omega_{\rm dec}$]. Its 2-dimensional pattern, $B(l,m)$, is a Fourier transformation of the weight distribution $W(u,v)$. With the normalization of $\iint W(u,v)dudv=1$,

[TABLE]

The area of $\Omega_{\rm syn}$ is the beam solid angle and is given as an integration over space,

[TABLE]

It depends only on $W(0,0)$, the weight at ($u$,$v$)=(0,0). Therefore, only the weight of the TP visibility matters, because the interferometer data have zero weight at ($u$,$v$)=(0,0). Therefore, $\Omega_{\rm syn}=W_{\rm TP}(0,0)$ when ALMA 12m, 7m, and TP are combined.

The convolution beam $\Omega_{\rm conv}$ is typically defined by a Gaussian fit to the central peak in the $B(l,m)$ pattern [Note: $\Omega_{\rm conv}=\Omega_{\rm clean}=\Omega_{\rm res}$]. Hence, it depends on the full weight distribution, $W(u,v)$; but in practice, the outer, long-baseline part of the visibility weight distribution is often the determining factor.

In summary, $\Omega_{\rm syn}$ is set at the center of the $uv$ space, while $\Omega_{\rm conv}$ is determined mainly by the outer part.

3.2 Beam Disparity and Flux Problem

For simplicity, we use Clean as an example here, but it applies to other deconvolvers as well. The Clean task iteratively finds peaks in a dirty map, whose brightness units are Jy/$\Omega_{\rm syn}$. In each iteration Clean identifies a peak in a dirty map, fits the 2-dimensional pattern of $\Omega_{\rm syn}$ to the peak and obtains its position and amplitude, subtracts it from the dirty map, and puts a $\delta$-function at that position in a clean-component map. The Clean task then convolves the clean-component map with $\Omega_{\rm conv}$ to produce a model map. The model map has brightness units of Jy/$\Omega_{\rm conv}$. Typically, Clean is not perfect and leaves emission in a residual map. The unit of the residual map remains Jy/$\Omega_{\rm syn}$. Clean naively adds the model and residual maps to produce a final map, whose unit is often denoted as Jy/beam. For the “beam” in the denominator, $\Omega_{\rm conv}$ is typically adopted. Therefore,

[TABLE]

where $F$, $M$, and $R$ stand for the final, model, and residual maps, respectively.

The areas and patterns of $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$ could be very different, as they depend on different parts of the weight distribution in the $uv$ space (Section 3.1). In particular, if their areas are not equal, $\Omega_{\rm syn}\neq\Omega_{\rm conv}$, the brightness of $R$ in eq. (3) is calculated incorrectly. For example, if $\Omega_{\rm syn}\gg\Omega_{\rm conv}$, the brightness in $R$ is overestimated by $\Omega_{\rm syn}/\Omega_{\rm conv}$ ($\gg 1$), leading to an overestimation of flux in $F$. This happens when the ALMA 12m, 7m, and TP visibilities are combined. In theory, the $R$ could be zero if we could clean down to zero flux, but in practice, it is nearly impossible to get there, as the map to be cleaned has noise.

In case of $\Omega_{\rm syn}\ll\Omega_{\rm conv}$, the brightness in $R$ is underestimated. The impact of this on $F$ may not be as noticeable as in the opposite case, since $R$ typically has much less flux than $M$. This happens when the 12m and 7m visibilities are imaged without TP, since $\Omega_{\rm syn}=W_{\rm TP}(0,0)=0$ (hence, $\ll\Omega_{\rm conv}$). This flux problem did not attract much attention in the past, as single-dish data were typically not included.

4 Weights

The optimal weighting scheme for TP visibilities is a topic of debate. Tp2vis provides some implementations of weighting schemes, as well as providing a framework for testing other schemes in the future. It has been very common to manipulate visibility weights in the history of radio interferometry. For example, the natural weighting scheme takes the sensitivities of individual visibilities measured during observations, and uses them as a weight distribution in the $uv$ space. The uniform weighting scheme takes a binary approach, setting the weights to 1 at the positions in the $uv$ space where visibilities exist and 0 for the rest, independent of sensitivities. The robust/briggs weighting connects the natural and uniform weightings and transforms one from the other seamlessly by adjusting the single “robust” parameter. Of course, one could come up with other arbitrary weighting schemes to accomplish particular scientific goals. No matter what weights are used, the final map should be usable for scientific applications, as long as the same weights are used consistently for both map and beam by deconvolvers.

Our main function, Tp2vis, sets the weights so that their sum represents the RMS noise of the original TP map. This is the default in Tp2vis. In case that users wish to manipulate the weights in other ways, the supplementary function, Tp2viswt, provides options to manipulate the weights. Currently, it has five modes:

• •

mode=’statistics’: output statistics of the current weights in MS(s).

• •

mode=’constant’, value=: set the weights of all visibilities to a constant specified by “value=”.

• •

mode=’multiply’, value=: multiply the current weights by a constant specified by ”value=”.

• •

mode=’rms’, value=: use the sensitivity-based weight for all visibilities (Section 4.1). An RMS value from TP map should be set with “value=”. This is the default weights in Tp2vis.

• •

mode=’beam’: use the beam size-based weight for all the TP visibilities (Section 4.2).

In what follows, we discuss two provisional approaches to optimizing the weight, one based on the RMS noise (Section 4.1) and another that matches the areas of synthesized and convolution beams (Section 4.2). This beam matching is important to ensure flux conservation in deconvolution, and we will revisit this in Section 5. Tp2vis adjusts the weight in the $uv$ space, but the same can be done in the image domain (Section 4.3).

4.1 The Sensitivity-Based Weight

One of the most intuitive weighting schemes is the one that represents the quality of TP data, namely the root-mean-square (RMS) noise of the TP map. This is the default weighting scheme used by the Tp2vis function, and the Tp2viswt function with the mode=’rms’ can set the weights to the sensitivity-based one. This section explains the calculation of the sensitivity-based weights in Tp2vis. Section 4.1.1 discusses the conversion of the RMS noise of a TP map into the weights of individual visibilities. For a joint deconvolution, the weights of 12m, 7m, and TP data should be set consistently. Hence, section 4.1.2 shows the corresponding weights of 12m and 7m data. [Note that CASA version 5 uses this definition of 12m and 7m weights by default, improving upon prior CASA versions.]

This sensitivity-based weighting is equivalent to the conventional “natural” weighting. For interferometer data, additional manipulations, such as “robust/briggs” or “uniform” weighting schemes, are often applied on top of the sensitivity-based weights. Such manipulations can also be applied to the sensitivity-based weights from Tp2vis. For simplicity, however, the discussions in this section will focus on the sensitivity-based (natural) weights.

4.1.1 Weights for TP Visibilities

With the natural weighting, the sensitivity of a map, or the RMS noise of the image ($\Delta S^{i}$), is a simple summation over all visibility sensitivities ($\Delta S^{v}_{k}$),

[TABLE]

where $N_{\rm vis}$ is the number of visibilities. All visibilities are assumed to have the same sensitivity ($\Delta S^{v}_{k}=\Delta S^{v}$), and a single RMS value is assumed for the map. Hence, the weight of each TP visibility is

[TABLE]

4.1.2 Corresponding Weights for 12m & 7m Visibilities

For consistency, the 12m and 7m data should also use the sensitivity-based weights. The sensitivity for a single visibility $k$ between antennas $i$ and $j$ is given by

[TABLE]

where $T_{{\rm sys},i}$, $T_{{\rm sys},j}$, and $B$, $t_{\rm vis}$ are the system temperatures of $i$ and $j$, bandwidth, and integration time of the visibility, respectively (Taylor et al., 1999; Thompson et al., 2007). For simplicity, we assume $T_{{\rm sys},i}=T_{{\rm sys},j}=T_{\rm sys}$. The coefficient $C_{ij}$ is

[TABLE]

where $\eta_{a,i}$, $\eta_{a,j}$, $A_{i}$, $A_{j}$, and $\eta_{q}$ are the aperture efficiencies and areas of antennas $i$ and $j$, and the quantum efficiency of correlator, respectively. With the natural weighting, the RMS sensitivity of final map is

[TABLE]

where $t_{\rm tot}$ is the total integration time $t_{\rm tot}=N_{\rm vis}t_{\rm vis}$. The weights of 12m and 7m visibilities should be set to $w_{k}=(1/\Delta S^{v}_{k})^{2}$ with eq. (6).

4.1.3 For Internal Consistency

This subsection is not necessary in practice, but is included for consistency among all data columns in the MS. It might become important in future versions of CASA. Tp2vis fills the WEIGHT and SIGMA columns. Only the WEIGHT column is used by Clean/Tclean in Casa – currently, other columns do not matter. However, one could also make the EXPOSURE column (i.e., integration time of each visibility) consistent.

The sensitivity equations for TP, corresponding to eqs. (6), (8), and (7), are

[TABLE]

and

[TABLE]

where

[TABLE]

and $\eta_{\rm mb}$ is the main beam efficiency of the TP antennas. $T_{\rm sys}$, $\Delta S^{i}$, and $B$ are determined by observation. With eq. (10) we can derive $t_{\rm tot}$, using $T_{\rm sys}$ and $\Delta S^{i}$ from MS and TP maps, respectively. Equations (5), (9), and (10) give $t_{\rm vis}=t_{\rm tot}/N_{\rm vis}$.

4.2 The Beam Size-based Weight

Section 3.2 discussed the problem in flux due to the disparity between the areas of $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$. The beam size-based weighting scheme discussed here is an attempt to adjust the weights of TP visibilities to equalize their areas. Note again that we use the symbols, $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$, to represent the areas of the beams as well as to refer the beams themselves.

4.2.1 The Synthesized and Convolution Beams from Combined INT+TP Visibilities

The pattern of the synthesized beam $\Omega_{\rm syn}$, denoted by $B(l,m)$, is given by eq. (1). The weight is normalized as $\iint W(u,v)dudv=1$. For interferometers (hereafter INT; e.g., ALMA 12m and 7m) and TP, the synthesized beam patterns, $B_{\rm INT}(u,v)$ and $B_{\rm TP}(u,v)$, are derived with their normalized weights, $W_{\rm INT}(u,v)$ and $W_{\rm TP}(u,v)$, respectively.

By defining a relative weight parameter $\beta$, a combined weight distribution of INT and TP can be expressed as

[TABLE]

The corresponding synthesized beam pattern is

[TABLE]

The area of the synthesized beam is derived from eqs. (2)(13),

[TABLE]

where we used $W_{\rm INT}(0,0)=0$ since INT has no contribution at zero spacing. Note that $W_{\rm TP}(0,0)$ is equal to the area of the synthesized beam of TP visibilities.

The convolution beam $\Omega_{\rm conv}$ is typically derived by a 2-d Gaussian fit to the central peak in $B(l,m)$. Using the major and minor axis diameters, $b_{\rm maj}$ and $b_{\rm min}$, from the fit, the area of $\Omega_{\rm conv}$ is

[TABLE]

4.2.2 Setting Up Weights for TP Visibilities

This weighting scheme attempts to achieve $\Omega_{\rm syn}=\Omega_{\rm conv}$. Equation (15) gives $\beta$ as

[TABLE]

$W_{\rm TP}(0,0)$ is equal to the area of the synthesized beam of TP visibilities alone, which, by construction, is the same as that of the Gaussian primary beam of the TP antennas. Hence,

[TABLE]

using the FWHM beam size of the TP antennas, $b_{\rm TP}$. The $\Omega_{\rm conv}$ should be calculated from INT+TP visibilities, but in practice it depends primarily –and almost solely– on the longest baselines. Hence, for simplicity, we derive $\Omega_{\rm conv}$ from INT visibilities alone. We can calculate $\beta$ from eqs. (17), (18), and (16) for INT.

Using $\beta$, we adjust the TP visibility weight distribution from $W_{\rm TP}$ to $\beta W_{\rm TP}$, by scaling the weights of individual visibilities ($w^{\rm INT}_{k}$ and $w^{\rm TP}_{k}$ for INT and TP visibilities in MSs, respectively).

For the following, we use the notation $\overline{W}$ for an unnormalized weight distribution. The ratio of the two terms in the numerator of eq. (13) is

[TABLE]

This transformation converts the normalized $W_{\rm INT}(u,v)$ into the unnormalized $\overline{W}_{\rm INT}(u,v)$, which is what is delivered from the ALMA observatory. The term to the right of the ratio mark gives the corresponding normalization for the weights of TP visibilities. When the natural weighting is used, the unnormalized $\overline{W}_{\rm TP}(u,v)$ and $w^{\rm TP}_{k}$ are related as $\iint\overline{W}_{\rm TP}(u,v)dudv=\sum_{k}w^{\rm TP}_{k}=N_{\rm vis}^{\rm TP}w^{\rm TP}_{k}$. Thus, for the normalized $W_{\rm TP}(u,v)$, we set $w^{\rm TP}_{k}=1/N^{\rm TP}_{\rm vis}$. Also, the term in the bracket is simply, $\iint\overline{W}_{\rm INT}(u,v)dudv=\sum_{k}w^{\rm INT}_{k}$. Therefore,

[TABLE]

This weight would satisfy $\Omega_{\rm syn}=\Omega_{\rm conv}$ when the INT and TP visibilities are combined.

4.3 Notes on Visibility and Weight Distributions

The weight distribution of TP visibilities can be controlled, either by changing the distribution of TP visibilities in $uv$ space, and/or by adjusting the weights of individual visibilities. Tp2vis takes the former approach: generating a Gaussian visibility distribution with a constant weight for all the visibilities.

The latter approach could be useful in future. For example, one could generates a grid distribution, instead of the Gaussian distribution, and change the visibility weights as a function of $uv$ distance so that the weight distribution follows a Gaussian. This way, one could avoid shot noise due to the randomly-distributed discrete visibilities. In the end, the Clean task would put the visibilities onto a grid. If the grid spacing is controlled consistently between Tp2vis and Clean, we can test grid-based approaches in future. Of course, the weight distribution can be realized on a grid without visibilities. An implementation of this would require an update of the Clean/Tclean task in Casa (U. Rau 2018, in private communication).

5 The Residual Scaling

The flux problem in Section 3.2 is due to the disparity between the areas of $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$ in eq. (3). The residual part of the final image from Clean is added incorrectly with inconsistent units of brightness. To circumvent this inconsistency, we could rescale the residual map by a factor of $\Omega_{\rm conv}/\Omega_{\rm syn}$ and add it to the model map to generate a new final map. The tp2vistweak function does this operation 333One should keep it in mind that the patterns of $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$ are still not the same, even though their areas are the same..

In theory, $\Omega_{\rm conv}/\Omega_{\rm syn}$ can be calculated with the images of $\Omega_{\rm syn}$ from Clean, but it turns out not to be that simple. The estimation of $\Omega_{\rm syn}=W(0,0)$ depends on how the visibilities are gridded in the imaging process. The integration of the synthesized beam, $\Omega_{\rm syn}=\int B(l,m)dldm$, requires an unrealistically large image size as its tiny responses at the far outskirts add up and their area is cumulatively large.

Instead, Jörsater & van Moorsel (1995) suggested a more practical approach by taking advantage of detected emissions. In parallel to eq. (3), we can write a similar equation for the dirty map $D$,

[TABLE]

where $I$ is a clean-component map, i.e., a map of the emissions that will be identified by Clean. We should note that the brightness unit is Jy/$\Omega_{\rm syn}$ for all terms, and that $M$ and $I$ in eqs. (3) and (21) represent exactly the same emissions in different units. The $D$, $F$, and $R$ are given by Clean. We can calculate the scaling factor as

[TABLE]

We multiply the $R$ map by this factor and add it to $M$ to construct a new $F$ as in eq. (3), but with a consistent brightness unit (consistent beam areas). This residual scaling method is included in the Aips task Imager and was used recently, e.g., by Walter et al. (2008) and Ianjamasimanana et al. (2017).

In principle, this method can be applied on a pixel-by-pixel basis, as long as the emissions are significantly detected across the map. In practice, this is not always the case. In addition, different spillovers of the beam patterns of $\Omega_{\rm syn}$ and $\Omega_{\rm conv}$ degrade the accuracy of the calculation on a pixel-by-pixel basis. Hence, Tp2vis integrates the part of map with significant detections and calculates a single beam ratio with eq. (22). Currently, Tp2vis uses a single ratio even for mosaics, which may be modified in the future as the synthesized beam may vary among mosaic fields.

One may wonder why Jörsater & van Moorsel (1995) could apply this technique to pure interferometer data without single-dish. In that case $\Omega_{\rm syn}=W(0,0)=0$, and the denominator of the left-hand side of eq. (22) is apparently zero. The answer lies in the fact that eq. (22) is the ratio of surface brightnesses in units of Jy/beam (see the right-hand side). Both the numerator (i.e., flux) and denominator (beam) of $D-R$ are zero, but their ratio, and hence the surface brightness, has a finite value. Hence, this method is applicable to pure interferometer data as well.

6 Imaging Tests

We test Tp2vis using model images. With their known emission distributions, the model images allow us to evaluate the accuracy of image deconvolution with Tp2vis. We use Casa version 5.4.0-68 and employ the Clean task for imaging. Examples with actual ALMA data are found in the Github repository of Tp2vis, but not in this paper, as they will evolve with future updates of Casa.

6.1 Setup

6.1.1 Model images

Two model images are used for tests (Figure 2). One is a single Gaussian emission distribution with a standard deviation of 1500 pixels in width (i.e., FWHM of $\sim 53\arcsec$ in the definition described below), and another is a pseudo sky image that mimics small to large structures within giant molecular clouds (GMCs), i.e., dense cores and extended emission (see Appendix A). The current Clean/Tclean tasks tend to cause flux divergence near the edge of a mosaic field of view. As a workaround, the emission around the edge is tapered off outward so that it reaches zero at around the $\sim 80\%$ level of the maximum primary beam coverage (see Section 6.1.4).

Both images have the same dimensions of $16,384^{2}$ and are configured to the same setup, i.e., an image center coordinate of (RA, DEC) = (12h00m00s, -35d00m00s), a cell size of $0.015\arcsec$ (hence the image size of $4.096\arcmin$), and a frequency of 115 GHz with a bandwidth of 2 GHz (a single channel). We run the Simobserve task to convert the images into the CASA format with these parameters. The $0.015\arcsec$ cell size is chosen to include the main structures in a 38-pointing mosaic (Section 6.1.2) and corresponds to baselines of $\sim 36$ km when observed at 115 GHz. Therefore, the images can be properly sampled by all the ALMA configurations including the one with the longest baseline ($\sim$16 km). The flux scale is arbitrary, and only the relative flux with respect to the peak is important in the tests. The peak flux is set to 1 for Figure 2.

6.1.2 Visibilities for 12m and 7m Arrays

We employ a 38-pointing mosaic pattern with a $0.4\arcmin$ spacing. This is a slightly finer spacing than the Nyquist sampling of the primary beam of the 12m array. The coordinates of the 38 mosaic pointings are determined by Simobserve. The same coordinates are used for the 12m, 7m, and TP arrays, and we explicitly set the $uv$ coverage to be exactly the same for all the 38 pointings so that the point spread function (PSF) shape does not vary across the mosaic field except near the edge (see Section 6.1.4).

Figure 2 shows the 38 pointing centers (crosses), the 40 and 80% levels of the 12m+7m-combined primary beam attenuation (green lines), and the 12m and 7m beam sizes (white circles) as references. To control the number of visibilities, we adjust total and individual integration times in Simobserve. The weights of visibilities are set separately, so these times affect nothing but the number of visibilities to be generated: the total integration time of each pointing per array configuration is set to 1,000 sec with a 10 sec integration. This results in 100 visibilities per baseline. We use one configuration for the 7m array (”aca.i.cfg”) and four configurations for the 12m array (”alma.cycle4.4.cfg”, ”alma.cycle4.3.cfg”, ”alma.cycle4.2.cfg”, and ”alma.cycle4.1.cfg”). The parameters of these array configurations are described in detail in the ALMA Proposer’s Guide of Cycles 4, 5, and 6.444https://almascience.nrao.edu/proposing/proposers-guide The tables of the antenna positions of these configurations are distributed with CASA.

6.1.3 Visibilities for TP Array

The TP array observes the sky as single-dish telescopes with a FWHM beam size of $56.7\arcsec$ at 115 GHz. TP data are delivered as a calibrated image, or cube, with a cell size of $5.6\arcsec$. Hence, we smooth the model images with a Gaussian with a FWHM of $56.7\arcsec$ and regrid them on a $5.6\arcsec$ cell size. We adopt the the same pointing coordinates as those of the 12m and 7m visibilities. Tp2vis is run on these images and generates 5,175 TP visibilities per pointing.

The model images do not have noise, and the visibility weights must be set arbitrarily. We apply the sensitivity-based weight (Section 4.1) with relative image RMS sensitivities of 1 for each of the four 12m configurations, 4.5 for the 7m configuration, and 7 for the TP visibilities. These relative sensitivities are chosen based on typical ALMA mosaic data we have at hand (e.g., of nearby galaxies, Galactic molecular clouds), Table A-2 in the ALMA Proposer’s Guide555https://almascience.nrao.edu/proposing/proposers-guide#section-57, and the ALMA sensitivity calculator. The visibility and sensitivity distributions of the TP, 7m, and 12m arrays are plotted in Figure 3.

6.1.4 Caveats

Extra efforts are needed to circumvent some of the current problems in Clean/Tclean in Casa. In general interferometer observations, the $uv$ coverages are different among mosaic pointings, and the PSF varies across the mosaic. The current version of Clean/Tclean (as of Casa 5.4) does not take this variation into account and uses a single PSF for all mosaic pointings in its minor clean cycles. This simplistic PSF deviates from the true PSF that reflects the $uv$ coverage at each position. The deviation is larger for a larger mosaic, which often leads to an artificial flux divergence. To mitigate this problem, we explicitly set the $uv$ coverage exactly the same for all the 38 pointings so that the PSF stays the same over the most part of the entire mosaic. This unique setup should not be a particular advantage for Tp2vis, since only the local PSF at the position of each emission affects clean results – the local PSF should be consistent with the local $uv$ coverage at that position, but it does not matter whether the PSF is the same at other parts of the mosaic.

However, we note that even with the unique $uv$ coverage for all pointings, the edge of a mosaic suffers from the PSF variation since it could be covered significantly by one array (e.g., the 7m, TP arrays), but not by the others (e.g., the 12m array), due to their different primary beam sizes. The PSF shape should be modified accordingly in such regions, as only a subset of the arrays contributes to the $uv$ coverage.666Such a spatially-variable PSF was implemented in the previous generation software Miriad, which has been successfully applied to ALMA mosaic data (Hirota et al., 2018; Sawada et al., 2018). There is a tendency for current Clean/Tclean results to show a flux divergence when significant emission exists near the edge. The cause has not been entirely isolated, but the simplistic PSF likely plays a part, since Clean converges when each of the 12m, 7m, and TP data are imaged separately. Therefore, for our tests, we taper off the emission near the edge of the model images. For the Gaussian model we subtract a small constant flux from the whole image so that the flux decreases to zero at around the 80% level of the peak of the 12m+7m+TP combined primary beam. The flux outside is set to zero (Figure 2a). For the pseudo sky image we first make a combined primary beam pattern image, $P$, with Clean. We then multiply the original image by “$(P^{2}-0.5)$”, so that the flux becomes zero at around the $\sim 70-80\%$ level of the peak of the 12m+7m+TP combined primary beam. The flux outside is set to zero (Figure 2b).

The Tclean task also suffers from the same problem and appears to have another problem in calculation of the PSF. These issues in Tclean are far beyond the scope of this paper, but their effects and magnitudes should be characterized in the future. Here we adopt Clean as it currently appears somewhat more stable than Tclean for mosaic imaging.

6.2 Deconvolution and Results

We use the Clean task for deconvolution. The termination of Clean iterations is typically controlled by a threshold from the RMS noise, or by the number of iterations. The model images do not have noise, and hence we use the number of iterations for termination and set it to $2\times 10^{5}$. The gain is set to 0.05. The region for clean is restricted to the area above $80\%$ of the peak of the 12m+7m+TP combined primary beam. We use the Briggs weighting with a “robust” parameter of +0.5. It results in a PSF of $1.59\arcsec\times 1.46\arcsec$ with a position angle of $87.8\deg$ for both test images.

Figures 4 and 5 compare the model and cleaned images: (a) and (c) the original model maps smoothed with the PSF, (b) and (d) the cleaned maps with 12m, 7m, and TP visibilities, (e) the residual maps after the residual rescaling, and (f) the accuracy maps. Only relative fluxes are important for comparisons, so the fluxes of all images are normalized by a single constant so that the peak of panel (a) becomes equal to unity. Figures 4 and 5 are normalized separately. The accuracy maps show deviations from the input models and are calculated as

[TABLE]

With this definition an accuracy map plus 1 makes a recovered flux ratio map at the $1.59\arcsec\times 1.46\arcsec$ resolution. The clean is applied over nearly three orders of magnitude in flux as the residuals after clean are $\lesssim 0.3\%$ of the peak flux (panel e). Over this wide dynamic range, the original flux is reproduced with errors mostly less than $\sim 5\%$ (panel f), though some very diffuse regions show errors as high as $\gtrsim$20%.

The purpose of this paper is to introduce the new Tp2vis package. Hence, these tests are performed with simple noise-free images. Further imaging simulations with more realistic setups would help evaluating systematic errors in actual ALMA observations.

7 Conclusions

We presented the Casa-based package “Total Power Map to Visibilities (Tp2vis)”, which converts a TP map into visibilities. The TP visibilities, along with 12m and 7m visibilities, can be jointly-deconvolved with standard deconvolvers (e.g., Clean). The Tp2vis package includes four functions: Tp2vis, Tp2viswt, Tp2vispl, and Tp2vistweak. The Tp2vis function generates a Gaussian visibility distribution, whose inversion, i.e., the synthesized beam, reproduces the primary beam of TP telescopes with natural weighting. The sensitivity-based weighing scheme is adopted as a default. The Tp2viswt function can manipulate the weights. The Tp2vispl function visualizes the TP weights with respect to those of 12m and 7m array data. After the joint deconvolution, the residual map should be re-scaled with the Tp2vistweak function. This rescaling is required due to the typical disparity between the areas of the synthesized beam and convolution beam.

The capability of Tp2vis was demonstrated using model maps. The joint deconvolution of 12m, 7m, and TP visibilities appeared very successful and reproduced the true emission distributions typically within 5% errors. The quality of images, of course, depends not only on the presence/absence of TP visibilities, but also on the distribution and contrast of the emission. Tests in wider parameter spaces with different deconvolvers are beyond the scope of this paper. Users are encouraged to do their own tests with model images that are similar to their objects of interest. We should note that our test was done under the idealized condition of no noise, but the 5% accuracy over the dynamic range of $\gtrsim 2$ orders of magnitude is quite encouraging. A manual of usage is presented in the Github repository (https://github.com/tp2vis/distribute).

We thank the anonymous referee for useful comments, the NRAO staff, in particular, Remy Indebetouw, Kumar Golap, Jennifer Donovan Meyer, Crystal Brogan, and John Carpenter for their help. We also thank Kazuki Tokuda, Fumi Egusa, Manuel Fernández, and Mercedes Vazzano for feedback on an early version of Tp2vis and Jim Barrett for comments on drafts. We acknowledge a research grant from the ALMA Cycle 4 Development Study program at NRAO. JK acknowledges supports from the National Astronomical Observatory of Japan and the Joint ALMA Observatory during his stay on sabbatical.

Appendix A GMC-like Model Map for Testing Interferometry Imaging

For tests of interferometer imaging procedures, we produce images of density/emission distributions similar to those in giant molecular clouds (GMCs). The method described below is similar to the one used for generating initial conditions for cosmological simulations of large scale structures or galaxy formation, and to the one presented by Dubinski et al. (1995) for molecular cloud-like density fields. We generate a gaussian random field that follows a density power spectrum (Section A.1) and manipulate it to enhance the contrast (Section A.2).

A.1 Model Density Field

We start from a representation of a density field with the density power spectrum $P(k)\propto\left<\widetilde{\delta}^{2}\right>\propto k^{-n}$. We assume a Gaussian random field (i.e., Fourier modes are uncorrelated) and the probability distribution function of its amplitude is drawn from a Gaussian distribution.

A.1.1 Definitions

Using the density at position $\mathbb{x}$, $\rho(\mathbb{x})$, and average density, $\rho_{0}$, the density contrast is defined as

[TABLE]

For the discrete Fourier transformation, we set the grid coordinates

[TABLE]

where $(l,m,n)$ and $(u,v,w)$ are integers, and $(L,M,N)$ are the numbers of grid points in $x$, $y$, and $z$ directions, and $D$ is the size of the image, respectively. We adopt the definition of the discrete Fourier transformation:

[TABLE]

This definition is adopted in most numerical packages (e.g., IDL, Python Numpy package, etc) and is convenient for implementation. The last terms, $\Delta u\Delta v\Delta w$ and $\Delta l\Delta m\Delta n$ (both = 1), are written explicitly because they need an evaluation later. Since the density field, $\delta(l,m,n)$, is a field of real numbers, each complex Fourier component has to satisfy

[TABLE]

A.1.2 Power Spectrum and Normalization

The density power spectrum is often defined in the form of a power law $P(k)=Ck^{-n}$ with a normalization constant of $C$. To avoid the singularity at $k=0$ the power spectrum is redefined here as

[TABLE]

where $k_{\rm min}$ is the minimum wavenumber ($k_{\rm min}=1/D$).

The normalization constant $C$ is calculated from the variance of the density contrast, $\sigma^{2}$ (input parameter):

[TABLE]

For the purpose of analytical evaluation, we can also obtain

[TABLE]

where the last equation is valid when $n>3$ (the integral diverges when $n\leq 3$). $\Gamma$ is the gamma function. [Note $\sigma^{2}\propto k_{\rm min}^{3-n}\propto D^{n-3}$, and the size-line width relation of GMCs indicates $n\sim 4$.] For a given set of $n$ ($>3$) and $\sigma$, the normalization coefficient of the power spectrum is

[TABLE]

The above prescription generates a 3-dimensional data cube and involves a Fourier transformation in 3-d. We can generate a 2-d image and apply a Fourier transformation in 2-d by using

[TABLE]

where

[TABLE]

A.1.3 Realizations

The power spectrum $P(k)$ and the fluctuation of density contrast $<\widetilde{\delta}^{2}(k)>$ are related as follows.

[TABLE]

The $\{...\}$ term is equal to $\left\{L\delta^{K}_{u,u^{\prime}}\cdot M\delta^{K}_{v,v^{\prime}}\cdot N\delta^{K}_{w,w^{\prime}}\right\}$, where $\delta^{K}$ is the Kronecker delta. The bracket terms at the end, $\left[\Delta l\Delta m\Delta n\right]$, $\left[\Delta u\Delta v\Delta w\right]$, etc, are expressed as $[...]$ to save space. From equations (A16) and (A17),

[TABLE]

or

[TABLE]

Note that $\left[\Delta u\Delta v\Delta w\right]=1$ and $\left[\Delta k_{u}\Delta k_{v}\Delta k_{w}\right]=1/D^{3}$.

The Fourier components $\widetilde{\delta}(k)$s are sampled using $<\widetilde{\delta}^{2}(k)>$ as the variance. A random number is generated for each $\mathbb{k}$ using a Gaussian distribution with a standard deviation of $\sqrt{<\widetilde{\delta}^{2}(k)>}=\sqrt{P(k)d^{3}\mathbb{k}}=\sqrt{P(k)/D^{3}}$. The phase of the complex number $\widetilde{\delta}(k)$ is also drawn from a uniform distribution in $[0,2\pi)$, and the symmetry condition, eq. (A5), is taken into account. Once all the Fourier components are generated, they are Fourier transformed to the image plane using equation (A3). We could also draw the random numbers without the symmetry condition and take the real or imaginary part after the Fourier transformation. In this case, the amplitude should be multiplied by $\sqrt{2}$.

Figures 6a,b show the model images generated with a power spectrum index of $n=4$, an average density fluctuation of $\sigma=0.3$, and an image size of $16,384^{2}$. Both the real and imaginary images are shown. The images have periodic boundaries, and we spatially shifted them so that significant features come near the image centers.

A.2 Enhance Density Contrasts

Figure 6a,b qualitatively resemble some parts of GMCs. However, by construction, they don’t have “dense cores”, whose environs are often observed with ALMA. We could develop dense cores by running hydrodynamics simulations with gravity, starting from these images as the initial condition, but this is obviously beyond the scope of this paper. Here, to mimic their high brightness, we simply convert the flux $f^{\rm old}$ (or density $\rho(\mathbb{x})$) into new flux $f^{\rm new}$ using the following equation:

[TABLE]

We adopt $\beta=0.7$. This conversion is applied three times to make the images resemble GMCs with some dense cores. The final images are displayed in Figure 6c,d in a log scale. They have larger dynamic ranges than the original ones, and thus, are more challenging to reproduce with interferometer imaging.