  The 3D velocities of M31 and M33 are important for understanding the
evolution and cosmological context of the Local Group. Their most massive stars
are detected by Gaia, and we use Data Release 2 (DR2) to determine the galaxy
proper motions (PMs). We select galaxy members based on, e.g., parallax, PM,
color-magnitude-diagram location, and local stellar density. The PM rotation of
both galaxies is confidently detected, consistent with the known line-of-sight
rotation curves: $V_{\rm rot} = -206\pm86$ km s$^{-1}$ (counter-clockwise) for
M31, and $V_{\rm rot} = 80\pm52$ km s$^{-1}$ (clockwise) for M33. We measure
the center-of-mass PM of each galaxy relative to surrounding background quasars
in DR2. This yields that $({\mu}_{\alpha*},{\mu}_{\delta})$ equals $(65 \pm 18
, -57 \pm 15)$ $\mu$as yr$^{-1}$ for M31, and $(31 \pm 19 , -29 \pm 16)$
$\mu$as yr$^{-1}$ for M33. In addition to the listed random errors, each
component has an additional residual systematic error of 16 $\mu$as yr$^{-1}$.
These results are consistent at 0.8$\sigma$ and 1.0$\sigma$ with the (2 and 3
times higher-accuracy) measurements already available from Hubble Space
Telescope (HST) optical imaging and VLBA water maser observations,
respectively. This lends confidence that all these measurements are robust. The
new results imply that the M31 orbit towards the Milky Way is somewhat less
radial than previously inferred, $V_{\rm tan, DR2+HST} = 57^{+35}_{-31}$ km
s$^{-1}$, and strengthen arguments that M33 may be on its first infall into
M31. The results highlight the future potential of Gaia for PM studies beyond
the Milky Way satellite system.
