Stable Langmuir solitons in plasma with diatomic ions

We study stable axially and spherically symmetric spatial solitons in plasma with diatomic ions. The stability of a soliton against the collapse is provided by the interaction of induced electric dipole moments of ions with rapidly oscillating electric field of a plasmoid. We derive the new cubic-quintic nonlinear Schrodinger equation which governs the soliton dynamics and numerically solve it. Then we discuss the possibility of implementation of such plasmoids in realistic atmospheric plasma. In particular, we suggest that spherically symmetric Langmuir solitons, described in the present work, can be excited at the formation stage of long-lived atmospheric plasma structures. The implication of our model for the interpretation of the results of experiments for the plasmoids generation is discussed.


Introduction
Stable spatial solitons are observed in the studies of optical phenomena (Segev, 1998), in solid states physics (Burger et al., 1999), and in the plasma research (Antipov et al., 1981). Typically a soliton can be stable owing to a certain nonlinearity. For example, spatial plasma solitons, described in frames of the classical electrodynamics, can exist due to the combined action of electron-ion and electron-electron nonlinear interactions. The former one was found to be focusing (Zakharov, 1972), whereas the latter interaction can be defocusing (Kuznetsov, 1976;Skorić and ter Haar, 1980;Davydova et al., 2005). Recently Haas and Shukla (2009) suggested that, if one accounts for the additional quantum pressure of electron gas, it may explain the appearance of two and three dimensional Langmuir solitons in dense plasmas.
In the present work we stall study the existence of stable Langmuir solitons in a plasma with ions possessing induced electric dipole moments (EDM). We shall demonstrate that the interaction of ion's EDM with a rapidly varying electric field of a plasma oscillation results in the arrest of the Langmuir wave collapse. Thus the existence of a spatial soliton becomes possible.
In our analysis we shall choose a plasma with diatomic nitrogen ions, which corresponds to a 1 realistic atmospheric plasma. Thus our results may be applied for the theoretical description of long-lived plasma structures observed in the atmosphere (Bychkov et al., 2010). It is interesting to notice that the role of EDM of charged particles for the explanation of the stability of atmospheric plasmoids was also discussed previously by Bergström (1973); Stakhanov (1973).
This work is organized as follows. In Sec. 2 we consider the general description of nonlinear waves in plasma. We introduce the new ponderomotive force, associated with EDM, which acts on the ion component of plasma. Then we derive a system of nonlinear equations for the electric field amplitude and the perturbation of the ion density. In Sec. 3 we reduce this system to a single nonlinear Schrödinger equation (NLSE), containing cubic and quintic terms, for the envelope of the electric field. This equation is analyzed numerically for the case of a radial plasma oscillation. We show that both axially and spherically symmetric stable solitons can exist. Then, in Sec. 4, we discuss a possible application of our model of stable spatial solitons to the description of long-lived atmospheric plasmoids as well as for the interpretation of some experimental results. Finally, in Sec. 5, we briefly summarize our results. The calculation of the permittivity of the ion component of plasma is presented in Appendix A.

Nonlinear Langmuir oscillations in plasma
In this section we shall derive the basic nonlinear equations for the description of Langmuir waves in plasma accounting for a ponderomotive force acting on nonpolar diatomic ions. If we study electrostatic plasma oscillations, i.e. when the magnetic field is zero, B = 0, the motion of the electron component of plasma obeys the following plasma hydrodynamics equations: where n e is the number density of electrons, v e is the electron velocity, m is the mass of an electron, e > 0 is the proton charge, E is the strength of the electric field, and p is the pressure. We should also consider Maxwell and Poisson equations for the electric field evolution, where n i is the ion number density and v i is the ion velocity.
In the zeroth approximation only electrons participate in a plasma oscillation, with the number density of ions being approximately constant: n i ≈ n 0 = const. Thus we may present the electric field in the form, where ω e = 4πe 2 n 0 /m is the Langmuir frequency for electrons and E 1 is the amplitude of the electric field. It should be noted that in the following (see, e.g., Sec. 3) we shall study axially and spherically symmetric plasma oscillations. In this case one can find a scalar potential φ 1 such as E 1 = −∇φ 1 in Eq. (3).
In a realistic situation ions will also participate in a plasma oscillation. Thus the ion density becomes n i = n 0 + n(r, t), where n is the perturbation of the ion density. It leads to the appearance of higher harmonics omitted in Eq. (3). The plasma hydrodynamic equations for the description of the ions evolution have the form (Škorić and ter Haar, 1980), where M is the ion mass, F i is the force acting of ions. The reason why one can omit the ion pressure term in Eq. (4) will be discussed at the end of this section.
Using the quasineutrality of plasma we can find that (Zakharov, 1972) where ϕ s is the slowly varying part of the electric potential, T e is the electron temperature, and U pm = |E 1 | 2 /(4πn 0 ) is the potential of the ponderomotive force which acts on a charged particle in a rapidly oscillating electric field given in Eq.
(3). Supposing that ions are mainly involved in the slow motion of plasma we get that F i = −e∇ϕ s . Finally, using Eqs. (1)-(5) one arrives to the system of equations of Zakharov (1972). More detailed derivation of the nonlinear plasma evolution equations which include electron-ion and electron-electron interactions can be found in the works by Kuznetsov (1976); Skorić and ter Haar (1980). It should be noted that Eq. (4) is derived under the assumption of ions having point like charges. However realistic atmospheric plasma contains mainly nitrogen or oxygen ions, which are diatomic (see also Sec. 4). In this case the simplified ion hydrodynamics Eq. (4) is incomplete since it does not take into account the internal structure of ions.
A diatomic molecule is nonpolar, i.e. it cannot have an intrinsic EDM because of the symmetry reasons. Nevertheless, this kind of molecules can acquire EDM, p i = α ij E j , in an external electric field. Here (α ij ) is the polarizability tensor. Hence the additional force, F pol = (p∇)E, will act on this particle placed in an external inhomogeneous electric field. Thus, if we study the plasma with diatomic ions, in Eq. (4) one should replace F i = −e∇ϕ s → −e∇ϕ s + f pol /n i , where f pol is the volume density of the ponderomotive force related to the matter polarization.
If an ion is diatomic and possesses an axial symmetry, one can always reduce the polarizability tensor to the diagonal form, (α ij ) = diag α ⊥ , α ⊥ , α , where α ⊥ and α are transversal and longitudinal polarizabilities. Now the expression for f pol can be obtained using Eq. (21) and the general technique developed by Tamm (1979), as where ε is the permittivity of the ion component of plasma, T i is the ion temperature, α = (2α ⊥ +α )/3 is the mean polarizability of an ion, and ∆α = α − α ⊥ . It should be noted that the general expression for the ponderomotive force f pol was derived by Tamm (1979) under the assumption of static fields with (∇ × E) = 0. However, as we mentioned above, we study electrostatic plasma oscillations with zero magnetic field (see also Sec. 3). Thus Eq. (6) remains valid.
(1)-(6) we get the following nonlinear coupled equations for the amplidute of the electric field, (7) and for the perturbation of the ion density, where r D = T e /4πe 2 n 0 is the Debye length and c s = T e /M is the sound velocity in plasma. To derive Eq. (8) we take into account that E 4 = 6|E 1 | 4 while averaging over the time interval ∼ 1/ω e . The first quadratic term in the rhs of Eq. (8) corresponds to the direct interaction of charged ions with the electric field whereas the second quartic term there, ∼ ∇ 2 |E 1 | 4 , is related to the induced EDM interaction. Hence the contribution of this second term to the Langmuir waves dynamics should be typically smaller. However, as we shall see in Sec. 3, in some cases it is the EDM term which arrests the collapse of Langmuir waves.
It should be noted that in Eq. (8) we neglect the contribution of the ion temperature to the sound velocity. Such a contribution would correspond to a nonzero ion pressure term in Eq. (4). Since we suppose that T i ≪ T e , we can omit the ion pressure. However we keep the ion temperature in the quartic nonlinear term in Eq. (8). In the rhs of Eq. (8) we also neglect term ∼ −n 0 α ∇ 2 |E 1 | 2 /M , which is small compared to the contribution of the Miller force. Indeed, the ratio of these terms is ∼ n 0 α . In Sec. 4 we shall use the following values of n 0 and α : n 0 ∼ 10 21 cm −3 and α ∼ 10 −24 cm −3 . For such a parameters, this ratio ∼ 10 −3 , that justifies the validity of Eq. (8).
Finally we notice that the quartic nonlinear term in Eq. (8) blows up if T i → 0. It means that our model is not valid at the very low ion temperature. However, in Sec. 4, where we shall discuss a possible application, we consider only reasonable, relatively high values of T i . Note that at extremely low temperatures quantum corrections to the ion and electron motion become important.

Cubic-quintic nonlinear Schrödinger equation
In this section we shall derive the nonlinear Schrödinger equation for the amplitude of the electric field. This equation will be analyzed in a particular case of a radial plasma oscillation. We shall numerically solve it and find the characteristics of Langmuir solitons. The soliton stability will be examined.
Let us suggest that the density variation in Eq. (8) is slow, i.e. ∂ 2 n/∂t 2 ≪ c 2 s ∇ 2 n. In this subsonic regime Eqs. (7) and (8) can be cast in a single NLSE, (9) which has both cubic and quintic nonlinear terms. NLSEs analogous to Eq. (9) were previously examined by Desyatnikov et al. (2000) in connection to the studies of the light bullet propagation in crystals. Note that in Eq. (9) we omit the index "1" in the amplitude of the electric field, i.e. E 1 ≡ E, in order not to encumber the formulas.
As we mentioned in Sec. 2, we shall examine axially or spherically symmetric plasma oscilla-tions, i.e. E = Ee r , where e r is a unit vector in radial direction and E is a scalar function. Introducing the following dimensionless variables: we can represent Eq. (9) in the form, which contains no dimensionless parameters.
Here d = 2, 3 is the dimension of space. One can check by the direct calculation that the plasmon number and the Hamiltonian are the integrals of Eq. (11).
Here Ω 2 = 2π and Ω 3 = 4π are the solid angles in two and three dimensions. We shall look for the solution of Eq. (11) as ψ(x, τ ) = e iλτ ψ 0 (x), where λ is a real number meaning the dimensionless frequency shift. By the proper choice of the phase we can always make the function ψ 0 to be real. The corresponding ordinary differential equation for the function ψ 0 is solved numerically using the MATLAB program.
Firstly, we analyze the stability of axially and spherically symmetric solitons by plotting N (λ) dependence. It is shown in Fig. 1(a) for 2D case and in Fig. 2(b) in 3D case. One should notice that in 2D case ∂N/∂λ > 0 in a quite broad range of λ. Thus according to Vakhitov and Kolokolov (1973) criterion (VKC), this kind of 2D solitons is stable. In 3D case the accuracy of calculations is significantly lower than that in 2D situation. To build a smooth N (λ) curve in Fig. 2(a) the least squares method was used since the points on this plot, obtained with numerical simulations in MATLAB, have rather big spread, especially at large λ. Meanwhile one can see that unstable and stable solitons coexist in 3D case. In Fig. 2(a) we get that ∂N/∂λ < 0 at λ 0.26. Applying VKC we conclude that this branch corresponds to unstable solitons. However, if λ 0.26, we obtain that ∂N/∂λ > 0, that corresponds to stable solitons.

Application
In this section we discuss the application of the results obtained in Secs. 2 and 3 for the studies of a natural atmospheric plasmoid called a ball lightning (BL).
BL is a glowing object appearing mainly during a thunderstorm (Keul, 2013). Despite the existence of numerous BL models (Bychkov et al., 2010), it is likely to be a plasma based phenomenon. Dvornikov and Dvornikov (2006);Dvornikov (2011Dvornikov ( , 2012aDvornikov ( ,b, 2013 put forward a hypothesis that BL is based on a radial oscillation of plasma and discussed both classical and quantum approaches to this problem. Note that analogous BL model was also considered by Fedele (1999);Shmatov (2003); Tennakone (2011).
Note that, in the BL model based on radial plasma oscillations, it is not required to have any special force which maintains the shape of a plasma structure. On the contrary, in case of a model based on the static distribution of electric charges inside a plasmoid one needs a self-consistent interaction, e.g., of the magnetic origin, which sustains a plasma object in equilibrium. However, in this case plasma will radiate electromagnetic waves and a plasmoid will lose its energy within the millisecond time scale.
We already mentioned in Sec. 2 that the magnetic field in the system of radially oscillating particles is equal to zero. This fact may explain the relative stability of a plasma structure since it does not lose energy by the electromagnetic radiation. Moreover, here we are using the concept of collisionless plasma. Later in this section we discuss how one can account for the visible radiation of BL described in frames of our model.
Let us discuss the conditions for the existence of a spherical nonlinear Langmuir soliton in the atmospheric plasma. Firstly, an atmospheric plasma can contain N + 2 (≈ 78%) and O + 2 (≈ 21%) ions. These ions are diatomic and nonpolar. We have already mentioned in Sec. 2 that diatomic ions cannot have an intrinsic EDM, whereas an induced EDM is well possible. Thus our results may be applied for the description of atmospheric plasmoids. Now let us evaluate the characteristics of the plasma structure. For the stability of a soliton its size should be greater than r D -otherwise thermal electrons will cross the soliton leading to the development of a turbulence. One can see in Figs. 1(b) and 2(b) that in dimensionless variables the typical soliton size can be up to 10. Using the values of the polarizabilities of a neutral nitrogen molecule (Alms et al., 1975), α ij ∼ 10 −24 cm −3 , and Eq. (10), we get that if a plasmoid becomes stable. The analysis of the soliton stability is based on the perturbation theory which requires that δω < ω e , where δω is the frequency shift. Using the parameter κ defined in Eq. (14) we get that this condition is equivalent to, κ > 4.7 × 10 3 , where we use Figs. 1(a) and 2(a) to obtain the typical value of λ ∼ 0.4. The ion density inside a plasmoid is lower than the ambient density (Zakharov, 1972). Integrating Eq. (8) in the subsonic regime and using Eq. (10) we should impose the constraint on the normalized ion density variation, where ν 0 (x) = |ψ 0 | 2 − |ψ 0 | 4 . Note that n is Eq. (16) is negative. The function ν 0 (x), corresponding to different numerical solutions of Eq. (11), is shown in Figs. 1(c) and 2(c) for axial and spherical plasma structures. As one can see in these plots, the maximal value of ν 0 is in the (0.05 − 0.25) range. In 2D case, the bigger the parameter λ is, the greater the maximum of ν 0 is. Thus, if we take λ 0.4 and κ > 6 × 10 3 , then using Eq. (16), we get that |n| < n 0 inside the whole plasmoid. In 3D case, one can see in Fig. 2(c) that the smallest value of the maximum of ν 0 (x) is archived for λ min ≈ 0.25. Note that λ min corresponds to the minimum of the function N (λ), cf. Fig. 2(a). Again considering that λ is close to λ min , one obtains that the ion number density is positive in a spherical plasma structure if κ > 6 × 10 3 .
The conditions of the plasma structure stability in Eqs. (14)-(16) are satisfied if we suppose that n 0 ∼ 10 21 cm −3 , T i = 300 K and T e = 10 6 K. Now let us estimate the typical size of a plasmoid requiring that R eff > r D . For the chosen electron density and temperatures we get the lower bound for the radius: R eff (10 −7 ÷ 10 −6 ) cm. Using the numerical solutions of Eq. (11), shown on Figs. 1(b) and 2(b), as well as obtained estimates for λ, we get the upper bound for the radius: R eff 10 −5 cm. Note that at n 0 ∼ 10 21 cm −3 and R eff ∼ 10 −5 cm, still there is a significant number of charged particles inside the plasmoid: 4 3 πR 3 eff n 0 ∼ 10 6 . It should be also noticed that the plasmoid radius is at least one order of magnitude greater than the intermolecular distance: It should be mentioned that the obtained small size of a plasma structure is very close to the estimates obtained by Dvornikov and Dvornikov (2006);Dvornikov (2012aDvornikov ( ,b, 2013, where radial nonlinear plasma oscillations were studied using the quantum approach. Note that a theoretical model of BL having a nano-sized kernel was developed by Alexeev (2008) using the methods of generalized quantum hydrodynamics. The existence of plasma structures having the size, which lies in the nanometer range, was theoretically suggested by Abrahamson and Dinniss (2000). Recently their hypothesis was experimentally tested by Dikhtyar and Jerby (2006); Paiva et al. (2007).
According to observations, the visible dimension of BL is about several centimeters (Doe, 2013). However its core should be much smaller. Otherwise it is difficult to explain how BL passes through small holes or cracks in dielectric materials without destroying them (Doe, 2013). Therefore the visible size of a plasmoid is likely to be caused by some auxiliary effects.
We can also assume that BL is a composite object. This assumption is based on the fact that sometimes BL can decay into pieces (Bychkov et al., 2010). The composite structure of plasma objects was also confirmed experimentally. The photographs taken by Paiva et al. (2007) show that the visible size of artificial plasmoids lies in the cm range. However, the analysis of materials which were in contact with these plasma objects demonstrates that plasmoids consist of multiple small kernels. The models of a composite BL were previously discussed by Meshcheryakov (2007);Dvornikov (2012b). It should be noted that the model of BL, consisting of multiple oscillating plasmoids, can explain the visible radiation of such a plasma structure. This radiation is due to the decay of plasma oscillations in some small fraction of oscillating kernels, whereas the majority of them are stable.
The results of the present work can be used to explain the formation of a plasmoid in natural conditions. Let us suppose that during a thunderstorm a high tension, related to a linear lightning strike, is applied to a small object with the appropriate size. It can be a natural spike, a piece of dust or soot etc. Then this object turns out to be embedded in plasma with n 0 ∼ (10 19 − 10 20 ) cm −3 . Moreover, this object will be a source of divergent waves, where the electron density can rise almost one order of magnitude (Gorbunov et al., 2010;Bulanov et al., 2012) and reach the value required to excite an axial or a spherical Langmuir soliton described in our work. Note that the plasma heating up to ∼ 10 6 K is also possible in the vicinity of this object. We remind that the typical electron temperature in a linear lightning bolt is (Rakov, 1998) T e ≈ 3 × 10 4 K. Note that at high electron temperature, atoms can be multiply ionized. It will also increase the electron density.
It should be mentioned that ions are unlikely to be involved in the divergent waves because of their low mobility. Thus the density of ions should keep the initial value, i.e. be lower than that of electrons, in some region near the object. It can explain the formation of a cavern (see Eq. (16) and the work by Zakharov (1972)).
Therefore we may consider a Langmuir soliton, which involves the interaction of induced EDM of ions with the oscillating electric field, as proto-BL, i.e. an atmospheric plasma structure at its initial stages of evolution. Under certain conditions, when other (maybe quantum) nonlinear effects become important, this proto-BL can then be transformed into a glowing object identified as BL.
We also mention that besides the explanation of the formation of natural plasmoids, the model of radial plasma oscillations, involving diatomic ions, can be used for the interpretation of the results of the experiments performed by Klimov et al. (1994);Dimitrov et al. (1994); Kirko et al. (1995), where spherical luminous structures were obtained in electric discharges in liquid nitrogen. Indeed, in case of a plasmoid in liquid nitrogen we can take that T i = 77 K and n 0 ≈ 3.44×10 21 cm −3 , which corresponds to 10% ionization. To satisfy conditions in Eqs. (14)-(16) one should require that T e ∼ 10 4 K. Such electron temperatures are achievable in laboratory plasmas. It makes rather plausible interpretation of the results of Klimov et al. (1994);Dimitrov et al. (1994); Kirko et al. (1995) in frames of our model.

Conclusion
In the present work we have studied stable spatial plasma structures possessing axial and spherical symmetry. The stability of solitonlike plasmoids against the collapse is provided by the defocusing interaction of the induced EDM of ions with the rapidly oscillating electric field. Note that an ion was supposed to be diatomic and nonpolar. In Sec. 2, starting from the complete set of plasma hydrodynamic equations we have derived the basic nonlinear equations for the envelope of the electric field and for the perturbation of the ion density (see Eqs. (7) and (8)). In Sec. 3 we have reduced these equations to a single cubic-quintic NLSE (11), written in the dimensionless variables. Then Eq. (11) was solved numerically. Applying VKC, we have found that mainly stable spatial solitons can exist in 2D case, whereas in 3D case both stable and unstable solitons are present.
It should be noted that VKC was originally formulated by Vakhitov and Kolokolov (1973) for a NLSE for a scalar "wave function" Ψ. In 2D or 3D cases, the kinetic term of such an equation contains the Laplace operator of Ψ, which reads Ψ ′′ + d−1 r Ψ ′ for the function depending only on the radial coordinate. The kinetic term of NLSE (11), derived in our work, has different structure. Nevertheless, as shown by Davydova et al. (2005), the application of VKC to the analysis of NLSE analogous to Eq. (11) is justified.
In Sec. 4 we have discussed a possible application of the obtained results to the studies of a long-lived atmospheric plasma structure, BL. It was found that in a realistic atmospheric plasma, mainly composed of N + 2 and O + 2 ions, one may expect that the initial stages of BL evolution can be described in frames of the plasmoids model developed in the present work. We have also suggested one of the possible ways of the creation of this kind of plasma structures in natural conditions. Besides the description of natural plasmoids, in Sec. 4 we have considered the application of our model for the interpretation of the results of experiments where plasma structures were generated in liquid nitrogen.
As we mentioned in Sec. 1, besides the mechanism proposed in the present work, there are other ways to arrest the collapse of a Langmuir wave. For example, it was shown by Kuznetsov (1976);Škorić and ter Haar (1980); Davydova et al. (2005) that the electronelectron nonlinear interaction prevents the soliton collapse. It should be noticed that in this case higher harmonics are generated in the system. However, a plasma structure, supported by this nonlinearity, was demonstrated by Dvornikov (2011) to exist only in the upper ionosphere, where the density is very low. We remind that the majority of the BL observations indicate that this phenomenon happens in the lower troposphere (Bychkov et al., 2010).
Note that in our work we have used polarizabilities of a neutral nitrogen molecule rather than of N + 2 , which in fact can differ significantly. Nowadays no measurements of N + 2 polarizabilities have been made. This fact can change the estimate of the critical background density, necessary for the plasmoid existence, towards its reduction.
A Polarization of a nonpolar diatomic gas in an external electric field The collective response of plasma results in the concept of "dressed" particles. It means that the Coulomb interaction of a charged particle is replaced by the Debye-Hückel potential, 1 r exp − r r D . If the plasma density is quite high, i.e. the Debye length is short, we may suppose that the electric charge of an ion is quite perfectly screened by surrounding electrons. In this case the Hamiltonian of a diatomic ion possessing an axial symmetry in an external electric field reads where M ⊥ is the vector of the angular momentum of an ion perpendicular to the ion axis, I is the moment of inertia of an ion, U = −(pE) = −E 2 α ⊥ sin 2 θ + α cos 2 θ is the potential energy of a polarized ion in an external electric field, and θ is the angle between the ion axis and the electric field direction. On the classical level the thermodynamic properties of a gas can be calculated on the basis of the canonical partition function, Z = z N , where N is the total number of ions. The reduced partition function z which includes the rotational degrees of freedom of an ion has the form, where dΩ = 2π sin θdθ is the solid angle differential.
Calculating the integral over the angular momentum components, we can express z rot in the following form: where ξ = E 2 ∆α/T i . The polarization of the gas in an external electric field can be calculated on the basis of the expression for the free energy F = −T i ln Z, as where V is the volume of a gas. Here we account for the decomposition of z ′ in Eq. (19). Basing on Eq. (20) we obtain the permittivity of the ion component of plasma ε = 1 + 8πn i α + 4 45 which was used in Sec. 2 to derive the ponderomotive force acting on ions. Using Eqs. (19) and (20), one can get the quartic contribution to the permittivity in Eq. (21). It has the form, The obtained correction to ε is small. However, using the fact that ∆ε in Eq. (22) is positive as well as the formalism developed in Secs. 2 and 3, we get that, if one accounts for this correction, it will defocuse Langmuir waves and further stabilize solitons.
Finally we mention that our calculations are based on the assumption of the constant electric field whereas in Secs. 2 and 3 the electric field is supposed to oscillate with the high frequency ∼ ω e . According to Alms et al. (1975), the typical frequency associated with the polarizability of a molecule is ∼ 10 15 Hz. This value is several orders of magnitude greater than plasma frequencies in realistic plasmas. Thus the approximation of the constant electric field is valid.