EDP Sciences
Free Access
Issue
A&A
Volume 517, July 2010
Article Number A12
Number of page(s) 21
Section Astronomical instrumentation
DOI https://doi.org/10.1051/0004-6361/200912873
Published online 26 July 2010

© ESO, 2010

1. Introduction

The instantaneous field of view of an interferometer is naturally limited by the primary beam size of the individual antennas. For the ALMA 12 m-antennas, this field of view is  ~9′′ at 690 GHz and  ~27′′ at 230 GHz. The astrophysical sources in the (sub)-millimeter domain are often much larger than this, but still structured on much smaller angular scales. Interferometric wide-field techniques enable us to fully image these sources at high angular resolution. These techniques first require an observing mode that in one way or another scans the sky on spatial scales larger than the primary beam. The most common observing mode in use today, known as stop-and-go mosaicking, consists in repeatedly observing sky positions typically separated by half the primary beam size. The improvement of the tracking behavior of modern antennas now leads astronomers to consider on-the-fly observations, with the antennas slewing continuously across the sky. The improvements in correlator and receiver technologies are also leading to techniques that could potentially sample the antenna focal planes with multi-beam receivers instead of the single-pixel receivers installed on current interferometers.

The ideal measurement equation of interferometric wide-field imaging is (1)where V is the visibility function of 1) up (the spatial frequency with respect to the fixed phase center) and 2) αs (the scanned sky angle), I is the sky brightness, and B the antenna power pattern or primary beam of an antenna of the interferometer (Thompson et al. 1986, Chap. 2). For simplicity, 1) we assume that the primary beam is independent of azimuth and elevation, and 2) we use one-dimensional notation without loss of generality. We do not deal with polarimetry (see e.g. Hamaker et al. 1996; Sault et al. 1996a, 1999) because it adds another level of complexity over our first goal here, i.e. wide-field considerations. Several aspects make Eq. (1) peculiar with respect to the ideal measurement equation for single-field observations. First, the visibility is a function not only of the uv spatial frequency (up) but also of the scanned sky coordinate (αs). Second, Eq. (1) is a mix between a Fourier transform and a convolution equation. It can be regarded, for example, as the Fourier transform along the αp dimension of the function, B(αp − αs)   I(αp), of the (αp, αs) variables. But Eq. (1) can also be written as the convolution: (2)where (3)and (4)For each up kept constant, V(up, αs) is the convolution of ℬ and ℐ. Indeed, ℐ(αp,up = 0) = I(αp), so we derive (5)i.e., the convolution equation for single-dish observations.

Ekers & Rots (1979) were the first to show that the measurement equation (Eq. (1)) enables us to recover spatial frequencies of the sky brightness at a much finer uv resolution than the uv resolution set by the diameter of the interferometer antennas. Interestingly enough, the goal of Ekers & Rots (1979) was “just” to find a way to produce the missing short spacings of a multiplying interferometer. However, Cornwell (1988) realized that Ekers & Rots’ scheme has a much stronger impact, because it explains why an interferometer is able to do wide-field imaging. Cornwell (1988) also demonstrated that on-the-fly scanning is not absolutely necessary to interferometric wide-field imaging. Indeed, the large-scale information can be retrieved in mosaics of single-field observations, provided that the sampling of the single fields follows the sky-plane Nyquist sampling theorem.

As a result, all the information about the sky brightness is coded in the visibility function. From a data-processing viewpoint, all the current radio-interferometric wide-field imaging methods (see, e.g., Gueth et al. 1995; Sault et al. 1996b; Cornwell et al. 1993; Bhatnagar & Cornwell 2004; Bhatnagar et al. 2008; Cotton & Uson 2008) belong to the mosaicking family1 pioneered by Cornwell (1988). In this family, the processing starts with Fourier transforming V(up, αs) along the up dimension (i.e. at constant αs) to produce a set of single-field dirty images before linearly combining them and forming a wide-field dirty image. In this paper, we propose an alternate processing, which starts with a Fourier transform of V(up, αs) along the αs dimension (i.e. at constant up). We show how this explicitely synthesizes the spatial frequencies needed to do wide-field imaging, which are linearly combined to form a “wide-field uv plane”, i.e., one uv-plane containing all the spatial frequency information measured during the wide-field observation. Inverse Fourier transform will produce a dirty image, which can then be deconvolved using standard methods. The existence of two different ways to extract the wide-field information from the visibility function raises several questions: are they equivalent? What are their relative merits?

We thus aim at revisiting the mathematical foundations of wide-field imaging and deconvolution. Sections 2 to 7 propose the new algorithm, which we call wide-field synthesis: Sect. 2 first defines the notations and it then lays out the basic concepts used throughout the paper. Section 3 states the steps needed to go beyond the Ekers & Rots scheme and explores the consequences of incomplete sampling of both the uv and sky planes. Section 4 discusses the effects of gridding by convolution and regular resampling. Section 5 describes how to influence the dirty beam shapes and thus the deconvolution. Section 6 states how to introduce short spacings measured either from a single-dish antenna or from heterogeneous interferometers. Section 7 compares the proposed wide-field synthesis algorithm with standard nonlinear mosaicking. Some detailed demonstrations are factored out in Appendix A to enable an easier reading of the main paper, while ensuring that interested readers can follow the demonstrations. Appendices B and C then explain how the wide-field synthesis algorithm can cope with non-ideal effects: Appendix B discusses how at least one standard way to cope with sky projection problems is compatible with the wide-field synthesis algorithm. Appendix C explores the consequences of using the on-the-fly observing mode. Finally, we assume good familiarity with single-field imaging in various places. We refer the reader to well-known references: e.g. Chap. 6 of Thompson et al. (1986) and Sramek & Schwab (1989).

2. Notations and basic concepts

2.1. Notations

In this paper, we use the Bracewell (2000)’s notation to display the relationship between a function I(α) and its direct Fourier transform (u), i.e., (6)where (α, u) is the couple of Fourier conjugate variables. We also use the following sign conventions for the direct and inverse Fourier transforms (7)and (8)As V is a function of two independent quantities (up and αs), the Fourier transform may be applied independently on each dimension, while the other dimension stays constant. Several additional conventions are used to express this. First, we introduce a specific notation to state that either the first or the second dimension stays constant: (9)and (10)Second, we use a bottom/top line to derive the notation of the Fourier transform along the first/second dimension from the notation of the original function. Third, on the Fourier transform sign (i.e.  ⊃ ), we explicitly state the dimension along which the Fourier transform is computed. For instance, if D is a function of (αp, αs), then the Fourier transform of D along the first dimension is expressed as (11)while the Fourier transform of D along the second dimension is expressed as (12)We also use a more compact notation when doing the Fourier transform on both dimensions simultaneously, i.e., (13)Finally, the convolution of two functions G and V is noted and defined as (14)For reference, Table 1 summarizes the definitions of the symbols used most throughout the paper. With the one-dimensional notation used throughout the paper, the number of planes quoted directly gives the number of associated dimensions of the symbols. Generalization to images would require a doubling of the number of planes/dimensions. Table 2 defines the uv and angular scales that are relevant to wide-field interferometric imaging, and Fig. 1 sketches the different angular scales. Each angular scale (θ) is related to a uv scale (d) through θ = 1/d, where θ and d are measured in radians and in units of λ (the wavelength of the observation). In the rest of the paper, we explicitely distinguish between θprim ≡ 1/dprim, the angular scale associated to the diameter of the interferometer antennas, and θfwhm, the full width at half maximum of the primary beam. The relation between θprim and θfwhm depends on the illumination of the receiver feed by the antenna optics. In radio astronomy, we typically have θfwhm ~ 1.2   θprim (see e.g. Goldsmith 1998, Chap. 6). Finally, the notion of anti-aliasing scale (θalias) is introduced and discussed in Sect. 4.2.

Table 1

Definition of the symbols used to expose the wide-field synthesis formalism.

Table 2

Definition of the uv and sky scales relevant to wide-field interferometric imaging.

thumbnail Fig. 1

Visualization of the different angular scales relevant to wide-field interferometric imaging. The notion of anti-aliasing scale (θalias) is introduced and discussed in Sect. 4.2.

Open with DEXTER

2.2. Basic concepts

Figure 2 illustrates the principles underlying 1) the setup to get interferometric wide-field observations and 2) our proposition to process them. For simplicity, we display the minimum possible complexity without loss of generality. The top row displays the sky plane. The middle row represents the 4-dimensional measurement space at different stages of the processing. As it is difficult to display a 4-dimensional space on a sheet of paper, the bottom row shows 2-dimensional cuts of the measurement space at the same processing stages.

thumbnail Fig. 2

Illustration of the principles of wide-field synthesis, which enables us to image wide-field interferometric observations. The top row displays the sky plane. The middle row displays the 4-dimensional visibility space and the bottom row shows 2-dimensional cuts of this space at different stages of the processing. In panels b) to d), the scanned dimensions (αs and us) are displayed in blue while the phased spatial scale dimensions (up) are displayed in red and the spatial scale dimensions (u) of the final wide-field uv plane are displayed in black. The grey zones of panels b.2) and c.2) show the regions of the visibility space without measurements (missing short-spacings). In detail, panel a) shows a possible scanning strategy of the sky to measure the unknown brightness distribution at high angular resolution: for simplicity it is here just a 7-field mosaic. Panels b.1) and b.2) sketch the space of measured visibilities: the uv plane at each of the 7 measured sky positions is displayed as a blue square box in panel b.1) and a blue vertical line in panel b.2). For simplicity, only 6 visibilities are plotted in panel b.1). Panels c.1) and c.2) sketch the space of synthesized visibilities after Fourier transform of the measured visibilities along the scanned coordinate (αs): at each measured spatial frequency up (displayed on the blue axes) is associated one space of synthesized wide-field spatial frequencies displayed as one of the red squares in panel c.1) and the red vertical lines in panel c.2). The wide-field spatial scales are synthesized 1) on a grid whose cell size is related to the total field of view of the observation and 2) only inside circles whose radius is the primary diameter of the interferometer antennas. Panels d.1) and d.2) display the final, wide-field uv plane. This plane is built by application of the shift-and-average operator along the black lines on panel c.2), lines that display the region of constant u spatial frequency in the (up, us) space. Standard inverse Fourier transform and deconvolution methods then produce a wide-field distribution of sky brightnesses as shown in panel e).

Open with DEXTER

2.2.1. Observation setup and measurement space

Panel a) displays the sky region for which we aim for estimating the sky brigthness, I(α). The field of view of an interferometer observing in a given direction of the sky has a typical size set by the primary beam shape. In our example, this is illustrated by any of the circles whose diameter is θprim. As we aim at observing a wider field of view, e.g. θfield, the interferometer needs to scan the targeted sky field. We assume that we scan through stop-and-go mosaicking, ending up with a 7-field mosaic.

After calibration, the output of the interferometer is a visibility function, V(up, αs), whose relation to the sky brightness is given by the measurement equation (Eq. (1)). Panel b.1) shows the measurement space as a mosaic of single-field uv planes: the uv plane coverage of each single-field observation is displayed as a blue sub-panel at the sky position where it has been measured and which is featured by the red axes. We assume 1) that the interferometer has only 3 antennas and 2) that only a single integration is observed per sky position. This implies only 6 visibilities per single-field uv plane. In panel b.2), the uv planes at constant αs are displayed as the blue vertical lines. The measured spatial frequencies belong to the  [−dmax, − dmin]  and  [ + dmin, + dmax]  ranges, where dmin and dmax are respectively the shortest and longest measured baseline length. dmin is related to the minimum tolerable distance between two antennas to avoid collision. Here, we chose dmin ~ 1.5   dprim. The grey zone between  − dmin, and  + dmin displays the missing short spacings.

2.2.2. Processing by explicit synthesis of the wide-field spatial frequencies

All the information about the sky brightness, I(α), is somehow coded in the visibility function, V(up, αs). The high spatial frequencies (from dmin to dmax) are clearly coded along the up dimension. The uncertainty relation between Fourier conjugate quantities also implies that the typical spatial frequency resolution along the up dimension is only dprim because the field of view of a single pointing has a typical size of θprim. However, wide-field imaging implies measuring all the spatial frequencies with a finer resolution, dfield = 1/θfield. The missing information must then be hidden in the αs dimension.

In Sect. 3, we show that Fourier transforming the measured visibilities along the αs dimension (i.e. at constant up) can synthesize the missing spatial frequencies, because the αs dimension is sampled from  − θfield/2 to  + θfield/2, implying a typical spatial-frequency resolution of the us dimension equal to dfield. Conversely, the αs dimension is probed by the primary beams with a typical angular resolution of θprim, implying that the us spatial frequencies will only be synthesized inside the  [ − dprim, + dprim]  range. Panels c.1) and c.2) illustrate the effects of the Fourier transform of V(up, us) along the αs dimension, in 4 and 2 dimensions, respectively. The red subpanels or vertical lines display the us spatial frequencies around each constant up spatial frequency.

In panels d.1) and d.2) (i.e. after the Fourier transform along the αs dimension), (up, us) contains all the measured information about the sky brightness in a spatial frequency space. However, the information is ordered in a strange and redundant way. Indeed, we show that (up, us) is linearly related to (up+ us). To first order, the information about a given spatial frequency u is stored in all the values of (up, us) which verifies u = up + us (black lines on panel c.2).

A shift operation will reorder the spatial scale information and averaging will compress the redundancy (illustrated by the halving of the number of the space dimensions). The use of a shift-and-average operator thus produces a final uv plane containing all the spatial scale information to image a wide field in an intuitive form. We thus call this space the wide-field uv plane. Panels d.1) and d.2) display this space, where the minimum relevant spatial frequency is related to the total field of view, while the maximum one is related to the interferometer resolution.

Sections 3 and 4 show that applying the shift-and-average operator to produces the Fourier transform of a dirty image, which is a local convolution of the sky brightness by a slowly varying dirty beam. As a result, inverse Fourier transform of and deconvolution methods will produce a wide-field distribution of sky brightness as shown in panel e) at the top right of Fig. 2.

3. Beyond the Ekers & Rots scheme

In the real world, the visibility function is not only sampled, but this sampling is incomplete for two main reasons. 1) The instrument has a finite spatial resolution, and the scanning of the sky is limited, implying that the sampling in both planes has a finite support. 2) The uv coverage and the sky-scanning coverage can have holes caused either by intrinsic limitations (e.g. lack of short spacings or small number of baselines) or by acquisition problems (implying data flagging). The incomplete sampling makes the mathematics on the general case complex. We thus start with the ideal case where we assume that the visibility function is continuously sampled along the up and αs dimension. We then look at the general case.

3.1. Ideal case: infinite, continuous sampling

Starting from the measurement Eq. (1), Ekers & Rots (1979) first demonstrated (see Sect. A.1) that2(15)For each constant up spatial frequency, the Fourier transform thus synthesizes a function, , which is simply related to (up+ us), the Fourier components of the sky brightness around up. (up, us) is only defined in the  [ − dprim, + dprim]  interval along the us dimension because ( − us) is itself only defined inside this interval, since ( − us) is the autocorrelation of the antenna illumination.

We search to derive a single estimate of the Fourier components of the sky brightness. Equation (15) indicates that the fraction gives us an estimate of (u) for each couple (up,us) that satisfies u = up + us. However, the information about is strangely ordered. There are two possible ways to look at this ordering. 1) Starting from the measurement space, the Ekers & Rots scheme synthesizes frequencies around each up measure inside the interval  [up − dprim,up + dprim]  at the dfield spatial frequency resolution. 2) Starting from our goal, we want to estimate at a given spatial frequency u with a dfield spatial frequency resolution. We thus search for all the couples (up,us) satisfying u = up + us, which are displayed in panel c.2) of Fig. 2 as the diagonal black lines. It immediately results that 1) there are several estimates of for each spatial frequency u and 2) the number of estimates varies with u. We can average them to get a better estimate of (u).

This last viewpoint thus suggests averaging in the (up, us) space along linepaths defined by u = up + us. Such an operator can mathematically be defined as (16)where F is the function to be averaged and W is a normalized weighting function. Using the properties of the Dirac function, we can reduce the double integral to (17)In this equation, we easily recognize a shift-and-average operator. The normalized weighting function plays a critical role in the following formalism, and we propose clever ways to define W in Sect. 5. In the ideal case studied here, W can be defined as

  • for us   in    [ − dprim, + dprim] ,

  • W(up,us) ≡ 0

    for other values of us.

In other words, we have just normalized the integral by the constant length of the averaging linepath.

3.1.1. Wide-field dirty image, dirty beam and image-plane measurement equation

Section 3.2 shows that the incomplete sky and uv sampling forbid us to apply the shift-and-average operator to the function. To guide us in this general case, we thus explore the consequences of applying this operator to in the ideal case. It is easy to demonstrate that the result is the Fourier transform of a dirty image, i.e., (18)Indeed, substituting with the help of Eqs. (17) and (15) and taking the inverse Fourier transform, we get (19)with (20)Here, Idirty conforms to the usual idea of dirty image, i.e., the convolution of a dirty beam by the sky brightness: (21)In contrast to the usual situation for single-field observations, the mix between a Fourier transform and a convolution of Eq. (1), associated with the specific processing3 changes the image-plane measurement equation from a convolution of a dirty beam with the product B   I to a convolution of a dirty beam with I. The dependency on the primary beam is still there. It is just transferred from a product of the sky brightness distribution into the definition of the dirty beam.

3.1.2. Summary and interpretation

In summary, a theoretical implementation of wide-field synthesis implies

  • 1.

    the possibility of Fourier transforming the visibility functionalong the αs dimension (i.e. at constant up), which gives us a set ofsynthesized uv planes;

  • 2.

    the possibility of shifting-and-averaging these synthesized uv planes to build the final, wide-field uv plane containing all the available information.

Using those tools, we are able to write the wide-field image-plane measurement equation as a convolution of a wide-field dirty beam (D) by the sky brightness (I), i.e., (22)We can write a convolution equation in this ideal case because the wide-field response of the instrument is shift-invariant; i.e., D only depends on differences of the sky coordinates.

It is well-known that for a single-field observation, the dirty beam is the inverse Fourier transform of the sampling function. The shape of this sampling function is due to the combination of aperture synthesis (the interferometer antennas give a limited number of independent baselines) and Earth-rotation synthesis (the rotation of the Earth changes the projection of the physical baselines onto the plane perpendicular to the instantaneous line of sight). By analyzing via a Fourier transform, the evolution of the visibility function with the sky position, the Ekers & Rots scheme synthesizes visibilities at spatial frequencies needed to image a larger field of view than the interferometer primary beam. We thus propose to call this specific processing: wide-field synthesis.

3.2. General case: incomplete sampling

Reality imposes limitations on the synthesis of spatial frequencies. Indeed, we have already stated that the visibility function is incompletely sampled both in the uv and sky planes. To take the sampling effects into account, we introduce the sampling function S(up, αs), which is a sum of Dirac functions at measured positions4. The sampling function cannot be factored into the product of two functions, each only acting on one plane. Indeed, the Earth rotation happening during the source scanning implies a coupling of both dimensions of the sampling function. In other words, the uv coverage will vary with the scanned sky coordinate. This leads us to a shift-dependent situation, precluding us from writing the wide-field image-plane measurement equation as a true convolution. We nevertheless search for a wide-field image-plane measurement equation as close as possible to a convolution because all the inversion methods devised in the past three decades in radioastronomy are tuned to deconvolve images. The simplest mathematical way to generalize Eq. (22) to a shift-dependent situation is to write it as (23)In this section, we show how the linear character of the imaging process allows us to do this. Section 3.2.1 derives the impact of incomplete sampling on the Ekers & Rots equation, and Sect. 3.2.2 derives the wide-field measurement equation in the uv plane. Section 3.2.3 interprets these results.

3.2.1. Effect on the Ekers & Rots equation

The sampled visibility function, SV, is defined as the product of S and V and its Fourier transform along αs, i.e., (24)and (25)Because SVup is the product of two functions of αs, we can use the convolution theorem to show that is the convolution of by , i.e., (26)By replacing with the help of the Ekers & Rots relation (Eq. (15)), we derive (27)As is bounded inside the  [ − dprim, + dprim]  interval, is a local average, weighted by , of around the up spatial frequency.

As expected, we recover Eq. (15) for the ideal case (i.e., infinite, continuous visibility function) because then . A more interesting case arises when the visibility function is continuously sampled over a limited sky field of view, i.e., After Fourier transform this gives (30)In this case, the local average of the sky brightness Fourier components happens on a typical uv scale equal to dfield. However, the sinc function is known to decay only slowly. Some observing strategy (e.g. quickly observing outside the edges of the targeted field of view to provide a bandguard) could be considered to apodize the sky-plane dependence of the sampling function, resulting in faster decaying functions, hence in less mixing of the wide-field spatial frequencies.

3.2.2. uv-plane wide-field measurement equation

Because we aim at estimating the Fourier component of , we introduce the following change of variables and , to derive (31)We then shift-and-average (up, us) to build the Fourier transform of a wide-field dirty image (32)Substituting the shift-and-average operator by its definition and using Eq. (31) to replace , we derive (33)This uv-plane wide-field measurement equation can be written as (34)if we enforce the following equality (35)This is one way to define , which is convenient though unusual. It is implicit in this definition that we need to make a change of variable (u′′ = u − u′) to derive (36)In the following, we use either one or the other definition of , depending on convenience.

3.2.3. Interpretation

Appendix A.2 demonstrates that the image and uv-plane wide-field measurement equations (Eqs. (23) and (34)) are equivalent if (37)The image-plane wide-field measurement equation (Eq. (23)) can be written as (38)Its interpretation is straightforward: the sky brightness distribution is convolved with a dirty beam, D(α′, α′′), which varies with the sky coordinate α′′. This raises the question of the rate of change of the dirty beam with the sky coordinate. This question is addressed in Sects. 4.2 and 5.

4. Gridding by convolution and regular resampling

We want to Fourier transform the raw visibilities along the sky dimension (αs) at some constant value in the up dimension. The raw data, however, is sampled on an irregular grid in both the uv and sky planes. We need to grid the measured visibilities in both the uv and the sky planes before Fourier transformation for different reasons. First, the gridding in the uv plane will handle the variation in the spatial frequency as the sky is scanned, i.e., the difficulty and perhaps the impossibility of Fourier-transforming at a completely constant up value. Second, the gridding along the sky dimension allows the use of Fast Fourier Transforms. As usual, we grid through convolution and regular resampling.

4.1. Convolution

4.1.1. Definitions

We first define a gridding kernel that depends on both dimensions, . This gridding kernel can be chosen as the product of two functions, simplifying the following demonstrations: (39)We then define the sampled visibility function gridded in both the uv and sky planes as Finally, when assessing the impact of the gridding on the measurement Eq. (34), a new function, (42)and its Fourier transforms naturally appear in the equations. Defining the following Fourier transform relationships (43)and (44)we easily derive (45)and (46)Using these notations, we have before gridding, (47)and (48)

4.1.2. Conservation of the wide-field measurement equation

Appendix A.3 demonstrates that the wide-field dirty image is here again the convolution of the sky brightness I by a wide-field dirty beam Dα or, in the Fourier plane, (49)with (50)where (51)We thus have equations that resemble those containing the sampling function alone, except for 1) the replacement of the generalized sampling function by its gridded version and 2) the way the variables are linked together both in the gridding of (i.e., Eq. (51)) and in the averaging of (i.e., Eq. (50)).

4.2. Regular resampling

It is well known that too low a resampling rate in one space implies power aliasing in the conjugate space (see e.g. Bracewell 2000; Press et al. 1992). Aliasing must be avoided as much as possible because it folds power outside the imaged region back into it. Table 3 defines the intervals of definition of the different functions we are dealing with (i.e., visibilities, primary beam, dirty image, and dirty beam), as well as the associated sampling rates needed to enforce Nyquist sampling. The boundary values of the definition intervals (|u|max and |α|max) are related to the sampling rates (∂α and ∂u, respectively) through (52)where nsamp is an integer characterizing the sampling. Nyquist sampling implies nsamp = 2. However, slight oversampling (e.g. nsamp = 3) is often recommended because the measures suffer from errors and the deconvolution is a nonlinear process. In this section, we examine the properties of the different functions to define their associated sampling rates.

Table 3

Interval ranges of definition and associated sampling rates for the used functions.

4.2.1. The αs sampling rate of the visibility function

When Fourier transforming the measurement Eq. (1) along the αs axis, we derive the Ekers & Rots Eq. (15). This equation implies that (up, us) is bounded inside the  [ − dprim, + dprim]  spatial frequency interval along the us axis. As a result, the visibility function needs to be regularly resampled at a rate of only 0.5/dprim to satisfy the Nyquist theorem. This was first pointed out by Cornwell (1988). This sampling rate is equal to θprim/2 or  ~θfwhm/2.4. The “usual, wrong” habit of sampling at θfwhm/2 is indeed undersampling with aliasing as a consequence. Mangum et al. (2007) discuss the consequences of undersampling in-depth in the framework of single-dish imaging.

4.2.2. The Up sampling rate of the visibility function

Now, the Fourier transform of the measurement Eq. (1) along the up axis gives (53)where (54)We use the tilde sign under V to denote the inverse Fourier transform of V along its first dimension. A well-known Fourier transform property implies that B has infinite support because is bounded. The resampling rate along the up axis therefore depends on the properties of the product of B(αp − αs) times I(αp) as a function of αp. While no unique answer exists, three facts help us to find the right sampling rate: 1) B falls off relatively quickly; 2) the result depends on the spatial distribution of the sky brightness and in particular on the dynamic range in brightness needed to accurately image it; 3) the measure of (αp, αs) has a limited accuracy owing to thermal noise, phase noise, and other possible systematics (e.g. pointing errors). For simplicity, we quantify the measurement accuracy by a single number, namely the maximum instrumental fidelity measured in the image plane as defined in Pety et al. (2001). There are two cases:

  • 1.

    the maximum instrumental fidelity limits the dynamic range inbrightness. For instance, Petyet al. (2001) showed that thefidelity of interferometric imaging at (sub)-millimeterwavelengths will be limited to a few hundred. In this case, (αp,αs) aliasing can be tolerated when the amplitude of B is less than afraction of the inverse of the maximum instrumental fidelity;

  • 2.

    the maximum instrumental fidelity is much greater than the image fidelity, as can be the case at centimeter wavelengths. In this case, (αp, αs) aliasing can only be tolerated when the amplitude of B is less than a fraction of the inverse of the dynamic range of the image.

The criterion derived in each case gives a typical image size (θalias), which can be converted into the desired up sampling rate. To be more quantitative, Fig. 3 models the normalized antenna power patterns of an antenna illuminated by a Gaussian beam of 12.5 dB edge taper and with a given blockage factor (ratio of the secondary-to-primary diameters). The top panel presents an ideal case without secondary miror, while the middle and bottom panels present simple models of the ALMA and PdBI antennas. The largest angular sizes at which the power patterns are less than a given value, , is a first-order estimate of θalias/2 to get a fidelity or dynamic range higher than . Table 4 gives the values of θalias/θfwhm as a function of the searched fidelity or dynamic range. This condition is sufficient but not necessary. Indeed, the aliasing properties also depend on the brightness distribution of the source.

thumbnail Fig. 3

Simple models of the antenna power patterns as a function of the sky angle in units of half the primary beam FWHM (θfwhm). In the 3 cases shown, the illumination is Gaussian with an edge taper of 12.5 dB but 3 different ratios of the secondary-to-primary diameters (i.e. fb, the antenna blockage factors) are considered (see e.g. Goldsmith 1998, Chap. 6). The middle and bottom panels respectively model ALMA and PdBI antennas. The red lines define the minimum angular sizes for which the antenna power pattern is less than a given fraction.

Open with DEXTER

Table 4

Minimum sizes of the dirty beam images to get an image fidelity or a dynamic range greater than a given value.

4.2.3. The u sampling rate of (u)

We have no garantee that the sky outside the targeted field of view is devoid of signal, so the only way to ensure a given dynamic range inside the targeted field of view is to choose the image size large enough so that the aliasing of potential outside sources is negligible. This means that the dirty image size must be equal to the field-of-view size plus the tolerable aliasing size (55)The conjugate uv distance and associated uv sampling then are (56)

4.2.4. The u′ and u′′ sampling rates of (u′,u′′)

The u′′ axis must thus be sampled at the same rate as the second dimension of the definition space of , i.e., as us. Moreover, u′ has in this equation a behavior similar to u ( = up + us). It must thus have the same sampling behavior as u. This sampling rate (u′ = dimage/nsamp) is quite high. Some deconvolution methods (see below) allow us to relax this sampling rate.

4.3. Absence of gridding “correction”

Imaging of single-field observations goes through the following steps: 1) convolution by a gridding kernel; 2) regular resampling; 3) fast Fourier transform; and 4) gridding “correction”. The so-called gridding “correction” is a division of the dirty beam and dirty image by the Fourier transform of the gridding kernel used in the initial convolution. This step is mandatory when imaging single-field observations to keep the image-plane measurement equation as a simple convolution equation (see e.g. Sramek & Schwab 1989). When imaging wide-field observations, as proposed here, the Fourier transform along the αs dimension, followed by the shift-and-average operation, freeze the convolution kernel into the dirty beam of the wide-field measurement equation. This is why the gridding “correction” step is irrelevant here.

5. Dirty beams, weighting, and deconvolution

In radioastronomy, the dirty beam is the response of the interferometer to a point source. In the wide-field synthesis framework, the response of the interferometer to a point source, D, a priori depends on the source position on the sky. D(α′, α′′) can thus be interpreted as a set of dirty beams, with each dirty beam referred to by its fixed α′′ sky coordinate. These simple facts raise several questions. What are the properties of the convolution kernel? Is it possible to modify these properties? How do we deconvolve the dirty image?

5.1. A set of wide-field dirty beams

With the wide-field synthesis framework proposed here, Appendix A.6 shows that (57)where (58)and (59)Δ(αp, αs) is the single-field dirty beam, associated with the uv sampling at the sky coordinate αs. And Ω(α′, α′′) will be called the image plane weighting function, while W(u′, u′′) is the uv plane weighting function. The set of wide-field dirty beams D is then the double convolution of the image plane weighting function and the single-field dirty beams, apodized by the primary beam at the current sky position αs.

While the shape of the single-field dirty beam is directly given by the Fourier transform of the sampling function, the shape of the wide-field dirty beam depends, directly or through Fourier transforms, on the sampling function (S), the primary beam shape (B), and the weighting function (W). Moreover, the wide-field dirty beam shape a priori varies slowly with the sky position, since it is basically constant over the primary beamwidth as stated in Sect. 4.2. It nevertheless varies, implying, for instance, a “slow” variation of the synthesized resolution over the whole field of view.

While the single-field and wide-field dirty beam expressions seem very different, they share the same property of expressing the way the interferometer is used to synthesize a telescope of larger diameter in the image plane. In other words, the sampling function for single-field imaging and for wide-field imaging express the sensitivity of the interferometer to a given spatial frequency. These uv space functions are called the transfer functions of the interferometer (Thompson et al. 1986, Chap. 5). Modifying the transfer function has a direct impact on the measured quantity. Once the interferometer is designed and the observations are done, the only way to change this transfer function is data weighting.

An ideal set of wide-field dirty beams, D(α′, α′′), would have the following properties. All the wide-field dirty beams should be identical (i.e., independent of the α′′ sky coordinate) and equal to a narrow Gaussian (its FWHM giving the image resolution). This would give the product of a wide Gaussian of uby a Dirac function of u′′, as the ideal wide-field transfer function, (u, u′′).

5.2. Dirty beam shapes and weighting

When imaging single-field observations, giving a multiplicative weight to each visibility sample is an easy way to modify the shape of the dirty beam and thus the properties of the dirty and deconvolved images. Natural weighting (which maximizes signal-to-noise ratio), robust weighting (which maximizes resolution), and tapering (which enhances brightness sensitivity at the cost of a lower resolution) are the most popular weighting techniques (see e.g. Sramek & Schwab 1989).

In the case of wide-field synthesis, a multiplicative weight can also be attributed to each visibility sample before any processing. However, the weighting is also at the heart of the wide-field synthesis because it is an essential part of the shift-and-average operation. No constraint has been set on the weighting function up to this point, which indicates that the weighting function (W) gives us a degree of freedom in the imaging process. We look in turn at both kinds of weighting. In both cases, an obvious issue is the definition of the optimum weighting functions. As in the case of single-field imaging, there is no single answer to this question. It depends on the conditions of the observation and on the imaging goals.

5.2.1. Weighting the measured visibilities

Natural weighting consists of slightly changing the definition of the sampling function. It is now set to a normalized natural weight where there is a measure and 0 elsewhere. The natural weight is usually defined as the inverse of the thermal noise variance, computed from the radiometric equation, i.e., from the system temperature, the frequency resolution, and the integration time. Using this weighting scheme before computing the first Fourier transform along the αs sky dimension makes sense because the observing conditions (and thus the noise) vary from visibility to visibility.

We propose to generalize this weighting scheme to other observing conditions than just the system noise. Indeed, critical limitations of interferometric wide-field imaging are pointing errors, tracking errors, atmospheric phase noise (in the (sub)-millimeter domain), etc. While techniques exist for coping with these problems (e.g., water vapor radiometer, direction-dependent gains: Bhatnagar et al. 2008), they are not perfect. The usual way to deal with the remaining problems is to flag the source data based on a priori knowledge of the problems, e.g., pointing measurement, tracking errors, rms phase noise on calibrators, etc. However, flagging involves the definition of thresholds, while reality is never black and white. It can thus be asked whether some weighting scheme could be devised to minimize the effect of pointing errors, tracking errors or phase noise on the resulting image. We propose to modulate natural weighting based on the a priori knowledge of the observing conditions.

5.2.2. Weighting the synthesized visibilities

Robust weighting or tapering the measured visibilities do not make sense in wide-field synthesis because the dirty image is made from the synthesized visibilities after the first Fourier transform along the αs sky dimension. A weighting function W then appears naturally as part of the shift-and-average operator. Its optimum value depends on the properties of the measured sampling function. Here are a few examples.

  • Infinite, continuous sampling.

    This is the ideal case studied inSect. 3.1. Knowing that the Ekers & RotsEq. (15) links the quantity we want to estimate, i.e.,, to many noisy5 measurements, (up, us), via a product by (assumed to be perfectly defined), we can invoke a simple least-squares argument (see e.g. Bevington & Robinson 2003) to demonstrate that the optimum weighting function is (60)with w(up, us) the weight computed from the inverse of the noise variance of (up, us). Using Eq. (20), it is then easy to demonstrate that , and then Idirty(α) = I(α). The dirty image is a direct estimate of the sky brightness; i.e., deconvolution is superfluous.

  • Complete sampling.

    The signal is Nyquist-sampled, but it has a finite support in both the uv and sky planes, implying a finite synthesized resolution and a finite field of view. In contrast to the previous case, this one may have practical applications, e.g., observations done with ALMA in its compact configuration. Indeed, the large number of independent baselines coupled to the design of the ALMA compact configuration ensure a complete, almost Gaussian, sampling for each snapshot. In this case, the best choice may be to choose the weighting function so that all the dirty beams are identical to the same Gaussian function. In this case, the deconvolution would also be superfluous.

  • This is the more general case studied in Sect. 3.2. The signal not only has a finite support but it also is undersampled (at least in the uv plane). The deconvolution is mandatory. The choice of the weighting function thus will depend on imaging goals. If the user needs the best signal-to-noise ratio, some kind of natural weighting will be needed. It is tempting to use Eq. (60) as a natural weighting scheme. However, the main condition for derivation of this weighting function, i.e., the Ekers & Rots Eq. (15), is not valid anymore, as the noisy measured quantity () is now linked to the quantity we want to estimate () by a local average (see Eq. (31)). This is why it was more appropriate to try to get a Gaussian dirty beam shape in the complete sampling case. If the signal-to-noise ratio is high enough, the user has two choices. Either he/she wants to maximize angular resolution power and needs some kind of robust weighting, or he/she wants to get the more homogeneous dirty beam shape over the whole field of view. This requirement cannot always be fully met. The Ekers & Rots scheme enables us to recover unmeasured spatial frequencies only in regions near to measured ones, because has a finite support.

5.3. Deconvolution

Writing the image-plane measurement equation in a convolution-like way is very interesting because all the deconvolution methods developed in the past 30 years are optimized to treat deconvolution problems (see e.g. Högbom 1974; Clark 1980; Schwab 1984; Narayan & Nityananda 1986). For instance, it should be possible to deconvolve Eq. (23) with just slight modifications to the standard CLEAN algorithms. Indeed, Eq. (23) can be interpreted as the convolution of the sky brightness by a set of dirty beams, so that the only change, once a CLEAN component is found, would be the need to find the right dirty beam in this set in order to remove the CLEAN component from the residual image.

Following Clark (1980) and Schwab (1984), most algorithms today deconvolve in alternate minor and major cycles. During a minor cycle, a solution of the deconvolution is sought with a simplified (hence approximate) dirty beam. During a major cycle, the current solution is subtracted either from the original dirty image using the exact dirty beam or from the measured visibilities, implying a new gridding step. In both cases, the major cycles result in greater accuracy. The iteration of minor and major cycles enables one to find an accurate solution with better computing efficiency. In our case, the approximate dirty beams used in the minor cycle could be 1) dirty beams of a much smaller size than the image; or 2) a reduced set of dirty beams (i.e., guessing that the typical variation sizescale of the dirty beams with the sky coordinate is much larger than the primary beamwidth); or 3) both simultaneously. The model would be subtracted from the original visibilities before re-imaging at each major cycle. The trade-off is between the memory space needed to store a full set of accurate dirty beams and the time needed to image the data at each major step. Some quantitative analysis is needed to know how far the dirty beams can be approximated in the minor cycle.

It is worth noting that the accuracy of the deconvolved image will be affected by edge effects. Indeed, the dirty brightness at the edges of the observed field of view is attenuated by the primary beam shape. When deconvolving these edges, the deconvolved brightness will be less precise, because the primary beam has a low amplitude there. This only affects the edges, because inside the field of view, every sky position should be observed a fraction of the time with a primary beam amplitude between 0.5 and 1. This edge effect is nevertheless expected to be much less troublesome than the inhomogeneous noise level resulting from standard mosaicking imaging (see Sect. 7.1).

6. Short spacings

6.1. The missing flux problem

Radio interferometers are bandpass instruments; i.e., they filter out not only the spacings longer than the largest baseline length but also the spacings shorter than the shortest baseline length, which is typically comparable to the diameter of the interferometer antennas. In particular, radio interferometers do not measure the visibility at the center of the uv plane (the so-called “zero spacing”), which is the total flux of the source in the measured field of view.

The lack of short baselines or short spacings has strong effects as soon as the size of the source is more than about 1/3 to 1/2 of the interferometer primary beam. Indeed, when the size of the source is small compared to the primary beam of the interferometer, the deconvolution algorithms use, in one way or another, the information of the flux at the lowest measured spatial frequencies for extrapolating the total flux of the source. The extreme case is a point source at the phase center for which the amplitude of all the visibilities is constant and equal to the total flux of the source: extrapolation is then exact. However, the larger the size of the source, the worse the extrapolation, which then underestimates the total source flux. This is the well-known problem of the missing flux that observers sometimes note when comparing a source flux measured by a mm interferometer with the flux observed with a single-dish antenna.

Wide-field synthesis does not recover the full short spacings. Let us assume that the visibility function is continuously sampled from dmin to dmax, with dmin ~ 1.5   dprim. The length of the averaging linepath6), L(u), can be interpreted as the number of measures that contribute to the estimation of . Figure 4 shows the variations of L(u) function when starting from a visibility function continuously defined in the  [dmin,dmax]  interval along the up dimension. We can expect to recover only inside the  [dmin − dprim,dmax + dprim]  interval. In particular, information on short spacings lower than dmin − dprim (e.g. the crucial zero spacing) cannot be recovered when using a homegeneous interferometer, and the short spacings in the interval  [dmin − dprim,dmin]  are recovered with increasing accuracy from dmin − dprim to dmin. Both effects imply the need for complementary instruments to accurately measure the missing short-spacings.

thumbnail Fig. 4

Length of the averaging linepaths displayed as black lines in panel c.2) of Fig. 2, as a function of the spatial scale in the final, wide-field uv plane. In the case of a continuous sampling of up between dmin and dmax, these quantities can be interpreted as the number of measures that contribute to the estimate of (u).

Open with DEXTER

6.2. Usual hardware and software solutions

To derive the correct result for larger source sizes, it is necessary to complement the interferometer data with additional data, which contain the missing short-spacing information. The IRAM-30 m single-dish telescope is used to complement the Plateau de Bure Interferometer. Short-spacing information can also be in part recovered with a secondary array of smaller antennas and shorter baselines (e.g. the CARMA interferometer). In the ALMA project, the short-spacing information will be derived by a combination of four 12 m-single-dish antennas and an interferometer of 12 antennas of 7 m called ACA (Atacama Compact Array).

From the software point-of-view, two main families of algorithms exist in the standard processing of mosaics. Either the short-spacing information is combined on the deconvolved image (i.e., the interferometer data is imaged and deconvolved separately) through a hybridization in the Fourier plane (see e.g. Pety et al. 2001), or the long and short-spacing information is imaged and/or deconvolved jointly. In this category, we find the pseudo-visibility technique, which produces interferometric-like visibilities from single-dish maps (see e.g. Pety et al. 2001; Rodríguez-Fernández et al. 2008, and references therein), and the multi-resolution deconvolution algorithms, which work on images containing different spatial frequency ranges.

In the next two sections, we show how wide-field synthesis naturally processes the short-spacing information either from single-dish or from heterogeneous arrays.

6.3. Processing short spacings from single-dish measurements

The single-dish measurement equation can be written as (61)where Isd is the measured single-dish intensity, Ssd the single-dish sampling function, and Bsd the single-dish antenna power pattern. As already stated in the introduction, the above integral is identical to the ideal measurement equation of interferometric wide-field imaging taken in up = 0. If we define a single-dish visibility function as (62)we can thus write the measured single-dish intensity as (63)The recognition that the single-dish measurement equation is a particular case of the interferometric wide-field measurement equation opens the way to treating both the single-dish and interferometric data sets through exactly the same processing steps. We just have to define a hybrid sampling function, Shyb, as the Fourier transform of the hybrid primary beam, , as and a hybrid weighting function, Whyb, as All the processing steps described in the previous sections (including a potential gridding step of single-dish, on-the-fly data) can then be directly applied to the hybrid data set. Using the wide-field synthesis formalism, we can easily write (70)with (71)and (72)We thus see that is a linear combination of the information measured by the single-dish () and by the interferometer (). There, Wsd(u) plays a particular role for two reasons. First, its dependency on the spatial frequency (u) enables us to filter out the highest spatial frequencies that are measured by the single-dish antenna with low accuracy. Second, it is well-known that the relative weight of the single-dish to interferometric data is a critical parameter in the processing of the short spacings from single-dish data (see e.g. Rodríguez-Fernández et al. 2008). This relative weight is a free parameter within the restrictions set by the noise level (i.e., we want the single-dish data to bring information and not just noise to the interferometric data), and a criterion must therefore be defined to adjust it to an optimal value. We refer the reader to the discussion of Sect. 5, which also applies here.

Table 5

Definition of the symbols used to expose the processing of the short spacings.

6.4. Processing short spacings from heterogeneous arrays

A heterogeneous array is an interferometer composed with antennas of different diameters. ALMA and CARMA are two such examples. The measurement equation for a heterogeneous array is (73)where bi and bj are the voltage reception patterns of the antenna pair that forms the ij baseline and the asterisk denotes the complex conjugate (Thompson et al. 1986, Chap. 3). The formalism developed in the previous sections holds as long as we redefine (74)A simple application of the correlation theorem implies that (75)The use of the baseline indices ij must be generalized throughout the equations because the knowledge of the antenna type must be attached to each individual data point (visibility). As a result, the wide-field synthesis formalism can be easily adapted to heterogeneous arrays at the price of additional bookkeeping.

6.5. Two textbook cases: IRAM-30 m + PdBI and ALMA + ACA

Figure 5 sketches why wide-field synthesis naturally handles the short spacings in two textbook cases. In the ideal case, the Fourier transform along the αs dimension produces visibilities, which are related to the wide-field spatial frequencies of the source brightness weighted by the transfer function of the interferometer. In this sense, Fig. 5 displays the natural weighting of the synthesized wide-field visibilities at the position of each measured visibility. Handling visibilities from antenna of different sizes just implies that the natural weighting function of the synthesized visibilities will have a different shape.

thumbnail Fig. 5

Sketches of the natural weighting of the synthesized wide-field visibilities. Each measured spatial frequency will produce wide-field spatial frequencies apodized by the transfer function () centered on the measured spatial frequency. The used transfer function depends on the telescopes used, explaining why wide-synthesis naturally handles the short spacing either from a single-dish antenna or from a heterogeneous array. The synthesized visibilities in the overlapping regions will then be averaged. Two textbook examples are illustrated: 1) the combination of data from the IRAM-30 m single-dish (red transfer function) and from the Plateau de Bure Interferometer (black transfer functions) at the top; and 2) the combination of data from ALMA 12 m-antennas used either in single-dish mode (red transfer function), in interferometric mode (black transfer functions) and of data from the ACA 7 m-antennas (blue transfer functions) at the bottom. The minimum uv distances measured by each interferometer were set from the minimum possible distance between antennas (24 m for PdBI, 15 m for ALMA and 9 m for ACA).

Open with DEXTER

The top panel of Fig. 5 displays how the IRAM-30 m single-dish is used to complement the Plateau de Bure interferometer visibilities. The bottom panel displays how ACA is used to produce the short spacing information for ALMA. The four 12 m-antennas will provide the single-dish information, while the 12 additional 7 m-antennas will form with ALMA a heterogeneous array. In the first design, ACA and ALMA form two independent interferometers; i.e., they are not cross-correlated. The single-dish antennas, ACA and ALMA, thus appear as three different instruments. It is thus possible to decompose the hybrid set of wide-field dirty beams obtained by processing the 3 sets of data together in 3 different sets of dirty beams (76)with (77)For a multiplying interferometer, (78)This implies that contributes at up = 0 in the sum over up in Eq. (77), contributes for 9   m < up ≲ 40   m and (u′, u′′) contributes for 15   m < up < 150   m in the most compact configuration of ALMA.

7. Comparison with standard nonlinear mosaicking

7.1. Mosaicking in a nutshell

Several excellent descriptions of the mosaicking imaging and deconvolution algorithms can be found (see e.g. Cornwell 1988; Cornwell et al. 1993; Sault et al. 1996b). Here, we summarize the approach implemented in the gildas/mapping software used to image and deconvolve the data from the Plateau de Bure Interferometer. This approach is based on original ideas by F. Viallefond in the early 90s (Gueth et al. 1995).

The basic ideas of nonlinear mosaicking are 1) imaging the different fields of the mosaic independently; 2) linearly adding the single-field dirty images into a dirty mosaic; and 3) jointly deconvolving the dirty mosaic.

7.1.1. Single-field imaging

For simplicity, we skip the gridding convolution in the following equations because the gridding step does not change the nature of the equations. Imaging the fields individually means that we will work at constant αs. We first define the single-field dirty image of the αs-field as (79)where the Fourier transform of the single-field dirty image is the product of the sampling function S(up,αs) and the visibility function V(up,αs): (80)From the previous equations, it is easily demonstrated that (81)where the single-field dirty beam is defined as (82)We can rewrite the previous equation as (83)meaning that the single-field dirty images can be written as a local convolution of BαsI and Δαs, the single-field dirty beam associated to the currently imaged field.

7.1.2. Mosaicking the dirty images

In gildas/mapping7, the single-field dirty images are formed on the same grid (in particular the same pixel size and the same image size covering about twice the mosaic field of view). These single-field dirty images are then linearly averaged as (84)where (85)and w(αs) is the sky plane weighting function, i.e., (86)In the previous equation, the αi are the positions of each sky-plane measurement, and σi is the rms noise associated with . Cornwell et al. (1993) demonstrates that the noise in the mosaic image naturally varies across the field as (87)In particular, it rises sharply at the edges of the mosaic.

7.1.3. Joint deconvolution

Standard algorithms of single-field deconvolution must be adapted to the mosaicking case because both the dirty beam and the noise vary across the mosaic field of view. We describe here the adaptations made in gildas/mapping of the simplest CLEAN deconvolution method, described in Högbom (1974). Adaptations of more evolved CLEAN deconvolution methods are also implemented following the same basic rules.

  • 1.

    We first initialize the residual and signal-to-noise maps from the dirty and noise maps(88)and (89)

  • 2.

    The kth CLEAN component is sought on the SNRk − 1 map instead of the Rk − 1 map to ensure that noise peaks at the edges of the mosaic are not confused with the true signal of the same magnitude.

  • 3.

    Using that the kth CLEAN component is a point source of intensity Ik at position αk, the residual and signal-to-noise maps are then upgraded as follows: (90)and (91)Here γ(~0.2) is the usual loop gain that ensures convergence of the CLEAN algorithms.

  • 4.

    Steps 2 and 3 are iterated as long as the stopping criterion is not met.

7.1.4. Wide-field measurement equation

To help the comparison between mosaicking and wide-field synthesis, we now go one step further than is usually done in the description of mosaicking; i.e., we write the image-plane measurement equation as a wide-field measurement equation of the same kind as Eq. (23). Substituting Eq. (81) into Eq. (84) and reordering the terms after inverting the order of the sum over αs and αp, one obtains (92)with (93)Taking the inverse Fourier transforms of Dmos, we get the mosaicking transfer function (94)with (95)

7.2. Comparison

While both mosaicking and wide-field synthesis produce image-plane measurement equations of the same kind (see Eqs. (23) and (92)), the comparison of the dirty beams (Eqs. (57) and (93)) and of the transfer functions (Eqs. (35) and (94)) immediately shows the different dependencies on the primary beams (B), the single-field dirty beams (Δ), the image-plane weighting functions (Ω), and their respective Fourier transforms (, and W). This means that mosaicking is not mathematically equivalent to wide-field synthesis, though both methods recover the sky brightness. These differences come directly from the differences in the processing. If we momentarily forget the gridding steps, mosaicking starts with a Fourier transform along the up dimension of the visibility function, and most of the processing thus happens in the sky plane, while wide-field synthesis starts with a Fourier transform along the αs dimension, and most of the processing thus happens in the uv plane.

Moreover, both methods are irreducible to each other. Wide-field synthesis gives a more complex dirty beam formulation in the image plane, which could give the impression that it is a generalization of mosaicking. Indeed, the wide-field image-plane weighting function can be chosen as the product of a Dirac function of α′ times a function ω of α′′, 91This implies a wide-field uv-plane weighting function independent of u′; i.e., W(u′, u′′) = (u′′). This choice is a clear limitation because it enables us to influence the transfer function only locally (around each measured up spatial frequency), while weighting is generally intended to globally influence the transfer function (see Sect. 5). Eitherway, in this case, the wide-field dirty beam can easily be simplified to (97)While this simplified formulation of the wide-field dirty beam is closer to the mosaicking formulation, they still differ in a major way: ω(α′′– αs) is a shift-invariant function contrary to Ωmos(α′′, αs). This is the shift-dependent property of Ωmos(α′′, αs), which implies the additional complexity (integral over us in addition to the integral over up) of the mosaicking transfer function (Eq. (94)) over the wide-field one (Eq. (35)).

One main difference between the two processing methods is that standard mosaicking prescribes a precise weighting function, while we argue that the wide-field weighting function should be defined according to the context (see Sect. 5). Another important difference is the treatment of the short spacings, which are naturally processed in the wide-field synthesis methods, but which needs a very specific treatment in mosaicking (see Sect. 6 and references therein). Finally, while mosaicking implies a gridding only of up dimension of the measured visibilities, wide-field synthesis naturally requires a gridding of both the up and αs dimensions. As the Nyquist sampling along the αs dimension is only 0.5/dprim, the gridding of the sky plane can result in a large reduction of the data storage space and cpu processing cost when processing on-the-fly and/or multi-beam observations.

8. Summary

Interferometric wide-field imaging implies scanning the sky in one way or another (e.g. stop-and-go mosaicking, on-the-fly scanning, sampling of the focal plane by multi-beams). This produces sampled visibilities SV, which depends both on the uv-plane and sky coordinates (e.g., up and αs).

Based on a basic idea by Ekers & Rots (1979), we proposed a new way to image the interferometric wide-field sampled visibilities: SV(up, αs). After gridding the measured visibilities both in the uv and sky planes, the gridded visibilities are Fourier-transformed along the αs sky dimension, yielding synthesized visibilities sampled on a uv grid whose cell size is related to the total field of view; i.e., it is much finer than the diameter of the interferometer antennas. We thus proposed calling this processing scheme “wide-field synthesis”.

The Fourier transform is performed for each constant up value. As many independent estimates of the uv plane are produced as independent values of up measured. A shift-and-average operator is then used to build a final, wide-field uv plane, which translates into a wide-field dirty image after inverse Fourier transform, i.e., (98)where W is a normalized weighting function. Using these tools, we demonstrated that:

  • 1.

    The dirty image () is a convolution of the sky brightness distribution (I) with a set of wide-field dirty beams () varying with the sky coordinate α, i.e., (99)Compared to single-field imaging, the dependency on the primary beam is transferred from a product of the sky brightness distribution into the definition of the set of wide-field dirty beams.

  • 2.

    The set of gridded dirty beams () can be computed from the ungridded sampling function (S), the transfer function (, the inverse Fourier transform of the primary beam), and the gridding convolution kernel (see Eqs. (42), (50) and (51)).

  • 3.

    The dependence of the wide-field dirty beams on the sky position is slowly-varying, with their shape varying on an angular scale typically larger than or equal to the primary beamwidth.

Adaptations of the existing deconvolution algorithms should be straightforward.

A comparison with standard nonlinear mosaicking shows that it is not mathematically equivalent to the wide-field synthesis proposed here, though both methods do recover the sky brightness. The main advantages of wide-field synthesis over standard nonlinear mosaicking are

  • 1.

    Weighting is at the heart of the wide-field synthesis because it isan essential part of the shift-and-average operation. Indeed, notonly can a multiplicative weight be attributed to each visibilitysample before any processing, but the uv-plane weighting function (W, see Eq. (98)) is also a degree of freedom, which should be set according to the conditions of the observation and the imaging goals, e.g. highest signal-to-noise ratio, highest resolution, or most uniform resolution over the field of view. The W weighting function thus enables us to modify the wide-field response of the instrument. On the other hand, mosaicking requires a precise weighting function in the image plane, which freezes the wide-field response of the interferometer.

  • 2.

    Wide-field synthesis naturally processes the short spacings from both single-dish antennas and heterogeneous arrays along with the long spacings. Both of them can then be jointly deconvolved.

  • 3.

    The gridding of the sky plane dimension of the measured visibilities, required by the wide-field synthesis, may potentially save large amounts of hard-disk space and cpu processing power relative to mosaicking when handling data sets acquired with the on-the-fly observing mode. Wide-field synthesis could thus be particularly well suited to process on-the-fly observations.

The wide-field synthesis algorithm is compatible with the uvw-unfaceting technique devised by Sault et al. (1996a) to deal with the celestial projection effect, known as non-coplanar baselines (see Appendix B). Finally, on-the-fly observations imply an elongation of the primary beam along the scanning direction. These effects can be decreased by an increase in the primary beam sampling rate. However, it may limit the dynamic range of the image brightness if the primary beam sampling rate is too coarse (see Appendix C).


1

In the rest of this paper, stop-and-go mosaicking refer to the observing mode, while mosaicking alone refer to the imaging family.

2

The convolution theorem, which states that the Fourier transform of the convolution of two functions is the product of the Fourier transform of both individual functions, is a special case for Eq. (15): it can be recovered by setting . Indeed, as already mentioned in the introduction, the ideal measurement Eq. (1) can be interpreted as a convolution with an additional phase term. By Fourier transforming along the αs dimension, the convolution translates into a product of Fourier transforms and , while the phase term translates into a shift of coordinates: .

3

I.e. direct Fourier transform along the αs dimension, shift-and-average to define a final wide-field plane, and inverse Fourier transform.

4

Loosely speaking, the sampling function can be thought as a function whose value is 1 where there is a measure and 0 elsewhere.

5

The noise is assumed to have a Gaussian probability distribution function.

6

The notion of averaging linepath has been introduced in Sect. 3.1 (see in particular Eq. (16)).

7

See http://www.iram.fr/IRAMFR/GILDAS for more information about the GILDAS software.

8

In contrast to the convention used in this paper, the dmax and dfield unit is meter instead of unit of λ in the second form of the criterion, in order to explicitely show the dependence on the wavelength.

Acknowledgments

This work has mainly been funded by the European FP6 “ALMA enhancement” grant. This work was also funded by grant ANR-09-BLAN-0231-01 from the French Agence Nationale de la Recherche as part of the SCHISM project. The authors thank F. Gueth for the management of the on-the-fly working package of the “ALMA enhancement” project. They also thank S. Guilloteau, R. Lucas and J. Uson for useful comments at various stages of the manuscript and D. Downes for editing their English. They finally thank the referee, B. Sault, for his insightful comments, which challenged us to try to write a better paper.

References

Appendix A: Demonstrations

A.1. Ekers & Rots scheme

Fourier-transforming the visibility function along the αs dimension at constant up, we derive with simple replacements We then use the following change of variables β ≡ αp − αs and dβ =  −dαs, to get

A.2. Incomplete sampling

We here demonstrate that Eqs. (23) and (34) are equivalent. To do this, we take the direct Fourier transform of Idirty(α) (A.6)and we replace I(α′) by its formulation as a function of its Fourier transform (A.7)We thus derive (A.8)Using the following change of variables α′′ ≡ α − α′, α′ = α − α′′ and dα′′ =  −dα′, the innermost integral can be written as In the last two steps, we have simply recognized two different steps of Fourier transforms of D. Finally, (A.12)

A.3. Gridding

The gridding kernel can be defined as the product of two functions, each one operating in its own dimension. We use this to study separately the effect of gridding in the uv and sky planes. We then use the intermediate results to get the effect of gridding simultaneously in both planes.

A.4. Gridding in the uv plane

We define the sampled visibility function gridded in the uv plane as (A.13)Using that the gridding is here applied on the up dimension, while the Fourier transform is applied on the αs dimension, it is easy to show that the gridding and Fourier-transform operations commute: Defining the Fourier transform of the uv gridded dirty image, we derive Using Eq. (31) to replace , we can write the Fourier transform of the uv gridded dirty image as (A.18)with (A.19)and (A.20)Using (A.21)and (A.22)we derive (A.23)or (A.24)Thus, is the uv gridded version of the generalized sampling function .

A.4.1. Gridding in the sky plane

We define the sampled visibility function gridded in the sky plane as Applying the convolution theorem on the Fourier transform along the αs dimension, we derive (A.27)Defining the Fourier transform of the sky-gridded dirty image, we derive Using Eq. (31) to replace , we can write the Fourier transform of the sky-gridded dirty image as (A.30)with (A.31)and (A.32)or, with the definition of (i.e., Eq. (45)), (A.33)Using (A.34)and the convolution theorem when taking the inverse Fourier transform of , we derive (A.35)Thus, is the sky gridded version of the generalized sampling function .

A.5. Gridding in both planes

Starting from the definition of (Eq. (41)), we Fourier-transform it along the sky dimension at constant up. Using that the gridding along the up dimension can be factored out of the Fourier transform, we derive (A.36)Using Eq. (A.27), we now replace in the previous equation to get (A.37)or (A.38)From this relation, it is easy to deduce that (A.39)Using the convolution theorem when taking the inverse Fourier transform of along the us dimension and replacing with Eq. (A.24), we finally derive (A.40)

A.6. Wide-field vs. single-field dirty beams

The notation (59) yields . Using this in Eq. (35) gives (A.41)Taking the inverse Fourier transform along the u′′ axis of Eq. (A.41) and reordering the integral to factor out the term independent of u′′, we can write (A.42)where (A.43)We now introduce the following definition (A.44)to derive Using the following change of variables v ≡ u′′ + u′ − up, u′′ = v − u′ + up and dv = du′′ on the innermost integral, we get Substituting this result into Eq. (A.42) and taking the inverse Fourier transform along the u′ axis, we can write (A.45)where Using the following change of variables v ≡ up − u′, u′ = up − v and dv = du′, we get (A.46)Substituting this result into Eq. (A.45) and re-ordering the terms, we can write (A.47)where

A simple application of the convolution theorem gives

where (A.48)Substituting this result into Eq. (A.47), we finally derive the desired expression, i.e., Eq. (57).

Appendix B: From the celestial sphere onto a single tangent plane

Equation (1) neglects projection effects, known as non-coplanar baselines. Any method which deals with interferometric wide-field imaging must take this problem into account. After a short introduction to the problem, we show how wide-field synthesis is compatible with at least one method, namely the uvw-unfaceting of Sault et al. (1996b). This method tries to build a final wide-field uv plane from different pieces, just as our wide-field synthesis approach does. Another promising method is the w-projection, based on original ideas of Frater & Docherty (1980) and first successfully implemented by Cornwell et al. (2008). We did not look yet at its compatibility with wide-field synthesis.

B.1. w-axis distortion

When projection effects are taken into account, the measurement equation reads (B.1)In this equation, we continue to work in 1 dimension for the sky cosine direction (αp), but we explicitly introduce the dependence along the direction perpendicular to the sky plane. This dependence appears in two ways, which is handled in very different ways. First, the factor can be absorbed into a generalized sky brightness function (B.2)After imaging and deconvolution, I(αp) can be easily restored from the deconvolved ℐ(αp) image. The second dependence appears as an additional phase, which is written as (B.3)Thompson et al. (1986, Chap. 4) shows that this additional phase can be neglected only if8(B.4)The first form of the criterion indicates that the approximation gets worse at high spatial dynamic range (i.e., θfield/θsyn ≪ 1) while the second form indicates that the approximation gets worse at long wavelengths.

B.2. uvw-unfaceting

For stop-and-go mosaicking, it is usual to delay-track at the center of the primary beam for each pointing/field of the mosaic. This phase center is also the natural center of projection of each pointing/field. Stop-and-go mosaicking thus naturally paves the celestial sphere with as many tangent planes as there are pointings/fields; i.e., this observing scheme is somehow enforcing a uvw-faceting scheme. In the framework of on-the-fly observations with ALMA, D’Addario & Emerson (2000) indicate that the phase center will be modified between each on-the-fly scan while it will stay constant during each on-the-fly scan. This is a compromise between loss of coherence and technical possibilities of the phase-locked loop. Using this hypothesis, the maximum sky area covered by the on-the-fly scan must take into account the maximum tolerable w-axis distortion.

The easiest way to deal with such data is to image each pointing/field around its phase center and then to reproject this image onto the mosaic tangent plane as displayed in Fig. 5 of Sault et al. (1996b). These authors point out that this scheme implies a typical w-axis distortion ϵ less than (B.5)where θcenter is the angle from the pointing/field center and θalias is the anti-aliasing scale defined in Sect. 4.2. In particular, ϵ is 0 at the phase center of each pointing/field. In other words, this scheme limits the magnitude of the w-axis distortion to its magnitude on a size equal to the anti-aliasing scale (i.e., a few time the primary beamwidth) instead of a size equal to the total mosaic field of view. This scheme thus solves the projection effect as long as the w-axis distortion is negligible at sizes smaller than or equal to the anti-aliasing scale. A natural name for this processing scheme is uvw-unfaceting because it is the combination of a faceting observing mode (i.e., regular change of phase center) and a linear transform of the uv coordinates to derive a single sine projection for the whole field of view.

Sault et al. (1996b) also demonstrate that the reprojection may be done much more easily and quickly in the uvw space before imaging the visibilities because it is then just a simple transformation of the uv coordinates, followed by a multiplication of the visibilities by a phase term. Finally, Sault et al. (1996b) note that it is the linear character of this uv coordinate transform which preserves the measurement Eq. (1). As the change of coordinates happens before any other processing, it also conserves all the equations derived in the previous sections to implement the wide-field synthesis.

Appendix C: On-the-fly observing mode and effective primary beam

Usual interferometric observing modes (including stop-and-go mosaicking) implies that the interferometer antennas observe a fixed point of the sky during the integration time. Conversely, the on-the-fly observing mode implies that the antennas slew on the sky during the integration time. This implies that the measurement Eq. (1) must be written as (Holdaway & Foster 1994; Rodríguez-Fernández et al. 2009): (C.1)where δt is the integration time and ûp and are the mean spatial frequency and direction cosine, defined as (C.2)In this section, we analyze the consequences of the antenna slewing on the accuracy of the wide-field synthesis.

C.1. Time averaging

In all interferometric observing modes, it is usual to adjust the integration time so that up(t) can be approximated as ûp. To do this, it is enough to ensure that up(t) always varies less than the uv distance associated with tolerable aliasing (dalias, see Sect. 4.2) during the integration time (δt) (C.3)where dmax is the maximum baseline length, ωearth is the angular velocity of a spatial frequency due to the Earth rotation (7.27 × 10-5 rad   s-1), θalias and θsyn are respectively the minimum field of view giving a tolerable aliasing and the synthesized beam angular values.

Table C.1

Definition of the symbols used to explore the influence of on-the-fly scanning on the measurement equation.

C.2. Effective primary beam

Assuming that condition (C.3) is ensured, we can write Eq. (C.1) with the same form as Eq. (1) by the introduction of an effective primary beam (Beff); i.e., (C.4)where (C.5)Using the following change of variables (C.6)we derive (C.7)with (C.8)and (C.9)In these equations, vslew(β) is the slew angular velocity of the telescope as a function of the sky position, δαs is the angular distance covered during δt, A is an apodizing function, and Π(β) is the usual rectangle function, which reproduces the finite character of the time integration.

C.3. Interpretation

The form of the measurement equation is conserved when averaging the visibility function over a finite integration time, as long as the true primary beam is replaced by an effective primary beam, which is the convolution of the true primary beam by an apodizing function. To go further, it is important to return to the two dimensional case. Indeed, the convolution must be done along the slewing direction, resulting in an effective primary beam elongated in a particular direction.

thumbnail Fig. C.1

Assessement of the relative error implied by the use of the true primary beam instead of the effective primary beam when analyzing interferometric on-the-fly data sets. Left: inverse Fourier transform of interferometer primary beam, (i.e. the autocorrelation of the antenna illumination). Right: relative error as a function of sampling rate of the primary beam. The curves of different colors show the results at different normalized uv distances (u/dprim) from the center of .

Open with DEXTER

In principle, the equations derived in Sect. 3 can be accommodated just by replacing the true primary beam by its effective associate. In practice, the probability to take into account the effective primary beam is low because its shape varies with time. Indeed, it is often assumed that the sky is slewed along a straight line at constant angular velocity. Even in this simplest case, it is advisable to slew along at least two perpendicular directions to average systematic errors, implying two different effective primary beams. However, practical reasons may/will lead to complex scanning patterns: 1) the limitation of the acceleration when trying to image a square region leads to spiral or Lissajous scanning patterns; 2) the probable absence of derotators in future multi-beam receivers (B. Lazareff, private communication) implies the need to take into account the Earth rotation in the scanning patterns of the off-axis pixels.

C.4. Approximation accuracy

In the following, we thus ask what is the trade-off accuracy of using the true primary beam instead of the effective primary beam. The first point to mention is that using different scanning patterns somehow helps because the averaging process then makes the bias less systematic. Following Holdaway & Foster (1994), we quantify the accuracy lost in the Fourier plane. Indeed, the Ekers & Rots scheme tries to estimate missing sky brightness Fourier components from their measurements apodized by the Fourier transform of the primary beam. In the Fourier space, the above convolution just translates into a product. The Fourier transform of the apodizing function thus degrades the sensitivity of the measured visibility, V(up, αs), to spatial frequencies at the edges of the interval  [up − dprim,up + dprim] . To guide us in our quantification of the accuracy lost, we now explore the simplest case of linear scanning at constant velocity, where vslew(β) is constant and δαs = vslew   δt. The Fourier transform of the apodizing function is then a sinc function: (C.10)The relative error implied by the use of the true primary beam instead of the effective primary beam is then (C.11)Figure C.1 shows this relative error as a function of the number of samples per primary beam FWHM in the image plane (i.e., θfwhm/δαs) for different uv distances (in units of dprim). We see that we derive a 1% accuracy at all u when we sample the image plane at a rate of 5 dumps per primary beam. However getting a 0.1% accuracy needs quite high sampling rates (about 15). This must be compared with the accuracy of knowledge of .

We note that if a better accuracy is needed than the one achievable with the highest sampling rate, it is in theory possible to replace in the correlator software the rectangle apodizing function by another function which falls more smoothly. To avoid the loss of sensitivity inherent to the use of such an apodizing function (by throwing away data at the edges of the time interval of integration), would require, for instance, to half-overlap the integration intervals. This would imply more book-keeping in the correlator software and some noise correlation between the measured visibilities.

All Tables

Table 1

Definition of the symbols used to expose the wide-field synthesis formalism.

Table 2

Definition of the uv and sky scales relevant to wide-field interferometric imaging.

Table 3

Interval ranges of definition and associated sampling rates for the used functions.

Table 4

Minimum sizes of the dirty beam images to get an image fidelity or a dynamic range greater than a given value.

Table 5

Definition of the symbols used to expose the processing of the short spacings.

Table C.1

Definition of the symbols used to explore the influence of on-the-fly scanning on the measurement equation.

All Figures

thumbnail Fig. 1

Visualization of the different angular scales relevant to wide-field interferometric imaging. The notion of anti-aliasing scale (θalias) is introduced and discussed in Sect. 4.2.

Open with DEXTER
In the text
thumbnail Fig. 2

Illustration of the principles of wide-field synthesis, which enables us to image wide-field interferometric observations. The top row displays the sky plane. The middle row displays the 4-dimensional visibility space and the bottom row shows 2-dimensional cuts of this space at different stages of the processing. In panels b) to d), the scanned dimensions (αs and us) are displayed in blue while the phased spatial scale dimensions (up) are displayed in red and the spatial scale dimensions (u) of the final wide-field uv plane are displayed in black. The grey zones of panels b.2) and c.2) show the regions of the visibility space without measurements (missing short-spacings). In detail, panel a) shows a possible scanning strategy of the sky to measure the unknown brightness distribution at high angular resolution: for simplicity it is here just a 7-field mosaic. Panels b.1) and b.2) sketch the space of measured visibilities: the uv plane at each of the 7 measured sky positions is displayed as a blue square box in panel b.1) and a blue vertical line in panel b.2). For simplicity, only 6 visibilities are plotted in panel b.1). Panels c.1) and c.2) sketch the space of synthesized visibilities after Fourier transform of the measured visibilities along the scanned coordinate (αs): at each measured spatial frequency up (displayed on the blue axes) is associated one space of synthesized wide-field spatial frequencies displayed as one of the red squares in panel c.1) and the red vertical lines in panel c.2). The wide-field spatial scales are synthesized 1) on a grid whose cell size is related to the total field of view of the observation and 2) only inside circles whose radius is the primary diameter of the interferometer antennas. Panels d.1) and d.2) display the final, wide-field uv plane. This plane is built by application of the shift-and-average operator along the black lines on panel c.2), lines that display the region of constant u spatial frequency in the (up, us) space. Standard inverse Fourier transform and deconvolution methods then produce a wide-field distribution of sky brightnesses as shown in panel e).

Open with DEXTER
In the text
thumbnail Fig. 3

Simple models of the antenna power patterns as a function of the sky angle in units of half the primary beam FWHM (θfwhm). In the 3 cases shown, the illumination is Gaussian with an edge taper of 12.5 dB but 3 different ratios of the secondary-to-primary diameters (i.e. fb, the antenna blockage factors) are considered (see e.g. Goldsmith 1998, Chap. 6). The middle and bottom panels respectively model ALMA and PdBI antennas. The red lines define the minimum angular sizes for which the antenna power pattern is less than a given fraction.

Open with DEXTER
In the text
thumbnail Fig. 4

Length of the averaging linepaths displayed as black lines in panel c.2) of Fig. 2, as a function of the spatial scale in the final, wide-field uv plane. In the case of a continuous sampling of up between dmin and dmax, these quantities can be interpreted as the number of measures that contribute to the estimate of (u).

Open with DEXTER
In the text
thumbnail Fig. 5

Sketches of the natural weighting of the synthesized wide-field visibilities. Each measured spatial frequency will produce wide-field spatial frequencies apodized by the transfer function () centered on the measured spatial frequency. The used transfer function depends on the telescopes used, explaining why wide-synthesis naturally handles the short spacing either from a single-dish antenna or from a heterogeneous array. The synthesized visibilities in the overlapping regions will then be averaged. Two textbook examples are illustrated: 1) the combination of data from the IRAM-30 m single-dish (red transfer function) and from the Plateau de Bure Interferometer (black transfer functions) at the top; and 2) the combination of data from ALMA 12 m-antennas used either in single-dish mode (red transfer function), in interferometric mode (black transfer functions) and of data from the ACA 7 m-antennas (blue transfer functions) at the bottom. The minimum uv distances measured by each interferometer were set from the minimum possible distance between antennas (24 m for PdBI, 15 m for ALMA and 9 m for ACA).

Open with DEXTER
In the text
thumbnail Fig. C.1

Assessement of the relative error implied by the use of the true primary beam instead of the effective primary beam when analyzing interferometric on-the-fly data sets. Left: inverse Fourier transform of interferometer primary beam, (i.e. the autocorrelation of the antenna illumination). Right: relative error as a function of sampling rate of the primary beam. The curves of different colors show the results at different normalized uv distances (u/dprim) from the center of .

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.