PRECISION ASTROMETRY WITH ADAPTIVE OPTICS

With AO, the image motion of the guide star is removed with a flat tip-tilt mirror. This stabilizes the image of the guide star with respect to the imager to high accuracy. Any residual tip-tilt error is removed in subsequent analysis by calculating only differential offsets between the target of astrometry (not necessarily the AO guide star) and the reference stars. However, the difference in the tilt component of turbulence along any two lines of sight (LOS) in the field of view (FOV) causes a correlated, stochastic change in their measured separation, known as differential atmospheric tilt jitter.

More specifically, in propagating through the atmosphere to reach the telescope aperture, light from the target star and light from a reference star at a finite angular offset traverse different columns of atmospheric turbulence that are sheared. Differential atmospheric tilt jitter arises from the decorrelation in the tilt component of the wavefront phase aberration arising from this shearing effect. This differential tilt leads to a random, achromatic, and anisotropic fluctuation in the relative displacement of the two objects. The three-term approximation to the parallel and perpendicular components of the variance arising from differential atmospheric tilt jitter, assuming Kolmogorov turbulence, is given by (Sasiela 1994)

$$\begin{eqnarray} {\sigma ^2_{\parallel,{\rm TJ}} \atopwithdelims []\sigma ^2_{\perp,{\rm TJ}}} &=& 2.67 \frac{\mu _2}{D^{1/3}}\left(\frac{\theta }{D}\right)^2{3 \atopwithdelims []1} - 3.68 \frac{\mu _4}{D^{1/3}}\left(\frac{\theta }{D}\right)^4{5 \atopwithdelims []1} \nonumber\\ &&+\, 2.35 \frac{\mu _{14/3}}{D^{1/3}}\left(\frac{\theta }{D}\right)^{14/3}{17/3 \atopwithdelims []1}. \end{eqnarray} \tag{ 1 }$$

In this equation, D is the telescope diameter and θ is the angular separation of the stars. The turbulence moments μ_m are defined as

$$\begin{equation} \mu _m = {\rm sec} ^{m+1}\xi \int _0^\infty dh C_n^2(h)h^m, \end{equation} \tag{ 2 }$$

where h is the altitude, ξ is the zenith angle, and C²_n(h) is the vertical strength of atmospheric turbulence. Typical C²_n(h) profiles yield σ_∥,TJ ≈ 20–30 mas for a 20'' binary when observed with a 5 m aperture. Note that the variance from differential tilt is a random error, and, thus, is also ∝τ_TJ/t, where τ_TJ is the tilt jitter timescale (on the order of the wind crossing time over the aperture; see Section 5) and t is the integration time.

The largest instrumental systematic that limits the accuracy of astrometry in any optical system is geometric distortion. These distortions can be stable—resulting from unavoidable errors in the shape or placement of optics—or dynamic—resulting from the flexure or replacement of optics.

If geometric distortions are stable, then a number of strategies can be employed to mitigate their effect. One method is to model the distortion to high accuracy; the most notable example is the calibration of the Hubble Space Telescope (HST; e.g., Anderson & King 2003a). This is particularly important for data sets obtained with multiple instruments or those that use the technique of dithering, since knowledge of the distortion is necessary to place stellar positions in a globally-correct reference frame. Alternatively, one could use a consistent optical prescription from epoch to epoch by using the same instrument and placing the field at the same location and orientation on the detector. Here we use both a distortion solution and a single, consistent dither position to achieve accurate astrometry.

Any changes in the geometric distortion must be tracked through routine, consistent calibration. The question of stability is particularly important at the Hale 200 inch telescope, since the AO system and the imaging camera (Palomar Adaptive Optics System (PALAO) and Palomar High Angular Resolution Observer (PHARO), respectively; see Section 3) are mounted at the Cassegrain focus, and PHARO undergoes a few warming/cooling cycles per month (see Section 3). The PHARO distortion solution¹ by Metchev (2006) accounts for changes in the orientation of the telescope (which are relatively small for our experimental design), but the overall stability of the system is best verified with on-sky data. One of the purposes of the data presented here is to track the system stability. We find that the combination of the Hale Telescope, PALAO, and PHARO is capable of delivering ≲100 μas astrometry.

Refraction by the Earth's atmosphere causes an angular deflection of light from a star, resulting in an apparent change in its position. The magnitude of this deflection depends on the wavelength and the atmospheric column depth encountered by an incoming ray. The former effect is chromatic, while the latter is achromatic. The error induced by differential chromatic refraction (DCR) has proven to be an important, and sometimes the dominant, astrometric limitation in ground-based efforts (e.g., Monet et al. 1992; Pravdo & Shaklan 1996; Anderson et al. 2006; Lazorenko 2006). These studies have shown that DCR can contribute ≈0.1–1 mas of error depending on the wavelength and strategy of the observations.

The observations presented here were conducted using a Br-γ filter at 2.166 μm with a narrow bandpass of 0.02 μm to suppress DCR. The increased S/N provided by AO allows sufficient reference stars to be detected even through such a narrow filter in a short exposure time. We reach K_s ≈ 15 mag in our 1.4 s exposures through this filter with the Hale 200 inch telescope (see Section 3). In addition, observations were acquired over a relatively narrow range of airmass (1.17–1.27) at each epoch to minimize the achromatic differential refraction.

In order to estimate the effect of atmospheric refraction on our data, we took the asterism in the core of M5 and refracted it to 37 and 32° elevation with the parallactic angles appropriate for the observations on 2007 May 28 using the slarefro function distributed with the STARLINK library (Gubler & Tytler 1998). The root mean square (rms) deviation in reference star positions between these two zenith angles was ≲250 μas and the shift in the guide star position with respect to the grid (see Section 4) was ≈10 μas. Thus, our consistent zenith angle of observations, narrowband filter, and observations in the NIR (where the refraction is more benign) make the contribution of this effect negligible for our purposes, and we make no effort to correct for it.

Performing a similar experiment using a broadband K filter with a field of ≈5000 K reference stars and a ≈3000 K target would lead to a systematic shift of ≈100 μas between zenith angles separated by 10°, which would be detectable by this experiment. Consequently, for observations where broadband filters are necessary, refraction effects must be considered and corrected.

In the case of a perfect optical system, a perfect detector, and no atmosphere, the astrometric precision is limited to one's ability to calculate stellar centers. The centering precision is determined by measurement noise, and we will use the two terms interchangeably. For a monopupil telescope, the uncertainty is

$$\begin{equation} \sigma _{\rm meas} = \frac{\lambda }{\pi D}\frac{1}{\rm SNR} = 284\,\mu {\rm as}\left(\frac{\lambda }{2.17\,\mu {\rm m}}\right) \left(\frac{5\,{\rm m}}{D}\right)\left(\frac{100}{\rm SNR}\right) \end{equation} \tag{ 3 }$$

(Lindegren 1978). AO allow us to achieve the diffraction limit even in the presence of the atmosphere and substantially boosts the SNR over the seeing-limited case, thereby decreasing measurement noise and improving astrometric precision.

In practice, the centering of a given stellar image is limited by spatial and temporal variations in the AO point-spread function (PSF). A great deal of time and effort has been spent determining the AO PSF and producing software packages to perform PSF fitting (e.g., Diolaiti et al. 2000; Britton 2006). However, any PSF-fitting software package is capable of calculating image positions at a ≲0.01 pixel level in a single image. For the observations considered here, this is ≲2 mas, a factor of 5–10 larger than the measurement noise in Equation (3), but this is much smaller than the tilt jitter mentioned in Section 2.1. As such, we have chosen to use simple and widely-available PSF centering software (DAOPHOT; Stetson 1987; see Section 3).

The fundamental quantity in differential astrometry is the measured angular offset between a pair of stars. We will denote the angular distance between two stars, i and j, as $\vec{d}_{ij}$. Since $\vec{d}_{ij}$ is measured from an image, we will denote its components in the Cartesian coordinate system of the detector, simply

$$\begin{equation} \vec{d}_{ij} = { x_j - x_i \atopwithdelims []y_j - y_i } \equiv { x_{ij} \atopwithdelims []y_{ij} }, \end{equation} \tag{ 4 }$$

where we have introduced the notation x_ij ≡ x_j − x_i and likewise for y. The variance in the angular separation between two stars is given by

$$\begin{equation} {\sigma _{\parallel }^2 \atopwithdelims []\sigma _{\perp }^2} = {\sigma _{\parallel,{\rm meas}}^2 \atopwithdelims []\sigma _{\perp,\rm meas}^2} + \frac{\tau _{\rm TJ}}{t} {\sigma _{\parallel,{\rm TJ}}^2 \atopwithdelims []\sigma _{\perp,{\rm TJ}}^2}, \end{equation} \tag{ 5 }$$

where σ²_∥,meas is the sum of the squares of the centering errors of each star parallel to the axis connecting the pair (and similarly for the perpendicular variance), and the remaining terms are as defined in Section 2.1.

Measurement of the offset between the target star (which we will denote with a subscript i = 0) and each of the N reference stars results in a set of N vectors, $\vec{d}_{0i}$. For simplicity, we will write these measured offsets as a single-column vector,

$$\begin{equation} {\bf d} = [ x_{01},\ldots,x_{0N},y_{01},\ldots,y_{0N}]^{\rm T}. \end{equation} \tag{ 6 }$$

The goal of differential astrometry is to use d to determine the position of the target star with respect to the reference grid of stars at each epoch.

There are many possible ways to construct the position of the astrometric target from a given d. Here we use the most general linear combination of the angular offsets, namely

$$\begin{equation} \vec{p} = {\bf Wd}, \end{equation} \tag{ 7 }$$

where W is the 2 × 2N weight matrix, given by

$$\begin{equation} {\bf W} = \left[\begin{array}{@{}cccccc@{}} w_{xx,01} & \cdots & w_{xx,0N} & w_{xy,01} & \cdots & w_{xy,0N} \\ w_{yx,01} & \cdots & w_{yx,0N} & w_{yy,01} & \cdots & w_{yy,0N} \end{array} \right]. \end{equation} \tag{ 8 }$$

We have used the notation w_xy,0i to denote the weighting of the offset from the target star to star i in the y-direction, used to determine the x component of the target's position, $\vec{p}$. For example, for a standard average of the x and y measurements to calculate $\vec{p}$, we would assign all the w_xx,0i = w_yy,0i = 1/N and w_xy,0i = w_yx,0i = 0.

In principle, we are free to assign weights in any manner we please. However, we find it convenient to choose the weights such that they satisfy

$$\begin{eqnarray} &&\sum _i w_{xx,0i} = 1,\qquad \sum _i w_{yy,0i} = 1, \nonumber\\&&[-12pt]\\ &&\sum _i w_{xy,0i} = 0,\qquad \sum _i w_{yx,0i} = 0.\nonumber \end{eqnarray} \tag{ 9 }$$

These constraints ensure that the components of $\vec{p}$ have physical units (e.g., pixels or arcseconds) and that its components are measured in the same coordinate system as d (presumably the detector coordinates). As a consequence, $\vec{p}$ represents the position of the target star in the sense that a proper motion of the target, $\vec{\epsilon }$, with respect to the fixed grid between two epochs will cause a change, $\vec{p} \rightarrow \vec{p} + \vec{\epsilon }$.

In order to determine if any change in $\vec{p}$ over time is meaningful, we must understand its statistical properties. Both differential tilt jitter and measurement errors are assumed to follow Gaussian statistics, so that each instance of target-reference grid offset measurements, d, is drawn from a multivariate normal probability distribution:

$$\begin{equation} P({\bf d}) = \frac{1}{\sqrt{2\pi \det {\boldmath \Sigma} _{\bf d}}} \exp \left(\frac{1}{2} [{\bf d}-\bar{{\bf d}}]^{\rm T} {\boldmath \Sigma} _{\bf d}^{-1} [{\bf d}-\bar{{\bf d}}]\right), \end{equation} \tag{ 10 }$$

where Σ_d is the covariance matrix and the bars above symbols denote using the average value of each matrix entry.

The statistics of $\vec{p}$ follow in a straightforward manner from Equation (10), given our choice in Equation (7). Since $\vec{p}$ is a linear function of d, each $\vec{p}$ is also drawn from a multivariate normal probability distribution with the covariance matrix

$$\begin{equation} {\boldmath \Sigma_{\bf p}} = {\bf W}^{\rm T}{\boldmath \Sigma} _{\bf d}{\bf W}, \end{equation} \tag{ 11 }$$

and the uncertainties of $\vec{p}$ are described by the eigenvectors and eigenvalues of Σ_p. Thus, our goal of optimally determining the target's position requires calculating the covariance matrix, Σ_d, from data or theory, and choosing W to minimize the eigenvalues of Σ_p.

We have chosen to measure positions and offsets in the Cartesian coordinates of the detector, so the form of the covariance matrix, given our above definitions, is

$$\begin{equation} {\fontsize{6}{8}{\hbox{${\boldmath \Sigma} _{\bf d} = \left(\begin{array}{@{}cccccc} \langle (\Delta x_{01})^2\rangle & \cdots & \langle (\Delta x_{01})(\Delta x_{0N})\rangle & \langle (\Delta x_{01})(\Delta y_{01})\rangle & \cdots & \langle (\Delta x_{01})(\Delta y_{0N})\rangle \\ & \ddots & \vdots & \vdots & \ddots & \vdots \\ & & \langle (\Delta x_{0N})^2\rangle & \langle (\Delta x_{0N})(\Delta y_{01})\rangle & \cdots & \langle (\Delta x_{0N})(\Delta y_{0N})\rangle \\ & & & \langle (\Delta y_{01})^2\rangle & \cdots & \langle (\Delta y_{01})(\Delta y_{0N})\rangle \\ & & & & \ddots & \vdots \\ {\rm symmetric} & & & & & \langle (\Delta y_{0N})^2\rangle \end{array} \right),$}}} \end{equation} \tag{ 12 }$$

where we have written as $\Delta x_{ij} \equiv (x_{ij} - \bar{x}_{ij})$ to simplify the notation (likewise for y).

The total covariance matrix has contributions from centering errors and differential atmospheric tilt jitter. Since these contributions are independent, the total covariance matrix can be written as Σ_d = Σ_meas + Σ_TJ, and each term can be derived separately.

4.2.1. The Covariance Matrix for Measurement Noise

In the absence of differential tilt jitter, it is straightforward to construct the covariance matrix for measurement noise alone, Σ_meas. The diagonal terms can be written as

$$\begin{equation} \big\langle \Delta x_{0i}^2\big\rangle \equiv \sigma _{x,0i}^2 = \sigma _{x,0}^2 +\sigma _{x,i}^2, \end{equation} \tag{ 13 }$$

where σ_x,i and σ_x,0 are the uncertainties in determining the x-position of star i and the target star, respectively. For the off-diagonal terms 〈Δx_0iΔx_0j〉, we can use the fact that

$$\begin{eqnarray} \langle \Delta x_{0i} \Delta x_{0j}\rangle &=& \case{1}{2}\big\langle \big\{ \Delta x_{0i}^2 + \Delta x_{0j}^2 - [\Delta x_{0i}- \Delta x_{0j}]^2\big\} \big\rangle \cr &=& \case{1}{2}\big\{ \big\langle \Delta x_{0i}^2\big\rangle + \big\langle \Delta x_{0j}^2\big\rangle - \langle [\Delta x_{0i}- \Delta x_{0j}]^2\rangle\!\big\} \cr &=& \case{1}{2}\big\{ \big\langle \Delta x_{0i}^2\big\rangle + \big\langle \Delta x_{0j}^2\big\rangle - \big\langle \Delta x_{ij}^2\big\rangle \big\} \cr &=& \case{1}{2} \big(\sigma _{x,0i}^2 + \sigma _{x,0j}^2 - \sigma _{x,ij}^2\big)\cr &=& \sigma _{x,0}^2, \end{eqnarray} \tag{ 14 }$$

where we have used only algebra and the above definitions. Equation (14) is the obvious result of the fact that the measurements of the target star's coordinates are common to all differential measurements, and so its uncertainty appears in all the off-diagonal covariance terms, 〈Δx_0iΔx_0j〉 and 〈Δy_0iΔy_0j〉. However, the cross-terms involving both x and y (e.g., 〈Δx_0iΔy_0j〉) vanish because σ_x,0 and σ_y,0 are uncorrelated for measurement noise alone.

4.2.2. The Covariance Matrix for Differential Tilt Jitter

The covariance matrix for differential atmospheric tilt jitter between a pair of stars is diagonal when written in an orthogonal coordinate system with one axis lying along the separation axis of the binary. From Equation (1), we see that it can be written as

$$\begin{eqnarray} {\boldmath\Sigma}_{\rm pair} &=& \left(\begin{array}{@{}cc@{}} \langle (d_\parallel - \bar{d}_\parallel)^2\rangle & \langle (d_\parallel - \bar{d}_\parallel)(d_\perp - \bar{d}_\perp)\rangle \nonumber\\ \langle (d_\parallel - \bar{d}_\parallel)(d_\perp - \bar{d}_\perp)\rangle & \langle (d_\perp - \bar{d}_\perp)^2\rangle \end{array} \right) \nonumber\\ &=& \left(\begin{array}{@{}cc@{}} \sigma ^2_{\parallel,{\rm TJ}} & 0 \\ 0 & \sigma ^2_{\perp,{\rm TJ}} \end{array} \right), \end{eqnarray} \tag{ 15 }$$

where d_∥ and d_⊥ are the angular offsets parallel and perpendicular to the axis connecting the pair of stars, respectively.

For a general field of N stars, no coordinate system exists that diagonalizes the full tilt jitter covariance matrix, Σ_TJ. But we can begin computing the entries by rotating Σ_pair into our x–y coordinates via R^TΣ_pairR, where

$$\begin{equation} {\bf R} = \left(\begin{array}{@{}cc@{}} \cos \phi & \sin \phi \\ -{\rm sin}\, \phi & \cos \phi \end{array} \right). \end{equation} \tag{ 16 }$$

The result is

$$\begin{eqnarray} &&{\bf R}^{\bf T}{\boldmath\Sigma}_{\rm pair}{\bf R} = \nonumber\\ &&\!\!\left(\begin{array}{@{}cc@{}} \sigma _{\parallel,0i}^2\cos ^2\phi _{0i} + \sigma _{\perp,0i}^2\sin ^2\phi _{0i} & \big(\sigma _{\parallel,0i}^2 - \sigma _{\perp,0i}^2\big)\cos \phi _{0i}\sin \phi _{0i} \\ \big(\sigma _{\parallel,0i}^2 - \sigma _{\perp,0i}^2\big)\cos \phi _{0i}\sin \phi _{0i} & \sigma _{\parallel,0i}^2\sin ^2\phi _{0i} + \sigma _{\perp,0i}^2\cos ^2\phi _{0i} \end{array} \right)\!,\nonumber\\ \end{eqnarray} \tag{ 17 }$$

where ϕ_0i is the angle between $\vec{d}_{0i}$ and our arbitrary Cartesian system measured counterclockwise from the x-axis, and we have introduced the notation that the uncertainty parallel to $\vec{d}_{ij}$ is σ_∥,ij and the uncertainty orthogonal to $\vec{d}_{ij}$ is σ_⊥,ij as calculated from Equation (1). Thus, we can identify the diagonal terms

$$\begin{eqnarray} &&\ \big\langle \Delta x_{0i}^2\big\rangle = \sigma _{\parallel,0i}^2\cos ^2\phi _{0i} + \sigma _{\perp,0i}^2\sin ^2\phi _{0i},\nonumber\\&&[-6pt]\\ &&\big\langle \Delta y_{0i}^2\big\rangle = \sigma _{\parallel,0i}^2\sin ^2\phi _{0i} + \sigma _{\perp,0i}^2\cos ^2\phi _{0i}\nonumber \end{eqnarray} \tag{ 18 }$$

and

$$\begin{equation} \langle \Delta x_{0i} \Delta y_{0i}\rangle = \big(\sigma _{\parallel,0i}^2 - \sigma _{\perp,0i}^2\big)\cos \phi _{0i}\sin \phi _{0i}. \end{equation} \tag{ 19 }$$

For the off-diagonal terms 〈Δx_0iΔx_0j〉, we notice that (as used in Equation 14)

$$\begin{eqnarray} \langle \Delta x_{0i} \Delta x_{0j}\rangle &=& \case{1}{2}\big\langle \big\{ \Delta x_{0i}^2 + \Delta x_{0j}^2 - [\Delta x_{0i}- \Delta x_{0j}]^2\big\} \big\rangle \cr &=& \case{1}{2}\big\{ \big\langle \Delta x_{0i}^2\big\rangle + \big\langle \Delta x_{0j}^2\big\rangle - \langle [\Delta x_{0i}- \Delta x_{0j}]^2\rangle\!\big\} \cr &=& \case{1}{2}\big(\sigma _{\parallel,0i}^2\cos ^2\phi _{0i} + \sigma _{\perp,0i}^2\sin ^2\phi _{0i} + \sigma _{\parallel,0j}^2\cos ^2\phi _{0j} \nonumber\\ && + \sigma _{\perp,0j}^2\sin ^2\phi _{0j} - \sigma _{\parallel,ij}^2\cos ^2\phi _{ij}\nonumber\\ && - \sigma _{\perp,ij}^2\sin ^2\phi _{ij}\big), \end{eqnarray} \tag{ 20 }$$

where, in the last step, we have used the fact that x_ij = x_0i − x_0j and the relations in Equations (18) and (19). The quantities 〈Δy_0iΔy_0j〉 can be obtained by interchanging sine and cosine in Equation (20).

For the remaining off-diagonal terms 〈Δx_0iΔy_0j〉, we can use the fact that

$$\begin{eqnarray} \langle \Delta x_{0i}\Delta y_{0j}\rangle &=&\langle \Delta x_{0i}[\Delta y_{0i}+\Delta y_{ij}]\rangle \cr &=&\langle \Delta x_{0i}\Delta y_{0i}\rangle + \langle \Delta x_{0i}\Delta y_{ij}]\rangle \cr &=&\langle \Delta x_{0i}\Delta y_{0i}\rangle + \langle [\Delta x_{0j}-\Delta x_{ij}]\Delta y_{ij}\rangle \cr &=&\langle \Delta x_{0i}\Delta y_{0i}\rangle + \langle \Delta x_{0j}\Delta y_{ij}\rangle - \langle \Delta x_{ij}\Delta y_{ij}\rangle \cr &=&\langle \Delta x_{0i}\Delta y_{0i}\rangle + \langle \Delta x_{0j}[\Delta y_{0j}-\Delta y_{0i}]\rangle - \langle \Delta x_{ij}\Delta y_{ij}\rangle \cr &=&\langle \Delta x_{0i}\Delta y_{0i}\rangle + \langle \Delta x_{0j}\Delta y_{0j}\rangle - \langle \Delta x_{ij}\Delta y_{ij}\rangle\nonumber\\ && - \langle \Delta x_{0j}\Delta y_{0i}\rangle. \end{eqnarray} \tag{ 21 }$$

Rearranging gives

$$\begin{eqnarray} &&\hspace{-1.8pc}\langle \Delta x_{0i}\Delta y_{0j}\rangle + \langle \Delta x_{0j}\Delta y_{0i}\rangle =\langle \Delta x_{0i}\Delta y_{0i}\rangle\nonumber\\ &&\hspace{-1pc}\quad +\, \langle \Delta x_{0j}\Delta y_{0j}\rangle - \langle \Delta x_{ij}\Delta y_{ij}\rangle. \end{eqnarray} \tag{ 22 }$$

All the terms on the right-hand side are known from Equation (19), and further investigation shows that the two terms on the left-hand side are equal. So, Equations (13), (14), (18)–(20), and (22) contain all the information required to construct the full covariance matrix, Σ_d.

The optimal choice of weights in Equation (3) are those that minimize the eigenvalues in Equation (11). For a 2 × 2 symmetric matrix, the sum of the eigenvalues is the trace of the matrix, so our problem becomes one of minimizing the trace of Σ_p subject to the constraints in Equation (9). Specifically, we will use the method of Lagrange multipliers (Betts 1980) to find the optimal weights, W', that minimize the quadratic equation

$$\begin{equation} {\rm Tr(}{\boldmath \Sigma} _{\bf p}) = \case{1}{2} {\bf W}^{\rm \prime T} \,{\bf SW}^{\prime }, \end{equation} \tag{ 23 }$$

where

$$\begin{equation} {\bf S} = \left[\begin{array}{@{}cc@{}} {\boldmath\Sigma} _{\bf d} & 0 \\ 0 & {\boldmath \Sigma} _{\bf d} \end{array} \right] \end{equation} \tag{ 24 }$$

is a 4N × 4N matrix, and

$$\begin{eqnarray} {\bf W}^\prime &=& [w_{xx,01},\ldots, w_{xx,0N}, \ w_{xy,01},\ldots,w_{xy,0N},\nonumber\\ && \ w_{yx,01},\ldots,w_{yx,0N},\ w_{yy,01},\ldots,w_{yy,0N}]^{\rm T} \end{eqnarray} \tag{ 25 }$$

is a vector of length 4N. Note that W' has identical entries as W in Equation (8); it is just written as a single vector to cast the minimization problem into a single Equation (23). We want to find the extrema of Equation (23) subject to the linear constraints in Equation (9), which can be written as

$$\begin{equation} {\bf CW}^\prime = {\bf V}, \end{equation} \tag{ 26 }$$

where we define the 4 × 4N symmetric matrix

$$\begin{equation} {\bf C} = \left[\begin{array}{@{}cccccccccccc@{}} 1 & \cdots & 1 & 0 & \cdots & 0 & 0 & \cdots & 0 & 0 & \cdots & 0\\ & & & 1 & \cdots & 1 & 0 & \cdots & 0 & 0 & \cdots & 0\\ & & & & & & 1 & \cdots & 1 & 0 & \cdots & 0\\ {\rm sym} & & & & & & & & & 1 & \cdots & 1 \end{array} \right], \end{equation} \tag{ 27 }$$

and

$$\begin{equation} {\bf V} = \left[\begin{array}{@{}c@{}} 1 \\ 0 \\ 0 \\ 1 \end{array} \right]. \end{equation} \tag{ 28 }$$

In this framework, the optimal weights are those that solve the system of linear equations:

$$\begin{equation} \left[\begin{array}{@{}cc@{}} {\bf S} & {{\bf C}^{\rm T}}\\ {\bf C} & {\bf 0} \end{array} \right] \left[\begin{array}{@{}c@{}} {{\bf W}^\prime} \\ {\lambda} \end{array} \right] = \left[\begin{array}{@{}c@{}} {\bf 0}\\ {\bf V} \end{array} \right]. \end{equation} \tag{ 29 }$$

Here, λ are the Lagrange multipliers, which will not be used further. Equation (29) can be solved via a matrix inversion.

Note that the general constraints on the weights we have written in Equations (9) and (26) have two somewhat unintuitive features. The first is that the y measurements are sometimes used to compute the x position and vice versa. The other property is that they allow for negative weights, meaning that in some cases, certain measurements will be subtracted in calculating the position of the astrometric target, $\vec{p}$. These two facts conspire to exploit the natural correlations inherent in the data. The flexible and possibly negative weights essentially allow the reference grid to be symmetrized, thereby using the known correlations to cancel noise so as to minimize the variance in $\vec{p}$.

As indicated in the above analysis, the single-epoch uncertainty in the location, $\vec{p}$, of the target relative to the grid of reference stars is represented by the eigenvalues and eigenvectors of the 2 × 2 matrix Σ_p (Equation 11). This matrix itself depends on the distribution of reference stars, the precision of centering measurements, and the degree of noise correlation due to differential tilt through the matrix Σ_d. In this way, the intrinsic precision of the measured value of $\vec{p}$ depends on these three factors.

To ascertain the behavior of Σ_p with the density of available reference stars, we performed a series of numerical simulations. In each simulation, N (2 ⩽ N ⩽ 100) stars were randomly distributed throughout a 25'' × 25'' FOV. We assumed that the target was a bright star in the middle of the field with a centering error of 0.5 mas and the reference stars were fainter, drawn from a Gaussian distribution with a mean centering error of 2 mas and a standard deviation of 1 mas (somewhat analogous to the situation for the guide star in M5; see Section 5). The full covariance matrix, Σ_d, was computed for each stellar configuration assuming these centering errors, the typical turbulence profile above Palomar Observatory, and a 1.4 s exposure time.

In the first simulation, Σ_d was contracted as in Equation (11) using standard averaging for W (w_xx,0i = w_yy,0i = 1/N, w_xy,0i = w_yx,0i = 0). For the second simulation, Σ_d was contracted using the optimal W as calculated by using the prescription in Section 4.3. In each case, the geometric mean of the two eigenvalues of the resulting matrix, Σ_p, was computed to form an estimate of the single-epoch measurement precision of $\vec{p}$. To average away random effects arising from the particular geometry of the random distribution of stars, each numerical simulation was repeated for 100 random distributions of stars for each value of N, and these were averaged to generate a mean value for the single-epoch measurement precision.

The resulting values for the single-epoch measurement precision of $\vec{p}$ are shown in Figure 2 as a function of the number of reference stars, along with the contributions of measurement noise and differential tilt jitter. In both simulations, the error due to measurement noise decreases as N^−0.3. However, in the limit of an infinite number of reference stars, this error asymptotes to the target star's measurement error. The rate at which the measurement noise decreases to this value depends on the distribution of reference star measurement errors.

Zoom In Zoom Out Reset image size

The important distinction between the two simulations is the contribution of tilt jitter to astrometric performance. In the simulation utilizing standard averaging, there is very little gain with increased stellar density (N^−0.15), and tilt jitter dominates the error budget. However, the optimal estimation algorithm rapidly (N^−0.7) eliminates the contribution of differential tilt by taking advantage of the correlations inherent in Σ_d and the flexibility to symmetrize the reference field through the choice of weights.

In order to test our expectation that tilt jitter dominates the astrometric error, we calculate the rms of the angular offsets for pairs of stars in the field (Figures 1 and 3). These results clearly show the characteristic signature of differential tilt. Namely, the rms separation along the axis connecting the two stars is larger than that of the perpendicular axis by a factor of $\approx \sqrt{3}$. However, the magnitude of the tilt jitter is smaller than the theoretical expectations, which suggests that some of the tilt jitter has been averaged away in the 1.4 s exposure time.

Zoom In Zoom Out Reset image size

We have no direct measurement of the wind speed profile over the telescope to calculate the expected tilt jitter timescale. Instead, we fit the observed σ²_ij and angular offsets using the model in Equation (5) with t = 1.4 s. The best-fit values are σ_meas,ij ≈ 2 mas and t/τ_TJ ≈ 7. This implies that the characteristic timescale for tilt jitter is ≈0.2 s, resulting in a wind crossing time of 25 m s⁻¹. Turbulence at higher altitudes contributes most to the differential atmospheric tilt jitter, and this velocity is typical of wind speeds in the upper atmosphere (Greenwood 1977). It is also clear from the figure that a number of stars have measurement noise that is much less than 2 mas; thus, this number should only be taken as characteristic of the faint stars.

The astrometric precision achieved in a single epoch is an important diagnostic of the measurement model. On a given night for a given star, we investigate the use of both the ≈500 images and the theory in Section 4.2 to calculate Σ_d, leading to the optimal weights. We then apply this weight matrix to the measured offsets to compute the target's position in each image, resulting in a time series in each component of $\vec{p}$ for each epoch. The properties of each time series are best explored by computing its Allan deviation (also known as the square root of the two-sample variance). The Allan deviation is calculated by dividing a time series into chunks, averaging each segment, and computing the rms of the resulting, shorter time series. If the time series is dominated by random errors, its Allan deviation will decrease as $1/\sqrt{t_{\rm avg}}$, where t_avg is the length of each chunk. It is also necessary to have sufficiently many segments so that an rms calculation is meaningful. Here, the longest timescale probed is ≈2–3 minutes for each 10–15 minute time series.

We compute the geometric mean of the Allan deviation in each dimension as a function of the averaging time for the AO guide star in Figure 4 after computing the covariance matrix from data. After 1.4 s, the guide star's positional precision is ≈600 μas. The precision subsequently improves as t^{−0.51±0.08} to ≈70 μas after 2 minutes, and has yet to hit a systematic floor. This suggests a precision of ≈30 μas for the full 10–15 minute data set, assuming that no systematic limit is reached in the interim.

Zoom In Zoom Out Reset image size

This level of precision is not limited to the AO guide star; similar performance is obtained on other stars in the core of M5. In Figure 5, we show the astrometric precision obtained on 2007 May 29 after 2 minutes for all detected stars as a function of their K_s magnitude. Precision below 100 μas is achieved on targets as faint as K_s ≈ 13 mag using a narrowband filter and 1.4 s individual exposures. This demonstrates the substantial S/N benefit afforded by AO.

Zoom In Zoom Out Reset image size

The astrometric precision shown in Figure 5 resulting from the theoretically-determined covariance matrix and optimal weights is ≈300 μas after 2 minutes for stars with K_s ≲ 13 mag. This level of precision is substantially better than the performance of simpler weighting schemes, but it is a factor of 2–4 worse than using the data to calculate the covariance matrix and weighting. There are several possible reasons for this reduction in precision. First, we have only used estimates of the measurement noise for each star used to calculate Σ_meas. Second, the turbulence profile used to construct Σ_TJ is estimated from the average C²_n(h) seen by the DIMM/MASS. This unit is located 300 m from the Hale telescope and uses Polaris to estimate the turbulence profile. As a consequence, there could be important differences between the measured atmospheric turbulence and that encountered by the light from M5. Finally, we have not attempted to capture the time variability of the turbulence, having used only the average values.

In Figure 6, we investigate the improvement of the AO guide star astrometry with the number of reference stars. We drew random subsets of the available grid stars, computed Σ_d from the data, calculated the optimal weights, and showed the geometric mean of the eigenvalues of Σ_p. To average over the geometry of a particular draw, we repeated this process ten times for each value of N and averaged the results. We see that the precision rapidly decreases as N^{−0.60±0.03}. This is slightly faster than what our simulations predict for 1.4 s of integration time. However, as noted above, our simulations are meant to approximate M5, but do not capture the true distribution of stellar measurement errors (which are difficult to decouple from tilt jitter) or any evolution in atmospheric turbulence during the observation.

Zoom In Zoom Out Reset image size

The goal of astrometry is to measure the position of the target star over many epochs. Astrometrically interesting timescales range from hours to years. Clearly, the optical systems must be stable over these spans for astrometry with AO to be viable. There are several obstacles that could render the single-epoch precision obtained in Section 5.2 meaningless. For example, PHARO is mounted at the Cassegrain focus, which results in flexure of the instrument as the telescope tracks, and undergoes warming and cooling cycles between observing periods (typically twice per month) that could cause small changes in the powered optics. Either of these facts could alter the geometric distortion and make astrometric measurements unrepeatable. In order to probe the system stability, we have designed our experiment to be as consistent as possible, and it has spanned many removal and reinstallations of PHARO over 2 months.

In order to investigate the accuracy of the M5 measurements, we first measured and corrected the small rotational (≲004) and plate scale (≲10⁻⁵) changes between the May 29 and July 22 data and the May 28 images. We also calculate the optimal weights for a given star on all three nights, and average them to create one weighting matrix to use for each epoch. This is not strictly optimal, since each night has different turbulence conditions for example, but it ensures that the scenario that $\vec{p} \rightarrow \vec{p} + \vec{\epsilon }$. In Figure 7, we see that the measured position of the AO guide star is accurate from epoch to epoch at ≈100 μas. The error ellipses are those estimated by continuing to extrapolate the precision found in Figure 4 by $1/\sqrt{t}$ to the full 10–15 minute time series. This is an impressive level of accuracy, but unfortunately is a factor of 3 worse than our expectation. It suggests that there is some instability, likely in the distortion, over the 2 months that limit the astrometric accuracy.

Zoom In Zoom Out Reset image size

The other stars in M5 show a similar level of astrometric accuracy (Figure 8) up to K_s ≈ 13 mag. This limit can be certainly pushed to be considerably fainter with increased integration time or a larger aperture. The achievement of such high levels of astrometric performance on faint targets, given the modest time investment, short integration time, and narrowband filters, illustrates the substantial S/N gain and potential for astrometry enabled by AO.

Zoom In Zoom Out Reset image size