THE M31 VELOCITY VECTOR. I. HUBBLE SPACE TELESCOPE PROPER-MOTION MEASUREMENTS

The data set we used for measuring the PM of M31 consists of HST images of three fields obtained in two separate epochs. The fields are the SPHEROID, OUTERDISK, and TIDALSTREAM fields of Brown et al. (2006) obtained under the science programs GO-9453 and GO-10265 (PI: T. Brown). Figure 1 shows the location of the three M31 target fields. The observation details for the first-epoch data are given in Brown et al. (2006); field coordinates and total exposure times are listed in their Table 1. In brief, the target fields were observed between 2002 December and 2005 January in two filters (F606W and F814W) with HST ACS/WFC to create color–magnitude diagrams (CMDs) that reach well below the M31 main-sequence turn-off. For each field, all exposures were obtained within 40 days so they can be safely treated as a single epoch of data. Total per-field per-filter exposure times ranged from 53 to 161 ks. Individual images had exposure times ranging from 1100 to 1300 s and were dithered such that no two exposures placed a star on the same pixel. In addition to the whole-pixel offsets, the dithering strategy also included sub-pixel shifts so that every object was observed at a range of locations relative to the pixel boundaries; such "pixel phase" coverage is critical for our program. The very deep first-epoch F606W data were used by us for constructing super-sampled stacked images, identifying stars and galaxies, building templates for the template-fitting method, and providing reference positions of stars and galaxies, as described in Section 3 below.

Zoom In Zoom Out Reset image size

Our second-epoch data were observed between 2010 January and 2010 August using ACS/WFC and WFC3/UVIS in F606W (HST program GO-11684, PI: van der Marel). We designed our observations to minimize random and systematic errors with the following strategies in mind. Pre-analysis of the first-epoch data revealed that the F606W filter provides a slightly better astrometric handle on the background galaxies than the F814W filter, so we chose to obtain second-epoch images only with F606W. Furthermore, we decided to observe the three target fields with two different instruments instead of just one. This yields comparable random error, but allows additional consistency checks on our PM results. Since the purpose of the second-epoch observations was astrometry, and not deep photometry, the total second-epoch exposure time was much less than in the first epoch. We observed each target field for two orbits (four half-orbit exposures) with WFC3/UVIS and one orbit (two half-orbit exposures) with ACS/WFC. Individual images had exposure times ranging from ∼1300 to 1500 s, slightly longer than in the first epoch, and were dithered in similar fashion. The baselines between the two epochs are in the range of 5–7.5 years.

For the ACS observations, we used the same orientations and coordinates as used by the first-epoch observations to maximize the overlapping area and so that any errors in the static distortion solution would naturally cancel out. Due to different sets of guide stars being used between the first- and second-epoch observations, the second-epoch ACS images were slightly offset and rotated with respect to the first-epoch images. However, the offsets are all within a negligible fraction of the ACS/WFC FOV, and the orientations are within 01 of each other.

For the WFC3 observations, we pointed the telescope so that the resulting images roughly share the same centers with the first epoch, but the WFC3 images were rotated by roughly 45° with respect to the ACS images. This was to place the parallel ACS/WFC fields overlapping with the first epoch parallel WFPC2 images, but the parallel fields will not be discussed in this paper as they are not useful for astrometry at the accuracy that we need. Because the FOV of WFC3/UVIS is 64% of that of the ACS/WFC, the observed WFC3 images are fully contained within the ACS images despite the ∼45° difference in their orientations. The details of our second-epoch data are summarized in Table 1.

Measurement of PMs involves measuring the displacement between the position of an object at one time and its position at another. If this displacement can be measured with respect to objects for which we know the absolute motion, then we can obtain absolute PMs. In our case, we will measure the displacement of M31 stars with respect to the background galaxies in the field to obtain an absolute PM. Our strategy is to first align the stars in the first- and second-epoch images, restricting the alignment to those stars confirmed to belong to M31. Then, we measure the average displacement of the background galaxies with respect to this moving frame of reference. The negative of this relative displacement yields the mean absolute PM of the M31 stars. This measurement of the mean does not require knowledge of the exact PM distribution of the M31 stars. This distribution can be complex, because different structural components of M31 contribute to each field. However, to transform the mean PM of the M31 stars in a field to an estimate of the M31 COM, one does need to construct a model for the internal kinematics of M31, which is addressed in Paper II of this series.

The overall PM derivation process is summarized below.

• 1.  
  We first create a high-resolution stacked image using the deep and well-dithered first-epoch images for each target field, using the distortion-corrected first exposure as the reference frame.
• 2.  
  Stars and galaxies are then identified from the stacked image for each field, and photometric measures and the CMD are used to specifically identify M31 reference stars.
• 3.  
  We construct a template for each star and galaxy from the high-resolution image.
• 4.  
  We then use the template to measure in a consistent fashion a position for each star/galaxy in each individual exposure of each epoch.
• 5.  
  We redefine the first-epoch reference frame using the template-based positions of the stars and determine the average first-epoch position of each galaxy in this frame.
• 6.  
  We determine the template-measured positions of the stars in each of the six second-epoch exposures for each field and use them to transform the template-measured positions of the galaxies into the first-epoch frame.
• 7.  
  We then take the difference between the second- and first-epoch positions of the galaxies to obtain the relative displacement of the galaxies with respect to the M31 stars.
• 8.  
  Finally, we multiply the relative displacements of the galaxies by −1 to obtain the mean absolute displacement of the M31 stars, since in reality the galaxies are stationary and the stars are moving. Division by the time baseline turns the displacements into PMs.

The following sections provide the details of each step.

We downloaded all the required data in the form of _flt.fits images from the STScI HST archive. These images are already bias-subtracted, dark-subtracted, and flat-fielded by the STScI data-reduction pipelines using the best available calibration data at the time of data inquiry. Each multi-extension _flt.fits file was converted to a single 4096 × 4096 image where the top chip scene is directly abutted to the bottom chip scene. This single-format image allows a more efficient way of handling the data, and all computer software we used throughout our analysis were specifically programmed to deal with this image format.

The HST detectors suffer from significant degradation in charge-transfer efficiency (CTE) as time progresses. This impacts the astrometry of objects by effectively shifting the centroid position of the flux distribution. This can be a particularly serious problem for our second-epoch ACS/WFC data, taken almost nine years after the installation of ACS. A pixel-based CTE correction routine has recently been developed by Anderson & Bedin (2010) for the ACS/WFC, and Brown et al. (2010) demonstrated that this correction works very well for astrometric purposes. All first-epoch ACS/WFC images were processed through this routine. For the second-epoch ACS/WFC images, we used a modified version of the CTE correction code that also included a correction for X-CTE, in addition to the publicly available correction for Y-CTE. We found that the X-CTE correction made very little difference to the positions measured in the images, compared to the other uncertainties in the results.

A CTE correction routine was not available for the WFC3/UVIS at the time of writing this paper, so no correction was made to the WFC3 images. However, since at the time of the observations, UVIS had spent only ∼10% as much time as ACS in the space environment, the CTE degradation of our UVIS images is expected to be quite small. Although UVIS appears to be losing CTE faster than ACS did, the studies that show this have focused on low-background images (e.g., Baggett et al. 2011). As our backgrounds are quite high, we expect CTE losses to be very low for our images. Hence, the CTE degradation should only have a minimal effect on our WFC3 PM measurements. A direct post-facto consistency check of this is provided by comparison of our ACS/WFC and WFC3/UVIS PMs derived for the same fields (see Section 4.2 below).

As the final step of the initial processing of individual images, we measured each star using the img2xym_WFC_9x10 program (Anderson & King 2006) on the ACS/WFC images, and using a similar program on the WFC3/UVIS images. Both programs utilize library PSFs to determine a position and a flux for each star. We will use these PSF-based star positions to provide our initial handle on the coordinate transformations from one image to another.

The positions of objects on the HST detectors must be corrected for geometric distortions before they can be used for any astrometric measurement. For the ACS/WFC detector, the distortion-correction solutions by Anderson (2006) have been used in several astrometric studies. These solutions are known to be better than 0.01 pixel or ∼0.5 mas and were shown to be stable between 2002 and 2007 (Anderson 2007; van der Marel et al. 2007). These available solutions were directly applied to the positions measured in our first-epoch ACS/WFC images.

While reducing data for another program (GO-11677, PI: H. Richer), one of us (J.A.) found that the ACS distortion solution appears to have changed by ∼0.005 pixel relative to the solution before the HST Service Mission 4 (SM4). This may be related to a similar variation of the PSF (J. Anderson et al., in preparation). It is unclear what caused these variations, but they appear to be stable over time. We used the GO-11677 data set to develop a correction to the pre-SM4 ACS/WFC distortion solution and applied it to the second-epoch data. For the WFC3/UVIS data, we used the distortion corrections provided by Bellini et al. (2011). These corrections are better than 0.008 pixel or ∼0.3 mas and appear to be stable over the UVIS lifetime so far.

The accuracies of the distortion corrections are comparable to our final measured motions (as will be shown later). However, geometric distortion affects both star and galaxy positions similarly, and hence drops out in a differential measurement to lowest order. This is true in particular if each galaxy position is measured with respect to that of nearby stars that fall on the same area of the detector. This is what we do in Section 3.8 below, to ensure that higher order geometric distortion correction residuals do not affect our final results.

The deep and well-dithered first-epoch ACS/WFC data are well-suited for constructing high-resolution stacked images. To do this, we first apply the distortion corrections (see Section 3.3) to the positions of stars measured in Section 3.2. We then adopt the first exposure of each field as the frame of reference, cross-identify stars in this exposure and stars in the other exposures, and use their distortion-corrected positions to construct a six-parameter least-squares linear transformation between the two frames. The six parameters involve x–y translation, scale, rotation, and two components of skew. The need for these linear transformations is discussed in Section 3.6.4 of Anderson & van der Marel (2010). We use these transformations to convert the star positions measured in the various images into the reference frame, giving us many estimates of the position of each star in the reference frame. We average these positions to refine the reference-frame positions for all the stars and use these new average positions to improve the transformations (which had initially been based on the positions as measured in the first frame itself). This procedure is iterated two to three times to improve the positions of both the stars and the transformations.

We use the star-based transformations to construct the stacked images. In order to get better sampling, we super-sample the stacked image by a factor of two relative to the native ACS/WFC pixel scale. The image-stacking process we used is similar to the commonly used Drizzle algorithm (Fruchter & Hook 2002) with a point kernel, except that we included an iterative procedure to regularize the sampling in a manner similar to iDrizzle (Fruchter 2011). The procedure involves no deconvolution, and as such the resulting image at every point simply represents the flux that an actual _flt image pixel would receive if it were centered at that location in the frame. It is this property that allows us to construct empirical templates (see Section 3.6.1) for our stars and galaxies. Figure 2 shows a 25'' × 25'' portion of the 2× super-sampled stacked image for the SPHEROID field with M31 stars and background galaxies identified, as described in the following sections.

Zoom In Zoom Out Reset image size

All of our target fields are dominated by M31 stars, but there are also other sources, such as foreground stars and background galaxies. Since our strategy for deriving the PM of M31 is to measure the displacement of the background galaxies relative to the co-moving frame of reference defined by the M31 stars, it is important to accurately identify both background galaxies and M31 stars before we proceed any further.

3.5.1. Identifying M31 Reference Stars

The selection of M31 stars was carried out following the procedure below. For each target field, we use the star list compiled in Section 3.4. Our initial star lists include only sources that are found independently in a large number of first-epoch F606W exposures (typically, >35%); as such, the list is almost entirely free of cosmic rays and image defects. There are, however, a large number of resolved sources, which could be either blended stars or extended background galaxies. To identify them, we made use of the mean quality-of-PSF-fit parameter (QFIT) reported by the img2xym_WFC_9x10 program (see Figure 3(a)). Selecting the objects with low values of QFIT allows us to filter out extended objects as well as stars that are too close to other stars to provide accurate position measurements. To ensure that all of the stars in the list are M31 members, we reduced the F814W images and constructed a CMD (see Figure 3(b)). From this, we selected stars that lie roughly on the M31 sequences. Although small in number, field stars in the Milky Way halo may be included in our selection, but we believe that most of them will be filtered out in our next step of selection. Finally, reference stars were selected based on their lack of motion with respect to the other M31 stars. This was done by aligning the second-epoch star positions with the first-epoch positions and iteratively rejecting the objects that have moved between the epochs (see Figure 3(c)). For any given ACS/WFC target field, the M31 stars should all be moving toward the same direction in space. Differential motions due to the internal kinematics of M31 within a single ACS/WFC field should be negligible compared to the observational uncertainties (see Paper II for details). We note that, in principle, better detection and photometry of fainter stars can be achieved by measuring stars directly from the stacked images, as has been demonstrated by Brown et al. (2006). However, as our main goal is doing astrometry, we are only interested in stars for which positions can be reliably measured in individual exposures. Our final lists of M31 reference stars include ∼10, 000 stars for the SPHEROID and OUTERDISK fields, and ∼5000 stars for the TIDALSTREAM field.

Zoom In Zoom Out Reset image size

3.5.2. Identifying Background Galaxies

Our goal is to measure accurate PMs. The most crucial part of our analysis is therefore to measure the positions of objects in the individual images as accurately as possible. Simple methods such as fitting the objects with a two-dimensional Gaussian function or finding the flux centroid are accurate only if the objects have well-sampled peaks. This is not the case for most of the background galaxies in our target fields. To measure the positions of objects, we therefore adopt the template-fitting method developed by Mahmud & Anderson (2008).

3.6.1. The Template-fitting Method

The details and general diagnostics of the template-fitting method are documented in Mahmud & Anderson (2008). Here, we summarize only the basic concepts behind the method. The distortion corrections in Section 3.3 and the linear transformations derived in Section 3.4 allow us to associate a position (x_r, y_r)_j in a given individual image j with a position (x_m, y_m) in the stacked image, and vice versa. We are thus able to build a model of what the galaxy or star should look like in an individual image by proper sampling of the stacked image.

The goal is to measure in consistent fashion a position for each object (star or galaxy) in each individual exposure. To do this, we construct a template model for each object, and use it to measure that object in every exposure. The location within the template that we will define to be its center is arbitrary, since in the end we care only about differences in position. So we define the center of the brightest pixel in the stacked image to be the object center (or "handle" in the parlance of Mahmud & Anderson). With the center so defined, we interpolate the stacked image at a super-sampled array of points, with the center at (0,0) but with a pixel orientation and spacing that corresponds to the transformed exposure's _flt coordinate system. The result is a template that can be interpolated at an array of locations to tell us how much flux we would expect in an array of _flt pixels for a presumed object center. The template is custom made for each individual exposure, based on the mapping of that exposure's pixel coordinate system into the master frame. The template centers are all the same.

For each object in each exposure, we use the previously established transformations to identify the location of the object in the exposure to within 1 pixel. We then evaluate an array of trial locations for the center, covering an area of 1 × 1 pixels with a spacing of 0.01 pixel. At each trial location, interpolation of the template with a bicubic interpolation algorithm tells us how much flux we should expect to measure in each pixel of the exposure. We compare this model to the actual observed pixel values to obtain a quality-of-fit quantity,

$$\begin{eqnarray} &&qfit_{ij} = {{\Sigma _{\rm pixels} \: | {\rm OBS} - {\rm MODEL}_{ij} | } \over {\Sigma _{\rm pixels} \: {\rm OBS}}}, \end{eqnarray} \tag{ 1 }$$

for each trial offset (i, j). The position of the galaxy or star then corresponds to the location of the template center that gives the minimum qfit_ij.

In Section 3.2, we used empirical "library" PSF fits to the stars to define the transformations for construction of the stacked images. However, for the final PM measurements, we use the template procedure for both galaxies and stars. This minimizes any potential for systematic errors due to differences in measurement techniques between different types of objects. The sums in Equation (1) were chosen to extend over the 5 × 5 pixel raster centered on the object's brightest pixel in the individual exposure. The size of this region was deliberately chosen to be small, since most of the positional information for the stars and compact galaxies of interest here is contained in the central core pixels.

3.6.2. First-epoch ACS Data

We applied the template-fitting method to all stars and galaxies in the first-epoch ACS/WFC data to determine their positions in each exposure. To illustrate the quality of the template-fitting method, we show in Figure 4 the best-fit qfit (as defined in Section 3.6.1) versus instrumental magnitude for the background galaxies in one of the first-epoch images. The top row of Figure 5 shows the data, template fits, and residuals for three galaxies of different brightness, and for one star.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.6.3. Second-epoch ACS Data

The stacked images constructed in Section 3.4 represent the astronomical scene convolved with the average effective PSF appropriate for the first-epoch data. The effective PSF represents the convolution of the instrumental PSF with the sensitivity profile of a pixel. The astronomical scene for any given object should not change between our epochs of data, but the effective PSF might change. This is certainly the case for the second-epoch WFC3/UVIS data, discussed below, since both the instrumental PSF and the pixelization are different than for the first-epoch ACS/WFC data. However, we also found that the effective PSF for the second-epoch ACS/WFC data was not the same as for first-epoch ACS/WFC data. Even though the pixelization remained the same, we found that a subtle change occurred in the ACS PSF since SM4 (see Section 3.3).

It is important that any change in the effective PSF be explicitly accounted for in the analysis. PSF changes can introduce centroid shifts, and these shifts can be different for point sources and for extended sources. Our technique for measuring the M31 PM relies on positional differences between stars and background galaxies. PSF changes can therefore produce spurious PM results.

In order to deal with the differences in the effective PSF between epochs, we model each object in a given second-epoch image as follows. First, we select bright and isolated stars from the M31 reference star list compiled in Section 3.5.1, and use them to derive a 7 × 7 pixel kernel that accounts for differences in effective PSF between the stacked image and the individual exposure. Next, we interpolate the first-epoch template onto the _flt pixel grid of the second-epoch image, as we did in Section 3.6.2. We then convolve the template with the kernel derived in the earlier stage and finally evaluate the goodness-of-fit quantity qfit to find the object position that provides the minimum qfit. The middle row of Figure 5 shows the data, template fits, and residuals for three galaxies of different brightness, and for one star, obtained with this procedure for one of the second-epoch ACS/WFC images (jb4404vuq).

The complexity in this procedure lies in the determination of the kernel. Each kernel has 49 free parameters. However, it is constrained by 5 × 5 pixel patches centered around all the objects in each image. This is an overdetermined problem that can be solved to find the best-fitting kernel. We cast the problem into the form of a least-squares matrix equation and solved for the kernel using the observed pixel values for several thousand bright and isolated stars distributed throughout the target fields. The solution was obtained through singular value decomposition (Press et al. 1992). To fit the kernels, we first determined the stellar positions of the bright and isolated stars without any kernel. Then, we kept those positions fixed and optimized the kernels.

We ended up using a single kernel for each second-epoch image. We experimented with using different kernels for different parts of each image (to deal with potential field-dependent PSF variations), but found that this did not change the main results. In deriving the kernel, we recognized that there is a degeneracy between shifting the kernel and shifting the astronomical scene. We removed this degeneracy by constructing the pixel residuals that went into the kernel with respect to the library-PSF-measured position (see Section 3.2) for each star. This ensured that the kernel would reproduce, on average, the PSF-based positions for the stars, and that it would properly account for changes in source shape due to known PSF variations between epochs. Either way, the degeneracy does not affect absolute PMs that are based on relative differences in position between stars and background galaxies.

The top panels of Figure 6 show the average of the derived normalized kernels for the second-epoch ACS/WFC images. Slight asymmetries are evident, primarily in the X-direction. Without our kernel-based corrections, such asymmetric PSF differences would bias astrometry.

Zoom In Zoom Out Reset image size

3.6.4. Second-epoch WFC3 Data

The processes outlined above lead to accurate positions of stars and galaxies in each individual exposure. We align the star positions between exposures using the same iterative approach described in Section 3.4, but this time using the newly derived positions based on the template-fitting method instead of the initial library-PSF-based positions. For each exposure for a given field, this procedure yields a six-parameter linear transformation with respect to the first exposure in the first epoch for that field. For all objects (stars and background galaxies) in all exposures, we apply the known geometric distortion solutions and these linear transformations, to obtain the position in the distortion-corrected frame of the first exposure (i.e., j8f801abq, j92c01b4q, and j92c28ccq for SPHEROID, OUTERDISK, and TIDALSTREAM, respectively). These distortion-corrected frames serve as our reference frames.

In the first-epoch data, we have many exposures. This yields multiple determinations for the position of each star or galaxy in the reference frame. The rms scatter among these determinations provides a measure of the random positional uncertainty in a single measurement for each object. Figure 7 shows the one-dimensional errors in positions for M31 stars (upper panels) and background galaxies (lower panels) as a function of instrumental magnitude (defined as m_instr = −2.5log [electrons]) in the three separate target fields. Brighter objects have a higher signal-to-noise ratio than fainter objects, and therefore have more accurately determined positions. Stars are more compact than background galaxies, and therefore also have more accurately determined positions. At fixed magnitude, the scatter in positional accuracy is larger for the galaxies than for the stars. This is because galaxies have a large variation in sizes and shapes, with the morphology also depending on wavelength. The dependence of astrometric accuracy on size and shape was discussed in Sohn et al. (2010), and the dependence on wavelength was discussed in Mahmud & Anderson (2008). At the bright end in Figure 7, the per-exposure accuracies level off at 0.01–0.02 WFC pixels for both stars and galaxies. When many measurements are averaged, the positional errors decrease as $1/\sqrt{N}$.

Zoom In Zoom Out Reset image size

Using the above procedure,we determine for each object its average first-epoch position in the reference frame. In a similar way, we measure template-based positions for each second-epoch exposure, correct these positions for distortion, and transform them into the reference frame using the template-measured star positions. We then compare these positions with the first-epoch averages to estimate the PM displacement for that object. Division by the time baseline of the observations then yields an estimate of the actual PM.

By construction, our method aligns the star fields between epochs. M31 stars will therefore have zero PM on average. Figures 8 and 9 show the residual motion of each star in X and Y (PM_X and PM_Y) as a function of chip coordinates for one of the second-epoch ACS/WFC images and for one of the second-epoch WFC3/UVIS images. These plots show that the "average star" (given our selection criteria and brightness limits) is measured with a precision of ∼0.05 pixel (but as shown in Figure 7, this is a strong function of brightness). When averaging over an ensemble of N stars in M exposures, the accuracy of the average PM is in principle smaller by a factor of $\sqrt{NM}$. The stellar PMs show that there remain some small systematic trends with position on the detector, at levels of ≲ 0.02 pixel. These trends may be related to small changes in the distortion solution between the first and second epochs, or other low-level systematic effects that are not accounted for in our analysis. Such systematic trends might affect PM results. It is therefore important that we correct for them. We note that we found similar residual trends in all three target fields. This indicates that whatever is causing the low-level systematic trends is related to the detector characteristics rather than the internal kinematics of M31.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The background galaxies were not used to define any transformations between the first-epoch and second-epoch data. Any motion between the galaxies and stars should therefore show up as a displacement between the first- and second-epoch galaxy positions in the reference frame. We measure these displacements in a two step process to remove any systematic PM trends related to the position on the detector, such as those visible in Figures 8 and 9.

First, we calculate the difference between the first- and second-epoch positions for each background galaxy in the reference frame. Then, for each galaxy, we compute the average displacement of stars in the vicinity of the galaxy. We subtract this displacement from the galaxy displacement. This "local correction" ensures that the displacement of each background galaxy is measured only with respect to the M31 stars that fall on the same part of the detector. This removes any remaining systematic PM residual associated with detector position. Each local correction is constructed using stars of similar brightness (±1 mag in the inner 5 × 5 pixels area) and within a 200 pixel region centered on the given background galaxy. The matching of the stars and background galaxies in brightness was motivated by the fact that some known detector effects (such as CTE) depend not only on detector position, but also on the source brightness.

In the top panels of Figures 10 and 11, we show the X and Y proper motions (PM_X and PM_Y) of each galaxy as a function of chip locations (X and Y), for the same second-epoch exposures as in Figures 8 and 9. There are many more stars than background galaxies in each frame, and the positions of the stars are generally more accurately determined than those of the galaxies (see Figure 7). The final PM error for each second-epoch exposure is therefore dominated entirely by the astrometric accuracy for the background galaxies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To determine an average PM of the M31 stars in each field, we start with the measured PM displacements of the background galaxies in the reference frame (including local corrections as described above). For each of the 18 different second-epoch images, we calculate the weighted average using the positional errors determined from the RMS among the first-epoch measurements. Outliers were rejected in this process with an iterative 3σ rejection scheme applied to the 2-d PM distribution. The uncertainty of the average was computed using the bootstrap method (Efron & Tibshirani 1993), with 10,000 bootstrap samples for each case. This provides a rigorous error estimate, which is ultimately based on the actual scatter between results inferred from different background galaxies. This naturally includes any random errors due to photon-counting statistics, as well as any potentially remaining systematic errors that lead to increased scatter (but any systematic errors that affect all sources equally would not be quantified by this).

We then transformed the results to average PMs (and their associate errors) of M31 stars in the directions north and west, in units of mas yr⁻¹. To do this, we used the orientation of the reference image (see Section 3.6.2) with respect to the sky (using the FITS header keyword ORIENTAT), the time baseline, and the fact we are measuring the relative displacement of the background galaxies with respect to the M31 stars in our target fields. In Table 2, we tabulate the PM¹(μ_W, μ_N) inferred from each second-epoch image and the corresponding error, along with the number of background galaxies used for the PM derivation. Figure 12 shows for each of the three target fields the PM estimates inferred from the six-independent second-epoch exposures (2 ACS/WFC + 4 WFC/UVIS).

Zoom In Zoom Out Reset image size

Table 2. Proper-motion Results for Individual Second-epoch Images

Field	Data Set	μW	μN	$\sigma _{\mu _{W}}$	$\sigma _{\mu _{N}}$	Nuseda
		(mas yr−1)	(mas yr−1)	(mas yr−1)	(mas yr−1)
	jb4404vsq	−0.0839	−0.0212	0.0512	0.0327	308
	jb4404vuq	0.0012	−0.0982	0.0297	0.0299	306
SPHEROID	ib4401rsq	−0.0085	0.0184	0.0606	0.0530	176
	ib4401ruq	−0.0564	−0.0379	0.0332	0.0442	176
	ib4401rxq	−0.0772	−0.0104	0.0660	0.0495	186
	ib4401s1q	−0.0895	−0.0126	0.0369	0.0339	172
	jb4405jgq	−0.0706	0.0085	0.0501	0.0504	310
	jb4405jiq	−0.1105	−0.0043	0.0547	0.0501	286
OUTERDISK	ib4402urq	−0.0123	0.0993	0.0679	0.0681	156
	ib4402uwq	−0.0668	−0.1609	0.0649	0.0759	152
	ib4402vqq	−0.0641	0.1353	0.0761	0.0757	161
	ib4402vuq	0.0237	−0.1112	0.0589	0.0553	152
	jb4406muq	0.0425	−0.0895	0.0574	0.0601	321
	jb4406mwq	0.0104	0.0098	0.0564	0.0518	317
TIDALSTREAM	ib4403enq	−0.1116	−0.0572	0.0670	0.0752	185
	ib4403epq	−0.0792	−0.0358	0.0766	0.0812	196
	ib4403esq	−0.0445	−0.0474	0.0883	0.0731	175
	ib4403ewq	0.0298	−0.0148	0.0794	0.0728	185

Note. ^aNumber of background galaxies used for deriving the average proper motion.

Download table as:  ASCIITypeset image

The final estimate of the average PM of the M31 stars in each field is calculated by taking the error-weighted mean of the six independent measurements listed in Table 2. The results are shown in red in Figure 12, and they are tabulated in Table 3.

Table 3. Final Proper-motion Results for the Three Target Fields

Field	μW	μN	$\sigma _{\mu _{W}}$	$\sigma _{\mu _{N}}$
	(mas yr−1)	(mas yr−1)	(mas yr−1)	(mas yr−1)
SPHEROID	−0.0458	−0.0376	0.0165	0.0154
OUTERDISK	−0.0533	−0.0104	0.0246	0.0244
TIDALSTREAM	−0.0179	−0.0357	0.0278	0.0272
Weighted averagea	−0.0422	−0.0309	0.0123	0.0117

Notes. ^aWeighted average of the results for the three-target fields. This is not an unbiased estimate of the PM of the M31 center of mass. It contains contributions also from the internal motions of stars in M31, which are modeled and corrected in Paper II.

Download table as:  ASCIITypeset image

Figure 13 compares the final PM results for the three different M31 fields. The weighted average of the results is shown in black. This weighted average is also listed in Table 3. The figure shows results in physical units of km s⁻¹, instead of mas yr⁻¹. To transform the units, we assumed a distance of 770 kpc (see references in van der Marel & Guhathakurta 2008), so that 0.1 mas yr⁻¹ = 365 km s⁻¹. Distance errors were not propagated in this conversion. The final weighted average PM differs from zero at the ∼4.3σ level. Therefore, we have actually measured a motion, and we have not merely put an upper limit on any motion.

Zoom In Zoom Out Reset image size

To have faith in the results, it is important to assess the internal consistency of the measurements. There are many checks available for this, since we have performed measurements using different exposures, with different instruments, and for different fields.

For a given field i (SPHEROID, OUTERDISK, or TIDALSTREAM), second-epoch instrument j (ACS/WFC or WFC3/UVIS), and coordinate direction k (West or North), we identify the set of l PM measurements μ_ijkl (either four or two measurements, depending on the instrument) with random errors Δμ_ijkl. There are a total of 36 measurements (see Table 2). We combine the different measurements l to obtain the 12 weighted averages ${\overline{\mu }}_{ijk}$ with random errors $\Delta {\overline{\mu }}_{ijk}$.

The quantity

$$\begin{equation} \chi ^2_{1} = \sum _{ijkl} \left({ { \mu _{ijkl} - {\overline{\mu }}_{ijk} } \over { \Delta \mu _{ijkl} } } \right)^2 \end{equation} \tag{ 2 }$$

provides a measure of the extent to which different measurements for the same field and with the same instrument agree to within the random errors. We find that χ²₁ = 26.2. In absence of systematic errors, one expects χ²₁ to follow a χ² probability distribution with N_DF = 36 − 12 = 24 degrees of freedom. The expectation value for such a distribution is N_DF, and the dispersion is ${\sim} \sqrt{2 N_{{\rm DF}}}$. Hence, the measurements for the same field with the same instrument are statistically consistent with each other. This is also visually evident from inspection of Figure 12, which shows, furthermore, that the agreement is least good for the OUTERDISK field.

The quantity

$$\begin{equation} \chi ^2_{2} = \sum _{ik} { { ({\overline{\mu }}_{i1k} - {\overline{\mu }}_{i2k})^2 } \over { \Delta {\overline{\mu }}_{i1k}^2 + \Delta {\overline{\mu }}_{i2k}^2 } } \end{equation} \tag{ 3 }$$

provides a measure of the extent to which measurements for the same field with different instruments agree to within the random errors. We find that χ²₂ = 8.5. In this case, N_DF = 6, so the measurements for the same field with the different instruments are also statistically consistent with each other.

Finally, the quantity

$$\begin{equation} \chi ^2_{3} = \sum _{ik} \left({ { {\overline{\mu }}_{ik} - {\overline{\mu }}_{k} } \over { {\Delta {\overline{\mu }_{ik}} } } } \right)^2 \end{equation} \tag{ 4 }$$

provides a measure of the extent to which measurements for different fields agree to within the random errors. Here, ${\overline{\mu }}_{ik}$ with random errors $\Delta {\overline{\mu }}_{ik}$ are the weighted averages over all exposures for given field (top three lines of Table 3). The ${\overline{\mu }}_{k}$ are the weighted averages over all fields (bottom line of Table 3). We find that χ²₃ = 1.9. In this case, N_DF = 6 − 2 = 4, so the measurements for different fields are also statistically consistent with each other. This is also visually evident from inspection of Figure 13.

The fact that all measurements are statistically consistent indicates that there is no additional scatter in the data that is unaccounted for by the random errors. This justifies the use of weighted averages in combining the different measurements (which propagates only random errors).

The final weighted average PM corresponds to a displacement of 0.0074 ACS/WFC pixel over the ∼7 year time baseline for the SPHEROID field. The per-coordinate error in the displacement is only 0.0024 pixel. The displacement errors are even lower for the other fields, which have shorter time baselines.

These displacements are very small, but they can nonetheless be accurately measured. This can be understood based on first principles. A single source can be centroided to an accuracy of the order of σ/(S/N), where σ is the Gaussian dispersion of the object size, and S/N is the signal-to-noise ratio of the observation. For a compact source with S/N ≳ 100, this yields uncertainties of 0.01–0.02 pixels, consistent with what we found in Figure 7. Therefore, even a single bright compact source in a single exposure can provide an accuracy comparable to the M31 PM displacement. Averaging over multiple sources in many exposures/fields reduces the random uncertainties to thousandths of a pixel and hence allows a solid measurement of the PM displacement.

The key to our robust measurement is, though, control of systematic errors, since those are not guaranteed to decrease by averaging multiple measurements. We have paid careful attention throughout our observation planning and subsequent analysis to identify possible sources of systematic errors and adjusted our methodology to minimize them as follows.

• 1.  
  The observations were obtained with the F606W filter, which has the best PSF and distortion models for both ACS/WFC and WFC3/UVIS (see Section 2.2).
• 2.  
  The ACS/WFC images were obtained using the same telescope orientation and pointing in both epochs (see Section 2.2). This allows any static astrometric residuals that depend on the position on the detector (e.g., static geometric distortion solution errors) to cancel out in the differential PM measurement.
• 3.  
  The ACS/WFC data were explicitly corrected for the effects of imperfect CTE, reducing any potential astrometric impact (see Section 3.2).
• 4.  
  We used different geometric distortion solutions for the different epochs of data, thus minimizing any impact due to potential time variations in the geometric distortion solutions (see Section 3.3).
• 5.  
  We carefully selected our samples of M31 stars and background galaxies to minimize any bias introduced by interlopers (see Section 3.5).
• 6.  
  We built a unique template for every individual source that was measured, thus avoiding any ad hoc assumptions about the properties of PSF or galaxy shapes (see Section 3.6.1).
• 7.  
  We fitted and accounted for PSF variations between different epochs, thereby minimizing any potential PM biases (see Section 3.6.3).
• 8.  
  We measured the relative PM of each background galaxy with respect only to its neighboring stars (the "local correction"), thus minimizing the impact of spurious spatially varying PM residuals (see Section 3.8).
• 9.  
  The above local measurement was made with respect only to stars of similar brightness, thus minimizing the impact of magnitude-dependent PM residuals, including CTE effects (see Section 3.8).
• 10.  
  We obtained the final PM for each field by averaging over background galaxies at many different locations on the detector. This minimizes the potential impact of local peculiarities in detector properties or the astronomical scene (see Section 4.1). Our approach is significantly advantageous over PM studies that use only a single background quasar in an image.
• 11.  
  We determined the random PM error for each exposure from the actual scatter between results from different background galaxies. The resulting errors include all sources of scatter, and not just those from photon-counting statistics (see Section 4.1).
• 12.  
  We used robust statistical measures with outlier rejection throughout our analysis, thus minimizing the potential influence of individual erroneous measurements (e.g., Section 4.1).
• 13.  
  We showed that PM measurements with different instruments, taken at different times and with different telescope orientations, as well as measurements of different fields, all yield statistically consistent results (see Section 4.2). This rules out a large range of possible residual systematic errors, including most possible astrometric residuals specific to a given instrument or detector.

In summary, there is strong reason to believe that the PM measurements and the quoted uncertainties are robust and free of the potential systematics considered in this section.

The results presented thus far do not directly measure the PM of the M31 COM. Instead, in every field, we measure the sum of the internal kinematics of the M31 stars and the COM motion. The rotation curve of M31 has an amplitude of ∼250 km s⁻¹ (Corbelli et al. 2010). The contributions from internal kinematics are therefore not necessarily negligible compared to the uncertainties in our measurements (see Figure 13).

In Paper II, we model the internal kinematics explicitly, and we derive an unbiased estimate for the COM PM. As it turns out, this estimate is not very different from the weighted average presented in Table 3. We also show in Paper II that this estimate is statistically consistent with the estimate for the transverse motion of M31 derived with the independent methods of van der Marel & Guhathakurta (2008), based on the kinematics of M31 and Local Group satellites. This is a further indication that any remaining systematic uncertainties in our measurements are likely to be small.

The observed PMs presented here are heliocentric motions, i.e., not corrected for the reflex solar motion in the Milky Way. This known reflex motion falls in the same quadrant of PM space as our measurements (see van der Marel & Guhathakurta 2008). The actual transverse motion of M31 with respect to the Milky Way is therefore closer to zero than what we show in Figure 13. This provides another useful consistency check on our measurements, since there are theoretical reasons to suspect that the transverse motion of M31 with respect to the Milky Way should be small (e.g., Peebles et al. 2001). In fact, we show in Paper II that the transverse motion implied by our data is consistent with zero, meaning that M31 is likely moving on a nearly direct radial orbit toward the Milky Way.