THIRD-EPOCH MAGELLANIC CLOUD PROPER MOTIONS. I. HUBBLE SPACE TELESCOPE/WFC3 DATA AND ORBIT IMPLICATIONS

Geha et al. (2003) identified a total of 54 QSOs behind the Clouds from their optical variability in the MACHO database. These QSOs provided the inertial reference frame against which the PM measurements were made over a 2 yr baseline in K1 and K2, based on observations obtained between 2002 August and 2005 June. For the most efficient use of HST resources we decided to observe in SNAPSHOT mode for both epochs. Our two SNAP programs of imaging with the HRC achieved an overall completion rate of 48% yielding a final sample of 21 QSOs behind the LMC and 5 behind the SMC.

We obtained an additional epoch of data during the period 2007 July to 2008 November. These observations executed with the HST WFPC2 camera, due to the failure of ACS at that time. We briefly describe these data in Appendix A, but due to their limited quality, we do not use them in the analysis presented here.

We obtained another epoch of data with the HST WFC3 camera during the period 2009 October to 2010 July. These observations also executed in SNAPSHOT mode, and we refer to them as the "third epoch." The final yield from this program was 11 observed QSO fields for the LMC and 4 observed QSO fields for the SMC. However, two QSO fields, one in each galaxy, proved unsuitable for PM determinations. As discussed in Appendix B, in one field the QSO showed a bright extended host galaxy, which complicates astrometry, and in the other field the nature of the QSO was not confirmed. These two fields are omitted from the subsequent discussion and analysis. Therefore, the final number of QSOs used for the three-epoch analysis is 10 for the LMC and 3 for the SMC.

The HRC data, which makes up the first two epochs of this study, is discussed in detail in K1. Briefly, for our main astrometry goals we chose the F606W filter which is a broad V filter with high throughput. Eight dithered exposures were taken in each epoch to minimize pixelization-related systematic errors, and the dither pattern was kept the same over the two epochs. Exposure times were chosen so as to achieve a signal-to-noise ratio (S/N) of at least 100 for the QSOs based on the known MACHO magnitudes at the time.

For the third epoch with WFC3/UVIS (PID 11730), we again used the F606W filter and used the four-point DITHER-BOX pattern which is made for optimal sampling of the point-spread function (PSF). We aimed at S/N ∼ 200 for the QSOs in order to match the astrometric quality of the ACS data. Even though the WFC3/UVIS pixels are bigger than ACS/HRC pixels, with adequate dithering this should not degrade the astrometry because the astrometric error is proportional to the FWHM/(S/N) of the target. With QSO brightnesses in the epoch 1 sample ranging from 16.5 ⩽ V ⩽ 22.0, and an average brightness of V = 19.5, this yielded total science exposure times for epoch 3 ranging from 2.6 to 17.7 minutes, with an average of 12 minutes. Tables 1 and 2 describe the combined three-epoch data set, including R.A./Decl. for the QSOs, whether or not they were observed with WFC3, the total WFC3 exposure times, and the corresponding time baseline for the PM measurement.

Figure 1 shows the locations of the QSOs behind the LMC and SMC. QSOs for which we obtained two epochs of data are marked with yellow diamonds. Their distribution behind the LMC is sparse but reasonably uniform. QSOs observed with the WFC3 are marked with red squares. Note that there is 1 "new" LMC field for which there is only first epoch ACS and third epoch WFC3 data.

Zoom In Zoom Out Reset image size

As in K1 and K2, we used the bias-subtracted, dark-subtracted, flat-fielded images (_flt.fits) provided by the STScI data reduction pipeline for our basic reduction purposes. Anderson & King (2004) found that these images are well-behaved astrometrically. They are not, however, corrected for geometric distortion. We did not use the geometrically corrected products created by the MultiDrizzle software (_drz.fits) since they involve a resampling that might degrade the astrometry. Instead of geometrically correcting the images, we geometrically corrected the positions that we measured on the images. The positions were measured by means of empirical PSFs constructed from other data sets (similar to the approach discussed by Anderson & King 2006 in ACS/ISR 06-01). The geometric corrections were performed according to Bellini et al. (2011). In their study, the geometric distortion is modeled with two components: a third-order polynomial and an empirically-derived look-up table for the residuals, one per chip and filter combination. The method is akin to that used in Anderson & King (2004) and is shown in Bellini et al. (2011) to be stable over time to better than 0.008 pixels (rms) in each coordinate.

Since we are performing relative astrometry we have freedom in constructing the masterframe (per QSO field) into which all relevant star positions are to be transformed. While PSF-fitting was performed on the flt images, we did use stars in the drz image as a starting point to defining our masterframe, because their WCS header information is the most accurate. Simple aperture photometry was performed on the drz image to obtain initial positions and magnitudes for real sources in the region of overlap of all four individual exposures with the orientation north up and east to the left. We then did a simple transformation of these initial positions to an (x, y)-frame of identical orientation, but with a chosen center and 25 mas pixels (for simplicity). A six-parameter linear transformation for the 4 individual flt WFC3 exposures into this frame was then calculated. Based on previous experience, we only use stars brighter than an instrumental magnitude, M_INSTR = −8 (S/N ∼ 40). The instrumental magnitude is defined as −2.5log (counts) in the F606W flt image (without application of a zero point). We then calculated a new masterframe position for each source by averaging the linearly transformed WFC3 positions. At this stage we iterated on the linear transformations, using only stars with measurement errors smaller than ∼0.1 pixels, that are common to all four WFC3 exposures. The positional errors are calculated as the rms scatter between multiple measurements of the same source.

In general, we can routinely center WFC3 stars to 0.02 pixels. Because our exposures are deeper and the field of view (FOV) is much larger than for ACS/HRC, we identify on average ∼200–800 sources in the WFC3 images, depending on the field. By contrast there are far fewer sources in the HRC fields.

In our original K1 analysis we argued, based on detailed calculations and consistency checks, that the degrading charge transfer efficiency (CTE) of the HRC was not expected to systematically affect our analysis. The main reason was that the SNAP nature of the program yielded random roll angles of the telescope for each QSO field. Given that we had N = 21 QSO fields with more or less random detector orientations on the sky, any CTE-induced astrometric shift along columns would average down to zero ∝N^−1/2. There is no explicit model for the underlying HRC CCD detector physics and charge trap properties, as there now is for other cameras (Anderson & Bedin 2010). The residual CTE effects, given the random roll angles of our exposures, would be expected to puff up the final rms of the measurements. We included this effect by using the rms scatter between fields as our final error estimates, instead of simply propagating the per-field weighted errors (which would have suggested a factor of two higher accuracy).

In K1, our final PM value was based on 13 fields that we deemed as high quality for various reasons, and because of this reduced number of fields, and the small number of total SMC fields, CTE issues remained a legitimate concern. P08 subsequently performed an independent reanalysis of our HRC data and used a simple correction for CTE that was a function of Y coordinate of the object, the time since installation, and ∝S/N^−0.42. This latter exponent came from the ACS Handbook. They found, as expected, that correcting for CTE did not have a large systematic effect on the analysis (they obtained a COM PM for the LMC within 1σ from that of K1). However, they did obtain better field-to-field agreement. This yielded smaller error bars on the COM motion, and the possibility to derive a rudimentary PM rotation curve.

For the present study we reanalyzed the HRC data using our own new prescription for CTE degradation as well as other improvements. The main differences in this reanalysis versus the analysis in K1 and K2 are: (1) we use new codes for the implementation of six-parameter linear transformations between the two epochs that employ error-weighting for the stars. In the previous analysis we only used rejection. (2) We employ a prescription to correct for CTE, the form for which is as follows:

$$\begin{eqnarray} y_{\rm NEW} &=& y_{\rm RAW} - \bigg(\alpha \times \frac{y}{1024.0} \times \frac{t}{365.25} \nonumber\\ &&\times \frac{{\rm MAX}[(M_{\rm INSTR}-M_{\rm INSTR,lim}),0]}{{\rm norm}}\bigg), \end{eqnarray} \tag{ 1 }$$

where y is the measured position of a star on the detector, M_INSTR is the instrumental magnitude, typically ranging from M_INSTR = −6 (faintest) to −12 (brightest) in our data, and t is the time since installation of ACS (Modified Julian Date = 52333). CTE affects faint stars more than bright stars, and stars brighter than M_{INSTR, lim} are assumed to be unaffected. Norm is simply a normalizing factor discussed below. The slope, α, was fit in the following fashion: we sequentially corrected averaged detector y positions by a range of slopes. After the application of the correction of a given slope, error-weighted transformations were calculated for stars from epoch 2 into the epoch 1 frame. The input error per epoch is the rms/$\sqrt{N}$ obtained from centering the star over N exposures in that epoch. We then calculated which slope gave the minimum total χ²/dof for these transformations summed over all 21 QSO fields, where dof stands for degrees of freedom, and is defined as follows: N_dof = N_data − N_param, where N_data is twice the number of fields and N_param is the number of free parameters optimized in the fit. The QSO is not included in these transformations.

Figure 2 shows the total χ²/dof for all the QSO fields as a function of α. We tried several M_{INSTR, lim} limits, M_{INSTR, lim} = (− 10, −11, −12), where −12 corresponds roughly to the brightest stars in the field. We normalized the magnitude correction such that it is the greater of zero, when a star is brighter than or equal to −12, or a fractional correction up to a size of approximately 1, for M_INSTR = −6 to −12. Therefore, for M_{INSTR, lim} = (− 10, −11, −12), the corresponding values for norm are (4, 5, 6). This range of M_{INSTR, lim} did not affect the final PM values by more than 1σ. In addition to fitting for CTE-slope, we tried a number of different faint end magnitude cutoffs for the stars used in the transformations, settling on only using stars brighter than M_INSTR = −8. We also employed rejection in the form of (∣residual/error∣ < 3), i.e., a cut based on the size of the residual of the transformation of a star divided by its measurement error.

Zoom In Zoom Out Reset image size

As a consistency check, we calculated the residual in pixels for each Magellanic Cloud star between epoch 2 and epoch 1. This PM should be zero on average, since differential streaming motions within each field are too small to be resolved by our measurements. We averaged the y-residuals over all stars in all QSO fields, and binned the results by M_INSTR. The QSO was again excluded. Figure 3 shows the results. Points in red correspond to α = 0 (i.e., no CTE-correction applied). A clear trend with magnitude in the y-residuals reveals CTE-induced astrometric shifts. By contrast the black points are for our best-fit CTE-correction parameters in Equation (1), which are α = 0.048 for M_{INSTR, lim} = −10 and norm = 4. There is no evidence for any remaining y-CTE trend larger than 0.003 pixels (i.e., smaller than our random errors). There is no indication that serial-transfer CTE in the x-direction might have significantly affected the analysis, so we do not apply any corrections for it.

Zoom In Zoom Out Reset image size

Henceforth we define the PM in the west (μ_W) and north (μ_N) directions in terms of the change in right ascension (α) and declination (δ) on the sky: μ_W ≡ −(dα/dt) cos(δ), μ_N ≡ dδ/dt. With the final ACS to ACS transformation terms in hand, we determined the PM for each source, including the QSO. A linear function with time was fit to the available positions. The best-fit slope gives the PM. The average PM of the stars in each field is zero by construction. However, the QSO appears to move due to its reflex motion with respect to the LMC/SMC stars in the foreground. Hence, the average absolute PM of the stars in the field is simply minus the PM of the QSO. The final error in this PM has two components. One component is the error in the average PM of all the stars: σ_〈PM〉.⁸ This error quantifies how accurately we were able to align the starfields between the epochs. This was added in quadrature to the formal error of the fit to the QSO's motion, δPM_QSO. So the error of a given field in either direction is δμ_W or $\delta \mu _N = \sqrt{\vphantom{A^A}\smash{\hbox{${\delta {\rm PM_{\rm QSO}}^{2}+ \sigma _{{\rm \langle PM \rangle }}^{2}}$}}}$.

Columns 7–10 of Tables 1 and 2 list the PMs thus inferred for all fields with two epochs of ACS data, and the associated errors, (δμ_W, δμ_N). These results are very consistent with those of P08 even though we have used different methodologies for PSF-fitting, linear transformations and CTE-correction. The average residual between our reanalysis and that of P08 is 0.08 mas yr⁻¹ lower in the west direction and 0.06 mas yr⁻¹ lower in the north direction, while the χ² summed over all 21 fields is 12 in the west direction and 8 in the north direction. This is acceptable compared to N_dof. The 1σ difference between our final values and those of P08 arises from the inclusion of error-weighting in the transformations, and not the specifics of the CTE correction. If we exclude error-weighting, we recover values closer to P08. Since there are some sparse fields in which only a small number of stars are used in the transformations, it is especially important to include error-weighting in the analysis.

The inclusion of the WFC3 data gives us 20 images per QSO field for 15 fields: 8 ACS first epoch (e1), 8 ACS second epoch (e2), and 4 WFC3 third epoch (e3). One LMC field has only e1 and e3 data, and hence 12 images. To calculate PMs we first correct the ACS positions for CTE as in Section 2.3. We then perform preliminary transformations for stars (brighter than M_INSTR = −8) in each of the 20 images into the masterframe defined in Section 2.2. All 20 resulting (x, y) centers of each star are averaged and the rms is calculated to make up the final masterframe and error (the QSO positions are not averaged). We then treat this as the initial target frame for transformations. As the input frame, we use the individual positions of the stars in each of the 20 images, with input error the rms per epoch. Again, as in Section 2.3, error-weighting, rejection, and iteration are used in the determination of the final transformation terms into the masterframe.

At this stage we lose the possible advantages of having such a large FOV and high number of sources in e3 because we are limited by the relatively sparse (and small) HRC fields. When we perform final transformations to put e1, e2, and e3 into the masterframe, we only use the sources that are common to the HRC. The number of sources used in the transformations, after all cuts, is listed in Column 6 of Tables 1 and 2. The variation in this number reflects the real variation in stellar density at different locations in the Clouds and is not an artifact of the analysis.

Figure 4 shows the x versus y (equivalent to W, N) positions in pixels of the QSO, after transformation into the masterframe, for four randomly chosen QSO fields. The fields are labeled according to Tables 1 and 2. Linear motion is clearly visible, as is the consistency across the three epochs. Each triangle represents 1 of the 20 positions of the QSO relative to the starfield. The scatter reflects the real error in centering the QSO within a given epoch. In most cases the quality (rms) of the WFC3 data is comparable to that of HRC as we expected based on the fact that we aimed for higher S/N in e3 to account for the larger pixels of the WFC3. This plot represents the reflex motion of the LMC/SMC. The actual PM and its random error were determined as described in Section 2.4 by fitting linear motion with time to all sources. The high quality of the results is visible by eye in Figure 4. For two of the LMC fields in Figure 4 we also show the x and y position in pixels of the QSO as a function of time in Figure 5. The corresponding linear fits are also shown.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Columns 11–14 of Tables 1 and 2 list the PMs thus inferred for all fields with both ACS and WFC3 data, and the associated errors, (δμ_W, δμ_N).

In Figures 6 and 7 (top) we show the measured QSO motions (colored symbols) in the (μ_W, μ_N)-plane compared to the star motions (open circles; centered around zero by construction) for the LMC and SMC, respectively. The stellar motions clearly separate from the QSO motions. The green filled circles show the results from linear fits to all three epochs of data, and the pink squares show results from Section 2.3 for e1 and e2 only. The fact that the averages are consistent with each other is visible by eye. The addition of a third epoch with a different HST instrument not only improves the accuracy of the measured PMs, but also demonstrates that there were no fundamental systematic problems with the earlier two-epoch analysis of the observations in K1, K2, and P08.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As discussed in K1, given the large size the LMC subtends on the sky and the fact that we know it possesses internal rotation, the motion observed for an individual QSO field must be written as: PM(field) = PM(COM) + PM_res(field). Here PM(COM) is the PM of the LMC COM and PM_res(field) is a field-dependent residual. The latter contains contributions from different effects: (1) variations as a function of position in the components of the three-dimensional (3D) velocity vector of the COM that are seen in the plane of the sky ("viewing perspective"), and (2) the internal rotation of the LMC (see formulae in van der Marel et al. 2002, hereafter vdM02). Both corrections make a similar contribution to the observed PMs. Given the location of our fields, the viewing perspective contributes ≲ 0.4 mas yr⁻¹ and the internal LMC rotation contributes ≲ 0.3 mas yr⁻¹. This is less than ∼20% of the overall PM, but that is significant given the accuracy that we are trying to achieve.

The field-dependent residuals depend on the COM PM and vice versa and hence we used an iterative procedure to evaluate them in K1. For the internal rotation model of the LMC we used the model of vdM02, based on the line-of-sight (LOS) velocity data for 1041 carbon stars. This was because our PMs themselves were not accurate enough to significantly constrain the LMC rotation model. P08 took a similar approach by leaving the LMC center and disk plane orientation fixed at the values found by vdM02. However, they did vary and fit the amplitude of the rotation curve.

With the addition of the WFC3 epoch, the per-field, per coordinate PM errors are now typically 0.03 mas yr⁻¹ (see Table 1, Columns 13 and 14), which is only ∼7 km s⁻¹. We show in Paper II that this provides enough information to independently constrain all parameters of the PM rotation field, without reference to prior knowledge from LOS velocity studies. Table 3 lists the results and compares them to the values from vdM02 which were used by K1 and P08. We also compare them to the values from a more recent study by Olsen et al. (2011), who analyzed the LOS velocities of ∼6000 massive red supergiants, oxygen-rich and carbon-rich asymptotic giant branch stars. A detailed discussion of the commonalities and differences between the results is presented in Paper II. We summarize here only the salient features relevant to the determination of the LMC COM motion.

The inferred orientation of the LMC disk, as measured by the inclination i and the position angle θ of the line of nodes, is consistent with the range of values previously reported using other methods (see reviews of van der Marel 2006 and van der Marel et al. 2009). The rotation curve amplitude inferred from the PM data falls between those determined from the LOS velocities of (old) carbon stars (vdM02; Olsen & Massey 2007) and (young) red supergiants (Olsen et al. 2011, see Table 3). This is as expected, since the (bright) stars in our HST fields contain a mix of older and younger stars (see the color-magnitude diagram in Figure 6 of K1). This confirms the accuracy of the new PM data, both in a random and a systematic sense. The rotation curve amplitude inferred from our PM data is not, however, consistent with the result of P08. They were able to infer a rotation curve from the two-epoch data, but their rotation velocity was surprisingly high: 120 ± 15 km s⁻¹. This is approximately 30–40 km s⁻¹ higher than the value derived from the LOS velocities of H I and red supergiants (Kim et al. 1998; Olsen et al. 2011). It would be particularly hard to understand how the stars in the LMC could be rotating significantly faster than the H I gas.

The most important difference between our new results from Paper II and the LOS study of vdM02 (used by K1 and P08) is in the position of the LMC center. In Paper II we find strong evidence that the center of the stellar PM rotation field is consistent with the position of the H I dynamical center. This makes theoretical sense, since the stars and gas are orbiting in the same gravitational potential. However, this result differs by 112 from the position advocated by vdM02, which agrees with the brightest part of the LMC bar. The use of a new rotation center affects the inferred PM of the LMC COM, (μ_{W, LMC}, μ_{N, LMC}), as discussed in Section 4.1.

In Table 1 we list for each field the estimate for the LMC COM that is obtained after subtraction of the contributions from viewing perspective and LMC rotation indicated by the best fit model. In Figure 6 (bottom) we show the estimates PM_est(COM) resulting from this procedure. The scatter is considerably reduced with respect to the top panel. The weighted average of the estimates is μ_{W, LMC} = −1.910 ± 0.008 mas yr⁻¹ and μ_{N, LMC} = 0.229 ± 0.008 mas yr⁻¹. The error bars on these values reflect only the propagation of the random errors in the data. We show in Paper II that there are also contributions from two other sources. First, there is excess scatter between measurements from different fields that is not accounted for by random errors, disk rotation, and viewing perspective. The best-fitting model has a reduced (χ²_min/N_dof)^1/2 = 1.8; here χ²_min is the χ² value of the best-fit model. Second, there are uncertainties in the geometry and rotation of the best-fitting LMC disk model that need to be propagated. The latter provide the dominant uncertainty in the final result. When these sources of uncertainty are taken into account as described in Paper II, the final estimate of the LMC COM PM becomes μ_{W, LMC} = −1.910 ± 0.020 mas yr⁻¹ and μ_{N, LMC} = 0.229 ± 0.047 mas yr⁻¹.

For the SMC we have PM measurements for only five QSO fields, of which only three fields have three-epoch data. This sparse coverage cannot provide much insight into the SMC geometry and rotation field. We therefore fit a relatively simple model to the PM data in order to determine the SMC COM PM (μ_{W, SMC}, μ_{N, SMC}).

To calculate a PM model prediction for each observed field, we keep the SMC center fixed at the H I kinematical center (α, δ) = (1625, −7242) (Stanimirović et al. 2004), the distance modulus fixed at m − M = 18.99 (Cioni et al. 2000), and the radial velocity fixed at v_sys = 145.6 km s⁻¹ (Harris & Zaritsky 2006). We account for the influence of viewing perspective as in vdM02. We also allow for the possibility of a single overall rotation of velocity V_rot in the plane of the sky (i.e., as though we were viewing a face-on disk). We treat the SMC COM PM as a free parameter that is optimized by minimizing the χ² of the model fit to the data.

To determine the uncertainties on the best-fit model parameters, we create pseudo-data in a Monte Carlo sense from the best-fitting model, with properties similar to the real data. For this we use the observational error bars, but multiplied by a factor (χ²_min/N_dof)^1/2 to account for the observed scatter between the measurements. Many different pseudo-data sets are created that are analyzed similarly to the real data set. The dispersions in the inferred model parameters are a measure of the 1σ random errors on the model parameters. In the Monte Carlo simulations we also propagate the uncertainties on the model parameters that are kept fixed. We use Δ(m − M) = 0.10, Δv_sys = 0.6 km s⁻¹, and an uncertainty of 02 per coordinate in the plane of the sky on the SMC center position.

LOS velocity studies provide some insight into the importance of rotation in the SMC. The old stellar population of the SMC, traced by red giants, is consistent with a pressure-supported (V < σ) spheroidal system with little evidence for rotation (V_rot < 17 km s⁻¹; Harris & Zaritsky 2006). This is consistent with the fact that many studies have found a large LOS depth in the SMC (e.g., Crowl et al. 2001). By contrast the H I gas in the SMC shows rotation with an amplitude of ∼40 km s⁻¹ (Stanimirović et al. 2004). This presumably indicates that the H I gas resides in a more disk-like distribution than the old stars. Young stars (spectral types O, B, A) show evidence for a rotation velocity gradient of a similar magnitude as the H I gas, but surprisingly, with a different major axis position angle (Evans & Howarth 2008).

The results of Harris & Zaritsky (2006) suggest that it is reasonable to assume that the SMC has limited rotation in the plane of the sky. If we analyze the PM data with fixed V_rot = 0, the inferred SMC COM PM is μ_{W, LMC} = −0.772 ± 0.033 mas yr⁻¹ and μ_{N, LMC} = −1.117 ± 0.043 mas yr⁻¹. This can be compared to the result of just taking the weighted average of the five SMC PM data points, which yields μ_{W, LMC} = −0.754 ± 0.013 mas yr⁻¹ and μ_{N, LMC} = −1.133 ± 0.016 mas yr⁻¹. These results are statistically consistent, despite the fact that the weighted average does not account for the influence of viewing perspective. In Table 2 we list the estimate for the SMC COM that is obtained after subtraction of the contribution from viewing perspective for each field. Taking the weighted average of results for individual fields yields artificially low error bars, because it does not account for the actual scatter between measurements for different fields, or for uncertainties in the geometry and rotation of the SMC.

To assess the importance of rotation, we have also performed model fits in which the rotation is treated as a free parameter to be determined from the PM measurements. This yields V_rot = 29.5 ± 24.3 km s⁻¹. The PM data are therefore reasonably consistent with the absence of rotation, but they do not provide a very useful constraint. The corresponding best-fit SMC COM PM estimate is μ_{W, LMC} = −0.694 ± 0.074 mas yr⁻¹ and μ_{N, LMC} = −1.055 ± 0.078 mas yr⁻¹. This is consistent with the estimate obtained assuming V_rot = 0 km s⁻¹. However, the uncertainties are now significantly increased, due to the allowed variations in the rotation model corresponding to an amplitude variation ΔV_rot = 24.3 km s⁻¹.

As our final SMC COM PM estimate we adopt the result for V_rot = 0 km s⁻¹, but allowing for a 1σ Gaussian uncertainty ΔV_rot. Based on the preceding discussion, the average star used in our PM analysis probably rotates significantly less than the H I. We therefore adopt somewhat arbitrarily ΔV_rot = 15 km s⁻¹. This yields for the SMC COM PM that μ_{W, LMC} = −0.772 ± 0.063 mas yr⁻¹ and μ_{N, LMC} = −1.117 ± 0.061 mas yr⁻¹. The uncertainties are almost twice those quoted above, obtained with ΔV_rot = 0 km s⁻¹. So as for the LMC, the quality of the data is now such that the uncertainties in the structure of the Cloud play an important role in how well we can establish the COM PM.

The value of (χ²_min/N_dof)^1/2 for the best-fit SMC model is ∼2.8, independent of whether or not V_rot is treated as a free parameter. By contrast, for the best-fit LMC model derived in Paper II, (χ²_min/N_dof)^1/2 = 1.8. Therefore, for both galaxies the measurements show somewhat more scatter than what can be accounted for by the models. This is not surprising, given the complexity of the data analysis, and the relatively idealized nature of the models. Additional scatter might be due, e.g., to subtle detector geometric distortion variations over time, residual astrometric effects due to imperfect CTE, or more complicated galaxy structures than are accounted for in the models. It is intriguing that (χ²_min/N_dof)^1/2 is larger for the SMC than for the LMC. This might indicate that the structure and kinematics of the SMC are more complex (or less in dynamical equilibrium) than that of the LMC.

PM measurements from various sources, including the present paper, are summarized in Table 4. For the LMC, K1 obtained from two-epochs of data that μ_{W, LMC} = −2.03 ± 0.08 mas yr⁻¹ and μ_{N, LMC} = 0.44 ± 0.05 mas yr⁻¹. P08 obtained from the same data that μ_{W, LMC} = −1.956 ± 0.036 mas yr⁻¹ and μ_{N, LMC} = 0.435 ± 0.036 mas yr⁻¹. If we analyze the new PM data with the same fixed LMC center, inclination angle, and line-of-nodes position angle from vdM02 used by K1 and P08, we obtain μ_{W, LMC} = −1.899 ± 0.017 mas yr⁻¹ and μ_{N, LMC} = 0.416 ± 0.017 mas yr⁻¹. This shows that the addition of our new third epoch of data, while keeping the LMC model the same, has the following effects: (1) the random errors in the final LMC COM PM components decrease, as expected based on the increased time baseline for about half the fields; (2) the value of μ_{W, LMC} decreases by about 1.6σ compared to the earlier estimates; and (3) the value of μ_{N, LMC} remains unchanged compared to the earlier estimates, within the uncertainties.

Despite this good consistency, our final LMC result of μ_{W, LMC} = −1.910 ± 0.020 mas yr⁻¹ and μ_{N, LMC} = 0.229 ± 0.047 mas yr⁻¹ deviates significantly from the final K1 and P08 estimates. This is because of the different LMC model used here, and in particular, the different position for the center. The difference is primarily evident as a decrease by ∼0.20 mas yr⁻¹ in μ_{N, LMC}. By contrast, our new value for μ_{W, LMC} agrees with that of P08 at the 1.1σ level. It is also worth noting that our final best-fit estimate does not have uncertainties that are much lower than those of P08, in particular in μ_{N, LMC}. However, the reason for this is that P08 and K1 underestimated the uncertainties by not propagating all relevant uncertainties in the LMC disk model. Therefore, our new values are in fact significantly more accurate (both in a random and systematic sense) than the older ones.

For the SMC, our results are more similar to those of P08 than K2. This is due primarily to the way in which we combined the results for different fields in K2. On a field-by-field basis, K2 and P08 actually obtained rather similar results. In K2 we did not explicitly correct for the astrometric effects of CTE. So we reported a value that was essentially an average of only two fields, because three of the four used QSO fields had been taken with the same roll angle at both epochs. We also increased our errors to reflect a "true" systematic plus random error, estimated from our analysis of the more numerous LMC fields. The K2 result therefore had large error bars. While our new three-epoch result agrees with K2 for μ_{N, SMC}, for μ_{W, SMC} it differs by 2.4σ. Again, this is mainly because of the way in which fields were combined in K2. The new three-epoch result, however, includes CTE corrections and a longer time baseline, and is therefore the more accurate one.

Our three-epoch result for the SMC agrees reasonably well with that of P08. The value for μ_{W, SMC} agree within the uncertainties, while for μ_{N, SMC} the results differ by 1.7σ. This agreement is similar as for the LMC, if one compares analyses with the same assumed LMC center. So overall, there is acceptable internal consistency between studies with different analysis methods, scientific instruments, and time baselines.

A number of older ground-based (and Hipparcos) analyses of the LMC and SMC PMs were discussed in K1 and K2. These studies had large uncertainties, but within these uncertainties the results were generally consistent with what we know now from HST data.

There are now more recent results from two groups that are worth discussing (see Table 4). Vieira et al. (2010) measured the PMs of the Clouds based on CCD and plate material from the Yale/San Juan Southern Proper Motion program, spanning a baseline of 40 years, and covering a large area in the inter-Cloud region as well. They ultimately tie their reference frame to the International Celestial Reference System. Within the uncertainties of ∼0.3 mas yr⁻¹ per coordinate, their LMC and SMC PM determinations agree with those presented here. They obtain a stronger constraint on the relative motion between the Clouds, than on the motion of either cloud individually: (μ_W, μ_N)_SMC-LMC = (0.91, −1.49) ± (0.16, 0.15) mas yr⁻¹. The corresponding difference between the Clouds' motions from our new HST results is (μ_W, μ_N)_SMC-LMC = (1.188, −1.383) ± (0.039, 0.064) mas yr⁻¹. These results agree for μ_N, but differ by 1.7σ for μ_W. Given the totally different methodologies, this is a satisfactory result. However, the uncertainties on the ground-based data are large enough that they do not really help to validate the HST data at the level of its own uncertainties.

Costa et al. performed two recent studies using the 2.5 m telescope at Las Campanas and centered on background QSOs. In Costa et al. (2009) they measure one fairly outlying QSO in the LMC over a five-year baseline. The inferred LMC COM PM depends significantly on the assumed LMC rotation velocity. In Table 4 we quote the average of their results for V_rot = 50 km s⁻¹ and 120 km s⁻¹ (which corresponds to a rotation velocity consistent with that derived in Olsen et al. 2011 and Paper II). The uncertainties are 0.13–0.15 mas yr⁻¹ per coordinate. Their PM value differs from our new HST result by 1.5σ in μ_{W, LMC} and 1.8σ in μ_{N, LMC}. This is more than expected based on random errors alone. However, we are unable to evaluate our HST results on the basis of this study, given that we have many more quasars and typically higher accuracy.

Costa et al. (2011) present an SMC COM PM estimate from five QSO fields over a seven-year baseline. This is similar to our study, but HST of course has better spatial resolution. Their PM value differs from our new HST result by 1.5σ in μ_{W, LMC}, but agrees in μ_{N, LMC} to within the uncertainties. Given the totally different methodologies compared to our work, this too is a satisfactory result. However, again the uncertainties on the ground-based data are large enough that they do not really help to validate the HST data at the level of its own uncertainties.

We explore a total of five models for the MW: three static models with total virial masses of 10¹², 1.5 × 10¹², and 2 × 10¹² M_☉; and two cases in which we allow for the mass evolution of the MW as is expected from the hierarchical buildup of mass in ΛCDM. While it is beyond the scope of this work to provide an exhaustive compilation of the observational determinations of MW mass, we do want to motivate our adopted range by pointing to the work of Gnedin et al. (2010), who obtain an MW virial mass of 1.6 × 10¹² M_☉ if the middle value for their data-based circular velocity determination is used. The uncertainty around this value is 20%.

For the static MW cases we follow B07 and model the MW as an axisymmetric three-component potential with a Navarro–Frenk–White (NFW) halo, Miyamoto–Nagai disk (Miyamoto & Nagai 1975) and a Hernquist bulge (Hernquist 1990). In all cases the bulge is modeled with a scale radius of 0.7 kpc and total mass of 10¹⁰ M_☉. The NFW halo is adiabatically contracted to account for the presence of the disk using the CONTRA code (Gnedin et al. 2004). The NFW density profile is also truncated at the virial radius (unlike in B07). While it is widely agreed that halos should be truncated, exactly how and where the truncation is done is arbitrary. The exact nature of this truncation does not affect our conclusions because we consider a one-parameter family of models in which MW halo mass is considered a free parameter.

Halo and disk parameters are listed in Table 6. In all cases the disk scale radius is kept fixed at 3.5 kpc while the total disk mass is allowed to vary in order to reproduce the correct circular velocity at the solar circle. The MW circular velocity is taken as the updated value of 239 km s⁻¹ at the solar radius of 8.29 kpc (McMillan 2011) rather than the standard IAU value in all models, unlike B07. Rotation curves for each MW model are presented in Figure 8. We follow the methodology outlined in B07 to capture the effects of dynamical friction acting on the Clouds owing to their passage through the MW's dark matter halo.

Zoom In Zoom Out Reset image size

In the cases where the LMC is not on a first passage but has had a past pericentric passage, the typical periods are quite large (≳ 5 Gyr). As such, we also want to investigate a simple model that accounts for the MW's mass evolution over the past 10 Gyr. In this case we model the MW only as an NFW halo (neglecting the disk, etc.). We model the mass evolution following the mean mass growth rate with redshift from Fakhouri et al. (2010). The expected 1σ scatter in this evolution at each redshift is roughly 20% (from Figure 1 of Boylan-Kolchin et al. 2010). We model the concentration evolution following Klypin et al. (2011), and the scaling of the virial radius with redshift and mass following relations in Maller & Bullock (2004).⁹ We consider current-day MW virial masses of 1.5 × 10¹² M_☉ and 2 × 10¹² M_☉ in this exercise. We only explore the higher MW mass range because we are interested in constraining whether the LMC is on a first infall (and as we will show below, it is always on a first infall in the low MW mass case).

Our expectation is that the LMC mass is the dominant uncertainty in its orbital history (over MW mass evolution), since dynamical friction will change the LMC's orbit on timescales shorter than the MW's mass evolution.¹⁰ Therefore, for each static MW mass model, a variety of LMC masses are explored, ranging from 3 × 10¹⁰ to 2.5 × 10¹¹ M_☉ (the motivation for which is discussed directly below). This approach is different from B07, where we only considered one LMC mass of 3 × 10¹⁰ M_☉.

Our knowledge of the mass of the LMC is limited by the fact that kinematic data are only available in the inner ∼9 kpc. The observational estimate of LMC total mass within this distance is 1.3 × 10¹⁰ M_☉ (van der Marel et al. 2009). For each mass case considered here, the LMC is modeled as a Plummer sphere, where the softening radius is chosen such that the total mass contained within 9 kpc is ∼1.3 × 10¹⁰ M_☉ as observed, and as listed in Table 6.

Saha et al. (2010) have detected stars out to 10 disk scale lengths in the LMC (at a scale length of 1.4 kpc, this amounts to stars as far out as 15 kpc). For a bound stellar component to exist at a minimum distance of 15 kpc from the LMC COM, the tidal radius must be at least this large. In a 2 × 10¹² M_☉ MW, the LMC's Roche radius at an MW–LMC separation distance of 49.5 kpc reaches 15.4 kpc if the LMC total mass is 3 × 10¹⁰ M_☉. This is the reason for our choice of lower bound in LMC mass. The highest total LMC mass explored is chosen to match the observed present-day LMC stellar mass and the relations for the expected infall mass from Guo et al. (2010). For the evolving MW mass case we only consider LMC mass = 5 × 10¹⁰ M_☉. For this LMC mass we expect to find solutions where the LMC has completed an orbit about the MW, and this thus serves as a fiducial case for testing the plausibility of a non-first infall scenario.

Whether the SMC is a binary companion to the LMC is primarily determined by the LMC mass, which provides further motivation for varying the LMC mass. If the SMC is a long-term companion of the LMC, it will have been tidally truncated fairly early on, maintaining a roughly constant mass within its tidal radius at late times. Therefore, we do not vary the SMC mass. The SMC is modeled as a Plummer sphere of total mass 3 × 10⁹ M_☉, with a Plummer softening of 1 kpc so that the mass within 3.5 kpc is 2.4 × 10⁹ M_☉ (Stanimirović et al. 2004). We do not account for dynamical friction acting on the SMC as it travels through the dark matter halo of the LMC. This is not significantly different than the approach used in B07.

In our study we use the new third epoch values and updated solar motion for both the LMC and SMC, and their corresponding errors (see lines 1 and 7 of Table 5).

As in B07 we follow the orbit of the LMC backward in time, using the present day velocities and positions and integrating the corresponding equations of motion.¹¹ For each combination of MW and LMC mass, we explore 10,000 Monte Carlo drawings from the LMC's velocity error distribution. We calculate the evolution of the LMC's Galactocentric radius, and keep track of how many pericentric passages it makes about the MW over a Hubble time, assuming the MW mass is static over this time period. In Figure 9 we show the mean number of pericentric passages that the LMC has made about the MW, as a function of LMC and MW mass. Note that in all cases the LMC has made at least one pericentric passage (N = 1): our data indicate that the LMC is currently just past pericenter, in agreement with previous works. This does not, however, imply that the LMC has made a complete orbit about the MW. The combination of LMC and MW masses that yield, on average, solutions where the LMC has completed at least one orbit are indicated by N ⩾ 2; the first infall scenarios are indicated by N = 1. N = 1.5 represent cases where there is roughly equal probability that the LMC is on a first infall or that it makes one orbit about the MW (50% of the solutions have N = 1 and 50% have N = 2). The size of the square indicates the rms error on the mean value, as listed in the legend.

Zoom In Zoom Out Reset image size

For the majority of the cases we consider here, the LMC makes no additional pericentric passages about the MW than the one that it is currently at (it is on a first infall: dark blue squares). As might be expected, the LMC makes more pericentric passages if the MW mass is higher (larger binding energy) and the LMC mass is lower (minimal dynamical friction). The dispersion in the mean values increases as the mass of the MW increases; there are more solutions with one or more complete orbits allowed within the error space within a Hubble time for higher mass MW models. The dispersion decreases as the LMC's mass increases; for high-mass LMC models (consistent with ΛCDM expectations, Guo et al. 2010) dynamical friction limits the number of closed orbit solutions within the error space. This forces the LMC on increasingly eccentric orbits with orbital periods longer than a Hubble time.

Note that we have ignored the mass evolution of the LMC owing to the MW's tidal field in this analysis. Because we have constrained all LMC models to match the total observed mass of the LMC within 9 kpc, higher mass LMC models will have more material at large radii and should thus be more affected by the omission of this physical effect than the lower mass LMC models. However, Figure 9 illustrates that if the LMC today is well-described by the high-mass models, it could not have made an earlier pericentric approach about the MW. As such, our argument is self-consistent.

Figure 9 does illustrate that solutions can be found where the LMC makes more than one complete orbit about the MW within a Hubble time. For these cases, Figure 10 illustrates the mean orbital period for the most recent orbit. Note that if the LMC has completed more than one orbit, the orbital period of the previous passage is expected to be larger than the values quoted for the most recent passage because of dynamical friction. For the 1 × 10¹² M_☉ MW, the periods are >9 Gyr. For the 1.5 × 10¹² M_☉ MW, the typical periods are >5 Gyr. In the most tightly bound scenario (LMC mass = 3 × 10¹⁰ M_☉, static MW mass ∼2 × 10¹² M_☉), the LMC completes on average 2 orbits about the MW where the most recent orbital period is 3–4 Gyr. This conclusion is no different than in B07 where we showed that the high MW model recaptures the isothermal sphere orbit from previous studies (albeit with larger orbital periods).

Zoom In Zoom Out Reset image size

For the cases where the LMC makes at least one complete orbit about the MW, i.e., LMC mass <5 × 10¹⁰ M_☉ and MW mass >1.5 × 10¹² M_☉, we look at the effect of a cosmologically motivated mass growth history for the MW on the orbital history of the LMC. In Figure 11 we show the LMC Galactocentric radius as a function of time in the past (t = 0 corresponds to today) for a current day MW mass of 1.5 × 10¹² M_☉(left) and 2 × 10¹² M_☉(right), and for LMC mass = 5 × 10¹⁰ M_☉ in both cases. In the first case, the period of the orbit is >8 Gyr for an evolving MW, which is ∼2 Gyr longer than for a static MW. The LMC goes well outside the virial radius of the MW, i.e., the orbit is dramatically more eccentric than in the static MW case. In the second case, it is shown that a 3 Gyr orbit would still on average increase by 1 Gyr if the mass evolution of the MW is taken into account. The period of the last complete orbit of the LMC about the MW is critical in figuring out the relative importance of MW tides in forming the Stream (discussed further in the next section).

Zoom In Zoom Out Reset image size

We are also interested in reevaluating the orbit of the Clouds about each other. A chance three-way encounter between the MW, LMC, and SMC at z = 0 has low probability, and therefore we seek to enumerate under what conditions, given the new velocities, the Clouds constitute a binary for a significant fraction of a Hubble time. Also, we have recently put forth a new Magellanic Stream model that relies on the past interactions between the Clouds when they are far from the MW virial radius, rather than on the influence of the MW on the Clouds-system (Besla et al. 2010; Diaz & Bekki 2011). The orbital history of the SMC about the LMC is critical to assess the viability of such a model.

We draw 10,000 cases from the LMC and SMC PM error distributions in Monte Carlo fashion, and orbits are computed for each combination of LMC and MW masses. The SMC is assumed to be tidally truncated and hence its mass is kept fixed at 3 × 10⁹ M_☉. At each time step, we compute the escape velocity of the SMC from the LMC, V_esc, and compare this to the relative velocity between the Clouds, V_LS. We do not require them to be bound today, but rather search for cases at some point in time where V_LS < V_esc, and keep track of the longest amount of time over which this condition remains true.

This analysis, therefore, accounts for capture events and the fact that the MW's tidal field might be disrupting the binary today. By also requiring that the escape velocity condition be satisfied, as opposed to searching for minima in the distance between the Clouds as many previous authors have done, we are choosing real binary configurations, rather than chance superpositions caused by cases in which the SMC is on its own elliptical orbit about the MW, and therefore formally gets closer and farther away from the LMC even though the LMC may be on a very different orbit (as illustrated in Rŭžička et al. 2010).

In all searched cases the SMC has made at least one close encounter with the LMC in the past. This is unsurprising, as it is a direct result of the relative orientation of their 3D velocity vectors. The SMC is currently ∼20 kpc from the LMC and moving away from it. Integrating the SMC's orbit backward in time will thus always yield a close encounter with the LMC within the past 500 Myr. This conclusion has been drawn by many previous authors (Rŭžička et al. 2010, GN96) and is the basis of the LMC–SMC collision theory put forth in Besla et al. (2012).

In Figure 12 we show the outcome of our search for binary orbits between the Clouds. The percentage of Monte Carlo drawings for which V_LS < V_esc for more than 2 Gyr in the past is shown as a function of LMC and MW mass. We have chosen the 2 Gyr duration since this is the estimated age of the Magellanic Stream (as discussed below), and so the LMC and SMC should be a binary for at least this amount of time in each LMC and MW mass combination. Colored squares indicate the corresponding percentages, as indicated in the legend.

Zoom In Zoom Out Reset image size

Unsurprisingly, the percentage of cases for which V_LS < V_esc increases as the mass of the LMC increases (increasing binding energy). However, we also note a dependence on MW mass. At an LMC mass of 3 × 10¹⁰ M_☉, the percentage of orbits that satisfies this criterion is 20% if the MW mass is 1 × 10¹² M_☉. This number drops to 4% if the MW mass is 2 × 10¹² M_☉, implying that bound configurations between the LMC and SMC are statistically unlikely for low LMC mass and high MW mass. As the LMC mass increases, these statistics improve, but for a MW mass of 2 × 10¹² M_☉, they only reach 20% once we get to an LMC mass of 8 × 10¹⁰ M_☉. The corresponding percentage for a 1 × 10¹² M_☉ MW is 50%.

The rms of the distribution of cases matching this escape velocity criteria increases as the LMC mass increases, representing a larger spread in the possible durations. As discussed above, as the MW mass increases there are much fewer binary LMC–SMC orbits, even for high LMC mass. The PMs indicate that the current relative velocity between the Clouds is high. This implies that the SMC must be on an eccentric orbit about the LMC. For larger MW mass models, the tidal field is more effective at disrupting such a wide binary pair.

Therefore, from our statistical analysis we conclude that large LMC masses (≳ 1 × 10¹¹ M_☉) are favored if the Clouds are to have been in a binary, and that further, an MW mass ≲ 1.5 × 10¹² M_☉ is also required. From Figure 9, orbits in this mass range almost always correspond to first infall scenarios. The assumption that the Magellanic Clouds constitute a long-lived binary pair thus implies that the Clouds are likely on their first infall about the MW.

Based on our searched parameter space, and the requirement that the LMC and SMC have been a long-lived binary, we adopt a canonical model with the highest MW mass and lowest LMC mass that will fulfill these criteria. This gives (LMC, MW) mass = 1.8 × 10¹¹ M_☉, 1.5 × 10¹² M_☉.

In Figure 13 we show an example of one of these long-lived binary states (N = 4–6) in the canonical model (note again that the SMC is assumed to be tidally truncated and therefore kept at fixed mass in our models). The left-hand panel shows the Galactocentric radius of the LMC, the Galactocentric radius of the SMC and the relative distance between the two Clouds as a function of time. The middle panel shows the past orbits of the Clouds in the Galactocentric (Y–Z)-plane. The right panel shows the past orbit as a function of (l, b). Interestingly, the new velocities now cause the past SMC orbit to cross over in the direction of the Stream, something that we did not achieve with the K2 velocities, but was advocated for in earlier works like GN96.

Zoom In Zoom Out Reset image size

We note that with the vdM02 center, the relative LMC–SMC velocity is 143 ± 31 km s⁻¹ (Table 5). With the new center adopted in this work, it is 128 ± 32 km s⁻¹. Since the relative velocity is higher in the former case, the orbit of the SMC would be less bound at fixed mass, and therefore we expected fewer bound cases. Surprisingly, we in fact find no binary LMC–SMC orbits when using the vdM02 velocities. We attribute this not just to the larger absolute value of relative velocity, but to the fact that the angle between the velocity vectors is larger. This is an interesting finding that argues in favor of the new center used here and derived in Paper II.

While dynamical friction owing to the passage of the Clouds through the MW's dark matter halo was accounted for, we did not explicitly account for dynamical friction acting on the SMC owing to its passage through the LMC's dark matter halo. The inclusion of this effect would place the SMC on an increasingly eccentric orbit in the past. We have recently advocated for a high-eccentricity SMC orbit about the LMC in our simulations of the Magellanic Stream; such eccentric orbits are needed to explain the large angular extent of the Stream and to prevent the Clouds from merging (Besla et al. 2010). However, such eccentric orbits are more easily disrupted by the MW's tidal field and thus work against maintaining long-lived binary LMC–SMC states. This further strengthens our argument that the preferred mass scale for the MW is less than 2 × 10¹² M_☉, since in reality there should be even fewer viable binary orbits in this model than we currently find.

Another way to approach these questions is to look in more detail at the LMC–SMC dynamics in LMC–MW mass combinations that give rise to more traditional orbital trajectories for the Clouds about the MW. Specifically, consider the case of a low-mass (tidally truncated) 3 × 10¹⁰ M_☉ LMC orbiting around a 2.0 × 10¹² M_☉ MW, with an orbital period of ∼4 Gyr. We draw 10,000 LMC and SMC PMs in Monte Carlo fashion from their error distributions, and plot the relative PM between the two Clouds in Figure 14. Each value is color-coded by the length of time for which the criterion V_LS < V_esc held true in the past. While it is clear that the large majority of pairs within the 1σ error ellipse do not produce a bound system for >2 Gyr, it is interesting that some orbits with a bound pair in excess of several Gyr do exist within this error space. A larger abundance of such orbits can be found at relative PMs that are more in the west direction than our measurements imply. So the combination of galaxy masses shown in this plot, while unlikely, cannot be ruled out strictly based on the requirement that the LMC and SMC must have been bound for at least 2 Gyr. However, it should be noted that this does not imply that traditional models for the Magellanic Stream (which have often used masses such as shown in this figure) are tenable, because such models produce orbital periods that are inconsistent with the age of the Stream (discussed more below).

Zoom In Zoom Out Reset image size

In summary, perhaps unsurprisingly, we find that to obtain long-lived binary LMC–SMC orbits the LMC mass must be relatively high, and the MW mass must be relatively low. The lower mass of the MW implies that the Clouds are on a first infall and, as such, the MW's tidal field would be insufficient to disrupt the binary pair. We are more readily able to find such long-lived orbits with the new LMC center derived here rather than that previously derived by vdM02. We also find that a close encounter between the Clouds in the recent past is a generic result.

Given that the Clouds are a binary pair, the more massive LMC sets the orbital period for the Clouds about the MW. As shown above, the LMC mass must be ≳ 1 × 10¹¹ M_☉ for the Clouds to have been a long-term binary system. The smallest orbital periods about the MW that we can obtain for such large LMC masses are in excess of 4 Gyr. This is problematic for models that rely on Milky Way tides to form the Stream, because age estimates of the Stream imply that the Stream is a young feature, much younger than 4 Gyr old. In fact, given the rate at which the Stream is being ablated away, as inferred from simulations of the observed anomalously high Hα emission in the Stream (Weiner & Williams 1996) by Bland-Hawthorn et al. (2007), and based on estimates of the survivability of high-velocity clouds in MW-type environments (Heitsch & Putman 2009; Kereš & Hernquist 2009), it is unlikely that the Stream could have survived about the MW for more than ∼1–2 Gyr.

The past orbit of the LMC and SMC, implied by the new velocities, is now better aligned with the location of the Stream HI than it was with the K1/K2 velocities. As discussed in B07, this is directly due to the lower value for μ_N, which in turn is obtained here directly from the dynamical center of the LMC having changed with respect to what was used in K1. However, we note that the LMC orbit still does not trace the Stream closely: there remains a sizeable offset. The new SMC PMs, however, do allow the past SMC orbit to cross over the Stream location, which is now easier to reconcile with models where the Stream is largely stripped from the SMC (GN96; Rŭžička et al. 2010; Diaz & Bekki 2011; Besla et al. 2012) than it was with the K2 SMC PMs.

As we learned from B07, the question of first infall is a model-dependent rather than simply a velocity-dependent one. We showed in that work that no matter whose velocities were used (K1/GN96/prior ground-based determinations), the Clouds were always on a first infall in a 1.2 × 10¹² M_☉ NFW MW halo, and that in order to recapture something resembling the "traditional" isothermal sphere orbit, we had to go to a 2 × 10¹² M_☉ NFW halo.

Given that the new Galactocentric velocities derived here are lower than those in K1, we want to critically reexamine the arguments for whether the LMC may be on a first passage, which we define as an orbital solution wherein the LMC has first entered the virial radius of the MW within the past 1–4 Gyr and has not completed an orbit in that time. We have already shown above that with our new velocities, a low LMC mass can make past pericentric passages about a high-mass MW. Therefore, arguments for a first passage, being still model-dependent in nature, can be cast as arguments against a low-mass LMC and a high-mass MW. Our arguments fall into three main categories: (1) orbital eccentricity and cosmological expectations, (2) LMC tidal radius, and (3) LMC–SMC binarity. We also discuss a fourth argument that attempts to reconcile the high gas content and active star formation of the Clouds with their current close proximity to the MW.

Orbital eccentricity and cosmological expectations: As shown in Figure 10 in all cases where the LMC does complete an orbit about the MW, the corresponding orbital periods are typically large (>4 Gyr). Such orbits take the LMC to apocentric distances of the order of the virial radius of the MW and thus imply large orbital eccentricities. The computed eccentricities for the 3 × 10¹⁰ M_☉ LMC and the 2 × 10¹⁰ M_☉ evolving MW model are 0.6–0.7. These eccentricities get higher for larger LMC masses. Comparing to Figure 5 of Boylan-Kolchin et al. (2011), only 20% of LMC analogs accreted at early times (>8 Gyr) have such orbits. Indeed, in the prevailing ΛCDM model of hierarchical structure formation, massive satellites that are accreted at early times are very rarely found on highly eccentric orbits about MW-type hosts at z = 0 (Boylan-Kolchin et al. 2011; Wetzel 2010; Stewart et al. 2008).

Generally, subhalos that are accreted on ∼radial orbits at early times are either preferentially destroyed or exist on more circularized orbits today (e.g., Benson 2005). While 50% of LMC analogs accreted more than 4 Gyr ago and 60% of LMC analogs accreted more than 8 Gyr ago have eccentricities <0.5 in the Millennium II sample, only 20% of LMC analogs accreted within the last 2–4 Gyr have such modest eccentricities (Boylan-Kolchin et al. 2011). A high MW mass, low LMC mass today is not ruled out from cosmological expectations, but the high eccentricities (0.6–0.7) calculated from our orbits are more typical of LMC analogs accreted within the past 4 Gyr. Therefore, from the eccentricities of our orbits, a first infall scenario for the Clouds is the favored orbital solution of ΛCDM theory.

We can also make a more general timing argument based on the fact that the Clouds and the MW must have been in close proximity at the time of the big bang (in analogy with similar arguments for, e.g., M31 and Leo I Li & White 2008). This implies that the Clouds had a pericentric approach with the MW at the time of the big bang. So if there has been more than one complete orbit since, then the period must be <T₀/2, where T₀ = 13.73 Gyr. Hence, orbits with a pericenter >6.9 Gyr ago are not physical, because the Clouds and MW were not together at the big bang. Of course, this is oversimplified for many reasons, e.g., the mass of the MW increases with time. But as we show in Figure 11, this typically increases the period by ∼2 Gyr. Therefore, finding an orbit that is consistent with the PM data with a pericenter at, e.g., 10 Gyr ago, does not rule out a first infall scenario, since such orbits are not physical based on such a timing argument.

LMC tidal radius: As discussed in Section 6.1 the LMC total mass must be at least 3 × 10¹⁰ M_☉ from the Saha et al. (2010) observations. Muñoz et al. (2006) have claimed a detection of LMC stars out to 20 kpc. If stars do indeed exist out to 20 kpc, the required LMC mass is ∼8 × 10¹⁰ M_☉. It should be pointed out that more traditional Magellanic Stream models (GN96; Diaz & Bekki 2012) assume a total LMC mass of 1 × 10¹⁰ M_☉ which is already ruled out by the observations. Their choice of LMC mass is important to point out since this is the reason their orbital periods are so different from ours.

LMC–SMC binarity: Studies of the star formation histories of the Clouds all agree on the fact that the star formation rate started to increase ∼4–6 Gyr ago, and was quiescent in both galaxies before this time (Harris & Zaritsky 2009; Cignoni et al. 2012). The coincident increase in star formation in both galaxies may imply that they were in a common envelope at the time, and corresponds to the epoch that the LMC captured the SMC in the Besla et al. (2012) model.

If we therefore take the stance that the two Clouds must have fallen into the MW as a bound system (although they need not necessarily be bound at present) we find that we are able to place limits on both the MW and LMC galaxy masses. We find that in general an LMC mass >1 × 10¹¹ M_☉ is needed for the SMC to have been bound to the LMC for >2 Gyr in the past, and that further an MW mass <1.5 × 10¹² M_☉ is required. Note that LMC orbits at this particular mass are almost always on first infall regardless of MW mass. The LMC mass required to keep the SMC bound to it goes up as the mass of the MW increases, in order to compensate for the increasing tidal field of the MW. But as the LMC mass increases the orbital eccentricity does as well, making a first infall scenario more likely based on the plausibility argument above that high eccentricities and early infall times are relatively rare in cosmological simulations.

It is harder to maintain an LMC–SMC binary for long periods of time in high-mass MW and low-mass LMC models: less than 5% of the 10,000 searched Monte Carlo cases result in binary configurations, meaning such models generally require that the LMC capture the SMC while in orbit about the MW—a statistically improbable event. This argument implies that the MW's mass is likely less than 2 × 10¹² M_☉.

The clouds and the morphology–distance relation: van den Bergh (2006) argued for a recent accretion of the Clouds, based on a comparison of the morphologies of MW and M31 satellites that showed that the Clouds are the only gas-rich dIrr galaxies at small Galactocentric distances, all other dIrrs being at large distances from their respective hosts. Some recent studies have looked at this issue from a larger statistical/cosmological point of view.

Tollerud et al. (2011) study a volume-limited spectroscopic sample of isolated galaxies in the Sloan Digital Sky Survey (SDSS) and find that bright satellite galaxies around MW-type hosts are significantly redder than typical galaxies in a similar luminosity range, and argue that this is indicative of environmental quenching. This is found to be in stark contrast to the LMC, which is anomalously blue in comparison to other LMC–MW analogs in SDSS. The authors attribute this to the fact that the LMC may be undergoing a triggered star formation event upon first infall.

Geha et al. (2012) have used the NASA-Sloan Atlas to demonstrate that dwarf galaxies in the field (with masses in the LMC range) all have active star formation (<0.06% of field-dwarfs in their sample have no star formation). By contrast, the majority of quenched galaxies are all within 2–4 virial radii of a massive host. Therefore, ending star formation in such dwarf galaxies appears to require the presence of a more massive neighbor. Had the Clouds been accreted early in the universe then, it is unlikely that these galaxies would currently have as much gas or star formation as they do.

Wetzel et al. (2012) used SDSS galaxy group/cluster catalogs together with N-body simulations to look at the quenching timescales of satellites at z ∼ 0. They find a quenching scenario in which satellite star formation histories are unaffected for 2–4 Gyr after infall but then quench rapidly (with an e-folding time of <0.8 Gyr). Interestingly, because of the time delay before quenching starts, satellites are found to experience significant stellar mass growth after infall, which the authors point to as a key reason for the success of the subhalo abundance matching technique: supporting the larger LMC masses that we favor in this study.

In conclusion, despite the revised lower estimates for the 3D velocities of the Magellanic Clouds presented in this work, taken together, Figures 9–12 make a strong case for a first infall scenario. This conclusion draws largely from a recognition that the fundamental change to our understanding of the orbital history of the Clouds comes not only from the PMs, but also from a cosmologically motivated understanding of the mass evolution and dark matter halo profile of our MW. This picture is thus consistent with the theory that the origin of the large scale gaseous structures of the Magellanic System (Stream, Bridge, and Leading Arm) and the internal structure and kinematics of the Clouds are a result of interactions between the LMC and SMC in a first infall scenario, rather than interactions with the MW (Besla et al. 2010, 2012).

The ACS reanalysis and the third epoch analysis give very consistent results at the level of the measured per-field PMs. In order to obtain COM PMs from these field PMs, the geometry and internal motions of the LMC/SMC need to be taken into account. The addition of the WFC3 epoch provides very small per-field PM errors of ∼0.03 mas yr⁻¹ (7 km s⁻¹, or 1.6% of the total PM). We are therefore able to independently constrain all parameters of the LMC PM rotation field, including the inclination, i, the position angle of the line-of-nodes, θ, the dynamical center, rotation velocity amplitude, and even the distance. This is a big improvement over K1. This procedure and the results are discussed in detail in Paper II. The features of import for the present study are that we obtain a rotation velocity that agrees very well with that obtained independently by Olsen et al. (2011) using a very large number of LOS velocities of stellar tracers. Moreover, we find strong evidence that the center of the LMC (as obtained from the PM fit to the rotation field) is consistent with the HI dynamical center, and not the center derived by vdM02, which agrees with the brightest part of the LMC bar.

For the SMC, the sparse coverage of QSO fields renders it difficult to gain much insight into its geometry. We therefore fit a relatively simple model to the SMC PM data, which allows for viewing perspective and a single overall rotation of the SMC in the plane of the sky. The old stellar population in the SMC shows little evidence for rotation (Harris & Zaritsky 2006) so we base our analysis of the SMC COM PM on the assumption that V_rot = 0 ± 15 km s⁻¹. Analyses in which V_rot is instead fit to the PM data do not yield a significantly different answer.

If we analyze the new LMC PM data with the same fixed geometric parameters as used in vdM02 we obtain a COM PM whose (1) value for μ_W is lower, but consistent within 1.6σ with the K1 result, (2) value for μ_N remains unchanged with respect to K1, (3) random errors are smaller by a factor of almost five than K1. If we analyze the PM data with the LMC model derived from the data itself, which is the preferred approach, then the value of μ_W and its random error remain largely unchanged. However, the value for μ_N changes significantly (by 4σ), and its random error increases by a factor of almost three. For the SMC, the three-epoch analysis agrees with K2 for μ_N, while for μ_W it differs by 2.4σ. This is largely due to how fields were combined in K2 and not due to intra-field differences. Our final SMC COM PM value is in rough agreement with the results of P08.

We attempt here to give as accurate an estimate as possible of the observational random errors in the COM PM, by propagating all unknowns in the geometry as well as the PM determination. This obviously increases errors with respect to previous studies which only propagated PM errors. Remarkably, the PM data are now no longer the dominant source of uncertainty, but rather uncertainties in the structure of the Clouds, dominate how well the COM PM can be established. Our final COM measurement errors, listed in line 1 of Table 4, do not represent a huge improvement over what was listed in previous HST works, despite the significantly improved accuracy of the present work. However, this is because those previous works underestimated the true random errors. Regardless, the uncertainty in the inferred transverse velocities of the Clouds are now dominated by their distance uncertainties, and not their COM PM uncertainties.

In order to turn the measured COM PMs into Galactocentric velocities, we need to know the solar velocity. Our understanding of the solar velocity has recently been revised upward, and this works to directly reduce the Galactocentric velocities of the Clouds by a comparable amount. Therefore, the Galactocentric v_tot = 321 ± 24 km s⁻¹ for the LMC presented here is 57 km s⁻¹ lower than that in K1, owing roughly equally to the decrease in μ_W described above and to the increase in solar velocity. There is a much smaller dependence on the new dynamical center derived here (see Table 5). The new center affects mainly the μ_N value, which determines the location of the orbit in projection on the sky rather than the tangential velocity (as described in B07). By contrast, if we use the same geometric model and solar velocity (IAU value) used in K1, v_tot decreases by 31 km s⁻¹ compared to the K1 value, which is comparable to the size of the 1σ error bar.

The SMC's inferred v_tot = 217 ± 26 km s⁻¹ is 85 km s⁻¹ lower than in K2. This is due mainly to the new PM derived here (which itself is due to how the fields were combined in K2, rather than differences in the derived per-field PMs). Roughly 20 km s⁻¹ of the decrease in v_tot is due to the revised solar velocity, and there is also a small contribution from the different SMC center used here, compared to the K2 analysis.

Given the new velocities obtained here, we have reevaluated the past orbital histories of the Clouds about the MW. We find that the dominant unknowns are the MW and LMC masses. However, some reasonable arguments allow us to narrow down the allowed ranges of MW and LMC masses. The results continue to make a strong case, as first argued in B07, that the Clouds are likely on their first infall into the MW.

From the search for bound orbits between the two Clouds (see Figure 12), we are able to identify a combination of masses that satisfy the criterion of long-lived binarity, and also produces orbits around the MW that are plausible from a cosmological point of view (Boylan-Kolchin et al. 2011; Busha et al. 2011a). This yields a preferred MW mass, ≲ 1.5 × 10¹² M_☉. The probability of a stable LMC– SMC binary configuration decreases as the MW mass increases. An LMC mass ≳ 1 × 10¹¹ M_☉ is needed in order to keep the SMC bound to the LMC for a reasonable fraction of a Hubble time. Taken together, this combination of MW and LMC masses imply that the LMC/the Clouds are on their first infall (see Figure 9).

From a cosmological point of view, these orbits do not represent a major alteration to the interpretation given in B07. Here we have explored a more full set of mass models for the LMC. The inferred LMC mass agrees with the expected total infall mass of the LMC, as calculated from the observed present-day stellar mass (vdM02) and the relations from halo occupation models (Guo et al. 2010; Sales et al. 2011). A generic outcome of our study is that regardless of what MW mass is used, large LMC masses are always on a first infall.

We have also implemented a simple model for the expected mass evolution of the MW over the past ∼10 Gyr. This, as might be expected, significantly increases the periods of LMC orbits that do make a previous pericentric passage (higher MW mass models), i.e., the orbits are highly eccentric. In a recent study, motivated by the discussion over the higher solar velocity, Zhang et al. (2012) have also revisited the Clouds' past orbital history, utilizing numerical simulations to explore the evolution of the MW. Like us here, they find that it is possible for the LMC to make past pericentric passages given the observational and theoretical error space. They do not vary LMC mass (it is kept fixed at 2 × 10¹⁰ M_☉), or investigate whether binarity imposes constraints, so we cannot directly compare our results with theirs, but were we to use such a low LMC mass, our studies would likely be consistent.

We put forth four arguments why we think orbits in which the LMC/SMC make a previous pericentric passage are implausible. One is the expected eccentricities of orbits from cosmological simulations along with a simple version of the timing argument (see Section 6.4), the second is the tidal radius of the LMC, and the third is the fact that long-lived SMC–LMC binary configurations have lower probability for such orbits, given the observed COM PMs of the Clouds. The fourth argument derives from the work of van den Bergh (2006) who conducted a morphological comparison of the satellites of the MW and M31, finding that the LMC/SMC are the only two gas-rich dIrrs at close Galactocentric distance to a host. He argues that they may be interlopers from a remote part of the Local Group rather than true satellites of the MW.

This line of argument is also consistent with the fact that LMC–SMC–MW analogs are relatively rare today in our local volume (Liu et al. 2011) and in cosmological simulations at z = 0 (Boylan-Kolchin et al. 2011). MW-type hosts are efficient at tidally disrupting such bound configurations. Furthermore, Tollerud et al. (2011) find that, compared to other R-band selected LMC–MW analogs in SDSS, the LMC is unusually blue in color. This fact is reconcilable in a first infall scenario where the LMC is accreted recently, because it has been able to retain enough gas to maintain a high star formation rate today (see also Wetzel et al. 2012). However, Tollerud et al. (2011) also find that the LMC's color is unusual compared to R-band selected LMC analogs in the field. Besla et al. (2012) present a theory in which the star formation rate of both the LMC and SMC has increased recently owing to interactions between these two galaxies. Given that LMC–SMC configurations are rare in general, it is thus unsurprising that the overall star formation rate of the LMC is higher than average isolated analogs.

Even in the case for the most long-lived LMC–MW configuration presented here (i.e., low LMC mass and high MW mass) the LMC makes at best ∼2 past pericentric passages (see Figure 9). The long implied orbital periods exceed the lifetime of the Magellanic Stream (1–2 Gyr; e.g., Heitsch & Putman 2009; Kereš & Hernquist 2009). This conclusion is consistent with our arguments for the need for a new formation mechanism for the Magellanic Stream that does not rely on a previous pericentric approach about the MW. Moreover, such a combination of LMC and MW masses make it difficult for the LMC and SMC to have been a long-lived binary system.

The existence of the LMC–SMC system as a binary for more than 2 Gyr in the past favors a combination of a high-mass LMC and an intermediate- or low-mass MW. Our work has found that such configurations always yield first infall orbital configurations. We show an example of an orbit that fulfills all the criteria set out above in Figure 13. This shows that viable orbital solutions can be found that are consistent with both the HST PMs and cosmological expectations. The new SMC PM now allows the past orbit of the SMC to cross-over the location of the HI in the Stream, which was very difficult to obtain with the old K2 PMs and is consistent with a picture in which the Stream is formed primarily by the removal of material from the SMC (GN96; Connors et al. 2006; Besla et al. 2010, 2012; Diaz & Bekki 2011, 2012).

Recently, the sample of QSOs behind the Magellanic Clouds has increased drastically. There are now well over 200 newly identified QSOs behind the LMC and 29 behind the SMC. The new QSOs were selected primarily based on their mid-IR colors (Kozłowski & Kochanek 2009). The majority of them were spectroscopically confirmed from follow-up data taken with 2DF/AAOMEGA (Kozłowski et al. 2011, 2012), while a smaller subset (around 40) were confirmed with Magellan/IMACS (by N. Kallivayalil & M. Geha). PM studies of fields around a carefully chosen subset of these QSOs could vastly improve our knowledge of the structure and dynamics of the Clouds.

With a larger sample of QSOs homogeneously distributed behind the SMC, we can investigate its rotation and structure, something we have been unable to tackle with only five QSO fields thus far (and with only three QSOs currently imaged with WFC3/UVIS). A larger number of QSOs will also improve our knowledge of its COM PM.

For the LMC, a large number of fields with PMs would allow us to better sample the rotation curve, including its inner slope. The increased coverage would also allow us to combine our PMs with radial velocity studies (e.g., Parisi et al. 2010; Olsen et al. 2011), thereby constituting a unique 3D data set. With a large number of QSO fields we may also be able to distinguish between the kinematics of different stellar populations within the LMC directly.

Finally, with exposures of sufficient depth, it is possible to measure absolute PMs with resolved background galaxies (Sohn et al. 2012), so we should not be confined to fields with QSOs. We are in the process of identifying fields in the Clouds with deep HST/ACS exposures in an earlier epoch, which should have numerous background galaxies that can be used to yield an accurate absolute PM of the stars in the field.