DEVELOPMENT OF A METHOD FOR THE OBSERVATION OF LIGHTNING IN PROTOPLANETARY DISKS USING ION LINES

We begin by estimating the dielectric strength of the Earth's atmosphere, in order to introduce the reader to the discharge models we later apply to the protoplanetary disk gas.

The dielectric strength of an insulating material is the maximum amplitude of the electric field that does not cause electric breakdown in the material. It is a physical property of central importance to discharge physics. Lightning on Earth is a discharge phenomenon in the air. However, it has been a longstanding mystery that lightning takes place under an electric field amplitude well below the dielectric strength of air measured in the laboratory.

The dielectric strength of the air at normal temperature and pressure (NTP; 20 °C and 1 atm) are well established from laboratory experiments (Rigden 1996):

$$\begin{eqnarray}&&{E}_{0}=30\;\mathrm{kV}\;{\mathrm{cm}}^{-1}.\end{eqnarray} \tag{ 1 }$$

The long-distance limit of Paschen's law states that the dielectric strength of gas depends linearly on the gas number density (Raizer 1997). However, in the case of the Earth's atmosphere the dependence of the dielectric strength on the number density is known to be steeper than linear. This is explained by the effects of electron loss via three-body interactions and also collisions with water vapor molecules. Empirical formulae are known (Phelps & Griffiths 1976; Takahashi 2009):

$$\begin{eqnarray}&&E={E}_{0}{\left(\displaystyle \frac{P}{{P}_{0}}\right)}^{1.5-1.65},\end{eqnarray} \tag{ 2 }$$

where E₀ and P₀ are the dielectric strength and the pressure of the air at ground level, respectively. The formula predicts the dielectric strength of air to be 17 kV cm⁻¹ at an altitude of 3 km and 10 kV cm⁻¹ at 6 km.

On the other hand, intracloud lightning is observed with an electric field amplitude of 140 V cm⁻¹ (French et al. 1996) to 150 V cm⁻¹ (Dye et al. 1986). Cloud-to-ground lightning is observed with an electric field amplitude of around 1 kV cm⁻¹ (Takahashi 1983) to 2 kV cm⁻¹ (Takahashi et al. 1999).

In this section we introduce three breakdown models that we are going to compare in this paper. Note that, in our formulation, the density functions of electrons and ions under a given electric field are well understood (see, e.g., Golant et al. 1980), and thus we use the same density function for all three breakdown models. What we do not understand well is the dielectric strength—the amplitude of the electric field at which breakdown takes place. The only difference among the three models is the assumed value of the dielectric strength.

• In the Townsend breakdown model (T) the critical electric field is such that an electron accelerated by the electric field over its mean free path gains kinetic energy large enough to ionize a neutral gas molecule. It has been widely used in the meteorological context, and also adopted into an astrophysical context, e.g., by Desch & Cuzzi (2000) and Muranushi (2010). This model explains laboratory gas discharge experiments (Equation (1)) well.
• Druyvesteyn & Penning (1940) derived the formulae for an equilibrium distribution of electrons under a constant electric field, neglecting the effects of inelastic collisions with atoms. When the electric field is weak, so that the work done by the electric field per mean free path eEl_mfp is much smaller than the electron kinetic energy, the equilibrium distribution is nearly isotropic. The distribution is expressed as the sum of the isotropic equilibrium and its first-order perturbation. The average energy and the average velocity of the mean motion, $\langle \epsilon \rangle$ and $\langle {v}_{z}\rangle ,$ satisfy $\langle \epsilon \rangle =0.43{{eEl}}_{{\rm{mfp}}}\sqrt{M/{m}_{{\rm{e}}}}$ and $\langle {v}_{z}\rangle =0.9\sqrt{{{eEl}}_{{\rm{mfp}}}}{({m}_{{\rm{e}}}M)}^{-1/4},$ respectively. Here, m_e and M are the masses of the electron and its collision partner, respectively.The Druyversteyn–Penning breakdown model (DP) assumes that breakdown takes place when $\langle \epsilon \rangle$ exceeds the ionization energy. Since the factor $\sqrt{M/{m}_{{\rm{e}}}}$ makes the average energy $\langle \epsilon \rangle$ in the DP breakdown model nearly 100 times larger than that in the Townsend breakdown model, the former model allows for breakdown under an electric field amplitude about 10⁻² times that of the latter. The model is introduced as a protoplanetary disk lightning model by Inutsuka & Sano (2005).
• Gurevich et al. (1992) have proposed the runaway breakdown model (R) and Gurevich & Zybin (2001) provided a detailed review of it. In this model, the equilibrium of the electrons with relativistic (∼1 MeV) kinetic energy, much larger than the average of the Maxwellian energy distribution, plays an important role. Because the ionization losses for electrons are inversely proportional to the kinetic energy in the non-relativistic regime, the mean free path for such fast electrons is much longer than that for thermal electrons. In the runaway breakdown model, the number of relativistic electrons grows exponentially, once the electric field is large enough to balance the ionization losses for a certain energy range of the relativistic electrons. (We define that the acceleration criterion is met for those electrons.) The ionization processes generate a spectrum of fast and slow electrons. Fast electrons that meet the acceleration criterion contribute to the exponential growth, while slow electrons are large in number and increase the degree of ionization of the matter, ultimately leadiing to its electric breakdown. Thus, runaway breakdown can take place in an electric field much weaker than that of a Townsend breakdown. The runaway breakdown model better explains observations of lightning in the Earth's atmosphere and is used as the discharge model in thunderstorm simulation studies, e.g., by Mansell et al. (2002).

In order to estimate the dielectric strength of gas, we need to compute the energy distribution of electrons. Since the interactions of electrons with even the simplest atoms and molecules have profound details (Itikawa 2000, 2001; Martienssen 2003), this requires difficult numerical computations (Chantry 1981). In this paper, we will instead resort to a simple calculation that reproduces the observed values from the discharge models.

First, we derive the dielectric strength of air at ground level from the Townsend model. Air consists of 78% N₂, 21% O₂, and 1% Ar (volume fractions). The number density of air at NTP is 2.504 × 10¹⁹ cm⁻³. The ionization energies of these chemical species are ${\rm{\Delta }}{W}_{{{\rm{N}}}_{2}}=15.6\;{\rm{eV}}$, ${\rm{\Delta }}{W}_{{{\rm{O}}}_{2}}=12.1\;{\rm{eV}},$ and ${\rm{\Delta }}{W}_{{\rm{Ar}}}=15.8\;{\rm{eV}},$ respectively. Of these ${\rm{\Delta }}{W}_{{{\rm{O}}}_{2}}\sim 12\;{\rm{eV}}$ is the smallest, so we estimate the electric field amplitude E_crit required to accelerate an electron to 12 eV, i.e., we solve 12 eV = eE_critl_mfp. The inelastic collisional cross sections (σ_inel) of N₂, O₂, Ar for 12 eV electrons are 0.8,  1.8,  0.0 × 10⁻¹⁶ cm⁻² (Itikawa 2000; Martienssen 2003). Therefore, the mean inelastic cross section of air for 12 eV electrons is 1.0 × 10⁻¹⁶ cm⁻². Therefore, ${l}_{{\rm{mfp}}}={({n}_{n}{\sigma }_{\mathrm{inel}})}^{-1}=4.0\times {10}^{-4\;}{\rm{cm}}.$ This gives

$$\begin{eqnarray}&&{E}_{{\rm{crit}}}=30\;\mathrm{kV}\;{\mathrm{cm}}^{-1},\end{eqnarray} \tag{ 3 }$$

which is in agreement with the dielectric strength of air at ground level (Equation (1)).

On the other hand, according to the DP model, the average kinetic energy of electrons under the electric field E is (Inutsuka & Sano 2005)

$$\begin{eqnarray}&&\langle \epsilon \rangle =0.43{{eEl}}_{{\rm{mfp}}}\sqrt{\displaystyle \frac{M}{{m}_{{\rm{e}}}}},\end{eqnarray} \tag{ 4 }$$

where M is the mass of the collision partner, and the dielectric strength E_crit is the solution of $\langle \epsilon \rangle ={\rm{\Delta }}W.$ In the case of air at NTP, since the mean molecular weight of air is 28.96 g mol⁻¹, M = 4.81 × 10⁻²³ g. Note that l_mfp in the DP model means the elastic mean free path ${l}_{{\rm{mfp}}}={({n}_{n}{\sigma }_{\mathrm{el}})}^{-1}=3.59\times {10}^{-5\;}{\rm{cm}}.$ l_mfp is calculated from the elastic cross sections of the elemental molecules at 12 eV (σ_el = 1.16 × 10⁻¹⁵ cm², 9.00 × 10⁻¹⁶ cm², 1.74 × 10⁻¹⁵ cm², respectively, for N₂, O₂, Ar; see Itikawa 2000; Martienssen 2003.) Therefore, E_crit = 3.38 kV cm⁻¹.

Finally, according to the runaway breakdown model the dielectric strength E_crit is the electric field amplitude where acceleration by the electric field balances the ionization loss for minimum ionizing electrons. Minimum ionizing electrons are electrons with kinetic energy ε such that for them the ionization loss is the smallest. The kinetic energy of the minimum ionizing electrons is about 1 MeV, where ionization loss is the dominant energy sink for the electrons (Gurevich & Zybin 2001). The ionization loss of an electron per unit time as a function of ε was formalized by Bethe (1930, 1932) and Bloch (1933). We use the following form of the Bethe formula from Longair (2010, chap. 5.5):

$$\begin{eqnarray}&&-\displaystyle \frac{d\varepsilon }{{dx}}=\displaystyle \frac{{e}^{4}\bar{Z}{n}_{n}}{8\pi {{\epsilon }_{0}}^{2}{m}_{{\rm{e}}}{c}^{2}}a(\gamma ),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\mathrm{where}\quad a(\gamma ) & = & {\left(\displaystyle \frac{c}{v}\right)}^{2}\left[\mathrm{ln}\displaystyle \frac{{\gamma }^{3}{{m}_{{\rm{e}}}}^{2}{v}^{4}}{2(1+\gamma ){\bar{I}}^{2}}-\left(\displaystyle \frac{2}{\gamma }-\displaystyle \frac{1}{{\gamma }^{2}}\right)\right.\\ & & \times \left.\mathrm{ln}2+\displaystyle \frac{1}{{\gamma }^{2}}+\displaystyle \frac{1}{8}{\left(1-\displaystyle \frac{1}{\gamma }\right)}^{2}\right].\end{eqnarray} \tag{ 6 }$$

Here, $\gamma ={(1-{v}^{2}/{c}^{2})}^{-1/2}$ is the Lorentz factor of the electron, $\varepsilon =(\gamma -1){m}_{{\rm{e}}}{c}^{2}$ is the electron kinetic energy, and $\bar{Z}{n}_{n}$ is the number density of ambient electrons of the matter. $\bar{I}$ is the mean excitation energy, a parameter to be fitted to laboratory experimental data. We use the value ${\bar{I}}_{{\rm{air}}}=86.3\;{\rm{eV}}$ from the ESTAR database (Berger et al. 2009).

For the case of air a(γ) takes its minimum a_min = 20.2 at γ = 3.89 or ε = 1.48 MeV. The dielectric strength E_crit is the solution of the following work-balance equation

$$\begin{eqnarray}&&{eE}-\displaystyle \frac{d\varepsilon }{{dx}}=0,\end{eqnarray} \tag{ 7 }$$

which is

$$\begin{eqnarray}{E}_{{\rm{crit}}} & = & \displaystyle \frac{{e}^{3}\bar{Z}{n}_{n}}{8\pi {{\epsilon }_{0}}^{2}{m}_{{\rm{e}}}{c}^{2}}{a}_{{\rm{min}}},\\ & = & 1.9\;\mathrm{kV}\;{\mathrm{cm}}^{-1}.\end{eqnarray} \tag{ 8 }$$

Here we summarize the three models. The dielectric strength of the gas is proportional to its number density. It is this proportional relation that leads to the constant ion velocity we present in this paper.

$$\begin{eqnarray}{E}_{{\rm{c,T}}} & = & \displaystyle \frac{{\rm{\Delta }}W}{e}{\sigma }_{\mathrm{inel}}{n}_{n}=30.1\;\mathrm{kV}\;{\mathrm{cm}}^{-1}\cdot {\left(\displaystyle \frac{{n}_{n}}{{n}_{0,\mathrm{air}}}\right)}^{1},\\ {E}_{{\rm{c,DP}}} & = & \displaystyle \frac{{\rm{\Delta }}W}{0.43}\sqrt{\displaystyle \frac{{m}_{{\rm{e}}}}{M}}{\sigma }_{\mathrm{el}}{n}_{n}=3.4\;\mathrm{kV}\;{\mathrm{cm}}^{-1}\cdot {\left(\displaystyle \frac{{n}_{n}}{{n}_{0,\mathrm{air}}}\right)}^{1},\\ {E}_{{\rm{c,R}}} & = & \displaystyle \frac{{e}^{3}{a}_{{\rm{min}}}\bar{Z}}{8\pi {\epsilon }_{0}{{mc}}^{2}}{n}_{n}=1.9\;\mathrm{kV}\;{\mathrm{cm}}^{-1}\cdot {\left(\displaystyle \frac{{n}_{n}}{{n}_{0,\mathrm{air}}}\right)}^{1}.\end{eqnarray} \tag{ 9 }$$

The minimum-mass solar nebula (MMSN) model (Hayashi 1981) has been widely used in studies of the protoplanetary disk, with fruitful results. Recent observations have contributed to sophistication of the disk models and have also reported qualitative values for the inner and outer edge radii of the disks (Kitamura et al. 2002; Andrews et al. 2009, 2010; Williams & Cieza 2011). However, such observational values for the geometry and mass of specific objects still contain uncertainty factors of 2–3 and are subjects of debate. See, e.g., Menu et al. (2014).

Some of the features common to the recent models are that the power-law index of the surface density distribution is close to 1 rather than 1.5, and that there is an exponential cut-off at the outer edge of the disk. Therefore, we use a simple model proposed by Akiyama et al. (2013) that captures these common features, and adopt the values of TW Hya reported by Calvet et al. (2002). Our disk model is as follows:

$$\begin{eqnarray}{\rm{\Sigma }}\left(r\right) & = & 6.4\times {10}^{2}{\left(\displaystyle \frac{r}{1\;\mathrm{AU}}\right)}^{-1}\\ & & \times \mathrm{exp}\left(-\displaystyle \frac{3r}{{r}_{{\rm{out}}}})\;\right.{\rm{g}}\;{\mathrm{cm}}^{-2}\;{\rm{for}}\;r\gt {r}_{{\rm{in}}},\\ {\rm{\Sigma }}\left(r\right) & = & 0\;{\rm{otherwise}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&T\left(r\right)=273{\left(\displaystyle \frac{r}{1\mathrm{AU}}\right)}^{-\displaystyle \frac{1}{2}}\;{\rm{K}}.\end{eqnarray} \tag{ 11 }$$

Here, r_out = 150 AU is the outer radius of our model disk. We also introduce an inner cut-off at r_in = 3.5 AU. The assumption of hydrostatic equilibrium leads to the vertical distribution of the gas:

$$\begin{eqnarray}\rho (r,z) & = & {\rho }_{0}(r)\mathrm{exp}\left(-\displaystyle \frac{{z}^{2}}{2{h}^{2}}\right)\\ & = & 5.08\times {10}^{-10}{\left(\displaystyle \frac{r}{1\;\mathrm{AU}}\right)}^{-\displaystyle \frac{3}{2}}\mathrm{exp}\left(-\displaystyle \frac{{z}^{2}}{2{h}^{2}})\;\right.{\rm{g}}\;{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}\mathrm{where}\quad h(r) & = & \displaystyle \frac{{c}_{s}}{{\rm{\Omega }}}\\ & = & 3.29\times {10}^{-2}{\left(\displaystyle \frac{r}{1\;\mathrm{AU}}\right)}^{\displaystyle \frac{5}{4}}\;\mathrm{AU},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{c}_{s}(r)=\sqrt{\displaystyle \frac{{k}_{B}T(r)}{\mu {m}_{p}}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{K}(r)=\sqrt{\displaystyle \frac{{{GM}}_{\odot }}{{r}^{3}}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{v}_{K}(r)=\sqrt{\displaystyle \frac{{{GM}}_{\odot }}{r}}.\end{eqnarray} \tag{ 16 }$$

Here μm_p is the mean molecular mass of the gas. Therefore the number density of H₂ is

$$\begin{eqnarray}&&{n}_{{{\rm{H}}}_{2}}(r,z)=1.52\times {10}^{14}\mathrm{exp}\left(-\displaystyle \frac{{z}^{2}}{2{h}^{2}}\right)\displaystyle \frac{{\rm{\Sigma }}(r)}{{\rm{\Sigma }}(1\;\mathrm{AU})}\;{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 17 }$$

Here we estimate the dielectric strength of the protoplanetary disk gas. In our disk model the gas density at the equatorial plane, r = 1 AU is n_{0, ppd} = 1.52 × 10¹⁴ cm⁻³. We assume that protoplanetary disk gas consists of H₂, He, CO, O₂ and their volume fractions are 0.92, 7.8 × 10⁻², 2.3 × 10⁻⁴, 1.3 × 10⁻⁴, respectively (Lodders et al. 2009, chap. 3.4.6). We used cross-section data for 15 eV electrons tabulated in Itikawa (2000) and Martienssen (2003; see Table 1), since ${\rm{\Delta }}{W}_{{{\rm{H}}}_{2}}=15.43$ eV and 15 eV is the closest table index that is found in the database.

Calculations similar to those in the previous section lead to the following values:

$$\begin{eqnarray}{E}_{{\rm{c,T}}} & = & \displaystyle \frac{{\rm{\Delta }}W}{e}({\sigma }_{\mathrm{inel}}){n}_{n}=3.0\times {10}^{-1}\left(\displaystyle \frac{{n}_{n}}{{n}_{0,\mathrm{ppd}}}\right){\rm{V}}\;{\mathrm{cm}}^{-1},\\ {E}_{{\rm{c,DP}}} & = & \displaystyle \frac{{\rm{\Delta }}W}{0.43}\sqrt{\displaystyle \frac{{m}_{{\rm{e}}}}{M}}{\sigma }_{\mathrm{el}}{n}_{n}=5.0\times {10}^{-2}\left(\displaystyle \frac{{n}_{n}}{{n}_{0,\mathrm{ppd}}}\right){\rm{V}}\;{\mathrm{cm}}^{-1},\\ {E}_{{\rm{c,R}}} & = & \displaystyle \frac{{e}^{3}{a}_{{\rm{min}}}\bar{Z}}{8\pi {\epsilon }_{0}{{mc}}^{2}}{n}_{n}=1.4\times {10}^{-3}\left(\displaystyle \frac{{n}_{n}}{{n}_{0,\mathrm{ppd}}}\right){\rm{V}}\;{\mathrm{cm}}^{-1}.\end{eqnarray} \tag{ 18 }$$

The goal of this section is to calculate the Doppler broadening of the molecular ion lines in the disk, which reflects the electric field strength in the protoplanetary disk. In order to establish the observational procedure, we calculate the collisional cross sections and the terminal velocities of the molecules. Then, we can estimate the optical depths and the spectral irradiances of the specific lines. We simulate the observational images using the calculated spectral irradiances. Finally, we establish a model discrimination procedure based on matched filtering.

We choose three ion species: HCO⁺, DCO⁺, and ${{\rm{N}}}_{2}{{\rm{H}}}^{+}$ lines, whose observations have been performed (Öberg et al. 2010, 2011). Such charged chemical species are accelerated to their respective terminal velocities by the electric field of the LMG. Let ${\varepsilon }_{I}$ be the kinetic energy of a particle of such an ion species I. We can calculate the value of ${\varepsilon }_{I}$ at the equilibrium by solving

$$\begin{eqnarray}&&{\kappa }_{I,n}{\varepsilon }_{I}={{eE}}_{{\rm{crit}}}{l}_{{\rm{mfp}},I},\end{eqnarray} \tag{ 19 }$$

Here, ${\kappa }_{I,n}=\displaystyle \frac{2{m}_{I}{m}_{n}}{{({m}_{I}+{m}_{n})}^{2}}$ is the fraction of ion energy lost per collision, and m_I and m_n are the masses of the ion and the neutral molecule, respectively (Golant et al. 1980).

As shown in Equations (9) and (18), the dielectric strength E_crit is proportional to the gas number density n_n. Let A be the proportionality factor and E_crit = An_n. Now, the mean free path ${l}_{{\rm{mfp}},I}=1/{\sigma }_{I}({\varepsilon }_{I}){n}_{n}$ is inversely proportional to the gas number density n_n. This means that the obtained kinetic energy ${\varepsilon }_{I}$ is independent of the gas number density:

$$\begin{eqnarray*}{\varepsilon }_{I} & = & {{eE}}_{{\rm{crit}}}{l}_{{\rm{mfp}},I}{({\kappa }_{I,n})}^{-1}\\ & = & \displaystyle \frac{{eA}{({m}_{I}+{m}_{n})}^{2}}{2{\sigma }_{I}({\varepsilon }_{I}){m}_{I}{m}_{n}}.\end{eqnarray*}$$

The value of A depends only on the lightning model, so it is universally the same in a protoplanetary disk. This feature is what we propose as a new signal of the lightning models in the disk. Recall that the proportionality factor A for the three lightning models is given by E_crit = An_n in Equations (18). The predicted ${\varepsilon }_{I}$ and the velocities of the ion species are shown in Table 2. In the Appendix we describe the detail of the cross-section model we have used in order to estimate the above cross sections.

Table 2.  The Terminal Velocities of the Molecular Ions for Cross-section Models XL, XM and XS

	Cross Section Model XL
	HCO+	DCO+	${{\rm{N}}}_{2}{{\rm{H}}}^{+}$
T	8.6 × 104 cm s−1	8.6 × 104 cm s−1	8.6 × 104 cm s−1
DP	2.5 × 104 cm s−1	2.5 × 104 cm s−1	2.5 × 104 cm s−1
R	2.2 × 103 cm s−1	2.2 × 103 cm s−1	2.2 × 103 cm s−1
	Cross Section Model XM
	HCO+	DCO+	${{\rm{N}}}_{2}{{\rm{H}}}^{+}$
T	4.2 × 105 cm s−1	4.2 × 105 cm s−1	4.2 × 105 cm s−1
DP	1.2 × 105 cm s−1	1.2 × 105 cm s−1	1.2 × 105 cm s−1
R	1.1 × 104 cm s−1	1.1 × 104 cm s−1	1.1 × 104 cm s−1
	Cross Section Model XS
	HCO+	DCO+	${{\rm{N}}}_{2}{{\rm{H}}}^{+}$
T	2.0 × 106 cm s−1	2.0 × 106 cm s−1	2.0 × 106 cm s−1
DP	5.9 × 105 cm s−1	5.9 × 105 cm s−1	5.9 × 105 cm s−1
R	5.2 × 104 cm s−1	5.2 × 104 cm s−1	5.2 × 104 cm s−1

Note. The terminal velocities are the characteristic features we use for observational measurement of the dielectric strength.

Download table as:  ASCIITypeset image

The column density corresponding to optical depth τ_ν0 = 1 is estimated as follows (see the appendix of Scoville et al. 1986):

$$\begin{eqnarray}{N}_{{\tau }_{\nu 0}=1} & = & \displaystyle \frac{3{\epsilon }_{0}{k}_{B}{T}_{{\rm{ex}}}{\rm{\Delta }}{v}_{{\rm{gas}}}}{2{\pi }^{2}B{\mu }^{2}\mathrm{cos}\theta }\displaystyle \frac{1}{(J+1)}\\ & & \times \mathrm{exp}\left(\displaystyle \frac{{hBJ}(J+1)}{{k}_{B}{T}_{{\rm{ex}}}}\right)\left(1-\mathrm{exp}\left(-\displaystyle \frac{h{\nu }_{0}}{{k}_{B}{T}_{{\rm{ex}}}}\right)\right),\end{eqnarray} \tag{ 20 }$$

where T_ex is the excitation temperature, ${\rm{\Delta }}{v}_{\mathrm{gas}}$ is the Doppler broadening of the target molecule, B is the rotational constant of the molecule, μ is its electric dipole matrix element, θ is the angle between the disk axis and the line of sight, hν₀ is the energy difference between the two levels, and J is the rotational quantum number of the lower state.

The optical depth of the disk with column density N is

$$\begin{eqnarray}&&{\tau }_{{\nu }_{0}}(N)=N/{N}_{{\tau }_{\nu 0}=1},\end{eqnarray} \tag{ 21 }$$

so that the intensity is

$$\begin{eqnarray}&&I({\nu }_{0})=B({\nu }_{0},T)\left(1-\mathrm{exp}(-{\tau }_{{\nu }_{0}}(N))\right).\end{eqnarray} \tag{ 22 }$$

Here $B({\nu }_{0},T)$ is the vacuum brightness of a blackbody at frequency ν₀ (see Lang 2006, chap. 2.7)).

The spectral irradiance of the disk $E(\nu )$ as a function of ν is

$$\begin{eqnarray}E(\nu ) & = & \displaystyle \frac{1}{{D}^{2}}\displaystyle \int \displaystyle \int I({\nu }_{0},r)\\ & & \times \mathrm{exp}(-\displaystyle \frac{{{mc}}^{2}d{(\nu ;{\nu }_{0},r)}^{2}}{2{k}_{B}T{{\nu }_{0}}^{2}})r\;{dr}\;d\varphi \;\mathrm{cos}\theta ,\\ \mathrm{where}\quad d(\nu ;{\nu }_{0},r) & = & \nu -{\nu }_{0}-\displaystyle \frac{{v}_{K}(r)}{c}\mathrm{cos}\varphi \;\mathrm{sin}\theta .\end{eqnarray} \tag{ 23 }$$

The integral is done in cylindrical coordinates (r, ϕ) and θ is the inclination angle of the disk.

We consider HCO⁺ 3–2, DCO⁺ 3–2, and N₂H⁺ 3–2 lines. Their frequencies are 267.56 GHz, 216.12 GHz, and 279.52 GHz, respectively. At 100 AU of the model disk T = 27 K.

For simplicity we assume that fractional abundances of the ion species are uniform within the disk. The population of the ionized species can be drastically increased or decreased as a result of lightning, but to estimate such a population is beyond the scope of this paper. We adopt the XR+UV-new chemical process model of Walsh et al. (2012), and use the abundance values at r = 100 AU, $z=3h,$ since the electric field is most likely to reach the critical amplitude at higher altitudes of the disk (Muranushi et al. 2012). Thus, we assume that the fractional abundance (relative to H₂) of N₂H⁺ is 5.3 × 10⁻¹⁰. We assume the fractional abundances of HCO⁺ and DCO⁺ to be 9.0 × 10⁻⁹ and 3.0 × 10⁻⁹, respectively, based on a paper by Mathews et al. (2013) that reports observation of enhancement in DCO abundance. Therefore, the column densities of HCO⁺, DCO⁺, and ${{\rm{N}}}_{2}{{\rm{H}}}^{+}$ are 1.8 × 10¹⁵ cm⁻², 6.0 × 10¹⁴ cm⁻², and 1.1 × 10¹⁴ cm⁻², respectively.

Assuming that there is no lightning and that the molecules are at their thermal velocities, ${N}_{{\tau }_{\nu 0}=1}$ for the three lines are 7.78 × 10⁹ cm⁻², 6.74 × 10⁹ cm⁻², and 1.39 × 10¹⁰ cm⁻², respectively. On the other hand, ${N}_{{\tau }_{\nu 0}=1}$ for the three lines are 6.15 × 10¹¹ cm⁻², 5.42 × 10¹¹ cm⁻², and 1.10 × 10¹² cm⁻², respectively, if the molecules are accelerated by the lightning electric field.

We can see that the disk is optically thick for all of the lines at 100 AU. However, all the three lines become two orders of magnitude more transparent under the effect of the critical electric field. This is a result of the molecular speed becoming faster. Consequently, observed line profiles are broadened. This Doppler broadening of the lines of the charged molecular species is the key observational feature to observe the characteristic speed of the molecules, and therefore the electric field strength in the protoplanetary disk.

We introduce seven disk models, as in Table 3. We calculate the line profiles for the three ion species with these seven disk models in order to study the ability to distinguish the lightning model from the line observations (Figure 1). The line profiles are obtained by performing the spectral irradiance integral (Equation (23)). We assume an isotropic distribution for the ion velocities, assuming that the electric field is turbulent. We simulate the channel maps using the spectral line radiation transfer code LIME by Brinch & Hogerheijde (2010). In Figure 2, we present the simulated channel maps of the HCO⁺ line for disks N, T25, and T50. We assumed that our model disk is located in the same way as TW Hya. That is, our model disk is at a distance of 54 pc and at an inclination angle of 7° (van Leeuwen 2007). Although we limit the disk parameters to this specific distance and inclination thoroughout this paper, our programs can be easily applied to other disk parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We apply the matched filtering method (North 1963) in order to distinguish the lightning model by ALMA. Matched filtering is the optimal method for discriminating models from noisy observations; it has been well studied and has wide applications not only in radio astronomy (Ellingson & Hampson 2003) but also in extrasolar planet astronomy (Jenkins & Doyle 1996; Doyle et al. 2000), in gravitational wave astronomy (Owen & Sathyaprakash 1999; Vaishnav et al. 2007; Hotokezaka et al. 2013), and even in ocean tomography (Munk & Wunsch 1979). We follow the treatment by Creighton & Anderson (2011).

Given that the noise levels for HCO⁺, DCO⁺, and ${{\rm{N}}}_{2}{{\rm{H}}}^{+}$ are 1.130 × 10⁻² Jy, 1.330 × 10⁻² Jy, and 1.800 × 10⁻² Jy, respectively, their noise spectrum power densities S_h per square arcsecond are 5.843 × 10⁻⁵ Jy² km s⁻¹, 8.094 × 10⁻⁵ Jy² km s⁻¹, and 1.483 × 10⁻⁴ Jy² km s⁻¹, respectively.

The measure-of-sensitivity ${\sigma }_{{\rm{mos}}}$ of the matched filter between two images h₁(x, y, v) and h₂(x, y, v) is

$$\begin{eqnarray}{\sigma }_{{\rm{mos}}} & = & 4{\displaystyle \int }_{{v}_{{\rm{min}}}}^{{v}_{{\rm{max}}}}{\displaystyle \int }_{(x,y)\in {\rm{image}}}\\ & & \times \displaystyle \frac{{| {h}_{1}(x,y,v)-{h}_{2}(x,y,v)| }^{2}}{{S}_{h}}\;{dx}\;{dy}\;{dv}.\end{eqnarray} \tag{ 24 }$$

Here, x and y are image coordinates in arcseconds and v is the velocity coordinate.

The measure-of-sensitivity values among the models using different lines are summarized in Table 4. The measure-of-sensitivity for any two different models is larger than 100, and the largest measure-of-sensitivity is greater than 1000. This is true even if we change the cross-section models (Table 5). Therefore an image like Figure 2 is not difficult to detect. However, no observation of a protoplanetary disk has been reported. Therefore, we can reject such forms of lightning models from observations. There are multiple alternative scenarios that observations suggest: (1) Protoplanetary disk lightning does not exist at all. (2) The probability of a protoplanetary disk with LMG is low, so that we have not observed one yet. (3) Protoplanetary disk LMG exists in the form of LMG clumps (protoplanetary "cumulonimbus clouds") much smaller than the size of the protoplanetary disks. (4) Protoplanetary disk LMG is scattered in many smaller clumps in the protoplanetary disk with a certain volume-filling factor, so that their total cross section covers a fraction of the disk image. This case is reduced to case (3) by considering the total cross section of the clumps.

Table 4.  The Measure-of-Sensitivity Values

HCO+	T25	DP25	R25
N	3729.7	3021.6	1277.1
T25		1508.4	3294.1
DP25			2608.7
HCO+	T50	DP50	R50
N	2488.4	2199.2	1277.4
T50		1418.7	2316.5
DP50			2125.5
DCO+	T25	DP25	R25
N	122.3	104.9	46.3
T25		44.7	115.3
DP25			99.3
DCO+	T50	DP50	R50
N	111.3	100.9	45.8
T50		42.8	95.3
DP50			81.6
${{\rm{N}}}_{2}{{\rm{H}}}^{+}$	T25	DP25	R25
N	5.8	5.2	2.3
T25		2.0	5.8
DP25			5.1
${{\rm{N}}}_{2}{{\rm{H}}}^{+}$	T50	DP50	R50
N	4.0	3.5	2.4
T50		2.3	3.8
DP50			3.4
3 Species	T25	DP25	R25
N	3857.9	3131.7	1325.6
T25		1555.2	3415.1
DP25			2713.0
3 Species	T50	DP50	R50
N	2603.7	2303.7	1325.6
T50		1463.9	2415.5
DP50			2210.4

Note. The measure-of-sensitivity was estimated among N, T25, DP25, and R25 models, and among N, T50, DP50, and R50 models, using either one or all three of our lines HCO⁺ 3–2, DCO⁺ 3−2, and N₂H⁺ 3−2.

Download table as:  ASCIITypeset image

We can put an upper limit on the size of such LMG clumps by thresholding the measure-of-sensitivity. For example, if the radii of LMG clumps are smaller than the values in Table 6, they are not detectable at 5.0σ . The matched filter studies show that the Townsend breakdown model is the easiest model to detect, the DP model being next, and runaway breakdown being most difficult. This tendency is explained as the wider the Doppler broadening is, the easier is the detection.