Modelling nonlinear electrostatic oscillations in plasmas

The nonlinear 1-D plasma electrostatic oscillation is formulated in an analytic framework that allows closed-form analytic solutions along the characteristics, and solved numerically in configuration space. Additionally, a novel iterative analytical form for the finite-amplitude oscillation solution is derived, which compares favourably with the other two techniques. A fresh insight into the evolution of the oscillation is gained, including defining the least achievable density in the nonlinear oscillation as half of the equilibrium value, and relating the associated maximum density achievable in terms of that minimum.


I. INTRODUCTION
The electrostatic plasma oscillation is arguably the defining characteristic of that medium: the unique balance between conduction and particle currents that produces the distinctive "ring" that can only happen in a plasma.][3][4] The main motivation for revisiting this classic problem is the unique context of cold plasma oscillations in pulsar crusts, where the magnetic field strength is so high that the associated material compression ensures that the positive ions can truly be considered to be stationary, while the abundant free electrons are constrained by the Landau levels to have momenta entirely aligned with the internal field direction. 5The wave-breaking of such oscillations is alleged to eject electrons from the metallic crust into the pulsar atmosphere immediately above it, thus, populating the environment with energetic electrons, the radiation from which can then create the electron-positron pairs, which are the defining characteristic of the pulsar envelope.
Given that nonlinear effects are crucial in this context, we present a reworked analytical framework in which we demonstrate a novel insight into the solution, motivated by the method of characteristics but extended by incorporating recursive solutions to reveal the full nonlinear evolution in closed form.While there have been several reformulations of the wavebreaking problem in both the classical and relativistic limits 4,[6][7][8][9] there is a focus on either driven oscillations (forced by the imposed laser pulse electric field profile) or by beam-plasma interactions (in which the beam profile is assumed); the analytical techniques used involve Lagrangian coordinates and sometimes phase mixing, where an imposed inhomogeneity in the underlying density is forced either as an analytical formula or as spatial fourier series.The approach in this paper is to formulate the equations using the theory of characteristics, and solve directly for the electric and velocity field, and the density, along the characteristic without making assumptions about the nature of the initial wave profile.The natural solutions that arise on the characteristics are explored fully, analytically and numerically.These solutions, arising from independent methods, are mapped back to configuration space to confirm the validity and accuracy of our analytical formulation, and together these approaches allow new insight into how such processes evolve.While there is a clear motivation in the pulsar context for revisiting this description, the analysis holds for all plasmas in which the approximations are sufficiently valid.

II. MODEL EQUATIONS
Consider a cold electron plasma, in the absence of a magnetic force.Let the uniform number density of positive ions be denoted n 0 , and for simplicity, assume that the equilibrium number density of electrons is also n 0 (that is, assume that the atoms are only singly ionized).Let the total number density of electrons be where the single spatial coordinate is x, and time is denoted by t.Further, let the advective derivative be denoted by D, so that where u is the (fluid) electron plasma velocity component in the x-direction.The relevant fluid equations for the electrons are then as follows: where E is the electric field, m is the electron rest mass, and 0 denotes @/@x.Noting that Du 0 ¼ ðDuÞ 0 À u 02 ; (6) a homogeneous equation for the evolving density can be formed where Eq. ( 3) has been used to eliminate u 0 , and where x p ¼ ½n 0 e 2 =ð 0 mÞ 1=2 is the usual plasma frequency, a constant here.
For the small-amplitude disturbance, the density satisfies the usual plasma oscillation but the role of density gradients in the non-linear evolution is clear, particularly if the density evolution equation is written in the form where x2 p ¼ ð1 þ ñ=n 0 Þx 2 p is the square of the effective plasma frequency, which tends to zero as ñ !Àn 0 ; this same limit produces a singularity in the right-hand side of Eq. (9).
Since the focus of this article is the plasma oscillation, the electrostatic condition can be assumed where we have set the constants e, m, and 0 to unity, for clarity (and so making x 2 p ¼ n 0 ).Combining Eq. ( 5) with Eq. ( 10) yields Using the same convention on constants in Eq. ( 4), we have a closed set of operator equations governing the electrostatic dynamics Applying D to the equations for u and E yields simple coupled equations in the operator notation Note that u and E are interchangeable in Eq. ( 14); moreover, there are no singularities in these equations.It is also notable that the velocity and electric field engage in oscillations at the original plasma frequency, along the characteristics; this is in keeping with other authors; 3,6 however, the plasma density equation is non-linear along the characteristic, as can be seen from Eq. ( 9), and its spatio-temporal behaviour is less straightforward than either of the electric or velocity fields.The density equation ( 9) does, in fact, have an explicit singularity, though only when the total density approaches zero, which is physically problematic in any event, and is in fact impossible, as the analysis in Sec.III will show.

III. THE DENSITY CALCULATION ON THE CHARACTERISTIC
The density fluctuation ñ can be calculated via the electric field gradient (in configuration space), or can be evaluated directly on the characteristic from the homogeneous equation given in Eq. ( 9).The beauty of the latter approach is that the extremal values of the density are identical in both co-ordinate systems, but it is easier to solve directly on the characteristic, even if the precise path of the characteristic (in configuration space) remains unspecified.Let y ¼ n/n 0 , and set x2 p ¼ 1. Taking the coordinate along the characteristic to be z, so that D can be replaced everywhere by d/dz, the normalised density equation can be written in the form which has the solution where Hence the normalised density solution along the characteristic takes the form for arbitrary constants c 1 , c 2 .Given the general condition that for all y.Hence in particular, at a critical turning point, where y 0 ¼ 0; y ¼ y c we have Consider the possible extremal values of the density profile: If n > 1, then both possible values for y c are negative; hence, the physically correct constraint must be n < 1, and consequently y > -1/2 everywhere along the characteristic.Since mapping the solution into configuration space cannot change the values of the solution, but only its gradient, then this analysis establishes as physically correct that the minimum plasma density is half of the background density, as noted by Davidson and Schram, 2 (see Fig. 1 for examples of the profile for various amplitudes), and apparent in a variety of subsequent treatments (for example, see Fig. 2 of Ref. 10, where the density profile is very similar to that given here in Fig. 1, and also suggested by Figs.3-6 in  Ref. 11).It is simple to infer that the extrema of the density profile given by Eq. ( 16) imply that the normalised density y must lie in the interval where À1/2 < y m < 0 is the minimum density depression produced by the disturbance.These extremal values have been calculated without requiring any prior assumption about the shape of the characteristic in configuration space: these results are, therefore, general.The extra insight into the extrema of density afforded by direct solution on the characteristic, rather than grappling with the evident nonlinearity of the Lagrangian transformation into configuration space, is an important novel aspect of this work.

IV. ITERATIVE ANALYTICAL SOLUTION IN CONFIGURATION SPACE
We now turn to solving the whole system in configuration space, by using the method of characteristics. 12,13The system to be solved is as follows: subject to the initial conditions uðx; 0Þ ¼ f ðxÞ; Eðx; 0Þ ¼ gðxÞ; for some initial spatial profiles f(x), g(x) of the disturbance.Analysing this system by the method of characteristics, we choose new variables n, s in place of x, and t, respectively.The characteristic equations are then where the initial conditions for both characteristics and dependent variables are xðn; 0Þ ¼ n; tðn; 0Þ ¼ 0; uðn; 0Þ ¼ f ðnÞ, and Eðn; 0Þ ¼ gðnÞ, and subscript s denotes @/@s.
Eliminating E from Eq. ( 26) yields which has general solution uðn; sÞ in which A and B are independent of s.A simple solution consistent with initial conditions is given by A ¼ f(n), B ¼ 0 (this minimises the algebra, but can be generalised later) so that the solutions for u and E along the characteristics can be given as uðn; sÞ ¼ f ðnÞ cosðx p sÞ; Eðn; sÞ ¼ Àx p f ðnÞ sinðx p sÞ; where The equations of the characteristic curves themselves can then be expressed as 2,14,15 xðn; sÞ In contrast to Section III, we have an expression for the shape of the characteristics in configuration space, which we can explore.

A. A full recursive analysis
For illustrative purposes, we can consider the simple case f ðnÞ ¼ u 0 sinðknÞ, which corresponds to a simple sinusoidal spatial disturbance in the velocity field at t ¼ 0. The characteristic equation is then where l ¼ u 0 sinðx p tÞ=x p , we can proceed as follows.Using the trigonometric product rule to simplify the representation for u, we can express the zeroth iterate u (0) as the case where n is simply replaced by x u ð0Þ ¼ ½u 0 cosðkn À x p tÞ n¼x ¼ u 0 cosðkx À x p tÞ: (34) The next step in the recursion u (1) is to replace the first occurrence of n with the expression x À l sinðknÞ as follows: The procedure can be repeated to yield nested formulae to the nth degree: for example, the fourth recursion is given by These recursive expressions quickly become difficult to interpret analytically, but are very revealing graphically.The first four velocity recursions after the zeroth are plotted in Fig. 2, for the particular value of x p t ¼ p=2, with u 0 ¼ k ¼ 1 and l ¼ 1.The characteristic steepening of the solution near x ¼ p and x ¼ 3p, and to a lesser extent near x ¼ 0 and x ¼ 2p is consistent with the onset of wave breaking.This is clearest in Fig. 3, where the profile of the solution becomes markedly steeper as time progresses; the parameters are consistent with Fig. 2, having extracted the time t from parameter l ¼ ðu 0 =x p Þ sinðx p tÞ by setting u 0 /x p ¼ 1.2 and evaluating for various choices of x p t ¼ q.Additional insight comes from an analysis of the gradient (in configuration space coordinate x) of the recursive solutions.The zeroth iterate u (0) has derivative ku 0 times the envelope function that defines the overall shape of the solution.For the gradient of u (1) , the factor d 1 multiplying the envelope solution is ku 0 ð1 À klÞ cosðkxÞ, and so we have for the ith iterate d i , i ¼ 1,…, 4 (38) (41) Â sinðkl cosðkl cosðkl cosðkxÞ À kxÞ À kxÞ À kxÞ À 1Þ: (42) Given that the graphs show the largest gradient in the vicinity of kx ¼ p/2, consider the gradient multiplier d n of the nth iterate evaluated there where D n is the additional multiplicative factor over the linear case (that is, the zeroth recursion) for the gradient of the nth recursion.Note that the gradient only converges with recursion number if kl < 1; hence a critical condition for wave steepening is ku 0 /x p ! 1, at which point the gradients cannot converge.In reality, the recursive solution shows a phase drift away from the simple assumption that the maximum gradient occurs at the same position in the higher iterates as for the zeroth (namely, at q ¼ p/2); however, Eq. ( 43) is an adequate guide to the envelope of the gradient, as shown in Fig. 4.Under the envelope approximation, it is clear that the gradient will not converge with iteration number unless jklj ¼ jku 0 sinðx p tÞj=x p < 1, that is, ku 0 /x p < 1, in which case the limit of the maximum gradient is given by FIG. 2. Graphs of recursive iterations 0 to 4 of the velocity solution, for the particular choice of Since the electric field is just the time derivative of the velocity field (from Eq. ( 30)), then the maximum electric field gradient can be approximated by for the case where the recursion process converges.Note that the gradient of the electric field is directly proportional to the number density fluctuation, from Poisson's equation.
In any event, it is clear that u 0 ¼ x p /k is the limiting velocity amplitude (since otherwise the gradient is unbounded), with E ¼ x 2 p =k as the limiting electric field magnitude-in agreement with the published literature, which offers the maximum electric field amplitude as n 0 e=ðk 0 Þ ¼ m e x 2 p =ðkeÞ, in dimensioned units.The iterative method, thus provides straightforward (if cumbersome) analytical expressions for the longitudinal electric field, velocity and number density, the accuracy of which can be assessed in comparison with the full numerical solutions presented in Sec.IV, but it is useful to compare the profiles generated in this way with results in the literature, particularly, 10 Fig. 2; our analytic expressions that are given in Fig. 5, agree closely with the advanced simulations in Ref. 10.
Given that u 0 < x p /k means the velocity gradient multiplier D 1 converges to ð1 À ku 0 =x p Þ À1 , it is still possible to consider a practical limit on the maximum amplitude of the velocity disturbance before unacceptably steep gradients develop.If h ¼ ku 0 /x p ! 1, then the gradient does not converge with recursion number, and wavebreaking must follow, since the slope of the velocity approaches vertical.However, the threshold value for the velocity gradient above which wave steepening is deemed to approach breaking point may well be reached for values of h < 1.For example, if h ¼ 0.8, then D 1 ¼ 5. Hence the practical limit to the magnitude u 0 of the velocity before wave-breaking sets in may be achieved with u 0 < x p /k.

V. DIRECT NUMERICAL SOLUTION OF THE CHARACTERISTIC EQUATIONS
In order to compare the iterated analytical solutions with the full, unapproximated ones, the differential equations for the characteristics themselves, and for the physical variables along those characteristics, can be solved numerically, with the solutions mapped back into the configuration space for comparison.The numerical packages in the computer algebra system Maxima were used in accordance with the characteristic solution method given by Whitham. 15umerical simulations of the velocity and electric fields are presented in this section, for the case x p ¼ 2, and for an FIG. 5. Plots of the longitudinal electric field (E), velocity (u) and electron number density (n) for the nonlinear case of Ref. 10, Fig. 2. FIG. 6. Surface plot of the velocity field from direct numerical simulation using the characteristic equations, for the case where x p ¼ 2, and the amplitude of the initial disturbance is B ¼ 0.5.Note that the surface is a spline interpolation to a uniform x, t grid from the characteristic curves.FIG. 7. The set of characteristics for the case B ¼ 1.5; the solutions for the velocity and electric fields are computed along these curves, and then interpolated onto a uniform x, t-grid.Note that the characteristics become very close together in the region around x ¼ 2, t ¼ 1, showing that gradients are becoming very steep here.FIG. 8. Left plot: Surface plot of the velocity field from direct numerical simulation using the characteristic equations, for the case where x p ¼ 2, and the amplitude of the initial disturbance is B ¼ 1.5.Note that the surface is a spline interpolation to a uniform x, t grid from the characteristic curves.Right plot: The recursive solution for the same parameters as in the left hand picture.FIG. 9. Left plot: Surface plot of the electric field from direct numerical simulation using the characteristic equations, for the case where x p ¼ 2, and the amplitude of the initial disturbance is B ¼ 1.5.Note that the surface is a spline interpolation to a uniform x, t grid from the characteristic curves.Right plot: The recursive solution for the same parameters as in the left hand picture.FIG. 10.Left: The normalised density plot in configuration space for the case B ¼ 1.5 and x p ¼ 2, calculated via numerical simulation along the characteristics of the electric field, which is first mapped onto the configuration space grid and smoothed by cubic splines before being differenced to yield the derivative (and hence the number density).Notice that the minimum density saturates at half of the equilibrium value, consistent with the analysis in Section III.Right: The normalised density plot in configuration space for the same parameters, but calculated using the derivative of the analytical iterative formula for the electric field; there is no smoothing applied.Although the correlation between the two solutions is not exact, the essential features are recovered.The iterative formulae are limited in Fourier harmonics, and this becomes more evident when differentiated.initial disturbance in the electric field of the form Eðx; t ¼ 0Þ ¼ B sinðxÞ, with the accompanying uðx; t ¼ 0Þ ¼ B cosðxÞ.The calculation of the velocity surface is shown in Fig. 6.Increasing the amplitude causes the evolution of sharp features: for example, changing from B ¼ 0.5 to B ¼ 1.5 produces a significantly distorted set of characteristic curves, shown in Fig. 7.The associated velocity and electric field surfaces are given in Figs. 8 and 9; in each case, the numerically integrated solution is shown on the left, and the analytic solution (for the same parameters) on the right.Given that the numerical solution is smoothed and interpolated from the characteristics, and that the analytic expressions are truncated at a finite iteration, the correspondence between the solution forms is encouraging.Indeed, the direct calculation of the density, which is given by the gradient of the electric field, shows as similar agreement between numerical and analytical forms, even under additional differencing operations; these are shown in Fig. 10.

VI. CONCLUSIONS
There are a number of novel aspects presented here.The operator form of the equations allows the direct, analytic solution of the density along the characteristic curve.The analytical solution Eq. ( 16) of the oscillation equations in characteristic form shows why the plasma density cannot fall below one half of its initial value anywhere: a restriction attributed initially to an artefact of the Lagrangian method 2 to avoid singularities, but now demonstrated to be an asymptotic limit to the density depression caused by an oscillation.The corollary is that the maximum value of the normalised density can be expressed in terms of its minimum, with an indicative maximum spatial gradient n 0 in configuration space given by n 0 $ 2n 0 kq 1 À q 1 À 2q ; where q ¼ jy m j is the absolute value of the normalised density minimum.Plotting Eq. ( 46) to give the maximum density gradient as a function of the minimum density produces the graph in Fig. 11; comparing the values here with those apparent in Fig. 5 for the density show the utility of the iterative analysis.Despite the underlying simplicity in the formalism, the iterative solution recovers the intrinsic behaviour of the full numerical solution, remarkably well, via the method of characteristics; Figures 8-10 show the similarities, albeit the limitations of a finite number of harmonics is evident in the electric field and density analytic solutions.
The iterative solution is a powerful, closed-form approximation to the full numerical solution for finite amplitude perturbations.No continuum solution technique can carry an analytical solution through the breaking point itself, since the equations change their character at such a point, but the 4th order iterated solution is a very good approximation to the full solution for both the velocity and the electric field.
At this point, we can return to the original motivation: the emission of electrons from the pulsar crust.From Eq. (45) it is clear that strong gradients in electric field can produce electrons with an energy gain DE of magnitude where the limiting wave velocity amplitude is u 0 ¼ x/k, at which point wave breaking occurs.Hence, the non-linear self-field in the oscillation dynamics would appear to be capable of delivering the requisite Fermi Energy 5 to enable electron escape from the crust near wave-breaking conditions.Finally, we note that the method of characteristics used here for numerical solution is formally a more general technique than the Lagrangian methods used in earlier literature to make progress: the latter are restricted to a physical interpretation of the Lagrangian trajectory, namely, that this constitutes a velocity streamline for the fluid, and so must remain unidirectional.There is no such constraint on the more abstract characteristics, since higher-order systems can have multiple characteristics passing through the same points in space and time. 15Indeed, we are preparing a more general analysis of the oscillation problem that relaxes the electrostatic condition imposed here, and gives rise to precisely such a scenario.

FIG. 1 .
FIG.1.Examples of the normalised density profile along the characteristic, for various choices of amplitude c 1 , as given in Eq. (17), with c 2 ¼ 0. Note the localised strong profile steepening at for z < p, and the flattening for p < z < 2p, where the solution does not decrease below À1/2.
(4).3.Graphs of u(4)as a function of x for various values of q ¼ xt showing the onset of wave steepening, keeping u 0 ¼ k ¼ 1, and u 0 /x p ¼ 1.0.FIG.4.Contour plot of @u (4) /@x as a function of kx (horizontal axis) and q ¼ x p t (vertical axis) with kl ¼ 1.2, and u 0 /x ¼ 1, showing the gradient extrema occurring around kx % p/2.