A multiwavelength analysis of the spiral arms in the protoplanetary disk around WaOph 6

[Full abstract in the paper] In recent years, protoplanetary disks with spiral structures have been detected in scattered light, millimeter continuum, and CO gas emission. The mechanisms causing these structures are still under debate. A popular scenario to drive the spiral arms is the one of a planet perturbing the material in the disk. However, if the disk is massive, gravitational instability is usually the favored explanation. Multiwavelength studies could be helpful to distinguish between the two scenarios. So far, only a handful of disks with spiral arms have been observed in both scattered light and millimeter continuum. We aim to perform an in-depth characterization of the protoplanetary disk morphology around WaOph 6 analyzing data obtained at different wavelengths, as well as to investigate the origin of the spiral features in the disk. We present the first near-infrared polarimetric observations of WaOph 6 obtained with SPHERE at the VLT and compare them to archival millimeter continuum ALMA observations. We traced the spiral features in both data sets and estimated the respective pitch angles. We discuss the different scenarios that can give rise to the spiral arms in WaOph 6. We tested the planetary perturber hypothesis by performing hydrodynamical and radiative transfer simulations to compare them with scattered light and millimeter continuum observations.


Introduction
The study of protoplanetary disks provides insight into our understanding of the formation and early evolution of planets. Modern observational techniques have enabled significant progress to be made in this task. Recent observations in both scattered light and millimeter continuum have shown the striking frequency with which these disks present structures, such as gaps, rings, or spirals (e.g., Avenhaus et al. 2018;Long et al. 2018;Andrews 2020;Cieza et al. 2020).
In particular, spiral arms have been observed in more than a dozen young disks, spanning from tens to hundreds of au. They have been found in scattered light (e.g., Muto et al. 2012;Wagner et al. 2015;Benisty et al. 2015;Stolker et al. 2016;Garufi et al. 2020;Muro-Arena et al. 2020), millimeter continuum (e.g., Pérez et al. 2016;Huang et al. 2018b;Rosotti et al. 2020), and more recently in CO gas emission (e.g., Tang et al. 2017;Kurtovic et al. 2018). However, for the disks that have been imaged at different wavelengths, some of these spirals are only present Member of the International Max-Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD), Germany in either scattered light, and not in the millimeter continuum, or vice versa. The mechanisms that drive spiral arms in these disks are probably manifold, and they are a matter of debate in many individual cases.
Considering that planets are born in protoplanetary disks, the observed structures have been frequently linked with the presence of planets forming within the disk (e.g., Muto et al. 2012;Pohl et al. 2015;Dong et al. 2018a;Calcino et al. 2020;Ren et al. 2020). These planets perturb the disks via gravitational interactions, and these perturbations can cause the formation of spirals. In these cases, the planets cause spiral arms both in the interior and exterior of their own orbit. In particular, when the planet is fairly massive, it can trigger a secondary spiral arm on the inside of its orbit. For a planet of tens of M Jup , the primary and secondary spiral arm reach approximately a m = 2 1 . From such planet-driven spirals, we can study the mass and location of the potential planets. Another mechanism that can drive spiral arms is gravitational instability (GI, e.g., Goldreich & Lynden-Bell 1965). This happens when the self-gravity perturbations in the disk dominate over the restoring forces of gas pressure and differential rotation (Toomre 1964), which can be translated to the disk being gas rich and relatively massive with respect to its host star (M disk /M * 0.1 Kratter & Lodato (2016)). Spiral arms generated by GI can allow one to better constrain the disk mass. A third possible cause for the presence of spiral structure is a flyby event of a close companion to the system (e.g., Cuello et al. 2019;Ménard et al. 2020) that perturbs the disk material causing over densities that form the spiral pattern. And finally, a combination of some of the physical processes described above has also been investigated (e.g., Pohl et al. 2015).
Distinguishing between the possible physical processes driving spiral structure in protoplanetary disks is not an easy task. In fact, measuring the disk gas mass alone is already challenging. Having a more comprehensive understanding of these disks can be useful to decipher between the physical processes that are taking place within them. One approach is to use multiwavelength observations to trace their different components. In particular, images obtained in scattered light are sensitive to micron-sized dust grains at the disk surface, which under typical disk conditions are very well coupled to the gas. Images at (sub)millimeter wavelengths, on the other hand, trace larger grains within the disk midplane, and they are less well coupled to the gas. Therefore, comparing images at these wavelengths with similar angular resolution can potentially reveal the different morphologies of different disk components. Actually, it may be possible to distinguish between planet-or GI-induced spirals by comparing scattered light and millimeter continuum observations, since dust trapping in spiral arms is likely to be more efficient in gravitationally unstable disks (Juhász et al. 2015;Dipierro et al. 2015). Previous studies have utilized observations in both near-infrared (NIR) and millimeter continuum to analyze the disks' morphology, as well as the gas and dust distributions within them (e.g., van Boekel et al. 2017;Rosotti et al. 2020), however, so far there are only a handful of such works.
In this paper, we present NIR polarimetric observations of WaOph 6 obtained using the VLT/SPHERE instrument (Beuzit et al. 2019). These observations have a spatial resolution of ∼51 mas (∼ 6 au), similar to that of the recent ALMA-DSHARP millimeter wavelength observations (Andrews et al. 2018, see Fig. 1). We report a NIR scattered light counterpart of the innermost spiral arms reported in Huang et al. (2018b). Additionally, we explore one of the scenarios that can give rise to the spiral pattern observed in the disk.
The paper is structured as follows: in Section 2 we introduce WaOph 6 and the previous studies on its disk. In Section 3 we describe our scattered light observations and data reduction procedure. Our results are presented in Section 4. We describe the modeling setup and the comparison between the simulations and the observations on Section 5. In Section 6 we discuss our results, and finally in Section 7 we summarize our findings and list our conclusions.

Stellar and disk properties
WaOph6 is a K6 star (Eisner et al. 2005), and member of the Ophiuchus moving group at a distance of 122.5 +0.3 −0.2 pc (Gaia Collaboration et al. 2020) located near the L162 dark cloud. It was first identified as a suspected T Tauri star by Henize (1976), and then confirmed by Walter (1986). Here we constrain the stellar mass and age based on the updated photometry and Gaia parallax. We retrieved the full spectral energy distribution (SED) from Vizier 2 and employed a Phoenix model of the stellar photosphere (Hauschildt et al. 1999) with effective temperature T eff = 4200 K (Eisner et al. 2005), surface gravity log(g)=-4.0, and an optical extinction A V = 2.8 ± 0.3 mag calculated from the V, R, and I photometric fluxes. We integrated the stellar model scaled to the average V magnitude and Gaia distance of 122.5 pc obtaining a stellar luminosity of L * = 1.91 +0.70 −0.51 L . Then, we placed the source on the HR diagram and constrain a stellar mass M * = 0.7±0.1 M and an age t = 0.6±0.3 Myr through different sets of PMS tracks (Parsec, MIST, Bressan et al. 2012;Choi et al. 2016;Baraffe et al. 2015) with error bars propagated from L * and T eff (± 100 K).
The disk around WaOph 6 has been a common target for millimeter continuum surveys looking to constrain and characterize the structure in protoplanetary disks (e.g., Andrews & Williams 2007;Andrews et al. 2009;Ricci et al. 2010;Andrews et al. 2018). Submillimeter Array (SMA) observations were used along with a parametric model to constrain density structure parameters (Andrews et al. 2009). The disk model that best fitted the thermal continuum data and spectral energy distribution (SED) was that of a flat cold disk with a total disk mass (gas + dust-to-gas ratio of 1:100) of 0.077 M . With observations from the Australia Telescope Compact Array (ATCA), Ricci et al. (2010) analyzed and modeled the SED of WaOph 6, adopting a distance of ∼130 pc and an outer radius (R out ) interval of 175 − 375 au, and they find dust mass estimates (M dust ) between 8 × 10 −5 M and 9.8 × 10 −5 M , depending on the assumed dust size distribution power-law index (q = 2.5 or q = 3). More recently, WaOph 6 was observed by ALMA within the DSHARP program (Disk Substructures at High Angular Resolution Project, Andrews et al. 2018). These millimeter continuum observations showed that the disk has a set of symmetric spiral arms that extend to ∼70 au, a gap at 79 au and a bright ring at 88 au (Huang et al. 2018b). The disk has an inclination (i) of 47.3 • and a position angle (PA) of 174.2 • obtained from ellipse fitting on the dust continuum emission (see Huang et al. 2018a, for more details), and gas observations have shown that it suffers from mild molecular cloud contamination (Reboussin et al. 2015). We summarize the stellar and disk physical parameters in Table 1, where we include the different values for the disk mass found in the literature, as well as our own M dust estimate obtained following the procedure described in Appendix A.
Up to now, only seven disks have been known to have spiral arms in millimeter continuum wavelengths: WaOph 6, Elias 27, IM Lup, HT Lup A, AS 205 N, MWC 758, and HD100453 (Huang et al. 2018b;Kurtovic et al. 2018;Dong et al. 2018a;Rosotti et al. 2020), and only the first three are single systems. Out of these three, only IM Lup has published polarized scattered light observations (Avenhaus et al. 2018), however, with no spiral arms visible at these wavelengths.

Observations and data reduction
WaOph 6 was observed with the VLT/SPHERE high-contrast instrument (Beuzit et al. 2019) Table 2). The observations were carried out with the IRDIS Dual-beam Polarimetric Imaging (DPI) mode (Langlois et al. 2014;de Boer et al. 2020;van Holstein et al. 2020) in H-band (λ c =1.625 µm; ∆λ=0.291 µm; where λ c denotes the  (Huang et al. 2018b). The beam size is shown in the lower left corner. Assuming a gas-to-dust ratio of 100:1 and based on SMA SED modeling. 2 Obtained with a power-law index for the grain size distribution q = 2.5.   Langlois et al. 2013, for more details on this mode). Two flux calibration frames (images of the target star without the coronagraph) were obtained with an exposure time of 2 s and a neutral density (ND1) filter to prevent saturation. We measure a point spread function (PSF) FWHM of ∼51 mas by fitting a Gaussian function to the flux frames. The weather conditions were stable during the observations with a seeing of ∼0.5", a coherence time (τ 0 ) of ∼4 ms, and wind speed of ∼3.8 m/s. The Strehl ratio was about 0.7, however, the low scattered light intensity resulted in a rather low signal-to-noise ratio (S/N).
For the data reduction, we used the IRDAP pipeline 3 version 1. 3.2. (van Holstein et al. 2020). First, the pipeline preprocesses the data by performing the usual sky background subtraction, flat fielding, bad-pixel identification and interpolation, and star centering corrections. Subsequently, polarimetric differential imaging (PDI) is performed by applying the double-sum and doubledifference method described in de Boer et al. (2020) to obtain a set of Stokes Q and U frames. Finally, the data are corrected for instrumental polarization and crosstalk effects by applying a detailed Mueller matrix model of the instrument (see van Holstein et al. 2020, for more details on the data reduction procedure), yielding the final Q and U images. The final PDI images are corrected for true north following the procedure established by Maire et al. (2016). IRDAP then obtains the linearly polarized intensity (PI) image using the final Q and U images, from (1) Next, the pipeline computes the azimuthal Stokes parameters following : where φ is the position angle (PA) measured east of north with respect to the position of the star. In the definition above, a positive signal in the Q φ image corresponds to a signal that is linearly polarized in the azimuthal direction, while a negative signal denotes radially polarized light in Q φ . U φ contains any signal polarized at ±45 • with respect to the radial direction. This means that for disks with low inclinations (i.e., i < 45 • ), almost all of the scattered light is expected to be included as a positive signal in Q φ , while the U φ image can be considered as an upper limit of the noise level. In the case of WaOph 6, we expect some physical signal in the U φ image due to the inclination of the disk (see Table 1). The resulting Q φ and U φ images are shown in Appendix B.

SPHERE/IRDIS-DPI observations
In Fig. 2 we show the final, processed Q φ and U φ images. Due to the low S/N, we sharpened the images by subtracting a version of them which was convolved by a Gaussian kernel with the size of 10 pixels, which removes low frequency structures, and then we convolved them with a Gaussian kernel of size 0.1 × FWHM to smooth the images in order to enhance the spiral features. We observe the launch of the spiral arms up to ∼0.3" (40 au), as seen in the Q φ image (Fig. 2, left). As mentioned in Section 3, the U φ image (Fig. 2, right) in this case contains almost no signal and can be used as an upper limit of the noise level. For a better visualization of the spiral features, we plotted the azimuthal profile by first deprojecting the filtered Q φ image, and taking the average flux within the ring between ∼27 and 36 au in azimuthal bins of 15 • . The distance range is chosen due to the presence of the coronagraph at lower radii, and the S/N decrease at higher radii. In Fig. 3, we show the smoothed azimuthal profile. The spiral arms are seen as the two peaks between 20 • − 100 • and 200 • − 310 • . To estimate the spiral arm intensity contrast between the spiral and inter-spiral regions, we measured the peak intensities of the spiral arms in radial bins spaced by 2 au. We estimated the inter-spiral region intensity by taking the minimum value in each bin. The contrast is then the ratio between the peak intensities and the inter-spiral intensities. On average, we find the spiral arm contrast to be 1.5.

Companion candidate
We detect a companion candidate (CC) in our data as shown in the total intensity image in the top panel of Fig. 4, where the CC is more visible than in the polarized light frames. The CC is located at a projected distance of ∼400 au (∼3") and has a brightness contrast of 10 −3 with respect to WaOph 6. We used archival HST data (from 1999-01-23) as an additional epoch to perform an astrometric analysis in order to verify if the CC is bound to the system. The resulting astrometry plot is shown in the bottom panel of Fig. 4, where the black curve traces the path a stationary background object would have followed relative to WaOph 6 between the two epochs, and the markers show the position of the CC at both the HST and the SPHERE epochs. Since the CC is located near the final position a background object would be located at, we conclude that the object is not bound to the system, and therefore could not be considered as an external perturber causing the spirals. We note that the CC is not reported in the Gaia EDR3, despite of its presence in the two data sets described above.

Comparison to ALMA Observations
WaOph 6 was observed within ALMA/DSHARP (Andrews et al. 2018) in Band 6, at a frequency of 239 GHz (1.3 mm). Observations at these wavelengths sample the millimeter-sized dust grains that are typically located in the disk midplane (see e.g., Villenave et al. 2020). On the other hand, our SPHERE observations trace the light scattered from submicron sized dust grains located at the disk surface which are typically well-coupled to the gas. In this section we do a first comparison of the two data sets.

15.26
Note: For the SPHERE data, N corresponds to the northern spiral and S to the southern spiral. In the case of the ALMA data, N1 corresponds to the northern inner spiral, S1 to the southern inner spiral, N2 to the northern outer spiral and S2 to the southern outer spiral. For the Archimedean fit, the pitch angles are estimated at 35 au.

Radial profiles
In Fig. 5 we plot the radial intensity profiles of the SPHERE/IRDIS-DPI H-band image (teal curve), and the ALMA image in Fig. 1 (crimson curve). The curves have been normalized to the maximum intensity value on each image for visualization purposes. In order to obtain a smooth profile, we apply a Gaussian kernel of size 0.1 × FWHM to the reduced Q φ image in Appendix B. We obtained the profiles by taking the azimuthal average of the image intensity in rings of radius 3 au. For this we considered the latest literature values for i and PA (listed on Table 1), and we use the aperture_photometry() function from the photutils python package, which allows to perform aperture photometry within elliptical annuli. This permits to get the radial profile without first deprojecting the image. As expected from the fact that the two images trace different dust sizes, the profiles do not perfectly overlap. A closer look between 70 au and 90 au shows that the substructures created by the gap and the ring described in Section 2, are present in both profiles, with a slight shift. The initial drop in the intensity of the SPHERE/IRDIS-DPI profile is due to the use of a coronagraph in these observations. The cut in the ALMA data profile can be attributed to the emission from the large grains being limited to the central ∼130 au of the disk.

Spiral arms
Next, we performed a spiral search on each data set to compare the location of the spiral arms at both wavelengths. In the case of the SPHERE data, the spiral features are better seen when the image is plotted in polar coordinates, where spiral arms appear as inclined lines. To obtain the polar plot we first deprojected the image and then converted to polar coordinates. We then used the python function peak_local_max from the skimage package to search for the peak emission points around the location of the spiral arms. To trace the spirals in the ALMA data we used two different images generated using the tclean task in CASA 5.4.1 (McMullin et al. 2007). We used uv-taper = ['0.010arcsec', '0.010arcsec', '10deg'], and two different robust values (-1.0 and 1.0) to generate an image with higher resolution in the central and in the outer regions, respectively. Then, we searched for peak emission points along the radial direction and took the 3σ emission ones as the spiral arms. We repeated this process  in both of the images described above. The resulting spirals are shown in Fig. 6, where we overplot the retrieved spiral arms on both the SPHERE (left) and the ALMA image (right). To obtain the SPHERE image we followed the same procedure described in Section 4 with a Gaussian kernel of size 0.2 × FWHM. We did this to match the scaling of the ALMA image for better comparison purposes. To generate the ALMA image, we used the frank (Jennings et al. 2020) tool to remove the azimuthally component of the emission of the disk and leave only the nonaxisymmetric features, as described in the Appendix C. We note that the spiral pattern is much more prominent in the mm continuum than in scattered light. Besides the low S/N of the SPHERE data and the fact that the disk has been reported to be cold, this could also be explained by the anisotropic scattering properties of the dust of different sizes. For a disk that is not edge-on, the larger the particles on the surface layer compared to the wavelength, the more they will scatter light into the disk in forward scattering. Therefore, the amount of light scattered in the line of sight direction would be smaller and would not be detected in scattered light images (e.g., Mulders et al. 2013). Additionally, we notice that there appears to be a break in the spiral arms on the ALMA image at ∼0.16" (∼20 au), which cannot be observed in our SPHERE data. In the following we treat this break as a separate set of spirals.
In order to characterize the spirals we considered two models. A logarithmic spiral given by: and an Archimedean spiral, defined as: where θ is the polar angle, r 0 is the radius for which θ=0, and b relates to the pitch angle (µ) of the spiral. The pitch angle is defined as the angle between the tangents to a spiral arm and a circle drawn from the center of the disk, it describes how tightly the spiral arms are wound. In the logarithmic case, the pitch angle is constant along all radii and it is given by µ = arctan(1/b), while for the Archimedean spiral, the pitch angle depends on the radius as µ = b/r. To test the symmetry of the spiral arms, the parameters r 0 and b were fitted separately, while we assumed i and PA to be fixed and equal to the literature values shown in Table 1. Therefore, we had four and eight free parameters for the SPHERE and ALMA data, respectively. To fit the data, we used the MCMC code based on emcee (Foreman-Mackey et al. 2013) and described in Kurtovic et al. (2018). A flat prior probability is used for the free parameters. For each fit, we used 250 walkers with two consecutive burning stages of 1000 and 500 steps, and then 1500 steps to sample the parameter space. The resulting pitch angles for each fit are given in Table 3, where for the SPHERE data, N corresponds to the northern spiral and S to the southern spiral; and in the case of the ALMA data, N1 corresponds to the northern inner spiral, S1 to the southern inner spiral, N2 to the northern outer spiral and S2 to the southern outer spiral. For the Archimedean fit, the pitch angles are estimated at 35 au. We find that there are some significant differences between the values of the inner and outer ALMA spirals for the logarithmic model, which leads us to conclude that the two sets might not be part of the same spiral arm. There seems to be additional structure in the region of the discontinuity, however, follow-up deeper observations would be needed to draw any conclusions of the origin of this break. Furthermore, we note that the values of the pitch angle for the corresponding spirals differ from one data set to the other, and that in the case of the Archimedean model, the scattered light pitch angles are slightly higher that those from submillimeter. This can be expected, since we are tracing different regions of the disk (the flared surface vs. the midplane). Moreover, we find that the pitch angles from the Archimedean model decrease with the distance from the star, in agreement with the results of Huang et al. (2018b). Figure 7 shows the spirals retrieved from both data sets in a polar map, along with the best-fit Archimedean model for each arm. We find a discontinuity in the spiral arms for the ALMA data at ∼50 au. This peculiarity has already been reported by Huang et al. (2018b), who also noted that this discontinuity appears to coincide with a region where there is additional bright emission between the main spiral arms, and comment that it could be explained by either the presence of a ring crossing this region, or "spurs" emerging from the main spirals.

Origin of the spirals
Large perturbations in the disk launch sound waves that result in a spiral shape due to the differential rotation of the disk. Theoretical models have shown that such large perturbations can be driven by the presence of an embedded planetary mass perturber (e.g., Boccaletti et al. 2013), GI (e.g., Goldreich & Tremaine 1979;Tomida et al. 2017), or a combination of both (e.g., Pohl et al. 2015). Two-arm spirals in disks can be driven by a massive, giant planetary companion ( 5M J ), that would typically be located at the tip of the primary arm. This scenario suggests, however, that these planets are fainter than predicted by "hot-start" evolutionary models (Dong et al. 2018b), since the number of detections is low.
On the other hand, one criteria to test the GI hypothesis is to use the so called Toomre (1964) parameter Q, given by where Ω k is the angular velocity, c s is the sound speed, and Σ is the surface density. These parameters vary with radius in the disk, resulting in Q being a function of the radius. The Toomre stability criteria states that a disk will be gravitationally stable if Q ≥ 1 and unstable if Q < 1.
In the case of WaOph 6, we estimated the Toomre parameter Q using equation 5 and assuming Σ ∝ 1/ √ r, where r is the distance from the star, throughout the disk (see Appendix D for the detailed calculation) in order to see whether the disk would be stable under these conditions. We obtained that Q varies from 2.2 − 33.2 from the outer part of the disk inward. This indicates that the disk is fairly stable according to the Toomre stability criterion, implying that a large perturbation driving the spiral arms should have come from a source other than GI. Nonetheless, we are aware that this is not an absolute proof of GI not taking place in the disk at some earlier evolutionary state. This analysis only shows that a disk with a surface density ∝ 1/ √ r around a 0.98 M star appears to be stable within our assumptions of the disk mass and the gas surface density, which should be taken with caution due to the considerable uncertainties in their calculation. In this context, we decide to test the planetary perturber hypothesis by performing hydrodynamical simulations and radiative transfer to compare with our observations.

Modeling
In order to test the hypothesis of a forming planet causing the spiral features, we performed 3D gas-only and 2D gas+dust hydrodynamical simulations in which a massive planet at a separations ≥ 90 au generates large-scale spiral arms interior to its orbit. To compare the results to our observations, we fed the resulting density distributions into a radiative transfer code to generate synthetic images.

Hydrodynamical models
To model the dust on the surface sampled by scattered light, we ran 3D simulations of the gas dynamics using the hydrodynamical code PLUTO (Mignone et al. 2007). The gas distribution was initially set following a vertically isothermal configuration at hydrostatic equilibrium (see, e.g., Fromang et al. 2011). The disk temperature depends on the cylindrical radius R from the center of the domain as T ∝ R −1/2 , whereas the gas pressure scale height is computed as H = c s /Ω K , where c s ∝ T 1/2 is the local sound speed and Ω K ∝ R −3/2 is the Keplerian angular velocity. With the chosen parameters, the disk aspect ratio depends on R as H/R ∝ R 1/4 , with H/R = 0.1 at R = 50 au. The gas volumetric density decreases with R as ρ ∝ R −7/4 , in such a way that the surface density varies as Σ ≈ √ 2πHρ ∝ R −1/2 . A locally isothermal equation of state is applied such that the initial temperature distribution is maintained through time.
The hydrodynamical equations were solved in spherical coordinates in a reference frame centered at the star-planet center of mass corotating with the system. The computational domain, given by the region (r, θ, φ) ∈ [8,210] au×[π/2−0.3, π/2+0.3]× [0, 2π], is discretized using a grid of resolution N r × N θ × N φ = 256 × 64 × 512 logarithmically spaced in the r-direction and uniformly spaced in the remaining ones. All computed fields are fixed to their initial values at the radial boundaries except for the density and the radial velocity, which are reflected and extrapolated in the ghost zones, respectively. On the vertical bound- aries, reflective conditions are applied. The gravitational potential is computed as the sum of the potentials of the star and the planet. To avoid divergences, the later is modified in the vicinity of the planet location following the prescription by Klahr & Kley (2006), employing a smoothing length equal to half the planet's Hill radius. For stability purposes, the planet mass is smoothly increased from 0 to its final value in a total time of 100 yr. We also included viscosity with constant α = 10 −3 . The resulting dust mass distribution was computed assuming a perfect coupling between the dust particles and the gas flow, with a uniform dust-to-gas mass ratio of 10 −2 .
To model the dust evolution in the midplane sampled mainly by the millimeter observations, we ran two-dimensional hydrodynamical simulations using the multi-fluid version of the code FARGO3D 4 (Benítez-Llambay & Masset 2016; Benítez-Llambay et al. 2019). It solves the Navier-Stokes equations of the gas and multiple dust species, each one modeled as a pressureless fluid that represents a specific grain size. We traced eight different dust species in our simulations. The initial gas temperature, gas surface density structure, equation of state and gas viscosity prescription of the 2D model are equivalent to our 3D simulations model, described above. The initial dust surface density in our simulation has the same structure as the gas surface density, while set by an initial dust-to-gas mass ratio of 10 −2 everywhere in the disk. We traced the dynamical evolution of eight dust fluids, that are logarithmically spaced in size, and followed a dust size distribution n(s) ∝ s −2.5 with minimum and maximum dust sizes of 10 µm and 100 µm. We set the dust intrinsic density to 2.0 g cm −3 . The dynamics of the dust fluids is dictated by its local Stokes number (dimensionless stopping time), defined as St = πa i ρ int /2Σ g , where a i is the grain size of the size-bin, ρ int the intrinsic grain density and Σ g the gas surface density. Dust diffusion was included in the simulation following the same prescriptions implemented in Weber et al. (2019), which are based on the results of Youdin & Lithwick (2007). Dust feedback onto the gas, dust growth, and dust fragmentation were not included in our simulations. The two-dimensional grid is linear in azimuth and logarithmic in radius, using 512 cells in φ covering 2π, and 256 cells in r covering from ∼ 8.4 au to ∼ 420 au. A planet was slowly introduced over 8 × 10 4 yr fixed at the given radius, driving spiral density waves in the disk. The planet's potential was smoothed by a length factor of 60% the disk scale height. For a more detailed description of the FARGO3D multi-fluid simulations see also Weber et al. (2019). Our test runs show that only dust fairly well coupled to the gas follows the spiral density waves (as shown by e.g., Veronesi et al. 2019;Sturm et al. 2020), with Stokes numbers below ∼ 10 −2 . Dust particles with larger Stokes number decouple from the gas and form axisymmetric rings. Fixing the disk gas surface density given the mass constraint from the observations, and the dust intrinsic density to a standard value, limiting the maximum dust size to 100 µm in our models is required to maintain the Stokes number of the dust observed at millimeter wavelengths below ∼ 10 −2 , therefore, the simulated dust particles trace the spiral arm structure. A summary of the parameters used in our simulations can be found in Table 4.
We sample the parameter space of a planet with masses between 2 − 15 M Jup and separations between 90 − 160 au (see Appendix E, Fig. E.1 for some of the resulting density maps). The lower limit in the mass range is chosen based on the results of Juhász et al. (2015), who concluded that observable spiral arms are formed for planets with M > 1M Jup . Tighter constraints on the lower limit of the planet mass can be obtained from spiral arm formation theory. For a planet with a mass larger than three thermal masses (M th ≡ c 3 s /ΩG = M (h/r) 3 p ), two spiral arms will form interior to its orbit (Bae & Zhu 2018).
Since we wanted to model m = 2 spirals, we chose to place the planets outside of the spiral arm (which extends up to 90 au in the millimeter continuum) based on the results of Dong et al. (2015). Assuming that the planet is outside ∼ 90 au and the disk aspect ratio of our model, we obtained that two spirals are formed for planet masses larger than ∼ 4.8 M Jup . Another criteria for the planet mass comes from the separation between the primary and secondary spiral arms (φ sep ). Fung & Dong (2015) obtained that this quantity scales with the planet mass, following φ sep = 102 • (q/0.001) 0.2 , where in this case q is the planet-tostar mass ratio. If we consider that the spiral arms have a separation range between 135 • and 180 • , we obtain that the planet mass should be between 4 and 17 M Jup . After a few test runs in our simulations, we realized that in order to observe the disk truncate at 90 au (consistent with the disk outer edge in the millimeter continuum), when increasing the separation, we should also increase the planet mass. Snapshots of the gas and total dust surface densities of our 3D and 2D hydro simulations for planets of 5, 10 and 15 M Jup at separations of 130, 140 and 160 au, respectively, are shown in the Appendix E.

Radiative transfer
In order to compare the results generated by the procedure described in the last section with our observations, we generated images in both polarized NIR and millimeter continuum using the radiative transfer code RADMC-3D (Dullemond et al. 2012). We obtained synthetic scattered light images using the dust mass distribution computed in the described 3D PLUTO simulations. Based on Ricci et al. (2010), we assumed a dust size distribution n(s) ∝ s −2.5 and model scattering by submicron particles with sizes ranging between 0.01 and 0.5 µm. To compute the dust mass in this range, we used the total dust mass estimated in this work (see Table 1) assuming maximum grain sizes in the mm, to obtain M dust,<0.5 µm = 10 −9 M . Opacities are computed assuming a dust composition of 60% astronomical silicates and 40% amorphous carbon grains, taking the optical constants respectively from Draine & Lee (1984) and Li & Greenberg (1997), and combining them following the Bruggeman mixing rule. Scattering matrices were computed assuming spherical dust grains using the BHMIE code (Bohren & Huffman 1983) for Mie scattering. For the scattered light computations, we approximated the grain size distribution using 5 size bins.
To smooth out oscillations in the polarization degree occurring when considering spheres of a single size (see, e.g., Keppler et al. 2018), we used a Gaussian size distribution within each bin with a FWHM of 20% of the corresponding grain size. The star was modeled as a point source located at the domain center emitting thermal radiation with characteristics summarized in Table 1. We used RADMC-3D to model anisotropic scattering with full treatment of polarization, using a total of 10 8 photon packages. The obtained Stokes Q and U frames were then convolved by a Gaussian PSF with a FWHM of 51 mas to reproduce the resolution of the VLT/SPHERE observations (see Section 3), after which we used equation (2) to obtain the resulting Q φ images.
To compare with the ALMA data, we computed radiative transfer predictions of the dust continuum, in this case, using the output of the dust and gas 2D simulation. We used the dust density field from the simulation as input for RADMC-3D. We expanded the two-dimensional surface density vertically, assuming a Gaussian shape, where the volumetric mass density for each dust bin follows: where H i indicates the pressure scale height of the dust bin. The vertical settling of the disk follows a standard diffusion model (Dubrulle et al. 1995): where H g is the gas pressure scale height, S t is the dust Stokes number, andα = α/S c z with α the α-viscosity value of the gas. S c z is the Schmidt-number, set to 1 S c z relates the dust diffusion coefficient with the gas viscosity D z = ν/S c z (see also Weber et al. (2019)). We used optool 5 to compute the dust asorption and scattering opacities of a mixture using standard Mie theory and Bruggeman rules. We assumed that the composition of the dust grains is a mixture of silicates (internal density of 3.2 g/cm 3 ), amorphous carbon (internal density of 2.3 g/cm 3 ), and vacuum. Assuming that the solids in the mixture are 60% silicates and 40% carbon, a volume fraction of 25% of vacuum in the mixture is required so its internal density is ∼ 2 g/cm 3 . The dust size distribution is equal to the values used for the simulation, set by the power law n(s) ∝ s −2.5 , with maximum and minimum dust sizes of 100 µm and 10 µm, respectively. The total dust mass in our models is ∼ 10 −4 M . We computed the dust temperature using the Monte Carlo method of Bjorkman & Wood (2001), and the continuum emission image via ray-tracing, taking into account absorption and scattering, assuming Henyey-Greenstein anisotropic scattering. We computed simulated ALMA observations from the radiative transfer synthetic continuum image using CASA (version 5.6) simobserve and tclean tasks. Following the observations setup from the DSHARP survey (Andrews et al. 2018), we simulated an 8 h integration in configuration C43-8 combined with a 15 min integration in C43-5. Finally, we cleaned the image using briggs weighting 1.0. We obtained a beam size of 55 × 53 mas and PA of ∼ −55 • , directly comparable to the ALMA observation.

Results and comparison to observations
All tested planets drive m = 2 spiral arms whose symmetry increases for larger planetary mass and have a low contrast in the dust surface density (as seen in the density plots shown in Fig. E.1). Given the asymmetry in the 5 M Jup case, we conclude that in case the spirals are caused by a planet, its mass should be at least of approximately 10 M Jup . In Fig. 8 we show resulting radiative transfer images for a 10 M Jup planet at a separation of 140 au, with spiral arms observable both in the scattered light and millimeter continuum observations. For a better comparison to the simulations, we apply a Gaussian kernel to the image on the left panel, similar to the one used for the SPHERE image in Fig. 6, left; and we subtract the azimuthal average flux on the image plane to enhance the spirals on the synthetic millimeter continuum image on the right, analogous to the procedure applied to Fig. 6, right. Additionally, we overplot the location of the observed spiral arms. The obtained images resemble the ones detected both in the scattered light and millimeter continuum observations, except for the fact that we are only able to fit either the inner or the outer spirals from the millimeter observations, but not both at the same time (see Fig. E.2, lower panel). This is likely due to missing physics in our simulations, as these models of spirals launched by a single planet are unable to reproduce the break in the spirals observed by ALMA, as well as the gap and the ring features (at 79 and 88 au, respectively) in the observations. We must also note that in order to see spiral arms induced by a planet in millimeter continuum, the dust must be fragmentation limited (e.g., Birnstiel et al. 2010) leading to a small dust maximum size, and therefore, to Stokes numbers small enough to follow the spirals. In protoplanetary disks, the maximum grain size is mainly set by radial drift or fragmentation of particles after collisions. The later depends on the disk viscosity and the threshold considered for the fragmentation velocity of the grains. Assuming low fragmentation velocities for ice grains (e.g., < 1 m/s, as suggested by recent laboratory experiments such as Musiolik & Wurm 2019;Steinpilz et al. 2019), and alpha=10 − 3 (as taken in the simulations), the maximum grain size in the entire disk is dominated by fragmentation, limiting the maximum size of 100µm (Pinilla et al. 2021). Kataoka et al. (2016) have found that dust with similar characteristics is traced by millimeter continuum observations of the similarly young disk HL Tau. These characteristics are not required to see spirals generated by GI, where the dust trapping in spirals is efficient for larger dust Stokes numbers (Rice et al. 2004). We also note that none of the parameter sets that we sample are able to reproduce the contrast nor the apparent break in the spiral arms shown in the ALMA data, which might be explained by additional physical processes occurring in the disk. However, more complex simulations including other effects (e.g., dust growth, fragmentation, dust feedback, gas temperature evolution) are beyond the scope of this work. Additionally, we note that a planet of 10 M Jup in such a young disk could have either formed via gravitational collapse when the disk was probably more massive and, therefore, gravitationally unstable (Boss 1997), or formed as a stellar companion from cloud fragmentation due to the planet/star mass ratio (∼1%, Reggiani et al. 2016). We would like to mention that this is a first attempt to find a plausible planetary model to explain the observed spiral pattern in the protoplanetary disk around WaOph 6 and that further, deeper observations would be needed to confirm or discard this scenario.
Since we employ an isothermal equation of state, the spirals produced in our simulations are induced solely by Lindblad resonances and not by buoyancy modes, which may be triggered when using finite cooling times. It is argued in Bae et al. (2021) that such modes cannot be observed in millimeter continuum observations, but could potentially be seen in scattered light. The pitch angles for buoyancy resonances shown in that work for up to 2 M Jup planets are generally below those seen in our SPHERE observations (see Table 3), which suggests that the observed spi-rals are likely not triggered by such modes. Future resolved CO line emission observations analyzing the disk kinematic structure could help discard or verify this hypothesis (Bae et al. 2021).

Observations in NIR and mm
Our SPHERE/IRDIS-DPI observations show the launch of an m = 2 spiral pattern in the disk around WaOph 6. This is a surprising finding, since so far, no spiral arms had been observed in scattered light in disks around K and/or M stars with ages < 1 Myr. Moreover, spiral arms have not been observed at these wavelengths in single T Tauri stars of any age . Disks with spiral arms detected in scattered light are thought to be older (with the caveat that stellar ages are highly uncertain), and with stellar hosts of spectral types from G to A (e.g., MWC 758, Dong et al. (2018a), HD 142527, Claudi et al. (2019, HD 100546, Pérez et al. (2020), AB Aur, Boccaletti et al. (2020), HD 100453, Benisty et al. (2017)). In the millimeter continuum, most of these disks show asymmetric morphologies, along with large cavities (e.g., Tang et al. 2017;Cazzoletti et al. 2018;Pineda et al. 2019). Further observations in both scattered light and millimeter continuum of K and M type stars with disks would be needed to determine whether spiral arms are a common feature in such young disks, as well as the possible implications that this might have in dust and gas evolutionary models. We also note that comparing observations at different wavelengths can contribute greatly to the understanding of the physical processes driving the different morphologies seen in protoplanetary disks.
From our hydrodynamical simulations, we observe that in order to obtain a spiral pattern that can be observed in the millimeter continuum data, the dust particles must have a limited maximum size. This has previously been observed in dust evolution simulations by Gerbig et al. (2019), and can be linked to the young age of the disk.

Upper limits on the brightness of point sources
We used the total intensity image derived from our SPHERE/IRDIS-DPI observations to obtain information on the detection limits for WaOph 6. We built the contrast curve in Fig. 9 by considering the contrast between WaOph 6 (the central brightest pixel) and a representative planetary signal in the total intensity image. We took the planetary signal to be three times the noise (root mean square) in 2 pixel wide annuli centered on the star, at different separations up to ∼6" (∼740 au). Additionally, we estimated the foreground extinction in the H-band toward WaOph 6. For this we first estimated the reddening by using the intrinsic J − H magnitude of a K6V star from Pecaut & Mamajek (2013), then using the values in Table  3 of Rieke & Lebofsky (1985), we obtained a visual extinction of A V = 5.08 mag, and an H-band extinction of A H = 0.88 mag.
To estimate the apparent magnitude of our proposed planet, we used two independent evolutionary model predictions. On one hand we considered the evolutionary models by Baraffe et al. (2003) for a 10 M Jup planet at 1 Myr. On the other hand, we used the evolutionary models proposed by Spiegel & Burrows (2012) for both a "hot" and "cold-start" scenarios, and we ex- nitude for WaOph 6, m H = 7.57 mag, we obtained the following contrasts: ∆mag= 7.70, ∆mag= 7.36, and ∆mag= 15.49 mag, respectively. Furthermore, we obtained the contrasts for our proposed planet in the "warm-start" scenario from the initial entropy values reported by Spiegel & Burrows (2012). And as an additional comparison, we used the Bern EXoplanet cooling curves (BEX, Mordasini et al. 2017) coupled with the COND atmospheric models (Allard et al. 2001) reproducing the cooling under "warm-start" initial conditions (Marleau et al. 2019), and thus denominated BEX-WARM model (see Asensio-Torres et al. 2021, and references therein for more details). As seen in Fig. 9, the detectability of our proposed planet strongly depends on the adopted formation model. In case of the "hot-start" scenario, the planet should have been observed, while for a large part of the "warm-start" and for the "cold-start" scenarios, the planet contrast is below our detection limits. Based on the planet mass and location, a "warm" to "cold" start model would be more plausible to explain its existence. Additional detection limits for WaOph 6 in the L'-band (λ 0 = 3.8 µm) have been recently reported by Jorquera et al. (2020). They do not detect any companion candidates to the star, but report detection probability maps obtained using the (Baraffe et al. 2003) models. From these they preliminary rule out the presence of companions with masses > 5 M Jup at separations > 100 au. However, they advise that these estimates might be optimistic, since they do not consider extinction effects, either toward WaOph 6, nor due to the disk dust. An additional caveat comes from the models, as they become very uncertain in accurately predicting the properties of very young planets. It is also important to note that our hydrodynamical simulations do not include additional physical processes that could be ongoing in the disk, coming from the fact that spiral arm formation by a planetary mass object is still not well understood. This could lead to an overestimation of the planet mass, which along with evolutionary models uncertainties, could explain our differing results.

Summary and conclusions
We have presented for the first time scattered light SPHERE/IRDIS-DPI observations of the disk around WaOph 6 in the H-band. We analyzed the disk morphology, and used archival ALMA data to compare with ours. We tested the planetary mass perturber hypothesis as the underlying cause for the spiral structure by performing hydrodynamical simulations and using radiative transfer. Our results are summarized below: 1. We observe the launch of a set of m = 2 spiral arms up to ∼0.3" (40 au) in our Q φ SPHERE/IRDIS-DPI images as seen in Fig. 2, left. These spirals were first detected using millimeter continuum observations from the ALMA/DSHARP survey. To our knowledge, WaOph 6 is the youngest disk to show spiral features in scattered light . We note that this might be of interest for dust and gas evolutionary models.
2. We observe a companion candidate at about 3" from the star in our data, as shown in the top panel of Fig. 4. After the astrometric analysis described in Section 4.2, we were able to determine that the CC is not bound to WaOph 6. With this we also discard the CC being a possible cause of the spiral structure.
3. Comparing our SPHERE observations with archival ALMA/DSHARP data, we find that both the gap and the ring features at 79 and 88 au, respectively, seem to be present in both data sets. We traced the spiral features in both observations as seen in Fig. 6. For the ALMA data, we notice a break in the spiral arms of the ALMA image at ∼0.16" (∼20 au), which is not observed in our SPHERE data. We treated this break as a separate set of spirals, however, its origin remains unknown. When plotting the spirals in polar coordinates (Fig. 7) we find a discontinuity in the spiral arms for the ALMA data at ∼50 au, already reported by Huang et al. (2018b).
4. To test the planetary mass perturber hypothesis we performed hydrodynamical simulations combined with radiative transfer to compare with the observations. We tested the parameter space of a planet with masses between 2 − 15 M Jup and separations between 90 − 160 au (i.e., outside of the spiral structure). All tested planets drive m = 2 spiral arms. However, none of the parameter sets that we sample are able to reproduce the contrast nor the apparent break in the spiral arms shown in the ALMA data, which may be due to additional physical processes occurring in the disk. Furthermore, the tested planets do not reproduce the gap nor the ring features at 79 and 88 au, respectively, these features need further investigation outside the scope of this work. Given the symmetry of the observed spirals, we find that, if these are caused by a planet, its mass is likely of at least 10 M Jup . This is a first attempt to explain the spiral structure seen in both data sets, and more data are needed to better constrain the underlying cause of the spiral features.

5.
To determine the sensitivity of our data to possible companions embedded in the disk, we generated the contrast curve in Fig. 9 from the total intensity image. With this we obtain contrast limits for a planetary/substellar companion forming inside the disk in polarized light. We estimate the contrast of our proposed planet using different evolutionary models, where the possibility of detection strongly depends on the formation scenario. A "warm" to "cold" starts would explain the nondetection of the planet in our SPHERE data.
In conclusion, the findings in this work highlight the still unknown complexity of WaOph 6. The striking presence of a spiral pattern in scattered light even in limited S/N data are worth further, deeper observations of this source. Which will additionally serve to confirm or discard a planetary perturber as a possible cause behind the spiral features.  Fig. E.2: Resulting radiative transfer images from our 3D (upper panels) and 2D (lower panels) hydrodynamical simulations for the planets described in Fig. E.1. The mass and separation increase from left to right. The images have been processed with the same techniques as the ones in Fig. 8 for a better comparison to the observations. For the lower panels, we show the images that fit the inner spirals.