Stable three-dimensional Langmuir vortex soliton
aa r X i v : . [ n li n . PS ] A p r Stable three-dimensional Langmuir vortex soliton
Volodymyr M. Lashkin
Institute for Nuclear Research, Pr. Nauki 47, Kyiv 03680, Ukraine ∗ We present a numerical solution in the form of a three-dimensional (3D) vortex soliton in un-magnetized plasma in the model of the generalized Zakharov equations with saturating exponentialnonlinearity. To find the solution with a high accuracy we use two-step numerical method combin-ing the Petviashvili iteration procedure and the Newton-Kantorovich method. The vortex solitonwith the topological charge m = 1 turns out to be stable provided the nonlinear frequency shiftexceeds a certain critical value. The stability predictions are verified by direct simulations of thefull dynamical equation. I. INTRODUCTION
A vortex soliton (spinning soliton) is the localized non-linear structure with embedded vorticity and ringlike inthe two-dimnsional (2D) or toroidal in 3D case field inten-sity distribution, with the dark hole at the center wherethe phase dislocation takes place: a phase circulationaround the azimuthal axis is equal to 2 mπ . An integer m is referred to as topological charge. The important inte-gral of motion associated with this type of solitary wave isthe angular momentum. Spinning solitons have attracteda great deal of attention primarily in such fields of nonlin-ear physics as nonlinear optics and Bose-Einstein conden-sates (BEC), where they are the subject of considerabletheoretical and experimental research (see a recent review[1] and also [2] and references therein). In models withcubic self-focusing nonlinearity both fundamental (i. e.nonspinning, m = 0) and spinning solitons collapse (theamplitudes grow to infinity in a finite time or propagationdistance) in 2D and 3D dimensions. Vortex solitons, un-like fundamental ones, in addition to the collapse-driveninstability, may undergo an even stronger azimuthal in-stability, which tends to break the axially symmetric ringor torus into fragments, each one being, roughly speak-ing, a fundamental soliton. Stabilization of 2D and 3Dvortex solitons both against the collapse and azimuthalinstability may be achieved by means of competing non-linearities, nonlocal nonlinearities or by trapping poten-tials (harmonic-oscillator and spatially-periodic ones) [1].In plasmas, unlike optics and BEC, relativelyfew works have addressed spinning solitons. Two-dimensional spinning solitons in underdense plasma withrelativistic saturating nonlinearity were found in Ref. [3].Those vortex solitons did not collapse, but were unstableagainst azimuthal symmetry-breaking perturbations thatcaused splitting of the rings into filaments which formstable fundamental solitons. For electron-positron plas-mas with the temperature asymmetry of plasma species,2D and 3D spinning solitons were studied in Ref. [4].The vanishing saturating nonlinearity in this case doesnot sustain solitonic solutions with sufficiently large am- ∗ Electronic address: [email protected] plitudes in contrast to the ordinary saturating nonlinear-ity [1]. Under this, 2D vortex solitons with amplitudes(or, equivalently, nonlinear frequency shifts) below andabove some critical value turn out to be stable while the3D ones split into filaments due to the azimuthal insta-bility. Stable 2D spinning solitons were also found inpartially ionized collision plasma with the so called ther-mal nonlinearity [5] and at the hybrid plasma resonance[6]. In both cases the stability was due to the nonlo-cal character of the nonlinearity that is when the non-linear response depends on the wave packet intensity atsome extensive spatial domain and the nonlinear termhas the integral form. In Ref. [5] the parameter of non-locality is related to the relative energy that the electrondelivers to the ion during single collision. It was shownthat the symmetry-breaking azimuthal instability of the2D vortex soliton with | m | = 1 is fully eliminated in ahighly nonlocal regime, while the multicharge vorticeswith | m | > II. MODEL EQUATION
In the simplest case of an unmagnetized plasma, dy-namics of nonlinear Langmuir waves is governed by theclassical Zakharov equations [16]. In the subsonic limit,when ω ≪ kv s , where ω and k are the characteristic fre-quency and wave number of wave motions respectively, v s = p T e /m is the ion sound velocity, T e and m areelectron temperature and the electron mass respectively,the Zakharov equations reduce to one equation∆ (cid:18) i ∂ϕ∂t + 32 ω p r D ∆ ϕ (cid:19) − ω p ∇ · (cid:18) δnn ∇ ϕ (cid:19) = 0 , (1)for the slow varying complex amplitude ϕ of the potentialof the electrostatic electric field EE = −
12 [ ∇ ϕ exp( − iω p t ) + c . c . ] (2)at the Langmuir frequency ω p = p πe n /m , where e isthe magnitude of the electron charge, n is the equilib-rium plasma density and r D = T e / πe n is the Debyelength. Here, δn is the plasma density perturbation δn = − |∇ ϕ | πn T e . (3)In the one-dimensional case Eqs. (1) and (3) for E = − ∂ϕ/∂x reduces to the well known nonlinear Schr¨odingerequation (NLSE). In 3D case, the cubic nonlinearity inEqs. (1) and (3) results in the wave collapse, i.e. the situ-ation where the wave amplitude becomes singular withina finite time [16]. It should be noted that Eqs. (1) and(3) in 2D and 3D cases can not be reduced to the mul-tidimensional scalar NLSE. These equations, unlike 2Dand 3D NLSE, contains the vector differential operators.This difference, in particular, results in the anisotropicform of the collapsing cavity, that is a region with neg-ative density δn and trapped Langmuir waves [10]. Fur-thermore, Eqs. (1) and (3) for E = − ∂ϕ/∂r in radiallysymmetric case contains an additional centrifugal term(( D − /r ) E , where D is the space dimension, and,therefore, the field E → r → r D where Landaudamping occurs. However, some experimental observa-tions [17–19] of 3D Langmuir collapse demonstrate sat-uration of the field amplitude on a spatial scale muchlarger compared to the Debye length r D . Under this,evolution of Langmuir wave packets shows slow dynam-ics with a characteristic time scale t ≫ ω − pi (subsonic regime), where ω pi is the ion Langmuir frequency. Fromthe theoretical point of view, the wave collapse may beprevented by including some extra effects such as higherorder nonlinearities, electron nonlinearities [20–22], sat-urating nonlinearity [23], nonlocal nonlinearity [6] etc..Then, the arrest of collapse can result in the formation ofstationary structures which turn out to be (quasi)stablein some regions of parameters. We consider the case ofsaturating nonlinearity when the characteristic times ofthe nonlinear processes to exceed significantly the timeof an ion passing through the cavity, and then both elec-trons in slow motions and ions can be considered to havea Boltzmann distribution [10, 24] δnn = exp (cid:18) − |∇ ϕ | πn T e (cid:19) − . (4)Substituting Eq. (4) into Eq. (1) we get∆ (cid:18) i ∂ϕ∂t + 32 ω p r D ∆ ϕ (cid:19) − ω p ∇ · (cid:26)(cid:20) exp (cid:18) − |∇ ϕ | πn T e (cid:19) − (cid:21) ∇ ϕ (cid:27) = 0 , (5)Introducing the variables r ′ , t ′ , n ′ and ϕ ′ by r = 32 r D r Mm r ′ , t = ω − pe r Mm t ′ ,δnn = 43 mM n ′ , ϕ = T e e √ ϕ ′ (6)we rewrite Eq. (5) in the dimensionless form (accentshave been omitted)∆ (cid:18) i ∂ϕ∂t + ∆ ϕ (cid:19) − ∇ · (cid:8) [exp( −|∇ ϕ | ) − ∇ ϕ (cid:9) = 0 . (7)Equation (7) can be written in Hamiltonian form i ∂∂t ∆ ϕ + δHδϕ ∗ = 0 , (8)where the Hamiltonian is H = Z (cid:2) | ∆ ϕ | − |∇ ϕ | − exp( −|∇ ϕ | ) + 1 (cid:3) d r (9)It follows immediately from (8) that the Hamiltonian H is conserved. Other integrals of motion are the plasmonnumber (a consequence of the gauge invariance) N = Z |∇ ϕ | d r , (10)the momentum (a consequence of the translational in-variance) P = Z p d r , (11)where the momentum density is p l = i ∇ k ϕ ∗ ∇ l ∇ k ϕ − c . c . ) , (12)and the angular momentum (a consequence of the rota-tion invariance) M = Z ([ r × p ] + i [ ∇ ϕ × ∇ ϕ ∗ ]) d r . (13)In the radially symmetric case stable three-dimensionalsoliton solution of Eq. (5) was predicted in Ref. [23].Since the Langmuir field is a longitudinal one we have E k = ( k /k ) E k , where E k is the space Fourier transformof the electric field E ( k is the wave vector, k ≡ | k | ) andnext consider the magnitude E instead of the potential ϕ . In the Fourier space equation (7) can be written as i ∂ϕ k ∂t − k ϕ k − k k · Z n k k ϕ k δ ( k − k − k ) d k d k = 0 , (14)where n k is the Fourier transform of the density pertur-bation determined by Eq. (4). On the other hand, theFourier transform of |∇ ϕ | in Eq. (4) is |∇ ϕ | k = − Z k · k ϕ k ϕ ∗ k δ ( k − k − k ) d k d k . (15)By introducing the operator ˆ L acting in the physicalspace as ( f ( r ) is the arbitrary function and ˆ f ( k ) is itsFourier transform)ˆ L f ( r ) = Z ˆ f ( k ) k k e − i k · r d k (16)and taking into account that E k = − ikϕ k , Eqs. (14) and(15) become i ∂E k ∂t − k E k + ˆ L k · Z n k ˆ L k E k δ ( k − k − k ) d k d k = 0 , (17)and |∇ ϕ | k = Z ˆ L k E k · ˆ L k E k δ ( k − k − k ) d k d k , (18)respectively. Then, in the physical space one can writeEq. (17) in the form i ∂E∂t + ∆ E + ˆ L · ( n ˆ L E ) = 0 , (19)where n = exp( −| ˆ L E | ) − . (20) III. VORTEX SOLITON SOLUTION
We look for a stationary solution of (19) in the form E ( r , t ) = A ( r ) exp( iλt ) , (21) where λ is the nonlinear frequency shift. Substituting(21) into (19) and (20), one can obtain − λA + ∆ A + ˆ L · ( n ˆ L A ) = 0 , (22)where n = exp( −| ˆ L A | ) − . (23)We are interested in the stationary solution in the formof a solitary vortex with axial symmetry A ( r ) = A ( r, z ) exp( imθ ) ≡ A ( r, z ) ( x ± iy ) | m | r | m | , (24)where r = p x + y and θ are the radial coordinate andthe azimuthal angle, respectively, in the cylyndrical co-ordinates ( r, θ, z ). The real function A ( r ) should satisfythe boundary conditions at the centre, A → r → A → r → ∞ . An integer m re-ferred as the topological charge. The signs ± correspondto ±| m | . Such solutions describe either the fundamen-tal soliton, when m = 0, or the vortex soliton with thetopological charge m = 0. In contrast to the case m = 0,one can see that due to the vector type of the nonlinear-ity it is difficult to write the equation for the function A ( r, z ) if m = 0. Thus, we must solve equations (22),(23) in cartesian coordinates ( x, y, z ) just for the function A ( r ) without specifying the vortex topology. The princi-pal difficulty for finding numerical solutions of the vortextype is the convergence of usual relaxation iteration pro-cedures to the ground state, i. e. to the fundamentalsoliton rather than the vortex soliton (or any other ex-cited states). To this end, we use two stage method whichcombines the Petviashvili iteration procedure (the firststep) [12, 13] and the Newton-Kantorovich method (thesecond step) [14, 15]. The latter, as it is known, belongsto a family of universally convergent iterative methodsand it can converge to any nonfundamental solution pro-vided that the initial condition is sufficiently close to thatsolution. In the k -space equations (22) and (23) can bewritten as − λA k − k A k = Z k · k kk n k A k δ ( k − k − k ) d k d k , (25)where n k is the Fourier transform of the plasma densityperturbation n determined by (23) and | ˆ L A | k = Z k · k k k A k A ∗ k δ ( k − k − k ) d k d k . (26)Note that the nonlinearity in (25), (26) has an explicitlyanisotropic character. To calculate nonlinear terms inphysical and Fourier spaces we use the identity ( f and g are arbitrary functions)( f g ) k = Z ˆ f k ˆ g k δ ( k − k − k ) d k d k . (27) |A| FIG. 1: The vortex soliton with λ = 0 .
2. Left column: iso-surface | A ( x, y, z ) | = 0 .
15; right column: the field | A | in the x − y plane (i.e. in vertical cross-section, z = 0). Equation (25) can be written in the form G k A k = B k , (28)where G k = − ( λ + k ) and B k accounts for the non-linear term. The Petviashvili iteration method for solv-ing Eq. (28) is presented in the Appendix. An initialguess is chosen in the form of the vortex soliton (24)with A ( r, z ) = √ λr | m | exp[ −√ λ ( r + z )]. Next we re-strict the case m = 1. The convergence was controlledby stopping the iteration when the value | s − | beganto increase, where s is determined by Eq. (A.2). Thisindicates that the iteration procedure jumps off the solu-tion corresponding the vortex soliton and begins to con-verge to the fundamental soliton. Under this, the foundapproximate vortex solution has a quite acceptable accu-racy: typically, depending on λ , one can reach the value | s − | ∼ − − − . Then, as the second step, thisobtained solution is used as an initial condition in theNewton-Kantorovich iterative method as described in theAppendix. The method reduces the corresponding non-linear problem into a sequence of linear equations (A.3)which are solved by the conjugate gradient method [32].The Newton-Kantorovich method has a quadratic rateof convergence and typically after several iterations weare able to find the vortex soliton solutions with a veryhigh accuracy (up to a machine one). An example of thevortex soliton solution with λ = 0 . N = Z | ˆ L A | d r (29)is plotted, as a function of the nonlinear frequency shift λ , in Fig. 2. The minimum of the function, N cr = 863,is located at λ cr = 0 . N λ FIG. 2: The plasmon number N of the 3D vortex soliton asa function of the nonlinear frequency shift λ .FIG. 3: Self-cleaning of a randomly perturbed stable vortexsoliton with m = 1 after the application of a random per-turbation with the parameter ν = 0 .
08. Panels (a), (b) and(c) display the shape of the perturbed vortex at the initialmoment, t = 0, at t = 6 and at t = 200. The unperturbedvortex has λ = 0 . N = 1347. any dimension is based on the Vakhitov-Kolokolov (VK)criterion [7, 25]. The VK criterion states that a station-ary solution ∼ exp( iλt ) with the energy N may be stableif ∂N/∂λ >
0, and is definitely unstable otherwise. Forthe fundamental solitons (i. e. ground states), this condi-tion is both a necessary and sufficient one. This criterion,however, does not apply to the equation (7) and more-over, provides a necessary, but generally, not sufficientcondition for the stability of vortex solitons. In the regionwhere the VK criterion holds, but the vortices in NLSEare unstable, they are vulnerable to the splitting insta-bility induced by perturbations breaking the azimuthalsymmetry [1].To investigate the vortex soliton stability within theframework of Eq. (7) we solved numerically the dynam-ical equations (19) and (20) initialized with our com-puted vortex solutions with added initial perturbation.The initial condition was taken in the form E ( r ,
0) = A ( r )[1 + νf ( r )], where A ( r ) is the numerically calcu-lated solution and f ( r ) is some function correspondingthe perturbation at the initial time t = 0. We consideredtwo forms of the function f ( r ). In the first case, f ( r ) isthe white Gaussian noise with the zero mean and variance σ = 1. In the second case, f ( r ) = sin x + i cos y . Theparameter of perturbation is ν = 0 . − .
1. In bothcases perturbations break the azimuthal symmetry. Anumerical simulation shows that the vortex soliton turnsout to be stable if λ > λ cr that is in the region formallypredicted by VK criterion. We could not see any evidenceof the splitting instability at least up to times t = 400that is much larger than the characteristic nonlinear time ∼ /λ . An example of stable evolution of the vortex soli-ton with λ = 0 . ν = 0 .
08. The vortex soliton cleans up itself fromthe noise at the characteristic nonlinear time ∼ /λ andthen survives over long time.Though we did not perform the linear stability analy-sis with the corresponding eigenvalue problem, it seemsthat the complex eigenvalues (if any) accounting for thesplitting instability are sufficiently small. In this connec-tion it is interesting to compare some results concerningthe stability of the vortex-soliton solutions obtained forthe generalized two-dimensional NLSE of the form i ∂ψ∂t + ∆ ⊥ ψ + f ( | ψ | ) ψ = 0 , (30)where f ( | ψ | ) is an arbitrary function of the intensity and∆ ⊥ = ∂ /∂x + ∂ /∂y . Following the approach origi-nally proposed in Refs. [26, 27], an approximate analyt-ical estimate for the growth rate of azimuthal instabilityagainst the azimuthal perturbations ∼ exp( i Ω t + iM θ )of the vortex soliton in the model (30) can be written as[28] Im Ω = M ¯ r m,λ Re " f ′ ( | ψ | ) | ψ | − M ¯ r m,λ / (31)Here the prime stands for the derivative, ¯ r m,λ is the meanvalue of the vortex radius defined in Ref. [26] and theamplitude | ψ | ≡ | ψ m,λ | is evaluated at r = ¯ r m,λ ( m isthe topological charge and λ is the nonlinear frequencyshift). For the competing nonlinearity f ( I ) = I − I in Eq. (30), where I = | ψ | , the function f ′ ( I ) maybe negative so that Im Ω = 0 and the vortex turns outto be stable for sufficiently large amplitudes [29]. Forthe saturating nonlinearity with f ( I ) = I/ (1 + I ) thefunction f ′ ( I ) is always positive and spinning solitons areunstable in agreement with the rigorous analysis [30]. Inthe case of the exponential saturating nonlinearity f ( I ) =1 − exp( − I ), which has not been previously studied, thefirst term in the square brackets in Eq. (31), thoughpositive, is exponentially small for large amplitudes andthe instability is either practically suppressed or absent.In our case, by analogy with the model (30) one couldalso expect strong decreasing of the growth rate or fullstability of the vortex soliton. It is important to emphasize that at the critical value λ cr the vortex soliton characteristic size ∼ λ − / cr is muchlarger than the electron Debye length so that the Landaudamping is negligible. The Landau damping is still smallenough for the solitons with λ . .
2, however for larger λ values one has to take into account the damping.While the fundamental soliton corresponds to theground state of the nonlinear eigenvalue problem, thevortex soliton may be regarded as an excited state. Un-der certain conditions, one can expect emergence of theexcited states (sometimes with a sufficiently long life-time) just as in linear problems. For example, dynamicalgeneration of the vortex soliton clusters from the initialcondition in the form of the fundamental soliton withsuperimposed discontinuous phases was demonstrated inthe two-dimensional NLSE with the parabolic trappingpotential [31]. IV. CONCLUSION
In conclusion, we have shown the possibility of exis-tence in an unmagnetized plasma of the stable 3D soli-tons with the nonzero angular momentum (vortex soli-tons). Such coherent structures, as well as fundamen-tal solitons, may be considered as ”elementary bricks” ofstrong Langmuir turbulence. The real picture of strongLangmuir turbulence – collapsing cavitons or quasista-tionary structures interacting with quasilinear waves –is apparently largely dependent on excitation conditions,amplitudes and spatial scales of initial perturbations.
Appendix: Petviashvili and Newton-Kantorovichiteration schemes
The Petviashvili iteration procedure [12, 13] for Eq.(28) at the n -th iteration is A ( n +1) k = sG − k B ( n ) k , (A.1)where s is the so called stabilizing factor defined by s = R | A ( n ) k | d k R A ∗ , ( n ) k G − k B ( n ) k d k ! α . (A.2)and α >
1. For the power nonlinearity, the fastest con-vergence is achieved for α = p/ ( p − p is thepower of nonlinearity [13]. For the exponential nonlin-earity in (23) we use, depending on λ , empirical value α = 1 . − . F ( A ) = 0with A = ( A, A ∗ )). Then Newton-Kantorovich iterationscheme [14, 15] is F ′ ( A ( n ) ) A ( n +1) = F ′ ( A ( n ) ) A ( n ) − F ( A ( n ) ) , (A.3)where the prime stands for the Frechet derivative of theoperator F ( A ) at the point A defined as F ′ ( A ) A = lim h → [ F ( A + h A ) − F ( A ))] . (A.4)Calculating the Frechet derivative reduces to linearizingthe corresponding nonlinear operator in h and F ′ ( A ( n ) )is determined by F ′ ( A ( n ) ) h = − λh + ∆ h + ˆ L · ( α ˆ L h − β ˆ L h ∗ ) , (A.5) where α = [exp( −| ˆ L A ( n ) | ) − − exp( −| ˆ L A ( n ) | ) | ˆ L A ( n ) | (A.6)and β = exp( −| ˆ L A ( n ) | )(ˆ L A ( n ) ) . (A.7) [1] B. A. Malomed, Physica D , 108 (2019).[2] Yu. S. Kivshar and G. P. Agrawal , Optical Solitons:From Fibers to Photonic Crystals (Academic, San Diego,2003).[3] V. I. Berezhiani, S. M. Mahajan, Z. Yoshida, and M.Pekker, Phys. Rev. E , 046415 (2002).[4] V. I. Berezhiani, S. M. Mahajan and N. L. Shatashvili,J. Plasma. Physics , 467 (2010).[5] A. I. Yakimenko, Y. A. Zaliznyak, and Y. Kivshar, Phys.Rev. E , 065603(R) (2005).[6] V. M. Lashkin, Phys. Plasmas ,102311 (2007).[7] V. I. Petviashvili and O. A. Pokhotelov, Solitary Wavesin Plasmas and in the Atmosphere (Gordon and Breach,Reading, PA, 1992).[8] W. Horton and Y.-H. Ichikawa,
Chaos and Structures inNonlinear Plasmas (World Scientific, Singapore, 1996).[9] V. M. Lashkin, Phys. Rev E , 032211 (2017).[10] E. A. Kuznetsov, A. M. Rubenchik, V. E. Zakharov,Phys. Rep. , 103 (1986).[11] A. Ruocco, G. Duchateau1, V. T. Tikhonchuk, S. H¨uller,Plasma Phys. Control. Fusion , 115009 (2019).[12] V. I. Petviashvili, Sov. J. Plasma Phys. , 257 (1976).[13] T. I. Lakoba, J. Yang, J. Comput. Phys. , 1668(2007).[14] L. V. Kantorovich, Acta Mathematica
63 (1939).[15] A. M. Ostrowski,
Solution of Equations in Euclidean andBanach Spaces (Academic Press, New York, 1973).[16] V. E. Zakharov, Sov. Phys. JETP , 908 (1972).[17] S. V. Antipov, M. V. Nezlin, A. S. Trubnikov, I. V. Kur-chatov, Physica D , 311 (1981). [18] A. Y. Wong, P. Y. Cheung, Phys. Rev. Lett. , 1222(1984).[19] P. Y. Cheung, A. Y. Wong, Phys. Fluids , 1538 (1985).[20] E. A. Kuznetsov, Sov. J. Plasma Phys. , 178 (1976).[21] M. M. ˆSkoriˆc, D. ter Haar, Physica , 211 (1980).[22] T. A. Davydova, A. I. Yakimenko, Yu. A. Zaliznyak,Phys. Lett. A , 46 (2005).[23] E. W. Laedke and K. H. Spatschek, Phys. Rev. Lett. ,279 (1984).[24] P. K. Kaw, G. Schmidt, and T. Wilcox, Phys. Fluids ,1522 (1973).[25] N. G. Vakhitov, A. A. Kolokolov, Radiophys. QuantumElectron.