The dimming of RW Auriga. Is dust accretion preceding an outburst?
Matías Gárate, Til Birnstiel, Sebastian Markus Stammler, Hans Moritz Günther
DDraft version January 9, 2019
Typeset using L A TEX twocolumn style in AASTeX62
The dimming of RW Auriga. Is dust accretion preceding an outburst?
Mat´ıas G´arate, Til Birnstiel, Sebastian Markus Stammler, and Hans Moritz G¨unther University Observatory, Faculty of Physics, Ludwig-Maximilians-Universit¨at M¨unchen, Scheinerstr. 1, 81679 Munich, Germany MIT, Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Avenue, Cambridge, MA 02139, USA (Received 2018 October 12; Revised 2018 November 23; Accepted 2018 November 24)
Submitted to ApJABSTRACTRW Aur A has experienced various dimming events in the last years, decreasing its brightness by ∼ M d = 6 × − M (cid:12) / yr and reaching dust-to-gas ratiosof (cid:15) ≈
5, preceding an increase in the gas accretion by a few years. This sudden rise of dust accretioncan provide the material required for the dimmings, although the question of how to put the dust intothe line of sight remains open to speculation.
Keywords: accretion, accretion disks — hydrodynamics — protoplanetary disks — stars: individual:RW Aur A INTRODUCTIONRW Aur A is a young star that in the last decade pre-sented unusual variations in its luminosity. The star hasabout a solar mass, it is part of a binary system, andis surrounded by a protoplanetary disk showing signa-tures of tidal interaction (Cabrit et al. 2006; Rodriguezet al. 2018). The star had an almost constant luminos-ity for around a century, interrupted only by a few shortand isolated dimmings (see Berdnikov et al. (2017) fora historical summary) until 2010, when its brightnesssuddenly dropped by 2 mag in the V band for 6 months(Rodriguez et al. 2013). Since 2010, a total of five dim-ming events have been recorded (see Rodriguez et al.2013, 2016; Petrov et al. 2015; Lamzin et al. 2017; Berd-nikov et al. 2017, among others). The dimmings can lastfrom a few months to two years, and reduce the bright-ness of the star up to 3 mag in the visual. Moreover,there is no obvious periodicity in their occurrence, and
Corresponding author: Mat´ıas G´[email protected] their origin is not yet clear (a summary of the eventscan be found in Rodriguez et al. 2018).1.1.
Observations of RW Aur Dimmings
Some observations in the recent years have shed lighton the nature of RW Aur A dimmings:During the event in 2014-2015 (Petrov et al. 2015), ob-servations by Shenavrin et al. (2015) show a increase inIR luminosity at bands L and M. The authors infer thathot dust from the inner regions is emitting the infraredexcess, while occluding the starlight and causing thedimming in the other bands.Observations by Antipin et al. (2015); Schneider et al.(2015) found that the absorption from optical to NIRwavelengths is gray, which could indicate the presenceof large particles causing the dimming ( (cid:38) µ m), andmeasured a dust column density of 2 × − g / cm ,although a similar absorption may be produced if anoptically thick disk of gas and small grains partiallycovers the star (Schneider et al. 2018).Then, the study of RW Aur A spectra by Facchini et al.(2016) also found that the inner accretion regions of thedisk are being occluded, and therefore the dimmings a r X i v : . [ a s t r o - ph . E P ] J a n G´arate et al. should come from perturbations at small radii.Finally, during the dimming in 2016, X-Ray observa-tions from G¨unther et al. (2018) indicate super-SolarFe abundances, along with a higher column density ofgas in the line of sight. The high N H /A V found by theauthors is interpreted as gas rich concentrations in theoccluding material, or as a sign of dust growth. Also,gray absorption was found again during this event.Given this information, several authors have discussedwhat mechanism would put the dust from the innermostregions of the protoplanetary disk, into the line of sight.Among the possible explanations are: a warped innerdisk (Facchini et al. 2016; Bozhinova et al. 2016), stellarwinds carrying the dust (Petrov et al. 2015; Shenavrinet al. 2015), planetesimal collision (G¨unther et al. 2018),and a puffed up inner disk rim (Facchini et al. 2016;G¨unther et al. 2018).Most of the proposed mechanisms rely on having enoughdust close to the star to cause the dimmings. So thefocus of this paper is to propose a mechanism that candeliver large amounts of dust to the inner regions of theprotoplanetary disk, by raising the dust accretion ratethrough the release of a dust trap.An increased dust concentration can explain someaspects of the dimmings, such as the high metallic-ity (G¨unther et al. 2018), the emission of hot dust(Shenavrin et al. 2015), and cause the dimmings pro-vided that another mechanism transports it into the lineof sight.1.2. A Fast Mechanism for Dust Accretion
A sudden rise in the dust accretion can occur in theearly stages of stellar evolution, when the stars are stillsurrounded by their protoplanetary disk composed ofgas and dust. The dynamics of the gas component ofthe disk is governed by the viscous evolution, whichdrives the accretion into the star (Lynden-Bell & Pringle1974), and the pressure support, that produces the sub-keplerian motion. On the other hand, the dust parti-cles are not affected by pressure forces, but suffer thedrag force from the gas. This interaction extracts an-gular momentum from the particles and causes themto drift towards the pressure maximum, with the driftrate depending on the coupling between the dust andgas motion (Whipple 1972; Weidenschilling 1977; Nak-agawa et al. 1986). This means that any bumps in thegas pressure act as concentration points for the dust. Inthese dust traps the grains accumulate, grow to largersizes, and reach high dust-to-gas ratios (Whipple 1972;Pinilla et al. 2012).One of the proposed mechanisms to generate a pressurebump is through a dead zone (Gammie 1996), a region with low turbulent viscosity (the main driver of accre-tion), due to a low ionization fraction which turns off themagneto-rotational instability (Balbus & Hawley 1991),that allows the gas to accumulate until a steady state isreached.The presence of a dead zone would allow the dust todrift and accumulate at its inner border (Kretke et al.2009), which can be located at the inner regions of theprotoplanetary disk ( r (cid:46) ust accretion preceding an outburst? MODEL DESCRIPTIONWe model a protoplanetary disk composed of gas anddust, solving the advection equation for both compo-nents. Our model includes multiple dust species andcoagulation as in Birnstiel et al. (2010), and the pres-ence of a dead zone in the inner disk. The advectionof gas and dust is solved using the mass conservationequation: ∂∂t ( r Σ g,d ) + ∂∂r ( r Σ g,d v g,d ) = 0 , (1)where r is the radial coordinate, Σ is the mass surfacedensity, v is the radial velocity, and the subindex ‘g’ and‘d’ refer to the gas and dust respectively.2.1. Viscous Evolution
The gas is accreted towards the star following the vis-cous evolution theory, for which the viscous velocity is: v ν = − g √ r ∂∂r ( ν Σ g √ r ) , (2)with ν the turbulent viscosity (Lynden-Bell & Pringle1974). If we neglect the effect of dust back-reaction wecan assume that the radial velocity of the gas is v g,r = v ν . We follow the Shakura & Sunyaev (1973) α modelfor the viscosity: ν = α c s Ω − K , (3)where the α parameter controls the intensity of the vis-cosity (and therefore the accretion), Ω K is the Keplerianangular velocity, and c s = (cid:112) γK b T /µm H is the adia-batic gas sound speed, with T the gas temperature, K b the Boltzmann constant, µ = 2 . γ = 1 . m H the hydro-gen mass.As a starting point for our model, we assume that thegas is in steady state, such that the accretion rate ˙ M g is constant over the disk radius following:3 π Σ g ν = ˙ M g . (4)2.2. Gas Orbital Velocity
The gas orbital motion is partly pressure supported,so for a negative pressure gradient it will orbit at sub-keplerian velocity. This point is particularly impor-tant for dust dynamics (Whipple 1972; Weidenschilling1977). The difference between the gas orbital velocityand the keplerian velocity, ∆ v g ,θ = v g ,θ − v K , is:∆ v g ,θ = − η v K , (5)with v K the keplerian speed and η = − (cid:0) ρ g r Ω K (cid:1) − ∂P∂r . (6) Here ρ g is the gas volume density in the disk’s midplane,and P = ρ g c s /γ is the pressure. In subsection 4.3, wewill consider the effects of the dust back-reaction on thegas orbital velocity.2.3. Dead Zone Model
Magneto-hydrodynamical models have predicted thepresence of a region with low turbulence at the innerregions of protoplanetary disks, commonly called “DeadZone”, caused when the ionization fraction is too lowfor the magneto-rotational instability (MRI) to operate(Gammie 1996).We parametrize the dead zone by using a variable α parameter over radius, while remaining agnostic to theunderlying physics, similar to previous research (Kretkeet al. 2009; Pinilla et al. 2016). Our profile is defined as: α ( r ) = α active − ∆ α · e rr − r < r α dead + ∆ α · e rr − r < r < r α dead + ∆ α · (1 − e − rr − ) r < r, (7)where α dead and α active are the characteristic valuesof the α parameter inside and outside the dead zone, r and r are its inner and outer edges, and ∆ α = α active − α dead . Although the dead zone shape is ratherarbitrary, it was chosen to have a smoother outer bor-der than the one of (Kretke et al. 2009), but retaininga sharp inner edge where the dust accumulates. Fig-ure 1 shows a diagram of the α profile, illustrating theshape of the different intervals and its main componentsto guide the reader with Eq. 7.In order to maintain the steady state from Eq. 4 we de-fine the initial surface density and temperature profilesas follow: Σ g ( r ) = Σ (cid:18) rr (cid:19) − α active α ( r ) , (8) T ( r ) = T (cid:18) rr (cid:19) − / , (9)where Σ and T are the values of the density and tem-perature at r = 1 au.2.4. Dead Zone Reactivation
While the dead zone is present the gas will remainin steady state and the dust will accumulate at its in-ner boundary. Yet, different processes can reactivatethe turbulence in the dead zone, allowing the accumu-lated material to flush towards the star. In the Gravo-Magneto instability (Martin & Lubow 2011) for exam-ple, the gas in the dead zone becomes gravitationally un-stable, raising the temperature to the point of triggering
G´arate et al. r r r deadactive Figure 1.
The diagram shows the shape of the α radialprofile for our dead zone model (in logarithmic scale). Thisconsists on three regions, the active inner zone limited by asharp decay at r , the dead zone with a smooth rise towardsits outer edge around r , and the outer active zone extendinguntil the outer boundary of the simulation. the MRI, and finally producing an accretion outburst.In our simulations we remain agnostic about the mech-anism that causes the reactivation, and only set the re-activation time t r arbitrarily, such that: α ( r, t > t r ) = α active . (10)2.5. Dust Dynamics
The dust dynamics are governed by the gas motionand the particle size. The time required for a particle ofsize a and material density ρ s to couple to the motionof the gas is called the stopping time, and is defined as: t stop = (cid:112) π ρ s ρ g ac s λ mfp /a ≥ / ρ s ρ g a ν mol λ mfp /a < / , (11)with the mean free path λ mfp = ( nσ H ) − (where n isthe number density, and σ H = 2 × − cm ), and themolecular viscosity ν mol = (cid:112) /πc s λ mfp (following thedefinitions in Birnstiel et al. 2010).A more useful quantity to describe the dust-gas couplingis the Stokes number (or dimensionless stopping time),which is defined as: St = t stop Ω K . (12)From this quantity we can quickly infer if a dust grainis coupled (St (cid:28)
1) or decoupled (St (cid:29)
1) to the gas.For the midplane this can be rewritten as:St = π aρ s Σ g λ mfp /a ≥ / π a ρ s λ mfp Σ g λ mfp /a < / . (13) The dust radial velocity is given in Nakagawa et al.(1986); Takeuchi & Lin (2002) as: v d = 11 + St v g ,r + 2 St St ∆ v g ,θ − D d Σ g Σ d ∂∂r ( Σ d Σ g ) . (14)Here, the first term of the dust velocity is responsiblefor small grains to move along with the gas, while thesecond term is responsible for the dust to drift towardsthe pressure maximum. The last term corresponds tothe dust diffusion contribution (see Birnstiel et al. 2010),with D d the dust diffusivity defined following Youdin &Lithwick (2007) as: D d = ν (1 + St ) . (15)The dust coagulation model is specified in Birnstiel et al.(2010). Since our model focus on the inner regions ofthe protoplanetary disks, the grain growth will alwaysbe limited by the fragmentation barrier (Brauer et al.2008), and the maximum grain size will be approxi-mately: St frag = 13 v αc s , (16)where v frag is the fragmentation velocity for dust parti-cles (Birnstiel et al. 2012). For silicates this correspondsto v frag ≈ / s (G¨uttler et al. 2010).2.6. Dust Back-reaction Effects
In protoplanetary disks, where the dust-to-gas ratio (cid:15) = ρ d /ρ g is assumed to be 0 .
01, the effect of the dustonto the gas is often neglected. However, in regions withhigher concentrations of solids, like in pressure bumps(Pinilla et al. 2012), the angular momentum transferedfrom the dust into the gas might be significant enoughto alter its dynamics.Many authors have already derived and included theeffect of back-reaction (coupled with viscous evolution)into numerical simulations, showing in which regimes itshould be considered (Tanaka et al. 2005; Garaud 2007;Kretke et al. 2009; Kanagawa et al. 2017; Taki et al.2016; Onishi & Sekiya 2017; Dipierro et al. 2018).In this paper we include the dust back-reaction in oneof our simulations to study its impact on our model. Todo so, the gas velocities v g ,r and ∆ v g ,θ are rewritten asfollows: v g ,r = Av ν + 2 Bηv K , (17)∆ v g ,θ = 12 Bv ν − Aηv K . (18)The back-reaction coefficients A and B measure the de-gree to which the gas dynamics are affected by the dust. ust accretion preceding an outburst? (cid:15) = 0 we obtain A = 1 and B = 0, recovering thetraditional velocities for the gas. When dust is present,the value of A decreases and B increases. Thus, fromEq. 17 we see that back-reaction slows down the viscousevolution with the term Av ν , and pushes the gas out-ward with the term 2 Bηv K .We find this contracted notation specially useful to sum-marize the back-reaction contribution to the gas dy-namic. SIMULATION SETUPIn this section we describe the observational constrainsrelevant for RW Aur A, the free parameters of ourmodel, and the setup of our 1D simulations using the twopoppy (Birnstiel et al. 2012) and
DustPy (Stamm-ler & Birnstiel, in prep.) codes.Our setup consists of three phases, the first phase sim-ulates the dust accumulation at the dead zone, using aglobal disk simulation over long timescales ( ∼ yr),but with a simplified and fast computational model forthe dust distribution using only two representative pop-ulations. As the first phase only tracks the evolutionof the surface density, in a second phase we recover thequasi-stationary particle size distribution at the innerdisk ( r ≤ Observational Constrains
RW Aur A is a young star with a stellar mass of M ∗ =1 . (cid:12) (Ghez et al. 1997; Woitas et al. 2001). The cir-cumstellar disk has an estimated mass around M disk ≈ × − M (cid:12) (Andrews & Williams 2005), presents ahigh accretion rate of ˙ M ≈ × − − × − M (cid:12) / yr DustPy is a new Python code that solves the diffusion-advection of gas and dust, and the coagulation-fragmentation ofdust, based on the Birnstiel et al. (2010) algorithm. (Hartigan et al. 1995; Ingleby et al. 2013; Facchini et al.2016), and extends from a distance of ∼ . T = 250 K, whichgives similar values to the Osterloh & Beckwith (1995)profile in the inner regions of the disk for our choice ofslope.Using these parameters, Eq. 3 and Eq. 4, we can con-straint the values for the density and turbulence. Fromthe disk accretion rate, mass and size we infer the valuefor the density Σ = 50 g / cm at r = 1 au, and the vis-cous turbulence α active = 0 .
1. These parameters yieldvalues of ˙ M g = 5 × − M (cid:12) / yr for the accretion rate,and M disk = 2 × − M (cid:12) for the disk mass (withoutconsidering the accumulation excess in the dead zone).The turbulence parameter α active used in our simula-tions is high, but necessary in order to account for thehigh accretion rates measured.3.2. Phase 1: Dust Concentration at the Dead Zone
In the first phase of our simulations we model the ac-cumulation of dust in a disk with a dead zone, to obtainthe dust-to-gas ratio radial profile.We use the TwoPopPy code to simulate a global pro-toplanetary disk with two representative populations ofthe dust species (details of the model can be found inBirnstiel et al. 2012). We initialize our simulations usingEq. 7, Eq. 8 and Eq. 9 for the α parameter, surface gasdensity and temperature profiles, with the values pro-vided by the observational constrains. For the dust-to-gas ratio we assume an uniform initial value of (cid:15) = 0 . r in =0 .
01 au to r out = 100 au, using n r = 500 radial gridcells with logarithmic spacing. In the fiducial model,the inner and outer boundaries of the dead zone are r = 0 .
51 au, r = 10 au, with a depth of α dead = 10 − .The simulation is evolved with this setup until the reac-tivation time t r = 10 yrs. Approximately at this pointthe dust reaches its maximum accumulation at the innerboundary of the dead zone, which will yield the maxi-mum dust accretion rate in the next phase. Since thegas is in steady state, we only evolve the dust in orderto minimize possible numerical errors. Inside the deadzone, the gas phase is (marginally) gravitationally sta-ble, with a Toomre parameter Q = c s Ω K / ( πG Σ g ) (cid:38) . G´arate et al. r (AU)10 ( g / c m ) Dust and Gas DensitiesGasDust ( t = 0 yrs)Dust ( t r = 10 yrs) Figure 2.
Gas and dust surface density obtained fromTwoPopPy, at the beginning and at the end of the dustconcentration phase. The gas (red) remains in steady stateduring this phase. The dust is initialized with a dust-to gasratio of (cid:15) = 0 .
01 (dashed blue line). The dust component isevolved for 10 yrs (solid blue line) in which the dust con-centrates at the inner disk, reaching (cid:15) ≈ .
24 at the innerboundary of the dead zone, and (cid:15) ≈ .
16 inside this region r < r = 0 .
51 au. dead zone inner edge reaching values of (cid:15) = 0 .
24, andconcentrating 110 M ⊕ between 0 .
51 - 0 . (cid:15) ≈ . Phase 2: Dust Size Distribution at the Inner Disk
In the second phase we want to recover the dust sizedistribution for multiple species, based on the dust-to-gas ratio and disk conditions obtained in the previoussection.We take the outcome of the TwoPopPy simulation asthe new initial conditions, and use the DustPy code tosolve the dust coagulation and fragmentation (followingthe study of Birnstiel et al. 2010) at the inner disk,while “freezing” the simulation exactly at the reactiva-tion time t = t r , while the dead zone is still present (i.e.still using Eq. 7).The mass grid consists of n m = 141 logarithmic-spacedcells, between m = 10 − − g, at every radius. Sincein this phase we only care about the inner disk, we adjustour simulation radial domain to be from r in = 0 .
05 auto r out = 5 au. The radial grid is defined as follow: •
25 linear-spaced grid cells at r = 0 . − .
09 au, •
120 logarithmic-spaced grid cells at r = 0 . − . •
20 logarithmic-spaced grid cells at r = 1 . − . S i z e ( c m ) Time: t r + 0.0 yr 5432101234 d u s t [ g / c m ] Figure 3.
Dust distribution in the inner region of the proto-planetary disk immediately before the dead zone reactivation( t = t r ). In the dead zone, where the turbulence and the col-lision speed of solids are lower, the dust particles can growto larger sizes ( a max ∼ a max ∼ µ m). The inner edge of the dead zone(marked by the white line) presents a high concentration oflarge dust grains. The innermost region is necessary to avoid numericalproblems with the inner boundary conditions. For opti-mization purposes we also turn off coagulation for r < .
09 au, since the growth and fragmentation timescalesare so short in this region that the simulation be-comes computationally unfeasible. Moreover, accord-ing to Akeson et al. (2005); Eisner et al. (2007) theinner boundary of RW Aur A disk should be around r ∼ . − . r = 0 . − . t = t r for everyradius. The dust distribution obtained at this phase isshown in Figure 3, where the grains adjust to the frag-mentation limit in the dead and active zones.3.4. Phase 3: Dead Zone Reactivation
For the final phase we simulate the evolution of dustand gas in the inner disk, after the reactivation of thedead zone ( t > t r ).Once again we use the DustPy code, this time to solvethe advection of gas and dust, along with the dustcoagulation-fragmentation. We start this phase from theconditions given at Section 3.3, using the same grid formass and radius, but now with the reactivated turbu- ust accretion preceding an outburst? Table 1.
Fiducial simulation parameters.Parameter ValueΣ
50 g / cm T
250 K r α active − α dead − r .
51 au r
10 au t r yrsBack-reaction Off Table 2.
Parameter variations.Simulation Parameter ChangedControl Simulation t r = 0 yrsShallow Dead zone α dead = 10 − Closer Inner Edge r = 0 .
25 auCloser Outer Edge r = 4 auBack-reaction On lence following Eq. 10. We let the simulation evolve for15 yrs, in which we expect that the material accumu-lated at the inner boundary of the dead zone will drifttowards the star. The results of this phase on the ac-cretion rate of gas and dust, as well as the final dustdistribution, will be shown in Section 4. Parameter Space
In Table 1 we summarize the parameters used for thedisk setup of our fiducial simulation. As the proper-ties of the dead zone are free parameters, chosen to bein a relevant range for the RW Aur dimming problem,we also require to explore (even briefly) the parameterspace for these properties, and see how they affect thefinal outcome of the simulations. We present five addi-tional simulations, changing one parameter of the fidu-cial model at a time, this way we explore the effect ofhaving: no initial dust accumulation at reactivation, dif-ferent dead zone properties, and the expected effects ofback-reaction in the final result. The parameter changesare described in Table 2. RESULTSIn this section we show the results obtained on thedust and gas dynamics after the “dead zone reactivation The simulation data files and a plotting script are available inzenodo: doi.org/10.5281/zenodo.1495061. ( g / c m ) DensityGasDust0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 d / g Dust-to-Gas RatioTime: t r + 15.0 yrTime: t r + 0 yr Figure 4.
The plots show the simulation state immediatelybefore the dead zone reactivation at t = t r (dashed lines), and15 yrs after it (solid lines). Top : Evolution of the gas (red)and dust (blue) surface densities. The initial state showsthe gas steady state profile and the accumulation of dust atthe inner boundary of the dead zone. After reactivation theaccumulation of dust flushes towards the star faster than thegas.
Bottom : Dust-to-gas ratio evolution. At the initial statethe inner region presents an already high solid concentrationthanks to mixing at the dead zone boundary. During theflushing the dust-to-gas ratio reaches values of (cid:15) = 5 at someof the regions where the dust concentration arrived beforethe gas. phase” (Section 3.4), including the final dust distribu-tion, the dust-to-gas ratio at the inner boundary, andaccretion rates of dust and gas. We also study the im-pact of the dead zone parameters on the final outcome,to see if these results follow the expected behavior.The gas and dust surface densities before and after thereactivation are shown in Figure 4. The initial sur-face density obtained from the first accumulation phaseshows a dust-to-gas ratio of (cid:15) = 0 .
24 at the dead zoneinner edge, and (cid:15) = 0 .
16 in the inner disk ( r (cid:46) . G´arate et al. S i z e ( c m ) Time: t r + 0.05 yr 5432101234 d u s t [ g / c m ] S i z e ( c m ) Time: t r + 15.0 yr 5432101234 d u s t [ g / c m ] Figure 5.
Dust distribution in the inner region of the proto-planetary disk after 0 .
05 and 15 yrs of the dead zone reactiva-tion.
Top : The dust that was accumulated at the dead zonediffuses to the inner region within ∼
10 collisional times, gen-erating high dust-to-gas concentrations. The original edge ofthe dead zone is marked in white.
Bottom : Afterwards, thedust drifts towards the inner disk regions ( r ∼ . − . ∼
15 yrs, while adjusting to the new fragmentationlimit.
After the reactivation ( t > t r ) the dust accumulated atthe dead zone is transported towards the inner regionsfaster than the gas, reaching the inner boundary of thedisk ( r in ∼ . − . (cid:15) > < − ), the only impactof dust back-reaction is to slow down the dust and gasevolution. Therefore no instabilities are generated and M ( M / y r ) Accretion rate at r: 0.15 AUGasDust
Figure 6.
Accretion rate of gas (red) and dust (blue) at r = 0 .
15 au for 15 yrs after the dead zone reactivation. Theinitial gas accretion rate is ˙ M g = 5 × − M (cid:12) / yr, in agree-ment with the observations of RW Aur A, while the dustaccretion rate is ˙ M d = 7 × − M (cid:12) / yr. Upon reactivationthe dust accretion rate increases faster than the gas, eventu-ally surpassing it after 10 yrs, and reaching a high value of˙ M d = 6 × − M (cid:12) / yr. we can proceed with our analysis.Upon entering the active zone the particles fragmentdue to the high turbulence and adjust to their fragmen-tation limit (see Eq. 16) in a few collisional timescales t coll , which can be approximated by: t coll = ( n d σ ∆ v turb ) − , (19)where n d is the number density of dust particles, σ ≈ πa is the collisional cross section, and ∆ v turb ≈√ α St c s is the turbulent collision speed (Ormel & Cuzzi2007).We find that during 10 collisional timescales ( t ∼ .
05 yrs) after the dead zone reactivation the dust grainsdiffuse inward faster than the gas, gaining a head startthat leads to high dust-to-gas ratio concentrations. Weattribute this feature to the sudden rise in the turbu-lence α at the dead zone inner edge, that increases thedust diffusivity and spreads the particles towards theinner regions (see Eq. 15 and Figure 5). After the dusthas adjusted to the new fragmentation limit, it driftsroughly with the viscous velocity of the gas v ν towardsthe inner boundary of the disk. The particles reach-ing the inner boundary reach maximum sizes between a max = 10 − µ m (see Figure 5).The accretion rate (measured at r = 0 .
15 au) of bothdust and gas increases after the reactivation of the deadzone, as the accumulated material arrives at the innerboundary of the disk (see Figure 6). Before the deadzone reactivation, the gas accretion rate is given by ust accretion preceding an outburst? M g = 5 × − M (cid:12) / yr,similar to the observational value of (Facchini et al.2016; Ingleby et al. 2013), and the dust accretion rate is˙ M d = 7 × − M (cid:12) / yr, this value comes from the dustdiffusing into the inner disk during the concentrationphase.After the dead zone reactivation the dust concentra-tion moves inwards, and the accretion rate at the innerboundary of the disk becomes dominated by the dust, tothe point of surpassing that of the gas. This high sup-ply of solid material, with ˙ M d = 6 × − M (cid:12) / yr couldcause hot dust and metallicity features of RW Aur A(Shenavrin et al. 2015; G¨unther et al. 2018), and providean ideal environment for the dimmings to occur (see Sec-tion 5.1). At this point we also note that the accretionrate of gas has increased up to ˙ M g = 10 − M (cid:12) / yr.In Section 5 we will discuss how the high accretion ofsolids could cause the dimmings in the context of previ-ous proposed mechanisms (dusty winds, puffed-up innerdisk rim, etc), and if we can expect future accretionsignatures from the gas. In the following subsections,we study the effect of the simulation parameters on thedust dynamics.4.1. Simulation without Dust Concentration
In our model we remained agnostic to the reactivationprocess of the dead zone, and allowed the dust to ac-cumulate for long enough time to reach concentrationsas high as (cid:15) = 0 .
24 at its inner edge. Depending on themechanism that reactivates turbulence, the flushing ofsolid material towards the star may occur earlier withlower dust concentrations, reducing the total accretionof solids. To model the limit case in which no dustconcentration occurs, we repeat our setup with a con-trol simulation, but now setting the reactivation timeto t r = 0 yrs.In Figure 7 we show the initial and final state of bothsimulations. Here the control simulation has an uni-form (cid:15) = 0 .
01 in the beginning, since no additional dustconcentration has occurred. Notice that we still assumethat the gas has reached the steady state density pro-file.After the reactivation the gas and dust excess at thedead zone drift towards the star. As in the fiducial case,the dust that was located at the dead zone inner edgearrives at the inner boundary of the disk before thegas, also in a time of ∼
15 yrs. The only difference isthat now the material being accreted has a dust-to-gasratio of (cid:15) = 0 .
28, which is still higher than the initial (cid:15) = 0 .
01, although not as extreme as the (cid:15) = 5 foundin the fiducial case. Here the dust accretion rate at theinner boundary of the disk ( r ≈ .
15 au) can reach up ( g / c m ) Control - Density (Time: t r )Gas (Fiducial)Dust (Fiducial)Gas (Control)Dust (Control)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 ( g / c m ) Density (Time: t r + 15.0 yr)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 d / g Dust-to-Gas Ratio (Time: t r + 15.0 yr)FiducialControl Figure 7.
Comparison between the “fiducial” simulation( t r = 10 yrs, solid lines) and the “control” simulation ( t r =0 yrs, dashed lines). Top : Gas and dust densities at thereactivation time ( t = t r ). In both cases the gas is in steadystate, but in the control simulation the dust did not had thetime to accumulate at the inner edge of the dead zone, inthis case the reactivation occurs with a uniform dust-to-gasratio (cid:15) = 0 . Mid : Gas and dust densities after the deadzone reactivation ( t = t r + 15 yrs). The dust that was at thedead zone in steady state still arrives faster than the gas tothe inner boundary of the disk. Bottom : Dust-to-gas ratioafter the reactivation. The maximum dust concentration inthe “control” simulation is now (cid:15) = 0 . G´arate et al. to ˙ M d = 3 × − M (cid:12) / yr.From here we learn that the dust arrival time at the in-ner boundary does not depend on the amount of solidsaccumulated at the dead zone inner edge, and that uponreactivation the accreted material will still carry a highconcentration of solids.4.2. Simulations for Different Dead Zone Properties
The dead zone shape is parametrized by its edges r and r , and the turbulence parameter α dead followingEq. 7, and altering these parameters also changes arrivaltime and dust-to-gas ratio of the accreted material.By shifting the inner edge of the dead zone to smallerradii ( r = 0 .
51 au → .
25 au) the dust concentrationduring the first phase will also be located closer to theinner boundary of the disk (see Figure 8). Now it onlytakes the dust between 3 − r = 10 au → (cid:15) = 2 . α dead of the dead zone. Our fiducial simulation consid-ered an α dead = 10 − , which in contrast with the ac-tive zone α active = 10 − leads to an accumulation ofmaterial in the dead zone with a factor of 1000 rela-tive to the steady state of a fully active disk, this ofcourse favors the accretion of massive amounts of gasand solids upon reactivation. By taking a shallower deadzone ( α dead = 10 − → − ) there is less gas and dustaccumulated, so upon reactivation the flushing of mate-rial is slower by a factor of a few (see Figure 10). Wealso find that for α dead = 10 − there is no significantconcentration at the dead zone inner edge, after 10 yrsthe dust-to-gas ratio rises only up to (cid:15) = 0 . (cid:15) = 0 .
06. Therefore we find that the turbulence param-eter α dead of the dead zone is the main determinant ofthe total amount of material being accreted, and thatdeeper dead zones are necessary to produce the dustdominated accretion seen in our fiducial simulation. ( g / c m ) CloseIn - Density (Time: t r )Gas (Fiducial)Dust (Fiducial)Gas (CloseIn)Dust (CloseIn)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 ( g / c m ) Density (Time: t r + 5.0 yr)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 d / g Dust-to-Gas Ratio (Time: t r + 5.0 yr)FiducialCloseIn Figure 8.
Same as Figure 7, but comparing the “fiducial”simulation ( r = 0 .
51 au) with the “closer inner edge” simu-lation ( r = 0 .
25 au).
Top : Initially, the dead zone is moreextended toward the inner boundary of the disk.
Mid : Afterreactivation, the dust concentration arrives in only 5 yrs tothe inner edge of the disk, also moving faster than the gas.
Bottom : The maximum dust-to-gas ratio for the “closer in-ner edge” is similar to the “fiducial” value with (cid:15) ≈ ust accretion preceding an outburst? ( g / c m ) CloseOut - Density (Time: t r )Gas (Fiducial)Dust (Fiducial)Gas (CloseOut)Dust (CloseOut)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 ( g / c m ) Density (Time: t r + 15.0 yr)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 d / g Dust-to-Gas Ratio (Time: t r + 15.0 yr)FiducialCloseOut Figure 9.
Same as Figure 7, but comparing the “fiducial”simulation ( r = 10 au) and the “closer outer edge” simula-tion ( r = 4 au). Top : Initially the dead zone is smaller andhas less material, leading to a lower dust concentration atits inner edge of (cid:15) = 0 . Mid : The dust excess once againarrives to the inner boundary of the disk before the gas in15 yrs, but in a lower concentration.
Bottom : The maxi-mum dust-to-gas ratio for the “closer outer edge” simulationis now (cid:15) ≈ ( g / c m ) Shallow - Density (Time: t r )Gas (Fiducial)Dust (Fiducial)Gas (Shallow)Dust (Shallow)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 ( g / c m ) Density (Time: t r + 15.0 yr)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 d / g Dust-to-Gas Ratio (Time: t r + 15.0 yr)FiducialShallow Figure 10.
Same as Figure 7, but comparing the “fiducial”simulation ( α dead = 10 − ) and the “shallower dead zone”simulation ( α dead = 10 − ). Top : There is no notable accu-mulation of material at the inner edge of the dead zone dur-ing the accumulation phase, with (cid:15) = 0 .
012 at most. Prob-ably because the particles are too small to get trapped, anddiffuse more easily to the inner region.
Mid : After 15 yrsthere is only a little dust excess traveling towards the innerregion, a small bump can still be appreciated in the dustsurface density profile.
Bottom : The final dust-to-gas ratioin the inner boundary is only (cid:15) = 0 .
06, just a factor of a fewabove the original (cid:15) = 0 .
01. This is because the lack of morematerial in the entire dead zone. G´arate et al.
Simulation with Dust Back-reaction
In all our results until this point we have neglectedthe back-reaction of dust to the gas, however for thedust-to-gas ratios presented during the dead zone reac-tivation ( (cid:15) (cid:38)
1) this effect should be relevant. In thissection we study its impact to see if our previous resultsremain valid.First, we should mention that in this setup we only con-sider the back-reaction in the “reactivation phase”, andassume that the dust will still accumulate at the inneredge of the dead zone even if back-reactions are con-sidered. This is still justified since our particles aretoo small (St < − ) to cause any perturbation be-yond slowing down the concentration process and wecan be infer it also by studying the single dust speciesscenario described in the Appendix A.1. Studies of On-ishi & Sekiya (2017) also showed that the back-reactionstill allows the dust to accumulate at pressure maxima,and that the dust traps do not self-destruct by this ef-fect when taking into account the vertical distributionof solids.For the reactivation phase we implement the gas ve-locities as described by equations Eq. 17 and Eq. 18.In the radial direction, the gas velocity now consist oftwo terms modulated by the back-reaction coefficients0 < A, B <
1, the term Av ν is slowing down the viscousevolution of the gas respect to the default value v ν , andthe term 2 Bηv K is pushing the gas in the direction op-posite to the pressure gradient.Since the vertical distribution of particles is not exactlythe same as the gas, the effect of the back-reaction is alsonot uniform in the vertical direction, the importance ofthis point is shown in Dipierro et al. (2018); Onishi &Sekiya (2017). To account for the vertical effect of back-reactions in our 1D simulations we take the verticallyaveraged velocity for the dust and gas, weighted by themass density to conserve the total flux. The details forthis implementation can also be found in the AppendixA.2.In Figure 11 we show a comparison between the “fidu-cial” and “back-reaction” simulations 15 years after thedead zone reactivation. When the back-reaction is con-sidered the most notable effect is the slowing down ofthe accretion of material by a factor of ∼
2. While thedust particles in the “fiducial” simulation take 15 yearsto reach the inner boundary of the disk, in the “back-reaction” simulation they need 30 years instead.The reason why no further effects are observed is thatthe high turbulence causes fragmentation of the parti-cles to smaller sizes (as seen in Figure 5), where theyare unable to “push” the gas backwards (i.e. the term B → ( g / c m ) Backreaction - Density (Time: t r )Gas (Fiducial)Dust (Fiducial)Gas (Backreaction)Dust (Backreaction)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 ( g / c m ) Density (Time: t r + 15.0 yr)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 d / g Dust-to-Gas Ratio (Time: t r + 15.0 yr)FiducialBackreaction Figure 11.
Same as Figure 7, but comparing the “fiducial”simulation and the “back-reaction” simulation.
Top : Bothsimulations start at the reactivation time with the same ini-tial conditions.
Mid : When back-reaction is considered, theevolution of the dust and gas component is slower than in thefiducial case. This is because the high dust concentrationsslow down the viscous evolution of the gas, which in turnalso slows down the drifting of the dust towards the innerdisk.
Bottom : Both the simulation with and without back-reaction present a high concentration of dust in the accretedmaterial, yet in the case with back-reaction the bulk of dustreaches a radii of only r = 0 .
22 au in 15 years, while the dustin the fiducial simulation is already at r = 0 .
15 au. ust accretion preceding an outburst? A ≤ (cid:15) = 2 . (cid:15) = 5 value). Thismeans that in all our previous results we should considerthat the accreted material will take more time (a factor ∼
2, if the dust-to-gas ratio is high enough) to reach theinner boundary of the protoplanetary disk, and as con-sequence that the accretion rate will also be reduced.An order of magnitude estimate for the back-reactioncoefficients can be obtained by approximating the parti-cle distribution to a single size population (see AppendixA.1), in which A ≈ ( (cid:15) +1) − and B ≈ St (cid:15) ( (cid:15) +1) − , how-ever to have an overall estimate for the entire disk thedust-to-gas distribution should be taken into account. DISCUSSIONWe have seen that through the reactivation of a deadzone, located in the inner regions of RW Aur A cir-cumstellar disk, large amounts of dust can be flushedtowards the star in timescales that can go between5 −
30 years. In this section we compare our resultson dust accretion with the properties of the dimmings,speculate upon the future accretion signatures of RWAur, and discuss which ingredients of our model may beimproved to better match the observations.5.1.
The Fast Accretion Mechanism in the Context ofRW Aur A Dimmings
Multiple observations during the dimmings of RW AurA reveal the presence of large amounts of dust at smallradii, in the line of sight, and in the material accretedby the star (see Section 1.1). Our model with dead zonereactivation provides a way to increase the dust at theinner rim of the protoplanetary disk by 2 orders of mag-nitude from its original value (from Σ d = 50 g / cm to5 × g / cm at 0 .
15 au). Although this alone does notexplain the dimmings, as a significant amount of mate-rial still needs to be moved towards the line of sight, itrelaxes the conditions of other mechanisms that can doso.In the case that the dimmings are caused by a dustywind coming from the inner regions (Petrov et al. 2015;Shenavrin et al. 2015), the dead zone reactivation in-creases the amount of dust entangled with the gas, andthen the particles would be dragged into the line of sightby the wind if they are small enough to be coupled tothe gas.Similarly, if the dimmings are caused by a puffed up in-ner rim of the disk (Facchini et al. 2016; G¨unther et al. 2018), the reactivation of the dead zone not only in-creases the amount of solid material in the line of sight(since our particles are small they should be well mixedwith the gas in the vertical direction), but also can in-crease the scale height of these regions with the rise intemperature due to accretion heating. If our model iscorrect, in the following decade(s) the gas accretion rateshould also increase by more than one order of magni-tude, up to ˙ M g = 10 − M (cid:12) / yr, as the gas excess fromthe dead zones adjusts to a new steady state.5.2. A Single Reactivation Event, or Multiple ShortReactivation Spikes?
The dimmings of RW Aur A last from a few monthsup to two years (see the list in Rodriguez et al. 2018),repeatedly moving from the dimmed state to the brightstate. The accretion process described in our model lastfrom several years to a few decades (depending on thesimulation parameters), which is clearly longer than thetypical dimming duration, and without presenting anydecrease in the dust accretion rate over the process. Yet,we can think of a few ways to reconcile the timescalesof our simulations with the dimmings.The first possibility is that the reactivation of the deadzone, and the subsequent increment in the accretion rate(see Figure 6), can also raise the temperature of the in-ner disk through viscous heating, resulting in a puffed-up inner rim. In this case, the dimmings will occur inthe time required to drag the dusty material into the lineof sight as the scale height of the inner disk increases,and will finish once the material settles back down or iscompletely accreted.In a similar way, if the dead zone reactivation is accom-panied by stronger stellar winds, these will determinethe timescale of the dimmings within the accretion eventdescribed by our model.Also, the accretion process described in our model doesnot need to increase smoothly, and can present vari-ability in shorter timescales. Instabilities that are notresolved by our model during the accretion process mayproduce a bumpy surface density profile in dust andgas, for example through the ring instability (W¨unschet al. 2005) for dead zones in layered disks. In this casethe dimming events would correspond to the local max-imums in the accretion of dust.Another possibility is that the dead zone reactivation isnot instantaneous. In the case that only inner edge be-comes active, a fraction of the accumulated dust willstart drifting towards the star, leaving most of thegaseous and solid material still trapped in the dead zone.If multiple reactivation events like this take place, thecorresponding dust excess will also arrive in intervals,4
G´arate et al. with a frequency depending on the reactivation mech-anism. In this case the dimmings would start when aspike of dust accretion arrives to the inner edge, andend once it decreases back to its steady state value. Toimprove our model we would need to resolve the ther-mal and gravitational instabilities that can reactivatethe dead zone (Martin & Lubow 2011).Finally, azimuthal asymmetries in the inner disk (suchas vortices) may also add variability to the accretionprocess, but these are not considered by our model.A further monitoring of the metallicity of the accretedmaterial would help understand the nature of the dustaccretion process and the mechanism that drags thedust into the line of sight. If the dust accretion contin-ues delivering material to the star, independently of thedimmings, then the metallicity should remain high evenafter the luminosity returns to the bright state. Oth-erwise, if the rise in metallicity keeps correlating withthe dimmings (as in G¨unther et al. 2018), the dimmingscould be an outcome of the sudden increase in the dustaccretion (although of course, correlation does not implycausality).5.3.
Validity of the Dead Zone Model
Our results showed that the exact values of dust ac-cretion rate and timescales depend sensitively on theparameters used for the dead zone. A dead zone inneredge closer to the inner boundary of the disk reducesthe timescale of the process, a closer outer edge reducesthe total mass of the dead zone, and the turbulence pa-rameter α dead and the reactivation time determine theamount of dust that can be trapped and flushed towardsthe star. In addition, the shape of the dead zone profilealso affects dust and gas surface density profiles, and aproper modeling of the gas turbulence would be requiredto obtain their final distributions after the dead zone re-activation.With this amount of free parameters, our simulationsprovide more a qualitative scenario than quantitativepredictions. Further constrains are necessary to deter-mine how relevant the reactivation of a dead zone canbe for RW Aur A dimmings. The first step is of courseto find the mechanism that puts the required amountsof dust in the line of sight, and obtain an estimate ofthe required dust surface density at the inner disk (notonly in the line of sight) to produce this phenomena.In parallel, any signature of enhanced gas accretion ratein the following years would speak in favor of the deadzone reactivation mechanism as one of the drivers of thedimmings.Additionally, constrains on the properties of the innerdisk would limit the parameter space described. In par- ticular, the mass in the inner 10 au of RW Aur A diskwould be useful to constrain the outer edge and turbu-lence parameter of the dead zone, and with them theamount of material that can be thrown to the star.5.3.1. The Dead Zone as an Accretion Reservoir
One point that speaks in favor of our model is its po-tential to sustain the large accretion rates of RW Aur forextended periods of time. Considering only the observedvalues for the disk mass and the accretion rate, the max-imum lifetime of RW Aur would be of M g / ˙ M g ∼ -10 yrs which is too short for a T Tauri star.Rosotti et al. (2017) defined the dimensionless accretionparameter: η acc = τ ∗ ˙ M /M disk , (20)with τ ∗ the age of the star, that indicates if the proper-ties of a disk are consistent the steady state accretion, inwhich case it follows η acc (cid:46)
1. The RW Aur A is around5 Myr old (Ghez et al. 1997) and presents an accretionparameter of η acc ≈
60, indicating either that the diskis not in steady state, or that the observed disk mass isunderestimated.In our model, the dead zone provides a reservoir of ma-terial able to sustain the detected accretion rates foraround 2Myrs (considering the parameters of our fidu-cial simulation), which is close to the estimated ageof the star. At the same time, at high densities thedusty material would be optically thick, remaining hid-den from the mm observations used to measure the diskmass.5.3.2.
Another Free Parameter for Dust Growth?
In our model we considered that the turbulent α that limits dust growth by fragmentation (Eq. 16), isthe same that drives the viscous evolution of the gas(Eq. 3), yet recent models explore the case with two in-dependent α values for each process (e.g., Carrera et al.2017). For our model, using a single α value means thatparticles reach bigger sizes while they remain in thedead zone, and fragment to smaller sizes in the activeregion.Allowing two independent alpha values would allow theformation of large particles in the active region. Theselarger particles drift faster, but also exert a strongerback-reaction on the gas. The “pushing” back-reactioncoefficient B is roughly proportional to the particle size(see Appendix A.1), and at the high dust-to-gas ratiosfound during the dead zone reactivation, it could bestrong enough to generate density bumps in the gas, oreven trigger the streaming instability for particles withlarge enough Stokes number (Youdin & Goodman 2005). ust accretion preceding an outburst? SUMMARYIn this work we studied a new mechanism that canincrease the concentration of solids in the inner regionsof a protoplanetary disk in timescales of ∼
10 years,through the reactivation of a dead zone.This study was motivated by the recent dimmings ofRW Aur A, which present a high concentration dustin the line of sight (Antipin et al. 2015; Schneideret al. 2015), and an increased emission from hot grainscoming from the inner regions of the protoplanetarydisk(Shenavrin et al. 2015), and subsequently observedsuper-solar metallicity of the accreted material(G¨untheret al. 2018).Using 1D simulations to model the circumstellar disk ofRW Aur A, we find that the dust grains accumulate atthe inner edge of the dead zone, which acts as a dusttrap, reaching concentrations of (cid:15) ≈ .
25. When the tur-bulence in this region is reactivated, the excess of gasand dust is released from the dead zone and advectedtowards the star. By effect of dust diffusion and gasdrag, the dust component can reach the inner boundaryof the protoplanetary disk before the gas component,producing high dust concentrations of (cid:15) ≈ M d =7 × − M (cid:12) / yr to 6 × − M (cid:12) / yr in only 15 years.This scenario can provide an ideal environment for othermechanisms, such as stellar winds (Petrov et al. 2015; Shenavrin et al. 2015) or a puffed up inner rim (e.g.,Facchini et al. 2016; G¨unther et al. 2018), to transportthe required amount of solid material into the line ofsight and cause the dimmings, although further studiesare required to link the surface density at the midplanewith the measured dust concentrations.Additionally, our simulations predict that in the follow-ing decade(s) the gas accretion rate should also rise byan order of magnitude, from ˙ M g = 5 × − M (cid:12) / yr to10 − M (cid:12) / yr if the dead zone reactivation is the mecha-nism transporting dust towards the disk inner region.ACKNOWLEDGMENTSWe would like to thank S. Facchini, and also theanonymous referee, for their useful comments that im-proved the extent of this work. M. G., T. B., andS. M. S. acknowledge funding from the European Re-search Council (ERC) under the European Unions Hori-zon 2020 research and innovation programme undergrant agreement No 714769 and funding by the DeutscheForschungsgemeinschaft (DFG, German Research Foun-dation) Ref no. FOR 2634/1. HMG was supportedby the National Aeronautics and Space Administrationthrough Chandra Award Number GO6-17021X issuedby the Chandra X-ray Observatory Center, which is op-erated by the Smithsonian Astrophysical Observatoryfor and on behalf of the National Aeronautics Space Ad-ministration under contract NAS8-03060.APPENDIX A. DETAILS OF DUST BACK-REACTIONThe gas and dust velocities incorporating the back-reaction of the dust into the gas is obtained from the momentumconservation equations shown in Nakagawa et al. (1986), and considering the force exerted by the multiple species ofdust as in (Tanaka et al. 2005). In this case the gas feels the drag force of multiple species, the pressure force, theviscous force, and the stellar gravity, while the dust only feels the stellar gravity and the drag force from the gas.The momentum equations for the gas and dust then are respectively:d v g d t = − (cid:90) ( v g − v d ) t stop ρ d ( m ) ρ g d m − GM ∗ r ˆr − ρ g ∂P∂r ˆ r + f ν ˆ θ, (A1)d v d d t = − ( v d − v g ) t stop − GM ∗ r ˆr , (A2)where the viscous force contribution in the azimuthal direction can be conveniently written as f ν = Ω K v ν / A , B are defined as: A = X + 1 Y + ( X + 1) , (A3) B = YY + ( X + 1) , (A4)6 G´arate et al. with X and Y following the notation of Okuzumi et al. (2012): X = (cid:90)
11 + St ρ d ( m ) ρ g d m, (A5) Y = (cid:90) St1 + St ρ d ( m ) ρ g d m. (A6)The X, Y integrals come from the contribution of the multiple dust species at the momentum conservation equationsdescribed by Tanaka et al. (2005). Here ρ d ( m ) is the dust volume density per mass, and the dust-to-gas ratio at eachmass bin (cid:15) ( m ) = ( ρ d ( m ) /ρ g ) dm determines the contribution of each particle species to the final result.To summarize the effect of back-reactions on the gas, we can understand the coefficient A as a “smoothing” factorand the coefficient B as a “pushing” factor. The “smoothing” factor A reduces the effect of the viscous force in theradial velocity (slowing down the accretion), and reduces the effect of the pressure gradient in the azimuthal velocity(making the gas orbital velocity more keplerian). The “pushing” factor B on the other hand, pushes the gas in thedirection opposite to the pressure gradient in the radial direction (with the term 2 Bηv K ), and slows down the orbitalvelocity (with the term Bv ν / Single Species Analysis
The expressions for A and B are difficult to study by eye, but we can consider the case with a single species ofparticles to simplify subsequent analysis. In this case the integrals X and Y become: X single = 11 + St (cid:15), (A7) Y single = St1 + St (cid:15). (A8)The back-reaction coefficients now have a simple expression that only depends on the Stokes number and the dust-to-gas ratio: A single = (cid:15) + 1 + St ( (cid:15) + 1) + St , (A9) B single = (cid:15) St( (cid:15) + 1) + St . (A10)From these equations we can start noticing some interesting values for A and B : • < A, B < • The limit without particles (cid:15) →
0, recovers the traditional gas velocities: – A → – B → • The limit with small particles, assuming St << – A ≈ ( (cid:15) + 1) − – B ≈ St (cid:15) ( (cid:15) + 1) − The last case is applicable to our simulations, since we have high dust-to-gas ratios of small particles only, in this casewe reach the limit of B ≈ A ≈ (1 + (cid:15) ) − , and therefore we only care about the slowing effect of the particles onthe gas. ust accretion preceding an outburst? Considering the Vertical Structure
The dust-to-gas ratio is not necessarily constant in the vertical direction. As shown by Okuzumi et al. (2012);Dipierro et al. (2018) this might have an impact in the effects of the back-reaction, since the velocities of the gas anddust will experience a different force in the different layers of the disk.We assume that both the gas and dust densities are distributed as a Gaussian in the vertical direction: ρ g ( z ) = Σ g √ πH g exp( − z H ) , (A11) ρ d ( m, z )d m = Σ d ( m ) √ πH d ( m ) exp( − z H ( m ) ) . (A12)Here ρ d ( m, z ) dm is the volume density of the particles in the mass bin m at a height z , and Σ d ( m ) is the surfacedensity of particles at the mass bin m . H g and H d ( m ) are the scale heights of the gas and the dust particles withmass m , respectively defined as: H g = c s Ω K , (A13) H d ( m ) = H g · min(1 , (cid:114) α min(St , / ) ) , (A14)the last equation coming from Birnstiel et al. (2010), that tells us that bigger particles are more concentrated towardsthe midplane, while small particles are distributed like the gas.With the vertical distribution of gas and dust we can obtain the mass weighted average radial velocity ¯ v g,d , defined toconserve the radial mass flux, as: Σ g,d ¯ v g,d = (cid:90) + ∞−∞ ρ g,d ( z ) v g,d ( z )d z. (A15)Now the densities, dust-to-gas ratio and Stokes number of the particles can be defined in the vertical direction, allowingus to have the vertical distribution of back-reaction coefficients and radial velocities.In terms of implementation, we perform the integral Eq. A15 over n z = 50 logarithmically spaced grid cells locallydefined between 10 − − H g . Our only assumption to simplify the numeric calculation is that the viscous velocity v ν and the sub-keplerian velocity ηv K are constant over the vertical direction. B. HIGHER RESOLUTION TESTTo validate our results we perform an additional simulation using the same parameters as in the fiducial set-up, butwith an increased radial resolution for the different phases.In the Phase 1 of dust concentration (Section 3.2) we increase to n r = 1000 radial grid cells between r = 0 . −
100 au.In the Phase 2 and 3 of build up and dead zone reactivation (Section 3.3, 3.4) we use: •
40 linear-spaced grid cells at r = 0 . − .
09 au, •
200 logarithmic-spaced grid cells at r = 0 . − . •
30 logarithmic-spaced grid cells at r = 1 . − . n m = 141 logarithmically spaced grid cells.In Figure 12 we see that the evolution of gas remains the same independent of the resolution used. The final dustsurface density profile and the dust-to-gas ratio are slightly lower for the higher resolution simulation over the entireinner region, the maximum concentration is reduced from (cid:15) = 5 to (cid:15) = 4.This decrease in the dust surface density apparently happens because, during the dust concentration phase, thematerial accumulated at the dead zone inner edge diffuses more efficiently towards the inner disk, slightly loweringthe dust-to-gas ratio in the dust trap. In the high resolution case, the accumulated dust mass between 0 . . ⊕ (in contrast with the 110 M ⊕ accumulated in the low resolution case).Other than the slight redistribution of dust, the high resolution simulation behaves in the same way as the standardcase, and our conclusions are maintained.We tested that the dust distribution obtained during the dust concentration phase (Phase 1) converges for resolutionsresolutions higher than n r = 1000 (tested up to n r = 3000), and therefore we do not expect any further changes inthe subsequent results after the dead zone reactivation (other than those already reported in this section).8 G´arate et al. ( g / c m ) High-Res - Density (Time: t r )Gas (Fiducial)Dust (Fiducial)Gas (High-Res)Dust (High-Res)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 ( g / c m ) Density (Time: t r + 15.0 yr)0.1 0.2 0.3 0.4 0.5 0.6 0.7r (AU)10 d / g Dust-to-Gas Ratio (Time: t r + 15.0 yr)FiducialHigh-Res Figure 12.
Comparison between the “fiducial” (solid line) simulation and the “high resolution” (dashed line) simulation.
Top : At higher resolutions, there is less dust accumulated at the dead zone inner edge, due to more effective diffusion of solidstowards the inner regions in the dust concentration phase.
Mid : After 15 yrs of evolution, the gas profile remains practicallyidentical with the increased resolution, the dust profile is slightly lower at all radii in the high resolution case.
Bottom : Thefinal dust-to-gas ratio in the high resolution case is slightly reduced, the maximum dust concentration now has a value of (cid:15) ≈ ust accretion preceding an outburst?19REFERENCES