On the statistical description of classical open systems with integer variables by the Lindblad equation
aa r X i v : . [ qu a n t - ph ] J u l On the statistical description of classical open systems with integer variables bythe Lindblad equation
E. D. Vol ∗ B. Verkin Institute for Low Temperature Physics and Engineering of the NationalAcademy of Sciences of Ukraine 47, Lenin Ave., Kharkov 61103, Ukraine. (Dated: July 2, 2018)We propose the consistent statistical approach to consider a wide class of classical open systemswhose states are specified by a set of positive integers(occupation numbers).Such systems are oftenencountered in physics, chemistry, ecology, economics and other sciences.Our statistical methodbased on ideas of quantum theory of open systems takes into account both discreteness of thesystem variables and their time fluctuations - two effects which are ignored in usual mean fielddynamical approach.The method let one to calculate the distribution function and (or)all momentsof the system of interest at any instant.As descriptive examples illustrating the effectiveness of themethod we consider some simple models:one relating to nonlinear mechanics,and others taken frompopulation biology .In all this examples the results obtained by the method for large occupationnumbers coincide with results of purely dynamical approach but for small numbers interestingdifferences and new effects arise.The possible observable effects connected with discreteness andfluctuations in such systems are discussed.
PACS numbers: 03.65.Ta, 05.40.-a
I. INTRODUCTION
Among the vast number of classical open dynamicalsystems under consideration in physics, chemistry, biol-ogy, economics and other sciences there are many suchwhose states in accordance with the sense of the problemare specified by a set of integer variables ( { n i } ) , where i = 0,1,2,..N (N-number of degrees of freedom). For exam-ple in physics n i - are occupation numbers of cell statesin phase space, in chemistry - numbers of molecules ofreactive elements, in ecology -numbers of individuals inpopulations which live in the area and interact with otherpopulations, in economics the number of companies oper-ating on the market. Since all these systems are classicaltheir dynamics as a rule is described by a system of dif-ferential equations of the form dn i dt = F i ( { n α } ) , (1)where F i ( { n α } ) -some nonlinear functions, depending onconcrete problem. Obviously notation Eq. (1) impliesthat variables n i in system Eq. (1) -are considered ascontinuous. In the case when all n i ≫ n i are occupation numbers averagedover some appropriate statistical ensemble. In the casewhen all n i &
1, dynamical description becomes inade-quate and the question naturally arises: is there consis-tent statistical approach which takes into account bothdiscreteness of variables n i and their time fluctuationsthat may be not small. In addition it is naturally to de- ∗ Electronic address: [email protected] mand that such approach gave the same results as dy-namical description in the large n i limit. In this paperwe propose such approach based on the ideas of quan-tum theory of open systems (QTOS) and consider someexamples that demonstrate its effectiveness. The rest ofthe article organized as follows. In the Sect.2 we brieflydescribe minimal information from (QTOS) which is nec-essary for understanding of the method used and presentmain steps of our method. In Sect.3 we consider simplemodel of nonlinear autonomous oscillator with soft excit-ing mode and give its statistical description on the basisof the method proposed. All the main features of themethod clearly come to light already in this representa-tive example. In Sect.4 with the help of our approachwe consider some problems from population dynamicsrelating to evolution of two interaction populations liv-ing in a certain area..We show that in the case of smallpopulations statistical description leads to a number ofdifferences from the ordinary dynamical picture. On theother hand in the case of large occupation numbers bothdescriptions are virtually identical. In conclusion we dis-cuss some generalizations of the method and its possibleexperimental verification. II. DESCRIPTION OF THE METHOD
In this section we briefly remind the main points of themethod proposed by author earlier [1] which allows oneto make the transition from known dynamical equationsof classical open system to the master equation for itsquantum analogue. The method based on the correspon-dence that can be set between quantum master equationin the Lindblad form and the Liouville equation for distri-bution function in phase space of the classical system ofinterest. This correspondence allows one using classicalequations of motion to restore the form of all operatorsinvolved in the Lindblad equation. Thereby we can ap-ply the procedure of quantization at least in semiclassicalapproximation in the case of a large class of nonhamil-tonian dynamical systems. In short (all details see in[1] ) the recipe of quantization proposed consists of threeconsecutive steps.
Step1:
The input classical dynamical equationsshould be presented in the form allowed the quantiza-tion (FAQ). For the purposes of present paper the mostconvenient form is complex representation of equationsof motion: dz i dt = − i · dHdz ∗ i + X α (cid:18) R α dR α dz ∗ i − R α dR α dz ∗ i (cid:19) , (2)where z i = x i + iy i √ , z ∗ i = x i − iy i √ are complex dynami-cal coordinates of the system of interest, H , R α , R α arefunctions of z i , z ∗ i ( H is real function, and R α , R α arecomplex, R α means function which conjugate to R α ). Itis necessary to emphasize that it is the most delicate stepof the method because it is difficult exactly to formalizethis point. Step2.
Having in hands representation Eq. (2) wecan use classical function H , R α , R α and with theirhelp determine their quantum analogues - operators b H , c R α , c R + . For this purpose the variables z i , z ∗ i mustbe replaced by the correspondence Bose operators b a i and b a + i with usual commutation rules: (cid:2)b a i , b a + j (cid:3) = δ ij ,[ b a i , b a j ] = (cid:2)b a + i , b a + j (cid:3) = 0. Step3.
The operators b H , c R α , c R + α found in this man-ner should be substituted into the quantum Lindbladequation for the evolution of density matrix of the sys-tem: d b ρdt = − i h b H, b ρ i + X α h c R α · b ρ, c R + α i + h c R α , b ρ · c R + α i . (3)The correspondence principle guarantees us that such ap-proach will give correct description of evolution of quan-tum open system at least with the accuracy up to thefirst order in h. III. REPRESENTATIVE MODEL.AUTONOMOUS NONLINEAR OSCILLATORWITH SELF EXCITING MODE
It is convenient to demonstrate on concrete exampleall the features of the approach proposed. We consideran oscillator with nonlinear damping in situation whenits equilibrium point loses stability and small fluctua-tions switch the system to the new stationary state thatcorresponds to the closed trajectory (limit cycle). Thefollowing system of equations gives the correct mathe-matical description of the behavior of the oscillator inthe vicinity of bifurcation point (see [2]): dxdt = ωy + µx − x (cid:0) x + y (cid:1) , (4) dydt = − ωx + µy − y (cid:0) x + y (cid:1) , where x , y are coordinates of oscillator in phase spaceand ω -its frequency. The system of the equations ofmotion Eq. (4) can be written in the complex form asone equation dzdt = − iωz + µz − z | z | , (5)where z = x + iy √ . It is easy to show that equation Eq.(5) can be represented in the FAC. For this purpose weintroduce the functions H = ωz ∗ z , R = √ µz ∗ and R = z . One can verify by direct checking that r.h.s. of Eq.(5) can be rewritten in the form dzdt = − i ∂H∂z ∗ + (cid:18) R ∂R ∂z ∗ − R ∂R ∂z ∗ (cid:19) + (cid:18) R ∂R ∂z ∗ − R ∂R ∂z ∗ (cid:19) . (6)According to the recipe of quantization the Lindbladequation for quantum analog of the system Eq. (4) canbe written in the form dρdt = − i [ H, ρ ] + (cid:2) R ρ, R +1 (cid:3) + (cid:2) R ρ, R +2 (cid:3) + h.c, (7)where b H = ω b a + b a , b R = √ µ b a + and b R = b a .We are interesting in the stationary solution of Eq. 7and from physical considerations we will seek it in theform: b ρ = b ρ ( b N ) = P | n i ρ n h n | , where | n i are the eigen-functions of the operator b N = b a + b a . It is convenient to in-troduce the generating function G ( u ) for the coefficients ρ n . By definition G ( u ) = n = ∞ P n =0 ρ n · u n . One can obtainfrom Eq. 7 in stationary case the following equation forthe function G ( u ):(1 + u ) d Gdu − µu dGdu − µG ( u ) = 0 . (8)The solution of Eq. 8 that satisfies all conditions of theproblem may be represented as G ( u ) = Φ(1 , µ, µ (1 + u ))Φ(1 , µ, µ ) , (9)where Φ( a, c, x ) is well known confluent hypergeometricfunction (see [3]) - Φ( a, c, x ) = 1 + axc + a ( a +1)2! c ( c +1) x + ... .Note that condition G ( u ) = 1 corresponds to normaliza-tion of c ρ st ( n ) − namely P n ρ n = 1. Having in hands theexpression Eq. 9 for the generation function one can findthe average values for any physical quantity of interestin the stationary state that is all moments of the distri-bution ρ ( n ). For example moments of first and secondorder are determined by relations: n = P n n · ρ n = dGdu , n − n = d Gdu . (all derivatives are taken at point u = 1).Let us consider now the behavior of our system in twolimiting cases: µ ≫ µ ≪
1. When µ ≫ a, c, x ) (see [3])we obtainfor G ( u ) the next expression G ( u ) ≈ e µ ( u − · (cid:18) u (cid:19) − µ , (10)The average number of quanta in stationary state n andtheir dispersion σ = n − n in this case are: n = µ +12 and σ = µ +14 . The relative fluctuation of quanta gen-erated in stationary state is √ σn = √ µ +1 µ +1 tends to zerowhen µ ≫
1. The obtained result testifies validity ofdeterministic approach in this case since in classical case n cl = | z | = µ when the system is moving along the limitcycle. Now let us consider the opposite case when µ ≪ a, c, x ) [3] we obtainedfor G ( u ) the next expression (up to the first order in µ ) G ( u ) ⋍ (2 + u )(1 + µ ) + µ (1 + u ) µ . (11)The relation Eq. 11 implies that G ( u ) tends to G ( u ) = u when µ tends to zero. The first two moments ofthe distribution ρ ( n ) are: n = and σ = . The relativefluctuation of n = √ σn is equal to √ n = 0 when µ tendsto zero moreover fluctuations of occupation number turnout to be large and very essential.It is necessary to note that in classical system Eq. 4in addition to variable | z | we have phase variable ϕ thatsatisfies to equation dϕdt = ω . However phase dynam-ics determines only velocity along limit cycle but doesnot influence on stationary state itself. In principle wecan turn ω to zero (which implies that H = 0) and allforegoing results do not change. It is the point whichexplains why the Lindblad equation can be successfullyapplied for statistical description of classical open sys-tems with integer variables. The reason is that in thecase when H = 0 in Eq. 4, h may be eliminated from theLindblad equation and wave properties of quantum sys-tem turn out to be irrelevant. After this crucial remarkwe can apply our approach to different classical systemswith integer variables in particular population dynamicsmodels. IV. STATISTICAL DESCRIPTION OFPOPULATION DYNAMICS MODELS.
We begin our consideration with well known Lotka-Volterra model (LVM) (see [4]) describing ”the interac-tion” between two populations: prays and predators (forexample hares and lynxes) living in the same territory. Let n ( t ) and n ( t ) the current numbers of preys andpredators correspondently. Then input dynamical equa-tions of the LVM can be written as dn dt = n ( a − n ) , (12) dn dt = − n ( b − n ) , where coefficients a , b are positive and have clear ecolog-ical meaning. Note that for convenience we have chosentime unit thus to put coefficient of the term n n in Eq.12 equal to 2. Now let us demonstrate that system Eq.12 can be represented in FAQ. For this purpose we intro-duce two auxiliary complex variables z and z so that n = | z | and n = | z | . Let us assume that evolu-tion z , z in time is governed by the following system ofequations: dz dt = λ z − z | z | , (13) dz dt = − λ z + z | z | . It is easy to see that system Eq. 13 implies the followingequations for n , n dn dt = 2 λ n − n n , (14) dn dt = − λ n + 2 n n , which is exactly coincides with Eq. 12 if one lets a = 2 λ , b = 2 λ .Now if we introduce three functions R = λ z ∗ , R = λ z , and R = z z ∗ it is easy to verify that Eq. 13 canbe represented in the form dz dt = X i =1 (cid:18) R i ∂R i ∂z ∗ − R i ∂R i ∂z ∗ (cid:19) , (15) dz dt = X i =1 (cid:18) R i ∂R i ∂z ∗ − R i ∂R i ∂z ∗ (cid:19) . Using the foregoing recipe of quantization we maintainthat description of LVM which takes into account thediscreteness of variables n and n and their time fluctu-ations is given by the next master equation dρdt = X i =1 h b R i ρ, b R + i i + h.c., (16)where b R = λ b a +1 , b R = λ b a and b R = b a b a +2 . Againwe will interesting in only the solutions of Eq. 16 whichhave the form b ρ = P n ,n | n n i ρ n , n h n n | . In this casefor the coefficients of the expansion ρ n n we obtain thenext general equation dρ n n dt = 2 λ (cid:2) n ρ n − ,n − ( n + 1) ρ n n (cid:3) +2 λ h ( n + 1) ρ n , n +1 − n ρ n n i +2 (cid:2) ( n + 1) n ρ n +1 ,n − − ( n + 1) n ρ n n (cid:3) (17). It is convenient to introduce the generating function G ( u, v, t ) = P n n ρ n n · u n · v n . Then after the simplealgebra we find that Eq. 17 implies the next equation for G ( u, v, t ): ∂G∂t = 2 λ ( u − ∂∂u ( uG )+2 λ (1 − v ) ∂G∂v +2( v − u ) ∂ ∂u∂v ( vG ) . (18)The Eq. 17 and Eq. 18 in principle give us all the nec-essary information about statistical behaviour of LVM.Now let us consider concrete results following from theseequations. Note that even in stationary case it is difficultto find exact analytical solutions of Eq. 18 for arbitrary λ and λ . But in the special case when λ = 1 + λ suchsolution easily can be found and have the form G ( u, v ) = (1 − κ ) (1 − κ u )(1 − κ v ) where κ = λ λ . Since G ( u, v ) = g ( u ) · g ( v )it is clear that n and n are independent variables inthis case and n = n = κ − κ = λ . The dispersion σ inthis case is equal to λ + λ and if we calculate relativefluctuation of n in stationary state we get the result: δn = √ σ n = q λ . Thereby we see this quantityis not small and can be easily measured. More detailedanalysis of statistical behaviour of LVM following fromEq. 17 for arbitrary λ and λ will be carried out else-where. Here we show only that using Eq. 18 one caneasy obtain a collection of explicit relations connectingdifferent moments of distribution ρ ( n , n ). For exampleif we differentiate stationary Eq. 18 with respect to u andafter that put u = v = 1 we obtain the simple relationbetween moments of first and second order which readsas: λ (1 + n ) − n − n n = 0 . (19)In a similar manner by differentiating of Eq. 18 withrespect to v we obtain the second relation: − λ n + n + n n = 0 (20)Relations Eq. 19 and Eq. 20 imply the helpful equationconnecting n and n namely n n = λ λ . It is worthto note that in classical LVM Eq. 14 similar relationexists: n n = λ λ , so we conclude that when n ≫ n ≫ n , n the difference between themmay be essential. To demonstrate this distinction andalso to compare the approach proposed in the presentpaper with usual Markovian description of such systemsproposed in well known article of Nicolis and Prigogine [5] and expanded in their later book [6] it is appropriateto consider the truncated case of LVM, with λ = λ = 0.It is obvious that total number of individuals N= n + n in this model will be constant .This fact greatly simplifiesfinding and analysis of solutions Eq. 18 which in this casetakes the form: ∂G∂t = ( v − u ) ∂ ∂u∂v ( vG ) (21)It is easy to see that Eq. 21 has solutions for any inte-ger N in the form of homogeneous polynomial in u and v of degree N , namely G N ( u, v, t ) = P k A k ( t ) u k v N − k where coefficients A k satisfy the normalization condition P k A k = 1. Thereby Eq. 21 is reduced to the linearsystem of equations of N+1order for coefficients A k ofthe form dA k dt = L km A m , where matrix elements L km can directly be found from Eq. 21 for any N. For exam-ple in simplest case when N=2 the matrix L km has theform − − .Let us consider now this case more detail. Let G ( u, v, t ) = au + buv + cv is generating function ofthe model. Then Eq. 21 implies the next system ofequations for evolution of coefficients a, b, c . dadt = − a, dbdt = 2 a − b, dcdt = 2 b. (22)Together with normalization condition a + b + c = 1 sys-tem Eq. 22 allows one to give statistical description ofthe model at any time. In particular Eq. 22 implies thatwhen t tends to infinity n tends to zero and n tendsto 2.This completely agrees with solutions of dynamicalequations in this case. From the other hand let us con-sider now the equation for generating function of thismodel obtained by Nicolis and Prigogine (see eq .10.66in their book [6]) which in our notation has the form dGdt = ( v − u ) v ∂ G∂u∂v . (23)In the case N = 2 Eq. 23 implies the system of equationsfor coefficients of expansion G = au + buv + cv whichdiffers from 22 namely dadt = 0 , dbdt = − b, dcdt = b. (24)Eq. (24) imply that when t tends to infinity both quan-tities n , n tend to nonzero values depending on initialconditions what obviously does not agree with dynami-cal equations. But when N ≫ n and n in such model can be represented of simple system ofequation of the form dn dt = an n − bn n , (25) dn dt = − an n + bn n . It is easy to see that total number of individuals in thismodel N= n + n conserves. But dynamical descrip-tion of the system with the help of Eq. (25) seems tobe oversimplified. In particular it implies that for any Nwhen a > b and t tends to infinity n tends to N , and n tends to zero. Now following the spirit of our methodwe will describe this system by the help of two operators c R = λ b a c a +2 and c R = λ c a +1 b a . Show that such statis-tical version is completely consistent. Actually acting asin previous examples we can write the master equationfor the density matrix of the system as d b ρdt = h c R b ρ, c R +1 i + h c R b ρ, c R +2 i + h.c.. (26)If again we are interesting in by solutions of the Eq.(26) of the form b ρ = P n n | n n i ρ n n n h n n | then forthe generating function of the problem G ( u, v, t ) = P n n ρ n n u n v n we obtain the equation dGdt = a ( u − v ) d dudv ( uG ) + b ( v − u ) d dudv ( vG ) , (27)where a = 2 λ and b = 2 λ . The equation Eq. (27)implies that generating function of stationary state ofthis model in the case when total number of individualsis equal to N can be represented in the form G st ( u, v ) = C N " ( bv ) N +1 − ( au ) N +1 bv − au = C N h ( bv ) N + ( bv ) N − ( au ) + ... i (28)where C N is the normalization factor, C N = (cid:0) b N + ab N − + ... (cid:1) − Having in hands expression Eq.(28) for the generation function one can find all statisti-cal characteristics of the model for any N. In particularfor N=2 for the average values of individuals in both kinswe obtain n = a + abb + ab + a , and n = b + abb + ab + a . Now we want to show that when N tends to infinity the proper-ties of our statistical model will be similar to ones of thedynamical model Eq. (25).Let us assume that b > a andlet κ = ab , ( κ < . Let us calculate now the ratio n n in this model. Using Eq. (28) we obtain n ←− n = N a N + ( N − a N − b + ...N b N + ( N − b N − a + ... = N κ N + ( N − κ N − + ... κ N + ( N − κ + .... κ N − (29). Expression Eq. (29) can be represented in the next convenient form n n = κ + κ ∂Lnf N ( κ ) ∂ κ N − κ ∂Lnf N ( κ ) ∂ κ , (30)where f N ( κ ) = − κ N − κ . It is easy to see that ∂ ( Lnf N ) ∂ κ = ( N − κ N − N κ N − +1(1 − κ )(1 − κ N ) . Obviously this expression tends to − κ when κ < n n tends to zero but this case is realizedonly for infinite population. From the other hand for anyfinite population of cannibals all its statistical propertiescan be well described by the approach proposed in our pa-per. In conclusion let us briefly summing up our consid-eration. We propose the statistical method of describingof various systems in physics,chemistry,ecology whosestates can be represented by integers. Although themethod proposed is based on quantum Lindblad equa-tion nevertheless we have shown that it can be success-fully used also for description of classical open systems with integer variables (at least in the cases when ”Hamil-tonian” of the open system is equal to zero) .All examplesconsidered in present paper demonstrate that approachproposed results in consistent conclusions and allows onein several cases to eliminate essential inaccuracies of pre-ceding considerations. V. ACKNOWLEDGEMENT
The author acknowledges I.V. Krive and S.N. Dolya forthe discussions of the results of the paper and valuablecomments. [1] E. D. Vol Phys. Rev. A73