Physics 728
Radio Astronomy:  Lecture #9
Prof. Dale E. Gary
NJIT

Calibration

Introduction

We have discussed how an interferometer array measures the complex visibilities of the source.  This is a good opportunity for a review of that relationship.  The basic need for calibration is to correct the measured visibilities V'(u,v) to approximate as closely as possible the true visibilities V(u,v).

The basic equation for a synthesis array is
 

V(u,v) = A(l,m) I(l,m) exp[-i2p(ul + vm)] dl dm .         (1)

where (l,m) are the direction cosines with respect to the phase center,
           (u,v) are the projected baseline coordinates in wavelengths, u = Bl,l, v = Bm,l,
           V(u,v) is the true visibility evaluated at u, v,
          A(l,m) is the normalized primary beam pattern (beam of a single antenna)
            I(l,m) is the brightness distribution of the source

Recall that the exponent comes from Bl. s = ul + vm + wn = ul + vm + tgn, so when we say "with respect to the phase center," we mean that the group delay is compensated so that tgn  = 0.

We can rewrite this to emphasize the discrete sampling with antenna pair i,j at times t by

Vij(t) = A(l,m) I(l,m) exp[-i2p(uij(t)l + vij(t)m)] dl dm .         (2)

The term ul + vm is the geometric phase difference Df, produced by the geometric path length difference between antenna i and antenna j from the part of the source at location (l,m) relative to the phase center.

Recall that

u
 1 
sin ho
cos ho
       0 
.. Bx
Bl. s v =

-sin do cos ho
  sin do sin ho
 
    cos do
By
  w  l
cos do cos ho
-cos do sin ho
 
    sin do
Bz

where the hour angle and declination (source coordinates) are ho, do, and the baseline lengths are Bx, By, Bz.  Then the geometric phase difference at the phase center (the w term, at l = m = 0) is

fg = 2p tgn  = (2p / l)[Bxcos do cos ho - Bycos do sin ho + Bzsin do].
We can see what can affect the geometric phase by taking the differential of this expression:
dfg = 2pn dtg  = (2p / l)[dBxcos do cos ho - dBycos do sin ho + dBzsin do
                                  + daocos do (Bxsin ho + Bycos ho)
                                    + ddo (-Bxcos ho sin do + Bysin ho sin do + Bzcos do)],           (3)
where we have used the relation between right ascension and hour angle: ho = LST - ao, so dho = dao.  Equation (3) shows how baseline errors (dBx, dBy, dBz) and source position errors (ao, do) will affect the error in group delay dtg (or yield an error in phase dfg).  Note that a clock error is equivalent to a source position error dao.

If we have a source whose position is known, we can use equation (3) to find the location of the antennas (this is called baseline determination).  The error in antenna position is largely independent of the baseline lengths.  For example, say that we can measure dfg to within 1o at 5 GHz (l = 6 cm).  Then we can measure dBx, dByand dBz to a precision of order

(1 / 360) 6 cm ~ 1 / 60 cm
even though B = (Bx2 + By2 + Bz2)1/2 = 5000 km or more (VLBI).

This has been used to verify and accurately measure motion of the Earth's crust due to plate tectonics.

Alternatively, once the baselines are accurately known, we can observe unknown sources and find their positions, (ao, do), to high accuracy.  This does require long baselines, because to minimize (dao, ddo) we have to have large (Bx, By, Bz).  So once again, VLBI is needed.  Determining accurate source positions is called astrometry, and VLBI is the most accurate method for establishing a highly accurate coordinate system for the sky, of order 1 milliarcsecond in absolute position.

Complex Gains
A properly designed synthesis array will preserve a linear relationship between the measured visibilities  V'(u,v) and the actual visibilities V(u,v), which in general can be written
Vij(t) = Gij(t)Vij(t) + eij(t) + hij(t)         (4)
where we have broken the constant term into an offset eij(t) and a stochastic noise term hij(t).  The proportionality factor Gij(t) is called the complex gain for the ij-th baseline.  All of the terms in this equation are "complex," but this is really just a convenient way of using one equation to represent the two outputs of each correlator -- the cosine and sine components (or real and imaginary).

The meaning of the complex gains can be seen by expressing them in terms of gains of individual antennas:

Gij(t) = gi(t) gj*(t) = ai(t) aj(t) exp[i(fi(t)- fj(t))]         (5)
where ai(t) is the multiplicative gain of the i-th antenna, and fi(t) is the phase of the i-th antenna. When the array is optimized, these numbers can be determined by measuring a point source at the phase center.  Since the expected amplitude of such a source is constant, and the expected phase is zero, any differences in the ai and any non-zero phase differences fi - fj must be due to the antennas themselves.  Note that for N antennas, we measure Gij on N (N - 1)/2 baselines, so for large N the problem is well over-determined.  In this case we can use a least squares algorithm to solve for the antenna-based complex gains gi, i = 1, 2, ... N.

However, first the array must be optimized, which requires some initial calibrations:

Antenna Pointing and Gain
Our equation for inverting the visibilities assumes that the normalized antenna primary beam distribution A(l,m) is independent of time and identical for each antenna.  The antenna tracking center must follow the same intended position for all antennas (generally the same position as the phase tracking center, although this is not essential).  Errors in tracking can cause reduced sensitivity and distortions in extended objects.  For point sources, the tracking accuracy should be better than 0.1qFWHM, where qFWHM = angular size of FWHM of the primary beam.  For extended objects (e.g. the Sun) that cover a large part of the primary beam, the pointing accuracy should be even better, say qFWHM/20.  For OVSA at, say, 10 GHz, where the 27 m primary beam is 4.6 arcmin, this would require about 14" accuracy.  For the 2 m antennas, the requirement is relaxed to about 3 arcmin.

The antenna pointing error is the difference between the actual pointing position and the desired one, and it has a complicated directional dependence due to

Often, rather than trying to describe axis and deformation errors based on engineering of the structure, an ad-hoc trigonometric dependence on the sky is determined using multiple measurements of sources around the sky.  For the OVSA 27 m antennas, we point one antenna at the source and move the other in a prescribed pattern. (J069 arcfile)

For the OVSA 2 m antennas, which are not sensitive enough for this, we use the solar disk (which is nearly a point source for the 2 m primary beam).  One problem is that this gives only a single declination track over a day, so it is not possible to optimize for pointing over the entire sky. (J067 arcfile)

Delay Calibration
We have already discussed delay calibration in some detail in Lecture 8, as well as the related ideas of bandpass calibration.  Recall that an error in delay of only a few nanoseconds will cause the different parts of the observing band to partially interfere (different frequencies in the band arrive with different phase), leading to a drop in the amplitude.  Obviously, we need to have the delays optimized.

Time and Place
The time of day and location of the antennas must be known to relatively high accuracy -- needed for determining the geometric delay.  A clock error of 1 s, or a baseline error of a few cm, will cause a serious phase shift of the source over, say, 10 minutes.  At OVRO, using a GPS clock and measuring baselines with cosmic source calibration (we will discuss this shortly), we get a time accuracy of << 1 ms, and baseline errors of about 3 mm.  Therefore, these effects are not serious over a short time interval, but may still be problematic over 8 hours.  This is one reason that we do phase calibration observations every ~ 2 hours.

Routine Corrections
Several additional corrections must be done, generally by the "on-line" computer in real time.  These are:

  1. System optimization
    1. level control, or gain control -- for VLA, this is done using an ALC (automatic leveling control) loop to keep the threshhold optimum on the digital correlator.  At OVRO, descrete attenuation is switched in as necessary to avoid saturation (and non-linearity).
    2. Round trip phase measurement (to stabilize for phase changes due to temperature effects in the LO reference frequency).
    3. Pointing control (adjusting the antenna pointing according to the pointing parameters).
  2. Timing corrections and calibrator fluxes
    1. The Earth does not rotate uniformly, due to wind, tide, and other global effects, precession and nutation.  These parameters are available in terms of timing corrections from government sources.
    2. The fluxes of calibrators can change and must be monitored by specially designed radio telescopes for absolute flux calibration.  We will see more about this shortly.
  3. Path length changes due to the troposphere and ionosphere.
    1. The refractive index of the atmosphere differs from unity by about 1 part in 3000, so there is an additional delay in traversing the atmosphere of about 8 ns, or an "extra path length" of about 23 m at the zenith.  This would not matter for a plane parallel atmosphere in which the antennas are alll at ground level, since only the relative phase matters.  However, if we refer the phase to a calibrator at some other sky location (different zenith angle) we have to take account of the extra path length.
      • DL = (Dh + ctgL sec z / ro) sec z
      where z is the zenith angle, L = 0.228 Ptot + 6.3 w is the zenith path length due to the atmosphere, Ptot is the total pressure at ground level, w is the vertical column water vapor content (cm), ro is the Earth radius, and Dh is the height difference of the antennas.  The second term is small, ~ 10-5 of baseline length, except for large z.  The wet component is generally a small fraction of the dry component, but it is also variable over the sky, so Df is in fact dominated by the wet component.
    2. The ionosphere causes a frequency-dependent refraction that is proportional to l2, so is worse at low frequencies.  The "extra" path length is
      1. Li = -40 n -2 Ne
      and is negative because the index of refraction n < 1 in an ionized medium.  Ne is the electron column density in units of 1018m-2 (typically near unity for Earth's ionosphere).  The differential path length is then
        DLi = ctgLi / (ro cos2z + 2H)
      where H is the height of the ionosphere (about 600 km).
  4. Absorption by troposphere and ionosphere
    1. In addition to refraction, the atmosphere also absorbs radiation, which has two effects: reduces flux by S = So exp(-ta sec z), increases noise (atmospheric emission) by Tatm [1 - exp(-ta sec z)]

    2. where ta = a0 + a1Pv = atmospheric opacity and Pv = partial pressure of water vapor in millibars.  ta can be measured using "tipping radiometers" which measure the sky brightness as a function of zenith angle.  The water vapor line near 22 GHz is a favorite part of the spectrum to use for this purpose.

Figure 1: The dependence of phase errors due to the ionosphere and troposphere
as a function of frequency.  Phase errors as small as about 1 degree near 1 GHz
can rise to 10-20 degrees at both high (20 GHz) and low (100 MHz) frequencies.
The phase errors can be much greater than this under certain conditions.

The Offset and Noise Terms
Before getting to the final calibration, it is worthwhile to look at the unwanted offset and noise terms in equation (4).  The offset term, eij(t), is generally kept small by design, and the amplitude part can be measured and subtracted by attenuating the signal.  This is analogous to "dark current."  The phase part is generally not a problem, since offsets on two antennas should not be correlated and it can be reduced to zero by using phase switching in the receiver.

The noise term, hij(t), is part of the system and generally vastly outweighs the signal, as we discussed in Lecture 5.  But it can be reduced by averaging over some length of time.  As we saw when discussing sensitivity, the noise is proportional to sqrt(t), where t is the integration time.  Over this same period, the signal remains constant, so S/N increases.

An interesting question is, what kind of average should we perform?  Often the data are represented as amplitude and phase, and one can either average amplitude and phase separately, or average as a vector average.  A vector average may seem the best choice, but in many (most) cases the phase variations are uncorrelated with the amplitude variations (e.g. due to weather, temperature, or cable length changes).  In this case, it is better to average amplitudes separately, and fit for any phase wind before determining phase).

Final Astrometric Calibration
We are fortunate that there are many point-like radio sources scattered around the sky and relatively well isolated from other confusing sources.  When we measure such sources with an interferometer array, we can use the known visibilities (constant amplitude and zero phase) to calibrate the system.  This is particularly useful because the signals traverse the entire system, including the ionosphere, troposphere, etc.  It is not perfect, however, because generally the calibrator sources are in different sky locations than the source of interest -- unless the source of interest is strong enough for self calibration.  We will discuss this in a moment.

As an example, let us take a look at the VLA calibrator manual.  As you can see, there are lots of point sources, but they may not be point sources at all frequencies and size scales.  Therefore, there are restrictions of baseline length (u,v spacing) and frequency.  Note again that these sources can be variable.  There is a nice database maintained by the VLA with a Java interface, to allow checking for variability.

Baseline-based Vs. Antenna-based Calibration
After suitable averaging, we are left with (for a point source)

Gij(t) =  Vij(t)/S         S = flux of the point source
where Gij is the baseline-based gain (complex -- includes both amplitude and phase components).  For a sufficient number of antennas, however, these are not all independent and we can do better.  Let us write the baseline gain factor as
Gij(t) = gi(t)gj*(t) gij(t)
where the last term is a residual term very close to unity that basically keeps track of errors (usually < 1%).  These complex numbers can be separated into amplitude and phase:
Aij(t) = ai(t)aj(t) aij(t)
Fij(t) = fi(t) - fj(t) + fij(t).
Let us write Vij(t) = aij(t) exp[ifij(t)] ;  Vij(t) = Aij(t) exp[ifij(t)] . Then for a point source the measured amplitudes and phases are
Aij(t) = ai(t)aj(t) aij(t) S
Fij(t) = fi(t) - fj(t) + fij(t)
which can be solved for ai and fi for all N antennas, provided gij is close to unity, using a least squares method (using logarithms for the amplitude terms).  The error terms can then be checked to show if there are any problems.  The errors are acceptable when
fij(t) = Fij(t) - fi(t) - fj(t)  < 1o
aij(t) = Aij(t) / [ai(t)aj(t) S]  within 1% of 1.00.
One advantage of antenna-based solutions is that they can be obtained even without all of the baselines.  This is especially important for partially resolved or confused calibrators, as we saw with the restrictions in the VLA calibrator list.  In this case one simply applies the solution on baselines that meet the criteria so that the source is a good point source.  Note that this approach uses the fact that we have observed a point source calibrator.  In fact, one can also do this on a more complicated source, as we will now describe in the context of Selfcal (self calibration).
Self Calibration
For our discussion of self calibration, we will use Tim Cornwell's NRAO Summer School lecture notes.