Nonideal MHD Simulation of HL Tau Disk: Formation of Rings

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 nonideal MHD simulations including effects from both ohmic resistivity and ambipolar diffusion (AD) to model the HL Tau disk. The nonideal 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 the 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 the snow line. Furthermore, AD leads to highly preferential accretion at the midplane, followed by magnetic reconnection. This results in a local zone of decretion that drains the mass in the field reconnection area, which leaves a gap and an adjacent ring just outside of it. Overall, under favorable conditions, both snow lines and nonideal MHD effects can lead to gaseous gaps and rings in protoplanetary disks.


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 ∼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 that 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., Dipierro et al. 2015;Dong et al. 2015), or a single planet opening multiple gaps .
Another possibility is the condensation fronts of various solids altering dust growth. The large dispersion of location and size of rings from recent observations show low correlation to snow lines or resonance orbits of planets 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., H O 2 and NH 3 . A more detailed study on dust coagulation by Okuzumi et al. (2016) shows that the 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 a low-turbulence inner disk, making a ringlike structure (Takahashi & Inutsuka 2016). Dust to gas back reaction and gas compression can cause a large gradient of dust radial drift velocity; this self-induced dust pile up is discussed in Gonzalez et al. (2017). If planets are formed early, dust rings can be the collisional outcome of planetesimals perturbed by planets and orbital resonances (Boley 2017). Magnetic selforganization is seen in the form of "zonal flows," resulting in radial variations in density and magnetic flux. These banded rings can be sustained as the pressure gradient is balanced by the Coriolis force (Johansen et al. 2009;Bai 2014). Recent MHD simulations suggest that the 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(Suriano et al. , 2018(Suriano et al. , 2019, when nonideal MHD diffusivities like ohmic resistivity and ambipolar diffusion (AD) 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 the 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 (submicron) grains, thus affecting the coupling between gas and magnetic fields (e.g., Wardle 2007;Bai 2011). This inspired the idea, which 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 variation of accretion rate and lead to axisymmetric structure in gas surface density.
In this study we first generated a nonideal MHD diffusivity profile from a dust coagulation model including the sintering effect. Then we implement the diffusion profile in a global 3D nonideal 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.

Nonideal MHD Diffusion
We solve the magnetohydrodynamic (MHD) equations using Athena++ (J. M. Stone et al. 2019, in preparation) with nonideal magnetic diffusion terms (Bai 2017a(Bai , 2017b. The basic nonideal induction equations are where v is gas velocity, B is magnetic field, | | = b B B is the unit vector representing field line direction. J is current density vector, withĴ as the current component perpendicular to the magnetic field. η O , η H , and η A are ohmic, Hall, and ambipolar diffusivities, and here for simplicity, we only consider ohmic resistivity and AD. The Hall diffusivity significantly complicates the picture, but they are more applicable to the inner disk (Béthune et al. 2016;Bai 2017b). Modeling HL Tau does not resolve the inner disk, so our simulation mainly focused on the outer part of the 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 min and a max (Equation (20) of Okuzumi et al. 2016). The minimum grain size a min is fixed to 0.1μm, whereas the maximum size a max and total dust surface density Σ d are taken from a snapshot of a dust evolution calculation that takes into account the low stickiness of CO 2 ice (Okuzumi & Tazaki 2019) and aggregate sintering (Okuzumi et al. 2016). The fragmentation threshold velocity of nonsintered aggregates is assumed to be 0.73ms −1 inside the H 2 O snow line, 7.3ms −1 between the H 2 O and CO 2 snow lines, and 0.73ms −1 outside the CO 2 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 max and Σ d are plotted in Figures 1(a) and (b). Because of sintering, dust particles behind the snow lines CO 2 , C 2 H 6 , and CH 4 experience enhanced collisional fragmentation and pile up there. The loss of H 2 O ice interior to the H 2 O snow line also enhances Figure 1. Radial profiles of the maximum grain size a min (panel (a)), total dust surface density Σ d (panel (b)), and midplane ionization rate ζ (panel (c)) used in the ionization calculation, as well as the profiles of the electron and ion abundances, n e /n and n i /n, obtained from the ionization calculation (panel (d)). The vertical lines mark the snow lines of, from left to right, H 2 O, CO 2 , C 2 H 6 , and CH 4 . For comparison, the ion and electron abundances in the dust-free case are also shown by the dotted curve in panel (d).
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 semianalytical ionization model of Okuzumi (2009). Our ionizing sources include the short-lived radionuclide 26 Al (Umebayashi & Nakano 2009) and stellar X-rays (Igea & Glassgold 1999;Bai & Goodman 2009) with an X-ray luminosity of 10 31 erg s −1 . The radial profile of the ionization rate per hydrogen molecule at the midplane, ζ, is plotted in Figure 1(c). At R15au, X-rays are able to penetrate down to the midplane and produce a high value of ζ. All positive ions are represented by HCO + , with its recombination rate coefficient given by 2.4×10 −7 (T/300 K) −0.69 cm 3 s −1 (Ganguli et al. 1988). Only for dust evolution and ionization calculations, we assume the gas disk model and T=310(R/1 au) −0.57 with M disk =0.2 M e and R c =150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 dustfree case is also shown in Figure 1(d)). At R15au, the effect of dust grains is relatively minor because of the elevated ionization rate. 8 The second step is to convert the midplane electron and ion abundances into the midplane ohmic and ambipolar diffusivities, η O (R, 0) and η 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 η O (R, 0) and η A (R, 0) for B=1mG. While η O is independent of B, η 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 ≈n e holds in most parts of the disk at R4au (see Figure 1(d)), we assume η A ∝B 2 when we compute η A (R, 0) for B¹1mG. Because of the B 2 dependence, the dimensionless ambipolar Elsässer number, is better suited to describe the strength of AD. Here pr = v B 4 a is the Alfvén velocity, and Ω 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 the water snow line. To ensure a reasonable timestep, we focus on the CO 2 and C 2 H 6 snow lines and remove any feature within the 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 to Bai (2017a). The strength of both diffusion effects drops quickly when getting away from disk portion, the difference is that Am remains constant within 2h, and η O ʼs decrease starts at midplane, ,0 0 . 5 91 t a n h 10 0.06 Here h is disk scale height at midplane. These profiles ensure the strongest diffusion in the disk midplane and drops gradually over the atmosphere. The AD Elsässer number Am is almost constant until ≈2h above disk midplane, making AD more dominant over ohmic resistivity at the disk surface, where density is lower.

Disk Setup
To minimize grid noise of the radially flowing disk wind, we perform our simulation in spherical polar coordinates (r, θ, f). In this paper, we transform some quantities to cylindrical coordinates (R, f, z) for analysis. Note that R=1 equals 10 au in physical units. For an accretion disk with the vertical hydrostatic equilibrium, the vertical density profile can be determined from temperature. We divided the disk vertically into two parts: the high density region within ≈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 (2017a), but we vary temperature on z/h instead of θ, 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 e to 1 M e : The radial grid in our simulation is set from r=0.1 to 100 in code units with logarithmic grid spacing. The θ grid extends from 0.05 to π−0.05. The azimuthal f grid spans from 0 to π/2. The initial density and temperature profile at midplane both follow a power-law function with index marked as p and ) and ambipolar diffusion ( ) at disk midplane (assuming constant plasma β=2×10 4 at midplane). In the radial direction, from inside out are H 2 O, CO 2 , C 2 H 6 , CH 4 , and CO snow lines. The diffusion profile produced by chemical calculation are plotted in dashed lines, and the solid lines are implemented in our MHD simulation. Note that we removed any feature inside the CO 2 snow line. The analytical functions of vertical diffusion profiles are described in Equation (3). q, respectively. Midplane density at R=1 is set to be unity in code units. The initial disk profile for density and temperature at midplane is 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 directions to determine v f on the 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): where r min is the radius of the spherical inner boundary. This density floor follows the same power-law slope (p) as the initial density profile at midplane, and gets smaller at disk atmosphere. We choose ρ 0,fl =3×10 −7 in code units. The initial magnetic field is set to be purely vertical. The initialization of the magnetic field is described by vector potential A so that ·  = B 0, same as Zhu & Stone (2018). The actual field strength is so we have B 2 ∝R p+ q , the same as midplane gas pressure, to maintain a constant plasma β at midplane. The actual disk parameters are set to match the disk model used in the 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 β=2×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×64×16 grid cells in r, θ, f directions at the root level. The grid size in the f 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 toward the midplane, so the disk scale height at r=1 is resolved by about nine grids at the finest level. We use the open boundary condition in the r direction while restricting mass flux from outside into the domain. Reflective and periodical boundaries are used for θ and f directions, respectively.

Results
We run our simulation for t=50T 0 , with T 0 as one orbital period at r=1 (10 au in real world). This is equivalent to 1581 orbits at the inner boundary. At this stage, our disk reaches a quasi-steady state within R∼3, so we can evaluate disk structure within the CO 2 and C H 2 6 snow lines (<30 au). The outer disk cannot reach steady state as we do not run for a sufficiently 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 nonideal MHD diffusion coefficients that are caused by CO 2 and C H 2 6 snow lines. We define the "disk portion" as the region below the wind launching location (where v θ 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 with larger radii as the disk is flared, but for our analysis we keep this value constant as our region of interest is between r=1 and r=3.
The most apparent structure is the gap/ring at r∼1.3 (CO 2 snow line). Comparing to a smooth disk profile, the density contrast ≈17%. This contrast is lower than the 37%∼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.

Radial Accretion Variation
For a viscous accretion disk, the gas radial velocity is Frank et al. 2002), i.e., gas drifts radially faster with a higher viscosity. If we calculate the Shakura-Sunyaev α, defined as the R-f stress to pressure ratio, , z−and z+ are upper and lower boundaries of the simulation domain. From the top right panel of Figure 4 we show α int ; the α 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 α 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 the disk surface can be characterized by the z-f stress: The total accretion rate is driven by a radial gradient of R-f stress and the difference of z-f stress between upper and lower disk surfaces: Here, up and bot are upper and bottom wind bases. The measured accretion rate is plotted by the solid red line in the bottom left panel of Figure 4.ṁ, calculated from stress based on Equation (11), is shown by the solid green line, which matches reasonably well with direct measurement. The drop and rise shape ofṁ match what is given by the α profile. Note the accretion driven by wind stress is much more dominant over the viscous stress in the midplane. Viscous stressṁ is swinging around zero across the disk. We can also justify this by plotting α that only considered the disk portion (Figure 4 upper right, orange line). The picture that nonideal MHD diffusion suppresses MRI at midplane and accretion is driven by wind is consistent with a previous study (Bai & Stone 2013).
A more direct point of view on the evolution of Σ g is measuring the change of the mass accretion rate (-dm dr) at each radius. This analysis is done in spherical polar coordinates so we can get the wind flux (ρv θ ) perpendicular to defined disk surface. The accretion rate is calculated by integrating ρv r instead of ρv R . From the bottom right panel of Figure 4 we find that the local wind loss and the radial change ofṁ disk show the same scale of contribution to Σ g . Since in our calculation the mass of disk mesh at r grows with positive-dm dr disk and declines with positive wind flux, the overlapping lines in this panel indicate that the disk Σ g is at quasi-steady state, as mass transfers between accretion and wind almost cancel 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 ≈0.1 M e (ALMA Partnership et al. 2015;Kwon et al. 2015), with an age of about 1-2 Myr. In our simulation, the mass accretion rate is in a reasonable range, ≈2×10 −7 m e yr −1 . The caveat is wind loss rate. By integrating "wind flux" in the bottom right panel of Figure 4, we find that the total wind loss rate can be up to ≈5.5×10 −6 m e yr −1 , 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.

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 CO 2 snow line) in Figure 5. In the top left panel we see that 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 profiles at R=1 and at the CO 2 snow line do not show significant differences. The same goes for the magnetic pressure to thermal pressure ratio, plasma β −1 , in the lower left panel. The upside down spike of β −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 θ−π/2=0.13 above the midplane, and becomes supersonic at θ−π/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 Figure 4. Radial profiles (from 1 to 100 au) of azimuthally and time averaged quantities (from t=49 to t=50 T 0 with Δt=0.1) over radius, with Am overplotted as a gray dashed line on the right y-axis in the left two panels. Top left: gas surface density Σ g , in gcm −2 . Top right: vertically integrated α, with the blue line as the total value, and the orange line containing only the disk portion. The green dashed line is the vertically averaged β within the disk portion, shown on the right y-axis. The three gray vertical lines mark the location of Σ g maxima ("rings") measured in the top left panel. Bottom left: measuredṁ vs. calculation from stress, with the red solid line as the measured midplaneṁ. The blue dashed line isṁ calculated from wind stress, the orange dashed line isṁ calculated from viscous stress, and the green solid line is the sum of these two. Bottom right panel shows quantities measured in spherical polar coordinates. The blue dashed line is Σ g but integrated in the θ direction within the disk portion (in gcm −2 , y-axis on the left). The blue solid line is the radial gradient of mass flux in the r direction, within the disk portion, and the orange solid line is local wind flux perpendicular to the disk surface (wind base). These two quantities are marked by the right y-axis, inṁ yr −1 au −1 . Note that in this plot, positive-dm dr disk means the disk gains mass at this radii, while positive wind flux means the disk loses mass. local sound speed) accretion happening at disk surface, which is similar to the flow pattern in global ideal MHD simulations (Zhu & Stone 2018). A similar trend can be seen in the radial mass flux (ρv r ) plot in the lower middle panel. The Rf and zf components of total stress are plotted in the rightmost two panels. The Rf stress at both radii are close to zero at midplane, as already shown by very low midplane α in Figure 4. Almost zero Rf stress at different radii also suggests that midplane accretion or decretion is not driven by the Rf component. As calculated in Section 3, the vertical (z) gradient of f T z also contributes toṁ. In the zf stress panel, maximizes at midplane, which is in agreement with the ρv r spike in the lower middle panel. At R=1.33 (13.3 au), the midplane decretion is driven by zf stress's negative z gradient. The surface accretion is related to the two positive z gradient parts around disk surface. The differences of f T z at these two radii also suggest very different magnetic field structure.
The origin of the differences in stress structure can be understood from the 2D plot on the z-R plane, showing how magnetic field is evolved in Figure 6. At the beginning of the simulation, the magnetic field only has a vertical (z) component, so B f , B R =0 everywhere. Apart from the dominant Keplerian shear (dv f /dR) at midplane, the disk is also differentially rotating vertically (dv f /dz), as midplane is the fastest and the atmosphere portion, which is mostly pressure gradient supported, rotates much slower. After the simulation starts, the field lines closer to midplane are dragged toward the rotation direction, making the B f positive below midplane and negative above midplane. The motion along the r direction is even more complicated. From t=0.5, at the 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 toward 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 the B f 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 f component, and B r , B θ varies slowly in the 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 zf /dz at midplane, and drives a highly concentrated accretion sheet. An alternative physical picture is the removal of angular momentum at midplane by Lorentz force. Since r and B f changes sign at midplane there is a infinite thin current sheet at the disk midplane. The Lorentz force ∝J r B θ 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, CO 2 snow line), the kinked field lines, together with field lines dragged by infalling gas toward midplane, eventually reconnect. The poloidal component B θ , or vertical component B z , that reconnects now has its direction reversed. This changes the sign of dT zf /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 shows that the magnetic field does not always change its polarity exactly at midplane. From R=2.2 to R=2.8, there is a large triangular structure of the B f anomaly. B f usually changes from positive to negative when moving from below to above midplane, but at the anomaly, the sequence of field polarity is positive-negativepositive-negative. This changes the direction of f T z and 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 (CO 2 snow line) is confined in a thin sheet, as plotted in the left panel in Figure 7. The formation of this magnetic anomaly is shown in Figure 6. In this area, there is a quick accretion flow above the disk surface. As mentioned before, B R changes sign in the accretion stream, and with Keplerian shear B f also changes sign. This overall B f structure sustains over time, though surface accretion is just a transient behavior in the early phase.

Parameter Study
In order to assess how robust the ring formation mechanism in the reference simulation is, especially how the diffusion structure contributes to ring build up, we performed several  additional simulations. The first one is a disk with a "featureless" diffusion profile, i.e., we use constant Am and η 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 η O =2.2×10 15 cm 2 s −1 (1.56×10 −5 in code units).
The results of the "featureless" model are shown as a solid green line in Figure 8 and in the leftmost panels in Figure 9. The rings and gaps at CO 2 and 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 f 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. (2018Suriano et al. ( , 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 f 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 Σ 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 causes the diffusion features at the snow line to give little contribution to the whole disk.
In the fourth test we simply make magnetic field 10 times stronger than the reference model, setting plasma β=2×10 3 . Note that we keep Am and η O unchanged so the strength of magnetic diffusion remains the same as that in the reference model. The result is significantly thicker rings and deeper gaps at nearly all locations (Figure 8). This is expected as magnetic stress is 10 times larger, as are its radial and vertical gradients. The "extra" density bump at R=2.2 does not show up in this case. The relevant B f anomaly and flow pattern are also missing in β=2×10 3 panels in Figure 9. With plasma β 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 f 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): where m is a parameter that specifies the degree at which poloidal fields bend, with m of infinity giving a pure vertical field. We chose m=0.5 the same as Bai (2017a). The initially bent field lines follow the gas flow better than the vertical field, and are not easily dragged by fast surface inflow when the simulation starts. This setup gets rid of the B f anomaly. In Figure 8, the surface density profile of the "poloidal field" almost overlaps with our reference run, with the absence of the "extra" density bump at R=2.2. It shows that an MHD simulation with moderate AD and magnetic field may start better with a poloidal field instead of pure vertical field.

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), and (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 f are 0.8 mG and 30 mG, respectively (Vlemmings et al. 2019). Chemistry modeling shows that 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 fields at z=h, 2h, and 5h above midplane. At r=50 au, z=5h gives B z =1 mG and B f =8 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 rates, we measuredṁ of the β=2×10 3 case. It is ≈1.5×10 −6 m e yr −1 , 7 times higher than the reference run. This gives˙| | µ M B 1.7 , and it is between the wind driven accretion measurement from Bai (2013; almost proportional to | | B z ), and the MRI driven accretion in Okuzumi et al. (2014; proportional to B z 2 ).

Summary
In this study, we have carried out a set of 3D nonideal MHD simulations of magnetically coupled disk-wind systems with Figure 8. Gas surface density variation (Σ g /Σ g,initial ) of six different models, divided into two bottom panels for better illustration. The top panel shows dust surface density adopted for ionization calculation. The vertical lines mark the location of CO 2 and C 2 H 6 snow lines. Note the "extra" density bump is not featured in cases of "poloidal field" and "β=2×10 3 . Also note the relative location of gaseous rings and dust rings. The deepest point of gap is located a bit inside the snow lines, as the region of decretion caused by magnetic reconnection extends a bit inward. different profiles of AD and ohmic resistivity. Our main conclusions are as follows.
(i) The snow line in HL Tau makes 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 AD 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 snow line locations and rings/ gaps can be complex. The dust model we adopted is based on a smooth gas disk, while the rings produced in our simulations can trap dust and change the overall profile of Σ d . In Figure 8 the snow lines 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 is being trapped at pressure maxima, the radius with the strongest diffusion (also Σ d maxima) can move outward, toward the gaseous ring. If the zone of magnetic reconnection relocates at the pressure maxima, it reduces both dust and gas surface density and 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 that the grains most concentrated at pressure bumps due to radial drift are larger dust particles, while in the ionization calculation it is the smallest grains that dominate the total electron absorbing cross section.
(iii) With very strong AD (Am∼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 AD, ohmic resistivity is almost irrelevant to ring formation at the outer part of the 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 the disk surface instead of the midplane in the AD simulation. One contribution factor could be that ohmic resistivity is only strong at the high density region. It declines rapidly from midplane to disk surface. A likely more important distinction between ohmic resistivity and AD is that 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 (β=2×10 4 ), AD-dominated nonideal MHD simulation can be sensitive to the initial profile of a 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 f anomaly at midplane, and a circular-like flow pattern is sustained. This flow pattern allows the formation of an "extra" ring between the CO 2 and C 2 H 6 snow lines. Stronger magnetic field will not be affected by initial surface velocity shear. A poloidal field that is bent away from the rotation axis follows the fast surface flow better than pure vertical field, and is also barely affected. We would like to point out the fact that the "extra" ring related to the specific initial setup does not rule out this ring formation mechanism in a 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 nonideal 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 future 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 snow line affects the nonideal MHD effects and eventually leads to gaseous gaps and rings. Our work also suggests that, besides these changes of the disk ionization structure, nonideal MHD effects in smooth disks can also lead to structure formation but these structures sensitively depends on the numerical setup and different nonideal MHD effects involved.