Non-Equilibrium ϕ4 theory for networks: towards memory formations with quantum brain dynamics

We investigate the time evolution of quantum fields in neutral scalar ϕ4 theory for open systems with the central region and the multiple reservoirs (networks) as a toy model of quantum field theory of the brain. First we investigate the Klein–Gordon (KG) equations and the Kadanoff–Baym (KB) equations in open systems in d + 1 dimensions. Next, we introduce the kinetic entropy current and provide the proof of the H-theorem for networks. Finally, we solve the KG and the KB equations numerically in spatially homogeneous systems in 1 + 1 dimensions. We find that decoherence, entropy saturation and chemical equilibration all occur during the time evolution in the networks. We also show how coherent field transfer takes place in the networks.


Introduction
What is memory and where is located in the human brain? We know that memory is rich and diverse, can be short-term and long-term but in both cases is imperfectly stable, and has a diffused nonlocal character in the brain [1][2][3]. Conventional neuroscience does not explain these properties in a satisfactory manner nor does it provide a clear mechanistic explanation of how and where memory is stored at a molecular level [4,5].
Quantum Field Theory (QFT) of the brain, also called Quantum Brain Dynamics (QBD) represents an attempt to provide a microscopic physical mechanism of memory formation and storage in the brain [6,7]. Its origin can be traced to the seminal work by Ricciardi and Umezawa in 1967 [8], which is based on the concept of spontaneous symmetry breaking (SSB) [9][10][11], equivalent to the macroscopic order of the system, triggered by external physical stimuli. In the 1970s, Quantum Brain Dynamics was further developed by Stuart et al [12,13]. According to this theory, the brain is a mixed physical system composed of classical neurons and quantum degrees of freedom named corticons and exchange bosons. Around the same time, Fröhlich proposed a theory of electric dipole moments of molecules in biological systems which are nonlinearly coupled and interact with the heat bath and phonon modes [14][15][16][17][18][19]. The electric dipole moments of biomolecules become dynamically ordered, and as a result, a coherent wave of dipole oscillations propagates through the biological system in the form of a Fröhlich condensate. Around the same time, in 1976, Davydov and Kislukha proposed a model in which a solitary wave can propagate along protein chains carrying energy quanta without dissipative losses. This form of biological energy transfer was called the Davydov soliton [20]. Both Fröhlich and Davydov theories describe nonlinear biological coherence phenomena and are derived as static and dynamic properties in the nonlinear Schrödinger equation stemming from the same quantum Hamiltonian, respectively [21]. In the 1980s, Del Giudice et al studied several properties of water dipole fields representing them using a QFT formalism adapted for biological systems [22][23][24][25]. In the 1990s, Jibu and Yasue offered a concrete representation of quantum degrees of freedom (referred to as corticons and exchange bosons) in the dynamics of the brain, namely water electric dipole fields and photon fields [6,[26][27][28][29]. In this representation, the Quantum Field Theory of the brain is a Quantum Electrodynamics (QED) description of water molecules' electric dipoles and their interactions. More specifically, this formulation adopted a superradiant phase for the coherent state of water dipoles and photons [30][31][32][33][34]. According to this theory, memory emerges as a coherent quantum state multi-energy-mode analysis in the non-equilibrium QFT. If a coherent quantum state can be shown to be robust against decoherence in multi-energy-mode analysis, we can then conclude that memory formation processes in open systems could be adequately described using this formalism. Hence, in case the information transfer between systems can be numerically simulated in QBD, this might offer a solid path toward a solution of the binding problem.
The aim of this work is to introduce a toy model of the QBD and describe its non-equilibrium multi-energymode dynamics using the Klein-Gordon equations and the Kadanoff-Baym equations [43][44][45] for open systems with the central region and the multiple reservoirs (networks) in the neutral scalar f 4 theory. We adopt quantum tunneling processes to describe the transfer of coherent fields and incoherent particles among systems engaged in information transfer. Our work is an extension to the case of open systems of the neutral scalar theory of an isolated system [46][47][48]. In this paper, we demonstrate that it is possible to extend the two-reservoir representation to the case of an open system with N res reservoirs. We describe decoherence, entropy production and chemical equilibration via numerical simulations in 1+1 dimensions. We also show how information is transferred within networks.
The present paper is organized as follows. In section 2 we derive time evolution equations for quantum fields in open systems with the central region and the multiple reservoirs. In section 3 we introduce a kinetic entropy current and show how it enters into the H-theorem. In section 4 we demonstrate numerical simulations performed. In section 5 we discuss our results. In section 6 we provide conclusions of this work.

Time evolution equations
In this section, we introduce the Klein-Gordon (KG) equations for coherent fields and the Kadanoff-Baym (KB) equations [43][44][45] for quantum fluctuations in open systems. In figure 1, we show the central region and the N res reservoirs. (This is the generalization to the multiple reservoirs from the two reservoirs described previously [49][50][51][52][53] with tunneling effects [54][55][56][57].) First, the Lagrangian density is given by, where m is the mass, λ C and λ α ʼs are coupling constants of the interaction, v α ʼs are tunneling coupling constants between the central region C and the reservoir α [49][50][51][52][53]. We shall adopt the closed time path formalism in figure 2, the Keldysh formalism [58,59], labelled by  to describe non-equilibrium phenomena, and use 2-Particle-Irreducible (2PI) effective action technique [60][61][62].   G  G  G G  G  ,  ,  2  Tr ln  2  Tr  1 Tr  á ñ º · ·(: the density matrix), d is the spatial dimension, and the Green functions G (x, y) are given in the matrix notation by, term represents all the 2PI diagrams.
Next we derive the KG equations and the KB equations by differentiating 2PI G by f and G, respectively. The KG equations are given by, where with with the self-energy G , which is given by, We neglect off-diagonal elements of the self-energy S, since they represent higher-order terms of tunneling coupling constants. Time evolution equations or solutions of each element of equation (15) are derived in a similar way to [53] by, where we have used the z-component of the Pauli matrix σ z =diag (1,−1) and the solutions g αα of time evolution equations given by, The product of Green functions or self-energy represents convolutions, for example Gg , ; , ; , 23 and,  , ; , 28 C y , ; , ; . 32 We need not derive F αβ and ρ αβ with (a b ¹ ) since they do not appear in the energy formula in appendix A. We can derive the KG equations for C f and f ā as, In Next-to-Leading-Order (NLO) approximation of the coupling expansion of λ C , it is possible to write selfenergy as, and, , ; , ; , 37 , ; , ; 1 12 , ; , ; , 38 , ; , ; . 39 And Σ F,αα , Σ ρ,αα and , S r aã are expressed by changing λ C , F CC and ρ CC by λ α , F αα and ρ αα in equations (37), (38) and (39).

Kinetic entropy current and the H-theorem
In this section, we introduce the kinetic entropy current in the first order of the gradient expansion [63][64][65], and give a proof of the H-theorem in the Next-to-Leading Order approximation of the coupling expansion and in the Leading-Order approximation in the tunneling coupling expansion. The variable of Green functions and self-energy is (X, p) with the center-of-mass coordinate X x y 2 º + and momentum p by the use of the Fourier transformation of the relative space-time x−y of (x, y) in this section. We set t 0  -¥ in figure 2.
In a similar way to [53], by use of the (C, C) component in the Kadanoff-Baym equations (15), we arrive at, The total entropy is given by, and s ¶ m m is given by,

Numerical simulations
In this section, we show the results of our numerical simulations of the KG and the KB equations in section 2 in 1+1 dimensions. We assume spatially homogeneous systems, which means that the condition of spatial homogeneity is satisfied for each system. The initial condition of the Green functions is the same as that in [53].
The initial temperature T ini is set to be the mass m. The (L, R) reservoirs are extended to 1, 2,L, N res . We prepare the coupling constants λ C =λ α =4.0, and set the tunneling coupling constants v α =0.2/N res , with N res =2, 3, 4, 5, 10 and 100. We set the memory time mt mem =35.25 and remove the information in earlier past. We use the momentum p x

Decoherence, entropy production and chemical equilibration
The initial condition for the background coherent fields is given by, At early times 0<mX 0 <10, the f ā is amplified, and at later times 10<mX 0 , the f ā is damped with subsequent oscillations and converges asymptotically to zero. The frequencies of C f are less dependent on N res at early times mX 0 <10, but the N res dependence in frequencies appears a little after the time when f a |¯| reaches its maximum value. Note that coherent fields disappear in the time evolution for any N res .
In figure 5, we show the time evolution of the total energy density (divided by N res +1) given in appendix A. The energy error is within 0.3 %. At early times 0<mX 0 <10, the values of the total energy density are oscillating around their initial values for any N res due to sudden switch-on of self-energy at the initial time. At 35.25<mX 0 , they tend to increase monotonically due to the removal of information in earlier past for any N res . The increase appears a little larger for larger N res . The values of the total energy become smaller for larger N res , since the initial energy density of the coherent fields is divided by N res +1. The time evolution of the energy density of the coherent fields is given in appendix A in figure 6. The initial energy density of the coherent fields at mX 0 =0 has the same value for any N res in this initial condition. The energy of the coherent fields is damped in the time evolution and converges to zero. We find that decoherence occurs in open systems. The difference of   the energy density of coherent fields for N res becomes gradually larger at larger mX 0 in the time evolution. The larger the N res is, the larger the damping of the energy density becomes. At around mX 0 ∼2, 7, 10, 13, and 17, the energy density of the coherent fields appears to increase a little for N res =10 and 100. For N res =2, 3, 4, and 5, the decrease of the energy density becomes smaller at mX 0 ∼4, 7 and 10.
In figure 7, we show the incoherent particle number distribution functions at mX 0 =0 and mX 0 =250 in N res =2, and mX 0 =250 in N res =5 in the central region C. The incoherent particle number distribution function n(X 0 , p) of the mode energy Ω(X 0 , p) is defined as, x y x y x y X 0 00 0 0 Here, statistical functions are labelled by CC or α α with α=1, 2, L, N res . The difference in the n(X 0 , p) in smaller Ω(X 0 , p)<2.0 for N res =5 and N res =2 seems explicit. The n(X 0 , p) in C seems to approach the Bose-Einstein distribution with the temperature T and the chemical potential μ. The temperature T is near m at around mX 0 =250. A typical difference appears in the values of the chemical potentials μ/m=0.157 for N res =5 and μ/m=0.234 for N res =2. Even at mX 0 =250, the values of the chemical potential remain nonzero. The larger N res is, the smaller the chemical potential becomes for N res =2, 3, 4, and 5. . Figure 8. Time evolution of the number density in the system C and the reservoir α.
In figure 8, we show the time evolution of the number density N X V n X p , dp 0 2 0 x ò = p ( ) ( )with the volume V. The number density in C increases at early times 0<mX 0 <10 due to field-particle conversion. The increase becomes larger for larger N res . At mX 0 ∼2, 4, and 7, the N(X 0 )/V in C decreases a little temporarily. At later times 10<mX 0 , the number density in C starts to decrease and converges to the constant values for N res =2, and 5. The number density in the reservoir α tends to increase in the time evolution and converges to the same constant values as those in the C. The constant values become smaller for larger N res for N res =2, and 5. We find that the time scales of the convergence to the constant values become larger for larger N res . For N res =10 and 100, the difference in N/V between C and the reservoir α appears even at mX 0 =250, but N/V in C tends to approach that in the reservoir α.
In figure 9, we show the time evolution of the entropy density in the system C and the reservoir α. We adopt the entropy density with the quasi-particle approximation given by, s dp n n n n 2 1 ln 1 ln . 52 We find the same tendency as that for the number density. The entropy density in the system C increases at early times 0<mX 0 <10 for N res =2 due to field-particle conversion, but decreases at later times 10<mX 0 due to incoherent particles being transferred from C to the reservoirs. The decrease of the entropy density in C and the increase of the entropy density in the reservoir α become slower for larger N res , since the effects of incoherent particle transfer become smaller. In figure 10, the time evolution of the total entropy density is depicted. We find that total entropy density tends to increase monotonically over the time evolution. At early times mX 0 <20 with the field-particle conversion, the increase is rapid compared with that at 20<mX 0 . The nonmonotonic behavior is due to higher  order corrections in the gradient expansion. At around mX 0 =2, 4, and 7, the total entropy density tends to decrease a little. The nonmonotonic behavior disappears at later times for mX 0 >40. The increase at the later times becomes larger for N res =100 than that for N res =2. The entropy density for N res =100 is smaller than that for N res =2 since the produced incoherent particles are spreading out in C and the N res reservoirs. We find that the values of the total entropy density divided by N res +1 for N res =2 and 5 at mX 0 =250 in figure 10 are equal to s 0 in C and α for N res =2 and 5 at mX 0 =250 in figure 9.

Information transfer in networks
In this section, we show how quantum coherence is transferred within networks of interacting domains. The initial condition of the coherent fields is given by,   In figure 11, the time evolution of the coherent field 1 f is depicted. The 1 f is damped with subsequent oscillations and converges to zero in the time evolution. When N res is larger, the damping of 1 f becomes slightly smaller due to smaller coherent field transfer from the first reservoir to the system C. The frequency is less dependent on N res . The time evolution in figure 11 seems similar to that in figure 3.
In figure 12, we show the time evolution of the coherent field C f . At early times 0<mX 0 <10, the amplitude increases, but it decreases at later times 10<mX 0 . The values of C f in the local minimum at around mX 0 =10.0 are −0.105, −0.072, −0.055, −0.044, −0.022, and −0.002 for N res =2, 3, 4, 5, 10, and 100, respectively. We find that the phase difference between 1 f in figure 11 and C f in figure 12 is π/2. The frequency is less dependent on N res .
In figure 13, we show the time evolution of the coherent field f b ( 1 b ¹ ). At early times 0<mX 0 <17.5,

Discussion
In this paper, we have derived the time evolution equations in open systems with the central region and the multiple reservoirs (forming networks), namely the Klein-Gordon (KG) equations for coherent fields and the Kadanoff-Baym (KB) equations for incoherent particles. (We have generalized our results in [53] to systems with multiple reservoirs.) We have investigated the relativistic f 4 theory to describe a similar dynamics situation to that with photons, since the time evolution equation of coherent photons is the relativistic KG equation. We have introduced the kinetic entropy current and shown the H-theorem for open systems with a central region and multiple reservoirs. We have found that decoherence, entropy saturation, and chemical equilibration occur in numerical simulations of the KG and KB equations in 1+1 dimensions. The decoherence occurs due to parametric resonance instability and field-induced processes in self-energy as elaborate on in section 2. Coherence is not robust in the f 4 theory with m 2 >0. The entropy production seems to be consistent with the proof of the H-theorem shown in section 3. The chemical equilibration occurs due to the particle numberchanging processes in KB equations and the tunneling processes between systems.
In section 4.1, we have prepared a nonzero value in C f at an initial condition. In the time evolution in figure 3, the coherent field C f is damped with subsequent oscillations. The damping of the coherent field C f seems to be slightly dependent on N res . We have shown the time evolution of f ā in figure 4. We find that the maximum amplitude is proportional to the tunneling coupling constant 0.2/N res . The frequency is less dependent on the tunneling coupling constant. In a similar way to [53], these results are explained by the forced oscillation in the Klein-Gordon equation for f ā . We find that the energy of coherent fields is damped over the time evolution (decoherence) in figure 6, although there are back reactions from particles to fields. Coherence is lost and not robust in neutral scalar 4 f theory. The damping of the energy of coherent fields is dependent on the number of the reservoirs N res . This is explained by incoherent particles being transferred between systems. As is shown in figure 8, the incoherent particle number density in the central region C increases at early times 0<mX 0 <10 due to field-particle conversion (decoherence). The field-particle conversion occurs mainly in C.
The maximum value of the number density in C becomes larger as the number of the reservoirs N res increases since the tunneling (with the tunneling coupling constant v α =0.2/N res ) of incoherent particles becomes smaller for larger N res . The larger the number density is, the larger the effects of the field-induced processes become. As has been shown by several authors [46][47][48], larger damping of coherent fields occurs at higher temperatures, which means that the number density in the system is larger. Since the number density in C becomes larger for larger N res , the damping of energy of coherent fields becomes larger. In figure 8, we also see that the number density in C and in the reservoirs converges to a constant value over the course of the time evolution for N res =2 and 5. This result means that chemical equilibration occurs. In this figure, the number density near the equilibrium state is smaller for larger N res , since the produced incoherent particles are spreading within networks and the number density in each system becomes smaller for larger N res . The time scales of convergence become larger for larger N res , because tunneling of incoherent particles (with the tunneling coupling constant v α =0.2/N res ) becomes smaller for larger N res . For N res =10 and 100, the chemical equilibrium state has not been realized even at mX 0 =250. In figure 9, we show the entropy density for each system. In this figure, we find similar results to the number density in figure 8. The entropy density in C increases at early times 0<mX 0 <10 due to field-particle conversion (decoherence), but decreases at later 10<mX 0 due to incoherent particles being transferred from C to the reservoirs. The entropy density in the reservoirs increases over the course of the time evolution. The entropy density in C and that in the reservoirs converge to the constant values for N res =2 and 5. Total entropy density in figure 10 tends to increase monotonically for N res =2, 5, 10, and 100. It seems that the monotonically-increasing behavior is common for any number of reservoirs N res . This result is consistent with the H-theorem stated in section 3. Nonmonotonic behavior is due to higher order terms in the gradient expansion. In the course of time evolution, the error in the total energy density is within 0.3 %. The error has the maximum values at around the initial time due to the sudden switchon of self-energy. The error of the total energy density at later times 35.25<mX 0 <250 with memory time mt mem =35.25 does not exceed the error at around the initial time. (For mX 0 >35.25, we have removed the information in earlier past.) This is because we take sufficiently large memory time mt mem =35.25 in time evolution of the KG and KB equations. The statistical function is approximately expressed as F x z p , , with the particle number distribution n(p), the damping factor γ p and the frequency Ω p . We can estimate the order as m 10 p 0 g= at mx 0 =35 in C for N res =2 and 100 in numerical simulations in section 4.1. Hence, we can remove the information in earlier past since m/γ p ∼10 is sufficiently smaller than mt mem =35.25. The quantity γ p is related to the spectral width or Σ ρ . The larger the number density is, the larger the spectral width and the γ p become. The reason why the error of the total energy density divided by N res +1 for N res =100 is a little larger than that for N res =2 at later times 35.25<mX 0 <250 might be explained by the difference of the number density in the reservoir α. Since m/γ p becomes larger (smaller damping) by the smaller number density in the reservoir α for N res =100 than that for N res =2, the information of Green functions in the reservoir α in earlier past, which is removed at mx mt 0 mem > , might be relatively significant for N res =100 compared with N res =2.
In coherent field transfer within networks in section 4.2, the damping of coherent field 1 f is smaller for larger N res since the tunneling coupling constant v α is proportional to 1/N res . The time evolution of 1 f in section 4.2 is similar to that in C f in section 4.1 because the tunneling coupling is small. The difference between these two gradually appears in the time evolution, since C f is connected to N res reservoirs, but 1 f is connected only to C.
The coherent field in one reservoir is transferred to C. The maximum amplitude of C f (the minimum value at mX 0 =10) seems to be proportional to the tunneling coupling constant v α =0.2/N res . The phase difference between 1 f and C f is π/2. These results are explained by the analysis of the forced oscillation in the KG equation in C in a similar way in [53]. We also find that the minimum value of f b ( 1 b ¹ ) at mX 0 =17.5in other reservoirs seems to be proportional to N 1 res 2 , which is explained by forced oscillations in the KG equation in the reservoir β. The frequency in f b is less dependent on the number of the reservoirs N res . This result corresponds to information transfer, in which the change of one system by external stimuli is transferred to the hub system C and spreading to other systems. Since the phases in f b are the same, this situation might represent a form of synchronization induced by external stimuli.
We set the positive mass squared m 2 >0 in this paper, although the negative mass squared m 2 <0 term is commonly adopted in discussing the Higgs mechanism. Even in m 2 >0, we can show the Higgs mechanism with the nonzero expectation value of charged bosons j by preparing nonzero charge density of fermions as in QED with charged fermions [66]. It is also possible to show the nonzero expectation value of charged bosons j by preparing nonzero charge density of incoherent charged bosons [67]. This is because the potential energy j F[¯] is given by m e A phase ) is equivalent to the chemical potential. The Higgs mechanism in m 2 >0 is due to a nonzero chemical potential. Then, the expectation values j becomes diverse due to the diverse charge density or chemical potential. This situation might be similar to the case in which the diquark condensate qq á ñ has different values due to a different chemical potential in the SU(2) lattice gauge theory [68,69]. We might be able to prepare diverse coherent states by preparing different charge density even in QBD. We set the same mass m in C and the reservoirs α to discuss the information transfer between systems with the quantum tunneling phenomena suggested in [29]. In section 3, we find that the momentum and the frequency of incoherent particles do not change in tunneling between C and α to the first order in the gradient expansion, since no convolution for the momentum and the frequency appears in g X p G X p , , ) in equation (48). In the quasiparticle approximation in which the spectral width is narrow enough in the small number density, the mass in C must be the same as that in α to describe the tunneling phenomena. Hence, we prepared the same mass in C and α. (We also set the same coupling constant λ C =λ α to avoid different thermal effects between systems.) In QBD, the excitation of the finite number of evanescent photons represents the effect of memory retrieval. By using quantum tunneling, we might be able to describe the associative memory in which incoherent particles are transferred between systems. Then, the difference between the masses of evanescent photons might classify each memory subdivided in each area in the brain, namely by external stimuli we might maintain only memory in networks with coherent states with similar values of the mass of the evanescent photons.
We can set the Planck constant to zero in the classical approximation. The decoherence occurs even in that case, in which the incoherent particle number distribution function in the equilibrium state becomes with the temperature T instead of the Bose-Einstein distribution, since the self-energy in open systems is the same as that in the isolated system. By comparing the quantum case with the classical one, we find that the difference between the frequencies of the coherent fields gradually appears in the time evolution due to the difference between the corresponding local self-energies given by momentum integration of the statistical function proportional to the Bose-Einstein distribution or T X p , near the equilibrium states. We have started the universal usage of the same spatial variable x in the introduction of the Lagrangian. By assuming spatial homogeneity, we have constructed a model which represents various areas connected with networks having a hub. Due to spatial homogeneity, the spatial difference between the spatial boundaries and others disappears. However, since we start with QFT, it is possible to describe a phase transition occurring in this model. We can describe macroscopic properties of matter with time-dependent order parameters based on QFT. It is also possible to describe quantum tunneling phenomena such as the Josephson effect. They are the key advantages of this model.
In this paper, we have studied equilibration processes in quantum networks possessing a central region and multiple reservoirs. The expectation value of the scalar coherent fields corresponds to the square root of the number of aligned dipoles and the mass of the evanescent photons in QED with water electric dipoles (QBD). When the mass of the evanescent photons remains in the equilibrium state in open systems (coherence is maintained), we can describe memory formation processes using QBD and trace the information transfer between systems in QBD within networks. In the future, we plan to describe equilibration processes in QBD in networks. This work in the f 4 theory applied to quantum networks can be extended to the case of QBD.

Conclusion
We can trace the time evolution of the quantum fields in the f 4 theory for open systems with the central region connected to multiple reservoirs using the Klein-Gordon and Kadanoff-Baym equations. Decoherence, entropy saturation and chemical equilibration occur over the course of the time evolution. The time scales of decoherence are smaller in the case with a larger number of the connected reservoirs N res , but the time scales of chemical equilibration are larger in the case of larger N res . The information transfer has been shown to occur between the interacting systems. In particular, the phase and frequency are less dependent on the number of the interacting reservoirs N res . and, The above results are derived by extending the result in systems with the two reservoirs to that in systems with the multiple reservoirs.

Appendix B. Numerical method
We show how to calculate time evolution of the Kadanoff-Baym equations in open systems in section 2, especially time integration using the Fortran programming language. We discretize the time variable as x na t 0 = and y ma t 0 = with integers n and m and time stepsize a t . Before the time evolution is implemented, the Green functions and self-energy at n, m=0, 1 are given by their initial conditions. First, we write the do-loop of n starting with 1 for time steps of evolution of x 0 . Inside the do-loop, we calculate n 1 )with the Klein-Gordon equations. Next, we calculate time evolution of the Kadanoff-Baym equations. We write the first do-loop of m (0mn+1) and momentum p, inside which we calculate in equation (23) ) if m<n+1 and the fixing ρ CC (n+1, n+1)=0, and calculate g F,αα in equation (25) and g ρ , α α in equation (26) with the fixing g ρ,αα (n+1, n+1)=0 in a similar way to the above equations or [70]. The do-loop of p is inside the do-loop of m. In case m=0, the sum of l on the right-hand side in the above equation represents n l F l m , , Next, we can calculate F α C (n+1, m) in equation (27), ρ α C (n+1, m) in equation (28), F C α (n+1, m) in equation (29) and ρ C α (n+1, m) in equation (30)   ) and ρ C α (n+1, n+1)=0 inside the same do-loop. Then end the second do-loop of m and p.
Next, we write the third do-loop of m (0mn+1) and momentum p, inside which we calculate F αα in equation (31) and ρ αα in equation (32) with g F, α α , g ρ,αα , F Cα and ρ C α in a similar way to the above calculations. Then end the third do-loop of m and p.
Finally, we calculate the energy density, the entropy density and the number density at x 0 , and calculate the local self-energy Σ loc,C (n+1) and Σ loc, α (n+1), and the nonlocal self-energy of Σ F,CC (n+1, m), n m 1,