Magnetic fields in gaps surrounding giant protoplanets
MMon. Not. R. Astron. Soc. , 1–14 (Year) Printed 6 May 2019 (MN L A TEX style file v2.2)
Magnetic fields in gaps surrounding giant protoplanets
Sarah L. Keith , (cid:63) and Mark Wardle Department of Physics & Astronomy and Research Centre in Astronomy, Astrophysics & Astrophotonics, Macquarie University,NSW 2109, Australia Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester, M13 9PL,United Kingdom
Accepted 2015 May 5. Received 2015 May 4; in original form 2015 February 4
ABSTRACT
Giant protoplanets evacuate a gap in their host protoplanetary disc, which gas mustcross before it can be accreted. A magnetic field is likely carried into the gap, po-tentially influencing the flow. Gap crossing has been simulated with varying degreesof attention to field evolution (pure hydrodynamical, ideal, and resistive MHD), butas yet there has been no detailed assessment of the role of the field accounting forall three key non-ideal MHD effects: Ohmic resistivity, ambipolar diffusion, and Halldrift.We present a detailed investigation of gap magnetic field structure as determinedby non-ideal effects. We assess susceptibility to turbulence induced by the magnetoro-tational instability, and angular momentum loss from large-scale fields.As full non-ideal simulations are computationally expensive, we take an a posteri-ori approach, estimating MHD quantities from the pure hydrodynamical gap crossingsimulation by Tanigawa et al. (2012). We calculate the ionisation fraction and estimatefield strength and geometry to determine the strength of non-ideal effects.We find that the protoplanetary disc field would be easily drawn into the gap andcircumplanetary disc. Hall drift dominates, so that much of the gap is conditionallyMRI unstable depending on the alignment of the field and disc rotation axes. Fieldalignment also influences the strong toroidal field component permeating the gap.Large-scale magnetic forces are small in the circumplanetary disc, indicating theycannot drive accretion there. However, turbulence will be key during satellite growthas it affects critical disc features, such as the location of the ice line.
Key words: accretion discs – magnetic fields – MHD – planets and satellites: for-mation
Giant planets capture their massive atmospheres from thesurrounding protoplanetary disc. As the planet grows, itsgravitational sphere of influence expands so that gas is cap-tured faster than the nebula can resupply it, and a gapis evacuated in the nebula around the planet (Lin & Pa-paloizou 1985; Artymowicz & Lubow 1996; Bryden et al.1999). Gaps have been seen in ∼ µ m scattered light aroundTW Hya and infrared polarimetry of surface dust aroundHD97048 (Debes et al. 2013; Garufi et al. 2014). Spiralwaves, likely launched by a protoplanet, have also been seenin HCO + line emission within a gap in the HD 142527 cir-cumstellar disc (Casassus et al. 2013).A circumplanetary disc encircles the planet following (cid:63) E-mail: [email protected];[email protected] the collapse and detachment of the protoplanet envelopefrom the nebula (Lunine & Stevenson 1982; Ayliffe & Bate2009a). Gas flowing across the gap supplies the circum-planetary disc and controls the protoplanetary accretionrate (Lubow & D’Angelo 2006). Gas with too much angu-lar momentum to reach the planet directly is captured bythe circumplanetary disc and its flow is controlled by an-gular momentum transport processes operating within thecircumplanetary disc (Fujii et al. 2011; Lubow & Martin2012; Rivier et al. 2012; Tanigawa et al. 2012, hereafterTOM12; Fujii et al. 2014; Keith & Wardle 2014, hereafterKW14; Szulágyi et al. 2014; Turner et al. 2014). The cir-cumplanetary flow pattern potentially comprises of high-altitude inflow, and midplane Keplerian and outflow com-ponents (Bate et al. 2003; D’Angelo et al. 2003, TOM12,Gressel et al. 2013). Although circumplanetary discs areyet to be observed, there are prospects for their detectionthrough infrared spectral energy distributions observed with c (cid:13) Year RAS a r X i v : . [ a s t r o - ph . E P ] M a y Sarah L. Keith and Mark Wardle the James Webb Space Telescope or the Atacama LargeMillimeter/submillimeter Array (ALMA; Isella et al. 2014;Dunhill 2015; Zhu 2015). This is strengthened by the recentobservation of discs around three brown dwarfs, in infraredcontinuum and CO line emission (Ricci et al. 2014).Magnetic fields likely play a role in the dynamics ofthe gap and circumplanetary disc system (Gressel et al.2013; Uribe et al. 2013, KW14, Turner et al. 2014). Thefield can have numerous effects, including generating hydro-magnetic turbulence via the magnetorotational instability(MRI; Balbus & Hawley 1991; Hawley et al. 1995), centrifu-gally driven disc winds and jets (Blandford & Payne 1982;Wardle & Königl 1993), and magnetic braking (Matsumoto& Tomisaka 2004).Magnetic forces are transmitted by charged particles,and so these mechanisms require a minimum ionisation frac-tion. Collisions between charged and neutral particles giverise to non-ideal effects which counteract flux freezing -Ohmic resistivity, Hall drift, and ambipolar diffusion. Suf-ficiently strong resistivity and ambipolar diffusion decouplethe motion of the gas and field (see Turner et al. 2014 andreferences within), and Hall drift can lead to complex fieldevolution sensitive to the global field orientation (Wardle1999; Wardle & Salmeron 2012; Kunz & Lesur 2013; Lesuret al. 2014, O’Keeffe & Downes 2014, Bai 2015).Gap-crossing dynamics is complex and so studies haverelied on numerical simulations to model the flow. These in-clude both grid-based (Bryden et al. 1999; Lubow et al. 1999;Bate et al. 2003; Lubow & D’Angelo 2006, TOM12) andsmoothed particle (Bryden et al. 1999; Ayliffe & Bate 2009b,2012) hydrodynamical simulations, and ideal magnetohydro-dynamical (MHD) simulations (Nelson & Papaloizou 2003;Papaloizou et al. 2004; Uribe et al. 2013). Global 3D MHDsimulations of the wider system have also been used to gen-erate synthetic observation maps for ALMA (Flock et al.2015).The recent inclusion of Ohmic resistivity adds an im-portant level of realism to gap-crossing modelling (Gresselet al. 2013). Ohmic resistivity is strong in the dense cir-cumplanetary disc, and can produce an MRI stable, deadzone with implications for planet growth and disc evolution.However, non-ideal effects are also at work in relatively lowdensity regions where the Hall effect and ambipolar diffusioncan be strong (Wardle 2007). Their inclusion may enhanceor suppress turbulence, affecting the flow and protoplanetgrowth rate.Indeed, the addition of ambipolar diffusion toresistive simulations radically alters the protoplanetary discflow between 1–5 au, producing laminar accretion poweredby a magnetocentrifugal wind, rather than the traditionalaccreting turbulent surface layers (Gressel et al. 2015).In this article we perform a comprehensive study of therelative importance of the three non-ideal effects on fluid dy-namics throughout the gap. Owing to the significant com-putational cost of including all three effects in a full, 3Dsimulation, we take an a posteriori, semi-analytic approach.We calculate MHD quantities using a snapshot from theTOM12 3D hydrodynamical simulation (described in Sec-tion 2.1). This allows us to calculate detailed ionisationmaps including ionisation from cosmic-rays, stellar X-raysand radioactive decay, accounting for grain charging (Sec-tions 2.2 and 3.1). We determine the strength of non-idealeffects (Sections 2.3 and 3.2) to estimate the field strength and geometry (Sections 2.4 and 3.3). We find that the fieldwould not alter the general fluid motion, justifying our useof an underlying hydrodynamical gap-crossing model (i.e.,TOM12). The implications of incorporating magnetic fieldsand non-ideal effects in gap crossing, such as the additionalforce from large-scale fields, are considered (Sections 4.1 and4.2). Finally, we present a summary and discussion of find-ings in Section 5.
In this section we outline the disc model used to describethe protoplanetary disc and gap. We give details of the ion-isation and diffusivity calculations, along with estimates forthe magnetic field strength and geometry.
We take a semi-analytic, a posteriori approach to the calcu-lation. We use a snapshot of a pure hydrodynamical simu-lation as the basis for estimating MHD quantities. The pro-toplanetary disc and gap model that we use is the TOM12three-dimensional, hydrodynamical simulation. The simula-tion models a protoplanetary disc surrounding a star of mass M ∗ , in which a gap has been carved out through gas cap-ture by an embedded protoplanet of mass M p , at fixed or-bital radius d p . The simulation corresponds to a protoplanetwith the present-day mass and orbital radius of Jupiter (i.e., M p = M J , and d p = d J , where d J ≡ . au). Gas is treatedas inviscid, and self-gravitational and magnetic forces areneglected.The planet is located at the origin of a local cartesianco-ordinate system ( x, y, z ) , which orbits the star at an or-bital angular frequency Ω p . The x -axis is oriented along theradial direction, ˆ r , which extends from the star to the planet;the y -axis is oriented in the azimuthal direction, ˆ φ , parallelto the the protoplanet orbit; and the z -axis is parallel tothe protoplanetary disc angular momentum vector. To cap-ture details of the fine flow structure near the protoplanetthe simulation has eleven levels of nested grids, each with n x × n y × n z = 64 × × data points. The total simu-lated portion of the disc extends over x ∈ [ − H p , H p ] , y ∈ [ − H p /d p , H p /d p ] , and z ∈ [0 , H p ] where H p is thescale height of the protoplanetary disc. The local cartesianapproximation requires that the Hill radius, R H = d (cid:18) M M ∗ (cid:19) ≈ . au (cid:18) d . au (cid:19) (cid:18) MM J (cid:19) (cid:18) M ∗ M (cid:12) (cid:19) − , (1)which is roughly the feeding zone of the protoplanet, satisfies R H (cid:28) d p . Here M (cid:12) is the mass of the Sun. The simulationdata used here was taken after 160.7 protoplanet orbits, al-lowing the simulation to approach a steady state. The gapreached a maximum density contrast relative to the unper-turbed disc of ∼ .The simulation used a non-dimensionalised form forthe MHD equations and so we must rescale the velocity v = ( v x , v y , v z ) and density ρ data for our calculations. c (cid:13) Year RAS, MNRAS000
We take a semi-analytic, a posteriori approach to the calcu-lation. We use a snapshot of a pure hydrodynamical simu-lation as the basis for estimating MHD quantities. The pro-toplanetary disc and gap model that we use is the TOM12three-dimensional, hydrodynamical simulation. The simula-tion models a protoplanetary disc surrounding a star of mass M ∗ , in which a gap has been carved out through gas cap-ture by an embedded protoplanet of mass M p , at fixed or-bital radius d p . The simulation corresponds to a protoplanetwith the present-day mass and orbital radius of Jupiter (i.e., M p = M J , and d p = d J , where d J ≡ . au). Gas is treatedas inviscid, and self-gravitational and magnetic forces areneglected.The planet is located at the origin of a local cartesianco-ordinate system ( x, y, z ) , which orbits the star at an or-bital angular frequency Ω p . The x -axis is oriented along theradial direction, ˆ r , which extends from the star to the planet;the y -axis is oriented in the azimuthal direction, ˆ φ , parallelto the the protoplanet orbit; and the z -axis is parallel tothe protoplanetary disc angular momentum vector. To cap-ture details of the fine flow structure near the protoplanetthe simulation has eleven levels of nested grids, each with n x × n y × n z = 64 × × data points. The total simu-lated portion of the disc extends over x ∈ [ − H p , H p ] , y ∈ [ − H p /d p , H p /d p ] , and z ∈ [0 , H p ] where H p is thescale height of the protoplanetary disc. The local cartesianapproximation requires that the Hill radius, R H = d (cid:18) M M ∗ (cid:19) ≈ . au (cid:18) d . au (cid:19) (cid:18) MM J (cid:19) (cid:18) M ∗ M (cid:12) (cid:19) − , (1)which is roughly the feeding zone of the protoplanet, satisfies R H (cid:28) d p . Here M (cid:12) is the mass of the Sun. The simulationdata used here was taken after 160.7 protoplanet orbits, al-lowing the simulation to approach a steady state. The gapreached a maximum density contrast relative to the unper-turbed disc of ∼ .The simulation used a non-dimensionalised form forthe MHD equations and so we must rescale the velocity v = ( v x , v y , v z ) and density ρ data for our calculations. c (cid:13) Year RAS, MNRAS000 , 1–14 agnetic fields in protoplanetary gaps We rescale the column density using the minimum mass So-lar nebula column density at the orbital radius of Jupiter Σ = 140 g cm − (Kitamura et al. 2002; Williams & Cieza2011). This model estimates the gas mass in the Solar neb-ula by augmenting the solid material in Solar System planetswith sufficient hydrogen and helium to bring the composi-tion up to Solar (Weidenschilling 1977; Hayashi 1981). Ob-servations of Class II Young Stellar Objects (i.e., those withdiscs and strong UV and H α emission indicating active ac-cretion) are currently limited to ∼ au resolution, howeverinferred disc column density profiles are broadly consistentwith this profile. As the gap gas density naturally drops overtime, we allow for simple reduction of the column density, Σ( x, y, z ) , though multiplication of a constant parameter f Σ ,where f Σ = 1 recovers the original TOM12 results.The simulation is isothermal, and we adopt the tem-perature of a black body in radiative equilibrium with So-lar luminosity, excluding any heating from inflow processes,at the orbital radius of Jupiter, T = 120 K (Wardle 2007).This is supported by interferometric CO images of T Tauriand dust temperatures from SED modelling (Guilloteau &Dutrey 1998; Andrews & Williams 2005; Piétu et al. 2007).We adopt Solar gas composition, with a mixture of80% molecular hydrogen and 20% atomic helium. This cor-responds to a mean molecular weight of m n = 2 . m p ,where m p is the mass of a proton. This allows us to calcu-late the neutral number density, n , isothermal sound speed c s = 0 . km s − , and aspect ratio: H p d = c s v k = 4 . × − (cid:18) dd J (cid:19) , (2)where v k = ( GM ∗ /d ) is the Keplerian velocity.The flow pattern within a gap transitions from orbit-ing the star [orbital radius (cid:112) ( x + d p ) + y ], to orbiting theplanet (orbital radius (cid:112) x + y ) within the Hill sphere. Inour analytic calculations we approximate the flow geometryby a composite Keplerian angular velocity functions: Ω = (cid:104) GM p (cid:0) x + y (cid:1) − (cid:105) for (cid:112) x + y ≤ R H , (cid:2) GM ∗ ( x + d p ) − (cid:3) otherwise . (3)TOM12 also simulated the circumplanetary disc, how-ever they note that as it is shear-dominated, artificial vis-cosity is strong and so the results are more reliable in thegap. Similarly we include the circumplanetary disc in ourcalculations with the caveat that the disc structure in thisregion is uncertain.The scale height in the circumplanetary disc, H c = H p (cid:32) (cid:112) x + y R H (cid:33) (4) ≈ . × − H p (cid:32) (cid:112) x + y . R H (cid:33) , (5)is significantly lower, and this needs to be accounted forwhen calculating gradients in the fluid flow (see Section 2.3).To that end we use a composite scale-height, H , which tran-sitions at the boundary of the circumplanetary disc hydro- static region: H = (cid:26) H c for (cid:112) x + y < . R H , z/H c < , and H p otherwise . (6) Charged particles transmit magnetic forces to the bulk neu-tral flow. In the gap, ionisation is largely caused by penetrat-ing cosmic-rays and stellar X-rays, however there is also aweak, pervasive contribution from decaying nuclides. We cal-culate the ionisation level by solving the coupled rate equa-tions for electron and ion number densities and mean graincharge subject to overall charge neutrality. We follow themethod outlined in KW14, but using a more detailed cal-culation of the grain charge. A summary of the method isgiven below.We calculate the grain number density, n g , for sphericalgrains with radius a g = 0 . µ m and bulk density 3 g cm − with a constant dust to gas mass ratio f dg ∼ − to ac-count for incorporation of solids into protoplanetary bodies(Pollack et al. 1994).Internal ionisation from radioactive decay is primarilyfrom the short-lived radionuclide Al, which contributes anionisation rate ζ R = 7 . × − s − ( f dg /f (cid:12) ) (Umebayashi& Nakano 2009). The ionisation rate has been scaled bythe ratio of the dust mass fraction to the Solar photosphericmetallicity, f (cid:12) = 1 . , so that the Al abundance is consis-tent (Asplund et al. 2009). The midplane dust fraction willbe enhanced by settling, however the consequent increase inthe ionisation fraction will be small as ionisation is predom-inantly by external radiation.External ionisation sources are attenuated by the col-umn density above and below each location in the simula-tion. Cosmic-ray ionisation occurs at the interstellar cosmic-ray ionisation rate, ζ CR = 10 − s − , attenuated by the at-tenuation depth Σ CR = 96 g cm − approaching from aboveand below the disc (Umebayashi & Nakano 1981, 2009).We account for ionisation from diffuse scattered stellar X-rays but neglect direct illumination, instead assuming thatshielding by the inner portion of the disc is effective. Ion-isation for solar luminosity stellar X-rays occurs at a rate ζ XR = 9 . × − s − ( d/d J ) − , with attenuation depth Σ XR = 8 g cm − (Igea & Glassgold 1999; Turner & Sano2008). The total ionisation rate, ζ , is the sum from the threesources ζ = ζ R + ζ CR + ζ XR .We solve for the electron, ion and grain charge numberdensities by calculating the grain charge needed for overallcharge neutrality in the steady state [see equations (26)–(29) in KW14]. We include charge focussing in the electronand ion capture rate coefficients for grains (Umebayashi &Nakano 1980): k ig = πa g (cid:114) k b Tπm i × exp (cid:16) − q | Z g | a g k b T (cid:17) Z g < (cid:16) q | Z g | a g k b T (cid:17) Z g > (7) k eg = πa g (cid:114) k b Tπm e × (cid:16) q | Z g | a g k b T (cid:17) Z g < (cid:16) − q | Z g | a g k b T (cid:17) Z g > (8)where Z g q is the mean grain charge, k b is Boltzmann’s con-stant, and m e is the electron mass. We model metals as asingle species, adopting the mass, m i , and abundance, x i , of c (cid:13) Year RAS, MNRAS , 1–14
Sarah L. Keith and Mark Wardle the most abundant metal - magnesium (Lide 2004; Asplundet al. 2009). The recombination rate for ions and electrons is k ei = 1 . × − ( T / K ) − . cm s − (McElroy et al.2013). If the grain charge is known these equations give theion and electron number densities: n e = n g k ig k ei (cid:34)(cid:18) k ei nζk eg k ig n g (cid:19) − (cid:35) (9) n i = n e k eg k ig . (10)Equations (7)–(10) are analogous to equations (30)–(33) in KW14 except for the inclusion of a non-zero graincharge for charge focussing. The protoplanetary disc hasa higher average ionisation fraction resulting in a highergrain charge which cannot be neglected (as it could bein the KW14 circumplanetary disc analysis). We deter-mine the equilibrium values of n e , n i , and Z g numericallyby solving equations (7)–(10) along with charge neutral-ity, for the grain charge. We use the Brent-Dekker method, gsl_root_fsolver_brent , for root-finding implemented inthe GNU Scientific Library to solve for the grain charge toan accuracy of 1% (Galassi et al. 2009). This method com-bines bisection and secant methods for rapid convergence.We also determine the effectiveness of collisional ionisa-tion produced in MRI turbulent regions. Currents generatedby an MRI field may be accelerated sufficiently to ionise neu-tral particles, providing additional ionisation (Muranushiet al. 2012; Okuzumi & Inutsuka 2015). We calculate thekinetic energy of accelerated particles using equations (37)–(38) in KW14 to determine if it exceeds the ionisation po-tential of neutral species (e.g., potassium has the lowest ion-isation energy of the abundant elements, with E ∼ . eV).We find that although electrons reach energies of up to 1 eVin the disc atmosphere ( z ∼ H p ) where collisions are infre-quent and the mean-free path is long, in the denser midplaneelectrons are not accelerated above . eV. This is too lowto ionise any atomic species and so the MRI does not con-tribute significantly to the ionisation. Thus, we neglect theeffect in our calculation of the ionisation fraction.Finally, we note that thermal ionisation is likely not ef-fective in the gap surrounding Jupiter. The protoplanetarydisc at d = d J is much cooler than the required T ∼ K,and further shock heating of accreting material at the cir-cumplanetary disc is weak. In the circumplanetary shock,shock heating flux F s ∼ σ (40 K ) (calculated for an accre-tion rate − M J / year, at r ∼ R J from the protoplanet;Cassen & Moosman 1981) is negligible compared with Solarirradiation, F (cid:12) = σ (120 K ) at this orbital radius, preclud-ing thermal ionisation. The evolution of the magnetic field and its role in the gasdynamics depends on the strength and nature of the mag-netic coupling between the field and gas. Magnetic couplingis manifest in the induction equation, ∂ B ∂t = ∇ × ( v × B ) − ∇ × [ η O ( ∇ × B ) + η H ( ∇ × B ) × ˆ B ] −∇ × [ η A ( ∇ × B ) ⊥ ] , (11) where η O , η H , and η A are the transport coefficients forOhmic resistivity, Hall drift, and ambipolar diffusion. Thesenon-ideal effects cause the magnetic field to slip through thegas, and can dominate the field evolution if they are toostrong. They are caused by neutral collisions disrupting the E × B drift of charged particles. Differences in the densitydependence of these effects leads to three different regimesin which Ohmic resistivity, Hall drift and ambipolar diffu-sion are most effective in high, moderate, and low densityregions, respectively. For a predominantly vertical magneticfield, Ohmic resistivity and Ambipolar diffusion act similarly(and distinctly from the Hall drift), and so they are oftencombined into the Pedersen diffusion η P = η O + η A . We cal-culate the strength of the three effects using the transportcoefficients from equations (53)–(55) in KW14 .How strong these effects must be to dominate the in-ductive term and decouple the motion of the field and gasdepends on the field geometry. Below we outline the cou-pling conditions for MRI turbulent, toroidal, and poloidalcomponents in turn.The MRI dynamo harnesses shear to generate a tur-bulent field. Diffusion restricts MRI turbulence by limitingbending of the smallest-scale field modes. The MRI oper-ates effectively if the fastest growing mode, with wavelength λ = 2 πv a / Ω , survives diffusion. Magnetic diffusion counter-acts dynamo growth of short wavelength fluctuations withwavelength λ < η/v a,z , for which the diffusion rate exceedsthe growth rate. Here v a,z = B z / √ πρ the vertical Alfvénspeed. Thus, the MRI operates in effectively ideal MHD con-ditions if Pedersen and Hall components are simultaneouslybelow the threshold (Sano & Stone 2002; Turner et al. 2014): η P (cid:46) v a,z Ω − and | η H | (cid:46) v a,z Ω − . (12)Here, the absolute value of the Hall term reflects its signednature. The MRI also requires the field must be weak enoughthat the fastest-growing mode is confined within the discscale height, so that for ideal MRI (Okuzumi & Hirose 2011): β z = 8 πρc s B z > π . (13)Unlike Ohmic resistivity and ambipolar diffusion, Halldrift is not diffusive and can change the behaviour of MRIturbulence, even in the presence of strong Ohmic and am-bipolar diffusion. Hall drift is antiparallel to the current den-sity and may cooperate with, or act against Keplerian shearaccording to the relative orientation of the protoplanetarydisc angular momentum vector and vertical component ofthe magnetic field, s = sign ( B · Ω ) (Wardle 1999; Balbus& Terquem 2001; Wardle & Salmeron 2012). For example, There is a typographical error in the expression for Ohmicresistivity in equation (53) of KW14. The scaling of thediffusivity with the neutral density should be η O ∝ ρ rather than the inverse proportionality given. Additionally, thegrain Hall Parameter in equation (51) of KW14 is a fac-tor of (cid:112) m g /m n too large, and should be written as β g =5 . × − Z g ( B/ G ) (cid:0) cm − /n (cid:1) (0 . µ m /a g ) (cid:112) K /T ,instead. This reduction has no impact on the KW14 results, as β g was already so small as to be negligible. Indeed, Kunz & Lesur (2013) found that the Hall effect can be anti-diffusive in unstratified shearing boxes.c (cid:13)
Year RAS, MNRAS000
Year RAS, MNRAS000 , 1–14 agnetic fields in protoplanetary gaps if the field is aligned with the rotation axis, the Hall ef-fect can destabilise the flow, enhancing turbulence. On theother hand, if the field and rotation axis are anti-aligned,Hall drift acts against the shear, suppressing turbulence. Inthe Hall-MRI regime, the Hall effect counteracts Pedersenstabilisation against field tangling if it dominates and thevertical field orientation is favourable (Wardle & Salmeron2012; Lesur et al. 2014). Hall-MRI requires | η H | > η P , | η H | > v a,z Ω − , and s = 1 , (14)along with a Hall-MRI weak-field criterion, which we showin Appendix A is the same as equation (13). Indeed, the Halleffect need not dominate diffusion to impact turbulence, aseven weak Hall drift modifies the wave number and growthrate of the fastest growing mode (Wardle & Salmeron 2012).A toroidal field surrounding the star or protoplanet,must also be maintained against diffusion which would tendto unwrap and straighten field lines. The impact of non-ideal effects on a toroidal filed are most important close tothe central object where gradients are strongest. A toroidalfield is preserved against diffusion if magnetic induction ex-ceeds non-ideal effects, as captured by the coupling threshold(c.f., Turner & Sano 2008; see Appendix B): η P (cid:46) Ω H , and | η H | (cid:46) Ω H . (15)A toroidal field can also be preserved against strong dif-fusion if Hall drift dominates. Hall drift is along the current,which is in the fluid rotation direction for both the alignedand anti-aligned cases (i.e., s = 1 and s = − ). So, either byenhancing the field wrapping (i.e., for s = 1 ) or by unwind-ing and rewrapping the field in the counter flow direction(i.e., for s = − ), the Hall effect produces a strong negativetoroidal field. See Appendix B for a discussion of the effect ofvertical field direction dependence for coupling of a toroidalfield. Simulations are needed to verify this behaviour, butwe highlight Hall-dominated regions with the anticipation itmay extend the reach of the toroidal field. We simply requirethat azimuthal winding from Hall drift exceeds diffusion: | η H | > η P , | η H | > Ω H . (16)Finally, the disc contains a poloidal field componentthroughout. Although magnetic coupling limits the mini-mum field bending radius of the poloidal component, an es-sentially vertical field will always couple, owing to its smallgradient. As we are not concerned with the exact geometryof the poloidal component, we simply assume that it is ap-proximately vertical, couples everywhere and so permeatesthe gap. Protoplanetary discs inherit a large-scale magnetic field fromtheir progenitor molecular cloud. Field measurements aredifficult as current observations are limited to ∼ au res-olution (Stephens et al. 2014). The disc field is certainly com-pressionally enhanced over the cloud field ( B ∼ – mG;Shu et al. 1987), while the equipartition field B eq = 8 πp pro-vides a maximum field strength the disc will support beforemagnetic forces exceed the thermal pressure, p = ρc s .Nevertheless, it is possible to gain a more precise esti-mate from the stellar accretion rate, ˙ M ∼ − M (cid:12) yr − ,as the field is thought to play a principle role in angular momentum transport (Hartmann et al. 1998; see referencesin McKee & Ostriker 2007). This can be made from the az-imuthal component of the axisymmetric momentum equa-tion, vertically integrated between the disc surfaces (Wardle2007), ρ [( v · ∇ ) v ] φ = ( B · ∇ B ) φ π . (17)This yields a minimum field strength for vertical, toroidal,and turbulent components in the disc (Wardle 2007): B z = 6 . × − G (cid:18) ˙ M − M (cid:12) yr − (cid:19) × (cid:18) Ω1 . × − s − (cid:19) (cid:18) dd p (cid:19) − , (18) B φ = 2 . × − G (cid:18) ˙ M − M (cid:12) yr − (cid:19) × (cid:18) Ω1 . × − s − (cid:19) (cid:18) H p . au (cid:19) − , and (19) B MRI = 5 . × − G (cid:18) ˙ M − M (cid:12) yr − (cid:19) × (cid:18) Ω1 . × − s − (cid:19) (cid:18) H p . au (cid:19) − (20)respectively. Components of the turbulent field are (Sanoet al. 2004): B MRI ,r = 0 . B MRI , (21) B MRI ,φ = 0 . B MRI , and (22) B MRI ,z = 0 . B MRI . (23)We now turn to estimating the magnetic field strengthin the gap. Turbulence would be continuously generatedthere so that equations (20)–(23) remain valid, however ourestimates of the vertical and toroidal fields are not sufficientin the gap.The vertical component would be drawn in with the flowand reduced or enhanced in keeping with the large densityvariation across the gap. We assume that the field is frozeninto the gas at the outer edge of the gap (at x g = 1 . au),drawn into the gap by the flow. The self-consistency of ouradoption of this assumption is verified by our results (seeSections 3.2, 4.2). The field would be pinned to the flow atthe lowest magnetically-coupled gas layer, reducing or en-hancing the field according to B ∝ Σ − . Nevertheless, forsimplicity, we assume the field is pinned at the midplane,with a field strength B z,g = 4 . × − G, given by equation(18) at x g = 1 . au, y = 0 where Σ g = 97 g cm − . Normal-ising to the field strength at the outer edge of the gap, thisleads to a vertical field strength profile of B z = B z,g (cid:18) ΣΣ g (cid:19) (24)in the gap.Determining the toroidal field strength is considerablymore difficult as the field would be wound up continuallyuntil field-line drift balances differential rotation. As thisis beyond the scope of our treatment, we simply assumethe constant ratio (cid:112) d/H p between the toroidal and verticalcomponents in the disc is preserved for our estimate in the c (cid:13) Year RAS, MNRAS , 1–14
Sarah L. Keith and Mark Wardle gap. This leads to a gap toroidal field strength profile of: B φ = B φ,g (cid:18) ΣΣ g (cid:19) , (25)where B φ,g = 2 . × − G is the toroidal field strength,likewise taken at x g = 1 . au, y = 0 . We cap the fieldstrength so that the associated magnetic pressure does notexceed the gas pressure: β φ = 8 πρc s B φ ≥ . (26)This allows the field to expand in response to strong mag-netic pressure at high altitude, correspondingly reducing thefield strength.Numerical simulations suggest that the toroidal com-ponent may be subject to periodic polarity reversals, asit is lost through magnetic buoyancy every ∼ orbits(Miller & Stone 2000; Shi et al. 2010; Hirose & Turner2011; Flock et al. 2012). The magnitude of the toroidalfield would also vary vertically between its surface value, B φ,s , and zero at the midplane, with the approximate scal-ing B φ ∼ B φ,s ( z/H p ) , appropriate for a thin, rotationally-supported disc that is well coupled to the magnetic field(Wardle & Königl 1993). For simplicity we take B φ con-stant with height, but account for vertical gradients (e.g., inAppendix B), and probe the effect of uncertainty by globalenhancement or reduction of the the poloidal and toroidalcomponents through multiplication by a constant parameter f B . The variation of diffusion and field gradients across the gapmeans that not all three field components are present ev-erywhere. Here we determine where diffusion permits eachfield component. We use the diffusivities to determine themagnetic structure based on the level of coupling for eachfield component (poloidal, toroidal, and turbulent), and onlyinclude the components where they couple. Our procedurefor calculating the field geometry is:(i) Poloidal field - we include a poloidal component at alllocations, calculated with equation (18) in the main proto-planetary disc (i.e., | x | > x g ), or using the flux-conservedform, equation (24), inside the gap. It is safe to include apoloidal component everywhere since a pure vertical fieldalways couples and that any radial bending will adjust tothe level of diffusion.(ii) Toroidal field - we calculate the toroidal componentstrength using equation (19) in the bulk protoplanetary discflow, and equation (25) in the gap. We cap the field strengthby the gas pressure through equation (26). We calculate η O , η H , and η A using the net field strength, B = (cid:0) B z + B φ (cid:1) / ,and determine if the field is coupled using equations (15)and (16). If the field is coupled the toroidal component iskept; otherwise it is set to zero.(iii) MRI turbulent field - finally, we determine if anMRI turbulent field is sustained. We calculate the tur-bulent field strength in both the protoplanetary discand the gap using equations (20)–(23). We calculate η O , η H , and η A using the net field strength, B = (cid:2) B MRI ,r + ( B φ + B MRI ,φ ) + ( B z + B MRI ,z ) (cid:3) / , and de-termine if the field is sustained using equations (12)–(13). Here we present the results of the calculations developed inSection 2, including the ionisation fraction, magnetic field,and strength of non-ideal MHD effects. Figures are shownfor a Jupiter-mass planet, orbiting a Solar mass star, at thepresent orbital distance of Jupiter. Unless otherwise stated,we take f dg = 10 − , f Σ = 1 , f B = 1 , and s = 1 . Fig. 1 shows the ionisation profiles within the gap and disc.The left and centre panels show 2D slices through the gap,with the colour scale showing the ionisation fraction. Theleft panel shows the top-down view in the x – y plane at themidplane, and the centre panel shows the edge-on profilein the x – z plane at y = 0 . The star is located beyond thesimulation boundary at x = − d J , y = 0 . Small white patchesin the top- and bottom-right corners of the edge-on profilesare beyond the simulation boundary. The right panel showsthe dependence of the vertical ionisation profile on the dustto gas mass ratio, at x = 0 . au, y = 0 .Shielding from the overlying gas reduces the ionisationfraction with height and in denser regions. X-rays penetrateto z ∼ H p , while cosmic rays reach everywhere except forthe circumplanetary disc (located at x, y (cid:46) . au), owing toits high column density. Instead, ionisation in the circum-planetary disc is from radioactive decay, which is a muchweaker ionising source. Consequently, the ionisation frac-tion in the circumplanetary disc is much weaker than in theprotoplanetary disc, as supported by previous studies (Fu-jii et al. 2011, 2014, KW14, Turner et al. 2014). The sharpdrop in the ionisation fraction profile at z ∼ H p occurs atthe shock boundary surrounding the circumplanetary disc,visible as a halo of poorly ionised gas in the edge-on viewshown in the centre panel of Fig. 1.Dust grains reduce the ionisation fraction by soaking upfree electrons. Grain charge capture is significant and mostefficient at the midplane where grains are abundant. It alsoreduces radioactive decay, influencing ionisation, and hencethe potential for magnetic coupling, in the circumplanetarydisc. Fig. 2 shows diffusivity profiles for the gap. The left andcentre panels show 2D slices through the gap, with thecolour scale showing the field perpendicular diffusivity η ⊥ =( η P + η H ) / , along with logarithmically spaced contour lev-els. The left panel shows the top-down view at the midplane,and the centre panel shows the edge-on profile at y = 0 . Theright panel shows the vertical profiles of the transport coef-ficients of Ohmic resistivity, the Hall effect and ambipolardiffusion, along with the coupling threshold for toroidal andMRI fields.Ohmic resistivity traces the high density structures suchas the spiral arms and circumplanetary disc and is lowestin the evacuated regions to the upper-right and lower-left of c (cid:13) Year RAS, MNRAS , 1–14 agnetic fields in protoplanetary gaps Figure 1.
Ionisation fraction profile slices through the gap, a posteriori computed for the hydrodynamical simulation of Tanigawa et al.(2012). Left and centre panels show 2D slices of the ionisation fraction log ( n e /n ) as contour plots with logarithmically spaced contoursfor y = 0 and z = 0 , respectively, and for f dg = 10 − . Right panel shows n e /n as a function of height for x = 0 . au , y = 0 for dust togas mass ratios f dg = 0 , − , − with solid, dashed, and dotted curves respectively. Figure 2.
As for Figure 1 but showing a contour plot of the field-perpendicular diffusivity, η ⊥ = (cid:113) η P + η H at the midplane (leftpanel), and y = 0 (centre panel). Ohmic resistivity (solid curve), Hall effect (dashed curve), and ambipolar diffusion (dotted curve) arealso shown as functions of height at x = 0 . au, y = 0 (right panel). The long-dashed and dot-dashed curves show the coupling thresholdsfor toroidal and MRI fields, Ω H and v a,z Ω − , respectively. the circumplanetary disc. It decreases with height in keepingwith the exponential density profile. Hall drift is relativelyconstant with height above the midplane and dominates upto z = 2 . H p . Ambipolar diffusion is strongest in the lowdensity atmosphere where neutrals are too sparse to influ-ence electron and ion motion.Non-ideal effects are orders of magnitude below the levelat which the toroidal field decouples, except in the circum-planetary disc. Here the field is uncoupled to the midplanegas flow, but is anchored to, and transported by, gas in theoverlying coupled region.Non-ideal effects are below the MRI suppression thresh-old in the disc atmosphere, z (cid:38) H p . The turbulent regionis limited by the weak field condition [equation (13)], ratherthan diffusion. Magnetic pressure is weaker than gas pres- sure at z (cid:46) H p , and the field is weak enough to bend.The strong field region reaches towards the midplane nearthe planet where field enhancement from flux-conservationis greatest. Fig. 3 shows the magnetic field strength and geometry. Thetop panels are contour plots, in the top-down (top-left panel)and edge-on (top-right panel) views, with colour scale show-ing log ( B ) , calculated for s = 1 . The bottom-left panelshows vertical profiles for the poloidal, toroidal and turbu-lent magnetic field components evaluated at x = 0 . au , y =0 . A turbulent field is only present in the active zone be-tween . (cid:46) z (cid:46) H p for s = − and for z (cid:46) H p for c (cid:13) Year RAS, MNRAS , 1–14
Sarah L. Keith and Mark Wardle
Figure 3.
Properties of the estimated magnetic field in the gap for the standard parameter set f dg = 10 − , f Σ = 1 , f B = 1 . Top-left and-right panels show contour plots of the field strength, B , at the midplane and y = 0 , respectively. Bottom-left panel shows the verticalprofile of field components at x = 0 . au, y = 0 . The vertical, toroidal, MRI field with s = 1 , and MRI field with s = − are shown as thedot-dashed, dashed, solid, and dotted curves respectively. Bottom-right panel shows the edge-on view of the field geometry, colour-codedaccording to the dominant field component for s = 1 : (a) poloidal field (yellow), (b) toroidal field (green), and (c) MRI turbulent field(blue). Orange, dark green, and dark blue regions indicate where the Hall effect sustains and influences the toroidal field. Hatched MRIunstable regions show where Hall effect dominates, so that the field is predominantly turbulent if s = 1 , and toroidal if s = − . s = 1 . The bottom-right panel is an edge-on view of the gap,colour coded according to the dominant field component.The gap is divided into the following regions: (a) poloidal(yellow), (b) toroidal (green), and (c) MRI turbulent (blue)field. Orange, dark green, and dark blue regions indicate thecorresponding regions where the Hall effect maintains cou-pling of the toroidal component, and so the field may becounter-wrapped when B z is aligned with the rotation axis.Hatched blue regions are Hall-MRI unstable and so they arepredominantly turbulent if s = 1 , or predominantly toroidalif s = − .The field strength is relatively uniform in the gap butflux conservation enhances the field in the circumplanetarydisc. A toroidal field would easily couple and permeate thegap, but is limited by magnetic pressure in the disc atmo-sphere, and so the field would be increasingly poloidal withheight. The gap would be turbulent between z ∼ . – H p ,except in the circumplanetary disc where the field is toostrong. The Hall effect extends the turbulent region to theprotoplanetary disc and circumplanetary disc midplanes if s = 1 , however the midplanes are MRI stable if s = − . Theupper boundary of the turbulent zone is limited by the weakfield requirement, rather than non-ideal effects. Fig. 4 shows the effect of varying parameters on thefield geometry: f dg = 0 , − (top- and bottom-left panels), f Σ = 0 . , (top-and bottom-centre panels), and f B = 0 . , (top- and bottom-right panels). In general, we find thatthe gap is either predominantly toroidal, ideal MRI, or Hall-MRI unstable. The upper MRI boundary only depends on β z through the column density and magnetic field strength,whereas the boundary of the Hall-MRI layer is controlled bythe height at which the Hall drift exceeds the MRI threshold, Ω v − a,z . Thus, the Hall-MRI region extends by lowering theionisation fraction, caused either through enhanced surfacedensity attenuating ionising radiation (e.g., for f Σ > ), orby grains soaking up free electrons (e.g., for f dg = 10 − ).Lowering the magnetic field strength reduces the Alfvénspeed, and consequently lowers the MRI threshold to belowthe Hall drift (e.g., for f B < ).On the other hand, the circumplanetary disc may hosta wide range of conditions. Here the Hall-MRI zone in thecircumplanetary disc is resilient to changes in f dg , but re-cedes if β z increases through enhanced density or reducefield strength, so that the field is predominantly toroidal.The field would be predominantly poloidal if Ohmic resis-tivity exceeds the Hall effect, either by enhanced resistivity c (cid:13) Year RAS, MNRAS , 1–14 agnetic fields in protoplanetary gaps through higher density, or reduced Hall EMF from a lowerfield strength. However, further analysis into the potentialfor MRI turbulence in the circumplanetary disc will rely onmodels which are not affected by numerical viscosity. In this section we consider the implications of including non-ideal effects in MHD modelling of gas flow. We determinethe extent to which diffusion limits small scale field gra-dients to calculate the minimum field bending radius. Wediscuss the implications of diffusion-limited field bending onthe resolution needed for simulations. We also estimate thestrength of large-scale magnetic forces, to compare with thehydrodynamical forces included in the TOM12 simulation.This allows us to determine where, if at all, magnetic forceshave the most influence on gas flow.
Non-ideal effects, particularly Ohmic resistivity and ambipo-lar diffusion, can wash out magnetic field gradients and limitthe field bending length-scale, L B . Here we calculate thesmallest length-scale that the field can bend on, L min , giventhe level of magnetic resistivity and diffusion throughout thegap as calculated in Section 3.2. The role of the Hall effecton magnetic field gradients is highly uncertain, being ableto enhance or resist gradients. Therefore we neglect its effectin this simple calculation, but note the potential key impactit can have on the large-scale structure we consider here.Just as equations (12), (14)–(15) specify the maximumdiffusivity which permits the gradients needed for a tur-bulent or toroidal field, we can invert this relationship tospecify the minimum bending length-scale given the mag-netic diffusivity. This is applicable to large-scale, quasi-staticmagnetic field structures (as opposed to rapidly fluctuatingturbulence), for which magnetic induction is balanced bymagnetic diffusion. We derive the relation between magneticdiffusion and a general field geometry using the inductionequation [equation (11)]. While it is generally necessary totreat each component of the induction equation separately,we can estimate the diffusion limit for the minimum gradientlength-scale, η P < vL B (cid:18) L − B + L − v L − B + L − η (cid:19) . (27)We calculate the velocity and diffusivity gradient length-scales, L v , and L η , using L f = f/ |∇ f | , (28)with f set to the relevant fluid velocity or Pedersen diffu-sivity, respectively. Note that as we use the velocity in theframe orbiting with the planet to calculate L v , we removethe Keplerian component with respect to the star, whichacts as a large constant offset in equation (28).We invert equation (27) to determine the minimum fieldgradient length-scale L min : L min = 2 η P v − ηvL η + (cid:115)(cid:18) − ηvL η (cid:19) + 4 ηvL v − . (29) Figure 5 shows a contour plot of log ( L min / au ) at y = 0 , with logarithmically spaced contour level. The mini-mum field bending length-scale follows L min ∝ η P /v every-where except the circumplanetary disc where the η P is largeand L v , L η are small. It increases toward the planet, andis smallest at z ∼ H p in the protoplanetary disc where thediffusivity is lowest.This can be used to gauge the minimum resolution re-quired for MHD simulations. For example, as the minimumfield bending length scale is large in the circumplanetarydisc, the resolution is determined by velocity and pressuregradients length-scales instead. In the protoplanetary disc,significantly higher resolution is needed to resolve turbu-lence which can be tangled on a very small length scales.The Hall effect can pose an additional challenge to simula-tions as it is dissipationless and so can introduce small-scalefield structure (e.g. Sano & Stone 2002, Kunz & Lesur 2013;Lesur et al. 2014; O’Keeffe & Downes 2014; Bai 2015). Large-scale magnetic fields have the potential to influencegas dynamics, rather than being merely passively drawnalong with the gas. We touched on this in Section 2.4when we compared the toroidal component of magnetic pres-sure with gas pressure. For example, circumplanetary jetshave been launched in ideal and resistive MHD simulations(Machida et al. 2006; Gressel et al. 2013).The TOM12 simulation is not magnetised; neverthelesswe can estimate the size of magnetic forces and comparethem with hydrodynamic forces in the simulation. This willindicate where, if at all, they have the greatest potential toalter the balance of forces and so influence the gas flow. Italso acts as a consistency check, validating (or otherwise) ouradoption of a hydrodynamic simulation as the underlyingmodel.The full momentum equation, including both hydrody-namic and magnetic forces is ρ ∂ v ∂t + ( v · ∇ ) v = − p ˆz × v − Ω p ˆz × ( ˆz × r ) − ∇ Φ − ρ ∇ p + 14 πρ ( ∇ × B ) × B , (30)where r = ( x, y, z ) and Φ is the combined gravitational po-tential from the star and planet: Φ = − GM ∗ | r − r ∗ | − GM p | r | . (31)The hydrodynamic terms in equation (30) are the timederivative and inertial force on the left hand side, and theCoriolis, centrifugal force, gravitational, and pressure gradi-ent forces on the right hand side. These are calculated in theframe co-orbiting with the protoplanet in the TOM12 sim-ulation. We neglect the time derivative because it is smallwith the simulation in almost steady state.The addition of the magnetic force [the final term onthe right hand side of equation (30)] will inevitably changethe balance of forces, however the impact will be negligibleunless the magnetic force is strong compared to the inertialforce. We calculated the inertial force F I = | ( v · ∇ ) v | fromthe TOM12 simulation and its associated gradient length c (cid:13) Year RAS, MNRAS , 1–14 Sarah L. Keith and Mark Wardle
Figure 4.
Same as bottom-right panel in Fig. 3, except varying the dust-to-gas mass ratio f dg = 0 , − (top- and bottom-left panels),column density enhancement factor f Σ = 0 . , (top-and bottom-centre panels), and magnetic field enhancement factor f B = 0 . , (top- and bottom-right panels). Figure 5.
Left panel shows contours of the minimum field bending-scale, L min , shown in the edge-on view at y = 0 . Centre and rightpanels shows the strength of magnetic forces relative to hydrodynamical forces. The centre panel shows the logarithm of the ratio of themagnetic force to the inertial force, log ( F B /F I ) at y = 0 . The dashed contour shows the location of F B = F I . The right panel showsvertical profiles of log ( F B /F I ) and log ( F B /F P ) , the logarithm of the ratio of the magnetic force to the inertial and pressure gradientforces at x = 0 . au, y = 0 as the solid and dashed curves, respectively. c (cid:13) Year RAS, MNRAS000
Left panel shows contours of the minimum field bending-scale, L min , shown in the edge-on view at y = 0 . Centre and rightpanels shows the strength of magnetic forces relative to hydrodynamical forces. The centre panel shows the logarithm of the ratio of themagnetic force to the inertial force, log ( F B /F I ) at y = 0 . The dashed contour shows the location of F B = F I . The right panel showsvertical profiles of log ( F B /F I ) and log ( F B /F P ) , the logarithm of the ratio of the magnetic force to the inertial and pressure gradientforces at x = 0 . au, y = 0 as the solid and dashed curves, respectively. c (cid:13) Year RAS, MNRAS000 , 1–14 agnetic fields in protoplanetary gaps scale L I = v | ( v · ∇ ) v | , (32)using second-order finite differencing, with good convergencecompared with first-order differencing.We estimate the magnetic force from large-scale non-turbulent fields using, F B ≡ πρ | ( ∇ × B ) × B | ∼ v a L B . (33)We do this by using the Alfvén speed calculated for thetotal, non-turbulent field strength. Our magnetic force esti-mate relies on a knowledge of the field gradient length-scale, L B . We expect that, in general, L B traces the fluid gradientlength-scales because either the field is dragged along withthe fluid or visa versa. If this is the case, we take L B = L I ,and it cancels from the ratio of the magnetic and hydrody-namic forces F B /F I . If, instead, the gas varies over a shorterlength-scale than diffusion permits the field to bend, we take L B = L min as calculated in Section 4.1. In this way, we com-pare the magnetic force over the same scale as the hydrody-namic forces, to the extent that diffusion permits. We havenot included the the force arising from effective MRI turbu-lent viscosity, as it is already the subject of extensive studywithin the literature (see Turner et al. 2014).We also consider the effect magnetic forces have on thescale height by comparing them with the pressure gradientforce, F P = | ( ∇ p ) /ρ | , with length-scale L P calculated usingequation (28).The centre-panel of Fig. 5 shows a contour plot of thelogarithm of the ratio of the magnetic force to the inertialforce, log ( F B /F I ) at y = 0 . The dashed contour denotes F I = F B . The right panel shows the vertical profile of theratio of the magnetic force to the inertial, and pressure gra-dient forces. Magnetic forces are strongest in the disc atmo-sphere, exceeding the inertial and pressure gradient forcesabove z (cid:38) . – . H p , where they are able to influence theprotoplanetary disc scale-height.Magnetic forces are weakest in the circumplanetary discwhere strong magnetic diffusion limits the field bendinglength scale to the diffusion length scale. Although artificialdiffusion lessens the reliability of the circumplanetary discstructure in the TOM12 simulation, KW14 shows that mag-netic diffusion is strong across a range of circumplanetarydisc models, owing to the column high density shielding thedisc from external ionising radiation. The strong diffusionincreases the minimum field bending length scale, therebyreducing the magnetic force. This indicates that large-scalemagnetic forces cannot yield any significant accretion in thecircumplanetary disc, which must come from the MRI orGravitational Instability-MRI limit-cycles (e.g., Martin &Lubow 2011; Lubow & Martin 2012). In this paper we examined the importance of non-ideal ef-fects in determining the magnetic field structure in a gapsurrounding a giant protoplanet. We modelled the gap us-ing a snapshot from the pure-hydrodynamical gap-crossingsimulation by Tanigawa et al. (2012). Our approach was to use this snapshot as a basis to a posteriori estimate key MHDquantities semi-analytically, which would otherwise be verychallenging to incorporate into simulations. We calculatedthe ionisation fraction produced by cosmic-rays, stellar X-rays and radioactive decay, including the effect of grains. Wecalculated Ohmic resistivity, ambipolar diffusion and Halldrift to determine whether an MRI field could be generatedand if a toroidal field could couple to the gas flow. We es-timated the magnetic field strength in the protoplanetarydisc from inferred accretion rates, and determined the gapfield strength from flux-freezing.We found that a magnetic field would be easily drawnfrom the protoplanetary disc into the gap and circumplane-tary disc. A toroidal field permeates the gap, but is weakenedby expansion from magnetic pressure above z > –3 H p . Thegap is MRI unstable at z > . H p with turbulence extend-ing down to the midplane if the vertical component of themagnetic field is along the rotation axis (i.e., if s = 1 ). Ifhowever, the vertical field and rotation axis are anti-aligned,the conditionally unstable regions are non-turbulent. As pro-toplanetary discs exhibit a range of field/rotation axis mis-alignment angles we expect the size of the MRI turbulentregion to differ between systems, with significant implica-tions for the flow dynamics (Hull et al. 2013; Krumholz et al.2013).This direction dependence of the MRI originates withthe Hall effect. The Hall effect is strong below z < H p andplays an important role in generating MRI instability. Wefound that it can also facilitate coupling of a toroidal field,despite strong magnetic Pedersen diffusion, and influencesthe orientation of the toroidal component. This may lead tocounter-wrapping of the toroidal field if s = 1 , where thefield drift opposes the Keplerian gas flow. Further simula-tions are needed to probe the role of the Hall effect in thissystem.By testing the sensitivity of the calculations to dust con-tent, column density, and magnetic field strength, we foundthat the gap was generally susceptible to the Hall- or ideal-MRI. On the other hand, the MRI in the circumplanetarydisc may be quite sensitive to disc conditions and could havea turbulent, toroidal or vertical field. As the circumplanetarydisc evolves, the disc may experience a range of different fieldconfigurations in keeping with the varying column densityand dust content. Understanding the evolution of a circum-planetary disc is key for satellite formation studies, whichis believed to be the formation site of moons. Bimodalityin circumplanetary disc dynamics caused by the Hall effectmay transfer to the growth of moons within the system. Forexample, turbulent heating will effect the location of thecircumplanetary ice-line.Finally, we have calculated the minimum magnetic fieldgradient length-scale limited by magnetic diffusion. Mag-netic diffusion resists field gradients, and so minimal fieldbending is permitted in the circumplanetary disc whereOhmic resistivity is strong. We find that large-scale (non-turbulent) magnetic forces are unable to drive accretion incircumplanetary discs. Turbulence aside, we found the largescale-flow features are well modelled by a hydrodynamicalfluid as large scale magnetic forces as small outside theprotoplanetary disc atmosphere. Here, magnetic coupling isstrong and so a strong field may be able to produce a jet,winds, or other variability. c (cid:13) Year RAS, MNRAS , 1–14 Sarah L. Keith and Mark Wardle
In summary, we find that the magnetic field surround-ing a giant protoplanet is mostly toroidal, with large bandsof ideal- and Hall-MRI turbulent zones. The magnetic fieldgeometry is dependent on the orientation of the vertical fieldcomponent, established during the collapse of the protostel-lar core. Non-ideal effects are important in the gap and needto be included in future MHD simulations of gap crossing.
ACKNOWLEDGEMENTS
We are grateful to Takayuki Tanigawa for the generous pro-vision of the simulation data, along with assistance in us-ing and interpreting the simulation. We also thank Pandey,and Oliver Gressel for valuable discussions, and PhilippaBrowning, Michael Keith, and the anonymous referee forhelpful comments on the manuscript. The authors thank theNiels Bohr International Academy for hospitality. This workwas supported in part by the Australian Research Coun-cil grant DP130104873. S.K. further acknowledges the sup-port of an Australian Postgraduate Award, funding from theMacquarie University Postgraduate Research Fund schemeand the International Space Science Institute. This researchhas made use of NASA’s Astrophysical Data System.
REFERENCES
Andrews S. M., Williams J. P., 2005, ApJ, 631, 1134Artymowicz P., Lubow S. H., 1996, ApJ, 467, L77Asplund M., Grevesse N., Sauval A. J., Scott P., 2009,ARA&A, 47, 481Ayliffe B. A., Bate M. R., 2009a, MNRAS, 397, 657Ayliffe B. A., Bate M. R., 2009b, MNRAS, 393, 49Ayliffe B. A., Bate M. R., 2012, MNRAS, 427, 2597Bai X.-N., 2015, ApJ, 798, 84Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Balbus S. A., Terquem C., 2001, ApJ, 552, 235Bate M. R., Lubow S. H., Ogilvie G. I., Miller K. A., 2003,MNRAS, 341, 213Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883Braiding C. R., Wardle M., 2012, MNRAS, 427, 3188Bryden G., Chen X., Lin D. N. C., Nelson R. P., PapaloizouJ. C. B., 1999, ApJ, 514, 344Casassus S., van der Plas G., M S. P., Dent W. R. F.,Fomalont E., Hagelberg J., Hales A., Jordán A., MawetD., Ménard F., Wootten A., Wilner D., Hughes e. a., 2013,Nature, 493, 191Cassen P., Moosman A., 1981, Icarus, 48, 353D’Angelo G., Kley W., Henning T., 2003, ApJ, 586, 540Debes J. H., Jang-Condell H., Weinberger A. J., RobergeA., Schneider G., 2013, ApJ, 771, 45Dunhill A. C., 2015, MNRAS, 448, L67Flock M., Dzyurkevich N., Klahr H., Turner N., HenningT., 2012, ApJ, 744, 144Flock M., Ruge J. P., Dzyurkevich N., Henning T., KlahrH., Wolf S., 2015, A&A, 574, A68Fujii Y. I., Okuzumi S., Inutsuka S.-I., 2011, ApJ, 743, 53Fujii Y. I., Okuzumi S., Tanigawa T., Inutsuka S.-I., 2014,ApJ, 785, 101Galassi M., Davies J., Theiler J., Gough B., Jungman G., Alken P., Booth M., Rossi F., 2009, GNU Scientific Li-brary Reference Manual, third edn. Network Theory Ltd.Garufi A., Quanz S. P., Schmid H. M., Avenhaus H., Buen-zli E., Wolf S., 2014, A&A, 568, A40Gressel O., Nelson R. P., Turner N. J., Ziegler U., 2013,ApJ, 779, 59Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015,ApJ, 801, 84Guilloteau S., Dutrey A., 1998, A&A, 339, 467Hartmann L., Calvet N., Gullbring E., D’Alessio P., 1998,ApJ, 495, 385Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440,742Hayashi C., 1981, Progress of Theoretical Physics Supple-ment, 70, 35Hirose S., Turner N. J., 2011, ApJ, 732, L30Hull C. L. H., Plambeck R. L., Bolatto A., et al., 2013,ApJ, 768, 159Igea J., Glassgold A. E., 1999, ApJ, 518, 848Isella A., Chandler C. J., Carpenter J. M., Pérez L. M.,Ricci L., 2014, ApJ, 788, 129Keith S. L., Wardle M., 2014, MNRAS, 440, 89Kitamura Y., Momose M., Yokogawa S., Kawabe R.,Tamura M., Ida S., 2002, ApJ, 581, 357Krumholz M. R., Crutcher R. M., Hull C. L. H., 2013, ApJ,767, L11Kunz M. W., Lesur G., 2013, MNRAS, 434, 2295Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56Lide D. R., 2004, CRC Handbook of Chemistry andPhysics, 85th Edition, 85 edn. CRC PressLin D. N. C., Papaloizou J., 1985, in Black D. C., MatthewsM. S., eds, Protostars and Planets II On the dynamicalorigin of the solar system. pp 981–1072Lubow S. H., D’Angelo G., 2006, ApJ, 641, 526Lubow S. H., Martin R. G., 2012, ApJ, 749, L37Lubow S. H., Seibert M., Artymowicz P., 1999, ApJ, 526,1001Lunine J. I., Stevenson D. J., 1982, Icarus, 52, 14Machida M. N., Inutsuka S.-I., Matsumoto T., 2006, ApJ,649, L129Martin R. G., Lubow S. H., 2011, ApJ, 740, L6Matsumoto T., Tomisaka K., 2004, ApJ, 616, 266McElroy D., Walsh C., Markwick A. J., Cordiner M. A.,Smith K., Millar T. J., 2013, A&A, 550, A36McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565Miller K. A., Stone J. M., 2000, ApJ, 534, 398Muranushi T., Okuzumi S., Inutsuka S.-I., 2012, ApJ, 760,56Nelson R. P., Papaloizou J. C. B., 2003, MNRAS, 339, 993O’Keeffe W., Downes T. P., 2014, MNRAS, 441, 571Okuzumi S., Hirose S., 2011, ApJ, 742, 65Okuzumi S., Inutsuka S.-i., 2015, ApJ, 800, 47Papaloizou J. C. B., Nelson R. P., Snellgrove M. D., 2004,MNRAS, 350, 829Piétu V., Dutrey A., Guilloteau S., 2007, A&A, 467, 163Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P.,Roush T., Fong W., 1994, ApJ, 421, 615Ricci L., Testi L., Natta A., Scholz A., de Gregorio-Monsalvo I., Isella A., 2014, ApJ, 791, 20Rivier G., Crida A., Morbidelli A., Brouet Y., 2012, A&A,548, A116 c (cid:13) Year RAS, MNRAS000
Andrews S. M., Williams J. P., 2005, ApJ, 631, 1134Artymowicz P., Lubow S. H., 1996, ApJ, 467, L77Asplund M., Grevesse N., Sauval A. J., Scott P., 2009,ARA&A, 47, 481Ayliffe B. A., Bate M. R., 2009a, MNRAS, 397, 657Ayliffe B. A., Bate M. R., 2009b, MNRAS, 393, 49Ayliffe B. A., Bate M. R., 2012, MNRAS, 427, 2597Bai X.-N., 2015, ApJ, 798, 84Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Balbus S. A., Terquem C., 2001, ApJ, 552, 235Bate M. R., Lubow S. H., Ogilvie G. I., Miller K. A., 2003,MNRAS, 341, 213Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883Braiding C. R., Wardle M., 2012, MNRAS, 427, 3188Bryden G., Chen X., Lin D. N. C., Nelson R. P., PapaloizouJ. C. B., 1999, ApJ, 514, 344Casassus S., van der Plas G., M S. P., Dent W. R. F.,Fomalont E., Hagelberg J., Hales A., Jordán A., MawetD., Ménard F., Wootten A., Wilner D., Hughes e. a., 2013,Nature, 493, 191Cassen P., Moosman A., 1981, Icarus, 48, 353D’Angelo G., Kley W., Henning T., 2003, ApJ, 586, 540Debes J. H., Jang-Condell H., Weinberger A. J., RobergeA., Schneider G., 2013, ApJ, 771, 45Dunhill A. C., 2015, MNRAS, 448, L67Flock M., Dzyurkevich N., Klahr H., Turner N., HenningT., 2012, ApJ, 744, 144Flock M., Ruge J. P., Dzyurkevich N., Henning T., KlahrH., Wolf S., 2015, A&A, 574, A68Fujii Y. I., Okuzumi S., Inutsuka S.-I., 2011, ApJ, 743, 53Fujii Y. I., Okuzumi S., Tanigawa T., Inutsuka S.-I., 2014,ApJ, 785, 101Galassi M., Davies J., Theiler J., Gough B., Jungman G., Alken P., Booth M., Rossi F., 2009, GNU Scientific Li-brary Reference Manual, third edn. Network Theory Ltd.Garufi A., Quanz S. P., Schmid H. M., Avenhaus H., Buen-zli E., Wolf S., 2014, A&A, 568, A40Gressel O., Nelson R. P., Turner N. J., Ziegler U., 2013,ApJ, 779, 59Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015,ApJ, 801, 84Guilloteau S., Dutrey A., 1998, A&A, 339, 467Hartmann L., Calvet N., Gullbring E., D’Alessio P., 1998,ApJ, 495, 385Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440,742Hayashi C., 1981, Progress of Theoretical Physics Supple-ment, 70, 35Hirose S., Turner N. J., 2011, ApJ, 732, L30Hull C. L. H., Plambeck R. L., Bolatto A., et al., 2013,ApJ, 768, 159Igea J., Glassgold A. E., 1999, ApJ, 518, 848Isella A., Chandler C. J., Carpenter J. M., Pérez L. M.,Ricci L., 2014, ApJ, 788, 129Keith S. L., Wardle M., 2014, MNRAS, 440, 89Kitamura Y., Momose M., Yokogawa S., Kawabe R.,Tamura M., Ida S., 2002, ApJ, 581, 357Krumholz M. R., Crutcher R. M., Hull C. L. H., 2013, ApJ,767, L11Kunz M. W., Lesur G., 2013, MNRAS, 434, 2295Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56Lide D. R., 2004, CRC Handbook of Chemistry andPhysics, 85th Edition, 85 edn. CRC PressLin D. N. C., Papaloizou J., 1985, in Black D. C., MatthewsM. S., eds, Protostars and Planets II On the dynamicalorigin of the solar system. pp 981–1072Lubow S. H., D’Angelo G., 2006, ApJ, 641, 526Lubow S. H., Martin R. G., 2012, ApJ, 749, L37Lubow S. H., Seibert M., Artymowicz P., 1999, ApJ, 526,1001Lunine J. I., Stevenson D. J., 1982, Icarus, 52, 14Machida M. N., Inutsuka S.-I., Matsumoto T., 2006, ApJ,649, L129Martin R. G., Lubow S. H., 2011, ApJ, 740, L6Matsumoto T., Tomisaka K., 2004, ApJ, 616, 266McElroy D., Walsh C., Markwick A. J., Cordiner M. A.,Smith K., Millar T. J., 2013, A&A, 550, A36McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565Miller K. A., Stone J. M., 2000, ApJ, 534, 398Muranushi T., Okuzumi S., Inutsuka S.-I., 2012, ApJ, 760,56Nelson R. P., Papaloizou J. C. B., 2003, MNRAS, 339, 993O’Keeffe W., Downes T. P., 2014, MNRAS, 441, 571Okuzumi S., Hirose S., 2011, ApJ, 742, 65Okuzumi S., Inutsuka S.-i., 2015, ApJ, 800, 47Papaloizou J. C. B., Nelson R. P., Snellgrove M. D., 2004,MNRAS, 350, 829Piétu V., Dutrey A., Guilloteau S., 2007, A&A, 467, 163Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P.,Roush T., Fong W., 1994, ApJ, 421, 615Ricci L., Testi L., Natta A., Scholz A., de Gregorio-Monsalvo I., Isella A., 2014, ApJ, 791, 20Rivier G., Crida A., Morbidelli A., Brouet Y., 2012, A&A,548, A116 c (cid:13) Year RAS, MNRAS000 , 1–14 agnetic fields in protoplanetary gaps Sano T., Inutsuka S.-I., Turner N. J., Stone J. M., 2004,ApJ, 605, 321Sano T., Stone J. M., 2002, ApJ, 577, 534Shi J., Krolik J. H., Hirose S., 2010, ApJ, 708, 1716Shu F. H., Adams F. C., Lizano S., 1987, ARA&A, 25, 23Stephens I. W., Looney L. W., Kwon W., Fernández-LópezM., Hughes A. M., Mundy L. G., Crutcher R. M., Li Z.-Y.,Rao R., 2014, Nature, 514, 597Szulágyi J., Morbidelli A., Crida A., Masset F., 2014, ApJ,782, 65Tanigawa T., Ohtsuki K., Machida M. N., 2012, ApJ, 747,47Turner N. J., Fromang S., Gammie C., Klahr H., Lesur G.,Wardle M., Bai X.-N., 2014, Protostars and Planets VI,pp 411–432Turner N. J., Lee M. H., Sano T., 2014, ApJ, 783, 14Turner N. J., Sano T., 2008, ApJ, 679, L131Umebayashi T., Nakano T., 1980, PASJ, 32, 405Umebayashi T., Nakano T., 1981, PASJ, 33, 617Umebayashi T., Nakano T., 2009, ApJ, 690, 69Uribe A. L., Klahr H., Henning T., 2013, ApJ, 769, 97Wardle M., 1999, MNRAS, 307, 849Wardle M., 2007, Ap&SS, 311, 35Wardle M., Königl A., 1993, ApJ, 410, 218Wardle M., Salmeron R., 2012, MNRAS, 422, 2737Weidenschilling S. J., 1977, Ap&SS, 51, 153Williams J. P., Cieza L. A., 2011, ARA&A, 49, 67Zhu Z., 2015, ApJ, 799, 16
APPENDIX A: WEAK-FIELD CONDITIONFOR HALL-MRI
Developing and sustaining MRI turbulence requires two con-ditions: (i) firstly, that field perturbations must grow fastenough, and (ii) secondly, the wavelength of the fastest-growing turbulent mode λ , must be contained within adisc scale-height (i.e., λ ≤ H ). Here we develop the weak-field condition for Hall-MRI. The wavelength of the fastest-growing Hall-MRI mode is [see equation (B14) of Wardle &Salmeron 2012] λ = πν (cid:20) sη H Ω − η P ν − v a,z + 10 v a,z Ω ν + Ω (cid:21) , (A1)where ν is the growth rate of the fastest-growing mode. Inthe Hall MHD limit, η P = 0 , the maximum growth rateattains the ideal rate ν = Ω − (Wardle & Salmeron 2012).Using this result in strong-Hall limit, η H Ω v − a,z (cid:29) , thewavelength of the fastest-growing mode is approximately λ ≈ π (cid:114) η H Ω , (A2)up to a constant factor of order unity.Ensuring that the two MRI conditions [(i), (ii) above]are met allows us to bracket the Hall-drift: v a,z Ω ≤ η H ≤ Ω (cid:18) H π (cid:19) , (A3)which, with the thin-disc approximation H = c s / Ω , leads tothe weak-field limit β z = 2 c s v a,z ≥ π . (A4) This is identical to that in the ideal and resistive MRIregimes [equation (13); Okuzumi & Hirose 2011]. APPENDIX B: NON-IDEAL EFFECTS ACTINGON A TOROIDAL FIELD
Here we determine the coupling threshold for a toroidalfield using the azimuthal component of the induction equa-tion. We consider a Keplerian disc with velocity, v = Ω r ˆ φ ,independent of height, and uniform transport coefficients η O , η H , η A for simplicity.Shear generates a toroidal component from a poloidalfield in the inductive term: [ ∇ × ( v × B )] φ = − Ω B r . Theresistive term, [ −∇ × ( η O ∇ × B )] φ = η O ( ∇ B φ − B φ /r ) ∼ η O B φ /H , is dominated by vertical gradients (see Section2.4; Wardle & Königl 1993), as are the Hall and ambipolarterms, which scale similarly. Comparing the two terms yieldsthe coupling criterion η < Ω H , (B1)where η = η P , η H .How will the field evolve once it decouples from the gasmotion? Recasting the induction equation to show the fieldline drift, V B , makes this clear (Wardle & Salmeron 2012): ∂ B ∂t = ∇ × (cid:110) ( v + V B ) × B − η O (cid:104) ( ∇ × B ) · ˆ B (cid:105) ˆ B (cid:111) , (B2)where V B = η P B ( ∇ × B ) ⊥ × ˆ B − η H B ( ∇ × B ) ⊥ (B3)The azimuthal component of the drift velocity, V B,φ , con-trols the wrapping of the field lines with the shear. Neglect-ing terms of O ( H/r ) , the azimuthal drift velocity is (Braid-ing & Wardle 2012): V B,φ = η P H B z B φ,s B − η H H B r,s B (B4) ≡ V P + V H , (B5)where we have separated the Pedersen and Hall drift com-ponents. Pedersen drift always resists winding, tending toreduce | B φ,s | and straighten field lines. On the other hand,the Hall component will enhance or resist winding accordingthe sign of B r,s .In the natural configuration for s = 1 (i.e., B r,s > ,and B φ,s < for z > ), V H < and V P < so that thecomponents cooperate in unwrapping the field. For domi-nant Hall, the unwrapping will overshoot, winding the fieldin the opposite direction so that B φ,s > . For the anti-aligned field, s = − , the natural configuration is B r,s < ,and B φ,s > above the midplane. Once again, Pedersen re-sists wrapping of the field lines and if it dominates the fieldwill tend towards vertical. On the other hand, now Hall drifttends to enhance field motion along ˆ v φ so that even if thefield is somehow wound against the flow (i.e., B φ,s > ) bothcomponents cooperate to restore B φ,s ≥ in equilibrium.This behaviour is summarised in Fig. B1, which shows This threshold is a factor of
10 (
H/r ) − ∼ × lower thanTurner & Sano (2008), as we allow for vertical gradients in thenon-ideal terms.c (cid:13) Year RAS, MNRAS , 1–14 Sarah L. Keith and Mark Wardle
Hall-modified ! toroidal field Diffusive ! vertical field ~Ideal-MHD ! toroidal field ⌘ H / (⌦ H ) ⌘ P / ( ⌦ H ) Figure B1.
Regions of η H – η P parameter space, shaded accord-ing to the toroidal field drift. Ideal MHD (light green) and Halldominated MHD (dark green) will sustain a toroidal field, al-though the orientation of the field varies with the vertical fieldorientation in the Hall-MHD region. Pedersen diffusion only per-mits a vertical field in the yellow region. how the resulting field geometry varies in the η H – η P pa-rameter space. Non-ideal effects are weak in the lower-leftquadrant (light green), with the toroidal field wound up bythe disc so that B z B φ < . Pedersen diffusion dominatesabove the diagonal (yellow region), so that the field tendstoward vertical. The Hall effect counteracts diffusion belowthe diagonal (dark green), so that the toroidal component isenhanced for s = − or counter-wrapped when s = 1 .This paper has been typeset from a TEX/ L A TEX file preparedby the author. c (cid:13) Year RAS, MNRAS000