Asymmetry, Gap Opening and High Accretion Rate on DM Tau:
A Hypothesis Based on Interaction of Magnetized Disk Wind with Planet

Abstract

Over two hundred protoplanetary disk systems have been resolved by ALMA, and the vast majority suggest the presence of planets. The dust gaps in transition disks are considered evidence of giant planets sculpting gas and dust under appropriate disk viscosity. However, the unusually high accretion rates in many T Tauri stars hosting transition disks challenge this theory. As the only disk currently observed with high turbulence, the high accretion rate ($\sim 10^{-8.3}M_{\odot}/yr$) observed in DM Tau indicates the presence of strong turbulence may within the system. Considering the recent theoretical advancements in magnetized disk winds is challenging the traditional gap-opening theories and viscosity-driven accretion models, our study presents a pioneering simulation incorporating a simplified magnetized disk wind model to explain the observed features in DM Tau. Employing multi-fluid simulations with an embedded medium mass planet, we successfully replicate the gap formation and asymmetric structures evident in ALMA Band 6 and the recently JVLA 7 mm observations. Our results suggest that when magnetized disk wind dominate the accretion mode of the system, it’s entirely possible for a planet with a medium mass to exist within the gap inside 20 au of DM Tau. This means that DM Tau may not be as turbulence as imagined. However, viscosity within the disk should also contribute a few turbulence to maintain disk stability.

UTF8gbsn

1 Introduction

The Atacama Large Millimeter/submillimeter Array (ALMA) has unveiled a wide array of intricate structures within mm/sub-mm observations of protoplanetary disks (PPDs), showcasing features such as cavities, gaps, rings, and vortices, as reported in numerous studies (e.g., ALMA Partnership et al., 2015; Long et al., 2018, 2019; Andrews et al., 2018; Huang et al., 2018; Lodato et al., 2019; Andrews, 2020). Notably, the pervasive presence of luminous rings alongside dim gaps stands out. Although some recently studies suggest that substructures can be generated by specific dust distributions without relying on planet-disk interactions (e.g., Pinilla et al., 2021; Kalyaan et al., 2021; Stadler et al., 2022; Jiang & Ormel, 2023). However, traditionally explained by the influence of planets embedded within the disk, altering the distributions of gas and dust (e.g., Paardekooper et al., 2022).

The interaction of giant planets with the disk leads to the generation and attenuation of density waves at resonances, creating noticeable depletions in the gas density nearby (Goldreich & Tremaine, 1980; Lin & Papaloizou, 1986, 1993). This interaction results in the formation of a pressure bump at the gas gap’s outer boundary, effectively halting the inward migration of dust particles poorly coupled to the gas (Paardekooper & Mellema, 2006; Zhu et al., 2011), hence the formation of pronounced dust rings. By adjusting parameters such as disk viscosity, thickness, and planet mass, one can replicate an assortment of dust emission patterns observed in the sub-mm ALMA observations (e.g., Rosotti et al., 2016; Zhang et al., 2018; Chen et al., 2021; Garrido-Deutelmoser et al., 2023; Wu et al., 2023a).

In the process of planet-disk interaction, the turbulent viscosity parameter $\alpha_{\rm v}$ (Shakura & Sunyaev, 1973) plays a role in maintaining the disk profile and preventing gap opening. However, current observational results suggest that the planet-forming regions in PPDs are not as turbulent as previously thought (Flaherty et al., 2017; Lodato et al., 2017; Dullemond et al., 2018; Rosotti et al., 2020; Doi & Kataoka, 2021). Although the methods used for these observations may vary, and moreover, some gas observations trace the upper disk layer in the outer parts of the disk. These latest observations can restrict the upper limit of the disk’s turbulence coefficient to $10^{-4}\sim 10^{-3}$, meaning that if PPDs are indeed weak turbulence or nearly inviscid, then a viable alternative explanation is needed for how disks accrete and evolve.

Recent theoretical work proposes that magnetically-driven disk winds might transport the angular momentum of gas at a faster rate than disk turbulence, dominating the accretion in PPDs (Bai, 2011, 2013, 2015; Armitage et al., 2013; Bai & Stone, 2013, 2014; Suzuki et al., 2016; Bai et al., 2016; Hasegawa et al., 2017; Lesur, 2021; Cui & Bai, 2021, 2022). This could be a promising alternative to reconcile the relationship between turbulence velocity and host stellar accretion rates.

Based on numerical simulation of planet-disk interactions in with MHD disk winds, Elbakyan et al. (2022) and Aoyama & Bai (2023) have shown that for Type-II planets, gap opening is easier under the same level of turbulence when magnetohydrodynamics (MHD) disk wind dominates accretion, compared to the traditional viscosity-driven accretion mode. And Wu et al. (2023b) has further discovered that when MHD disk wind dominates in the angular momentum transfer within the disk, different degrees of asymmetric structures can be generated depending on the strength of the wind. These structures, distinct from those produced in viscosity-dominated disks (Zhang et al., 2018; Long et al., 2022), can be directly imaged by ALMA.

In this work, we use a simplified MHD disk wind prescription for the first time to simulate a case study of DM Tau—the well studied transition disk around the T Tauri star (e.g., Chung et al., 2024). We focus on reproduce its substructures, and keep it with a high accretion rate. By embedding a medium mass ($\sim 53$ times Earth mass) planet and using multi-fluid in simulations, we reproduce the gap opening and asymmetric structures caused by planet-disk interaction, and make synthetic images to compare with ALMA Band 6 (230 GHz) and Karl G. Jansky Very Large Array (JVLA) Q Band (44 GHz) observations. It need be noticed that the hypothesis of multi-planet cannot be ruled out in DM Tau system (Xu & Wang, 2024). However, for the sake of simplicity, we only consider the case of single-planet in this exploratory work.

The paper is organized as follows: In Sec.2, we introduce the characteristics of DM Tau. In Sec.3, we introduce the numerical setup and parameters for implementing hydrodynamic simulations and making synthetic maps. In Sec.4, we describe the comparison between simulated observational images obtained using the MHD disk wind model and actual ALMA and JVLA observations. The discussions and conclusions are detailed in Sec.5 and Sec.6.

2 The DM Tau Disk System

DM Tau is a T Tauri star with a high active accretion level $\sim 10^{-8.3}M_{\odot}/yr$, located at 145 pc far from us with a mass of 0.53 $M_{\odot}$ (Gaia Collaboration et al., 2023). It hosts a transition disk inclined at i $\sim$ $36.1^{\circ}$ and extending to 120 au in dust emission, which radius encloses 90% of the dust disk emission.

Recent observational studies have confirmed the dim gap at approximately 4 to 21 au from the center of DM Tau, with a bright ring at 21 au that features some asymmetric clumpy structures (Kudo et al., 2018; Hashimoto et al., 2021). This suggests that there might be forming planets within the 4-21 au region. This type of structure, also observed in some other disks. In general, it might be generated by a super-Earth-mass planet residing in a system with low viscosity (Dong et al., 2017, 2018; Pérez et al., 2019; Francis & van der Marel, 2020).

However, if the high accretion rate implied by observations is generated by MRI viscosity, it would mean that this source is not a low-viscosity but a high-viscosity system (Delage et al., 2022). This would require a sufficiently massive planet (may be nearly 10 times the mass of Jupiter) to produce such structures (Pinilla et al., 2012; Francis et al., 2022). And the massive planets will carve out eccentric gaps and generate streamer structures (Chen et al., 2024), thus the hypothesis of a single-massive-planet would have difficulty generating the observed asymmetric structures in DM Tau.

At the same time, recent observations have shown that the two asymmetric structures on the DM Tau ring exhibit different observational characteristics at different wavelengths. Specifically, in previous high-resolution ALMA Band 6 observations (Kudo et al., 2018; Hashimoto et al., 2021), the asymmetric structure on the west side of DM Tau was brighter than the one on the south side. However, in the latest JVLA 7 mm observations (Liu et al., 2024), the asymmetric structure on the south side of DM Tau is much brighter, showing a major asymmetry, while the west side asymmetric structure has become very dim. These observational results suggest that the formation of substructures in DM Tau requires a new theoretical explanation.

3 Numerical Setup

To use MHD disk wind model to fit the observable signatures on DM Tau, we use the hydrodynamic grid-basics code FARGO3D (Benítez-Llambay & Masset, 2016) with its multi-fluid version (Benítez-Llambay et al., 2019).

3.1 Hydrodynamic Model

Our work consider a 2D cylindrical polar coordinate ($r$, $\varphi$), locally isothermal protoplanetary disk system with an embedded planet, composed of gas and dust. We neglect the disk self-gravity, planet migration and accretion in all of our simulations. The disk model in our simulations is described by the Navier-Stokes equations:

$$\frac{\partial\Sigma_{\rm g}}{\partial t}+\nabla\cdot\left(\Sigma_{\rm g}\textbf{\emph{V}}\right)=0\;,$$ (1)

$$\frac{\partial\bm{V}}{\partial t}+\left(\bm{V}\cdot\nabla\right)\bm{V}=-\nabla\Phi-\frac{\nabla P}{\Sigma_{\rm g}}-\frac{\sum\limits_{i}\Sigma_{\rm d,\it i}\bm{F}_{i}}{\Sigma_{\rm g}}-\nabla\cdot\tau+\bm{S}\;,$$ (2)

$$\frac{\partial\Sigma_{\rm d,\it i}}{\partial t}+\nabla\cdot\left(\Sigma_{\rm g}\bm{W}_{i}+\bm{j_{i}}\right)=0\;,$$ (3)

$$\frac{\partial\bm{W}_{i}}{\partial t}+\left(\bm{W}_{i}\cdot\nabla\right)\bm{W}_{i}=-\nabla\Phi+\bm{F}_{i}\;,$$ (4)

where $\Sigma_{\rm g}$, V and $P$ are the surface density, velocity and pressure of the gas, respectively. The $\Phi$ is the gravitational potential from the central stellar and the planet. The $\Sigma_{\rm d,\it i}$, $\bm{W}_{i}$, $\bm{j_{i}}$ and $\bm{F}_{i}$ is the surface density, velocity, dust diffusion flux and the drag force from gas acting on the dust of unit mass for $i$-th dust species, respectively. And the $\tau$ is viscous stress tensor.

In general, dust is turbulently diffused in viscous disks, with its mass flux following the gradient of gas density to replenish the gap. However, in MHD windy disks, the wind solely extracts angular momentum, leading to a consistently inward gas velocity. Based on the prescription given by Elbakyan et al. (2022), which is an extension from 1D prescription Tabone et al. (2022a, b), we introduced the wind accretion term $\bm{S}$ in Eq. 2, to simulate the role of the MHD disk wind in the angular momentum transfer process within the disk:

$$\bm{S}=\Gamma\frac{\hat{\bm{r}}}{r},$$ (5)

here ${\hat{\bm{r}}}$ is set to the unit vector in radial direction. The $\Gamma$ is a distinct torque, that serves to reduce the rotational angular momentum of the gas:

$$\Gamma=\frac{1}{2}\sqrt{\frac{GM_{*}}{r}}V_{\rm dw}\;.$$ (6)

For DM Tau, we set $M_{\star}=0.53~{}M_{\odot}$ of the host star in our simulations (Gaia Collaboration et al., 2023), and $V_{\rm dw}$ is the characteristic radial inward velocity of gas due to MHD disk winds corresponding to this torque follow Tabone et al. (2022a):

$$V_{\rm dw}=-\frac{3}{2}\alpha_{\rm dw}h^{2}V_{\rm K}\;,$$ (7)

here $\alpha_{\rm dw}$ is a dimensionless parameter to describe the strength of the disk wind, defined in a similar manner to the turbulent viscosity parameter $\alpha_{\rm v}$. Therefore, the total accretion level in the disk is $\alpha_{\rm acc}=\alpha_{\rm dw}+\alpha_{\rm v}$, see Tabone et al. (2022a) for more details. And $h$ is the disk aspect ratio, $V_{\rm K}$ is the Keplerian velocity at $r$.

To summarize, the impact of wind torque (Eq.6) translates to an augmented radial velocity in the gas (Eq.7), consistently directed towards the host star. This wind-induced radial gas motion can match or exceed the viscous velocity of the gas, contingent on the magnitude of $\alpha_{\rm dw}$. This effect, in turn, significantly affects the radial drift velocity of dust particles, facilitating the transport of pebbles across any gap (Wu et al., 2023b).

3.2 Numerical Parameters

The simulation spans a computational domain from 0.3 $r_{0}$ to 5.0 $r_{0}$, with $r_{0}=20$ au denoting the unit length in code, the planet is also fixed at this radial position. The domain resolution includes 512 radial cells and 656 azimuthal cells. This resolution comes from previous simulations of 2D MHD disk winds using FARGO3D (Kimmig et al., 2020; Elbakyan et al., 2022; Wu et al., 2023b).

We set the initial gas surface density profile to

$$\Sigma_{\rm g}(r)=\Sigma_{0}({r}/{r_{0}})^{-1}\;,$$ (8)

where $\Sigma_{0}=4\times 10^{-4}M_{\star}/r_{0}^{2}$ in code units, this means that the disk mass is about 0.005 $M_{\odot}$ within a range of $\sim$100 au, which is a little lower than the value of recently estimates for DM Tau (Trapman et al., 2022), given by a novel way using N₂H⁺ and C¹⁸O.

In our model featuring a local isothermal temperature distribution, the aspect ratio is determined as follows:

$$h(r)=c_{s}/V_{\rm K}=H/r=h_{0}(r/r_{0})^{0.25},$$ (9)

here we set $h_{0}=0.03$, which means a flatter disk compared to some other disks is considered in our simulations.

In our model, we introduce five distinct dust populations, with particle sizes ranging logarithmically and uniformly from 30 $\mu$m to 3 mm. The dust feedback is included in our simulations. The size distribution at the onset follows a power law of $s^{-3.5}$, as suggested by (Mathis et al., 1977). We set the initial dust-to-gas mass ratio to $\epsilon=1\%$, and assume a uniform internal density of $2.3~{}g/cm^{3}$ for the dust grains.

In terms of boundary conditions, we implement wave-damping boundary conditions for the gas component, where the gas density and velocity are gradually adjusted back to their original radial profiles, to maintain simulation stability. For the dust components, we employ open boundary conditions at the inner edge of the grid, specifically for their radial velocities. This is to prevent dust from gathering on the inner region, thus avoiding the influence of boundary effects on subsequent radiation transfer processing and simobserve.

3.3 Simobserve

We utilized the publicly software package RADMC-3D code (Dullemond et al., 2012) to conduct our dust radiative transfer calculations. Then to process our FARGO3D simulation results for use as inputs in RADMC-3D, we used the publicly available code fargo2radmc3d (Baruteau et al., 2019). The physical details of the radiative transfer calculations are described in the Appendix A.

To calculate the dust absorption and scattering opacities, we apply Mie theory, utilizing OpTool (Dominik et al., 2021) for the computations. The dust composition follows that used in the DSHARP project (Birnstiel et al., 2018; Zhang et al., 2018), encompassing both absorption and scattering total opacities for specified wavelengths and grain sizes. For the determination of the dust continuum emission’s specific intensity, we deploy ray tracing with $10^{8}$ photons for scattering.

Based on the known observation data of DM Tau, we assume the disk distance is 145 parsecs far from us (Piétu et al., 2007; Gaia Collaboration et al., 2021, 2023). And the star radius is set to 0.60 solar radii with a stellar effective temperature of 3705 K (Gaia Collaboration et al., 2023). For the geometry of the disk, we set the disk inclination of DM Tau is $35.2^{\circ}$ and the disk position angle is $158^{\circ}$ (Kudo et al., 2018).

To mock ALMA and JVLA observations, we used the SIMOBSERVE and SIMANALYZE functions within the CASA 6.4.9 package (McMullin et al., 2007), converting synthetic images into Fourier domain visibility data sets. We utilized the alma.out28.cfg and vla.a.cfg files provided within the CASA package, which include the longest baselines for ALMA (extending up to 16 km) and JVLA (extending up to 21 km). We adjusted the “inbright” and “integration” parameter to make The root-mean-square (RMS) noise level produced by CASA as close as possible to the values generated by realistic observations (Hashimoto et al., 2021; Liu et al., 2024).

4 Results

Wu et al. (2023b) explored a wide range of $\alpha_{\rm v}-\alpha_{\rm dw}$ parameter space based on a fixed planet mass. In this exploratory work, we use the same planet-to-stellar mass ratio $q=0.0003$ as in that parameter space. Which corresponds to approximately 53 times Earth mass in DM Tau. As introduced in Sec.3.2, we placed this planet in a fixed orbit at a distance of 20 au from the stellar.

We note that some recent studies (e.g., McNally et al., 2019; Wu et al., 2023a) suggest that planet migration may lead to the formation of different substructures within the disk, and the migration behavior of planets in a windy-disk differs from that in a traditional laminar disk (Kimmig et al., 2020). However, we have overlooked the impact of planet migration in this work like most similar works in the past (e.g., Zhang et al., 2018; Ricci et al., 2018; Wu et al., 2023b, 2024c). In next work, we will explore planet migration within MHD disk wind (Wu & Chen 2024 in prep).

All of our simulations were run for 1500 orbital periods at $r_{0}$ (corresponding to $\sim 0.2$ Myr for physical units), which is a reasonable value for the age of DM Tau (Simon et al., 2000; Luhman et al., 2010; Guilloteau et al., 2014).

Combining the parameter space study from Wu et al. (2023b) with the morphological characteristics of DM Tau in dust continuum emission, we conclude that the roles of MHD disk wind and viscosity in the angular momentum transfer process of DM Tau should be well-matched in strength, i.e. $\alpha_{\rm v}\sim h\alpha_{\rm dw}$ (see Figure 2 in Wu et al., 2023b). Based this, we obtained our fiducial model (as shown in Fig.1 and Fig.2) by adjusting $\alpha_{\rm dw}$ and $\alpha_{\rm v}$. The numerical parameters used in the fiducial model are listed in Tab.1.

Figure 1: Left panel: Previously published ALMA Band 6 (230 GHz, $\lambda=1.3$ mm) continuum image on DM Tau (Hashimoto et al., 2021) from ALMA data $\#$2018.1.01755.S (the project code marked at the upper-left corner), with the $0.027\times 0.018$ arcsec beam size (beam position angle is $33^{\circ}$). Right panel: Simulated ALMA Band 6 continuum images of our fiducial model (the values of $\alpha_{\rm dw}$ and $\alpha_{\rm v}$ are written in the upper-left corner) with a 53 times Earth mass planet at 20 au (marked by white “x”) from the central star, with the $0.02\times 0.02$ arcsec beam size. We mark the beam size at the left-lower corner for each panel. The RMS noise in each panel is at the same level, i.e. RMS = $8.4\mu Jy/beam$. Since the inner boundary of our simulation domain is at 6 au, to prevent the influence of boundary effects, we’ve applied a mask with a black ellipse to each panel. Additionally, we’ve marked the position of the stellar with a white “+” for reference. In each panel, the approximate outline of the bright ring around 21 au is delineated with black solid line contours, and the high-brightness clumps on the ring, representing the asymmetric structures, are indicated with white dashed line contours.

Figure 2: Similar to the two top panels in Fig.1, but for JVLA Q Band (40$\sim$48 GHz, $\lambda=7$ mm) continuum image (Liu et al., 2024, the left panel), and simulated JVLA continuum image (the right panel) with our fiducial model.

Table 1: The parameters used in fiducial model for DM Tau. Where $r_{p}$ and $M_{p}$ are the planet position and planet mass, respectively. $M_{\oplus}$ is the Earth mass and planet to stellar mass ratio is $q$.

Parameter	Value
$r_{p}$	20 au
$M_{p}$	53 $M_{\oplus}$
$M_{\star}$	0.53 $M_{\odot}$
$q$	$3\times 10^{-4}$
$\Sigma_{0}$	$4.7g/cm^{2}$
$h_{0}$	0.03
$\alpha_{\rm dw}$	0.076 ($\sim 97\%$ of $\alpha_{\rm acc}$)
$\alpha_{\rm v}$	$0.002$ ($\sim 3\%$ of $\alpha_{\rm acc}$)

4.1 Observation at 1.3 mm and 7 mm

As shown in Fig.1, when $\alpha_{\rm dw}=0.076$ and $\alpha_{\rm v}=0.002$, our fiducial model matches the observations very well. At this point, the MHD disk wind contributes $\sim 97\%$ of the accretion rate in the system, while viscosity accounts for the remaining $\sim 3\%$. Under this scenario, some of the characteristic substructures reflected in ALMA observations of DM Tau, the uneven bright ring at 21 au, the high-brightness clumps on the ring causing asymmetric structures, and the wide and deep gap spanning several au, are all represented in our fiducial model.

What is even more remarkable is that our fiducial model continues to show a disk morphology that matches the observed results at longer wavelengths. As shown in Fig.2, the high-contrast clump on the ring observed in the recently JVLA 7 mm observations (Liu et al., 2024), which causes asymmetry, is also represented in our fiducial model. It also accurately reproduces the uneven distribution of emission on DM Tau’s dust ring. The situation where some clumps with high contrast at 1.3 mm observations appear dimmer in the 7 mm observations is also reflected in the fiducial model.

The substructures in the disk displayed in the fiducial model originate from planet-disk interactions as well as the roles of wind-viscosity in the angular momentum transfer process. We will provide a brief discussion based on the simulation results in Sec.5.2.

However, we noted that compared to our simulation, the latest 7 mm observations show the gap within 20 au of DM Tau to be a bit more shallow. The reason for this observation discrepancy is that our radiative transfer calculations at 1.3 mm and 7 mm wavelengths considered only dust thermal emission. Yet, Liu et al. (2024) concluded that at 7 mm wavelengths for DM Tau, flux densities are contributed by both dust thermal emission and time-varying free-free emission of gas (e.g., Zapata et al., 2017; Liu et al., 2021). The differences between observations and simulations in this regard precisely validate this conclusions about varying emissions contributed at different wavelengths in DM Tau.

We also remind readers that due to a lack of detailed parameters for DM Tau, our fiducial model currently cannot produce synthetic images at 1.3 mm wavelength that match the actual observations. This might be because we currently do not have a sufficient understanding of the dust species composition in DM Tau. This will require enhancements based on further observational data and analysis of DM Tau in the future.

4.2 Pure Wind or Pure Viscosity Cases

We also plotted pure viscosity cases and pure MHD disk wind cases as reference models. As shown in Fig.3, when $\alpha_{\rm v}=\alpha=0.078$, a planet with 53 $M_{\oplus}$ is unable to carve out a gap that is sufficiently deep and wide in the disk to allow for the formation of a bright dust ring behind the planet. Conversely, when $\alpha_{\rm dw}=\alpha=0.078$, the strong wind interacting with the planet causes dust to accumulate on the ring in a crescent-shaped structure, creating a very pronounced and intense asymmetric structure that contrasts sharply with the surrounding ring, which does not match current observations.

The results of the reference models indicate that the accretion in DM Tau must stem from the interaction between wind and viscosity, where the wind contributes the majority of the accretion rate ($\sim 97\%$) to achieve the formation of asymmetric structures and gap opening. Although the contribution of viscosity to accretion is only a small fraction ($\sim 3\%$), it plays an important role in maintaining the stability of the disk structure.

Figure 3: Top panels: Similar to Fig.1, but for simulated ALMA continuum images with pure viscosity case (the left panel) and pure wind case (the right panel). Lower panels: The sky model processed by RADMC-3D for reference cases.

5 Discussions

Compared multi-fluid simulations to the latest ALMA Band 6 and JVLA Q Band observations, we found that the accretion in DM Tau should stem from the combined action of MHD disk wind and viscosity. In terms of contribution to turbulence, MHD disk wind will dominate. Yet in shaping the morphology and substructures of the disk, the roles of MHD disk wind and viscosity are equally important, despite viscosity contributing only a minor part to the accretion rate.

For regions within  30 au of PPDs, this possibility exists. Through 2D MHD simulations in spherical polar coordinates with the Hall effect and ambipolar diffusion, Bai & Stone (2017) found that when a poloidal field is aligned with disk rotation, the contributions of viscosity and wind to angular momentum transfer in the midplane region are comparable.

The substructures within DM Tau should also depend on planet-disk interaction and do not require planets of larger than Jupiter. Observations suggest that within the 20 au region of DM Tau (Flaherty et al., 2020), there exists sufficient material for the formation of planets up to 100 times the mass of Earth. Our study believes that planets of this mass are already capable of producing substructures consistent with observations.

Without considering the complexities of multi-planet systems (Xu & Wang, 2024), the currently known models show that only the MHD disk wind model can simultaneously explain the asymmetry, gap opening, and high accretion rate observed in DM Tau. Other models, such as the gravitational instability model (Xu, 2022; Xu et al., 2023), can reproduce asymmetric clumps and high accretion rates but fail to account for the gap opening issue. And the dusty Rossby wave instability (Liu & Bai, 2023) can explain the asymmetry and gap opening in the ring but fails to account for the high accretion rate observed in DM Tau.

5.1 Gap Opening

Traditionally, the required minimum gap opening planet to stellar mass ratio is very sensitive to disk aspect ratio $H/r$ and viscosity $\alpha_{\rm v}$. Crida et al. (2006) defined a parameter $\mathcal{P}$

$$\mathcal{P}=\frac{3H}{4r_{\rm H}}+\frac{50\alpha_{\rm v}}{q}\left({\frac{H}{r}}\right)^{2}$$ (10)

to quantify the relationship between gap opening and planet to stellar mass ratio. At $\mathcal{P}\textless 1$ the planet migrates in the type II regime, with the potential to open a gap with a factor of 10 depression in the gas surface density. Here $r_{\rm H}=r(M_{\rm p}/3M_{\star})^{1/3}$ is the Hill radius.

When the contribution of MHD disk wind to turbulence is taken into account, Elbakyan et al. (2022) derived Eq.10 to

$$\mathcal{P}=\frac{3H}{4r_{\rm H}}+\frac{50\alpha_{\rm v}}{q}\left({\frac{H}{r}}\right)^{2}+\frac{70\alpha_{\rm dw}^{3/2}}{q}\left({\frac{H}{r}}\right)^{3}.$$ (11)

Obviously, when $\alpha_{\rm dw}\gg\alpha_{\rm v}$ , the value of $\mathcal{P}$ will be dominated by $\alpha_{\rm dw}$. This also means that the extent of gap opening at this time will be determined by the strength of the MHD disk wind.

5.2 Asymmetry in Ring

According to the toy model which given by Wu et al. (2023b), the asymmetric structures on the ring should originate from vortices. One popular hypothesis is that vortices arise due to the Rossby Wave Instability (RWI, Lovelace et al., 1999; Li et al., 2001), triggered by sharp gradients in the gas surface density at the boundaries of planetary gaps or other perturbation. These vortices can capture dust into clumps, creating asymmetries observable in brightness maps. Typically, sharper gas surface density gradients promote RWI development. Although both wind and viscosity work to reduce the contrast in gas surface density across a planetary gap (Elbakyan et al., 2022), their impact on sustaining vortices appears to be contradictory.

The MHD disk wind transfers angular momentum at the disk’s mid-plane in a manner akin to inward radial drift, Wu et al. (2023b) believed that during the process of angular momentum transfer, the wind pushes the fluid inward as much as possible. Due to planet-disk interactions, the fluid then accumulates just behind the planet, creating a vortensity gradient (as shown in Fig.4, see the orange (cyan) line for 1 (3) mm dust). According to RWI, this process makes it easier for vortices to form on the ring compared to a viscosity-dominated disk, manifesting in observations as asymmetric structures with high signal (as shown in the pure wind case in Fig.3). This is also why in the 7mm JVLA observations, the asymmetric structures appear more significant compared to the observations in ALMA Band 6.

It is important to note that, unlike the traditional vortices which lack bright dust rings (i.e., the pure-wind case in Fig.3), our fiducial model still retains complete dust ring structures. Upon incorporating appropriate viscosity, the interaction between wind and viscosity tends to break down large clumps into smaller, numerous clumps distributed across the ring. This accounts for the source of asymmetric structures we observe in Fig.1. Additionally, appropriate viscosity and a moderate planet mass generally maintain structural symmetry and stability in the dust rings and gaps. This is something that massive-planet models are unable to achieve (Chen et al., 2024).

Similar vortex structures also seen in 3D MHD simulations (Aoyama & Bai, 2023; Hu et al., 2024), which suggest that the tendency for MHD disk wind to facilitate the formation of vortices/asymmetric structures is not merely an artifact of our wind prescription or the 2D geometric structure used in our simulations.

Figure 4: The averaged radial profiles for each fluid surface density in fiducial model with a 53 $M_{\oplus}$ planet at 20 au from the central star.

5.3 Caveats and Outlook

Unlike real 3D MHD simulations, our simplified 2D MHD disk wind prescription allows for longer run times and a wider range of parameter exploration, and the dust dynamics can be considered in simulations. However, this also means our hydrodynamic simulations are limited by the 2D geometry, with the MHD wind only added as a torque term driving radial gas inflow. This approximation cannot capture feedback effects on the initial magnetic field structure or wind profiles during gap opening.

A recent study conducted a detailed examination using three-dimensional magnetohydrodynamic simulations of star-disk interactions (Aoyama & Bai, 2023). They explored how MHD winds transport gas through planetary gaps and examined the migration torques exerted on the embedded planets by the disk. The conclusion was that due to magnetic flux concentration near the planet, the wind torque within the gap is increased by 3-5 times. This implies that $\alpha_{\rm dw}$ should be stronger in the gap region, naturally resulting in stronger turbulence. However, while these 3D MHD simulations are complex and capture additional physical features, they are computationally expensive and thus cannot run for the long timescales related to ALMA disk ages. They also cannot incorporate dust dynamics to generate simulated observational images for comparison with ALMA observations. This means that in the future, we should look to optimize our 2D MHD disk wind prescription based on the latest 3D simulation results for more refined outcomes.

Additionally, despite many observations of DM Tau, some key parameters like disk mass and scale height remain unclear, yet these parameters are crucial for observational structures. To better model and analyze the evolution of DM Tau theoretically, more critical parameters from observations are needed for better constraints.

It’s important to note that our simulations do not cover the inner disk within 4 au of DM Tau, Hashimoto et al. (2021) have found substructures and the potential presence of a type I planet. If there is indeed high turbulence in DM Tau, it’s necessary to consider possible planet migration (McNally et al., 2020; Kimmig et al., 2020; Ziampras et al., 2024; Turpin & Nelson, 2024) and strong magnetic field problems, especially in the inner region. Wu et al. (2024a) indicates that high turbulence can cause chaotic migration of type I planets, and studies show that planet migration can also affect the formation of substructures (Wu et al., 2023a). Moreover, under appropriate magnetic field conditions, the dust dynamics in the inner region might differ from the usual scenario (Krapp et al., 2018; Wu et al., 2024b), which future theoretical research needs to consider.

Furthermore, recent $[$O I$]$ line spectral mapping observations have identified evidence of wind in TW Hya (Fang et al., 2023), PDS 70 (Campbell-White et al., 2023) and T Cha (Bajaj et al., 2024). Though there is still debate over whether it’s MHD disk wind or photoevaporative wind (Lin et al., 2024; Sellek et al., 2024). However, our study suggests DM Tau is undoubtedly a source suitable for future related observations.

6 Conclusions

In this work, we present the results of multi-fluid hydrodynamic simulations that incorporate the 2D prescription of MHD disk wind. Through planet-disk interaction and wind-viscosity interaction, we successfully fit some observational characteristics of DM Tau in different wavelength (ALMA Band 6 and JVLA Q Band), thereby investigating the source of high accretion rate in this system. Our main findings are as follows:

• We find that the high accretion rate in DM Tau could arise from the combined action of MHD disk wind and viscosity, where MHD disk wind should account for about $97\%$ of the contribution, with the remaining approximately $3\%$ coming from viscosity. The high-contrast clumps that produce asymmetric structures in the bright dust ring around 21 au in DM Tau, as well as gap opening, could be mainly result from the interaction between the wind and a 53 $M_{\oplus}$ planet, affecting the distribution of dust. While the maintenance of the ring and gap morphology comes from appropriate viscosity. Based on suggestions from MHD simulations (Bai & Stone, 2017), this may imply that the potential Hall effect in DM Tau cannot be ignored.

• Our study indicates that a planet fitting the material mass of the DM Tau system can indeed exist within the 4-21 au gap in DM Tau, and the formation of a gap consistent with observational structures does not require the presence of a gas giant planet as large as 1 $M_{J}$ or even 10 $M_{J}$.

• Our research also supports recent conclusions from sub-cm wavelength observations of DM Tau (Liu et al., 2024), namely that at wavelengths $\textless$ 1.3 mm, the flux density is dominated by optically thick dust thermal emission, while at sub-mm and even longer wavelengths, the flux density should difference.

Although Flaherty et al. (2020) calculated an alpha value of $\alpha_{\rm v}\sim 0.078$ by detecting nonthermal gas motions in DM Tau, indicating that DM Tau is a highly turbulent system, our results do not actually conflict with this. This is because the turbulence level measured by Flaherty et al. (2020) comes from regions beyond 100 au from DM Tau and above the midplane ($\sim 2h$), whereas our focus is on the midplane and the region within $\sim$ 20 au. Moreover, MHD simulations also suggest that the turbulence level at the midplane should be lower than that in higher regions (Simon et al., 2015).

Overall, as a system with such a high accretion rate and varied substructures, DM Tau warrants further research effort in the future, especially with the help of upcoming observational facilities like SKA and ngVLA in longer wavelength (Ricci et al., 2018; Ilee et al., 2020; Ricci et al., 2021; Wu et al., 2024c). These will complement the existing capabilities of ALMA and JVLA in terms of wavelength or resolution, thus better constraining key parameters of the disk (e.g. disk mass and temperature). Meanwhile, MHD disk wind, as an important addition to existing accretion and planet formation theories, deserves to be applied across a broader range of systems.

Aknowledgement

Y.W. thanks the anonymous referee for highly constructive suggestions, which significantly enhanced the quality of the manuscript. Y.W. thanks Jun Hashimoto and Hauyu Baobab Liu for their discussions and share observational data. Y.W. also thanks Xue-Ning Bai, Kan Chen, Yi-Xian Chen and Wenrui Xu for very useful suggestions.

This research used DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk).

Data availability

The data obtained in our simulations can be made available on reasonable request to the corresponding author.

Appendix A Radiative Transfer Calculations

During the process of performing radiative transfer calculations on the results of hydrodynamic simulations using RADMC-3D, we first model the vertical expansion of the 2D dust density map into 3D by adopting follow Gaussian distribution for the dust densities along the vertical axis with 40 grids:

$$\rho_{\rm d}(r,z)=\Sigma(r)\exp(-z^{2}/2h_{\rm d}^{2}))/\sqrt{2\pi}.$$ (A1)

To process the dust scale height, $h_{d}$, a common method is using the formula

$$h_{d}\approx\sqrt{D_{\rm g}/(\Omega{\rm St})},$$ (A2)

as detailed by Dubrulle et al. (1995); Youdin & Lithwick (2007), here $D_{\rm g}$ is the vertical turbulent diffusion coefficient at the midplane.

As shown in Eqs. 1 - 4, since our code is multi-fluid, the wind torque does not directly interact with the dust in the Navier-Stokes equations. That means our wind torque does it affect the diffusion evolution of the dust itself in simulations. Therefore, we have ignored all contributions of MHD winds to dust diffusion.

The coefficient $D_{\rm g}$, is often equated to the viscosity $\nu=\alpha_{\rm v}c_{s}h_{\rm g}$ under the assumption of isotropic turbulence (Birnstiel et al., 2016; Baruteau et al., 2021). Yet, findings from MHD simulations introduce a layer of complexity (Zhu et al., 2015; Bai et al., 2016), highlighting the potential inaccuracies in this assumption. The presence of wind can elevate dust to higher altitudes, however, this typically applies to micron-sized dust with very small Stokes numbers, which have $h_{d}$ intrinsically reaching the wind base at $\sim h$ (Armitage et al., 2013), rather than to millimeter-sized dust.

Our analysis, particularly focused on disks observed face-on, has shown that variations in the vertical diffusion coefficient minimally impact the observed results. Thus, for simplicity and to concentrate on analyzing data from relatively face-on perspectives, we adopt a uniform dust scale height $h_{\rm d}=0.1~{}h_{\rm g}$ across all dust components, aligning with the approach of Miranda et al. (2017). This adjustment equates to a diffusivity coefficient $D_{\rm g}=10^{-4}h_{g}^{2}\Omega$ for pebbles with a Stokes number of 0.01. Simulations incorporating both ideal and non-ideal MHD effects (Bai & Stone, 2014; Baehr & Zhu, 2021; Zhou et al., 2022) suggesting possible variations in dust diffusivity along vertical and radial dimensions.

However, some works indicate that the vertical setting of particles can affect several observables, such as the total millimeter flux and shape of the substructures, even if the disk inclination is low (e.g., Pinilla et al., 2021). This implies that we should explore this aspect further in future research.

We apply the DSHARP opacities version without ice, it assume our dust species comprised of $41.1\%$ astro-silicates, 49.6% CHON organics and $9.3\%$ FeS, with internal density $\rho_{\rm d}=2.11~{}g/cm^{3}$ consistent with our simulations.