Non-ideal MHD simulation of HL Tau disk: formation of rings

TL;DR

This study uses 3D non-ideal MHD simulations to explore how snow lines and magnetic effects can produce rings and gaps in the HL Tau disk, aligning with recent high-resolution observations.

Contribution

It introduces a comprehensive model combining dust evolution, non-ideal MHD effects, and snow line physics to explain ring formation in protoplanetary disks.

Findings

Gaps and rings form rapidly from varying accretion rates across snow lines.

Ambipolar diffusion causes preferential midplane accretion and magnetic reconnection.

Magnetic effects and snow lines together can produce observable disk structures.

Abstract

Recent high resolution observations unveil ring structures in circumstellar disks. The origin of these rings has been widely investigated under various theoretical scenarios. In this work we perform global 3D non-ideal MHD simulations including effects from both Ohmic resistivity and ambipolar diffusion (AD) to model the HL Tau disk. The non-ideal MHD diffusion profiles are calculated based on the global dust evolution calculation including sintering effects. Disk ionization structure changes dramatically across the snow line due to the change of dust size distribution close to snow line of major volatiles. We find that accretion is mainly driven by disk wind. Gaps and rings can be quickly produced from different accretion rates across snow line. Furthermore, ambipolar diffusion (AD) leads to highly preferential accretion at midplane, followed by magnetic reconnection. This results a local zone of decretion that drains of mass in the field reconnection area, which leaves a gap and an adjacent ring just outside it. Overall, under the favorable condition, both snow lines and non-ideal MHD effects can lead to gaseous gaps and rings in protoplanetary disks.

Equations

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{v}\times \mathbf{B}\right)-\frac{4\pi }{c}\nabla \times \left({\eta }_{\mathrm{O}}\mathbf{J}+{\eta }_{\mathrm{H}}\mathbf{J}\times \mathbf{b}+{\eta }_{\mathrm{A}}{\mathbf{J}}_{\perp }\right),$$

$${\eta }_O\left(R,z\right)$$

$$Am(R,z)$$

$$T(R,z)=\{\begin{array}{ll}T(R,0)& \text{if }|z|<2h\\ T(R,0)+49.5T(R,0)\times & \\ (1+\sin (\frac{|z|/h-2}{10}\pi -\frac{\pi }{2}))& \text{if }2h\le |z|\le 12h\\ 100T(R,0);& \text{if }|z|>12h\end{array}$$

$$\rho (R,z=0)$$

$$T(R,z=0)$$

$${\rho }_{fl}(R,z)=\{\begin{array}{ll}{\rho }_{0,fl}(\frac{r_{min}}{R_0}{)}^p(\frac{r_{min}}{z}{)}^2\frac{r_{min}}{\sqrt{R^2+z^2}}& \text{if }R\le r_{min}\\ {\rho }_{0,fl}(\frac{R}{R_0}{)}^p(\frac{r_{min}}{z}{)}^2\frac{r_{min}}{\sqrt{R^2+z^2}}& \text{if }R>r_{min}\end{array}$$

$$B(R,z)=Bz(R)=\{\begin{array}{ll}{B}_0(\frac{r_{min}}{R_0}{)}^{\frac{p+q}{2}}& \text{if }R\le r_{min}\\ B_0(\frac{R}{R_0}{)}^{\frac{p+q}{2}}& \text{if }R>r_{min}\end{array}$$

$${\alpha }_{int}={\int }_{z-}^{z+}{\mathbf{T}}_{R\varphi }dz/{\int }_{z-}^{z+}{P}_{\mathrm{gas}}dz$$

$${\mathbf{T}}_{z\varphi }=\left(\rho v_z\left(v_{\varphi }-\overline{v_{\varphi }}\right)-B_z{B}_{\varphi }\right)$$

$$\frac{\dot{m}_{acc}v_{K}}{4\pi}=\frac{d}{dr}\left(R^{2}\int_{\rm bot}^{\rm up}{\bm{T}}_{R\phi}dz\right)+R^{2}\Big{|}_{\rm bot}^{\rm up}{\bm{T}}_{z\phi}$$

$$A_{\varphi }(r,\theta )=\frac{2B_{z0}{R}_0}{3-\alpha }{\left(\frac{R}{R_0}\right)}^{-\frac{\alpha -1}{2}}[1+(m\tan \theta {)}^{-2}{]}^{-\frac{5}{8}}$$

Full text

accretion, accretion disks; protoplanetary disk; magnetohydrodynamics (MHD)

1 Introduction

Accretion disks play an important role in forming many astrophysical objects. Protoplanetary disks, while supplying mass to central stars, are also birth places of planets. HL Tau is a Class I-II protostellar source surrounded by a circumstellar disk with a size of $\sim$100 au. The 2014 Long Baseline Campaign of the Atacama Large Millimeter/submillimeter Array (ALMA) provided unprecedented details of the HL Tau disk(ALMA Partnership et al., 2015). Patterns of multiple bright and dark rings that are symmetric to the central star are clearly seen. Subsequent ALMA observations find rings in dust continuum emission are very common in protoplanetary disks (e.g., Huang et al., 2018; Long et al., 2018).

Axisymmetric dust rings in a protoplanetary disk can be produced in several ways. With dust radial drift under the action of gas drag, a major mechanism of creating dust radial concentration is dust trapping at local gas pressure maxima (Whipple, 1972). Planet-disk interaction has been the most popular way to carve gaps in protoplanetary disks. The filtration effect on larger dust particles can make dust gaps more prominent than gas structure (e.g., Zhu et al., 2012; Zhang et al., 2018). The gap feature in HL Tau can be either explained by multi planets with one of them in each gap(e.g., Dong et al., 2015; Dipierro et al., 2015), or a single planet opening multiple gaps(Bae et al., 2018).

Another possibility is the condensation fronts of varies solids altering dust growth. The large dispersion of location and size of rings in from recent observations show low correlation to snow lines or resonance orbits of planets (Huang et al., 2018; Long et al., 2018; van der Marel et al., 2019). Though snow lines are not likely to be the common origin of rings and gaps, it still does a good job in the case of HL Tau. It was first proposed by Zhang et al. (2015) that smaller dust can be depleted by rapid pebble growth at snow lines. The location of dark rings in the ALMA image are close to snow lines of some major volatiles, e.g.,$\mathrm{H_{2}O}$ and $\mathrm{NH_{3}}$. A more detailed study on dust coagulation by Okuzumi et al. (2016) shows sintering effect between dust aggregates can change dust sizes near snow lines even with low volatile abundance.

Other theoretical models include secular gravitation instability triggered from low turbulence inner disk making ring-like structure (Takahashi & Inutsuka, 2016). Dust to gas back reaction and gas compression can cause large gradient of dust radial drift velocity, this self induced dust pile-ups is discussed in Gonzalez et al. (2017). If planets are formed early, dust rings can be the collision outcome of planetesimals perturbed by planets and orbital resonances (Boley, 2017). Magnetic self-organization is seen in form of ‘zonal flows’, resulting radial variations in density and magnetic flux. These banded rings can sustain as the pressure gradient is balanced by the Coriolis force (Johansen et al., 2009; Bai & Stone, 2014). Recent MHD simulations suggest formation of rings and gaps can be achieved through redistribution of magnetic flux, direct feeding of avalanche accretion and midplane magnetic reconnection(Suriano et al., 2017, 2018, 2019), when non-ideal MHD diffusivities like Ohmic resistivity and ambipolar diffusion are considered.

The abundance of dust is not just a disk profile tracer in radio continuum observation, but can also play a crucial role constraining disk accretion. In the traditional scenario of viscous driven accretion, gas surface density varies at radial viscosity jumps such as transition from magnetorotational instability (MRI) active zone to dead zone (e.g., Gammie, 1996; Flock et al., 2015; Hu et al., 2016). Disk ionization degree strongly depends on the abundance of small (sub-micron) grains, thus affecting the coupling between gas and magnetic fields (e.g., Wardle, 2007; Bai, 2011). This inspired the idea, that will be presented in this paper, to incorporate snow line induced dust distribution into MHD disk simulations. In the sintering operating regime, easier fragmentation produces more smaller grains (Okuzumi et al., 2016), reduces ionization rate and increases magnetic diffusion. This could generate radial variantion of acccretion rate, leads to axisymmetric structure in gas surface density.

In this study we first generated non-ideal MHD diffusivity profile from dust coagulation model including sintering effect. Then we implement the diffusion profile in global 3D non ideal MHD simulation to analyze disk accretion and density structure. In Section 2 we briefly introduce basic physics and numerical setups. In Section 3 we present diagnostics of disk radial and vertical structure. In Section 4 we explore the impact of variation field strength, shape and magnetic diffusion. We will summarize the paper in section 5.

2 Numerical setup

2.1 Non ideal MHD diffusion

We solve the magnetohydrodynamic (MHD) equations using Athena++ (Stone et al. 2019, in preparation) with non-ideal magnetic diffusion terms (Bai & Stone, 2017; Bai, 2017). The basic non-ideal induction equations are:

[TABLE]

where $\bm{v}$ is gas velocity, $\bm{B}$ is magnetic field, ${\bm{b}}={\bm{B}}/|B|$ is the unit vector representing field line direction. $\bm{J}$ is current density vector, with $\bm{J}_{\bot}$ as the current component perpendicular to the magnetic field. $\eta_{\mathrm{O}}$, $\eta_{\mathrm{H}}$ and $\eta_{\mathrm{A}}$ are Ohmic, Hall and ambipolar diffusivities, and here for simplicity, we only consider Ohmic resistivity and ambibolar diffusion (AD). The Hall diffusivity significantly complicates the picture, but they are more applicable to the inner disk (Béthune et al., 2016; Bai, 2017). Modeling HL Tau does not to resolve the inner disk,so our simulation mainly focused on outer part disk where most snow lines are.

We compute the radial and vertical profiles of the magnetic diffusivities in three steps. The first step is to calculate the radial ionization structure at the midplane by considering the balance among ionization, recombination in the gas, and grain charging. The grain size distribution is given by a power law with minimum and maximum grain sizes $a_{\rm min}$ and $a_{\rm max}$ (Equation 20 of Okuzumi et al. 2016). The minimum grain size $a_{\rm min}$ is fixed to $0.1~{}\micron$, whereas the maximum size $a_{\rm max}$ and total dust surface density $\Sigma_{d}$ are taken from a snapshot of a dust evolution calculation that takes into account the low stickiness of CO2 ice (Okuzumi & Tazaki, 2019) and aggregate sintering (Okuzumi et al., 2016). The fragmentation threshold velocity of nonsintered aggregates is assumed to be $0.73~{}\rm m~{}s^{-1}$ inside the H2O snow line, $7.3~{}\rm m~{}s^{-1}$ between the H2O and CO2 snow lines, and $0.73~{}\rm m~{}s^{-1}$ outside the CO2 snow line. In regions where aggregate sintering takes place, the fragmentation threshold is taken to be $40\%$ of the threshold for nonsintered aggregates. The dust evolution calculation assumes weak gas turbulence with a velocity dispersion of 1% of the sound speed. The snapshot is taken at the time at which the computed total millimeter fluxes from the dust disk match the ALMA observations (ALMA Partnership et al., 2015). The adopted radial profiles of $a_{\rm max}$ and $\Sigma_{d}$ are plotted in Figures 1(a) and (b). Because of sintering, dust particles behind the snow lines CO2, C2H6, and CH4 experience enhanced collisional fragmentation and pile up there. The loss of H2O ice interior to the H2O snow line also enhances collisional fragmentation and results in another traffic jam of dust in this inner region.

To compute the ionization balance for grains of various size distribution, we use the semi-analytical ionization model of Okuzumi (2009). Our ionizing sources include the short-lived radionuclide ${}^{26}{\rm Al}$ (Umebayashi & Nakano, 2009) and stellar X-rays (Igea & Glassgold, 1999; Bai & Goodman, 2009) with an X-ray luminosity of $10^{31}~{}\rm erg~{}s^{-1}$. The radial profile of the ionization rate per hydrogen molecule at the midplane, $\zeta$, is plotted in Figure 1(c). At $R\gtrsim 15~{}\rm au$, X-rays are able to penetrate down to the midplane and produce a high value of $\zeta$. All positive ions are represented by ${\rm HCO}^{+}$, with its recombination rate coefficient given by $2.4\times 10^{-7}(T/300~{}\rm K)^{-0.69}~{}\rm cm^{3}~{}s^{-1}$ (Ganguli et al., 1988). Only for dust evolution and ionization calculations, we assume the gas disk model $\Sigma_{g}=M_{\rm disk}\exp(-R/R_{c})/(2\pi R_{c}R)$ and $T=310(R/1~{}\rm au)^{-0.57}$ with $M_{\rm disk}=0.2M_{\sun}$ and $R_{c}=150~{}\rm au$ (Okuzumi et al., 2016), where $R$ is the radial distance. The electron and ion abundances obtained from our ionization calculation are shown in Figure 1(d). Reduction of the electron and ion abundances by dust grains is visible in the regions where small dust particles pile up (for comparison, $n_{e}/n$ ($=n_{i}/n$) in the dust-free case is also shown in Figure 1(d)). At $R\gtrsim 15~{}\rm au$, the effect of dust grains is relatively minor because of the elevated ionization rate.111It can be shown analytically that $n_{e},n_{i}\propto\zeta$ in the limit where grain charging dominates over gas-phase recombination, and $n_{e},n_{i}\propto\zeta^{1/2}$ in the opposite, dust-free limit (e.g., Okuzumi, 2009). These indicate that the relative importance of grain charging decreases as $\zeta$ increases.

The second step is to convert the midplane electron and ion abundances into the midplane Ohmic and ambipolar diffusivities, $\eta_{O}(R,0)$ and $\eta_{A}(R,0)$. We use the exact expressions of the diffusivities derived from the conductivity tensor (e.g., Equations (13) and (15) of Salmeron & Wardle, 2008) to compute $\eta_{O}(R,0)$ and $\eta_{A}(R,0)$ for $B=1~{}\rm mG$. While $\eta_{O}$ is independent of $B$, $\eta_{A}$ depends on $B$. However, as long as $n_{i}=n_{e}$, the ambipolar diffusivity is simply proportional to $B^{2}$ (Balbus & Terquem, 2001). Because the condition $n_{i}\approx n_{e}$ holds in most part of the disk at $R\gtrsim 4~{}\rm au$ (see Figure 1(d)), we assume $\eta_{A}\propto B^{2}$ when we compute $\eta_{A}(R,0)$ for $B\not=1~{}\rm mG$. Because of the $B^{2}$ dependence,

the dimensionless ambipolar Elsasser number, $Am=v_{a}^{2}/\left(\eta_{\mathrm{A}}\Omega_{K}\right)$ is better suited to describe the strength of ambipolar diffusion. Here $v_{a}=B/\sqrt{4\pi\rho}$ is the Alfven velocity, $\Omega$ is orbital angular velocity.

The obtained diffusivities are shown in the left panel of Figure 2 by the dashed lines. The diffusion gets exceptionally strong inside water snow line. To ensure a reasonable timestep, we focus on the $\mathrm{CO_{2}}$ and $\mathrm{C_{2}H_{6}}$ snow lines, and removes any feature within water snow line. The solid lines are actual profiles we used in our MHD simulations.

The final step is to extend the midplane diffusivities to the diffusivities at arbitrary $z$. The vertical variation of gas ionization state is more complicated when considering both thermal dynamics and chemical evolution (Wang & Goodman, 2017). For simplicity we adopt an analytical profile similar as Bai & Stone (2017). The strength of both diffusion effects drops quickly when getting away from disk portion, the difference is $Am$ remains constant within $2h$, and $\eta_{O}$’s decrease starts at midplane.

[TABLE]

Here $h$ is disk scale height at midplane. These profiles ensure strongest diffusion in the disk midplane and drops gradually over the atmosphere. The AD Elsasser number $Am$ is almost constant till $\approx 2h$ above disk midplane, making AD is more dominant over Ohmic resistivity at disk surface, where density is lower.

2.2 disk setup

To minimize grid noise of the radially flowing disk wind, we perform our simulation in spherical polar coordinates $(r,\theta,\phi)$. In this paper, we transform some quantities to cylindrical coordinates $(R,\phi,z)$ for analysis. Note that $R=1$ equals 10 $au$ in physical unit. For an accretion disk with the vertical hydrostatic equilibrium, the vertical density profile can be determined from temperature. We divided disk vertically into two parts: the high density region within $\approx 2h$ above midplane (disk portion, defined later), and the hotter low density atmosphere. In between we use a sinusoidal function to smooth out the temperature transition. The overall profile is similar to Bai & Stone (2017), but we vary temperature on $z/h$ instead of $\theta$, and the actual value at disk portion and atmosphere is close to Owen et al. (2010), with central star mass extended from 0.8 $M_{\odot}$ to 1 $M_{\odot}$:

[TABLE]

The radial grid in our simulation is set from r = 0.1 to 100 in code units with logarithmic grid spacing. The $\theta$ grid extends from 0.05 to $\pi-0.05$.The azimuthal $\phi$ grids spans from [math] to $\pi/2$. The initial density and temperature profile at midplane both follow power law function with index marked as $p$ and $q$, respectively. Midplane density at $R=1$ is set to be unity in code unit. The initial disk profile for density and temperature at midplane is:

[TABLE]

The initial vertical density profile is calculated by assuming hydrostatic equilibrium at $R-z$ plane, i.e,$v_{R}=v_{z}=0$. Following Nelson et al. (2013), we use the force balance equations at vertical and radial direction to determine $v_{\phi}$ on $R-z$ plane. The gas density in the atmosphere can be orders of magnitude lower than the midplane. In order to maintain a reasonable timestep size, a density floor that varies with location is applied (similar to Zhu & Stone (2018), with minor modification close to the pole region):

[TABLE]

where $r_{min}$ is the radius of spherical inner boundary. This density floor follows same power law slope ($p$) as the initial density profile at midplane, and gets smaller at disk atmosphere. We choose $\rho_{0,fl}=3\times 10^{-7}$ in code unit.

The initial magnetic field is set to be pure vertical. The initialization of magnetic field is described by vector potential $\bm{A}$ so that $\divergence B=0$, same as Zhu & Stone (2018). The actual field strength is:

[TABLE]

so we have $B^{2}\propto R^{p+q}$, the same as midplane gas pressure, to maintain a constant plasma $\beta$ at midplane.

The actual disk parameters are set to match the disk model used in magnetic diffusivity calculation, as p=-2.2218, q=-0.57 (see also, Okuzumi et al., 2016). Unlike that model, we do not apply exponential taper for disk density at the outer region to avoid very small density at the outer boundary. We use the same magnetic field strength as Hasegawa et al. (2017), i.e, a constant plasma $\beta=2\times 10^{4}$ at disk midplane at initialization. Following Okuzumi et al. (2016), the disk’s aspect ratio $H/R$ is 0.058 at $R=R_{0}$. The simulation domain has $128\times 64\times 16$ grid cells in $r,\theta,\phi$ directions at the root level. The grid size in $\phi$ direction is twice the value of the other two directions, as the initial setup is axisymmetric and we mainly focus on radial and vertical structure evolution. We use three levels of static mesh refinement towards the midplane, so the disk scale height at r=1 is resolved by about 9 grids at the finest level. We use open boundary condition in $r$ direction while restricting mass flux from outside into the domain. Reflective and periodical boundaries are used for $\theta$ and $\phi$ direction, respectively.

3 Results

We run our simulation for t=50$T_{0}$, with $T_{0}$ as one orbital period at r=1 ( 10 au in real world). This is equivalent of 1581 orbits at inner boundary. At this stage, our disk reaches quasi-steady state within R$\sim$ 3, so we can evaluate disk structure within the $\mathrm{CO_{2}}$ and $\mathrm{C_{2}H_{6}}$ snow lines ($<$ 30 au). The outer disk cannot reach steady state as we do not run sufficient long time and no material is fed at the outer boundary. The first thing we notice from our simulations is rings and gaps in gas surface density. In Figure 3 we see axisymmetric gaps overlay with the jumps of non-ideal MHD diffusion coefficients that are caused by $\mathrm{CO_{2}}$ and $\mathrm{C_{2}H_{6}}$ snow lines. We define the ”disk portion” as the region below the wind launching location (where $v_{\theta}$ rises significantly from zero), which is roughly 0.13 radians above the midplane at r=1.0, slightly above two scale heights. The value gets higher to larger radii as the disk is flared, but for our analysis we’ll keep this value constant as our mostly interested region is between r=1 and r=3.

The most apparent structure is the gap/ring at $r\sim 1.3$ ($\mathrm{CO_{2}}$ snow line). Comparing to a smooth disk profile, the density contrast $\approx 17\%$. This contrast is lower than $37\%\sim 78\%$ contrast we saw in surface brightness (Zhang et al., 2015). Considering dust is usually more concentrated at pressure bump, the actual dust density profile could be closer to observations.

3.1 Radial accretion variation

For a viscous accretion disk, the gas radial velocity is $v_{r,g}=-3\nu/(2r\sqrt{1-r_{*}/r})$ (Frank et al., 2002), i.e, gas drifts radially faster with a higher viscosity. If we calculate the Shakura-Sunyaev $\alpha$, defined as the $R-\phi$ stress to pressure ratio:

[TABLE]

where ${\bm{T}}=\rho v_{R}\left(v_{\phi}-\bar{v_{\phi}}\right)-B_{R}B_{\phi}$, $z-$ and $z+$ are upper and lower boundary of the simulation domain. From top right panel of Figure 4 we show $\alpha_{int}$, the $\alpha$ drops are consistent with $Am$ dips. As gas flows inward, the $v_{r,g}$ drops quickly so gas piles up, creating axisymmetric rings; then the viscosity jumps up very sharply within the $\alpha$ dip, just like the $Am$ profile, so $v_{r,g}$ increases and drains the material behind the ring, then a gap is formed.

In a magnetically coupled disk-wind system, more field diffusivity makes the disk more laminar, and wind can play a key role in driving accretion. Angular momentum transport through wind generated from disk surface can be characterized by the $z-\phi$ stress:

[TABLE]

The total accretion rate is driven by radial gradient of $R-\phi$ stress and the difference of $z-\phi$ stress between upper and lower disk surface:

[TABLE]

Here $up$ and $bot$ are upper and bottom wind base. The measured accretion rate is plotted in solid red line in the bottom left panel of Figure 4. $\dot{m}$ that is calculated from stress based on Equation 11 is shown in solid green line, which matches reasonably well with direct measurement. The drop and rise shape of $\dot{m}$ matches what is given by the $\alpha$ profile.

Note the accretion driven by wind stress is much more dominant over the viscous stress in the midplane. Viscous stress $\dot{m}$ is swinging around zero across the disk. We can also justify this by plotting $\alpha$ that only considered the disk portion (Figure 4 upper right, orange line). The picture that non-ideal MHD diffusion suppresses MRI at midplane and accretion is driven by wind is consistent with previous study (Bai & Stone, 2013).

A more direct point of view on the evolution of $\Sigma_{g}$ is measuring the change of the mass accretion rate ($-d\dot{m}/dr$) at each radius. This analysis is done in spherical polar coordinates so we can get the wind flux ($\rho v_{\theta}$) perpendicular to defined disk surface. The accretion rate is calculated by integrating $\rho v_{r}$ instead of $\rho v_{R}$. From the bottom right panel of Figure 4 we find the local wind loss and the radial change of $\dot{m_{disk}}$ show same scale of contribution to $\Sigma_{g}$. Since in our calculation, the mass of disk mesh at $r$ grows with positive $-d\dot{m}_{disk}/dr$ while declines with positive wind flux, the overlapping lines in this panel means the disk $\Sigma_{g}$ is at quasi-steady state, as mass transfer between accretion and wind almost cancels out each other at $t=50T_{0}$.

Figure 4 also raises the question of total mass budget. Modeling of previous millimeter data suggest an HL Tau disk mass of $M_{d}\approx 0.1M_{\odot}$ (ALMA Partnership et al., 2015; Kwon et al., 2015), with age about 1 to 2 Myr. In our simulation, the mass accretion rate is in a reasonable range, $\approx 2\times 10^{-7}m_{\odot}\ yr^{-1}$. The caveat is wind loss rate. By integrating ”wind flux” in bottom right panel of Figure4, we find the total wind loss rate can be up to $\approx 5.5\times 10^{-6}m_{\odot}/yr$, this is more than an order of magnitude over accretion rate. This particularly high wind flux could be due to the overly dense disk atmosphere.

3.2 Inside midplane: magnetic reconnection

To understand the accretion beneath the wind launching point, we plot the vertical cut of azimuthally averaged quantities at R=1 (blue lines) and R=1.33 (orange lines, location of $\mathrm{CO_{2}}$ snow line) in Figure 5. In the top left panel we see gas density peaks at midplane and drops sharply within the disk portion. Beyond the disk the decrease is slower as temperature gets higher. The density profile at R=1 and at $\mathrm{CO_{2}}$ snow line do not show significant difference. Same is for the magnetic pressure to thermal pressure ratio, plasma $\beta^{-1}$, at lower left panel. The upside down spike of $\beta^{-1}$ at midplane indicates strong field stretch. The top middle panel shows gas radial ( in spherical $r$ ) velocity over local sound speed ( $v_{r}/c_{s}$ ). At R=1.0, the accretion is concentrated at midplane. Wind starts at $\theta-\pi/2=0.13$ above the midplane, and became supersonic at $\theta-\pi/2=0.57$. With similar $v_{r}$ structure in the wind region, $v_{r}$ in disk portion is quite different at R=1.33. At midplane we see decretion, and fast ( up to 0.1 local sound speed) accretion happening at disk surface, which is similar to flow pattern in global ideal MHD simulations Zhu & Stone (2018). The similar trend can be seen in radial mass flux ($\rho v_{r}$) plot at the the lower middle panel. The $R\phi$ and $z\phi$ components of total stress is plotted in rightmost two panels. The $R\phi$ stress at both radius are close to zero at midplane, as already shown by very low midplane $\alpha$ in Figure 4. Almost zero $R\phi$ stress at different radius also suggesting that midplane accretion or decretion is not driven by $R\phi$ component. As calculated in Section 3, vertical ($z$) gradient of ${\bm{T}}_{z\phi}$ also contributes to $\dot{m}$. In the $z\phi$ stress panel, $d{\bm{T}}_{z\phi}/dz$ maximizes at midplane, which is in agreement with the $\rho v_{r}$ spike in the lower middle panel. At R=1.33 (13.3 au), the midplane decretion is driven by $z\phi$ stress’s negative $z$ gradient. The surface accretion is related to the two positive $z$ gradient parts around disk surface. The differences of ${\bm{T}}_{z\phi}$ at these two radius also suggest very different magnetic field structure.

The origin of the differences in stress structure can be understood from the 2D plot on $z-R$ plane, showing how magnetic field is evolved in Figure 6. At the beginning of simulation, the magnetic field only has vertical ($z$) component, so $B_{\phi},B_{R}=0$ everywhere. Apart from the dominant Keplerian shear ($dv_{\phi}/dR$) at midplane, the disk is also differentially rotating vertically ($dv_{\phi}/dz$), as midplane being the fastest, and the atmosphere portion is mostly pressure gradient supported, rotates much slower. After simulation starts, the field lines closer to midplane are dragged towards the rotation direction, making the $B_{\phi}$ positive below midplane and negative above midplane. The motion along the $r$ direction is even more complicated. From t=0.5, at $\mathrm{CO_{2}}$ snow line, the velocity shear between wind and surface accretion forms a vortex like flow pattern, that drives gas infall to the midplane. This gas motion drags field lines towards midplane, making a kinked form of magnetic field. Accretion at midplane also drags field lines inward, amplifying horizontal field component $B_{R}$, and the polarity of $B_{R}$ changes sign in the accretion stream. Close to midplane, Keplerian shear stretches horizontal field lines, just like the growth of azimuthal fields in the linear phase of MRI. The process makes the field kink at midplane even sharper. The vertical gradient of $B_{\phi}$ component tends to get steeper with AD (Brandenburg & Zweibel, 1994). The strength of AD is proportional to $B^{2}$, so it minimizes at magnetic null (for $B_{\phi}$ component, and $B_{r}$, $B_{\theta}$ varies slowly in $z$ direction at midplane) where magnetic field changes sign at midplane, stopping AD from smoothing out the sharp field structure.

The steep magnetic gradient leads to large $dT_{z\phi}/dz$ at midplane, drives a highly concentrated accretion sheet. An alternative physical picture is the the removal of angular momentum at midplane by Lorentz force. Since $J_{r}\sim\partial(B_{\phi}\sin\theta)/\partial\theta$, and $B_{\phi}$ changes sign at midplane, there is a infinite thin current sheet at the disk midplane. The Lorentz force $\propto J_{r}B_{\theta}$ then maximizes at midplane and drives accretion at its highest rate. The accretion flow now drags field lines inward, pinching the magnetic field radially. At R=1.33 (13.3 au, $\mathrm{CO_{2}}$ snow line), the kinked field lines, together with field lines dragged by infalling gas towards midplane, eventually reconnect. The poloidal component $B_{\theta}$, or vertical component $B_{z}$ that reconnects now have their direction reversed. This changes sign of $dT_{z\phi}/dz$ at midplane, driving mass flux outward. This outward moving mass can also be seen in the mass flux panel (left) in Figure 7. The outward mass flux concentrates gas just outside the field reconnection area, forming a gaseous ring.

The middle panel of Figure 7 also shows that the magnetic field does not always changes its polarity exactly at midplane. From R=2.2 to R=2.8, there is an large triangular shaped structure of $B_{\phi}$ anomaly. $B_{\phi}$ usually changes from positive to negative when moving from below to above midplane, but at the anomaly, the sequence of field polarity is positive-negative-positive-negative. This changes direction of ${\bm{T}}_{z\phi}$ thus alters the direction of mass flux. This zone of decretion covers a relatively large area close to midplane, while the outward mass flow in the field loop of magnetic reconnection ($\mathrm{CO_{2}}$ snow line) is confined in a thin sheet, as plotted in left panel in Figure 7. The formation of this magnetic anomaly is shown in Figure 6. In this area, there’s a quick accretion flow above the disk surface. As mentioned before, $B_{R}$ changes sign in accretion stream , with Keplerian shear, $B_{\phi}$ also changes sign. This overall $B_{\phi}$ structure sustains over time, though surface accretion is just a transient behavior in the early phase.

4 Parameter study

In order to assess how robust the ring formation mechanism in the reference simulation is, especially how the diffusion structure contribute to ring build up, we performed several additional simulations. The first one is a disk with ”featureless” diffusion profile, i.e, we use constant $Am$ and $\eta_{\mathrm{O}}$ at midplane for the whole disk. The exact values are set to match the diffusion strength at R=1.16 in the reference run, with $Am=0.23$ and $\eta_{\mathrm{O}}=2.2\times 10^{15}cm^{2}s^{-1}$ ($1.56\times 10^{-5}$ in code units).

The results of ”featureless” model are shown as solid green line in Figure 8 and leftmost panels in Figure 9. The rings and gaps at $\mathrm{CO_{2}}$ and $\mathrm{C_{2}H_{6}}$ snow lines disappeared with smooth diffusion profiles. There is however a density bump at R=1.8. Judging from Figure 9, we see very similar $B_{\phi}$ and mass flow pattern as R=2.2 to R=2.8 in the reference run. This density bump has the same origin as the one at R=2.2 in the reference model. Another significant difference in Figure 9 is the absence of magnetic reconnection for the featureless run. This is similar to Suriano et al. (2018, 2019), as reconnection only happens when AD is strong enough at midplane. The field reconnection in the reference model is very likely to prevent surface accretion from dragging field lines inwards as much as in the featureless run. This moves the $B_{\phi}$ anomaly outward, thus the ”extra” density bump changes from R=1.8 to R=2.2 in the reference model.

The second and third additional simulations are ”AD only” and ”Ohmic only” runs. We use the same setup as the reference run but only include Ambipolar diffusion or Ohmic resistivity in these simulations. The ”AD only” model is almost identical to the reference model in both 1D $\Sigma_{g}$ profile and 2D field and flow pattern. The ”Ohmic only” case shows barely any significant feature in disk surface density and a vastly different structure in Figure 9. Accretion is dominated by fast flows at disk surface and disk midplane is filled with chaotic patterns of both accretion and decretion. The majority of accretion not being at midplane makes the diffusion features at snow line give little contribution to the whole disk.

In the 4th test we simply make magnetic field $\sqrt{10}$ times stronger than the reference model, setting plasma $\beta=2\times 10^{3}$. Note we keep $Am$ and $\eta_{\mathrm{O}}$ unchanged so the strength of magnetic diffusion remains the same as reference model. The result is significant thicker rings and deeper gaps at nearly all locations ( Figure 8). This is expected as magnetic stress is 10 times larger and so its radial and vertical gradients. The ”extra” density bump at R=2.2 does not show up in this case. The relevant $B_{\phi}$ anomaly and flow pattern are also missing in $\beta=2\times 10^{3}$ panels in Figure 9. With plasma $\beta$ 10 times smaller, the magnetic field is so strong that the initial surface accretion is not able to drag and bend it to make the $B_{\phi}$ anomaly.

In the last test, ”poloidal field”, we do not start with a pure vertical field. Poloidal magnetic fields are initialized with vector potential generalized from Zanni et al. (2007):

[TABLE]

where $m$ is a parameter that specifies the degree that poloidal fields bend, with $m$ goes infinity giving a pure vertical field. We chose $m=0.5$ the same as Bai & Stone (2017). The initially bent field lines follows the gas flow better than vertical field, and is not easily dragged by fast surface inflow when the simulation starts. This setup gets rid of the $B_{\phi}$ anomaly. In Figure 8, the surface density profile of ”poloidal field” almost overlaps with our reference run, with the absence of the ”extra” density bump at $R=2.2$. It advises that MHD simulation with moderate AD and magnetic field may better start with a poloidal field instead of pure vertical field.

4.1 Observation signature

Observational tests for the derived disk structure in our simulations include: 1) studying dust and gas surface density across the gaps (Figure 8), 2) probing the azimuthally averaged gas velocity structure (Figure 7) using the method in Teague et al. (2018), 3) measuring the disk magnetic field orientation and strength directly.

One exciting recent observational constraint is from using the molecular lines’ Zeeman effect to study disk magnetic fields directly. Recent ALMA circular polarization measurements for CN lines in TW Hya suggest that, at 42 au, the upper limit of $B_{z}$ and $B_{\phi}$ are $0.8\ \mathrm{mG}$ and $30\ \mathrm{mG}$, respectively (Vlemmings et al., 2019). Chemistry modelling shows CN are most abundant at z/R between 0.3 and 0.4 beyond 10 au (Cazzoletti et al., 2018). In Figure 10 we measured poloidal (vertical) and toroidal magnetic field at z=h, 2h and 5h above midplane. At r=50 au, z=5h gives $B_{z}=1\ \mathrm{mG}$ and $B_{\phi}=8\ \mathrm{mG}$, which are roughly consistent with the derived upper limit. This also suggests that ALMA Zeeman measurements start to constrain MHD disk models. The future circular polarization measurements for molecular lines will provide stringent tests to disk accretion models.

We have shown the field strength in our reference run. To scale the field strength to systems with other accretion rate, we measured $\dot{m}$ of $\beta=2\times 10^{3}$ case. It is $\approx 1.5\times 10^{-6}m_{\odot}\ yr^{-1}$, 7 times higher than the reference run. This gives $\dot{M}\propto|B|^{1.7}$, and it is between the wind driven accretion measurement from Bai (2013) (almost proportioanl to $|B_{z}|$), and the MRI driven accretion in Okuzumi et al. (2014) (proportional to $B_{z}^{2}$).

5 Summary

In this study, we have carried out a set of 3D non-ideal MHD simulations of magnetically coupled disk-wind systems with different profiles of ambipolar diffusion and ohmic resistivity. Our main conclusions are as follows.

(i) snow line in HL Tau make dust aggregates pile up in the sintering zones due to the combined effect of radial drift and sintering-induced fragmentation. The concentration of smaller dust leads to stronger ambipolar diffusion and Ohmic resistivity in these areas. By varying magnetic diffusion coefficients we can generate gaseous rings and gaps in a weakly ionized gas disk with moderate magnetic field. Generally speaking, varying diffusion strength changes accretion structure, making material pile up radially. This is a new dust feedback mechanism. Comparing to aerodynamic drag, it requires much less dust to take effect.

(ii) The correlation between snowline locations and rings/gaps can be complex. The dust model we adopted is based on smooth gas disk, while the rings produced in our simulations can trap dust and change the overall profile of $\Sigma_{d}$. In Figure 8 the snowlines are located a bit outside the deepest point in the gap. The width of gaseous rings in the reference run is overall narrower than rings in the dust model. As dust are being trapped at pressure maxima, the radius with strongest diffusion (also $\Sigma_{d}$ maxima) can move outward, towards the gaseous ring. If the zone of magnetic reconnection relocates at the pressure maxima, it reduces both dust and gas surface density, slows down ring formation. The exact strength of this gas to dust back reaction still needs further investigation. One factor that complicates this scenario is the ones are most concentrated at pressure bumps due to radial drift are larger dust particles, while in ionization calculation it is the smallest grains that dominates the total electron absorbing cross section.

(iii)With very strong ambipolar diffusion ($Am\sim 0.01$), the gaseous rings and gaps can be generated from magnetic reconnection at midplane. High AD strength confines accretion at midplane and the vertical differentiation of radial velocity drags field lines so close that they reconnect. The reconnected field loop alters the direction of midplane mass flow. This drains mass inside the reconnection zone and carves gaps in the disk. The mass accumulated outside the field loop form dense rings.

(iv) Compared to ambipolar diffusion, Ohmic resistivity is almost irrelevant to ring formation at outer part of disk. In the ”Ohmic only” case, the disk is better magnetically coupled and behaves more like ideal MHD. The accretion is dominated by mass flow at disk surface instead of midplane in the AD simulation. One contribution factor could be that Ohmic resistivity is only strong at high density region. It declines rapidly from midplane to disk surface. A probably more important distinction between Ohmic resistivity and ambipolar diffusion is the former does not rely on the mutual angle between current and magnetic field, so the diffusion from Ohmic is isotropic. Ohmic resistivity is less effective in modifying magnetic field and accretion structure.

(v) With moderate magnetic field strength ($\beta=2\times 10^{4}$), AD dominated non-ideal MHD simulation can be sensitive to initial profile of magnetic field. Fast surface flow drags field lines radially, and then the radial pinch is twisted by vertically differentiated rotation. Eventually this process leaves a $B_{\phi}$ anomaly at midplane, and a circular like flow pattern is sustained. This flow pattern allows the formation of an ”extra” ring between the $\mathrm{CO_{2}}$ and $\mathrm{C_{2}H_{6}}$ snow line. Stronger magnetic field will not be affected by initial surface velocity shear. A poloidal field that is bent away from rotation axis follows the fast surface flow better than pure vertical field, and is also barely affected. We’d like to point out the fact that the ”extra” ring is related to specific initial setup does not rule out this ring formation mechanism in real astrophysical process. Formation of certain rings may not be related to specific diffusion structure.

(vi) The wind loss rate in our simulations is exceptionally high. The vertical structure of the disk, including both temperature and non-ideal diffusion, are simple prescribed profiles. Wang et al. (2019) presented first MHD disk simulation with consistent thermochemistry. Temperature at disk atmosphere and magnetic diffusivities can be obtained directly in real time. This would be a nice direction of improvement for furture study.

Overall, our work suggests that changes in the disk ionization structure can robustly lead to the change of the disk accretion structure, which leads to gaseous gaps and rings. In our particular example, the change of the dust size distribution across the snowline affects the non-ideal MHD effects and eventually leads to gaseous gaps and rings. Our work also suggests that, besides these changes of the disk ionization structure, non-ideal MHD effects in smooth disks can also lead to structure formation but these structures sensitively depends on the numerical setup and different non-ideal MHD effects involved.

Acknowledgements

We thank Mordecai Mac-Low, Phil Armitage, Jeremy Goodman and Zhenyu Wang for helpful discussions, and Jiming Shi for help on the code. X. H. and Z. Z. acknowledge support from the National Aeronautics and Space Administration through the Astrophysics Theory Program with Grant No. NNX17AK40G and Sloan Research Fellowship. Simulations are carried out with the support from the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG-AST130002. S.O. is supported by JSPS KAKENHI Grant Numbers JP16K17661, JP18H05438, and JP19K03926. X.-N.B acknowledges support from the Youth Thousand Talent Program of China.K.T. is supported by JSPS grants JP16H05998, JP16K13786 and JP17KK0091. K.T. acknowledges support by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan as Exploratory Challenge on Post-K computer (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar System).