Observation of sharply peaked solitons in dusty plasma simulations

We report on observations of cusp-like sharply peaked dust acoustic solitons in a one-dimensional fluid simulation of a dusty plasma system. These nonlinear singular structures, dithering at the wave breaking point, are found to be long-lived entities that emerge spontaneously from a variety of large-amplitude initial perturbations. Such excitations can play a significant role in wave breaking and wave-particle dynamics in dusty plasmas and may explain some recent experimental observations.


Introduction
A dusty plasma medium consists of micron-sized particles immersed in an ionized gaseous medium of freely moving ions and electrons. The massive dust particles become heavily charged by the flux of ions and electrons constantly impinging on them [1]. Dusty plasmas occur naturally in many astrophysical situations ranging from cometary tails, through planetary rings 1 Authors to whom any correspondence should be addressed. 2 to noctilucent clouds occurring in the upper layers of the atmosphere [2,3]. They are also found in many industrial and laboratory devices such as plasma torches and magnetic fusion machines such as tokamaks [4]. Like ordinary gaseous plasmas, dusty plasmas also exhibit a variety of collective behaviours with the dust component enriching the dynamics.
The linear and nonlinear characteristics of collective oscillations in a dusty plasma system have been the subject of many theoretical and experimental studies in recent years ( [5] and references therein). The dust acoustic wave (DAW), in particular, has received a great deal of attention. This low-frequency longitudinal mode is an analogue of the ion acoustic mode found in normal electron-ion plasmas with the massive dust particles now providing inertia and the pressure contributions for sustaining the wave coming from both electrons and ions [6]. Its linear propagation characteristics are well understood theoretically and have also been widely confirmed experimentally ( [7] and references therein; [8]). There is also an extensive body of theoretical literature on the subject of nonlinear evolution of DAWs, mostly centred on the idea of excitation of solitons [6,[9][10][11][12]. These localized one-dimensional (1D) pulse structures belong to a class of exact solutions to integrable nonlinear partial differential equations (PDEs), such as the Korteweg-DeVries (KdV) equation (and its generalizations), and have been applied to a wide variety of physical systems. The 1D idealization holds good in many practical situations where the time scale for the breakup of the soliton due to bending instabilities in the perpendicular direction is quite long compared to the propagation time of the soliton over an experimental length. Likewise, collisional effects and other damping mechanisms can be effectively controlled by varying the plasma parameters, ensuring that the solitons have a long lifetime. Consequently, KdV solitons and other similar 1D solitons have been experimentally studied by many researchers (for a detailed review see [13]) ever since the first laboratory observations of Ikezi et al [14]. More recently, such pulses were also observed in dusty plasmas [12]. Other nonlinear structures that have been studied in the context of DAWs are shock waves and 2D vortices.
An important class of nonlinear solutions that has not received much attention in a dusty plasma is that of sharply peaked solitons. These solutions occur near the wave breaking amplitude and have a spatial structure that is distinctly different from the smooth pulse soliton solutions. Mathematically, these sharply peaked solitons, dithering at the wave breaking amplitude, are singular in nature and may exhibit singularities in the field and/or its higher derivatives at its maximum amplitude [15]. They are similar to singular structures observed in other situations such as slowly propagating envelope solitons in laser-plasma interactions. The physical mechanism responsible for the excitation of such structures is the incident laser pulseproduced ponderomotive pressure driving the ion waves to near wave breaking amplitudes [16]. Similar nonlinear states involving upper hybrid waves have also been reported for theoretical studies of electron-ion magnetized plasmas [17][18][19] and in numerical simulations [20] of the Jaulent-Miodek equations [21]. More recently, such singular structures were observed in experiments with dusty plasmas. Schwabe et al [22] observed cusp-like structures on the surface of voids formed in a dusty plasma system created in a microgravity environment of the outer space. Teng et al [23] carried out wave breaking experiments in the laboratory and saw cusp-like structures in the bulk of a compressible dusty plasma fluid. This paper is motivated to some extent by these recent experimental observations and by a desire to theoretically explore the accessibility and evolution of sharply peaked solitons for DAWs. To achieve this goal we carried out large-scale 1D numerical simulations of bulk longitudinal oscillations in a dusty plasma using a standard three-fluid model that has been 3 widely used in the literature for solitonic studies. Our simulations show the spontaneous excitation of a single sharply peaked propagating structure or a chain of sharply peaked pulses when large-amplitude initial perturbations in the form of a single pulse or a wave train are imposed on the system. These structures are quite long lived (over several dust acoustic frequency periods) and are accessible from a range of initial conditions. Our simulation results show a remarkable resemblance to the cusp-like structures observed in the above-mentioned recent experiments and may provide important insights into the role of such structures in the wave breaking and wave-particle interaction processes of dusty plasma systems.

Governing equations
Our simulation model equations describing the dusty plasma medium in 1D are ( [6,7] and references therein) where n and u are the normalized density and velocity of the dust fluid component, Z d is the dust charge number, σ = T i /T e is the ratio of the electron and ion temperatures, T d is the dust temperature, µ e = n e0 /Z d n d0 and µ i = n i0 /Z d n d0 . The subscripts d, e, i refer to dust, electron and ion species, respectively, and the subscript '0' refers to equilibrium values. The various quantities are normalized as follows: φ = eφ/k B T i , n = n d /n d0 , n i = n i /n i0 , n e = n e /n e0 , t = ω pd t, u = u d /L D ω pd and x = x/L D . Here ω pd = 4π(Z d e) 2 n d0 /M d 1/2 is the dust plasma frequency with e, M d , k B denoting the electronic charge, dust mass and Boltzmann constant, respectively. Further, L D = k B T i /4π Z d n d0 e 2 1/2 is a convenient length scale that makes the coefficient of the ∂φ/∂ x term in (2) unity. Note that the charge neutrality condition requires the equilibrium densities to satisfy the condition n i0 = Z d n d0 + n e0 . This standard model for DAW assumes that the electrons and ions obey a Boltzmann distribution and the dust pressure is modelled by a simple equation of state P = k B nT d . Equations (1) and (2) are numerically solved by the flux corrected scheme of Boris et al [24]. At each time step the scalar potential φ is determined from the Poisson equation (3). The latter is solved by employing a successive relaxation scheme. Our numerical code has been appropriately benchmarked in the small-amplitude limit to accurately agree with the linear dispersion relation of the DAW.

Results and discussion
When the amplitude of the initial perturbation is gradually increased, we observe the excitation and evolution of the usual KdV soliton solutions that have been analytically predicted before. For a further increase in the amplitude we observe the excitation of both a smooth soliton and a sharply peaked soliton. This is shown in figure 1 where an initial large-amplitude Gaussian pulse (figure 1(a)) is seen to split into two oppositely propagating pulses (figure 1(b)) and as we follow the evolution of the right propagating pulse it is seen to separate into a smooth small-amplitude soliton and a large-amplitude sharply peaked solution (figures 1(c) and (d) show this at various stages in time). In figure 2 we have shown the spatial structures of the density n (solid line), the velocity u (dashed dot line) and the potential φ (dash line) of the soliton solution. The plots clearly reveal the development of a singular structure. For these simulations the dust species was assumed to be cold, i.e. α = v 2 thd = 0.0. The zero-temperature case is, however, a mathematical idealization. Keeping this in view we have also carried out simulations of an initial Gaussian pulse for the case of finite dust temperature. We show the results of the finite-temperature case in figure 3. In this case the singular characteristics, namely a cusp formation in density profile, can also be observed clearly. We note from our simulations that the profiles and the maximum value of the fields u and φ of the sharp structures that finally form are insensitive to any refinements in the grid resolution for both the zero and the finite-temperature cases. The density profile and its maximum value are, however, insensitive to grid resolution only for the finite-temperature case. When the dust temperature is chosen as zero the maximum of the density is dependent on the choice of the grid resolution. We would show that in this case one expects the theoretical value of the density to shoot off to infinity, which obviously can never be captured, no matter how much the grid is refined. The sharply peaked soliton and the regular soliton move apart in time as their propagation velocities are different because of the difference in their respective amplitudes. We also observe from our simulations that the maximum value of the dust fluid speed u = u sim max (which occurs at the cusp point) is close to (β − v thd ), where β is the sharply peaked soliton propagation speed and v thd is the dust thermal velocity. This is seen from data gathered and consolidated from several simulations and summarized in table 1.
To gain a better understanding of our numerical simulation results we carry out a stationary wave analysis of the full set of equations (1)- (3). Assuming that all dependent variables are functions of ξ = x − βt only, we can reduce the full set of equations to the following ordinary differential equation (ODE) in n as where Here α = v 2 thd is a parameter. The analytical estimate for β is obtained from the upper bound condition for the formation of the localized structure of the pseudopotential V (n) [25]. From the expression it is clear that V (n) blows up when the denominator of the first bracket goes to zero. This occurs when n = n max = β/v thd and here clearly the first derivative in density ∂n/∂ξ blows up from equation (4). Also from the continuity equation we have n = β/(β − u); thus at this point we would have u = u max = β − v thd , as observed in our simulations studies. Similarly, an expression for φ max can be obtained from equation (6). The analytical values corresponding to these expressions for n max , u max and φ max have been tabulated for various cases in table 1. We also show the maximum values of the density, velocity and potential fields observed in the simulation, for the sharp pulse that eventually forms spontaneously from a given large-amplitude Gaussian initial condition. We observe that for u max and φ max there is good agreement between the analytical estimation and observations obtained from simulations. The estimates and the observed values of the density in simulations show slight differences. For the cold dust (T d = 0), the density can be infinite for the cusp solutions, which as stated earlier cannot be captured in simulations. In any case the simulations for T d = 0 is only a mathematical idealization. We also wish to state that the choice of σ = T i /T e = 1.0 made in our simulations, although unrealistic, does not influence the cusp behaviour, which is the main theme of this paper. The cusp occurs at the point where V (n) blows up, i.e. when the denominator of the first factor in equation (5) (which does not depend on σ ) goes to zero. Furthermore, basically the realistic condition of T i T e would simply make ion shielding more important. Charge density on the dust changes the electron density, influencing the electron Debye length, and again enters through shielding effects only.
We also choose to consider different initial conditions for our simulations. For instance, we have chosen the fields corresponding to the regular large-amplitude solitons and the cusp solutions obtained from equation (4) as initial conditions. The regular soliton solutions and the cusp structure for a particular set of parameters have been shown in figures 4(a) and (b), respectively. In figure 5(a), the evolution of regular soliton structure of figure 4(a) is shown. As expected the undistorted propagation of the soliton is captured well by our simulation. The choice of the cusp structure of figure 4(b) as our initial condition again leads to a stable propagation. The evolution of the profiles of density n and velocity u are shown in figures 5(b) and (c), respectively. It should be noted that we had specifically chosen to show here the pathological case of T d = 0 for which ideally the density solution should blow off to infinity. Clearly, the density profile for this case cannot be accurately represented in the numerical solution chosen as the initial condition and/or its simulation. The evolution, therefore, shows fluctuations in the amplitude of density. However, we would like to point out that even for this case, the velocity evolution shows a very stable propagation, as has been demonstrated in the plots of figure 5(c). For those simulations for which the dust temperature is finite even the evolution of density shows no perceptible change in its amplitude. It appears that when a sharply peaked structure forms, the gradients in density and velocity become large close to  Figure 7. The evolution for sinusoidal perturbation for the case of the finitetemperature dusty plasma system α = 0.1 has been shown. The various panels (a)-(d) correspond to t = 0, t = 1.227, t = 3.067 and t = 60.132, respectively. The other parameters are the same as that of figure 6. In this case, the development of shocks initially turning into cusp structure can also be clearly observed.
the peak. These are smoothed by the numerical integration schemes at the grid scale which may be interpreted as viscosity/hyperviscosity operating on the soliton peak. This keeps the amplitude dithering close to the wave breaking amplitude. In reality, these numerical effects will be replaced by collisional effects or wave-particle interaction effects in collisionless plasmas.
To further test the accessibility and excitation conditions of the cusp solitons, we have next changed the initial conditions to be in the form of a sinusoidal perturbation. This is also close to the experimental conditions of [23]. Again for very small amplitudes the wave train remains sinusoidal and propagates with a phase velocity given by the linear dispersion relation of the DAW. When the amplitude of the initial sinusoidal perturbation is chosen to be large, we observe that in the beginning the sinusoidal perturbations get steepened to form a periodic train of shocks. These shock structures then eventually evolve into a chain of cusp-like sharply peaked structures which are stable and survive over hundreds of DAW periods. These stages of evolution have been illustrated in figures 6 and 7 for a cold and finite-temperature dusty plasma system, respectively. The evolutionary stages of a sinusoidal perturbation leading to the emergence of a train of sharply peaked structures bear a remarkable resemblance to the experimental findings of Teng et al [23], where a self-excited DAW is found to grow in amplitude till it approaches the wave breaking condition. The measured dust density profiles in that stage appear as a series of propagating cusp profiles. These steepened structures are seen to then significantly influence the dust micro-dynamics, leading to particle trapping and disordered motion that brings about a phase transition from a liquid state of the dust fluid to a gaseous state. Our simple model simulations, based on a fluid model, cannot reproduce the wave particle dynamics observed in the experiment but does seem to capture well the emergence of propagating cusp structures that are close to wave breaking conditions. One further shortcoming of our model is that it is valid in the weakly coupled regime of a dusty plasma where correlation effects arising from strong dust-dust interactions have been neglected. However, the basic nonlinear wave propagation phenomena displayed by the present fluid model may not be significantly altered by these correlations in the framework of a hydrodynamic description.

Summary and conclusion
To conclude, we have presented strong evidence for the existence of cusp solitons in a dusty plasma medium on the basis of 1D numerical solutions using a fluid model. These nonlinear states of the DAW are shown to be long-lived structures and are accessible from a variety of initial conditions. The conditions for their onset and existence are obtained from a travelling wave analysis of the evolution equations and shown to be in good agreement with the numerical simulations. Our results assume significance in the light of the recent observation of similar structures in laboratory experiments and their potential role in influencing the micro-state of dusty plasmas.