Density operator approach to turbulent flows in plasma and atmospheric fluids

We formulate a statistical wave-mechanical approach to describe dissipation and instabilities in two-dimensional turbulent flows of magnetized plasmas and atmospheric fluids, such as drift and Rossby waves. This is made possible by the existence of Hilbert space, associated with the electric potential of plasma or stream function of atmospheric fluid. We therefore regard such turbulent flows as macroscopic wave-mechanical phenomena, driven by the non-Hermitian Hamiltonian operator we derive, whose anti-Hermitian component is attributed to an effect of the environment. Introducing a wave-mechanical density operator for the statistical ensembles of waves, we formulate master equations and define observables: such as the enstrophy and energy of both the waves and zonal flow as statistical averages. We establish that our open system can generally follow two types of time evolution, depending on whether the environment hinders or assists the system's stability and integrity. We also consider a phase-space formulation of the theory, including the geometrical-optic limit and beyond, and study the conservation laws of physical observables. It is thus shown that the approach predicts various mechanisms of energy and enstrophy exchange between drift waves and zonal flow, which were hitherto overlooked in models based on wave kinetic equations.

The first attempts to describe turbulence in plasma were made as early as the 60's, in works by Vedenov, Drummond and Pines, Kadomtsev, and others [16]. About a decade later, Hasegawa and Mima proposed a new approach to low-frequency turbulence in nonuniform strongly magnetized plasma [17]. In their approach, the fluid approximation and a perturbation approach were applied, resulting in the nonlinear Hasegawa-Mima equation (HME), which has since become a very popular model in the theory of zonal flows, drift waves and their interactions [18][19][20].
In a different branch of physics, the theory of atmospheric fluids, RW-related turbulence occurs as a result of fluid rotation, one example being the turbulence in Earth's atmosphere caused by the planet's rotation. Its evolution equations turn out to be very similar to the * Electronic address: https://orcid.org/0000-0002-9960-2874 Hasegawa-Mima equation [11]. In such systems, the role of drift waves is played by the Rossby waves, which allows direct analogies between plasmas and atmospheric flows [21][22][23].
It is known that turbulence transfers energy from large scales to smaller ones, thus causing it to dissipate. For instance, it is this kind of dissipation which leads to difficulties with magnetic fusion confinement. Therefore, it is important to study dissipative processes associated with drift or Rossby waves and zonal flows. There are many approaches to describing dissipative (open) systems, each with different underlying assumptions and various degree of rigor.
In our case, it is of great help that a sufficiently general DW/RW system can be mapped onto a wave-mechanical system described by a Hamiltonian operator, akin to a Schrödinger equation in quantum mechanics, the only difference being that this operator turns out to be non-Hermitian. This non-Hermiticity reaffirms the fact that we are dealing with an open system [24]. By virtue of this flow-Schrödinger analogy, one can apply a quantumstatistical technique based on master equations with non-Hermitian Hamiltonians (NH), developed relatively recently [25][26][27][28][29][30][31]. This approach is proven to be robust for a large set of physical systems and phenomena, to mention just recent literature . The advantage of this approach is that it makes the whole theory of dissipative DW/RW systems a subset of the modern theory of open systems; which is statistical mechanical by nature, and describes both microscopical and macroscopic phenomena, as well as interplays between them from basic quantum-mechanical principles [64].
The outline of the paper is as follows. In Section 2, we give a brief introduction to the generalized Hasegawa-Mima model, and the underlying assumptions thereof. In Section 3, we formulate a mapping between the HME flow equations and wave equations of a Schrödinger type, which allows us to describe wave-mechanical phenomena in zonal flows. In Section 4, we generalize the wavemechanical approach from state vectors to their statistical ensembles by introducing density operator and master equations. We find that two types of time evolution are possible, described by non-normalized and normalized density operators, and describe differences between the two. In Sections 5 and 6, we formulate a phase space approach for finding solutions of master equations in, respectively, non-normalized and normalized density operator cases. Section 7 summarizes our findings and presents our conclusions.

HASEGAWA-MIMA MODEL
Let us derive a generalized Hasegawa-Mima equation as an example, which is widely used to describe electrostatic two-dimensional turbulent flows [18,19]. These flows can occur in magnetized plasma, which exhibits the drift-wave turbulence. The plasma itself is usually assumed to be collisionless and quasi-neutral, its magnetic field B being uniform, ion temperature being much smaller than the electron temperature, and electrons following Boltzmann distribution [65,66].
As mentioned in our introduction, such flows can also occur in atmospheric fluids on a rotating planet, where the role of drift waves is played by Rossby waves. We are going to describe both types of physical systems using the same approach, the only difference being the values of some parameters, therefore our results should be applicable to both types of waves.

Hasegawa-Mima equation for plasma
We assume a conventional geophysical coordinate system, where x = (x, y) are coordinates on a twodimensional plane, such that the x axis lies in the zonal flow direction, and the y-axis points in the direction of the local gradient of the plasma density (in atmospheric fluids, it would be the Coriolis force parameter).
Let us derive a generalized Hasegawa-Mima equation in the case of plasma being entrapped in the magnetic field, which is directed along z-axis: where B = B 0 is a constant magnitude and e z is a unit vector normal to the plane. Plasma is assumed to be confined to a slab, which is perpendicular to the magnetic field, so that perturbations of density and potential can propagate only in the (x, y) plane. It is inhomogeneous, with density gradient pointing along x-axis.
From now on, we assume the cold ions approximation T i ≪ T e . Electron and ion components' temperatures are assumed to be constant, therefore the electron density is n e = n 0 exp (eφ/T e ), where φ = φ 0 +φ being electric potential, n 0 is unperturbed density and T e is the electron temperature. If eφ ≪ T e , then n e ≈ n 0 (1 + eφ/T e ). Keeping perturbations of electron density up to a first order only, n e = n 0 +ñ e , we obtaiñ therefore, in the formula φ = φ 0 +φ the first term vanishes.
In a slab plasma with v z ≪ v ⊥ , with the electrostatic field being E = −∇φ, where ∇ = e x ∂ x + e y ∂ y , the equation of motion for an ion is where m i and q are ion's mass and charge, respectively. Applying a perturbation theory, we can write Eq. (3) as where ω ci = qB/m i is the ion's cyclotron frequency. To a zeroth order of perturbation theory, we thus obtain the equation: whose solution describes rotary motion, so that we have v (0) In the first order, one can write whose solution is given by v (1) where we denoted the convective derivative as where by ∂ t we denote a partial time derivative. Ignoring higher-order corrections, we obtain By substituting n i = n 0 +ñ i into the continuity equation for the ion component's fluid, and using the quasi-neutrality conditionñ i =ñ e , we obtain the following equation which transforms, using Eqs. (2) and (8), and neglecting terms which are higher than the first order, into the equation In this equation, the term (∇φ × e z ) · ∇φ vanishes due to orthogonality, while the terms with (eφ/T e )∇· can be neglected as being higher-order small. Multiplying Eq. (9) by −T e /e and recalling an expression for the ion's cyclotron (gyro)radius, where c s = T e /m i being an ion's characteristic sound speed, we obtain In terms of dimensionless valuest = ω ci t, (x,ý) = (x/R s , y/R s ) andφ = eφ/T e , we thus obtain where Q(x, t) is an added external forces' and dissipation term, β ≡ V * /c s is a positive constant parameter, V * ≡ −(T e /|e|B) ∂x ln n 0 , and w = (∇ 2 − 1)φ is a vorticity defined as a measure of local angular velocity, see Appendix A for details.
Placing both plasma and atmospheric flow cases on the same footing, one can write the generalized Hasegawa-Mima equation in the final form: where the generalized vorticity w(x, t) can be written in the form: where v = e z × ∇φ is the fluid velocity on the (x, y) plane,α is a preselected operator whose value equals the unit operatorÎ (DW case) or to zero (RW case), φ(x, t) is the electric potential or stream function, L D is the plasma sound radius or deformation radius. Here, and below, we work with the dimensionless values, but omit primes for brevity.

Fluctuations
Let us define zonal averages according to the formula where L x is the system's extent along x axis. From now on we assume that a tilde and bar refer to, respectively, fluctuation and zonal-averaged values. We therefore start with the following expansion [23]: φ =φ(x, y, t) + φ(y, t), Note that if we want to restore the parameter L D , we would obtaiñ Furthermore, Eq. (14) can be written in the quasilinear approximation as because eddy-eddy interactions can be omitted [13]. Fluid velocity can be written in terms of its components as v = −e x ∂ y φ + e y ∂ x φ, therefore, its fluctuating and zonal-averaged components are, respectively: is the xth component of a ZF velocity v = e x U . With the use of Eqs. (19), (23) and due to independence of φ from x, the second term of Eq. (21) can be rewritten as where we used the notation U k ≡ ∂ k U/∂y k = −∂ k+1 φ/∂y k+1 . A third term in Eq. (21) can be written as Furthermore, a second term in (22) yields where we used integration by parts in the last line. Substituting Eqs. (25), (26) and (27) into Eqs. (21) and (22), we obtain: whereζ is an external dissipation source with zero zonal average, µ dw and µ zf are the constant parameters of the terms describing the simplest kinds of dissipation of, respectively, driftons and zonal flows, caused by the environment. Using Eq. (19), these equations can be also rewritten in terms of the functionφ: whereξ is an external dissipative potential with zero zonal average [12].

Observables
The important observables of the model are enstrophy and energy, which can be defined as the following integrals: where subscripts indicate a corresponding component, either a drift-wave or a zonal-flow one. The total enstrophy and energy are, respectively which are usually expected to be conserved values in those physical systems, and which allow description in terms of wave kinetic equations (WKE).

WAVE-MECHANICAL ANALOGY IN ZONAL FLOWS
Building on the ideas presented in works [5,22,23], let us formulate a formal mapping of Eq. (30) to an effective Schrödinger-like equation which is somewhat analogous to a quantum-mechanical description. For reasons, which will be later specified, around Eq. (58), we shall refer to this analogy as wave-mechanical.
Judging by the form of Eq. (30), one can expect that its solutions form the Hilbert space of normalizable smooth functions, whose inner product is defined as an integral over the whole (x, y) plane, d 2 x. Postponing discussion of the physical implications of this until the end of this section, we introduce, in Dirac's bra-ket notations, a state vector |φ , such that in the coordinate representation one obtainsφ(x, t) = x|φ , similarly for ξ(x, t).
Using the resolution of identity, d 2 x|x x| =Î, we obtain where ℵ(t) is a norm of a state vector |φ . This state vector is almost what we are looking for, but the proper Schrödinger analogy requires a normalized state vector (ray), which spans a projective Hilbert space. We thus introduce so that Ψ|Ψ = 1. Furthermore, using an operator of spatial translations in the (x, y) plane,p = −i∇ = −i(e xpx + e ypy ), we can define the Hermitian In terms of these operators, Eq. (30) can be rewritten in the Schrödingerlike form: or, in a general basis where the Hamiltonian operator is given bŷ where we introduced a similarity transform ofÛ : where dot means an ordinary derivative with respect to time.
It should be noticed that the term i‫2(/א‬ℵ) occurs in the Hamiltonian, caused by the transition from |w to the proper state vector (38), which indicates that dissipation exists even if µ (0) dw is zero. In general, this term is time-dependent, except in two cases: when ℵ is constant (in which case the term vanishes) or when ℵ depends on time exponentially (then the term becomes a constant representing a rate of loss or gain). Note that because the canonical Schrödinger picture per se presumes no explicit time dependence of a Hamiltonian operator, time evolution of our system at the general function ℵ(t) must be described, not by a Schrödinger equation, but by a master equation as will be discussed later, in Section 4.1.
Note also, that in Ref. [23] the flow-Schrödinger analogy was proposed, but using wavefunctionω(x, t) instead. From Eq. (19), it is clear thatω is a surjective function ofφ (due to the partial differentiation thereof), which should rather be interpreted as a source density. All this suggests that the Hilbert space spanned by state vectors |φ , or by projective rays |Ψ , is the underlying one. Therefore, it is not surprising that by using theφ-associated Hilbert space, one obtains more general Hamiltonian, as will be demonstrated below, and also decreases the number of inverse differential operators in the resulting evolution equations, both for the state vector and function U .
Furthermore, in terms of the normalized state vector, some of the integral values of Section 2.3 can be written in a more convenient form to use in the next section. We thus obtain in Dirac's notations where we used Eq. (38). Furthermore, one can see that the Hermitian adjoint of an operator (41), does not coincide with Eq. (41). Therefore, our model belongs to a class of theories with non-Hermitian Hamiltonians, usually referred as NH theories. Furthermore, operator (41) can be decomposed into its Hermitian and anti-Hermitian parts: is a Hermitian operator, often referred as the decay rate operator. In our case, we obtainĤ and { , } and [ , ] are the anticommutator and commutator, respectively. Using these, Eq. (54) can be written in the form which is more convenient for our calculations. Let us turn our attention to the evolution equation for function U , which must be also rewritten in a wavemechanical form. From Eqs. (31) and (38) we obtain where an overscore denotes a zonal average per usual.
To conclude this section, we showed that drift waves can be described not only as waves but also as "quanta" (states in a Hilbert space) which is somewhat similar to the de Broglie's wave-particle duality in quantum mechanics. However, one must emphasize that this analogy between fluctuation equations in turbulent flows and a Schrödinger-type equation is not an exactly quantummechanical one. The reason for this is that the flow equation does not contain a Planck constant, but only an analogue thereof, which is determined by pertinent scales of length, time and mass. Up to a dimensionless coefficient, we can define this effective Planck constant as where we used magnetized plasma's characteristic scales of length R s , time ω −1 ci , and mass m i , taken from Section 2.1. In the case of atmospheric fluids one could use, respectively, the planet's average radius, angular speed of planet's rotation, and mass of a characteristic fluid parcel of planet's atmosphere. For the sake of brevity, the value of eff can be set to one, by rescaling and working in properly chosen units.
Nevertheless, this is yet another example of macroscopic phenomena which have quantum-like features and can thus be described by virtue of notions and methods originating from quantum mechanics. One of such formalisms is an analogue of quantum-statistical master equation approach in the theory of open quantum systems, which will be described in the next section. While Eqs. (39) and (40) are simply a way of rewriting the original flow equations, the next section's formalism is actually a generalization, which ushers in a new physics.

DENSITY OPERATOR FORMALISM
In this section, we describe statistical approach based on the density operator, which is analogous to the von Neumann approach in quantum mechanics of mixed states, i.e., probabilistic mixtures of pure states [25][26][27][28][29][30]. In the matrix representation of a density operator, pure states are described by diagonal elements of a density matrix, whereas mixed states -by off-diagonal components thereof. It is the latter which are regarded in quantum theory as an essentially quantum-mechanical effect, which does not usually have a classical counterpart.
When it comes to NH systems, the density operator formalism is more general than the state-vector one, because it allows the description of not only pure states, or their superpositions, but also mixed states. This generalization is especially important, because dissipative effects are known to evolve some pure states into mixed ones [29]. In other words, restricting ourselves to state vectors ab initio would pose an implicit assumption of forbidding the transition from pure to mixed states. This assumption does not seem to be realistic when it comes to open systems.
Besides, if one were dealing with an equation for state vectors driven by a non-Hermitian Hamilton operator, such as Eq. (40), then one would arrive at complex eigenvalues of energy, which seem somewhat contradictory to the notion of energy per se. On the contrary, in a density operator approach, no complex-valued energies occur: the role of anti-Hermitian part of the Hamiltonian is to describe the time evolution of energy eigenvalues of its Hermitian part.

Master equation: Non-sustainable evolution
Let us regard the operator |Ψ Ψ| as a special case of the reduced density operatorŴ . From Eq. (40) and its adjoint, one can obtain an equation forŴ , which is thus called a master equation: is a dissipator operator, which is self-adjoint by construction. Here, the arrow stands for "for pure states tends to", which indicates thatD is a mixed-state generalization of the operator |ξ Ψ| + |Ψ ξ |. In that case, one should use the hybrid approach, which deals with systems which are driven by both Liouvillian (or Lindblad) and non-Hermitian Hamiltonian parts of their master equations [26].
If we neglect the dissipator term for the sake of simplicity, we arrive at a canonical NH master equation where the NH operators are given by Eqs. (51)- (53). Physical observables in this approach are defined as statistical averages: for a given observable's operatorÂ, and similarly for correlation functions [27]. Let us recall Eq. (57), which is supposed to supplement our master equation. Unlike the master equation, it is not an operator equation, therefore U can be considered a c-number. Altogether we obtain where we denotedŴ ℵ = ℵŴ . In the last equation, we introduced the partial trace operation 'tr', which is an algebraic matrix tracing or summation over all states of the type |Ψ k Ψ k |, but without averaging (integrating) over configuration space.
One can see that the ℵ-dependent term in Eq. (52) can be removed by rescaling the density operator. Note that the averages are still computed with respect toŴ , therefore Â W = ℵ −1 Tr(ÂŴ ℵ ) = ℵ −1 Â W ℵ , according to Eq. (62). Recalling that the typical averages in our case, such as (44) and (46), have a coefficient ℵ, we can transform our density operatorŴ →Ŵ ℵ = ℵŴ , Â W → Â W ℵ = ℵ Â W , then drop the subscript 'ℵ', and assume ℵ → 1 from now on. This contributes to the remarks made after Eq. (43): the ℵ-dependent term in the Hamiltonian disappears when working within the framework of a more general approach.
We thus finally obtain the following set of evolution equations while the averages are defined according to Eq. (62). For example, using Eqs. (44)-(47), we can define enstrophy and entropy as the following statistical averages: while the total values are given by Eqs. (36), per usual. Finally, taking trace of Eq. (65), one can obtain the following equatioṅ supplemented with the initial condition T W (t = t 0 ) = 1, where T W (t) ≡ TrŴ = Î W , and t 0 is the moment of time at which the environment switches on. This equation clearly indicates that the trace of the operatorŴ is not necessarily conserved during evolution. This means that the drifton (sub)system experiences drain or loss of its degrees of freedom, which can result in its total decay or critical instability. While this can indeed be the case in some systems, it is not compulsory for all open systems. In fact, next we are going to consider the possibility of other type of evolution.

Master equation: Sustainable evolution
Dynamical systems, which can be described by non-Hermitian Hamiltonians, can follow two types of evolution; which can be referred as non-sustainable and sustainable, by analogy with some photobiological systems, where this difference becomes striking [33]. The former type is the one described by the densityŴ obeying the master equation (61). The sustainable type is the one described by the normalized density operator so that Tr(ρ) = 1 at all times, which automatically removes the probability gain/loss problem [25]. This density would be a solution of the equation which is both nonlocal and nonlinear with respect to a density operator. Fortunately, for practical computations, this equation can be transformed into a linear equation of type (65) by using the ansatz (72). For the sustainable type of evolution, physical observables are defined as for a given observable's operatorÂ, and similarly for correlation functions [27]. Using the ansatz (72), one can relate the two types of averages: which is often convenient for optimizing computations. Furthermore, recalling the mean-value origin of the last term in Eq. (29), the operatorŴ must be replaced witĥ ρ in the equation for the function U : according to remarks preceding Eq. (74). Similarly to Eqs. (67)-(70), one can define "normalized" enstrophy and entropy as the following statistical averages:Z where U is now a solution of Eq. (76). While these definitions look very similar to their non-sustainable evolution analogues (67)- (70), one should notice the differences; such as the coefficients, functions of time, which occur both in the definitions themselves and in Eq. (76). To summarize Sections 4.1 and 4.2: in the case of non-Hermitian Hamiltonian systems, one can have two types of evolution, each described by its own sets of equations. In the case of non-sustainable evolution, this set consists of Eqs. (65) and (66), whereas in the case of sustainable evolution -of Eqs. (73) and (76). The choice between these types of NH evolution, or even the switch between them, is a process which is external to the subsystem itself (it should also be remembered that we are dealing with reduced density operators here). In other words, this choice must be considered a special case of the effect induced by the environment, which can also include the measuring apparatus itself, as well as any pre-and post-selection measurement protocols. For example, in quantum optics, the calibrating and resetting of fields inside optical fibers, to prevent them from damage or to prepare for the next run, would be another example of sustainable-type evolution.

Time evolution of averages
Master equations are differential equations for operators, which makes solving them a technically challenging task. However, some physical information can be extracted without actually finding a density operator explicitly. One of methods of such an extraction are equations for averages, which can be derived from a master equation.
Let us consider the non-sustainable type of NH evolution, described by Eqs. (65) and (66). Multiplying Eq. (65) with various operators and tracing, one can obtain time evolution equations for corresponding averages as ordinary differential equations. For example: where we used identities mentioned after Eq. (41) and the identity {ÂB,Ĉ} =Â[B,Ĉ]+ [Â,Ĉ]B. If some of the resulting differential equations form a nontrivial closed subset, they can then be solved for the corresponding averages as functions of time. Otherwise, one has to apply approximations or iteration methods to the full set.

PHASE-SPACE FORMULATION: NON-SUSTAINABLE EVOLUTION
Master equations of the types described in the previous section are differential equations for operators. This creates a formidable technical problem when it comes to solving them. As a result, there are a number of methods which allow us to solve them, each with a different degree of accuracy or completeness. Solving the equations for averages is one of the approaches -it allows us to compute averages as solutions of a set of differential equations for functions, instead of dealing with operators. However, in the case of an infinite-dimensional Hilbert space, the number of such equations is not necessarily finite or reducible to something simple, but dependent on other symmetries of the problem. The second method, which is to be discussed in this section, is based on mapping between functions in the quantum phase space formulation and Hilbert space operators in the Schrödinger representation [67][68][69][70], which was adapted for density operators for various NH systems [71][72][73][74], including zonal-flow models specifically [5,22,23].
There also exists a subtle point when it comes to the completeness of the phase space description of density operators for NH systems including zonal-flow ones. As mentioned at the beginning of Section 4, it is important to preserve and extract, as fully as possible, an essentially non-classical piece of information about time evolution of a probabilistic mixture of states of NH-driven systems [29]. This is straightforward within the framework of the Hilbert space matrix formulation of a density operator from Section 4, where the structure is clear: this information is encoded mostly in the off-diagonal elements of a density matrix, as we know it from quantum mechanics of mixed states. It is far less obvious how one can separate "off-diagonal" effects from "diagonal" ones in the phase space formulation of a density operator, in presence of additional approximations which are technically inevitable to obtain definite results.

Wigner-Weyl transform
Let us introduce the Weyl symbol of the density oper-atorŴ as the integral [69]: from which the following property can be easily deduced in the case of real-valued wavefunctions.
Correspondingly, Weyl symbols for operators in the previous section will turn into functions (c-numbers). For example, using the resolution of identity, we obtain for the operatorp x : for the operatorp −2 D : for the similarity transform operator and its combinations, using Eqs. (54)-(56): where {{, }} and {{{, }}} are, respectively, sine and cosine Moyal brackets defined in Appendix B. Furthermore, using the Moyal product rule (A16) we obtain for operator products, such as W(p xp Other properties of the Moyal product, which will be used in what follows, such as the associativity (A19), are listed in Appendix B.
Performing the Weyl transform of Eq. (65), we obtain the equation where where we took into account in the last step that the zonalaveraged Wigner function W does not depend on x and satisfies the condition (87). Finally, the integrals (67)- (70) are transformed into the form and the total values are as defined in Eq. (36).

Eikonal approximation
Equation (95) remains difficult to solve, except in a few cases when the Moyal products series can be truncated. In our case, it seems that no truncation is possible, therefore one has to resort to the eikonal or geometrical-optics approximation (A25), which is analogous to a leadingorder WKB approximation with respect to the effective Planck constant (58). Then Eqs. (93) and (94) become approximately and the whole system (95)-(98) simplifies to a system of integro-differential equations where where {, } c is a canonical Poisson bracket as defined in Appendix B.
Notice that in the formula for H, all derivatives of U are even, while for G -all are odd, which is similar to even-odd function classification in quantum mechanics. In this case, this is caused by a discrete Z 2 symmetry with respect to the mirror transformation x → x, y → −y, which leaves invariant both H and G. This symmetry thus supplements the parity-time reversal symmetry x → −x, t → −t, of the system (105)-(108), which occurs when µ By zonal averaging of Eqs. (105) and (106), we obtain the following equations where we also used Eqs. (43) and (121). The former of these equations belongs to a class of generalized WKE models, some other examples to be found in [23]. Finally, in the eikonal approximation, the integrals (99)-(102) can be simplified to while the total values are as defined in Eq. (36), as before. The rates of these values are where we used evolution equations (109), (110), and formulae from Appendix B.
Using these expressions for rates, one can immediately see that total enstrophy and total energy (36) evolve in our model according to the formulae therefore, depending on values of µ's, our system's enstrophy and energy can evolve in different ways. Let us consider the following special cases: • Overall conservation occurs when In this case, both energy and enstrophy can flow from the drift-wave component to the zonal-flow one, and back, but do not leave the system: • Overall exponential gain or loss occurs when where λ is a real-valued constant. In this case, the total energy and enstrophy of our system both increase (decrease) if λ is negative (positive), at an exponential rate: which ultimately leads to either critical instability or to the complete depletion of the system. Thus, our system has a finite lifetime τ = 1/λ in this case.
To conclude, an eikonal approximation reveals a number of important features of dissipative processes induced by drift-wave turbulence, including the behaviour of average values of energy and vorticity of its drifton and zonal-flow components.

PHASE-SPACE FORMULATION: SUSTAINABLE EVOLUTION
In this section, we adopt a Wigner-Weyl formalism in the case of the normalized density operatorρ. In the computations that follow, an ansatz (72) will be of great assistance, because it offers us a shortcut for our calculations, then use the results related to non-sustainable evolution.

Wigner-Weyl transform
Let us introduce the Weyl symbol of the normalized density operatorρ as the ratio [75]: where and W(Ŵ ) is given by Eq. (86). Therefore, an ansatz (72) can be applied in the form whereẀ (x, p, t) is a solution of the equation wherè wherè where functionÙ , according to Eqs. (76) and (129), is a solution of the equation whereẀ is a zonal average of the solution of Eq. (131). Furthermore, averages (77)-(80) take the form where we used the ansatz (129).

Eikonal approximation
By analogy with Section 5.2, we obtain that Eqs. (134) and (135) become in a leading-order approximatioǹ therefore, the system (131)-(136) simplifies to a system of integro-differential equations Therefore, performing zonal averaging of Eqs. (143), we obtain the zonal-averaged master equation while Eq. (144) stands as is.
It is useful to compute the trace of the auxiliary density operator in this approximation. Using Eqs. (71), (129), (130), (146) and (147), we obtaiṅ where is an average of the truncated decay rate operator, dwÎ , with respect to the auxiliary (non-normalized) density operator.
Finally, in the eikonal approximation, averages (137)-(140) can be simplified tò therefore, rates of these values can be computed, using the evolution equations (144), (147), (148), and formulae from Appendix B, as where we used G ρ (t) ≡ G W /TẀ to denote an average of the truncated decay rate operator with respect to the main (normalized) density operator.
Using these formulae, one can see that "normalized" total enstrophy and total energy (36) evolve in our model according to the formulae where we used the total values defined by Eqs. (36) as analogy, per usual. Notice that, unlike their analogues from Section 5.2, these rates do not depend on µ (0) dw .

CONCLUSION
In this paper, we have presented a statistical mechanical approach to describing dissipative phenomena in zonal flows of plasmas and atmospheric fluids, such as drift waves and Rossby waves, which is based on Landau-von Neumann's density operator in a theory of open quantum systems. This became possible due to an occurrence of Hilbert space associated with an electric potential or stream function, which is regarded as the fundamental Hilbert space of the theory. This results in a formal mapping between flow and wave equations, which allows us to describe zonal flows and associated phenomena as macroscopic wave-mechanical effects. As a consequence of this mapping, flow equations can be rewritten as Schrödinger-like equations, with two important differences: for dimensionality purposes, one uses an effective Planck constant whose value is not necessarily equal to the quantum-mechanical Planck constant, and the resulting Hamiltonian operator is not necessarily Hermitian.
The second feature requires us to treat the entire theory from within the framework of the non-Hermitian Hamiltonian approach, where anti-Hermitian parts of Hamiltonian operators usually describe an effect of the environment. Fortunately, such an approach has been already developed and applied to various open systems. According to this formalism, one has to generalize from state vectors to density operators, because pure states do not necessarily stay pure during the time evolution in the presence of dissipation and noise. Moreover, the density operator approach prevents the occurrence of complexvalued energies. In this approach, energies are always real-valued, and refer to the subsystem itself, whereas any anti-Hermitian components describe the effects of the environment upon this subsystem (decay rates, et cetera).
Thus, after deriving a Hamiltonian operator and evolution equations for state vectors in our DW/RW models, we made a transition from state vectors to density operators, introduced NH master equations and defined observables; including the enstrophy and energy of both the waves and zonal flow. The upshot is that two types of density operator's evolution exist, these can be referred to as non-sustainable and sustainable, by analogy with some photobiological systems where they can be visualized.
The non-sustainable type is the one described by the non-normalized density operator. During such evolution, an open (sub)system experiences drain or loss of its degrees of freedom, which can result in its total decay or critical instability. While this can indeed be the case in some systems, it is not a compulsory feature of all open systems.
The sustainable type is the one described by the density operator which is normalized at all-times. This automatically removes the problem of probability's gain or loss in NH systems, thus enhancing their stability. An example of such stability would be an environment-assisted stability in photobiological systems. Consequently, various observables, including those related to enstrophy and energy, behave in a way different from the nonsustainable case.
Because we are dealing with reduced density operators; a selection of one or another type of NH evolution, or even the switch between them at some point in time, is a process external to the subsystem itself. In fact, this process is neither unitary nor continuous. It can be considered a special case of the effect induced by the environment, which can also include the measuring apparatus itself.
Furthermore, to establish a method for solving master equations for a given DW/RW model, we considered a phase-space formulation of the theory. We introduced relevant Weyl-Wigner transforms and rewrote evolution equations and observables in the Moyal form. We also studied a leading-order approximation of the Wigner approach, for both types of evolution, which is analogous to the eikonal or geometrical approximation in optics or WKB approximation in quantum mechanics.
As it turns out, the statistical-mechanical density operator formalism shows itself to be an approach to dissipative phenomena related to zonal flows, which has clear foundations and notions emanating from quantum mechanics. It also allows us to take into account the wave-mechanical effects in the above-mentioned systems, thus producing more realistic descriptions thereof. As such, the formalism can also be extended to other dissipative systems in plasma and atmospheric physics, which allow the wave-mechanical analogy and non-Hermitian Hamiltonian operator description.
When assuming an eikonal or leading-order WKB approximation, these brackets have the following properties where the effective Planck constant is defined in Eq. (58). Finally, let us give some useful formulae for our case. If configuration space is two-dimensional and phase-space functions A and B do not depend on a coordinate x, we obtain where D(∂) are total derivative terms with respect to arguments listed in braces. Such terms can be omitted when working inside reduced phase space integrals: d 2 p dy D(∂) = 0, assuming that physical values vanish at phase space borders. For instance, for the functions we obtain the following identity: where n being integer, hence because neither U nor W depend on x.