Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.mrao.cam.ac.uk/~bn204/mk2/publications/2009/ALMAMemo587.pdf Äàòà èçìåíåíèÿ: Wed Nov 25 23:11:25 2009 Äàòà èíäåêñèðîâàíèÿ: Thu Apr 8 11:55:58 2010 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: m 94

ALMA Memo 587, 1­16; (2009)

Printed 19 March 2009

ALMA Memo #587 Inference of Coefficients for Use in Phase Correction I
B. Nikolic
Mullard Radio Astronomy Observatory, Cavendish Laboratory, Cambridge CB3 0HE, UK email: b. nikolic@ mrao. cam. ac. uk http: // www. mrao. cam. ac. uk/ ~ bn204/

5 March 2009

ABSTRACT

We present a Bayesian approach to calculating the coefficients that convert the outputs of ALMA 183 GHz water-vapour radiometers into estimates of path fluctuations which can then be used to correct the observed interferometric visibilities. The key features of the approach are a simple, thin-layer, three-parameter model of the atmosphere; using the absolute measurements from the radiometers to constrain the model; priors to incorporate physical constraints and ancillary information; and a Markov Chain Monte Carlo characterisation of the posterior distribution including full distributions for the phase correction coefficients. The outcomes of the procedure are therefore estimates of the coefficients and their confidence intervals. We illustrate the technique with simulations showing some degeneracies that can arise and the importance of priors in tackling them. We then apply the technique to an hour-long test observation at the Sub-Millimetre Array and find that the technique is stable and that, in this case, its performance is close to optimal. The modelling is described in detail in the appendices and all of the implementation source code is made publicly available under the GPL.
1 INTRODUCTION The performance of millimetre and sub-millimetre wave interferometers is often limited by the fluctuation of the properties of the Earth's troposphere along the lines of sight of each of the elements of the interferometer. If not corrected, these fluctuations lead to a fluctuating delay to each element and subsequent loss of correlation (and therefore sensitivity) and a limit on the maximum usable baseline length. Some recent simulations of the effect of these fluctuations on ALMA science have been presented by us (Nikolic et al. 2008) and other authors (e.g., Asaki et al. 2005). ALMA plans to correct for these fluctuations by a combination of two techniques: (i) Fast-switching, that is interleaving science observations with observations of near-by phase calibrators that allow antenna phase errors to be solved for (ii) Water-vapour radiometry, that is observing the emission of atmospheric water vapour along the line of sight of each element of the array, and inferring and correcting from these observations the fluctuating path to each element. The current plan for ALMA is that fast-switching will be used with a cycle period of between around 10 and 200 seconds while fluctuations on timescales from 1 second up to the fast-switching timescale will be corrected by the water vapour radiometry technique. The actual fast-switching time that will be used will depend on the weather conditions and the scientific requirements of each project as well as the achieved accuracy of phase correction with the water-vapour radiometers. Furthermore we expect to be able to use the WVRs to correct for the phase transfer errors, that is the
c 2009 The Authors

errors that arises due to the different lines of sight to the phase calibration and science sources. One of the key requirements for radiometric phase correction is a prescription for converting the observed sky brightnesses in the neighbourhood of the 183.3 GHz water vapour line, as measured by the WVRs, into a path delay that can be used to correct the observed visibilities. Any such prescription is complicated by the significant pressure-broadening and the high cross section of this line (it can be close to saturated even in the dry conditions of the ALMA site) which means that the optimum conversion strategy is quite a sensitive function of the prevalent atmospheric conditions. In this paper we assume that over the short timescales of interest (i.e., less than approximately 200 seconds) the differential path fluctuation between two telescopes can be predicted by a constant times the differential fluctuation of the sky brightness between the same telescopes. The constant of proportionality is the phase correction "coefficient" and we give in this paper an initial prescription for calculating these coefficients (Section 2). We analyse this prescription with simulations (Section 3) and then apply it to a test observation at the Sub-Millimetre Array (Section 4).

2

METHOD

Our aim is to find estimates of, and confidence intervals on, the coefficients to be used for phase correction. We denote the coefficients as dL/dTB,i where L is the excess path to the telescope and TB,i are the sky brightnesses as measured by the four channels of the WVRs. We use this notation in general although this differential only makes mathematical sense when there is a known and invertible function which connect the cause of fluctuation in L with TB,i .


2
In this paper we assume that the fluctuation in path are caused by fluctuations of the total water vapour column only so that the fluctuation water vapour column is also the sole cause of fluctuations in TB,i . The approach we take in this paper is to construct a physical model of the atmosphere with a number of unknown parameters and use some observables and the basic physical considerations to constrain the possible values the parameters. Once the distribution of possible values of model parameters are known, we can use the same model to compute the distribution of phase correction coefficients. The physical model employed in these initial studies is extremely simple: · We assume water vapour is the only cause of opacity and path fluctuation · We assume all the water is concentrated in a thin layer at a single pressure and temperature · We make the plane parallel approximation for computing effects due to changes in elevation of the antenna This means that the model has only three unknown parameters, namely total water vapour column (c), temperature (T ) and pressure (P). The reason for using such a simple model is that it allows very extensive computational analysis while containing the basic ingredients which we know must influence phase correction. As we gain experience with simulations and observations at the ALMA site, we intend to extend this model in directions which show an actual improvement in the inference of the phase correction coefficients. For the time being we will assume that the only true observables that we have available to constrain the model are the four absolute sky brightness temperatures measured by the WVRs. In additional to the direct observables, we have some constraints on the possible values the model parameters can take; these are called priors. In principle the priors are a joint probability function of all of the model parameters, but in the present paper we simplify by taking them to be separable into a product of functions of one parameter only, i.e., p(cT P) = p(c) p(T ) p(P). The method of solution of this problem that we describe here is the standard Bayesian Markov Chain Monte Carlo approach. As usual, the inference problem is described by the Bayes equation (see, for example Jaynes 2003; Sivia & Skilling 2006): p( |D) = p(D| ) p( ) . p(D) (1)
250

200

150 Tb (K) 100 50 0 175 177.5 180 182.5 (GHz) 185 187.5 190

(a) Precipitable water vapour from 0.6 mm to 1.3 mm
250

200

150 Tb (K) 100 50 0 175 177.5 180 182.5 (GHz) 185 187.5 190

(b) Temperature from 230 to 300 K
250

200

150 Tb (K) 100 50 0 175 177.5 180 182.5 (GHz) 185 187.5 190

where the symbols have following meaning: is the vector of model parameters (in this case {c, T , P}) D is the observed data (in this case the sky brightness temperatures observed by the WVR) p( ) is the prior information (in this case constraints on model parameters as mentioned above) p(D| ) is the likelihood, i.e., the probability of observing the data we have given some model parameters p(D) is the so-called Bayesian evidence, that is a measure of how good our model is in describing the data p( |D) is the posterior distribution of model parameters

(c) Pressure from 400 to 750 mBar Figure 1. Variation of the shape of the 183 GHz water vapour line with changes in water vapour column, temperature and pressure

2.1 Likelihood The computation of the likelihood of an observation given model parameters can be factored into three stages:

(i) Calculation of the sky brightness temperature as a function of frequency (ii) Calculation of the temperatures recorded by the WVRs given a sky temperature. Here we model frequency response of the units and the coupling to the sky. (iii) Calculating the probability of observed data given the model temperatures obtained in the previous step The sky brightness is computed using a simplified model, in which we assume that the only relevant contributors to the atmo2009 ALMA Memo 587


3
250

200

150 Tb (K)

100

50

0 175 177.5 180 182.5 (GHz) 185 187.5 190

(a) Prototype radiometers
250

200

150 Tb (K)

100

50

0 175 177.5 180 182.5 (GHz) 185 187.5 190

(b) Production radiometers Figure 2. Filter design centres and bandwidths of the prototype (top) and production (bottom) radiometers with the 183 GHz water vapour line also shown (in red). The heights of the rectangles representing the filters are inversely proportional to the square root of bandwidth, and are therefore an indication of the relative sensitivity of the filters.

spheric opacity in this band are the water vapour line at 183 GHz and the water-vapour pseudo-continuum. We assume the water vapour line has a Gross line shape and parameters derived from the HITRAN database entry (Rothman et al. 2005) and suitable correction for pressure and temperature as detailed in Appendix A. The parametrisation that we use for the pseudo-continuum follows closely that in the program am by Paine (2004) and is also given in Appendix A. With the opacity calculation, the sky brightness can be calculated using simple radiative transfer. The resulting sky brightness temperature is shown for a range of conditions in Figure 1 illustrating the change in the shape of the water vapour line with pressure, temperature and total water vapour column. The three model parameters are the total water vapour column (c), temperature (T ) and pressure (P). The water vapour column parameter is taken to refer to the zenith column; for the cases when we investigate non-zenith measurements, we assume that this parameter scales according to the plane-parallel approximation while the other parameters remain unchanged with a change in elevation of the telescope. With the sky brightness known, the temperature seen in each of the WVR filters can be calculated. In the present study we assume the WVRs are double-sideband (like the production units on ALMA will be) and that their filter set corresponds to the prototype
2009 ALMA Memo 587

filter set. The reason for using the prototype filter set definition is that later we will make an analysis of a sample test data set collected with the prototype radiometers at the SMA. The prototype filter set is shown in Figure 2, which also shows for contrast the filter set of the production ALMA WVRs. In this study we assume the filters are perfectly sharp and at their design centre frequencies and bandwidths. The average sky brightness across each filter sideband is calculated using the five point Gauss-Legendre Quadrature (e.g., Abramowitz & Stegun 1964) which provides reasonable (but, since it is not an adaptive algorithm, a non-uniform) accuracy and is extremely efficient since the sky brightness needs to be calculated at only five frequencies per filter sideband. The second effect which needs to be taken into account at this stage is the non-perfect coupling of the radiometer in the sky. Our initial analyses of sky-dip measurements during this test campaign (to be published subsequently) suggest that the coupling was about 0.91 and that the termination temperature of the parts of the beam that did not reach the sky was about 265 K. We therefore use these values in the analysis of the data from the SMA later (Section 4); however for the simulations shown in Section 3 for simplicity we assume perfect coupling. The results at stage of computation are the four temperatures, TB,i , seen by the WVR as functions of {c, T , P}. Given these parameters , what is the probability of observing a set of temper atures TB,i ? This is the probability that we denote by the likelihood function L and it is governed by the instrumental effects within the radiometer. For the present analysis, the dominant source of error is the uncertainty in the absolute calibration of the radiometers. The production radiometers specification require this error to be smaller than 2 K at all times. If the underlying errors were normally-distributed (which may be a reasonable approximation) that would require an underlying of less than 1 K. We do not however know at this time precisely what the final distribution will be and further it is likely that the calibration errors will be correlated between the filter outputs. Nevertheless, for simplicity we presently take the errors to be normally and independently distributed so that the likelihood function takes the well known form: log p(D| ) = log L = -


i

TB,i - TB,i

2



T ,i

(2)

where we take T ,i = 1 K . When we better understand the calibration uncertainties in calibration of the WVR units, it will be important to incorporate this information into the above equation. (By comparison the thermal noise in one second of integrating time will be below 0.1 K.)

2.2

Priors

Besides the sky brightness observed by the WVRs, we have (and will have) some constraints on the model parameters that are the results of physical considerations, or derived from independent past observations. For the purposes of this paper we will consider three different priors, shown in Table 1, with the aim of illustrating the effects they have on the final results. The priors are specific to the Mauna Kea site rather than ALMA site since they will also be used for analysis of testing data collected at the Mauna Kea. The conclusions derived from them are expected, however, to transfer directly to the ALMA site.


4
1 2 3 Basic Reasonable Pressure constraint p(c) = p(c) = p(c) = 1 0 1 0 1 0 0 mm < c < 5 mm otherwise 0 mm < c < 5 mm otherwise 0 mm < c < 5 mm otherwise p(T ) = p(T ) = p(T ) = 1 0 1 0 1 0 200 K < T < 320 K otherwise 260 K < T < 280 K otherwise 260 K < T < 280 K otherwise p(P) = p(P) = p(P) = 1 0 1 0 1 0 100 mB < P < 650 mB otherwise. 530 mB < P < 610 mB otherwise. 570 mB < P < 590 mB otherwise.

Table 1. Priors used for the simulations and the analysis of SMA data

The first prior we consider (number 1 in Table 1) is an extremely relaxed prior which puts very loose constraints on the pressure and temperature of the water vapour layer: (i) Uniform probability that the pressure is between 100 and 650 mBar. Since the mean pressure at the peak of Mauna Kea is 605 mBar this prior only assumes that the water vapour is not extremely high in the atmosphere (ii) Uniform probability that the temperature is between 200 and 320 K. Clearly this is a much wider range than typical tropospheric temperatures at altitudes where there is significant water vapour. (iii) Uniform probability that the zenith water vapour column is between 0 and 5 mm. This uninformative constraint is used for all of the other priors too. Since this prior is less informational then the constraints we will have even without any ancillary measurements at the site, it is used to illustrate the degeneracies present if no priors at all are present. The second prior we consider (number 2 in Table 1) is designed to be representative of information on water vapour have might have with basic understanding of the site but without any sophisticated ancillary measurements. In this prior, we assume we know the temperature of the water vapour layer within 20 K and the pressure of the layer to within 80 mBar. This prior is used to analyse the test data from Mauna Kea. Finally, the third prior we consider (number 3 in Table 1) is used as an illustration of the improvement in accuracy that can be obtained with a tight prior on one of the parameters. In this case we still assume that we know the temperature to within 20 K but that we know the pressure of the water vapour layer to within 20 mBar instead of 80 mBar assumed in prior 2. Such a constraint on pressure of the water vapour layer might be derived from, for example, determination of its height using one of the techniques described below.

T = 0.2 K and P = 0.5 mBar. We use the implementation of the MCMC algorithm in the BNMin1 library (Nikolic 2009). By construction of the Metropolis algorithm, the density of points of the chain in a volume of parameter space is an estimate of the posterior probability p( |D). This straight-forward approach does not however allow estimation of p(D) on its own, and so it is not possible to analyse the relative benefits of different models. There are however techniques available for estimation of p(D) and we intend to implement these in the future to allow proper model comparison. Marginalisation and calculating the dL/dTB

2.4

2.3 Markov Chain Monte Carlo With the expressions for the likelihood and the priors that we have described above, we have the necessary information to compute the Bayes equation (Eq. 1). In general, this is computationally expensive (see for example, MacKay 2003) because of the large volume of parameter space that must be characterised in order to determine the maximum of p(D| ) p( ) and the numerical value p(D) = d p(D| ) p( ) . The approach we take is standard Markov Chain Monte Carlo (MCMC) using the Metropolis et al. (1953) algorithm (for a tutorial review, see also Neal 1993). In this approach a chain of points in the parameter space is calculated such that the next point in the chain is found by proposing a new point by random displacement from the current point and calculating the relative likelihood of the two points. If the new point is more likely it is accepted onto the chain; if it is less likely, it accepted with probability determined by the ratio of the likelihoods of the current and new points. The proposal distributions we use in this paper are Gaussian with c = 0.001 mm,

We use the information contained in the Markov Chains in two ways: we marginalise and histogram the points to make estimates of the model parameters; and we, for each point in the chain, calculate the phase correction coefficients dL/dTB,i . The marginalisation of Markov Chains in directions parallel to parameter axes is trivially accomplished by simply ignoring those parameters. Subsequent computation of histograms is also easily done (we use the numpy.histogram routine for this). Computation of the coefficients dL/dTB,i is more involved as new physics must be introduced to compute the delay introduced by the water vapour layer. In the interest of simplicity, we split this calculation into two parts: calculation of the non-dispersive and the dispersive path delay. We compute the non-dispersive path delay using the Smith-Weintraub equation as described in Appendix B1 taking into account the temperature and the pressure of the water vapour. As we are interested here in the effects of water vapour only, we do not here consider the effect a change in temperature will have on the refractive index of the dry air. The dispersive delay calculation is more complicated as well as dependent on the observing frequency (unlike the rest of the discussion presented in this memo). However, at the frequencies of relevance to the SMA test data presented later, the dispersive effect is relatively small, i.e., around 5%. We therefore only calculate an adjustment using the ATM by Pardo et al. (2001) which is then used to scale the non-dispersive path. The details of this calculation are presented in Appendix B2 and the conclusion is that we scale up the non-dispersive phase coefficients by a factor of 1.05 to take into account the dispersion at 230 GHz. With the path delay calculated we calculate the coefficients dL/dTB,i by making a small perturbation to the parameter c, i.e., the quantity of water vapour, and computing the differential as a finite difference.

3

SIMULATION

In this section we present analysis of a simulated single observation of the sky brightness with a WVR. The goals for this section are to illustrate the outputs of the technique we have described above, to show the effects and the importance of the priors, and to illustrate
2009 ALMA Memo 587


5
the approximate accuracy with which it will be possible to predict the correction coefficients dL/dTB,i from the inputs consisting of the sky brightness only. The atmosphere from which we simulate our data point has 1 mm of water vapour toward the zenith at a temperature of 270 K and pressure of 580 mBar. We assume the observation is toward the zenith. As described above, we simulate the sky brightness measured by the filter set of prototype WVRs. We find in this case that the simulated temperatures are: 194.8, 142.6, 90.7 and 47.5 K for the inner to outer channels respectively. These simulated sky temperatures are then used for subsequent inference as the observable temperatures (TB,i ). The first inference we present is with the weak priors from row 1 of Table 1. Recall that in this case we assume weak constraints on the model parameters, including allowing the pressure to be higher than the ground level pressure; hence, this case is an illustration of the inference when essentially no prior information is supplied. The posterior distribution of the model parameters for this case is shown in Figure 3. Like in the other figures of the model parameter posterior distribution, we present these results as the marginalised posterior distributions of each of c, T and P; and as the joint distribution of c vs P, c vs T and c vs T with the remaining (third) parameter marginalised. The marginalised probabilities in the top row of Figure 3 show that in this case we can place relatively weak constraint on the amount of water vapour (to about 5%) and very poor constraints on both the pressure and the temperature of the water vapour. It can also be clearly noted that the posterior distributions are not well approximated by a Gaussian distribution; for example, the posterior distribution p(c) shows a long tail toward higher water columns. The reason for the poor inference can be understood from the second row of Figure 3 which shows the joint distributions of each combination of the model parameters. Considering first the plot of the joint probability p(T P) in the right panel of the lower row, Figure 3, we can see that the retrieval of the pressure and temperature are almost exactly degenerate. In other words if a certain combination of pressure and temperature explain the observed data well, then a higher temperature and a proportionally higher pressure describes the observation also sufficiently well. This degeneracy between pressure and temperature also affects the accuracy of the retrieval of the water vapour column as shown in the left and middle panels of the lower row of Figure 3. There we see that the extreme values of pressure and temperature that are permissible due to the degeneracy give rise to a tail of likelihood to higher values of the retrieved water vapour column. The next result we describe combines the same simulated data point with priors in the second row of Table 1. These are the priors that we might have without significant ancillary information, i.e., that we know the temperature of the water vapour layer to within 20 K and its pressure to within 80 mBar. The posterior distribution from this inference is shown in Figure 4. It can be seen from the top row of this figure that these results are qualitatively different from the inference with very noninformative priors. In this case, the inference of the water vapour is well approximated by Gaussian with a full-width-half-maximum of about 0.012 mm and the entire distribution is within 0.02 mm of the model value. The inferences of the temperature and pressure are still poor however; in fact, it can be seen that their distributions fill almost entirely the space allowed by their priors, indicating that the priors in this case are providing important information. The lower row of Figure 4 again provides an explanation for the marginalised distributions of the model parameters. The
2009 ALMA Memo 587

pressure-temperature joint distribution again shows the degeneracy which explains the poor retrieval of each of those model parameters individually. The water column-pressure and water columntemperature distributions still show spreads but with two important differences: (i) The priors mean that the range over which pressure and temperature can vary is much smaller, therefore leading to smaller errors in the retrieved water column (ii) The joint distributions are highly elongated along the vertical axes, which means a relatively large change in temperate or pressure is required to cause an error in the retrieved water column The condition (ii) is somewhat specific to the simulated point in parameter space, and will not be as true for a general combination of filter centres/widths and observed sky brightnesses. Since the inference shown in Figure 4 is constrained to a reasonable volume of the parameter space we can compute the dL/dTB,i to find how well we can predict the phase correction coefficients. The results are shown in Figure 5 as the marginalised distribution of each of the coefficients (upper part of the figure) and also the joint distribution of each pair of coefficients (lower part). It can be seen in the upper part of this figure that the posterior distributions are non-Gaussian, in fact almost square, and of width of about 10%. This large spread is probably dominated by the uncertainty in the retrieved temperature of the water vapour which influences its refractive index (see Equation B3). The joint distributions plotted in the lower part of Figure 5 show that in this case the errors on inference of the coefficients are very highly correlated, with similar correlation for each pair of the channels. This means that it is unlikely that making use of all of the radiometer channels simultaneously would reduce the error in phase correction due to the uncertainty in the inference of the coefficients. Under different conditions and perhaps with additional observational data, this may be however be possible. The last result that we show in this section is an inference with a tight prior on the pressure (or, equivalently, the height) of the water vapour layer, i.e., 570 mBar < P < 590 mBar, but with the prior on temperature as before. The posterior distribution of the model parameters for this case is shown in Figure 6. We again find that the results of the inference are qualitatively different, primarily in that the posterior distribution of the temperature of the water vapour is now approximately Gaussian and with a full-width-half-maximum of about 8 K, significantly less than its prior range of 20 K. The posterior distribution for the pressure however is, as expected, completely dominated by the prior and simply flat over its prior range. The implication is that a tighter a-priori range of the pressure allows a much better inference of the temperature of the water vapour layer. The posterior distributions of the dL/dTB,i coefficients for this posterior distribution are shown in Figure 7. The qualitatively different inference of the temperature has an important impact on the quality of the inference of these coefficients too, as can be seen from the approximately Gaussian shape of coefficients for channels 1 to 3 and by their significantly tighter ranges. For example. the distribution of coefficient for channel 2 for the inference with prior 2, (Figure 5) is essentially flat with a range of about 0.008 mm/K. With a tight prior on the pressure however, the distribution of this coefficient is close to Gaussian with a full-width-half-maximum of 0.003 mm/K, clearly substantially better.


6
0.1 0.08 0.06
0.01 0.02 0.025 0.02 0.015 f f 0.01 0.005 0.005

0.015

f 0.04 0.02 0 0.95

0

0 200 220 240 T (K)
5000

1

1.05 c (mm)

1.1

1.15

260

280

300

300

400

500 P (mbar)

600

700

5000 600 280 4000 550 P (mbar) P (mbar) 3000 260 3000 550 4000 600

8000

6000

500

T (K)

500

4000

450

2000 240 1000 220

2000

450 2000

400

1000

400

350 0 1 1.025 1.05 c (mm) 1.075 1.1 1.125 0 1 1 1.025 1.05 c (mm) 1.075 1.1 1.125 0 0 1

350 0 220 240 260 T (K) 280 0 1

Figure 3. Posterior distribution of model parameters from a retrieval of simulated measurement at {c = 1.0 mm, T = 270 K, P = 580 mBar} with very weak priors (row 1 of Table 1). The top row shows the marginalised distributions of each of the model parameters, while the bottom row shows the joint distribution of combinations of two parameters with the third marginalised.
0.04
0.0125 0.01 0.0075 0.015 0.0125 0.01 0.0075 0.005 0.0025 0 260 265 270 T (K)
1000 600 590 580 P (mbar) 600 570 560 550 200 540 0 0.98 0.99 1 c (mm) 1.01 1.02 0 1
0.98 0.99 1 c (mm) 1.01 1.02 262.5 0 0 1 T (K) 270 400 267.5 265 277.5 800

0.03

0.02

f

f

0.005

0.01

0.0025

0 0.97

0

0.98

0.99

1 c (mm)

1.01

1.02

1.03

275

280

f

520

540

560 P (mbar)

580

600

620

1250 600 590
600

800

275 272.5

1000

580 P (mbar) 750 570 560 550 250 540 0 265 270 T (K) 275 0 1 500

400

200

Figure 4. Posterior distribution of model parameters derived from four absolute sky brightness temperatures only, i.e., like Figure 3, but now including reasonable flat priors on all three of the c, T and P model parameters.

3.1 Discussion One of the main aims of this section was to illustrate the outputs of an analysis based on the methods described in Section 2, i.e., viewing the problem of estimating the phase correction coefficients as Bayesian inference problem. The main outcome of such an analysis are the posterior distributions for the phase correction coefficients such as presented in Figures 5 and 7. The posterior distributions allow us both to pick a specific coefficient to use for each channel of the radiometer and give us confidence intervals for the accuracy of those coefficients. Obtaining such confidence intervals is important since the combi-

nation of dynamic scheduling, wide range of ALMA configurations and a large number of projects will mean that each science observation is likely to be in conditions which are just marginal for its requirements. Therefore, if phase correction doesn't work as well as expected, it is likely to seriously impact the science aims. We should note however that the confidence intervals calculated using the model we have specified and therefore do not properly capture the probability that simply our model is not very good. Within the Bayesian framework this can be done through the evidence, or p(D), value and we intend to implement this functionality in the near future.
2009 ALMA Memo 587


7
0.02
0.015 0.0125 0.01

0.015

0.01

0.0075 0.005 0.0025

f

0.005

0 0.06

f

0.065
dL d TB,1

0.07 (mm/K)

0.075

0.08

0 0.0625

0.065

0.0675
dL d TB,2

0.07 (mm/K)

0.0725

0.075

0.015 0.0125 0.01 0.0075 0.005

0.02

0.015

0.01

f

f

0.005
0.0025 0 0.084

0.086

0.088
dL d TB,3

0.09 (mm/K)

0.092

0.094

0.096

0 0.145

0.15

0.155
dL d TB,4

0.16

0.165

(mm/K)

(a) Marginalised distributions of the phase correction coefficients
1.5 · 10 0.072
4

0.094

5000

0.16

2000

1.25 · 10

4

0.092 (mm/K)

4000
(mm/K)

0.1575 1500

(mm/K)

0.07

1 · 104 7.5 · 10
3

0.09

3000

0.155 1000 0.1525

dL d TB,2

dL d TB,3

0.068 5 · 103 0.066 2.5 · 10
3

2000 0.088 1000 0.086 0 0.065 0.07
dL d TB,1

dL d TB,4

0.15

500

0.1475

0 0 1 0.065 0.07
dL d TB,1

0

0.075

0.075

0

1

0.065

0.07
dL d TB,1

0.075

0

1

(mm/K) 4000

(mm/K)

(mm/K) 3000

0.094

0.16 1500

0.16

2500 0.1575 2000 (mm/K)

0.092 (mm/K)

3000
(mm/K)

0.1575

0.09

0.155 1000 0.1525 500

0.155 1500 0.1525 1000

dL d TB,3

2000

dL d TB,4

0.088 1000 0.086 0 0.066 0.068
dL d TB,2

0.15

dL d TB,4

0.15 500

0.1475 0

0.1475 0 0 1 0.086 0.088
dL d TB,3

0.07 (mm/K)

0.072

0

1

0.066

0.068
dL d TB,2

0.07 (mm/K)

0.072

0.09 (mm/K)

0.092

0.094

0

1

(b) Join distributions of the phase correction coefficiens Figure 5. Posterior distributions of the coefficients dL/dTB,i used to convert brightness fluctuations into path fluctuations for the model parameters posterior shown in Figure 4.

Beside the distributions of the phase correction coefficients, the outcome of the analysis is also the full joint distribution of the model parameters, which in this case are the water vapour column, and its pressure and temperature. It is the availability of this full posterior distribution that makes reliable estimates of the phase co-

efficients possible. This becomes particularly important when more parameters are added to the problem, as will no doubt be necessary in our case. With more parameters it becomes less and less feasible to pick a single representative point in the parameter space to

2009 ALMA Memo 587


8
0.04
0.03 0.025 0.02
0.0075 0.0125 0.01

0.03

0.02

0.015 0.01

f

f

f 0.005 0.0025

0.01
0.005

0 0.97
590

0

0

0.98

0.99

1 c (mm)

1.01

1.02

1.03

260

265

270 T (K)

275

280
590
2000

570

575

580 P (mbar)

585

590

277.5

800
1500

585

600

275 272.5

585 600 P (mbar) 580 400

P (mbar)

580

400

T (K)

270 267.5

1000

575

200
265 262.5

500

575

200

570 0.98 0.99 1 c (mm) 1.01 1.02

0 0 1
0.98 0.99 1 c (mm) 1.01 1.02

0 0 1

570 265 270 T (K) 275

0 0 1

Figure 6. Model parameter posterior distribution with tight prior on pressure.
0.03 0.025 0.02 0.015 0.01

0.04

0.03

0.02

f

f

0.01
0.005 0 0.0625

0.065

0.0675
dL d TB,1

0.07 (mm/K)

0.0725

0.075

0.0775

0 0.064

0.066

0.068
dL d TB,2

0.07 (mm/K)

0.072

0.074

0.025 0.02 0.015 f 0.01
f

0.015 0.0125 0.01 0.0075 0.005

0.005

0.0025 0 0.149

0 0.086

0.088
dL d TB,3

0.09 (mm/K)

0.092

0.094

0.15

0.151

0.152
dL d TB,4

0.153

0.154

0.155

(mm/K)

Figure 7. The posterior distribution of dL/dTB,i corresponding to Figure 6.

calculate the phase correction coefficients at and making use of the full distribution becomes increasingly important. The model we have been using in this paper is fairly simplistic in that it assumes that the only observables we have are the four absolute sky brightness temperatures observed by the radiometers. What we find is that if we have no further constraints at all, we can estimate the water vapour column to an accuracy of about 5%. We

can not however constrain the temperature of the water vapour at all because it is almost exactly degenerate with the pressure and therefore we can not make any estimate of the phase correction coefficients (since they depend on the temperature, see Equation B3). This shortcoming can be improved on by adding even a fairly loose constraint on the pressure and temperature of the water, such as can be derived from historical distributions of water vapour in
2009 ALMA Memo 587


9
the atmosphere and its temperature. Such loose constraints should be able to provide estimates of the phase correction coefficients in the 10% range. Even tighter a-priori constraints can provide further improvements as shown in Figures 6 and 7. With the model as simple as the one we have presented here and no data from the site of ALMA itself it would be somewhat premature to discuss in detail already the improvements in accuracy particular constraints can make. We however note that we will have a substantial amount of ancillary information such as: · Ground level air temperature, pressure and relative humidity and a number of points at the site · Inference of the vertical temperature profile of the atmosphere at the centre of the array from a commercial O2 line sounder · Library of vertical profiles of atmospheric parameters from past radiosonde launches · Meso-scale meteorological forecasts · Inference of atmospheric parameters from specialised telescope observations such as sky-dips and crossed beam observations Inferences from these measurements can be used as a-priori probabilities on model parameters, or, indeed, in some cases, as further observables which are analysed simultaneously with the sky brightness measurement. It should be noted however that further and better prior information on the temperature and pressure parameters will leave two other sources of uncertainty: errors due to the inaccuracies in the model that we use and limitations in estimating the water vapour column which arises due to calibration accuracy of the WVRs. The best way of tackling these uncertainties may turn out to be to use the observed correlation between changes in the excess path to the telescopes ( L) and the fluctuations of the temperature observed by the radiometers ( TB,i ) as an additional observable that can constrain the models. We will discuss this approach in the next memo in this series. Finally, we now consider ways in which the model presented above may be improved. Firstly, the model may fairly easily be re-parametrised in terms of the height of the water vapour layer, the temperature lapse rate of the atmosphere and the ground level pressures and temperatures instead of the current parametrisation in terms of pressure and temperature of the layer itself. This would allow us to easily take into account the measured groundlevel temperature and pressure and the further information that we will have on the temperature lapse rate (through historical radiosonde measurements and oxygen line profiling) and water vapour height (through historical radiosonde data and specialised telescope scans). The second improvement is to consider the effect of a thick layer of water vapour, perhaps with an exponential fall off in water vapour content as a function of height. Assessing the benefits of such improvements to the model will require computation of the evidence value and, ideally, test data from the ALMA site.
1.25 0.05 1.225 0.04 1.2 c (mm) 0.03 1.175 1.15 1.125 0.02

0.01

0 0 50 100 150 200 250 0 1 time (sample number)

Figure 9. Retrieval of the zenith water vapour column as a function of time from the data recorded by one of the radiometers during the February 17th observation. The priors used in the retrieval are the priors from row 2 of Table 1. The retrieval was calculated once every 25 seconds for the hourlong observation.

4 APPLICATION TO TEST DATA COLLECTED AT SMA In this section we analyse an observation taken with the two prototype ALMA water vapour radiometers (Hills et al. 2001) at the Submillimetre Array (SMA) on Mauna Kea. The observation was made on 17 February 2006 and consisted of an approximately one hour long track of a bright quasar. The interferometric visibility between the two SMA antennas with the water-vapour radiometers was recorded together with the sky
2009 ALMA Memo 587

brightnesses seen by the radiometers. This observation was part of the testing campaign of the prototype radiometers at the SMA which will be described in detail in a separate paper. The effective LO frequency of the observation was 235 GHz and both the upper and the lower sideband were recorded. In the present analysis we use only the upper side band with the on-sky frequency of 240 GHz. The water-vapour radiometers were read out with an integration time of 1 second while the interferometer was read out at 2.5 seconds. In order to match the two data sets, the radiometer data ware re-sampled to 2.5 second time resolution, and a small adjustment to the time-stamps was made to correct for a known timing drift problem. The radiometer data were recorded in already calibrated format and required no further processing. The interferometric data were converted to a text format by M. Reid of the SMA. Subsequently the variation in the observed visibility was transformed to the path fluctuation between the two antennas. The total sky brightness temperatures observed by four channels of two radiometers during this observation are shown in Figure 8. It can be seen that observed brightness is increasing during the observation which is a consequence of the decreasing elevation of the source and therefore increasing airmass as the observation progressed. It can also be seen that the innermost two channels (channels 1 and 2, top row of Figure 8) were almost saturated, indicating significant water vapour along the line of sight. The overall levels of the blue and red curves in this plot allow us to make inference about the atmospheric conditions at the time of the observation while the relative fluctuations between the two curves, once multiplied by the phase correction coefficients, should correlate closely with the path fluctuation measured by the interferometric observation. A retrieval of atmospheric parameters can be made from every one second integration of each of the two radiometers without significant loss in accuracy since, as mentioned before, the thermal noise in one second is much smaller than the expected calibration error of the radiometers. In practice however we expect retrievals will be made rather less often since it is likely that variation in atmospheric conditions will be fairly slow at the ALMA site. In this study we will present retrievals at three different time-resolutions:


10
275 280

270 TB,1 (K) TB,2 (K) 17.25 17.5 17.75

260

265

240

260

220

255 16.75 250

17

18

200 16.75 180

17

17.25

17.5

17.75

18

t (hours UT)

t (hours UT)

225 TB,3 (K) TB,4 (K) 17.25 17.5 17.75

160

200

140

175

120

150 16.75

17

18

100 16.75

17

17.25

17.5

17.75

18

t (hours UT)

t (hours UT)

Figure 8. The sky brightness temperatures observed by the two prototype WVRs while at tests at the SMA, which was tracking a quasar at the time. Blue is one of the radiometers and red is the other. The four panels represent the four radiometer channels.

(i) Calculation of the marginalised zenith water vapour at 300 time points during the observation (ii) A detailed analysis of a single data point in the middle of the observation (iii) Calculation of the marginalised phase correction coefficients at three points of observation We first present a sequence of 300 retrievals from sky brightness temperatures measured at intervals spaced by 25 seconds during the observation. As with the other retrievals in this section, these were made with the priors as shown in row 2 of Table 1. We plotted these retrievals in Figure 9 by marginalising them to obtain the posterior distribution of the zenith water vapour column and plotting these histograms in colour scale as a function of time. Therefore in this plot, time runs in the horizontal direction, the zenith water vapour column parameter along the vertical and the colour represents the posterior probability. The first point to note is that because we retrieve for the zenith water vapour column, we do not expect to see a large increase in the parameter c as the air-mass increases during the observation. The overall range of the water vapour measured in the retrievals is about 5% indicating that at least approximately the plane-parallel approximation holds and the referencing to zenith column is a reasonable approach. Secondly it can be noted that although fluctuations in the water vapour column can be seen, adjacent retrievals (which are separated by 25 seconds in time) generally show similar values indicating a high degree of stability in the retrieved posterior distribution. The magnitude of the fluctuations seen in the retrieved water column is about 50 µm of water vapour which corresponds approximately to a 300 µm of path fluctuation (using the rough con-

version of 1 mm water 6 mm path). As we will present later, e.g., Figure 13, this roughly corresponds to the fluctuations seen in the path by the interferometric measurements. Therefore, over the hour of the observation and a large change in airmass there is no obvious evidence of instability in the water vapour column retrieval. A more detailed analysis of the retrieval from one set of sky brightness measurements in the middle of the observation is presented in Figure 10 in the same format as the previous plots of the posterior distributions of model parameters. As the line of sight water vapour during this observations is significantly higher than the simulations shown in Section 3, some qualitatively differently different results are seen. The first noticeable feature is that in this retrieval the temperature of the water vapour layer is in fact quite well constrained, to a range of around 4 Kelvin, in contrast the results at 1 mm water vapour line of sight column where the posterior distribution filled the entire prior range of 20 K. This is a direct consequence of the almost complete saturation of the innermost channel of the radiometer which means that it is in effect measuring the temperature of the water vapour rather than line of sight column. The saturation of the inner channel, however, now also causes a degeneracy between the pressure and the water vapour column parameter leading to a non-Gaussian posterior distribution for c. The reason for this is that at lower pressure, the line opacity becomes highly peaked, but this can not be detected because of the saturation. The outcome therefore is that at this point in parameter space the water vapour column is somewhat less well constrained and the temperature is better constrained. It should again be noted that the inference errors we present are derived from the model itself and therefore do not account for inaccuracies of the model.
2009 ALMA Memo 587


11
0.05 0.04 0.03
0.02 0.04
0.03 0.025 0.02 0.015 0.01

0.03

f

f

0.02 0.01
0.01
0.005

0 1.15
610

0

f

0

1.175

1.2

1.225 c (mm)

1.25

1.275

1.3

272

274

276 T (K)

278

280
610

525

550

575 P (mbar)

600

625

3000 4000 279 2500 3000 278 2000 P (mbar) 277 T (K)

600

600

2000

590 P (mbar)

590

1500

580 2000 570

580 1000 570

276 275

1500

1000 560 500

560

1000 274 500

550 0 1.175 1.2 1.225 c (mm) 1.25 1.275 0 1

273 0 1.175 1.2 1.225 c (mm) 1.25 1.275 0 1

550 0 274 276 T (K) 278 0 1

Figure 10. Example posterior distribution of model parameters, from a real observation at the SMA at a low elevation (28 degrees) and with priors. Note that c represents the zenith water vapour column.
0.05 0.04 0.03
0.02

0.04

0.03

f

0.02 0.01
0.01

0 0.75

f

0

1

1.25
dL d TB,1

1.5 (mm/K)

1.75

2

0.3

0.325

0.35

0.375
dL d TB,2

0.4

0.425

0.45

(mm/K)

0.04

0.04

0.03

0.03

0.02

0.02

f

0.01

f 0.22 0.23
dL d TB,3

0.01

0 0.21

0.24 (mm/K)

0.25

0.26

0 0.24

0.25
dL d TB,4

0.26 (mm/K)

0.27

0.28

Figure 11. The posterior distribution of dL/dTB,i corresponding to Figure 10.

We next consider the posterior distribution of the phase correction coefficients for this detailed retrieval, as shown in Figure 11. Because of the high line of sight water vapour column, it can be seen that the coefficient for channel 1 is very high, i.e., around 1.1 mm/K. This means that thermal noise of 0.1 K within the radiometer would be sufficient to produce a high random path fluctuation
2009 ALMA Memo 587

of 110 µm. Furthermore, other effects which are not modelled are likely to cause large errors in the path fluctuations calculated from this channel. Therefore, in this case we do not expect this channel to provide useful data for the phase correction itself. It can also be seen in Figure 11 that although the posterior distributions of the phase correction coefficients are reasonable cen-


12
trally peaked, they do all have a tail to higher values. This tail is a result of the non-Guassian inference of the total water vapour column, as seen in the top-left panel of Figure 10. The errors on these coefficients are not dominated by the error on inferred water vapour temperature because, as discussed above, the saturated inner channel in allows for this error to be minimised in this case. Finally it should be noted that the confidence interval of the inferred phase correction coefficients is about 10 to 15% in this case. The final part set of inferences we made was to repeat the prediction of the phase correction coefficients at the beginning and at the end observation so that we can examine how they change during the course of the observation. The results of this analysis are shown in Figure 12, with each retrieval on a separate row. It can be seen in that the values of all of the coefficients increase as the observation progresses. This is of course due to the increase in airmass. The greatest increase is seen in channel 1, which starts at a median value of about 0.45 and finishes at a median value of about 2.25. 4.1 Phase correction We next examine how well we would have been able to do phase correction at the SMA during this observations if we had used the phase correction coefficients retrievals as shown previously. As mentioned previously the SMA interferometer made it possible to infer the fluctuation of path between the two antennas of the array that the WVRs mounted on them. We can also difference the signals between two radiometers to calculate the differential fluctuation of sky brightness as seen from the two antennas. Those two signals are compared (for each of the four channels of the WVRs) in Figure 13. Of course we hope that when multiplied by the phase correction coefficient, the sky brightness fluctuations will look very similar to the path fluctuation. The comparison shown in Figure 13 can immediately confirm that channel 1 is unlikely to yield a useful phase correction, while channel 2 show some correlation with path, and channels 3 and 4 show very high correlation. We next quantify the performance that the phase correction technique would have produced based on the phase correction coefficients inferred in this section. For this we consider both the observation as a whole and three sub-sets of the observation, namely the first, the middle and the last thirds. The reason for considering the subsets is to investigate if the change in the phase correction coefficients inferred during the observation (Figure 12) is reflected in improved correction performance. Before processing the data further we remove from both the interferometric path and the sky temperature measurements a threeminute running mean. This trend removal should approximately model the effect of the further offset calibration scheme that ALMA will use in practice. Although we do not present them here, the data without the running mean removal lead to generally the same conclusions as below. We next asses the potential performance of the phase correction. For these tests we consider phase correction using the fluctuations observed in only one channel at the time. In principle, using a combinations of channels rather than a single one should give better performance because of the reduced effect of thermal noise and possibly the averaging out of errors in the inferred phase correction coefficients. We do pursue this multiple-channel approach further in this paper but we expect it will be used in practice with ALMA. We calculated the corrected path fluctuation using the inferred phase correction coefficients and we also find the `best-fit' phase coefficient, i.e., the coefficient which would have produced smallest corrected path fluctuation given the data that we have. The results are present in Table 2 with following columns: (i) The root-mean-square of the uncorrected path fluctuation as seen by the interferometer, raw (ii) The approximate median inferences for the phase correction coefficients for each channel and for each observation subset dd^L,i TB as shown in Figure 12 (iii) The root-mean-square path fluctuation after correction using radiometric data from each channel, r,i (iv) The phase correction coefficient that would have given the optimum reduction of path fluctuation, Opt( ddL,i ) TB (v) The resultant root-mean-square path fluctuation that would have been obtained if the optimum correction coefficient were used, opt,i As can be seen from the table, channel 1 can not be used in these conditions used for phase correction, as expected from the results already presented. This is evident in phase correction with the `optimum', i.e., best-fit, phase correction coefficient which makes a negligible improvement to the path fluctuation, reducing it from 157 micron to 156 µm. If the phase correction were done with the inferred coefficient of 1.25 mm/K (this was calculated for the middle of the observation) then the path fluctuation would have been drastically increased. This is due to the high effect of thermal noise in this channel when the line is essentially saturated and other effects such as temperature variations of the water vapour. The results for the other three channels show that significant improvements in corrected versus uncorrected path fluctuations can be made based on the radiometer data. Furthermore, it can be seen path corrected with our inferred coefficients dL/dTB,i is close to the path corrected with the best-fit coefficient, i.e., is fairly close to the limit that be achieved with the present data. The difference between using the inferred and best-fit coefficients results in changes of as little as 1 µm rms in path for some of the data sections up to 7 µm path rms in the worst combination of section and filter. Some trends can however be noticed in the data presented in the table. Firstly, excluding channel 1 which is saturated, the inferred phase correction coefficients are systematically underestimating the best-fit coefficients. The only exception is for channel 2 in the last third of the observation, which may however be also affected by saturation. Secondly there is evidence that at lower elevations the gap between correction using inferred and best-fit coefficients widens, suggesting that refinements to the model may be necessary.

4.2

Discussion

The analysis presented in Table 2 shows that the Bayesian inference approach that we described in Section 2 can produce very useful phase correction coefficients even with a very simple atmospheric model. We have shown that the inferred coefficients give a corrected path fluctuation to within about five percent of what the optimal coefficient would have given. There are however a number of interesting points that these data (and the simulations) raise. Firstly, the confidence interval for the inference from our model analysis (Figure 12) show a range of about 10% and furthermore the optimal coefficients found empirically are about 15% higher than the median inferred values. This needs to be compared against the specification for ALMA phase
2009 ALMA Memo 587


13
Channel #1
0.04 0.04 0.03 0.03

Channel #2
0.04 0.03

Channel #3
0.04 0.03

Channel #4

0.02

0.02

0.02

0.02

f

f

f

0.01

0.01

0.01

f 0.165 0.17
dL d TB,3

0.01

0 0.35

0.4

0.45
dL d TB,1

0.5 (mm/K)

0.55

0.6

0.65

0 0.19

0.2

0.21
dL d TB,2

0.22 (mm/K)

0.23

0.24

0 0.16

0.175 (mm/K)

0.18

0.185

0 0.21

0.215

0.22
dL d TB,4

0.225 (mm/K)

0.23

0.235

0.05 0.04 0.03

0.04

0.04

0.04

0.03

0.03

0.03

0.02

0.02

0.02

f

f

f

0.02 0.01
0.01 0.01 0.01

0 0.75

0

1

1.25
dL d TB,1

1.5 (mm/K)

1.75

2

0.3

0.325

0.35

0.375
dL d TB,2

0.4

0.425

0.45

0 0.21

f 0.22 0.23
dL d TB,3

0.24 (mm/K)

0.25

0.26

0 0.24

0.25
dL d TB,4

0.26 (mm/K)

0.27

0.28

(mm/K) 0.04

0.05 0.04 0.03

0.04

0.04

0.03

0.03

0.03

0.02

0.02

0.02

f

f

0.02 0.01
0.01 0.01

f

f

0.01

0 1.75

0

2

2.25

2.5
dL d TB,1

2.75 (mm/K)

3

3.25

3.5

0.5

0.525

0.55

0.575
dL d TB,2

0.6

0.625

0.65

0 0.27

0.28

0.29
dL d TB,3

0.3 (mm/K)

0.31

0.32

0 0.27

0.275

0.28

0.285
dL d TB,4

0.29

0.295

0.3

(mm/K)

(mm/K)

Figure 12. Change in the posterior distribution of the four phase-correction coefficients during the one-hour long scan tracking a quasar. Top row is beg ginning of the scan (UT 17.1), middle row is the middle of the scan (UT 17.5) and bottom row is the end of the scan (UT 17.8).
d^ L dTB,1

UT Range (hours) 17.1­17.9 17.1­17.4 17.4­17.7 17.7­17.9



raw



r,1

Opt(

dL dTB,1

)



opt,1

d^ L dTB,2



r,2

Opt(

dL dTB,2

)



opt,2

(µm) 157 101 197 162

(µm/K) 1260 499 1260 2415

(µm) 420 175 421 783

(µm/K) 64 64 96 20 (a) Part I

(µm) 156 99 195 162

(µm/K) 383 221 383 588

(µm) 95 61 99 119

(µm/K) 404 274 469 552

(µm) 95 59 94 119

UT Range (hours) 17.1­17.9 17.1­17.4 17.4­17.7 17.7­17.9

raw (µm) 157 101 197 162

d^ L dTB,3



r,3

Opt(

dL dTB,3

)



opt,3

d^ L dTB,4



r,4

Opt(

dL dTB,4

)



opt,4

(µm/K) 242 179 242 310

(µm) 74 53 77 86

(µm/K) 286 220 290 403

(µm) 71 51 70 79

(µm/K) 268 234 268 297

(µm) 64 47 64 79

(µm/K) 304 259 303 361

(µm) 61 46 61 75

(b) Part II Table 2. Results of analysis of SMA data from February 17th. The part I of the table shows results for channels 1 and 2; the second part for channels 3 and 4. The UT range column indicates the time range of the data analysed: the first row corresponds to the whole set and the other three rows to respective thirds of d the data. The column raw is the RMS of phase fluctuations without correction. The column dT^L is the best estimate from retrieval of the correction coefficient B,i for the ith channel and for that section in time (see also Figure). The column r,i is the RMS phase fluctuation after correction using the estimated coefficient. dL The column Opt( dT ) is the coefficient which would give the minimum RMS give the radiometer and interferometer data. The column opt,i is that optimum B,i RMS.

correction which is: Lcorrected c 1+ 10 µm + 0.02 â Lraw . 1 mm (3)

rection coefficients to better to 2% accuracy, but the present analysis of the model and the data suggest that we will do significantly worse than that. There are a number of improvements that we are planning that will allow better inference of the phase correction coefficients. One of these is the use of the observed correlation between the change in path and sky temperature which is likely to allow a significant improvement of the inference accuracy and will be presented in one of the next papers in this series. Secondly, it should be noted

If we interpret the first term on the right hand side as the budget for the thermal noise contribution to corrected phase fluctuations then the second term, which is proportional to the raw phase fluctuation, among other effects (e.g., Nikolic et al. 2007) must account for the imperfect phase correction due to errors on the inferred phase correction coefficients. Therefore we need to be able to infer phase cor2009 ALMA Memo 587


14
-2 -2.5 -3 -3 TB,1 (K) -3.5 -4 -5 -4.5 -5 1000 -6 1000 TB,2 (K) L (µ m) 17.25 17.5 17.75 -4 -2

500

500

L (µ m)

0

0

-500

-500

-1000 16.75

17

18

-1000 16.75

17

17.25

17.5

17.75

18

t (hours UT)

t (hours UT)

(a) Channel 1 (innermost channel)
-3 4

(b) Channel 2

-4 TB,3 (K) TB,4 (K) L (µ m) 17.25 17.5 17.75

3

2

-5

1

-6

0 -1 1000

-7 1000

500

500

L (µ m)

0

0

-500

-500

-1000 16.75

17

18

-1000 16.75

17

17.25

17.5

17.75

18

t (hours UT)

t (hours UT)

(c) Channel 3

(d) Channel 4 (outermost channel)

Figure 13. Comparison of the differences between each of the channels (parts a-d) on the two radiometers at SMA test with the fluctuation of path observed by the interferometer.

that in this case one of the limiting factors on the inference is the absolute calibration of the radiometers. Although we do not think that the a-priori calibration of the radiometers when mounted on the ALMA telescopes will be better than what we assumed here, it will be probably the case that the major sources of calibration inac-

curacy will be fixed in time. Therefore we may be able to improve the accuracy of these devices over time by empirical corrections to their calibration. We have not discussed in detail in this paper the dispersive path delay, in particular how it depends on the parameters other
2009 ALMA Memo 587


15
than the water vapour column. This can be investigated by more fully featured atmospheric models along the lines of the simple calculations performed in Appendix B2. We have also not considered the `dry' path fluctuations, that is path fluctuation due to changes in refractive index of dry air. It is not presently clear how significant in practice these will be but data to be collected during the commissioning of ALMA should provide constraints on this. ¨ A. J. Stirling, R. Williamson, V. Belitsky, R. Booth, M. Hagstrom, L. Helldner, M. Pantaleev, L. E. Pettersson, T. R. Hunter, S. Paine, A. Peck, M. A. Reid, A. Schinckel and K. Young. We would like to thank S. Paine for making the code for am publicly available, which allowed us to match our model of the 183 GHz line to his. We would also like to thank J. Pardo for making the code for ATM publicly available. We used ATM for computation of the dispersive phase adjustment. We thank Robert Laing for a careful reading of the manuscript and helpful comments.

5 CONCLUSIONS We have described an approach for calculation of the phase correction coefficients based on: · A very simple, essentially minimal, model of the atmosphere with three parameters: water vapour column, pressure and temperature · Observables which are the four absolute sky brightness temperatures measured by the water vapour radiometers · Priors on the model parameters based either on basic considerations or ancillary measurements · Bayesian inference The approach produces the full probability distributions of the phase correction coefficients making it easier to asses in a timely way the accuracy of phase correction that can be expected. The simulations we presented in Section 3 show that although there are inherent degeneracies in retrieval of atmospheric paramours from sky brightness only, with the appropriate use of prior information fairly good retrievals may be obtained. We have then applied this approach to one test observation taken during the testing of the ALMA prototype radiometers at the SMA. We found that the inferred water vapour column during this observation is stable and its fluctuations are not greater than expected given the observed path fluctuations (Figure 9). In tests of applying the phase correction, we found that good phase corrections can be made by taking representative inferences of the phase correction coefficients and applying them to correct the observed path fluctuation: in fact the typical corrected path fluctuation is within about 5% of the best correction that could be achieved with these data, assuming only one channel is used at a time. In order for it to be possible to meet the very demanding ALMA specifications, it is clear to we will need to make use of further data and possibly more complicated models. The approach we are currently working on is to incorporate information from the occasional phase calibration observations of the observed correlation between path and sky brightness fluctuations into the inference process. The source code for all of the algorithms presented in this paper is available for public download under the GNU Public License at http://www.mrao.cam.ac.uk/~bn204/alma/ memo- infer.html. The observational data from the SMA is also available at the same internet address.

REFERENCES Abramowitz M., Stegun I. A., 1964, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, ninth dover printing, tenth gpo printing edn. Dover, New York Asaki A., Saito M., Kawabe R., Morita K., Tamura Y., Vila-Vilaro B., 2005, Simulation series of a phase calibration scheme with water vapor radiometers for the atacama compact array. ALMA Memo Series 535, The ALMA Project Hills R. E., Gibson H., Richer J. S., Smith H., Belitsky V., Booth R., Urbain D., 2001, Design and development of 183ghz water vapour radiometers. ALMA Memo Series 352, The ALMA Project Jaynes E. T., 2003, Probability Theory : The Logic of Science. Cambridge University Press MacKay D. J. C., 2003, Information Theory, Inference, and Learning Algorithms. Cambridge University Press Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., 1953, The Journal of Chemical Physics, 21, 1087 Neal R. M., 1993, Probabilistic inference using markov chain monte carlo methods. Tech. Rep. CRG-TR-93-1, University of Toronto, http://www.cs.toronto.edu/~radford/ftp/ review.pdf Nikolic B., 2009, BNMin1 Minimisation/Inference library. http: //www.mrao.cam.ac.uk/~bn204/oof/bnmin1.html Nikolic B., Hills R. E., Richer J. S., 2007, Limits on phase correction performance due to differences between astronomical and water-vapour radiometer beams. ALMA Memo Series 573, The ALMA Project Nikolic B., Richer J. S., Hills R. E., 2008, Simulating atmospheric phase errors, phase correction and the impact on alma science. ALMA Memo Series 582, The ALMA Project Paine S., 2004, The am atmospheric model. Tech. rep., SMA Technical Memo, revision 3 Pardo J. R., Cernicharo J., Serabyn E., 2001, IEEE Trans. on Antennas and Propagation, 49 Rothman L., et al., 2005, JQSRT, 96 Sivia D. S., Skilling J., 2006, Data Analysis­A Bayesian Tutorial, 2nd edn. Oxford Science Publications Stirling A., Richer J., Hills R., Lock A., 2008, Turbulence simulations of dry and wet phase fluctuations at chajnantor. ALMA Memo Series 517.1, The ALMA Project, (Original version published 2005)

ACKNOWLEDGEMENTS This work is a part of the ALMA Enhancement programme funded by the European Union Framework Programme 6. The Cambridge effort within this programme has been lead by R. E. Hills and subsequently J. S. Richer. We thank the team involved in construction and the testing of the ALMA prototype water vapour radiometers: P. G. Anathasubramanian, R. E. Hills, K. G.Isaak, M. Owen, J. S. Richer, H. Smith,
2009 ALMA Memo 587


16
APPENDIX A: CALCULATION OF ATMOSPHERIC OPACITY A1 The 183 GHz water vapour line The water line at 183 GHz is assumed to have the Gross line shape: ( ; S, 0 , ) = 4S 2 , 2 - 2 2 + 4 2 2 0 (A1) Of the remaining two terms in Equation B1, the dominant term is the second and it can be transformed to be in terms of water vapour density (see Stirling et al. 2008): dn R dw Mv T (B2)

where is the frequency and the set are {S, 0 , } are the parameters of the line. In the present model we are assuming that water vapour is a relatively small constituent of air and so that the volume mixing ratio for water is negligible. In this case, the line parameters can be computed as follow: 0 = F + P â air nair T = air â P â ref T S = H âe
El (T -Tref ) T Tref

where Mv = 18.02 g mol-1 is the molecular weight of water vapour and w is the density of the water vapour. Assuming a thin layer and by integrating along the line of sight we get: dl R 1741 K dc dc. Mv T T (B3)

B2

Dispersive path delay

(A2) (A3) (A4)

1 - e-h /T Q(Tref ) . 1 - e-h /Tref Q(T )

Here Tref is the reference temperature of the catalogue, which in this case is 296 K; T and P are the temperature and pressure of the air containing the water vapour, and Q(T ) is the partition function of water vapour which we take from a look-up table. The remaining parameters are taken from the HITRAN entry for the 183.3 GHz water vapour line and are shown in Table A1.

The dispersive delay adjustment was calculated using the ATM code Pardo et al. (2001). The calculation was made using a simple program dispersive that takes a frequency and ground-level relative humidity and computes the dispersive and non-dispersive delay due to the water vapour and the dry atmosphere at the supplied frequency. The program dispersive, and the code for the library ATM, are available at http://www.mrao.cam.ac.uk/~bn204/ alma/atmomodel.html, packaged as aatm. All of the code is made available under the GPL license. For this paper we used aatm-0.06, i.e., version 0.06.

A2 Continuum The continuum cross section is : C( ; C0 , T0 , m) = C0 T0 T
m

â N â 2,

(A5)

where N is number density of air, and the parameter set {C0 , T0 , m} describe the continuum. For water, the values we assume are: Param Value Unit C0 T0 m 6.1 â 10-48 300.0 2.6 cm5 /GHz2 K none

APPENDIX B: EXCESS PATH LENGTH B1 Non-dispersive path delay The non-dispersive excess path length due to water is computed using the Smith-Weintraub equation which gives the following expression for the excess refractive index (Stirling et al. 2008): n - 1 = 10-6 Pd Pw Pw + + 2 , T T T (B1)

where Pd is the partial pressure of the dry air and Pw is the partial pressure of the water vapour and the constants have the following values: · = 77.6 â 10-2 K Pa-1 · = 64.8 â 10-2 K Pa-1 · = 3.776 â 103 K2 Pa-1 . We are presently concerned with the excess path introduced with the water vapour only, so we can drop the first term of Equation B1.
2009 ALMA Memo 587


17
Parameter H F air El air self nair Units cm2 â GHz GHz GHz â mbar-1 K GHz â mbar-1 GHz â mbar-1 none Value 2.333884 â 10-21 183.310107 -5.0298 â 10-5 195.908398 2.9114 â 10-3 1.6149 â 10-2 0.77 Comment Line intensity Frequency Air-induced freq. shift Lower state energy level Air broadening Self-broadening Temperature dep. of air broadening

Table A1. The HITRAN parameters and values for the 183.3 GHz line (Rothman et al. 2005).

2009 ALMA Memo 587