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

ALMA Memo 588, 1­15; (2009)

Printed 24 July 2009

ALMA Memo #588 Inference of Coefficients for Use in Phase Correction II: Using the Observed Correlation Between Phase and Sky Brightness Fluctuations
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/

24th July 2009

ABSTRACT

By observing bright and compact astronomical sources while also taking data with the 183 GHz Water Vapour Radiometers, ALMA will be able to measure the `empirical' relationship between fluctuations in the phase of the astronomical signal and the fluctuations of sky brightness around 183 GHz. Simulations of such measurements assuming only thermal noise in the astronomical and WVR receivers are presented and it is shown that accurate determination of the empirical relationship should be possible in a relatively short time. It is then proposed that the best way of using these empirical coefficients is to include them as a constraint on a physical model of the atmosphere ­ this allows them to be used for longer period of time, increasing the efficiency of observing. This approach fits naturally into the analysis framework presented in the previous memo, which has now been extended to implement it. The technique is illustrated via simulations and on a short data set collected at the SMA.
1 INTRODUCTION The Atacama Large (sub-)Millimetre Array (ALMA) will calibrate and correct the effects of the Earth's atmosphere on the phase of the observed radiation by a combination of two techniques: (i) Frequent (every 20­200 s, depending on the conditions and the scientific programme) fast-switching observation of quasars (ii) Water Vapour Radiometry (WVR) using dedicated 183 GHz radiometers to be installed on all 12 m-diameter antennas The first technique requires suspension of science observing for a few seconds but allows a direct observation of the excess path due to the atmosphere between any two antennas. The second technique will be applied continuously, probably even during small changes in observing direction, but requires a model which will translate WVR outputs to estimates of excess path. A model for calculating the coefficients that translate WVR outputs to estimates of path fluctuation can be constructed using the physical properties of the atmosphere and its constituents. For example, the previous memo in this series describes an extremely simple model (Nikolic 2009b) of just a single, thin, layer of water vapour. Whatever the model is that is chosen, there will however be some errors which will limit the accuracy with which the phase correction coefficients are estimated. These errors can be divided in two categories: · Uncertainties in the parameters of the model, e.g., height and temperature of the water vapour in the atmosphere · Incorrect physical assumption or omission of important effects,
c 2009 The Authors

e.g., in a thin water-vapour layer model, it may be the assumption of thinness that introduces a significant error When the commissioning of ALMA begins, it will be possible to asses the accuracy of such models by observing quasars and examining how well the predicted path fluctuations correspond to those measured from the visibilities of the quasar. Clearly, such testing this will allow for better models with most of the relevant physics to be selected. In this memo, I consider a way of taking this a step further and using such test observations to constrain the model parameters as well model section. There are two reasons why this may be practical and useful with ALMA: (i) ALMA will have high sensitivity and has been designed to efficiently do the fast-switching observations. This means that observations of quasars which allow both the interferometric path and WVR sky brightness to be recorded can be done with only a small interruption to science observing. Therefore the relevant `test' data can be obtained cheaply and continuously and should be used to best advantage (ii) Such observations of quasars allow empirical measurement of the relationship between outputs of the WVRs and interferometric phase. Clearly this information must be very relevant for any model which tries to predict this same relationship The main part of the memo begins with an analysis (Section 3) of how well the relationship between fluctuations of WVR outputs and interferometric path can be estimated in practice, with both a single-base line system and a full 50 antenna ALMA array. Then, I will show (Section 4) that such empirical estimates can be used to


2
better constrain the parameters of an atmospheric model, taking in this case the simple model previously proposed. Finally the procedure will be illustrated (Section 5) on some test data collected at the Submillimeter Array with prototype ALMA WVRs.
dL/dT (µ m/K) 600 500 400 300 200 100 0 0 1 2 c (mm) 3 4 5

2 LINEARISATION OF FLUCTUATIONS For the application of phase correction to an array such as ALMA, we are primarily interested in fluctuations (over timescales from one second to about 200 seconds) of path to each of the antenna rather then the absolute value of the path. The reasons for this are: · An interferometer such as ALMA is not sensitive to absolute changes in path, only the relative paths between its antennas · The fast-switching phase calibration can easily remove any phase errors on time scales longer than about 200 s · Instrumental drifts in both the interferometer and the WVRs may become dominant sources of phase error in ALMA timescales longer than about 200 s The fluctuations in path L will be estimated from the fluctuations of the WVR outputs, TB,i according to: L

Figure 1. Model values of the phase correction coefficients, dL/dTB,i , for a wide range of line-of-sight water vapour columns (c), calculated for the ALMA production radiometer. Colours represent the four channels of the WVRs: red is channel 1 (the inner-most channel), yellow-green is channel 2, turquoise is channel 3 and blue channel 4 (the outer-most channel).

(iii) Number and size of antennas in the array In the regime when path errors are small, they are closely related to the estimated flux sensitivity of the calibration observation as follows. The error in the estimated phase from each visibility is approximately S/S where S is the noise in measured flux, S is the calibration source strength and is in radians. The sensitivity of radio interferometers is given, e.g., by Wilson et al. (2009): S = 2M kTsys Ae 2N t (2)


i

dL TB wi dTB,i

,i

(1)

where wi is the weighting of the channels and dL/dTB,i are the `phase-correction coefficients'. The coefficients are predicted by atmospheric models and, when observing a strong source of known structure, can be estimated observationally (as described in the next section). Inherent in this expression is the assumption that L and TB,i are linearly related. This certainly will not be the case for large changes in conditions as the 183 GHz is close to saturation even in the dry conditions prevalent at the ALMA site. However for small fluctuations and short time scales, this should be a reasonable approximation. Throughout this memo, it is assumed that the timescales and baselines are short enough that this linear relationship holds. This means one can consider dL/dTB,i as constant during intervals of interest. If the elevation of the antennas are reasonably constant, it is expected that this will be true on timescale a few minutes or so, which is sufficiently long for the technique described here. If the astronomical source is setting or rising quickly however, this may not be true, and care should be taken to take into account second-order effects.

Where the symbols have following meaning: N Number of correlations, n(n - 1)/2 for a full array of n antenna M Digitisation, correlation, efficiency factor, 1 Tsys System temperature, including atmospheric absorption Ae Effective antenna area, 90 m2 for the 12 m ALMA antennas t Integration time Bandwidth and final sensitivity is for a single polarisation. Since we are interested in measurement of antenna phases, the relevant number of correlations is N = n - 1, i.e., 49 for the case of the main ALMA array. If we take roughly representative parameters for phase calibration at 90 GHz, i.e., M = 1, N = 49,Tsys = 50 K, Ae = 90 m2 , t = 1 s, = 16 GHz (to account for the fact that both polarisations will be measured) we find S = 1.2 mJy. Therefore to achieve 1 degree phase solution accuracy, approximately a 70 mJy calibration source needs to be observed. In the remainder of this paper, we will assume that the strength of the calibration sources are indeed observed and that the resulting path accuracy is 10 µm. Such phase error would lead to only a small contribution to de-correlation even at the shortest operational wavelengths for ALMA. For initial tests of the phase correction techniques at ALMA, it is likely only a single baseline will be available. In this case the error on the phase measurement will be approximately 35 times greater and so sources of around 2.5 Jy should be observed for the results here to be applicable. The intrinsic, thermal-like, errors on the measurements by the WVRs will be in the range of 0.07 K to 0.1 K root-mean-square in one second of integration time, depending on the WVR channel.
2009 ALMA Memo 588

3 HOW ACCURATELY CAN WE ESTIMATE PHASE CORRECTION COEFFICIENTS OBSERVATIONALLY? There are two unavoidable (but also, easy to estimate) sources of error in estimating the phase correction coefficients directly from the correlation between the observed fluctuations in path and sky brightness. These are the thermal-like errors arising in the mixers and front-end amplifiers of both the astronomical and of the WVR receiving systems. The noise in the astronomical system will produce an error in the observed visibility between each pair of antennas, and therefore in the inferred path fluctuation between them. The magnitude of this error will depend on: (i) Receiver noise temperature (ii) Strength of the calibration source


3
1

TB /dL/dT (mm)

0.5

fluctuations, and the assumed noise on path and sky brightness temperatures. We will consider simulations of two types of observations: (i) Calibration observations of that are interleaved with science observations. We assume calibration observation are 1 second long and the intervals between them are in the range of 2 to 200 seconds, with total length of sequence of 1000 seconds (ii) Dedicated calibration observations of length up to 200 seconds Finally, in order to estimate the accuracy with which observations can be used to constrain the phase correction coefficients, we repeat the simulations 100 times to create a Monte-Carlo simulation. 3.2 Finding the best-fitting slope

0

-0.5

-1

-1

-0.75

-0.5

-0.25 l (mm)

0

0.25

0.5

Figure 3. Example of simulated data, with errors in both coordinates. Data are on an approximately 300 m baseline and sampled every one second for 1000 seconds. Error on path is assumed to be 10 µm and error on sky brightness equivalent to 50 µm. The phase screen used is shown in Figure 2.

Although the goal of this memo is to find the constraints on the phase correction coefficients, simple approaches (e.g., Nikolic 2009b) can easily yield estimates which are good enough to translate the above errors in terms of Kelvin into corresponding path errors, at the frequency at which the calibration source is observed. Using the simple single layer model presented previously (Nikolic 2009b), phase correction coefficients as a function of watervapour column only are shown in Figure 1. This figure therefore shows approximately how the thermal errors in the WVRs will translate to path errors for a range of conditions at the ALMA site. In the median conditions at the site, which correspond to water vapour column of 1 mm, the expected error is in the range of less than 10 µm to about 30 µm, depending on the WVR channel. Therefore, in the results presented above, we will do the analysis with assumed errors of 10 µm, 20 µm and 50 µm. There are likely to be other sources in error which have been neglected here, for example: (i) Instrumental drift of interferometer path (ii) Drift in the WVRs sky temperature (specified to be less than 0.1 K in 10 minutes) Although we neglect these in the present work, they may turn out to be a significant limitation in practice. We note however that we limit ourselves to observations around of length of ten minutes or shorter which should restrict the instrumental drifts. 3.1 Simulating fluctuations Beside the intrinsic errors on measurements of the interferometric path fluctuation and the sky brightness temperature, the accuracy with which the correction coefficients can be estimated will also depend on the magnitude of the fluctuations and the length and number of data in the observation. In order to quantify the typical magnitude of fluctuations, we have made simulations using the methods presented in detail in earlier memos (Nikolic et al. 2007, 2008). The phase screen generated for this study is shown in Figure 2 and corresponds to a 100 m thick turbulent layer in the atmosphere. Fluctuations are scaled so that they are 200 µm on a 300 m baseline, i.e., they correspond to the median conditions at the ALMA site. A typical set of simulated data is shown in Figure 3, illustrating the distribution of points due to the statistics of atmospheric
2009 ALMA Memo 588

Since we are assuming the approximate linear relationship between fluctuation in path and the sky brightness, we are only interested in finding the slope of the best fitting line which connects the two quantities. By design of the ALMA WVR system, the errors on the observed path and the observed sky brightness are approximately the same magnitude. As a result, the computation of the best-fitting line is slightly more involved than the usual case in which the only significant errors are in one (`y') coordinate. We adopt the following notation: y and x represent the true ^ ^ values of path and temperature and they are assumed to be related exactly by: y = ax + b ^ ^ (3)

where a is the slope coefficient that we seek to find. The observed values y, x differ from the true values by random errors that are normally distribution with known variance, i.e., y = y + y where ^ y N (0, y ) and similarly for x, x = x + x where x N (0, x ). ^ The values of x and y can in practice estimated a-priori using the calculations from the beginning of this section. A relatively simple, maximum-likelihood, procedure for finding the best fitting line is given by Press et al. (1992). This involves calculating the approximate negative log-likelihood of the observation given model parameters: F (a, b) =
2 i (yi - axi - b) 2 + a2 2 y x

(4)

and then simply minimising function F (a, b) with respect to the two parameters. We implement this using the Levenberg-Marquardt algorithm which produces an estimate of a and of the formal error a although this will only be approximately correct. The implementation is available in BNMin1, version 1.6 (Nikolic 2009a). (A complete solution to the fitting is given by Gull (1989), but has not been used here). The formal error derived from the fitting procedure, a , is used for combining estimates from different baselines in the array (see below), but the final quoted accuracy is always determined by a Monte-Carlo simulations from 100 realisation of the simulated data sets. 3.3 Single baseline

As mentioned previously, in the initial stages of ALMA commissioning and testing of the WVR-based phase correction, it is likely that we will be working with a single baseline configuration. For this the reason, we have performed an analysis of how well the phase


4
North-South (pixels) 200 10 100 0 -10 0 0 200 400 600 800 East-West (pixels) 1000 1200 1400 0 1

Figure 2. A one-eight subsection of the turbulent phase screen used in the simulation of empirically determining the correlation between temperature brightness and phase fluctuations. The scale is in arbitrary units, as the screen is later re-scaled so that the fluctuation on a 300 m baseline is 200 µm.

0.1 0.05 Normalised error 0.02 0.01 0.005 0.002 0.001 1 10 100 1000 Calibration interval (s)
Figure 4. Single baseline measurement: Normalised error on the estimated slope of the correlation between phase and sky-brightness fluctuations as a function of the time interval between calibration observations. It is assumed the interferometric path measurement has a normally-distributed error of 10 µm and the radiometric sky brightness measurement has a normally distributed error of 10 µm (red line), 20 µm (green) and 50 µm (blue).

0.1 Normalised error

The total observation time in all cases is 1000 seconds. As described above, the errors due thermal noise in the WVRs will depend on the conditions of observation and channel of the WVRs that is analysed, so we only show three representative cases by the three coloured lines on the plot. As can be seen in the figure, for the single baseline configuration we can expect errors in the range of 1% to 5% if the interval between calibration observations is in the range 20­100 s. One possibility therefore is to carry out experiments to simultaneously track the phase of a few quasars at a range of elevations. This is more efficient than doing these observations sequentially since allowing for time to elapse between taking data points results in greater phase changes and therefore higher accuracy. The results shown in Figure 5 are for a dedicated calibration scan, during which the phase calibration source is continuously observed. The vertical axis again represents the normalised error on the measured slope, while the horizontal axis represents the length of the calibration scan. The three lines are again calculated for representative values of errors due to thermal noise in the WVRs. It can be seen that for a range of parameters, the errors in inferred phase correction in this procedure are larger than in the procedure using interleaved observations of calibrators. The reason is again that interleaved observations allow greater changes in path to occur. Never the less, for scans longer than about 60 seconds, errors of the order of 5% or less can be expected.

0.01

3.4

Full array

0.001 10 20 50 100 200 500 1000 Length of scan (s)
Figure 5. Single baseline measurement: Normalised error on the estimated slope of the correlation between phase and sky-brightness fluctuations as a function of the length of the specialised scan for estimating it. It is assumed the interferometric path measurement has a normally-distributed error of 10 µm and the radiometric sky brightness measurement has a normally distributed error of 10 µm (red line), 20 µm (green) and 50 µm (blue).

correction coefficients can be estimated observationally with just a single baseline between two 12 m diameter antennas. The results are shown in Figures 4 and 5. The first figure (5) is for the observing strategy which consists of interleaved observations of calibrator and science targets. The vertical axis represents the normalised error in the estimated phase correction coefficient and the horizontal axis the interval between the phase calibration scans.

We next consider a more typical situation for ALMA, in which the whole bi-lateral array of 50 12 m-diameter antennas is used. In this situation, an estimate of the phase correction coefficient can be derived from each baseline, but with the caveat that determinations from individual baselines will depend on baseline length and will not be statistically independent. This can be tackled easily within the simulation already developed for the single baseline as follows. First, it should be noted that when calculating the antennabased complex gain solution from quasar calibration observations there is a degree a freedom that is always unconstrained by the observation, i.e., the "common-mode". This missing information is usually regularised by either assuming that the phase errors on one of the antennas are always zero, or that the mean phase error of all antennas is zero. The measurements from the WVRs are however intrinsically absolute measurements in which this common-mode is still present. Therefore correlating complex gain solution of an antenna with the WVR signal could give quite misleading results, if this common mode signal is significant. This is likely to be the case even in the extended configurations of ALMA, because of the tendency of atmospheric turbulence phenomena to increase in magnitude on longer length-scales.
2009 ALMA Memo 588


5
10
-1

10

-1

Normalised error

10

-2

Normalised error 1 10 100 1000

10

-2

10

-3

10

-3

10

-4

10 Calibration interval (s)

-4

1

10

100

1000

Scan duration (s)
Figure 7. As Figure 6, but for single continuous calibration observations of varying length. This figure is also for the full 50 antenna ALMA array.

Figure 6. Full 50 antenna ALMA array measurement: Normalised error on the estimated slope of the correlation between phase and sky-brightness fluctuations as a function of the interval between consecutive calibration scans. The total scan is assumed to be 1000 s long.It is assumed the interferometric path measurement has a normally-distributed error of 10 µm and the radiometric sky brightness measurement has a normally distributed error of 10 µm (red line), 20 µm (green) and 50 µm (blue). The whole array of 50 ALMA antennas is simulated.

better than 1% should be easily achievable unless significant errors in addition to those introduced by thermal noise become apparent.

One can get around this problem by considering the difference in signals between pairs of radiometers (thereby removing the common-mode) and correlating it with the difference in phase between the gains of the two antennas. This procedure is most naturally done by skipping the step of calculating the antenna phase gain solutions completely, since individual visibilities are in the present case nothing more than direct measurements of difference in phase of gains of pairs of antennas. Therefore in the calculation below we work directly with visibilities recorded on each baseline and correlate with the differences in the radiometer signals. The simulation method was as follows. Antenna path fluctuation were simulated for each of the antennas for the duration of the observation. The estimated noise is then added to the simulated antenna paths, and to the simulated radiometer outputs. Each of the baseline is then formed and slope of correlation estimated individually, together with the formal error on the slope, a . Finally estimates on each baseline are combined by: a= ¯ 1 j 1/
2 a, j j

4 4.1

USING `EMPIRICAL' PHASE CORRECTION COEFFICIENTS Direct use of estimated phase correction coefficients

The simplest way of using the empirically determined phase correction coefficients is simply to apply them directly to correct the phases for science observations, possibly interpolating between adjoining determinations. If the empirical estimate of the coefficients is made close in time and elevation to the science target, these should provide good correction. In practice this direct approach may not be desirable because: · Changes in elevation (possibly even azimuth, if it is a large change) between calibration and science observations, or during science observation, will cause an inaccuracy in the phase correction coefficients · More time than necessary may be spent on calibration · It may be difficult to predict the dispersive part of the phase since model parameters will be poorly constrained Beside the elevation change, it is difficult at this time to predict the limiting factors of this direct approach, but it should be possible to test this soon as ALMA starts to be commissioned at the high site.



a j /

2 a, j

(5)

where the sum over j runs over all baselines. It is essential to weight the estimates by their formal errors, because the errors on short baselines are much larger than on the long baselines. Note here we made an inherent assumption that the coefficients on all antennas are the same. There are two reasons we can foresee now that this will may not be the case: (i) Different atmospheric conditions at different antennas, especially in the most extended configurations when antennas will be separated by up to 15 km and have significantly different altitudes (ii) Different filter responses of WVR units If these differences are found to be significant in practice they can be included in the modelling although at a cost of a significant increase in complexity. The results for the full array simulations are shown in Figures 6 and 7, for the interleaved and continuous quasar observations respectively. With the full array, the accuracy of the measurements is of course much higher than on a single baselines and accuracies of
2009 ALMA Memo 588

4.2

Constraining the model parameters

A more efficient way of using the empirically determined dL/dTB,i values is to use them as an additional observational input to the Bayesian analysis procedure described by Nikolic (2009b). In this approach, the observed quantities would be:
(i) The four observed absolute sky brightness temperatures, TB,i and their associated error distribution (ii) The four empirically determined correlations between path and temperature fluctuations, dL/dTB,i , and the associated error which due to the simplified fitting procedure described in Section 3.2 will be parametrised by standard deviation only, a,i

Since an atmospheric model can directly predict both the absolute sky brightness and the correction coefficients, we can write down


6
the joint likelihood function for all these observables: 2 TB,i - TB,i dL/dTB,i - dL/dTB log L = - - T ,i a,i i there is a degeneracy between water vapour and the temperature, which prevents strong further constraints on the latter Overall, the result illustrates that even poorly estimated empirical phase correction coefficients can significantly improve the constraints on the model parameters that are not already well determined by the absolute measurements of the WVRs, in particular by removing some of the degeneracies present in the models. The significance of this increases in more complex physical models with more parameters, in which further degeneracies are bound to appear. The final figure in this section (Figure 10) shows the marginalised distributions of the phase correction coefficients corresponding to the model parameter distributions shown in Figures 8 and 9. Since observations of the empirical parameters are used in the inference of the model parameters, it can be expected that when the errors on the empirical coefficients are small enough, they completely dominate the marginalised distribution too. This can indeed be seen for errors of 0.5% and 2% ­ the distributions are close to the Gaussian distribution of empirical coefficients which were put into the problem as observations. On the other hand when errors on empirical coefficients are of the order of 10%, the constraints from the absolute brightness and the empirical coefficients combine to produce an error on phase correction coefficients which is smaller than the input 10%.

,i

2



(6) where the index i runs from 1 to 4, corresponding to the four channels of the WVRs and I've again assumed that the errors on the absolute sky brightness are normally distributed. This likelihood function, in conjunction with any relevant priors, can then be analysed using the Markov Chain Monte Carlo (MCMC) procedure as described in the previous memo. The first question to consider is: how much do the empirical phase correction coefficients constrain the model parameters in the simple model of a single, thin, water vapour layer? As in the previous memo, this can be analysed using the calculated distributions of model parameters using a simulated data point as the input. In this case, the simulated model parameters are c = 1.0 mm, T = 270 K, P = 580 mBar which are then used to calculate both four sky brightness temperatures and the four phase correction coefficients. As before, I then assign a 1 K normally distributed random error on the absolute sky brightness and a range of errors, from 100% to 0.5%, on the phase correction coefficients. The MCMC procedure is the used to estimate the confidence intervals on the model parameters {c, T , P}. The results of this experiment are shown in Figure 8, which shows the inferred distribution of model parameter for errors on empirical coefficients which are 100%, 10%, 2% and 0.5% of their values. By comparing the two top rows of this figure, it can be seen that even relatively noisy empirical estimates of the dL/dTB significantly constrain the model parameters. For example, in the second row of Figure 8, which corresponds to 10% error on empirical estimates, the temperature of the water vapour layer is already constrained to a range of about 35 K and the pressure is constrained to a range of about 100 mBar. This is contrast to the top row, which corresponds to essentially no useful information from the empirical coefficients. Here the possible ranges of temperature and pressure are 100 K and 300 mBar respectively. It is also useful to look in detail at the effect adding information from empirical coefficient has on the distribution of the model parameter representing the water vapour column. Even with the weak constraints provided by empirical coefficients with 10% relative error, it can be seen that the non-Gaussian tail from the distribution to be eliminated as a possibility. However, tighter constrains from empirical coefficients with smaller errors do not significantly reduce the uncertainty of the retrieved water vapour column while they moderately improve the errors on the pressure and temperature parameters. The reason for this can be understood from Figure 9, which shows the change in the joint distribution of parameters representing water vapour column and temperature, as the uncertainty on empirical coefficients is reduced: · The improvement in the inference of the water vapour column between the 100% and 10% error on coefficients (i.e., the first and second rows) is significant and is in turn due to the better constraints on temperature of water vapour which is provided by the empirical coefficient measurement · There is little further improvement in the water vapour column. This is because this parameter, when there are no degeneracies, is already strongly constrained by the absolute brightnesses measured by the WVRs · Over a smaller temperature range (in this case 266 K­275 K),

4.3

Transfer to different elevation

One of the key reasons for combining the information from the empirical coefficients with a physical model of the atmosphere is that it provides a way of extrapolating the results to slightly different conditions. Probably the most important example of this is a change in the elevation of the antennas, which will usually be necessary when moving from the calibration source to the science source and also while tracking the science source as it rises or sets. By making the usual approximation of a plane parallel atmosphere we have a direct prescription on how to modify the parameters of the model. In the case of the present, very simple, model this is simply to scale the water vapour column as sec( ) while leaving the other parameters unchanged. This transformation can be carried out on the full posterior distribution of the parameters as obtained from the inference procedure, leading to a new estimate for the distributions of the correction coefficients at the new elevation. An example of the procedure of transferring the results of the inference to a different elevation are shown in Figure 11. In this example a simulated observation at zenith with an assumed 2% error on empirical coefficients is analysed to provide the posterior distribution of the model parameters. These are then transferred to an elevation of 70 degrees (zenith angle of 20 degrees) and the new distribution of phase correction coefficients are computed. As can be seen the distributions of the coefficients are shifted between the two elevations corresponding to different levels of saturation of the 183 GHz water vapour line. The overall shapes are still largely the same however, and still approximately normally distributed, and there are no undesirable tails. This of course is simply the consequence of the good constraints on temperature and pressure we obtained at the zenith which we know (in this model, and probably in practice) will not change with elevation of the antennas. More advanced extrapolations are also possible. In particular, the sky toward the science source will have an intrinsically different amount of water vapour compared to the calibration source direction, over and above what is predicted by the plane parallel atmospheric
2009 ALMA Memo 588


7
0.08 0.06 0.04 0.02
0.02 0.015 0.01 0.005 0.02 0.015 0.01 0.005

f

f

0 0.95

0

f

0 200 220 240 T (K) 260 280 300 300 400 500 P (mbar) 600 700

1

1.05 c (mm)

1.1

1.15

1.2

(a) 100% error on dL/dTB
0.04 0.03 0.02 0.01
0.025 0.02 0.015 0.025 0.02 0.015 f 0.01 0.005 0 0.01 0.005 0 240 250 260 270 T (K) 280 290 300 500 525 550 575 P (mbar) 600 625 650

f

0 0.95

0.975

1 c (mm)

1.025

1.05

f

(b) 10% error on dL/dTB
0.04 0.03 0.02 0.01

0.04 0.03 0.02 0.01

0.03

0.02 f 0.01 0 260 265 270 T (K) 275 280 550 575 600 P (mbar) 625 650

f

f
0.975 1 c (mm) 1.025 1.05

0 0.95

0

(c) 2% error on dL/dTB
0.04 0.03 0.02 0.01

0.04 0.03 0.02 0.01

0.04 0.03 0.02 0.01

f

f

0 0.95

0
0.975 1 c (mm) 1.025 1.05

f 265 267.5 270 T (K) 272.5 275

0 580 590 600 P (mbar) 610 620

(d) 0.5% error on dL/dTB Figure 8. Marginalised distributions of model parameters when the retrieval is made simultaneously from the absolute brightness temperatures (assuming 1 K Gaussian error on each of the four channels) and from the correlation between path and brightness fluctuations with Gaussian errors as shown below each of the panels. A weak prior on the possible temperature and pressure of the water vapour is applied to these retrievals, i.e., the constraint that the temperature is between 200 K and 320 K and the pressure between 100 and 650 miliBar (this is the row 1 of Table 1 in Nikolic 2009b).

model and the change in elevation intrinsic change can be constrained ness measurements by the WVRs source as an additional input to the

between the two sources. That by using the absolute sky brightin the direction of the science retrieval.

A possible implementation of this approach is to have two variables representing water vapour column, one in the direction of the calibration source at time of empirical coefficients were determined and one in the direction of the science target at the current time.
2009 ALMA Memo 588

The two variables can then be tied together with a prior obtained through experience of how quickly the conditions at the site change. The observables in this case are the absolute sky brightness and empirical coefficients in the direction of the calibrator and the current absolute sky brightness in the direction of the science target. Such an approach would be able to maximise the time during which the information obtained from empirical coefficients is useful.


8
100% error on dL/dTB
290 280 3000 280 260 T (K) T (K) 2000 270 1000 1500

10% error on dL/dTB

240 1000 220 250 0 1 1.05 c (mm) 1.1 1.15 0 1 0.97 0.98 0.99 1 c (mm) 1.01 1.02 1.03 0 0 1 260 500

2% error on dL/dTB
277.5 2500 275 272.5 T (K) 270 267.5 265 262.5 0.97 0.98 0.99 1 c (mm) 1.01 1.02 1.03 2000

0.5% error on dL/dTB
274 4000 272 3000 T (K) 270 2000
1000

1500

268 1000
500

266
0 0 1

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

Figure 9. Joint distribution of the model parameters representing the water vapour column (horizontal axis) and the temperature of the water vapour layer (vertical axis) based on inference with the information from empirical phase correction coefficients included. The four panels are for the four different assumed errors on the observed coefficients. Note that the scales changes between the panels.

5 ILLUSTRATION WITH DATA FROM THE SMA In this Section I make use of an one-hour long stretch of test data obtained with the prototype ALMA water-vapour radiometers at the Submillimetre Array (SMA) to illustrate the above concepts of obtaining the empirical phase correction coefficients and using them in the retrieval of atmospheric parameters. This data set is the same data set as used in the previous memo of this series (Nikolic 2009b). The original data obtained at the SMA are shown in the topleft panel of Figure 12. The top plot in this panel is the difference between the sky brightness temperatures recorded at each integration time by the WVRs on the two antennas used in this test. Only the data from the outermost channel of the WVR is shown (as discussed in the previous memo, the data from the innermost channels are not very useful for correction in this case as the water vapour line is almost optically thick). The bottom plot is the interferometric phase recorded between the antennas and converted to a path-length fluctuation. Like most of the data that was collected at the SMA, the integration time in this experiment was 2.5 s. Also shown in Figure 12 are the trivial simulations of the type of phase measurement that may be obtained by fast-switching observations with ALMA. These were made by simply decimating the data collected at the SMA so that the remaining data resemble fast-switching observations with cycles of 25, 60 and 120 seconds. It is clear from these plots that even with a 120 second cycle time there are correlations between the measured phase and sky temperature

Interval (s) 2.5 25 60 125

Ch 1 (µm/K) -161 -181 -199 -275

Ch 2 (µm/K) 94 88 7 34

Ch 3 (µm/K) 388 397 423 394 (0.02) (0.22) (0.8) (1.2)

Ch 4 (µm/K) 351 345 361 335 (0.01) (0.14) (0.44) (0.7)

Table 1. Summary of inferred linear coefficient between phase and sky brightness fluctuations for the full (top row) and the decimated data from the SMA. The values for channels 1 and 2 were calculated using the standard least-squares straight line fitting since the correlation is so poor that the algorithm of Section 3.2 produces results with extremely large formal errors. The number in the parenthesis for channels 3 and 4 is the estimated formal error from the line fitting procedure.

difference which contain useful information about the empirical phase correction coefficients. By again making the linear assumption (Section 2), these data can be used to estimate the empirical phase correction coefficients. For these data (as discussed in the previous memo) channels 1 and 2 have a high optical depth and are not useful for phase correction. For the same reason, application of the line fitting with errors in both coordinates (Section 3.2) produces formal errors much larger than the actual coefficients. Therefore the line fitting with errors in both coordinates is only applied to channels 3 and 4 while for channels 1
2009 ALMA Memo 588


9
0.03 0.03 0.03

0.02 f f

0.02 f 0.01

0.02

0.01

0.01

0 0.05

0.075
dL d TB,2

0.1 (mm/K)

0.125

0.15

0 0.08

0.1
dL d TB,3

0.12 (mm/K)

0.14

0.16

0 0.14

0.16

0.18
dL d TB,4

0.2 (mm/K)

0.22

0.24

(a) 100% error on dL/dTB
0.03

0.03

0.03

0.02 f

0.02 f
f

0.02

0.01

0.01

0.01

0 0.06

0.065

0.07
dL d TB,2

0.075 (mm/K)

0.08

0.085

0 0.08

0.085
dL d TB,3

0.09 (mm/K)

0.095

0.1

0 0.14

0.145

0.15
dL d TB,4

0.155 (mm/K)

0.16

0.165

(b) 10% error on dL/dTB
0.04 0.03 0.02 0.01 0 0.064 0.04 0.03 0.02 0.01 0 0.084 0.04 0.03 0.02 0.01 0 0.144

f

f

0.066

0.068
dL d TB,2

0.07

0.072

0.086

0.088
dL d TB,3

0.09

0.092

f

0.146

0.148
dL d TB,4

0.15 (mm/K)

0.152

0.154

(mm/K)

(mm/K)

(c) 2% error on dL/dTB
0.04 0.03 0.02 0.01 0 0.0675 0.06 0.04 0.03 0.02 0.01 0 0.145

0.04 f f 0.02 0 0.087 0.0875 0.088
dL d TB,3

f

0.068

0.0685
dL d TB,2

0.069

0.0695

0.0885 (mm/K)

0.089

0.0895

0.1475
dL d TB,4

0.15 (mm/K)

0.1525

(mm/K)

(d) 0.5% error on dL/dTB Figure 10. Marginalised distributions of the phase correction coefficients corresponding to the parameter distributions in Figure 8 (to make the panels more readable, only three of four coefficients are shown).

and 2 the simple traditional fitting with error in one coordinate only is used. The data from channels 1 and 2 are not however used in any other way. The results of this analysis is shown in Table 1. Concentrating on channels 3 and 4 which are relevant in this case, it should first be noted that decimation of the data from the full data to one data point every 125 seconds does not drastically change the value of the empirical phase correction coefficients, although some variation is seen and the formal errors on the values increase. This is consistent with the results of the simulations (Sec2009 ALMA Memo 588

tion 3.3) which showed that increasing the fast-switching cycle time decreases the accuracy of empirical phase coefficients in only a gradual manner. Secondly, it can be seen that the empirical coefficients for Channels 3 and 4 are somewhat different to the `best-fitting' coefficients obtained in Section 4 and Table 2 of Nikolic (2009b). The reason is that best-fitting coefficients of Nikolic (2009b) are found by minimising the sums of squares of the residual phase errors while in this paper I try to measure the true correlation between the measured


10
Zenith
0.04 0.03 0.02 0.01 0 0.064 0.04 0.03 0.02 0.01 0 0.064 0.04 0.03 0.02 0.01 0 0.084 0.04 0.03 0.02 0.01 0 0.144

20 degree zenith-angle

0.04 0.03 0.02 0.01 0 0.0675

f

f

0.066

0.068
dL d TB,1

0.07 (mm/K)

0.072

0.074

0.07

0.0725
dL d TB,1

0.075 (mm/K)

0.0775

0.08

0.04 0.03 0.02 0.01 0 0.0675

f

0.066

0.068
dL d TB,2

0.07

0.072

f

0.07

0.0725
dL d TB,2

0.075

0.0775

(mm/K) 0.04 0.03 0.02 0.01 0 0.086

(mm/K)

f

0.086

0.088
dL d TB,3

0.09

0.092

f

0.088
dL d TB,3

0.09 (mm/K)

0.092

0.094

(mm/K)
0.04 0.03 0.02 0.01 0 0.145

f

0.146

0.148
dL d TB,4

0.15 (mm/K)

0.152

0.154

f

0.1475

0.15
dL d TB,4

0.1525 (mm/K)

0.155

0.1575

Figure 11. Transfer of inferred phase coefficients from zenith to elevation 20 degrees away.

2009 ALMA Memo 588


11
4 4

3

3 TB (K) L (µ m) 17.2 17.4 17.6 17.8

TB (K)

2

2

1

0 -1 1000

1

0 1000

500

500

L (µ m)

0

0

-500

-500

-1000 17

-1000 18 17

17.2

17.4

17.6

17.8

18

t (hours UT)

t (hours UT)

(a) Every integration, i.e., every 2.5 s
3 2.5 2 TB (K) 1.5 1 0.5 0 400 1 0.5 500 TB (K) 2 1.5 3 2.5

(b) Every 10th integration, i.e., every 25 s

200

250

L (µ m)

L (µ m)

0 -200 -400 -600 17

0

-250

17.2

17.4

17.6

17.8

-500 18 17

17.2

17.4

17.6

17.8

18

t (hours UT)

t (hours UT)

(c) Every 24th integration, i.e., every 1 minute

(d) Every 50th integration, i.e., every 2 minutes

Figure 12. Illustration of decimation of single-baseline data from the SMA to simulate inference of correlation between sky brightness and path fluctuation using fast-switching calibration observations. In each panel, the upper plot shows the difference in radiometric measurements (TB , using channel four of the WVRs) and the bottom plot shows inferred path between the antennas (L).

2009 ALMA Memo 588


12
phase and the sky brightness fluctuation. Because there are sources of error in both of these quantities, simply minimising the squares of the residuals will lead to an underestimate of the correction coefficients. There is a good discussion of this effect by Stirling et al. (2005). The final part of sections shows the results of the joint analysis of the absolute sky brightness and the empirical phase correction coefficients (as discussed in 4.2) to place constraints on the model parameters. For this analysis I took the absolute sky brightness observed in the middle of the SMA observation shown in Figure 12. The elevation of this time was 22 degrees and the source was setting rapidly. For this reason, I computed the empirical coefficients for a shorter time around the middle of the scan only. As described above, the results for channels 1 and 2 were not usable in these conditions, while the values of the coefficients for channels 3 and 4 were 313 and 324 µm/K respectively. Because in this memo I again do not take into consideration the dispersive effects of the water vapour, these figures need to be adjusted for the expected dispersive phase at the frequency of observation (240 GHz in this case). As in Nikolic (2009b) I take this correction factor to be 1.05, resulting in coefficients of 298 and 309 µm/K for Channels 3 and 4 respectively. These empirical coefficients were then analysed together with the four absolute brightness which in this case were 268.3, 247.4, 200.4, 135.9 K for channels 1 to 4 respectively. The errors assigned on the absolute sky brightness were again 1 K while I assigned error of 4% on the empirical coefficients. The marginalised distributions of the model parameters produced by this analysis are shown in Figure 13 side-by-side with the results of an analysis based on the absolute sky brightness only (essentially the same as shown in the previous memo). From this comparison, it can be seen that significantly different distributions are inferred for the all three of the model parameters. In particular, even the inferred water vapour column is moved to a significantly higher value, from about 1.2 to 1.3 mm of water vapour at zenith. Also the pressure distribution is moved from being mostly constrained by the high-end of the prior range to being mostly constrained by the low-end of the prior range. These changes in model parameters when the observations of the empirical coefficients are introduced into the analysis should of course be interpreted as the necessary moves to reconcile the predicted coefficients with those actually observed. The fact that the two model parameter distributions do not overlap does, however, suggest that the inclusion of the observed empirical coefficients is also correcting for a modelling and/or calibration error and not just for degeneracies present in the problem as can be seen in Figure 8. Finally, the marginalised distributions of the predicted phase correction coefficients are shown in Figure 14, again in a side-byside format with the retrieval based just on absolute sky brightnesses on the left and the retrieval based on the empirical coefficient observations and the absolute sky brightness on the right. As expected, the distributions of the phase correction coefficients are narrower (and closer to Gaussian) when the observed coefficients are included the analysis, both for the channels with observational constraints (channels 3 and 4) and for the other channels (1 and 2). The second point that should be noted is that the marginalised distribution of channels 3 and 4 are not centred on the observed values that were entered into the analysis. The reason is that the assigned 4% error allows for variation in the retrieves values while the observed absolute sky brightness temperatures clearly prefer a lower value of the coefficients. 6 SUMMARY AND DISCUSSION

The results presented above may be divided into two main parts. Firstly, in Section 3, the limits due to thermal noise on observational determination of the phase correction coefficients are calculated. Because only thermal noise in the receivers is considered, this analysis should be regarded a best-case scenario; in practice effects such as electronic and/or mechanical phase drifts will always be become important on a long enough timescale. These effects are by their nature very much less predictable than thermal noise in receivers and realistic estimates of their magnitude will be possible only once the commissioning of ALMA at the array operations site begins. In terms of the technique of estimating the empirical phase correction coefficient from observed data two points from this section should be noted: (i) Best-fitting linear relationship should be computed taking into account errors on both of the coordinates in the fit (i.e., both on the measured radiometric fluctuation and measured phase). (ii) When combining estimates from baselines in a twodimensional, array, they must be correctly weighted by their error estimates The last of these is important because errors on estimates on short baselines are much greater than the estimates from long baselines. The main implications of results in Section 3 is that it should be possible to accurately and relatively quickly estimate the empirical phase correction coefficients. Even with a single baseline array, it should be possible to get to about 1% accuracy in about 200 seconds. With full ALMA in a medium configurations accuracies of order of 0.1% could be reachable if non-thermal effects do not become important. These predicted accuracies should be compared to results of apriori modelling of the phase correction coefficients as for example given in the previous memo (Nikolic 2009b). There, I estimated the error on the predicted phase correction coefficients due to the uncertainties in the input parameters only, i.e., intrinsic error on the measurement of absolute sky brightness with the WVRs and the poorly known pressure/temperature of the water vapour in the atmosphere. The error estimates in Nikolic (2009b) do not therefore include modelling error. I found that if we do not have a good constraint on temperature and or height of the water vapour, the range of calculated coefficients spans about 8%, while if there are good constraints, the range spans around 3-4%. Therefore, according to the above calculations, relatively short empirical measurements are likely to have smaller errors than the a-priori predictions. For this reason empirical observations should be very useful for phase correction, at least until we learn how to estimate the atmospheric parameters which are not well constrained by absolute measurements from the WVRs (the parameter that is well constrained by the WVRs is of course the total water vapour column). The second part of the memo (Sections 4 and 5) concerns the way empirically determined phase correction coefficients are used in practice. The thesis put forward in this part of the memo is that the best way of using these coefficients is to regard them as additional observational measurement which are used to constrain a physical model of the atmosphere, such as the simplified model presented in the first memo of this series (Nikolic 2009b). The advantages of this approach are that: (i) It is possible to easily and accurately transfer measured empir2009 ALMA Memo 588


13
Without empirical coefficients
0.05 0.04 0.03 f 0.02 0.01
f 0.05 0.04 0.03 0.02 0.01

With empirical coefficients

0 1.15
0.04

0

1.175

1.2

1.225 c (mm)

1.25

1.275

1.3
0.04

1.2

1.225

1.25

1.275 c (mm)

1.3

1.325

1.35

0.03

0.03

0.02

0.02

f

0.01

f 272 274 276 T (K) 278 280

0.01

0

0 268 0.1 0.08 0.06 f 0.04 270 272 T (K) 274 276 278

0.03 0.025 0.02 0.015 0.01 0.005 0 525 550 575 P (mbar) 600 625 f

0.02

0 530 540 550 560 P (mbar) 570 580 590

Figure 13. Illustration of the retrieval of model parameters using as input only the absolute sky brightness (left column) and the absolute sky brightness and the empirical coefficients for channels 3 and 4 (right column) for SMA observation on 17 February 2006. The absolute sky brightness was taken from the middle of the observation scan, while empirical coefficients were estimated for a few minutes around the middle of the scan. In this case I assigned an error of 4% to the measurements of the empirical coefficients, and the error on the absolute sky measurement is as before assumed to be 1 K.

ical coefficients to different conditions, such as different elevation or a part of the sky with slightly different pwv column (ii) Constraints in the underlying physical model can improve the accuracy of the empirical measurements (iii) The same constraints can also flag implausible measurements of the empirical coefficients (iv) It is possible to do quantitative model selection, using all of the available information The existing algorithms presented in the previous memo can naturally be extended for such an analysis and this has been already implemented in the libAIR library.
2009 ALMA Memo 588

Finally here are some caveats and directions for future work. One of the topics that has not been addressed in this memo are the effects of dispersion. If they can be modelled accurately, there is no reason why dispersive effects can not be taken into account in the calculation of the predicted dL/dTB and the subsequent analysis can proceed as before. A second topic which may be important is that of the so-called `dry' fluctuations, that is, fluctuations in phase due to changes in the density of dry air rather than quantity of water vapour. If these are completely independent of the `wet' fluctuations then it will certainly not be possible to correct for them, and they will manifest


14
Without empirical coefficients
0.05 0.04 0.03 f 0.02 0.01
f 0.05 0.04 0.03 0.02 0.01

With empirical coefficients

0 0.75

1

1.25
dL d TB,1

1.5 (mm/K)

1.75

2

0 1.25

1.5

1.75
dL d TB,1

2 (mm/K)

2.25

2.5

2.75

0.04

0.05 0.04 0.03

0.03

0.02

f

f 0.02 0.01

0.01

0 0.3 0.325 0.35 0.375
dL d TB,2

0.4

0.425

0.45

0 0.375

0.4

0.425
dL d TB,2

0.45 (mm/K)

0.475

0.5

0.525

(mm/K) 0.05 0.04 0.03

0.04

0.03

0.02

f

f 0.02 0.01 0.22 0.23
dL d TB,3

0.01

0 0.21

0.24 (mm/K)

0.25

0.26

0 0.23

0.24

0.25
dL d TB,3

0.26 (mm/K)

0.27

0.28

0.04

0.06 0.05

0.03 0.04 0.02 0.03 0.02 0.01 0.01 0 0.24 0 0.25 f f 0.25
dL d TB,4

0.26 (mm/K)

0.27

0.28

0.26
dL d TB,4

0.27 (mm/K)

0.28

0.29

Figure 14. Inferred phase correction coefficients corresponding to Figure 13, i.e., based on only the absolute sky brightness (left column) and on the absolute sky brightness and the empirical coefficients for channels 3 and 4 (right column).

2009 ALMA Memo 588


15
themselves as an extra source of error when calculating the empirical correction coefficients. There errors assigned to the path-fluctuation coordinate used in the line-fitting procedure should take into account this extra source of error. It is likely however that the `dry' fluctuations are partially correlated with the `wet' fluctuations, as for example predicted by the Large Eddy Model simulations of Stirling et al. (2008). In this case the observed phase correction coefficients will be scaled by an unknown factor from their predicted values. If this effect is significant, it will likely be counter productive to let the inference procedure try to match the observed coefficients. Rather, the scaling should be introduced into the model as a further parameter and Bayesian evidence calculation then used to determine the significance of the dry fluctuations.
Internal revision control information: 207(2009-07-24 16:33:15 +0100)

ACKNOWLEDGEMENTS I would like to thank Robert Laing for a careful reading of the manuscript and suggestions for improving it and making it clearer. This work is a part of the ALMA Enhancement programme at the University of Cambridge funded by the European Union Framework Programme 6. More information on the programme is available at http://www.mrao.cam.ac.uk/projects/alma/ fp6/. I would like to thank the team involved in construction and the testing at the SMA of the ALMA prototype water vapour radiometers: P. G. Anathasubramanian, R. E. Hills, K. G.Isaak, M. Owen, J. S. Richer, H. Smith, 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.

REFERENCES Gull S. F., 1989, in Maximum Entropy and Bayesian Methods, Skilling J., ed. Nikolic B., 2009a, BNMin1 Minimisation/Inference library. http: //www.mrao.cam.ac.uk/~bn204/oof/bnmin1.html --, 2009b, Inference of Coefficients for Use in Phase Correction I. ALMA Memo Series 587, The ALMA Project 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 Press W., Teukolsky S., Vetterling W., Flannery B., 1992, Numerical Recipes in C, 2nd edn. Cambridge University Press, Cambridge, UK Stirling A., Holdaway M., Hills R. E., Richer J. S., 2005, Calculation of Integration Times for WVR . ALMA Memo Series 515, The ALMA Project 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) Wilson T. L., Rohlfs K., Hnttemeister S., 2009, Tools of Radio Astronomy. Springer
2009 ALMA Memo 588