A Decade of MWC 758 Disk Images: Where Are the Spiral-arm-driving Planets?

Planets form in gaseous and dusty protoplanetary disks around young stars that are a few million years old. Forming planets gravitationally interact with the host disk, producing structures such as gaps, spiral arms, and vortices (Kley & Nelson 2012). By comparing observations with theoretical models, spatially resolved disk structures may yield rich information about the properties of embedded planets, such as their orbits, and dynamical constraints on their masses.

In the past decade, near-infrared imaging of disks with high spatial resolution has discovered spiral arms at tens of astronomical units in a few systems (e.g., SAO 206462: Muto et al. 2012; Garufi et al. 2013; Stolker et al. 2016; LkHα 330: Akiyama et al. 2016; MWC 758: Grady et al. 2013; Benisty et al. 2015; HD 100453: Wagner et al. 2015; Benisty et al. 2017; and HD 141569 A: Mouillet et al. 2001; Clampin et al. 2003; Konishi et al. 2016). Hydrodynamical and radiative transfer simulations have suggested two mechanisms for reproducing such structure: gravitational instability (Lodato & Rice 2005; Dong et al. 2015a), which occurs in disks with sufficient mass (Kratter & Lodato 2016), and companion-disk interaction (Dong et al. 2015b; Zhu et al. 2015; Bae et al. 2016). Because the host disks in these few systems are probably not massive enough to trigger the gravitational instability (e.g., Andrews et al. 2011), the latter scenario is more likely.

Detailed numerical simulations have quantified the dependencies of arm separation and contrast on the companion mass and disk properties (Fung & Dong 2015; Dong & Fung 2017). These relations have been used to infer the parameters of hypothesized arm-driving companions. A proof of concept of this mechanism has recently been provided by the HD 100453 system, where both the arms and the companion have been found, with their physical connections numerically supported (Dong et al. 2016b; Wagner et al. 2018). Extensive direct imaging observations have been carried out to look for the predicted arm-driving companions in a few other systems. Assuming hot-start planet formation models (e.g., Baraffe et al. 2015), they have generally ruled out all but planetary mass objects of a few Jupiter masses or less (e.g., Maire et al. 2017).

Companion-driven arms corotate with their driver. Therefore, by measuring their pattern speed, the orbital period, and thus semimajor axis, of their companion can be constrained (e.g., Lomax et al. 2016). We perform such an exercise for the spiral-arm system MWC 758, taking advantage of observations of the arms over a decade-long baseline established by a 2005 HST/NICMOS observation and 2015/2017 VLT/SPHERE and Keck/NIRC2 observations.

MWC 758 is a Herbig Ae star located at ${151}_{-9}^{+8}$ pc (Gaia Collaboration et al. 2016) with an age of 3.5 ± 2.0 Myr (Meeus et al. 2012), and mass of ∼2.0 M_⊙.¹⁵ The disk has a low inclination of ∼20° (Isella et al. 2010). Its two prominent, roughly 180° rotationally symmetric arms were first discovered with Subaru/HiCIAO (Grady et al. 2013), and subsequently characterized in detail with VLT/SPHERE (Benisty et al. 2015), with a third arm and point-source candidate at ∼011 (17 au) recently reported in Reggiani et al. (2017) using Keck/NIRC2. Numerical simulations by Dong et al. (2015a) suggested that both arms can be produced by a multi-Jupiter-mass planet at ∼06 from the star.

In this section, we describe the observations and data reduction for our 2005 HST and 2017 Keck programs.

2.1. NICMOS

The HST/NICMOS coronagraphic instrument observed the MWC 758 system in total intensity with the F110W filter (λ_cen = 1.1 μm) on 2005 January 7 (Proposal ID: 10177, PI: G. Schneider), and the unresolved disk morphology was presented in Grady et al. (2013). To retrieve the morphology of the spiral arms, we obtain calibrated NICMOS images of MWC 758 and another 814 reference star exposures, i.e., point-spread functions (PSFs), from the Archival Legacy Investigations for Circumstellar Environments (ALICE) project (PI: R. Soummer; Choquet et al. 2014; Hagan et al. 2018). We align the observations for better astrometry by employing a Radon-transform-based technique (Pueyo et al. 2015; Poteet et al. 2018), which focuses on the diffraction spikes in each exposure. To minimize color mismatch, telescope breathing, and cold mask alignment, we select the 81 closest PSFs in the L²-norm sense, and perform PSF subtraction with the Non-negative Matrix Factorization (NMF) method (Ren et al. 2018), which is shown to preserve the morphology of circumstellar disks better, especially in reference differential imaging scenarios.

In Figure 1, we present the reduction results of the NICMOS observations at two telescope orientations (three exposures each) and their signal-to-noise (S/N) maps. We argue the physical existence of the detection since the spiral pattern is (1) consistent within the same telescope orientation, as shown from the S/N maps, which are calculated from dividing the combined result by the pixel-wise standard deviation of the ones constituing them; (2) independent of telescope orientation (30° separation), which excludes the scenario of unsuccessful noise removal; (3) not depending on the number of NMF components, reducing the possibility of underfitting and overfitting; (4) not resembling the null detections in the ALICE archive, as well as a reduction consistency using a principal-component-analysis-based reduction method (Soummer et al. 2012).

Zoom In Zoom Out Reset image size

The integrated flux for 03 < r < 05 is 2.0 ± 0.8 mJy at 1σ level, consistent with the upper limit reported in Grady et al. (2013). We notice flux variations between the two telescope orientations; however, we do not address the origin of this difference in this letter, but focus on the morphology only.

2.2. NIRC2

2.3. SPHERE

To measure the pattern speed of the arms, we first scale the surface brightnesses of the reduced NICMOS, NIRC2, and SPHERE images by the distance-dependent factor r². The results in Cartesian and polar coordinates are presented in Figure 2.

Zoom In Zoom Out Reset image size

There are three main differences among our observational data sets: (1) the pixel size of the NICMOS instrument is ∼8 times larger than the other two; (2) the NICMOS and NIRC2 observations measure the total intensity, while SPHERE traces the polarized light, and (3) the NICMOS and SPHERE observations are at ∼1 μm, while the NIRC2 observation is at ∼3.8 μm. For (1), we interpolate the NICMOS image to match the pixel scales of the others. For (2), we reduced the SPHERE ADI total intensity observation, compared it with the polarized map, and found no discernible discrepancy; this agreement is also endorsed by the simulation in Dong et al. (2016a). For (3), we compare observations at roughly the same central wavelength (λ_cen).

3.1. Measurement of Rotation of Spirals

3.1.1. NIRC2: 2015 versus 2017 (1.28 years apart)

To mitigate any systematic offset between instruments, and provide an initial constraint on the rotation, we compare two observations from the same Keck/NIRC2 instrument: our 2017 February 02 observation, and the 2015 October 24 observation (PI: E. Serabyn, Reggiani et al. 2017), which is aligned and reduced with an identical procedure.

We quantify the rotation of the spiral arms as follows: in polar coordinates, we first fit Gaussian profiles to the brightness of the spiral arms at a fixed radial separation; then for each arm, we perform the weighted Least Square Dummy Variable (LSDV, Appendix) analysis to fit the same morphological profile in both epochs to obtain their relative rotation. For the southern and northern primary arms (marked with "1" and "3" in Figure 3), we obtain a rotation of ${\rm{\Delta }}{\theta }_{S}^{(1.28\,\mathrm{year})}=0\buildrel{\circ}\over{.} 77\pm 10\buildrel{\circ}\over{.} 65$,¹⁶ and ${\rm{\Delta }}{\theta }_{N}^{(1.28\,\mathrm{year})}=-0\buildrel{\circ}\over{.} 70\pm 6\buildrel{\circ}\over{.} 78$, respectively. Since spiral arms in disks are trailing patterns, the MWC 758 arms are expected to rotate in a clockwise direction, i.e., Δθ ≥ 0, we therefore adopt the constraints from the southern primary arm, $\dot{\theta }={\rm{\Delta }}{\theta }_{S}^{(1.28\,\mathrm{year})}/1.28\,\mathrm{year}=0\buildrel{\circ}\over{.} 6\pm 8\buildrel{\circ}\over{.} 3$ yr⁻¹, as the rotation of the two.

Zoom In Zoom Out Reset image size

For a consistency check, we measure the rotation with another method: in polar coordinates, we obtain the cross-correlation maps (Tonry & Davis 1979), and measure $\dot{\theta }\,=0\buildrel{\circ}\over{.} 7\pm 56\buildrel{\circ}\over{.} 2$ yr⁻¹. We adopt the value from the LSDV method, since it is less biased by the nonspiral structures in the entire field of view as the cross-correlation method, and the best-fit values agree within 01.

3.1.2. 2005 NICMOS versus 2015 SPHERE (10.17 years apart)

By analyzing the rotation between the NICMOS and SPHERE images, we narrow down the uncertainty determined from the NIRC2 observations here. We adopt the NIRC2 best-fit and uncertainty values, rotate the SPHERE result back to the NICMOS epoch, then the two images are expected to have no azimuthal shift. By fitting identical profiles and offsets for the rotated SPHERE and original NICMOS observations, and given the northern arm is blended with its secondary arm (marked by "4" in Figure 3) but with a smaller uncertainty than the southern arm, we adopt the results from the southern arm, obtaining a conservative measurement of ${\rm{\Delta }}{\theta }^{(10.17\mathrm{year})}=6\buildrel{\circ}\over{.} 1\pm 29\buildrel{\circ}\over{.} 4$. This corresponds to a statistical uncertainty for the angular speed: ${(\delta \dot{\theta })}_{\mathrm{statistical}}=29\buildrel{\circ}\over{.} 4/10.17\,\mathrm{year}=2\buildrel{\circ}\over{.} 89$ yr⁻¹.

3.2. Additional Systematics

From the previous analyses of statistical and systematical uncertainties, we obtain a total uncertainty in the rotation of the arms at 3σ:

$$\begin{eqnarray*}\delta \dot{\theta } & = & \sqrt{(\delta \dot{\theta }{)}_{\mathrm{statistical}}^{2}+(\delta \dot{\theta }{)}_{\mathrm{systematic}}^{2}}\\ & = & \sqrt{(\delta \dot{\theta }{)}_{\mathrm{statistical}}^{2}+[(\delta \dot{\theta }{)}_{\mathrm{CD}}^{2}+(\delta \dot{\theta }{)}_{\mathrm{FPM}}^{2}+(\delta \dot{\theta }{)}_{\mathrm{pixel}}^{2}]}\\ & = & \sqrt{{2.89}^{2}+({0.06}^{2}+{1.08}^{2}+{1.18}^{2})}=3\buildrel{\circ}\over{.} 31\,{\mathrm{yr}}^{-1},\end{eqnarray*}$$

where the pixel uncertainty is accounted for the NICMOS pixel size of 12° at r ∼ 06 (i.e., ($\delta {\dot{\theta })}_{\mathrm{pixel}}=12^\circ /10.17\,\mathrm{year}\,=1\buildrel{\circ}\over{.} 18$ yr⁻¹). Together with the best-fit value, we obtain

$$\begin{eqnarray}&&\dot{\theta }=0\buildrel{\circ}\over{.} {6}_{-0\buildrel{\circ}\over{.} \,6}^{+3\buildrel{\circ}\over{.} \,3}\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 1 }$$

where the lower limit is physically constrained from the clockwise rotation of the MWC 758 arms.

For hypothesized arm-driving planet(s) on a circular orbit (eccentricity e = 0), the best-fit pattern speed corresponds to a period of T = 598 years, or a radial separation of r_best = 89 au; and the 3σ upper limit leads to T = 92 years and r_3σ = 26 au; see Figure 4 for graphical representations.

Zoom In Zoom Out Reset image size

For e > 0, if the planet reaches its apogee in the current epoch, the limit on the arm motion translates into a stellocentric separation ${r}_{\min }{(e)=(1-e)}^{1/3}r$ au. For giant planets of several Jupiter mass forming in and interacting with a gaseous disk, their eccentricities are not expected to grow beyond ∼0.1 (Dunhill et al. 2013; Duffell & Chiang 2015). Furthermore, density waves excited by Jovian planets with e ≳ 0.2 start to deviate from their usual morphology, as the waves launched at different epochs interact with each other (H. Li & S. Li 2018, private commnication), which provide poor fits to the arms around MWC 758. In this case, e = 0.2 leads to a 7% decrease in the minimum stellocentric separation, and the hypothetical arm-driving planet(s) should be located at least 25 au from the star in 2017 at 3σ.

One might consider another scenario in which the two spiral arms exchanged positions between 2005 and 2015, i.e., rotated ∼180°. However, this means the major arms should have rotated an additional 226 between the 2015 and 2017 NIRC2 observations, this is ruled out at 6σ with our constraints.

We present reprocessed 2005 HST/NICMOS observations of the MWC 758 disk, and successfully retrieve the two spiral arms in the system revealed by ground-based high-contrast imaging facilities. Thanks to the HST image, we are able to establish a 10 year baseline in observations to constrain the pattern speed of the major arms. Together with a 2015 VLT/SPHERE data set, and two Keck/NIRC2 observations in 2015 and 2017, we obtain a rotation speed of $0\buildrel{\circ}\over{.} {6}_{-0\buildrel{\circ}\over{.} \,6}^{+3\buildrel{\circ}\over{.} \,3}\,{\mathrm{yr}}^{-1}$ at 3σ for the two major spiral arms. The results correspond to a best-fit value of 89 au (059), and a 3σ lower limit of 26 au (017), for the orbital distance of the hypothesized arm-driving perturber on a circular orbit.

Our constraint is consistent with the Dong et al. (2015b) prediction of the arm-driving planet at ∼90 au (∼06). In addition, we rule out at a 6σ level the scenario that the companion candidate at 011 (17 au) reported by Reggiani et al. (2017) is driving the major two spiral arms, assuming the candidate is on a circular orbit coplanar with the arms. This is further supported by Bae & Zhu (2017), that a Jovian planet can drive only one external arm assuming a reasonable disk scale height. For new arm(s) and planet candidate(s) (e.g., Reggiani et al. 2017), more observations are needed to confirm their existence and dynamical connections.

The possible arm-driving planets in the MWC 758 system are excellent targets for future observations in direct imaging observations both from the ground and with the James Webb Space Telescope, and in millimeter observations using ALMA to search for evidence of a circumplanetary disk (e.g., Eisner 2015; Perez et al. 2015; Zhu 2015; Szulágyi et al. 2018).

We are grateful to the anonymous referee for constructive suggestions that improved our Letter, and Sean Brittain for insightful discussions. B.R. acknowledges the computational resources from the Maryland Advanced Research Computing Center (MARCC), which is funded by a State of Maryland grant to Johns Hopkins University through the Institute for Data Intensive Engineering and Science (IDIES). E.C. acknowledges support from NASA through Hubble Fellowship grant HST-HF2-51355 awarded by STScI, operated by AURA, Inc. under contract NAS5-26555, and support from HST-AR-12652, for research carried out at the Jet Propulsion Laboratory, California Institute of Technology. T.E. was supported in part by NASA Grants NNX15AD95G/NEXSS, NNX15AC89G, and NSF AST-1518332. This research has made use of data reprocessed as part of the ALICE program, which was supported by NASA through grants HST-AR-12652 (PI: R. Soummer), HST-GO-11136 (PI: D. Golimowski), HST-GO-13855 (PI: E. Choquet), HST-GO-13331 (PI: L. Pueyo), and STScI Director's Discretionary Research funds, and was conducted at STScI, which is operated by AURA under NASA contract NAS5-26555. The input images to ALICE processing are from the recalibrated NICMOS data products produced by the Legacy Archive project, "A Legacy Archive PSF Library And Circumstellar Environments (LAPLACE) Investigation," (HST-AR-11279, PI: G. Schneider). This work benefited from NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea, host to the W. M. Keck Observatory, has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

When the morphological patterns of the spiral arms do not change among different epochs, the only difference is their relative azimuthal offset. This is known as the "fixed effect" in statistics, which has been extensively studied with the classical LSDV method. LSDV generalizes the ordinary least square (OLS) method using dummy variables. In this section, we first describe the classical LSDV method, then introduce our generalization of the method by taking into account the uncertainty from input data.

A.1. The Classical LSDV Method

In polar coordinates, the location of spiral arms can be represented by $({\theta }_{\mathrm{PA}}\pm \delta {\theta }_{\mathrm{PA}},r)$ pairs, which represents the position angle (with its uncertainty) and radial separation. For a given radial separation, its position angle can be estimated from fitting Gaussian profiles at different azimuthal directions. Assume there are E epochs, each has n data pairs, we can use a Taylor polynomial of degree p to represent the arm morphology (e.g., Grady et al. 2013; Benisty et al. 2015; Reggiani et al. 2017). The classical LSDV method finds for all the data the best fit of

$$\begin{eqnarray}&&{\theta }_{\mathrm{PA},i}=f({r}_{i})=\displaystyle \sum _{j=1}^{p}{c}_{j}{r}_{i}^{j}+\displaystyle \sum _{k=1}^{E}{d}_{k}{D}_{k}(i),\end{eqnarray} \tag{ 2 }$$

where the dummy variables D_k(i) = 1 only when the $({\theta }_{\mathrm{PA},i},{r}_{i})$ pair is obtained from epoch k, and 0 otherwise. The coefficients d are then the position angles of the spiral arms when r = 0.

Let set ${{\mathbb{R}}}^{s\times t}$ contain s-by-t real-valued matrices, if we denote the θ_PA's by ${\boldsymbol{\Theta }}\in {{\mathbb{R}}}^{{nE}\times 1}$, with ${{\boldsymbol{\Theta }}}_{i}={\theta }_{\mathrm{PA},i};$ the r's and dummy variables by ${\boldsymbol{R}}\in {{\mathbb{R}}}^{{nE}\times (p+E)}$, with ${{\boldsymbol{R}}}_{i(\cdot )}\,=[{r}_{i},{r}_{i}^{2},\cdots ,\,{r}_{i}^{p},{D}_{1}(i),{D}_{2}(i),\cdots ,\,{D}_{E}(i)];$ and the coefficients ${\boldsymbol{\beta }}={[{c}_{1},\cdots ,{c}_{p},{d}_{1},\cdots ,{d}_{E}]}^{T}\in {{\mathbb{R}}}^{(p+E)\times 1}$. We now write Equation (2) in a matrix OLS form:

$$\begin{eqnarray}&&{\boldsymbol{\Theta }}={\boldsymbol{R}}{\boldsymbol{\beta }}+{\boldsymbol{\epsilon }},\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{\epsilon }}\in {{\mathbb{R}}}^{{nE}\times 1}$ is the residual. Its cost function,

$$\begin{eqnarray}C({\boldsymbol{\Theta }},{\boldsymbol{R}};{\boldsymbol{\beta }}) & = & {{\boldsymbol{\epsilon }}}^{T}{\boldsymbol{\epsilon }}\\ & = & {({\boldsymbol{\Theta }}-{\boldsymbol{R}}{\boldsymbol{\beta }})}^{T}({\boldsymbol{\Theta }}-{\boldsymbol{R}}{\boldsymbol{\beta }})\\ & = & \displaystyle \sum _{i=1}^{{nE}}{\left[{\theta }_{\mathrm{PA},i}-\left(\displaystyle \sum _{j=1}^{p}{c}_{j}{r}_{i}^{j}+\displaystyle \sum _{k=1}^{{nE}}{d}_{k}{D}_{k}(i)\right)\right]}^{2},\end{eqnarray} \tag{ 4 }$$

is minimized by

$$\begin{eqnarray}&&\hat{{\boldsymbol{\beta }}}={({{\boldsymbol{R}}}^{T}{\boldsymbol{R}})}^{-1}{{\boldsymbol{R}}}^{T}{\boldsymbol{\Theta }},\end{eqnarray} \tag{ 5 }$$

where ^T and ⁻¹ stand for matrix transpose and inverse. The standard deviations of $\hat{{\boldsymbol{\beta }}}$ are calculated from the element-wise square root of the diagonal elements in the variance-covariance matrix:

$$\begin{eqnarray}\delta \hat{{\boldsymbol{\beta }}} & = & \sqrt{{\delta }^{2}\hat{{\boldsymbol{\beta }}}}=\sqrt{\mathrm{diag}\{E[{(\hat{{\boldsymbol{\beta }}}-{\boldsymbol{\beta }})(\hat{{\boldsymbol{\beta }}}-{\boldsymbol{\beta }})}^{T}]\}}\\ & = & \sqrt{\mathrm{diag}\{{\hat{\sigma }}^{2}{({{\boldsymbol{R}}}^{T}{\boldsymbol{R}})}^{-1}\}},\end{eqnarray} \tag{ 6 }$$

where ${\hat{\sigma }}^{2}=\tfrac{{({\boldsymbol{\Theta }}-{\boldsymbol{R}}\hat{{\boldsymbol{\beta }}})}^{T}({\boldsymbol{\Theta }}-{\boldsymbol{R}}\hat{{\boldsymbol{\beta }}})}{{nE}-(p+E)}.$

A.2. The Weighted LSDV Method

To take into account the measurement uncertainty in our study, we generalize the classical LSDV method into a weighted form. The weighted LSDV method minimizes the chi-squared statistic:

$$\begin{eqnarray}{\chi }^{2}({\boldsymbol{\Theta }},{\boldsymbol{R}};{\boldsymbol{\beta }}) & = & {\left(\displaystyle \frac{{\boldsymbol{\Theta }}-{\boldsymbol{R}}{\boldsymbol{\beta }}}{\delta {\boldsymbol{\Theta }}}\right)}^{T}\left(\displaystyle \frac{{\boldsymbol{\Theta }}-{\boldsymbol{R}}{\boldsymbol{\beta }}}{\delta {\boldsymbol{\Theta }}}\right)\\ & = & \displaystyle \sum _{i=1}^{{nE}}{\left[\displaystyle \frac{{\theta }_{\mathrm{PA},i}}{{\delta }_{{\theta }_{\mathrm{PA},i}}}-\left(\displaystyle \sum _{j=1}^{p}{c}_{j}\displaystyle \frac{{r}_{i}^{j}}{{\delta }_{{\theta }_{\mathrm{PA},i}}}+\displaystyle \sum _{k=1}^{{nE}}{d}_{k}\displaystyle \frac{{D}_{k}(i)}{{\delta }_{{\theta }_{\mathrm{PA},i}}}\right)\right]}^{2},\end{eqnarray} \tag{ 7 }$$

where the division operation is element-wise; and $\delta {\boldsymbol{\Theta }}\in {{\mathbb{R}}}^{{nE}\times 1}$ stores the uncertainty for ${\boldsymbol{\Theta }}$. With substitution

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\theta }_{\mathrm{PA},i}^{{\prime} }=\displaystyle \frac{{\theta }_{\mathrm{PA},i}}{{\delta }_{{\theta }_{\mathrm{PA},i}}},\\ r{{\prime} }_{i}^{j}=\displaystyle \frac{{r}_{i}^{j}}{{\delta }_{{\theta }_{\mathrm{PA},i}}},\\ D{{\prime} }_{k}(i)=\displaystyle \frac{{D}_{k}(i)}{{\delta }_{{\theta }_{\mathrm{PA},i}}},\end{array}\right.\end{eqnarray} \tag{ 8 }$$

where ' denotes the (element-wise) division of $\delta {\theta }_{\mathrm{PA},i}$, we have a matrix form of

$$\begin{eqnarray}{\chi }^{2}({\boldsymbol{\Theta }},{\boldsymbol{R}};{\boldsymbol{\beta }}) & = & \displaystyle \sum _{i=1}^{{nE}}{\left[\theta {{\prime} }_{\mathrm{PA},i}-\left(\displaystyle \sum _{j=1}^{p}{c}_{j}r{{\prime} }_{i}^{j}+\displaystyle \sum _{k=1}^{{nE}}{d}_{k}D{{\prime} }_{k}(i)\right)\right]}^{2}\\ & = & {({\boldsymbol{\Theta }}^{\prime} -{\boldsymbol{R}}^{\prime} {\boldsymbol{\beta }})}^{T}({\boldsymbol{\Theta }}^{\prime} -{\boldsymbol{R}}^{\prime} {\boldsymbol{\beta }})\\ & = & C({\boldsymbol{\Theta }}^{\prime} ,{\boldsymbol{R}}^{\prime} ;{\boldsymbol{\beta }}),\end{eqnarray} \tag{ 9 }$$

as in Equation (4), whose best-fit values and standard deviations can thus be obtained from Equations (5) and (6).

With two epochs of observations in our fitting, we have E = 2 and focus only on the difference of the last two terms in $\hat{{\boldsymbol{\beta }}}$, i.e., ${\rm{\Delta }}d={\hat{d}}_{1}-{\hat{d}}_{2}$, as well as the uncertainty $\delta ({\rm{\Delta }}d)\,=\sqrt{{\delta }^{2}{\hat{d}}_{1}+{\delta }^{2}{\hat{d}}_{2}-2\mathrm{Cov}({\hat{d}}_{1},{\hat{d}}_{2})}$. In our fitting efforts, we investigated Taylor polynomials up to p = 3 degrees to study the different morphological pattern of the spiral arms; however, no significant difference was observed, we therefore only report the linear results in this Letter.