MRI-active inner regions of protoplanetary discs. I. A detailed model of disc structure
Marija R. Jankovic, James E. Owen, Subhanjoy Mohanty, Jonathan C. Tan
MMNRAS , 1–21 (2020) Preprint 26 February 2021 Compiled using MNRAS L A TEX style file v3.0
MRI-active inner regions of protoplanetary discs. I. A detailedmodel of disc structure
Marija R. Jankovic, ‹ James E. Owen, Subhanjoy Mohanty and Jonathan C. Tan , Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK Dept. of Astronomy, University of Virginia, Charlottesville, Virginia 22904, USA Dept. of Space, Earth & Environment, Chalmers University of Technology, Gothenburg, Sweden
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Short-period super-Earth-sized planets are common. Explaining how they form near theirpresent orbits requires understanding the structure of the inner regions of protoplanetarydiscs. Previous studies have argued that the hot inner protoplanetary disc is unstable to themagneto-rotational instability (MRI) due to thermal ionization of potassium, and that a localgas pressure maximum forms at the outer edge of this MRI-active zone. Here we present asteady-state model for inner discs accreting viscously, primarily due to the MRI. The structureand MRI-viscosity of the inner disc are fully coupled in our model; moreover, we accountfor many processes omitted in previous such models, including disc heating by both accretionand stellar irradiation, vertical energy transport, realistic dust opacities, dust effects on discionization and non-thermal sources of ionization. For a disc around a solar-mass star with astandard gas accretion rate ( 𝑀 „ ´ M d yr ´ ) and small dust grains, we find that the innerdisc is optically thick, and the accretion heat is primarily released near the midplane. As aresult, both the disc midplane temperature and the location of the pressure maximum are onlymarginally affected by stellar irradiation, and the inner disc is also convectively unstable. Aspreviously suggested, the inner disc is primarily ionized through thermionic and potassium ionemission from dust grains, which, at high temperatures, counteract adsorption of free chargesonto grains. Our results show that the location of the pressure maximum is determined by thethreshold temperature above which thermionic and ion emission become efficient. Key words: planets and satellites: formation – protoplanetary discs
Close-in super-Earths – planets with radii 1–4 R ‘ and orbital pe-riods shorter than „
100 days – are common around solar-type andlower mass stars (Fressin et al. 2013; Dressing & Charbonneau 2013,2015; Mulders et al. 2018; Hsu et al. 2019; Zink et al. 2019). Yet howthese planets form is still an open problem. One theory posits thatthey are born further away from the star, and subsequently migratethrough the protoplanetary disc to their present orbits (Terquem &Papaloizou 2007; Ogihara & Ida 2009; McNeil & Nelson 2010;Cossou et al. 2014; Izidoro et al. 2017, 2019; Bitsch et al. 2019).However, this hypothesis predicts that the planets should be water-rich, in contrast to the water-poor composition inferred from theirobserved radius distribution combined with atmospheric evolutionmodels (Owen & Wu 2017; Van Eylen et al. 2018; Wu 2019).An alternative proposal is that these super-Earths form at ornear their present orbits, out of solid material which migrates to theinner disc prior to planet formation. These solids are expected to ‹ E-mail: [email protected] arrive in the form of pebbles, which radially drift inwards from theouter disc due to gas drag (Hansen & Murray 2013; Boley & Ford2013; Chatterjee & Tan 2014; Hu et al. 2018; Jankovic et al. 2019).However, these pebbles need to be trapped, i.e., their radial driftmust be halted, in order for them to coalesce into planets instead ofdrifting into the star. Such trapping may occur as follows.In the inner protoplanetary disc, at the short orbital periodsat which close-in super-Earths are observed, the gaseous disc isthought to accrete via the magneto-rotational instability (MRI; Bal-bus & Hawley 1991). The MRI leads to turbulence in the disc,which drives viscous accretion, but its magnitude is sensitive to thedisc’s ionization state. The MRI is expected to be efficient in the hot,innermost parts of the disc (where it can be activated by thermalionization of trace alkali elements, at temperatures Á © a r X i v : . [ a s t r o - ph . E P ] F e b disk (interior to the pressure maximum) limits the grain size there byinducing high collisional velocities between grains, causing themto fragment. A decrease in the size of grains reduces their radialdrift, providing an alternative way to accumulate the grains arrivingfrom the outer disc (Jankovic et al. 2019). In general, therefore, thestructure of the inner disc, undergoing MRI-induced accretion, islikely to play a key role in the formation of close-in super-Earths.In previous work (Mohanty et al. 2018), we presented a modelof the inner disc in which the disc structure and accretion due tothe MRI are treated self-consistently. We obtained an inner discstructure in line with the expectations sketched above, and inferreda location for the pressure maximum consistent with the orbitaldistances of close-in super-Earths. However, this model includes anumber of simplifying assumptions about the disc’s physical andchemical structure. First, the disc is considered to be verticallyisothermal, and second, heating by stellar irradiation is neglected.In reality, the disc temperature will vary vertically, with a profile thatis particularly non-trivial when both accretion heating and stellarirradiation are accounted for (e.g. D’Alessio et al. 1998). Third,the model assumes a constant dust opacity, when in fact the latterdepends on the temperature and on grain properties and abundance.Fourth, Mohanty et al. (2018) assumed that the only source ofionization in the disc is thermal (collisional) ionization of potas-sium, and only considered gas-phase interactions. Potassium is in-deed a good representative of thermally-ionized species in the innerdisc due to its low ionization potential and high abundance (Desch& Turner 2015). However, ions and free electrons are also adsorbedonto the surfaces of dust grains in the disc, where they quickly re-combine. Dust can thus reduce the disc ionization level and suppressthe MRI (Sano et al. 2000; Ilgner & Nelson 2006; Wardle 2007;Salmeron & Wardle 2008; Bai & Goodman 2009; Mohanty et al.2013). Conversely, hot grains ( Á
500 K) can also emit electronsand ions into the gas phase (for electrons, this process is known asthermionic emission). Such temperatures are easily attainable in theinner disc, and thermionic and ion emission have been hypothesizedto be important sources of ionization there (Desch & Turner 2015).Neither of these grain effects are treated by Mohanty et al. (2018).Finally, while Mohanty et al. (2018) showed a posteriori that, in theMRI-accreting inner disc, X-ray ionization of molecular hydrogenmay be competitive with thermal ionization of potassium in supply-ing free electrons (due to the low gas surface density in this region),they did not actually include X-ray ionization in their model.In this paper we present a new model of an MRI-accreting,steady-state inner disc that addresses each of the shortcomings ofthe Mohanty et al. (2018) model listed above. The vertical structurein our model is calculated self-consistently from viscous dissipation(due to the MRI-induced viscosity), stellar irradiation, and radiativeand convective cooling, with realistic opacities due to dust grains.Ionization in the disc is determined by thermal ionization of potas-sium, thermionic and ion emission from dust grains, and ionizationof molecular hydrogen by stellar X-rays, cosmic rays and radionu-clides. In §2 we detail all the components of our model and themethods used to find self-consistent steady-state solutions for thedisc structure. We present our results in §3, discuss the relative im-portance of the various physical and chemical processes in §4, andsummarise our findings in §5.The aim of this paper is two-fold: to present the methodology ofour calculations, and to discuss the detailed physics of the inner discin the context of fiducial disc and stellar parameters. In a companionpaper (Paper II: Jankovic et al. in prep. ), we investigate how the innerdisc structure varies as a function of these parameters, and discussthe implications for the formation of super-Earths in the inner disc.
We consider a disc that is viscously accreting and in steady-state,i.e., has a constant mass accretion rate. Our model of the disc struc-ture is described in §2.1. The disc structure depends on the disc’sradiative properties, i.e. opacities, and on the viscosity. Our cal-culation of opacities is summarised in §2.2, and the prescriptionfor MRI-driven viscosity given in §2.3. The latter viscosity is afunction of the disc’s ionization state, calculated using a chemicalnetwork described in §2.4. The disc structure, opacities, ionizationand viscosity are calculated self-consistently at every point in thedisc, using numerical methods supplied in §2.5.The key parameters of our model are the steady-state gas ac-cretion rate 𝑀 through the disc, stellar mass 𝑀 ˚ , stellar radius 𝑅 ˚ ,stellar effective temperature 𝑇 ˚ , and value of the viscosity in the ab-sence of the MRI (the dead-zone viscosity). Additionally, the discopacities and ionization state, and thus the disc structure, dependon the properties of the dust; most importantly, on the dust-to-gasratio 𝑓 dg and the maximum dust grain size 𝑎 max .We assume that the disc only accretes viscously, due to theMRI. Disc accretion may additionally be partially driven by mag-netic winds, which could also affect the inner disc structure (e.g.,Suzuki et al. 2016). Therefore, inclusion of wind-driven accretion isan important issue; however, it is one which we do not tackle here. Our disc model largely follows the work of D’Alessio et al. (1998,1999). We consider a thin, axisymmetric, Keplerian, steady-statedisc that is viscously accreting. We assume that the disc is in verti-cal hydrostatic equilibrium, heated by viscous dissipation and stellarirradiation, and that energy is transported by radiation and convec-tion. Since the disc is vertically thin, we neglect energy transport inthe radial direction. Furthermore, at a given disc radius, our viscos-ity depends only on local conditions and the vertical mass column(see §2.3). As such, the disc structures at different radii are onlycoupled by the stellar irradiation, as it penetrates the disc along theline-of-sight to the central star.
In a thin Keplerian disc in vertical hydrostatic equilibrium, the gaspressure profile at any given radius follows from 𝑑𝑃𝑑𝑧 “ ´ 𝜌 Ω 𝑧, (1)where 𝑃 is the gas pressure, 𝜌 the gas volume density, Ω the Keple-rian angular velocity, and 𝑧 the height above the disc midplane. Weadopt the ideal gas law. The disc is heated by the viscosity that drives the accretion. Theviscosity 𝜈 is parametrized by 𝜈 “ 𝛼𝑐 { Ω , where 𝛼 is the viscosityparameter and 𝑐 s the isothermal sound speed (Shakura & Sunyaev1973). The local viscous dissipation rate at any location in the discis given by Γ acc “ 𝛼𝑃 Ω . (2) MNRAS , 1–21 (2020)
RI-active inner disc The flux generated by viscous dissipation that is radiated throughone side of the disc at any disc radius 𝑟 is 𝐹 acc “ 𝜋 𝑀 𝑓 𝑟 Ω , (3)where 𝑓 𝑟 “ ´ a 𝑅 in { 𝑟 comes from the thin boundary-layer condi-tion at the inner edge of the disc, 𝑅 in “ 𝑅 ˚ is the radius of the inneredge of the disc, and we are assuming a zero-torque inner boundarycondition (e.g. Frank et al. 2002).We also consider heating due to stellar irradiation. Stellar fluxpropagates spherically outwards from the star, and the resultantheating at any disc location depends on the attenuation of this fluxalong the line-of-sight to the star. Accounting for this, however,while simultaneously neglecting scattering of starlight and radialenergy transport within the disc (as we do in our 1+1D modelhere), leads to multiple equilibrium solutions or none at all (e.g.Chiang et al. 2001). Such behaviour does not appear in 2D discmodels (e.g. Dullemond 2002), which are nevertheless too complexfor our purposes here. Instead, we treat the irradiation heating ateach disc radius in isolation, by considering heating only due tothe stellar flux that impinges on the disc surface at that radius(at some grazing angle 𝜙 calculated self-consistently; see furtherbelow), and then propagates vertically towards the midplane (Calvetet al. 1992; Chiang & Goldreich 1997). In this framework, theattenuation, i.e. the optical depth to the stellar flux, is approximatedas 𝜏 irr « 𝜏 irr , z { 𝜇 , where 𝜏 irr , z is the optical depth in the verticaldirection, and 𝜇 ” sin 𝜙 . Local heating due to stellar irradiation isthen given by Γ irr “ 𝜅 ˚ P 𝜌 𝐹 irr 𝜇 𝑒 ´ 𝜏 irr , z { 𝜇 , (4)where 𝜅 ˚ P is the disc Planck opacity to stellar irradiation (see section2.2) and 𝐹 irr is the total incident stellar flux at the specified discradius. In an optically-thick disc, the latter is given by 𝐹 irr “ 𝜎 𝑆𝐵 𝑇 ˚ ˆ 𝑅 ˚ 𝑠 ˙ 𝜇, (5)where 𝑠 is the distance to the star (spherical radius) from the discsurface at the specified radius. Finally, the grazing angle 𝜙 is given by 𝜙 “ sin ´ 𝜋 𝑅 ˚ 𝑟 ` tan ´ 𝑑 log 𝑧 irr 𝑑 log 𝑟 𝑧 irr 𝑟 ´ tan ´ 𝑧 irr 𝑟 . (6)Here the first term is the value of 𝜙 for a flat disc and comes fromthe finite size of the stellar disc, and the other two terms are due todisc flaring. 𝑧 irr p 𝑟 q is the height above the disc midplane at whichthe stellar flux is absorbed; we take 𝑧 irr p 𝑟 q to be the height wherethe optical depth to the stellar irradiation is 𝜏 irr “ {
3. Specifically,in calculating the height 𝑧 irr , the optical depth 𝜏 irr is obtained byintegrating along the spherical radius 𝑠 to the star, as opposed to theapproximation used in the local heating term. We determine 𝑧 irr and 𝜙 self-consistently with the disc structure following the procedureby D’Alessio et al. (1999), as outlined in section 2.5. When calculating the factor 𝑓 𝑟 in eq. (3) we have assumed that the innerdisc edge is at the stellar radius, which would imply that the disc surfacecan only see a half of the stellar disc. If the disc surface can only see a halfof the stellar disc, the expression for the total absorbed stellar flux, eq. (5),should include an additional factor of 1 { In a thin disc, the optical depth to the disc’s own radiation is muchsmaller in the vertical direction than in the radial. Thus we expectradiative energy transport to be primarily vertical, and that is theonly direction we consider.The frequency-integrated moments of the radiative transferequation in the Eddington approximation (i.e. assuming that theradiation is isotropic, as is valid in the optically thick regime) andthe equation for energy balance are then 𝑑𝐹𝑑𝑧 “ Γ acc ` Γ irr , (7) 𝑑𝐽𝑑𝑧 “ ´ 𝜌𝜅 R 𝜋 𝐹, (8)4 𝜌𝜅 P p 𝜎 SB 𝑇 ´ 𝜋𝐽 q “ Γ acc ` Γ irr , (9)where 𝐹 and 𝐽 are the radiative flux and mean intensity respectively.We have also assumed here that the 𝐽 and 𝐹 weighted opacities canbe approximated by the Planck mean opacity 𝜅 P and the Rosselandmean opacity 𝜅 R respectively (following e.g. Hubeny 1990, see§2.2).Together with the ideal gas equation of state, equations (1, 7-9) form a closed set in 𝑃 , 𝐹 , 𝐽 and 𝑇 . Together with appropriateboundary conditions, they determine the disc vertical structure. Oneboundary condition is imposed at the disc midplane, where the flux 𝐹 p q “ 𝑧 surf above the midplane. Theboundary condition for the flux 𝐹 there is obtained by integratingeq. (7) from 𝑧 “ 𝑧 surf , from which it follows that 𝐹 p 𝑧 surf q “ 𝐹 acc ` 𝐹 irr . The boundary condition for the mean intensity 𝐽 isgiven by 𝐽 p 𝑧 surf q “ 𝜋 𝐹 p 𝑧 surf q . Finally, we assume that the gaspressure at the top of the disc has a small constant value, 𝑃 p 𝑧 surf q “ ´ dyn cm ´ , which is another boundary condition. The precisevalue of 𝑃 p 𝑧 surf q is arbitrary, but does not affect our results as longas it is sufficiently small. Overall, then, we have four boundaryconditions on three differential equations. Note that the temperature 𝑇 p 𝑧 surf q at the disc surface follows from the algebraic eq. (9), theboundary conditions on 𝐽 and 𝑃 , and the ideal gas law. If radiative energy transport yields a thermal structure such that thegradient ∇ “ 𝑑 ln 𝑇𝑑 ln 𝑃 is greater than the adiabatic gradient ∇ ad “p 𝛾 ´ q{ 𝛾 , then the gas is unstable to convection. In disc regionswhere this is the case, we assume that energy transport by convectionis efficient and the gas is vertically isentropic, so that ∇ “ ∇ ad atsuch locations (e.g. Shu 1992; Rafikov 2007; Garaud & Lin 2007).We adopt 𝛾 “ .
4, valid for an H dominated disc. Equations (8)and (9) are then replaced by 𝑑𝑇𝑑𝑧 “ ´ ∇ ad 𝑇𝑃 𝜌 p 𝑇, 𝑃 q Ω 𝑧 “ ´ ∇ ad 𝜇𝑚 H 𝑘 B Ω 𝑧. (10) Radiative transport is controlled by the Rosseland-mean opacity 𝜅 R in optically-thick regions (eq. 8), and by the Planck-mean opacity 𝜅 P in optically-thin regions (eq. 9). Additionally, the absorptioncoefficient for the stellar flux is a Planck-mean opacity 𝜅 ˚ P at thestellar effective temperature. We assume that the only source ofthese opacities are dust grains. Gas opacities are important in thevery innermost regions of discs, where most dust species have sub-limated. However, these regions are not of particular significance MNRAS000
4, valid for an H dominated disc. Equations (8)and (9) are then replaced by 𝑑𝑇𝑑𝑧 “ ´ ∇ ad 𝑇𝑃 𝜌 p 𝑇, 𝑃 q Ω 𝑧 “ ´ ∇ ad 𝜇𝑚 H 𝑘 B Ω 𝑧. (10) Radiative transport is controlled by the Rosseland-mean opacity 𝜅 R in optically-thick regions (eq. 8), and by the Planck-mean opacity 𝜅 P in optically-thin regions (eq. 9). Additionally, the absorptioncoefficient for the stellar flux is a Planck-mean opacity 𝜅 ˚ P at thestellar effective temperature. We assume that the only source ofthese opacities are dust grains. Gas opacities are important in thevery innermost regions of discs, where most dust species have sub-limated. However, these regions are not of particular significance MNRAS000 , 1–21 (2020) for the early stages of planet formation that we are interested in,as dust is required to form solid planet cores. Beyond the silicatesublimation line, gas opacities may still be important in hot, opti-cally thin regions (Malygin et al. 2014). However, as we find, theinner disc is significantly optically thick and including gas opacitieswould only alter the structure of the hot disc atmosphere. Therefore,we completely ignore the contribution of gas to the opacities.Furthermore, we assume that the only dust species present aresilicate grains. Other species that can be comparable in abundanceto silicates are water ice and carbonaceous grains (e.g. organics,Pollack et al. 1994). However, due to their low sublimation temper-atures, neither water ice nor carbonaceous grains are expected inthe hot MRI-active regions, and we generally limit our calculationsto the inner 1 AU of the disc.To calculate the opacities, we adopt the optical constants of“astronomical silicates” from Draine (2003). We use the miescatmodule of the python wrapper radmc3dPy for RADMC3D (Dulle-mond et al. 2012) to obtain dust absorption and scattering coef-ficients as functions of radiation wavelength and grain size. Foreach grain size 𝑎 , the coefficients are averaged within a size bin ofwidth Δ ln 𝑎 “ .
02. Next, we assume that the grain bulk densityis 𝜌 gr “ . ´ , and that the grain sizes are described by theMRN distribution, wherein the number density of grains in the sizerange r 𝑎, 𝑎 ` 𝑑𝑎 s is given by 𝑛 p 𝑎 q 𝑑𝑎 𝑎 ´ 𝑞 𝑑𝑎 (Mathis et al. 1977),with 𝑞 “ .
5. We adopt a minimum grain size of 𝑎 min “ . 𝜇 m,and a variable maximum grain size 𝑎 max . The size-dependent ab-sorption and scattering coefficients are then weighted by grain massand averaged over the grain size distribution. Finally, the absorp-tion coefficient is integrated over frequency to obtain the Planck-mean opacity 𝜅 P p 𝑇 q , and the total extinction coefficient yields theRosseland-mean opacity 𝜅 R p 𝑇 q . Following D’Alessio et al. (1998),we calculate the mean absorption coefficient for the stellar flux as afrequency-integrated absorption coefficient weighted by the Planckfunction at the stellar effective temperature: 𝜅 ˚ P ” 𝜅 P p 𝑇 ˚ q .Figure 1 shows the opacities per unit mass of gas, assum-ing a dust-to-gas ratio of 𝑓 dg “ ´ , a maximum grain size of 𝑎 max “ 𝜇 m, and 𝑇 ˚ “ 𝜅 ˚ P is a constant, since the wavelength-dependent dust absorption coefficient does not depend on the localtemperature. The Planck-mean opacity 𝜅 P is in general expectedto increase with increasing temperature. This is because the wave-length at which the Planck function peaks is inversely proportionalto the temperature, and for grains smaller than the wavelength ofpeak emission (and small grains contribute to the opacities most) ab-sorption is expected to increase with decreasing wavelength. How-ever, due to the silicate absorption feature at 10 𝜇 m, 𝜅 P decreaseswith temperature in the range „
500 K – 1000 K.
Our model of the MRI-driven viscosity closely follows that of Mo-hanty et al. (2018) (see also Bai 2011). Here we only summarizethe main points. A well-ionized circumstellar disc which followsthe laws of ideal magnetohydrodynamics (MHD) is susceptible tothe MRI (Balbus & Hawley 1991). The MRI leads to turbulence,producing an accretion stress and acting as a source of viscosity.The resulting Shakura-Sunyaev viscosity parameter 𝛼 is 𝛼 AZ “ 𝛽 , (11) Temperature [K]10 O p a c i t y [ c m / g ] PR*P
Figure 1.
Planck-mean opacity 𝜅 P , Rosseland-mean opacity 𝜅 R , and Planck-mean opacity at the stellar effective temperature 𝜅 ˚ P , as functions of disctemperature, assuming dust-to-gas ratio 𝑓 dg “ ´ and a maximum dustgrain size 𝑎 max “ 𝜇 m. where the subscript ‘AZ’ denotes an MRI-active zone, 𝛽 ” 𝑃 { 𝑃 𝐵 is the plasma parameter, and 𝑃 𝐵 “ 𝐵 { 𝜋 is the magnetic fieldpressure (Sano et al. 2004) .However, even in the inner regions of protoplanetary discs,non-ideal MHD effects can quench the MRI, leading to so-calleddead zones. The non-ideal effects considered here are Ohmic andambipolar diffusion. Ohmic diffusion will not suppress the MRI if(Sano & Stone 2002) Λ “ 𝑣 A 𝑧 𝜂 O Ω ą , (12)where Λ is the Ohmic Elsasser number, 𝑣 A 𝑧 ” 𝐵 𝑧 { a 𝜋𝜌 is thevertical component of the local Alfven velocity and 𝜂 O is the Ohmicresistivity. Here we utilize the relationship between the strength ofthe vertical component of the magnetic field, 𝐵 𝑧 , and the strength ofthe r.m.s. field, 𝐵 : 𝐵 𝑧 „ 𝐵 {
25 (Sano et al. 2004). We also assumethat 𝐵 is vertically constant. Our method of determining the valueof 𝐵 is described in section 2.5.Similarly, the ambipolar Elsasser number is defined by 𝐴𝑚 “ 𝑣 A 𝜂 A Ω , (13)where 𝜂 A is the ambipolar magnetic resistivity. However, in thestrong-coupling limit, valid in protoplanetary discs, the MRI can beactive even if 𝐴𝑚 ă
1, as long as the magnetic field is sufficientlyweak (Bai & Stone 2011). The criterion for active MRI in thepresence of ambipolar diffusion is 𝛽 { 𝛽 min ą , (14)where the minimum value of 𝛽 neccessary to sustain the MRI is afunction of the ambipolar Elsasser number: 𝛽 min p 𝐴𝑚 q “ «ˆ 𝐴𝑚 . ˙ ` ˆ 𝐴𝑚 . ` ˙ ff { . (15)Thus, whether the MRI is active or not depends on the magneticresistivities, 𝜂 O and 𝜂 A , as well as on the magnetic field strength The numerical factor is indeed 1 { 𝛼 parameterin eq. (11); see discussion in Appendix B of Mohanty et al. (2018).MNRAS , 1–21 (2020) RI-active inner disc 𝐵 . The magnetic resistivities, calculated following Wardle (1999),express the coupling between the gas and the magnetic field, whichis principally determined by the degree of ionization of the gas.If either of the two criteria given by eqns. (12) and (14) is notfulfilled, the MRI is not active. In such MRI-dead zones, we assumethere is a small residual viscosity 𝛼 DZ , driven either by propagationof turbulence from the MRI-active zone, or by hydrodynamic insta-bilities (for more discussion, see Mohanty et al. 2018). In this work,we also impose a smooth transition between the active and the deadzones, necessary to ensure numerical stability in the integration ofthe equations of disc structure (such smoothing was not employedby Mohanty et al. 2018). To cover both non-ideal effects that lead todead zones, we define 𝐶 ” min p Λ , 𝛽 { 𝛽 min q . Then, at any locationin the disc where | 𝐶 ´ | ă . 𝛼 given by 𝛼 “ 𝛼 DZ ` 𝛼 AZ ´ 𝛼 DZ ` exp ´ ´ 𝐶 ´ Δ ¯ , (16)where Δ “ ´ . We implement the same simple chemical network adopted by Desch& Turner (2015). This network tracks the number densities of fivespecies: free electrons ( 𝑛 e ), potassium ions ( 𝑛 K ` ), neutral potas-sium atoms ( 𝑛 K ), potassium atoms adsorbed (condensed) onto dustgrains ( 𝑛 K , cond ), and atomic ions ( 𝑛 i ; i.e., ions of atomic speciesother than potassium; see below).In the gas-phase, potassium atoms can be thermally ionized:collisions of neutral potassium atoms with H molecules producepotassium ions and free electrons. The potassium ions and freeelectrons can also recombine in the gas phase, either radiativelyor via three-body collisions with H molecules (the latter processdominates at the high densities prevalent in the inner disc).Furthermore, non-thermal sources (we consider stellar X-rays,cosmic rays and radionuclides; see section 2.4.3) can ionize H (Glassgold et al. 1997; Ercolano & Glassgold 2013). The chargeis quickly transferred from the ionized hydrogen to other abundantgas species through collisions, producing molecular and atomicions (e.g. HCO ` , Mg ` ). Notably, in application to the MRI, theexact composition of the gas that this leads to is unimportant in thepresence of dust, and simple chemical networks reproduce the gasionization levels well (Ilgner & Nelson 2006). Thus, it is assumedthat the ionization of molecular hydrogen at a rate 𝜁 directly pro-duces atomic ions and free electrons at a volumetric rate 𝜁𝑛 H . Theatomic ion species in this chemical network may thus be understoodas a representative of the various chemical species abundant in thegas-phase, whose mass is taken to be that of magnesium. It is as-sumed that the number density of molecular hydrogen is constant,which is valid for low ionization rates. Just as potassium, the atomicions also recombine radiatively and in three-body recombinations.Importantly, all gas-phase species (electrons, ions and neutralatoms) collide with and are adsorbed onto dust grains at a rate R 𝑗 , coll “ 𝑛 𝑗 𝑛 gr 𝜋𝑎 ˆ 𝑘 B 𝑇𝜋𝑚 𝑗 ˙ { ˜ 𝐽 𝑗 𝑆 𝑗 , (17)where 𝑛 𝑗 is the number density of the gas-phase species, 𝑛 gr isthe number density of the grains, 𝑎 gr is the grain size, 𝑚 𝑗 is thegas-phase species mass, ˜ 𝐽 𝑗 is the modification of the collisional cross-sections for charged species due to dust grain charge (Draine& Sutin 1987), and 𝑆 𝑗 is the sticking coefficient. It is assumed thatall grains have the same charge; this is valid since the dispersionin the distribution of charge states is generally found to be small(Draine & Sutin 1987). It is further assumed that potassium ionsrapidly recombine on the grain surface to form condensed potassiumatoms. Thus the ions are effectively destroyed upon adsorption.At high temperatures electrons on the dust grains have a finiteprobability of leaving the grain, producing so-called thermionicemission. The emission depends on the energy required for theelectron to escape the grain. For a neutral grain this is the workfunction 𝑊 , a property of the material out of which the grainsare made. The rate at which free electrons are produced throughthermionic emission is R therm “ 𝑛 gr 𝜋𝑎 𝜆 R 𝜋𝑚 e p 𝑘 B 𝑇 q ℎ exp ˆ ´ 𝑊 eff 𝑘 B 𝑇 ˙ , (18)where 𝑊 eff “ 𝑊 ` 𝑍𝑒 𝑎 gr (19)is the effective work function due to grain charge 𝑍𝑒 .Potassium atoms will also evaporate from the grains only athigh temperatures. The vaporization rate of condensed potassiumatoms is given by R K , evap “ 𝑛 K , cond 𝜈 exp ˆ ´ 𝐸 a 𝑘 B 𝑇 ˙ , (20)where 𝜈 is the vibration frequency of potassium atoms on the dustgrain surface lattice, and 𝐸 a “ .
26 eV is the binding energy, whosevalue is chosen to reproduce the condensation temperature of potas-sium (1006 K, Lodders 2003). These potassium atoms may be emit-ted into the gas phase as both neutral atoms and ions, contributingfurther to the gas’ ionization state. The ratio of ions to neutralsamong the emitted particles is given by 𝑛 ` K 𝑛 K “ 𝑔 ` 𝑔 exp ˆ ` 𝑊 eff ´ IP 𝑘 B 𝑇 ˙ , (21)where 𝑔 ` 𝑔 is the ratio of statistical weights of the ionized and neutralstate of potassium, and IP the ionization potential of potassium. Wedenote by 𝑓 ` the fraction of all emitted particles that leave the grainas ions.It is important to note here that there are two sources of theneutral potassium atoms condensed on grains, whose evaporation isdescribed above: first, gas-phase potassium ions that are adsorbedonto grain surfaces and recombine there into neutral atoms; andsecond, gas-phase neutral potassium atoms that are adsorbed di-rectly onto the grains. The reionisation and subsequent ejectionfrom grains of the former simply returns ions to the gas-phase thatwere originally adsorbed from it, and thus clearly cannot increasethe ion fraction in the gas-phase beyond what it would be in theabsence of grains. The ionisation on the grain surface, and ejectionas ions, of particles that were originally neutral in the gas-phase,however, can increase the gas-phase ion fraction beyond what itwould be without grains; it is this channel that makes ion emissionfrom grains such a crucial effect.Clearly, the contribution of the dust grains to gas ionization lev-els depends on the work function 𝑊 of the grain material. We adoptthe fiducial value of Desch & Turner (2015), 𝑊 “ “ .
34 eV, indicating that thermionic emission is
MNRAS000
MNRAS000 , 1–21 (2020) important for the production of free electrons in the same tempera-ture range as thermal ionization of potassium. Importantly, for thisgiven value of the work function, grains become negatively chargedat high temperatures as a large fraction of potassium evaporatingfrom the grains is in ionized state. This results in a reduction of theeffective work function 𝑊 eff , since, for negatively charged grainsthermionic emission is higher (eq. 19).In this work, we assume that the abundances of hydrogen andpotassium atoms are 𝑥 H “ . ˆ ´ and 𝑥 K “ . ˆ ´ ,respectively, relative to the total number density of all atomic parti-cles (Keith & Wardle 2014), and that the mean molecular weight is 𝜇 “ . 𝑚 H3 . The grain material density is 𝜌 gr “ . ´ , thesame as in our calculation of the dust opacities. The input for thechemical network are temperature 𝑇 , pressure 𝑃 , hydrogen ioniza-tion rate 𝜁 , dust-to-gas ratio 𝑓 dg and dust grain size 𝑎 gr . All otherkinetic rates, parameters and coefficients are the same as in Desch& Turner (2015). We note that the chosen value of the stickingcoefficient for electrons ( 𝑆 𝑒 “ .
6) is compatible with the detailedcalculation by Bai (2011), who consider work function values of1 eV and 3 eV. His results suggest that for a work function of 5 eV, at1000 K, 𝑆 𝑒 is indeed a few times 0.1 for neutral grains, and increasesfurther for negatively-charged grains. For the sticking coefficientsfor the ions we adopt 𝑆 𝑗 “
1, and similarly for the neutral atomswe adopt a sticking coefficient of 𝑆 “ 𝑓 dg and grain size 𝑎 gr , we pre-calculate and tabulate the equilibrium number densities of electronsand ions, and the average grain charge, as functions of temperature 𝑇 ,pressure 𝑃 and hydrogen ionization rate 𝜁 . We find the equilibriumsolution following the same method as Desch & Turner (2015).Time derivatives of all number densities are set to zero, and sorate equations yield an algebraic system of equations. For a givenaverage grain charge 𝑍 , this system of equations is solved iterativelyto find number densities of all five species. The grain charge is thenfound by solving the equation of charge neutrality. The described chemical network incorporates only one grain sizepopulation. Ideally, we would consider a number of grain size pop-ulations, with the same size distribution used in our calculation ofdust opacities. However, this would greatly enhance the computa-tional complexity of the problem. At the same time, it is clear thatdust grains of different size contribute differently to the equilib-rium ionization levels. To the lowest order of approximation, alldust-related reaction rates are regulated by the total grain surfacearea. Thus, it can be expected that the ionization levels are mostsensitive to the smallest grains. Bai & Goodman (2009) consideredthe effects of dust on the ionization levels in the cold regions ofprotoplanetary discs, in application to the onset of the MRI dueto non-thermal sources of ionization. They considered chemicalnetworks with two grain size populations and found that the grainpopulations behave independently, as the charge transfer betweenthe grains is negligible. They further found that the ionization levelsare largely controlled by a quantity 𝑓 dg { 𝑎 𝑝 gr , where the exponent 𝑝 varies between 𝑝 “ 𝑝 “ The total number densities of molecular hydrogen and potassium arethen related to the gas density as 𝑛 H “ 𝑥 H {p ´ 𝑥 H q 𝜌 { 𝜇 and 𝑛 K “ 𝑥 K {p ´ 𝑥 H q 𝜌 { 𝜇 , respectively. Note that the fiducial value of the hydrogen ionization rate used by Desch& Turner (2015) is in fact 𝜁 “ . ˆ ´ s ´ (S. Desch, priv. comm.). We repeat a similar exercise for the above chemical network,suited to the hot inner regions of protoplanetary discs. For a set ofvalues of the exponent 𝑝 = 1, 1.25 and 1.5, we vary the grain size 𝑎 gr and dust-to-gas mass ratio 𝑓 dg while keeping 𝑓 dg { 𝑎 𝑝 gr constant.We calculate the ionization levels as a function of temperature,and for different sets of pressure and hydrogen ionization rates, soas to probe different conditions in different regions of the innerdisc. The results are shown in Fig. 2. The grain surface area ( 𝑝 “
1) controls the equilibrium ionization levels when the hydrogenionization dominates, but does not determine the temperature atwhich the ionization levels rise due to thermionic and ion emission.On the other hand, for 𝑝 “ .
5, this temperature depends veryweakly on the grain size. That is, regardless of the actual dust grainsize, a quantity 𝑓 dg { 𝑎 . regulates the temperature at which denseinterior of the disc becomes ionized.Therefore, for the chemical network calculations, we use a sin-gle grain size of 𝑎 gr “ ´ cm, but employ an effective dust-to-gasratio 𝑓 dg , eff that satisfies 𝑓 dg , eff 𝑎 ´ 𝑝 gr “ ş 𝑎 max 𝑎 min 𝑑𝑛 p 𝑎 q 𝑚 p 𝑎 q{ 𝜌 g 𝑎 ´ 𝑝 ,where 𝑛 p 𝑎 q is the same grain size distribution used to calculate dustopacities. Since we find that thermionic and ion emission are moreimportant than the non-thermal ionization of hydrogen in the innerdisc, we use 𝑝 “ .
5. From the right panel of Fig. 2, it appears thatthis choice will lead to a large error in the gas ionization state wederive in the low-density non-thermally ionized disc regions, due tothe large variation in such regions in the ionization state producedby different grain sizes . However, in reality the error will be muchsmaller than implied by Fig. 2, since the majority of grains areskewed towards small sizes in a realistic grain size distribution.
In our calculation of the MRI-driven viscosity, we consider molec-ular hydrogen ionization rate due to radionuclides, cosmic rays andstellar X-rays. The ionization rate of molecular hydrogen due toshort-lived and long-lived radionuclides is 𝜁 R “ . ˆ ´ s ´ , (22)which predominantly comes from decay of Al (Umebayashi &Nakano 2009). The ionization rate of molecular hydrogen due tointerstellar cosmic rays is (Umebayashi & Nakano 2009) 𝜁 CR p 𝑧 q “ 𝜁 CR , ISM 𝑒 ´ Σ p 𝑧 q 𝜆 CR ˜ ` ˆ Σ p 𝑧 q 𝜆 CR ˙ ¸ ´ , (23)where 𝜁 CR , ISM “ ´ s ´ is the interstellar cosmic ray ionizationrate, Σ p 𝑧 q is the integrated density column from the top of the discto the height 𝑧 above disc midplane, and 𝜆 CR “
96 g cm ´ is theattenuation length for cosmic rays (Umebayashi & Nakano 1981).For the ionization rate of molecular hydrogen due to stellarX-rays we use Bai & Goodman (2009) fits to the Igea & Glassgold(1999) Monte Carlo simulations, 𝜁 X p 𝑧 q “ 𝐿 X erg s ´ ´ 𝑟 ¯ ´ . p 𝜁 𝑒 ´p Σ p 𝑧 q{ 𝜆 q 𝑐 ` (24) 𝜁 𝑒 ´p Σ p 𝑧 q{ 𝜆 q 𝑐 q , where 𝐿 X is stellar X-ray luminosity, 𝜁 “ ˆ ´ s ´ , 𝜆 “ . ˆ ´ g cm ´ , and 𝑐 “ . 𝜁 “ ´ s ´ , 𝜆 “ . ´ , and 𝑐 “ . MNRAS , 1–21 (2020)
RI-active inner disc terms of column densities of hydrogen nucleus into the surface den-sity lengths using the hydrogen abundance given above. We adoptthe saturated relationship 𝐿 X “ ´ . 𝐿 bol (e.g. Wright et al.2011). For both cosmic rays and X-rays, we ignore the contributioncoming through the other side of the disc. This is valid since we findthat the gas surface densities are mostly larger than the attenuationlengths of the ionizing particles and, even at low gas surface densi-ties, this can only increase the ionization rates by at most a factorof 2. At a given orbital radius, magnetic field strength and grazing angle 𝜙 , the disc’s vertical structure is determined as a solution to theboundary value problem given by eqns. (1), (7)-(9) and the idealgas law, and, where convectively unstable, by eqns. (1), (7) and (10).This boundary value problem is solved using the shooting method(Press et al. 2002). The equations are integrated from the top ofthe disc ( 𝑧 “ 𝑧 surf ) to the disc midplane ( 𝑧 “ 𝑧 surf is then found, such that 𝐹 p q “
0. This root-findingproblem is solved using the Ridders’ method, with an exit criterionthat | 𝐹 p q| ă ´ 𝐹 p 𝑧 surf q .Equations of the vertical disc structure are integrated on afixed vertical grid (i.e., in the 𝑧 -direction), with points uniformin the polar angle. Since eq. (9) is an algebraic equation, and fornumerical stability, we use a fully implicit integration method. Forthe radiative-transfer problem we use a Runge-Kutta method ofthe second order, i.e. the trapezoidal method. In the trapezoidalmethod the equations (1), (7)-(9) are discretised as a system ofnonlinear equations to be solved in every integration step (a systemof nonlinear equations in 𝑃 𝑛 ` , 𝐹 𝑛 ` , 𝐽 𝑛 ` , 𝑇 𝑛 ` to be solved bya root-finding algorithm), 𝑃 𝑛 ` “ 𝑃 𝑛 ` ℎ Ω p´ 𝜌 𝑛 𝑧 𝑛 ´ 𝜌 𝑛 ` 𝑧 𝑛 ` q ,𝐹 𝑛 ` “ 𝐹 𝑛 ` ℎ p Γ 𝑛 ` Γ 𝑛 ` q ,𝐽 𝑛 ` “ 𝐽 𝑛 ` ℎ p´ 𝜋 qp 𝜌 𝑛 𝜅 R p 𝑇 𝑛 q 𝐹 𝑛 ` 𝜌 𝑛 ` 𝜅 R p 𝑇 𝑛 ` q 𝐹 𝑛 ` q , “ 𝜌 𝑛 ` 𝜅 𝑃 p 𝑇 𝑛 ` qp 𝜎 SB 𝑇 𝑛 ` ´ 𝜋𝐽 𝑛 ` q ´ Γ 𝑛 ` , where ℎ is the integration step, 𝜌 “ 𝜌 p 𝑇, 𝑃 q is given by the ideal gaslaw and Γ “ Γ acc ` Γ irr . Here, Γ irr ,𝑛 “ Γ irr ,𝑛 p 𝑇 𝑛 , 𝑃 𝑛 , 𝜏 irr ,𝑧,𝑛 q , withthe optical depth to stellar irradiation obtained using 𝜏 irr ,𝑧,𝑛 ` “ 𝜏 irr ,𝑧,𝑛 ` ℎ 𝜅 ˚ P p 𝜌 𝑛 ` 𝜌 𝑛 ` q . Furthermore, viscous dissipation isa function of the MRI-driven viscosity 𝛼 and thus depends onthe local ionization levels. The latter are a function of the lo-cal temperature, pressure and the hydrogen ionization rate whichdepends on the column density from the top of the disc. Thus, Γ acc ,𝑛 “ Γ acc ,𝑛 p 𝑇 𝑛 , 𝑃 𝑛 , 𝑁 𝑛 q , with the column density given by 𝑁 𝑛 ` “ 𝑁 𝑛 ` ℎ 𝑘 ´ p 𝑃 𝑛 { 𝑇 𝑛 ` 𝑃 𝑛 ` { 𝑇 𝑛 ` q .The equation for pressure 𝑃 𝑛 ` can be rearranged into anexplicit form 𝑃 𝑛 ` “ 𝑃 𝑛 ´ ℎ Ω 𝜌 𝑛 𝑧 𝑛 ` ℎ Ω 𝜇𝑚 H 𝑘 B 𝑧 𝑛 ` 𝑇 𝑛 ` . Then, the above system of equations is equivalent to a single non-linear equation in 𝑇 𝑛 ` , greatly simplifying the problem. In everyintegration step we use the Ridders’ method to solve this equationfor the temperature 𝑇 𝑛 ` (down to a relative precision of 10 ´ ) and consequently for all other quantities. This includes the MRI-drivenviscosity 𝛼 , which is thus calculated self-consistently at each stepof integration. At every step and in every iteration of the root-solver opacities are interpolated from pre-calculated tables usingcubic splines, and the ionization levels (e.g., free electron numberdensity) using tri-linear interpolation.Additionally, at each integration step we check if the result-ing temperature gradient is unstable to convection, and if so, thetemperature 𝑇 𝑛 ` is obtained analytically using 𝑇 𝑛 ` “ 𝑇 𝑛 ` ℎ ∇ ad 𝜇𝑚 H 𝑘 B Ω p´ 𝑧 𝑛 ´ 𝑧 𝑛 ` q . With 𝑇 𝑛 ` known, all other quantities follow same as above. A disccolumn can, in principle, become convectively stable again at someheight above disc midplane. To calculate the mean intensity 𝐽 ata boundary between a convective and a radiative zone, we use theenergy balance equation (9).For some model parameters there can be a range of orbitalradii and values of the magnetic field strength for which there aremultiple solutions for the disc height 𝑧 surf (i.e. multiple solutionsfor the equilibrium vertical disc structure). This happens when acomplex ionization structure leads to strong variations in the vis-cosity 𝛼 as a function of height above disc midplane, making thetotal produced viscous dissipation a non-monotonous function of 𝑧 surf . When there are multiple solutions, we choose the solution withsmallest 𝑧 surf . It is likely that at least some of the additional solutionsare thermally unstable and/or unphysical, as the strong variationsin both the levels of turbulence and the levels of ionization shouldbe removed by turbulent mixing (see section 3.2.3). Note that themultiple solutions in our model, when they exist, correspond to thesame, radially-constant input accretion rate. Steady accretion is thusassured regardless of the choice of the solution. At a given orbital radius, for a fixed grazing angle, 𝜙 , and magneticfield strength, 𝐵 , the above procedure yields an equilibrium verticaldisc structure characterised by a vertically-averaged viscosity¯ 𝛼 “ ş 𝑧 surf 𝛼𝑃𝑑𝑧 ş 𝑧 surf 𝑃𝑑𝑧 .
As in the vertically-isothermal model (Mohanty et al. 2018), fora sufficiently small and a sufficiently large magnetic field strength 𝐵 the MRI is suppressed in the entire disc column and ¯ 𝛼 “ 𝛼 DZ .There can be an intermediate range of magnetic fields strengthsfor which the MRI is active, and the vertically-averaged viscosity¯ 𝛼 peaks at some value of 𝐵 . At every orbital radius, we choose 𝐵 such that ¯ 𝛼 is maximized. The underlying assumption here is thatthe magnetic fields are strengthened by the MRI-driven turbulence.It is assumed that the initial magnetic field configuration was suf-ficient to start the instability. The induced magnetohydrodynamicturbulence then amplifies the field strength. We assume the equi-librium configuration is one in which the turbulence is maximized,where the turbulence is parametrized by ¯ 𝛼 . Similar arguments wereemployed by Bai (2011) and Mohanty et al. (2013).To maximize ¯ 𝛼 p 𝐵 q , we use the Brent method with a target This is indeed necessary. An iterative method in which the disc thermalstructure is decoupled from the density structure and the heating terms(e.g. Dullemond et al. 2002) does not converge to a solution in the case ofMRI-driven viscosity.MNRAS000
As in the vertically-isothermal model (Mohanty et al. 2018), fora sufficiently small and a sufficiently large magnetic field strength 𝐵 the MRI is suppressed in the entire disc column and ¯ 𝛼 “ 𝛼 DZ .There can be an intermediate range of magnetic fields strengthsfor which the MRI is active, and the vertically-averaged viscosity¯ 𝛼 peaks at some value of 𝐵 . At every orbital radius, we choose 𝐵 such that ¯ 𝛼 is maximized. The underlying assumption here is thatthe magnetic fields are strengthened by the MRI-driven turbulence.It is assumed that the initial magnetic field configuration was suf-ficient to start the instability. The induced magnetohydrodynamicturbulence then amplifies the field strength. We assume the equi-librium configuration is one in which the turbulence is maximized,where the turbulence is parametrized by ¯ 𝛼 . Similar arguments wereemployed by Bai (2011) and Mohanty et al. (2013).To maximize ¯ 𝛼 p 𝐵 q , we use the Brent method with a target This is indeed necessary. An iterative method in which the disc thermalstructure is decoupled from the density structure and the heating terms(e.g. Dullemond et al. 2002) does not converge to a solution in the case ofMRI-driven viscosity.MNRAS000 , 1–21 (2020) Temperature [K]10 n e / n H f dg / a gr = const Temperature [K]10 n e / n H f dg / a = const Temperature [K]10 n e / n H f dg / a = const Figure 2.
Ionization fraction ( 𝑛 e { 𝑛 H ) as a function of temperature. Different colors correspond to different combinations of pressure and hydrogen ionizationrates: 𝑃 “ dyn cm ´ and 𝜁 “ ´ s ´ (blue), 𝑃 “ ´ dyn cm ´ and 𝜁 “ ´ s ´ (orange), 𝑃 “ ´ dyn cm ´ and 𝜁 “ ´ s ´ (green).Different linestyles correspond to different grain sizes, 𝑎 gr “ ´ cm (solid), 𝑎 gr “ ´ cm (dashed), 𝑎 gr “ ´ cm (dotted), each with a differentdust-to-gas ratio such that the ratio 𝑓 dg { 𝑎 𝑝 gr remains constant and as evaluated for 𝑎 gr “ ´ cm, 𝑓 dg “ .
01. Different panels show the calculations fordifferent values of the exponent 𝑝 , as indicated in each panel title. Exponent 𝑝 “ 𝑝 “ . 𝑎 gr . See Section2.4.2. absolute precision of 10 ´ in log 𝐵 . There can be a range of orbitalradii where there are multiple local maxima in ¯ 𝛼 as a function of 𝐵 .This is essentially for the same reasons that cause multiple solutionsin disc height 𝑧 surf at a fixed value of 𝐵 . In general we choose 𝐵 corresponding to the global maximum in ¯ 𝛼 . However, in some cases,this is a function of the grazing angle 𝜙 at a fixed orbital radius,and the procedure to determine the grazing angle (described below)does not converge. There we choose a local maximum with a largestmagnetic field strength. At any given orbital radius, the angle between the incident stellarradiation and the disc surface is determined self-consistently withthe disc structure following D’Alessio et al. (1999). A self-consistentdisc structure is found by iteratively updating the grazing angle andre-calculating the entire disc structure. We use a logarithmic grid fororbital radius. First, we calculate 𝜙 using eq. (6) by assuming that 𝑧 irr “ 𝜏 irr p 𝑟, 𝑧 q , whichyields an updated 𝑧 irr at each radius.Critically, to calculate the updated value of the grazing angle 𝜙 at each radius, the derivative 𝑑 log 𝑧 irr { 𝑑 log 𝑟 is approximated byassuming that 𝑧 irr is a power-law, 𝑧 irr 𝑟 𝑏 , within a radius bincentered at the given radius. So, at each radius we fit for the slope 𝑏 using 𝑧 irr at that radius and at a number of radial grid points interiorand exterior to it. Then, the value of the grazing angle is updatedand the vertical disc structure re-calculated at all radii.This procedure is repeated until the grazing angle has con-verged at every radius, i.e., until the relative difference in 𝜙 betweentwo consecutive iterations is less than 0.5% at all radii. For the firstradial point we always assume the flat-disc approximation ( 𝑧 irr “ 𝑑 log 𝑧 irr { 𝑑 log 𝑟 . Using our self-consistent model we can now investigate the struc-ture of the inner regions of a protoplanetary disc that is viscouslyaccreting with a constant accretion rate.In section 3.1 we consider a model in which the vertical struc-ture of the disc is calculated self-consistently from viscous heating,heating by stellar irradiation, radiative and convective energy trans-port, and self-consistent radiative properties of dust, while the disc’sionization state is calculated by only considering the thermal (col-lisional) ionization of potassium, using the Saha equation.In section 3.2 we present results for our full model that includesadditional chemical species, including dust grains, in the chemicalnetwork, as well as non-thermal sources of ionization.Unless otherwise stated, throughout this paper we assume asolar-mass star, 𝑀 ˚ “ d , with a stellar radius 𝑅 ˚ “ d ,effective stellar temperature 𝑇 ˚ “ 𝑀 “ ´ M d yr ´ , and viscosity in the MRI-dead zone 𝛼 DZ “ ´ .Our adopted gas accretion rate is the median from observations ofsolar-mass classical T Tauri stars (e.g., Hartmann et al. 1998; Ma-nara et al. 2016, 2017). In reality, there is more than an order ofmagnitude scatter in the data around this median; we explore theeffects of varying 𝑀 , as well as other parameters, in a companionpaper, as noted below. The stellar parameters are from the evolu-tionary models of Baraffe et al. (2015), for a solar-mass star at anage of 5 ˆ yr. At later times, the stellar luminosity decreases, asthe star contracts towards the zero-age main-sequence. Our choiceof a relatively young classical T Tauri star thus maximises the stel-lar luminosity, allowing us to examine the largest possible effect ofstellar irradiation on the inner disc. Lastly, we adopt a standard ISMdust-to-gas ratio of 𝑓 dg “ ´ , and a maximum dust grain size of 𝑎 max “ 𝜇 m.The above set of parameters represents our fiducial model . Inthis work we investigate in detail the main physical and chemicalprocesses that shape the structure of the inner disc for this model. Ina companion paper (Jankovic et al. in prep. ), we examine how theinner disc structure changes as a function of the model parameters,and build a picture of how this structure may affect planet formation. MNRAS , 1–21 (2020)
RI-active inner disc In this section we explore how thermal processes shape the inner discstructure. The ionization state of the disc is set here only by the ther-mal ionization and recombination of potassium, which is assumedto be entirely in the gas phase. Note that there is no other chemistryin this model: in the absence of grain-related reactions (grains hereonly contribute to the opacity), and without non-thermal sourcesto ionize H , gas-phase thermal ionization and recombination ofpotassium are the only reaction pathways available.Mohanty et al. (2018) presented a model in which the discis vertically isothermal, with a constant opacity (=10 cm g ´ ).Throughout this section, we shall compare this simple model withthree models of increasing complexity within our framework: (i) amodel with the same constant opacity ( 𝜅 R “ 𝜅 P “
10 cm g ´ ),but with the disc vertical structure calculated self-consistently fromviscous heating and cooling by radiation and convection; (ii) a modelwith the same heating and cooling processes as (i), but now withopacities also determined self-consistently; and (iii) similar to (ii),but now also including heating by stellar irradiation.To make it easier to interpret the effects of these thermal pro-cesses, we first focus on how the disc’s temperature changes whenthe increasingly complex models are used. We then discuss howthis leads to differences in the disc’s ionization state, the locationswhere the MRI is active, and the radial surface density and midplanepressure profiles. Finally we provide a short summary of the keyfindings at the end of this section. Figure 3 shows the temperature as a function of disc radius andheight above the midplane for our three models, increasing in com-plexity from the left to right panels.Noticeably, the temperature profiles deviate from being ver-tically isothermal. The dashed lines here show the disc pressurescale height, and the solid lines show the disc photosphere for out-going radiation (i.e., where the Rosseland-mean optical depth is 𝜏 R “ { general feature of any active, optically thick disc where the ac-cretion heat is released near the midplane (see also Garaud & Lin2007). We show this analytically in section 4.The model that includes heating by stellar irradiation, shown inthe right panels of Figs. 3 and 4, features a temperature inversion inthe disc upper layers. This inversion has been discussed in detail byD’Alessio et al. (1998), and is a simple consequence of the fact thatthe photosphere for the outgoing radiation (solid line) lies belowthe photosphere for the incoming stellar radiation (dotted line),where the latter corresponds to the irradiation surface 𝑧 irr at which 𝜏 irr “ {
3. Above this surface the disc upper layers are heatedby stellar irradiation. Below, the disc becomes optically thick tostellar radiation (i.e., stellar photons do not penetrate here) and thetemperature drops. Going deeper still, below the disc photosphereto the outgoing radiation, the disc becomes optically thick to its ownradiation and the temperature rises again. In the irradiated disc, close to the star the disc midplane isas hot as the disc upper layers, and further away the midplaneis significantly hotter than the upper layers. It would thus appearthat accretion heating dominates in the inner disc. However, at agiven radius, the total flux of stellar radiation absorbed by a verticalcolumn in the disc is in fact at least an order of magnitude greaterthan the total flux generated by viscous dissipation ( 𝐹 irr Á 𝐹 acc ),throughout the inner disc (see top panel of Fig. 5). In spite of this,the disc temperature near the midplane is only weakly affected byirradiation. We discuss this result in section 4.2.The ratio 𝐹 irr { 𝐹 acc varies non-monotonically, following thegrazing angle 𝜙 , shown in the bottom panel of Fig. 5. At the inneredge of our calculation domain, the first „
10 points (shown in gray)are affected by boundary effects. This is a well-known problem indisc models that account for stellar irradiation using the grazingangle prescription (Chiang et al. 2001). Importantly, far from theinner edge, the value of the grazing angle and the disc structure donot depend on the disc structure at the inner edge.
The ionization structure in our models with a self-consistent verticalstructure is qualitatively different from that in the vertically isother-mal model of Mohanty et al. (2018). In the latter, the temperature isconstant with height above the midplane while the density decreases,and thus the ionization fraction ( 𝑛 e { 𝑛 H ) increases with height (bythe Saha equation). In our models without stellar irradiation heat-ing, the temperature declines monotonically with height, and hencethe ionization fraction decreases as well (Fig. 6: middle row, solidcurves). When we include irradiation heating, on the other hand,the fractional ionization does increase in the uppermost irradiatedlayers, due to the rising temperature there (dashed curves).As a consequence of the above, both ambipolar and Ohmicresistivities increase with height in our non-irradiated discs. In otherwords, in these models, both ambipolar and Ohmic diffusion quenchthe MRI from above. In the irradiated disc, the Ohmic resistivityfalls in the uppermost layers; however, ambipolar diffusion stillincreases in these layers, again stifling the MRI there. Consequently,the vertical extent of the MRI-active region (i.e., where the local 𝛼 > 𝛼 DZ ) at any given radius in our irradiated disc model is nearlyidentical to that in the non-irradiated discs (Fig. 6, top row). The differences in the ionization structure of our vertically self-consistent models and the vertically isothermal model lead to dif-ferences in where the MRI is active in the disc. A plot of the viscosityparameter 𝛼 as a function of disc radius and height above midplane,Fig. 7, shows that the MRI-active zone (i.e., where 𝛼 ą 𝛼 DZ ) isconfined to the vicinity of the midplane in our models. This is qual-itatively different from the vertically isothermal case, where theMRI-active zone occurs around the midplane in the innermost re-gion, but rises into the upper layers at larger radii (see Mohanty et al.2018). The latter configuration, where the active zone is sandwichedbetween a dead zone around the midplane below and an ambipolar-dead zone above, emerges as a consequence of the ionization in-creasing with height above the disc midplane in the isothermal case,as discussed above. Moreover, comparing our non-irradiated discs(first two panels of Fig. 7) with the irradiated one (last panel), wesee that heating by stellar irradiation makes little difference to theextent of the MRI-active zone, since it only weakly affects the mid- MNRAS , 1–21 (2020) Radius [AU]10 H e i g h t [ A U ] T e m p e r a t u r e [ K ] Radius [AU]10 H e i g h t [ A U ] T e m p e r a t u r e [ K ] Radius [AU]10 H e i g h t [ A U ] T e m p e r a t u r e [ K ] Figure 3.
Temperature as a function of location in the disc for a viscously-heated constant-opacity model (left), viscously-heated model with realistic opacities(middle), and a viscously and irradiation-heated model with realistic opacities (right). In each panel the solid line shows the disc photosphere ( 𝜏 R “ {
3) andthe dashed line shows the pressure scale height ( 𝑃 “ 𝑒 ´ { 𝑃 mid ). The dotted line in the right-hand-side panel shows the surface at which 𝜏 irr “ {
3. Notethat the inclusion of heating by stellar irradiation does not strongly affect the disc midplane temperature. See Section 3.1.1. Radius [AU]10 H e i g h t [ A U ] radiativeconvective 10 Radius [AU]10 H e i g h t [ A U ] Radius [AU]10 H e i g h t [ A U ] Figure 4.
Radiative and convective zones (light and dark red, respectively) for a viscously-heated constant-opacity model (left), viscously-heated model withrealistic opacities (middle), and a viscously and irradiation-heated model with realistic opacities (right). In each panel the dashed line shows the pressure scaleheight ( 𝑃 “ 𝑒 ´ { 𝑃 mid ). In all three models the disc is convectively unstable within few scale heights. See Section 3.1.1. plane temperature, and the hot uppermost layers are dead due toambipolar diffusion.Finally, in Fig. 8 we compare radial profiles of the vertically-averaged viscosity parameter ¯ 𝛼 and the MRI-generated magneticfield strength 𝐵 from the four models (the vertically isothermaland constant opacity disc of Mohanty et al. 2018, plus our threenon-isothermal discs).In line with theoretical expectations and the results of Mohantyet al. (2018), the vertically-averaged viscosity parameter ¯ 𝛼 decreasesas a function of orbital radius. At some distance from the star theMRI is completely quenched and the average viscosity parameterreaches the minimum value, ¯ 𝛼 “ 𝛼 DZ . Interestingly, the ¯ 𝛼 radialprofile is both qualitatively and quantitatively similar in all fourmodels. However, the radial profile of the magnetic field strength 𝐵 reveals qualitative differences. In the vertically-isothermal model, 𝐵 p 𝑟 q features a sharp drop at „ „ All four steady-state models discussed above feature a maximum inboth the local gas pressure and the surface density (Fig. 9), at theradial location where the vertically-averaged viscosity parameter ¯ 𝛼 falls to its minimum value (see Fig. 8). Inwards of this location, theMRI-driven accretion efficiency (i.e., ¯ 𝛼 ) increases; for a radially-constant (i.e., steady-state) gas accretion rate, this leads to a decreasein both the surface density and the midplane pressure. Since the ra-dial profile of ¯ 𝛼 is similar in the four models, their surface densityand midplane pressure profiles are similar too. Within our threeself-consistent models, the addition of heating by stellar irradiationmoves the pressure maximum outwards by „ „ In this section we have considered how the details of the disc thermalstructure model affect the disc’s ionization state and where in thedisc the MRI is active. We have also examined the importance of thevarious physical effects on the location of the MRI-induced pressuremaximum in the inner disc.
MNRAS , 1–21 (2020)
RI-active inner disc F i rr / F a cc Radius [AU]0.0250.0500.0750.1000.1250.150 I n c i d e n t a n g l e () [ r a d ] Figure 5.
Top:
Ratio of the total irradiation heating to the total viscous dis-sipation ( 𝐹 irr { 𝐹 acc ) as a function of radius. The total (vertically-integrated)absorbed stellar flux 𝐹 irr is at least an order of magnitude larger than thetotal viscous dissipation 𝐹 acc at any given radius. Bottom:
Grazing angle ofstellar irradiation ( 𝜙 ) as a function of radius. Note that the 𝐹 irr { 𝐹 acc curvefollows the trend in 𝜙 . In both panels, the gray portion of the curve indicatesthe region affected by the inner boundary condition on 𝜙 . See Section 3.1.1. Our results show that when the disc vertical structure and theMRI-driven viscosity are considered self-consistently, the MRI isactive around the disc midplane. In previous work using a simplifiedvertically-isothermal model, a different configuration arises in thevicinity of the pressure maximum, where the MRI is active in a layerabove (and below) the dead midplane. Though both approachesyield a similar location for the pressure maximum for our fiducialparameters, the difference in the physical behaviour of the disccautions against the use of the simplified model.Furthermore, for the fiducial parameters considered here, wefind that the midplane temperature is primarily determined by vis-cous dissipation, and not stellar irradiation. As a result, the locationof the pressure maximum, where the midplane temperature and theionization levels fall below a critical value, is negligibly affected bystellar irradiation. Finally, we find that the inner disc is convectivelyunstable, even if the disc opacity is constant. Hence, even thoughthe location of the pressure maximum is similar to that in previouswork, the disc’s structure is qualitatively different. These findingsare further discussed in sections 4.2 and 4.3.
In this section we build upon our models from the previous section.We first explore the effects of dust on the disc’s ionization state and on the MRI. We then construct our complete model of the inner discby also considering non-thermal sources of ionization.
Here we consider two models. The first is our most complex modelfrom the previous section, which incorporates a self-consistent ver-tical structure and realistic opacities, and includes stellar irradiation,but accounts for only gas-phase thermal ionization and recombina-tion of potassium (solid lines in Figs. 10 and 11). The second is amodel including all these processes, as well as dust effects on theionization (dashed lines in Figs. 10 and 11); the latter effects in-clude adsorption of neutral atoms and free charges onto dust grains,recombinations on the grain surfaces, and thermionic and ion emis-sion from the grains.Figure 10 depicts the radial profiles of the vertically-averagedviscosity ¯ 𝛼 , magnetic field strength, midplane temperature and themidplane fractional ionization ( 𝑛 e { 𝑛 H ) for the two models. Wesee that the radial profile of ¯ 𝛼 , and the radial extent of the MRI-active zone (i.e., where ¯ 𝛼 ą 𝛼 DZ ), are similar in the two models.Consequently, their surface density and midplane pressure profiles,plotted in Fig. 11, are also similar. Interestingly, Fig. 10 also showsthat when the influence of dust on the ionization is included, both¯ 𝛼 and the ionization fraction are somewhat higher at a given radiuswithin the active region, while the midplane temperature is lower.We discuss these effects further in section 4.1.Figure 12 illustrates the vertical structure of the model withdust at three different disc radii. The top row shows the local vis-cosity parameter 𝛼 , the middle row shows the number densities offree electrons and ions, and the bottom row shows the various con-tributions to the free electron production rate per unit volume. Asin the thermally ionized disc (Fig. 6), the MRI is active in the hotdisc midplane. The ionization levels decrease with height above themidplane as the temperature decreases, and increase again in thedisc atmosphere heated by stellar irradiation. In the hot upper layersambipolar diffusion quenches the MRI.The plots of free electron and ion production rates show thatthermionic and ion emission are the dominant ionization sources.Thus, Fig. 10 is misleading in the sense that, while the differencesin the global structure are very small when dust effects are added, it is not because dust effects are minor . Clearly, dust dominates thechemistry in the inner disc. Rather, for the parameters assumed here,the ionization levels as a function of temperature and density, whendust is included, are similar to the levels obtained from gas-phasethermal ionization only, due to the similarity between the ionizationpotential of potassium and the grain work function . In this section we present our full model of the inner disc which,in addition to the above dust effects, also includes non-thermalsources of ionization: stellar X-rays, cosmic rays and radionuclidesall ionize H , ultimately producing metallic ions and free electrons.The resulting radial profiles of ¯ 𝛼 , magnetic field strength, midplanetemperature, density and ionization levels are shown in Fig. 10 asdotted lines. The radial profile of ¯ 𝛼 seems to overlap with the modelwith only dust effects (dashed line). However, in the region wherethe MRI is dead in the dust-only model, ¯ 𝛼 is slightly higher than 𝛼 DZ in the new model including H ionization. The magnetic fieldstrength also does not drop to zero in this model, further revealingthat the MRI remains active here. MNRAS000
Here we consider two models. The first is our most complex modelfrom the previous section, which incorporates a self-consistent ver-tical structure and realistic opacities, and includes stellar irradiation,but accounts for only gas-phase thermal ionization and recombina-tion of potassium (solid lines in Figs. 10 and 11). The second is amodel including all these processes, as well as dust effects on theionization (dashed lines in Figs. 10 and 11); the latter effects in-clude adsorption of neutral atoms and free charges onto dust grains,recombinations on the grain surfaces, and thermionic and ion emis-sion from the grains.Figure 10 depicts the radial profiles of the vertically-averagedviscosity ¯ 𝛼 , magnetic field strength, midplane temperature and themidplane fractional ionization ( 𝑛 e { 𝑛 H ) for the two models. Wesee that the radial profile of ¯ 𝛼 , and the radial extent of the MRI-active zone (i.e., where ¯ 𝛼 ą 𝛼 DZ ), are similar in the two models.Consequently, their surface density and midplane pressure profiles,plotted in Fig. 11, are also similar. Interestingly, Fig. 10 also showsthat when the influence of dust on the ionization is included, both¯ 𝛼 and the ionization fraction are somewhat higher at a given radiuswithin the active region, while the midplane temperature is lower.We discuss these effects further in section 4.1.Figure 12 illustrates the vertical structure of the model withdust at three different disc radii. The top row shows the local vis-cosity parameter 𝛼 , the middle row shows the number densities offree electrons and ions, and the bottom row shows the various con-tributions to the free electron production rate per unit volume. Asin the thermally ionized disc (Fig. 6), the MRI is active in the hotdisc midplane. The ionization levels decrease with height above themidplane as the temperature decreases, and increase again in thedisc atmosphere heated by stellar irradiation. In the hot upper layersambipolar diffusion quenches the MRI.The plots of free electron and ion production rates show thatthermionic and ion emission are the dominant ionization sources.Thus, Fig. 10 is misleading in the sense that, while the differencesin the global structure are very small when dust effects are added, it is not because dust effects are minor . Clearly, dust dominates thechemistry in the inner disc. Rather, for the parameters assumed here,the ionization levels as a function of temperature and density, whendust is included, are similar to the levels obtained from gas-phasethermal ionization only, due to the similarity between the ionizationpotential of potassium and the grain work function . In this section we present our full model of the inner disc which,in addition to the above dust effects, also includes non-thermalsources of ionization: stellar X-rays, cosmic rays and radionuclidesall ionize H , ultimately producing metallic ions and free electrons.The resulting radial profiles of ¯ 𝛼 , magnetic field strength, midplanetemperature, density and ionization levels are shown in Fig. 10 asdotted lines. The radial profile of ¯ 𝛼 seems to overlap with the modelwith only dust effects (dashed line). However, in the region wherethe MRI is dead in the dust-only model, ¯ 𝛼 is slightly higher than 𝛼 DZ in the new model including H ionization. The magnetic fieldstrength also does not drop to zero in this model, further revealingthat the MRI remains active here. MNRAS000 , 1–21 (2020) L o c a l v i s c o u s no irradiationwith irradiation10 n e / n H no irradiationwith irradiation10 R e s i s t i v i t i e s [ c m / s ] Height [AU] 10 , no irr. A , no irr. O , w. irr. A , w. irr. Figure 6.
Local viscous parameter 𝛼 (top), ionization fraction 𝑛 e { 𝑛 H (middle) and magnetic resistivities (bottom) at three different disc radii (as indicatedin the panel titles), in two different models: one with realistic opacities but no irradiation (solid lines), and one that also includes irradiation (dashed lines).Despite the high ionization fraction in the irradiation-heated disc upper layers, 𝜂 𝐴 is still very high in these layers, i.e., ambipolar diffusion quenches the MRIthere. See Section 3.1.2. Radius [AU]10 H e i g h t [ A U ] L o c a l v i s c o u s Radius [AU]10 H e i g h t [ A U ] L o c a l v i s c o u s Radius [AU]10 H e i g h t [ A U ] L o c a l v i s c o u s Figure 7.
Local viscous parameter 𝛼 as a function of location in the disc for a viscously-heated constant-opacity model (left), viscously-heated model withrealistic opacities (middle), and a viscously and irradiation-heated model with realistic opacities (right). In each panel the solid line shows the disc photosphere( 𝜏 R “ {
3) and the dashed line shows the pressure scale height ( 𝑃 “ 𝑒 ´ { 𝑃 mid ). The dotted line in the right-hand-side panel shows the surface at which 𝜏 irr “ {
3. Note that the heating by stellar irradiation has a very weak effect on the extent of the MRI-active region. See Section 3.1.3.MNRAS , 1–21 (2020)
RI-active inner disc V e r t i c a ll y - a v e r a g e d v i s c o u s Radius [AU]0.50.00.51.01.52.0 l o g B [ G a u ss ] real. opac., acc + irr real. opac., acc const. opac., acc const. opac., acc , vert. isoth. Figure 8.
Vertically-averaged viscosity parameter ( ¯ 𝛼 ) and the MRI-generated magnetic field strength ( 𝐵 ) as functions of radius, in a vertically-isothermal model as well as in self-consistent models of varying complexity(viscously heated constant-opacity model, viscously heated model with re-alistic opacities, and viscously and irradiation-heated model with realisticopacities), as indicated in plot legend. The vertically-averaged viscosity pa-rameter profile is similar in all four models, but the magnetic field strengthprofile is qualitatively different in the vertically-isothermal model (in whichthe field strength decreases as a function of radius outwards from „ Figure 13 shows the vertical structure of this disc at three differ-ent radii. It reveals that, at large radii where the MRI was quenchedin previous models, the MRI now remains active at large heights(see also Fig. 14). The situation is reminiscent of the appearanceof an MRI-active layer near the disc surface in the thermally ion-ized, vertically isothermal model of Mohanty et al. (2018), but thephysical reason is very different: in the latter model, it is due to theisothermal assumption, which leads to increasing fractional ioniza-tion with height; in our present non-isothermal model, it is due toadditional, non-thermal sources of ionization (mainly X-rays; seebelow) which elevate the fractional ionization near the disc surface.The plot of ionization levels in Fig. 13 shows that the fractionalabundances of free electrons and metal ions increase strongly to-wards the disc surface, where non-thermal ionization dominates. Asa result, potassium ions are depleted in the upper layers via recom-bination with the abundant electrons. Furthermore, the values of themetal ion abundance near the surface are themselves noteworthy:they greatly exceed the solar abundance of Mg, and even of C andO. The reason is as follows. In our chemical network, the ionizationof an H molecule by a non-thermal source effectively produces afree electron and a metal ion, the latter through (implicitly included) S u r f a c e d e n s i t y [ g c m ] Radius [AU]10 M i dp l a n e p r e ss u r e [ d y n c m ] real. opac., acc + irr real. opac., acc const. opac., acc const. opac., acc , vert. isoth. Figure 9.
Surface density (top) and midplane pressure (bottom) as functionsof radius in a vertically-isothermal and self-consistent models of varyingcomplexity (same as in Fig. 8). In each model the density and the midplanepressure maximum correspond to the point at which the viscous ¯ 𝛼 reachesa minimum value. See Section 3.1.3. rapid charge exchange between a metal atom and the H ion. This isvalid as long as the number of H ions remains lower than the totalnumber of metal atoms; in this case, the precise total abundance ofmetals is unimportant, and is not accounted for in our calculations.When the H ion abundance exceeds that of metal atoms, however,our network fails: it yields a spuriously high metal ion abundance,when in reality H ions dominate (since there are no remainingmetal atoms to transfer their charge to). Additionally, when the H ion abundance becomes very high (e.g., when 𝑛 H ` { 𝑛 H exceeds „ ´ ), our assumption that H is effectively a neutral species dueto charge transfer also breaks down.Clearly, above this level our chemical model is not applicable,as the ionized hydrogen would become an important species andour assumption of a constant hydrogen number density would beinvalid. Nevertheless, this is above the MRI active zone at all orbitalradii (see the black solid line in Fig. 14, and also the blue solid linein Fig. 13), and is thus not germane to our conclusions.Figure 15 compares the contributions from stellar X-rays andcosmic rays to the hydrogen ionization rate. As expected, the un-attenuated ionization rate due to X-rays is higher at the disc surface,but cosmic rays can penetrate deeper in the disc. Nevertheless, theMRI-active region in the upper disc layers is mostly ionized by X-rays, in agreement with previous work (e.g. Glassgold et al. 1997;Ercolano & Glassgold 2013).Finally, note that in this outer region, where the MRI is ac- MNRAS000
Surface density (top) and midplane pressure (bottom) as functionsof radius in a vertically-isothermal and self-consistent models of varyingcomplexity (same as in Fig. 8). In each model the density and the midplanepressure maximum correspond to the point at which the viscous ¯ 𝛼 reachesa minimum value. See Section 3.1.3. rapid charge exchange between a metal atom and the H ion. This isvalid as long as the number of H ions remains lower than the totalnumber of metal atoms; in this case, the precise total abundance ofmetals is unimportant, and is not accounted for in our calculations.When the H ion abundance exceeds that of metal atoms, however,our network fails: it yields a spuriously high metal ion abundance,when in reality H ions dominate (since there are no remainingmetal atoms to transfer their charge to). Additionally, when the H ion abundance becomes very high (e.g., when 𝑛 H ` { 𝑛 H exceeds „ ´ ), our assumption that H is effectively a neutral species dueto charge transfer also breaks down.Clearly, above this level our chemical model is not applicable,as the ionized hydrogen would become an important species andour assumption of a constant hydrogen number density would beinvalid. Nevertheless, this is above the MRI active zone at all orbitalradii (see the black solid line in Fig. 14, and also the blue solid linein Fig. 13), and is thus not germane to our conclusions.Figure 15 compares the contributions from stellar X-rays andcosmic rays to the hydrogen ionization rate. As expected, the un-attenuated ionization rate due to X-rays is higher at the disc surface,but cosmic rays can penetrate deeper in the disc. Nevertheless, theMRI-active region in the upper disc layers is mostly ionized by X-rays, in agreement with previous work (e.g. Glassgold et al. 1997;Ercolano & Glassgold 2013).Finally, note that in this outer region, where the MRI is ac- MNRAS000 , 1–21 (2020) V e r t i c a ll y - a v e r a g e d v i s c o u s thermalthermal+dustthermal+dust+H ionization0.50.00.51.01.52.0 l o g B [ G a u ss ] M i dp l a n e t e m p e r a t u r e [ K ] Radius [AU]10 M i dp l a n e n e / n H Figure 10.
Vertically-averaged viscous parameter ( ¯ 𝛼 ), MRI-generated mag-netic field strength ( 𝐵 ), midplane temperature and the midplane free elec-tron fraction ( 𝑛 e { 𝑛 H ) as functions of radius, for a model with thermalionization only, a model with thermal ionization and dust effects, and amodel which also includes realistic non-thermal sources of ionization of H ( 𝜁 “ 𝜁 R ` 𝜁 CR ` 𝜁 X ). The disc structure is quantitatively similar in allthree models; however, the main sources of ionization in the models withdust are thermionic and ion emission. See Section 3.2.1. S u r f a c e d e n s i t y [ g c m ] thermalthermal+dustthermal+dust+H ioniz.10 Radius [AU]10 M i dp l a n e p r e ss u r e [ d y n c m ] Figure 11.
Surface density (top) and midplane pressure (bottom) as func-tions of radius for a model with thermal ionization only, a model with thermalionization and dust effects, and a model which also includes realistic non-thermal sources of ionization of H ( 𝜁 “ 𝜁 R ` 𝜁 CR ` 𝜁 X ). The disc structureis similar in all three. See Section 3.2.1. tive in the upper disc, the gas accretes primarily through the deadzone (driven by the dead-zone viscosity 𝛼 DZ ), since the density atthe dead disc midplane is much higher than in the X-ray-ionizedMRI-active layer. This is why the vertically-averaged viscosity pa-rameter is close to the dead-zone value, ¯ 𝛼 „ 𝛼 DZ , outwards fromthe pressure maximum. In the results presented so far, for our fiducial choice of disc and dustparameters, solutions for the vertical disc structure (i.e., solutionsin disc height 𝑧 surf ) appear to be unique. In general, there is also asingle peak in the vertically-averaged viscosity ¯ 𝛼 as a function ofmagnetic field strength 𝐵 (which determines our choice for 𝐵 ; seesection 2.5.2). The exception is the vicinity of the orbital radius atwhich the MRI is quenched at the disc midplane (see Fig. 14). There,¯ 𝛼 p 𝐵 q has two peaks, one corresponding to the solution where theMRI is active at disc midplane, and the other to the solution wherethe MRI is active in the upper disc layers, mostly ionized by stellarX-rays. As discussed in section 2.5, we choose 𝐵 such that ¯ 𝛼 ismaximized in this case as well.Note that there could be, in principle, an MRI-active layerhigh up in the disc at shorter orbital radii as well, in addition tothe active layer at midplane. Here, this does not appear due to ourassumption that the magnetic field strength 𝐵 is vertically constant.At the high 𝐵 necessary to drive efficient accretion at midplane, MNRAS , 1–21 (2020)
RI-active inner disc L o c a l v i s c o u s n / n H n e / n H n K + / n H I o n i z a t i o n r a t e [ c m s ] Height [AU] 10 thermionicion emissionthermal Figure 12.
Local viscous parameter 𝛼 (top), fractional abundance of charged species 𝑛 { 𝑛 H (middle) and ionization rates (bottom; thermionic R therm , ionemission R K , evap 𝑓 ` , thermal 𝑘 𝑛 H 𝑛 K ) at three different radii (as indicated in panel titles) for the model with thermal ionization and dust effects. Thermionicemission is the primary source of free electrons in the MRI-active regions. See Section 3.2.1. the MRI is quenched in the low-density disc atmosphere due toambipolar diffusion (since the ambipolar criterion for active MRIalso encapsulates the requirement that the magnetic pressure shouldbe less than the thermal pressure). Therefore, it is only when thetemperature drops and high-temperature ionization effects (thermalionisation of potassium, thermionic and ion emission from grains)can no longer drive the accretion efficiently in the midplane that ourmodel features the active layer (generated by X-ray ionisation) highabove the disc midplane.For a different choice of parameters, e.g. if the maximum dustgrain size is 𝑎 max “ 𝜇 m, there may exist an additional rangeof orbital radii where there are multiple peaks in ¯ 𝛼 p 𝐵 q and alsomultiple solutions for the disc vertical structure (for 𝑧 surf ) at a fixedvalue of the magnetic field strength 𝐵 . Similarly to the case above,this issue arises due to the competing effects of high-temperaturesources of ionization and X-rays. As these sources of free electronshave different dependencies on the disc structure (density, tempera-ture, column density), their combination leads to non-monotonousvariations in the electron number density as a function of heightabove the disc’s midplane. Since the viscous dissipation due to theMRI is sensitive to the ionization fraction, the total dissipation canbe a non-monotonous function of 𝑧 surf . Since the solution for the vertical disc structure is determined by an equilibrium between aninput and an output total heat, this can lead to multiple solutions in 𝑧 surf .To illustrate this issue, we show in Fig. 16 an example ofthree thermally-stable solutions for the vertical disc structure at afixed value of magnetic field strength 𝐵 that appears in our modelfor a maximum grain size 𝑎 max “ 𝜇 m (we ignore thermallyunstable solutions). Note that the dependence of the overall discstructure and the location of the pressure maximum on the dust grainsize is presented and thoroughly discussed in our companion paper(Jankovic et al. in prep. ). Here we only discuss how we deal withthe multiple solutions. Fig. 16 shows the viscosity 𝛼 as a functionof height in the top panel and the ratio 𝑛 e { 𝑛 H in the bottom.Evidently, small variations in the free electron number densitycorrespond to large variations in the viscosity 𝛼 , all at heights belowone disc pressure scale height (indicated by gray lines). This impliesthat the difference between these solutions is likely unphysical fortwo reasons. First, the viscosity 𝛼 is in reality driven by turbulence,and turbulent motions should not abruptly change over length-scalesmuch smaller than a single pressure scale height. Second, chemicalspecies can also be expected to be spatially mixed by turbulence, andso such vertical variations in the ionization fraction as we obtain here MNRAS000
Local viscous parameter 𝛼 (top), fractional abundance of charged species 𝑛 { 𝑛 H (middle) and ionization rates (bottom; thermionic R therm , ionemission R K , evap 𝑓 ` , thermal 𝑘 𝑛 H 𝑛 K ) at three different radii (as indicated in panel titles) for the model with thermal ionization and dust effects. Thermionicemission is the primary source of free electrons in the MRI-active regions. See Section 3.2.1. the MRI is quenched in the low-density disc atmosphere due toambipolar diffusion (since the ambipolar criterion for active MRIalso encapsulates the requirement that the magnetic pressure shouldbe less than the thermal pressure). Therefore, it is only when thetemperature drops and high-temperature ionization effects (thermalionisation of potassium, thermionic and ion emission from grains)can no longer drive the accretion efficiently in the midplane that ourmodel features the active layer (generated by X-ray ionisation) highabove the disc midplane.For a different choice of parameters, e.g. if the maximum dustgrain size is 𝑎 max “ 𝜇 m, there may exist an additional rangeof orbital radii where there are multiple peaks in ¯ 𝛼 p 𝐵 q and alsomultiple solutions for the disc vertical structure (for 𝑧 surf ) at a fixedvalue of the magnetic field strength 𝐵 . Similarly to the case above,this issue arises due to the competing effects of high-temperaturesources of ionization and X-rays. As these sources of free electronshave different dependencies on the disc structure (density, tempera-ture, column density), their combination leads to non-monotonousvariations in the electron number density as a function of heightabove the disc’s midplane. Since the viscous dissipation due to theMRI is sensitive to the ionization fraction, the total dissipation canbe a non-monotonous function of 𝑧 surf . Since the solution for the vertical disc structure is determined by an equilibrium between aninput and an output total heat, this can lead to multiple solutions in 𝑧 surf .To illustrate this issue, we show in Fig. 16 an example ofthree thermally-stable solutions for the vertical disc structure at afixed value of magnetic field strength 𝐵 that appears in our modelfor a maximum grain size 𝑎 max “ 𝜇 m (we ignore thermallyunstable solutions). Note that the dependence of the overall discstructure and the location of the pressure maximum on the dust grainsize is presented and thoroughly discussed in our companion paper(Jankovic et al. in prep. ). Here we only discuss how we deal withthe multiple solutions. Fig. 16 shows the viscosity 𝛼 as a functionof height in the top panel and the ratio 𝑛 e { 𝑛 H in the bottom.Evidently, small variations in the free electron number densitycorrespond to large variations in the viscosity 𝛼 , all at heights belowone disc pressure scale height (indicated by gray lines). This impliesthat the difference between these solutions is likely unphysical fortwo reasons. First, the viscosity 𝛼 is in reality driven by turbulence,and turbulent motions should not abruptly change over length-scalesmuch smaller than a single pressure scale height. Second, chemicalspecies can also be expected to be spatially mixed by turbulence, andso such vertical variations in the ionization fraction as we obtain here MNRAS000 , 1–21 (2020) L o c a l v i s c o u s n / n H n e / n H n K + / n H n i / n H I o n i z a t i o n r a t e [ c m s ] Height [AU] 10 thermionicion emissionthermalionization of H Figure 13.
Local viscous parameter 𝛼 (top), ionization fraction (middle; electrons: 𝑛 e { 𝑛 H , potassium ions: 𝑛 K ` { 𝑛 H , metal ions: 𝑛 i { 𝑛 H ) and ionizationrates (bottom; thermionic R therm , ion emission R K , evap 𝑓 ` , thermal 𝑘 𝑛 H 𝑛 K and non-thermal 𝜁 𝑛 H ) at three different radii (as indicated in panel titles) forthe model with all sources of thermal and non-thermal ionization. Non-thermal ionization produces an MRI-active region high above disc midplane at largerradii (see the top right panel). The blue line in the right panel on the second row indicates the upper boundary of this MRI-active region at the given radius,and the ionization fraction at which it occurs. Note that the metal ion fraction becomes unrealistically large only to the right of the blue line, i.e., only abovethe high-altitude MRI-active layer. See Section 3.2.2. might be smoothed over in reality. Since resolving these issues isbeyond the scope of our models, we simply always choose a solutionwith minimum 𝑧 surf , which also appears to always correspond to amaximum ¯ 𝛼 at the given magnetic field strength. As previously suggested by Desch & Turner (2015), the inner discis primarily ionized through thermionic and potassium ion emis-sion from dust grains. These processes counteract adsorption offree charges onto dust grains at the high temperatures present in theinner disc. We show that, for our fiducial parameters, the introduc-tion of dust effects on disc ionization has very little effect on thelocation of the pressure maximum. This is because thermionic andion emission become efficient above roughly the same thresholdtemperature as thermal ionization of potassium. Additionally, weshow that non-thermal sources of ionization are unimportant for thefiducial parameters considered here. Thus, similar to our results inSection 3.1, while the inclusion of new physics in our model has not fundamentally changed the position of the pressure maximum,the physics setting its position is very different from that describedin Mohanty et al. (2018).
A key feature of the inner disc is that the MRI drives high viscosityin the innermost regions, close to the star, but becomes largelysuppressed at larger orbital distances (in the so-called dead zone).This leads to the formation of a local gas pressure maximum thatmay play a key role in planet formation at short orbital distances(Chatterjee & Tan 2014). This decrease in viscosity is expectedto arise because the innermost regions are hot enough (>1000 K)to thermally ionize potassium (coupling the gas to the magneticfield), but further out temperature and ionization levels decreasesubstantially (Gammie 1996). In a previous study, we showed thatin a thermally-ionized disc coupled self-consistently to an MRI
MNRAS , 1–21 (2020)
RI-active inner disc Radius [AU]10 H e i g h t [ A U ] L o c a l v i s c o u s Figure 14.
Local viscous parameter 𝛼 as a function of location in the discfor the model with all sources of thermal and non-thermal ionization. In theinnermost disc, thermionic and ion emission ionize the dense regions aroundthe disc midplane, producing the high MRI-driven 𝛼 there. At larger radii,the MRI is active in a low-density layer high above the disc midplane, dom-inated by non-thermal sources of ionization. The solid black line indicatesthe surface in the disc above which the ionization fraction 𝑛 e { 𝑛 H ą ´ ,i.e., above which the assumptions of our simple chemical network break;this surface is above the MRI-active region at all radii. See Section 3.2.2. viscosity, the inner edge of the dead zone lies at a few tenths of anAU (Mohanty et al. 2018).One of the key differences between this work and that of Mo-hanty et al. (2018) is that here we also take into account the effectsof dust on the disc’s ionization state. Small dust grains present inthe disc sweep up free electrons and ions from the gas, and theserecombine quickly on the grain surfaces. In the bulk of the proto-planetary disc dust grains therefore efficiently lower the ionizationfraction, decoupling the magnetic field from the gas and suppress-ing the MRI (Sano et al. 2000; Ilgner & Nelson 2006; Wardle 2007;Salmeron & Wardle 2008; Bai & Goodman 2009; Mohanty et al.2013). However, in the inner regions of protoplanetary discs dustgrains also act to increase the ionization levels, as at high temper-atures they can also emit electrons and ions into the gas (Desch &Turner 2015). The balance between thermal ionization and theseprocesses then determines how well ionized the inner disc is, andthus the extent of the high-viscosity region and the location of thedead zone inner edge.The top panel of Fig. 10 shows that addition of dust only weaklyaffects the vertically-averaged viscosity ¯ 𝛼 in the inner disc. At agiven orbital radius, ¯ 𝛼 is even slightly higher than in the model withno dust, implying that thermionic and ion emission are importantsources of ionization. In fact, as Desch & Turner (2015) showed,thermionic and ion emission can become the main source of freeelectrons at high temperatures. For our disc model this can be seenin the bottom panel of Fig. 13, which shows that at the hot discmidplane thermionic and ion emission dominate over other sourcesof ionization. Clearly then, for the chosen parameters, the adsorptionof charges onto grains is more than offset by the expulsion of chargesfrom hot grain surfaces.Similarity in the resulting disc structure in the models withand without dust grains can be explained by the similar dependencyon temperature that thermal and thermionic/ion emission have (andwhich can also be deduced from the bottom panel of Fig. 13). As dis-cussed by Desch & Turner (2015), thermal ionization of potassium becomes efficient at temperatures above „ “ .
34 eV. The temperature at whichthermionic and ion emission become important is determined by thework function 𝑊 of the material out which dust grains are made,and for silicates 𝑊 „ „ 𝑊 „ IP; see eq. (21)).Dust grains then become negatively charged, which reduces theeffective potential that electrons need to overcome for thermionicemission (the effective work function 𝑊 eff ).The above results could change significantly as dust grainsgrow or as they accumulate in the inner disc (as needed for theformation of solid planet cores). Dust adsorption of free charges,for example, becomes much less efficient for larger grains, since thetotal grain surface area decreases (Sano et al. 2000; Ilgner & Nelson2006). We investigate how dust growth and varying dust-to-gas ratioaffect the inner disc structure in our companion paper (Jankovic etal. in prep. ). We find that the absorbed flux of stellar irradiation is many timeshigher than the heat flux generated by accretion at any given radius,yet irradiation has a very small effect on the disc midplane tem-perature. Consequently, the ionization levels and the MRI-drivenviscosity are similar in the models with and without stellar irradia-tion. Why does the irradiation have such a marginal effect?Essentially, it is because the stellar irradiation heats the disc’soptically-thin regions, from which heat escapes easily. Accretionheat is generated deep in the disc, where the optical depth is muchhigher. In the absence of stellar irradiation, the midplane temper-ature in the optically thick disc is 𝜎 SB 𝑇 „ 𝐹 acc 𝜏 mid (Hubeny1990, though in our work the midplane temperature is somewhatlower due to convection). If on top of a viscously-heated layerof optical thickness 𝜏 mid " 𝜏 upper , it follows from eq. (8) and (9) that 𝜎 SB 𝑇 „ 𝐹 acc 𝜏 mid ` 𝐹 irr 𝜏 upper (again, neglecting convection).Here 𝜏 upper is the optical depth of the disc to its own radiation downto a height at which the disc becomes optically thick to stellar irra-diation. Then, if 𝜏 mid is sufficiently larger than 𝜏 upper , the midplanetemperature is determined by viscous dissipation.Our results are consistent with those of D’Alessio et al. (1998),who also found that models with and without stellar irradiation yieldroughly the same midplane temperatures in the optically-thick innerdisc. Similarly, Flock et al. (2019) considered 2D static radiation-hydrodynamics models of the inner disc heated by stellar irradiationonly, but found that the midplane temperature (and the orbital radiusof the dead zone inner edge) would increase appreciably if accre-tion heat were included, on the condition that the accretion heat isreleased near the optically-thick midplane.Note, additionally, that we have not considered the details ofthe inner disc edge or the dust sublimation line. If the inner rim ofthe disc is puffed-up, it would throw a shadow over a portion ofthe inner disc (Natta et al. 2001; Dullemond et al. 2001), furtherreducing the importance of stellar irradiation. MNRAS , 1–21 (2020) X , C R , R [ s ] Height [AU] XCRR
Figure 15.
Ionization rates of molecular hydrogen due to stellar X-rays ( 𝜁 X ), cosmic rays ( 𝜁 CR ), and radionuclides ( 𝜁 R ) for the model with all sources ofthermal and non-thermal ionization. Stellar X-rays are the dominant source of ionization in the disc upper layers. See Section 3.2.2. L o c a l v i s c o u s Height [AU]10 n e / n H Figure 16.
Degeneracy in the vertical disc structure at a fixed magnetic fieldstrength in a model where the maximum dust grain size is 𝑎 max “ 𝜇 m:local viscous parameter 𝛼 as a function of height (top), and the fractionalelectron number density ( 𝑛 e { 𝑛 H ; bottom) as functions of height in threedifferent equilibrium solutions. The different solutions arise from verticalvariations in the viscous 𝛼 ; the lengthscales of these variations are muchsmaller than the disc pressure scale height (shown for each solution by thevertical gray lines). The solutions shown here are at a radius of 0.22 AU, forfield strength log 𝐵 𝑧 “ .
226 and grazing angle 𝜙 “ .
06. The vertically-averaged viscous ¯ 𝛼 for the solid, dashed and dotted lines are, respectively,2 ˆ ´ , 1 . ˆ ´ , and 1 . ˆ ´ . See Section 3.2.3. In section 3.1 we showed that a large region of the inner disc isconvectively unstable. We find this to be the case even when theopacities are constant, i.e., a super-linear growth of the opacity withtemperature (Lin & Papaloizou 1980) is not needed. Here, the hightemperature gradient is established because the heat is depositeddeep within the optically-thick disc. In the presented models mostof the viscous dissipation happens near the midplane, where the MRIis active, but the same is also true for a vertically-constant viscosity 𝛼 . This result can also be confirmed analytically. We consider asimplified problem of radiative transfer in the optically thick limit,where the temperature is given by 𝜎𝑇 “ 𝜏𝐹 p 𝑧 surf q (Hubeny1990). Assuming a constant disc opacity, the equation of hydro-static equilibrium can be re-written as 𝑑𝑃𝑑𝜏 “ Ω 𝑧𝜅 R . The temperaturegradient is then given by ∇ “ 𝑑 ln 𝑇𝑑 ln 𝑃 “ 𝜅 R 𝑃 Ω 𝑧𝜏 . The appropriate upper boundary condition for this problem is thedisc photosphere ( 𝜏 “ { 𝑃 surf “ Ω 𝑧 surf 𝜏 surf { 𝜅 R (assuming that the disc is verticallyisothermal above the photosphere, Papaloizou & Terquem 1999).At the photosphere, given the chosen boundary condition, we have ∇ “ {
4. Near the midplane the optical depth is 𝜏 MID “ 𝜅 R Σ ,and we may estimate the midplane pressure as 𝑃 MID “ 𝜌 MID 𝑐 , MID “ Σ 𝐻 𝑐 , MID “ Ω Σ 𝐻, where the disc’s scale height 𝐻 is related to the midplane tempera-ture through hydrostatic equilibrium. Substituting 𝜏 MID and 𝑃 MID into the expression for the temperature gradient, we have ∇ MID “ 𝐻𝑧 , implying that such a disc should become convectively unstablea bit below one scale height. However, we can further estimatethe gradient ∇ near one scale height, by assuming that there 𝑃 𝐻 „ 𝑃 MID 𝑒 ´ { „ . 𝑃 MID , and 𝜏 𝐻 „ . 𝜏 MID , as followsfrom vertically isothermal, Gaussian profiles of pressure and den-sity. Thus, near 𝑧 “ 𝐻 , we have ∇ 𝐻 “ 𝐻𝑧 , showing that the disc should become convectively unstable aboveone scale height.
MNRAS , 1–21 (2020)
RI-active inner disc Furthermore, in the convectively unstable regions we use asimple approximation that convection is efficient and the tempera-ture gradient is isentropic. More detailed calculations would yieldan answer in which the temperature gradient lies between the isen-tropic one and the gradient given by the radiative transport. As itturns out, the result would not differ much from what is obtainedhere. In the optically-thick limit considered above, the difference inthe temperature profile when the entire flux is transported by radia-tion and when the entire flux is transported by convection is a veryweak function of optical depth, and remains small for rather largeoptical depths (Cassen 1993).Limiting the temperature gradient is the only role of convec-tion in our simple viscous model. In real discs, convection mightinteract with the MHD turbulence induced by the MRI. For exam-ple, in simulations of ideal MHD, Bodo et al. (2013) and Hiroseet al. (2014) find that convection can increase the angular momen-tum transport driven by the MRI by increasing the magnetic fieldstrength, although it appears that this is only the case when convec-tion is particularly strong (Hirose 2015). Concurrently, it is foundthat the relationship between the induced stress and the magneticfield strength are not modified. The consequences for the non-idealMHD regime, relevant in protoplanetary discs, are not clear. In oursolutions the value of the vertically-averaged MRI-driven viscositywould decrease with both a decrease and an increase in the magneticfield strength due to non-ideal effects, as discussed in section 2.5.2.Convection itself is not expected to drive the angular momentumtransport at a level comparable to the MRI (e.g. Lesur & Ogilvie2010; Held & Latter 2018), and in any case it is not self-sustainable,i.e., it requires an additional source of heat near disc midplane toestablish the high temperature gradient.
In our model, we have not considered the possibility that the tur-bulent elements driving the angular momentum transport may alsotransport energy (Ruediger et al. 1988). Such turbulent energy trans-port flux would be analogous to convection, transporting energydown the entropy gradient, the difference being that turbulent ele-ments may persist at sub-adiabatic temperature gradients (as theyare driven by other instabilities) in which case they transport en-ergy from cooler to hotter regions (Balbus 2000). We do not expectincluding this mode of energy transport in our calculations to ap-preciably change any of our results. D’Alessio et al. (1998) foundthat turbulent energy transport accounts for less then 20 % of thetotal energy flux at any given orbital radius. In our work, the con-vectively stable upper layers of the disc are MRI-dead, and thusthe thermal diffusivity due to turbulence is likely very low. In theregions of the disc where the radiative flux alone would yield super-adiabatic temperature gradient, we already assume that convectionefficiently establishes the adiabat. The proportion of turbulent en-ergy flux could be higher than the convective energy flux in suchregions, but the temperature would not change significantly. Hence,we have not included this effect in our analysis.
The criterion for ambipolar diffusion to quench the MRI that weemploy is valid in the strong-coupling regime (Bai & Stone 2011).Strong coupling requires that ionization equilibrium be achieved ona timescale 𝑡 chem shorter than the dynamical timescale 𝑡 dyn “ 𝜋 { Ω .Previously, Mohanty et al. (2018) reported that this condition is not Radius [AU]10 H e i g h t [ A U ] t c h e m / t d y n Figure 17.
Ratio of the shortest recombination timescale ( 𝑡 chem ) to thedynamical timescale ( 𝑡 dyn ) as a function of the location in the disc. Thesolid line indicates where 𝑡 chem { 𝑡 dyn “
1. Everywhere below this line thechemical equilibrium timescale is shorter than the dynamical timescale,justifying our assumption of the strong-coupling regime. See Section 4.5. fulfilled in most of the inner disc, as slow radiative recombinationsmake the chemical equilibrium timescale long. However, even in theabsence of dust, three-body recombination (recombinations throughcollisions with the abundant molecular hydrogen) are much fasterthan radiative recombinations, and adsorption onto grains is evenfaster (Desch & Turner 2015).Since we solve directly for the equilibrium ionization state, wedo not have access to the timescale 𝑡 chem . However, we can estimateit as 𝑡 chem “ 𝑛 e { R , where R is the fastest of the above threerecombination rates (in general, but not always, that is adsorptiononto dust grains). Fig. 17 shows, for our fully self-consistent modelwith our full chemical network, that 𝑡 chem { 𝑡 dyn ă We present a steady-state model of the inner protoplanetary discwhich accretes viscously, primarily due to the MRI. In this model,the disc is heated by viscous dissipation and stellar irradiation, andcools radiatively and convectively. The disc is ionized by thermalionization, thermionic and ion emission from dust grains and bystellar X-rays, cosmic rays and radionuclides, and we also accountfor adsorption of charges onto dust grains. The disc’s structure (den-sity, temperature), viscosity due to the MRI, opacity and ionizationstate are calculated self-consistently everywhere in the disc (both asa function of radius and height). To the best of our knowledge thisis the first model that self-consistently couples all these processesto describe the structure of the inner regions of a steadily accretingprotoplanetary disc.We investigate how these various processes affect the structureof the inner disc and the extent to which the MRI can drive efficientaccretion, i.e., the locations of the inner edge of the dead zone andthe gas pressure maximum. For the fiducial parameters consideredin this work (stellar parameters: 𝑀 ˚ “ d , 𝑅 ˚ “ d , 𝑇 ˚ “ 𝐿 X “ ´ . 𝐿 bol ; disk parameters: gasaccretion rate 𝑀 “ ´ M d yr ´ , viscosity in the MRI-dead zone MNRAS000
1. Everywhere below this line thechemical equilibrium timescale is shorter than the dynamical timescale,justifying our assumption of the strong-coupling regime. See Section 4.5. fulfilled in most of the inner disc, as slow radiative recombinationsmake the chemical equilibrium timescale long. However, even in theabsence of dust, three-body recombination (recombinations throughcollisions with the abundant molecular hydrogen) are much fasterthan radiative recombinations, and adsorption onto grains is evenfaster (Desch & Turner 2015).Since we solve directly for the equilibrium ionization state, wedo not have access to the timescale 𝑡 chem . However, we can estimateit as 𝑡 chem “ 𝑛 e { R , where R is the fastest of the above threerecombination rates (in general, but not always, that is adsorptiononto dust grains). Fig. 17 shows, for our fully self-consistent modelwith our full chemical network, that 𝑡 chem { 𝑡 dyn ă We present a steady-state model of the inner protoplanetary discwhich accretes viscously, primarily due to the MRI. In this model,the disc is heated by viscous dissipation and stellar irradiation, andcools radiatively and convectively. The disc is ionized by thermalionization, thermionic and ion emission from dust grains and bystellar X-rays, cosmic rays and radionuclides, and we also accountfor adsorption of charges onto dust grains. The disc’s structure (den-sity, temperature), viscosity due to the MRI, opacity and ionizationstate are calculated self-consistently everywhere in the disc (both asa function of radius and height). To the best of our knowledge thisis the first model that self-consistently couples all these processesto describe the structure of the inner regions of a steadily accretingprotoplanetary disc.We investigate how these various processes affect the structureof the inner disc and the extent to which the MRI can drive efficientaccretion, i.e., the locations of the inner edge of the dead zone andthe gas pressure maximum. For the fiducial parameters consideredin this work (stellar parameters: 𝑀 ˚ “ d , 𝑅 ˚ “ d , 𝑇 ˚ “ 𝐿 X “ ´ . 𝐿 bol ; disk parameters: gasaccretion rate 𝑀 “ ´ M d yr ´ , viscosity in the MRI-dead zone MNRAS000 , 1–21 (2020) 𝛼 DZ “ ´ ; dust parameters: dust-to-gas mass ratio 𝑓 dg “ ´ ,maximum grain size 𝑎 max “ 𝜇 m), we find that:(i) Inwards of the pressure maximum, the MRI is active onlyaround the disc midplane. This differs from the predictions ofvertically-isothermal models, and is possibly important for the evo-lution of dust grains in the inner disc.(ii) Since the inner disc is optically thick, stellar irradiation onlyweakly influences the midplane temperature, and thus also onlyweakly affects the location of the dead zone inner edge.(iii) Most of the inner disc is convectively unstable, which weshow is a property of any optically thick disc in which viscousheating occurs near the midplane. This motivates further work intoa coupled MRI-convective instability in the limit of non-ideal MHD.(iv) As suggested by Desch & Turner (2015), dust controls theionization state of the inner disc, and thus the onset of the MRI.Thermal ionization plays a secondary role, as thermionic and ionemission from dust grains ionize the hot dense regions.(v) High above disc midplane stellar X-rays produce an MRI-active layer. However, the X-rays barely change the overall viscosityat short orbital distances, or the location of the pressure maximum.(vi) The pressure maximum resides at „ ´ and small dust grains ( 𝑎 max “ 𝜇 m), whichmay be expected in the early stages of dust evolution in the disc.How these results depend on the model parameters, including dustgrain size and dust-to-gas ratio, is explored in a companion paper(Jankovic et al. in prep. ), where we also use our new inner discmodel to speculate on possible formation pathways for close-insuper-Earths. ACKNOWLEDGEMENTS
We thank the reviewer for helpful suggestions that improved themanuscript. We thank Richard Booth, Thomas Haworth, ZhaohuanZhu, Steven Desch, Neal Turner and Colin McNally for helpful dis-cussions. M.R.J. acknowledges support from the President’s PhDscholarship of the Imperial College London and the UK Science andTechnology research Council (STFC) via the consolidated grantST/S000623/1. JEO is supported by a Royal Society UniversityResearch Fellowship. This project has received funding from theEuropean Research Council (ERC) under the European Union’sHorizon 2020 research and innovation programme (Grant agree-ment No. 853022, ERC-STG-2019 grant, PEVAP). JCT acknowl-edges support from NASA ATP grant Inside-Out Planet Formation(80NSSC19K0010).
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Bai X.-N., 2011, ApJ, 739, 50 Bai X.-N., Goodman J., 2009, ApJ, 701, 737Bai X.-N., Stone J. M., 2011, ApJ, 736, 144Balbus S. A., 2000, ApJ, 534, 420Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577, A42Bitsch B., Raymond S. N., Izidoro A., 2019, A&A, 624, A109Bodo G., Cattaneo F., Mignone A., Rossi P., 2013, ApJ, 771, L23Boley A. C., Ford E. B., 2013, arXiv e-prints, p. arXiv:1306.0566Calvet N., Magris G. C., Patino A., D’Alessio P., 1992, Rev. Mex. Astron.Astrofis., 24, 27Cassen P., 1993, in Lunar and Planetary Science Conference. Lunar andPlanetary Science Conference. p. 261Chatterjee S., Tan J. C., 2014, ApJ, 780, 53Chiang E. I., Goldreich P., 1997, ApJ, 490, 368Chiang E. I., Joung M. K., Creech-Eakman M. J., Qi C., Kessler J. E., BlakeG. A., van Dishoeck E. F., 2001, ApJ, 547, 1077Cossou C., Raymond S. N., Hersant F., Pierens A., 2014, A&A, 569, A56D’Alessio P., Cantö J., Calvet N., Lizano S., 1998, ApJ, 500, 411D’Alessio P., Calvet N., Hartmann L., Lizano S., Cantó J., 1999, ApJ, 527,893Desch S. J., Turner N. J., 2015, ApJ, 811, 156Draine B. T., 2003, ApJ, 598, 1017Draine B. T., Sutin B., 1987, ApJ, 320, 803Dressing C. D., Charbonneau D., 2013, ApJ, 767, 95Dressing C. D., Charbonneau D., 2015, ApJ, 807, 45Dullemond C. P., 2002, A&A, 395, 853Dullemond C. P., Dominik C., Natta A., 2001, ApJ, 560, 957Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Com-mercon B., Flock M., 2012, RADMC-3D: A multi-purpose radiativetransfer tool (ascl:1202.015)Ercolano B., Glassgold A. E., 2013, MNRAS, 436, 3446Estrada P. R., Cuzzi J. N., Morgan D. A., 2016, ApJ, 818, 200Flock M., Turner N. J., Mulders G. D., Hasegawa Y., Nelson R. P., BitschB., 2019, A&A, 630, A147Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics: ThirdEditionFressin F., et al., 2013, ApJ, 766, 81Gammie C. F., 1996, ApJ, 457, 355Garaud P., Lin D. N. C., 2007, ApJ, 654, 606Glassgold A. E., Najita J., Igea J., 1997, ApJ, 480, 344Hansen B. M. S., Murray N., 2013, ApJ, 775, 53Hartmann L., Calvet N., Gullbring E., D’Alessio P., 1998, ApJ, 495, 385Held L. E., Latter H. N., 2018, MNRAS, 480, 4797Hirose S., 2015, MNRAS, 448, 3105Hirose S., Blaes O., Krolik J. H., Coleman M. S. B., Sano T., 2014, ApJ,787, 1Hsu D. C., Ford E. B., Ragozzine D., Ashby K., 2019, AJ, 158, 109Hu X., Tan J. C., Zhu Z., Chatterjee S., Birnstiel T., Youdin A. N., MohantyS., 2018, ApJ, 857, 20Hubeny I., 1990, ApJ, 351, 632Igea J., Glassgold A. E., 1999, ApJ, 518, 848Ilgner M., Nelson R. P., 2006, A&A, 445, 205Izidoro A., Ogihara M., Raymond S. N., Morbidelli A., Pierens A., BitschB., Cossou C., Hersant F., 2017, MNRAS, 470, 1750Izidoro A., Bitsch B., Raymond S. N., Johansen A., Morbidelli A., Lam-brechts M., Jacobson S. A., 2019, arXiv e-prints, p. arXiv:1902.08772Jankovic M. R., Owen J. E., Mohanty S., 2019, MNRAS, 484, 2296Keith S. L., Wardle M., 2014, MNRAS, 440, 89Lesur G., Ogilvie G. I., 2010, MNRAS, 404, L64Lin D. N. C., Papaloizou J., 1980, MNRAS, 191, 37Lodders K., 2003, ApJ, 591, 1220Malygin M. G., Kuiper R., Klahr H., Dullemond C. P., Henning T., 2014,A&A, 568, A91Manara C. F., et al., 2016, A&A, 591, L3Manara C. F., et al., 2017, A&A, 604, A127Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425McNeil D. S., Nelson R. P., 2010, MNRAS, 401, 1691Mohanty S., Ercolano B., Turner N. J., 2013, ApJ, 764, 65MNRAS , 1–21 (2020)
RI-active inner disc Mohanty S., Jankovic M. R., Tan J. C., Owen J. E., 2018, ApJ, 861, 144Mulders G. D., Pascucci I., Apai D., Ciesla F. J., 2018, AJ, 156, 24Natta A., Prusti T., Neri R., Wooden D., Grinin V. P., Mannings V., 2001,A&A, 371, 186Ogihara M., Ida S., 2009, ApJ, 699, 824Owen J. E., Wu Y., 2017, ApJ, 847, 29Papaloizou J. C. B., Terquem C., 1999, ApJ, 521, 823Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush T., FongW., 1994, ApJ, 421, 615Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 2002, Nu-merical recipes in C++ : the art of scientific computingRafikov R. R., 2007, ApJ, 662, 642Ruediger G., Elstner D., Tschaepe R., 1988, Acta Astron., 38, 299Salmeron R., Wardle M., 2008, MNRAS, 388, 1223Sano T., Stone J. M., 2002, ApJ, 577, 534Sano T., Miyama S. M., Umebayashi T., Nakano T., 2000, ApJ, 543, 486Sano T., Inutsuka S.-i., Turner N. J., Stone J. M., 2004, ApJ, 605, 321Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33Shu F. H., 1992, The physics of astrophysics. Volume II: Gas dynamics.Suzuki T. K., Ogihara M., Morbidelli A., Crida A., Guillot T., 2016, A&A,596, A74Terquem C., Papaloizou J. C. B., 2007, ApJ, 654, 1110Umebayashi T., Nakano T., 1981, PASJ, 33, 617Umebayashi T., Nakano T., 2009, ApJ, 690, 69Van Eylen V., Agentoft C., Lundkvist M. S., Kjeldsen H., Owen J. E., FultonB. J., Petigura E., Snellen I., 2018, MNRAS, 479, 4786Wardle M., 1999, MNRAS, 307, 849Wardle M., 2007, Ap&SS, 311, 35Wright N. J., Drake J. J., Mamajek E. E., Henry G. W., 2011, ApJ, 743, 48Wu Y., 2019, ApJ, 874, 91Zink J. K., Christiansen J. L., Hansen B. M. S., 2019, MNRAS, 483, 4479This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000