© ESO, 2017

1. Introduction

Tiny spherules that crystallised from molten rock are abundant in primitive chondrite meteorites from the asteroid belt (Krot et al. 2009). These chondrules of typical sizes between 0.1 and 1 mm comprise approximately 60–80% of the mass of the ordinary and enstatite chondrite classes, mixed with smaller metal grains and microscopic matrix particles (Dodd 1976). The chondrule abundance is more varied in the carbonaceous chondrites, ranging from approximately 50% in CO and CV chondrites to 0% in the CI chondrites (Scott & Krot 2003). The lack of chondrules in the CI chondrites could nevertheless be the result of chondrule destruction during extensive aqueous alteration of the parent body (Endress et al. 1996).

The sheer abundance of chondrules implies that these particles played a key role in the formation and/or growth of planetesimals in the solar protoplanetary disc. Chondrule-sized particles can concentrate into dense clumps in the gaseous protoplanetary disc either by preferential concentration between turbulent eddies near the dissipation scale of the turbulent gas (Cuzzi et al. 2008) or by the streaming instability that concentrates particles into dense filaments through their collective drag force on the gas (Youdin & Goodman 2005; Johansen & Youdin 2007; Bai & Stone 2010). The streaming instability concentrates particles as small as chondrules at the location of the asteroid belt if the turbulence in the protoplanetary disc is sufficiently weak to allow mm-sized particles to sediment (Carrera et al. 2015; Yang et al. 2017). Chondrules also have ideal sizes to become accreted onto young planetesimals and protoplanets. Johansen et al. (2015) showed that the observed size distribution of asteroids (particularly the steeply declining differential size distribution for asteroids larger than 100 km in diameter and the transition to a shallower size distribution above 400 km in diameter) can arise due to the accretion of chondrules with sizes similar to those measured in ordinary chondrites. Chondrule accretion also drives the continued growth to Mars-sized protoplanets, both in the asteroid belt and in the terrestrial planet region. Therefore understanding the origin of chondrules is imperative for understanding the formation of asteroids and terrestrial planets.

The identification of radiogenic ²⁶Mg in the CAI component (calcium-aluminium-rich inclusions) of chondrites, produced through the decay of the short-lived radionuclide ²⁶Al (Lee et al. 1977), provides a strong connection between the solar protoplanetary disc and nearby supernovae in the Sun’s birth cluster. The ratio of ²⁶Al to its stable sibling ²⁷Al is narrowly peaked at a value of approximately 5 × 10^-5 in the CAIs found in the CV chondrites (Amelin et al. 2002; Connelly et al. 2012). Age determination of chondrules based on radiogenic ²⁶Mg yields chondrule ages that are systematically several million years younger than the CAIs (Villeneuve et al. 2009). However, direct measurement of ages of individual chondrules, based on the long-lived Pb-Pb system, indicates that the amount of ²⁶Al was originally lower in the chondrule-forming region compared to the region where CAIs condensed out (Connelly et al. 2012; Larsen et al. 2016). This could reflect a real heterogeneous distribution of ²⁶Al in the solar protoplanetary disc or that the dust grain carrier of ²⁶Al was not incorporated into chondrules (M. Bizzarro, priv. comm.). The abundance of ²⁶Al in the asteroid belt may have been as much as five times lower than the value inferred from CAIs.

Decay of ²⁶Al is believed to have heated some meteorite parent bodies to differentiate completely while others (the parent bodies of the chondrites) acquired an onion layering of internal heating degree below the melting temperature. The measured systematic magnetisation of some chondrites (Weiss & Elkins-Tanton 2013), supported by theoretical models of chondrule accretion on existing planetesimals (Johansen et al. 2015), indicates that some differentiated bodies have outer layers of primitive chondrules. A reduction in the primordial amount of ²⁶Al in the asteroid belt bodies, as indicated by the older Pb-Pb ages of chondrules (Connelly et al. 2012), strongly reduces the degree of internal heating and hence the expected relative frequency of fully differentiated, partially differentiated and primitive bodies in the asteroid belt (Larsen et al. 2016).

Chondrules appear to have crystallised after flash-heating of the unknown chondrule precursor particles (likely fluffy dust aggregates of approximately mm sizes). An additional component of small dust grains in meteorites, the matrix, has been argued to have condensed out of a solar composition gas in connection with chondrule formation, mainly based on the similar cooling times needed to explain both chondrule and matrix crystal growth (Scott & Krot 2005). The mechanism that led to widespread melting and partial evaporation of chondrule precursors is one of the most troubling unsolved problems of meteoritics and planet formation. Chondrule precursors must have been either heated and cooled very quickly or heated in the close vicinity of other chondrule precursors, in order to build up the saturation vapour pressure needed to stabilise the liquid rock phase, a value much in excess of the ambient pressure of the solar protoplanetary disc. Heating a dust aggregate in the nominal pressure of the solar protoplanetary disc will result in complete sublimation before the melting temperature is reached, unless the heating and cooling were extremely rapid. However, such rapid heating and cooling yield the wrong texture of the synthetic chondrules produced in experiments (Hewins et al. 2005).

The requirement for chondrule precursors to heat to the liquid phase without vanishing by sublimation, combined with the observed lack of fractionation of volatiles in chondrules, gives a number density constraint of minimum 10 chondrule precursors per m³, or a minimum mass density of approximately 10^-5 kg m^-3 (Cuzzi & Alexander 2006). Retention of the volatile element sodium during chondrule formation, inferred from the lack of internal sodium gradients within large olivine crystals, implies even higher densities, between 10^-3 and 1 kg m^-3 (Alexander et al. 2008). The solar protoplanetary disc had a characteristic density of 10^-7 kg m^-3 at the distance of the asteroid belt in its pristine stages and likely an order of magnitude lower after a few million years of evolution (Bitsch et al. 2015). The Roche density – above which particle clumps are gravitationally bound against the tidal force from the central star – is ρ_R ≈ 2.8 × 10^-5 kg m^-3 at a distance of 2.5 AU from a solar-mass star. Chondrules therefore, appear to have heated in an environment of 100–10 000 times the Roche density. Alexander & Ebel (2012) showed that matching the abundance of compound chondrules, fusioned during cooling at the short-lived viscous stage, also requires such ultrahigh densities. These conditions are hard to reconcile with any nominal environment of the protoplanetary disc, except perhaps very close to the star (0.1–0.5 AU) where the Roche density is 2–4 orders of magnitude higher than in the asteroid belt.

A plethora of theories have been formulated to explain the formation of chondrules. The most popular is the shock model where chondrule precursors are heated in shock waves that originate in the protoplanetary disc itself or from planetesimals that plow through the gas on eccentric orbits. These models can match the right heating and cooling times (Desch & Connolly 2002; Ciesla et al. 2004), but do not explain the high precursor densities needed to sustain the liquid phase and maintain volatile abundances. Current sheets driven by the magnetic field in the protoplanetary disc provide another mechanism for heating chondrule precursors, but require a very high starting temperature of more than 800 K to reach the relevant temperature for melting chondrule precursors (McNally et al. 2013).

The planetesimal collision model posits that chondrules are secondary in nature and arise as debris from planetesimal collisions. The collision model comes in two flavours: collisions between pre-melted planetesimals (Sanders & Taylor 2005) and collisions between solid planetary embryos (Johnson et al. 2015). The first actually taps into the decay energy of ²⁶Al, as the planetesimals melted through this energy input. The splash ejecta is then envisioned to crystallise as mm-sized chondrules (Asphaug et al. 2011). In the embryo model the melting happens in jets that are ejected from the impact between solid embryos. Collisions can in principle explain the high densities of chondrule precursors as the collision debris crystallises, as well as the retention of volatiles (Dullemond et al. 2014). Producing enough chondrules to form all the chondritic bodies that must have existed in the primordial asteroid belt is nevertheless a challenge. In the model of Johnson et al. (2015), creating even 1% of the primordial asteroid belt in chondrules required a total planetesimal mass in excess of three times the mass of the primordial asteroid belt, as inferred from extrapolating between the solid components of the terrestrial planets and the giant planets (Hayashi 1981). Most of the chondrules in the simulations of Johnson et al. (2015) were in fact produced in the terrestrial planet formation region where collision frequencies are higher.

Lightning as a mechanism for heating chondrule precursors was proposed by Whipple (1966). Lightning occurs in terrestrial volcanic ash plumes and causes large ash particles to melt to form chondrule-like spherules (Genareau et al. 2015). Dust aggregates hit directly by an electric discharge will nevertheless break rather than melt (Güttler et al. 2008), but exposure to the photons from the lightning results in melt spherules with size distributions similar to chondrules (Eisenhour & Buseck 1995; Poppe et al. 2010). In fact, terrestrial lightning processes approximately 70% of its total energy through photons (Eisenhour et al. 1994). Radiation is attractive for chondrule melting since it naturally introduces a smallest chondrule size, in agreement with the observed log-normal or Weibull size distributions of chondrules (Dodd 1976; Teitler et al. 2010), from the inefficient photon absorption by small dust aggregates.

Obtaining lightning discharge in protoplanetary discs requires (a) a charging mechanism; (b) a charge separation mechanism and (c) low enough gas conductivity to allow the electric field to build up to a value where the current collisionally ionises the hydrogen molecule. Charging and charge separation are natural consequences of vertical or radial sedimentation, since large pebbles are charged triboelectrically when colliding with dust grains (Desch & Cuzzi 2000; Muranushi 2010). Pilipp et al. (1998) nevertheless found that vertical sedimentation of pebbles is many orders of magnitude too slow to compete with the neutralising current. Desch & Cuzzi (2000) considered charge separation as chondrule-sized particles concentrate at km scales in strong disc turbulence. They found that the gas ionisation by decay of ²⁶Al must be reduced by several orders of magnitude in order to build up the break-down electric field.

The elevated densities needed for chondrule formation make an intriguing link to planetesimal formation by particle concentration by the streaming instability and gravitational collapse (Johansen et al. 2014). In this paper we therefore, explore lightning discharge and chondrule formation in dense pebble environments. We show that secondary electrons released as positrons from the decay of ²⁶Al hammer through dust surfaces leads to positive charging of the pebbles. Charge separation by radial drift is nevertheless very slow inside dense pebble filaments and we find that the pebble current cannot compete with the neutralising current, even if abundance of ²⁶Al, and hence the gas conductivity, is decreased by an order of magnitude. Therefore we consider instead charge separation driven by the current of positrons released by the decay of ²⁶Al. Positrons stream out of dense pebble structures and cause a large-scale negative charging of both circumplanetesimal discs and the sedimented mid-plane layer of pebbles in the protoplanetary disc. This charging is related to a direct-charging nuclear battery where a radioactive source charges a capacitor. The electric field can build up values where the electrons ionise the hydrogen molecule and cause lightning discharge. We construct simple models of the circumplanetesimal pebble disc and the protoplanetary disc environments. We show that the efficiency of converting the energy of the decay of ²⁶Al to photons released in lightning discharge can be high in both these settings, although many uncertainties remain in the modelling of circumplanetesimal discs and the charging of filamentary pebble structures with density variation on many scales.

The paper is organised as follows. In Sect. 2 we review the physics of lightning discharge and show by simple analytical arguments that charge separation driven by differential drift of the pebbles relative to the dust and gas is orders of magnitude slower than the neutralising current. This conclusion is supported by thorough calculations of the equilibrium charge on gas, dust and pebbles presented in Appendices A–C, where we consider pebble charging by emission and absorption of positrons from the decay of ²⁶Al and release of secondary electrons after positron absorption (this mechanism is elaborated on in Appendix D). In Sect. 3 we analyse the positron current emanating from pebble filaments and sheets. We show that the positron current is higher than the neutralising current when the dust-to-gas ratio is above approximately 50. We discuss in Sect. 4 the operation of the nuclear battery effect in circumplanetesimal discs and in the sedimented mid-plane layer of pebbles. In the following Sect. 5 we demonstrate that the cooling rates of chondrules formed in dense pebble regions agrees with the long cooling times inferred from crystal growth experiments. Finally we summarise our results in Sect. 6.

2. The problem with lightning discharge

In this section we review the physics of lightning discharge and demonstrate how charge separation by pebble sedimentation or differential radial drift between pebbles charged oppositely of dust grains and gas drives currents in protoplanetary discs that are many orders of magnitudes lower than the neutralising gas current.

2.1. Lightning discharge

Lightning discharge requires the build up of an electric field E strong enough that the accelerated electrons ionise the hydrogen molecule. This leads to a cascade effect where the additional electrons increase the conductivity and hence the charge separation is neutralised in thin discharge channels. The breakdown electric field is given in Desch & Cuzzi (2000) as (1)Here, we have normalised to a nominal gas density value ρ_g for the asteroid belt region in the primordial state of the solar protoplanetary disc and to an ionisation potential of e_ion = 60 eV. The ionisation energy of H₂ is 15.6 eV, but Desch & Cuzzi (2000) argued that impact energies below 60 eV lead to dissociation of the hydrogen molecule. Okuzumi & Inutsuka (2015) used a distribution function of electron energies and found that the breakdown of the hydrogen molecule occurs already at electron energies around 5–10 eV; hence we multiplied Eq. (1) by a factor f_BD that we will assume in this paper to be 0.1.

The breakdown electric field can only be built up if the neutralising current, driven by growing electric field, remains smaller than the current driven by the motion of the charged solid particles. The gas has an intrinsic conductivity that is a function of the density and mobility of ions and electrons. Ionisation in the dense mid-plane layer is dominated by the decay of short-lived radionuclides, most importantly by release of positrons in the decay of ²⁶Al (Desch & Cuzzi 2000). External sources of ionisation, such as cosmic rays, can be ignored in the dense mid-plane in the asteroid formation region.

2.2. Ion and electron densities

The gas conductivity depends strongly on the number densities of ions and electrons, which are in turn dictated by adsorption onto dust grains. We consider throughout this paper for simplicity a two-component model for the solid particles. Component 1 consists of microscopic dust grains that are the main absorbers of electrons and ions, while component 2 consists of macroscopic (dust aggregate) pebbles that can form chondrules after flash heating. The number densities of electrons, ions, positrons, small grains and pebbles are n_e, n_i, n_β, n₁, and n₂. The two species of solid particles have masses m₁ and m₂, radii R₁ and R₂, and charges (in units of the elementary charge e) q₁ and q₂. The gas is characterised by a neutral density n_n = ρ_g/m_n, where m_n is the neutral (hydrogen molecule) mass. Cross sections are denoted σ, with single subscript used to indicate the target species or double subscript to indicate that both collision partners matter for the cross section.

Positrons released by the decay of ²⁶Al lose their energy to ionisation of the hydrogen molecule (with ionisation potential E_ion) and to stopping in solid particles. The positron kinetic energy has an average value of T_β = 0.66 MeV and we take the energy loss E_ion = 37 eV per ionisation (Glassgold 1995), yielding in the absence of dust around 18 000 ionisations before the positron finally annihilates with a bound electron to create a γ-ray. This γ-ray, in turn, has a very long stopping length and its effects on the mid-plane ionisation can be ignored. A more typical end to the positron trajectory is penetration in a solid particle. The stopping column density of positrons in silicon is Σ_β ≈ 3 kg m^-2 and the stopping length is χ_s ≈ 1.3 mm (see Sect. 3.1). We use the terminology , , and to denote the effective radii and masses that are available for absorbing and releasing positrons. These quantities will only be different from R₁, R₂, m₁ and m₂ if the particles are larger than χ_s. The effective radius can be written and the effective mass , where we subtract the inactive mass present deeper than χ_s in the solid.

We can calculate the positron production rate per kilogram of rocky material from the expression (2)Here \begin{lxirformule}$b=0.82$\end{lxirformule} is the branching ratio for positron emission, x₂₆ = 3 × 10^-10 is the nominal ratio of ²⁶Al nuclei to hydrogen molecules, τ is the decay time of ²⁶Al (1.03 Myr) and Z_sil is the dust mass fraction in rock (we use the approximate solar composition value of 0.005, following Lodders 2003). This gives r_β = 5 × 10⁵ kg^-1 s^-1.

In regions of the protoplanetary disc where the number density of dust is very high, the number densities of ions and electrons are given by a balance between ionisation, at rate ζ per neutral molecule, and absorption by dust grains (e.g. Okuzumi 2009), The index d introduced here implies summation over dust species. The ion and electron speeds, v_i and v_e, are discussed further in Sect. 2.3. The ionisation rate by positrons, ζ, is given by the expression (5)Here σ_nβ is the ionisation cross section of molecular hydrogen to positrons (6)where we have approximated the stopping column density in molecular hydrogen by that in silicon, Σ_β. The number density n_β of positrons in the gas, moving at the speed v_β, follows directly from the production equation (7)The first term represents the emission rate, with . The second term represents the absorption in solid particles and the third term the absorption in gas molecules. The equilibrium state with ṅ_β = 0 yields (8)We now rewrite the denominator in terms of the stopping column density Σ_β = ρ_•χ_s = ρ_gχ_s,g, where we introduced the material density ρ_• and the stopping length in the gas χ_s,g = (n_nσ_nβ)^-1(T_β/E_ion) and assumed that the stopping column is independent of the material (molecular hydrogen or silicate rock). That yields now, for solid particles with radii much smaller than χ_s, the approximate expression (9)Combining Eqs. (9) and (5) yields (10)with ζ₀ = 0.005r_βΣ_βσ_nβ ≈ 1.5 × 10^-19 s^-1 denoting the ionisation rate at the unsedimented value Z = 0.005. Here, Z = Z₁ + Z₂ is the total mass loading of solid particles in the gas. The ionisation rate thus rises from ζ₀ for unsedimented particles well-mixed with the gas in the entire vertical column to ζ_∞ = 200ζ₀ ≈ 3 × 10^-17 s^-1 when the particle mass loading is high, Z ≫ 1, e.g. in a dense mid-plane layer of sedimented pebbles.

2.3. Ion and electron currents

Knowing the ion and electric densities from Eqs. (3) and (4) allows us to calculate the neutralising current of ions and electrons as a function of the electric field strength E through The mean ion and electron velocities, ⟨ v_i ⟩ and ⟨ v_e ⟩ follow from the expression (Mori & Okuzumi 2016) (13)Here, α denotes either ions or electrons. The mean time between collisions is given for electrons as (14)and for ions as (15)The constant in the denominator expresses that the interaction cross section between ions and neutrals actually decreases with increasing ion speed, yielding a constant product of cross section and speed (Wannier 1953). The current densities in Eqs. (11) and (12) are proportional to E (through Eq. (13)) only if n_i and n_e are independent of E. In fact, as first pointed out by Okuzumi & Inutsuka (2015), n_i and n_e can decrease with increasing E when the electric field is sufficiently strong. If E is above a threshold (see Eqs. (19) and (20) for the thresholds for ion and electron heating, respectively), the electrons and ions obtain kinetic energies above the mean thermal energy of the neutrals through acceleration by the electric field. In the presence of such a strong electric field, the plasma particles collide with and adsorb onto dust grains more frequently as the field strength increases, and hence their number density decreases with increasing E. For the ion component, the added energy resides mainly in the mean motion because of the larger inertia of the ion molecules (e.g., HCO⁺) than that of neutrals (mainly H₂ and He). Hence, for electrically heated ions, the magnitude of the mean velocity ⟨ v_i ⟩ in Eq. (11) becomes equal to the mean speed ⟨ v_i ⟩ in Eq. (3). Using this fact together with Eqs. (3) and (11), one can show that the neutralising ion current approaches the constant expression (Okuzumi & Inutsuka 2015) (16)Inserting the ionisation rate from Eq. (10), we find the limiting values of the ion current, for low and high total particle mass loading, as We here assumed that the adsorption rate of ions is dominated by the dust component, i.e. that adsorption onto pebbles can be ignored, and scaled the result to dust grains of radius R₁ = 1 μm and material density ρ_• = 10³ kg m^-3. In Eq. (17) we scaled the ratio of the total mass loading of particles, Z = Z₁ + Z₂, to the mass loading in dust, Z₁, to 10. The ions experience significant heating, and obtain their plateau current density value, above the critical electric field strength of (Okuzumi & Inutsuka 2015) (19)This critical field strength is 1–2 orders of magnitude smaller than the breakdown strength (Eq. (1)) and has the same linear scaling with gas density. Hence the ions are well in their heated regime at the electric field strengths relevant for lightning discharge.

Electrons have larger inertia than neutrals and hence the electron kinetic energy is dominated by the random motion. The critical electric field for electron heating E_crit is (Okuzumi & Inutsuka 2015) (20)Here ℓ is the mean free path of the electrons to collide with neutrals. The limiting values of the electron current at strong electric field (21)Hence the electron current density plateaus at a value approximately 60 times higher lower than the ion current density. We can therefore neglect the electron current contribution to the neutralising current at the breakdown value of the electric field.

2.4. Sedimentation current

Charge separation can be driven by aerodynamical size sorting, if small and large solid particles are charged oppositely. Such opposite charging can e.g. be obtained by triboelectric charging in collisions between small dust and large pebbles. Desch & Cuzzi (2000) showed that triboelectric effects can charge pebbles up to approximately 10⁴ electron charges. The charge separation in particle components can then be transferred to a spatial charge separation by the differential motion of the particles through the gas. Charged pebbles sedimenting towards the mid-plane carry a current this way. The sedimentation current is (22)Here n₂ = Z₂ρ_g/m₂ is the pebble number density, v₂ is the terminal speed of the pebbles and q₂ is the pebble charge. The terminal speed at a height above the mid-plane of one scale-height, z = H, is given by (23)Here Ω denotes the Keplerian frequency at the distance r from the star and Ω²H is the gravitational acceleration towards the mid-plane. This expression simplifies to (24)That yields a sedimentation current density of (25)This is more than eight orders of magnitude lower than the neutralising ion current for Z₁ = 0.005 (Eq. (17)). The sedimentation current would increase at a higher pebble mass loading Z₂, but this is not compatible with the assumption that pebbles fall at the terminal speed at one scale-height above the mid-plane. Increasing the pebble charge or decreasing the pebble size and semi-major axis would also boost the current, but overall it seems hopeless to gain the eight orders of magnitude needed to fight the neutralising gas current.

2.5. Radial drift current

Particles also separate aerodynamically by radial drift caused by the radial pressure support of the gas (Weidenschilling 1977; Nakagawa et al. 1986). Tanaka et al. (2005) provide general solutions for the equilibrium velocities of gas and multiple particle components that couple via drag forces. The radial velocity of pebbles relative to the smaller dust grains embedded in the gas becomes (26)Here St₂ = R₂ρ_•/ (Hρ_g) is the pebble Stokes number and Δv is the sub-Keplerian speed parameter (Youdin & Johansen 2007). Tanaka et al. (2005) introduced the inverse Stokes number Γ = St^-1 and we derived Eq. (26) in the limit Γ₁ ≫ 1, Γ₂ ≫ 2, Γ₁ ≫ Γ₂ and Z₂ ≫ 1. The relative drift speed at high Z is dominated by the outwards motion of gas and dust pushed through the almost stationary pebble component, but given charge neutrality q₂n₂ = −(q₁n₁ + n_i−n_e) we can use q₂n₂ to describe the charge density of the current. Inserting typical parameters for the pebbles and the gas we obtain (27)The low mass loading limit (Z₂ ≪ 1) yields the same value and scaling, but without the dependence on Z₂. The drift current density for high mass loading is then (28)The maximum radial drift current is thus more than 20 times higher than the maximum sedimentation current (Eq. (25)) at Z₂ ≫ 1, but only becomes comparable to the neutralising gas current at the breakdown electric field (Eq. (18)) for Z₁ ~ 10⁶.

2.6. Full model

The derivations presented in the previous subsections were made under a set of simplifying assumptions, e.g. that the adsorption rate of electrons and ions on dust particles are not affected by the charge of the dust. We also rather handwavingly referred to the triboelectric charging model of Desch & Cuzzi (2000) when setting the characteristic charge of pebbles to q₂ = 10⁴e. In Appendices A and B we therefore present a fully numerical solution to the complex charge equation equilibrium. We include the release of secondary electrons from penetrating positrons and show that this leads to a large positive charging of the pebbles, at the expense of a negative charge on the gas. However, the results confirm that aerodynamical charge separation is much too slow to compete with the neutralising current. Therefore we consider in the next section a current driven by the positrons themselves.

3. Nuclear battery effect

We now describe a charging mechanism akin to a nuclear battery. The simplest form of a direct charging nuclear battery has a radioactive isotope on one plate of a capacitor. The nuclear decay products are absorbed by the opposite plate. If the emitted particles are charged (positrons, electrons, alpha particles or fission fragments), then they will charge up both plates and build up a voltage difference over the capacitor. In a protoplanetary disc, dense pebble structures emit positrons from the decay of ²⁶Al. This leads to a net negative charging of the pebbles – provided that the positrons stream into regions of significant gas that stops the positrons but does not emit any positrons in the opposite direction. Importantly, there will only be a net current if the positrons are stopped by gas, as the dust component is in itself in complete balance between emitting and absorbing positrons.

3.1. Characteristic current

The characteristic positron current density arising from a pebble region, dense enough to be optically thick to its own positrons, is (29)Here e is the positron charge, r_β is the emission rate of positrons per unit dust mass and unit time, ρ_p is the local mass density of solids and ℓ_β is the stopping length of the positrons. We divided the expression by a factor six to take into account the random direction of the positrons. We can express the stopping length in terms of the stopping column Σ_β = ρ_pℓ_β, yielding the positron current expression (30)Where the activity of ²⁶Al relatively to the value at the formation conditions of CAIs, f₂₆, is defined as (31)We used here Σ_β ≈ 3 kg/m² in silicon¹. This current density is orders of magnitude higher than the currents carried by pebble sedimentation (Eq. (25)) or radial drift (Eq. (28)) and is larger than the neutralising ion current for dust mass-loadings Z₁> 50 (Eq. (18)). Such mass loading of small dust in the gas is hard to achieve under nominal conditions in the protoplanetary disc. However, in Sect. 4 we apply the nuclear battery effect to dense filaments in this mid-plane layer and to circumplanetesimal pebble discs. Dust is produced continuously in collisions between pebbles and hence the dust density is expected to follow the density of larger pebbles. In the coagulation-fragmentation simulations of Birnstiel et al. (2011) small dust constitutes approximately 1–10% of the mass in pebbles of mm-cm sizes. We generally assume in this paper that this dust fraction is valid also in overdense pebble regions. The pebble-pebble collision timescale can be estimated as (32)The first lightning discharge will thus be delayed until small dust has been produced to reduce the conductivity. However, the dust production time-scale to produce small dust fragments is much shorter than the decay time of ²⁶Al, within the range of likely physical parameters in Eq. (32). We assume that the dust within dense pebble filaments undergoes continuous release in pebble-pebble collisions and recoagulation to maintain a high dust mass loading for millions of years.

As we discuss further in the next subsection, the neutralising gas current in dense pebble structures may be even lower than the expression in Eq. (18), if the gradient in the particle density is high. In that case positrons streaming out of the region are not replaced by positrons streaming in and hence the ionisation rate inside of the structure is reduced.

3.2. Optimal conditions for the nuclear battery

Consider a dense pebble region emitting positrons that are absorbed in the surrounding gas. The charging of the pebbles builds up an electric field that points from the surroundings into the negatively charged pebble region. Initially the positrons are stopped by collisions with gas molecules, but as the pebble region becomes increasingly negatively charged, the positrons lose more and more of their kinetic energy to work against the electric field. The width of the pebble structure is now key to obtaining neutralisation through lightning discharge. If the pebble structure has a too large width, W, then the positrons are trapped in the pebble region before the breakdown electric field value is reached. The critical width W^† is given through the relation (33)Here T_β = 0.66 MeV is the kinetic energy of the positrons, e is the positron charge and E_BD is the breakdown value of the electric field. Using the expression for the breakdown electric field from Eq. (1) we obtain (34)Here f_BD defines the ratio of the actual breakdown field to the nominal value given in Desch & Cuzzi (2000). A one-dimensional pebble filament should have a width comparable to W for efficient conversion of positron energy into lightning discharge; broader filaments do not reach the breakdown field strength at all and narrower filaments reach the breakdown field strength over such short length-scales that positrons only deposit a minor fraction of their energy in the electric field. Two-dimensional pebble sheets, on the other hand, build up an electric field that has a constant value up to a height comparable to the extent of the sheet itself. Hence a pebble sheet loads energy into the electric field efficiently independently of the width, for any width less than W^†. Reaching the breakdown electric field, ions and electrons accelerated by the electric field to neutralise the pebble region must now release an energy similar to the original kinetic energy of the positrons. This way the energy of the positron emission becomes converted to kinetic energy and heat in the lightning discharge current.

The density of the pebble structure does not affect the critical width found in Eq. (34). However, the highest positron current is obtained when the column density of the pebble structure is similar to the stopping column of the positrons. Structures of lower column density exhibit a lower current density than the characteristic value given in Eq. (30), while structures of higher column density self-absorb most of the emitted positrons and therefore have a low efficiency in converting the positron emission to charge separation. The highest positron current combined with the highest charging efficiency are therefore obtained when the pebble density obeys (35)That gives a pebble density of (36)Higher pebble densities yield inefficient charging and lower pebble densities result in a reduction of the positron current. Combining with the critical width for conversion of positron kinetic energy to electric potential, from Eq. (34), yields the optimal solids-to-gas ratio for a pebble filament and the minimum solids-to-gas ratio for a pebble sheet as (37)Lower solids-to-gas ratios can still yield efficient conversion of positron energy to heating, but the current emanating from pebble structures of lower densities will be lower than the characteristic current. The characteristic current is nevertheless quite high, at J_β = 3.6 × 10^-14 C m^-2 s^-1, so even a factor 10 reduction yields currents in excess of the neutralising current at Z₁ = 500.

Fig. 1 The charging rate by positron release in an axisymmetric pebble filament of half-width 20 km. The filament has a peak density of 1000 times the gas density and the background dust density is 0.1 times the gas density. The colored contours show the negative charging rate and the white contour lines show the dust-to-gas ratio. The charging rate has been obtained by a Monte Carlo calculation of the emission and absorption of 107 positrons. A hundred randomly picked positron trajectories are shown in green, released at the green dot and absorbed at the red dot (which happens mostly outside of the plotted region). The filament charges rapidly negative as most of the released positrons are absorbed in the surrounding gas.

3.3. Simulations of nuclear battery effect

The calculations presented in the previous subsection do not capture the details of the positron emission and absorption process. To get a better handle on realistic values of the positron current and the charging efficiency we perform a Monte Carlo simulation of positron emission and absorption in a dense particle filament. We use a Gaussian profile for the density of particles, with peak mass loading ϵ₂ relative to the gas, embedded in a constant background particle density with mass loading ϵ₁, (38)Note that ϵ₁ and ϵ₂ do not represent respectively small and large particles (and hence are unrelated to Z₁ and Z₂ employed in Sect. 2); from the view of positron emission the particle size is insignificant, as long as the particles are smaller than approximately 1 mm, so that the positron can freely escape from the interior. We create positrons at random positions that are weighted with the ratio of the local particle density ϵ₁ + ϵ₂ to the total dust mass. The positrons are given random 3D directions and are absorbed after interacting with the characteristic stopping column of Σ_β in gas and particles. In Fig. 1 we present the results of a simulation with ρ_g = 10^-7 kg m^-3, ϵ₁ = 0.1, ϵ₂ = 1000 and W = 20 km. The center of the filament charges rapidly negative, at 100 electron charges per cubic meter per second. The measured current density streaming out of the filament is shown in Fig. 2, as a function of the solids-to-gas ratio (a measure of the distance from the centre of the filament). The current density is above 3 × 10^-14 C m^-2 s^-1 in most of the region and peaks at a value of 4.5 × 10^-14 C m^-2 s^-1, near the characteristic value found in Eq. (30). The rate of change of the electric field is related to the current density through the simple expression , with vacuum permittivity ϵ₀. The right axis of Fig. 2 shows that the breakdown value of the electric field is reached within a characteristic time-scale of only 10³ s. The filament thus quickly builds up an electric field strength near the breakdown value and the continued charging by positron release is neutralised by lightning, keeping the electric field around the breakdown value.

Fig. 2 Current density of the positrons as a function of the solids-to-gas ratio in the particle filament shown in Fig. 1. The current density is above 2 × 10-14 C m-2 s-1 through the entire filament. The right axis shows the rate of change of the electric field. The breakdown electric field will be reached in the filament after a characteristic time-scale of only approximately 103 s.

The efficiency of the nuclear battery effect can be calculated from the net charging. We perform here a series of experiments where we cover a pebble structure and its surroundings by a fixed grid and monitor the emission and absorption of positrons in each cell. We consider a gas density of ρ_g = 10^-7 kg m^-3 and a box size of 10⁵ km to capture the positron stopping length in the gas of approximately 3 × 10⁴ km. The dust-to-gas ratio was chosen to be ϵ₁ = 0.1 and the pebble-to-gas ratio varied from ϵ₂ = 1 to ϵ₂ = 10³. The pebble structure scale is fixed at 3 × 10³ km, larger than the example shown in Fig. 1 because we must be careful to capture here the large length-scales of positron absorption in the gas. The transition between being thin and thick to positron emission happens then at a pebble mass loading of ϵ₂ ≈ 10.

Fig. 3 Charging efficiency (solid lines) and current density (dashed lines) as a function of the mass loading in pebbles, for both 1D filaments (blue lines) and 2D sheets (red lines), for a characteristic structure width of 3000 km. The charging efficiency is high at low pebble densities when the positrons can escape freely from the structure, but falls for high pebble densities due to self-absorption of positrons within the structure. The current behaves oppositely: at low pebble densities the current is low, since the structure is thin to its own positron emission, while the current density plateaus around the characteristic value of 3.6 × 10-14 C m-2 s-1 for high pebble densities.

Releasing a total of N_β positrons inside the pebble structure and measuring a net charge of Q_j on the grid, where j is the grid index, yields an efficiency of (39)Here we sum only over regions inside of the pebble structure, defined as ρ₂>ϵ₁ρ_g. This exclusion is necessary since positrons released near the boundary of the computational domain are not replaced with a counterstreaming flux and hence appear to have high charging efficiency. A uniform particle density would yield η_c ≈ 0, since all cells would have equally many emitted and absorbed positrons. The efficiency of the nuclear battery effect is shown in Fig. 3 as a function the pebble mass loading ϵ₂. We consider both 1D filaments and 2D sheets. The most efficient charging is obtained when the pebble density is low, as most positrons can escape such dilute pebble structures. The current is only 10^-14 C m^-2 s^-1 for low values of the pebble density, but plateaus around the characteristic value of 3.6 × 10^-14 C m^-2 s^-1 when the pebble structure transitions from being thin to thick to its own positron emission (at ϵ₂ ≈ 10).

4. Application to circumplanetesimal discs and protoplanetary discs

The nuclear battery effect identified and described in Sect. 3 provides a pathway to convert the decay energy of ²⁶Al to lightning discharge that can flash heat solids. The total energy present in positrons released by ²⁶Al – of decay rate r_β = 5 × 10⁵ kg^-1 s^-1, energy T_β = 0.66 MeV and decay constant τ₂₆ = 1.03 Myr – is (40)The heating and partial melting of chondrule precursors requires approximately 1.6 × 10⁶ J kg^-1 (Hevey & Sanders 2006). The positrons thus release, over the life-time of ²⁶Al of a few million years, approximately the energy needed to melt all rock in the protoplanetary disc². Transferring the positron energy to the electric field and then to lightning discharge can therefore lead to flash heating of a significant fraction of the solids in the disc, if the energy transfer and heating efficiencies are both high. We discuss here two settings where the nuclear battery effect could operate to thermally process solids: in pebble discs orbiting young planetesimals and in the sedimented mid-plane layer of pebbles in the main protoplanetary disc.

4.1. Nuclear battery in a circumplanetesimal pebble disc

The high abundance of sodium within chondrules can be understood if chondrule precursors were heated in an environment that was several orders of magnitude denser than the Roche density at the distance of the asteroid belt from the Sun (Alexander et al. 2008). Such dense regions should be undergoing gravitational collapse. The gravitational collapse phase is short, on the order of a few orbital periods (Nesvorný et al. 2010; Wahlberg Jansson & Johansen 2014). How can chondrules in chondrites reflect formation conditions that were prevalent for only the last tiny fraction of the free life of a pebble in a protoplanetary disc before it was incorporated in a planetesimal? Actually, densities higher the Roche limit do not necessarily imply gravitational collapse, if the collapsing region is rotationally supported. Johansen & Lacerda (2010) demonstrated how pebbles accreted onto protoplanets enter a prograde particle disc in stable orbit around the central protoplanet. The Roche density in such a particle disc will be orders of magnitude higher than the Roche density that results from the tidal force of the central star alone. Could the high particle densities required for chondrule heating have been sustained for hundreds of thousands of years in thin pebble discs orbiting around planetesimals and protoplanets?

Fig. 4 Conditions for the nuclear battery effect in circumplanetesimal discs. The top left plot shows the pebble column density as a function of the distance from the planetesimal, for two values of the planetesimal radius (Rp = 100 km and Rp = 1000 km) and three values of the random pebble speed c. The vertical dot-dashed lines indicate the Hill radius at 2.5 AU and 25 AU, respectively, and the horizontal dashed line indicates the stopping column of positrons. The pebble column density is below the stopping column of positrons in the outer region of the circumplanetesimal disc, unless the random pebble motion is very high, so positron emission leads to an efficient charging of the circumplanetesimal disc. The top right plot shows the integrated disc mass relative to the mass of the central planetesimal. The small planetesimal has a relatively more massive disc than the larger protoplanet and hence can convert a higher fraction of the central mass to chondrules. The bottom left plot shows the scale-height of the pebble disc. The stopping length of positrons at the breakdown value of the electric field is indicated for the nominal breakdown current at 2.5 AU. The bottom right plot shows the mid-plane density of pebbles compared to the minimum pebble density required for lack of isotopic fractionation and stability of the liquid phase by Cuzzi & Alexander (2006) and the density range implied by the observed high abundance of Na within chondrules by Alexander et al. (2008). The pebble density in the circumplanetesimal disc environment lies perfectly within these experimental constraints.

A circumplanetesimal particle disc will not evolve in isolation from its surroundings. The sub-Keplerian gas acts as a friction on the outermost pebbles, unbinding them from the gravity of the planetesimal so that they mix with other pebbles in the protoplanetary disc and participate in new generations of planetesimal formation and pebble accretion. Planetesimals that enter eccentric orbits, e.g. perturbed by the growing protoplanets (Johansen et al. 2015), experience even stronger gas headwinds that could erode off pebbles even more efficiently. A large fraction of the chondrules that are flash-heated in a circumplanetesimal disc are nevertheless expected to become accreted onto the central planetesimal. Radiometric dating of individual chondrules, by the long-lived Pb-Pb chronometer, shows that chondrules in the same chondrite display a wide range of ages (Connelly et al. 2012). This age difference could reflect either that the parent bodies of those meteorites accumulated late from a mixture of chondrules that formed in other circumplanetesimal discs or that the chondrules that are accreted through a circumplanetesimal disc experience efficient mixing between chondrules formed at various epochs.

We can estimate the column density of a circumplanetesimal disc from the assumption that the disc evolves towards a state that is marginally unstable to self-gravity. In this view the circumplanetesimal disc accretes mass onto the central planetesimal by large-scale gravitational stresses when onfalling pebbles lower the Safronov-Toomre Q parameter to below unity (Safronov 1960; Toomre 1964). The marginal case of Q = 1 gives a column density of (41)Here c is the particle velocity dispersion and Ω_pla(r) is the Keplerian frequency profile of the circumplanetesimal disc. The velocity dispersion c is generally not known. Completely elastic particle collisions cause an exponential transfer of the relative shearing motion of the particles to random particle motion. The evolution of the random motion stops when the coefficient of restitution drops at higher collision speeds (Goldreich & Tremaine 1978). Estimates of the velocity dispersion in Saturn’s rings motivates setting this transition in the elasticity of the collisions at 0.001 m s^-1. However, the fluffy silicate pebbles that form in the inner regions of protoplanetary discs may have very different energy dissipation in inelastic collisions than the ice particles of Saturn’s rings and hence we keep c as a free parameter in the circumplanetesimal disc model. The profiles of particle column density, particle scale-height and mid-plane particle density, all calculated for Q = 1, come out as In these expressions we have normalised the radial distance from the planetesimal, r, by 100 times the planetesimal radius, R_pla. We plot the radial profiles of pebble column density, ratio of disc mass to planetesimal mass, pebble scale-height and pebble mid-plane density for a planetesimal of radius 100 km and a protoplanet of radius 1000 km in Fig. 4. The conditions for the nuclear battery depend strongly on the adopted value of the random pebble speed c as well as on the planetesimal size. Discs with random pebble speed in the approximate range between 0.01 m s^-1 to 0.1 m s^-1 have column densities close to the stopping column of positrons between 100 and 1000 planetesimals radii where most of the disc mass resides.

The ratio of the circumplanetesimal disc mass to the mass of the central planetesimal determines the maximum production efficiency of chondrules. The top right panel of Fig. 4 shows that the protoplanet of radius 1000 km has only a small mass fraction in the circumplanetesimal disc and hence can only produce a relatively small amount of chondrules relative to its mass. This may nevertheless not necessarily be in conflict with the high mass fraction of chondrules in chondrites, if these chondrules were produced in pebble discs around larger protoplanets. As the outer edge of the disc expands by angular momentum exchange with the inner disc, the chondrules formed in the circumplanetesimal disc are recycled into the protoplanetary disc and could contribute significantly to the mass reservoir from which the smaller chondrite parent bodies formed, even if the production efficiency relative to the protoplanet mass is low.

Importantly, the bottom right panel of Fig. 4 illustrates how the pebble density inside of the Hill radius of the planetesimal reaches several orders of magnitude above the Roche density in the protoplanetary disc. Thus the mid-plane conditions around both the planetesimal and the protoplanet are right in the middle of the constraints inferred by Alexander et al. (2008) and Alexander & Ebel (2012). The circumplanetesimal discs that formed in the simulations by Johansen & Lacerda (2010) had sizes around 10% of the Hill radius, or 100 times the planetesimal radius. These discs would nevertheless not form around small planetesimals, due to the headwind of the gas. Therefore, it is doubtful whether small planetesimals can host extensive circumplanetesimal discs, unless the gravitational collapse in itself leads to the formation of a compact, rotationally supported pebble clump that subsequently expands by viscous stresses into its Hill radius. Such a dense disc expanding outwards could have enough density to withstand the friction from the sub-Keplerian gas. This model is nevertheless, very speculative and this illustrates well the need for computer simulations to understand the formation and dynamics of circumplanetesimal discs better.

Fig. 5 Conditions for the nuclear battery effect in the protoplanetary disc mid-plane. The top left plot shows the pebble column density for three values of gas depletion fg and the stopping column of positrons. The mid-plane layer transitions from optically thick to optically thin to positron stopping at a gas mass of 10% of the nominal value in the asteroid belt. The top right plot shows the breakdown value of the electric field. The bottom left plot shows the scale-height of the mid-plane layer (dotted line) together with the stopping length of the positrons at the breakdown value of the electric field. The mid-plane layer is clearly much thicker than the optimal thickness for transferring the positron energy to the electric field, but substructures of 1% of the mid-plane thickness are more in agreement with the optimal width (dot-dashed line). The bottom right plot shows the mid-plane density of gas, the Roche density and the conditions for chondrule formation. The lack of isotope fractionation and stabilisation of melted silicate can be achieved at solids-to-gas ratios of 10 (at 1 AU) to 100 (at 2.5 AU), but the stricter density conditions imposed by the observed high abundance of Na (Alexander et al. 2008) can only be achieved interior of 1 AU where the Roche density is high.

Pebbles in the circumplanetesimal disc will spend most of their time in the outer regions of the disc where the accretion time-scale is longest. The viscous time-scale of a circumplanetesimal particle disc of optical depth τ (=Σ_pσ_p/m_p) is (see e.g. Charnoz et al. 2009) (45)The viscous time-scale in the limit τ ≫ 1 can now be expressed as (46)This is clearly a very long time-scale for a protoplanet of radius R_p = 10³ km, which validates our argument that accretion will happen by gravitational stresses and thus will be regulated by the continuous infall of new pebble material onto the circumplanetesimal disc. A planetesimal of radius R_p = 100 km has a much shorter accretion time-scale, unless the random pebble motion is less than 0.1 m s^-1.

Johansen et al. (2015) showed that planetesimals can grow to protoplanets of Mars mass, initially by planetesimal-planetesimal collisions and later mainly by pebble accretion, in the asteroid belt after a few million years. The pebble accretion process would likely endow these protoplanets with extensive circumplanetesimal pebble discs. Drazkowska et al. (2016) demonstrated that the pile up of particles in the terrestrial planet region can trigger planetesimal formation already after a few hundred thousand years evolution of the protoplanetary disc. Therefore protoplanets may have populated the terrestrial planet region even earlier than the asteroid belt region and hence chondrules could have been produced around these protoplanets and subsequently transported to the asteroid formation region to become incorporated into the asteroids and protoplanets that formed there.

4.2. Nuclear battery in the protoplanetary disc

The circumplanetesimal disc model presented in the previous subsection is appealing because the mid-plane pebble density reaches orders of magnitude higher than the Roche density, as required for heating chondrule precursors to the melting temperature and avoiding isotopic fractionation and maintaining the high abundance of Na. This model nevertheless implies that the chondrules found in chondrites were first processed in the environment near another planetesimal – and hence we must evoke a generation of planetesimals and protoplanets that formed prior to the asteroids in the asteroid belt and produced extensive amounts of chondrules that became spread throughout the protoplanetary disc.

Here, we explore an alternative setting for lightning discharge and thermal processing of solids, namely the sedimented mid-plane layer in the protoplanetary disc. This mid-plane layer can not support particle densities in excess of the Roche density for long time-scales and thus the conditions are at odds with the formation of chondrules in the asteroid belt (we show below that the required densities can be reached inside of 1 AU). The mid-plane layer is nevertheless prone to charging by positron emission and the photons from the resulting lightning discharge could flash-heat and sublimate material in its path. These silicate vapours could in turn recondense as tiny matrix particles. Some authors have actually suggested that mineral growth in chondrules implies free propagation of the cooling radiation into a region of relatively low pebble density (Miura & Yamamoto 2014). Such conditions would be in agreement with the mid-plane layer of the protoplanetary disc as a chondrule formation site.

At the nominal ionisation efficiency of ²⁶Al, we require a mass loading in small dust, Z₁, of at least 50 in the mid-plane layer, to suppress the neutralising current. The pebble mid-plane layer may nevertheless experience reduced positron ionisation, since the positrons streaming out from the mid-plane layer are not replaced by positrons streaming in from below and above the mid-plane layer, since those regions are underdense in pebbles. Hence the conductivity will be reduced significantly compared to thick structures. Thus even a dust-to-gas ratio as low as 10 in the mid-plane layer may yield low enough conductivities to drive charging of the mid-plane layer faster than the neutralisation.

We show in Fig. 5 the conditions for the nuclear battery effect in the sedimented mid-plane layer of pebbles in the protoplanetary disc. We consider a protoplanetary disc model with gas density ρ_g = f_g × 10^-7 kg m^-3 at 2.5 AU, where f_g measures the depletion of gas due to accretion, and a column density profile that falls with distance from the star r as Σ_g ∝ r^-1. The mid-plane layer scale-height is set to H_p = 0.01 H, a typical value when stirred by the streaming instability (Carrera et al. 2015). We consider three values of f_g – 1, 0.1 and 0.01 – to parameterise the time evolution of the protoplanetary disc.

The column density of pebbles is thick to positron stopping in the interior of the disc and thin in the exterior. The transition point is at approximately 2.5 AU for a gas density that is 10% of the nominal value of 10^-7 kg m^-3 at 2.5 AU. A mid-plane layer thin to positron stopping is ideal for the nuclear battery effect, as the positrons from the entire mid-plane layer, irrespective of any substructure, will be absorbed by the gas outside of the mid-plane. The mid-plane layer is nevertheless much wider, by two orders of magnitude at 2.5 AU, than the optimal width for loading positron energy into the electric field. Therefore substructure of width 1% of the mid-plane layer width must be present to reach the breakdown value of the electric field before the positrons are stopped by the electric field. This implies that the mid-plane first charges up by sending positrons to the gas above and below the mid-plane layer, but that the subsequent trapping of positrons in the mid-plane layer charges up smaller and smaller structures until the breakdown value of the electric field is reached. Simulations of the streaming instability show that filaments form on many scales within the mid-plane layer (Johansen & Youdin 2007; Bai & Stone 2010; Johansen et al. 2015). It remains to be seen whether this density spectrum is prone to charging to the breakdown value of the electric field; we plan to present such calculations in a future publication.

The Roche density in the protoplanetary disc is higher inside of 3 AU than the constraint proposed by Cuzzi & Alexander (2006) to suppress isotopic fractionation in the heating and stabilise the molten chondrules against sublimation. The stricter conditions proposed by Alexander et al. (2008) to main the high abundance of Na within the chondrules can only be achieved below the Roche density well inside of 1 AU in the disc.

5. Heating by lightning discharge in a dense pebble region

We consider in this section the heating effect of lightning discharge into a dense region of pebbles and dust grains. The solid particles act like a thick, insulating blanket that facilitates high temperatures and long cooling times for the relatively moderate discharge powers expected in protoplanetary disc conditions. The discharge is assumed to have a total length h into a cylindrical pebble filament of width R_∗. The cylindrical symmetry simplifies the radiative transfer equations significantly. Radiative transfer within a disc-like geometry is more complex, but we expect that heating by lightning penetration into a circumplanetesimal disc will be similar to the cylindrical geometry, as long as the region heated by lightning maintains approximately cylindrical shape.

5.1. Radiative equilibrium

We develop the analytical and numerical heating model around a discrete radiative-transfer approach where the cylinder is divided into annuli of width ΔR corresponding to an optical depth interval of unity. The optical depth with respect to radiation is dominated by the small dust component. The mean free path of photons, with wavelengths comparable to or smaller than the dust so that we can consider geometric absorption, is (47)The outer edge of the ith cylindrical annulus is denoted R_i = iΔR. We assume that each cylindrical ring radiates like a black body whose radiation is absorbed entirely in its two neighbouring annuli. That gives the temperature evolution equation for ring i as (48)Here ρ_p = ρ₁ + ρ₂ is the total density of solid particles, c_p is the specific heat capacity of the particles, V_i is the volume of annulus i, σ_SB is the Stefan-Boltzmann constant and T_i is the temperature of annulus i. We then introduce the source function . The energy flux over the edge of annulus i is (49)In equilibrium this equals the energy flux release by the lightning, ℱ_i = IEh. Here I is the current carried by the lightning and E the breakdown electric field. The continuous equation for the energy flux is (50)We define the optical depth τ measured inwards from the edge of the pebble region, with a total optical depth of τ_⋆. The optical depth thus relates to the radial coordinate R through R = ΔR(τ_⋆−τ). Rewriting Eq. (50) in terms of optical depth τ we get (51)Rearranging we find (52)Integration yields (53)The constant C is fixed by requiring that all the released energy is radiated outwards at the τ = 0 surface, corresponding to the requirement τ_⋆S(1/2) = IE/ (2πΔR) in the outermost annulus with τ = 1/2. This gives the final source function expression (54)The temperature at optical depth τ = τ_⋆−1/2, the approximate optical depth of the first annulus, in the limit τ_⋆ ≫ 1, is (55)We used here a typical current carried by terrestrial lightning and landed fortuitously in the ballpark temperature for chondrule heating. The characteristic density of protoplanetary discs is many orders of magnitude lower than the density of the terrestrial atmosphere. The discharge process may nevertheless be taken to follow the so-called Townsend scalings (Ebert et al. 2010), which yield a discharge length-scale proportional to the number density of neutrals n, a discharge time-scale proportional to 1 /n and a discharge current independent of n. Hood & Kring (1996) used the 100 ms duration of terrestrial lightning to argue that protoplanetary disc lightning has a characteristic time-scale of 10⁴–10⁵ s. However, terrestrial lightning strokes consist of several flashes of 30μs duration, which would yield 10² s time-scales in the asteroid belt region of the hydrogen-dominated solar protoplanetary discs with ρ_g = 10^-7 kg m^-3.

Note that the central temperature in Eq. (55) has a very weak scaling with τ_⋆ – even a dramatic increase in the optical depth will not lead to a significant temperature increase in the core of the filament. This is a result of the radial thinning of the radiation at large distances from the lightning. The temperature at the optical at the optical surface (τ = 1/2) relates to the deep temperature T_⋆ through (56)That gives T(1/2) = 366 K for T_⋆ = 1754 K and τ_⋆ = 100. The total energy content of the pebble region at the equilibrium temperature is obtained by integrating over the energy contents of the cylindrical annuli, (57)Inserting Eq. (54) we obtain the expression (58)We find empirically that the term in hard brackets within the integral has a mean value, weighted by τ, of approximately unity for a wide range of values of τ_⋆. Hence the integral is well approximated by . That gives the energy content per length (59)Here we assumed that the density of solid particles is dominated by the pebble component, ρ_p ≈ ρ₂. The time to heat the dense pebble region is (60)

5.2. Heating and cooling rates

Chondrule precursors appear to have been heated to near-liquidus temperatures, reaching the range 1500–2200 K (Hewins & Radomsky 1990). This constraint is mainly based on what appears to be relict, unmelted grains inside chondrules. If these grains are instead interpreted as interlopers that penetrated a fully molten spherule (Connolly & Hewins 1995), then chondrule precursors could have been heated to much higher temperatures. The cooling rates of chondrules can be constrained by studying crystal growth in furnace experiments of chondrule analogues or from the retention of volatiles such as Na in chondrules. Hewins & Radomsky (1990) found that cooling rates in the range 10²–10³ K h^-1 are consistent with experimental data. The upper range of these cooling rates corresponds to what is inferred for matrix grains (Scott & Krot 2005). On the other hand, Yurimoto & Wasson (2002) and Miura & Yamamoto (2014) inferred much higher rates, in the range 10⁵–10⁷ K h^-1. The upper range of those rates are comparable to cooling rate of isolated molten droplets, (61)Much lower cooling rates can be achieved in an optically thick region of chondrules. We measure the heating and cooling rates of pebbles exposed to radiation from lightning discharge by solving Eq. (48) numerically starting at a background temperature of T = 200 K. Absorption and release of latent heat ℒ by melting and crystallisation is included using a simplified approach where we translate the internal energy e to temperature T through the relations Here T_s is the solidus temperature (taken to be 1500 K), T_l is the liquidus temperature (taken to be 2000 K), and the internal energies at solidus and liquidus, respectively, are defined as Absorption and release of latent heat does not change the equilibrium temperature profile, as given in Eq. (54), but delays both the initial heating and the subsequent cooling.

Fig. 6 Temperature as a function of distance from the centre of the lightning discharge channel for a current I = 105 A, gas density ρg = 3 × 10-7 kg m-3, dust density ρ1 = 10-3 kg m-3, pebble density ρ2 = 10-2 kg m-3 and optical depth τ⋆ = 300. The discharge channel has a width of 100 m, indicated by the vertical dotted line. The discharge is assumed to last 2 × 103 s, followed by a cooling phase that is evolved up to a time of 104 s. Heating is rapid and uniform within the discharge channel, while the outer regions only heat after radiation diffuses out from the discharge region. The equilibrium temperature follows the analytical solution outside of the discharge channel.

The width of the lightning discharge channel must be scaled properly to the protoplanetary disc conditions. Horányi & Robertson (1996) argued that the discharge channel has a width of 1000–10 000 times the mean free path of the electrons. This scaling relation yields (67)Here f_dis denotes the discharge width in units of the mean free path of electrons. We use here f_dis = 10³ and note that the equilibrium temperature solution of Eq. (54) is still valid outside the discharge channel, while the temperature within the discharge channel is lower than the equilibrium solution as the decreased heat flux necessitates a lower temperature gradient.

In Fig. 6 we show the evolution of the temperature profile of the pebble filament for selected times. We assume that the lightning discharge heats for a time-scale of t_heat = 2 × 10³ s. The temperature climbs steadily during this time, but only approaches the analytical equilibrium solution (Eq. (54)) outside the discharge channel. The temperature is relatively constant within the discharge channel itself and significantly below the analytical solution due to the lower energy flux. Cooling occurs relatively isothermally over the pebble filament. In Fig. 7 we show the heating and cooling rate at selected regions in the pebble filament as a function of time. Heating is initially rapid, but slows down as the equilibrium temperature is approached. The cooling rate lies at a few times 10³ K h^-1 in the first few hundred seconds after heating terminates and subsequently falls to a few times 10² K h^-1 as the temperature drops to 1000 K. These cooling times lie towards the upper end of the values inferred from chondrule textures (Hewins & Radomsky 1990). The cooling time of the filament should broadly follow the heating time, so following Eq. (60) lower cooling rates could be achieved in filaments of large optical depth or lower gas density.

Fig. 7 Heating rate (red) and cooling rate (blue) of the pebble filament heated by the lightning discharge shown in Fig. 6. The curves show, from top to bottom, the heating and cooling rate at optical depth τ = 0, 50, 100, 150, 200. The cooing rates lie at a few times 103 K h-1 in the first 100 seconds after heating terminates at theat = 2 × 103 s and subsequently fall to below 103 K h-1 over the next 1000 s. These cooling rates are in good agreement with the upper range of cooling rates inferred from chondrule textures (Hewins & Radomsky 1990).

6. Summary

We have demonstrated in this paper that the decay energy of the short-lived radionuclide ²⁶Al can be harvested to drive lightning discharge in protoplanetary discs. Positrons released from dense pebble structures load their energy into the electric field that builds up by the negative charging of these structures. This process is akin to a direct-charging nuclear battery where a radioactive isotope charges a capacitor. The total energy available in the kinetic energy of the released positrons is comparable to the energy required to melt all the solid material in a protoplanetary disc. Such an immense energy released in lightning discharge may be observable in nearby protoplanetary discs with e.g. the ALMA telescopes (Muranushi et al. 2015).

The neutralising gas current, driven with an efficiency (conductivity) that is primarily set by positron ionisation of hydrogen molecules, can no longer compete with the positron current when the dust-to-gas ratio of small grains is higher than a few tens. This leaves lightning discharge as the natural path to charge neutralisation, releasing the energy stored in the electric field in thin channels. Eisenhour et al. (1994) demonstrated that heating by infrared laser light yields melting primarily of pebbles larger than 0.1 mm, due to the inability of smaller pebble aggregates in absorbing infrared photons. In this picture, the largest chondrules represent the largest available chondrule precursors pebbles, while the lack of small chondrules in most chondrites reflects inefficient heating of small pebbles by the heat front propagating from the lightning channel.

We find that the optimal conditions for heating pebbles by lightning discharge arise in dense particle discs that orbit around planetesimals and protoplanets. The high dust content of such discs implies a low conductivity and the rapid transition from dense pebble environment to the surrounding protoplanetary disc dominated by gas is ideal for driving an outwards directed current of positrons. Circumplanetesimal discs are appealing from a cosmochemical point of view also, since the prevalent conditions there match the high dust densities required for chondrule melting and high abundance of volatiles (Cuzzi & Alexander 2006; Alexander et al. 2008). The formation and dynamics of such circumplanetesimal discs is nevertheless not well understood and should be the focus of future research.

The sedimented mid-plane layer of pebbles in the protoplanetary disc can also become charged by the emission of positrons, if the mass loading in small dust grains is high enough in the mid-plane (larger than 10–50) to reduce the neutralising current. The densities of solids in the mid-plane layer are too low in the region of the asteroid belt to allow heating to the liquid phase and the solids are expected to sublimate completely upon flash heating. The subsequent cooling and recondensation could form small particles similar to the matrix grains found between chondrules in chondrite meteorites, although the long cooling times inferred for some matrix grains may imply that matrix and chondrules formed under elevated densities (Scott & Krot 2005). The density in the mid-plane is much higher in sub-AU distances than in the asteroid belt, so the conditions close to the star are more in agreement with the cosmochemical constraints on the elevated density of chondrule precursors.

The solar protoplanetary disc appears to have been injected with a nominal value of ²⁶Al in its birth cluster. Stellar systems are born with a range of abundances from ten times below to ten times above the solar abundance (Vasileiadis et al. 2013; Parker et al. 2014; Lichtenberg et al. 2016). Our results show that the energy budget for chondrule melting could have been supplied by the decay of ²⁶Al. The pebbles in the rare systems with ten times the solar abundance of ²⁶Al must have undergone widespread and repeated melting that could even affect the water contents, and hence the habitability, of terrestrial planets, if ices are continuously driven off the pebbles in large regions of the protoplanetary disc. Therefore, it is imperative to understand the heating mechanism that processes solids in protoplanetary discs to form chondrules and matrix. We have demonstrated here that the decay of ²⁶Al – already generally accepted to be the energy source for melting and differentiating asteroids – may also have provided the energy needed for the widespread thermal processing and melting of dust and pebbles in protoplanetary discs.

Acknowledgments

A.J. thanks the Knut and Alice Wallenberg Foundation (grants 2012.0150, 2014.0017, 2014.0048), the Swedish Research Council (grant 2014-5775) and the European Research Council (ERC Starting Grant 278675-PEBBLE2PLANET and ERC Consolidator Grant 724687-PLANETESYS) for their financial support. S.O. is supported by the JSPS Grants-in-Aid for Scientific Research (Grant Number JP17K18812). The authors are thankful to the referee, Steven Desch, for a thorough referee report that led to many improvements in the presentation of the paper.

References

Appendix A: Charge equations

In this appendix we formulate the equations for production of electrons and ions by high-energy positrons and the charging of dust particles. In the following appendix, Appendix B, we will find the equilibrium state for electron and ion densities and the charge of dust particles and calculate the conductivity of the gas and the neutralising current as a function of the electric field and compare this to the current driven by the radial drift of the pebbles.

Appendix A.1: Ionisation rates

The evolution equations for the electron and ion densities are The first term appearing in both Eqs. (A.1) and (A.2) represents gas ionisation by the passing positron (here σ_nβ and v_β are the cross sections of ionisation and the speed of the positrons). The second term represents dust absorption of electrons and ions. The index d implies summation over species. The adsorption rates of electrons and ions onto the dust grains (absorption cross section of dust times speed of the charged species) are denoted K_de and K_di, respectively. These rates are complicated functions of the charge of the dust; we derive expressions for K_de and K_di in Appendix C.

The third term of Eq. (A.1) represents the emission of N_sd secondary electrons as a positron passes the surface of a particles. In order to obtain the values of N_sd, averaged over impact parameter, we have performed a Monte Carlo simulation of the penetration of a positron through a solid particle of finite radius. The results are discussed in Appendix D. We find N_s1 = 0.04 for micrometer-sized particles and N_s2 = 0.08 for millimeter-sized pebbles. Secondary electrons leave silicon targets with typical energies of approximately 10 eV (Dapor et al. 2008); this energy is matched by the electric potential energy from a positive charge of +10⁷e on a millimeter-sized pebble. Therefore we can safely assume for all pebble charges found in Appendix B that secondary electrons can freely leave the surface of the pebbles.

The third term of Eq. (A.2) represents the absorption of positrons by a gas molecule as the positron comes to rest, to create a γ-ray when it combines with an electron. The rate of this event equals the ionisation rate by the positron times the ratio of the ionisation energy of the hydrogen molecule, taken to be E_ion = 37 eV, to the kinetic energy of the positron (T_β = 0.66 MeV for the decay of ²⁶Al).

Appendix A.2: Dust charging

The charging of small dust grains and macroscopic pebbles follows the evolution equations The first two terms represent collisions with electrons and ions, the third term the production of secondary electrons, the fourth term the charging by release of positrons (r_β is the positron production rate per unit mass) and the fifth term is the absorption of positrons (χ_s ≈ 1.3 mm is the stopping length of the positrons in the solid particles).

The last term present in both Eqs. (A.3) and (A.4) represents the transfer of charge from small dust grains to pebbles by collisions (Okuzumi 2009). We adopt here the closure that the number density of dust and pebbles remains constant, under the assumption that small dust grains transfer all their charge to a pebble following a bouncing collision. The adsorption rates K_de, K_di and K₁₂ are described in Appendix C.

Fig. A.1 The ionisation rate ζ versus the dust-to-gas ratio. The initial rise of the ionisation rate is due to the increased density of 26Al, while the plateau occurs when the rise is balanced by an increased absorption rate of positrons in dust grains. For the production of secondary electrons, which occurs in collisions with dust, we have converted to an effective rate on the gas particles for comparison. Secondary electrons come to dominate after Z ≈ 1000.

Appendix A.3: Positron density and charge equilibrium

The positron number density follows the equation (A.5)The first term is the emission rate, with . The second term is the absorption in solid particles and the third term the absorption in gas molecules. All these terms are balanced in corresponding terms in the evolution of the number densities of electrons and ions and the charging of dust. We see that the charge is conserved, since (A.6)The integral of the equation is Q, the free charge density, which we set to be zero. This means that any of the four dynamical equations can be written as a sum over the other three, so that the number of independent dynamical equations is only three.

The ionisation rate by positrons is equivalent to (A.7)We plot the effective ionisation rate ζ of positrons in Fig. A.1 as a function of the mass loading of dust Z = ρ_p/ρ_g, calculated from the equilibrium state of Eq. (A.5). The ionisation rate increases initially with Z, as the concentration of ²⁶Al increases relative to the gas. However, above Z = 1 the ionisation rate flattens out. At those high particle densities the increase in positron emission rate is balanced by the stopping of the positrons in the solid particles. We also show in Fig. A.1 the production rate of secondary electrons, calculated as an effective ionisation rate of the gas. The contribution of secondary electrons to the electron density dominates at mass loading above Z ~ 1000.

Appendix A.4: Electron and ion heating

The speed of the electrons, of mass m_e, and the ions, of mass m_i, can be obtained from their temperatures T_e and T_i through Okuzumi & Inutsuka (2015) showed that the electrons and ions obtain energies above the mean thermal energy of the molecules in the presence of strong electric fields. This is due to transfer of energy from the mean motion along the electric field to random motion after collisions with neutrals. The increased random motion in term decreases both the number densities and mean current speed of electrons and ions. Following Mori & Okuzumi (2016) we write the electron temperature as (A.10)The critical electric field for electron heating E_crit is (A.11)Here m_e and m_n are the electron and neutral masses, respectively. The mean free path of the electrons is simply ℓ_e = (n_nσ_en)^-1 with σ_en = 10^-19 m².

The ions are heated at higher values of the electric field. Mori & Okuzumi (2016) give an expression for the ion temperature, (A.12)The ion temperature increases as the square of the electric field, while the electron temperature only increases linearly with the electric field in the limit of high field strength.

Appendix B: Charge equilibrium

The direct temporal integration to equilibrium of the charging equations presented in Appendix A is made very expensive by the large range of time-scales present in the problem. Particularly at high particle mass loading and strong electric fields, relevant for the heating of chondrules, the time-scale for collisions between electrons and dust particles is very low. Therefore we search for the equilibrium condition by solving directly for . Unfortunately, none of these dynamical equations are analytically solvable, due to the complicated expressions for the adsorption rates K_de and K_di (see Appendix C).

Fig. B.1 The electron and ion densities (top left; blue lines indicate electrons and red lines ions), the pebble charge (top right, blue lines indicate negative charge, red lines positive charge, black lines the charge relative to a model that neglects charge transfer in dust collisions), the electron and ion currents (bottom left) and the total current (bottom right). The electron and ion critical fields are marked in blue and red dotted lines, while the break down electric field is shown with black dotted lines (we show both the nominal value and a value ten times lower relevant when taking into account the distribution of electron energies). The positive pebble charging is mainly due to the emission of secondary electrons as positrons penetrate dust surfaces. The equilibrium charge increases with increasing field strength, since the field accelerates electrons and ions to high speeds, thus increasing the collision rate with the charged pebbles.

Appendix B.1: Numerical solution method

One of the four variables in the solution vector (n_e,n_i,q₁,q₂) can be eliminated by using the charge conservation condition from Eq. (A.6). Our algorithm for finding the equilibrium state consists of stepping through n_i from 10^-30n_n to 10^-5n_n. For each value of n_i we find n_e that satisfies ṅ_e = 0. This requires additional knowledge of q₁ and q₂. The value of q₁, given n_e and n_i, is calculated from and q₂ is calculated straightforwardly from the charge conservation equation. We then minimise the value of ṅ_i along the curve ṅ_e = 0. All zero points are found by bisection to a tolerance level of 10^-20. We implemented the algorithm in Fortran 90 and compiled with quadruple precision in order to avoid underflows in the very large terms that are present in the dynamical equations.

We consider electric field strengths E between 10^-4 V m^-1 and 10³ V m^-1 and a density of small particles relative to the gas density Z₁ between 10^-3 and 10⁴. The density of pebbles is fixed at a 10 times higher value, Z₂ = 10Z₁. The radius of small dust is set to 1 μm and the radius of pebbles to 1 mm; both components are assumed spherical with an internal density of ρ_• = 10³ kg m^-3. For each value of E and Z₁ we solve the charge equilibrium equations to obtain n_e, n_i, q₁ and q₂.

Appendix B.2: Results

The equilibrium state is plotted in Fig. B.1. We show the electron and ion densities, the charge on the macroscopic pebbles and the current densities of electrons and ions, all as functions of the electric field and the density of small dust grains. We also indicate the critical electric fields for electron heating (blue dotted line) and ion heating (red dotted line) as well as the break down electric field to trigger lightning discharge (black dashed lines), for the nominal value of Desch & Cuzzi (2000) as well as for a 10 times lower value that may reflect the actual breakdown field strength (Eq. (1)). Both ion and electron densities fall after reaching their respective critical electric field strengths, due to the increase in collision rates with dust. The ion density falls faster than the electron density as the ions interact weakly with neutrals and hence are accelerated to very high speeds that result in high collision rates with the dust particles. The electron density falls much more slowly when the mass loading reaches Z₁ ~ 100 (corresponding to a total metallicity of Z ~ 10³, see Fig. A.1). This is due to the release of secondary electrons as positrons pass through lots of dust grains at those densities.

Appendix B.3: Limiting values for electron and ion currents

Okuzumi & Inutsuka (2015) showed that the heating of ions causes the ion current density to plateau out at (B.1)Here ζ = 1.7 × 10^-17 s^-1 is the ionisation rate by positrons from ²⁶Al, valid when Z ≫ 1 (see Fig. A.1). The limiting ion current density is in good agreement with Fig. B.1. For the electrons we obtain the limiting expression (B.2)This is also in excellent agreement with Fig. B.1 until Z ≈ 100. Beyond this mass loading the electron production rate rises due to release of secondary electrons. The electron current is nevertheless less important than the ion current in setting the overall current density at all relevant values of the mass loading. Importantly, the limiting values of the ion and electron currents are independent of the gas density and hence the results presented in Fig. B.1 are relatively independent of the distance from the star, except for a small effect of the disc temperature on the conductivities.

Appendix B.4: Charge flux

Figure B.1 shows that the charge on large pebbles, q₂, is positive for high mass loading. This is due to the emission of secondary electrons as a positron passes through a particle surface. Pebbles are large enough to completely stop positrons. Therefore a positron will on the average pass through λ₂/λ₁ small dust grains before it is absorbed by a pebble, with λ₁ and λ₂ indicating the mean free path of the positron relative to dust and pebbles. The relative production rate of secondary electrons is (B.3)Here A₁/A₂ is the total surface area of dust grains relative to the total surface area of pebbles within a given volume. The total adsorption rate of secondary electrons by the two dust components is also proportional to A_i, Therefore, secondary electron release only yields net charging if N_s1 is not equal to N_s2. In Appendix D we show that N_s1 = 0.04 and N_s2 = 0.08. This asymmetry in the secondary electron production yields net positive charging of the pebbles at the expense of the dust grains. The resulting pebble charge can be very high, up to 10⁵ electron charges for strong electric fields. The charging of pebbles takes much higher values when the charge flux from collisions with small dust is ignored (as indicated by the black contour lines in the upper right plot of Fig. B.1).

In Fig. B.2 we show the charge flux for the six components of the charging model (including the neutrals), normalised to the charge flux to the dust and pebble component by positron ionisation, for Z₁ = 100 and E = 10 V m^-1. Charging by decay of ²⁶Al (DE) and stopping of positrons (ST) are in close balance for both dust and pebbles, with just a small fraction of positrons stopped by the gas (AN). The production rate of secondary electrons (SE) by the dust is less than 100 times (i.e., the ratio of total dust surface to total pebble surface) that of the pebbles, due to the lower efficiency of secondary electron production of the dust. This asymmetry of secondary electron production is indeed the source of pebble charging, as the pebbles must charge positively to pull in additional electrons that can compensate for the loss of secondary electrons. Charge transfer in dust collisions (CO) is another source of pebble flux on the pebbles; as is evident in Fig. B.1 this dust flux reduces the pebble charge by approximately a factor ten at (E,Z₁) = (10 V m^-1,100).

Fig. B.2 Diagram of the charge fluxes between the five components neutrals, electrons, ions, dust, pebbles and positrons for E = 10 V m-1 and Z1 = 100. The fluxes are normalised by the charge flux to dust and pebbles by ionisation. The abbreviations are: IO = ionisation, AD = adsorption, SE = secondary electrons, CO = dust collisions, DE = nuclear decay and positron emission, ST = stopping of positrons in solids, AN = annihilation of positrons in gas.

Fig. B.3 Neutralising gas current J (black contours) and drift current Jdri (red for positive current and blue for negative), as a function of electric field E and dust metallicity Z1. The total current of the solids is separating to the left of the thick yellow zero-line and neutralising to the right of the zero-line. The drift current is at best 2–3 orders of magnitude smaller than the neutralising gas current.

Appendix B.5: Drift current

The large, positive charging of the pebbles implies that spatially separating pebbles from gas and dust will lead to the build-up of an electric field. The question is whether the charge separation is continuously neutralised by the electron/ion current or whether it builds up rapidly enough to necessitate neutralisation by lightning discharge. We compare in Fig. B.3 the drift current to the gas current. We calculate the drift current from (B.4)where v_1x, v_2x and u_x are the pebble, dust and gas velocities in the radial direction. These velocities result from both radial drift, as the pressure gradient on the gas forces large pebbles to drift outwards at the cost of small dust and gas drifting inwards, as well as the electric acceleration of the dust and pebbles. The equilibrium drift velocities are found by finding the dynamical equilibrium of the gas velocity (u_x,u_y) and the particle velocities (v_1x,v_1y), (v_2x,v_2y) in the local frame corotating with the protoplanetary disc at the Keplerian frequency Ω. This equilibrium is given by the condition Here τ₁ and τ₂ are the friction times of dust and pebbles, respectively, and Δv denotes the sub-Keplerian speed resulting from the radial pressure gradient (Nakagawa et al. 1986). In Fig. B.3 we overplot contours of the drift current on the neutralising gas current. The drift current is separating (opposite sign as the electric field) to the left of the zero-line neutralising (same sign as the electric field) to the right of the zero-line. The zero-line of the drift current is reached well before the breakdown value of the electric field at high dust mass loadings. This implies that the mobility of the solids leads to neutralisation of the charge separation well before lightning discharge can occur.

Fig. B.4 Neutralising gas current and drift current for a model with reduced 26Al abundance (10% of nominal value) and decreased gas density (ρg = 10-8 kg m-3). Both the neutralising gas current and the breakdown value of the electric field are a factor 10 lower than in Fig. B.3, but the drift current is still at least a factor 10 too small to compete with the neutralising current.

In Fig. B.4 we present a model where the abundance of ²⁶Al has been reduced by a factor 10 relative to the nominal value. This leads to a corresponding decrease in the gas conductivity (see Eq. (B.1)). We also reduced the gas density by a factor 10 to mimic conditions at the late evolution stages of a protoplanetary disc after a few Myr of gas accretion onto the star. The pebble charge is actually not reduced by the decrease of the positron production, since the equilibrium is set by the relative efficiency of secondary electron production between small dust grains and macroscopic pebbles. The drift current is also relatively unaffected by the change of gas density, since the reduction in the number density of the solids is balanced by their higher mobility relative to the gas. The lowered gas density decreases both the critical electric field strengths for electron and ion heating as well as the breakdown electric field. The results nevertheless show that the zero-drift current remains relatively unchanged, although the lowered breakdown value of the electric field has approached the zero-drift current by almost a factor 10. A further reduction in the gas conductivity and the breakdown electric field by an order of magnitude could open a pathway to separating charges faster than neutralisation. However, such low gas densities would imply particle densities much lower than those required for chondrule formation.

In conclusion, it is very hard to obtain charge-separating radial drift currents that are higher than the neutralising current. An additional concern about chondrule formation by lightning in overdense filaments is the time-scale for which these high densities can be sustained. Simulations of planetesimal formation by the streaming instability demonstrate that the particle densities can reach very high values of several thousand times the gas density during the collapse (Johansen et al. 2015). This phase is nevertheless short, on the order of a few orbital periods, and hence can not lead to widespread thermal processing of the particles involved in the collapse. Therefore we consider in Sect. 3 a potentially much power powerful source of currents in protoplanetary discs, namely the flux of positrons emanating from dense pebble regions.

Appendix C: Adsorption rates

The cross section σ_dj for collisions between a charged species j and a dust grain of type d with charge q_d, radius R_d and surface potential φ_d = (4πϵ₀)^-1q_d/R_d is given by (see e.g. Shukla & Mamun 2002) (C.1)Here m_j and v_j are the mass and speed of the charged species and the ± symbol refers to ions of charge +e and electrons of charge −e, respectively. The index j can now be used to sum over a velocity distribution. We denote the electron temperature by T_e and the ion temperature by T_i. These temperatures can obtain values different from the thermal temperature T due to heating by the electric field, see discussion in Appendix A.4. For the electrons we assume that the mean speed is negligible compared to the thermal speed (Okuzumi & Inutsuka 2015; Mori & Okuzumi 2016). The absorption rates of electrons are therefore given by expressions that are independent of the mean speed, The adsorption rate of ions must be calculated by considering the separate flux contributions from the ion temperature T_i and the ion speed along the electric field, ⟨ v_i ⟩. The adsorption rate is calculated through the integral (C.4)We then insert a Maxwellian distribution function that includes a mean velocity ⟨ v_i ⟩, (C.5)Here v_i = | v_i | is the magnitude of v_i and θ is the angle between v_i and ⟨ v_i ⟩. We denote by v₀ the collision speed below which the effective collision cross section is zero. That gives For the repulsion case with φ_d> 0 we obtain (C.8)For the attraction case with φ_d< 0 we obtain an expression already given in Shukla & Mamun (2002), (C.9)Note that the repulsion and accretion expressions for K_di are the same for φ_d = v₀ = 0.

Appendix C.1: Dust adsorption

The adsorption of small dust grains by large pebbles is given by the expression (C.10)Here v₁₂ is the relative drift speed between dust grains and pebbles, calculated from the drag force equilibrium of Eqs. (B.5)–(B.10), and E₁ is the kinetic energy (C.11)

Appendix D: Secondary electrons

Positrons released by ²⁶Al in regions of very high particle density are primarily stopped by other particles. The passage of a positron through the surface of a solid particle is accompanied by the emission of N_s secondary electrons. This leads to a net negative charging of dust particles, since a small fraction of the positrons are absorbed by the gas. The secondary electrons also contribute to the conductivity of the gas. In this appendix we calculate the yield of secondary electrons as a function of the size of the particle that the positron passes through. The secondary electron yield is an important parameter in the charge equilibrium model presented in the main paper.

Fig. D.1 Secondary electron production versus the incident angle μ = cos(θ) for particles of radii μm, mm and cm. The secondary electron production rises with decreasing μ and experiences a jump at the angle when the positron penetrates all the way through the solid. Finally, the secondary electron production falls when the pathway through the particle is short at very low μ (dotted line). The integrated yield (dashed lines) is about 0.08 for mm and cm solids and 0.04 for μm solids.

The production of secondary electrons by a positron experiencing energy loss dE_β in the interval between χ and χ + dχ along its penetration path through a solid particle can be written according to Lin & Joy (2005) as (D.1)Here ϵ is the ionisation potential of electrons inside the material, α is the inverse stopping length of the released electrons and χ_⊥ is the shortest distance to the surface of the particle. The particle sizes we consider are much larger than the stopping length of the electrons and hence the value of α obtained from experiments with flat surfaces can be applied to curved particles as well. The perpendicular distance χ_⊥ can be written as function of χ for a positron with angle of incidence θ relative to the normal of the surface as (D.2)We use the parameters ϵ = 90 eV and λ = 1 /α = 2.7 nm relevant for pure silicon (Lin & Joy 2005). The energy loss rate is calculated from the Bethe-Bloch equation modified for positrons (e.g. Olive et al. 2014, chap. 31), (D.3)Here the constant K = , with N_A denoting Avogadro’s number, r_e the classical electron radius, m_e the electron mass and c the speed of light. The material is characterised by the density ρ_i, nuclear charge Z_i and atomic number A_i. The positron energy is defined through β = v/c and . The function F is specific to positrons and is defined as (D.4)We calculate the number of secondary electrons produced as a function of μ = cos(θ) for particles of μm, mm and cm sizes. Secondary electrons are mainly produced within a few times the mean free path λ from the surface. We therefore only integrate Eq. (D.1) at distances within 10λ from the surface, in order to speed up the integration. We use an energy of E₀ = 0.66 MeV for the positron. The results are shown in Fig. D.1. The yield of secondary electrons increases initially with decreasing μ, as the positron moves closer to the surface. The yield experiences a sudden increase at an angle where the positron passes all the way

through the particle. Finally, the yield falls along a characteristic curve (D.5)for very low values of μ. Here 2Rμ is the passage length through the particle and χ_s is the total stopping length in the material (1 mm for silicate particles). The integrated yield gives the total number of secondary electrons through (D.6)This gives an integrated yield of N_sec = 0.08 for particles of mm and cm sizes and N_sec = 0.04 for μm-sized particles.

All Figures

Fig. 1 The charging rate by positron release in an axisymmetric pebble filament of half-width 20 km. The filament has a peak density of 1000 times the gas density and the background dust density is 0.1 times the gas density. The colored contours show the negative charging rate and the white contour lines show the dust-to-gas ratio. The charging rate has been obtained by a Monte Carlo calculation of the emission and absorption of 107 positrons. A hundred randomly picked positron trajectories are shown in green, released at the green dot and absorbed at the red dot (which happens mostly outside of the plotted region). The filament charges rapidly negative as most of the released positrons are absorbed in the surrounding gas.

In the text

	Fig. 2 Current density of the positrons as a function of the solids-to-gas ratio in the particle filament shown in Fig. 1. The current density is above 2 × 10-14 C m-2 s-1 through the entire filament. The right axis shows the rate of change of the electric field. The breakdown electric field will be reached in the filament after a characteristic time-scale of only approximately 103 s.
In the text

Fig. 3 Charging efficiency (solid lines) and current density (dashed lines) as a function of the mass loading in pebbles, for both 1D filaments (blue lines) and 2D sheets (red lines), for a characteristic structure width of 3000 km. The charging efficiency is high at low pebble densities when the positrons can escape freely from the structure, but falls for high pebble densities due to self-absorption of positrons within the structure. The current behaves oppositely: at low pebble densities the current is low, since the structure is thin to its own positron emission, while the current density plateaus around the characteristic value of 3.6 × 10-14 C m-2 s-1 for high pebble densities.

In the text

Fig. 4 Conditions for the nuclear battery effect in circumplanetesimal discs. The top left plot shows the pebble column density as a function of the distance from the planetesimal, for two values of the planetesimal radius (Rp = 100 km and Rp = 1000 km) and three values of the random pebble speed c. The vertical dot-dashed lines indicate the Hill radius at 2.5 AU and 25 AU, respectively, and the horizontal dashed line indicates the stopping column of positrons. The pebble column density is below the stopping column of positrons in the outer region of the circumplanetesimal disc, unless the random pebble motion is very high, so positron emission leads to an efficient charging of the circumplanetesimal disc. The top right plot shows the integrated disc mass relative to the mass of the central planetesimal. The small planetesimal has a relatively more massive disc than the larger protoplanet and hence can convert a higher fraction of the central mass to chondrules. The bottom left plot shows the scale-height of the pebble disc. The stopping length of positrons at the breakdown value of the electric field is indicated for the nominal breakdown current at 2.5 AU. The bottom right plot shows the mid-plane density of pebbles compared to the minimum pebble density required for lack of isotopic fractionation and stability of the liquid phase by Cuzzi & Alexander (2006) and the density range implied by the observed high abundance of Na within chondrules by Alexander et al. (2008). The pebble density in the circumplanetesimal disc environment lies perfectly within these experimental constraints.

In the text

Fig. 5 Conditions for the nuclear battery effect in the protoplanetary disc mid-plane. The top left plot shows the pebble column density for three values of gas depletion fg and the stopping column of positrons. The mid-plane layer transitions from optically thick to optically thin to positron stopping at a gas mass of 10% of the nominal value in the asteroid belt. The top right plot shows the breakdown value of the electric field. The bottom left plot shows the scale-height of the mid-plane layer (dotted line) together with the stopping length of the positrons at the breakdown value of the electric field. The mid-plane layer is clearly much thicker than the optimal thickness for transferring the positron energy to the electric field, but substructures of 1% of the mid-plane thickness are more in agreement with the optimal width (dot-dashed line). The bottom right plot shows the mid-plane density of gas, the Roche density and the conditions for chondrule formation. The lack of isotope fractionation and stabilisation of melted silicate can be achieved at solids-to-gas ratios of 10 (at 1 AU) to 100 (at 2.5 AU), but the stricter density conditions imposed by the observed high abundance of Na (Alexander et al. 2008) can only be achieved interior of 1 AU where the Roche density is high.

In the text

Fig. 6 Temperature as a function of distance from the centre of the lightning discharge channel for a current I = 105 A, gas density ρg = 3 × 10-7 kg m-3, dust density ρ1 = 10-3 kg m-3, pebble density ρ2 = 10-2 kg m-3 and optical depth τ⋆ = 300. The discharge channel has a width of 100 m, indicated by the vertical dotted line. The discharge is assumed to last 2 × 103 s, followed by a cooling phase that is evolved up to a time of 104 s. Heating is rapid and uniform within the discharge channel, while the outer regions only heat after radiation diffuses out from the discharge region. The equilibrium temperature follows the analytical solution outside of the discharge channel.

In the text

Fig. 7 Heating rate (red) and cooling rate (blue) of the pebble filament heated by the lightning discharge shown in Fig. 6. The curves show, from top to bottom, the heating and cooling rate at optical depth τ = 0, 50, 100, 150, 200. The cooing rates lie at a few times 103 K h-1 in the first 100 seconds after heating terminates at theat = 2 × 103 s and subsequently fall to below 103 K h-1 over the next 1000 s. These cooling rates are in good agreement with the upper range of cooling rates inferred from chondrule textures (Hewins & Radomsky 1990).

In the text

Fig. A.1 The ionisation rate ζ versus the dust-to-gas ratio. The initial rise of the ionisation rate is due to the increased density of 26Al, while the plateau occurs when the rise is balanced by an increased absorption rate of positrons in dust grains. For the production of secondary electrons, which occurs in collisions with dust, we have converted to an effective rate on the gas particles for comparison. Secondary electrons come to dominate after Z ≈ 1000.

In the text

Fig. B.1 The electron and ion densities (top left; blue lines indicate electrons and red lines ions), the pebble charge (top right, blue lines indicate negative charge, red lines positive charge, black lines the charge relative to a model that neglects charge transfer in dust collisions), the electron and ion currents (bottom left) and the total current (bottom right). The electron and ion critical fields are marked in blue and red dotted lines, while the break down electric field is shown with black dotted lines (we show both the nominal value and a value ten times lower relevant when taking into account the distribution of electron energies). The positive pebble charging is mainly due to the emission of secondary electrons as positrons penetrate dust surfaces. The equilibrium charge increases with increasing field strength, since the field accelerates electrons and ions to high speeds, thus increasing the collision rate with the charged pebbles.

In the text

Fig. B.2 Diagram of the charge fluxes between the five components neutrals, electrons, ions, dust, pebbles and positrons for E = 10 V m-1 and Z1 = 100. The fluxes are normalised by the charge flux to dust and pebbles by ionisation. The abbreviations are: IO = ionisation, AD = adsorption, SE = secondary electrons, CO = dust collisions, DE = nuclear decay and positron emission, ST = stopping of positrons in solids, AN = annihilation of positrons in gas.

In the text

Fig. B.3 Neutralising gas current J (black contours) and drift current Jdri (red for positive current and blue for negative), as a function of electric field E and dust metallicity Z1. The total current of the solids is separating to the left of the thick yellow zero-line and neutralising to the right of the zero-line. The drift current is at best 2–3 orders of magnitude smaller than the neutralising gas current.

In the text

	Fig. B.4 Neutralising gas current and drift current for a model with reduced 26Al abundance (10% of nominal value) and decreased gas density (ρg = 10-8 kg m-3). Both the neutralising gas current and the breakdown value of the electric field are a factor 10 lower than in Fig. B.3, but the drift current is still at least a factor 10 too small to compete with the neutralising current.
In the text

Fig. D.1 Secondary electron production versus the incident angle μ = cos(θ) for particles of radii μm, mm and cm. The secondary electron production rises with decreasing μ and experiences a jump at the angle when the positron penetrates all the way through the solid. Finally, the secondary electron production falls when the pathway through the particle is short at very low μ (dotted line). The integrated yield (dashed lines) is about 0.08 for mm and cm solids and 0.04 for μm solids.

In the text