Weakly nonlinear ion waves in striated electron temperatures

The existence of low-frequency waveguide modes of electrostatic ion acoustic waves is demonstrated in magnetized plasmas for cases where the electron temperature is striated along magnetic field lines. For low frequencies, the temperature striation acts as waveguide that supports a trapped mode. For conditions where the ion cyclotron frequency is below the ion plasma frequency we find a dispersion relation having also a radiative frequency band, where waves can escape from the striation. Arguments for the formation and propagation of an equivalent of electrostatic shocks are presented and demonstrated numerically for these conditions. The shock represents here a balance between an external energy input maintained by ion injection and a dissipation mechanism in the form of energy leakage of the harmonics generated by nonlinear wave steepening. This is a reversible form for energy loss that can replace the time-irreversible losses in a standard Burgers equation.


I. INTRODUCTION
We study the propagation of weakly nonlinear lowfrequency ion waves in magnetized inhomogeneous plasmas.The ion temperature is homogeneous, but the electron temperature T e varies in the direction perpendicular to an externally imposed homogeneous magnetic field, so T e is locally enhanced in a magnetic flux tube.Conditions like these occur often in nature for plasmas out of equilibrium [1] and can also be produced in laboratory plasmas [2].Wave propagation in such conditions has been studied experimentally [3] and numerically [4].
Low-frequency ion waves have different characteristics for the cases where ci / pi < 1 and ci / pi > 1, with ci and pi being the ion cyclotron and the ion plasma frequencies, respectively.The present study is concerned with the case where ci < pi .For low-frequency waves we can have an interesting situation where the fundamental mode has a frequency below the ion cyclotron frequency, while some or all of its harmonics are above.We demonstrate by analytical and numerical methods how these conditions have importance for weakly nonlinear low-frequency ion waves.For low-frequency electrostatic plasma waves in fluid models it is often argued that the Korteweg-de Vries (KdV) equation is the most relevant one for describing weakly nonlinear waves [5].For the particular problem analyzed in the present work we find that the basic physical properties are best described by a model that includes the basic elements of the Burgers equation [6], but incorporating a time-reversible process for the effective dissipation mechanism.
The numerical part of the study is based on a particle-in-cell (PIC) code and is thus capable of resolving also kinetic plasma effects.Since the analysis is based on a fluid model, we use large electron-ion temperature ratios T e /T i 1 in order to Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
make the results consistent so that, for instance, linear Landau damping by ions can be ignored.These high-temperature ratios are seldom found in nature but can be obtained in laboratory experiments.While the linear fluid dispersion relation is well known for all parameter combinations, comparably less attention has been given to the phase and group velocities derived from it.These are central for the analysis presented in this work and therefore described in the relevant detail.
Sections II and III give the basic linear results.Sections IV and V extend the linear analysis with results to give arguments for a generalized Burgers equation.The analysis distinguishes the spatial variation in the directions parallel and perpendicular to the magnetic field B and thereby the axis of the electron temperature striation.We systematically keep only the dominant terms in series expansions of terms that depend on the B-perpendicular variable.The so-called Fubini solution is mentioned here explicitly for later reference.Results of numerical simulations are presented in Sec.VI, for which preliminary results were reported previously [4,7], while Sec.VII presents a summary and our conclusions.Appendix A discusses a particularly simple top-hat model for an electron temperature enhancement, where relatively simple analytical results can be obtained for the linear wave dynamics.Appendix B discusses a simplified model accounting for the wave radiation pattern escaping from the enhanced electron temperature striation.

II. HOMOGENEOUS MAGNETIZED PLASMAS
First we consider small-amplitude low-frequency ion waves in a homogeneous magnetized plasma.The basic fluid equations for this case are the linearized ion continuity equation and the momentum equation with the assumption of a simple equilibrium state, having a constant plasma density n 0 and homogeneous magnetic fields B 0 .These basic equations then reduce to the linear continuity equation for compressional perturbations and the linear ion momentum equation of the form where M is the ion mass and T i is the ion temperature, taking the ion charge to be q, allowing for multicharged ions q = Ne.The adiabatic exponent for the ion dynamics is γ i ≡ C p /C V in standard notation.For most plasma applications we can use the ideal gas value γ i = 5/3.The electrons are assumed to be locally Boltzmann distributed at all times.Physically, we assume that the electrons can flow along magnetic field lines, from wave crest to trough, to maintain this equilibrium.In the case in which the wave vector is exactly perpendicular to B 0 , this flow is not possible, so we have to exclude a small region of the dispersion diagram around k z = 0 from the analysis, but this is a minor restriction, which will not be elaborated here.Its resolution requires the electron inertia to be retained in the electron momentum equation.In the following we will explicitly assume that we are dealing with electrostatic waves and let the magnetic field be constant, as appropriate for low-β plasmas.
We use the ideal gas law for the ion equation of state p i = nT i and let a tilde denote perturbations of equilibrium.Boltzmann's constant is included in the temperature to simplify the notation.We also assume quasineutrality ñe ≈ ñi = ñ, i.e., considering perturbations where all wavelengths are much longer than the Debye length λ De .From the electron momentum equation we obtain ñ/n 0 = e ψ/T e for isothermal electrons, ignoring electron inertia.Introducing C 2 s ≡ (T e + γ i T i )/M and using plane-wave assumptions such as ñ = n exp(ik • r − iωt), we find after standard calculations the well-known dispersion relation of the form In the limit k z ≈ 0 we find the dispersion relation for electrostatic ion cyclotron waves The effective ion Larmor radius a i ≡ C s / ci enters as a convenient length scale indicating the wavelength where ion cyclotron waves become strongly dispersive.The limit k z → 0 has to be interpreted in a restrictive sense, as already mentioned, since the phase velocity component ω/k z parallel to B 0 has to be smaller than the electron thermal velocity u T e , i.e., |ω/k z | u T e , so that the electrons can maintain a local isothermal Boltzmann equilibrium by flowing along magnetic-field lines.In the limit k ⊥ → 0, the dispersion relation reduces to (ω 2 − 2 ci )(ω 2 − k 2 z C 2 s ) = 0 containing ion sound waves and the electrostatic ion cyclotron resonance.
More generally, the dispersion relation has two branches, one for ω 2 < 2 ci and one for ω 2 > 2 ci .The full dispersion relations are illustrated in Figs. 1 and 2. The group velocity has in general a direction different from the wave-vector direction.The wave properties of the two branches are very different, as illustrated best by the angle between the group velocity and the wave vector.For very low frequencies ω ci , these two vectors are almost perpendicular, while they are close to parallel when ω ci .The point (ω,k ⊥ ,k z ) = ( ci ,0,a i ) is singular, due to the merging of two dispersion branches.The figures also show the direction and relative magnitude of the group velocity in the (k ⊥ ,k z ) plane.We note that by the quasineutrality assumption, the ion plasma frequency does not enter the expressions.It therefore does not have any meaning to ask for the ratio ci / pi in Figs. 1 and 2. The assumption of quasineutrality in general restricts the applicability of results like (3) to relatively low frequencies as compared to the ion plasma frequency ω pi .For low frequencies ω ci , we find the dispersion relation ω = ±C s k z , giving the group-velocity vector u g = ±C s ẑ, where ẑ is the unit vector along B 0 .In this limit the dispersion relation of the wave is similar to that of shear Alfvén waves and can be considered as an electrostatic analog of this wave type.
To obtain a more generally applicable but also more complicated result, we have to include Poisson's equation explicitly, which gives the relation n/n 0 = (eψ/T e )[1 + (kλ De ) 2 ] for the perturbations in ion density, with the electron Debye length being λ De .By inspection we find that with the replacement in (3), we can relax the assumption of quasineutrality and find When ci < pi we find that the low-frequency branch (ω < ci ) is modified only slightly by the finite Debye length corrections.For ω ≈ pi there is a significant dispersive effect.

III. INHOMOGENEOUS MAGNETIZED PLASMAS
We now relax the assumption of homogeneous plasmas and allow the electron temperature to vary in the direction perpendicular to B. For the following analysis it is essential that the ion cyclotron frequency is taken to be smaller than the ion plasma frequency, i.e., ci < pi .The relevant frequencies are assumed to be so low that an inertialess electron component can be taken to be in local Boltzmann equilibrium at all times.We assume quasineutrality n e ≈ n i .For a linearized fluid model we readily derive a basic equation in the form where ψ is the electrostatic potential, related to the relative density perturbations as eψ/T e = η ≡ n 1 /n 0 .Since we are here only interested in cases where T e /T i 1, we took T i = 0 in (5).For the case T e = const, a linear dispersion relation is readily obtained from (5) by a Fourier transform with respect to t and x,z as discussed in Sec.II.
If we let T e = T e (x), with z being along the magnetic field B and x in the transverse direction, we can still Fourier transform with respect to time and along z.We denote the Fourier transformed normalized potentials by ψ.Normalizing frequencies and x positions so that ≡ ω/ ci and ξ ≡ x ci /C s , respectively, we readily obtain the expression where we introduced a normalized propagation speed γ through γ 2 ≡ (ω/k z ) 2 (M/T 0 ), where T 0 is a reference electron temperature that enters also the definition of the sound speed C s ≡ √ T 0 /M.Thus, γ measures the ratio between the phase velocity and a sound speed, so γ = const would correspond to the nondispersive wave propagation.In all cases analyzed here we use T 0 ≡ 1 2 [T e (0) + T e (|∞|)].The case with T e (−∞) = T e (∞) requires a special analysis not considered here.
The expression (6) has the form of an eigenvalue equation, with 1/γ 2 being the eigenvalue.Approximate methods [8] can be used to find qualitative solutions of (6) for the general case, while some simple models can be solved analytically (see Appendix A).In the following we will show results obtained numerically.To obtain the full space-time variation of the potential we have ψ(z,x,t) = A ψ(x) exp(ik z z − iωt), where we reintroduced x = ξC s / ci .The amplitude coefficient A has physical dimensions of potential, while we have the eigenfunction ψ(x) dimensionless as before.
We mentioned before the analogy between the dispersion relations of the present low-frequency electrostatic waves and shear Alfvén waves.For inhomogeneous conditions the analogy can be extended [4] by noting the existence of an electrostatic wave continuum, similar to the Alfvén continuum [9].

A. Small-amplitude waveguide modes
We present numerical solutions in Fig. 3 for < 1.For illustration we choose a temperature variation in the form T e (ξ )/T 0 = 1 − 1 2 D + D exp(−ξ 2 /W 2 ), with width W = 4 and D = 24/13 to give a significant spatial variation of the electron temperature.We find that in the low-frequency limit ω < ci the waves are confined to the electron temperature striation (here referred to as waveguide modes), corresponding to a discrete set of eigenfunctions ψm , with mode number m.For instance, for Gaussian variations of T e (x) we have that the mode number m corresponds to the number of zero crossing of ψm (x) [4].For the fully-three-dimensional problem we would FIG. 3. Illustrative examples of a numerical solutions of (6), assuming here an electron temperature profile given by T e (ξ )/T 0 = 1 − 1 2 D + D exp(−ξ 2 /W 2 ), with W = 4 and D = 24/13, as shown in the upper panel.The waveguide mode solutions are given for = 0, 1 4 , 1 2 , 3 4 .We have ψ0 (ξ = 0) = 1 in all cases.Note the large temperature ratio T e (ξ = 0)/T e (|ξ | → ∞) used here.In the upper panel T e (|ξ | → ∞) = 0, although its numerical value is small.have two indices ψkm corresponding to the two directions perpendicular to B. From ψm we can obtain the corresponding eigenmodes for the B-parallel velocity u z .The value of γ depends on and we find, for instance, for = 0, 1 4 , 1 2 , 3 4 , that γ 2 0 = 1.3942, 1.3805, 1.3311, and 1.201, i.e., a relatively weak variation of γ 0 with .In Fig. 4 we show a more detailed illustration for the variation of the eigenvalue.In Appendix A we give further illustrative results for a different and simpler boxlike electron temperature variation.This case is easier to solve, but has no direct relevance for a physically realizable situation.The toy model serves to illustrate, for instance, that there can be cases where only one waveguide mode can be found.
By inspection of Fig. 3 we note that the eigenfunctions change only little in spite of the large change in and the corresponding change in the eigenvalues is also small, as listed.Only for close to unity, say, around 0.9 or larger, do we see significant variations in ψ0 (x) and γ 0 .For shallow temperature FIG. 4. Illustration of the variation of the eigenvalue for varying .Note that for fixed plasma parameters we can approximate γ 0 by an average constant value to a good accuracy as long as < 1.Note the restricted range on the vertical axis.The relative variation of γ 2 0 is less than ±7%.6), for the same conditions as in Fig. 3.The antisymmetric modes are here obtained corresponding to the γ 1 values given in the figure.We have

FIG. 5. Examples of numerical solutions of (
variations and narrow temperature ducts, we have only the lowest-order mode ψ 0 (x).For γ smaller than the minimum value of T e (x)/T 0 we have a continuum of eigenvalues with corresponding eigenfunctions.
While Fig. 3 illustrates localized guided modes, we note that also free modes with a propagation component in the x direction exist.These reduce to obliquely propagating plane waves in the case in which the local electron temperature enhancement disappears.The localized waveguide modes disappear when the temperature striation disappears.By superposition of free modes we can construct an initially localized solution, just as a superposition of plane waves can give a pulse.This solution will however disperse with increasing time, while the waveguide mode will remain localized in the direction perpendicular to B at all times.Such a solution will be discussed in the following.
Physically, we can explain the waveguide modes by considering Fig. 1.If we want a wave packet to escape from the electron temperature striation, the wave vector has to point out of the striation.The group velocity associated with the wave packet will however be pointing along B 0 when ω < ci < pi and the wave packet will not leave the striation.
The lowest-order antisymmetric trapped modes are shown in Fig. 5 for the conditions used in Fig. 3. Symmetric and antisymmetric higher-order modes are illustrated in Fig. 6.Only the lowest-order symmetric mode ψ0 (ξ ) with mode number m = 0 is positive everywhere for all γ .This mode has also the highest propagation velocity (see also Appendix A).Physically, it can be argued that m = 0 is the fastest mode since it is maximized where the local sound speed is maximum and decays rapidly in the region where the local sound speed is reduced.For instance, the continuum modes or the higher-order m modes are different in this respect.
The foregoing analysis demonstrated that we have a discrete set of eigenvalues γ m (i.e., characteristic propagation velocities) with corresponding eigenfunctions ψ m .As an illustration we show in Fig. 7 the low-frequency limit of the dispersion relations for the modes given in Figs. 3, 5, and 6.There are here only five dispersion relations.The fastest mode corresponds to the mode number m = 0, which is also faster than the continuum modes.FIG. 6. Examples of higher-order modes (m = 2, 3, and 4) obtained by numerical solutions of ( 6), for the same temperature profile as in Fig. 3  Solution of ( 6) gives the eigenfunctions ψj (ξ ) for normalized potential.The corresponding eigenfunction for normalized density is T e (ξ ) , found by separating the variables in the linear expression for Boltzmann distributed electrons.This gives eφ(z,t) where the constant is dimensionless.Similarly, we find the normalized velocity eigenfunction to be FIG. 7. Low-frequency-long-wavelength limit of the dispersion relations for the trapped waveguide modes shown in Figs. 3, 5, and 6 obtained by the γ values given there.The larger the number of zero crossings of the eigenmode, the smaller the phase velocity.For illustration we inserted an open circle for some (ω 0 ,k z0 ) on one of the slower branches and illustrated how this wave can decay, with the decay products indicated by colored circles.The nonlinear wave steepening gives rise to higher harmonics, here illustrated by a filled circle for the second harmonic (2ω 0 ,2k z0 ). or ûj (ξ ) = ηj (ξ ) − d 2 ψj (ξ )/dξ 2 , obtained from the linear ion continuity equation where n(z,t)/n 0 = eψ(z,t)/T 0 is used and C s is derived by T 0 .Note that d 2 ψj (x)/dx 2 is given by the right-hand side of ( 6), there in dimensionless units.For the case in which d 2 ψj (ξ )/dξ 2 is small, e.g., for a wide electron temperature duct, we will have ûj (ξ ) ≈ ηj (ξ ).In Fig. 8 we illustrate the three lowest-order eigenfunctions, noting that they appear similar.

B. Free modes
For ω > ci the right-hand side of (6) changes sign, and the nature of the eigenmodes changes as well, to become free modes as shown in Fig. 9.If we let the electron temperature striation vanish to have a uniform T e , then the free modes degenerate to two obliquely propagating plane waves.As already mentioned, we have free modes also for ω < ci where γ 2 < min[T e (ξ )]/T 0 .The relations between the normalized phase velocities of the waveguide and free modes of propagation are discussed in Appendix A for the particular model for T e (ξ ) analyzed there.In general, we find that γ is larger for waveguide modes than for free modes FIG. 9. Examples of numerical solutions of (6), assuming here an electron temperature profile given as in Fig. 3.Here we show continuum solutions.The free modes are obtained here with = having the same mode number in the central region of the electron temperature striation.

IV. WEAKLY NONLINEAR WAVES FOR ω < ci
The linear analysis of the previous section will be generalized to include the lowest-order nonlinear effects.We will not use a strict reductive perturbation method, but rather build the analysis on the basic physical features of the problem.
To simplify the problem we assume a slab geometry for the electron temperature T e (x) with x ⊥ B ẑ and ignore variations with the spatial coordinate y.This model is relevant for numerical simulations.
By inspection of the dispersion relation for the waveguide modes (see, for instance, Fig. 7) we find that for any mode with mode number m 1, the wave number and frequency matching conditions for wave-number components along the z direction can be satisfied for wave decay [10] for some η a , η b , and η c with a = b = c.The present study will be restricted to the highest phase velocity m = 0 mode, where the three wave decay conditions cannot be fulfilled for forward propagation.
For waves excited at a boundary (as in our case), only forward decay is relevant.An example of a forward decay is illustrated in Fig. 7.For modes with m > 0 the nonlinearity can mix different m values to give wave decay where wave number and frequency matching conditions are fulfilled for forwardpropagating modes, i.e., k 0 = k 1 + k 2 and ω 0 = ω 1 + ω 2 as shown in Fig. 7.At later times, when the amplitudes of the (ω 1,2 ,k 1,2 ) modes reach significant amplitudes, also these will generate harmonics.By the construction it is readily seen that the highest phase velocity dispersion branch cannot have a decay that fulfills wave-number and frequency matching.The conditions are similar to those found for long-wavelength electrostatic modes in magnetized plasma waveguides with immobile ions [11,12].One difference is that for the latter problem only waveguide modes are found.For the highest phase-velocity mode with density eigenmode η0 (x) we expect wave steepening to be the dominant nonlinearity [6,[13][14][15].This process will be studied here.
In general, the nonlinear terms couple the various modes to give products of x modes.These can be expanded as, for instance, ηi (x) ηj (x) = q ζ qij ηq (x), where we assume that the set ηm is complete and orthogonal [11].The eigenmode continuum has to be included here, but for our applications the dominant contribution will come from the discrete set, so only this will be written explicitly.For the lowest-order nonlinear waveguide mode η0 (x) we have in particular For the Gaussian variations of T e (x) studied here, we will have ζ 000 > ζ j 00 for all j 1, since η0 is the only eigenfunction that is positive everywhere.
The electron profile can also be expressed in terms of an expansion of eigenfunctions ηm , for instance, in the combination ζ Tj ≡ ∞ −∞ ηj (ξ ) η0 (ξ )T e (x)dξ/T 0 .For a large class of relevant electron temperature profiles we can ignore all higher-order modes and retain only the ζ T 0 contribution.A model for the temperature variation is solved in Appendix A. The precise values of the constants ζ j 00 , etc., are of minor importance for the basic physics of the problem.

A. Weakly nonlinear simple waves
In the ion continuity equation we have basically two different velocity contributions entering.
The u ⊥ velocity component has contributions from the E × B/B 2 velocity and from the polarization drift [16,17] and we can write The ∇ψ × B drift does not contribute to the divergence in ( 7) since (∇ ⊥ ψ × B) • ∇ ⊥ n = 0 for the η 0 mode; this conclusion will be true also if, for instance, a cylindrical temperature enhancement in three dimensions is considered, with the lowest-order waveguide mode η00 , here with a double index.The polarization drift, which is a small correction, on the other hand gives rise to a B-transverse pulsation of the plasma in the electron temperature enhancement.Since these polarization fields affect the wave motion perpendicular to B, they have little influence on the wave steepening in the B 0 -parallel direction.Polarization fields will be important for B 0 -transverse particle accelerations, which can be important in different contexts [18,19].With these arguments we ignore the second term in (7) and retain the last nonlinear term.With these approximations, the ion continuity equation reduces to for the low-frequency dynamics.The ẑ component of the momentum equation for cold ions becomes in the same limit Both density and velocity perturbations depend on the two spatial variables, i.e., n = n(x,z,t) and u z = u z (x,z,t).Using the orthogonality and completeness of the set of eigenfunctions the product nu z in (7) can be decomposed to become a sum of z-varying products where the dominant term will correspond to the m = 0 mode shown in Fig. 3.The coefficient of the dominant term in the sum will be denoted by ζ 0 .
Following the arguments outlined before, we assume that initially and then at all subsequent times a disturbance is to retain its m = 0 structure.The problem will be to solve for the dynamic nonlinear equation for the z-varying magnetic-field-aligned structure.Accordingly, we assume n(x,z,t) → n 0 + n(z,t) η0 (x), ( 10) It turns out to be an advantage to use η0 as the basis here.With this notation we write (8) as η0 (x) where now η0 (x) cancels.To simplify the notation we wrote ζ 000 → ζ 0 .Similarly, we rewrite (9) as where ζ T originates from an expansion of η0 (ξ )T e (ξ )/T 0 .The expression for the Boltzmann distributed electrons is here written in the form eψ(z,t) which can be argued by using the series expansion of ln(1 + x) and approximating the dominant terms.
The basic set of equations ( 13)-( 15) represents the basis for the following analysis.The equations are not exact, but contain the most important physics.They will fail for studying higher-order modes by not including the possibility of wave decay of waveguide modes.For the lowest-order mode such a decay is not possible, as mentioned.

B. Simple wave approximation
We have that the space-time variation of the plasma density is through the variation of the potential, i.e., n = n(ψ(z,t)), by the assumption of Boltzmann distributed electrons together with the assumption of quasineutrality.We make now the simple wave assumption that the space-time dependence of potential is ψ = ψ(u z (z,t)).The potential is hereby assumed to depend on space and time implicitly through the velocity.The basic idea is identical to the one used for deriving simple waves in fluids and plasmas [6,17].The assumption is a natural extension of the linear theory, where n(z,t) and u z (z,t) are linearly proportional.With this assumption, we write the part of the ion continuity equation retained here as where we used and similarly for ∂n/∂z.
From the ion momentum equation we find The set of equations is closed by the assumption of isothermally Boltzmann distributed electrons.Combining ( 16) and ( 17) we find ) is known by (15).We finally obtain By the model of locally Boltzmann distributed isothermal electrons in (15) we find By this we obtain the well-known standard form of the simple wave equation [6,[13][14][15]] with u z = u z (z,t).

C. Fubini solution
The simple wave equation ( 20) can be applied for initialand boundary-value problems.The initial-value problem is usually studied in a frame of reference moving with the velocity C s0 .Assuming that the initial value is specified as u z (0,z) = G(z) in terms of some known function G(z), we can find the nonlinear solution in the form u z (z,t) = G(z − ζ 0 u z (z,t)t).This solution exhibits the well-known steepening and ultimate unphysical breaking of the waveform [6].The basic equation (20) can also be solved for a boundary-value problem.In this case we assume that u(0,t) = F (t) is given.This formulation of the problem has meaningful solutions only in the rest frame so we use the full expression (20).The solution for this case is u z (z,t) = F (t − z/(ζ 0 u z ± C s0 )).We can introduce an inverse function F −1 (t) to find In the standard reference case it is usually assumed that a signal in the form of a harmonic waveform with frequency ω 0 is excited at the boundary at z = 0.It is then found [14,20,21] that higher harmonics develop at increasing distances from the excitation as the waveform steepens with increasing z.This latter phenomenon is best illustrated by what is known as Fubini's solution, given by given here in terms of the normalized distance ρ ≡ z(ω 0 /C s0 )(u 0 /C s0 ), where the Bessel functions J n were introduced.The amplitude of the fundamental frequency decreases slowly with increasing z, while the amplitudes of the harmonics at the same time increase slowly (Bessel also seems to deserve credit for obtaining a result similar to the Fubini solution, although by using a different method [20]).This problem was originally analyzed for sound waves in neutral gases excited by a harmonically oscillating piston, where the displacement was expressed as with u 0 /C s0 1.The result given by ( 21) is approximate.It accounts qualitatively for the weakly nonlinear evolution of a variety of linearly nondispersive waves excited at a plane boundary as well as some plasma waves with acousticlike linear dispersion relations.

D. The KdV and Burgers equations
As shown by (20), the low-frequency ion waves in an electron temperature duct will steepen and eventually break like nonlinear one-dimensional simple waves.The limiting case with multivalued ion velocities is unphysical and mechanisms to stop the wave breaking have to be sought.In classical fluid mechanics we have two candidates for correction terms that arrest this unphysical wave breaking: viscosity (giving the Burgers equation) and dispersion (giving the KdV equation).The equations are here reproduced in a frame moving with the sound speed with a constant ν entering as a positive diffusion coefficient (in the original fluid version of the equation it is the kinematic viscosity [6]) and α determining the wave dispersion in the KdV equation.The Burgers equation has shock solutions, where dissipation balances wave steepening, while the KdV equation has the soliton solutions [17,22] where dispersion balances wave steepening.In classical fluid models of lowfrequency electrostatic plasma dynamics, e.g., ion sound waves, the KdV model has found applications [5].For the present case we note that in the wave-number range of interest, the linear wave has no dispersion.It will be argued that a generalization of the Burgers model is the most relevant one for the present problem.For later reference we here rewrite it in a spatially Fourier transformed versions as introducing the notation for the convolution used also later.To obtain the Fourier presentation in this form we first rewrite the nonlinear term as 1 2 ∂u 2 (z,t)/∂z.For small νk 2 z u z (k z ) the right-hand side of ( 23) is negligible and the wave steepening due to the nonlinear term [6] due to the nonlinear term [6] evolves uninhibited, thereby cascading wave energy to larger k z .At a some time we can no longer ignore νk 2 z u z (k z ) and the evolution saturates.In the case of a finite energy pulse, we have u z (k z ,t → ∞) → 0 for all k z .For a constant energy input, for instance at the boundary, an equilibrium state will develop in the form of a shock where energy is dissipated at the same rate that it is injected.A generalization of this latter scenario will be relevant for our study.

V. ELECTROSTATIC SHOCKS
The solutions of ( 20) have the well-known steepening of the initial condition.The characteristic time for wave breaking is approximately L/ max{u z (t = 0)}, where L is the characteristic scale length of the initial perturbation along B and max{u z (t = 0)} is the maximum value of the initial velocity perturbation.The model equation ( 20) assumes η0 (x) and γ 0 being used also when > 0, but this approximation is acceptable for at least 0 < 0.75, as can be seen from Fig. 3.For the basic ideas outlined in the present work, this restriction is of minor consequence.Polarization drifts become increasingly larger as is increased, but these do not affect the dynamics parallel to B, which is covered by (20).Also, a boundary-value problem can be addressed with the well-known Fubini solution for harmonic frequency generation with increasing distance from the source [14].
The formation of shocks as described by the Burgers equation [6] can be understood as a balance between the energy input by an external source (a moving piston, for instance) and viscous dissipation, where the kinematic viscosity coefficient is ν.In one spatial dimension, this particular nonlinear problem can be solved analytically by linearizing the equation through a Cole-Hopf transformation [6,17] and demonstrate, for instance, that the shock thickness varies as ∼ν/U with the basic parameters of the problem, where U is here the injection velocity.In principle, the Burgers equation can apply for any continuous viscous fluid media, also plasmas.Experiments performed in the strongly magnetized plasma of the Risø Q machine [23] demonstrated that for moderate electron to ion temperature ratios T e /T i , the strong ion Landau damping prohibited the formation of shock.For large temperature ratios, the ion Landau damping is reduced and there is a possibility for forming steady-state nonlinear shocklike forms, propagating at a constant speed [23,24].
In the present study we present a novel mechanism of an effective energy dissipation: selective radiation or leakage of short-wavelength ion sound waves.By numerical methods we also demonstrate that electrostatic shocks can form as a balance between these losses and the standard nonlinear wave steepening as described by the nonlinear term in the simple wave equation ∂u/∂t + u∂u/∂z = 0 [6,14].Studies in two spatial dimensions are sufficient for illustrating the basic ideas and the numerical analysis of the present paper is restricted to this case.
If we initialize the system with characteristic wavelengths λ corresponding to frequencies ω ci , i.e., λ C s / ci , the short-time evolution will be governed by ( 20) and we will have shorter and shorter scales developing as for the usual breaking of waves [6,14].This process is however arrested when the characteristic length scales become of the order of the effective ion Larmor radius C s / ci , where the modes become radiating, and are no longer confined to the waveguide.We propose a phenomenological expression for the process, best written in Fourier space in a frame moving with C s0 .We retain the dominant nonlinearity to be wave steepening, corresponding to harmonic generation for an initially plane wave as illustrated in Fig. 7.This part is accounted for by the nonlinear part of (20).In our model, this steepening can continue uninhibited until the harmonic wave numbers generated by the process exceed ci /C s0 .At this time the wave energy is radiated out of the temperature striation.Analytically, we can account for the transition wave number by introducing the Heaviside step function in a loss term in the wave-number representation, together with a characteristic time scale for the energy loss, where the symbol ⊗ denotes convolution as in ( 23), H is the Heaviside step function, and T characterizes the time it takes for the energy of the ω > ci waves to be lost from the waveguide.The right-hand side of ( 24) has the same role as −νk 2 z u z (k z ) in (23), while the left-hand sides of ( 23) and ( 24) are identical apart from a numerical constant.We have T = T (k z ), but it will depend also on parameters such as the waveguide width and the other plasma parameters.The term −H(|k z | − ci /C s0 )/T (k z ) has the role of −α in the Burgers equation (22).We expect that loss of wave energy is delayed when the width of the waveguide is increased, resulting in larger T .The wave steepening can thus continue for a longer time, giving decreasing shock widths T s for increasing T .The argument can be quantified by arguing that a characteristic time of energy flow in the direction perpendicular to B 0 is the width W of the electron temperature striation divided by the B 0 -perpendicular component of the group velocity taken at k ⊥ ∼ 1/W.As an order of magnitude we find that u g⊥ ≈ k ⊥ C 2 s / ci .Within the present model the waveform will steepen uninhibited until the shock width becomes of the order of C s / ci , at which time the harmonic frequencies will exceed ci to become radiative and the high-frequency wave energy will be lost from the waveguide.Consequently, the B 0 -parallel scale of the shock is controlled by a quantity referring to the B-perpendicular dynamics, which is an unexpected feature.

VI. NUMERICAL SIMULATIONS
The propagation of weakly nonlinear low-frequency electrostatic waves is studied numerically by a 2 1  2 -dimensional PIC code described elsewhere [4].The electron dynamics are modeled by an isothermally Boltzmann distributed fluid.The steady-state background electron temperature striation T e (x) is assumed to have a Gaussian variation so that T e (|∞|)/T i = 1 FIG.10.Illustration of the z variation for a propagating wave with the parameters ω = 0.1π pi and ci = pi /2 for different applied amplitudes.For this case the applied frequency is below ci , while some of the higher harmonics are above.The applied amplitude is increased in units of 0.1δn/n 0 .in the reference case, while T e (0)/T i > 1, where the ion temperature T i = const.Several values of T e (0)/T i were investigated.For the reference case, the spatial width of the electron striation is W = 28λ Di = 3.9C s0 / ci .Most of the theoretical analysis in the foregoing sections assumes quasineutrality, but the numerical code is general by retaining Poisson's equation, thereby resolving also the Debye length scales.

A. Small-amplitude waves
To illustrate some of the nonlinear features of the waves at ω < ci < pi we show in Fig. 10 a case where the wave frequency is chosen so that its fundamental is below the cyclotron frequency ci , while some of its harmonics are above.For small amplitudes the wave propagates as a guided mode in the electron temperature duct, but when the amplitude is increased it generates harmonics at frequencies exceeding ci .These waves are radiated out of the duct and the wave gets damped due to this loss of energy.Close to the excitation at the boundary we note that a Fubini-type [14,20] steepening waveform develops [see also (21)].For larger distances the waveforms become distorted.
For comparison we show in Fig. 11  The electron temperature enhancement is modeled by a localized Gaussian variation in the x direction.We used T e /T i = 50 at x = 0 so the ion Landau damping is negligible there.Only a part of the spatial simulation domain is shown.In both cases we show the spatial potential variation for a fixed time (x,z,τ = 80 −1 pi ).
as in Fig. 10.We note that the distortion of the waveform is considerably reduced.Although it is not a proof, the observation is consistent with a model where a large-amplitude waveform is distorted by leakage of higher harmonics.The steepening given analytically by ( 21) is now more pronounced at large distances from the excitation at z = 0. Results illustrating the waveguide and the free modes are shown in Fig. 12.The properties are clearly different and consistent with the interpretation given before.The waveguide mode remains inside the electron striation, while the highfrequency free modes are dispersing or leaking, consistent also with laboratory experimental results [3].The amplitudes are small in this simulation, with imposed density oscillations δn/n 0 = 0.1 at the boundary.The decrease of wave amplitude seen for ω > ci in the left panel of Fig. 12 is not a sign of damping in the sense of energy loss, but rather a spatial dispersion of wave energy emitted form a localized source.Several similar observations are reported for different wave types, some supported by a ray tracing analysis [25].This apparent damping will be found also in a fluid model without loss terms.The observed effective damping is caused by wave energy dispersing in space and dissipated by linear ion Landau damping outside the electron temperature striations where T e (|∞|)/T i = 1 in all cases considered.In order to emphasize the physical effects we discuss here, we consider only high-temperature ratios T e (0)/T i 25 in order to reduce the effects of linear and nonlinear ion Landau damping.Such high-temperature ratios (even as large as T e /T i = 100) can be obtained in discharge plasmas under laboratory conditions [26].

B. Shock solutions
In Fig. 13 we show an example of the formation and propagation of a shock.The density perturbation is chosen so that only the trapped waveguide mode is excited.The initial density perturbation has an error-function-type spatial variation in the B-field aligned direction.Ions are continuously injected at the boundary at y = 0 to maintain the energy input.We have T e /T i = 25, ci / pi = 1/2, and δn i /n 0 ≈ 0.22, where δn i refers to the actual detected density perturbation at the time where the shock is fully formed and not to the initial imposed perturbation.
The expressions (i)-(iv) that we have used for modeling the shock are summarized in Table I.Note that model (iii) scales differently from the others by being a solution to Burgers equation.It is possible to define a shock thickness T s = (n mx − n mn )/(dn/dx)(x s ) in all cases for comparison between models.The corresponding analytical expressions are also listed in Table I.When applied to our simulation results, they all give approximately (within the uncertainty) the same shock thickness T s .
When it comes to different fitting method, we have tried optimization methods like trust region reflective and optimize T s = 4s l the same objective function χ 2 and they give basically the same results.To illustrate the robustness of the shock fitting we use in Fig. 13 the model (iii) in Table I, while similar results were illustrated in [7], but there using model (ii) from Table I.
In Fig. 14 we show the shock thickness T s and the normalized shock velocity U for various combinations of the parameters T e /T i , ci / pi , and δn i /n 0 (see also the summary in Table II).The thickness T s and the shock position z s are determined by fitting various functions, listed in Table I, to the density variations at each time step.The results are robust  in the sense that they show no significant deviation between the various methods of analysis.The shock velocity U is obtained by following the shock position as function of time.
For small or moderate values of δn i /n 0 we find a close to linear relationship between δn i and U .As δn i → 0 (shown with the dashed line) we have U approaching the appropriate value for γ 0 , apart small corrections due to finite T i ignored in (6).Note here that γ 0 differs slightly from the value C s0 used for velocity normalization.For ci / pi = 1 2 we find that an approximate fivefold increase in δn i /n 0 corresponds to approximately 50% reduction T s .
To illustrate the radiative mode escaping from uide as the shock is steepening we show in Fig. 15 the (x,τ ) variation for a fixed z position.The channel of enhanced electron temperature is in the center of the figure.Here max T e = 75T i and min T e = 5T i , so ion sound waves are only weakly ion Landau damped outside the electron temperature enhancement.We have also studied cases with T e = T i outside the electron temperature striation.With these conditions the waves are strongly damped by the ion Landau damping, but an enhancement of wave potential is noticeable nonetheless.A simplified analysis explaining the basic features of Fig. 15 is presented in Appendix B.
In Fig. 16 we show ion phase space at a selected time t = 99 pi .For the high-temperature ratios used in this study, we have the ion sound speed to be much larger than the ion thermal velocity.Consequently, the number of ions reflected by the shock is negligible.
We studied also a weak-magnetic-field limit, with ci / pi = 1 4 and a narrow striation (see case 3 in Table II and Fig. 14).Finite ion Larmor radius effects that were ignored in deriving (5) will begin to be important here, since the striation is now only ∼10 ion Larmor radii wide.For this limit we find that the thickness T s is almost constant.Even weaker magnetic fields will require a fully kinetic analysis to account for the effects of collisionless ion viscosity [27].

VII. CONCLUSION
In the present study we analyzed the propagation of weakly nonlinear low-frequency waves in magnetized plasmas with inhomogeneous electron temperatures.The case where ci < pi was found to have interesting features concerning trapped and propagating wave modes.
We reported arguments for the formation of shocks in electron striations in magnetized plasmas when ci < pi and T e T i .By PIC simulations we demonstrated the formation and propagation of electrostatic shocks under these conditions for a wide combination of parameters and verified that the solutions have the basic characteristics of classical shocks, in spite of the system being dissipationless, i.e., fully time reversible.The decreasing shock thickness for increasing amplitudes is found also for classical shocks as described by the Burgers equation, but in our case the dissipation mechanism is leakage from the electron temperature striation of the short-scale lengths generated by the nonlinear wave steepening.This steeping due to the u∂u/∂z term in, for instance, (22) cascades wave energy to shorter wavelengths in the direction along B, where the energy at short wavelengths eventually gets lost by radiation out of the electron temperature striation as soon as k z > ci /C s .The nonlinear dynamics discussed in the present paper differs from other models applied for weakly nonlinear sound waves [28,29].It is often found that the basic nonlinear features of acoustic-type waves can be manifested in one spatial dimension, along a magnetic field for instance.In our case the weak magnetization is essential and the results depend on the model being developed in two or three spatial dimensions with one of them including the magnetic-field direction.
By numerical simulations, we demonstrated that the shock thickness T s scales inversely proportional to the Bperpendicular width of the electron temperature striation.This observation is consistent with expectations based on the model equation (24).The shock velocity increases approximately linearly with the strength of the plasma density perturbation.We verified that changes in the initial B 0 -parallel scale length of the initial density profile do not change the time-asymptotic saturated shock thickness T s for any δn i /n 0 .
For nonlinear waves described by a Korteweg-de Vries equation, a shocklike structure is followed by Airy-type ripples, originating from the dispersion term in the KdV equation.These ripples are absent in our results (see Figs. 13 and 16), although they develop in a later stage.Likewise, at the hightemperature ratio used here (with C s0 u th,i , the ion thermal velocity), we find no ions being reflected by the shock.A backward-propagating rarefaction wave is of no concern here.
The difference between the wave characteristics for ω < ci and ω > ci , respectively, has implications for observation in the ionospheric and magnetospheric plasmas with large electron-ion temperature ratios.For striated conditions, as discussed here, we would expect low-frequency waves to be localized to the striations, while for ω > ci the waves should be observed to be dispersed over all space.This feature should be observable also for moderate electron-ion temperature ratios.
The results may have general applications in magnetized plasmas with large temperature ratios T e T i where there are inhomogeneities in the electron temperature.In nature this could be due to electron precipitation.We find it more interesting, however, that the model introduces a concept for a time-reversible mechanism that has the effects of dissipation, in our case shock formation.We anticipate that similar effects can have relevance for waves in inhomogeneous fluid flows, for instance.The plasma conditions studied in the present work are physically realistic and also relevant for conditions found in nature, but they admittedly represent a special and limiting case.We find it interesting, however, that physically realizable conditions can be found where the reversible wave phenomenon has a role similar to dissipation for nonlinear wave propagation in a way that can be described solely by time-reversible fluid equations, i.e., by a model completely different from the Vlasov equation.
The assumption of Boltzmann distributed electrons used in this work is restrictive.For low-frequency phenomena it is easily justified, but it might fail for higher frequencies where a dynamic electron model is needed.The instability of drift waves in electron temperature gradients has been investigated experimentally [30].The growth rate of these instabilities are small however, unless the temperature gradients are steep.

APPENDIX B: ION CYCLOTRON WAVE RADIATION FROM A MOVING SOURCE
Due to the dispersive nature of the ion cyclotron waves (see Fig. 2), the radiation pattern in Fig. 15 differs from the standard Mach cone found for sound radiation from supersonic sources.The aim of this Appendix is to illustrate the consequences of the dispersive nature of the ion cyclotron waves being radiated by a moving source.The spatial variation of the sound speed will be ignored and the propagating supersonic shock will be modeled by a moving point source.
The qualitative features of the radiation pattern in Fig. 15 can be explained by the following arguments.The dispersion relation (3) with constant electron temperature corresponds to the linear differential equation for two spatial dimensions in an (x,z) plane  axis.The radiation pattern is stationary in the frame moving with U , assuming that transient effects have vanished.By the substitution ∂/∂t → −U∂/∂z we find the form in the frame of reference moving with the source.Division by U 2 and two integrations with respect to z gives Two integration constants are set to zero since we require that no signal has arrived in front of the moving disturbance at z → ∞ and set ñ = 0 there.With U > C s it is here an advantage to introduce the normalized variables ξ ≡ x( ci /C s ) 1 − (C s /U ) 2 and ζ ≡ z( ci /U ) to give which can be solved for given boundary conditions at the position of the source at ζ = 0.The solution corresponds to a source in the electron temperature enhancement for −W < x < W as in, e.g., Fig. 12, moving at a speed exceeding the sound speed in the region outside the electron temperature striation |x| > W. Simple physical arguments for qualitatively estimating the radiation pattern from sources moving at constant speed are given elsewhere [17,31].Here we present results from numerical solutions in Fig. 20.The figure gives the full (x,z) variation of the stationary wave pattern.When sampled at a fixed z position, this wave field is seen moving with constant velocity U and it will give a signal varying with (z,t) very similar to the observation found in Fig. 15.The simple model (B4) does not include any damping mechanism to represent ion Landau damping or the effects of spatially varying electron temperatures.Complete agreement with observations therefore cannot be expected.Nevertheless, we find the basic characteristics of the waves for |x| > W trailing the moving source in Fig. 15 to be well represented.

FIG. 1 .
FIG. 1. Illustration of the low-frequency branch for electrostatic ion waves in magnetized plasmas.The frequencies are normalized by ci and the wave numbers by the effective ion Larmor radius a i ≡ C s / ci .The arrows in (b) indicate the relative magnitude and direction of the group-velocity vectors in the (k ⊥ ,k z ) plane.

FIG. 2 .
FIG. 2. High-frequency branch with ω 2 > 2 ci for the electrostatic dispersion relation for ion waves in magnetized plasmas.Frequencies are normalized by ci and wave numbers with a i ≡ C s / ci as in Fig. 1.The arrows in (b) indicate the relative magnitude and direction of the group-velocity vectors in the (k ⊥ ,k z ) plane.

but with = 1 4 .
FIG.6.Examples of higher-order modes (m = 2, 3, and 4) obtained by numerical solutions of (6), for the same temperature profile as in Fig.3but with = 1 4 .The corresponding γ m values are given in the figure.

5 FIG. 11 .
FIG. 11.Spatial variation of the wave amplitude for the case with ci = pi and applied frequency ω = 0.1π pi as in Fig.10.The perturbation amplitudes are increased in units of 0.1δn/n 0 .

FIG. 13 .
FIG.13.Spatial variation of the potential ψ and ion density n i during the formation and propagation of an electrostatic shock.The results are obtained by a numerical simulation.For the present case we have T e /T i = 25 and ci / pi = 1/2.Here δn i is the actual detected density perturbation at the time where the shock is fully developed.The background density n 0 is normalized to unity.(b) shows a sample of shock fitting for the ion density using model (iii) in TableI.The result is δn i /n 0 = 0.24 ± 0.002, z s = 824 ± 1λ Di , and T s = 140 ± 4.4λ Di .
FIG. 14.(a) Shock velocity U in units of C s0 and (b) the saturated shock thickness T s in units of C s0 / ci , shown for different combination of the parameters T e /T i and ci / pi , as summarized in Table II, and varying δn i /n 0 .One single point at δn i /n 0 = 0.22 having the highest shock velocity and the most narrow width corresponds to Fig. 15.

FIG. 15 .
FIG. 15.Illustration of wave radiation from the propagating shock.The figure is obtained for the (x,τ ) variation of the electrostatic potential at a fixed position z = 900λ Di .For this illustrative case we have a width W = 85λ Di of the electron temperature enhancement.We show the x variation of the electrostatic potential for a fixed z position and varying time (x,z = 900λ Di ,τ ).
FIG.19.Variation of γ with varying for W = 2.5 and = 1/4.See also Fig.4.Note the restricted scale on the vertical axis.

ξFIG. 20 .
FIG. 20.Radiation pattern obtained by solving (B4) for a localized boundary condition at z = 0 and a periodic boundary condition in the x direction.The figure uses normalized axes as defined in the main text.

TABLE I .
Four expressions for n(x) have been used for fitting models to the shocks, with four parameters n mn , n mx , x s , and s l , representing minimum and maximum density, shock position, and width, respectively.

TABLE II .
Parameters used for studies of shock velocity and thickness in Fig.14, with case numbers given there.