Zubarev ’ s method of a nonequilibrium statistical operator and some challenges in the theory of irreversible processes

Once Gibbs announced his method of statistical ensembles for equilibrium many-particle systems, a serious effort was mounted to extend this method to nonequilibrium situations. Nevertheless, it was a long time before such an extension was developed in a form general enough to study virtually any class of nonequilibrium processes. One of the most compact and elegant approaches of this kind is Zubarev’s method of a nonequilibrium statistical operator (NSOM). This method is detailed in books [1,2] (see also Zubarev’s article [3] where the history of NSOM is presented). The Zubarev formulation of nonequilibrium statistical mechanics turned out to be very convenient for concrete applications and at the moment there exists extensive literature which shows NSOM in action. Here we will not dwell on technical aspects and concrete applications of Zubarev’s method. The reader may wish to refer to papers in this issue and in special literature. On the other hand, since nonequilibrium statistical mechanics is on a stage of active development, it makes sense to take a quick look at more general questions related to NSOM, which may be considered as challenges and call for further investigation.


Introduction
Once Gibbs announced his method of statistical ensembles for equilibrium many-particle systems, a serious effort was mounted to extend this method to nonequilibrium situations.Nevertheless, it was a long time before such an extension was developed in a form general enough to study virtually any class of nonequilibrium processes.One of the most compact and elegant approaches of this kind is Zubarev's method of a nonequilibrium statistical operator (NSOM).This method is detailed in books [1,2] (see also Zubarev's article [3] where the history of NSOM is presented).The Zubarev formulation of nonequilibrium statistical mechanics turned out to be very convenient for concrete applications and at the moment there exists extensive literature which shows NSOM in action.Here we will not dwell on technical aspects and concrete applications of Zubarev's method.The reader may wish to refer to papers in this issue and in special literature.On the other hand, since nonequilibrium statistical mechanics is on a stage of active development, it makes sense to take a quick look at more general questions related to NSOM, which may be considered as challenges and call for further investigation.Some of these questions were discussed with D. N. Zubarev during our work on books [2,4].

Nonequilibrium statistical ensembles and irreversible thermodynamics
It is a matter of common knowledge that statistical mechanics must provide microscopic foundations for the thermodynamic description of many-particle systems.For equilibrium, the situation seems to be wholly satisfactory since Gibbsian statistical ensembles give rigorous formulations of thermodynamic quantities and relations.Following this line of reasoning, it is hoped that an extension of Gibbs' method to the nonequilibrium case would justify the basic postulates of irreversible thermodynamics.On this way, however, one meets a number of fundamental problems.Some of them are yet to be solved.Here we will briefly discuss this point in the context of NSOM.
For orientation, we first recall the basic ideas of NSOM which are required in further considerations.In Zubarev's approach, the starting point for the construction of the nonequilibrium statistical operator1 ̺(t) is the Liouville equation with an infinitesimally small source term that breaks the time reversal symmetry of the equation [1,2]: where L is the Liouville operator determined by the quantum or the classical Poisson bracket with the Hamiltonian Ĥ: iL̺(t) = [̺(t), Ĥ]/i (the quantum case) and iL̺(t) = {̺(t), Ĥ} (the classical case).The source term on the right-hand side of equation (2.1) selects a retarded solution of the Liouville equation, which coincides with an auxiliary statistical operator ̺ rel (t) in a distant past.This relevant statistical operator (another frequently used term is a quasi-equilibrium statistical operator ) is intimately related to the thermodynamic description of the system.
Here we would like to emphasize the trivial but sometimes overlooked fact that the Liouville equation, being a differential equation, does not uniquely determine the statistical operator, so that one has to introduce some initial or boundary conditions describing the physical situation.
A standard procedure for constructing ̺ rel (t) is based on the extremum of the information entropy for the given average values of some relevant dynamical variables { Pm } describing the macroscopic state of the system [5].In many concrete problems one deals with local variables Pi (r) which depend on space.In such cases, the composite index m = {i, r} contains spatial coordinates of the point, and the summation over m includes integration over r.
The extremum condition for the information entropy leads to the following explicit expression for the relevant statistical operator: where F m (t) are Lagrange multipliers determined by the self-consistency conditions for the given values P m (t).The quantity is called the Massieu-Planck function (functional).
Starting from equation (2.1) with the relevant statistical operator given by equation ( 2.3), one can derive a formally closed set of the evolution equations for the observables P m (t).The corresponding mathematics is presented, for example, in book [2].
Formulation of the boundary condition for ̺(t) in terms of the relevant statistical operator is just the point where a thermodynamic description is introduced.It can easily be shown that the Massieu-Planck function (2.5) can be regarded as a nonequilibrium thermodynamic potential in the variables {F m (t)}.Indeed, varying the both sides of equation (2.5) with respect to F m (t) and then taking the self-consistency conditions (2.4) into account, we obtain the thermodynamic relations We can also go to the description in the variables {P m (t)} by means of a nonequilibrium generalization of the Legendre transformation.The potential conjugated to Φ(t) is the nonequilibrium thermodynamic entropy S(t) which is defined as the information entropy in the relevant statistical ensemble, i.e.

S(t)
where again self-consistency conditions have been used.This equation is well known from irreversible thermodynamics [6].It gives the entropy production in terms of the thermodynamic forces F m (t) and the macroscopic fluxes ∂P m (t)/∂t.The crucial point in the above formalism is the choice of the set of relevant variables { Pm } which is adequate for the description of a nonequilibrium system.One frequently used guide-line for choosing relevant variables is the so-called "hierarchy" of widely spaced relaxation times in many-particle systems.By this is meant that nonequilibrium processes typically proceed on different stages, each characterized by a proper time scale.Different time scales can be determined by the characteristic time intervals ∆t at which the change of the system state occurs.The larger is the characteristic time interval, the smaller set of slowly varying observables can be used.This idea of a reduced description of nonequilibrium systems was first formulated and developed by Bogoliubov [7] in the context of the classical kinetic theory and now is, in fact, the heart of all the methods in nonequilibrium statistical mechanics.
In Zubarev's and many other approaches, the level of the reduced description is determined by the set of dynamical variables { Pm } in the relevant statistical operator.For example, the usual hydrodynamic description implies that the variables Pm ≡ Pj (r) correspond to locally conserved quantities, such as the densities of energy, mass, and momentum.It should be emphasized, however, that, in principle, nothing prevents us from using variables which are not "slow" variables on a chosen time scale.An example is the Grad method in the classical kinetic theory [8], where the fluxes of hydrodynamic quantities are treated as independent variables in constructing normal solutions of the Boltzmann equation.
Note that, depending on the choice of relevant variables, we have in fact different thermodynamic descriptions of the same nonequilibrium system.However, this "versatility" of the formalism brings up some questions.First, what is the physical meaning of thermodynamic parameters F m for different sets of relevant variables ?Second, what is the relation between the evolution equations for different sets of relevant variables and how these equations are related to the observable quantities ?Third, is there a general rule for the choice of Pm if one uses slow variables as well as fast variables ?
At present we know enough to solve the above questions only for linear irreversible processes when the state of a system is close to total equilibrium.It can be shown (for a detailed discussion see, e.g., [4]) that in the limit ε → +0 all sets of relevant variables are equivalent in the sense that they give the same results for observable quantities (transport coefficients and susceptibilities), if one deals with exact expressions.On the other hand, using the same approximation for dif-ferent sets of relevant variables, say the Born approximation in time correlation functions, different results for transport coefficients and susceptibilities can be obtained.Roughly speaking, the more dynamical variables { Pm } are included into the relevant statistical operator, the closer we get to the exact results for transport coefficients and susceptibilities.Thus it appears reasonable to consider different sets { Pm } as trial sets of variables within some variational principle.This interpretation can be justified rigorously in the special case of stationary processes [4].The physical meaning of the variational principle is very simple: as long as the parameter ε in equation (2.1) is kept finite, the exact values of transport coefficients and susceptibilities correspond to the maximum of the entropy production in the system.This resembles the situation in the kinetic theory of rarefied gases [9] where the exact solution of the Chapman-Enskog equations gives transport coefficients which correspond to the maximum of the entropy production (the Kohler variational principle).For non-stationary linear irreversible processes, the relation between NSOM and variational principles calls for further investigation.
It goes without saying that in far-from-equilibrium situations the problem of the choice of proper relevant variables becomes much more complicated.At the phenomenological level, different sets of state variables are introduced in the socalled extended irreversible thermodynamics [10,11] which now is established as a formal closed theory.Recently a serious effort was made to justify the extended irreversible thermodynamics on the basis of general principles of statistical mechanics and, in particular, on the basis of NSOM (for a detailed discussion and references see, e.g., [12,13]).
A rather general scheme for the selection of relevant variables can be formulated [12,14] in the case when the total Hamiltonian admits a decomposition into two parts, Ĥ = Ĥ0 + Ĥ′ , where the secular part Ĥ0 includes kinetic energies of the particles and, may be, some strong interactions associated with very short relaxation times on the chosen time scale.The other part, Ĥ′ , is considered as a small perturbation describing slow relaxation processes in the system.Then it is assumed that appropriate variables Pm satisfy a closure condition (in the quantum case) where Ω mn are some c-numbers.This condition plays an important role in the kinetic theory (see, e.g., [1,2,15,16]) where Ĥ0 describes noninteracting particles (or quasiparticles), while Ĥ′ is an interaction term in the Hamiltonian.In the context of the kinetic theory, the relevant dynamical variables, Pm = a † l a l ′ , are bilinear in the creation and annihilation operators corresponding to some oneparticle quantum states |l .
If the variables Pm and (or) the secular part Ĥ0 of the Hamiltonian are not bilinear in the creation and annihilation operators, condition (2.11) produces, in general, an infinite chain of variables coupled through their commutator with Ĥ0 .Therefore, in concrete calculations one is led to apply a truncation procedure [14] which resembles the variational methods of the kinetic theory but is difficult to formulate explicitly.Note also that in many cases of interest the decomposition of the Hamiltonian into the secular part and the interaction term is not unique since this decomposition can depend on the order of the approximation (see, e.g., [17]).
In summary one may conclude that now the choice of relevant variables for highly nonequilibrium systems is based essentially on heuristic arguments supported by explicit calculations for some concrete models.The most natural way to a more general scheme is to look for a variational principle applicable to a wide range of far-from-equilibrium situations.We believe that the works along this line are of great importance for the development of NSOM.

Correlations in the kinetic theory
The kinetic theory is a traditional field of application of the nonequilibrium statistical mechanics.In a broad sense, any problem where the evolution of the single-particle distribution f 1 (t) (the classical distribution function or the quantum density matrix) is of interest may be assigned to the kinetic theory.Over many years it has been agreed that the main objective of the kinetic theory is to derive a closed kinetic equation for f 1 (t).In the context of NSOM this means that the dynamical variable f1 determined by the relation f 1 (t) = f1 t is the only variable Pm in the relevant statistical operator (2.3).Physically, the corresponding boundary condition for the Liouville equation is a complete uncoupling of all many-particle distribution functions to one-particle functions in the distant past.This is nothing but Bogoliubov's boundary condition in the kinetic theory of dilute gases [7].It is clear that for moderately dense gases and liquids such a boundary condition is not adequate since the relevant statistical operator describing an ideal gas differs greatly from the true nonequilibrium statistical operator which must involve many-particle correlations.These correlations are associated, for instance, with the bound states of particles, conservation laws, etc.The shortcoming of Bogoliubov's boundary condition for correlated systems is manifested in strong memory effects in the kinetic equation.Several modifications of Bogoliubov's boundary conditions have been proposed to include relevant correlations, thus preserving the Markovian character of the kinetic equation.Unfortunately, for the most part these modifications are based on heuristic considerations adapted to special cases.It should be noted, however, that in changing Bogoliubov's boundary condition it would be reasonable to proceed from a physical principle applicable to as many real systems as possible.Here we shall discuss the arguments put forward in [18].
The existence of long-lived correlations in a system is connected with the fact that there are some collective dynamical variables Ĉj which vary slowly in time.The locally conserved hydrodynamic variables may serve as an example.It is significant that Ĉj cannot all be expressed exactly in terms of a one-particle distribution (say, the energy density).Thus, in addition to f1 , one has from the beginning to include Ĉj in the set of relevant variables.This simple reasoning is conceptually in line with the Bogoliubov principle of the reduced description of nonequilibrium systems [7].The only new aspect is that we consider a sequence of decreasing time scales rather than a sequence of increasing ones, as it was done by Bogoliubov.
To clarify the meaning of this statement, we assume that the system of interest possesses a hierarchy of relaxation times τ i , (i = 0, 1, . ..), and τ i ≪ τ j+1 .The shortest relaxation time, τ 0 , defines the so-called dynamic stage of evolution when no reduced description can be used (for gases τ 0 can be estimated as the collision time).Let M (i) = { P (i) m } (i 1) be a set of relevant dynamical variables which are slowly varying variables on the time scale with the characteristic time interval ∆t i such that τ i ≪ ∆t i ≪ τ i+1 .Then the Bogoliubov formulation of the principle of a reduced description is that M (i+1) is a subset of M (i) , and that M (1) is sufficient for the most detailed description of the system.Stated more precisely, it is assumed that the equations of motion for P (1) m t can be taken as basic evolution equations for all the time scales.This is the case in the kinetic theory of dilute gases, where P m t may be interpreted as a total set of moments of the single-particle distribution function f 1 (t), and the Boltzmann equation plays the role of the basic evolution equation.Then, if the characteristic time interval ∆t is taken to be much larger than the mean free time of a particle, one can use only a few moments of f 1 (t) (the hydrodynamic variables) to describe the system.Unfortunately, even for moderate gases we do not know a closed kinetic equation which is applicable to all the time scales, so that we have to start directly from the N-particle Liouville equation and the above scheme of a reduced description becomes ineffective.Let us, however, consider the characteristic time scales in the reverse order.The largest time scale is determined by the characteristic interval ∆t τ eq , where τ eq is the time necessary for the relaxation of the system to total equilibrium.On this scale, only the average values of the integrals of motion, Ĉj , are needed to describe the system.The hydrodynamic time scale is associated with the characteristic interval such that τ loc ≪ ∆t ≪ τ eq , where τ loc is the time for the establishment of local equilibrium in a volume that is small but still contains a large number of particles.Now the state of the system can be described by the averaged densities of the integrals of motion (or by their spatial Fourier components Ĉj,k t ).Note that a new set of the relevant dynamical variables { Ĉj,k } includes all the variables Ĉj = Ĉj,k=0 of the previous set.This rule is quite general since any variable which is slowly varying on some scale of time is also slowly varying on shorter scales.
Generally speaking, the set of slowly varying variables on the hydrodynamic scale is not exhausted with Ĉj,k , because an arbitrary high order product, like Ĉj 1 ,k 1 Ĉj 2 ,k 2 • • • Ĉjs,ks , is also a slowly varying variable.Therefore, to obtain the Markovian equations of evolution, hydrodynamic fluctuations should be included into the set of relevant variables.In doing so, one can construct the thermodynamics of fluctuations [19] and derive the hydrodynamic Fokker-Planck equation for the distribution functional [20].If there exist other slow relaxation processes in the system (for instance, the formation of long-lived bound states or the relaxation of internal degrees of freedom of the particles), the corresponding dynamical variables should be incorporated along with the hydrodynamic variables.As one goes to shorter time scales, some new dynamical variables can also be taken to construct the relevant statistical operator determining the boundary condition for the true nonequilibrium statistical operator of the system.
The above procedure has no restriction to its generality and does not require that a basic kinetic equation like the Boltzmann equation be known, as it is necessary when one follows the original Bogoliubov scheme of the reduced description.The remaining problem, of course, is the choice of a proper set of relevant variables on small time scales.
Turning back to the kinetic theory, we note that the characteristic relaxation time for the single-particle distribution is distinctly smaller than the hydrodynamic relaxation time τ loc .Thus, strictly speaking, kinetic processes must be considered together with the evolution of the locally conserved quantities and other "slow variables".This is not clearly seen in the case of dilute gases since all the quantities of interest can be expressed approximately terms of the single-particle distribution.For dense and strongly correlated systems, however, the interplay between the kinetic and long-lived collective modes is of the utmost significance.Linear irreversible processes in classical fluids [21][22][23] and plasmas [24] are the examples.
The inclusion of long-lived correlations in the kinetic theory fits naturally into the scheme based on NSOM.Having constructed the relevant statistical operator with the kinetic and collective dynamical variables, one can use equation (2.1) to obtain the corresponding boundary conditions for the classical or quantum hierarchy.For a more detailed discussion of this point see our paper in this issue [25] and [2,4].

Nonequilibrium statistical mechanics of relativistic systems
Recent experimental studies of the ultra-relativistic heavy-ion collisions have attracted considerable attention to the theory of nonequilibrium processes in relativistic quantum systems.It is notable in this connection that one of the first Zubarev's papers on the method of the nonequilibrium statistical operator was devoted to the relativistic transport theory [26] (see also [1]).In this short paper Zubarev derived Green-Kubo formulae for transport coefficients and paved the way for nonequilibrium statistical mechanics of quantum fields.Later, a covariant formulation of the method was developed [27,28] and transport coefficients were calculated for some models [29].
The existing version of the relativistic NSOM works well in the case when disturbances are slowly varying on space and time scales associated with the microscopic correlation length and correlation time (for low-density systems these are the mean free path and the mean collision time).The Markovian evolution equations derived within the statistical approach coincide with the equations of the phenomenological relativistic hydrodynamics [30].If one attempts to study the processes involving steep gradients and rapid variations in time, the effects of nonlocality and memory should be taken into account.Note that in the relativistic context these effects have a fundamental significance.The point is that the customary (Markovian) equations of relativistic hydrodynamics are parabolic differential equations and, therefore, predict infinite speed of propagation of thermal disturbances in contradiction with the principle of causality.Formally, the generalized transport equations derived within NSOM [26,27] or other equivalent approaches [31] are nonlocal in space and time, but the structure of the transport coefficients is very complicated to be used in concrete problems.So, the task is to construct a suitable perturbation scheme for evaluating non-Markovian corrections to the hydrodynamic equations.
Another possible way to improve the customary relativistic hydrodynamics on shorter time scales is to use an extended set of variables describing the system.An analogous phenomenological version of irreversible relativistic thermodynamics was proposed by Israel [32].It is based on the assumption that the entropy density contains extra terms which depend on the heat flux and viscous stresses.Such a theory is in line with the molecular hydrodynamics developed for classical fluids [21] and predicts a wave-like propagation of thermal modes, which does not violate causality.Israel's extended relativistic thermodynamics was then confirmed by a generalized Grad method applied to the relativistic Boltzmann equation [33,34].Note, however, that the relativistic Boltzmann equation is inadequate for many of the field models, so that it might be good to use a more fundamental evolution equation.In the context of NSOM, one again faces the problem of the choice of a proper set of relevant variables on space and time scales, where relaxation processes are considered on the same footing as hydrodynamic ones.
We see that, setting aside questions specific to the Lorentz symmetry, etc., further development of the nonequilibrium relativistic statistical mechanics is closely related to the solution of the problems discussed in section 2. What is needed is a sufficiently justified procedure for constructing relevant statistical operators on short space and time scales.

Connection between NSOM and some other methods
As a final topic of our brief review, we compare NSOM with some other methods based on "first principles" of statistical mechanics.
At the level of a linear response of equilibrium systems to mechanical perturbations, the first general theory was developed in the pioneering work by Kubo [35].In Kubo's theory, it is assumed that in the distant past the system is in equilibrium with a heat bath and then its evolution is described by the Liouville equation with the Hamiltonian involving external mechanical perturbations.Within the framework of NSOM, Kubo's response formulae are obtained by setting ̺ rel (t) = ̺ eq , where ̺ eq is an equilibrium statistical operator [4].In a certain sense Kubo's scheme may be considered as a particular case of NSOM when the set of relevant variables contains only the integrals of motion.It should be emphasized that Kubo's theory describes the response of an isolated system to external mechanical perturbations.In general, this response need not be equal to the response of a system being in contact with a heat bath during the evolution.That is why Kubo's theory meets with some difficulties in the case of static perturbations [36].The formalism based on NSOM seems to be more satisfactory, because the use of the relevant statistical operator ̺ rel (t) to formulate the boundary condition for the Liouville equation can be treated as violating a complete isolation of the system.Although the static (thermodynamic) susceptibilities can be obtained from Kubo's formulae if the limits ω → 0 and ε → +0 are taken in the correct order [4], these formulae are not always convenient in concrete problems since they require a partial summation of the perturbation series.On the other hand, with a proper set of relevant variables, NSOM gives a reasonable approximation for susceptibilities and kinetic coefficients using rather crude approximations for correlation functions appearing in the response equations.
A linear response to thermal (internal) perturbations is conveniently studied in the memory function formalism (it is also called the Mori-Zwanzig formalism [37,38]) based on the projection method.The basic equations of the memory function formalism follow from the response equations derived by NSOM, when external mechanical perturbations are absent [4,39].The advantage of NSOM is that it can also be used to study a linear response to mechanical perturbations and the cross effects.
The projection formalism was generalized to far-from-equilibrium situations by Zwanzig [40], Robertson [41], and Kawasaki and Gunton [42].In the context of NSOM, Zwanzig's approach corresponds to the following representation of the relevant statistical operator: ̺ rel (t) = P̺(t), where P is a time-independent projector.In the original Zwanzig formulation, P selects the diagonal part of ̺(t) with respect to eigenstates of the main (secular) part H 0 of the Hamiltonian.Note, however, that Zwanzig's projector can also be constructed for the cases where the reduced description of nonequilibrium states is not limited to using the diagonal part of the statistical operator.
Robertson's method [41] is based on the fact that the time derivative of the relevant statistical operator (2.3) can be represented in the form ∂̺ rel (t)/∂t = P R (t) ∂̺(t)/∂t, where P R (t) is some time-dependent projector.This property allows one to derive a formally closed evolution equation for ̺ rel (t) (the generalized master equation), which, in turn, leads to a closed set of equations for the observables Pm t = Pm t rel .Robertson's scheme was improved by Kawasaki and Gunton [42], who introduced a new projector P rel (t) which has the property P rel (t) ̺(t) = ̺ rel (t).This property is analogous to the property of Zwanzig's projector P.
It can be shown (see, e.g., [2]) that all the above approaches are, in fact, different formulations of NSOM, except that sometimes one uses the initial condition for the nonequilibrium statistical operator instead of the boundary condition in a distant past.
The boundary conditions for the Liouville equation can be formulated as the so-called ergodic conditions which are a far-reaching generalization of Bogoliubov's boundary condition of the weakening of initial correlations [2,43].They reflect the fact that during the evolution of a macroscopic system its initial statistical operator (say, ̺ rel (t)) must tend to the statistical operator ̺(t) which is an integral of the Liouville equation.In [15,16] the ergodic conditions were used as the basis for constructing nonequilibrium statistical operators.Zubarev and Kalashnikov [43] (see also chapter 2 in [2]) showed that ergodic conditions lead to the same expression for the nonequilibrium statistical operator as the Liouville equation with the broken time-reversal symmetry, equation (2.1).In other words, the method of ergodic conditions is equivalent to NSOM.
Finally, we would like to mention a further scheme for constructing the nonequilibrium statistical operator.It is based on the principle of the extremum of the information entropy, including memory effects.The supplementary conditions have the form for −∞ < t ′ t.These conditions include information about the evolution of the system, whereas self-consistency conditions (2.4) include information only about the state of the system at the given time t.The conditional extremum of the information entropy with memory effects leads to the nonequilibrium statistical operator in the so-called "exponential form", which coincides with the statistical operator obtained as a retarded solution of the Liouville equation for the operator σ(t) = − ln ̺(t) [1,44].The above derivation is not without appeal since it shows a connection with irreversible thermodynamics, including memory effects.
Recently an extremum-entropy approach to NSOM was used to construct a family of nonequilibrium statistical operators which differ in the representation of the Lagrange parameters [14,45].Each of the representations corresponds to a special choice of the infinitesimal source term in the Liouville equation and, therefore, the question arises as to whether all these representations are equivalent in the sense that they give the same evolution equations for the observables.This seems to be the case because the source term is designed only to select the retarded solution of the Liouville equation and after all this the term tends to zero.Nevertheless, it is not inconceivable that there are exceptions to this rule.As far as we know, the role of the form of the source term has never been investigated.To summarize, one can say that, physically, the central ideas of the existing general approaches to nonequilibrium statistical mechanics are closely related to one another and fit naturally into the scheme of NSOM.It should be emphasized, however, that some of the fundamental problems concerning all these approaches are still open.Little is known, for instance, about the behaviour of nonequilibrium averages and transport coefficients in the thermodynamic limit.Another important problem is to construct nonequilibrium ensembles representing the state of an open system in a strong contact with other nonequilibrium systems, accompanied by fluxes of matter, energy, and entropy.In the case of weakly coupled systems we have a well-developed theory based on NSOM or other similar methods, but many of the most intriguing self-organization phenomena in open systems are beyond this weak-coupling limit.
We have not discussed the connection between NSOM and the nonequilibrium Green function formalism (see, e.g., [46,47]) which is often considered an alterna-tive to the methods based on the Liouville equation for the statistical operator.This point is discussed in our paper [25].Suffice it to note here that, in fact, these approaches complement each other and their combination appears to have a considerable promise.