© ESO 2022

1 Introduction

The limited lifetime of a few million years (Haisch et al. 2001; Mamajek et al. 2004; Ribas et al. 2015) of protoplanetary disks has sparked interest in various mechanisms that transport and extract angular momentum from the disk. One possibility to drive accretion onto the central star is viscous turbulence, usually referred to as α-viscosity in the context of disks (Shakura & Sunyaev 1973). Turbulence might be generated by a multitude of physical instabilities in gas and dust. The well studied magnetorotational instability (MRI; Balbus & Hawley 1991) present in an ionized, differentially rotating disk threaded by a weak magnetic field can lead to values of α roughly ranging from 10⁻³ to 10⁻² and it was studied in a local shearing box (Fromang & Papaloizou 2007; Fromang et al. 2007; Bai 2011) as well as in global simulations (Dzyurkevich et al. 2010; Flock et al. 2011).

The MRI can only operate if the ionization fraction of gas is sufficiently high, which is not necessarily the case in the midplane of the disk (Igea & Glassgold 1999), where nonideal magnetohydrodynamic (MHD) effects such as ohmic diffusion dominate. These so-called dead zones prevent magnetically driven accretion in the midplane, whereas layered accretion in the ionized upper part of the disk atmosphere is still possible (Gammie 1996). In the more dilute outer regions of the disk, the ionization fraction rises again and the edge between the inner dead zone and the outer MRI-active part might cause jumps in the density profile and therewith enable dust or pebble traps (Lyra et al. 2009; Dzyurkevich et al. 2010, 2013). Purely hydrodynamic instabilities such as the global baroclinic instability (Klahr & Bodenheimer 2003) and the convective overstability (Klahr & Hubbard 2014), on the other hand, are alternative turbulence driving mechanisms that can also operate in the weakly ionized parts of the disk. In the presence of short gas cooling times, the vertical shear instability (VSI) causes vertical oscillations which are able to transport small dust grains towards the upper layers of the disk (Stoll & Kley 2014, 2016; Flock et al. 2020; Blanco et al. 2021). Vortices triggered by the VSI in combination with the Rossby wave instability (Lovelace et al. 1999) may act as long-lived dust traps (Manger & Klahr 2018; Pfeil & Klahr 2021).

Since recent observation hint towards small viscosity and thus small values of α on the order of 10⁻⁴ to 10⁻³ (Flaherty et al. 2015, 2017), the effect of disk winds might be a possible explanation for the finite disk lifetimes. In the context of photoevaporation, ionizing radiation either from the central star or the from the stellar environment heats the surface layers of the disk and drives winds due to the thermal pressure gradient. Depending on the dominant type of radiation, the wind mass loss rates are on the order of 10⁻¹⁰ M_⊙ yr⁻¹ with EUV radiation (Hollenbach et al. 1994; Font et al. 2004), and 10⁻⁸ M_⊙ yr⁻¹ to 10⁻⁷ M_⊙ yr⁻¹ with FUV radiation (Adams et al. 2004; Gorti & Hollenbach 2009) and X-ray radiation (Ercolano et al. 2009, 2021; Owen et al. 2010, 2012; Picogna et al. 2019, 2021).

An alternative way to produce strong disk winds would be by considering the effect of magnetic fields. In the presence of a global magnetic field the differential shearing motion of the disk generates a magnetic pressure gradient, accelerating matter and thereby driving the outflow. In the strong field limit, the gas follows the magnetic field lines and is magneto-centrifugally accelerated, as described in Blandford & Payne (1982). Contrary to photoevaporation, magnetically driven winds extract angular momentum from the disk and drive accretion. The accretion flow has to pass by the mostly vertically aligned magnetic field in the disk which only allows stationary solutions if the coupling between gas and magnetic field is weak, thus if significant diffusive nonideal MHD effects are present (Wardle & Koenigl 1993; Ferreira & Pelletier 1993, 1995). The magnetic launching mechanism is considered to be responsible for high-velocity jet outflows and has been studied extensively in global resistive simulations (Casse & Keppens 2002; Zanni et al. 2007; Tzeferacos et al. 2009; Sheikhnezami et al. 2012; Fendt & Sheikhnezami 2013). Stepanovs & Fendt (2014, 2016) determined the transition between magnetocentrifugally driven and magnetic pressure-driven jets at β ≈ 100 where β is the ratio of thermal over magnetic pressure.

With the growing interest in winds in the context of protoplanetary disks, numerous studies have been performed over the lastdecade, first in the form of local shearing box simulations (Suzuki & Inutsuka 2009; Suzuki et al. 2010; Fromang et al. 2013; Bai & Stone 2013; Bai 2013). Global nonideal MHD simulations have been carried out by Gressel et al. (2015) including ohmic and ambipolar diffusion finding MRI-channel modes to develop in the upper layers of the disk. They furthermore stress the importance of ambipolar diffusion enabling a more laminar wind flow and decreased mass loss rates. Béthune et al. (2017) find asymmetric wind flows, nonaccreting configurations and self-organizing structures such as rings. The results of Suriano et al. (2018) corroborate the ones by Béthune et al. (2017) such that rings emerge in ambipolar diffusion dominated disks by reconnection events. A combination of both photoevaporation and magnetically driven winds, also called magnetothermal winds were studied by Wang et al. (2019) and Rodenkirch et al. (2020) in the regime of magnetic pressure supported winds. With increasing magnetic field strengths, photoevaporative winds smoothly transition into magnetic winds driven by the magnetic pressure gradient in the upper layers of the disk (Rodenkirch et al. 2020). In contrast to the results of Bai (2017); Suriano et al. (2018); Wang et al. (2019) and Rodenkirch et al. (2020), simulations in Gressel et al. (2020) demonstrate magneto-centrifugal winds mainly driven by magnetic tension forces.

Dust entrainment in protoplanetary disk winds has been mostly examined in photoevaporative winds. Low gas densities in the wind limit the maximum grain size to a few μm. For EUV-photoevaporation Owen et al. (2011) determined a grain size limit of ≈ 2.2 μm. Hutchison et al. (2016b) performed small scale smoothed particle hydrodynamics simulations of a dusty EUV-driven wind and formulated an analytical expression for the maximum entrainable dust grain size. They conclude that typically the limiting size is less than 4 μm for typical T Tauri stars and they also argue that dust settling may lower this limit further. Franz et al. (2020) simulated dust entrainment in postprocessed X-ray and extreme ultraviolet (XEUV)-wind flows of Picogna et al. (2019) and found a larger maximum dust size of 11 μm. Dust entrainment in magnetically driven winds was tested by vertical 1D models in the work of Miyake et al. (2016), concluding that grains in the range of 25 μm to 45 μm can float up to 4 gas scale heights above the disk mid plane. Giacalone et al. (2019) used a semianalytical MHD-wind description to compute the dust transport and report a maximum dust size of 1 μm. Considering the effect of radiation pressure, submicron dust grains can be blown away from the disk and support clearing out remaining dust (Owen & Kollmeier 2019).

In this paper we present fully dynamic, multifluid, global, azimuthally symmetric, nonideal MHD simulations of protoplanetary disk winds including XEUV-heating to model photoevaporative flows. We aim to study dust entrainment in magnetically and thermally driven disk winds in a dynamic fashion that also considers the vertical dust distribution due to turbulent diffusion of dust grains. We thus extend the previous study in Rodenkirch et al. (2020) by including dust as a pressureless fluid. We furthermore postprocess the dust density maps to examine observational signatures of such dusty winds.

In Sect. 2, we describe the theoretical and computational concepts of the model. The results are presented in Sect. 3 and Sect. 4 discusses limitations of the model as well as the comparison to other studies. Finally, in Sect. 5 we give concluding remarks.

2 Model

The simulations were carried out with the FARGO3D code Benitez-Llambay & Masset (2016); Krapp et al. (2018); Benitez-Llambay et al. (2019) in a 2.5D axisymmetric setup using a spherical mesh with the coordinates (R, θ, φ). Throughout this work the cylindrical radius will be referred to as r, whereas the spherical radius will be denoted as $R=\sqrt{r^2+z^2}$. Observational inclinations, flow angles and viewing angles will be expressed in terms of latitude equivalent to θ ranging from 90 ° to −90 °.

2.1 Basic equations

The FARGO3D code solves the following set of equations: $${\partial }_{\text{t}}{\rho }_{\text{g}}+\nabla \cdot \left({\rho }_{\text{g}}{v}_{\text{g}}\right)=0,{\rho }_{\text{g}}\left({\partial }_{\text{t}}{v}_{\text{g}}+v_{\text{g}}\cdot \nabla v_{\text{g}}\right)=-\nabla P+\nabla \cdot \Pi -{\rho }_{\text{g}}\nabla \Phi +\frac{1}{4\pi }\left(\nabla \times B\right)\times B,{\partial }_{\text{t}}{\rho }_{\text{d},\text{i}}+\nabla \cdot \left({\rho }_{\text{d},\text{i}}{v}_{\text{d},\text{i}}+j_{\text{i}}\right)=0,{\rho }_{\text{d},\text{i}}\left[{\partial }_{\text{t}}{v}_{\text{d},\text{i}}+\left(v_{\text{d},\text{i}}\cdot \nabla \right)v_{\text{d},\text{i}}\right]=-{\rho }_{\text{d},\text{i}}\nabla \Phi +{\sum }_{\text{i}}{\rho }_{\text{i}}{\text{f}}_{\text{i}},$$

where Eqs. (1) and (4) describe the mass conservation of gas and dust with ρ_g and ρ_d,i being the gas and dust density, respectively. The suffix i denotes the corresponding dust species. Equations (2) and (5) represent the momentum equations of gas and dust. We note that we did not consider the dust backreaction onto the gas in Eq. (2) since the dust-to-gas mass ratio is assumed to be low. The symbol P refers to the thermal gas pressure, B the magnetic flux density vector, v the velocity vectors, $\Phi =-\frac{G{M}_{*}}{R}$ the gravitational potential and f_i the drag forces between dust and gas.

Turbulent mixing of dust grains caused by the underlying turbulent viscosity in the gas was considered by defining the diffusion flux j_i (Weber et al. 2019) $$j_{\text{i}}=-D_{\text{i}}\left({\rho }_{\text{g}}+{\rho }_{\text{d},\text{i}}\right)\nabla \left(\frac{{\rho }_{\text{d},\text{i}}}{{\rho }_{\text{g}}+{\rho }_{\text{d},\text{i}}}\right),$$(6)

where the diffusion coefficient D_i depends on the kinematic viscosity ν and the Stokes number St (Youdin & Lithwick 2007): $$D_{\text{i}}=\nu \frac{1+{\text{St}}^2}{{\left(1+{\text{St}}^2\right)}^2}.$$(7)

The dimensionless Stokes number can be expressed in terms of the stopping time t_s as follows: $$\text{St}={\Omega }_{\text{K}}{t}_{\text{s}}={\Omega }_{\text{K}}\frac{{\rho }_{\text{mat}}{a}_{\text{i}}}{{\rho }_{\text{g}}{v}_{\text{th}}},$$(8)

with the Keplerian angular frequency ${\Omega }_{\text{K}}=\sqrt{G{M}_{*}/r^3}$, the material density of the dust grain ρ_math and the mean thermal velocity of the gas: $$v_{\text{th}}=\sqrt{\frac{8k_{\text{B}}{T}_{\text{g}}}{\mu m_{\text{p}}\pi }}.$$(9)

The Boltzmann constant is expressed as k_B, the gas temperature as T, the mean molecular weight as μ and the proton mass as m_p. The viscosity stress tensor is given by $$\Pi =\rho \nu \left[\nabla v+{(\nabla v)}^T-\frac{2}{3}(\nabla \cdot v)I\right],$$(10)

where I is the unit tensor, and the energy equation can be stated as $$\frac{\partial e}{\partial t}+\nabla \cdot (ev)=-P\nabla \cdot v.$$(11)

The adiabatic equation of state is given by $$P=(\gamma -1)e,$$(12)

where γ = 5∕3 is the adiabatic index. The induction equation including ohmic and ambipolar diffusion reads as follows (Béthune et al. 2017): $$\frac{\partial B}{\partial t}-\nabla \times \left(v\times B-\frac{4\pi }{c}{\eta }_{\text{ohm}}\cdot J+{\eta }_{\text{am}}\cdot J\times b\times b\right)=0,$$(13)

where b = B∕|B| the normalized magnetic field vector, J the electric current density vector, η_ohm the ohmic diffusion coefficient and η_am the ambipolar diffusion coefficient.

2.2 Disk model

2.2.1 Initial conditions of the gas in the disk

The vertically integrated surface density profile Σ(r) was assumed to be a power law $\Sigma (r)={\Sigma }_0{\left(\frac{r}{r_0}\right)}^{-\text{p}}$ with the reference surface density Σ₀ at the characteristic radius r₀ of the system and the volumetric gas density takes the form $$\rho (r,z)={\Theta }_{\text{H}}(r-r_{\text{c}})\frac{{\Sigma }_0}{\sqrt{2\pi }H(r)}{\left(\frac{r}{r_0}\right)}^{\text{p}}\text{exp}\left[\frac{1}{h^2}\left(\frac{r}{\sqrt{r^2+z^2}}-1\right)\right],$$(14)

where h = H∕r is the local aspect ratio and H = c_s∕Ω_K the gas scale height with the local isothermal sound speed $c_{\text{s}}=\sqrt{k_{\text{B}}T(r)/(\mu m_{\text{p}})}$. Within a certain cutoff radius r_c, we set the density to zero with a sharp transition and we thus mimiced a disk with an inner cavity of size r_c. The transition is initially taken care of by the Heaviside function Θ_H(r − r_c) in Eq. (14). The sharp transition quickly softens and relaxes to the shape of a smooth inner rim during the simulation run after a few orbits. The cavity in the model was applied for both gas and dust.

In a simple blackbody disk being irradiated by stellar light with a grazing angle ϕ onto its surface, the temperature scales with radius as (Chiang & Goldreich 1997) $$T_{\text{g}}(r)={\left(\frac{\varphi }{2}\right)}^{\frac{1}{4}}{\left(\frac{R_{*}}{r}\right)}^{\frac{1}{2}}{T}_{*}={\left(\frac{\varphi L_{*}}{8\pi r^2{\sigma }_{\text{SB}}}\right)}^{1/4},$$(15)

where R_*, T_* and L_* are the stellar radius, temperature and luminosity, respectively. The Stefan-Boltzmann constant is denoted by σ_SB. Thus, we assumed a radial power law profile for the temperature with T(r) = T₀ r^−q and the gas scale height scales with radius as $H(r)\propto \sqrt{T(r)}\propto r^{(3-q)/2}$.

The initial values of the azimuthal gas velocity are slightly subkeplerian (Nelson et al. 2013): $$v_{\varphi }(r,z)={\Omega }_{\text{K}}r\sqrt{h^2(p+q)+(q+1)-\frac{qr}{\sqrt{r^2+z^2}}}.$$(16)

We assumed no initial velocities in the radial and polar direction, that is v_r = v_θ = 0.

2.2.2 Coronal gas structure

The hydrostatically stable coronal gas structure with density ρ_c (R) and pressure P_c(R) serves as a density floor throughout the simulation and prevents excessively low initial densities which prevents numerical instabilities. We used the same prescription as used in the previous models in Rodenkirch et al. (2020) based on Sheikhnezami et al. (2012).

The density was set to $${\rho }_c(R)={\rho }_{c0}(r_i){\left(\frac{r_i}{R}\right)}^{\frac{1}{\gamma -1}},$$(17)

where r_i is the radius at the inner boundary of the simulation domain and γ the adiabatic index. The corresponding pressure profile is given by: $$P_c(R)={\rho }_{c0}(r_i)\frac{\gamma -1}{\gamma }\frac{G{M}_{\odot }}{r_i}{\left(\frac{r_i}{R}\right)}^{\frac{\gamma }{\gamma -1}}.$$(18)

The density at the inner boundary ρ_c0(r_i) is related to the disk’s density at the midplane and r_i with the density contrast $\delta ={\rho }_{c0}(r_i)/\left[\Sigma (r_{\text{i}})/\left(\sqrt{2\pi }H\right)\right]$. Depending on Σ₀, we set the density contrast to values ranging from δ = 10⁻⁷ to 5 × 10⁻⁷. In absence of X-ray/EUV heating, the gas temperature was directly set to the coronal temperature.

2.2.3 Dust structure

Vertical dust settling towards the midplane becomes more effective in the upper layers of the disk and the overall vertical thickness of the dust layer is over-estimated in models of a dust scale height uniquely depending on midplane values as in Dubrulle et al. (1995) (Dullemond & Dominik 2004). We therefore used the following vertical dust distribution $${\rho }_{\text{d}}(r,z)={ϵ}_{\text{dg}}\rho (r,z)\left\{\frac{\text{St}(r,z=0)}{\alpha }\left[\text{exp}\left(\frac{z^2}{2H^2}\right)-1\right]-\frac{z^2}{2H^2}\right\}$$(19)

given in Fromang & Nelson (2009). The dust-to-gas mass ratio of the dust fluids is denoted by ϵ_dg and the dimensionless diffusion coefficient follows the α-viscosity prescription ν = αc_sH (Shakura & Sunyaev 1973) with the kinematic viscosity ν. We assumed a constant value for the Schmidt number Sc = 1.

The initial azimuthal velocities were set to $v_{\text{d},\varphi }={\Omega }_{\text{K}}\sqrt{r^2+z^2}$. Similar to the gas, the initial radial and polar velocities of the dust fluids were set to zero.

2.2.4 Ionization rate

In all of the following models X-ray radiation was considered including both a radial and scattered component. The parametrization was based on the prescription of Igea & Glassgold (1999) and Bai & Goodman (2009): $${\zeta }_{\text{X}}=\text{\&}{\zeta }_{\text{X},\text{sca}}\left[\exp {\left(-\frac{{\Sigma }_{\text{top}}}{{\Sigma }_1}\right)}^{0.65}+\exp {\left(-\frac{{\Sigma }_{\text{bot}}}{{\Sigma }_1}\right)}^{0.65}\right]{\left(\frac{R}{\text{au}}\right)}^{-2}\text{nn\&}+{\zeta }_{\text{X},\text{rad}}\exp {\left(\frac{{\Sigma }_{\text{rad}}}{{\Sigma }_2}\right)}^{0.4}{\left(\frac{R}{\text{au}}\right)}^{-2}.$$(20)

Both coefficients were multiplied by a factor of 20 compared to the cited values to model a fiducial X-ray luminosity of 2 × 10³⁰ erg s⁻¹ and were set to ζ_X,rad = 1.2 × 10⁻¹⁰ s⁻¹ and ζ_X,sca = 2 × 10⁻¹⁴ s⁻¹. The critical column densities accounted to Σ₁ = 7 × 10²³cm⁻² and Σ₂ = 1.5 × 10²¹ cm⁻².

The luminosity increase was chosen to be compatible with the photoevaporation model of Ercolano et al. (2008); Picogna et al. (2019) where the luminosity L_X = 2 × 10³⁰ erg s⁻¹ was used as the fiducial value. Taking into account the mass-luminosity relation presented in Flaischlen et al. (2021), Fig. 1 based on Getman et al. (2005), a luminosity of L_X = 1 × 10²⁹ erg s⁻¹ would correspond in average to a M_*≈ 0.15 M_⊙ star. The increased value is thus better suited for the T Tauri star with M_* = 0.7 M_⊙ assumed in our model.

Additionally, the cosmic ray ionization rate ζ_cr was modeled by the parametrization of Umebayashi & Nakano (2009) where a characteristic column density of 96 g cm⁻² was chosen. The symbols Σ_bot and Σ_top refer to the column densities integrated along the polar direction, starting from the lower boundary (suffix bot)] and the upper boundary (suffix top). The radial column densities Σ_rad were computed starting from the radial inner boundary. $${\zeta }_{\text{cr}}=\text{\&}5\times {10}^{-18} {\text{s}}^{-1}\exp \left(-\frac{{\Sigma }_{\text{top}}}{{\Sigma }_{\text{cr}}}\right){\left[1+{\left(\frac{{\Sigma }_{\text{top}}}{{\Sigma }_{\text{cr}}}\right)}^{3/4}\right]}^{-4/3}\text{nn\&}+\exp \left(-\frac{{\Sigma }_{\text{bot}}}{{\Sigma }_{\text{cr}}}\right){\left[1+{\left(\frac{{\Sigma }_{\text{bot}}}{{\Sigma }_{\text{cr}}}\right)}^{3/4}\right]}^{-4/3}.$$(21)

To account for short-lived radio nuclides in the disk, an ionization rate of ζ_nuc = 7 × 10⁻¹⁹ s⁻¹ caused by the ²⁶Al isotope (Umebayashi & Nakano 2009) acts as an ionization rate floor value. The resulting total ionization rate then simply sums up to ζ = ζ_X + ζ_cr + ζ_nuc.

2.2.5 Ionization fraction

We computed the ionization fraction of the gas via a semianalytical model for fractal dust aggregates elaborated by Okuzumi (2009). The main assumption here is an approximately Gaussian charge distribution in the dust grains. This approximation is reasonably well fulfilled if the dust aggregate consists of at least 400 monomers Dzyurkevich et al. (2013). In this approach, the ionization equilibrium was described by the following equations Okuzumi (2009): $$\frac{1}{1+\Gamma }-\left[\frac{s_{\text{i}}{u}_{\text{i}}}{s_{\text{e}}{u}_{\text{e}}}\text{exp}(\Gamma )+\frac{1}{\Theta }\frac{\Gamma g(\Gamma )}{\sqrt{1+2g(\Gamma )}-1}\right]=0,$$(22)

with the parameter Θ: $$\Theta =\frac{\zeta n_{\text{g}}{e}^2}{s_{\text{i}}{u}_{\text{i}}\overline{\sigma }\overline{a}{n}_{\text{d}}^2{k}_{\text{B}}T}.$$(23)

The function g(Γ) is evaluated as: $$g(\Gamma )=\frac{2\tilde{\alpha }\zeta n_{\text{g}}}{s_{\text{i}}{u}_{\text{i}}{s}_{\text{e}}{u}_{\text{e}}{(\overline{\sigma }{n}_{\text{d}})}^2}\frac{\text{exp}(\Gamma )}{1+\Gamma }.$$(24)

After solving Eq. (22) numerically for Γ with a simple Newton-Raphson scheme, finally the ion and electron densities n_i and n_e can be calculated: $$n_{\text{i}}=\frac{\zeta n_{\text{g}}}{s_{\text{i}}{u}_{\text{i}}\overline{\sigma }{n}_{\text{d}}}\frac{\sqrt{1+2g(\Gamma )}-1}{(1+\Gamma )g(\Gamma )},$$(25) $$n_{\text{e}}=\frac{\zeta n_{\text{g}}}{s_{\text{e}}{u}_{\text{e}}\overline{\sigma }{n}_{\text{d}}}\frac{\sqrt{1+2g(\Gamma )}-1}{\text{exp}(-\Gamma )g(\Gamma )}.$$(26)

The stickingcoefficients s_i and s_e were set to 1 and 0.3 respectively, as described in Okuzumi (2009). The thermal ion and electron velocities u_i,e follow: $$u_{\text{i},\text{e}}=\sqrt{\frac{8k_{\text{B}}T}{\pi m_{\text{i},\text{e}}}},$$(27)

where m_e describes the electron mass and the ion mass is set to m_i = 24 m_p, corresponding to magnesium.

With the assumption of fractal dust aggregate with a dimension of D = 2, the dust aggregate size becomes ã = a₀N^1∕2 with the number of monomers N of size a₀. Similarly, the geometric cross section $\tilde{\sigma }$ of these aggregates can be expressed as $\tilde{\sigma }\approx \pi a_0^{2}N$. The dust number density n_d can be calculated as $n_{\text{d}}={ϵ}_{\text{ion}}{n}_{\text{g}}/(\frac{4\pi }{3}{a}_0^{3}N)$ with the dust-to-gas ratio ϵ_ion used in the ionization model.

In the limit of a vanishing dust-to-gas mass ratio the ionization fraction x_e approaches $$x_{\text{e}}=\frac{n_{\text{e}}}{n_{\text{g}}}=\sqrt{\frac{\zeta }{\tilde{\alpha }n}},$$(28)

where ζ is the total ionization rate, $\tilde{\alpha }=3\times {10}^{-6}{\text{cm}}^3{\text{s}}^{-1}/\sqrt{T}$ the dissociative recombination rate coefficient and n the neutral number density of the molecular species (Oppenheimer & Dalgarno 1974; Ilgner & Nelson 2006).

In the same manner as described in Bai (2017), we increased the ionization fraction in the wind region as a proxy for the effect of FUV-radiation: $$x_{\text{fuv}}=2\times {10}^{-5}\text{exp}\left[-{\left(\frac{{\Sigma }_{\text{rad}}}{{\Sigma }_{\text{fuv}}}\right)}^4+\frac{0.3{\Sigma }_{\text{fuv}}}{{\Sigma }_{\text{rad}}+0.03{\Sigma }_{\text{fuv}}}\right].$$(29)

2.2.6 Magnetic field

The magnetic field was set such that the initial plasma beta β at the midplane is independent of the radius r. The poloidal field was initialized from the vector potential A_ϕ $$A_{\varphi }(r,\theta )=\sqrt{\frac{8\pi P}{\beta }}{\left(\frac{r}{r_0}\right)}^{\frac{2p+q+1}{4}}{\left[1+m\cdot \text{tan}{\left(\theta \right)}^{-2}\right]}^{-\frac{5}{8}}$$(30)

described in Bai (2017) and exhibits an outward-bent, hourglass shape depending on the parameter m where m = 0 would correspond to a completely vertical field. In order to avoid excessively small time steps in the low density regions close to the rotation axis, we limited the local Alfvén velocity $v_{\text{A}}=B/\sqrt{4\pi {\rho }_{\text{g}}}$ to 18v_K(r₀) by increasing the density in the cell accordingly, similar to Riols & Lesur (2018). This procedure does not affect the bulk of the wind flow.

2.2.7 Nonideal MHD

Using the electron fraction x_e, the local ohmic diffusion coefficient η_ohm can be written as (Blaes & Balbus 1994): $${\eta }_{\text{ohm}}=\frac{c^2{m}_n{m}_e}{4\pi e^2}\frac{{\langle \sigma v\rangle }_{\text{e}}}{x_{\text{e}}}\approx \frac{234T^{\frac{1}{2}}}{x_{\text{e}}}\frac{{\text{cm}}^3}{\text{s}},$$(31)

where m_n and m_e are the neutral and electron mass respectively. The electron-neutral collision frequency was set to ${\langle \sigma v\rangle }_{\text{e}}=8.28\times {10}^{-10}{T}^{-1/2}{\text{cm}}^3{\text{s}}^{-1}$ (Draine et al. 1983). The ambipolar diffusion coefficient η_am becomes: $${\eta }_{\text{am}}=\frac{B^2}{4\pi {\langle \sigma v\rangle }_{\text{i}}{n}_{\text{g}}^2{x}_{\text{i}}}.$$(32)

The ion-neutral collision rate was set to ${\langle \sigma v\rangle }_{\text{i}}=1.9\times {10}^{-9}{\text{cm}}^3{\text{s}}^{-1}$ (Draine et al. 1983). Since large diffusivities severely limited the simulation time step a super time stepping scheme was implemented and more details are given in Appendix C. We furthermore limited the magnetic Reynolds number for ambipolar diffusion Rm = c_s H∕η_am to Rm ≥ 0.05, similar to Gressel et al. (2015), confirming that this approximation does not significantly alter the simulation outcome.

2.2.8 Heating and cooling

We employed the same temperature presciption as described in Rodenkirch et al. (2020) based on Picogna et al. (2019) to mimic the outflow due to photoevaporation in the simulation runs modeling a warm (magneto-) photoevaporative wind. According to their model, XEUV heating is significant up to radial column densities of Σ_r,crit = 2.5 × 10²² cm⁻². In these regions, the temperatures T_photo were updated to the values given by Eq. (1) in Picogna et al. (2019). Beyond the critical column density Σ_r,crit the temperature was cooled to the original gas temperature T₀ = T_g.

In the case of the cold magnetic wind model b4c2np (see Table 1), the temperature in the corona was only slightly increased by 20% compared to the bulk temperature in the disk. We used a simple β-cooling recipe in the whole simulation domain to exponentially damp the temperatures to T₀ on a time scale β_cool $$T(t^{\text{n}}+\Delta t)=T(t^{\text{n}})+\left(T(t^{\text{n}})-T_0\right)\text{exp}\left(\frac{\Delta t}{{\beta }_{\text{cool}}}\right).$$(33)

2.3 Boundary conditions

We employed simple symmetric boundary conditions for most of the variables and boundaries where the values of the active cells are copied to the ghost cells. At the inner radial boundary, the radial velocity followed an outflow boundary condition where v_r = 0 in the ghost zone if v_r > 0 in the active cell. In the contrary case, v_r was copied to the ghost cell. The same procedure was applied to the outer radial boundary with the opposite sign. With this formulation mass influx from outside of the domain was avoided.

At the θ_min and θ_max boundaries v_ϕ, v_θ, B_ϕ, B_θ and the EMF in radial direction change sign in the ghost zone to mimic a polar boundary. In MHD disk models the radial inner boundary condition can lead to spurious effects in the simulation domain. The same applies to the outer radial boundary where, if ill defined, artificial collimation affects the flow (Ustyugova et al. 1999). To avoid these difficulties we added a zone with additional artificial ohmic diffusion, so that the coupling between the gas and the magnetic field is weak enough that spurious effects were damped and the simulation domain remained largely unaffected. A similar approach was used at the inner boundary in Cui & Bai (2021) to stabilize the inner region of the domain.

Dimensionless values of the artificial ohmic resistivity ware η_ohm,in = 2 × 10⁻⁴ and η_ohm,out = 1 × 10¹ at the inner and outer damping zone, respectively. The extent of the damping zones is 10% of the inner and outer radius. With the choice of these values we found that the simulation domain and the timestep are largely unaffected by the boundary conditions. Adding an inner cavity with greatly reduced amount of gas and dust furthermore ensures negligible impact of the inner boundary condition on the wind flow.

2.4 Radiative transfer model and post-processing

We post-processed the dust density outputs with the radiative transfer code RADMC-3D (Dullemond et al. 2012). In order to obtain three dimensional dust density distributions the axisymmetric simulation outputs were extended in the azimuthal direction by repeating the densities in 300 cells in ϕ. For the thermal Monte-Carlo step and the image reconstruction with scattered light n_phot = 10⁶ and n_{phot_scat} = 10⁸ photon packages were used, respectively.

The thermal Monte-Carlo calculation was based on a modified version of the Bjorkman–Wood method (Bjorkman & Wood 2001) and the resulting temperatures were used in the synthetic images instead of the temperature profile prescribed in the MHD simulations. Every synthetic observation was created with only one dust species since the limited number of fluid would lead to artificial discrete step-like features in the image. Only considering one dust species at a time may lead to unrealistic temperature profiles but after checking the output of the thermal calculation of RADMC-3D the temperatures agreed with the dust temperature estimated computed in Appendix A in the relevant regions. We further emphasize that the synthetic images were produced in the H-band at a wavelength of 1.6 μm where scattered light effects dominate and thermal emission from the dust grains at these low temperatures is negligible.

The dust opacities were computed with optool (Dominik et al. 2021) and are similar to the DIANA opacities (Woitke et al. 2016) with a material mixture of 87% pyroxene, 13% carbon (by mass) but with an increased porosity of 64% to be consistent with a dust material density of ρ_mat = 1 g cm⁻³ in the simulations. Concerning the dust grain sizes, we computed the opacities for a₁ = 0.1 μm, a₂ = 1 μm and a₃ = 10 μm corresponding to the dust sizes in the simulations. Anisotropic scattering effects were included by using the Henyey-Greenstein phase function approximation (Henyey & Greenstein 1941). A classical T Tauri star (CTT) was assumed with a temperature of T_* = 4000 K and a radial extent of R_* = 2.55 R_⊙ at a distance of 150 pc which is in the range of typical properties of observed CTTs (Alcalá et al. 2021). The resulting images were convolved with a gaussian beam size of 75 mas to account for typical observational resolution limits using VLT/SPHERE (Beuzit et al. 2019; Avenhaus et al. 2018).

3 Results

Section 3.1 describes the relevant parameters and simulation runs. In the following subsections we present the simulation results starting with a description of the wind flow structure emerging from the disk in Sect. 3.2 and continue with the dust dynamics in Sect. 3.3. The effectiveness of dust entrainment and the outflow angle of the grains is addressed in Sect. 3.4, Sect. 3.5 and Sect. 3.6. Finally, in Sect. 3.7 we compare synthetic observations computed from the simulation output with existing observations.

3.1 Parameters and normalization

Table 1 summarizes the key parameters of the different simulation runs presented in this study. The first part of the label describes wether a magnetic wind (for instance b4 with β = 10⁴) or a photoevaporative wind without magnetic fields was modeled (with label ph). Subsequently, the letter “c” refers to the cavity location (such as c2 for a cavity at 2 au). The sublabels np, a3 and eps2 refer to “no photoevaporation”, α = 10⁻³ and ϵ_ion = 10⁻², respectively.

In the code all quantities were treated in dimensionless units that are scaled as follows: $$v_0=r_0{\Omega }_K(r_0),{\rho }_0=M_{\odot }/r_0^3,P_0={\rho }_0{v}_0^2,$$

where the length unit was set to r₀ = 5.2 au, the semi-major axis of Jupiter as a natural scale of stellar systems. The initial density profile was scaled to reach a surface density of Σ = 200 g cm⁻² at r = 1 au and decreases radially with power law slope of p = 1. The central star mass was set to 0.7 M_⊙ and the disk follows a moderate aspect ratio of $H/r_0=0.055{(r/r_0)}^{1/4}$ which corresponds to a flared disk with a temperature slope of q = 1∕2. These parameters correspond to a stellar luminosity of L_* = 1.5 L_⊙ with the assumption of a grazing angle of ϕ = 0.02. The luminosity is in line with the stellar parameters of T_* = 4000 K and R_* = 2.55 R_⊙ mentioned in Sect. 2.4.

A temperature floor of T_min = 10 K and a cooling time scale of β_cool = 10⁻⁴ was applied. The fiducial value of the inner cavity size was chosen to be r_c = 2 au. Further parameters are the adiabatic index of γ = 5∕3 corresponding to an ideal atomic gas and the mean molecular weight μ = 1.37125, consistent with the value chosen in Picogna et al. (2019) for the composition of the ISM (Hildebrand 1983).

In Sect. 2.2.5, the model necessitates an assumption of the dust content in the disk. In this work the ionization model was parametrized with a static value of the dust-to-gas ratio ϵ_ion = 10⁻³ and a monomer size of a₀ = 5 μm that constitute aggregates with 400 particles, resulting in a grain size of 100 μm. Increasing the dust-to-gas ratio by one order of magnitude did not alter the ionization fraction profile by a great amount, as described in Appendix B. In order to simplify the model, the dynamically simulated dust fluids were not coupled to the ionization model and the relative insensitivity to the chosen parameters in the midplane supports this approximation.

All simulation runs were carried out with three dust species of the sized a = 0.1 μm, 1 μm and 10 μm, intended to represent the three cases of full entrainment, slow decoupling, and fast decoupling from the gas in the wind. The material density of the grains was set to ρ_mat = 1 g cm⁻³.

The computational grid adopts a spherical symmetry and contains 500 logarithmically spaced cells in radial direction whereas 276 cells were used in the polar (θ) direction. Inorder to increase the vertical resolution closer to the disk mid plane, the polar cells were arranged in a stretched grid, with a cell size atthe polar boundaries that was twice as larger compared to the mid plane. With this approach a resolution of 7 cells per H in polar direction at r = r₀, z = 0 was reached. To avoid numerical issues at the polar boundaries the grid limits in this direction were chosen such that θ_min ≈ 87 ° and θ_max ≈−87 °. The radial domain spans from 1 au to 200 au as indicated in Table 1.

3.2 Gas flow structure

Figure 1 visualizes the gas wind flow for the simulation runs phc2, b5c2 and b4c2. The latter two include magnetic fields whereas the former only includes the photoevaporation recipe. Comparing both magnetic winds in Fig. 1 the magnetic field structure consists mostly of straight field lines for β =10⁴, whereas the weaker field shown in the central panel exposes a more irregular field structure towards the inner region of the disk. For this field strength and configuration, the wind is in the transition between a thermally dominated and a magnetically dominated wind. We note that these are simulation snapshots and not averaged wind flows.

In terms of wind total mass loss rates all of these three runs display similar results, as shown in Fig. 2. Generally, the wind converges to a steady-state flow after a time scale about 200 orbits at 1 au and the simulation time chosen for these runs is thus sufficient to explore relevant wind dynamics in the inner region.

The sharp cutoff of the disk at r = 2 au quickly relaxes towards a stable equilibrium after a few orbits. The choice of the initial shape of the inner rim is thus not relevant for the windflow. We observe accretion of gas through the inner cavity when a magnetic field is present. This effect, which can be best observed in the right panel of Fig. 1, is more pronounced for larger magnetic field strengths and is caused by wind-driven accretion.

The photoevaporation gas mass loss rates reach ≈2.83 × 10⁻⁸ M_⊙ yr⁻¹, in good agreement with Picogna et al. (2019). The mass loss rates in the simulation runs b5c2 and b4c2 are ≈ 3.13 × 10⁻⁸ M_⊙ yr⁻¹ and ≈ 3.03 × 10⁻⁸ M_⊙ yr⁻¹, respectively. In the case of apurely magnetically driven wind without any heating by ionizing radiation of model b4c2np the wind mass loss rates are similar with a value of 3.19 × 10⁻⁸ M_⊙ yr⁻¹. The mass loss rates were computed at 95% of the outer radial simulation boundary. The region within five gas scale heights of the disk was excluded from the mass loss rate to avoid the corruption of the results from circulation within the disk.

Fig. 1 Gas density maps of the simulations phc2, b5c2 and b4c2 after 500 orbits at 1 au. Magnetic field lines are represented by the white lines in the corresponding panels. The color map scaled logarithmically annotated by the exponent to a base of ten and represents the volume density of the gas in g cm−3.

Fig. 2 Total wind mass loss rates of gas and dust for the simulation runs phc2, b5c2, b4c2 and b4c2np. Results from 10 μm grains are excluded due to the negligible outflow and the lack of entrainment in the wind for all runs shown above.

3.3 Dust flow structure

In Fig. 3 the dust flows are shown for the grain sizes 0.1 μm, 1 μm and 10 μm. Generally, 10 μm sized grains are hardly entrained with the wind flow since the drag force is too weak to counter act the gravitational pull. 0.1 μm grains are well coupled to the gas flow and can be considered as tracers of the wind. Throughout the analysis we refer to the flow angle or inclination angle, meaning the angle enclosed between the dust flow direction and this disk mid plane axis at z = 0. The inclinations θ_d of the dust flows are smaller in terms of latitude for the thermally driven wind shown on the left-hand side compared to the magnetically driven wind on the right-hand side. The magnetic field is well coupled to the gas in the wind region where the ionization fraction x_e,i surpasses 10⁻⁴ and the field topology forces the wind flow to be more inclined compared to the thermally driven wind.

The same phenomenon also appears for 1 μm grains shown in the second row of Fig. 3. Due to the weaker coupling and the lack of pressure support the dust flow is less inclined with respect to the gaseous wind compared to smaller grains. The dust flow exhibits wave-like patterns in the case of β = 10⁵ and becomes more stable for larger magnetic field strengths. These wave structures coincide with the irregular magnetic field patterns displayed in Fig. 1. We furthermore point out that the dust flow consisting of 1 μm grains is efficiently launched in the inner part of the disk, whereas beyond several au the wind significantly weakens. The dusty wind cannot be sustained in the lower density regions further out and up in the disk since the Stokes number significantly increases.

An example of a dust flow in a cold magnetically driven wind is given in Fig. 4 for an initial β = 10⁴. The inclination with respect to the mid plane is smaller compared to the warm magnetothermal or photoevaporative winds. Grains with a size of 1 μm are entrained at a rather shallow inclination angle and similar to the warm winds, dust entrainment is negligible considering 10 μm dust particles. While material in the winds including photoevaporation is ejected with a velocity of ≈ 15 km s⁻¹ the outflow velocity of the cold magnetic wind in b4c2np only reaches values of ≈ 5 km s⁻¹. The decreased wind speed causes a smaller drag force between gas and thus lowers the dust entrainment efficiency. Given that the mass loss rate is approximately equal to the one of the heated winds, the gas and dust densities in the wind region are significantly higher.

In Fig. 5, the cumulative mass loss rate depending on the cylindrical radius of the foot point of the wind is shown. The wind streamlines were traced backwards from the cells close to the outer radial boundary and the foot point was registered when the streamlines crosses the surface at five gas pressure scale heights. Only streamlines with foot point r_f ≥ 3 au were considered since the wind flow structure and the limited resolution did not in every case ensure a successful construction of traceable streamlines close to the inner cavity. The mass loss rates were obtained from a velocity and density field averaged over 50 orbits at 1 au starting from 600 orbits.

Warm magnetic winds lead to more dust entrainment and a stronger magnetic field increases the dust mass loss rate significantly. The dashed lines in Fig. 5 mark the radius within from within 99% of the mass loss occurs, which illustrate less efficient dust entrainment for larger grains further out in the disk. Numerical values are 32.3 au, 34.5 au, 43 au and 88.8 au for the simulations phc2, b5c2, b4c2, b4c2a3 and 0.1 μm grains, respectively. Analogously, the limits for 1 μm grains are 14.6 au, 12.7 au, 14.6 au and 44.3 au. The general trend is that the dust grains are preferably entrained starting from the inner regions of the disk. The limiting radius decreases with increasing dust size and hence we expect a much larger ejection region of dust for the smallest grain size. Considering the results of the simulation b4c2a3 with an increased viscosity of one order of magnitude (α = 10⁻³) the limiting radius lies significantly further outward. Increased turbulent diffusion effectively transports dust grains towards the wind launching front at larger radii.

Going back to Fig. 2 the mass loss rate for the two relevant dust species stabilizes over a few hundred orbits in the case of a thermal wind. In the case of magnetically driven winds however, dust mass loss rates reach a quasi-steady state but are slowly increasing throughout the simulation. In the presence of a cold magnetic wind, the mass loss rates of larger micron-sized dust grains are lower compared to the ionized winds since the flow angle is rather shallow and the entrainment efficiency suffers from the decreased wind speed.

Fig. 3 Dust density maps for the simulations phc2, b5c2 and b4c2 after 500 orbits at 1 au. Dust velocity streamlines are annotated by the white arrowed lines. The upper, center and lower three panels visualize the flow of grains with size 0.1 μm, 1 μm and 10 μm, respectively.

Fig. 4 Dust density maps for the simulation b4c2np after 500 orbits at 1 au. Dust velocity streamlines are annotated by the white arrowed lines.

Fig. 5 Cumulative mass loss rates of gas and dust for the simulation runs phc2, b5c2, b4c2 and b4c2a3 starting from 3 au. The cylindrical radius or foot point of the loss is obtained by tracing backwards the corresponding streamline from the outer radial boundary down to 5 scale heights above the disk mid plane. The flow field is averaged over 50 orbits at 1 au starting at 600 orbits. The vertical dashed lines mark the limit of the region where 99% of the mass loss occurred.

Fig. 6 Numerically integrated dust trajectories for snapshots after 500 T0 of the simulation runs phc2, b5c2 and b4c2. The orange color map represents the gas density and the thin arrowed white lines denote the gas velocity streamlines. The thicker colored lines visualize the dust trajectories for grain sizes ranging from 0.1 μm to 10 μm.

3.4 Maximum grain size and flow angle

In order to quantify the maximum entrained grain size and the entrainment angle, we chose the approach to numerically integrate streamlines of various dust species based on the gas velocity field from the simulations. The following system of equation was solved numerically with the given velocity field of the corresponding simulation output: $$v_{\text{d},\text{r}}=v_{\text{d},\varphi }-\frac{G{M}_{*}}{{\left(r^2+z^2\right)}^{\text{n}rac32}}r-\frac{v_{\text{d},\text{r}}-v_{\text{g},\text{r}}}{t_{\text{s}}},v_{\text{d}],\text{z}}=-\frac{G{M}_{*}}{{\left(r^2+z^2\right)}^{\frac{3}{2}}}z-\frac{v_{\text{d},\text{z}}-v_{\text{g},\text{z}}}{t_{\text{s}}}.$$

The numerical integration was carried out with the VODE solver Brown et al. (1989). We verified that the streamlines were comparable to the ones in the dust velocity field output from the simulation. No significant deviations were visible.

The results are shown in Fig. 6 for 20 logarithmically distributed grain sizes ranging from 0.1 μm to 10 μm. The resulting trajectories match the corresponding dust streamlines from the actual simulation. As the starting point we chose r = 8 au, z = 2.5 au which is situated well above the wind launching front. As expected from the simulation outputs 10 μm grains fall back towards the disk surface in the case of winds including photoevaporation (models phc2, b5c2 and b4c2). Generally, grains are lifted more efficiently in the case of warm, magnetic winds. Grains larger than ≈ 3 μm to 4 μm are lifted further but eventually fall back onto the disk at larger radii in the presence of a thermal wind.

For warm, magnetically driven winds, the limiting grain size shifts towards ≈ 6 μm. In a cold, magnetic wind the picture changes and only submicron dust grains are lifted far from the disk surface and particles larger than 1 μm fall back onto the disk or are not entrained at all. The decreased wind speed is responsible for the lower dust entrainment efficiency as mentioned in Sect. 3.3.

Since the dust flow inclination angle is strongly dependent on the grain size, we quantified the inclination by fitting linear functions to the dust trajectories presented in Fig. 6. We only considered the final 25% part of the cylindrical radius of the trajectory for the fits. The corresponding slopes and inclination angles based on the simulations phc2, b5c2, b4c2 and b4c2np are plotted in Fig. 7. For small grains and warm winds the difference in the inclination angle is rather small. The inclination of the flow of 1 μm grains ranges from roughly 67° to 69 °. Considering 3 μm grains however, the inclination angle is 36.1° assuming a thermal wind compared magnetic winds with 44.4° and 46.5° for β = 10⁵ and β = 10⁴, respectively.Generally, the inclination angle is larger for warm, magnetic winds, especially considering micron sized grains.

A significant difference in inclination is however observed in a cold magnetic wind. While 0.1 μm-sized dust grains are entrained by an angle of 59.1°, the inclination drops down rapidly to roughly 25° for micron-sized grains, as expected.

Fig. 7 Asymptotic dust flow inclination depending on the grain size for the simulation runs phc2, b5c2, b4c2 and b4c2np.

Fig. 8 Normalized vertical density slices of the phc2 at different radial locations in the disk. The violet lines represent the normalized gas density while the turquoise lines display the dust density of 10 μm grains. The gray lines denote the respective Stokes number profile.

Fig. 9 Normalized vertical density slices and Stokes number profiles of the phc2 similar to Fig. 8 for the grain sizes 0.1 μm (straight lines), 1 μm (dashed lines), 10 μm (dotted lines) at 5 au.

3.5 Wind launching surface

Since the dust entrainment efficiency seems to be dependent on the dust reservoir at the location of the foot point of the wind, we therefore provide a more detailed analysis of the radial dependencies of the dust entrainment efficiency in this subsection.

We only focus on winds including photoevaporation where the wind launching surface is marked by a sharp transition in the gas density. Figures 8 and 9 visualize vertical slices of normalized gas and dust densities as well as the respective Stokes number profile. As expected, in Fig. 8 the dust scale height of 10 μm sized grains narrows with increasing radius as the Stokes number within the bulk of the disk increases due to the negative radial density slope. The kink at roughly 4.5 gas scale heights represents the wind launching surface. Clearly, the dust scale height is too small to lift a considerable amount of grains to the launching surface, especially at larger radii.

A comparison between the different dust species at 5 au is given in Fig. 9. The trend of larger dust scale heights with smaller Stokes numbers is readily visible. Weaker wind flows by several orders of magnitude appear for 1 μm sized grains compared to 0.1 μm particles. Only the vertical distribution of the smallest 0.1 μm grains that are well coupled to the gas appear to correspond the gas density distribution.

In the following, we aim to examine the dust entrainment efficiency depending on the radial location in the disk. We therefore traced various quantities, such as dust-to-gas ratio, Stokes number, temperature and gas number density along the wind launching surface. Numerically we defined this surface to be located about 0.3 gas scale heights above the kink in the vertical gas density profile in order to decrease the impact of fluctuations close to the surface of the disk. The kink location was determined by numerically evaluating local extrema in the second derivative of the vertical gas density profile which we find to be a robust method to find the launching surface.

The corresponding quantities at these positions are displayed in Fig. 10 where thermally and magnetically driven winds with 0.1 μm and 1 μm-sized grains are taken into account. Clearly, the dust entrainment efficiency, defined by the dust to gas ratio close to the wind launching front shown in the upper most panel of Fig. 10, decreases significantly with increasing radius. At the inner part of the flow close to the cavity at 2 au, 0.1 μm- sized grains are basically perfectly entrained with the wind for both photoevaporative and magnetically driven outflows since the initial dust-to-gas ratio is ϵ_dg = 10⁻². In the outer part of the disk, the dust-to-gas ratio in the wind decreased by one to two orders of magnitude. For 1 μm-sized grains the effect is more pronounced with a decline of six orders of magnitude. It becomes apparent that generally the magnetically driven wind with β = 10⁴, represented by the dotted lines, leads to a higher dust-to-gas ratio in the wind flow compared to weaker field strengths or a purely thermally driven wind.

One could argue that dust entrainment becomes less effective at larger radii because of the more diluted gas in the wind and thereby larger Stokes numbers. This is however not the case as shown in the second panel of Fig. 10. In the inner part of the disk up to ≈8 au to 9 au the Stokes number close to the launching surface increases until it stagnates and remains constant at the outer part of the disk. The decrease in dust entrainment efficiency can thus not be attributed to a weaker coupling between gas and dust in these regions. The gas density in the wind close to the launching surface indeed decreases with radius as shown in the bottom panel of Fig. 10. On the other hand, the gas temperature heated by the ionizing radiation from the central star in turn also decreases with radius. The peak in temperature is reached at ≈ 8 au to 9 au, agreeing with the stagnation point of the Stokes number profile. By linear fitting the quantities at the wind launching surface in logarithmic space we determine the power law slopes. Beyond 8 au the thermal velocity v_th decreases ∝ r^−0.13 whereas the gas density profile follows the power law ∝ r^−1.41. Looking at the definition of the Stokes number in Eq. (8), we expect an approximately constant relation St ∝ Ω_K∕(ρ_g ⋅ v_th) ∝ r^0.04 with a Keplerian angular frequency following Ω_K ∝ r^−1.5.

Fig. 10 Physical quantities at the wind launching surface of the simulation runs phc2 (straight lines), b5c2 (dashed lines) and b4c2 (dotted lines). In the top two panels the blue lines represent 0.1 μm -sized grains, whereas the green lines denote 1 μm-sized grains.

Fig. 11 Streamlines starting from r = 8 au, z = 2.5 au. Shown are poloidal velocities of gas and dust along the streamline (left two panels) and the respective Stokes numbers St as well as the Alfvén Mach number MA (right two panels).

3.6 Gas and dust streamlines

In Fig. 11 velocities, Stokes number and Alfvén Mach number are plotted for a representative streamline starting at r= 8 au. Since the flow is close to a steady state we can consider the streamline as a proxy for the actual trajectory of a dust and gas parcel moving in the wind. The radial velocities in the left-most panel are almost constant with the cylindrical radius whereas the vertical velocities in z-direction significantly increase up to values on the order of 15 km s⁻¹. The models including magnetic effects generally lead to slightly faster wind flows increasing with the magnetic field strength. Depending on the dust grain size the velocity can be significantly less compared to the gas velocity, differing by roughly 5 kms⁻¹. Dust velocities are also consistently faster in the magnetic wind case. Looking at the Stokes numbers of the two smallest dust species, the values up to unity are reached for micron-sized grains far from the disk. The right-most panel demonstrates that the wind is onlysubalfvénic for a small part close to the wind launching front. Even in the simulation run b4c2 with an initial β =10⁴ the wind quickly reaches superalfvénic speeds which in consequence represents a small magnetic lever arm. This picture is in line with recent models of magnetothermal winds such as Gressel et al. (2015); Bai (2017); Rodenkirch et al. (2020).

3.7 Synthetic observations

Given that the wind flow is rather thin and dust is not perfectly entrained for larger grain sizes and distances from the central star, the observability of such dusty winds deserves a closer look. As a first simple estimate we computed the radial optical depth depending on the latitude θ. The results are displayed in Fig. 12 for three dust species ranging from 0.1 μm to 10 μm in the H-band. The optical depth τ_obs,i of dust species i was calculated by $${\tau }_{\text{obs},\text{i}}={\displaystyle {\int }_{r_{\text{min}}}^{r_{\text{max}}}{\kappa }_{\text{sca},\text{i}}{\rho }_{\text{d},\text{i}}}$$(39)

with the scattering opacity κ_sca,i. The inner and outer radial boundaries are denoted by r_min and r_max, respectively.In Fig. 12, the bulk of the disk appears to be optically thick for all dust species in the H-band. The wind itself however is optically thin throughout the whole flow and its optical depth reaches values of τ_obs≈ 10⁻¹ in the direction of the wind base. Naturally, the optical depth decreases with larger viewing angles due to the thinner dust and gas density in these regions.

For 1 μm-sized grains the optical depth is comparable to the one for 0.1 μm grains but tapers off at a flatter angle of 70° due to the less inclined dust flow. Finally, the larger 10 μm grains are almost irrelevant since dust entrainment is negligible and the densities are therefore approaching the dust density floor value. In Figs. 13 and 14, synthetic observations are shown with inclination angles of 20 ° and 0° in latitude, respectively. In comparison to the disk brightness the wind structure can be clearly identified as an hourglass-shaped signature. Smaller 0.1 μm grains are confined to a narrower cone with diffuse emission towards the outer regions. A rather sharp transition between the region devoid of dust close to the rotation axis and the inner most part of the dusty wind is visible.

The wind is not completely illuminated by the stellar light because of the decreasing dust content towards the outer regions of the disk. Comparing the morphology of the scattered light signatures from 0.1 μm grains to emission of 1 μm grains, the larger dust particles lead to a conical shape with a wider opening angle with respect to the rotation axis. This feature is expected from Fig. 7 where the increased grain size leads to weaker gas-dust coupling and therefore causes a dust flow with a smaller inclination. In addition, larger grains are entrained less efficiently in the outer regions of the disk compared to smaller grains, due to the stronger dust settling towards the disk mid plane. In combination with the larger opening angle of the flow we therefore detect a less diffuse, conical shape with a sharp transition in brightness towards the outer regions of the disk.

It is also visible that magnetic winds with stronger fields confine the dust flow more towards the rotation axis and exhibit a larger inclination angle of the dust flow. The cold MHD wind model exhibits a brighter wind region compared to the models including photoevaporation. As expected, the opening angle of the 0.1 μm wind flow is larger compared to the warm wind models and resembles the emission of 1 μm if photoevaporative heating is included. In the cold wind model 1 μm-sized grains are entrained in a shallow inclination angle and could be identified as a disk with a larger aspect ratio than the thermodynamics would otherwise allow in the system. At a viewing angle of 20 ° the backside of the disk almost shadows the dusty wind, making it difficult to identify such a dusty outflow on the basis of a low resolution or a low signal-to-noise ratio. Generally, we find the brightness contrast between wind and disk to be less in the case of 1 μm grains if the viewing angle is 20°. In the edge-on view shown in Fig. 14, the wind signature brightness is comparable to the disk brightness for both grain sizes.

Fig. 12 H-band radial optical depth τobs originating from the central star depending on the inclination expressed in latitude for three dust species with sizes 0.1 μm, 1 μm and 10 μm indicated by the color scheme. The continuous lines refer to the photoevaporation simulation phc2, the dashed lines to the magnetic wind b5c2 with β = 105 and the dotted lines to the run b4c2 with β = 104.

4 Discussion

In the following, we discuss limitations of our wind model and compare the key results with other works in the literature. An advantage of the fluid description over a Lagrangian particle formulation of the dust dynamics in the simulations, is the ability to resolve low density regions if a large density contrast between the midplane and the corona as well as the wind is present, which is the case in the disk-wind system shown here. The disadvantage is the lack of proper modeling of a dust grain size distribution. A single-size dust population is unlikely due to dust growth and fragmentation. The synthetic observations thus only serve as a hint towards the detectability of dusty winds and the expected shape of the flow assuming the corresponding predominant dust size. Including multiple dust fluids in the synthetic observations would lead to misleading structures in the emission because of the rather sharp transition between the inner dust free cone and the wind region.

Deducing the maximum dust grain size from the wind flow angle seems to be rather difficult since the underlying parameter space is degenerate. On the one hand, the opening angle of the dusty wind increases with the grain size. On the other hand, the same effect occurs with colder magnetically driven winds. A distinction could be made by probing the wind speed, which is smaller in the cold wind case, albeit being a specific example.

If the wind is optically thin for the ionizing radiation of the central star, the difference between photoevaporative winds and warm magnetically driven winds is minor. In future works, the parameter space could be extended by varying the luminosity of the central star and by using a more sophisticated photoevaporation model as presented in Ercolano et al. (2021); Picogna et al. (2021).

Within the simulated time frame we do not observe any significant asymmetry between the upper and lower wind. Such asymmetric flows were observed in simulations by Béthune et al. (2017), Suriano et al. (2018) and Riols & Lesur (2018). Ring formation by magnetically driven winds does not occur in our simulation runs compared to Suriano et al. (2018) and Riols & Lesur (2018) where ambipolar diffusion is prescribed to be the dominant nonideal MHD effect. In their work however, larger field strength in the regime of β = 10³ were probed. Riols & Lesur (2018) also observed ring formation for β = 10⁴ and included fluid dust focusing on the dust dynamics in the mid plane.

Given that the same photoevaporation recipe is used in this paper as the one in Franz et al. (2020) our results are compatible in this regime. They conclude a maximum entrainable grain size of 11 μm. As shown in Fig. 6 our photoevaporative models allow a short-term lifting of 10 μm dust grains which eventually fall back onto the disk. In the complete dust fluid picture, the small dust grain scale height prevents any dust entrainment starting from the wind launching front. More recently, Franz et al. (2022) studied the observability of dust entrainment in photoevaporative winds based on their models in Franz et al. (2020). They conclude that the wind signature in scattered light is too faint to be observed with SPHERE IRDIS, but might be detectable with the JWST NIRCam under optical conditions. In their synthetic observations a similar “chimney”-shaped wind emission is visible compared to our photoevaporative models.

Giacalone et al. (2019) found a positive correlation between gas temperature and maximum dust grain size in their semi-analytical magneto-centrifugal wind model which corresponds best to our cold magnetically driven wind model b4c2np. Our results corroborate their findings of a maximum entrainable dust size in the submicron regime in a cold wind around a typical T Tauri star. Our models including photoevaporation furthermore agree with Hutchison et al. (2016a) where the maximum grain size was determined to be < 10 μm depending on the disk radius for a typical T Tauri star. They additionally argue that dust settling most certainly limits this maximum size further, which we confirm in our models, where the diminishing dust entrainment efficiency can be attributed to the dust settling towards the mid plane.

In recent scattered light observations of RY Tau (Garufi et al. 2019), a broad dusty outflow obstructing the underlying disk was detected. Both lobes are inclined by roughly 45 °. The presence of a jet in this system advocates a magnetic wind launching mechanism that causes the features observed in the near-infrared. We do not observe an optically thick wind that would be able to completely obstruct the disk, since the dust densities are too low, especially in the models including photoevaporation. A possibility might hence be a slow, cold, magnetic wind at higher field strengths. The large opening angle in the observation would suggest a wind flow containing dust grains of the size of several μm if photoevaporation would be the significant wind mechanism. If a cold magnetically driven wind was the dominant factor, smaller submicron grains would be the solution compatible with our models.

In this work, we neglected the effect of radiation pressure that may be able to blow out small dust grains. Franz et al. (2020) argued that a photoevaporative wind can entrain grains of sizes roughly 20 times larger compared to a compatible radiation pressure model from Owen & Kollmeier (2019). Since we use similar parameters for the photoevaporation model, the effect of radiation pressure should not significantly affect the dust dynamics in the studied parameter space.

In the photoevaporation models, the gaseous part of the wind reaches temperatures of up to 10⁴ K. We verified thatthe dust temperature stays below the dust sublimation temperature in Appendix A.

Fig. 13 Radiative transfer images in the H-band (λ = 1.6 μm) of simulation runs with photoevaporation (left most column), β = 105 (center left column), β = 104 (center right column) and a cold magnetic wind with β = 104 (right most column). The upper panels are based on dust grains with a = 0.1 μm whereas the lower panels demonstrate synthetic images with a = 1 μm. The object is modeled at a distance of 150 pc with a solar type star at an inclination of 20°.

Fig. 14 Radiative transfer images in the H-band (λ = 1.6 μm) of the same simulation runs as in Fig. 13 but with an edge-on view at an inclination of 0°.

5 Conclusion

We presented a suite of fully dynamic, global, multifluid simulations that model dust entrainment in photoevaporative as well as magnetically driven disk wind using the FARGO3D code. We demonstrated that both types of winds are able to efficiently transport dust in the wind flow and addressed the observability of such dusty winds with the help of the RADMC-3D radiative transfer code. In the following we summarize the main results:

• The maximum entrainable dust grain size a_max both depends on the type of the wind launching mechanism and the turbulent diffusivity of the disk. Whereas magnetic winds including XEUV-photoevaporation only show minor differences in a_max, ranging from 3 μm to 6 μm, cold magnetically driven winds at comparable field strengths only entrain submicron sized grains efficiently.

• With increasing radius the dust to gas ratio in the wind drops rapidly, mainly due to the smaller dust scale height in comparison to the gas scale height. Dust grains are unable to reach the wind launching surface in the outer regions of the disk and the dust content in the wind starting from these locations hence decreases. In the case of warm, ionized winds the Stokes number stays approximately constant with radius beyond 8°. Generally, magnetically driven winds including photoevaporation are slightly more efficient in entraining dust compared to pure photoevaporative flows.

• The dust flow angle θ_d decreases with increasing grain size to the weaker coupling between gas and dust. In warm, ionized winds the flow angle reaches values of 67° ≤ θ_d ≤ 69° for a = 1 μm sized grains. At a grain size of a = 3 μm a photoevaporative dust flow provides θ_d ≈ 36° whereas a warm magnetic wind leads to a steeper dust flow with θ_d ≈ 45°. In a cold magnetic wind the flow angle is significantly smaller, being θ_d ≈ 59° for a = 0.1 μm and θ_d ≈ 25° for micron-sized grains.

• Radiative transfer images in the H-band suggest a rather weak, optically thin emission for photoevaporative winds. A more pronounced, conical shape is visible considering 1 μm grains. The larger dust flow angle θ_d leads to a narrower emission region in the case of warm, ionized magnetic winds. Cold magnetic winds appear with a significantly larger opening angle with comparable grain sizes. The emission region of the wind appears brighter and more easily detectable due to the increased dust density of the slower magnetic wind.

Our results presented in this paper might help to constrain the magnetic field strength as well as the dominating wind mechanism if such dust signatures will be detected in future observations. The dust grain size could be determined via the inclination angle of the flow. Although the parameter space is degenerate, probing the wind speeds might help to differentiate between thermally and cold magnetically driven winds.

Acknowledgements

We would like to thank R. Franz for helpful discussions. We furthermore thank the anonymous referee for helpful comments that further improved the manuscript. Authors Rodenkirch and Dullemond acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 325594231 under grant DU 414/23-1. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 37/935-1 FUGG. Plots in this paper were made with the Python library matplotlib (Hunter 2007).

Appendix A Dust temperature estimate

The relatively high temperatures in the ionized wind models heat the dust grains that are entrained in the gas flow. We show in the following that dust particles do not reach the dust sublimation threshold and would thus remain solid along their trajectories. The expression of the heating rate per unit surface q_d is described in Gombosi et al. (1986), taking into account both contact and friction heating. We follow the approach of Stammler & Dullemond (2014) in this regard and evaluate the subsequent expressions: $$q_{\text{d}}={\rho }_{\text{g}}{C}_{\text{H},\text{d}}\left(T_{\text{rec}}-T_{\text{d}}\right)|v_{\text{g}}-v_{\text{d}}|,$$(A.1)

where the heat transfer coefficient C_H,d is expressed as $$C_{\text{H},\text{d}}=\frac{\gamma +1}{\gamma -1}\frac{k_{\text{B}}}{8\mu m_{\text{p}}{s}^2}\left[\frac{s}{\sqrt{\pi }}\text{exp}\left(-s^2\right)+\left(\frac{1}{2}+s^2\right)\text{erf}\left(s\right)\right],$$(A.2)

with the error function erf and the normalized difference between dust and gas velocities: $$s=\frac{|v_{\text{g}}-v_{\text{d}}|}{\sqrt{2k_{\text{B}}{T}_{\text{g}}/(\mu m_{\text{p}})}}.$$(A.3)

The recovery temperature T_rec is given by the expression (Gombosi et al. 1986): $$T_{\text{rec}}=\frac{T_{\text{g}}}{\gamma -1}\left[2\gamma +2(\gamma -1)s^2-\frac{\gamma -1}{\frac{1}{2}+s^2+\frac{s}{\sqrt{\pi }}\text{exp}\left(-s^2\right){\text{erf}}^{-1}\left(s\right)}\right].$$(A.4)

The dust equilibrium temperature is then finally computed by solving $$q_{\text{d}}+\frac{R_{*}}{r^2}{ϵ}_{\text{d}}{\sigma }_{\text{SB}}{T}_{*}^4-4{\sigma }_{\text{SB}}{T}_{\text{d}}^4=0$$(A.5)

with the thermal cooling efficiency factor $${ϵ}_{\text{d}}=\frac{{\kappa }_{\text{P}}(T_{*})}{{\kappa }_{\text{P}}(T_{\text{d}})},$$(A.6)

using the Planck mean opacity κ_P with the corresponding temperature. Fig. A.1 provides the solutions of the equilibrium dust temperatures in comparison with the gas temperatures in a photoevaporative flow including 0.1 μm grains. We assume a sun-like star in the computation and find maximum dust temperatures of T_d ≈ 345 K at a radial distance of r = 2 au. Due to the thin gas in the wind region the ambient heating has an almost negligible impact on the dust temperature. Without the flux of the central star the dust grains would reach T_d ≈ 40 K only being heated by the ionized gas. In the wind flows presented here dust temperatures are clearly below the dust sublimation threshold which would be ≈ 980 K assuming pyroxene as the grain material embedded in a thin medium of the density ρ_g = 1 × 10⁻¹⁶ g cm⁻³ (Pollack et al. 1994).

Fig. A.1 Vertical slices of gas and dust temperatures at different radial location in the disk and wind region of the photoevaporation simulation phc2. The gray lines represent the gas temperatures whereas the red lines display the equilibrium dust temperatures of the 0.1 μm grains.

Appendix B Ionization fraction

Since the ionization fraction depends on the dust-to-gas ratio ϵ_ion in the disk a comparison between ϵ_ion = 10⁻³ and ϵ_ion = 10⁻² of the simulation runs b4c2 and b4c2eps2 is shown in Fig. B.1. The ionization is relatively insensitive to the increased dust-to-gas ratio and only decreases by a small amount in the weakly ionized disk mid plane. In the dust density slice the effects are negligible.

Fig. B.1 Normalized dust density (blue lines) and ionization fraction (gray lines) in the simulations b4c2 and b4c2eps2 with a dust-to-gas ratio of ϵ = 10−3 and ϵ =10−2, respectively. The snapshot is taken after 1000 orbits at 1 au.

Appendix C Super time stepping

This sectiondescribes the implementation of the Super Time-Stepping technique (STS) for parabolic or mixed advection-diffusion problems. Introduced by Alexiades et al. (1996), the first-order method divides one explicit super step ΔT in multiple substeps τ₁, τ₂, …τ_N. The algorithm ensures stability while maximizing $\Delta T={\sum }_{\text{i}}^{\text{N}}{\tau }_{\text{i}}$ The optimized length of each substep is based on Chevychev polynomials: $${\tau }_{\text{j}}=\Delta t_{\text{exp}}\left[(\nu -1)\text{cos}\left(\frac{2j-1}{N}\frac{\pi }{2}\right)+\nu +1\right],$$(C.1)

where Δt_exp is the explicit time step. The sum of all substeps then gives: $$\Delta T=\Delta t_{\text{exp}}\frac{N}{2\sqrt{\nu }}\left[\frac{{(1+\sqrt{\nu })}^{2N}-{(1-\sqrt{\nu })}^{2N}}{{(1+\sqrt{\nu })}^{2N}+{(1-\sqrt{\nu })}^{2N}}\right].$$(C.2)

The parameter ν adopts values between 0 and 1. In the limit of ν → 0 the method is N times faster than the explicit integration. Lower values of ν however can lead to oscillations and eventually to an unstable solution. In FARGO3D the magnetic field is updated with the method of characteristics and constrained transport (MOCCT) (Hawley & Stone 1995) whereby the constrained transport method described in Evans & Hawley (1988) is used (Benitez-Llambay & Masset 2016). In this method the evolution of the magnetic field is determined by computing the electromotive forces (EMFs) Benitez-Llambay & Masset (2016): $$E=v\times B.$$(C.3)

We compute the contribution to the total EMFs of each substep j with duration τ_j. The result is added to the EMFs, weighted by a factor of τ_j∕ΔT. The magnetic field updates at every substep are only applied on a temporary buffer field and the main B-field is then updated with the total EMFs, containing contributions from nonideal terms and the ordinary MHD dynamics.

In order to test the correctness of the implementation, a decaying standing Alfvén wave test for ambipolar diffusion as described in Choi et al. (2009) Sec. 4.2 is used. The time dependence of the first normal mode is given by Choi et al. (2009): $$h(t)=h_0\text{sin}({\omega }_{R}t)e^{{\omega }_{I}t}$$(C.4)

with ω = ω_R + iω_I. The terms ω_R and ω_I are the real and imaginary parts of the angular frequency of the Alfvén wave. The dispersion relation is given by Choi et al. (2009): $${\omega }^2+i\frac{c_{\text{A}}^2{k}^2}{{\langle \sigma v\rangle }_{\text{i}}{n}_{\text{i}}}\omega +c_{\text{A}}^2{k}^2=0,$$(C.5)

where k is the real wave number and c_A the characteristic Alfén speed.

In the test according to Choi et al. (2009) the magnetic field is set to $B=B_0\widehat{\text{x}}$ with B₀ = 1. The density is uniformly set to 1. The wave is initialized in the $\widehat{x}$ direction with an initial velocity of: $$v_{\text{y}}=v_0{c}_{\text{A}}\text{sin}(kx).$$(C.6)

The amplitude of the perturbation is set to v₀ = 0.1, the characteristic Alfén velocity to $1/\sqrt{2}$, the ion density to ρ_i = 0.1 and the wave number to k = 2π∕L where L = 1 is the domain size. The collision rate coefficient ${\langle \sigma v\rangle }_{\text{i}}$ is varied from 100 to 1000. Instead of an oblique wave test, the simulations are run in a one-dimensional configuration.

In Fig. C.1 and Fig. C.2 the test results are shown and the numerical results fit the analytical solution well. Stability is found to be fulfilled up to five substeps which is sufficient for the simulation runs presented in this paper. An overall loss of accuracy in comparison to the solution without STS is also anticipated since this approach is only first-order accurate.

Fig. C.1 Decaying Alfvén wave test for ${\langle \sigma v\rangle }_{\text{i}}=1000$. The dots correspond to the simulation results whereas the dashed line is the analytical solution.

Fig. C.2 Decaying Alfvén wave test for ${\langle \sigma v\rangle }_{\text{i}}=100$ similar to Fig. C.1

References

All Tables

All Figures

	Fig. 1 Gas density maps of the simulations phc2, b5c2 and b4c2 after 500 orbits at 1 au. Magnetic field lines are represented by the white lines in the corresponding panels. The color map scaled logarithmically annotated by the exponent to a base of ten and represents the volume density of the gas in g cm−3.
In the text

	Fig. 2 Total wind mass loss rates of gas and dust for the simulation runs phc2, b5c2, b4c2 and b4c2np. Results from 10 μm grains are excluded due to the negligible outflow and the lack of entrainment in the wind for all runs shown above.
In the text

	Fig. 3 Dust density maps for the simulations phc2, b5c2 and b4c2 after 500 orbits at 1 au. Dust velocity streamlines are annotated by the white arrowed lines. The upper, center and lower three panels visualize the flow of grains with size 0.1 μm, 1 μm and 10 μm, respectively.
In the text

	Fig. 4 Dust density maps for the simulation b4c2np after 500 orbits at 1 au. Dust velocity streamlines are annotated by the white arrowed lines.
In the text

Fig. 5 Cumulative mass loss rates of gas and dust for the simulation runs phc2, b5c2, b4c2 and b4c2a3 starting from 3 au. The cylindrical radius or foot point of the loss is obtained by tracing backwards the corresponding streamline from the outer radial boundary down to 5 scale heights above the disk mid plane. The flow field is averaged over 50 orbits at 1 au starting at 600 orbits. The vertical dashed lines mark the limit of the region where 99% of the mass loss occurred.

In the text

	Fig. 6 Numerically integrated dust trajectories for snapshots after 500 T0 of the simulation runs phc2, b5c2 and b4c2. The orange color map represents the gas density and the thin arrowed white lines denote the gas velocity streamlines. The thicker colored lines visualize the dust trajectories for grain sizes ranging from 0.1 μm to 10 μm.
In the text

	Fig. 7 Asymptotic dust flow inclination depending on the grain size for the simulation runs phc2, b5c2, b4c2 and b4c2np.
In the text

	Fig. 8 Normalized vertical density slices of the phc2 at different radial locations in the disk. The violet lines represent the normalized gas density while the turquoise lines display the dust density of 10 μm grains. The gray lines denote the respective Stokes number profile.
In the text

	Fig. 9 Normalized vertical density slices and Stokes number profiles of the phc2 similar to Fig. 8 for the grain sizes 0.1 μm (straight lines), 1 μm (dashed lines), 10 μm (dotted lines) at 5 au.
In the text

	Fig. 10 Physical quantities at the wind launching surface of the simulation runs phc2 (straight lines), b5c2 (dashed lines) and b4c2 (dotted lines). In the top two panels the blue lines represent 0.1 μm -sized grains, whereas the green lines denote 1 μm-sized grains.
In the text

	Fig. 11 Streamlines starting from r = 8 au, z = 2.5 au. Shown are poloidal velocities of gas and dust along the streamline (left two panels) and the respective Stokes numbers St as well as the Alfvén Mach number MA (right two panels).
In the text

	Fig. 12 H-band radial optical depth τobs originating from the central star depending on the inclination expressed in latitude for three dust species with sizes 0.1 μm, 1 μm and 10 μm indicated by the color scheme. The continuous lines refer to the photoevaporation simulation phc2, the dashed lines to the magnetic wind b5c2 with β = 105 and the dotted lines to the run b4c2 with β = 104.
In the text

Fig. 13 Radiative transfer images in the H-band (λ = 1.6 μm) of simulation runs with photoevaporation (left most column), β = 105 (center left column), β = 104 (center right column) and a cold magnetic wind with β = 104 (right most column). The upper panels are based on dust grains with a = 0.1 μm whereas the lower panels demonstrate synthetic images with a = 1 μm. The object is modeled at a distance of 150 pc with a solar type star at an inclination of 20°.

In the text

	Fig. 14 Radiative transfer images in the H-band (λ = 1.6 μm) of the same simulation runs as in Fig. 13 but with an edge-on view at an inclination of 0°.
In the text

	Fig. A.1 Vertical slices of gas and dust temperatures at different radial location in the disk and wind region of the photoevaporation simulation phc2. The gray lines represent the gas temperatures whereas the red lines display the equilibrium dust temperatures of the 0.1 μm grains.
In the text

	Fig. B.1 Normalized dust density (blue lines) and ionization fraction (gray lines) in the simulations b4c2 and b4c2eps2 with a dust-to-gas ratio of ϵ = 10−3 and ϵ =10−2, respectively. The snapshot is taken after 1000 orbits at 1 au.
In the text

	Fig. C.1 Decaying Alfvén wave test for ${\langle \sigma v\rangle }_{\text{i}}=1000$. The dots correspond to the simulation results whereas the dashed line is the analytical solution.
In the text

	Fig. C.2 Decaying Alfvén wave test for ${\langle \sigma v\rangle }_{\text{i}}=100$ similar to Fig. C.1
In the text