Wave beams, packets and pulses in inhomogeneous non-Hermitian media with dispersive gain or damping

Wave beams, packets or pulses are known to be subject to a drift if the properties of the medium change across their extension. This effect is often analyzed considering the dispersive properties of the oscillation, related to the real part of the dispersion relation. The evolution of Gaussian beams/packets/pulses in nonuniform media in the presence of gain or damping is investigated in detail, with particular emphasis on the role of dispersion on both the real and the imaginary part of the dispersion relation. In the paraxial limit, the influence of a non-Hermitian medium on the evolution of the wave can be treated employing the equations derived by Graefe and Schubert in the frame of non-Hermitian quantum mechanics (Phys. Rev. A 83 060101(R)). Analytic solutions of the corresponding paraxial equations are obtained here for a one-dimensional complex dispersion relation characterized by a linear or quadratic dependence on the transverse coordinate (a space coordinate for beams and packets, the time in the co-moving frame for a pulse). In the presence of a constant gradient in both the real and the imaginary part of the dispersion relation, the contribution of the latter can lead to a faster or slower propagation with respect to the Hermitian case, depending on the parameters. In focusing media, a constant gain can counteract dispersive or inhomogeneous damping producing packets of asymptotically constant shape. The analytic formulas derived in this paper offer a way to predict or control the properties of beams/packets/pulses depending on their initial conditions and on the characteristics of the medium.


Introduction
The propagation of wave beams, packets and pulses in dispersive media is a fundamental problem in wave theory, of relevance both in classical and in quantum physics [1,2].Although these three terms denote in principle three different phenomena, they all refer to waves exhibiting some sort of localization, which can be described in terms of an 'evolution' coordinate for their maximum, or center-of-mass, and one or more 'transverse' coordinates related to their extent.In this context, with the word packet we refer to a wave characterized by an envelope localized in some region of space and evolving in time.With pulse we describe the somehow symmetric case of a structure evolving in space (e.g.along an optical fiber) and characterized by a finite extent with respect to the time measured in a frame moving with the group velocity v g (so that time and space coordinates are swapped as compared to a packet).Finally, a beam is the result of a continuous excitation of the wave by an antenna and both coordinates are space coordinates.The commonality in the mathematical description lies in the fact that a first-order derivative with respect to the evolution coordinate is retained in the basic equation, while a second-order derivative with respect to the transverse coordinate accounts for the 'wave effects' .One often refers to the underlying equation as the parabolic equation in the case of a pulse and as paraxial equation in the case of spatial localization, especially for beams [3] (in quantum mechanics, this is of course the Schrödinger equation).In this paper, for the sake of simplicity we will generally use mainly the word 'packet' and in the equations choose the time t as the evolution coordinate, while for the transverse coordinate we choose a spatial coordinate x.Correspondingly, dispersion effects will appear through a wavevector k.We will refer to the approach consisting in retaining second-order derivatives in the 'transverse' direction only as the paraxial approximation.The case of a pulse is readily obtained by substituting the time with an evolution coordinate z, the transverse coordinate x with the delayed time and the wavevector k with the frequency ω.The analogy between the propagation of a diffracting beam and the evolution of a pulse/packet in a dispersive medium has been known for a long time [4], and has been employed recently to investigate with optical means a variety of quantum-mechanical phenomena [5].Although these considerations are not new, in this paper we make an effort to facilitate the 'translation' of our results across these different fields, and we briefly discuss this analogy and the corresponding notation further in section 2.
This paper addresses the evolution of localized wave structures as described above in a medium with following properties: (i) inhomogeneity, i.e. the properties of the medium change in general both across the extent of the packet and in the direction of propagation3 ; (ii) finite loss or gain, i.e. the energy of the packet is not conserved during its evolution; (iii) dispersiveness, i.e. the medium responds differently to different frequencies and/or wavevectors.Problems characterized by the properties listed above range from non-Hermitian quantum mechanics [7] to classical-physics applications like non-Hermitian optics and photonics [8], seismic waves [9], waves in solids [10], plasma waves [11], etc.The present study is originally motivated by plasma-physics applications.In controlled nuclear fusion experiments, the injection of collimated electromagnetic beams as typically employed for electron cyclotron heating [12] or electromagnetic wave packets arising in radial domains prone to small-scale instabilities [13] can be described by paraxial methods.How powerful paraxial techniques are in heating applications has been demonstrated in a number of papers, e.g.[14][15][16][17][18][19], and in [20] for the specific case of geodesic acoustic oscillations [21,22] under the assumption of a Hermitian dielectric response of the plasma.Another important field of application of paraxial waves in media with loss and gain is the study of analytic solutions to the problem of the propagation of a beam in a complex optical 'potential' , e.g.employing ABCD matrices [3,23,24], see e.g.[25][26][27] and references therein.It is stressed, however, that in the present paper no particular assumption on the specific physical system under consideration is made.It is simply supposed that the response of the medium is formulated in terms of a corresponding (generally complex) dispersion relation.The aim of this paper is to present both the mathematical model and the particular solutions derived later in general terms.In this sense, different physical systems 'just' differ in the concrete physical meaning of the coefficients of the dispersion relation.
The impact of the three properties listed above (inhomogeneity, loss/gain, dispersiveness) are not always treated on the same footing.In classical applications of the WKB (semiclassical) limit like geometrical optics [28,29], the phase (or eikonal) is supposed to be real and describes the slow evolution of the wavevector as the wave propagates.The wave 'trajectory' (ray) is determined by the solution of Hamiltonian equations.In this approach, the corresponding Hamiltonian (dispersion relation) is generally assumed to be real and possible damping or gain have no influence on this 'trajectory' .Besides the gradients of the dispersive properties, however, a spatial variation of the damping (or gain) also contributes to the displacement of the wave.This can be expected, since the wave has a larger amplitude in regions of weaker damping (or stronger gain), so that its center-of-mass is shifted towards these regions.If the dispersive properties and the gain/damping can be properly tuned, the combination of both can be used to attain a better control of the wave and thus new functionalities, like in PT -symmetric optics and, more generally, non-Hermitian optics [30][31][32][33].
The treatment we apply in this paper to quantify the effect of gain/damping on the propagation of paraxial packets was developed by Graefe and Schubert [34] in the context of non-Hermitian quantum mechanics.These authors employ a paraxial expansion of the underlying wave-kinetic equation [35], derived from a non-Hermitian Hamiltonian operator.The solution is given in terms of the parameters characterizing a Gaussian Wigner function, i.e. the position and the momentum (wavevector) at the center of the packet, a matrix characterizing its Gaussian shape, and its amplitude.In this approach, the trajectory of the packet center is the sum of two contributions.One stems from the real part of the Hamiltonian and appears in the form of Hamilton's equations of motion, like in standard geometric optics.The imaginary part of the dispersion relation, which is associated to the damping (or growth) of the packet, generates an additional contribution (a gradient flow with time-dependent metric), in which the Gaussian profile of the packet determines the motion of its center together with the properties of the medium.An application of these paraxial equations to the case of electromagnetic beams propagating in a medium with a complex, PT -symmetric potential was given in [36].For quadratic Hamiltonians, this approach was demonstrated to be equivalent to the solution of the paraxial equations for the corresponding Gaussian field [37].In this case the equations are solved for a complex center of the packet, and an appropriate procedure to project the complex phase-space position to physical (real) phase-space is given.
A particular focus of the present paper is on dispersive damping (wavevector or frequency dependent as discussed above and in section 2).Dispersive absorption of electromagnetic fields is a common situation [11] in space [38], astrophysical [39] and laboratory [40] plasma, but it can be considered as a general behavior of dispersive systems (sound waves in solids, etc).In hot plasmas, in particular, damping, but also resonant destabilization (gain), are often determined by collisionless resonant (inverse) damping processes like Landau damping, cyclotron damping.The resonance condition is determined by the position in both real space and wavevector space.The particular mechanism responsible for loss (collisional, collisionless, etc) or gain (population inversion, anisotropic velocity distribution, equilibrium gradients, etc), is not essential in the description adopted here, as long as it can be encoded in the (imaginary part of the) dispersion relation.By the same token, the details of the particle distribution in phase space, which can be essential for the determination of the complex dispersion relation for a particular wave, do not need to be further specified here.Some examples are discussed in more detail in section 6.The effect of dispersive absorption on Gaussian beams [41] and packets [42] in nuclear fusion applications has been discussed recently.Depending on the parameters, the strength of the damping can vary significantly, and the non-dispersive part of the damping can be smaller or larger than the dispersive part.In [42], analytic solutions for the case of a uniform medium were derived and the results applied to the description of geodesic acoustic packets.It was shown among other things that a damping rate exhibiting a quadratic dependence on the wavevector leads to a diffusive broadening of the packet, which adds to the dispersive broadening.This behavior, which is expected from the theory [3], was confirmed by physically comprehensive kinetic simulations.Moreover, dispersive damping leads to an additional contribution to the (Gouy) phase.A question not addressed in [42] is the role of damping on the transverse displacement of the packet when the medium is inhomogeneous.This is discussed in detail in the present paper, in particular on the basis of analytic solutions obtained for the case of a dispersion relation depending linearly or quadratically on the space coordinate and quadratically on the wavevector.The theory discussed here, however, is not limited to quadratic Hamiltonians.
A quadratic Hamiltonian represents a simple and physically relevant model, for which analytic solutions of the Schrödinger equation are possible, also in the case of a non-Hermitian Hamiltonian, see e.g.[43], and can be found in terms of Gaussian packets.In section 3, the paraxial Equations of Graefe and Schubert are derived for a general pseudo-differential Hamiltonian (under the assumption that the paraxial limit is satisfied) and for a real packet center, as in [34].In this paper, however, we start from the (Schrödinger) equation for the field, rather than from the corresponding wave kinetic equation.This derivation, that produces the same equations for the packet center as in [34], has the advantage that it also leads to a direct derivation of an equation for the phase.Section 4 reports the derivation of analytic solutions of the paraxial equations for the evolution of a one-dimensional packet (or for a two-dimensional wave beam).Section 4.1 summarizes previous work for reference, while a linear and quadratic dependence of the dispersion relation on the transverse coordinate are investigated in sections 4.2 and 4.3), respectively.Examples of the solutions derived in section 4 are discussed in section 5.As expected, the impact of the imaginary part of the dispersion relation on the propagation depends also on the packet size, so that the drift of the packet can be either accelerated or decelerated with respect to the Hermitian case.Dispersive damping is shown to act in a similar way as a quadratic dependence of the damping coefficient on the transverse coordinate.Summary and concluding remarks follow in section 6.

Dispersion relations for beams and packets and notation
As stated in the introduction, the analogy between a diffracting beam and a packet or pulse evolving in the presence of dispersion is well known.Nevertheless, it is useful to briefly compare two examples of dispersion relation in the paraxial limit, one for a wave beam and one for a packet/pulse, to show explicitly how they can be described by a formally identical Hamiltonian function, which just expresses the dispersion relation in the mathematical description presented later.For a wave beam, the frequency of the injected wave is known and the dispersion relation is usually expressed as k = K(ω) with k = |k| or, multiplying both sides by c/ω, as n = N(ω).For a wave packet, we assume that a localized displacement in the medium from its equilibrium at some initial time leads to an oscillation at a frequency ω = Ω(k), determined by the local dispersion relation in terms of the wavelength(s) characterizing the displacement.For a pulse, the space and time coordinates are swapped with respect to a packet, and so are wavevector and frequency.The letter β is often used [3,44] to cast the dispersion relation in the form b = β(ω), where b is the propagation constant.
We assume that an electromagnetic wave beam satisfies a Helmholtz equation in the (r, z)-plane, with a possibly complex refractive index N 2 (ω; r, z) assumed to be weakly (in the WKB sense) varying in space.In the paraxial limit, we can assume propagation e.g.mainly in z-direction and expand the dispersion relation r .To choose a concrete example, we take a refractive index with a parabolic dependence on the transverse coordinate r, N 2 (r, z) = N 2 0 (z) + q(z)r 2 , but the present considerations are obviously not limited to this case.In sections 4 and 5, we will consider both a linear and a quadratic scaling of the dispersion relation with r.A Taylor expansion of the square root yields [45] where the dependence of N on ω is left implicit to simplify the notation.In the paraxial WKB approach adopted in this paper, the dynamics of the wave is described by a set of equations involving a Hamiltonian function (and derivatives thereof).For the present case we can write e.g. with Analogous considerations can be applied to a dispersion relation of the form which is typical of a number of physical systems, like an array of coupled pendula, de Broglie matter waves or several kind of waves and oscillations in plasma physics, like Langmuir waves, kinetic Alfvén waves, dispersive geodesic acoustic oscillations, etc.Also in this case we can expand the previous formula assuming weak dispersion and a non-dispersive part of the frequency of the form F(r, t) = F 2 0 (t) + q(t)r 2 , with the result The corresponding Hamiltonian can be written as with Comparing equation ( 2) with equation ( 5), we see that the role of the 'axial' coordinate z in the former case is played by the time t in the latter.Similarly, the role of the 'axial' wavevector n z is taken by the frequency ω.For a pulse, the Hamiltonian is found by swapping the time and space (and the corresponding Fourier-transformed) coordinates in equation (5).In this case, the Taylor expansion of the dispersion relation is performed around the frequency midband and the coefficient G of the dispersive term is often denoted by β ′ ′ or β 2 and called group velocity dispersion, cf table 1.
Apart from this redefinition of the 'axial' , or 'propagation' coordinate (and its dual coordinate in Fourier space), the Hamiltonians in equation ( 2) and in equation ( 5) are formally the same.Thus, the dynamics shown in figure 1 can be attributed to (i) the time evolution of a packet initialized at t = 0 and evolving in time, as indicated by the label on the x-axis, as well as (ii) to the propagation of a beam injected from the left and propagating in the (z, r)-plane, where the z-coordinate takes the place of the time t, or (iii) to a pulse in the (z, η)-plane, with the delayed time η = t − z/v g representing the 'transverse' coordinate.It should be noted that figure 1 represents the limit of vanishing group velocity, see discussion in section 5.1.
We stress that the coefficients F and G in equations ( 2) and ( 5) are in general complex in the present treatment.As shown by the examples above, the quantities F and G containing the information about the properties of the medium can depend (weakly) on the evolution coordinate, i.e. on time t, resp.on the propagation coordinate z.Although the derivation presented in section 3 includes this case, in the examples considered later in this paper this dependence is not investigated.The only non-uniformity retained in sections 4 and 5 is along the 'transverse' coordinate r.
At the end of this section, a last remark on the notation and nomenclature adopted for a Gaussian envelope should be made.We will be dealing with wave structures with a transverse Gaussian profile, which we will write in the form (taking again a packet as an example) where q(t) denotes the position of the center of the packet, the complex coefficient s(t) = s(t) + iϕ(t) describes the Gaussian envelope and we assume ϕ(t) > 0. It is evident that ϕ = 2/w 2 (t) is related to the transverse width w of the field at 1/e level.The parameter s describing the relative phase of different parts of the packet can be both positive and negative, s = 0 corresponding to an in-phase oscillation of the whole packet.In beam optics, this parameter corresponds to the curvature of the phase front, while in fiber optics, where the pulse is localized in time (in the co-moving frame), the parameter s is connected to the frequency (rather than wavevector) spectrum and expresses the chirping of the pulse4 .

Derivation of the paraxial equations for a complex Hamiltonian
We are interested in a paraxial solution of the Schrödinger equation Later in this paper, we derive and discuss analytic results obtained for the Hamiltonian operator which leads to a dispersion relation of the form The derivation reported in this section, however, is valid for a more general, pseudo-differential Hamiltonian operator of the form (the time dependence of the field ψ is omitted for simplicity).Appendix shows how the standard differential operator in equation ( 8) can be derived as a special case of equation (10).We derive the paraxial equations for a field ψ of the form Here, the eikonal S(x) is a real function and the Gaussian profile is singled out later from the amplitude a(x).
Introducing the new variable x − x ′ = s, equation (10) becomes The first derivative of S is related to the packet wavevector, while the second derivative expresses the 'slow' (in terms of the WKB limit) variation of the wavevector as the packet propagates, and we can expand In the paraxial WKB ordering, the medium varies on a scale L much larger than the transverse dimensions of the packet w.The amplitude a is also expanded to the second order in s, while in the dispersion relation H only the first-order term needs to be retained.With the notation H x = ∂H/∂x and omitting the arguments of the various terms for simplicity, we can hence rewrite equation (12) as where the term s 2 H x a ′ /2 is of higher order and can be omitted.Now we exploit the fact that Integrating by parts, the derivatives with respect to k are applied to H, so that where the δ-function arising from the integration over s implies that in the second step H and its derivatives must be evaluated for k = S ′ .Since the displacement s does not appear in the equations any longer, in order to keep consistency with the notation employed in previous work, the letter s is used from now on to denote the real part of the complex Gaussian parameter s(t) = s(t) + iϕ(t).The real quantities s and ϕ have been introduced in section 2 and are defined explicitly in equation ( 16) below.Now we select the field ψ(x, t) to take the shape of a Gaussian packet centered around the real position q(t): where p is real and expresses the wavevector on the trajectory of the packet center, while A is a complex amplitude.Calculating the derivatives of a and S and substituting them into equation ( 15) leads to where the derivatives of H have the same arguments as the first term on the right-hand side.The final step is to expand H and its derivatives around the packet center (q, p).Expansion around p yields while expansion of x = q + (x − q) up to second order in x − q gives finally where again in the last two equations the derivatives of H have the same arguments as the first term.The left-hand side of the Schrödinger equation ( 7) with the packet ansatz ( 11) and ( 16) becomes Equating equal powers of x − q in equations ( 20) and ( 19) yields If the damping is ordered small (γ/ω ∼ λ/L), as in standard paraxial theory [15], the imaginary part of H drops out of all the derivatives and survives only in the last equation.In this limit, the propagation of the packet is determined by the real part of the Hamiltonian only.In the general case, however, H k and H x in equation ( 21) are complex and the evolution equations for the center of the packet can be found splitting the real and imaginary part.With H = H − iΓ (here we follow the convention [34,37,42] that a positive value of γ corresponds to damping), the result is This set of equations for the packet center has been obtained originally in [34] from a paraxial expansion of the wave kinetic equation, see also equations ( 45) and( 46) below.Note that through the imaginary part of H the propagation acquires a dependence on the shape of the packet described by s(t) and ϕ(t).
The equation (22) for s is the standard Riccati equation of paraxial theory [15,47], but the derivatives of the Hamiltonian are in general complex-valued.Finally, the equations for amplitude and phase can be found by separating real and imaginary part of equation ( 23): We note that p and q in the first term on the right-hand side of the equation for the phase are real and can be calculated either from equation (24) or applying the projection (42) to the solution of the equations for the complex trajectory (41) discussed in section 4.2, see also [48].The corresponding term on the right-hand side of equation (B9) of [42] is to be interpreted in the same way (in that reference, the evolution of the phase was actually discussed only under conditions for which q = 0).Equation ( 7) can be solved in terms of the parameters determined through equations ( 22), ( 24) and (25) as long as the paraxial approximation is satisfied.In terms of characteristic lengths of the problem, this requires that the packet width w is intermediate between the wavelength λ and the inhomogeneity scale of the medium L, i.e. λ ≪ w ≪ L [14].It is noted in passing that these are precisely the conditions under which a paraxial (Schrödinger) equation can be derived from a standard wave equation involving second-order derivatives, see e.g.[45] for a recent discussion.While a Gaussian is an exact solution of a paraxial equation in the case of a quadratic Hamiltonian, in general the solution given in terms of the equations derived above is only an approximation to the exact result.On the other hand, the simplicity of the ordinary differential equation equations ( 22), ( 24) and ( 25) makes them very attractive for practical applications (the method has proven to be quite robust also in critical cases [49]).
A standard method in optics, which is related to the approach employed here, is based on the so-called complex ABCD matrices [3,24,25].The method is designed to model a cascade of optical elements, e.g.multielement optical resonators.Smoothly nonuniform media are described as a set of successive complex elements.The method adopted here treats this situation naturally and seems to be more powerful in the treatment of complicated trajectories and arbitrary dispersion relations.A detailed comparison of these approaches is, however, outside the scope of this work.

Analytic solutions for a quadratic Hamiltonian
In this section, we derive analytic solutions for the propagation of a wave packet according to equations ( 22), ( 24) and ( 25), considering initially a uniform medium (this part summarizes the results of [42]), then including first a linear dependence and then a quadratic dependence of the Hamiltonian on the space variable.For the sake of clarity, we adopt here the time t as a coordinate describing the propagation, as appropriate for the evolution of a wave packet, and we will use r to denote the (transverse) space variable.According to the discussion in section 2, the results derived in the following can be applied to a wave beam or pulse through an appropriate change of the coordinates.

Homogeneous medium
The solution of the paraxial equations including damping for a uniform medium has been derived in [42] (some of these results are known in different contexts, see e.g.[3]).The Hamiltonian takes the form with constant H 0 and H kk .In this case the solution of equation ( 24) for the center of the packet (we assume p 0 = p(t = 0) = 0) are simply p = 0, q = q 0 = q(t = 0).The Gaussian parameters s and ϕ can be determined solving equation (22).Since H rr = 0 and H kr = 0 if the medium is uniform, the solution of equation ( 22) involves only H kk and Γ kk and reads with Here we can identify a typical time scale characterizing the dispersion effects (corresponding to the Rayleigh range for beams and the dispersion length for pulses), while the wavevector-dependent damping introduced by Γ kk leads to a second time scale which characterizes a diffusive broadening (for positive Γ kk ) of the packet due to the selective absorption of the high spectral components [3].We will see below that t D characterizes also a slowing-down rate of the packet propagation in the presence of gradients.We observe incidentally that equation ( 28) suggests two further (inverse) time scales H kk s 0 and Γ kk s 0 .The packets considered in this paper, however, have been calculated with s 0 = 0 and we do not discuss the time scales related to s 0 any further.With s 0 = 0, according to equation (27), the width of the packet can be written The same solution can be expressed in terms of the matrix G expressing the Gaussian shape of the Wigner function W(z, t) related to the field ψ, where z = (r, k) is a phase-space variable, z 0 = (q, p), and the symmetric matrix G is related to the Gaussian parameters s and ϕ used above by and for H rr = Γ rr = 0, H kr = Γ kr = 0 has the solution which can be shown to coincide with equation (27).
For the amplitude and the phase of the packet, only the term H kk s on the right-hand side of equation ( 25) contributes and the solution becomes for the amplitude, and which is a generalization of the Gouy phase in the presence of dispersive damping.

Inhomogeneous medium: linear profiles
We consider now the complex model Hamiltonian H = H − iΓ with (complex) constants H 0 , H r and H kk .With this Hamiltonian, the solution of equation ( 22) for the Gaussian envelope is the same as in the homogeneous case reported in the previous Section, since again H rr = 0 and H kr = 0.The presence of a complex gradient H r , however, implies that the position of the center of the packet changes in time, unlike in the uniform case.
A straightforward solution of the equations of motion for the packet center can be computed employing the paraxial equations for the complex packet center (q, p) derived in [37] (see also [42]), i.e.
which are the standard Hamilton equations with a complex Hamiltonian.The physical (real) values of position and wavevector at the center of the packet can be obtained through the projection [37] where the elements of the G-matrix are related to s = s + iϕ through equation (33).Equation (41) with the model Hamiltonian ( 40) can be solved straightforwardly with the result Applying the projection operator defined in equation ( 42), the physical phase-space coordinates of the center of the packet are found to be (for simplicity, we choose as initial conditions q 0 = 0 and p 0 = 0) The same result can be achieved from the integration of equation (24), which in terms of the matrix G read As stated above, for the model Hamiltonian (40) the elements of the matrix G are the same as in the homogeneous case and are expressed by equations ( 35)- (37).Taking into account that Γ k = Γ kk k, equation ( 45) is Multiplying both sides by 1 + Γ kk G rr,0 t and integrating the left-hand side by parts one finds immediately the expression for p(t) reported in equation (44).Integration of equation ( 46) with the above solution for p is lengthy but straightforward, involving only integrals of rational functions.The final result is which can be verified to coincide with the solution for q in equation ( 44), once equations ( 36) and ( 37) are substituted.
Concerning the trajectory given by equation ( 48), one can immediately note that the gradient flow due to the non-uniform imaginary part of the dispersion relation (finite Γ r ) forces the trajectory into the region of smaller damping (or higher gain), as discussed in section 1.In the absence of non-uniform damping, the trajectory follows the standard t 2 -law, with a constant acceleration due the constant force −H r .This can be modified if dispersive damping is included (finite Γ kk ).Thus, even in the absence of non-uniform damping, a constant dispersive damping (Γ kk > 0) leads to an increase of the denominator of equation ( 48) and hence to a slowing-down of the motion of the packet center.A dispersive gain Γ kk < 0 leads to the opposite effect.
Different time scales can be identified in equation (48).In the denominator, the term Γ kk G rr,0 coincides (assuming s 0 = 0 and hence G rr,0 = ϕ 0 ) with the inverse time scale 1/t D defined in equation (30).A short t D will lead to a stronger impact of the denominator on the evolution of the trajectory.Normalizing the position q to some typical scale of the system L, and assuming again s 0 = 0, so that also G kr,0 = 0, five further time scales can be identified in equation (48) examining the term within curly brackets.The first time scale is connected with the term linear in time and reads t 1 ∼ L/Γ r G kk,0 ∼ L/Γ r w 2 0 .This time scale becomes shorter with increasing initial packet width.The two time scales related to the quadratic term in t are t 2H = (L/H r H kk ) 1/2 and t 2Γ ∼ (L/Γ r Γ kk ) 1/2 .The former is akin to the acceleration of a particle in a non-dissipative Hamiltonian system mentioned above, while the later expresses a similar effect of the combined inhomogeneous and dispersive damping.Finally, there are the two scales , which become short for small initial packet width.Obviously, the cubic term in t always dominates for t → ∞, but such long times might or might not be relevant in practical applications depending on the parameter of both the medium and the packet.

Inhomogeneous medium: parabolic profiles
We start deriving the solution for the Hamiltonian (model without dispersive damping, Γ kk = 0).This Hamiltonian is PT -symmetric if Γ 0 = 0 and the related packet dynamics has been discussed in [36].Here we obtain the same result as in [36] using the complex paraxial equations, i.e. equation ( 41) for the complex trajectory and equation ( 22) for the envelope.
Assuming H rr H kk > 0 and defining the real frequency Ω = √ H rr H kk , the solution of equation ( 41) is H rr (50) with α and β constants related to the initial conditions to be determined below.For the Hamiltonian under consideration, equation ( 49), the evolution of s is the same as for a lossless lens-like medium [50], cf also appendix of [51]: (here we retain s 0 ̸ = 0 to facilitate the comparison with [36]).Splitting the real and imaginary part yields G kr = −s/ϕ and G kk = 1/ϕ, which can be inserted in equation ( 42) to find the real-space projection of the complex trajectory (50).Choosing real α and β in equation ( 50) yields and hence which coincides (apart from the different sign of Γ r due to the convention followed here) with the solution derived in [36] using the equations for the real packet center.If dispersive damping is included (Γ kk ̸ = 0), the Hamiltonian becomes i.e. the coefficient of the last term is now complex, and the solutions obtained above are replaced by H rr (55) with a complex frequency Ω = H rr H kk and This result can be easily checked to reproduce equation (51) for Γ kk = 0. Splitting the real and imaginary part of equation ( 56), and consequently projecting the complex trajectory (55) to real space, becomes rather cumbersome in this case.However, it is easy to check that q → 0 for long times.Moreover, the Gaussian envelope attains asymptotically a constant value independent of the initial conditions, s → i H rr /H kk .This result remains true if also the damping (or gain) coefficient is supposed to have a parabolic dependence on the spatial coordinate (Γ r = 0, Γ rr > 0).In this case Ω = H rr H kk and s → i H rr /H kk .These findings generalize those obtained by Kogelnik [25] for non-dispersive damping (Γ rr ̸ = 0, Γ kk = 0).It can be noticed that the damping term −Re[H kk s]/2 in equation ( 25), which is constant for long times, can be compensated in principle by an equal constant gain Γ 0 < 0 (with the convention adopted here for the sign of the imaginary part of the Hamiltonian), leading to a packet with asymptotically constant width and amplitude.

Discussion of the packet dynamics
The analytic solutions derived in section 4 are discussed in this Section on the basis of some examples, obtained from numerical integration of equations ( 41), ( 22) and ( 25), using equation (42) to project the complex trajectory, solution of the complex Hamiltonian equations equation (41), to real space.It is remarked that the equations derived in section 3 for the parameters of a Gaussian packet are a set of ordinary differential equations, which can be solved straightforwardly using standard numerical algorithms.In this section, the equations are solved using a fourth-order Runge-Kutta method [52][53][54].These numerical results have been checked to reproduce the analytic solutions derived in section 4. As an illustrative example, we compare for one case the solution of these Equations for the packet parameters with a direct solution of the paraxial equation ( 7) produced by a pseudo-spectral code [54,55] at a much more significant computational cost.
We start again from the case of a homogeneous medium, as a reference for the following cases of a linear and a parabolic medium.In all plots, H 0 = Re[H 0 ] = 1, so that the packet oscillates at the angular frequency ω = 2π.In general, we select p 0 = 0 as initial condition (no 'oblique injection').A deviation from a straight 'trajectory' , however, develops naturally in the presence of medium non-uniformity, as shown in the examples below.For an analytic treatment, which includes a comparison with exact solutions for oblique beam launch and high frequency waves in a linear medium, the reader is referred to [41,49].

Homogeneous medium
The propagation of a Gaussian beam/packet in a homogeneous, lossless medium is well known [2] and is characterized by a diffractive/dispersive broadening occurring on the Rayleigh scale.This is visible in the left panel of figure 1, which compares a case without dispersive damping (Γ kk = 0) with a case with finite dispersive damping (right panel).In both cases, the dispersion is anomalous (H kk = 1 × 10 −4 > 0) and the initial Gaussian width (1/e-level of the field) is w 0 = 0.1, so that t R = 50.
Concave phase fronts can be observed in the absence of damping, while dispersive damping (right panel) not only reduces the field amplitude, but also leads to flatter phase fronts, i.e. the whole packet tends to oscillate at the same frequency across its whole extent.For the parameters considered in the right panel (H kk = Γ kk ), the 'diffusive' time scale introduced in equation ( 30) equals the Rayleigh time.It is recalled [42] that when H kk ̸ = 0 and Γ kk → 0 (i.e.t D → ∞), the width of the packet follows the well-known Gaussian-optics formula, i.e. it scales linearly with time for t ≫ t R .This is easily seen from the corresponding limit of equation (31).In this case, the amplitude scales as 1/ w(t) and the energy of the packet is conserved.In the opposite limit Γ kk ̸ = 0, H kk → 0 (i.e.t R → ∞), the packet evolves like in a diffusion problem, so that w ∼ t 1/2 for t ≫ t D , the amplitude is proportional to 1/w, and energy is not conserved.

Inhomogeneous medium: linear profiles
Figure 2 shows the evolution of a packet, excluding first any damping effects, but allowing for a linear spatial variation of its frequency H 0 through H r > 0 in equation (40).A case with normal dispersion (H kk < 0, left) and one with anomalous dispersion (right, for the same value of H kk > 0 as in figure 1) are compared.The center of the packet follows a parabolic trajectory, while the wave fronts become tilted with respect to the central trajectory, as particularly evident in the case of anomalous dispersion.It can be noticed that in the left plot of figure 2 the group velocity is positive, while the phase velocity is negative, whereas both velocities are negative in the plot on the right.
Introducing damping changes this picture significantly.In figure 3, the evolution of the packet is computed for the same dispersion parameters as in the right panel of figure 2, but now with w 0 = 0.3, Γ kk = 1 × 10 −4 and Γ r = −0.2/3.Since we take Γ 0 = 0, Γ r < 0, the medium leads to a damping in the lower  half-plane (r < 0) and a gain in the upper half plain.In figure 3, besides the solution of the paraxial equations, the Schrödinger equation (7) with the Hamiltonian of equation ( 40) is also solved numerically.The two solutions show an excellent agreement, as expected, given the fact that a Gaussian is an exact solution of the Schrödinger equation for a quadratic Hamiltonian.It can be seen that the packet moves upwards, i.e. into the region of higher gain and its phase fronts become nearly flat.
The different contributions to the trajectory of the packet center for the same complex Hamiltonian as in figure 3 are detailed in figure 4 for two different initial widths, namely w 0 = 0.3 (left) as in figure 3 and w 0 = 0.08 (right).In the legend, 'term 1' refers to the first term within brackets in equation (48), which scales as t/(1 + t/t D ), while 'term2, H' shows the contribution of the term proportional to H r H kk t 2 /(1 + t/t D ), 'term2, Γ' that proportional to Γ r Γ kk t 2 /(1 + t/t D ), and 'term 3' denotes the term proportional to t 3 /(1 + t/t D ).The parameters H r , Γ r , H kk and Γ kk are selected in such a way that the terms proportional to t 2 in the brackets on the right-hand side of equation ( 48) cancel each other.As expected from the discussion in section 4.2, term 1 is seen to dominate when the initial width of the packet is larger (w 0 = 0.3, left panel), while term 3 is dominant for smaller initial widths (w 0 = 0.08, right panel).The diffusive time scale t D decreases from t D = 450 for w 0 = 0.3 to t D = 32 for w 0 = 0.08.In the former case the denominator on the right-hand side of equation ( 48) leads only to a small deviation with respect to a straight trajectory q ∝ t for the time range considered in the plots, where the packet is followed up to t = 100.In contrast, for w 0 = 0.08, i.e. t D = 32, the denominator in equation (48) reduces significantly the drift of the trajectory, which follows a t 2 -law at long times rather than a t 3 -scaling.

Inhomogeneous medium: parabolic profiles
As it is well known, for a Hermitian medium with a parabolic potential well (lens-like medium in optics), the width of the beam/packet oscillates as a consequence of the competition between diffractive/dispersive broadening and focusing properties of the medium.A constant width results from taking s 0 = 0 and the  width of the fundamental eigenmode as initial width, w 2 0 = 2(H kk /H rr ) 1/2 [17,25,56].If the center of the packet is displaced with respect to the minimum of the potential, it also oscillates around the equilibrium position.The (PT -symmetric) case of a medium with a linear damping/gain profile in space was discussed in detail in [36], where it was shown that even if the packet starts with its center at the potential minimum, it undergoes an oscillation driven by the balance between the drift towards higher-gain (smaller damping) regions and the focusing effect of the medium.
As discussed in section 4.3, including dispersive absorption, or an absorption coefficient with a parabolic profile, leads to a damping of both oscillations mentioned above, i.e. that of the Gaussian envelope (which tends to a constant independent of the initial conditions) and that of the packet center.An example is shown in figure 5, where the evolution of the packet has been computed for H rr = −100, Γ rr = 30, is added in such a way that the amplitude of the packet tends to a constant value as well.Figure 6 compares in more detail the long-time behavior of the width w(t) and the amplitude A(t).According to equation (25), the latter evolves for long times as A ∼ exp(−Re[H kk s]/2) if no extra gain is included.
The last example illustrates a case with stronger dispersive damping (Γ kk = 3 × 10 −4 in figure 7 as compared to Γ kk = 3 × 10 −5 in figure 5), showing that for long time the packet retains also a constant phase-front curvature.This is in agreement with the asymptotic behavior of the Gaussian envelop as predicted by equation ( 56), i.e. s → i H rr /H kk for long times.If H rr and H kk are predominantly real, as in figure 5, the real part s of s, which describes the curvature of the phase front, is small and the wave front is flat.A larger value of Γ kk leads to a larger real part of s at large times and hence to asymptotically curved wave fronts.We note that in figures 5 and 7, both Γ kk and Γ rr are finite, but a qualitatively similar behavior can be reproduced also when either Γ kk or Γ rr vanish.

Summary and concluding remarks
In this paper, the evolution of a wave packet, or the propagation of a wave beam or pulse, in a medium characterized by the presence of dispersive gain or damping, has been addressed in the paraxial limit.The equations for the parameters characterizing a Gaussian packet (beam center, beam envelope, amplitude and phase) have been derived starting from the underlying paraxial / parabolic / Schrödinger equation for a general pseudo-differential Hamiltonian.The result of this derivation is a combination of the equations for the real beam center, as in the corresponding original derivation for the related Wigner function [34], an equation for the complex Gaussian envelope, which is a Riccati equation whose coefficients involve second-order derivatives of the complex Hamiltonian [37], and an equation for the amplitude and the phase.The last quantity is not directly accessible in a treatment based on the Wigner function (but see also [48]).For the case of a quadratic Hamiltonian, the trajectory can be also calculated assuming complex phase-space coordinates of the packet center.These satisfy standard Hamiltonian equations with a complex Hamiltonian.The physical phase-space coordinates of the packet center are the found by projecting the complex solution onto real space using the procedure introduced in [37] and summarized in section 4.2.In the derivation of the equation and their analytic solutions, the 'propagation coordinate' has been chosen to be the time t and the 'transverse coordinate' to be a spatial direction r, having in mind the application of these equations to packets.However, the same results can be applied to a beam or a pulse according to the discussion in section 2 with an appropriate change of coordinates, as summarized in table 1.
Analytic solutions of the paraxial equations have been computed for the case of a Hamiltonian in which both the real and imaginary part depend on the wavevector, i.e. in the presence of both dispersion and dispersive damping (or gain).In the presence of a constant gradient, i.e. for a 'linear' medium, the dispersive damping (in our convention, Γ kk > 0) is shown in general to slow down the drift of the beam/packet due to the gradients (while a dispersive gain Γ kk < 0 accelerates it).This effect becomes stronger with decreasing packet width.On the other hand, a gradient of the (non-wavevector-dependent) damping/gain coefficient leads to additional contributions to this drift.In the limit Γ kk → 0, we find a contribution that scales linearly with time and is large for large packet widths, while another contribution evolves as t 3 and is larger for smaller widths.In the case of a 'parabolic' medium, it is shown that a k 2 -dependent absorption produces a similar effect as a quadratic dependence of the damping rate with respect to space [25], i.e. it leads to a constant packet width for long times.Compensating the loss through a corresponding constant gain leads to a packet with constant shape and amplitude at large times.
As mentioned in section 1, the present work was originally motivated by applications to electromagnetic wave packets that can form radially in tokamak plasmas, with specific reference to geodesic acoustic packets.They exhibit interesting properties, like different dispersion regimes (sign of Re[G] in equation ( 9)) depending on the ratio between electron and ion temperature, see [57] and references therein.The real and the imaginary part of the dispersion relation exhibit different dependencies on the plasma parameters, and hence their magnitude and gradient can largely differ.For example, in the real part of the dispersion relation the role of dispersion is often just that of a small correction to a k-independent oscillation frequency, while the dispersive damping can be much larger than the k-independent damping.Moreover, the geodesic-acoustic-mode (GAM) frequency depends only weakly on the tokamak safety factor q (measure of the winding of the magnetic field lines), while the damping rate strongly depends on it.An application of the treatment developed in this paper to GAM packets is planned for the near future.Also including a source should allow to expand the range of applications of the method.
Another typical application of paraxial methods in plasma physics is the propagation and absorption of high-frequency collimated beams, typically in the electron-cyclotron range of frequencies.Also in this case, the real and the imaginary part of the dispersion relation have in general quite different dependencies on the plasma parameters and the present approach could prove itself useful.In this paper, we considered only a scalar wave equation.In an anisotropic plasmas, different modes of propagation exist, which in theoretical treatments are usually separated by ordering the anti-Hermitian part of the dielectric tensor as small, so that the dispersion tensor is Hermitian to lowest order and the eigenvalues are real.A generalization of the present approach to vector wave equations is required in order for it to be applicable in this case.
In the case of fusion plasmas, the dielectric properties of the medium can be influenced only to a limited amount.In many wave applications, the plasma configuration is treated as given and the evolution of beams or packets is calculated correspondingly.In optics, on the other hand, it is in general rather the opposite, and the properties of the medium can be tuned, at least to some extent, according to the desired propagation characteristics.In this sense, the approach employed here and first derived in [34] offers both a powerful technique to treat the propagation of paraxial waves with rather arbitrary complex dispersion relation and provides an explicit demonstration of how the properties of the waves can be predicted and controlled at the level of analytical formulas.

Figure 6 .
Figure6.Width (left) and amplitude (right) of a packet propagating in a lens-like medium with the same parameter as in figure5(apart for q0 = 0).In the right plot, the amplitude obtained by including a constant gain Γ0 = −|ImΩ|/2, which compensates the damping due to Γ kk and Γrr, is also shown.The value of Γ0 does not alter the evolution of the width in the left plot.The thin red lines display the long-time trends predicted analytically.

Figure 7 .
Figure 7. Absolute value of the field for the same parameters as in figure 5, except for Γ kk = 3 × 10 −4 .A constant gain Γ0 = −|ImΩ|/2 is included.In this case, the wave front is asymptotically curved.

Table 1 .
Analogy between the description of wave packets, pulses and beam.