Dust-Dust Collisional Charging and Lightning in Protoplanetary Discs

We study the role of dust-dust collisional charging in protoplanetary discs. We show that dust-dust collisional charging becomes an important process in determining the charge state of dust and gas, if there is dust enhancement and/or dust is fluffy, so that dust surface area per disc volume is locally increased. We solve the charge equilibrium equations for various disc environments and dust number density $\eta$, using general purpose graphic processors (GPGPU) and {\sc cuda} programming language. We found that as dust number density $\eta$ increases, the charge distribution experience four phases. In one of these phases the electrostatic field $E$ caused by dust motion increases as $E \propto \eta^4$. As a result, macroscopic electric discharge takes place, for example at $\eta = 70$ (in units of minimum-mass solar nebula (MMSN) values, considering two groups of fluffy dust with radii $10^{-2}\unit{cm}$, $10^{2}\unit{cm}$). We present a model that describes the charge exchange processes in the discs as an electric circuit. We derive analytical formulae of critical dust number density for lightning, as functions of dust parameters. We estimate the total energy, intensity and event ratio of such discharges (`lightning'). We discuss the possibility of observing lightning and sprite discharges in protoplanetary discs by Astronomically Low Frequency ({\em ALF}) waves, {\em IR} images, {\em UV} lines, and high energy gamma rays. We also discuss the effects of lightning on chondrule heating, planetesimal growth and magnetorotational instability of the disc.


INTRODUCTION
Planets are formed in protoplanetary discs from interstellar dust. The electric charge state of the dust aggregates in the protoplanetary discs is one of the key parameters in understanding a number of aspects of protoplanetary discs and protoplanetary formation.
The inner structure of the dust aggregates, relative velocity, and electric charge are key parameters that determine the growth and migration of dust aggregates. Dust relative velocity (Brauer et al. 2008) includes random motion caused by turbulence (Ormel & Cuzzi 2007) and Brownian motion (Blum et al. 1996), and bulk motion caused by vertical sedimentation (Dullemond & Dominik 2004) and radial migration (Weidenschilling 1977). The collision velocity governs the growth rate (Suyama et al. 2008), compactification (Weidling et al. 2009), and disruption (Wada et al. 2008a), of the dust. Okuzumi (2009) considered the charge state of the dust aggregates in protoplanetary discs. They assumed that the charge state is determined by absorption equilibrium of ions and free electrons. Since electrons have much larger thermal velocity compared to positive ions, plasma absorption makes all dust to charge weakly negative. The repulsive Coulomb force may suppress dust-dust collisional growth for all but the heaviest dust species who can overcome the Coulomb barrier.
It is also possible that the charge state of the dust is affected by dust-dust collision. The effect has been simply ignored in most research due to the fact that in protoplanetary discs, dust has low number density, and is surrounded by weakly ionised plasma. We give quantitative estimates of the effect of dust-dust collision as a function of dust size, fractal dimension, and number density, and show for the first time that the dust-dust collision is actually an important factor in high dust number density regions of protoplanetary discs.
One of the possible dust-dust collisional charging mechanisms is known as the triboelectric process, where two bodies exchange electrons and sometimes molecular ions when they come into contact (e.g. Sickafoose et al. 2001). Another mechanism is possible for materials with spontaneous surface charge, such as H2O ice crystals (e.g. Kudin & Car 2008). In this mechanism the surface matter within typical depth ∼ 1.0 × 10 −4 cm (Mason & Dash 2000) is exchanged together with contained charge.
Surface space charge due to electron spill-out is widely known among metals and semiconductors (Somorjai 1994), the charge separation being ∼ 10 −7 cm deep for metals and ∼ 10 −5 cm deep for semiconductors. H2O is unique in that molecular ions OH − and H3O + holds the charge, and that proton exchange between the molecules (Grotthus mechanism; c.f. Agmon 1995) allows charge diffusion much faster than molecular ion diffusion. Thus surface charge separation develops as deep as ∼ 2.0 × 10 −4 cm (Dash et al. 2001). For example, ammonia lacks the mechanism (Goncalves et al. 1999). It is important that charge separation layer is deeper than exchange depth, because if the entire charge separation layer is exchanged, charge transport is neutral and collisional charging do not take place.
The dust-dust collisional charging due to the exchange of this spontaneous surface charge of ice crystals, is an established model in the context of meteorology (e.g. Takahashi 1978;Baker et al. 1987;Dash et al. 2001) that explains lightning on earth. When two ice dust of different surface states collide, they exchange their surface charge, producing charged dust. When the charged particles within nonconducting gas are separated by some external force, electric field grows between them. At the point the electric field is larger than the dielectric field strength of the gas, rapid ionisation of the gas occurs, converting the electrostatic energy into kinetic energy of the electrons and ions. This is electric discharge. Lightning in the earth's atmosphere is one of the most prominent, and well studied examples of electric discharge phenomena; in thunderclouds, typically 3.0 × 10 10 esu, or 1.0 × 10 1 C of electric charge is repeatedly separated and neutralized with typical length scales 1.0 × 10 5 cm (Koshak & Krider 1989).
In protoplanetary discs, lightning is one of the candidate mechanisms for chondrule heating, although compared to other models e.g. heating by shock wave (e.g. Miura et al. 2008), some difficulties have been pointed out (Weidnschilling 1997). For example, electric field cannot grow large enough to cause electrostatic breakdown in standard discs (Gibbard et al. 1997). Moreover, when mm-sized silicate aggregates made of µm-sized monomers are subject to electric discharge, they generally fragment without being thermally processed (Güttler et al. 2008).
Lightning in protoplanetary discs is strongly related to turbulence. The relative random velocity between the charged dust species that sets the dust to collide, results from the turbulence. Also the difference of the bulk velocities between the charged dust species that leads to macroscopic charge separation results from the turbulence.
The turbulent state of the accretion discs is often expressed in terms of viscous α parameter introduced by Shakura & Sunyaev (1973). Since the specific angular momentum increases outward in Keplerian discs, they satisfy Rayleigh's hydrodynamical stability criterion, and there are no clear mechanism for hydrodynamic turbulence in protoplanetary discs (Sano et al. 2004). On the other hand, the angular velocity decreases outward in Keplerian discs, they satisfy criterion for magnetorotational instability (MRI).
Therefore, if a protoplanetary disc is ionised enough to sustain magnetic field, MHD turbulence is excited and α parameter can be as large as 1.0 × 10 −3 ∼ 1.0 × 10 −1 (Sano et al. 1998). If the ionisation is suppressed, on the other hand, α ≃ 1.0 × 10 −5 . For a typical protoplanetary disc it is believed that so-called 'dead zones' form between 0.1 AU and 10 AU where instabilities are damped and gas flow is almost laminar (e.g. Gammie 1996). But it is possible that MRI is active in the whole disc, if sufficient ionisation degree is maintained, for example by turbulent mixing (Turner et al. 2007) or by self-sustained ionisation (Inutsuka & Sano 2005). Thus ionisation state of the protoplanetary discs is critical in determining α and understanding the fate of planetesimals and protoplanetary discs (e.g. Kretke & Lin 2007;Brauer et al. 2008).
The purpose of this paper is twofold: One is to solve the local charge exchange equilibrium of gas and dust numerically, for various dust parameters such as radii, fractal dimensions and dust number density, with dust-dust collisional charging taken into consideration; Given the results, the other goal is to determine the critical dust number density ηcrit under which lightning to take place, as analytical functions of other dust parameters such as radii, fractal dimensions and disc environment parameters such as temperature and gas number density. This paper is organized as follows. We define the terms we use in Table 1, and we list the symbols we frequently use in Table 2. In §2 we introduce the dynamic charge exchange equations and its equilibrium solution in schematic forms. We introduce circuit diagram to depict them (Fig. 1). In §3 we examine the processes in protoplanetary discs that set the parameters for the charge equilibrium equations. Crucial parameters are dust number density, the amount of charge exchange in single dust-dust collision, and relative velocity. In §4 we estimate the electrostatic field strength, and define the critical number density ηcrit for lightning in the protoplanetary discs. At this point all the equations are specified, and we solve them numerically. In §5 we show the results of the simulations. We describe four distinct phases of the charge distribution and explain the results using circuit diagrams. We also give analytical estimates for electric field strength in protoplanetary discs and critical number density η for lightning to occur. In §6 we discuss the possibility of various phenomena caused by the highly charged dust and lightning in protoplanetary discs, and their observations. particle gas dust plasma neutral gas ion electron smaller dust larger dust | | | | × | neutral cation anion cationic anionic Table 1. The terminology we use in this paper. 'Particle' is generic term for all components in the protoplanetary discs. Solid components are 'dust,' and the others are 'gas.' 'Gas' components are further subdivided into 'neutral gas,' and charged components, or 'plasma.' Finally, 'plasma' consists of 'electron,' the negative charge carrier, and various molecular 'ion,' the positive charge carrier. On the right side of the table, 'dust' is classified by their size as 'smaller dust' and 'larger dust.' Either can be 'anionic' or 'cationic' dust, depending on the material they consist of. We also use the one-letter symbols 'g', 'e', 'i', 'S', and 'L' for neutral gas, electron, ion, Smaller and Larger dust. The symbols for 'Cationic' and 'Anionic' dust are 'C' and 'A'. We use variable I to represent one of these symbols.

MODEL DESCRIPTION
In this section we describe our models. In 2.1 we model the disc and the dust at the unperturbed state, then introduce the models for dust number density. In 2.2 we model the charge density and charge separation processes.

Disc Model
Unless otherwise mentioned, we focus on a local, uniform box at certain orbital radius r near the equatorial plane of the protoplanetary disc. We model the protoplanetary disc based on the minimum-mass solar nebula (MMSN) model (Hayashi 1981). The gas surface density Σ ⊲⊳ g (r), disc scale height h ⊲⊳ (r) , and the temperature T ⊲⊳ (r) of the disc are where r is the distance from the central star. This leads to gas density distribution The dust-to-gas ratio in MMSN is approximately 1.0×10 −2 . We use the model by Cuzzi & Zahnle (2004), and introduce two species of dust, the smaller dust and the larger dust (see Table 1.) We further assume that surface density of the larger dust is 10 per cent of the total dust surface density. These two species are also either 'cationic' and 'anionic.' The 'cationic' species receives the positive electric charge through dust-dust collision. See Appendix A for the justification of this two-dust model. We can also represent the role of various molecular ions by one abstract ion species 'i,' according to Okuzumi (2009).
The motivation for this two-dust model is twofold. First, the two dust model is the simplest model that can handle the dust-dust collisional charge separation and the macroscopic relative velocity between the dust species. Second, the charge tendency of the dust and their size are strongly correlated. In one scenario, older dust are larger and also anionic. In another scenario, dust made of ice is larger and also cationic compared to dust made of silicate. (see §3.4 for the details.) Therefore, we expect that instead of considering four (cationic smaller dust, cationic larger dust, anionic smaller dust, and anionic larger dust) species of dust, we can correlate the two size species with the two charge tendency species, (Table 1), although both correspondences (smaller dust is cationic / larger dust is cationic) are possible.
To summarise, we define the reference density of the smaller dust ρ ⊲⊳ S (r) and the density of the larger dust ρ ⊲⊳ L (r) as We further assume that within a local condensation region, density for each component of the disc are multiplied. Figure 1. The circuit diagram of the charge exchange process in dust plasma. Each arrow represents 'current' density J, which has the unit esu cm −3 s −1 , the amount of charge passed from one component to the other per unit disc volume per unit time. The arrow points from the component that receives negative charge to the component that receives positive charge. In this figure, i and e are ions and electrons created from ionising neutral disc gas. C and A are cationic and anionic dust defined in §A.
J ei represents gas ionisation as 'current' from e vertex to i vertex; J iA , J iC , J Ae , and J Ce are ion and electron absorption to dust; J AC is dust-dust collisional charge separation and J (n) CA is neutralization current of charged dust-dust absorption.
Alternatively, we can think of protoplanetary discs with different gas or dust density than MMSN. We denote the ratio of the density of gas, smaller dust, and larger dust by ηg, η S , η L , respectively. Then the density of gas, smaller dust and larger dust is given by Mass of the smaller dust and the larger dust are m S and m L , respectively. The number density is density divided by dust mass: We estimate the mass as a function of the dust radius and the fractal dimension in §3.1.

Charge exchange equations
There are four species of charge carrier in our modelions, electrons, cationic, and anionic dust (Table 1). Charge exchange processes between these species are ionisation, plasma absorption, and dust-dust collision. The ionisation of the neutral gas molecules generates the ions and the electrons. Plasma absorption decreases the number of plasma particles and passes the lost charge to the dust aggregates. The dust aggregates also get charged by dust-dust collision. We label the particle species with letter I. The charge density carried by species I is QI (the unit is esu cm −3 ), and the charge transferred from species I to species I ′ is J I,I ′ (the unit is esu cm −3 s −1 ).
The charge density QI of a species I is the product of their number density nI and their average charge per particle qI. For dust species, we assume that nI is known from num-ber density model while qI is unknown; for ion and electrons we know qI but do not know nI. This constitutes the four dynamical equations for four unknown variables q A , q C , ni, ne : The current terms J I,I ′ are where we have included neutral gas ionisation Je,i, dustplasma absorption JA,i, JC,i, JA,e, JC,e, dust-dust collisional charge-up JA,C , and dust-dust collisional neutralization J (n) C,A terms. Here, vi and ve are the thermal velocity of the ions and the electrons, ng is the number density of the neutral gas, ζ is the ionisation rate, which is dominated by cosmic ray ionisation near equatorial, r = 2.7 AU of MMSN (Umebayashi & Nakano 2009). The exact value for these terms are given in §3. We have neglected, for example, the gas-phase recombination.
We want to solve the equilibrium equations for the dynamic equations (12-15): together with charge neutrality equation: We use circuit diagram ( Fig. 1) to depict the dynamical equations (12-15), and to interpret the numerical equilibrium solutions (23-27) in §5. The circuit diagram represents charge-exchange processes; each vertex represents the species of charge reservoir and each arrow represents the charge exchange process. The size of the vertex circles represents the amount of charge QI. The thickness of the arrows represents the amount of charge transfer J I,I ′ . We define the direction of the arrows so that the arrows point to the positive charge receivers.
In the system of equations depicted by a circuit diagram, charge density of each vertex QI corresponds to an unknown quantities. Therefore, the number of unknown quantities is equal to the number of vertexes NV . On the other hand, at the equilibrium, sum of the current flowing into each vertex is required to be zero (Kirchhoff's Laws); this gives us NV equations but only NV − 1 of them are independent. Charge neutrality gives us 1 equation. Thus we have NV equations for NV unknown values.

CHARGE EQUILIBRIUM OF GAS AND DUST
In this section we specify the current terms of the dynamic equations (16-22), especially the dust-dust collisional charging terms JA,C −J

Fluffy dust model
We use model of dust aggregates by Wada et al. (2008b). We consider dust aggregates composed of a large number of spherical monomers with radius rm = 0.1 µm. Each dust species I has its mass mI, the number of monomers that constitute the dust NI, and representative radius rI. We define the fractal dimension of the fluffy dust DI in the following simple manner: The dust mass is expressed in terms of monomer mass mm as follows: Wada et al. (2008b) studies the collision of the fluffy dust of the radii 1.0 × 10 −5 ∼ 9.1 × 10 −4 cm. The effect of offset collisions, collision between dust of much different sizes, and dust much larger than 9.1 × 10 −4 cm are yet to be confirmed. Therefore we make the following assumptions on smaller dust-larger dust collision.
• If the smaller dust graze at the larger dust, i.e. if the line that passes the gravitational centre of the smaller dust and is parallel to the relative velocity vector do not intersect with the larger dust, the two dust aggregates do not stick to each other. Therefore the grazing cross section is of the order of r S r L . In this case they separate ∆q A , C of charge, which is the product of charge surface density σ ch and contact surface area S kiss . This contributes to the dust-dust charging current, JC,A.
• If the smaller dust bump into the larger dust, i.e. if the line that passes the gravitational centre of the smaller dust and is parallel to the relative velocity vector do intersect with the larger dust, the smaller dust do not penetrate the larger dust but becomes a part of the larger dust. The cross section is of the order r L 2 . In this case all the charges the smaller dust have are removed from the smaller dust charge density and added up to the larger dust charge density. This contributes to the dust-dust neutralization current, J (n) C,A .

Collisional cross section of charged spherical object
In this section, we estimate collisional cross sections for dust. The collisional cross sections for two electrically charged Figure 2. A grazing collision between a smaller dust (the blue solid sphere) and a larger dust (the red wire-frame sphere). The smaller dust creates a trench on the larger dust (the black cylinder).
spherical particle is given by where q,q ′ is each particle's charge, T is the temperature of their relative motion and πa 2 is the geometric cross section (e.g. Spitzer 1941). Equation (31) represents the effect of Coulomb focusing: particles of the opposite charge attract each other and collide more often than when they are neutral. On the limit qq ′ a −1 > > kBT we can approximate the cross section as σcou(q) ≃ −πaqq ′ (kBT ) −1 , which is bi-linear on q and q ′ . On the other hand, cross section (30) represents the effect of Coulomb repulsion: for the collision between particles of the same charge only a portion of particles that belongs to the long tail of Boltzmann's distribution for temperature T can overcome the Coulomb barrier and collide. On the limit qq ′ a −1 > > kBT the cross section vanishes quickly, but never reaches 0.
We use Coulomb cross sections (30), (31) to estimate the event rate of gas-dust collision and dust-dust collision.

Collisional cross section and contact surface of fluffy dust
The amount of charge exchanged in a collision, ∆q A , C , is product of area of contact S kiss , upper limit of charge exchanged per unit surface area of contact σ ch , and the nondimensional efficiency factor η ch . We leave the detailed argument to determine η ch σ ch to §3.4. Here we assume that η ch σ ch is known and describe how to estimate contact surface area S kiss . Since it requires another detailed simulation to estimate S kiss qualitatively, we resort to an order-of-magnitude estimate for this part of the work.
We illustrate the collision between a smaller dust and a larger dust in Fig. 2. The smaller dust grazes the larger dust, pushes away the monomers that belong to the larger dust and creates a trench on the larger dust. The trench is a portion of the black cylinder in the figure. The radius and the length of the cylinder is r S and (r S r L ) 1/2 , respectively. Therefore, the surface area of the trench SC is of order and the number of monomers NC required to fill the surface of the trench is Their total surface area is also of the order of SC . However, S kiss ≃ r S 3/2 r L 1/2 overestimates the actual contact surface area if the large dust is so fluffy that there is not enough monomers in the trenched volume to fill the trench surface.
From the definition of the fractal dimension (28), the number density of monomers within the larger dust material is On the other hand the volume of the trench is Therfore, the number of particle contained in the trench is and their total surface area is If NF < NC , the surface of the trench is only partially covered by the monomers, and we estimate S kiss ≃ SF . On the other hand, if NF > NC , NF monomers are crushed onto the trenched surface, and since they overlap, about NC monomers will take part in the charge exchange. In this case we estimate S kiss ≃ SC. To summarize, we assume that S kiss is the smaller of (32) or (37):

Charge separation processes
There are generally two classes of possible charge separation processes in protoplanetary discs. One is surface charge exchange, where each dust has some kind of spontaneous charge separation (Kudin & Car 2008), so at the initial condition each dust charge is zero as a whole (globally neutral), but there are charge separation within the dust particles (locally charged). For example, water ice crystals tend to gather negative charge at its surface and positive charge inside. When two dust aggregates with different charge collide and melt partially, they exchange molten material and the charge included in the molten material. As a result each dust gets globally charged.
The other charge separation mechanism may be triboelectric processes (e.g. Desch & Cuzzi 2000). In this case, at the initial condition each dust is both globally and locally neutral. When two dust aggregates made of materials with different electron affinity collide, the surface electrons move from one material to the other. As a result each dust gets globally charged.

Surface charge exchange I -larger dust is anionic
The mechanism we consider the most plausible for the dustdust collisional charge separation is surface charge exchange between ice dust. For the dust aggregate of ice mantled silicate, Cuzzi & Zahnle (2004) proposed a condensation scenario, that at the snow line ice larger dust drifting inward dissociate and many smaller dust form.
There are established models on charge separation caused by ice-ice dust collision in the context of thundercloud meteorology (for review, see e.g. Dash et al. (2001)). We will carefully import them as a charge separation model in protoplanetary discs. The essential steps to cause lightning on earth are (1) spontaneous charge separations on ice crystal surfaces, (2) existence of different dust species with different spontaneous charge separation per surface area, (3) collisions between the different dust that leads to global charging of each dust and (4) relative motion between the globally charged dust to create electrostatic field.
For (1), we argue that the charge separation per surface area is quantitatively the same as the values measured in laboratory experiments. For (2), dominating dust species in charge separation process in protoplanetary discs is uncertain, and we discuss two possibilities (c.f. §3.4.1, §3.4.2 ) in this work. For (3) and (4), we make simple estimations for the collision rate and relative velocity in protoplanetary discs.
Ice crystal surface is intrinsically charge-separated. Ice is negatively charged near the surface, and the inside is positive. The typical charge surface density for stable ice surface is σ ch ≃ 3.0 esu cm −2 or σ ch ≃ 6.2×10 9 e cm −2 and the typical skin depth of the charged layer is d ch ≃ 2.0 × 10 −4 cm, though charge surface density for fast-growing ice surfaces are larger and shallower (Dash et al. 2001). This charge separation has a general explanation as a result of interaction between hydroxide(OH − ) and hydronium (H3O + ) ions and a hydrophobic surface (Kudin & Car 2008), and the above value of typical charge surface density is observed at liquid water-air surfaces as well as at ice crystal-air surfaces (Takahashi 2005). Therefore we use the value for ice-vacuum surfaces as well.
In the thundercloud, there are varieties of ice crystals with different surface charge densities, depending on the surface history of the ice crystals. Newly formed surfaces have larger charge surface density than old surfaces, because they have higher fractal dimension and deeper amorphous layers.
We now consider how surface charge exchange works in the model of Cuzzi & Zahnle (2004). Larger dust that migrate towards the snow line has old surface and has less negative charge surface density, while smaller dust formed at the snow line have new surface and larger negative charge surface density, as in meteorological case. Note that before collision each dust is globally neutral.
At the collision, the surface of the dust aggregates melts and the surface charge density is exchanged, and averaged. The larger dust, having less surface charge density than the smaller dust, receives more negative charge than it gives. Therefore the larger dust becomes anionic, smaller dust becomes cationic.
Laboratory experiments (Takahashi 1978), in-situ observations and meteorological estimates (Gaskell et al. 1978;Christian et al. 1980) suggest that for mm-size ice crystals, at least 10 per cent of the total surface charge within contact surface is exchanged in a single collision; experiments by Mason & Dash (2000); Dash et al. (2001) suggests almost η ch = 1.0. As a conservative estimate, we use η ch = 0.1 unless mentioned otherwise.

Surface charge exchange II -larger dust is cationic
It may be possible that charge separation processes occurring in protoplanetary discs are different from those occurring in the terrestrial thunderclouds. The collision time-scale in the protoplanetary discs is much longer than that in a thundercloud, so long that sintering may take place (Sirono 1999). As a result, The surface state of old ice larger dust and young ice smaller dust might resemble each other. If they are identical, some random charge exchange by collision is still possible, but they do not exchange charge on average.
However, compared to thundercloud, protoplanetary discs are more dirty and fine-grained; they contain much dust made of materials other than ice such as silicates, and the monomer size is 0.1 µm rather than 1 mm. Since the monomer size is smaller than typical skin depth of the charge separation d ch ≃ 2.0 × 10 −4 cm mentioned above, it is possible that ice smaller dust and silicate smaller dust with thin ice mantles formed at the snow line is inefficient in separating charge. There may be silicate aggregates with no surface charge separation. Meanwhile old larger dust that have travelled from the far end of the protoplanetary disc have undergone sintering and have developed thick mantles with full surface charge separation.
In such scenario, the larger dust has more surface charge separation than the smaller dust. Therefore, collision between a larger dust and a smaller dust still leads to charge separation but the larger dust becomes cationic, and the smaller dust is anionic in this case. We assume that η ch = 0.1 and σ ch ≃ −3.0 esu cm −2 in this case (The charge exchange rate has the same magnitude but the opposite sign compared to that of §3.4.1.) Both scenarios, the larger dust is anionic and the larger dust is cationic are plausible. They may even take place in the different parts of the same disc simultaneously. Therefore, we have decided to take both scenarios into consideration. To that end, we treat the concept of cationic and anionic dust separately from the size of the dust.

Triboelectric charge separation
Desch & Cuzzi (2000) have proposed that collision between large silicate grains and fine iron metal grains leads to triboelectric charge separation. For instance, silicate dust of radius 3.0 × 10 −2 cm will gain 5.4 × 10 3 e charges per dust. The process can be built into our model in the same manner as we treat surface charge exchange processes.

Relative velocity
When a cloud of positively and negatively charged dust is separated much larger than plasma Debye length the electrostatic field between them become observable. In order to cause such macroscopic charge separation, there must be a significant relative bulk motion between anionic and cationic dust. Inward migration of large dust is a source of this bulk motion. The sedimentation may act in the same way. Also Desch & Cuzzi (2000) have proposed that largest eddies in turbulence of protoplanetary discs cause bulk motion between smaller dust and larger dust. Such effects on the relative velocity between dust species in MMSN has been studied (see Brauer et al. (2008) and references therein).
Here, we simply assume that the largest contribution to the smaller dust-larger dust relative velocity is the bulk motion of the larger dust, and the velocity is ∆v L , S ≡ u L ≡ 3.4 × 10 3 cm s −1 , the catastrophic collision velocity of the ice dust aggregates of 9.1 × 10 −4 cm size dust (Wada et al. 2008b). Note that the non-sticking velocity threshold decrease as the monomer size increase (Blum & Wurm 2000). We also check our analytic formulae with smaller values of ∆v L , S and u L assumed.
Dust migration speed are comparable to this value at some stages of the dust growth. On the other hand, turbulent motion is faster than the value for most of our parameter range (c.f. Table 3). Turbulent mode that is larger than the scale of interest can be treated as bulk motion, and can be used to explain the charge separations of the scale. The scale can be as large as of order of disc scaleheight (Balbus & Hawley 1991).

The charge equilibrium equations
By substituting the results of analyses up to here into (12-15) we have the following dynamic equation for charge transport: where the current density terms (16-22) become: Ji,e = −eζng.
In (44), the amount of current exchange ∆q A , C is product of contact surface area S kiss and surface charge density σ ch , each described in §3.3 and §3.4. The contact surface area S kiss is the function of dust radii and dust fractal dimensions; see equation (38). The surface charge density σ ch depends on the dust material. The relative velocity of the larger dust and the smaller dust is ∆v L , S = 3.4 × 10 3 cms −1 , as we have discussed in §3.5. The cross section term σcou is the Coulomb cross section introduced in §3.2. We assume vi and ve to be thermal velocities of ions and electrons. For ionisation in MMSN at r = 2.7 AU, cosmic ray ionisation is the main contributor and ζ ≃ 10 −18 (Umebayashi & Nakano 2009). We introduce the nondimensional dust number density η (dust number density in unit of MMSN values), so that in equations (7-9), ηg = 1, and η S = η L = η. From those density term, the number density terms ng, n S , n L are given as ρg/mg, ρ S /m S , ρ L /m L . The masses of dust aggregates m S , m L are function of their radii and fractal dimensions; see equation (29).
All the variables that appear in the current density terms (44-49) are controlled by five parameters; radii of the dust aggregates (r S , r L ), their fractal dimension (D S , D L ), and the nondimensional dust number density η.
The equilibrium equations (23-27) become: J L ,e + J S ,e + Ji,e = 0 (53) Again note that, out of four Kirchhoff's Laws (50-53) only three of them are independent, and the charge neutrality condition (54) is necessary.

CRITICAL DUST NUMBER DENSITY FOR LIGHTNING
In this section we derive the strength of electric field generated by the relative motion of the large and small dust, and set conditions for macroscopic electric discharge events, or lightning.
Lightning occurs when the maximum electric field in the plasma Emax exceeds the critical value E dis . The critical electric field E dis is determined by the condition that an electron accelerated by the field has kinetic energy large enough to ionise a neutral gas molecule. Let l mfp be the mean free path for electron. Then an electron accelerated in electric field of strength E receive the energy of order e E l mfp . The ionisation potentials ∆Wion for H, H2, and He molecules are 13.6 eV, 15.4 eV, and 24.6 eV respectively (Duley & Williams 1984). We use ∆Wion = 15.4 eV in this work. Therefore the critical value E dis of electric field for the lightning satisfies: Next we derive the value of Emax. When the differential motion between the oppositely charged dust species continues much longer than the plasma Debye length, it can be interpreted as current carried by the dust jD generating electrostatic field, and the plasma counter-current jp is induced in the neutralizing direction . We consider that jp is carried by electrons, and neglect current carried by positive ions because it is at most the same order as that by electrons. Moreover, even if positive ions are accelerated to ∆Wion and ionise other molecules, they increase the electron number density only linearly, not exponentially.
The dust current jD is estimated simply, by the product of dust charge density Q L and macroscopic motion u L , as: On the other hand the particle current jp is determined by the Ohm's law: where ν is the electric conductivity, Emax is determined at the equilibrium of these two currents jD and jp: By substituting (57), (58), and (59) into (60), we obtain Now that we know both Emax and E dis , the condition for electric discharge is By substituting (56) and (61) into (62) , we have the following form of the condition for electric discharge: Within our parameter range of interest, the behaviour of the left hand side of (63) as we increase η is that it first keeps values much smaller than the right hand side and then it monotonically increases (c.f. Figure. 3, 4). Thus there is a unique value of η at which the equality for (63) holds. We define this value to be ηcrit, the critical dust number density at which lightning takes place. Note that the condition doesn't depend on the detail of the electron stopping processes because we can eliminate l mfp from the condition.

RESULTS
We have performed two sets of numerical experiments. In the first set of experiments, we fixed the set of parameters, r S , r L , D S , and D L to some typical values. We varied the dust number density η, and calculated charge density for each species of particles at the equilibrium.
In the second set of numerical experiments, we varied the set of input parameters, r S , r L , D S , and D L , and for each set of input parameters we calculated the dust number density required to cause electric discharge ηcrit.
For all these simulations we assumed the environment at the equatorial plane and the snowline of the MMSN model; r = 2.7 AU, T ⊲⊳ = 1.7 × 10 2 K, ρ ⊲⊳ g = 1.6 × 10 −10 g cm −3 , ρ ⊲⊳ S = 1.6 × 10 −12 g cm −3 , ρ ⊲⊳ L = 1.6 × 10 −13 g cm −3 . The results of the first set of experiments are in §5.1. We found that the dust-plasma system experience four phases as we increase η. We interpret this result in §5.2. The results of the second set of experiments are in §5.3. We derive the analytic formula for ηcrit in §5.4.

Equilibrium charge density of particles as a function of dust number density
We found that as we increase η while keeping other dust parameters constant, the equilibrium charge densities QI = qInI experience four phases (Table 4). Fig. 3 and Fig. 4 shows the typical four phases behaviour. In this and the next sections, we explain the origin of the four phases, using the circuit diagrams (Fig. 5) as a great help. The four-phase behaviour we describe here is independent of most of the details of charge exchange processes. In fact Fig. 3 model and Fig. 4 model have the opposite sign for dust-dust collisional charge exchange, but the evolutions are almost similar. The rest of the discussion in following sections is based on the former case, which we consider is most plausible (see §3.4.1). The discussion is easily generalized to the other case.
To analyse the result, we first identify the dominant processes by comparing the competitive current in circuit diagram, then write down all the unknown values in simple polynomials of η. Fig. 5 illustrates the transition of dominant process in the circuit as dust number density η increases. The two particles with the largest charge density is marked by larger circle. There are always two of them, one carrying most of the system's positive charge and the other negative, thus charge neutrality holds. The arrows and their line width represents direction and amount of currents. Labels for dominant currents are marked with thick rectangle, sub-dominant currents with thin rectangle, negligible cur- rents with dashed rectangle. The names and conditions for each phase is listed in Table 4. There are two major consequences of the size difference. Larger dust is much fewer in number density. So in the fewer dust limit (η < < 1) the larger dust carries much less charge density than smaller dust do. Since larger dust is the fewer, one larger dust collides with smaller dust much more often than one smaller dust does with larger dust. Therefore larger dust are the species that experience the quick charge density raise in (c)charge-up phase. The main role of the smaller dust is to absorb plasma and keep the charge neutrality.

Four phases of charge separation as a function
of dust number density

Ion-electron plasma phase
In ion-electron plasma phase ( Fig. 5 (a)), the dominant path of charge transfer is the next-dominant path is Therefore, we have following current hierarchy: The amount of current for path (64) is constrained by edge e − → i + ; since we have assumed that ζ and ng is independent of η, so is Je,i.
From charge neutrality (27), Qe = Qi and therefore ne = ni. So equation Ji, C ≃ J C ,e is satisfied by setting, in equations (19) and (21), Equation (68) tells us that σcou(q C , e)/σcou(q C , −e) is constant of η. This means q C ∝ η 0 because the only η-dependent term in σcou is q C . By definition of dust number density factor η, n C ∝ η 1 , so Q C ∝ η 1 . By similar argument we can deduce In other hand, to satisfy Ji, In this phase, ions and electrons are the major carriers of positive and negative charge. Equation (68) also tells us that σcou(q C , e)/σcou(q C , −e) ≃ ve/vi > > 1. This is interpreted as follows: Since thermal velocity of electron is much faster than that of molecular ions, electron is more rapidly absorbed to neutral dust than ions. Therefore dust continues to acquire negative charge, until its negative charge is enough to repulse most of the electrons inflow to attain a current equilibrium. Both cationic and anionic dust are forced to charge negative to hold back the overwhelming electron absorption.
To summarise,

Ion-dust plasma phase
The system enters ion-dust plasma phase when the negative charge in dust Q C become comparable to that in plasma Qe.
Charge neutrality (27) requires free electrons to decrease. So the Coulomb barrier of dust species become weaker until Coulomb cross section approximates geometric cross section σcou(q C , e) ≃ σcou(q C , e) ≃ πa C 2 ∝ η 0 where electrons and ions are equally absorbed to the dust.
In ion-dust plasma phase ( Fig. 5 (b)), the dominant path is still and the next-dominant path is still and the same current hierarchy holds: However, now that σcou(q C , e) ≃ σcou(q C , e), equation Ji, C ≃ J C ,e is satisfied by setting, in equations (19) and (21), So the ratio ni/ne is kept constant to ve/vi = 6.1 × 10 1 . Still, in order to have Ji, C ∝ η 0 and ≃ J C ,e ∝ η 0 we need ni, ne ∝ η −1 . Since qi, qe ∝ η 0 , we have Qi, Qe ∝ η −1 .
In this phase the cationic dust carry most of the negative charge while ions carry most of the positive charge of the system. Therefore, the charge neutrality equation (27) is dominated by these two components, and Q C ∝ Qi ∝ η −1 .
In this phase anionic dust also feels the same environment as cationic dust, so Q A ∝ η −1 . However as η approaches to (c)charge-up phase, dust-dust collisional charge separation J A , C gradually comes into play and Q A increases. Therefore in Fig. 3 we can see the power law Q A ∝ η −1 only at the beginning of (b)ion-dust plasma phase.
In ion-electron plasma phase and ion-dust plasma phase the dust-dust collisional charging is ineffective. So we can understand these two phase without dust-dust collisional charging (see Okuzumi (2009) and references therein.)

Charge-up phase
The system enters (c)charge-up phase when J A , C becomes larger than J A ,e. Now anionic dust has their own negative charge supply from dust-dust collision, their negative charge grow quickly, and σcou(q A , −e) become rapidly small. At this point, the circuit switches one of its current path.
In charge-up phase (Fig. 5 (c)), the dominant path is still but the next-dominant path is The amount of current for path (82) is constrained by edge e − → i + ; since we have assumed that ζ and ng is independent of η, so is Je,i. The amount of current for path (83) is constrained by edge A → C (16); since we have assumed that ∆q A , C is independent of η, J A , C ∝ η 2 . Therefore, we have following hierarchy: The path (82) is as same in ion-dust plasma phase, leading to Qi, Qe ∝ η −1 , and charge neutrality requires Q C ∝ η −1 .
In dust charge-up phase, however, anionic dust has so much charge that electrostatic potential for electron and ion at the surface of larger dust is larger than their thermal energy; this is qq ′ a −1 > > kBT limit of the Coulomb cross section (30), (30). Thus σcou(q A , −e) → 0 and σcou(q A , e) ∝ q A in (18). Substituting ni ∝ η −1 and n A ∝ η 1 into Ji, A ∝ η 2 , we have q A ∝ η 2 and Q A ∝ η 3 .
To summarise, At this phase, by substituting equations (86) (87) into equation (61) we have The Emax has the dependency of η 4 in this phase, instead of E ∝ η 2 dependence used, for example, in Gibbard et al. (1997). Moreover, at the end of dust charge-up phase there is a steep increase in Q L and steep decrease in Qe. These means that the electric discharge condition (63) meets at smaller value of η.

Dust phase
The system enters (d)dust phase when J A , C becomes larger than Ji, C . Now the charge states of both anionic and cationic dust is governed by dust-dust collision, and the plasma component is sub-dominant to the dust. In dust phase (Fig. 5 (d)), the dominant path is the dust-dust collision is now short-circuiting. The nextdominant path is The amount of current for path (90) is constrained by edge A → C (16); since we have assumed that ∆q A , C is independent of η, J A , C ∝ η 2 . The amount of current for path (91) is constrained by edge e − → i + ; since we have assumed that ζ and ng is independent of η, so is Je,i.
Therefore, we have following hierarchy: A , C . Therefore only η dependent term q C must satisfy q C ∝ η 0 , leading to Q C ∝ η 1 . Charge neutrality leads to Q A ∝ η 1 .
The path (91) gives us Qi, Qe ∝ η −1 , same as in iondust plasma phase and in dust charge-up phase.
At the boundary of (c)charge-up phase and (d)dust phase there is a jump of dust charge. This is because when η cross the boundary dust charge grows until dust-dust collisional neutralization can compensate dust-dust charge separation.
To summarise,

Critical dust number density as function of dust parameters
We now explain the details of the second numerical experiments, where we varied the set of input parameters, r S , r L , D S , and D L , and for each set of input parameters we calculated the critical dust number density ηcrit at which the lightning strikes. The numerical results strongly suggest that the parameter space (r S , r L , D S , D L ) is subdivided into several regions, at each of which ηcrit is a simple analytic function of parameters (r S , r L , D S , D L ). The parameter ranges are with additional constraints 1.0 < η < 1.0 × 10 6 .
Constraint (101) requires that the smaller dust is smaller than the larger dust. Constraint (102) comes from empirical fact that larger dust aggregates have experienced more compactification, and have higher fractal dimension (Suyama et al. 2008;Wada et al. 2008a). Constraint (103) is cutoff value of our computation. We visualize the four-dimensional field ηcrit(r S , r L , D S , D L ) in the figures at the last of this paper, by choosing some representative points and presenting several 2-dimensional sections that passes the point. As the mass, the radius, and the fractal dimension of a dust is related by equation (29), we have some freedom of choosing the direction of 2-dimensional section. We keep mI constant when we vary DI (the dust puff up with constant mass); we keep DI constant when we vary rI (the dust mass increase with constant fractal dimension).
The third set of Fig. 8, is the ηcrit averaged over the parameters that do not appear in the axes, to show the tendency of overall dependence on the parameters, and to demonstrate the precision of the analytic formulae.
The fourth set of Fig. 9 uses the same representative dust as in Fig. 6, but is the result of another simulations, where we are now extremely pessimistic and assume that the charge exchange is four orders of magnitude inefficient (η ch = 1.0 × 10 −5 instead of η ch = 1.0 × 10 −1 ). Even though, the number density ηcrit required for lightning has raised only by two order of one magnitude. The critical number density is ηcrit = 6.59×10 3 for the representative parameter.
The fifth set of Fig. 10 shows the averaged ηcrit for the pessimistic case η ch = 1.0 × 10 −5 . We later examine the accuracy of our formulae by fitting Fig. 10 with the formulae using correction factors determined by Fig. 8 data.

Analytic formulae for lightning conditions
In this section we derive the analytic form, of ηcrit and lightning conditions. Numerical results obtained in §5.3 are of great help in deriving these analytic formulae. We show at the end of §5.4.4 that by our analytic formulae we can fit 364325 numerically-obtained points distributed among six decades with 21 per cent precision. Moreover, the formulae 'predicts' results of another simulation with 59 per cent precision, where charge exchange is 10 4 times inefficient. These agreements are good evidences for correctness of both numerical and analytical results.
We made plots like Fig. 3 and Fig. 4 for many points within our parameter space, and found that η = ηcrit is met at the boundary of (c)charge-up phase and (d)dust phase in most cases, and sometimes in (d)dust phase or (c)charge-up phase. Therefore we derive analytic form of η corresponding to these three cases in §5.4.1 to §5.4.3, and combine them in §5.4.4

Analytic formulae for charge-up phase / dust phase boundary
We first calculate η (cd) , the value of η corresponding to the (c)charge-up phase / (d)dust phase boundary. The boundary satisfies J L , S = Ji, S (Table 4). We use the approximations in (c)charge-up phase to find the break point of the phase. Then because we can ignore the neutralization current J (n) L , S , approximate σ L , S = πr L 2 , σ S ,i = πr S 2 (geometric cross sections) and σ L ,i = πr L 2 |eQ L | (r L kBT n L ) −1 (qq ′ a −1 > > kBT limit of Coulomb cross section (31)).
In dust charge-up phase, both ions and electrons are mainly absorbed by smaller dust, so from (52) and (53) we have and the absorption cross sections are geometric : σ S ,i = σ S ,e = πr S 2 . By substituting (106) into (105) cf. Fig. 3(c) and in equation (84).
We have come to a simple result, that ηcrit satisfies Substituting equations (104) , (108), together with n S = η (cd) n ⊲⊳ S and n L = η (cd) n ⊲⊳ L into equation (109) and solving for dust number density η, we have We have introduced a nondimensional correction factor α (cd) , a constant that does not depend on r S , r L , D S , D L . We need this to compensate the error arising from using the formulae in (c)charge-up phase to find the break point of itself. The actual value for α (cd) is in §5.4.4.

Analytic formula for ηcrit in dust phase
Next, we derive the analytic formula of the critical density η Imposing J L , S = 0 in (44), and by approximating the charge neutrality (54) with Q S + Q L = 0, we have By approximating equation (53) with J S ,e + Ji,e = 0, we have where 1 + χ = 1 + Q S e n S r S kBT ; the factor (1 + χ) comes from the Coulomb cross section (31). Substituting Q L and Qe, the equality for the lightning condition (63) We have introduced another nondimensional correction constant α (d) as we did in §5.4.1.

Analytic formula for ηcrit in charge-up phase
Finally, we derive the analytic formula of the critical density η (c) crit , where the condition for electric discharge (63) is met in (c)charge-up phase (c.f. §5.2.3, Figure 5(c)).
By approximating equations (52) and (54) with J S ,i = Ji,e and Q S + Qi = 0, we have In equation (50), we can ignore J L ,e and further ignoring the second term in equation (44), we have here we used the qq ′ a −1 > > kBT limit of Coulomb cross section (31). Solving this for Q L , we have And from equation (53) we have By substituting these Q L and Qe to the equality for the lightning condition (63), replacing n S = η We have introduced a third nondimensional correction constant α (c) as we did in previous sections.
We have plotted these analytic solutions (124-126) combined with the condition (122) in solid-line contours from Fig. 6 to Fig. 8. The red thin contour represents the parameter ranges where η (cd) contributes. The blue thick contours represents the parameter ranges where η (d) contributes, where blue solid contour means χ < 1 and blue dashed contour χ > 1. The thick yellow-sleeved red contours represents the parameter ranges where η (c) contributes. The numerical solutions, on the other hand, are plotted in colour maps and the black dashed contours.

CONCLUSIONS AND DISCUSSIONS
We have shown that as dust number density η increase, the charge density distribution experience four phases: (a)ionelectron plasma phase, (b)ion-dust plasma phase, (c)chargeup phase and (d)dust phase. The former two phases are studied in detail by Okuzumi (2009), while the latter two phases are unique results of taking dust-dust collision into consideration. We have calculated the dust number density ηcrit at which lightning strikes, as function of dust radius r S , r L and fractal dimension D S , D L numerically. Using the numerical results we have derived the analytical formulae for ηcrit: equations (122), (124)(125)(126). Because the generated electrostatic field Emax(η) grows more rapidly than estimate by Gibbard et al. (1997) in (c)charge-up phase and (d)dust phase, lightning in protoplanetary discs are possible with smaller dust number densities. We discuss the consequences in this section.

Energetics and direct observations
We estimate the total energy of a lightning event in a protoplanetary disc at r = 2.7 AU. For MMSN, the number density of the gas is 4.7 × 10 13 cm −3 in the region. The typical electron mean free path at this site is l mf p ≃ 1.2 × 10 2 cm. By equation (56) we know the critical electric field E dis ≃ 4.3 × 10 −4 G. The sphere with radius of the disc scale-height h ≃ 2.4 × 10 12 cm contains the electric energy W ≡ E dis 2 /8π × 4πh 3 /3 ≃ 4.3 × 10 29 erg. When the lightning strikes, the energy is concentrated into lightning bolt of radius w and length h, where w is related to l mf p by w ≃ 5000 l mf p ≃ 6.0 × 10 5 cm (Pilipp et al. 1992). If all the energy is used to heat the gas within the lightning bolt, the gas can be heated to 1.6 × 10 7 K.
The ultimate energy source for this electric discharge event is the gravitational energy of the accreting matter. In our model the mass accretion ratio of uncondensed larger dust isṀ = 2πrΣ ⊲⊳ L u L ≃ 3.3 × 10 17 g sec −1 . The gravitational energy released within condensation region h is L ≡ GM⊙Ṁ hr −2 ≃ 6.6 × 10 28 erg sec −1 . For the largest energy event W = 4.3 × 10 29 erg, The upper limit of the event rate is 1.5 × 10 −1 sec −1 .

Astronomically Low Frequency (ALF) Waves
The change density evolution, electromagnetic pulse, and electromagnetic waves accompanying lightning in terrestrial thunderclouds are observed (e.g. Koshak & Krider 1989;Lin et al. 1979). The typical wavelength of the electromagnetic waves are similar to the scale height of the thundercloud. These are called extremely low frequency waves. The electromagnetic waves from lightning can be basically modelled as solutions of Maxwell equations, including lightning current as a source term (e.g. Rakov & Uman 1998). When we apply these models to the protoplanetary discs, the electromagnetic wave spectrum is extend between the event duration and light crossing time of the system, or 9.6 × 10 −5 ∼ 1.2 × 10 −2 Hz. This frequency range is at least two orders of magnitude lower than any frequencies with established observational methods. It is difficult to make a fair choice for the successor to the frequency list 'very low frequency (VLF),' 'ultra low frequency(ULF),' 'super low frequency (SLF),' and 'extremely low frequency (ELF).' We opt for Astronomically Low Frequency (ALF) waves and hope that the reader will forgive us! Anyway the frequency is so low that we will need an astronomical budget to build an astronomically large detector to receive it, considering its wavelength of order of an astronomical unit.

Infrared (IR) Observations
The energy of the lightning contributes to the local heating of the protoplanetary discs, which might be resolved by advanced telescopes such as Atacama Large Millimetre Array (ALMA). The most possible observational evidence is excess of heating near the snowline. To distinguish the cause of the heating with other heating model candidates, the variability or correlation function of the heating might be useful. This is because lightning propagates at the speed of ionised electrons, which is much faster than the speed of sound.

Ultraviolet (UV) Observations
The ionisation electrons of the lightning excite various electron levels in gas molecules and dust. There is possibility of observing fluorescence photons from such excited molecules. Although the disc gas is generally expected to be thick for ultraviolet photons, there are categories of lightning that extends toward thin regions of the gas, known as sprites and elves (e.g. Williams 2001). The sprites and elves are phenomena similar to lightning observed in the mesosphere of the earth, possibly caused by electric fields induced by the thunderclouds. Fluorescence lines from such regions can be observed by future ultraviolet missions like THEIA (Spergel et al. 2009). Also, some observational results on protostellar and protoplanetary systems today have difficulties in explaining either lack or excess of UV (e.g. Nomura & Millar 2005;Chapillon et al. 2008;Pérez et al. 2008;Herczeg & Hillenbrand 2008). If excess of UV photons is observed compared to the model, it might be from the sprite discharges and elves from the surface of the protoplanetary discs; on the other hand if the chemical composition model require more UV photons than is observed, lightning hidden in the disc mid-plane might be providing them.

High Energy Gamma Rays
Detection of burst-like gamma-ray is reported from terrestrial thunder clouds. The burst precedes a cloud-to-ground lightning, lasts for ∼ 40 seconds, extends to 10 MeV. The spectrum can be interpreted as consisting of bremsstrahlung photons from relativistic electrons (Tsuchiya et al. 2007;Enoto et al. 2008). These relativistic electrons are secondary electrons generated by cosmic rays, and accelerated by the electric fields through process known as avalanche amplification (Roussel-Dupré & Gurevich 1996). If a charged particle is accelerated by the protoplanetary thundercloud field, through similar process, its kinetic energy reaches eEh ≃ 3.1 × 10 11 eV.

Chondrule heating by lightning
Chondrule heating by lightning scenario is now considered unlikely (Weidnschilling 1997;Gibbard et al. 1997;Güttler et al. 2008). The reasons that prohibit the scenario can be summarized as following three problems.

Energetics problem
The ultimate energy source (gravitational potential of the protoplanetary disc) is sufficient to melt the chondrules; but most of the energy earned by ingoing larger dust go to the outgoing gas by angular momentum exchange (Weidnschilling 1997); little contribute to the random motion, the energy source for the lightning.

Neutralization problem
Unlike the earth atmosphere, the protoplanetary discs are filled with weakly ionised plasma which rapidly responds to electric field. Neutralization effect can be further subdivided to microscopic neutralization of individual dust and macroscopic neutralization of large-scale electric field necessary to cause lightning. If a dust get charged by dust-dust collision, the dust absorbs plasma of opposite polarity in ∼ 10 sec and returns to equilibrium charge state. Moreover, even if there is charged dust and bulk motion between the oppositely charged dust, the electric field caused by the dust induces Ohmic current in the plasma. The current will quickly neutralize the electric field.

Destruction problem
After all, there is an experimental evidence by Güttler et al. (2008) that lightning destroys the dust aggregates rather than melting them.

Solution to the problems
This work can provide answer for the first and second problem. energetics problem, the larger dust and the gas (containing smaller dust that are coupled to the gas) is now 'harnessed' by electric field. Outgoing gas is not free in carrying the gravitational energy away; instead the gas converts its gravitational energy into electric field energy, fully contributing to lightning. For the neutralization problem, we have shown in this work that with reasonably high dust number density η, the dust-dust charge separation can dominate over the plasma neutralization, and the electrostatic field can grow up to critical value.
For the third problem, we point out that in Güttler et al. (2008)'s experiment, either the electron mean free path is many orders of magnitude shorter, or the electron kinetic energy is much larger compared to the protoplanetary-disc environment. They used air at pressures between 10 and 10 5 Pa. Air consists of 78 per cent nitrogen, 21 per cent oxygen, and 1 per cent argon. Their molecular van der Waals radii are 1.6 × 10 −8 cm, 1.5 × 10 −8 cm, and 1.9 × 10 −8 cm, respectively (Bondi 1964).
Therefore, the electron mean free path and the electron kinetic energy, We = e E l mfp , was l mfp ∼ 4.8 × 10 −1 cm, We = 1.6×10 4 eV for 10 Pa case, and l mfp ∼ 4.8×10 −5 cm , We = 1.6 eV for 10 5 Pa case, respectively. On the other hand in protoplanetary discs, typical mean free path and electron kinetic energy are l mfp = 1.2 × 10 2 cm, We = 15.4 eV.
It might be possible that protoplanetary-disc lightning is effective in melting dust aggregates, although experimental lightning is ineffective in heating and led to disruption of the dust, due to shorter mean free path or higher energy electron. The minimum size of the structures that electron can form is of order of its mean free path. If the electron mean free path is much shorter than the dust aggregates, as in 10 5 Pa case, the electron current may concentrate on the most conductive part of the dust aggregate, leading to partial heating and explosion of the dust. On the other hand if the electron is much more energetic, as in 10 Pa case, it may react differently on dust monomers.
To reproduce the mean free path and electron energy simultaneously, one must reproduce the electric field strength E = 4.3 × 10 −4 G of protoplanetary discs; while the electric field used in the experiment E = 1.1 × 10 2 G was much stronger. This much stronger electric field itself, might be the cause of dust aggregate dissociation, due to much stronger electric force exerted on electron-absorbed dust monomers. Also the discharge time-scale in the experiment was much smaller than that in the protoplanetary discs, which might have led to the catastrophic results.
We think that the effect of lighting on dust aggregate in protoplanetary-disc environment is yet to be confirmed in future experiments and simulations.

Effects on magnetorotational instability (MRI) and disc environment
The dust-dust collisional charging and lightning is not a sideeffect of some other processes, but is one of the key processes in protoplanetary discs that affects each other. The lightning is powered by gravitational energy of the migrating larger dust. The migration of the larger dust as well as the long term evolution of the gas disc is governed by the disc viscosity. The best candidate for providing the disc viscosity is MRI. And MRI is controlled by gas ionisation degree, which in turn is controlled by the dust charge state and lightning. Even the longest estimate for time-scale of the lightning 1.0 × 10 4 sec is much smaller than the time-scale of MRI, which is at least of the order of Kepler timescales. Lightning occur in low-ionisation regions where MRI is prohibited (dead zones), and even if the lightning instantly raise the ionisation rate, the free electrons and ions will quickly be absorbed by the dust. Therefore we expect that MRI and lightning cannot co-exist. However lot of profound phenomena are possible. Just for an example let us think of a twolayer dead-active zone model but with dust-dust collisional charging. The dead-zone is filled with lightning, inducing sprite discharges towards active zones, which sustains the ionisation rate and MRI. The MRI in turn shovels the dust into the dead-zone.
Such global models are beyond the reach of this paper. Nevertheless we conclude this paper by stating that the dust-dust collisional charging is a necessary component for understanding the planetesimal formation and global behaviour of the protoplanetary discs.           Figure 10. Value of η crit as function of r S , r L , D S , and D L . Parameters do not appear in x-axis or y-axis are uniformly averaged over the parameter range as in Fig. 8, and we assume η ch = 1.0 × 10 −5 as in Fig. 9. Numerical results are in colour maps and black dashed contours; analytical values in coloured solid contours (c.f. §5.4.4 for the details of the plots.)

APPENDIX B: SIMULATIONS
In this section, we briefly describe our numerical methods. We need to solve the equilibrium equations (50-54), for various environmental parameters. Especially we vary η for each set of other parameters. Then we know the minimum η that satisfies the electric discharge condition (63), for the each set of other parameters. This kind of problem, a massive parameter parallelism, is typically suitable for massively parallel computing hardware (e.g. Ford 2009), such as general purpose graphic processors (GPGPUs) or GRAPE-DR (Makino 2008). We describe the CPU and GPGPU based programmes we used in this research to solve equations (40-43) in this section.

B1 Direct integral solver
The most straightforward means of finding the equilibrium solutions (40-43) is to directly integrate the dynamic equations (40-43). Nevertheless, constant-time-step direct integral cannot solve (40-43) correctly for some of the parameter range. This is because the current densities QI differ many orders of magnitude for such parameters. We must choose the integration time-step dti carefully. This leads us to the use of a adaptive time step.
The adaptive-time-step direct-integral solver follows the dynamic equations (40-43) in terms of discretized time ti where time ti is incremented by dynamic time-step dti: Our choice of the dynamical time-step dti is as follows: The direct integral solver is reliable, in sense that it is less prone to implementation mistakes because it almost straightforwardly reflects the equations (40-43), and that out of possible many equilibrium solutions (40-43) the solver will always find the desired equilibrium.
However, as η become much larger or much smaller than unity, we have found that charge distribution in the system get unbalanced. As we try to update the species I with little charge but large current, the dynamic time-step dti (B4) become the time-scale the equilibrium is reached, and the simulations becomes time consuming. Use of higher-order integral schemes are futile because we cannot take timestep much larger than dti. In addition to that, floating point numbers mainly available on GPU are single-precision floating point numbers. Computation of double-precision floating point numbers are either not supported or order of magnitude slower on common GPU.

B2 Binary search solver
So, we need to find an alternative method to solve the equilibrium equations (50-54) without directly integrating the dynamic equation, avoiding the addition between numbers of different magnitude as long as possible.
Binary search is a common method for finding zero point of a one-parameter function f ; solving equation f (x) = 0 for x. To solve the system of equations (50-54), we divide the problem into set of one-parameter problems, and conquer by recursive use of binary search. In doing so, we have to be careful in choosing which of equations (50-54) we use to find zero point of which freedom QI. Wrong choice leads to wrong result.
First, Qi and Qe can be analytically expressed in terms of Q S and Q L as follows: Qe = −eζng (n L σ L ,e + n S σ S ,e) ve .
Where we made abbreviations σ L ,i = σcou (q L , e, r L , kBT ) , σ L ,e = σcou (q L , −e, r L , kBT ) and so on. Further eliminations of freedoms is possible but complicated, because of complicated and sign-sensitive form of the Coulomb cross sections (30-31). Instead we resort to numerical methods to find out the equilibrium solution Q S and Q L for each η, and then find ηcrit, the minimum η that satisfies the electric discharge condition (63). We now describe how to solve the system of equations (50),(51),(54), and to find ηcrit, provided that for any oneparameter f we can solve f (x) = 0.
Let us name the left-hand-sides of equations (50), (51), (54) as f L , f S , and fΣ. We eliminate Qi and Qe from these using (B5) and (B6), and regard them as functions of η, Q L , Q S as follows: fΣ (η, Q L , Q S ) ≡ Q L + Q S + Qi (η, Q L , Q S ) We have also defined fcrit according to (63). For each fixed set of η and Q L , fΣ(η, Q L , Q S ) is a oneparameter function of Q S . According to the assumption we can solve fΣ(η, Q L , Q S ) = 0 for Q S . We define Q 0 S (η, Q L ) to denote the solution, so that fΣ η, Q L , Q 0 S (η, Q L ) = 0 (B13) holds.
Then fcrit(η, Q 0 L (η), Q 0 S (η, Q 0 L (η))) is a one-parameter function of η. According to the assumption we can solve fcrit(η, Q 0 L (η), Q 0 S (η, Q 0 L (η))) = 0 for η. We define η 0 to de-note the solution, so that holds, which is the ηcrit we are looking for. With this method, whenever a solver approaches the zero point of one of f , the f will consist of two or more terms of same magnitude, and of other terms of smaller magnitude. Smaller terms are irrelevant to the equilibrium. So we will always be comparing the terms of same magnitude. Because of this, the method gives sufficiently precise solutions even with single precision floating point numbers.

B3 Binary search solver on GPGPU
Graphic processing units (GPUs) are processors specialized for computer graphics, widely used in personal computers, workstations, and video game devices. But as more and more realistic computer graphics have been demanded, GPUs became capable of more and more types of calculations, and finally evolved into general purpose graphic processing units (GPGPUs), who are programmable for general computation, not limited to graphic processing. Due to the nature of graphic processing tasks, GPUs' parallel computation capacities are are one or two orders of magnitude larger compared to that of CPUs. On the other hand marketplace competition and mass production keep the GPUs' price low. Although parallel programming has been a hard task for programmers, the parallelism found in nature, together with GPGPU's power and price makes it very alluring as next generation computational platform for computational astrophysics, and computational natural science. Use of GPGPU have already started in several fields of astronomy and astrophysics, such as signal processing (e.g. Harris et al. 2008;Wayth et al. 2009), N-body simulations of gravity (e.g. Hamada & Iitaka 2007;Belleman et al. 2008;Moore & Quillen 2008), gravitational lensing (e.g. Thompson et al. 2009), orbital dynamics (e.g. Ford 2009), radiation-transfer (e.g. Jonsson & Primack 2009), and also in various other branches of science (e.g. van Meel et al. 2007;Andrecut 2008;Barros et al. 2008;Januszewski & Kostur 2009).
With GPGPU we can challenge problems that had been computationally formidable. To begin this challenge, we have constructed Tengu, (Tenmon-GPGPU cluster; GPGPU cluster for astrophysical purposes. It consists of 10 computer nodes, each equipped with NVIDIA's GPGPU. We use the programming language cuda to write codes for GPGPUs. cuda is compatible with c++, so we benefit both from expressive power of c++ and computational power of GPGPUs.
Thanks to this, we organize our c++ and cuda codes in the following way. We made c++ classes representing the protoplanetary disc, dust plasma, problem initial conditions and solutions, and the numerical solvers. Each solver inherit from an abstract solver class. Thus the users of the solvers, namely the programme parts that carries out tests and numerical experiments can use any solver they prefer, without noticing what algorithm the solvers use nor on what hardware they run. We write most of the code in c++ and compile them by gcc. The GPGPU related details are separated in several .cu files by means of pimpl idiom. We compile .cu files and link the object files using nvcc, cuda compiler provided by NVIDIA.
Another example of such benefit is thrust (http://code.google.com/p/thrust/), a cuda counterpart of what standard template library (STL) is in c++. With thrust, for example, device and host memory management is automated. Memories are allocated and freed automatically in the constructor and destructor of the container classes. Copying data between host memory and device memory are simply expressed by substitution = operators.

B4 Testing
We choose the Test-Driven Development style for this study. We have tested that the charge conservation and the current conservation conditions are held, for each equilibrium solution that each solver give. We have also tested that the value of the currents satisfy equations by comparing them with the simplest implementation.
Why do we test our codes? We need tests in numerical physics because the codes must compile correctly, the codes must translate the algorithms correctly, and algorithms represent the physical concepts correctly.
In order to assure these, we are accustomed to perform various tests during a code development, by examining the internal states and outputs of the programme. Furthermore, we want to make sure that criteria once tested always meet thereafter, and that the tests cover all the important aspects of the code. As the code grow, it becomes more and more effective to build up a system of test rather than to perform tests manually. This technique is known as Test-Driven Development (e.g. Erdogmus et al. 2005).
The programme is divided into many functional modules, and we test that these modules give expected output for given input. This is compared to code-reading style of tests, where testers finds faults in the code by reading them. Although code reading is capable of finding more mistakes (Basili & Selby 1987) it is only effective when the code is short and it is possible to trace the comprehensive behaviour of the code line-by-line. Unit tests, on the other hand cares only on the input and output. It is effective even if we are trying new languages or hardware, or we cannot debug trace on them. We use googletest, Google's framework for writing automated c++ tests (http://code.google.com/p/googletest/). We must also consider the time cost of the test. Because we do computationally heavy tasks, examining the entire behaviour of the programme is not practical. With systematized tests, we can ensure the equivalence of the codes as we optimize them, or as we transplant it onto another language or hardware. We further construct a 'test ladder,' an analog of distance ladder in cosmology. We develop a chain of successively faster algorithms, and feed them with randomly generated inputs and check if their response is same up to required precision. At the one end of the ladder is a code that almost directly trace the equations, correctness of whose implementation is self-evident. At the other end of the ladder is optimized, massively parallel code running on GPGPU.
Finally, we evaluate the computational optimization achieved by measuring the speed of each solvers in terms  Table B1. The name of the computation runs, the size of initial conditions sets, and the wall clock time needed to solve η crit for all initial conditions in seconds. The speed of the codes are also listed, in terms of number of solved cases per time, relative to the first case.
(1) We used the direct integral solver ( §B1), running it in parallel on 40 CPU cores.
(2) We used the binary search solver ( §B2), on a Core 2 Quad 9300 CPU in single thread.
(3) We used the cuda version of binary search solver, ( §B3), on single GTX280 GPGPU. (4) We used the same programme ( §B3) and the same GPU for actual numerical experiment. The run generates a set of data that corresponds to η crit as function of dust radius r S , r L and fractal dimension D S , D L . Or it corresponds to one page of the result figure, e.g. Fig. 6.
of how many problems they solve per wall clock time. See Table B1 for optimization results.