arXiv:1308.4139v1 [astro-ph.SR] 19 Aug 2013 To appear in ApJ Reexamination of Induction Heating of Primitive Bodies in Protoplanetary Disks Raymond L. Menzel1 and Wayne G. Roberge1 New York Center for Astrobiology and Department of Physics, Applied Physics and Astronomy, Rensselaer Polytechnic Institute, 110 8th Street, Troy, NY 12180, USA menzer@rpi.edu, roberw@rpi.edu ABSTRACT We reexamine the unipolar induction mechanism for heating asteroids originally proposed in a classic series of papers by Sonett and collaborators. As originally conceived, induction heating is caused by the “motional electric field” which appears in the frame of an asteroid immersed in a fully-ionized, magnetized solar wind and drives currents through its interior. However we point out that classical induction heating contains a subtle conceptual error, in consequence of which the electric field inside the asteroid was calculated incorrectly. The problem is that the motional electric field used by Sonett et al. is the electric field in the freely streaming plasma far from the asteroid; in fact the motional field vanishes at the asteroid surface for realistic assumptions about the plasma density. In this paper we revisit and improve the induction heating scenario by: (1) correcting the conceptual error by self consistently calculating the electric field in and around the boundary layer at the asteroid-plasma interface; (2) considering weakly-ionized plasmas consistent with current ideas about protoplanetary disks; and (3) considering more realistic scenarios which do not require a fully ionized, powerful T Tauri wind in the disk midplane. We present exemplary solutions for two highly idealized flows which show that the interior electric field can either vanish or be comparable to the fields predicted by classical induction depending on the flow geometry. We term the heating driven by these flows “electrodynamic heating”, calculate its upper limits, and compare them to heating produced by short-lived radionuclides. Subject headings: astrobiology — MHD — minor planets, asteroids: general — protoplanetary disks – 2 – 1. Introduction The mineralogy of asteroids inferred from studies of meteorites and asteroid spectroscopy implies that the asteroids were heated during the first few millions years of their lifetimes (e.g., Ghosh et al. 2006 and references therein). The degree of heating depended strongly on heliocentric distance. Igneous-type asteroids at the inner edge of the asteroid belt show signs of heating to > 1200 K and are thought to be pieces of the exposed metallic cores of differentiated asteroids (Keil 2000). At the outer edge of the asteroid belt no heating occurred and the asteroids there are thought to be composed of unaltered, primitive solar system materials (McKinnon 1989). Near the center of the asteroid belt, metamorphic type asteroids show mineralogical evidence of aqueous alteration at temperatures of 300– 450 K and thermal metamorphism at temperatures up to 1200 K (Keil 2000 and references therein). Some of the asteroids in this region produce fragments which reach Earth in the form of carbonaceous chondrites. The latter may contain innumerable organic compounds, including small amounts of amino acids, the distribution of which depends on the severity of thermal alteration inside the parent body (Glavin et al. 2011). In addition, excesses of L-type amino acids have been found in aqueously altered meteorites (Cronin & Pizzarello 1997; Glavin & Dworkin 2009). Because asteroids may have provided significant amounts of organic material to the early Earth (Chyba & Sagan 1992), understanding how asteroids were heated may be central to understanding the origin of life. Two theories have been proposed to explain asteroid heating: the decay of short-lived radionuclides (SLRs, Urey 1955) and unipolar induction heating (Sonett et al. 1970). Heating by SLRs such as 26Al is a widely accepted scenario. The discovery of an excess abundance of the daughter nuclide 26Mg in calcium-aluminum rich inclusions (CAIs) from the Allende meteorite implies an 26Al abundance large enough to cause significant heating or even melting (Wasserburg et al. 1977; Lee et al. 1977). However the gradient of heating across the asteroid belt requires a finely tuned dependence of asteroid accretion times on heliocentric distance. For example Grimm & McSween (1993) were able to reproduce the observed gradient by assuming that all bodies at each specific heliocentric distance R inside the asteroid belt accreted instantaneously, regardless of size, at a time τac relative to CAI formation determined by the equation τac ∝ R n , (1) where n is an adjustable parameter ranging from 1.5 to 3. Other models have attempted to relax this assumption; however the later models predict that most of the mass of the asteroid belt would be contained in small bodies that would not achieve melting or thermal metamorphism (McSween et al. 2002). This would imply an unobserved excess of small, unheated asteroids in all parts of the asteroid belt. Furthermore, although 26Al heating may – 3 – have been important in our solar system, the SLR scenario is unlikely to occur elsewhere: Ouellette et al. (2010) predict that the probability of another protoplanetary disk receiving the same concentration of SLRs is only ∼ 10−3–10−2 . For these reasons we revisit the unipolar induction heating mechanism described in Sonett et al. (1970). In their scenario an asteroid or similarly unmagnetized body is immersed in a uniform, fully-ionized solar wind with velocity v0, magnetic field B = B0, and electric field E = 0 in the wind’s rest frame. Sonett et al. observed that in the frame of the asteroid there appears a motional electric field, Em = − v0 c × B0 ≡ E0. (2) If the solid body is not a perfect insulator, a nonzero electric field inside it will drive currents which will generate heat via Ohmic dissipation. Sonett et al. calculated the electric field inside the body by treating the latter as a dielectric sphere immersed in a uniform electric field E0. They assumed that the interior currents are able to form a closed circuit through the plasma and estimated the distortion of the ambient magnetic field B0 by secondary magnetic fields generated by currents in- and around the body. However no attempt was made to account for distortions in the velocity field, v, which was assumed to be v0 everywhere. Unfortunately the classical induction scenario is based on a subtle misconception. Expression (2) is a local relation, in the sense that Em is the electric field, as measured in the body frame, at the location where the velocity is equal to v0. Consequently, E0 is the body frame electric field in the freely streaming plasma. If the velocity v depends on position, x, then the motional field is also position dependent. Thus, Em(x) = − v(x) c ×B0 (3) is the motional field at position x. Strictly speaking, unipolar induction describes the heating of bodies that do not perturb the velocity field of the surrounding plasma. This is a good approximation in the solar system today,1 where the collision mean free path exceeds the size of any solid body by many orders of magnitude. However the collision mean free path was of order meters in the solar nebula so the distortion of the velocity field cannot be ignored in calculations of asteroid heating. In general, friction between the plasma and body surface will cause the formation of a shear layer, in which v decreases systematically to zero at the body surface. It follows that the motional electric field also vanishes at the body surface. 1The unipolar induction mechanism was originally conceived to explain the magnetic fields of the Moon and other bodies in the solar system today (Sonett & Colburn 1968; Schwartz et al. 1969) and later applied to asteroids. – 4 – In the following sections we show that, while the motional electric field vanishes on and inside a large body, the total electric field may not. In Section 2 we give the equations which govern the velocity and electromagnetic fields in and around a body immersed in a flowing plasma. We adopt a multifluid, magnetohydrodynamic description appropriate for weakly ionized protoplanetary disks. In this paper we make no assumptions about the origin of the flow, which could be due to the body’s orbital motion (Weidenschilling 1977a; Morris et al. 2012) or passing shock waves (e.g., Desch & Connolly 2002; Miura & Nakamoto 2006; Hood et al. 2009). Section 3 describes our disk model and the method used to calculate the abundances of charged particles and transport coefficients in the disk midplane. In Section 4 we present two simple examples which show that, although the motional electric field vanishes at the body surface, a non-zero total electric field can exist due to magnetic field gradients in the flow set up by Ohmic dissipation, the Hall effect, and ambipolar diffusion. The effects of small dust grains and possible relevance of our results to asteroid heating are discussed in Section 5 and our results are summarized in Section 6. 2. Multifluid, Magnetohydrodynamic Flow Past an Arbitrary Body 2.1. Governing Equations for the Shear Flow Consider the flow of plasma past a large, unmagnetized body in a weakly ionized protoplanetary disk. Close to the body surface the plasma is distorted by shearing motions. The dynamics of the shear layer are described by the equations for mass conservation, ∂ρ ∂t + ∇·(ρv) = 0, (4) momentum conservation, ρ  ∂v ∂t + (v·∇) v  = − ∇  B2 8π  + 1 4π (B·∇) B − ∇P + α∇2v, (5) and energy conservation, ∂U ∂t + ∇· [(U + P) v] = Γ − Λ, (6) plus the induction equation for a weakly-ionized multifluid plasma, ∂B ∂t = ∇× (v×B) − ∇× [ηo∇×B] − ∇× h ηh (∇×B) ×Bˆ i − ∇× [ηa (∇×B)⊥ ] (7) (Wardle 2007; Balbus 2011), where ρ is mass density, v is fluid velocity, P is thermal pressure, and U is the density of kinetic plus internal energy. The quantity Γ − Λ is the net rate of – 5 – heating (Γ) minus cooling (Λ) per unit volume due to the absorption and emission of photons, x-ray absorption, cosmic-ray ionization, etc. Strictly speaking, v is the velocity of the neutral particles which form the bulk of the weakly-ionized plasma. In Equation (7), Bˆ is the unit vector parallel to B and (∇×B)⊥ is the component of ∇×B perpendicular to B. The plasma is described by four transport coefficients: the shear viscosity, α, and the diffusivities ηo, ηh, and ηa associated with Ohmic dissipation, the Hall effect, and ambipolar diffusion, respectively. The solution must satisfy boundary conditions at infinity and also at the plasma-body interface. We work in a frame whose origin lies somewhere inside the body and moves with it. Then the boundary conditions at infinity are lim |x|−→∞        ρ(x, t) P(x, t) U(x, t) v(x, t) B(x, t)        =        ρ0 P0 U0 v0 B0        , (8) where the subscript “0” denotes values in the undisturbed plasma. The free-stream velocity v0 and ambient magnetic field B0 are parameters. The ambient density, pressure, and energy density must be obtained from a model of the protoplanetary disk. At the plasma-body interface the velocity satisfies the no-slip condition, v = 0 [plasma-body interface], (9) and the magnetic induction satisfies nˆ× (Bp/µp − Bb/µb) = 4π c K [plasma-body interface] (10) and nˆ·(Bp − Bb) = 0 [plasma-body interface] (11) (e.g., Jackson 1975), where µ is the magnetic permeability and K is the surface current. The unit normal, nˆ, points away from the body and quantities with subscripts b (“body”) and p (“plasma”) are evaluated just in- and outside of the body, respectively, in Equations (10)– (11). Dimensional analysis of Equations (4)–(7) yields the fundamental length- and time scales for the shear flow. The thickness of the shear layer is predicted to be of order Lsf, where Lsf ≡ (ηα) 1/2 B0 ∼ 10  η 1015 cm2 s −1 1/2  α 10−5 Poise 1/2  B0 0.1 G−1 km, (12) – 6 – and2 η ≡ ηo + ηa + ηh. (13) It is important to note that the scaling for η in Equation 12 is appropriate for a dustfree plasma. If dust is present in the disk, the diffusivities and thus Lsf will change (See Sections 3 and 5.1). There are three characteristic time scales: the time scale for viscous forces to accelerate the plasma, τacc ≡ ρLsf 2 α ∼ 106  nH 1013 cm−3   η 1015 cm2 s −1   B0 0.1 G−2 s, (14) the advection time scale, τadv ≡ Lsf v ∼ 10  η 1015 cm2 s −1 1/2  α 10−5 Poise 1/2  B0 0.1 G−1  v0 km s −1 s, (15) and the magnetic diffusion time τdif ≡ Lsf 2 η ∼ 10−3  α 10−5 Poise   B0 0.1 G−2 s. (16) The scaling in Equations (12), (15), and (16) is appropriate if α is the molecular viscosity of the gas (e.g., Schaefer 2010); however the effective viscosity could be much larger if the flow is turbulent. Although the transition from laminar to turbulent flow is not well understood for MHD flow over bodies, this transition has been widely studied for MHD channel flow (Hartmann 1937; Thess et al. 2007 and references therein) and is generally determined by the ratio F ≡ Re Ha (17) of the Reynold number, Re ≡ ρvL α , (18) to the Hartmann number, Ha ≡ BLr σo α , (19) where v is the velocity at the center of the channel, L is the channel width, and σo is the Ohmic conductivity of the plasma (See Equation [81]). For MHD channel flow with a 2 It’s important to note that in some cases the Hall diffusivity ηh can be negative (Wardle 2007), in which case the definition given in Equation (13) would not be useful. However we adopt this definition because in our models ηh is positive in the midplane of the disk for all radii of interest. – 7 – uniform magnetic field perpendicular to the walls, the flow is found to remain laminar if the parameter F remains less than a critical value Fc ∼ 100 determined by both experiments and numerical simulations (Moresco & Alboussi`ere 2004; Thess et al. 2007). Although we obviously are not considering MHD channel flow, it seems plausible that the idealized flows over planar body surfaces (See Section 4) considered in this paper may behave similarly. For realistic bodies, Fc will probably depend greatly on the geometry of the flow; however here we will assume that Fc is equal to the above value. If α is the molecular viscosity of the plasma, we find that Re = ρv0Lsf α ∼ 105  nH 1013 cm−3   v0 km s−1   η 1015 cm2 s −1 1/2  B0 0.1 G−1  α 10−5 Poise −1/2 , (20) Ha = B0Lsfr σp α ∼ 1010  η 1015 cm2 s −1 1/2  σo 106 s −1 1/2 , (21) and hence that F = Re Ha ∼ 10−5  nH 1013 cm−3   v0 km s−1   σo 106 s −1 −1/2  B0 0.1 G−1  α 10−5 Poise −1/2 , (22) which shows that F ≪ Fc, suggesting that the flows we study are indeed laminar. If these flows are in fact turbulent, a larger effective viscosity could be implied. A probable upper limit on α is set by assuming that the hydromagnetic turbulence which likely dominates the large-scale dynamics of a protostellar disk extends all the way down to scales of order the body size, Lbody. If this were true, α would be comparable to αt = κtρcsH (23) (Shakura & Sunyaev 1973), where κt is a dimensionless number, cs is the sound speed, H = cs/Ω is the local disk scale height, and Ω is the local angular velocity. Putting in numbers appropriate for the asteroid belt gives αt ∼ 104  κt 0.01   nH 1013 cm−3   T 100 K1/2  H 1012 cm Poise, (24) so that αt ∼ 109α and Lsf ≫ Lbody for bodies of asteroidal size. However even in this highly speculative scenario the time scales in Equations (14)–(16) would all be . 100 yr., and therefore steady flow can still be assumed. The electric field inside the body is coupled to the electric field in the shear layer via boundary conditions at the body-plasma interface (see Equations [31]–[32]); thus it is necessary to calculate the electric field Ep in the plasma. Equation (7), Faraday’s Law, – 8 – Helmholtz’ Theorem3 , and macroscopic charge neutrality together imply that in a steady flow Ep = − v c ×B + ηo c ∇×B + ηh c h (∇×B) ×Bˆ i + ηa c (∇×B)⊥ . (25) Equation (25) says that in general the electric field at each point in the shear flow is the sum of four contributions, Ep = Em + Eo + Eh + Ea, (26) of which the first is the motional E-field. We refer to the other contributions as the Ohm field, Eo ≡ ηo c ∇×B, (27) the Hall field, Eh ≡ ηh c h (∇×B) ×Bˆ i , (28) and the ambipolar field, Ea ≡ ηa c (∇×B)⊥ , (29) because they are the electric fields required to maintain the magnetic field gradients set up by Ohmic dissipation, etc. The boundary condition on the electric field at infinity is lim |x|−→∞ Ep = E0. (30) The boundary conditions at the body-plasma interface are nˆ× (Ep − Eb) = 0, (31) and nˆ·(Dp − Db) = 4πΣ, (32) where Σ is the surface charge density, D = ǫE is the dielectric displacement and we assume, for the present, that the dielectric constant, ǫ, is a scalar. 2.2. The Body Interior To calculate the electric and magnetic fields inside the body one must solve Maxwell’s equations there. The light-crossing time for an asteroid of size Lbody, τℓc = 0.3  Lbody 100 km ms, (33) 3A vector field is uniquely determined if its divergence and curl are both known. – 9 – is utterly negligible compared to the flow time scales (Section 2.1), so the fields inside the body are well approximated by Faraday’s Law, ∇×Eb = − 1 c ∂Bb ∂t , (34) and Amp`ere’s Law with zero displacement current, ∇×Bb = 4πσbµb c Eb + 1 µb ∇µb×Bb. (35) We have assumed that the macroscopic charge density vanishes inside the body and that Ohm’s Law holds with a scalar conductivity σb. However Equations (34)–(35) make no assumptions about the position dependence of the material properties µb, σb, and dielectric constant ǫb. Equations (34)–(35) describe the dissipation of electromagnetic field energy by Ohmic heating of the body material. This can be seen by momentarily neglecting the position dependence of µb. Then Equations (34) and (35) reduce to diffusion equations for the electric field, ∂Eb ∂t = c 2 4πσbµb ∇ 2Eb, (36) and magnetic field, ∂Bb ∂t = c 2 4πσbµb ∇2Bb. (37) In an isolated body, Ohmic dissipation would eventually reduce the electromagnetic field energy in the body to zero. However the electric and magnetic fields in the shear layer and surrounding plasma constitute an energy reservoir which tends to replenish losses inside the body, via the diffusion of E and B from the plasma into the body. The time to reach a steady state in which diffusive gains balance Ohmic losses is just the diffusion time, τEB = 4πσbµbL 2 body c 2 . (38) Estimates of τEB depend on the very uncertain electrical conductivities of primitive solar system materials. For example, estimates of σb for plausible asteroidal constituents span ∼ 16 orders of magnitude (Ip & Herbert 1983), and the situation is further complicated if the body contains a spatially connected metallic phase (H. Watson, private communication). Direct measurements of σb have been carried out for a few chondritic meteorites (Schwerer et al. 1971; Brecher 1973). They are consistent with an Arrhenius-law temperature dependence, σb = σ0 exp (−Ea/kBT), (39) – 10 – with σ0 ∼ 1012 s −1 and Ea ≈ 0.1 eV. In the absence of more extensive measurements we will adopt these values, giving τEB ∼ 10–103 µb  Lbody 100 km2 s (40) over the temperature range 100–200 K. Given the shortness of the times in Equation (40), it seems reasonable to neglect the time dependence of Eb and Bb, uncertainties in σb notwithstanding. Then Faraday’s Law reduces to ∇×Eb = 0, (41) i.e., finding Eb inside the body reduces to a boundary-value problem in electrostatics. Once Eb is known, Amp`ere’s Law reduces to a set of coupled ordinary differential equations for Bb. In Sections 4.1–4.2 we show how this all works out for two extremely simple models. 2.3. Self Heating Dissipative processes inside the shear flow transform ordered kinetic energy into heat. Since the transport coefficients depend on temperature, it is important to know whether the heating is large enough to raise the temperature above ambient. The friction associated with viscosity heats the plasma at a rate Γvisc = τ : ∇v, (42) where τ is the viscous stress tensor and ∇v is the velocity gradient tensor. In Cartesian coordinates τij = α  ∂vi ∂xj + ∂vj ∂xi  i, j = x, y, z (43) and Γvisc = X i,j τij ∂vi ∂xj . (44) Taking v ∼ v0 and ∂vi/∂xj ∼ v0/Lsf gives the order-of-magnitude estimate Γvisc ∼ 10−7  v0 km s−1 2  B0 0.1 G2  η 1015 cm2 s −1 −1 erg cm−3 s −1 . (45) Notice that the RHS of Equation (45) is independent of the viscosity. The shear flow is also heated by the friction associated with streaming motions between charged and neutral particles; however electron-neutral scattering can generally be neglected – 11 – (Balbus 2011). Let ni be the number density of ions with mass mi and charge Zie and nn be the number density of neutral particles with mass mn. If E or B is nonzero, the ions will be accelerated by the Lorentz force. On a very short time scale they reach a terminal velocity with respect to the neutrals, such that the Lorentz force is balanced by the drag force associated with elastic ion-neutral scattering. The terminal drift velocity is vi − v = c " βi Bp·E ′ p
Bp B3 p + β 2 i 1 + β 2 i E ′ p×Bp B2 p + βi 1 + β 2 i Bp× E ′ p×Bp
B3 p # (46) (Wardle 1998), where vi is the ion velocity and E ′ p = Ep + v c ×Bp (47) is the electric field in the frame of the neutral particles (=the bulk plasma). The drift velocity depends on the dimensionless ion Hall parameter, βi ≡ Ωi τin, (48) where Ωi = ZieB mic (49) is the ion cyclotron frequency, τin ≡ 1 + mi/mn nnhσvi in (50) is the time scale for the ions to reach their terminal drift speed, and hσvi in is the momentum transfer rate coefficient for elastic ion-neutral scattering. Ion-neutral streaming heats the neutral gas at a rate Γin ≈ ni nn hσvi in µin |vi − v| 2 (51) (Draine et al. 1983; Chernoff 1987), where µin is the ion-neutral reduced mass. 3. Transport Coefficients In Section 4 we present analytic solutions which describe shear flow past an infinite slab. These solutions are cast in dimensionless units which depend only on certain combinations of the transport coefficients. However to express the solutions in terms of ordinary units we require numerical values of α, ηo, ηh, and ηa. We estimate these as follows. – 12 – 3.1. Protoplanetary Disk Model We use the density and temperature predicted by the minimum mass solar nebula (MMSN) model (Weidenschilling 1977b; Hayashi 1981) as described by Sano et al. (2000). The gas temperature is T(R) = 280  R AU−1/2 K, (52) where R is distance from the central star. The mass density is ρg(R, Z) = ρg(R, 0) exp " −  Z H 2 # , (53) where Z is vertical distance from the midplane. The midplane density of the MMSN is ρg(R, 0) = 1.4 × 10−9  R AU−11/4  M∗ M⊙ 1/2  µg 2.34 1/2 g cm−3 , (54) where M∗ is the mass of the central star, µgmH is the mean mass per neutral particle, and mH is the hydrogen mass. We set M∗ = 1 M⊙ and µg = 2.33, the latter corresponding to a gas composed of molecular hydrogen and helium with number densities 0.5nH and 0.1nH, respectively, where nH is the number density of hydrogen nuclei. Then the gas number density at the midplane is ng(R, 0) = 0.6 nH(R, 0) = 3.5 × 1014  R AU−11/4 cm−3 . (55) The midplane temperature and number density are plotted versus R in Figure 1. We also require the ambient magnetic field, B0. It is not yet possible to directly measure the magnetic fields in protoplanetary disks; however several lines of evidence suggest that B0 ∼ 0.1–1 G, including simulations of star formation (Desch & Mouschovias 2001) and constraints implied by protostellar mass accretions rates (Wardle 2007; Bai & Goodman 2009). We will adopt the value B0 = 0.3 G for our calculations discussed in Section 5 which is approximately the geometric mean of the extremes. 3.2. Ionization Equilibrium The magnetic diffusivities depend on the number densities of ions, electrons, and charged dust grains (Section 3.3). To estimate these we use the semianalytical model of ionization – 13 – equilibrium developed by Okuzumi (2009). Okuzumi’s model calculates the number densities of electrons, dust grains in different charge states, plus the total ion number density, ni ≡ X k n (k) i , (56) where the sum is over different ionic species k = Mg+ , HCO+ , H + 3 , . . .
. Okuzumi’s results are in excellent agreement with full numerical calculations in which different ionic species are treated explicitly.4 The inputs to Okuzumi’s model are listed in Table 1, where βr ≡ 1 ni X k n (k) i β (k) r (57) is the average recombination rate coefficient of the ions, ui ≡ 1 ni X k n (k) i u (k) i (58) is the average of the mean ion thermal speed, and ζ ≡ 1 ng X j n (j) g ζ (j) , (59) where ζng is the number of ionizations per unit volume per unit time and the sum is over all gas-phase species j. More detailed treatments of ionization equilibrium in disks (e.g., Sano et al. 2000) show that the abundance of atomic ions exceeds the abundance of molecular ions by 2–3 orders of magnitude. We will take Mg+ to be a proxy for the ions and set ui =  8kBT πµimH 1/2 , (60) where µi = 24, and βr = 2.8 × 10−12 (T/300 K)−0.7 , (61) which is the rate coefficient for radiative recombination of Mg+ (McElroy et al. 2013). The ionization rate, ζ, is crucial but uncertain. In a protoplanetary disk it has contributions from cosmic rays, x rays, and live radionuclides so that ζ = ζcr + ζxr + ζra. (62) 4See Okuzumi (2009) Figure 3 but note that the figure has an error: the code which plotted Figure 3b multiplied every quantity by a factor of ≈ 0.36 (S. Okuzumi, private communication). – 14 – The cosmic ray ionization rate is ζcr = ζcr,0 Fcr(R, Z), (63) where ζcr,0 is the rate far from the midplane. The factor Fcr(R, Z) ≡ 1 2  exp −Σ + g (R, Z)/Σcr h 1 + Σ + g /Σcr
3/4 i−4/3 + exp −Σ − g (R, Z)/Σcr h 1 + Σ − g /Σcr
3/4 i−4/3  (64) describes attenuation by the gaseous disk (Okuzumi 2009), with Σ + g (R, Z) ≡ Z +∞ Z ρg (R, Z′ ) dZ′ , (65) Σ − g (R, Z) ≡ Z Z −∞ ρg (R, Z′ ) dZ′ , (66) and Σcr = 96 g cm−2 . Here we are interested only in points with Z = 0 so that Σ + g (R, Z) = Σ− g (R, Z) = 1 2 Σg(R), (67) where Σg(R) = 1700  R AU−3/2 g cm−2 (68) is the total surface density. Since the cosmic ray ionization rates of H2 and He are the same to within about 20%, we will take ζcr,0 to be the rate for molecular hydrogen. Values of the latter inferred from observations of molecular clouds span about two orders of magnitude, from 3 × 10−18 s −1 to 4 × 10−16 s −1 (Padovani et al. 2009 and references therein). We will adopt the value suggested by Dalgarno (2006), ζcr,0 = 5 × 10−17 s −1 , (69) which is approximately the geometric mean of the extremes. Our models also include ionization by stellar x rays and radionuclides. For the former we adopt the fit of Turner & Sano (2008) to calculations of x-ray ionization by Igea & Glassgold (1999): ζxr(R, Z) = ζxr,0 R AU
−2  Lxr 2×1030 erg s−1  ×  exp −Σ + g (R, Z) /Σxr + exp −Σ − g (R, Z) /Σxr , (70) – 15 – where Σxr = 8.0 g cm−2 , ζxr,0 = 2.6 × 10−15 s −1 , and the stellar x-ray luminosity, Lxr, is a parameter (Table 1). The rate of ionization by radioactivity depends strongly on the presence or absence of short-lived radionuclides (SLRs) like 26Al. It increases from ζra = 1.4×10−22 s −1 if SLRs are absent to ζra = 7.6 × 10−19 s −1 if SLRs are present with the abundances inferred for the early solar nebula (Umebayashi & Nakano 2009). The fact that chondrules are a few Myr younger than the first solids to condense from the solar nebula (the CAIs) suggests that the solar nebula was a few Myr old when the asteroids formed, at least for the parent bodies of chondritic meteorites. Since the half life of 26Al is 0.7 Myr, we will neglect ionization by SLRs. The midplane ionization rates are plotted vs. R in Figure 2. The ionization balance and magnetic diffusivities are sensitive to the presence of small (with radii a less than a few µm) dust grains. Dust tends to sweep up electrons, the only particles that are well coupled to B at the densities of interest here (Section 3.3). This reduces the coupling between gas and magnetic field with a consequent increase in the diffusivities. For example Wardle (2007) showed that a standard interstellar population of grains with a = 0.1 µm would increase the diffusivities enormously, by ∼ 8 orders of magnitude in the midplane. However grain growth can reduce the total dust surface area5 to a point where the effect of grains on the diffusivities becomes negligible; Wardle (2007) showed that (for nonfractal grains) this occurs when a exceeds a few µm. In this paper we are concerned with disk ages of a few Myr, which is longer than the time scale for grain growth (Dullemond & Dominik 2005). The expectation that grain growth is significant at t ∼ Myr is borne out by millimeter observations of disks, which suggest that most of the dust is locked up in millimeter and larger-size particles at about t = 2 Myr (Williams 2012 and references therein). However Spitzer observations of the infrared features produced by amorphous (at 10 µm) and crystalline (23 µm) silicates imply the presence of µm-sized grains in numerous disks (e.g., Furlan et al. 2011). Because these disks are optically thick the observed small grains lie in the disk atmospheres, i.e., not in the midplane. Nevertheless, it is important to know whether enough µm-sized dust might reside in the midplane to affect the diffusivities. To estimate the density of small dust at the midplane we exploited the calculations of Birnstiel et al. (2011), who described the equilibrium size distribution established by coagulation and fragmentation. Birnstiel et al. modeled a grain of radius a as a constantporosity sphere of mass m(a) = 4 3 π ρs a 3 , (71) 5 If the grains grow as spheres with fractal dimension Df = 3. However the initial stages of grain growth generally produce fractal structures with Df < 3 (Blum & Wurm 2008 and references therein). To be consistent with the grain model of Birnstiel et al. (2011) considered below [see Equation (71)] we will set Df = 3. – 16 – where ρs is the bulk density of the porous solid. At any given time the mass distribution extends from the monomer mass, m0, to some upper limit, mc, determined by the growth process. For a range of models with different assumptions about the underlying physics they found that the mass distribution can be approximated by a power law6 , n(m) = A m−γ m0 < m < mc, (72) where n(m) dm is the number density of grains with masses in (m, m + dm) and γ lies in the range from about 0.5 to 2. We set m0 = m(a0) and mc = m(ac), where a0 = 0.1 µm and ac = 1 mm, corresponding to the coagulation of interstellar dust to sizes consistent with the millimeter observations of disks. The total abundance of dust with all sizes is characterized by χd, where χdρg is the total dust mass density and χd ∼ 0.01. As a measure of how much small dust is present we define χsd(a)ρg to be the mass density of dust with radii less than a. Then it is easy to show that Equation (72) implies χsd(a) = " (a/a0) 3(2−γ) − 1 (ac/a0) 3(2−γ) − 1 # χd if γ 6= 2. (73) Figure 3 plots χsd for a few plausible values of γ. It is apparent that the theory allows a very broad range in the mass fraction of micron-size and smaller grains, from ∼ 10−14 to 10−2 . Figure 4 plots the abundances of ions, electrons, χd+, and χd− relative to nH calculated using Okuzumi’s model as a function of R in the disk, assuming χsd = 10−4 . Here χd+ and χd− describe the total amount of positive and negative charge contained by dust grains, and are defined as χd+ = X Z>0 nd(Z) nH |Z| (74) and χd− = X Z<0 nd="" z="" nh="" 75="" where="" is="" the="" charge="" of="" grains="" in="" units="" elementary="" we="" assume="" that="" all="" dust="" are="" single-sized="" with="" radius="" a="" and="" consider="" two="" cases:="" m="" when="" find="" positively="" negative="" charged="" dominant="" carriers="" plasma="" for="" r="" 2="" au="" 4="" negatively="" 6strictly="" speaking="" distribution="" given="" equation="" 72="" does="" not="" account="" effect="" settling="" disk="" however="" numerical="" simulations="" performed="" by="" birnstiel="" et="" al="" 2011="" including="" effects="" turbulent="" mixing="" vertically="" integrated="" can="" be="" approximated="" power="" laws="" depending="" on="" whether="" or="" greater="" than="" certain="" size="" see="" their="" section="" 3="" figure="" 17="" still="" dominate="" over="" electrons="" while="" ions="" become="" carrier="" positive="" species="" if="" grain="" increased="" to="" observe="" abundances="" exceed="" those="" 1="" 5="" more="" abundant="" continue="" particular="" interest="" flows="" considered="" this="" paper="" electron="" abundance="" xe="" ne="" 76="" since="" only="" particles="" strongly="" tied="" magnetic="" field="" midplane="" shown="" figures="" 6="" calculations="" were="" carried="" out="" values="" sd="10−4" toward="" high="" end="" plausible="" range="" assumed="" single-size="" dust-free="" also="" comparison="" confirm="" very="" sensitive="" smallest="" latter="" smaller="" few="" wardle="" 2007="" example="" reduced="" about="" three="" orders="" magnitude="" at="" once="" small="" falls="" unless="" 0="" will="" calculate="" diffusivities="" it="" should="" kept="" mind="" power-law="" steeper="" our="" underestimate="" viscosity="" shear="" depends="" temperature="" adopt="" calculated="" schaefer="" 2010="" which="" well="" 26="" t="" k="" poise="" 77="" both="" ortho-="" para-h2="" 300="" obtained="" from="" particle="" as="" follows="" :="" o="ec" 78="" h="" 79="" 18="" p="" 80="" b="" x="" j="" nj="" zj="" 81="" ohmic="" conductivity="" 82="" hall="" 83="" pedersen="" q="" sums="" mass="" mj="" number="" density="" zje="" include="" various="" states="" within="" four="" standard="" deviations="" mean="" okuzumi="" 2009="" addition="" each="" depend="" its="" dimensionless="" parameter="" turn="" drag="" time="" jn="" defined="" force="" produced="" friction="" neutral="" gas:="" fjn="" vj="" vn="" 84="" ion="" 50="" vi="" 10="" 9="" cm3="" s="" use="" en="σg" ng="" 85="" 128kbt="" me="" 86="" g="" cm2="" draine="" 1983="" expression="" salpeter="" 1979="" limit="" subsonic="" gas-grain="" drift="" implies="" dn="√" vth="" 87="" 2kbt="" gmh="" 88="" 19="" coupling="" between="" characterized="" then="" coupled="" lines="" moves="" them="" frozen="" into="" recover="" ideal="" mhd="" parameters="" plotted="" vs="" 7="" 8="" different="" model="" adopted="" here="" predicts="" e="" i="" region="" corresponding="" asteroid="" belt="" regime="" move="" but="" do="" notice="" extremely="" so="" terms="" conductivities="" negligible="" hold="" present="" free="" fliers="" scenario="" nevertheless="" affects="" profoundly="" reducing="" gas-phase="" no="" another="" ionization="" illustrate="" dependence="" specific="" now="" compare="" diffusivity="" described="" recent="" work="" dzyurkevich="" 2013="" although="" used="" makes="" similar="" assumptions="" physical="" conditions="" they="" make="" significantly="" regarding="" rates="" properties="" most="" important="" difference="" models="" fractal="" dimension="" df="" set="" equal="" modeled="" fluffy="" aggregates="" rate="" an="" order="" introduce="" concept="" metal="" line="" beyond="" freeze="" molecular="" ionic="" quantitative="" function="" these="" listed="" table="" note="" recreate="" results="" required="" total="" magnesium="" relative="" hydrogen="" nuclei="" 100="" times="" larger="" value="" quoted="" appendix="" 11="" 12="" show="" contained="" left="" panel="" 20="" positively-="" negatively-charged="" completeness="" plot="" mg="" hco="" right="" demonstrates="" comparing="" qualitatively="" fractional="" predicted="" exceeds="" expected="" leads="" large="" largest="" general="" main="" just="" increase="" thickness="" layer="" steady="" flow="" past="" infinite="" slab="" complexity="" governing="" equations="" seems="" inevitable="" around="" bodies="" realistic="" shapes="" have="" numerically="" highly="" idealized="" problem="" 13="" studied="" analytically="" geometry="" unrealistic="" solutions="" how="" electric="" inside="" body="" couples="" self="" heating="" insignificant="" point="" equals="" undisturbed="" energy="" omitted="" validity="" assumption="" checked="checked" conserved="" identically="" necessary="" solve="" momentum="" induction="" v="[vx(z)," bp="[0," ep="" plus="" eb="" amp="" ere="" law="" bb="" parallel="" first="" simplest="" possible="" case="" ambient="" surfaces:="" b0="" y="" 89="" 90="" 21="" pressure="" constant="" under="" satisfied="" reduces="" single="" ode="" vx="" d="" dz2="0," 91="" boundary="" limz="" 92="" w="" 93="" symmetry="" requires="" even="" solution="" explicitly=""> 0. The electric field inside the body can be calculated exactly. Outside the body the magnetic field is simply advected by the flow. Since there is no compression or bending of the field lines, the magnetic field is uniform with Bp = B0 everywhere. It follows that the Ohmic, ambipolar, and Hall electric fields all vanish in the plasma and hence that Ep equals the motional field. Now the no-slip boundary condition implies that Ep = 0 just outside each body surface. We find the electric field just inside each surface by applying the boundary conditions across the plasma sheath, whose thickness is neglected because it is of order meters and thus much less than the dimensions of the shear layer. Using Equations (31) and (32) with Ep = 0, we find Ebs± = ∓ 4πΣ ǫb ˆz (94) where the upper and lower signs correspond respectively to the upper and lower slab surfaces. The two electric fields in Equation (94) are just the electrostatic fields associated with electric charge in the plasma sheaths at the upper and lower surfaces; their existence and potential consequences for heating are not the subject of this investigation. In any case, the sum of the “sheath fields” vanishes for the geometry considered here: In order to satisfy ∇×Eb = 0 and ∇·Eb = 0 everywhere inside the body each sheath field must be uniform (for an infinite slab), so that Equation (94) is actually valid throughout the body interior. But the sheath fields sum to zero so the total electric field inside the body is Eb = 0. We conclude that it is possible for induction heating to vanish when one accounts for the interaction of a body with the flow around it. Of course the fact that Eb = 0 exactly is a consequence of the extreme symmetry of an infinite slab. For realistic bodies of finite extent the electric field will not be zero, but instead attain some value dependent on the nonzero magnetic field gradients associated with departures from planar symmetry. The electric field in the plasma must vary continuously with z, from E0 at infinity to zero at the slab surfaces. One could calculate Ep(z) in principle from the velocity solution – 22 – by evaluating Ep(z) = − v(z) c ×B0. (95) However it is well known that the solution for v does not exist for an infinite plane: the solution to Equation (91) cannot satisfy boundary conditions (92) and (93) simultaneously because the shear layer is infinitely thick. Since the magnetic forces vanish for this geometry, a rough picture of the variation can be obtained by considering viscous flow past a semiinfinite flat plate, for which the shear layer has finite thickness δ. Then the boundary condition at infinity can be replaced by vx(W + δ) = v0. (96) From scaling arguments, Blasius (1908) inferred that δ = δB(x) ≈  αx ρv0 1/2 (97) for flow past a semi-infinite plate, where x is the distance downstream from the leading edge of the plate. The thickness of the shear layer varies with x but the rate of change, dδB dx = 1 2  α ρv0x 1/2 , (98) goes to zero as x−→∞. Thus, δ becomes approximately independent of x for large x and the flow approximates flow past an infinite slab. We will interpret δ in Equation (96) to be δB(x0), where x0 is a large but unspecified value of x for which the approximation holds to some desired precision. Then the solution to Equation (91) subject to the approximate boundary conditions (93) and (96) is vx(z) ≈  z − W δ − W  v0 (99) and the electric field in the plasma is Ep(z) ≈ −  z − W δ − W  E0zˆ. (100) That the electric field in Equation (100) has nonzero divergence shows that the approximation for δ is unphysical. However the nonzero macroscopic charge density implied by Gauss’s Law is e(ne − ni) = v0B0 4πc(δ − W) , (101) which goes to zero as δ → ∞. – 23 – 4.2. Perpendicular Field Geometry Now consider the case where the ambient magnetic field is perpendicular to the slab faces, with B0 = B0 zˆ. Because the slab is infinite all variables depend only on z and the electric field can be calculated trivially. Far from the body, where magnetic field gradients vanish, the electric field in the body frame is just the motional field with components7 E0 = [0, v0B0/c, 0] . (102) For steady flow the electric field at an arbitrary point in the plasma8 must also satisfy ∇×Ep = 0 which, together with the fact that Ep only depends on z, implies that Ep = E0 everywhere in the plasma. Now the electric field inside the body is also uniform because ∇·Eb = 0 and ∇×Eb = 0. Noting that Ep is tangent to the slab surfaces, and that the tangential component of E does not change across the surfaces, we see that Eb = E0 everywhere inside the body, and hence that E = E0 everywhere. This simple example shows that it is possible to have electric fields ∼E0 inside the body. In the following subsections we show how this is physically possible by calculating the motional, Ohm, Hall, and ambipolar contributions to Ep. 4.2.1. The Body Interior For infinite planar geometry the condition ∇·B = 0 implies that Bz is constant, so Bz = B0 − ∞ < z < +∞. (103) We assume a solution of the form Bb = [B1x(z), B1y(z), B0], where the notation indicates that B1x and B1y are regarded as small perturbations: B1x, B1y ≪ B0. This is necessary because the x- and y-components of the magnetic field are small outside of the slab (Section 4.2.2) and it would not be possible to satisfy the boundary conditions on B at the slab surfaces, in general, if the components of the magnetic field were small outside and large inside. Because Eb is known, B1x and B1y follow from Amp`ere’s Law. Writing out the x and y components of Equation (35) gives dB1x dz = 4πµbσbEy c (104) 7Looking at the velocity and magnetic field solutions in Equations (130)–(133), we note that the motional electric field in the plasma also contains an unphysical non-zero z-component (Epz ) in this case. However Epz is found to be of order ∼ v0B1x,y c , which is ≪ E0 and is therefore neglected. 8To avoid numerous disclaimers we imply henceforth that “in the plasma” means points in the plasma outside of plasma sheaths. – 24 – and dB1y dz = − 4πµbσbEx c , (105) where it has been assumed for convenience that µb is constant and σb is the electrical conductivity inside the body. It is apparent from the symmetry of the problem that B1x(z) and B1y(z) must either be odd or even functions. However even functions are excluded: B1 is odd outside the body (Section 4.2.2) and it would not be possible to satisfy the boundary conditions at the body surfaces in general if B1 were even inside and odd outside. The boundary conditions on Equations (104)–(105) are therefore B1x(0) = 0 (106) and B1y(0) = 0. (107) The solution of Equations (104)–(107) is B1x(z) = 4πµbE0 c Z z 0 σb (z ′ ) dz′ (108) and B1y(z) = 0, (109) where we have exploited the fact that E = E0 everywhere. In practice Equation (108) sets an upper limit on the conductivity, above which our treatment of B1x as a perturbation would break down. The conductivity appears inside the integral in Equation (108) to allow for the sensitive dependence of σb on the temperature (Equation [39]), which varies with z if the body is heated significantly. 4.2.2. The Shear Flow Now consider the velocity and magnetic field in the plasma. We seek solutions to Equations (5) and (7) of the form v = [vx(z), vy(z), 0] , (110) Bp = B0 + [B1x(z), B1y(z), 0] , (111) and P = P(z), where we have anticipated that the magnetic field is deflected only slightly and regard B1x and B1y as infinitesimal perturbations. Substituting the assumed solutions into Equations (5) and (7), setting the time derivatives to zero, and retaining terms up to – 25 – first order in B1x and B1y yields six differential equations for the components of v and B. However the z component of the linearized induction equation is satisfied identically. The z component of the linearized momentum equation reduces to an equation for the pressure gradient required to maintain hydrostatic equilibrium in the z direction; it decouples from the other equations when the temperature dependence of the transport coefficients is neglected. This leaves four ODEs for vx, vy, B1x, and B1y: (ηo + ηa) d 2B1x dz2 + ηh d 2B1y dz2 + B0 dvx dz = 0, (112) − ηh d 2B1x dz2 + (ηo + ηa) d 2B1y dz2 + B0 dvy dz = 0, (113) B0 4π dB1x dz + α d 2 vx dz2 = 0, (114) and B0 4π dB1y dz + α d 2 vy dz2 = 0. (115) It follows from Equations (112)–(115) that vx and vy are even functions of z and that B1x and B1y are odd. Thus it is only necessary to calculate the flow explicitly for z > 0. Equations (112)–(115) contain the second derivatives of vx, vy, B1x, and B1y so that eight boundary conditions are required. At the body surface the velocity satisfies the no-slip condition, vx(W) = 0, (116) vy(W) = 0, (117) and the tangential components of the magnetic field satisfy Equation (10): B1x (W−) µb = B1x (W+) µp , (118) and B1y (W−) µb = B1y (W+) µp . (119) Far from the body the velocity approaches the free-stream velocity, lim z→+∞ vx = v0, (120) lim z→+∞ vy = 0, (121) and the magnetic field perturbations approach constant (infinitesimal) values: lim z→+∞ B1x = const., (122) – 26 – and lim z→+∞ B1y = const. (123) For a body of finite size B1x and B1y would vanish at infinity. However our infinite planar slab acts as an infinite current sheet; thus it produces a magnetic field whose magnitude is constant everywhere outside the body. We introduce the dimensionless position ¯z ≡ z/Lsf, velocity, u = v/v0, and magnetic induction, b ≡ B1/B∗, where B∗ ≡ v0 (α/η) 1/2 (124) is the natural unit of magnetic induction. In dimensionless units Equations (112)–(115) become (1 − λh) b ′′ x + λhb ′′ y + u ′ x = 0, (125) − λhb ′′ x + (1 − λh)b ′′ y + u ′ y = 0, (126) b ′ x + 4πu′′ x = 0, (127) and b ′ y + 4πu′′ y = 0, (128) where primes indicate derivatives with respect to ¯z. Apart from the boundary conditions the exterior solution depends on just one dimensionless parameter, λh ≡ ηh ηh + ηo + ηa , (129) which varies from zero (no Hall effect) to unity (Hall effect only, no Ohmic dissipation or ambipolar diffusion). Figures 14–15 show the value of λh as a function of R for two cases: the dust-free case, and the case where χsd = 10−4 and the dust is assumed to be single-sized with a = 1 µm. We point out that in both cases the shear flows are “Hall dominated” in the sense that λh > 1 2 for most radii of interest.9 Equations (125)–(128) are solved in Appendix A, where we show that vx = h 1 − cos kIξ e−kRξ i v0, (130) vy = sin kIξ e−kRξ v0, (131) 9Although the Hall effect dominates the dynamics of the flows studied here, this may not be the case for other flows, particularly in the context of the magnetorotational instability (MRI; Balbus & Hawley 1991). For example, 3D numerical MHD simulations done by Sano & Stone (2002) show that the Hall effect may be neglected when determining the condition for MRI turbulence, even though the value of the Hall diffusivity may exceed that of the Ohmic and ambipolar diffusivities. – 27 – B1x = B1x W+
+ r 4π D h cos θ/2 − cos (kIξ − θ/2) e −kRξ i B∗, (132) and B1y = r 4π D h sin θ/2 + sin (kIξ − θ/2) e −kRξ i B∗, (133) where ξ ≡ z − W Lsf (134) is dimensionless height above the body surface and B1x W+
= 4πµpE0 c Z W 0 σb (z ′ ) dz′ . (135) The solution depends on λh via the dimensionless quantitities kR (λh) ≡ cos θ/2 √ 4πD (136) and kI (λh) ≡ sin θ/2 √ 4πD (137) (Figure 16), where θ ≡ arctan  λh 1 − λh  (138) and D ≡ q λh 2 + (1 − λh) 2 . (139) Notice that 0 ≤ θ ≤ π/2 and kR, kI ≥ 0. Now that the flow velocity and magnetic field solutions are known, the motional, Ohm, Hall, and ambipolar components of the electric field can be found by evaluating Equations (26)–(29). These fields will be discussed in the next section. 5. Discussion 5.1. Electric Fields in Multifluid Shear Flows As demonstrated in the previous section, the solutions for the velocity, magnetic, and electric fields depend only on the transport coefficients α, ηh, ηo, and ηa. Although the shear viscosity α does not depend on the amount of small dust in the disk, dust may profoundly affect the diffusivities (see Section 3.3, Figures 9–10). In order to illustrate the effects of – 28 – dust on the components of the electric field we will consider two cases: (1) the limit where all the dust is locked into millimeter or larger grains and their is no small dust present in the disk (χsd = 0) and (2) a case where the abundance of small grains is taken to be χsd = 10−4 and all the small dust is the same size with a radius of a = 1 µm. Figures 17–18 describe the electric field in a perpendicular shear flow when no small dust (χsd=0) is present in the disk. These electric fields were computed using our disk model (Section 3) where the diffusivities were calculated at R = 3 AU assuming B0 = 0.3 G. The plots are independent of v0 because E is plotted in units of E0. The motional, Ohm+ambipolar10 and Hall electric fields are indicated and the total field Ep is also plotted. The figures show the inevitable result that Ep = E0. What is interesting about the figure is how this result comes about. The Ohm and ambipolar fields are relatively small. The motional field dominates (Em ≈ E0) far from the body and the Hall field dominates (Eh ≈ E0) close to the body. For a body of realistic shape the four contributions to the electric field would obviously not “conspire” to sum to E0 as in Figures 17–18. However since the magnitudes of the various magnetic field gradients should not be very sensitive to body shape (for large bodies), it seems clear that real shear flows will have electric fields ∼ E0. If small dust is present in the disk, the diffusivities and thus the electric field will be affected. Figures 19–20 describe the electric field for the case analogous to Figures 17–18 except with χsd = 10−4 and a = 1 µm. Comparing Figures 19–20 with Figures 17–18, it is evident that even with small dust with abundances as high as χsd = 10−4 present in the disk, the flow remains Hall-dominated at R = 3 AU with a total electric field Ep = E0. However dust does produce a dramatic effect on the characteristic length scale of the shear layer Lsf, which grows from ∼ 10 km in the absence of small dust to ∼ 100 km when small dust is present. This effect is important to note because our approximation of the asteroid as an infinite plane is only valid when Lbody ≫ Lsf. Taking Lbody ∼ 100 km for a typical large asteroid, we find that our solutions should provide a decent approximation in a dust-free plasma, however their validity should be questioned when µm-sized dust is present in the abundances discussed above. In the extreme case that ∼ 0.1 µm dust grains are present in the midplane, the magnetic diffusivities will be even greater corresponding to an even larger characteristic length scale of shear layer and rendering our approximation of an asteroid as an infinite plane completely unrealistic. In this case numerical calculations of the flow over a finite body would be required. 10The sum of the Ohm and ambipolar E fields is plotted because, for perpendicular geometry, they are proportional to one another. – 29 – 5.2. The Importance of Self Heating In Section 4.2 we solved the equations of motion for a multifluid shear flow on the assumption that heating of the flow by viscous dissipation and ion-neutral scattering is negligible. This will be a good assumption if the heating rates are ≪ Λ, where Λ is the radiative cooling rate of the plasma. 5.2.1. Radiative Cooling In the absence of dust, radiative cooling is mainly due to rotational and vibrational transitions of H2, CO, and H2O molecules. We calculate the cooling rate due to molecular line emission by interpolating the results given by Neufeld & Kaufman (1993) at several values of R. These cooling rates depend on the velocity gradients in the plasma, which determine the probability that an emitted photon has of escaping (See Neufeld & Kaufman 1993). Since the viscous dissipation and ion-neutral scattering heating rates are greatest at the asteroid surface (See Equations [148] and [152] below), we use the velocity gradients in the shear flow dv dz = s dvx dz 2 +  dvy dz 2 (140) at the body surface (ξ = 0). The radiative cooling rates also depend strongly on the abundance of water vapor in the gas, which may freeze out. Since the location of the snow line in protoplanetary disks is still a topic of active research (e.g., Martin & Livio 2012 and references therein), we assume that the abundance of water nH2O/nH is constant for all R and calculate the cooling rate for two cases: nH2O/nH = 1.1×10−4 (Figures 21–22) corresponding to a standard abundance water in the gas phase, and nH2O/nH = 0 (Figures 23–24) assuming that all the water has frozen out. If dust is present, its thermal emission will also contribute to the cooling of the plasma. It is easy to show that the shear layer is optically thin to the dust thermal emission, and thus the dust radiative cooling rate is Λdr = 4πnda 2 hQabsiσsbT 4 d , (141) (Draine 2011) where σsb is the Stefan-Boltzmann constant, Td is the dust temperature, and hQabsi is the dust grains’ Planck-averaged emission efficiency with hQabsi ≈ 0.13  a 1 µm   T 100 K2 (silicates). (142) – 30 – The cooling rate due to dust thermal emission is plotted as a function of R in the disk in Figure 25, assuming that the dust grains are a single-size distribution with a = 1 µm and χsd = 10−4 . 5.2.2. Self Heating Rates The heating rate due to viscous dissipation can be calculated by evaluating Equation 44 using the velocity solution given in Equations 130–131. For the shear flow considered in Section 4.2, the components of the viscous stress tensor become τxy = τyx = 0, (143) τyz = τzy = α dvy dz , (144) τzx = τxz = α dvx dz , (145) and τxx = τyy = τzz = 0. (146) Substituting these stresses into Equation (44) gives Γvisc = α " dvx dz 2 +  dvy dz 2 # . (147) Using the solutions found in Section 4.2 to calculate the derivatives of vx and vy and substituting them into the equation above gives Γvisc = v 2 0B2 0 4πDη e −2kRξ . (148) In addition to the viscous dissipation heating rate, the heating rate in the shear flow due to friction between the ion and neutral fluids can also be calculated. Substituting the assumed forms of the plasma velocity, magnetic, and electric fields described in Section 4.2 into Equation (46), the ion-neutral velocity differences are found to be vix − vx = βi 1 + β 2 i (βiv0 − βivx + vy), (149) viy − vy = βi 1 + β 2 i (v0 − vx − βivy), (150) – 31 – and viz − vz ≈ 0. (151) The ion-neutral frictional heating rate is Γin ≈ ninnhσviinµinv 2 0β 2 i 1 + β 2 i e −2kRξ . (152) Upon inspection of Equations (148) and (152), it is clear that the heating rates due to viscous dissipation and ion-neutral scattering depend strongly on the free-stream velocity of the plasma v0 and are largest at the asteroid surface (ξ = 0). In order to assess the importance of self heating, we consider flows with v0 = 1 km/s and v0 = 5 km/s, which could correspond to flows driven by the relative motions between the asteroid and the gas on an eccentric orbit (Morris et al. 2012) or a passing shock wave in the disk. As done in the previous section, we present two cases: the dust-free case, and the case where χsd = 10−4 with a = 1 µm. Figure 26 describes the heating rates due to viscous dissipation and ion-neutral scattering at the body surface ξ = 0 as a function of R in a dust-free plasma. The left and right panels describe flows with v0 = 1 km/s and v0 = 5 km/s respectively, and assume B0 = 0.3 G. For both flows we note that the heating rate due to viscous dissipation dominates that of ion-neutral scattering at all distances R of interest. Comparing the heating rates (Figure 26) to the radiative cooling rates when gas-phase water is present in the abundance nH2O/nH = 1.1 × 10−4 (Figure 21), we observe that for flows with v0 = 1 km/s the cooling rate exceeds the heating rate for R . 4 AU. For 6 . R . 8 AU, the heating and cooling rates become comparable. If the free-stream flow velocity is increased to v0 = 5 km/s, the cooling rate still dominates the heating rate for R . 2 AU. Once R ≈ 3 AU, the heating and cooling rates become comparable. For R & 4 AU, we find that the heating rates start to exceed the cooling rate. Thus we conclude that self heating of the plasma is insignificant for flows with v0 = 1 km/s at R . 6 AU. However for faster flows with free-stream velocities & 5 km/s, significant self heating of the plasma to temperatures above ambient (Equation 52) may become important, particularly at R & 3 AU. If gas-phase water is absent in the disk, the heating rate is found to exceed the cooling rate (Figure 23) at all R for both v0 = 1 and 5 km/s, implying that significant self heating of the plasma would occur. Figure 27 shows the analogous heating rates in the case were single-sized small dust grains are present in the disk with an abundance χsd = 10−4 and radius a = 1 µm. Once again we note that viscous dissipation dominates the heating at all considered values of R. Comparing the heating rates (Figure 27) to the cooling rates by dust thermal emission (Figure 25) and molecular line emission (Figures 22) when a standard abundance of gasphase water is present in the disk, we find that the cooling rate exceeds the heating rate at – 32 – all R for v0 = 1 km/s. Under these conditions, dust emission dominates the total radiative cooling rate of the plasma for R . 4 AU, but becomes insignificant for R & 8 AU. For faster flows with v0 = 5 km/s, the cooling rate is again found to exceed the heating rate for R . 8 AU. In this case dust thermal emission again dominates the cooling for R . 3 AU, but becomes smaller than the rate of cooling by molecular line emission for R & 4 AU. If gas-phase water is absent, the cooling rate becomes dominated by thermal dust emission at all R of interest for both v0 = 1 and 5 km/s. Comparing the heating rates to the cooling rate due to dust emission, we find that the cooling rate exceeds the heating rate everywhere for v0 = 1 km/s. For faster flows with v0 = 5 km/s, the cooling rate still exceeds the heating rate for R . 4 AU. However for R & 6 AU, the heating rate is found to exceed the cooling rate, implying that significant heating of the plasma would occur. 5.3. Chiral Asymmetry and the Hall Effect As noted by Wardle (2007), the Hall effect is capable of imparting a handedness to multifluid MHD. Using our velocity solutions calculated in Section 4.2, we demonstrate this feature and the resulting chiral asymmetry in the plasma flow. Outside the shear layer the neutral velocity is just the free-stream velocity v0. Inside the shear layer, however, the velocity vector rotates. As shown in Figure 28, if the unperturbed magnetic field B0 points in the +z-direction the velocity rotates counterclockwise in the x-y plane and acquires a component vy, which is positive near the body surface. However, if B0 points in the −zdirection, the velocity rotates clockwise through the same angle in the x-y plane and acquires a component vy, which is negative close to the body. As is clearly seen in Figure 28, these rotations define a chiral asymmetry in the neutral velocity profile which depends on the relative orientation of v0 and B0. The observed chiral behavior of the flow can be explained by considering the motions of the neutral and charged particles. Since electromagnetic forces do not act directly on the neutral (bulk) plasma and the electron-neutral collisional force is neglected, the only force that can be responsible for producing the neutral velocity component vy is drag between ions and neutral particles. In the absence of the Hall effect, the ions and electrons flow together in the x-direction. However if the Hall effect is present, the Lorentz force causes the ions and electrons to move relative to one another in the x-direction. Assuming macroscopic charge neutrality (ni = ne) in the flow, this relative motion between the ions and electrons defines a component of the current density in the x-direction given by Jx ≡ eni (vix − vex ). (153) This component of the current density is associated with a gradient in the y-component of – 33 – the magnetic field through Amp`ere’s Law Jx = − c 4π dBy dz , (154) which produces a y-component of the magnetic tension force  1 4π (B · ∇) B  y ≈ B0 4π dBy dz . (155) For weakly ionized plasmas the ion-neutral drag force can generally be written as ρρiγi (vi − v) ≈ J c × B, (156) (e.g., Balbus 2011), where γi ≡ 1/ (ρi τin). Substituting Equation (154) into Equation (156) and looking at the y-component of the equation gives B0 4π dBy dz − ρρiγi (viy − vy) ≈ 0, (157) which shows that as the magnetic tension force pushes the ions in the y-direction, they drag the neutral particles with them. Thus, when the free-stream velocity of the neutral particles (v0) points in the +x-direction and the undisturbed magnetic field (B0) points in the +z-direction, magnetic tension due to the Hall effect forces the ions to move and drag the neutral particles with them in the +y-direction, causing the flow to rotate in the counterclockwise direction; while if the direction of B0 is reversed and now points in the −zdirection, magnetic tension due again to the Hall effect forces the ions and neutral particles to move in the −y-direction, causing the flow to rotate in the clockwise direction. As shown by the solutions in Section 4.2, this component of the neutral velocity outside of the plane defined by the directions of v0 and B0 can be large and thus cannot necessarily be treated as a perturbation in MHD planar shear flow analysis. 5.4. Electrodynamic Heating Although the motional electric field decreases to zero at the surface of a large body, Section 4.2 shows that it is nevertheless possible to have interior electric fields ∼ E0, i.e., electric fields comparable to those invoked in classical induction heating. It is therefore interesting to inquire whether significant heating may ensue. To emphasize the fact that the interior electric field is due ultimately to the dynamical interaction of the body with surrounding plasma, we will refer to the heating process as “electrodynamic heating.” – 34 – The benchmark for electrodynamic heating is the rate of heating by short-lived radionuclides (SLRs). We will use the rate of heating by 26Al (Urey 1955), which is known to be the dominant species (e.g., Ghosh et al. 2006). Recent work on 60Fe suggests it may also contribute to within an order of magnitude of 26Al (Castillo-Rogez et al. 2007); however we ignore its contribution here because we are only interested in making order of magnitude estimates. The volumetric heating rate is generally given by Γ26(t) = Γ26,0 exp (−0.966t/Myr), (158) where t is the time elapsed since CAI formation and Γ26,0 = 5.33 × 10−4  ρs 3000 kg m−3   xAl 1 wt%   f26,0 5 × 10−5  W m−3 , (159) where ρs is the density of the solid material, xAl is the total (including all isotopes) aluminum abundance, and f26,0 is fraction of all aluminum in 26Al at the time of CAI formation (Lee et al. 1976; MacPherson et al. 1995, 2010). This expression assumes that the half-life of 26Al is 7.17 × 105 yr and that each decay delivers 3.12 MeV of thermal energy to the solid (Castillo-Rogez et al. 2009). Also of interest is the total thermal energy deposited in the body per unit volume, ∆E26, obtained by integrating expression (158) over time. If a body is assumed to accrete instantaneously at time t = t0, then ∆E26 = 2 × 1010  ρs 3000 kg m−3   xAl 1 wt%   f26,0 5 × 10−5  exp (−0.966t0/Myr) J m−3 . (160) In comparing different heating mechanisms, ∆E26 is probably the more useful benchmark, because the time scale for 26Al to completely decay is much shorter than the time scale for thermal energy to diffuse through the body. For a body of size Lbody the diffusion time is τdif ≡ L 2 body κ > 300  Lbody 100 km Myr, (161) where κ < 10−6 m2 s −1 is the thermal diffusivity and the numerical value describes a variety of chondritic materials (Yomogida & Matsui 1983). Of course, thermal diffusion is always important in a thin layer just below the body surface, however the thickness of this layer is ≪ Lbody for times ≪ τdif (e.g., see Figure 4 of Ghosh & McSween 1998). For comparison, the rate of electrodynamic heating is Γed < σb E 2 0 , (162) – 35 – where the inequality means that this is an upper limit for reasons discussed below. If we assume that the temperature dependence of the conductivity obeys the Arrhenius law given in Equation (39), then Γed < 1.11 × 10−2  σ0 1012 s −1   v0 km s−1 2  B0 0.1 G2 exp (−1161 K/T) W m−3 . (163) Equation (163) assumes an activation energy, Ea = 0.1 eV, and σ0 consistent with measurements for a few chondritic meteorites (Schwerer et al. 1971; Brecher 1973); these are plausible values but highly uncertain. The total energy deposited in the body by electrodynamic heating in a time ∆t is just ∆Eed < 4 × 1011  σ0 1012 s −1   v0 km s−1 2  B0 0.1 G2  ∆t Myr exp (−1161 K/T) J m−3 . (164) The heating rates due to both the decay of 26Al and electrodynamic heating are plotted in Figure 29, with the unperturbed magnetic field strength once again taken to be B0 = 0.3 G. Several electrodynamic heating curves are shown corresponding to plasma flows driven by the body’s orbital motions (v0 ∼ 0.1 km/s), passing shock waves (v0 ∼ 10 km/s), and an intermediate value (v0 ∼ 1 km/s). Although larger electrodynamic heating rates are observed for faster flows such as the ones driven by shock waves, it is unlikely that these flows could produce heating inside the body comparable to 26Al because they would need to operate for time periods of order ∆t & 103 yr, A more likely scenario is that electrodynamic heating comparable to 26Al is driven by the orbital motions of the asteroids relative to the gas because this mechanism should operate for a much longer time period of order ∆t ∼ 10 Myr. This scenario can be tested quantitatively by comparing the total energy deposited by orbital motion driven electrodynamic heating to 26Al heating in a simple model asteroid. Assuming that the model asteroid body accreted instantaneously at a time t0 = 2.85 Myr after the the formation of the CAIs (Ghosh & McSween 1998) and has the density, xAl, and f26,0 values given in Equation (159), the total energy deposited by each mechanism is plotted as a function of time relative to the t0 in Figures 30–31. We consider two values of the free-stream velocity: v0 = 0.1 km corresponding to the relative velocity calculated by Weidenschilling (1997a), and v0 = 1 km/s corresponding to relative motion between the gas and a body on an eccentric orbit (Morris et al. 2012). For simplicity we neglect the temperature dependence of the body’s electrical conductivity and instead use a range of constant conductivity values, the largest of which is determined by substituting the ambient plasma temperature at R = 3 AU (Equation [52]) into Equation (39). As the plot shows, electrodynamic heating is only competitive with 26Al when the body’s electrical conductivity is of order ∼ 109 s −1 for v0 = 0.1 km/s, which is a – 36 – high but plausible value for some asteroids. For v0 = 1 km/s, comparable electrodynamic heating may be possible if the electrical conductivity of the body is of order ∼ 107 s −1 , which is probably more plausible, but still on the high end. For bodies with smaller electrical conductivities, electrodynamic heating will not be significant. It is very important to stress that the above analysis only describes possible upper limits on electrodynamic heating. Since the current which drives electrodynamic heating must leave the body, pass through a magnetized plasma sheath at each body surface, and close through the bulk plasma, the electrodynamic heating rate will not generally be given by Equation (162), but instead depend on the entire current circuit. In order to determine the limiting effects of the plasma on the current, the dynamics of the ions and electrons in the magnetic sheaths and the resistive properties of the specific path taken by the current through the bulk plasma must be taken into account. The dynamics of charged particles in magnetic sheaths have been investigated numerically using two-fluid (ions and electrons) models for simple geometries (Tskhakaya et al. 2005; Pandey et al. 2008 and references therein), and similar calculations for either side of the body are required in order to establish the effects these sheaths would have on the total current density flowing through the body. Limitations on the current due to both sheath and bulk plasma effects in a solar wind have been calculated by Srnka (1975) for the induction heating mechanism described in Sonett et al. (1970). However we leave as future work the analogous calculation for weakly ionized protoplanetary disks. Once the entire current circuit can be completely described, more accurate heating rates and temperature profiles produced inside the body by MHD shear flows can be calculated. 6. Summary In this paper we reexamined the “unipolar induction” heating mechanism developed by Sonett et al. for primitive bodies in weakly ionized protoplanetary disks. We were motivated by the fact that some asteroids once experienced thermal conditions that were conducive to a rich prebiotic chemistry, including the production of amino acids with chiral asymmetries of the type observed in terrestrial proteins. Understanding whether and how primitive bodies are/were heated is therefore an important issue for astrobiology. According to current wisdom, asteroids in the solar nebula were heated by the radioactive decay of short-lived radionuclides (SLRs), principally 26Al. However the SLR model requires a finely tuned dependence of asteroid accretion times on heliocentric distance, which have not yet been shown to produce both the correct mass distribution and heating gradient across the asteroid belt. Furthermore, the existence of other heating mechanisms may – 37 – be crucial for the viability of prebiotic chemistry in other disks, where the probability of the SLR scenario is only ∼ 10−3 –10−2 (Ouellette et al. 2010). For these reasons it seems appropriate to explore other heating mechanisms. Our principal results are: 1. We pointed out a subtle conceptual error in the induction heating mechanism as originally conceived, in consequence of which the electric field inside the body was calculated incorrectly. 2. We described the steps required to calculate the electric field correctly for bodies of arbitrary shape moving through weakly ionized plasmas of the type expected in protoplanetary disks including dust. 3. We presented a highly idealized example which demonstrated that it is possible for the electric field inside the body to vanish. Under these circumstances there is no heating. 4. We presented another highly idealized example which demonstrated that large electric fields are possible, in the sense that E is comparable to the field predicted by classical induction. 5. We demonstrated that the main effect of micron-sized dust grains on the flow and electric fields is to increase the thickness of the shear layer from ∼ 10 to ∼ 100 km. 6. We assessed the possible importance of heating by viscous dissipation and ion-neutral scattering in the shear flow. If there is no small dust in the disk, we find that heating may be significant for flows with v0 & 5 km/s at R & 3 AU when gas-phase water is present in the disk, and v0 & 1 km/s at all R values considered when water is absent. When micron-sized dust is present in the disk with an abundance of χsd = 10−4 , we find that heating is insignificant for v0 . 5 km/s at R . 8 AU when gas-phase water is present and at R . 6 AU when gas-phase water is absent. 7. We demonstrated an interesting property of shear flows around primitive bodies in protoplanetary disks, namely, that the velocity field is chirally asymmetric. As pointed out by Wardle (2007), the existence of the asymmetry is due to the Hall effect and its sense depends on the orientation of the ambient magnetic field relative to the primitive body’s velocity. The significance of the asymmetry, if any exists, remains to be determined. 8. We discovered a new “electrodynamic heating” mechanism and quantitatively compared its heating rate with the rate of heating produced by the decay of 26Al. We – 38 – found that electrodynamic heating can only be competitive with 26Al heating in asteroids with relatively high but plausible electrical conductivities of order ∼ 107 -109 s −1 depending on the value of v0; however we stress that this is an upper bound on the heating and listed a series of problems which must be solved in order to assert this unambiguously. The authors are grateful to the anonymous referee for a thorough, constructive report which greatly improved the paper. We also acknowledge helpful comments by Mike Gaffey, Satoshi Okuzumi, Heather Watson, and Stuart Weidenschilling. This work was supported by the New York Center for Astrobiology, a member of the NASA Astrobiology Institute, under grant #NNA09DA80A. A. Solution for Perpendicular Field Geometry A.1. Formal Solution Equations (125)–(128) are a set of coupled, linear ODEs for the velocity and magnetic field derivatives. In matrix form they become y ′ = A·y, (A1) where y ≡ b ′ x , b′ y , u′ x , u′ y t , (A2) is the vector of unknowns, A ≡   0 0 h11 h12 0 0 h21 h22 −ǫ 0 0 0 0 −ǫ 0 0   , (A3) is the coupling matrix, ǫ ≡ 1/4π, h ≡ 1 D2  − (1 − λh) λh −λh − (1 − λh)  , (A4) and D 2 = (1 − λh) 2 + λh 2 . (A5) – 39 – The solution of Equation (A1) is y = X m Cm φm e kmz¯ , (A6) where {km} are the four eigenvalues of A, {φm} are the corresponding eigenvectors, and the expansion coefficients {Cm} are determined by the boundary conditions. Two of the eigenvalues have positive real parts and are excluded by the boundary conditions at z = +∞. The two physically admissible eigenvalues are k− = −kR − ikI (A7) and k+ = −kR + ikI (A8) where kR and kI are functions of λh defined in eqs. (136)–(137). The corresponding eigenvectors are φ− = [−4πk−, 4πik−, 1, −i] t (A9) and φ+ = [−4πk+, −4πik+, 1, i] t . (A10) The formal solution therefore reduces to y = C− φ− e k−z¯ + C+ φ+ e k+z¯ . (A11) A.2. Boundary Conditions To apply the boundary conditions it is useful to rewrite Equation (A11) in the equivalent form y = A− φ− exp (k−ξ) + A+ φ+ exp (k+ξ) (A12) in terms of the redefined expansion coefficients A− and A+ and the dimensionless distance ξ from the body surface [see Equation (134)]. Writing out the third and fourth components gives u ′ x = A− exp (k−ξ) + A+ exp (k+ξ) (A13) and u ′ y = −iA− exp (k−ξ) + iA+ exp (k+ξ). (A14) Integrating the last two equations and applying the boundary conditions at ξ = ∞ yields the velocity solution, ux = 1 + A− k− exp (k−ξ) + A+ k+ exp (k+ξ) (A15) – 40 – and uy = − iA− k− exp (k−ξ) + iA+ k+ exp (k+ξ). (A16) The constants A− and A+ are determined by the boundary conditions at ξ = 0 with the result ux = 1 − 1 2 exp (k−ξ) − 1 2 exp (k+ξ) (A17) uy = i 2 exp (k−ξ) − i 2 exp (k+ξ). (A18) Writing out k± in terms of kR and kI and taking the real parts gives eqs. (130)–(131). With ux and uy known, the magnetic field solution can be found by integrating b ′ x = 4πu′′ x (A19) and b ′ y = 4πu′′ y (A20) and applying the boundary conditions, (118) and (119). Expressing the result in dimensional units gives Equations (132) and (133). REFERENCES Bai, X., & Goodman, J. 2009, ApJ, 701, 737 Balbus, S.A. 2011, in Physical Processes in Circumstellar Disks around Young Stars, ed. P.J.V. Garcia (Chicago: Univ. Chicago) 237 Balbus, S.A., & Hawley, J.F. 1991, ApJ, 376, 214 Blasius, H. 1908, Zeitschrift f¨ur angewandte Mathematik und Physik, 56, 1 Birnstiel, T., Ormel, C.W., & Dullemond, C.P. 2011, A&A, 525, A11 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 Brecher, A. 1973, Meteoritics, 8, 17 Castillo-Rogez, J.C., Matson, D.L., Sotin, C., Johnson, T.V., Lunine, J.I., & Thomas, P.C. 2007, Icarus, 190, 179 Castillo-Rogez, J., Johnson, T.V., Lee, M.H., Turner, N.J., Matson, D.L., & Lunine, J. 2009, Icarus, 204, 658 – 41 – Chernoff, D.F. 1987, ApJ, 312, 143 Chyba, C., & Sagan, C. 1992, Nature, 355, 125 Cronin, J.R., & Pizzarello, S. 1997, Science, 275, 951 Dalgarno, A. 2006, PNAS, 103, 12269 Desch, S.J., & Connolly, H.C. Jr. 2002, Meteoritics & Planetary Science, 37, 183 Desch, S.J., & Mouschovias, T.Ch. 2001, ApJ, 550, 314 Draine, B.T. 2011, Physics of the Interstellar and Intergalactic Medium, (Princeton: Princeton University Press) Draine, B.T., Roberge, W.G., & Dalgarno, A. 1983, ApJ, 264, 485 Draine, B.T., & Salpeter, E.E. 1979, ApJ, 231, 77 Dullemond, C.P., & Dominik, C. 2005, A&A, 434, 971 Dzyurkevich, N., Turner, N.J., Henning, T., & Kley, W. 2013, ApJ, 765, 114 Furlan, E., et al. 2011 ApJS, 195, 3 Ghosh, A., & McSween, H.Y. Jr. 1998, Icarus, 134, 187 Ghosh, A., Weidenschilling, S.J., McSween, H.Y. Jr., & Rubin, A. 2006, in Meteorites and the Early Solar System II, ed. D.S. Lauretta & H.Y. McSween Jr. (Tucson: Univ. Arizona), 555 Glavin, D.P., Callahan, M.P., Dworkin, J.P., & Elsila, J.E. 2011, Meteoritics & Planetary Science, 45, 1948 Glavin, D.P., & Dworkin, J.P. 2009, Meteoritics & Planetary Science Supplement, 5009 Grimm, R.E., & McSween, H.Y., Jr. 1993, Science, 259, 653 Hartmann, J. 1937, Kongelige Danske Videnskabernes Selskab Matematisk-Fysiske Meddelelser, 15, 1 Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35 Hood, L.L., Ciesla, F.J., Artemieva, N.A., Marzari, F., & Weidenschilling, S.J. 2009, Meteoritics & Planetary Science, 44, 327 – 42 – Igea, J., & Glassgold, A.E. 1999, ApJ, 518, 848 Ip, W.-H., & Herbert, F. 1983, Moon & Planets, 28, 43 Jackson, J.D. 1975, Classical Electrodynamics, 2nd Ed. (New York: Wiley) Keil, K. 2000, Planet. Space Sci., 48, 887 Lee, T., Papanastassiou, D.A., & Wasserburg, G.J. 1976, Geophys. Res. Lett., 3, 41 Lee, T., Papanastassiou, D.A., & Wasserburg, G.J. 1977, ApJ, 211, L107 MacPherson, G.J., Bullock, E.S., Janney, P.E., Kita, N.T., Ushikubo, T., Davis, A.M., Wadhwa, M., & Krot, A.N. 2010, ApJ, 711, L117 MacPherson, G.J., Davis, A.M., & Zinner, E.K. 1995, Meteoritics, 30, 365 Martin, R.G., & Livio, M. 2012, MNRAS, 425, L6 McElroy, D., Walsh, C., Markwick, A.J., Cordiner, M.A., Smith, K., & and Millar, T.J. 2013, A&A, 550, A36 McKinnon, W.B. 1989, Nature, 340, 343 McSween, H.Y., Jr., Ghosh, A., & Weidenschilling, S.J. 2002, Meteoritics & Planetary Science, 37, A98 Miura, H., & Nakamoto, T. 2006, ApJ, 651, 1272 Moresco, P., & Alboussi`ere, T. 2004, Journal of Fluid Mechanics, 504, 167 Morris, M.A., Boley, A.C., Desch, S.J., & Athanassiadou, T. 2012, ApJ, 752, 27 Morris, M.A., Desch, S.J., & Ciesla, F.J. 2009, ApJ, 691, 320 Neufeld, D.A., & Kaufman, M.J. 1993, ApJ, 418, 263 Okuzumi, S. 2009, ApJ, 698, 1122 Ouellette, N., Desch, S.J., & Hester, J.J. 2010, ApJ, 711, 597 Padovani, M., Galli, D., & Glassgold, A.E. 2009, A&A, 501, 619 Pandey, B.P., Samarian, A., & Vladimirov, S.V. 2008, Plasma Physics and Controlled Fusion, 50, 055003 – 43 – Sano, T., Miyama, S.M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 Sano, T., & Stone, J.M. 2002, ApJ, 577, 534 Schaefer, J. 2010, Chem Phys, 374, 15 Schwartz, K., Sonett, C.P, & Colburn, D.S. 1969, The Moon, 1, 70 Schwerer, F.C., Nagata, T., & Fisher, R.M. 1971, The Moon, 2, 408 Shakura, N.I., & Sunyaev, R.A. 1973, A&A, 24, 337 Sonett, C.P., & Colburn, D.S. 1968, Phys. Earth Planet. Interiors, 1, 326 Sonett, C.P., Colburn, D.S., Schwartz, K., & Keil, K. 1970, Ap&SS, 7, 446 Srnka, L.J. 1975, Ap&SS, 36, 177 Thess, A., Krasnov, D., Boeck, T., Zienicke, E., Zikanov, O., Moresco, P., & Alboussi`ere, T. 2007, Gesellschaft f¨ur Angewandte Mathematik und Mechanik Mitteilungen, 30, 125 Tskhakaya, D.D., Shukla, P.K., Eliasson, B., & Kuhn, S. 2005, Physics of Plasmas, 12, 103503 Turner, N.J., & Sano, T. 2008, ApJ, 679, L131 Umebayashi, T. & Nakano, T. 2009, ApJ, 690, 69 Urey, H.C. 1955, PNAS, 41, 127 Wardle, M. 1998, MNRAS, 298, 507 Wardle, M. 2007, Ap&SS, 311, 35 Wasserburg, G.J., Lee, T., & Papanastassiou, D.A. 1977, Meteoritics, 12, 377 Weidenschilling, S.J. 1977a, MNRAS, 180, 57 Weidenschilling, S.J. 1977b, Ap&SS, 51, 153 Williams, J.P. 2012, Meteoritics & Planetary Science, 47, 1915 Yomogida, K., & Matsui, T. 1983, J. Geophys. Res., 88, 9513 This preprint was prepared with the AAS LATEX macros v5.2. – 44 – Table 1. Parameters in Transport Coefficient Calculation Symbol Definition Value(s) ui Average†of mean ion thermal speed Equation 60 βr Average† recombination rate coefficient Equation 61 ζ Average††ionization rate See text. Lxr Stellar x-ray luminosity 2 × 1030 erg s−1 ζra Rate of ionization by radionuclides 1.4 × 10−22 s −1 nH Number density of H nuclei From disk model T Gas temperature From disk model χsd Mass fraction, small grains ≤ 10−4 Df Dust fractal dimension 3 µi Ion mass/H atom mass 24 Si Sticking efficiency, ion-grain collision 1 Se Sticking efficiency, e−-grain collision 0.3 ρs Bulk density of dust material 1.4 g cm−3 †Over all ionic species. See text. ††Over all gas-phase species. See text. – 45 – Table 2. Parameters in Dzyurkevich et al. (2013) Model Symbol Definition Value(s) χsd Abundance of small dust 10−4 ui Average†of mean ion thermal speed Equation 60 βrMg+ Mg+ recombination rate coefficient 3 × 10−11/ √ T cm3 s −1 βrHCO+ HCO+ recombination rate coefficient 3 × 10−6/ √ T cm3 s −1 Lxr Stellar x-ray luminosity 1029 erg s−1 ζcr Rate of ionization by cosmic rays See Dzyurkevich et al. (2013) ζxr Rate of ionization by x rays See Dzyurkevich et al. (2013) ζra Rate of ionization by radionuclides 7 × 10−19 (χsd/10−2 ) s−1 nH Number density of H nuclei From disk model T Gas temperature From disk model a0 Monomer radius 0.1 µm N Number of monomers/aggregate 400 Df Dust fractal dimension 2 µmg Mg mass/H atom mass 24 µHCO HCO mass/H atom mass 29 Si Sticking efficiency, ion-grain collision 1 Se Sticking efficiency, e−-grain collision 0.3 ρs Bulk density of dust material 1.4 g cm−3 nMg Total number density of Mg nuclei/nH 2 × 10−9†† B Magnetic field strength See Dzyurkevich et al. (2013) †Over Mg+ and HCO+. See text. ††100 times greater than the abundance quoted in Appendix B of Dzyurkevich et al. (2013). – 46 – 2 4 6 8 1013 1014 R(AU) ng !cm −3 " 2 4 6 8 100 150 200 250 R(AU) T(K) Fig. 1.— The midplane number density of neutral particles (left) and temperature (right) plotted vs. distance R from the central star. – 47 – 2 4 6 8 10−21 10−20 10−19 10−18 10−17 R(AU) ζ ! s −1 " Fig. 2.— The midplane ionization rate per gas particle as a function of R showing the contributions of cosmic rays (filled circles), radioactivity (dashed), x rays (dotted) and the total rate (solid). – 48 – 2 4 6 8 10 10−16 10−14 10−12 10−10 10−8 10−6 10−4 a(µm) χsd a0 = 0.1 µm ac = 1 mm χd = 0.01 γ = 0.5 γ = 1.0 γ = 1.5 γ = 1.75 Fig. 3.— Plot of χsd vs. grain radius, a (see Equation [73]) for a few values of the exponent in Equation (72). – 49 – 2 4 6 8 10−18 10−16 10−14 10−12 10−10 R(AU) Abundance 2 4 6 8 10−18 10−16 10−14 10−12 10−10 R(AU) Abundance χsd = 10−4 a = 0.1 µm χsd = 10−4 a = 1 µm Fig. 4.— The abundances of ions (filled circles), electrons (solid curve), and charge contained by single-sized dust grains (dotted curve = postively charged, dashed curve = negatively charged; See Equations [74]–[75]) relative to nH in the disk midplane are plotted versus distance from the central star. The assumed abundance of small dust and grain radii are indicated for each case. – 50 – 2 4 6 8 10−14 10−13 10−12 10−11 10−10 10−9 R(AU) Electron Abundance no dust 10 µm 1 µm 0.1 µm χsd = 10−4 Fig. 5.— Variation of the electron abundance, xe, with R including single-size dust grains with three different radii indicated in the legend. The abundance of small dust is χsd = 10−4 . – 51 – 2 4 6 8 10−14 10−13 10−12 10−11 10−10 10−9 R(AU) Electron Abundance no dust 10 µm 1 µm 0.1 µm χsd = 10−6 Fig. 6.— As in Figure 5 but for χsd = 10−6 . – 52 – 2 4 6 8 10−8 10−6 10−4 10−2 100 102 R(AU) Hall Parameter B = 0.1 G e e*i i d Fig. 7.— Hall parameters as a function of R for the ions (i) and electrons (e) and a dust grain with a = 1 µm and Z = ±1 (d). The curve labeled i*e is the product of the ion and electron Hall parameters. Values are for the disk midplane. – 53 – 2 4 6 8 10−8 10−6 10−4 10−2 100 102 R(AU) Hall Parameter B = 1 G e e*i i Fig. 8.— As in Figure 7 but for B = 1 G. – 54 – 2 4 6 8 1013 1014 1015 1016 R(AU) Di ffusivity !cm2 s −1 " 2 4 6 8 1013 1014 1015 1016 R(AU) Di ffusivity !cm2 s −1 " χsd = 0 B = 0.1 G χsd = 0 B = 1 G Fig. 9.— Variation of the Hall (dashed curve), Ohmic (solid), and ambipolar (dash-dotted) diffusivities with R for a dust-free plasma. Results are plotted for two values of the magnetic field B. Values are for the disk midplane. – 55 – 2 4 6 8 1015 1016 1017 1018 1019 R(AU) Di ffusivity !cm2 s −1 " 2 4 6 8 1015 1016 1017 1018 1019 R(AU) Di ffusivity !cm2 s −1 " χsd = 10−4 B = 0.1 G χsd = 10−4 B = 1 G Fig. 10.— As in Figure 9 but for χsd = 10−4 . The dust is a single-size distribution with a = 1 µm. – 56 – 2 4 6 8 10−18 10−16 10−14 10−12 R(AU) Abundance 2 4 6 8 10−18 10−16 10−14 10−12 R(AU) Abundance Mg+ HCO+ Total Fig. 11.— Left panel: As in Figure 4, but for the parameters described in Table 2. Right panel: The number densities of Mg+, HCO+, and their sum relative to the number density of hydrogen nuclei in the disk. Both panels are for the model described in Dzyurkevich et al. (2013). – 57 – 2 4 6 8 1017 1018 1019 1020 1021 1022 1023 R(AU) Di ffusivity (cm2 s −1 ) Fig. 12.— As in Figure 9, but for the model described in Dzyurkevich et al. (2013). – 58 – xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx z x x x x x x 2W v0 Fig. 13.— Flow past an infinite slab of thickness 2W. At the upper body surface the outward normal is +zˆ. Far from the body the velocity is v0 = v0 xˆ. We consider two cases, where the undisturbed magnetic field far from the body is parallel (B0 = B0 yˆ, “parallel field geometry”) and perpendicular (B0 = B0 ˆz, “perpendicular field geometry”) to the body surfaces. – 59 – 2 4 6 8 0 0.25 0.5 0.75 1 R (AU) λ h B0 = 0.1 G B0 = 1 G Fig. 14.— The parameter λh defined in Equation (129) for a dust-free plasma is plotted vs distance from the central star. The two curves correspond to different B0 values as indicated. – 60 – 2 4 6 8 0.25 0.5 0.75 1 R(AU) λ h B0 = 0.1 G B0 = 1 G a0 = 1 µm χsd = 10−4 Fig. 15.— As in Figure 14, but with χsd = 10−4 and a = 1 µm. – 61 – 0 0.25 0.5 0.75 1 0.1 0.2 0.3 λh k R,kI kR kI Fig. 16.— Dimensionless quantities kR and kI defined in eqs. (136) and (137) plotted vs. the dimensionless parameter λh. – 62 – 0 20 40 60 80 100 −0.5 0 0.5 z (km) Ex/E0 Hall motional sum R = 3 AU B0 = 0.3 G χsd = 0 Ohm+ambipolar Fig. 17.— The x-component of the electric field, Ep, in units of E0 as a function of height z above the body surface. The values of R and B0 are indicated. The motional, Ohm+ambipolar, and Hall fields are indicated. The curve labeled “sum” is the total electric field. – 63 – 0 20 40 60 80 100 −0.5 0 0.5 1 z (km) Ey/E0 sum Hall motional R = 3 AU B0 = 0.3 G χsd = 0 Ohm+ambipolar Fig. 18.— As in Figure 17 but for the y component of the electric field. – 64 – 0 500 1000 1500 2000 −0.5 0 0.5 z (km) Ex/E0 R = 3 AU B0 = 0.3 G χsd = 10−4 a = 1 µm Hall motional sum Ohm+ambipolar Fig. 19.— As in Figure 17 but for χsd = 10−4 . The dust is a single-size distribution with a = 1 µm. – 65 – 0 500 1000 1500 2000 0 0.5 1 z (km) Ey/E0 R = 3 AU B0 = 0.3 G χsd = 10−4 a = 1 µm sum Hall motional Ohm+ambipolar Fig. 20.— As in Figure 18 but for χsd = 10−4 . The dust is a single-size distribution with a = 1 µm. – 66 – 100 200 300 400 500 10−7 10−6 10−5 10−4 T (K) Λmr (erg cm −3 s −1 ) 100 200 300 400 500 10−7 10−6 10−5 10−4 10−3 T (K) Λmr (erg cm −3 s −1 ) v0 = 1 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 1.1 × 10−4 χsd = 0 v0 = 5 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 1.1 × 10−4 χsd = 0 Fig. 21.— The rate of radiative cooling by H2, CO, and H2O molecules plotted vs. the gas temperature, computed using the cooling rate function of Neufeld & Kaufman (1993). Each curve represents the cooling rate a specific distance from the central star as follows: R = 2 AU (solid), R = 3 AU (dashed), R = 4 AU (dashed-dotted), R = 6 AU (solid with open circles), and R = 8 AU (dashed with closed circles). The free-stream flow velocity, ambient magnetic field, molecular abundances, and small grain abundance are indicated. – 67 – 100 200 300 400 500 10−8 10−7 10−6 10−5 10−4 T (K) Λmr (erg cm −3 s −1 ) 100 200 300 400 500 10−7 10−6 10−5 10−4 T (K) Λmr (erg cm −3 s −1 ) v0 = 1 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 1.1 × 10−4 χsd = 10−4 a = 1 µm v0 = 5 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 1.1 × 10−4 χsd = 10−4 a = 1 µm Fig. 22.— As in Figure 21 but for χsd = 10−4 . The dust is a single-size distribution with a = 1 µm. – 68 – 100 200 300 400 500 10−11 10−10 10−9 10−8 10−7 10−6 T (K) Λmr (erg cm −3 s −1 ) 100 200 300 400 500 10−11 10−10 10−9 10−8 10−7 10−6 T (K) Λmr (erg cm −3 s −1 ) v0 = 1 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 0 χsd = 0 v0 = 5 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 0 χsd = 0 Fig. 23.— As in Figure 21, but with nH2O/nH = 0. – 69 – 100 200 300 400 500 10−11 10−10 10−9 10−8 10−7 10−6 T (K) Λmr (erg cm −3 s −1 ) 100 200 300 400 500 10−11 10−10 10−9 10−8 10−7 10−6 T (K) Λmr (erg cm −3 s −1 ) v0 = 1 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 0 χsd = 10−4 a = 1 µm v0 = 5 km/s B0 = 0.3 G nCO /nH = 7 × 10−5 nH2O/nH = 0 χsd = 10−4 a = 1 µm Fig. 24.— As in Figure 22, but with nH2O/nH = 0. – 70 – 2 4 6 8 10−9 10−8 10−7 10−6 10−5 10−4 10−3 R(AU) Λdr (erg cm −3 s −1 ) a0 = 1 µm χsd = 10−4 Fig. 25.— The cooling rate due to thermal emission from dust grains is plotted vs. distance R from the central star. The dust is a single-size distribution with a = 1 µm and χsd = 10−4 . – 71 – 2 4 6 8 10−9 10−8 10−7 R(AU) Γ !erg cm −3 s −1 " 2 4 6 8 10−8 10−7 10−6 R(AU) Γ !erg cm −3 s −1 " B0 = 0.3 G v0 = 5 km s−1 χsd = 0 B0 = 0.3 G v0 = 1 km s−1 χsd = 0 Fig. 26.— The volumetric heating rate Γ due to viscous dissipation (solid curve) and ionneutral scattering (dashed curve) in the plasma at the body surface (ξ = 0) are plotted versus distance from the central star. The values of the free-stream velocity, ambient magnetic field, and small dust abundance are indicated. – 72 – 2 4 6 8 10−10 10−9 R(AU) Γ !erg cm −3 s −1 " 2 4 6 8 10−9 10−8 10−7 R(AU) Γ !erg cm −3 s −1 " B0 = 0.3 G v0 = 5 km s−1 χsd = 10−4 a = 1 µm B0 = 0.3 G v0 = 1 km s−1 χsd = 10−4 a = 1 µm Fig. 27.— As in Figure 26 but for χsd = 10−4 . The dust is a single-size distribution with a = 1 µm. – 73 – 0 0.2 0.4 0.6 0.8 1 1.2 −0.4 −0.2 0 0.2 0.4 vx/v0 v y/v 0 B0 in the +z-direction B0 in the −z-direction B0 = 0.3 G v0 = 1 km/s R = 3 AU Fig. 28.— The neutral bulk plasma velocity profile for the flow calculated in Section 4.2 is plotted in the x-y plane for two different orientations of the magnetic field with v0 = v0 xˆ and the indicated parameters. The solid and dot dashed curves correspond to B0 = B0zˆ and B0 = −B0ˆz respectively. Each point on the curves represents the velocity at a certain height above the body surface, with the black circle at the point (1,0) indicating the edge of the shear layer (ξ = ±∞, v = v0), and the intercept at (0,0) indicating the body surface (ξ = 0, v = 0). – 74 – 2 4 6 8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 R(AU) Γ !W m −3 " 2 Myr 4 Myr 6 Myr t = 0 10 km/s 1 km/s 0.1 km/s Fig. 29.— Heating rates vs. heliocentric distance. Dashed curves: 26Al heating, where the labels indicate the time since CAI formation. Solid curves: electrodynamic heating, where the labels indicate v0. The electrodynamic heating rate depends on heliocentric distance via the temperature dependence of the conductivity (see text). – 75 – 0 2 4 6 8 10 x 106 100 102 104 106 108 1010 1012 ∆t (yr) ∆ E (J/m − 3 ) 26Al Heating Electrodynamic Heating (σb = 7.63 × 108 s −1 ) Electrodynamic Heating (σb = 106 s −1 ) Electrodynamic Heating (σb = 104 s −1 ) B0 = 0.3 G v0 = 0.1 km/s R = 3 AU Fig. 30.— The total amount of energy deposited in our model asteroid body per volume by both 26Al and electrodynamic heating driven by the orbital motions (v0 = 0.1 km/s, B0 = 0.3 G) of the asteroid is plotted versus time relative to the time of accretion. Curves are shown for several different plausible values of the asteroid’s electrical conductivity. – 76 – 0 2 4 6 8 10 x 106 102 104 106 108 1010 1012 1014 ∆t (yr) ∆ E (J/m − 3 ) 26Al Heating Electrodynamic Heating (σb = 7.63 × 108 s −1 ) Electrodynamic Heating (σb = 106 s −1 ) Electrodynamic Heating (σb = 104 s −1 ) B0 = 0.3 G v0 = 1 km/s R = 3 AU Fig. 31.— As in Figure 30 but with v0 = 1 km/s.