# Dynamical Evidence of a Spiral Arm-driving Planet in the MWC 758 Protoplanetary Disk

Bin Ren, Ruobing Dong, Rob van Holstein, Jean-Baptiste Ruffio, Benjamin Calvin, Julien Girard, Myriam Benisty, Anthony Boccaletti, Thomas Esposito, Elodie Choquet, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Bin Ren, Ruobing Dong, Rob van Holstein, Jean-Baptiste Ruffio, Benjamin Calvin, et al., Dynamical Evidence of a Spiral Arm-driving Planet in the MWC 758 Protoplanetary Disk. The Astrophysical Journal Letters, 2020, 898 (2), pp.L38. 10.3847/2041-8213/aba43e hal-03083224

HAL Id: hal-03083224

https://hal.science/hal-03083224v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Dynamical Evidence of a Spiral Arm-driving Planet in the MWC 758 Protoplanetary Disk

Bin Ren (任彬)<fnref n="1" />, Robbing Dong (董若冰)<fnref n="2" />, Rob G. van Holstein<fnref n="3" />,<fnref n="4" />, Jean-Baptiste Ruffio<fnref n="1" />, Benjamin A. Calvin<fnref n="1" />, Julien H. Girard<fnref n="5" />, Myriam Benisty<fnref n="6" />,<fnref n="7" />, Anthony Boccaletti<fnref n="8" />, Thomas M. Esposito<fnref n="9" />, Élodie Choquet<fnref n="10" />, Dimitri Mawet<fnref n="1" />, Laurent Pueyo<fnref n="5" />, Tomas Stolker<fnref n="11" />, Eugene Chiang<fnref n="9" />, Jozua de Boer<fnref n="3" />, John H. Debes<fnref n="5" />, Antonio Garufi<fnref n="12" />, Carol A. Grady<fnref n="13" />, Dean C. Hines<fnref n="5" />, Anne-Lise Maire<fnref n="14" />, François Ménard<fnref n="6" />, Maxwell A. Millar-Blanchaer<fnref n="1" />, Marshall D. Perrin<fnref n="5" />, Charles A. Poteet<fnref n="5" />, and Glenn Schneider<fnref n="15" />

Received 2020 May 8; revised 2020 July 8; accepted 2020 July 8; published 2020 July 29

## Abstract

More than a dozen young stars host spiral arms in their surrounding protoplanetary disks. The excitation mechanisms of such arms are under debate. The two leading hypotheses—companion-disk interaction and gravitational instability (GI)—predict distinct motion for spirals. By imaging the MWC 758 spiral arm system at two epochs spanning ~5 yr using the SPHERE instrument on the Very Large Telescope, we test the two hypotheses for the first time. We find that the pattern speeds of the spirals are not consistent with the GI origin. Our measurements further evince the existence of a faint “missing planet” driving the disk arms. The average spiral pattern speed is $0^{\circ}22 \pm 0^{\circ}03 \text{ yr}^{-1}$, pointing to a driver at $172_{-14}^{+18}$ au around a $1.9 M_{\oplus}$ central star if it is on a circular orbit. In addition, we witness time-varying shadowing effects on a global scale that are likely originating from an inner disk.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Coronagraphic imaging (313); Planetary system formation (1257); Orbital motion (1179)*

*Supporting material: data behind figure*

## 1. Introduction

Spiral arms, spanning from tens to hundreds of astronomical units, are found in more than a dozen protoplanetary disks in visible to near-infrared light with high-contrast imaging (e.g., Grady et al. 1999, 2013; Muto et al. 2012; Wagner et al. 2015; Monnier et al. 2019; Garufi et al. 2020; Ménard et al. 2020; Muro-Arena et al. 2020). Their origin has profound implications for both planet formation and disk evolution (Dong et al. 2018; Brittain et al. 2020). In the companion-disk interaction scenario (Kley & Nelson 2012; Dong et al. 2015b; Zhu et al. 2015; Bae et al. 2016), the masses and locations of the drivers can be inferred (Fung & Dong 2015; Dong & Fung 2017), while in the gravitational instability (GI) scenario (Lodato & Rice 2005; Dong et al. 2015a; Kratter & Lodato 2016), the disk masses can be constrained. To test the two hypotheses, great effort has been expended to search for faint companions in disks, and to accurately measure disk masses. However, both approaches are notoriously difficult. As a result, only the spiral

arms in two systems have been confirmed to be induced by stellar companions (HD 100453: Rosotti et al. 2020; UXTau: Ménard et al. 2020).

The pattern speed of the spirals provides an exciting route to test the hypotheses: the arms corotate with the driver in the companion scenario, and undergo local Keplerian motion on timescales much shorter than the dynamical timescale in the GI scenario. Long temporal baseline and high-resolution imaging are needed to accurately assess the motion of the spirals (Ren et al. 2018, hereafter R18). Thus far, there has been no observational attempt to dynamically test the two arm formation, and thus motion, hypotheses. In this Letter, we image the MWC 758 protoplanetary disk with a 4.71 yr baseline to investigate the change in the brightness and the motion of the spiral arms between two observations, and dynamically test the two arm motion hypotheses.

Located at $160.2 \pm 1.7$ pc (Gaia Collaboration et al. 2018), the A8V Herbig star MWC 758 (Vieira et al. 2003) is $10.9_{-1.0}^{+12.0}$ Myr

<footnote n="1">¹Department of Astronomy, California Institute of Technology, 1216 East California Boulevard, Pasadena, CA 91125, USA; ren@caltech.edu</footnote>

<footnote n="2">²Department of Physics & Astronomy, University of Victoria, Victoria, BC V8P 1A1, Canada</footnote>

<footnote n="3">³Leiden Observatory, Universiteit Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4">⁴European Southern Observatory, Alonso de Córdova 3107, Casilla 19001, Vitacura, Santiago, Chile</footnote>

<footnote n="5">⁵Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="6">⁶Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="7">⁷Unidad Mixta Internacional Franco-Chilena de Astronomía (CNRS, UMI 3386), Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="8">⁸LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC, Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="9">⁹Department of Astronomy, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="10">¹⁰Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="11">¹¹Institute for Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland</footnote>

<footnote n="12">¹²INAF, Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="13">¹³Exoplanets and Stellar Astrophysics Laboratory, Code 667, Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="14">¹⁴Space Sciences, Technologies, and Astrophysics Research (STAR) Institute, Université de Liège, Liège, Belgium</footnote>

<footnote n="15">¹⁵Steward Observatory, The University of Arizona, Tucson, AZ 85721, USA</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** VLT/SPHERE observations of MWC 758 spiral arms; (a) and (b) are the 2015 and 2019 scaled surface brightness Stokes Q<sub>φ</sub> images; (c) is the surface brightness ratio from dividing (b) by (a). In 2019, most of the southern region is more than twice the corresponding brightness of 2015; the northern region is less than half except for a brightened clump. The 184 mas diameter coronagraph used in 2019 blocks light in the central circular regions.</figcaption></figure>

(The data used to create this figure are available.)

old with an estimated mass of 1.9 ± 0.2 M<sub>☉</sub> (Garufi et al. 2018)<sup>16</sup>. It hosts a protoplanetary disk with two prominent spiral arms in near-infrared light (Grady et al. 2013; Benisty et al. 2015; Reggiani et al. 2018; R18; Wagner et al. 2019). The spiral arms have been hypothesized to be driven by GI (Dong et al. 2015a) or a planetary companion (Dong et al. 2015b; Baruteau et al. 2019). High-contrast imaging searches have identified two candidates in the disk (Reggiani et al. 2018; Wagner et al. 2019), and their nature is still under investigation. Under the planet driver assumption, the inner candidate in Reggiani et al. (2018) has been ruled out as the arm driver by the motion measurements in R18, unless it is on a highly eccentric orbit (Calcino et al. 2020).

## 2. Observation and Data Reduction

We observe MWC 758 at two epochs using the infrared dual-band imager and spectrograph on SPHERE (Beuzit et al. 2019) with the differential polarimetric imaging mode at Y-band (1.04 µm). The first epoch is on 2015 March 4 under ESO program 60.A-9389(A) (PI: M. Benisty; Benisty et al. 2015). The second epoch is on 2019 November 18 under ESO program 104.C-0472(A) (PI: B. Ren).

In both epochs, the detector integration time is 32 s per frame. We obtain the 2015 data in field-tracking mode using the apodized Lyot coronagraph with apodizer APO2, which is optimized for 5.2λ/D focal masks, and the 145 mas diameter Lyot mask ALC1 (coronagraph combination name: N_ALC_Y, inner working angle: IWA = 72.5 mas, 1 pixel is 12.25 mas; Maire et al. 2016). We have four polarimetric cycles. In each cycle, the half-wave plate (HWP) cycles through switch angles 0°, 22°, 5°, 45° and 67°, 5° to measure Stokes Q and U. At each HWP position, there are six integrations. The atmospheric seeing (as measured by the differential image motion monitor) is 1.′′08 ± 0.′′17<sup>17</sup> and the coherence time is 3 ms. We obtain the

2019 data in pupil-tracking mode to clean and stabilize the diffraction pattern using the apodized Lyot coronagraph using apodizer APO1, which is optimized for 4λ/D focal masks, and the 185 mas diameter Lyot mask ALC2 (coronagraph combination name: N_ALC_YJH_S, IWA = 92.5 mas). We have 10 polarimetric cycles with HWP switch angles 0°, 45°, 22°, 5° and 67°, 5°. There are two integrations at each HWP position. The seeing and coherence time are 0.′′63 ± 0.′′06 and 5 ms, respectively. In our 1 hr observation blocks, we total 3072 s on-source integration time in 2015 and 2569 s in 2019, with 128 s and 32 s on-sky time at the end of each observation, respectively.

We reduce the two data sets using the IRDAP data-reduction pipeline (van Holstein et al. 2017, 2020) which employs a fully validated Mueller matrix model to minimize reduction bias. We use the images with star polarization subtracted for our analysis. Specifically, we use the Q<sub>φ</sub> images that show the light polarized parallel or perpendicular to the radial direction from the star, and trace the dust particles on the surface of a disk (Monnier et al. 2019). We measure that the flux in the central star’s point-spread function halo in the total intensity images of the 2015 data is 90% ± 2% of that in the 2019 data. We therefore divide the 2015 Q<sub>φ</sub> image by 0.9 to minimize effects from the central star illumination and/or observation conditions.

We scale the surface brightness distribution for the two Q<sub>φ</sub> images for comparison and analysis. First, we deproject the images to face-on view assuming an inclination of 21° and a position angle of 62° for the disk (Isella et al. 2010; Boehler et al. 2018). Next, we compute the mid-plane stellocentric distances (r) for all pixels. We multiply the value at each pixel by (r/r<sub>0</sub>)<sup>2</sup>, where r<sub>0</sub> = 0.′′5, to enhance the visibility of features at large distances. We present the resulting surface brightness maps at the two epochs in Figure 1, as well as the ratio between the two. We mark the star location with a white cross, whose longer axis is aligned to the disk major axis.

## 3. Analysis

### 3.1. Shadowing Effects

We identify time-varying illumination patterns on a global scale in Figure 1. Since scattered light probes the disk surface,

<footnote n="16"><sup>16</sup> The previous age estimate of 3.5 ± 2.0 Myr (Vieira et al. 2003) has been updated after Gaia DR2 by, e.g., Garufi et al. (2018) and Vioque et al. (2018) who also calculated the star mass. Here we adopt the Garufi et al. (2018) values since their star mass is consistent with the CO line emission measurement by Isella et al. (2010).</footnote>

<footnote n="17"><sup>17</sup> The uncertainties in this Letter are 1σ unless otherwise specified.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Arm location pairs in polar coordinates. The error bars are the peak positions of spiral arms in the form of $(\theta, r)$ pairs in the deprojected version of Figure 1, where $\theta$ is the counter-clockwise deviation from the northeast semimajor axis. The colored lines are the best-fit companion-driven model to the observations. The angular offset between the colored lines is $1^{\circ}04 \pm 0^{\circ}14$. The observation does not follow the gravitational instability prediction for a $1.9 \pm 0.2 M_{\odot}$ central star.</figcaption></figure>

we interpret the variations as moving shadows. Similar large-scale shadowing effects have only been observed for TW Hya in Debes et al. (2017). For MWC 758, the southern region is twice as bright in 2019 than in 2015, while the northern arm on the opposite side (~1 o'clock) has dimmed by a factor of two.

In scattered light, the brightness variation at a location may be caused by a change in the shadow casting inner disk structure at the same azimuth, resulting from mechanisms such as the precessing of an inner disk behind the coronagraphic mask (Nealon et al. 2019) or fluctuations in the dust structure in the inner disk arising from dust dynamics (Stolker et al. 2017). We detect a central source degree of linear polarization of $0.50\% \pm 0.06\%$ and an angle of linear polarization of $126^\circ \pm 3^\circ$ using IRDAP, which could originate from an inclined inner disk (van Holstein et al. 2020). A clump right outside the north edge of the coronagraph (marked in panel (c) in Figure 1) has brightened by a factor of two from 2015 to 2019. Given that its brightening coincides with the fainting of the spiral arm tip in the northern region at a larger radii, we hypothesize that the clump may be shadowing the outer disk.

## 3.2. Arm Motion

We deproject the scaled $Q_{\phi}$ images to face-on views to measure the location of the spiral arms. For each angle $\theta$, which is defined as the counter-clockwise deviation from the northeast semimajor axis of the disk, we fit a Gaussian profile to its corresponding radial profile to obtain the peak location $r$ with error $\delta r$ using `scipy.optimize.curve_fit` (Virtanen et al. 2020). We obtain the $(\theta, r)$ pairs with a $1^\circ$ step, and present the measurements in Figure 2.

We constrain the morphology and quantify the angular offset between the two epochs for each arm under the two hypotheses. Under different motion mechanisms, a $(\theta, r)$ pair in the first epoch will advance to $(\theta + \Delta\theta, r)$ in the second epoch, where $\Delta\theta$ is the angular offset between the epochs. On one hand, in the GI scenario, each part of the arm moves

roughly on a circular orbit at the local Keplerian velocity on a short timescale, $\Delta\theta \propto r^{-3/2}$, and the arms wind up with time as the local pattern speed decreases with increasing stellocentric distance (see, e.g., Pfalzner 2003 for the winding up of spiral arms). On the other hand, in the companion scenario, an entire arm corotates around the star as a rigid body with its driver, and the angular offset between epochs is radius independent. We fit *p*-degree polynomials to the $(\theta, r)$ pairs in both epochs with predicted locations to simultaneously constrain arm morphology and obtain the motion between different epochs; see the Appendix for the mathematical formalism.

In the GI-induced scenario, if the two arms are undergoing rotation at the local Keplerian speed, the fitted pattern speed is $(0^{\circ}058 \pm 0^{\circ}009) \times (0^{\prime\prime}5/r)^{3/2} \text{ yr}^{-1}$. To take into account of the $0^{\circ}08$ true north uncertainty of SPHERE (Maire et al. 2016), which affects the position angle measurement toward the same direction within each epoch, we first propagate the $0^{\circ}009 \text{ yr}^{-1}$ measurement uncertainty using the temporal separation between the two epochs, then we combine it with the instrumental true north uncertainty for two observations assuming no correlated noise. Finally, we obtain an uncertainty of $\sqrt{(0.009 \times 4.71)^2 + 2 \times 0.08^2}/4.71 = 0^{\circ}03 \text{ yr}^{-1}$, thus the updated motion rate is $(0^{\circ}06 \pm 0^{\circ}03) \times (0^{\prime\prime}5/r)^{3/2} \text{ yr}^{-1}$. This rate corresponds to a central star mass of $0.014_{-0.010}^{+0.018} M_{\odot}$, two orders of magnitude smaller than the current estimate of $1.56_{-0.08}^{+0.11} M_{\odot}$ or $1.9 \pm 0.2 M_{\odot}$ (Garufi et al. 2018; Vioque et al. 2018). We thus rule out the GI origin of the spirals at $>5\sigma$ levels. Furthermore, the symmetric two-arm morphology in a GI disk in scattered light suggests a disk-to-star mass ratio of $\gtrsim 0.25$. (Dong et al. 2015a), which corresponds to a high accretion rate. Therefore, the disk would have been dissipated given the age of MWC 758. For illustration, we use the constrained morphology of 2015 spiral arms to predict their locations in 2019, see Figure 2.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Dynamically constrained, assumed circular orbit of a single arm driver. For a 1.9 $M_\odot$ central star, the best-fit $\pm 1\sigma$ orbit has a radius of $1''07_{-0''09}^{+0''11}$ ($172_{-14}^{+18}$ au).</figcaption></figure>

In the companion-driven scenario, the prominence and the symmetry of the two arms in the MWC 758 system suggest that they are produced by one companion of at least a few Jupiter masses (Fung & Dong 2015; Dong & Fung 2017). We thus fit the same pattern speed to both arms. The two spirals can be well fit by rigid body rotation at a rate of $0.216 \pm 0.016 \text{ yr}^{-1}$. Taking into account of SPHERE's true north uncertainty, the updated motion rate is

$$ \omega = 0^{\circ}22 \pm 0^{\circ}03 \text{ yr}^{-1}. \qquad (1) $$

This pattern speed points to a driver located at $172_{-14}^{+18}$ au, or $1''07_{-0''09}^{+0''11}$, from the 1.9 $M_\odot$ central star (Figure 3).

Our best-fit measurement of the companion-driven spiral pattern speed is consistent with the R18 measurement within $3\sigma$,<fnref n="18" /> while our derived uncertainty is ~40 times smaller, thanks to the use of the same instrument and the $Q_\phi$ maps that are the least biased by postprocessing methods. In numerical simulations, an approximately five Jupiter mass arm driver located at ~0.9 has been proposed by Baruteau et al. (2019), which is within $2\sigma$ from our best-fit companion location assuming a 1.9 $M_\odot$ central star, or $1\sigma$ assuming a 1.56 $M_\odot$ central star.

# 4. Discussion

## 4.1. Direct Imaging Constraints

We obtain the direct imaging constraints on the mass of the putative planet orbiting MWC 758 with hot-start evolutionary models (i.e., Sonora, Bobcat; M. Marley et al. 2020, in preparation) using 6.73 h of Keck/NIRC2 L'-band archival observations: 2025 s on 2015 October 24 (Program ID: C220N2, PI: E. Serabyn, Reggiani et al. 2018), 11,160 s on 2016 February 12 (Program ID: U131N2, PI: E. Chiang), 3200 s on 2016 October 24 (Program ID: C221N2, PI: G. Ruane, Reggiani et al. 2018), and 7830 s on 2017

February 2 (Program ID: U072, PI: E. Chiang). Following the method described in Ruffio et al. (2018) and taking into account of the orbital period uncertainty for the driver, we obtain an apparent L-band magnitude of 18 at 99.9% confidence level (i.e., a $3\sigma$ equivalent upper limit). Adopting an age of 10.9 Myr (Garufi et al. 2018), this corresponds to a mass of $5 M_{\text{Jupiter}}$.

Using the contrast curve of Keck/NIRC2 in the $M_\alpha$-band (Mawet et al. 2019) and assuming Gaussian noise, we expect that a $5 M_{\text{Jupiter}}$ planetary driver, whose mass has been predicted by Dong et al. (2015b) and Baruteau et al. (2019) with the semimajor axis updated in this study, can be detected at $5\sigma$ level if it is a hot-start planet using four half-nights of NIRC2 $M_\alpha$-band high contrast imaging observations (or at $3.5\sigma$ with the Spiegel & Burrows 2012 cold-start model using the same observation).

## 4.2. Eccentric Driver

If the driver has a non-zero eccentricity $e$, its semimajor axis cannot be uniquely determined from an assessment of the instantaneous angular frequency due to the unknown orbital phase. Giant planets of several Jupiter masses interacting with a gaseous disk are expected to have their eccentricities quickly damped to below the disk aspect ratio, here about 20% (Dunhill et al. 2013; Duffell & Chiang 2015). Simulations have also shown that density waves excited by planets with $e \gtrsim 0.2$ develop wiggles and bifurcations, as the waves launched at different phases interact (Li et al. 2019; Muley et al. 2019), which provide poor fits to the arms around MWC 758. A modest eccentricity introduces an uncertainty in the inferred planet location comparable to the uncertainty from pattern speed measurements—for example, $e = 0.2$ translates to a range of possible companion locations from 155 to 190 au.

## 4.3. Diverse Motion

Noticing that the two arms could be excited by different companions or by different mechanisms (e.g., Forgan et al. 2018), we investigate their motion separately.

In the GI-induced scenario, the northern arm alone rotates at an angular speed of $0.105 \pm 0.013 \times (0''5/r)^{3/2} \text{ yr}^{-1}$, and the southern arm at $0.025 \pm 0.011 \times (0''5/r)^{3/2} \text{ yr}^{-1}$. A Keplerian disk around a 1.56 $M_\odot$ central star (Vioque et al. 2018) would be rotating at an angular speed of $0.63 \times (0''5/r)^{3/2} \text{ yr}^{-1}$, which is a factor of >5 faster than the measurement and inconsistent with observations for both arms at >5$\sigma$ levels. Therefore, an even faster rotating Keplerian disk around a 1.9 $M_\odot$ central star (Garufi et al. 2018) is inconsistent with the motion rates.

In the companion-driven scenario, if the arms are driven by different companions, we measure that the northern arm rotates by $0.211 \pm 0.019 \text{ yr}^{-1}$, and the southern arm rotates by $0.228 \pm 0.027 \text{ yr}^{-1}$. The two rates are within $1\sigma$ from each other, consistent with the expectation that the two arm are corotating and driven by the same driver.

## 4.4. Model Selection

From a statistical approach, our fitting results have a $\chi^2$ value of 3772 in the single-companion-induced scenario, and 3812 in the global GI-induced scenario. Given that the two mechanisms are applied to the same number of data points (i.e., location pairs) and have the same number of variables (i.e., 1 rotation speed variable and $2p + 2 = 8$ polynomial coefficient

<footnote n="18">¹⁸ Only $3\sigma$ uncertainties are well-constrained in R18.</footnote>
variables), the Schwarz information criterion (SIC; Schwarz 1978) difference is then $\Delta\text{SIC} = \Delta\chi^2 = 40$, which is greater than $\Delta\text{SIC} = 10$ threshold for “decisive” evidence (Kass & Raftery 1995) for model selection, making a single-planet driver the preferred mechanism.

From another approach, assuming identical arm morphology between the two epochs, we can focus on the marginalized distribution for the speed parameter $\omega$ to quantify the difference. In this way, the above $\Delta\chi^2 = 40$ difference corresponds to a confidence level of $\sqrt{40}\sigma = 6.3\sigma$ (Chapter 15.6 of Press et al. 1992), which makes the single-planet driver mechanism more consistent with our observations. Similarly, we apply the above analysis to each individual arm and find that the planet-driven scenario is more consistent with the observations.

### 4.5. Robustness Estimation

In our analysis, we have investigated polynomials up to $p=7$ degrees. We obtain the lowest SIC that penalizes excessive use of parameters at the cubic form when $p=3$. When $p \ge 3$, we observe no discernible best-fit angular speeds, and thus we use the cubic description of the spiral arms in our analysis. To robustly obtain the best fit and uncertainty for these parameters, we have investigated the impact from different chi-squared minimization methods, including `scipy CURVE Fit` (Virtanen et al. 2020) and the orthogonal least-squares-fitting code `scipy.odr` (Boggs et al. 1989), and no discernible difference was obtained. We present in this study the results from `scipyড` optimize ``curve_fit``.

The flaring of the disk (e.g., Stolker et al. 2016; Rosotti et al. 2020) does not bias our estimation. We use `diskmap` (Stolker et al. 2016) to deproject the disk images with various flaring exponents (i.e., 0, 0.5, 1, 1.1, 1.2) and repeat the motion measurement; the results are all consistent within $1\sigma$. In addition, we randomly varied the inclination and position angle for the disk within $\pm 5^\circ$ $10^3$ times, and the motion rates for the single-planet driver are $0.22 \pm 0.06 \text{ yr}^{-1}$ and consistent with the original estimate within $2\sigma$.

Our measurements are not biased by star centering uncertainties in two aspects. First, we use the same pipeline (i.e., IRDAP) with identical reduction parameters to minimize systematic offset. Second, even if there were offsets, the inclination of the disk would impact in the individual arm rotation rates by returning different angular speed measurements in the planet-driven scenario, which are indistinguishable since our measurements are within $1\sigma$.

The impact from individual location pairs is negligible. We experiment by randomly discarding up to 25% of the pairs and repeating the speed measurement procedure $10^4$ times. The best-fit rotation rate for a single driver is found to be $0.22 \pm 0.02 \text{ yr}^{-1}$, consistent with our initial measurement to within $1\sigma$.

The morphology of the spiral arms is consistent with being circular when the stellocentric separation is less than 40 au ($0.25$, Figure 2), and such regions have a seemingly outward motion in Figure 2. We argue here that this does not bias our results. When we ignore these location pairs and repeat our fitting, the results do not change by more than $1\sigma$. In fact, since the number of data points in these regions is less than 25% of the total number of data points, this scenario has been investigated in the above procedure of random location pair rejection.

### 4.6. Possible Systematics

The 2015 March data are taken in field-tracking mode with an non-ideal HWP control law,<fnref n="19" /> which was rectified in 2015 late April. The non-ideal control law causes the polarization direction to rotate on the detector during the observations, which we correct for using the Mueller matrix model of IRDAP. We confirm the proper correction of the images with IRDAP by comparing the uncorrected and corrected polari-metric images cubes for both the 2015 and 2019 data.

The 2015 data are taken with a sub-optimal order of HWP switch angles ($0^\circ, 22.5^\circ, 45^\circ, 67.5^\circ$ instead of $0^\circ, 45^\circ, 22.5^\circ, 67.5^\circ$ for the 2019 data) and with a high number of integrations per HWP position (six instead of two for the 2019 data). As a result, for the 2015 data a measurement of $Q$ or $U$ lasts approximately 10 min 10 s, compared to 2 min 25 s for the 2019 data. Because the polarization direction rotates on the detector during the 2015 measurements, there may be a global impact on the final $Q_\phi$ image. The 2019 $Q_\phi$ image may also be slightly affected, because in the reduction of pupil-tracking data the images are derotated after computing the double difference (van Holstein et al. 2017). Given that the disk has a non-zero inclination and position angle, these effects could bias the speed measurement when we deproject the image to face-on views. Therefore, it should be reflected in the measured individual rotation speeds under the planet-driver scenario that calculates global offsets. Nevertheless, since the individual arm rotation rates are consistent within $1\sigma$, we do not expect the global impact from observation strategy to bias our results at more than the $1\sigma$ level.

There are caveats in our measurements. Finite inclinations are known to produce distortions in images that are hard to correct in deprojection (Dong et al. 2016). However, this effect usually becomes prominent only at inclinations larger than ~20°, and we do not expect strong morphology distortions in the MWC 758 disk whose inclination is ~20°. In addition, perturbations from a theorized inner companion in the disk (Baruteau et al. 2019) may cause slight changes to the shape of the spiral arms that are unrelated to their primary driver, thus affecting our pattern speed measurements. Furthermore, the change in illumination may slightly change the observed features on disk surface (Montesinos et al. 2016). Future multi-epoch observations of MWC 758 are necessary to quantify such effects.

## 5. Summary

We have established a 5 yr baseline and obtained the most accurate pattern speed measurement of spiral arms in a protoplanetary disk to date. For the two prominent spiral arms surrounding MWC 758, we witness global scale shadowing effects and measure the motion between the two epochs to test their formation and motion mechanisms.

We found that the measured motion of spirals disfavors their GI origin. This is the first time that this has been shown for any protoplanetary disk. Meanwhile, our motion analysis suggests a single planet driving both spiral arms. For a $1.9 M_\odot$ central star, our measurement pinpoints a semimajor axis of $172_{-14}^{+18}$ au for the planet driver if its orbit is circular. Using archival Keck/NIRC2 L'-band observations totaling 6.73 hr, we obtain a $3\sigma$-equivalent upper limit of $5 M_{\text{Jupiter}}$ for the location of this driver

<footnote n="19">¹⁹ See SPHERE User Manual at http://www.eso.org/sci/facilities/paranal/instruments/sphere/doc.html.</footnote>
using hot-start planet formation models. The inferred spiral arm driver in the MWC 758 system is ideal for Keck/NIRC2, Very Large Telescope/Enhanced Resolution Imager and Spectrograph, and James Webb Space Telescope direct detections in longer wavelengths, and for Atacama Large Millimeter/ submillimeter Array circumplanetary disk exploration.

We thank the anonymous referee for comments that improved the clarity of this Letter, and Cassandra Hall for useful discussions. This research is partially supported by NASA ROSES XRP, award 80NSSC19K0294. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 060.A-9389(A) and 104.C-0472(A). B.R. thanks Christian Ginski for discussions on shadowing effects, Rémi Soummer for initiating the Archival Legacy Investigations of Circumstellar Environments (ALICE) project that set up the stage for R18 and this Letter. T.E. was supported in part by NASA Grants NNX15AD95G/NEXSS, NNX15AC89G, and NSF AST-1518332. R.D. acknowledges the support from the Natural Sciences and Engineering Research Council of Canada. A.L.M. acknowledges the financial support of the F.R.S.-FNRS through a postdoctoral researcher grant. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

Facilities: Very Large Telescope (SPHERE), Keck:II (NIRC2).

Software: diskmap (Stolker et al. 2016), IRDAP (van Holstein et al. 2017, 2020), scipy (Virtanen et al. 2020).

## Appendix

### Spiral Arm Motion

To constrain spiral arm morphology and motion, we note that $\theta$ was expressed as function of $r$ to allow for matrix inversion using linear algebra in R18. After inspecting the stability of the high-resolution $Q_\phi$ images in this study, here we switch their relationship in order to allow for precise measurement of $r$ as a function of $\theta$. For a rotating spiral arm, a $(\theta_i, r_i)$ pair will be updated to $(\theta_i + \omega_{\text{model},i}, t, r_i)$ at a new epoch, where $t$ is the temporal separation between the two observations, and $\omega_{\text{model},i}$ the angular speed in the scenario that is either companion-driven (“comp”) or gravity instability-induced (“GI”). For a total of $E$ epochs, we describe the location pairs using polynomials with $(E-1)$ dummy variables,

$$r_i^{(\text{model})}(\theta_i) = \sum_{j=0}^{p} c_j \left( \theta_i + \omega_{\text{model},i} \sum_{k=2}^{E} t_k D_k \right)^j, \quad (\text{A1})$$

where $p \in \mathbb{N}$ describes the degree of the polynomial at the first epoch when $r_i(\theta_i) = \sum_{j=0}^{p} c_j \theta_i^j$ and $c_j \in \mathbb{R}$ is the coefficient for the $j$th power term, $t_k \in \mathbb{R}$ is temporal separation between epochs $k$ and 1, and $D_k \in \{0, 1\}$ is a dummy variable that equals 1 only when the $(\theta_i, r_i)$ pair is obtained at epoch $k$. We note that the above equation is to describe the location pairs

using $E$ polynomials that are mutually related through angular offsets.

For a total of $m$ measured location pairs, we minimize the following chi-squared statistic:

$$\chi^2 = \sum_{i=1}^{m} \left( \frac{r_i - r_i^{(\text{model})}(\theta_i)}{\delta r_i} \right)^2, \quad (\text{A2})$$

to obtain the motion rate. In this way, we can simultaneously constrain the morphological parameters and motion rate using all available location pairs.

#### A.1. Companion Driven

If a spiral arm is driven by a companion on a circular orbit (“comp”; Kley & Nelson 2012; Dong et al. 2015b; Zhu et al. 2015; Bae et al. 2016) that is located at a stellocentric position of $r_{\text{comp}}$, the entire arm corotates with the companion at the Keplerian angular speed of the companion, $\omega_{\text{comp}}$. In this way, an arm observed at different epochs is shifted in the azimuthal direction while maintaining its shape in the disk plane. A $(\theta_i, r_i)$ pair will be updated to $(\theta_i + \omega_{\text{comp}}t, r_i)$ at a new epoch. Equation (A1) then becomes

$$r_i^{(\text{comp})}(\theta_i) = \sum_{j=0}^{p} c_j \left( \theta_i + \omega_{\text{comp}} \sum_{k=2}^{E} t_k D_k \right)^j. \quad (\text{A3})$$

We note that the physical meaning of the above equation is to fit offset identical polynomials to the data.

In the companion-driven scenario, we substitute Equation (A3) into Equation (A2), i.e.,

$$\chi^{2(\text{comp})} = \sum_{i=1}^{m} \left[ \frac{r_i - \sum_{j=0}^{p} c_j \left( \theta_i + \omega_{\text{comp}} \sum_{k=2}^{E} t_k D_k \right)^j}{\delta r_i} \right]^2, \quad (\text{A4})$$

to obtain the motion rate when the arms co-move.

When there are a total of $s$ spiral arms, we denote their arm location measurements with $(r_{i,l}, \theta_{i,l})$ for $l \in \{1, \dots, s\}$. These spiral arms rotate at the same rate if all of them are driven by the same driver; then Equation (A3) becomes

$$r_{i,l}^{(\text{comp})}(\theta_{i,l}) = \sum_{l=1}^{s} \sum_{j=0}^{p} c_{j,l} \left( \theta_{i,l} + \omega_{\text{comp}} \sum_{k=2}^{E} t_k D_k \right)^j D_l, \quad (\text{A5})$$

where $D_l \in \{0, 1\}$ is a dummy variable that equals 1 only when the $(\theta_{i,l}, r_{i,l})$ pair is obtained from spiral arm $l$. The corresponding $\chi^2$ minimization formula is obtained by substituting the $r$ expressions in Equation (A5) to the $\chi^2$ expression in Equation (A2). The physical meaning of the above equation is to fit multiple arms using Equation (A3) but with the constraint that their motion rates are identical.

#### A.2. GI Induced

If a spiral arm is excited by GI (Lodato & Rice 2005; Dong et al. 2015a; Kratter & Lodato 2016), each part of the arm rotates at its local Keplerian angular speeds in the disk plane on a timescale much smaller than the local dynamical timescale (spiral arms disappear and reemerge on a longer timescale). For
any $(\theta_i, r_i)$ location pair, its location at a new epoch will be $(\theta_i + \frac{r_0^{3/2}}{r_i^{3/2}}\omega_0 t, r_i)$, where $\omega_0$ is the Keplerian angular speed at stellocentric separation $r_0$. In this scenario, Equations (A1) and (A2) have a power-law attenuation in their angular speed terms, i.e.,

$$r_i^{(\text{GI})}(\theta_i) = \sum_{j=0}^{p} c_j \left( \theta_i + \omega_0 \frac{r_0^{3/2}}{r_i^{3/2}} \sum_{k=2}^{E} t_k D_k \right)^j, \quad (\text{A6})$$

and

$$\chi^{2(\text{GI})} = \sum_{i=1}^{m} \left[ \frac{r_i - \sum_{j=0}^{p} c_j \left( \theta_i + \omega_0 \frac{r_0^{3/2}}{r_i^{3/2}} \sum_{k=2}^{E} t_k D_k \right)^j}{\delta r_i} \right]^2, \quad (\text{A7})$$

respectively. We note that the physical meaning of the above two equations is to fit lines that have $r^{-3/2}$-dependent angular offsets.

Similarly, for a total of *s* spiral arms undergoing the same local Keplerian motion, the corresponding power law attenuation form of Equation (A5) is

$$r_{i,l}^{(\text{GI})}(\theta_{i,l}) = \sum_{l=1}^{s} \sum_{j=0}^{p} c_{j,l} \left( \theta_{i,l} + \omega_0 \frac{r_0^{3/2}}{r_i^{3/2}} \sum_{k=2}^{E} t_k D_k \right)^j D_l. \quad (\text{A8})$$

The corresponding $\chi^2$ minimization formula is obtained by substituting the *r* expression in Equation (A8) to the $\chi^2$ expression in Equation (A2). The physical meaning of the above equation is similar to Equation (A5) but for the GI-induced mechanism.

In this study, we have $E = 2$ epochs with a temporal separation of $t = 4.71$ yr. We report the motion with polynomial degree $p = 3$ for both scenarios.

## ORCID iDs

Bin Ren (任彬) https://orcid.org/0000-0003-1698-9696

Ruobing Dong (董若冰) https://orcid.org/0000-0001-9290-7846

Rob G. van Holstein https://orcid.org/0000-0003-1520-8405

Jean-Baptiste Ruffio https://orcid.org/0000-0003-2233-4821

Julien H. Girard https://orcid.org/0000-0001-8627-0404

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Thomas M. Esposito https://orcid.org/0000-0002-0792-3719

Élodie Choquet https://orcid.org/0000-0002-9173-0740

Dimitri Mawet https://orcid.org/0000-0002-8895-4735

Tomas Stolker https://orcid.org/0000-0002-5823-3072

Eugene Chiang https://orcid.org/0000-0002-6246-2310

John H. Debes https://orcid.org/0000-0002-1783-8817

Dean C. Hines https://orcid.org/0000-0003-4653-6161

Anne-Lise Maire https://orcid.org/0000-0002-2591-4138

François Ménard https://orcid.org/0000-0002-1637-7393

Maxwell A. Millar-Blancher https://orcid.org/0000-0001-6205-9233

Marshall D. Perrin https://orcid.org/0000-0002-3191-8151

Charles A. Poteet https://orcid.org/0000-0003-4845-7483

Glenn Schneider https://orcid.org/0000-0002-4511-5966

## References

Bae, J., Zhu, Z., & Hartmann, L. 2016, ApJ, 819, 134

Baruteau, C., Barraza, M., Pérez, S., et al. 2019, MNRAS, 486, 304

Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Boehler, Y., Ricci, L., Weaver, E., et al. 2018, ApJ, 853, 162

Boggs, P. T., Donaldson, J. R., Byrd, R. H., & Schnabel, R. B. 1989, TOMS, 15, 348

Brittain, S. D., Najita, J. R., Dong, R., & Zhu, Z. 2020, ApJ, 895, 48

Calcino, J., Christiaens, V., Price, D. J., et al. 2020, arXiv:2007.06155

Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205

Dong, R., & Fung, J. 2017, ApJ, 835, 38

Dong, R., Fung, J., & Chiang, E. 2016, ApJ, 826, 75

Dong, R., Hall, C., Rice, K., & Chiang, E. 2015a, ApJL, 812, L32

Dong, R., Najita, J. R., & Brittain, S. 2018, ApJ, 862, 103

Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015b, ApJL, 809, L5

Duffell, P. C., & Chiang, E. 2015, ApJ, 812, 94

Dunhill, A. C., Alexander, R. D., & Armitage, P. J. 2013, MNRAS, 428, 3072

Forgan, D. H., Ilee, J. D., & Meru, F. 2018, ApJL, 860, L5

Fung, J., & Dong, R. 2015, ApJL, 815, L21

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Garufi, A., Avenhaus, H., Pérez, S., et al. 2020, A&A, 633, A82

Garufi, A., Benisty, M., Pinilla, P., et al. 2018, A&A, 620, A94

Grady, C. A., Muto, T., Hashimoto, J., et al. 2013, ApJ, 762, 48

Grady, C. A., Woodgate, B., Bruhweiler, F. C., et al. 1999, ApJL, 523, L151

Isella, A., Natta, A., Wilner, D., et al. 2010, ApJ, 725, 1735

Kass, R. E., & Raftery, A. E. 1995, JASA, 90, 773

Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211

Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

Li, Y.-P., Li, H., Li, S., & Lin, D. N. C. 2019, ApJ, 886, 62

Lodato, G., & Rice, W. K. M. 2005, MNRAS, 358, 1489

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, Proc. SPIE, 9908, 990834

Mawet, D., Hirsch, L., Lee, E. J., et al. 2019, AJ, 157, 33

Ménard, F., Cuello, N., Ginski, C., et al. 2020, A&A, 639, L1

Monnier, J. D., Harries, T. J., Bae, J., et al. 2019, ApJ, 872, 122

Montesinos, M., Perez, S., Casassus, S., et al. 2016, ApJL, 823, L8

Muley, D., Fung, J., & van der Marel, N. 2019, ApJL, 879, L2

Muro-Arena, G. A., Ginski, C., Dominik, C., et al. 2020, A&A, 636, L4

Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJL, 748, L22

Nealon, R., Pinte, C., Alexander, R., et al. 2019, MNRAS, 484, 4951

Pfalzner, S. 2003, ApJ, 592, 986

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C. The Art of Scientific Computing (2nd ed.; Cambridge: Cambridge Univ. Press)

Reggiani, M., Christiaens, V., Absil, O., et al. 2018, A&A, 611, A74

Ren, B., Dong, R., Esposito, T. M., et al. 2018, ApJL, 857, L9

Rosotti, G. P., Benisty, M., Juhász, A., et al. 2020, MNRAS, 491, 1335

Ruffio, J.-B., Mawet, D., Czekala, I., et al. 2018, AJ, 156, 196

Schwarz, G. 1978, AnSta, 6, 461

Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174

Stolker, T., Dominik, C., Min, M., et al. 2016, A&A, 596, A70

Stolker, T., Sitko, M., Lazareff, B., et al. 2017, ApJ, 849, 143

van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020, A&A, 633, A64

van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, Proc. SPIE, 10400, 1040015

Vieira, S. L. A., Corradi, W. J. B., Alencar, S. H. P., et al. 2003, AJ, 126, 2971

Vioque, M., Oudmaijer, R. D., Baines, D., et al. 2018, A&A, 620, A128

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJL, 813, L2

Wagner, K., Stone, J. M., Spalding, E., et al. 2019, ApJ, 882, 20

Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88