© The Authors 2023

1 Introduction

Dust grains in protoplanetary disks play a crucial role in disk evolution and planet formation (see, e.g., reviews by Lesur et al. 2023, Drazkowska et al. 2023). The dust grains in the midplane tend to pick up free electrons, which changes the character of the chemistry and has been shown to be connected to the degree of ionization (Bai & Stone 2011; Ivlev et al. 2016). Thus, grain charges are important for dust dynamics and instabilities, such as the magneto-rotational instability (MRI; see, e.g., Salmeron & Wardle 2008; Perez-Becker & Chiang 2011; Bai 2011; Thi et al. 2019). Additionally, grain charges and MRI can play an important role in and can be used for explaining gap formations in a disk, as shown in Dullemond & Penzlin (2018), Regály et al. (2021, 2023).

Dust grains accumulate toward the midplane via gravitational settling, enhancing the concentrations of the large grains in particular (see, e.g., Dubrulle et al. 1995; Dullemond & Dominik 2004; Riols & Lesur 2018). Dust grains undergo growth as well, where the charge of the grains can play an important role. Dust grains tend to uniformly charge negatively in the shielded midplane regions, which can create an electrostatic repulsion between the grains large enough to prevent collisions and hence growth. This effect is known as the “charge barrier” (see, e.g., Okuzumi 2009; Okuzumi et al. 2011, and Akimkin et al. 2020).

In this work, we investigate if charged dust grains can trigger discharge phenomena similar to lightning on Earth. Lightning is a universal phenomenon that occurs in the atmospheres of several different types of astrophysical objects. Most prominently, lightning occurs in Earth’s atmosphere, and due to the high temperatures (Orville 1968a,b,c) associated with it, lightning has several implications on the Earth’s atmospheric chemistry. Lightning also plays a central role in the abiotic fixation of nitrogen (Noxon 1976; Hill et al. 1980), which is the creation of several nitrogen bearing species, such as NO_x (Price et al. 1997), from elemental nitrogen N₂. Nitrogen fixation via lightning is also central in the prebiotic fixation of nitrogen (Navarro-González et al. 2001), a process which made nitrogen available to the first lifeforms on Earth. Further implications of this process are the abundance of ozone in Earth’s atmospheres, as ozone is a product of reactions between NO_x and elemental oxygen (Wild 2007; Murray et al. 2013). Lightning is not just a prominent feature of Earth, and it has been observed on other planets and moons in the Solar System. The Pioneer Venus Orbiter missions revealed lightning discharges on Venus from clouds in the ionosphere (Russell et al. 2006, 2007, 2008), and the potential chemical implications have been discussed by Delitsky & Baines (2015). Cassini revealed that lightning also occurs on Saturn (Baines et al. 2009; Dyudina et al. 2010).

There has also been a vast amount of work on the theoretical front, including discussions on the existence of lightning on objects such as brown dwarfs (Helling et al. 2016; Hodosán et al. 2016) and exoplanets (Méndez Harper et al. 2018; Helling & Rimmer 2019). The findings of these discussions have revealed that lightning is a phenomenon that is ubiquitous. There has also been a significant amount of work contributed to answer the question of whether lightning is feasible in protoplanetary disks. Lightning induced by triboelectric or collisional charging could be a reasonable mechanism in forming chondrules, as shown by Desch & Cuzzi (2000). The triboelectric charging they implemented in their work is not considered in this work. Similar work has been done in Muranushi (2010), who calculated a critical dust number density at which lightning should occur. Compared to our present work, they are missing an explicit turbulence analysis. In the realm of lightning in disks, one should also mention the works of Inutsuka & Sano (2005), Okuzumi & Inutsuka (2015), and Okuzumi et al. (2019). In particular, the latter two investigate how dust grains charge in weakly ionized plasma and additionally study the heating of dust grains by electric fields, which are induced via magneto hydrodynamics (MHD). This work differs in the way electric fields are generated. They investigated how MRI could induce electric fields that heat the dust grains, while we investigate Kolmogorov turbulence. In addition, no MHD effects are considered in this work. We find that we should also mention Johansen & Okuzumi (2017), as they argue that lightning discharge can be caused by positron emission from pebbles that contain radioactive Al²⁶. In our work, the main source of ionization is cosmic rays, and we do not consider radioactive decay. In order for lightning to emerge, certain conditions have to be met. Charges need to build up and then have to be separated, and this separation and buildup have to result in a large enough electric field for an electron cascade to occur (Williams 1985; Gurevich & Zybin 2005; Saunders 2008). In clouds on the Earth, the charges are built-up by ice crystals becoming charged via collisions with each other (Takahashi 1978; Jayaratne et al. 1983; Saunders 1993; Saunders & Peck 1998). During these collisions, faster growing and therefore smaller particles tend to charge positively, whereas the larger particles charge negatively (Baker et al. 1987; Dash et al. 2001). For charge separation in clouds, an updraft is generally required, and in this updraft, a lower positive layer and an upper negative layer form (Stolzenburg et al. 1998; Bruning et al. 2007). Another way for lightning to emerge on Earth is during volcanic eruptions (Brook et al. 1974; Gaudin et al. 2018). While the fundamentals are similar to lightning in clouds, some (important) details are different. In volcanic eruptions, dust particles in the volcanic plumes are the particles that become triboelectrically charged (Houghton et al. 2013). Triboelectric charging leads to the smaller particles becoming charged more negatively and larger particles being charged more positively (Lacks & Levandovsky 2007; Forward et al. 2009).

Lastly, one should also mention the new experimental findings in triboelectric charging by Jungmann & Wurm (2021) and Wurm et al. (2022). In the first publication, Jungmann and Wurm state that triboelectric charging could be instrumental in overcoming the bouncing barrier, and in the second, the authors investigate how triboelectric charging could play a role in ionizing the gas phase. The former explanations might lead the reader to ask whether the conditions in protoplanetary disks can lead to a similar emergence of lightning. We think so, as the collection of dust in high densities found in the midplane of protoplanetary disks due to gravitational settling resembles conditions found on Earth and other planets, moons, and even brown dwarfs where lightning has been shown to emerge. In addition to this, turbulent motions, which could result in charge separation, are found in protoplanetary disks.

In Sect. 2, we explain our model, which includes explanations on the chemistry of our code, all the numerical basics and improvements we have done to our code, and how we implemented cosmic rays. In Sect. 3, we present the results of our simulations. In particular, we investigate the role electrons play in our simulations, how the dust grains become charged, and whether specific regions of the disk can be favorable for the emergence of lightning. We furthermore investigate how turbulent eddies could induce a charge separation in these areas and if electric fields that emerge from this separation could be large enough for lightning to emerge. In Sect. 4, we summarize our results. We provide a glossary of the different symbols used in this work in Table A.1.

2 The thermo-chemical disk model

To simulate the size-dependent charging of dust grains in a protoplanetary disk, we used the Protoplanetary Disk Model (PRODIMO) developed by Woitke et al. (2009, 2016). PRODIMO is a 2D thermo-chemical disk modeling code that combines detailed continuum and line radiative transfer with the solution of a chemical network, the determination of the non-local thermal equilibrium (LTE) population of atomic and molecular states, and the heating-cooling balance for both gas and dust.

The model uses a simple parametric density setup chosen to represent a 2 Myrs old T Tauri star. All further assumptions about the star; the disk shape; the irradiation of the disk with FUV photons, X-rays, and cosmic rays; the dust material composition and opacities; and dust settling are explained in detail in Woitke et al. (2016, see their Table 3). Here, in this publication, we only explain the changes made to the PRODIMO code that differ from the code used in Woitke et al. (2016), and differences to the setup from Woitke et al. (2016). The selected values for the stellar, disk shape, and dust parameters are summarized in Table 3.

2.1 Chemical setup

We choose the so-called large DIANA-standard chemical network with 235 species as our base network for this work. It has been described in detail by Kamp et al. (2017). Most UV, cosmic ray, two-body, and three-body reaction rates are from the UMIST 2012 database (Woodall et al. 2007; McElroy et al. 2013). In addition, there are rates for the H2 formation on grains, reactions for electronically excited molecular hydrogen ${\text{H}}_2^{\star }$, polyaromatic hydrocarbon (PAH) charge chemistry (Thi et al. 2019), and a simple freeze-out chemistry for the ice phases of all neutral molecules (see Woitke et al. 2009; Kamp et al. 2017). The large DIANA-standard network also includes X-ray processes with doubly ionized species from Aresu (2012).

We added a number of dust species to this base chemical network to represent the different charging states q of grains of different sizes a. By solving the chemical rate network, we determined the dust charge distribution function for every given size and charge f(a, q) from which we could then compute, for example, the mean charge of grains of different sizes, including the feedback of the charged grains on all other results of the chemistry.

In contrast to Woitke et al. (2016), we have chosen to omit the PAH charging states in the chemistry in this paper, but PAH heating is still taken into account. If present, the PAHs can introduce similar effects as charged grains. Since we want to study the chemical feedback of charged grains in disks in particular, we wanted to avoid any confusion between PAH charging and dust charging. In addition, preliminary tests showed that the impact of PAHs in the midplane is negligible, due to the PAHs freezing out. Therefore, we also omitted them to lower the numeric effort.

2.2 Dust size distribution and settling

In PRODIMO, the dust component enters into the modeling (i) as an opacity source for the continuum, (ii) its effect on line radiative transfer, and (iii) as a set of chemical species that can pick up and release charges in the chemical rate network. The dust grains, however, are neither created nor destroyed, nor do they grow in the chemistry module.

Our assumptions about dust settling and opacity are explained in Woitke et al. (2016). The dust size distribution function before settling [cm⁻¹] is assumed to be

$$f_0\left(a\right)\propto a^{-a_{\text{pow}}},$$(1)

which is the same everywhere in the disk, using N_size log-equidistant size grid points between the minimum grain radius a_min and the maximum radius a_max. The proportionality constant in Eq. (1) is derived from the local gas density ρ(r, z), the dust material density ρ_gr, and the global gas-to-dust mass ratio, assumed to be 100 in this paper. In each vertical column at radius r, all grains of size a are vertically redistributed with a scale height that is smaller than the gas scale height. Here, in this publication, we used the settling description by Dubrulle et al. (1995) with the settling parameter α_settle = 0.01, as outlined in Woitke et al. (2016). After applying the settling, the dust size distribution function f(a, r, z) is numerically available on the N_size size grid points at each point (r, z) in the model. The term f(a, r, z) is the basis for the opacity calculations.

2.3 Dust size bins in the chemistry

We used a much smaller number of dust size bins N_bin (between two and nine) to represent the position-dependent dust size distribution function after settling f(a, r, z) in the chemistry. The idea in this section is to introduce two power-law indices, κ and ζ, to make the sparse representation exact with regard to two selectable dust size moments. We then argue which of the dust size moments should be selected in order to minimize the numerical errors in the chemistry introduced by the sparse resolution of f(a, r, z) in size space. We start by computing the following moment of the dust size distribution function

$$\langle a^{\zeta }\rangle \left(r,z\right)={\displaystyle {\int }_{a_{\min }}^{a_{\max }}f\left(a,r,z\right)a^{\zeta }da,}$$(2)

where f(a, r, z) [cm⁻¹] is normalized as ${\displaystyle {\int }_{a_{\min }}^{a_{\max }}f\left(a,r,z\right)da=1}$. For example, ζ = 2, n_d 4π 〈a²〉 is the total dust surface per cubic centimeter [cm² cm⁻³]. Next, we divide 〈a^ζ〉 by N_bin and numerically adjust the bin boundaries a_j such that each size bin contains the same portion of 〈a^ζ〉.

$${\langle a^{\zeta }\rangle }_j={\displaystyle {\int }_{a_{j-1}}^{a_j}f\left(a,r,z\right)a^{\zeta }da=\frac{\langle a^{\zeta }\rangle \left(r,z\right)}{N_{\text{bin}}}.}$$(3)

Here, a_j terms are the size bin boundary values, where a₀ = a_min and $a_{N_{\text{bin}}}=a_{\text{max}}$. Once the a_j are fixed this way, we compute similar quantities

$${\langle a^{\kappa }\rangle }_j={\displaystyle {\int }_{a_{j-1}}^{a_j}f\left(a,r,z\right)a^{\kappa }da}$$(4)

and define the average size of our grains [cm] in the dust bin j as

$${\overline{a}}_j={\left(\frac{{\langle a^{\zeta }\rangle }_j}{{\langle a^{\kappa }\rangle }_j}\right)}^{1/\left(\zeta -\kappa \right)}.$$(5)

Finally, the dust particle densities [cm⁻³] in the size bin j are calculated as

$$n_{\text{d},j}=n_{\text{d}}\left(r,z\right)\frac{{\langle a^{\kappa }\rangle }_j}{{\left({\overline{a}}_j\right)}^{\kappa }}.$$(6)

From Eqs. (3) to (6), it follows that

$$\begin{array}{c}{n}_{\text{d,}j}{\left({\overline{a}}_j\right)}^{\kappa }=n_{\text{d}}\left(r,z\right){\displaystyle {\int }_{a_{j-1}}^{a_j}f\left(a,r,z\right)a^{\kappa }da}\\ n_{\text{d,}j}{\left({\overline{a}}_j\right)}^{\zeta }=n_{\text{d}}\left(r,z\right){\displaystyle {\int }_{a_{j-1}}^{a_j}f\left(a,r,z\right)a^{\zeta }da,}\end{array}$$

which means that our dust bins {(n_d,j ā_j) | j = 1, … , N_bin) exactly represent two dust size moments where the dust size distribution function is weighted with a^κ and a^ζ, respectively (see Table 1).

If κ = 0 is chosen, the quantity 〈a^κ〉_j becomes the fraction of dust particles having sizes between a_j−1 and a_j, and consequently, the total dust number density $n_{\text{d,1}}+n_{\text{d,2}}+\dots +n_{\text{d,}{N}_{\text{bin}}}=n_{\text{d}}$ becomes exact. In a similar way, choosing ζ = 2 causes the representation of the total dust surface to be exact. All other dust size moments, however, show deviations from their true integral values. Increasing the number of dust bins N_bin would reduce these problems.

Thus, with κ and ζ, we can choose two dust size moments that are represented exactly by the bins and are most relevant for the problem at hand. Since the effects of the charged grains in the chemistry, via their collisional and photoionization rates, scale with their total cross-section, a choice of ζ = 2 seems appropriate. That scaling, however, is not exact because there are second-order dependencies of the rates on grain charge over radius q/a (see Sect. 2.8). In addition, the conservation of the total charge is an important principle. The total number of elementary charges on the grains per volume, assuming q/a = const is valid (see Sect. 2.8), is

$$Q=n_{\text{d}}{\displaystyle {\int }_{a_{\min }}^{a_{\max }}q\left(a\right)f\left(a\right)da}\approx n_{\text{d}}\left(\frac{q}{a}\right){\displaystyle {\int }_{a_{\min }}^{a_{\max }}af\left(a\right)da.}$$(7)

Therefore, it seems important to get the first dust size moment correct as well. In contrast, there is no direct link between the total dust particle density and the chemistry nor to any dependencies of the chemistry on the total dust volume. Therefore, our default choice is κ = 1 and ζ = 2.

Table 1 shows that in this case, we only need to represent the smaller grains up to a few microns in size. A larger value of ζ would also come with an increased level of computational problems and computing time, for example, because we would need to extend the maximum number of charging states q_max to many thousands for large grains (see Sect. 2.8).

2.4 Dust charging reactions

The charging of dust grains is part of the chemistry of PRODIMO. We focus only on the most relevant reactions here. A full description of all considered dust grain charging processes can be found in Thi et al. (2019). The three most important reaction types are electron attachment, photoionization, and charge exchange between dust grains and molecular ions.

Photoionization. Photoionization describes the process of stripping away electrons from the dust grains via photons with sufficiently high energies.

$$\text{Z}+hv\to {\text{Z}}^{+}+e^{-},$$(8)

where Z and Z⁺ stand for grains with charges q and q + 1, respectively. Here, q can have any integer value greater than zero.

For positively charged and neutral grains, this process is called photoejection, and for negatively charged grains, this process is called photodetachment. According to Weingartner & Draine (2001), the energy threshold for photoejection is

$$h{v}_{\text{pe}}=\text{IP}=W_0+W_c,$$(9)

where IP is the ionization potential of the dust grain and W₀ the work function of silicate grains, assumed to be 8 eV by Cuzzi et al. (2001) after measurements of Feuerbacher & Fitton (1972). The term W_c is an additional term that is dependent on the charge of the grain, and it increases the ionization potential for positive grains, making it harder to ionize the grain or decreasing the ionization potential and making it easier to ionize the grain. This additional term, W_c, is defined as

$$W_c=\left(q+\frac{1}{2}\right)\frac{e^2}{a},$$(10)

where q is the charge of the grain normalized by the elementary charge, e is the elementary charge, and a is the grain radius. We note that the factor 1/2 is still open for discussion according to Wong et al. (2003), where a factor of 3/8 was proposed. For photodetachment, the threshold photon energy hv_pd is equal to the grain electron affinity EA plus the minimum energy E_min at which the tunneling probability becomes relevant.

$$h{v}_{\text{pd}}=\text{EA}\left(q+1,a\right)+E_{\min }\left(q,a\right)$$(11)

$$\text{EA}\left(q,a\right)=W_0-E_{\text{bg}}+\left(q+\frac{1}{2}\right)\frac{e^2}{a}$$(12)

$$E_{\min }\left(q,a\right)=-\left(q+1\right)\frac{e^2}{a}{\left[1+{\left(\frac{27\text{Å}}{a}\right)}^{0.75}\right]}^{-1},$$(13)

where E_bg is the band gap, assumed to be 5 eV in Weingartner & Draine (2001); a is the dust grain radius; and q is the grain charge in units of the electron charge. For the total rate coefficient of photoionization, one has to combine both photodetachment and photoejection. The rate coefficient therefore takes the form of

$$k_{\text{ph}}\left(q\right)=\pi a^2{\displaystyle {\int }_{v_{\text{pe}}}^{v_{\max }}{\eta }_{\text{eff}}}{Q}_{\text{abs}}{J}_{v}dv+\pi a^2{\displaystyle {\int }_{v_{\text{pd}}\left(q\right)}^{v_{\max }}{\eta }_{\text{pd}}}{Q}_{\text{abs}}{J}_{v}dv,$$(14)

where J_ν is the direction-averaged flux of photons per area and second [cm⁻²s⁻¹]; ν_pe is the threshold frequency of photoejection; ν_pd is the threshold frequency for photodetachment; ν_max is the maximum frequency of photons, being equivalent to an energy of 13.6eV in PRODIMO; Q_abs is the frequency-dependent absorption efficiency; η_eff and η_pd are the yields of the photoejection and photodetachment (Weingartner & Draine 2001) of silicate; and a is the dust grain radius.

Electron attachment. The rate coefficient for electron attachment,

$$Z+e^{-}\to {\text{Z}}^{-},$$(15)

is derived by averaging over the Maxwellian velocity distribution of the impinging electrons. Here Z and Z⁻ stand for grains with charges q and q − 1, respectively, where q can have any integer value less than or equal to zero. This results in a reaction coefficient k_e(q) of the form

$$k_e\left(q\right)=n_e{S}_e\sqrt{\frac{8k_{b}T}{\pi m_e}}{\sigma }_Z{f}_q,$$(16)

where n_e is the electron density, S_e is a sticking coefficient assumed to be larger than 0.3 (Umebayashi & Nakano 1980) and set to 0.5 in our simulations, T is the gas temperature, m_e is the electron mass, and k_b is the Boltzmann constant. The dust grain cross-section is represented by σ_Z = πa², and a is the dust grain radius. The term f_q is a charge-dependent factor that is derived by taking the Coulomb interaction between the incoming electron and charged grain into account. For neutral grains, we have f_q = 1. For positively charged grains (q > 0), there is an additional attraction that enlarges the rate, and for negatively charged grains (q < 0), there is Coulomb repulsion

$$f_q=\{\begin{array}{ll}1+\frac{W_c}{kT}\hfill & W_c>0\hfill \\ 1\hfill & W_c=0\hfill \\ \exp \frac{W_c}{kT}\hfill & W_c<0\hfill \end{array}.$$(17)

Charge exchange. Dust grains can get charged via charge exchanges between dust grains and molecular ions. The most prevalent charge exchange for this study is the exchange of a negative charge between a negatively charged dust grain and a positively charged molecular ion,

$${\text{Z}}^{-}+{\text{M}}^{+}\to \text{Z+\hspace{0.17em}M.}$$(18)

The general form of the rate coefficient is similar to the electron attachment. Since we only allowed for exothermal charge exchange reactions, the rate coefficient takes the form

$$k_{\text{ex}}^M\left(q\right)=n_M{S}_M\sqrt{\frac{8kT}{\pi m_M}}{\sigma }_Z{f}_q^M,$$(19)

where, analog to Eq. (17),

$$f_q^M=\{\begin{array}{ll}1-\frac{q{q}^M{e}^2}{akT}\hfill & q{q}^{\text{M}}<0\hfill \\ 1\hfill & q{q}^{\text{M}}=0\hfill \\ \exp \left(-\frac{q{q}^M{e}^2}{akT}\right)\hfill & q{q}^{\text{M}}>0\hfill \end{array},$$(20)

where n_M is the density of the molecules, q^M is the charge of the impinging molecule (usually +1), S_M is a sticking coefficient set to one (Thi et al. 2019), n_M is the number density of the molecule, and mM is the molecular mass.

We note that we excluded passive ion attachment in our network. Meaning, in our network, molecules and dust grains cannot attach to each other via electromagnetic forces only. We discuss the potential implications of including passive ion attachment in Sect. 5.

We excluded endotherm charge exchange reactions between protonated molecules and dust grains, where the work function minus the band gap plus the proton affinity of the molecule exceeds 13.6 eV (see Appendix D and Table 2). In order to identify these endotherm charge exchange reactions, we needed to know the heat of formation of negatively charged dust grains $H_f^0\left(Z^q\right)$, which is given by the work function minus the band gap plus small correction terms of order e²/a (see Appendix D). In our model, we consider porous silicate grains internally mixed with amorphous carbon (see Table 3), but these grains can also be ice coated. For pure silicate grains, values for the work function minus the bad gap lie between about 3 eV and 8 eV (Weingartner & Draine 2001). Yet, following the explanations of Helling et al. (2011) and looking at the compiled values of work functions of different dust components from Desch & Cuzzi (2000) and Kopnin et al. (2004), a value between 2 eV and 6 eV seems appropriate. In this paper, we use $H_f^0\left(Z^q\right)=q\times 5.89\text{\hspace{0.17em}eV}$, where q < 0 for negative grains, 0 for neutral grains and q > 1 for positive grains, in accordance with Thi et al. (2019). This is supported by Rosenberg (2001), who stated that a mixture of dust grain material tends to lower the work function overall. According to the data collected by Thi et al. (2019), water ice (coated) grains should have higher work functions. We can assume that our dust grains are water ice coated according to D’Angelo et al. (2019), who showed that at temperatures of 300 – 500 K, a water layer should be found on dust grains. This is further supported by Thi et al. (2020), who showed, with a PRODIMO model similar to ours, that phyllosilicates can be found in the midplane of a protoplanetary disk.

2.5 Chemical network adjustment for charge and proton exchange reactions

Particular to this work is that the code has the ability to automatically generate all possible charge exchange reactions between molecules and dust grains and proton exchange reactions between molecules and then add them to the network if they meet certain criteria. In addition, we allowed the code to add endothermic reactions to the network if they meet a predetermined requirement. We tested the implications of the different networks this creates in Appendix B.

Proton exchange reactions. Proton reactions come in the following form:

$${\text{M}}_1+{\text{M}}_2{\text{H}}^{+}\leftrightarrow {\text{M}}_1{\text{H}}^{+}+{\text{M}}_2,$$(21)

where M₁ and M₂ are neutral molecules while M₁H⁺ and M₂H⁺ are the corresponding protonated molecules. We started our procedure by first identifying pairs of molecules and their protonated counterpart (M₁, M₂H⁺). If a reaction for the pair already exists in the network, we simply moved on. However, if a reaction between the pair was not found in the network, we added an approximated reaction to the network with certain assumptions.

The aim was to create a reaction rate with the Arrhenius equation,

$$k_{\text{ex}}\left(T\right)=\alpha {\left(\frac{T}{T_0}\right)}^{\beta }{e}^{-\gamma /T},$$(22)

where T₀ refers to the reference temperature of 293 K. In order to achieve this, we had to provide the different factors α, β, and γ. For α and β, we find their values with the following method. For each unprotonated molecule, M₁, we counted how many protonation reactions for M₁ already exist in the network. We summed up the different α_i and β_i parameters of the existing reactions. We averaged the summed-up parameters by simply dividing them by the number of reactions found in the network, n_reac:

$${\alpha }_{\text{mean}}=\frac{1}{n_{\text{reac}}}{\displaystyle \sum {\alpha }_i},$$(23)

$${\beta }_{\text{mean}}=\frac{1}{n_{\text{reac}}}{\displaystyle \sum {\beta }_i}.$$(24)

If we only found one reaction, we assumed α_mean to be 5 × 10⁻¹⁰ and β_mean to be -0.5. We then checked pairs of unprotonated and protonated molecules (M₁, M₂H⁺) and created reaction rates with the calculated α_mean and β_mean and the corresponding product (M₁H⁺, M₂). Lastly, we needed to calculate γ given by the reaction enthalpy ∆H_r (Eq. (25)). This was done by taking the difference of the heats of formation of the products and reactants $H_f^0$, with data taken from measurements at 0 K from Millar et al. (1997).

$$\Delta H_r=H_f^0\left(M_1{H}^{+}\right)+H_f^0\left(M_2\right)-H_f^0\left(M_1\right)-H_f^0\left(M_2{H}^{+}\right).$$(25)

Here, two cases can arise. In case of an exothermic reaction (∆H_r < 0), we can just add the reaction and set γ = 0. In case of an endothermic reaction (∆H_r > 0), the reaction is only added when ∆H_r/k_b is smaller than a predetermined barrier. In cases where ∆H_r is small enough, γ is set to ∆H_rand the reaction is added anyway. In our standard case, we exclude endothermic reactions (i.e., we set the barrier to zero). In the appendix we investigate the impact of such endothermic reactions by using a barrier value of 5000 K.

Charge exchange reactions – molecules. For charge exchange reaction, we assume two types: non-dissociative and dissociative. We followed the same approach as with the proton exchanges, meaning we had to supply the different parameters of the Arrhenius equation. Our approach for γ is the same as with the proton exchange reactions. We calculated the difference of the formation enthalpies dH_f and added it as γ if the reaction was exothermic or below the predetermined barrier. For α and β we simply assumed the values of 5 × 10⁻¹⁰ and -0.5, respectively.

Charge exchange reactions – Dust. If dust grains were involved in a charge exchange reaction, we adjusted the rate coefficient according to

$$k_{j,\text{ex}}\left(q\right)=A\frac{\sqrt{\frac{T}{300K}}}{n_{d,\text{ref}}}{\left(\frac{a}{a_{\text{ref}}}\right)}^2.$$(26)

Here, A is the first factor in the Arrhenius equation and taken from the results from Leung et al. (1984), who calculated these rates for the reference dust particle density n_d,ref = 2.64 × 10⁻¹² cm⁻³ and reference grain radius a_ref = 0.1 μm.

2.6 Dust charge moments

In principle, each charge state q of the grains in size bin j could be included in the chemical rate network. However, since micron-sized grains can already collect thousands of elementary charges (Stark et al. 2015; Tazaki et al. 2020), this method would mean including a couple thousand species, which is computationally very challenging. In order to avoid this problem, the dust grains and all their charge states are represented by three charge moments, ${\text{Z}}_{\text{m,}j}^{-}$, ${\text{Z}}_{\text{m,}j}$, and ${\text{Z}}_{\text{m,}j}^{+}$, in each size bin j.

$$\left[Z_{\text{m,}j}^{+}\right]=n_{\text{d,}j}{\displaystyle \sum _{q=1}^{q_{\max }}\left(+q\right)f_j\left(q\right)}$$(27)

$$\left[Z_{\text{m,}j}^{-}\right]=n_{\text{d,}j}{\displaystyle \sum _{q=-q_{\max }}^{-1}\left(-q\right)f_j\left(q\right)}$$(28)

$$\left[Z_{\text{m,}j}\right]=n_{\text{d,}j}{q}_{\max }-\left[Z_{\text{m,}j}^{-}\right]-\left[Z_{\text{m,}j}^{+}\right],$$(29)

where q_max is the maximum number of elementary charges on the dust grains, both positive and negative, and n_d,j is the dust number density of bin j. The moments ${\text{Z}}_{\text{m,}j}^{+}$ and ${\text{Z}}_{\text{m,}j}^{-}$ express the total number of positive and negative charges on the dust grains in size bin j per volume, respectively. The neutral moment Z_m,j is defined in such a way that it becomes zero when all grains have a maximum charge, q_max or −q_max, and stays positive for any other charge distribution function f_j(q). The charge distribution function was normalized to one (i.e., ∑ f_j(q) = 1). According to Eq. (29), we can identify a constant quantity,

$${ϵ}_{Z_j}=\frac{q_{\max }{n}_{\text{d,}j}}{n_{\langle \text{H}\rangle }}=\frac{\left[Z_{\text{m,j}}^{+}\right]+\left[Z_{\text{m}\text{.j}}\right]+\left[Z_{\text{m,j}}^{-}\right]}{n_{\langle \text{H}\rangle }},$$(30)

which we used in PRODIMO to define a quasi-element abundance for the dust grains in size bin j. Here, n_〈H〉 refers to the total number density of hydrogen nuclei. The mean charge of the grains in size bin j is

$$\langle q_j\rangle ={\displaystyle \sum _{-q{}_{\max }}^{-1}q{f}_i\left(q\right)+{\displaystyle \sum _1^{q_{\max }}q{f}_i\left(q\right)=\frac{\left[Z_{\text{m,}j}^{+}\right]-\left[Z_{\text{m,}j}^{-}\right]}{n_{\text{d,}j}}.}}$$(31)

2.7 Solving for the charge distribution function

To determine the discrete charge distribution function f_j(q) for all charge states q in every size bin j, we employed a method called linear chains of reactions. We note that we assumed the charge distribution function to be in a steady state. The general idea is that in our model, the charge states q of a grain can only change by one increment, either forward to q + 1 or the reverse to q − 1 . This is illustrated in Fig. 1.

The processes that add a charge, q → q + 1, are photoionization (called photoejection for neutral grains and photodetachment for negative grains) and charge exchange between a dust grain and a molecular cation. The forward rates [1/ s] were calculated in the following way:

$${\text{forward}}_j\left(q\right)=k_{j,\text{ph}}\left(q\right)+{\displaystyle \sum _{\text{cations\hspace{0.05em}}M}{n}_M{k}_{j,\text{ex}}^M\left(q\right).}$$(32)

Here, n_M is the volume density of the molecule M that interacts in a charge exchange reaction with the dust grain Z.

Processes that decrease the charge, q + 1 → q, include the recombination with electrons and charge exchange reactions with molecular anions. The reverse rates were calculated as

$${\text{reverse}}_j\left(q\right)=n_e{k}_{j,\text{e}}\left(q\right)+{\displaystyle \sum _{\text{anions\hspace{0.05em}}M}{n}_M{k}_{j,\text{ex}}^M\left(q\right).}$$(33)

After determining the processes of the forward and reverse rates, the charge distribution function f_j(q) could be determined, separately in each size bin, in a step-down approach. We set f_j(q_max) = 1 and performed a sequence of steps where we calculated f_j(q) from f_j(q +1) until q = −q_max was reached. During each step, we used

$$f_j\left(q\right)=f_j\left(q+1\right)\frac{{\text{reverse}}_j\left(q+1\right)}{{\text{forward}}_j\left(q\right)}.$$(34)

Eventually, we calculated the normalization constant as

$$f_{j,\text{norm}}={\displaystyle \sum _{-q_{\max }}^{q_{\max }}{f}_j\left(q\right)}$$(35)

and normalized it as f_j(q) → f_j(q)/f_j,norm.

The results only depend on the UV to optical photon fluxes, the electron density ne, and the molecular ion densities M. In Fig. 2, we show the resulting f_j(q) for six dust size bins at two positions in the disk. The left side of the figure shows a point in the midplane at r = 0.1 au and z = 0, and in the right side, a point in the upper areas of the disk also with r = 0.1 au but with z = 0.29 au is shown.

In the midplane, all grains charge up negatively because there are no UV-photons here, and the rate coefficients for electron attachment are larger than the rate coefficients for charge exchange with molecular ions. In addition, grains charge up negatively because free electrons have a higher mean thermal speed than molecular ions, that is, $\sqrt{\frac{{\text{8k}}_{\text{b}}\text{T}}{\pi m_e}}>>\sqrt{\frac{{\text{8k}}_{\text{b}}\text{T}}{\pi m_M}}$. The lower part of the figure shows that the resulting equilibrium charge scales about linear with grain size because it is the electric potential q/a that enters into the rate coefficients. For example, all grains continue to charge up negatively until the rate for electron attachment, with decreases exponentially with q/a, balances the rates for charge exchange with molecular ions, which also depend on q/a. As a result, we get similar q/a values for all grain sizes. This behavior can also be seen and has been discussed in Okuzumi (2009). There, the authors showed that the shape of the charge distribution function is approximately Gaussian, with its width scaling with the square root of the grain size (compare to our Fig. 2). An overview of the behavior of the charge distribution function for the whole simulation, instead of examples at single points in the simulation (as seen in Fig. 2), can be seen in Fig. 3.

2.8 Effective rates in the chemical rate network

In order to solve the chemical network with the dust charge moments as species, as introduced in Sect. 2.6, we needed to derive the effective rate coefficients for the moments. We only demonstrate the derivations for the most relevant processes here. For a full description of all reactions that concern dust grains and their moment representation, we refer to Appendix C of Thi et al. (2019). For photoejection and photodetachment,

$$Z_{\text{m,}j}+hv\to Z_{\text{m,}j}^{+}+e^{-}$$(36)

$$Z_{\text{m},j}^{-}+hv\to Z_{\text{m,}j}+e^{-},$$(37)

we have

$$\begin{array}{c}\frac{d\left[Z_{\text{m,}j}^{+}\right]}{dt}=k_{\text{m,}j,\text{ph}}\left[Z_{\text{m,}j}\right]=n_{\text{d.}j}{\displaystyle \sum _0^{q_{\max }-1}{k}_{j,\text{ph}}\left(q\right)f_j\left(q\right)}\\ \frac{d\left[Z_{\text{m,}j}^{+}\right]}{dt}=k_{\text{m,}j,\text{ph}}^{-}\left[Z_{\text{m,}j}^{-}\right]=n_{\text{d.}j}{\displaystyle \sum _{-q_{\max }}^1{k}_{j,\text{ph}}\left(q\right)f_j\left(q\right)}.\end{array}$$

We note that these effective rates represent all photoreactions from charge state (q) to charge state (q + 1 ), not just for the neutral and negatively charged grains. To gain the effective rate coefficient, we divided by [Z_m,j] and $\left[{\text{Z}}_{\text{m,}j}^{-}\right]$, respectively:

$$k_{\text{m,}j,\text{ph}}=\frac{n_{\text{d,}j}}{\left[Z_{\text{m,}j}\right]}{\displaystyle \sum _0^{q_{\max }-1}{k}_{j,\text{ph}}\left(q\right)f_j\left(q\right)}$$(38)

$$k_{\text{m,}j,\text{ph}}^{-}=\frac{n_{\text{d,}j}}{\left[Z_{\text{m,}j}^{-}\right]}{\displaystyle \sum _{-q_{\max }}^{-1}{k}_{j,\text{ph}}\left(q\right)f_j\left(q\right).}$$(39)

In the code, once f_j(q) has been determined as explained in Sect. 2.7, we computed $\left[{\text{Z}}_{\text{m,}j}^{+}\right]$, $\left[{\text{Z}}_{\text{m,}j}\right]$, and $\left[{\text{Z}}_{\text{m,}j}^{-}\right]$ according to Eqs. (27) to (29) and then determined the effective photo rates from f_j(q) by Eqs. (38) and (39). The rate coefficients are hence fully determined by the discrete charge distribution functions f_j(q).

For the electron attachment, two new rate coefficients had to be derived. The derivation is analogous to the photo processes, but one has to consider the different charge states of the grains:

$$Z_{\text{m},j}+e^{-}\to Z_{\text{m},j}^{-}$$(40)

$$Z_{\text{m},j}^{+}+e^{-}\to Z_{\text{m},j}.$$(41)

As the positively charged grains enhance the electron attachment, we needed to describe their behavior with the rate coefficient $k_{\text{m,j,e}}^{+}$. Neutral grains do not show this kind of enhancement; therefore, we defined another rate coefficient k_m,j,e. We also included the contribution of negatively charged grains to the rate for the neutral grains, even though their contribution to the overall rate coefficient is low. For positive charged grains and neutral grains, with the added contribution of negatively charged grains, we could define

$$\begin{array}{l}\frac{d\left[Z_{\text{m},j}^{-}\right]}{dt}=k_{\text{m},j,\text{e}}\left[Z_{\text{m},j}\right]=n_{\text{d},j}\left(k_{j,\text{e}}\left(0\right)f_i\left(0\right)+{\displaystyle \sum _{-q_{\max }+1}^{-1}{k}_{j,\text{e}}\left(q\right)f_i\left(q\right)}\right)\hfill \\ \frac{d\left[Z_{\text{m},j}\right]}{dt}=k_{\text{m},j,\text{e}}^{+}\left[Z_{\text{m},j}^{+}\right]=n_{\text{d},j}{\displaystyle \sum _1^{q_{\max }}{k}_{j,\text{e}}\left(q\right)f_i\left(q\right)}.\hfill \end{array}$$

Here, k_j,_e(q) are the rate coefficients for the individual charge states defined by Eq. (16). The effective rate coefficients for electron attachment are hence

$$k_{m,j,\text{e}}=\frac{n_{\text{d},j}}{\left[Z_{m,j}\right]}\left(k_{j,\text{e}}\left(0\right)f_j\left(0\right)+{\displaystyle \sum _{-q_{\max }+1}^{-1}{k}_{j,e}\left(q\right)f_i\left(q\right)}\right)$$(42)

$$k_{m,j,\text{e}}^{+}=\frac{n_{\text{d},j}}{\left[Z_{m,j}\right]}{\displaystyle \sum _1^{q_{\max }}{k}_{j,e}\left(q\right)f_i\left(q\right)}.$$(43)

The charge exchange reactions between dust grains and molecules also had to be adjusted. We considered the rates for neutral grains and negative grains separately:

$$Z_{\text{m},j}+M^{+}\to Z_{\text{m},j}^{+}+M$$(44)

$$Z_{\text{m},j}^{-}+M^{+}\to Z_{\text{m},j}+M.$$(45)

This resulted in the following charge exchange rates:

$$\frac{d\left[Z_{\text{m},j}^{+}\right]}{dt}=k_{\text{m},j,\text{ex}}^M\left[Z_{\text{m},j}\right]\left[M^{+}\right]=n_{\text{d},j}\left[M^{+}\right]{\displaystyle \sum _0^{q_{\max }-1}{k}_{j,\text{ex}}^M\left(q\right)f_i\left(q\right)}$$(46)

$$\frac{d\left[Z_{\text{m},j}\right]}{dt}=k_{\text{m},j,\text{ex}}^{M,-}\left[Z_{\text{m},j}^{-}\right]\left[M^{+}\right]=n_{\text{d},j}\left[M^{+}\right]{\displaystyle \sum _{-q_{\max }}^{-1}{k}_{j,\text{ex}}^M\left(q\right)f_i\left(q\right)},$$(47)

from which we derived the rate coefficients for charge exchange with molecule M⁺:

$$k_{\text{m,}j,\text{ex}}^M=\frac{n_{\text{d,}j}}{\left[Z_{\text{m,}j}\right]}{\displaystyle \sum _0^{q_{\max }-1}{k}_{j,\text{ex}}^M\left(q\right)f_i\left(q\right)}$$(48)

$$k_{\text{m,}j,\text{ex}}^{M{,}_{-}}=\frac{n_{\text{d,}j}}{\left[Z_{\text{m,}j}^{-}\right]}{\displaystyle \sum _{-q_{\max }}^{-1}{k}_{j,\text{ex}}^M\left(q\right)f_i\left(q\right)}.$$(49)

With these effective rate coefficients for the dust charge moments, we could either solve our rate network system for the time-independent solution, or we could advance the ordinary system of first-order differential equations (ODE system) in time from an initial vector of particle densities. Appendix C explains how to deal with the problem that arises because f_j(q) depends on the electron and molecular ion densities and the rate coefficients depend on particle densities.

2.9 Cosmic ray implementation

As explained in Sect. 3, we expected lightning to occur, if at all, in the dense and shielded midplane regions of the disk close to the star, where cosmic rays are the most relevant source for ionization. It was therefore important to use a model for the cosmic ray penetration that is as realistic as possible.

We used the description of cosmic ray attenuation and ionization following Padovani et al. (2009) and as implemented by Rab et al. (2017, Rab et al. 2018). According to this model, the H₂ cosmic ray ionization rate ζ_cr s⁻¹ is given by

$${\zeta }_{\text{cr}}=\{\begin{array}{cc}{\zeta }_{\text{low}}& ,N_{\langle \text{H}\rangle }<{10}^{19}{\text{cm}}^{-2}\\ \frac{{\zeta }_{\text{low}}{\zeta }_{\text{high}}}{{\zeta }_{\text{high}}{\left(\frac{N_{\langle H\rangle }}{{10}^{20}{\text{cm}}^{-2}}\right)}^{\text{a}}{\zeta }_{\text{low}}\left(\exp \frac{\sum }{{\sum }_0}-1\right)}& \text{,otherwise}\end{array},$$(50)

where N_〈H〉 [cm⁻²] is the vertical hydrogen nuclei column density and Σ ≈ (1.4 amu) N_〈H〉 [g cm⁻²] the vertical mass column density. The fitting parameters are represented by ζ_low, ζ_high, Σ₀, and a. These fitting parameters originate from the fitting of cosmic ray spectra and account for different attenuations at different column densities. The term ζ_1ow is used to fit the spectra for low column densities smaller than 10¹⁹ cm⁻², where the attenuation can be described as a power law, and ζ_high is used for column densities higher than 10¹⁹ cm⁻², where the attenuation can be described with an exponential term. Additionally, Eq. (50) is applied twice in our models (once to account for cosmic rays impacting the disk from above and once to account for cosmic rays impacting the disk from below) by adjusting N_〈H〉 = 2N_〈H〉(r, 0) − N_〈H〉(r, z), where N_〈H〉(r) is the total vertical hydrogen nuclei column density at point r and N_〈H〉 (r, z) is the local hydrogen nuclei column density at the point (r, z) for which the cosmic ray ionization rate has to be determined. In our main simulation, we made use of the cosmic ray spectra named Solar Max in Cleeves et al. (2013). We do not include stellar energetic particles (SEPs) in this work.

3 Resulting dust charge and ionization properties

This section summarizes our results concerning the mean charge of the dust grains, how the grain charge distribution function depends on the position in the disk and on dust size, and how the dust charge is linked to the degree of ionization and molecular ions in the gas, with particular emphasis on the midplane regions. We first illustrate our standard case, where we only generated exchange reactions but no new protonation reactions and only considered exothermic reactions. A full discussion can be found in Appendix B. These results are used later to discuss whether charge separation and lightning can occur in disks (see Sect. 4). An overview of the hydrogen nuclei density and resulting radial and vertical extinction can be seen in Fig. 4. Additionally, an overview of the temperature structure can be seen in Fig. 5.

3.1 Electron concentration

Figure 6 shows the resulting electron concentration n_e/n_〈H〉 as a function of position (r, z) in our disk model. On the basis of these results, we introduced six different disk regions, A-F (see Fig. 6), where the physical and chemical processes leading to ionization and dust charge are qualitatively different in each case.

3.2 Regions F, E, and D – the plasma regions

Regions D–F in Fig. 6 are radiation dominated regions, where the UV and X-ray photons create a large degree of ionization along with positive dust charges. We found that region F is dominated by the ionization of H and H₂, region E is dominated by the ionization of atomic carbon, and region D by the ionization of atomic sulfur. According to the assumed element abundances of H, C, and S in our disk model, the degree of ionization in regions F, E, and D is about 1, 10⁻⁴, and 10⁻⁷, respectively. The grains charge up positively via photo-effect, Z + hv → Z⁺ + e⁻, and their charge is balanced by electron recombination, Z⁺ + e⁻ → Z. The total number of charges on the grains, however, is insignificant in comparison to the number of free electrons and positive ions in the gas phase. Any dynamical displacement of the charged grains would easily be balanced out by slight motions of the free electrons in the gas, making charge separations very unlikely to occur in these plasma regions (see further discussion in Sect. 4).

3.3 Region A – the dust dominated region

We define region A as the disk region where the number of negative charges on the dust grains is balanced by the abundance of molecular cations. Free electrons are very rare in this region (degree of ionization 10⁻¹⁸… 10⁻¹⁴) and unimportant for the charge balance. In our disk model, region A is the area that stretches out from just behind the inner rim to about r ≈ 1 au and z/r ≲ 0.1. The outer boundary coincides with the location where water and ammonia ice emerge (snowline), and the upper boundary is roughly given by a vertical visual extinction of $A_V^{\text{ver}}=10$, which makes sure that UV photons cannot penetrate into region A.

Figure 7 shows that we found ${\text{NH}}_4^{+}$ to be by far the most important molecular cation in region A. The chemical processes that lead to this kind of charge balance in region A are multi-staged and illustrated in Fig. 8.

Since region A is entirely shielded from UV photons and X-rays, cosmic rays are found to be the only relevant ionization source, in particular

$${\text{H}}_2+\text{CR}\to {\text{H}}_2^{+}+{\text{e}}^{-}.$$(51)

While the free electrons are quickly picked up by the dust grains via Z + e⁻ → Z⁻, the H₂⁺ molecules, after fast reactions with H₂, form the molecular cation H₃⁺, which has a relatively low proton affinity (see Table 2). Therefore, the surplus protons in H₃⁺ are quickly passed on to other abundant molecules, creating more complex molecular cations, such as H₃O⁺, HCNFT, and NH₄⁺. Further proton exchange reactions with abundant neutrals tend to increase the abundances of the protonated molecules, which have the highest proton affinities.

The NH₄⁺ molecule has an extremely high proton affinity of 8.9 eV, which means that it cannot recombine and dissociate on grain surfaces, that is, the reaction ${\text{Z}}^{-}+{\text{NH}}_4^{+}\to {\text{Z+NH}}_3+$ Η is energetically forbidden. This creates a dead end in which the concentration of ${\text{NH}}_4^{+}$ is enhanced until an equilibrium is established with the dissociative recombination reaction ${\text{NH}}_4^{+}+{\text{e}}^{-}\to {\text{NH}}_3+\text{H}$ with the extremely rare free electrons.

However, the other protonated molecules below the threshold shown in Table 2, denoted by MH⁺ in Fig. 8, can dissociatively recombine on the surface of the negatively charged dust grains, which creates a charging balance for the dust grains via the following three reactions:

$$\text{Z}+{\text{e}}^{-}\to {\text{Z}}^{-},$$(52)

$${\text{Z}}^{-}+{\text{MH}}^{+}\to \text{Z}+\text{M}+\text{H},$$(53)

$${\text{Z}}^{-}+{\text{A}}^{+}\to \text{Z}+\text{A},$$(54)

where M is the neutral abundant molecules and MH⁺ is their protonated counterparts, in particular H₃O⁺ and HCNH⁺ in region A. There are also some simple charge exchange reactions with atomic ions A⁺, such as Mg⁺ and Fe⁺, whereas the electron affinity of Na⁺ of 5.14 eV is too small to detach an electron from a silicate grain on impact (work function 8 eV). We therefore eliminated the reaction ${\text{Z}}^{-}+{\text{Na}}^{+}\to \text{Z}+\text{Na}$ Na from our databases. None of the chemical reaction rates visualized in Fig. 8 involve activation barriers, which is important to remember when using this reaction scheme to discuss how the charge balance between Ζ⁻, MH⁺, NH₄⁺, and e⁻ reacts to any changes in density, CRI rate, or temperature. However, reaction (52) depends on the number of negative charges already collected by the grains, which is hence a result of the model.

Figure 7 shows that NH₄⁺ remains the most abundant molecular cation up to a radial distance of about 0.7 au, which coincides with the formation of ammonia ice in the midplane of our disk model. At this point, NH₃ is no longer an abundant molecule, and protonated silicon monoxide SiOH⁺ takes over the role of proton keeper from NH₄⁺. The proton affinity of SiOH⁺ is 8.1 eV, which means that it also cannot recombine dissociatively on negatively charged dust grain surfaces. Eventually, SiO freezes out as well, at 0.9 au, and we transit into region B.

3.4 Region Β – the intermediate region

We identify region Β as the disk region where the concentrations of (1) negative charges on dust grains, (2) free electrons, and (3) molecular cations are all of the same order, that is, ≈ 10⁻¹³. This happens in a transition region between about 1 au and 3 au in our disk model (see Fig. 9). Again, we required $A_V^{\text{ver}}\lesssim 10$, ruling out photoionization, which corresponds to z/r ≲ 0.1 in our disk model.

Figure 10 shows the stepwise freeze-out of oxygen, nitrogen, and carbon into ices with an outward falling temperature, which happens mostly in region B. The freeze-out starts with water ice, then ammonia ice around 1 au, and eventually several hydrocarbon ice phases, such as C₃H₂#, C₂Hs#, and C₂H₄# farther out. At the end of this ice formation zone, at a distance of about 3 au, the midplane is virtually devoid of any molecules other than H₂, He, noble gases, and some sulfur molecules, such as H₂S, CS, and H₂CS.

Therefore, the transition region Β is carbon rich, with oxygen and nitrogen already being strongly depleted by ice formation. Among the remaining hydro-carbon molecules in region B, we found Cyclopropenylidene C₃H₂ to have the highest proton affinity (see Table 2). Consequently, the reaction diagram in Fig. 8 would need two slight modifications for to be valid for Region Β as well. Firstly, one needs to replace NH₄⁺ by C₃H₃⁺, because C₃H₃⁺ takes over its role as proton keeper. Secondly, we would replace M with H₂S,CS and H₂CS.

3.5 Region C – the metal-poor region

In region C (Fig. 11), there are almost no molecules other than H₂ left in the gas phase because nearly all the oxygen, carbon, nitrogen, and sulfur is in the ice. Since the various electron recombination rates scale with n² but the cosmic ray ionization rate scales with n, the electron density steadily increases toward larger radii, whereas the concentration of negative charges on dust grains Ζ⁻ stays about constant. In fact, we used n_e = [Z⁻] to set the boundary between regions Β and C, which happens at about 3 au in our model.

Consequently, we observed a charge equilibrium between free electrons and H₃⁺ in region C, with some traces of F⁺, whereas the negatively charged dust grains became increasingly unimportant for the charge equilibrium. Toward the outer edge of region C, the disk becomes vertically transparent again, some interstellar UV and scattered UV starlight reaches the midplane, and we found increasing amounts of ionized sulfur in the gas phase. At a radial distance of about 100 au in our model, F⁺ and S⁺ become more abundant than H₃⁺, and we enter region D.

3.6 Dependence of grain charge on physical parameters

In order to investigate the charging behavior of dust grains in more detail, we wanted to investigate the trends seen in Fig. 2 further and see how they react to the changing of parameters that influence dust, gas, and electron abundance. In order to do this, we changed five specific parameters and studied how the affected species reacted to these changes. We changed the cosmic ray intensity ζ, the minimum dust grain size a_min, the total gas density ρ, the dust-to-gas ratio, and the dust temperature T_d. We did this by changing the parameters compared to our standard case (Table 3) and solved the chemistry only at the point in our simulation with r = 0.1 au and z/r = 0 in the midplane. An overview of the parameter space for each changed parameter can be found in Table 4. We only ever changed one parameter at a given time, and the other unchanged parameters remained equivalent to our standard case.

The main results of this section can be found in Fig. 12. Through the plots in the figure, we investigate the abundances of the three species that are the most influential to the charge balance at the investigated point and how the charge distribution function, as seen in the lower parts of Fig. 2, compared to q_j/a, changes when varying the different parameters. In order to see changes in the charge distribution functions compared to q_j/a, we created these plots for every dust bin of the different simulations. In order to make a comparison feasible, we took an average of the peaks of these charge distribution functions since Fig. 2 already revealed a very uniform behavior between the different dust bins. In addition, we calculated the standard deviations of the different q_j/a values σ_j. These standard deviations were averaged as well, and this average can be seen in Fig. 12 as the gray shaded areas. In the following paragraphs, we describe the different subplots of Fig. 12 for each parameter in more detail.

Minimal dust size. For testing the influence of the minimum dust size on the dust charging behavior, we tested five different minimum dust sizes, starting with a_min = 0.05 μm, which is the default value of our standard simulation, and increasing by a factor of ten at each step until we reached a_min = 500 μm. The first big trend we observed was the decrease of the abundance of the negatively charged dust grains. This is due to the smaller surface area that a few large dust grains have, compared to many smaller ones, since one of the main contributors to the abundance of negatively charged dust grains is the total surface area. This resulted in lower rate coefficients and therefore less negatively charged dust grains, as they react less with electrons. This also explains the two other trends we observed: an increase in electron abundance and a decrease in ${\text{NH}}_4^{+}$ abundance until a saturation is reached at a_min = 50 μm. As explained before, this is due to the lower amount of dust grains. Dust grains are the main reaction partner of electrons, resulting in the normally high amount of negatively charged dust grains. If fewer dust grains react with electrons, one gets more free electrons.

These free electrons are also part of the main destruction path for the ${\text{NH}}_4^{+}$, which explains the decrease of the ${\text{NH}}_4^{+}$ molecules. Both the electrons and the ${\text{NH}}_4^{+}$ abundances converge toward each other because the dust grain abundance becomes fully negligible with a sufficiently high minimum dust grain size.

Dust-to-gas ratio. We tested the influence of different dust-to-gas ratios by simulating five different dust-to-gas ratios, starting with 10⁻⁴ and again increasing at each step by a factor of ten up to a dust-to-gas ratio of one. The clear trend that can be seen is the decrease of the electron abundance with increasing dust-to-gas ratios. To understand these results, it is central to remember the role dust grains play in the midplane chemistry. Dust grains are the main reaction partner of free electrons. Therefore, increasing the amount of dust grains by increasing the dust-to-gas ratio decreased the amount of electrons. Secondly, as mentioned in the last paragraph, the main destruction pathway of ${\text{NH}}_4^{+}$ molecules is via reaction with free electrons. Hence, a reduction of electrons via having more dust grains resulted in an increase of the ${\text{NH}}_4^{+}$ molecules as well.

Cosmic rays. To investigate the impact different cosmic ray ionization rates have on our results, we chose to simulate six different cosmic ray ionization rates from ζ = 10⁻¹⁶ to ζ = 10⁻²². We note that we assumed a constant cosmic ray ionization rate, in contrast to the method we chose for our large simulations. This was done to make it easier to compare results. At the point we investigated, the cosmic ray ionization rate for our standard case is 5 × 10⁻²⁰.

The strongest trend one can see is the strong relation between the abundance of the electrons and the cosmic ray ionization rate. The trend is that the lower the cosmic ray ionization rate is, the lower the amount of electrons. This is caused by the fact that the main source of electrons in these highly shielded regions is the ionization of molecular hydrogen via cosmic rays.

We also saw no change in the abundance of ${\text{NH}}_4^{+}$ and negative charged dust grains. This is the case because the cosmic ray ionization rate also influences how many ${\text{H}}_3^{+}$ molecules are available. A decrease in cosmic ray ionization rate also results in less molecular hydrogen ionization. This results in a constant abundance of both negative dust grains and ${\text{NH}}_4^{+}$ since there is a decrease of the species that are influential for both the creation and destruction of negatively charged dust grains and ${\text{NH}}_4^{+}$. For negatively charged dust grains, the important creation and destruction species are electrons for the creation and protonated molecules for the destruction. For ${\text{NH}}_4^{+}$ it is the opposite.

Gas density. To see if changing the gas density changes our result, we simulated different models with a differing disk mass. As the amount of gas in our model cannot be directly increased or decreased, we had to change the disk mass instead, which indirectly increases or decreases the gas density. We simulated five different disk masses, varying from 10⁻⁴ M_⊙ up to 1 M_⊙. We note, however, that we plotted against gas density and not disk mass. The results are very similar to what we got from the runs where we varied the cosmic ray ionization rate, and this is because increasing the gas density also increases the shielding and hence decreases the amount of cosmic ray ionization and therefore free electrons. Thus, the effects mentioned in the previous section mostly hold true for this part as well.

Dust temperature. We further wanted to investigate how changing the dust temperature could influence the charging behavior of the grains. Thus, we changed the dust temperature from 100 K to 1800 K in steps of 100 K. We note that we also changed the gas temperature at the investigated spot to be equal to the dust temperature at the beginning. The findings can be summarized into four different parts.

At the very first point at 100 K, we observed one clear deviation from the standard picture in the concentration of the ${\text{NH}}_4^{+}$ molecules. The molecules are less abundant by two orders of magnitude compared to the normal case. This is due to the fact that at these lower temperatures, we see a similar case to region B in our larger simulations, where nitrogen bearing species are frozen out, therefore reducing nitrogen abundances in the gas. Also similar to region B is that ${\text{C}}_3{\text{H}}_3^{+}$ is now the most abundant positive species, as carbon has not been frozen out yet and is therefore abundant in the gas phase.

For temperatures between 200 K to 1300 K, the concentrations of the different species resemble the standard case but with two exceptions. First, we found that the ability of the dust grains to gather electrons scales linearly with their temperature. Second, due to the dust grains being able to carry more negative charges, we found a slight increase in the concentrations of the negatively charged dust grains and a decrease in the electron concentration. This also resulted in an increase in the concentration of the ${\text{NH}}_4^{+}$ molecules.

From 1400 K to 1800 K, we observed Na⁺ become a dominant species. This is due to the creation reactions of Na⁺ becoming more efficient. The main creation reaction that becomes more efficient is the collisional ionization of Na via H₂ :

$$\text{Na}+{\text{H}}_2\to {\text{Na}}^{+}+{\text{H}}_2+{\text{e}}^{-}.$$(55)

This reaction is endothermic in nature, but the barrier of dH_f / k_b ≈ 60000 is still low enough to be overcome more often at temperatures greater than 1400 K. This results in Na⁺ becoming the most abundant positive species after 1400 K, and after 1600 K, it has a similar abundance as the dust grains. Additionally, we observed an increase in the electron concentration, as one of the products of reaction in Eq. (55) are electrons. This increase of electrons results in an overall reduction of the ${\text{NH}}_4^{+}$ molecules, as the main destruction reaction for ${\text{NH}}_4^{+}$ is the recombination with electrons.

4 A simple model for turbulent charge separation

4.1 Turbulence induced electric fields

In this section, we develop a simple model to estimate the maximum electric field E that can arise when a mixture of negatively charged grains and molecular cations is shaken by turbulence. Figure 13 shows a physical sketch of the situation. We used the Richardson-Kolmogorov approach for turbulence (Richardson 1926; Kolmogorov 1941). Therefore, we assume a superposition of turbulent eddies k with various spatial scales ℓ_k, timescales τ_k, and characteristic velocities v_k = ℓ_k/τ_k. In the inertial range

$$\begin{array}{cc}{\upsilon }_k\propto {(\dot{ϵ}{\ell }_k)}^{1/3}& ({\ell }_{\eta }<{\ell }_k<{\ell }_L)\end{array},$$(56)

is valid, where $\dot{ϵ}$ is the energy dissipation rate toward smaller and smaller turbulent scales and ℓ_η and ℓ_L are the smallest (thermalization) and largest (driving) turbulent length scales, respectively. In the co-moving frame of a curved gas flow related to the eddy k, the dust grains are accelerated by the centrifugal force $m{\upsilon }_k^2\text{}/{\ell }_k$ until an equilibrium with the frictional force F_fric in the gas is established. This problem is well known for constant gravity g (see, e.g., Woitke & Helling 2003.) After an initial acceleration timescale τ_acc, the grains of mass m and size a reach a constant drift velocity ${\stackrel{°}{\upsilon }}_{\text{dr}}$ also known as the final fall speed. We took the results from Woitke and Helling for the case of a subsonic flow and large Knudsen numbers, also known as the Epstein regime. After replacing g by ${\upsilon }_k^2\text{}/{\ell }_k$ the results read:

$$\begin{array}{lll}{\tau }_{\text{acc}}=\frac{a{\rho }_{\text{m}}}{\rho {\upsilon }_{th}},\hfill & F_{\text{fric}}=\frac{m{\upsilon }_{\text{dr}}}{{\tau }_{\text{acc}}},\hfill & {\stackrel{°}{\upsilon }}_{\text{dr}}=\frac{{\upsilon }_k^2}{{\ell }_k}{\tau }_{\text{acc}},\hfill \end{array}$$(57)

where τ_acc is the acceleration (or stopping) time, ρ_m is the dust material density ≈ 2 g cm⁻³, and ${\upsilon }_{\text{th}}=\sqrt{8kT\text{}/\left(\pi \overline{m}\right)}$ is the thermal velocity, with $\overline{m}$ being the mean molecular weight of the gas particles.

The drift of the charged grains leads to a charge separation, which causes an electric field. How exactly this field will look is complicated and depends on the geometry and the superposition of the electric currents caused by the drifting grains in the various turbulent eddies. However, we can estimate the maximum electric field that can be generated by a single long-lived eddy by considering the case $\left({\tau }_E,{\tau }_{\text{acc}}\right)<\text{t}<t_k$ where τ_Ε is the electric field buildup timescale and t_k is the eddy turnover time. We assumed the electric field builds up in a coherent manner, from larger to smaller eddies, but as shown later by Eqs. (61) and (71), the main contribution comes from smaller eddies.

In this case, after the dust grains become accelerated by the eddies (Fig. 13, panel 2: acceleration), an electric field will continue to build up (Fig. 13, panel 3: intermediate) until the molecular cations in the gas follow the drifting grains because of their mobility in the E field created by the grains. We call this equilibration (Fig. 13, panel 4: equilibration), and it is characterized by a vanishing electric current:

$$j_{\text{e}1}={\displaystyle \sum _j([z_{\text{m},j}^{+}]-[z_{\text{m},j}^{-}])}{\stackrel{°}{\upsilon }}_{\text{dr},j}-n_{\text{e}}{\upsilon }_{\text{e}}+n_{\text{I}}{\upsilon }_{\text{I}}=0.$$(58)

As before, j is an index for the dust size bins; $\left[Z_{m,j}^{+}\right]$, $\left[Z_{m,j}^{-}\right]$ are the concentrations of the dust charge moments [charges cm⁻³] (see Eqs. (27) and (28)); and n_e and n_I are the electron and molecular cation densities [charges cm⁻³].

The drift velocities of the molecular cations and electrons are given by

$${\upsilon }_{\text{I}}=-{\mu }_{\text{I}}E$$(59)

$${\upsilon }_{\text{e}}={\mu }_{\text{e}}E,$$(60)

where μ_I and μ_e are the mobilities of the molecular cations and free electrons, respectively, in units of cm² s⁻¹ V⁻¹. The E field is measured in V cm⁻¹. The different signs in Eqs. (59) and (60) are because of the opposite directions of the drift velocities of molecular cations and electrons in a given electric field.

If the turbulent eddy under consideration lives long enough for the drifting grains to build up the maximum, "equilibrated" electric field, it follows from Eq. (58) that

$$E=\frac{{\displaystyle {\sum }_j\left(\left[Z_{\text{m,}j}^{+}\right]-\left[Z_{\text{m,}j}^{-}\right]\right)}{\stackrel{°}{\upsilon }}_{{}_{\text{dr},j}}}{n_{\text{I}}{\mu }_{\text{I}}+n_{\text{e}}{\mu }_{\text{e}}}.$$(61)

From Eq. (61) we can immediately see which turbulent flows can potentially produce large E fields. Normal plasmas have n_e ≈ n_I » [Z⁺], [Z⁻], in which case we do not expect large turbulence-induced electric fields. Dusty plasmas, like the one we expect in region A in the disk, have [Z⁻] ≈ n_I » n_e » [Z⁺]. The induced electric fields in this case may be large when ${\stackrel{°}{\upsilon }}_{\text{dr}}$ is large, μ_Ι is small, and n_e is sufficiently small to play no role, in which case n_I and $\left[Z^{-}\right]={\displaystyle {\sum }_j\left[Z_{\text{m,}j}^{-}\right]}$ cancel, and we find

$$E=\frac{\langle {\stackrel{°}{\upsilon }}_{\text{dr}}\rangle }{{\mu }_{\text{I}}},$$(62)

where $\langle {\stackrel{°}{\upsilon }}_{\text{dr}}\rangle ={\displaystyle {\sum }_j\left[Z_{\text{m,}j}^{-}\right]{\upsilon }_{\text{dr,}j}/{\displaystyle {\sum }_j\left[Z_{\text{m,}j}^{-}\right]}}$ is the charge-mean grain drift velocity. There is another case when [Z⁻] ≈ [Z⁺] » n_I, n_e. In that case, which we call “region X” (not found in our present disk models), large E fields might be produced when many small grains are strongly charged negatively and many large grains are strongly charged positively, or vice versa, such that the nominator in Eq. (61) remains large despite charge neutrality in an isolating gas with few charged gas particles. We note that tribolelectric charging, which is not included in our current disk models, might be able to turn parts of region A into a region X.

We mention that at this point, there are alternative approaches to how we consider ion and electron mobilities. In Okuzumi & Inutsuka (2015), they present analytic expressions for ion and electron drift velocities in an H₂ gas.

In order to estimate the duration of the equilibration phase in Fig. 13, henceforth called the E field buildup timescale τ_Ε, we considered the travel time of the grains

$${\tau }_{\text{E}}=\frac{d}{{\stackrel{°}{\upsilon }}_{\text{dr}}},$$(63)

where d is the charge separation distance, that is, the distance between the moving grains and the moving molecular ions causing the electric field E. To estimate that distance, we considered the most simple geometry of a parallel plate capacitor (with equations in SI units):

$$E_{\text{SI}}=\frac{Q_{\text{SI}}}{{ϵ}_{0}A}$$(64)

$$Q_{\text{SI}}=e{\left[Z^{-}\right]}_{\text{SI}}A{d}_{\text{SI}}$$(65)

$$\Rightarrow d_{\text{SI}}=\frac{{ϵ}_0{E}_{\text{SI}}}{e{[Z^{-}]}_{\text{SI}}},$$(66)

where ϵ₀ is the electric vacuum permittivity and e is the electron charge in SI units; E_SI = 10² E is the maximum electric field in [Vm⁻¹]; A is a horizontal area [m²] that cancels; [Z⁻]_SI = 10⁶ [Z⁻] is the concentration of negative charges on dust grains in $\left[{\text{m}}^{-3}\right];d_{\text{SI}}={10}^{-2}$ d is the charge separation distance in [m]; and Q_SI is the total negative charge on the grains in volume d_SI × A in [C]. Equations (57), (63), and (66) provide an estimate for the time required to clear a volume d × A from dust grains of size a by centrifugal forces in a turbulent eddy k, which results in an overpopulation of negative charges at the lower boundary and missing negative charges at the upper boundary, which together build up the maximum electric field.

Thus, we have three different timescales to consider that may depend on the considered particle size a and the selected turbulent eddy k. In most cases, we found

$${\tau }_k\gg {\tau }_{\text{acc}}(a)\gg {\tau }_{\text{E}}(a,k)$$(67)

to be valid, that is, the turbulent eddy lives long enough to allow the grains to accelerate quickly and then drift for long enough distances to cause the maximum electric field. If relation (67) does not hold, Eq. (61) is not valid.

For the ion mobility μ_I we considered the NH₄⁺ cation, as we found this molecule to be the most abundant gas charge carrier in region A. Abedi et al. (2014) presented measurements for the mobility of ${\text{NH}}_4{}^{+}$ in a N₂ gas at standard pressure and temperature, P₀ = 760 Torr and T₀ = 273 K, which we denote by ${\mu }_{{\text{NH}}_4^{+},{\text{N}}_2}^0=2.22{\text{cm}}^2{\text{s}}^{-1}{\text{V}}^{-1}$ As the most abundant species in the disks is H₂, we needed to adjust this value. This was done by multiplying with the ratio of the reduced masses of N₂ and H₂ with NH₄⁺, denoted with ${\overline{m}}_{{\text{N}}_2}$ and ${\overline{m}}_{{\text{H}}_2}$, respectively, and scaling with the gas particle density (see Eq. (4) in Abedi et al. 2014):

$${\mu }_{\text{I}}={\mu }_{{\text{NH}}_4^{+},{\text{N}}_2}^0{\left(\frac{{\overline{m}}_{{\text{N}}_2}}{{\overline{m}}_{{\text{H}}_2}}\right)}^{1/2}\frac{P_0}{P}\frac{T}{T_0}.$$(68)

The electron mobility μ_e was calculated from the drift velocity measurements for electrons of Ryzko (1965) and calculated as

$${\mu }_{e,0}=\frac{v_{\text{dr},\text{e}}}{E},$$(69)

with υ_dr,e measured in dry air to be

$${\upsilon }_{\text{dr},\text{e}}=256\times {10}^3\frac{E}{E_0}\frac{P_0}{P}\frac{T}{T_0}[\text{cm\hspace{0.17em}}{\text{s}}^{-1}],$$(70)

where P₀ and T₀ are the same as for the ion mobility and E₀ = 75 Vcm⁻¹. This formula results in an electron mobility that is about a factor of 200 larger than the mobility of NH₄ ⁺. As the original measurements were done in dry air, we would have to adjust our mobility calculations in the same way we did for ion mobility, but due to only minor differences in reduced masses when considering the electron mass, we omitted making a similar correction.

According to our charge separation model, the electric fields that can potentially be caused by turbulence are proportional to the centrifugal accelerations ${\upsilon }_k^2/{\ell }_k$ that the turbulent eddies can provide. Following Eq. (56), we found

$$\frac{{\upsilon }_k^2}{l_k}\propto l_k^{-1/3},$$(71)

that is, the smallest eddies cause the largest centrifugal forces. We therefore considered the Kolmogorov timescale t_η, that is, the turnover timescale of the smallest eddy in the inertial subrange. Following Ormel & Cuzzi (2007), we determined τ_η by

$${\tau }_{\eta }=\frac{{\tau }_L}{\sqrt{ℛe}},$$(72)

where τ_L is the turnover timescale of the largest eddy in our system and ℛe is the Reynolds number defined by

$$ℛe=\frac{{\upsilon }_L{\ell }_L}{v_{\text{kin}}},$$(73)

where υ_L is the velocity of the largest eddy, ℓ_L is its size, and v_kin is the kinematic viscosity. The kinematic viscosity is from Woitke & Helling (2003). Using their formula for the dynamic viscosity, we got ${\upsilon }_{\text{kin}}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.{\upsilon }_{\text{th}}\overline{\ell }$. For the size and velocity of the largest eddy, we assumed ${\ell }_L=\sqrt{\alpha }{H}_p$ and ${\upsilon }_L=\sqrt{\alpha }{c}_s$ (see Ormel et al. 2007), where H_p is the pressure scale height, which is set in our simulations as

$$H_p(r)=10\text{\hspace{0.17em}au\hspace{0.17em}}{\left(\frac{r}{100\text{\hspace{0.17em}au}}\right)}^{1.15}.$$(74)

The resulting size of the smallest turbulent eddy ℓ_η is considerably larger than the gas mean free path $\overline{\ell }$ by about four orders of magnitude, but it is also much smaller in scale height by about eight orders of magnitude.

Since we wanted to maximize the effect of the turbulence, we had to choose the τ_k that maximizes the centrifugal acceleration. Additionally we had to choose a time that allows the electrical fields to develop. Therefore, we followed the principle that we choose the shortest time that is still long enough for the fields to develop. This is illustrated by Fig. 14. Following this principle of smallest time that still lets the fields develop, we can see that in most cases we can choose τ_k = τ_η. In cases with larger grains, a > 10⁻¹cm, we would choose τ_k = τ_acc to ensure the grains have enough time to reach their drift velocity ${\stackrel{\circ }{\upsilon }}_{\text{dr}}$.

With these results, we could then calculate what electric fields would result from the interaction of our dust and gas mixture in a turbulent eddy and answer the question of whether lightning could emerge with this mixture already. Nonetheless, we still had to consider how large these electric fields could be. For lightning to emerge in any medium, the electric fields generated have to be large enough for an electron cascade to occur. We followed the description and Eq. (56) from Muranushi & Tomiyasu (2009). With this we were able to define a critical electric field, E_crit, that is the electric field at which an electron cascade can occur because the electrons are accelerated enough to ionize the neutral components of the gas:

$$E_{\text{crit}}=\frac{\text{Δ}W}{e\lambda }.$$(75)

Here, ∆W is the ionization energy required to ionize the neutral components, e is the coulomb constant, and λ is the mean free path. For ∆W, we took the ionization energy of H₂ with 15.4 eV. To determine, λ we followed Eq. (10) of Woitke & Helling (2003).

4.2 Electric field magnitudes

In this section, we discuss our investigation of the magnitude of the previously discussed turbulence-induced electric fields. We took the data from our simulations that we made for investigating the dust charge behavior from Sect. 3.6 to perform this investigation. The resulting electric fields are shown in Fig. 15.

Minimal dust size. We found that if one changes the minimal dust size, one first sees an increase in electric field strength due to the increase in dust size and in drift velocity and therefore in electric field strength as well. However, with a minimum dust grain size of more than 0.5μm, a decrease in electric field strength occurs. This is due to the fact that larger dust grains are less abundant in our model. Regarding the way we set up the simulations, only very large dust grains were considered, and therefore the dust abundance drops. According to Eq. (61), this of course has a negative impact on the total electric field strength.

Dust-to-gas ratio. For increasing dust-to-gas ratios, we observed a steady increase in electric field strength until a saturation took place at our standard value of 10⁻². The reason for the low electric field values at lower dust-to-gas ratios is again the fact that fewer dust grains can contribute positively toward the creation of electric fields.

Cosmic rays. The influence of cosmic rays on the electric fields is not very surprising given their strong link to electron abundance. High cosmic ray ionization rates result in high electron abundances and therefore low electric fields, due to the attenuating effect electrons have on the buildup of electric fields. The reverse is true for low cosmic ray ionization rates, however, and this effect stagnates at around 10⁻¹⁹[s⁻¹].

Gas density. One could expect that the gas density has a similar effect on the electric fields as the cosmic ray ionization rate, but this does not fully hold. While the general trend is the same, that is, where high gas densities reduce the electron abundance due to the larger shielding, there are more effects that must be mentioned. Compared to the other tests where the critical field is kept constant, we found that the critical field rises with increasing gas densities. This is also the case with higher gas density: the mean free path increases, which increases the required energy for an electron avalanche to occur. We also found that for gas densities lower than our standard value, one can observe an enhancing factor for the electric fields. This is because of the reduced drag on the dust grains in lower gas densities. This in turn increases the drift velocities of the dust grains and can therefore enhance the charge separation and the electric fields.

However, one general trend holds true for all of the studied combination of parameters discussed above, namely, the resulting electric fields are several orders of magnitude smaller than the required critical field. Nonetheless, we can still use this study to state which conditions are best for lightning to emerge. We therefore state that the most likely parameter conditions for lightning to emerge in a protoplanetary disk include having large and abundant dust grains in an area with a low cosmic ray ionization rate and low gas density.

Dust temperature. To calculate the electric field strength for the simulations where we varied the dust temperature, we had to adjust our calculations of the electric field, as the most abundant positive ion for certain temperatures is Na⁺ and not NH₄⁺. We adjusted Eq. (61) to account for the concentration of Na⁺, $n_{{\text{Na}}^{+}}$, and their mobility ${\mu }_{{\text{Na}}^{+}}$ :

$$E=\frac{{\displaystyle {\sum }_j\left(\left[{\text{Z}}_{\text{m,}j}^{+}\right]-\left[{\text{Z}}_{\text{m,}j}^{-}\right]\right)}{\upsilon }_{\text{dr},j}}{(n_{\text{I}}{\mu }_{\text{I}}+n_{{\text{Na}}^{+}}{\mu }_{{\text{Na}}^{+}})+n_{\text{e}}{\mu }_{\text{e}}}.$$(76)

We calculated the mobility of Na⁺ in a manner similar to the ${\text{NH}}_4^{+}$^:

$${\mu }_{{\text{Na}}^{+}}=K_{0,{\text{Na}}^{+}}\frac{P_0}{P}\frac{T}{T_0}.$$(77)

We took 17.5 [cm² V⁻¹ s⁻¹] for $K_{0,{\text{Na}}^{+}}$, according to Loeb (1931), who measured the mobility of Na⁺ in H₂. We chose the two normalization factors [P₀] and [T₀] to be the same as for our ${\text{NH}}_4^{+}$ calculations.

With these additional calculations in place, we could analyze the behavior of the electrical field at different temperatures. We note that in our analysis, we omitted the first point at 100 K, as including it would have required taking ${\text{C}}_3{\text{H}}_3^{+}$ into account with its mobility. We identified three different trends in the electric field curve. First, we identified a steady increase in electric field strength from 200 K up to 1300 K. This can be explained by the steady increase we observed in the concentration of negatively charged dust grains and the ${\text{NH}}_4^{+}$ relative to the electron concentration.

From 1400 K to 1600 K, we first observed a drop by nearly a factor of two and then a plateau. This can be explained by the Na⁺ ions having become the most relevant positive species. As their mobility is larger by nearly a factor of ten, it is harder for the dust grains to obtain high amounts of separation compared to the case with ${\text{NH}}_4^{+}$.

Lastly, we observed a continuous drop in electric field strength until 1800 K. This can be explained by the electron concentration becoming relevant again. As shown in Sect. 3.6, the dominant reaction creating Na⁺ also creates free electrons. Due to the very high mobility of the free electrons, they can strongly inhibit the buildup of the electrical field. Their mobility is about two orders of magnitude higher than the ones we found for our positive molecules. Hence, their effect can already be seen when their concentration is still at only about 1% of the negative dust grains and positive molecules.

In conclusion, the picture that presents itself when we modify the dust temperature is more complicated than for our other examples. In isolation, where we would just observe a system of dust grains and one other positive species, one could conclude that higher temperatures show a positive impact on electrical field strength. But as shown in our example, this is an oversimplification, and one should always consider the other effects an increase in temperature can have.

5 Discussion of the charge balance in region A under consideration of ion attachment

In this section, we discuss the charge balance in a case where a physical attachment between NH₄⁺ and negative dust grains would be possible. Although the dissociative recombination reaction Z⁻ + NH₄₊ → Z + NH₃ + H is energetically forbidden, the NH₄⁺ ion could be trapped in the electrostatic potential of a strongly negatively charged grain, either after a collision with its surface or after an inelastic collision with a gas particle close to the surface. This could occur when the electrostatic potential $e^2\frac{q}{a}$ is much larger than the molecule’s thermal energy kT. At a temperature of 100 K, this corresponds to a negative dust charge of |q|/a ≫ 6 μm⁻¹. In such cases, the NH₄⁺ molecule would be expected to stick to the surface or to stay in close vicinity of the surface.

Such an attachment reaction could be formally written as

$${\text{Z}}^{-}+{\text{NH}}_4^{+}\leftrightarrow ({\text{Z}}^{-}-{\text{NH}}_4^{+}).$$(78)

The compound $Z^{-}-{\text{NH}}_4^{+}$ would hold both the negative charges of the accumulated electrons and the positive charges of the accumulated molecular ions and would therefore be considerably less negative in the far field due to shielding by the NH₄⁺.

We assumed that the overall effects of such a mechanism would not significantly change the charge balance we found in region A but could result in an increase of the non-shielded negative dust grain charge q/a. This would introduce a new destruction path for NH₄⁺ with such a reaction, therefore decreasing the amount of NH₄⁺ in the gas. This could increase the electron concentration in the charge balance, but the $Z^{-}-{\text{NH}}_4^{+}$ compound is less negatively charged than a dust grain Z⁻, and it would result in more free electrons attaching to the compound compared to a dust grain. This process would increase the non-shielded grain charge q/a until the shielded negative grain charge becomes similar to the q/a calculated in this paper. In conclusion, such a process would act as a sink of both ${\text{NH}}_4^{+}$ and free electrons, while it would increase the non-shielded negative q/a.

However, in order to implement these effects into this work, one would need more data about the microphysical details of these processes, as such a reaction would not only need to be considered for NH₄⁺ but also for every other positive molecule with a proton affinity higher than about 8 eV (Table 2), and one would need a kinetic formulation of the back reactions. Therefore, inclusion of this mechanism goes beyond the scope of this paper, but it could be the subject of a follow-up study.

6 Summary and conclusion

We identified six different regions in our disk simulations with significantly different charging behaviors of dust and gas. The three upper regions (D, E, and F) are all photon dominated, the dust grains are positively charged, and the charge balance of the grains is characterized by an equilibrium between photoionization and electron attachment. The free electrons are generated by photoionization, leading to H⁺ in region F, C⁺ in region E, and S⁺ in region D.

In the lower disk regions (A, B, and C), we found a more complex picture with negatively charged grains. In these regions, the electron concentration is significantly lower due to the almost complete shielding from UV and X-rays. This leaves cosmic ray ionization of molecular hydrogen as the most significant mechanism for electron production. The grain charge balance is characterized by an equilibrium between electron attachment and charge exchange reactions with molecular cations.

Region C is metal poor due to ice formation, and H₃⁺ is the most abundant molecular cation. Region B has a very complex hydrocarbon chemistry, as nitrogen and oxygen are already frozen out, but carbon is not. The most abundant molecular cation was found to be C₃H₃⁺, due to its large proton affinity. We found that negatively charged dust grains Z⁻ become more abundant than free electrons inward of about r = 3 au in our standard T Tauri disk model.

In region A, we found very low concentrations of free electrons between about 10⁻¹³ and 10⁻¹⁸. The low electron concentration is because of the free electrons attaching themselves very quickly to dust grains at the high gas densities in region A. A charge balance becomes established between negatively charged dust grains Z⁻ and positively charged molecular ions, whereas the free electrons do not play an important role in the charge balance in region A. We found NH₄⁺ to be the most abundant molecular cation by far, due to its very high proton affinity, and since region A is too warm for NH₃ ice to form, NH₄⁺ is available.

Additionally, we found that the charge of a dust grain q depends linearly on its size a in almost all cases and in all disk regions (i.e., q/a = const). In disk region A, we found q/a ≈ −350 charges per micrometer to be a typical value. This linear dependence is because of the electric grain potential oc q/a, which enters the various photo-electron attachment and charge exchange rates.

We studied how the quantities characterizing the charging behavior of gas and dust depend on certain gas, dust, and cosmic ray properties in region A. These properties include the electron density n_e, the molecular cation density $n_{{\text{ion}}^{+}}$, the density of negative charges on grains Z⁻, and q/a. In particular, we found that the electron concentration depends linearly on the cosmic ray ionization rate in region A, while this rate has no effect on Z⁻and $n_{{\text{ion}}^{+}}$ . We found that at a cosmic ray ionization rate as low as 10⁻²² s⁻¹, the electron concentrations can drop to < 10⁻¹⁹. The gas density has the opposite effect, as lower gas densities cause higher electron concentrations.

The minimum dust size and the dust-to-gas ratio also have pronounced effects on n_e, $n_{{\text{ion}}^{+}}$, and Z⁻. In particular, higher dust/gas ratios cause lower electron concentrations, an effect that seems essential to understanding ring formation in early disks via disk instabilities (see Regaly et al. 2023).

Higher gas temperatures more strongly favor negatively charged grains, as the higher thermal velocities of the free electrons allow them to still collide with the grains, even if they are already strongly charged negatively. However, at temperatures roughly over 1300 K, thermal ionization of Na leads to a sudden increase of the electron concentration.

We find region A to be the most likely site for lightning to occur in disks. This is because the electron concentration in region A is so low that the gas effectively becomes an insulator, and the dynamics of the charged grains can cause an electrification.

We developed a simple model for the electrification due to turbulence in which a charge separation between the positive NH₄⁺ molecules and the negative dust grains occurs due to the centrifugal force in turbulent eddies. While this mechanism does not generate electric fields large enough to overcome the critical electric field for an electron cascade to occur in our standard T Tauri disk model, we still used it to explore areas in the disk that present the best conditions for lightning to emerge. According to this electrification model, we favor low cosmic ray ionization rates, high dust-to-gas ratios, large dust grains, and small gas densities. The model also favors warm gas temperatures, but no higher than about 1300 K, where the thermal ionization of Na inhibits the buildup of electric fields by creating too many charged gas particles that enhance the gas conductivity.

So far, we have not yet included triboelectric charging rates in our disk models. This effect, which generates negatively charged small grains and positively charged large grains by dust-dust collisions, could lead to a situation where both Z⁺ and Z⁻ are larger by orders of magnitude compared to our current results, whereas the concentrations of free electrons and molecular cations can remain very low. In that case (i.e., region X), our turbulent electrification model (Eq. (61)) predicts much larger electric fields that could potentially reach the critical E field for electron cascade (i.e., lightning).

In preliminary studies, we noticed that the grain charge chemistry is very sensitive to the heat of formation of charged dust grains $H_f^0$(Z^q) connected to the work function and the band gap. As alluded to at the end of Sect. 2.4, the underlying physics here are complicated, and the value is rather uncertain, depending on the material mixture, surface properties, and ice coating. We would like to investigate this question further in an upcoming study, where we would test different chemical potentials and show the effects this value might have on the midplane chemistry and the potential emergence of lightning.

Acknowledgements

T.B., P.W., and U.G.J. acknowledge funding from the European Union H2020-MSCA-ITN-2019 under Grant Agreement no. 860470 (CHAMELEON). U.G.J. further acknowledges funding from the Novo Nordisk Foundation Interdisciplinary Synergy Programme grant no. NNF19OC0057374.

Appendix A Additional table

Appendix B On the impact of automatically generated charge exchange and protonation reactions

In this section, we investigate the effect that automatically generated charge exchange, protonation, and endothermic reactions can have on our model. For this, we incorporated three different models (Tab. B.1): one where we do not allow exchange, protonation, or endothermic reactions to be generated and only use the network as provided, from now on called “noEX”; one where we allow charge exchange reactions to be generated, which is the case for our standard case discussed before; and one where we allow charge exchange, protonation, and endothermic reactions tobegenerated, fromnow oncalled”EX.”Wenotethatfor endothermicreactions, wesetthebarrieroftheactivationenergy for the reactions to be at 5000 Kk_b, which equals to 0.431 eV.

To illustrate the differences between the three models, we wanted to show a list of the different reactions considered for the different setups. However, such a list would consider thousands of reactions since our noEX setup consists of 5530 reactions, the standard setup of 7009, and the EX of 8672, and we would have to show the differences with example reactions, meaning reactions of a specific type not considered by a certain setup. Therefore, we have chosen instead to present a curated list of reactions that show what kind of reactions are additionally considered for the different setups (Tab. B.2).

Appendix B.1: noEx model

For the noEx model in region A, we found that, firstly, the whole region A starts a bit farther in, due to $C_3{H}_4^{+}$ dipping at a later point. Secondly, we found that C₄H⁺ is more abundant much farther in. For region B, we again found the trend that the region emerges more inward than in our standard case. Overall, we found the chemistry to be very comparable. We however found that the species H2S+ and H2CS+ are more abundant and overcome our set threshold, while SiOH+ and H3CO+ do not overcome this threshold and are therefore no longer abundant. For region C and outward, we found no notable differences.

Appendix B.2: EX model

For region A, we mainly found that compared to our standard model, Na⁺ becomes the second most abundant positive charge carrier and replaces NaH2⁺. In addition, we found that Fe⁺ is no longer abundant enough to meet our set threshold. However, we found protonated methanol $C{H}_{3}O{H}_2^{+}$ to be abundant in the EX model. For region B, we found a much less complex chemistry. In particular, C₄H⁺ and HCNH⁺, which are abundant in our standard model, do not show up anymore. We only found the sulfur bearing species, which are also abundant in the standard model, dominate as the second and third most abundant cations after ${\text{C}}_3{\text{H}}_3^{+}$. For region C, we again found no significant differences.

Appendix C Solving the chemical rate network with dust charge moments

The charge and size distribution of the dust grains are represented in our chemical models by three charge moments, ${\text{Z}}_{\text{m,}j}^{+}$, Z_m,j, and ${\text{Z}}_{\text{m,}j}^{-}$, in each dust size bin, using a small number of dust size bins j є {1, …, J}. In Sect. 2.8, we introduced the effective chemical rates for these moments for photoionization, electron attachment, and charge exchange reactions, which properly describe the overall effect of all grains on the gas and ice chemistry. On one side, the charge distribution function f_j(q) must be known to calculate the effective rate coefficients. However, on the other side, as explained in Sect. 2.7, the determination of f_j(q) requires knowing the photon flux and the concentrations of all chemical species, including electrons and molecular anions and cations.

Our formulation with dust charge moments and effective rates still conforms with the basic shape of a coupled system of first-order ODE system,

$$\frac{d{n}_i}{dt}=P_i-L_i=F_i\left(n\right),$$(C.1)

which enables us to use numerical ODE solvers to advance the chemistry. The term n_i is the particle density of chemical species i, n = {n₁, …,n_I} is the solution vector, P_i and L_i are the total production and destruction rates, F_i represents the components ofthe right-hand side (RHS) vector, and I is the numberofchemical species that include 3 × J dust moments. From n, we could in principle update f_j(q) in each call of the RHS vector, then update the effective rate coefficients, and calculate all F_i. However, this would result in an extremely slow computational method because the calculation of f_j(q) requires the computation of tens of thousands of rates for all individual charging states. In addition, we would need implicit solvers, which require the Jacobian derivative ∂F_i/∂n_j, that is, for each integration step one would need to carry out the determination of f_j(q) a couple of 100 times. Normally in astrochemistry, one can calculate the rate coefficients only once, prior to calling the ODE solver, because the rate coefficients do not depend on the particle densities, and then compute F_i and ∂F_i/∂n_j from these rate coefficients in very efficient, quick ways. In contrast, according to our formulation of the problem, the rate coefficients here depend on the particle densities in complicated ways.

We figured that there is a much more efficient way to solve our chemical equations by means of an iterative approach. To obtain a time-independent solution of our chemical equations, as used in this paper, we followed the following algorithm:

(1) Provide an initial guess of the particle densities n based on the results obtained for the previously solved grid point.

(2) Compute f_j(q) from n.

(3) Compute all rate coefficients, including the effective rate coefficients, from f_j(q).

(4) Solve F(n^new) = 0 for given rate coefficients to obtain the new solution vector n^new.

(5) Compute the relative changes of electron density n_e and mean grain charges <q_j> in all size bins j between iterations.

(6) If these relative changes are <10⁻⁵ and <10⁻², respectively, then take n^new as solution and stop the iteration.

(7) Provide a better guess for n and jump back to (2).

For step (7), we applied the following strategy. Each iteration was assigned a quality defined as Q^new = - n_e. At first instance, we took n ← n^new as the improved guess in step (7), that is, we performed a Λ iteration. In all subsequent calls, we used the previous and current electron densities and Q values to provide a better guess of the solution vector,

$$0=Q^{\text{new}}+X\left(Q-Q^{\text{new}}\right)$$(C.2)

$$n\leftarrow n^{\text{new}}+X\left(n-n^{\text{new}}\right),$$(C.3)

by determining the value for X that nullifies Q in case Q is a linear function of ne, that is, we performed a Newton step. Precisely speaking, we put $0=Q^{\text{new}}+\frac{\partial Q}{\partial n_{\text{e}}}\left(n_{\text{e}}^0-n_{\text{e}}^{\text{new}}\right)$ where we found $\frac{\partial Q}{\partial n_{\text{e}}}=\left(Q-Q^{\text{new}}\right)/\left(n_{\text{e}}-n_{\text{e}}^{\text{new}}\right)$ from the two previous iterations to determine the next improved guess of the electron density $n_{\text{e}}^0$ related to $X=\left(n_{\text{e}}^0-n_{\text{e}}^{\text{new}}\right)/\left(n_{\text{e}}-n_{\text{e}}^{\text{new}}\right)$. If the magnitude of the new quality Q^new was much smaller than the magnitude of the old quality Q, X was close to zero, and we essentially performed another Λ iteration. The reader may verify that the resulting electron density according Eq. (C.3) is indeed $n_{\text{e}}^0$. This numerical method typically converges within a couple of iterations. In a few cases, however, the method started to oscillate. We therefore monitored the Q values and saved the n vectors that were encountered during the iteration. Once we found the negative and positive entries for Q, we bracketed the solution and used the two iterations that produced the smallest positive and largest negative Q values to provide a better next guess for n in step (7).

In case we required a time-dependent solution of Eq. (C.1), not performed in this paper, we used a different strategy. The problem in this case is that with constant rate coefficients, the solution may start to alternate between two charging states with every time step taken. For example, highly negatively charged grains do not take any more electrons (i.e., the k_m,j,e become small), which means that during the next time step, the grains charge up less negatively, and so on. Our solution here was to monitor the dependencies of the effective rate coefficients k_r on electron density n_e and use the correction power law

$$k_r=k_r^{\text{ref}}{\left(\frac{n_{\text{e}}}{n_{\text{e}}^{\text{ref}}}\right)}^{p_r}$$(C.4)

whenever the ODE solver required the RHS or the Jacobian. The term $n_{\text{e}}^{\text{ref}}$ is a reference electron density from which $k_r^{\text{ref}}$ was calculated, and n_e is the actual electron density occurring during the integration. The term p_r is a rate coefficient correction power law index. The initial conditions were $n=n\left(t=0\right),n_{\text{e}}^{\text{ref}}=n_{\text{e}}\left(t=0\right)$, and p_r = 0.

The whole time-dependent approach can be summarized in eight steps:

(1) Compute f_j(q) from n and rate coefficients $k_r^{\text{ref}}$ from f_j(q).

(2) Advance the chemistry for time step ∆t from n with rate coefficients according to Eq. (C.4), resulting in n^new(t + ∆t).

(3) Compute new f_j(q) and new rate coefficients $k_r^{\text{new}}$ from n^new.

(4) Update the rate coefficient correction indices as p_r = $\left(\log k_r^{\text{new}}-\log k_r\right)/\left(\log n_{\text{e}}^{\text{new}}-\log n_{\text{e}}\right)$.

(5) Assess the uncertainty ∆n based on the difference between k and k^new and ∆t.

(6) Limit the increase of the next time step ∆t when large relative uncertainties ∆n_i/n_i are found.

(7) If t + ∆t reaches the requested total integration time, take n^new as the solution and stop the iteration.

(8) Update $t=t+\nabla t,n=n^{\text{new}},n_{\text{e}}^{\text{ref}}=n_{\text{e}}^{\text{new}},\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}}{k}_r^{\text{ref}}=k_r^{\text{new}}$ and jump back to (2).

This algorithm can slow down the computation of a t-dependent solution when strong oscillations would occur, resulting in large ∆n_i where smaller time steps are indeed required. However, the algorithm manages to sufficiently damp these oscillations and quickly re-increases the time step once the conditions have changed and the oscillations disappear.

Appendix D Dissociative molecular ion attachment reactions

Consider a reaction between dust grains and NH₄⁺ molecules, in which the NH₄+ molecule dissociates,

$$Z^{-q}+{\text{NH}}_4{}^{+}\to Z^{-\left(q-1\right)}+{\text{NH}}_3+\text{H}$$(D.1)

where z^−(q−1) is a dust grain that is less negatively charged than Z^−q by a factor of one. The reaction enthalpy, which tells us whether a reaction is exothermic or endotherm, is given by

$$\text{Δ}{H}_r=H_f\left(Z^{-\left(q-1\right)}\right)+H_f\left({\text{NH}}_3\right)+H_f\left(\text{H}\right)-H_f\left(Z^{-q}\right)-H_f\left({\text{NH}}_4^{+}\right),$$(D.2)

where H_f is the heat of formation of a substance from a chemical standard state, which we set here as purely atomic in the absence of electric fields. For neutral molecules, we used the atomization energies $H_f^0$ derived from the standard enthalpies at 0K from UMIST (Woodall et al. 2007), but for H_f (NH₄⁺), we needed to take into account the electric potential in the presence of the electric field of the dust grain at its surface

$$H_f\left({\text{NH}}_3\right)=H_f^0\left({\text{NH}}_3\right)$$(D.3)

$$H_f\left(\text{H}\right)=H_f^0\left(\text{H}\right)$$(D.4)

$$H_f\left({\text{NH}}_4{}^{+}\right)=H_f^0\left({\text{NH}}_4{}^{+}\right)-\frac{q{e}^2}{a}.$$(D.5)

The heat of formation of a negatively charged dust grains is

$$H_f\left(Z^{-q}\right)=H_f\left(\text{atoms}\to Z^0\right)+{\displaystyle \sum _{i=1}^q{H}_f}\left(Z^{-\left(i-1\right)}\to Z^{-i}\right),$$(D.6)

where H_f (atoms → Z⁰) is the heat of formation of a dust grain from the atoms it is composed of. The term H_f (Z⁻⁽ⁱ⁻¹⁾ → Z⁻ⁱ) describes the energy liberated when moving a single electron from infinity to the grain’s surface (takes energy) and then attaching it to the surface (liberates energy E₀)

$$H_f\left(Z^{-\left(i-1\right)}\to Z^{-i}\right)=-E_0+\frac{\left(i-1\right)e^2}{a}.$$(D.7)

Using the Gaussian sum formula, the total formation enthalpy of a negatively charged grain is hence

$$H_f\left(Z^{-q}\right)=H_f\left(\text{atoms}\to {\text{Z}}^0\right)-q{E}_0+\frac{e^2}{a}\frac{q}{2}\left(q-1\right).$$(D.8)

From Eq. (12), we could identify the constant E₀ = W₀ − E_bg, where W₀ is the work function and E_bg the band gap, and we assumed both to be independent of q.

The reaction enthalpy according to Eq. (D.2) is then given by

$$\begin{array}{ll}\text{Δ}{H}_r\hfill & =E_0+\frac{e^2}{a}+H_f^0\left({\text{NH}}_3\right)+H_f^0\left(\text{H}\right)-H_f^0\left({\text{NH}}_4{}^{+}\right)\hfill \\ \hfill & =E_0+\frac{e^2}{a}+P_A\left({\text{NH}}_4{}^{+}\right)-13.6\text{\hspace{0.17em}eV}\hfill \end{array}$$(D.9)

where the proton affinity is $P_A\left({\text{NH}}_4{}^{+}\right)=H_f^0\left({\text{NH}}_3\right)+H_f^0\left({\text{H}}^{\text{+}}\right)-H_f^0\left({\text{NH}}_4{}^{+}\right)$, and $H_f^0\left(\text{H}\right)-H_f^0\left({\text{H}}^{+}\right)=-13.6\text{eV}$ is the ionization energy of hydrogen. Since the term e²/a is negligible even for grains as small as 0.1 μm, Eq. (D.9) allowed us to state that a dissociative attachment reaction is exothermal when the work function minus the band gap plus the proton affinity of the respective molecule is smaller than 13.6 eV, independent of q.

There is, however, a small correction term of order e²/a that complicates things. In our derivation, this term arises from the removal of one electron and one positive molecular ion from the electric field at the grain surface by the reaction, and this term is size dependent. Similar correction terms appear in Eqs. (12) and (13). In the current code version, we used E₀ = 5.89 eV, as explained in more detail at the end of Section 2.4.

References

All Tables

All Figures