**1. Introduction**

Numerous attempts to understand memory in a brain have been made over one hundred years starting at the end of 19th century. Nevertheless, the concrete mechanism of memory still remains an open question in conventional neuroscience [1–3]. Conventional neuroscience is based on classical mechanics with neurons connected by synapses. However, we still cannot answer how limited connections between neurons describe mass excitations in a brain in classical neuron doctrine.

Quantum field theory (QFT) of the brain or quantum brain dynamics (QBD), is one of the hypotheses expected to describe the mechanism of memory in the brain [4–6]. Experimentally, several properties of memory, namely the diversity, the long-term but imperfect stability and nonlocality (Memory is diffused and non-localized in several domains in a brain. It does not disappear due to the destruction in a particular local domain. The term 'nonlocality' does not indicate nonlocality in entanglement in quantum mechanics.), are suggested in [7–9]. The QBD can describe these properties by adopting infinitely physically or unitarily inequivalent vacua in QFT, distinguished from quantum mechanics which cannot describe unitarily inequivalence. Unitarily inequivalence represents the emergence of the diversity of phases and allows the possibility of spontaneous symmetry breaking (SSB) [10–13]. The vacua or the ground states appearing in SSB describe the stability of the states. Furthermore, the QFT can describe both microscopic degrees of freedom and macroscopic matter [10]. To describe stored information, we can adopt the macroscopic ordered states in QFT with SSB involving long-range correlation via Nambu–Goldstone (NG) quanta. In 1967, Ricciardi and Umezawa proposed a quantum field theoretical approach to describe memory in a brain [14]. They adopted the SSB with long-range correlations mediated by NG quanta in QFT. Stuart et al. developed QBD by assuming a brain as a mixed system of classical neurons and quantum degrees of freedom, namely corticons and exchange bosons [15,16]. The vacua appearing in SSB, the macroscopic order, are interpreted as the memory storage in QBD. The finite number of excitations of NG modes represents the memory retrieval. Around the same time, Fröhlich proposed the application of a theory of electric dipoles to the study of biological systems [17–22]. He suggested a theory of the emergence of a giant dipole in open systems with breakdown of rotational symmetry of dipoles where dipoles are aligned in the same direction (the ordered states with coherent wave propagation of dipole oscillation in the Fröhrich condensate). In 1976, Davydov and Kislukha studied a theory of solitary wave propagation in protein chains, called the Davydov soliton [23]. It is found that the theory by Fröhlich and that by Davydov represent static and dynamical properties in the nonlinear Schödinger equation with an equivalent quantum Hamiltonian, respectively [24]. In the 1980s, Del Giudice et al. applied a theory of water electric dipoles to biological systems [25–28]. In particular, the derivation of laser-like behavior is a suggestive study. In the 1990s, Jibu and Yasue gave a concrete picture of corticons and exchange bosons, namely water electric dipole fields and photon fields [4,29–32]. The QBD is nothing but quantum electrodynamics (QED) with water electric dipole fields. When electric dipoles are aligned in the same directions coherently, the polaritons, NG bosons in SSB of rotational symmetry, emerge. The dynamical order in the vacua in SSB is maintained by long-range correlation of the massless NG bosons. In QED, the NG bosons are absorbed by photons and then photons acquire mass due to the Higgs mechanism and can stay in coherent domains. The massive photons are called evanescent photons. The size of a coherent domain is in the order of 50 μm. Furthermore, two quantum mechanisms of information transfer and integration among coherent domains are suggested. The first one is to use the super-radiance and the self-induced transparency via microtubules connecting two coherent domains [31]. Super-radiance is the phenomenon indicating coherent photon emission with correlation among not only photons but also atoms (or dipoles) [33–37]. The atoms (or dipoles) cooperatively decay in a short time interval due to correlation; coherent photons with intensity proportional to the square of the number of atoms (or dipoles) are emitted. The pulse wave photons in super-radiance propagate through microtubules without decay. Then the self-induced transparency appears, since microtubules are perfectly transparent in the propagation. The second one is to use the quantum tunneling effect among coherent domains surrounded by incoherent domains [32]. The effect is essentially equivalent to the Josephson effect between two superconducting domains separated by a normal domain. Del Giudice et al. studied this effect in biological systems [28]. In 1995, Vitiello has shown that a huge memory capacity can be realized by regarding a brain as an open dissipative system and doubling the degrees of freedom with mathematical techniques in thermo-field dynamics [38]. In dissipative model of a brain, each memory state evolves in classical deterministic trajectory like a chaos [39]. The overlap among distinct memory states is zero at any time in the infinite volume limit. However, finite volume effects allow states to overlap one another, which might represent association of memories [6]. In 2003, exclusion zone (EZ) water was discovered experimentally [40]. The properties of EZ water correspond to those of coherent water [41].

However, we have never seen the dynamical memory formations based on QBD at the physiological temperature in the presence of thermal effects written by quantum fluctuations. Hence, there are still criticisms related with the decoherence phenomena (We should use the mass of polaritons in estimating the critical temperature of ordered states, not that of water molecules themselves.)in memory formations in QBD [42]. So, we need to derive time evolution equations of coherent fields and quantum fluctuations and show numerical simulations of memory formation processes in non-equilibrium situations to check whether or not memory in QBD is robust against thermal effects. Futhermore, in 2012 Craddock et al. suggested the mechanism of memory coding in microtubules with phosphorylation by Ca2<sup>+</sup> calmodulin kinase II [43]. It will be an interesting topic to investigate how water electric dipoles and evanescent photons are affected by phosphorylated microtubules.

The aim of this paper is to derive time evolution equations, namely the Schrödinger-like equations for coherent dipole fields, the Klein–Gordon equations for coherent photon fields, the Kadanoff–Baym equations for quantum fluctuations [44–46], with the two-particle-irreducible effective action technique with Keldysh formalism [47–51]. We derive both the equilibration for quantum fluctuations and the super-radiance for background coherent fields from the single Lagrangian in quantum electrodynamics (QED) with electric dipole fields. We arrive at the Maxwell–Bloch equations for the super-radiance by starting with QED with electric dipole fields in 2 + 1 dimensions. When we consider electric fields in super-radiance, we only need two spatial dimensions, one axis for the amplitude and another axis for the propagation. Hence we have discussed the case in 2 + 1 dimensions in this paper. By using our equations for super-radiance in this paper, we can describe information transfer via microtubules. Then, microtubule-associated proteins can make an important contribution to information transfer with interconnections among microtubules. We also derive the Higgs mechanism and the tachyonic instability for coherent fields in the Klein–Gordon equation for coherent electric fields. In two energy level approximation for electric dipole fields, namely with the ground state and the first excited states, the Higgs mechanism appears in normal population in which the probability amplitude in the ground state is larger than that in the first excited states. The penetrating length in the Meissner effect due to the Higgs mechanism is 6.3 μm derived by using coefficients in 2 + 1 dimensions and the number density of liquid water molecules in 3 + 1 dimensions. On the other hand, the tachyonic instability appears in inverted population in which the probability amplitudes in the first excited states are larger than that in the ground state. Then the electric field increases exponentially while the system is in inverted population. The increase stops at times when normal population is realized. Our analysis also contains the dynamics of quantum fluctuations in non-equilibrium cases. We also derive the Kadanoff–Baym equations for quantum fluctuations with the leading-order self-energy in the coupling expansion. The Kadanoff–Baym equations describe the entropy producing dynamics during equilibration as shown in the proof of the H-theorem. Entropy production stops when the Bose–Einstein distribution is realized. By combining time evolution equations (the Klein–Gordon equations for coherent electric fields and the Schrödinger-like equations for coherent electric dipole fields) and the Kadanoff–Baym equations for quantum fluctuations, we can describe the dynamical behavior of dipoles with thermal effects written by quantum fluctuations. Our analysis will be applied to memory formation processes in QBD. In particular, by extending our method to the case in open systems (networks), we can also trace dynamical memory recalling processes with excitations of particles in coherent domains via quantum tunneling processes, which are described by the Kadanoff–Baym equations. We can perform the simulations of the dynamical recalling processes in QBD with our equations to understand our thinking processes.

This paper is organized as follows. In Section 2, we introduce the two-particle-irreducible effective action in the closed-time path contour to describe non-equilibrium phenomena and derive time evolution equations. In Section 3, we introduce a kinetic entropy current in the first order of the gradient expansion, and show the H-theorem in the leading-order approximation of the coupling expansion. In Section 4, we show the time evolution equations, the conserved total energy and the potential energy in spatially homogeneous systems in an isolated system. In Section 5, we derive the super-radiance by analyzing the time evolution equations for coherent fields. In Section 6, we discuss our results. In Section 7, we provide the concluding remarks. In the Appendix A, we show how quantum fluctuations appear as additional terms in the Klein–Gordon equations. In this paper, the labels *i*, *j* = 1 and 2 represent *x* and *y* directions in space, the labels *a*, *b*, *c*, *d* = 1, 2 represent two contours in the closed-time path, the labels *α* = −1, 1 represent the angular momentum of electric dipoles. The speed of light, the Planck constant divided by 2 *π* and the Boltzmann constant are set to be 1 in this paper. We adopt the metric tensor *ημν* = diag(1, −1, −<sup>1</sup>) with *μ*, *ν* = 0, 1, 2.

#### **2. The Two-Particle-Irreducible Effective Action and Time Evolution Equations**

We begin with the following Lagrangian density to describe quantum electrodynamics (QED) with electric dipoles in 2 + 1 dimensions in the background field method [52–55],

$$\begin{split} \mathcal{L}[\mathbf{Y}^\*(\mathbf{x}, \boldsymbol{\theta}), \mathbf{Y}(\mathbf{x}, \boldsymbol{\theta}), A(\mathbf{x}), a(\mathbf{x})] &= -\frac{1}{4} F^{\mu\nu}[A+a] F\_{\mu\nu}[A+a] - \frac{(\partial^{\mu} a\_{\mu})^2}{2\alpha\_1} \\ &+ \int\_0^{2\pi} d\theta \left[ \mathbf{Y}^\* i \frac{\partial}{\partial x^0} \boldsymbol{\Psi} + \frac{1}{2m} \mathbf{Y}^\* \boldsymbol{\nabla}\_i^2 \mathbf{Y} \right. \\ &\left. + \frac{1}{2I} \mathbf{Y}^\* \frac{\partial^2}{\partial \theta^2} \mathbf{Y} - 2cd\_c \mathbf{Y}^\* u^i \boldsymbol{u}^j \boldsymbol{\Psi} F^{0i} [A+a] \right], \end{split} \tag{1}$$

where *A* is the background coherent photon fields, *a* is the quantum fluctuations of photon fields, *Fμν*[*A*] = *∂μAν* − *∂νAμ* is the field strength, the *α*1 is a gauge fixing parameter, the *m* is the mass of a dipole, the *I* is the moment of inertia, *ui* = (cos *θ*, sin *θ*) is the direction of dipoles and 2*ede* is the absolute value of dipole vector. The variable *θ* represents the degrees of freedom of rotation of dipoles in 2 + 1 dimensions. The dipole–photon interaction term −2*ede*Ψ<sup>∗</sup>*ui*Ψ*F*0*<sup>i</sup>*[*<sup>A</sup>* + *a*] has the similar form to that in [27]. We shall expand the electric dipole fields Ψ and Ψ∗ by the angular momentum and consider only the ground state and the first excited states in energy-levels. Then we can write them as,

$$\begin{aligned} \Psi(\mathbf{x},\boldsymbol{\theta}) &= \frac{1}{\sqrt{2\pi}} \left( \psi\_0(\mathbf{x}) + \psi\_1(\mathbf{x})e^{i\boldsymbol{\theta}} + \psi\_{-1}(\mathbf{x})e^{-i\boldsymbol{\theta}} \right), \\ \Psi^\*(\mathbf{x},\boldsymbol{\theta}) &= \frac{1}{\sqrt{2\pi}} \left( \psi\_0^\*(\mathbf{x}) + \psi\_1^\*(\mathbf{x})e^{-i\boldsymbol{\theta}} + \psi\_{-1}^\*(\mathbf{x})e^{i\boldsymbol{\theta}} \right), \end{aligned} \tag{2}$$

in 2 + 1 dimensions. (In 3 + 1 dimensions, we might expand Ψ and Ψ∗ by spherical harmonics.) We can rewrite the terms in the above Lagrangian as,

$$\int d\theta \Psi^\*(\mathbf{x}, \theta) i \frac{\partial}{\partial \mathbf{x}^0} \Psi(\mathbf{x}, \theta) \quad = \underbrace{\Psi\_0^\* i}\_{\cdot} \frac{\partial}{\partial \mathbf{x}^0} \Psi\_0 + \Psi\_1^\* i \frac{\partial}{\partial \mathbf{x}^0} \Psi\_1 + \Psi\_{-1}^\* i \frac{\partial}{\partial \mathbf{x}^0} \Psi\_{-1} \tag{3}$$

$$\int d\theta \frac{1}{2m} \Psi^\* \nabla\_i^2 \Psi^\* = \frac{1}{2m} \left[ \psi\_0^\* \nabla\_i^2 \psi\_0 + \psi\_1^\* \nabla\_i^2 \psi\_1 + \psi\_{-1}^\* \nabla\_i^2 \psi\_{-1} \right],\tag{4}$$

$$\int d\theta \frac{1}{2I} \Psi^\* \frac{\partial^2}{\partial \theta^2} \Psi^\* = -\frac{-1}{2I} \left[ \psi\_1^\* \psi\_1 + \psi\_{-1}^\* \psi\_{-1} \right]. \tag{5}$$

We also write the dipole–photon interaction term with electric fields *F*0*<sup>i</sup>* = −*Ei* by,

$$\begin{split} \int d\theta 2cd\_{\varepsilon} \Psi^\* u^i \Psi E\_i &= \left. \epsilon d\_{\varepsilon} \int d\theta \left[ (E\_1 - iE\_2) \Psi^\* \epsilon^{i\theta} \Psi + (E\_1 + iE\_2) \Psi^\* \epsilon^{-i\theta} \Psi \right] \right| \\ &= \left. \epsilon d\_{\varepsilon} \left[ (E\_1 - iE\_2) (\psi\_0^\* \psi\_{-1} + \psi\_1^\* \psi\_0) + (E\_1 + iE\_2) (\psi\_0^\* \psi\_1 + \psi\_{-1}^\* \psi\_0) \right] \right|, \end{split} \tag{6}$$

with the direction of dipoles *ui* = (cos *θ*, sin *<sup>θ</sup>*).

Next, we show two-particle-irreducible (2PI) effective action [47–49] for electric dipole fields and photon fields. Starting with the above Lagrangian density, we write the generating functional with the gauge fixing condition for quantum fluctuation,

$$\text{gauge fixing } \therefore a^0 = 0,\tag{7}$$

and perform the Legendre transformations. Then we arrive at,

$$\begin{split} \Gamma\_{2\text{Tr}}[A, d^i \vec{\psi}, \vec{\Psi}^\*] &= \int\_{\mathcal{C}} d^{d+1} \mathbf{x} \left[ -\frac{1}{4} F^{\mu \nu} [A + \vec{a}] F\_{\mu \nu} [A + \vec{a}] + i \vec{\psi}\_0^\* \frac{\partial}{\partial \mathbf{x}\_0} \vec{\psi}\_0 + \sum\_{a = -1, 1} i \vec{\psi}\_a^\* \frac{\partial}{\partial \mathbf{x}\_0} \vec{\psi}\_a \\ &+ \frac{1}{2m} \left( \vec{\psi}\_0^\* \nabla\_i^2 \vec{\Psi}\_0 + \sum\_{a = -1, 1} \vec{\psi}\_a^\* \nabla\_i^2 \vec{\Psi}\_a \right) - \frac{1}{2I} \sum\_{a = -1, 1} \vec{\psi}\_a^\* \vec{\Psi}\_a \\ &+ e d\_{\mathcal{C}} \sum\_{a = -1, 1} \left[ (E\_1 + i a E\_2) (\vec{\psi}\_0^\* \vec{\Psi}\_a + \vec{\Psi}\_{-a}^\* \vec{\Psi}\_0) \right] \Bigg] \\ &+ i \text{Tr} \ln \mathbf{\Delta}^{-1} + i \text{Tr} \Delta\_0^{-1} \Delta + \frac{i}{2} \text{Tr} \ln D^{-1} + \frac{i}{2} \text{Tr} D\_0^{-1} D + \frac{\Gamma\_2[\Delta, D]}{2}, \tag{8} \end{split} \tag{9}$$

where the C represents the Keldysh contour [50,51] shown in Figure 1, the spatial dimension *d* = 2, the bar represents the expectation value · with the density matrix. The 3 × 3 matrix *i***Δ**−<sup>1</sup> 0 (*<sup>x</sup>*, *y*) is defined as follows,

$$\begin{array}{rcl} i\Delta\_{0}^{-1}(\mathbf{x},y) & \equiv & \left. \frac{\delta^{2} \int\_{\mathbf{x}} \mathcal{L}}{\delta\psi^{\*}(y)\delta\psi(\mathbf{x})} \right|\_{\mathbf{z}=0} \\ &=& \left[ \begin{array}{cccc} i\frac{\partial}{\partial\mathbf{x}^{0}} + \frac{\nabla\_{i}^{2}}{2\mathbf{m}} - \frac{1}{2\mathsf{I}} & ed\_{\mathsf{e}}(E\_{1}+iE\_{2}) & 0 \\\\ ed\_{\mathbf{e}}(E\_{1}-iE\_{2}) & i\frac{\partial}{\partial\mathbf{x}^{0}} + \frac{\nabla\_{i}^{2}}{2\mathbf{m}} & ed\_{\mathsf{e}}(E\_{1}+iE\_{2}) \\ 0 & ed\_{\mathsf{f}}(E\_{1}-iE\_{2}) & i\frac{\partial}{\partial\mathbf{x}^{0}} + \frac{\nabla\_{i}^{2}}{2\mathbf{m}} - \frac{1}{2\mathsf{I}} \end{array} \right] \delta\_{\mathcal{C}}^{d+1}(\mathbf{x}-\mathbf{y}), \end{array} \tag{9}$$

for −1, 0 and 1, and the *iD*−<sup>1</sup> 0,*ij*(*<sup>x</sup>*, *y*) is written by,

$$\begin{array}{rcl} iD\_{0jj}^{-1}(\mathbf{x}, y) & \equiv & \frac{\delta^2 \int\_{\mathcal{X}} \mathcal{L}}{\delta a^i(\mathbf{x}) \delta a^j(y)} \\ & = & -\delta\_{lj} \partial\_{\mathbf{x}}^2 \delta\_{\mathcal{C}}^{d+1}(\mathbf{x} - y), \end{array} \tag{10}$$

where *i* and *j* run over spatial components 1, ···, *d* = 2 in 2 + 1 dimensions. The 3 × 3 matrix <sup>Δ</sup>(*<sup>x</sup>*, *y*) is,

$$
\Delta(\mathbf{x}, \mathbf{y}) \quad = \begin{bmatrix}
\Delta\_{-1-1}(\mathbf{x}, \mathbf{y}) & \Delta\_{-10}(\mathbf{x}, \mathbf{y}) & \Delta\_{-11}(\mathbf{x}, \mathbf{y}) \\
\Delta\_{0-1}(\mathbf{x}, \mathbf{y}) & \Delta\_{00}(\mathbf{x}, \mathbf{y}) & \Delta\_{01}(\mathbf{x}, \mathbf{y}) \\
\Delta\_{1-1}(\mathbf{x}, \mathbf{y}) & \Delta\_{10}(\mathbf{x}, \mathbf{y}) & \Delta\_{11}(\mathbf{x}, \mathbf{y})
\end{bmatrix}, \tag{11}
$$

where <sup>Δ</sup>−<sup>10</sup>(*<sup>x</sup>*, *y*) = <sup>T</sup>C *δψ*−<sup>1</sup>(*x*)*δψ*<sup>∗</sup><sup>0</sup> (*y*) with time-ordered product TC in the closed-time path contour. The Green's function of dipole fields <sup>Δ</sup>−<sup>10</sup>(*<sup>x</sup>*, *y*) is also written by the 2 × 2 matrix Δ*ab*−<sup>10</sup>(*<sup>x</sup>*, *y*) with *a*, *b* = 1, 2 in the contour. The Green's function for photon fields *Dij*(*<sup>x</sup>*, *y*) represents,

$$D\_{ij}(\mathbf{x}, \mathbf{y}) = \langle \mathbf{T}\_{\mathcal{C}} a\_i(\mathbf{x}) a\_j(\mathbf{y}) \rangle. \tag{12}$$

**Figure 1.** Closed-time path contour C. The label "1" represents the path from *t*0 to ∞ and the label "2" represents the path from ∞ to *t*0.

*Entropy* **2019**, *21*, 1066

Finally we write time evolution equations for coherent fields and quantum fluctuations. The 2PI effective action satisfies the following equations,

$$\left.\frac{\delta\Gamma\_{\rm 2PI}}{\delta\Delta}\right|\_{\mathfrak{A}=0} = \left.0\right|\_{\mathfrak{A}=0} \tag{13}$$

$$
\left.\frac{\delta\Gamma\_{2\text{PI}}}{\delta D}\right|\_{\mathfrak{A}=0} = \left.0\right|\tag{14}
$$

$$
\left.\frac{\delta\Gamma\_{\rm 2PI}}{\delta\delta^i}\right|\_{\mathfrak{A}=0} = \left.\frac{\delta\Gamma\_{\rm 2PI}}{\delta A^i}\right|\_{\mathfrak{A}=0} = 0,\tag{15}
$$

$$\left.\frac{\delta\Gamma\_{2\text{Pl}}}{\delta\bar{\Psi}\_{-1,0,1}^{(\*)}}\right|\_{\mathfrak{d}=0} = \left.0\right|\tag{16}$$

due to the Legendre transformation of the generating functional. Equation (13) is written by,

$$i\Delta\_0^{-1} - i\Delta^{-1} - i\Sigma = 0,\tag{17}$$

with *i***Σ** ≡ −12 *δ*Γ2 *δ***Δ** . The matrix of self-energy **Σ** can be written by diagonal elements,

$$\Delta = \text{diag}(\Sigma\_{-1-1}, \Sigma\_{00}, \Sigma\_{11})\_\prime \tag{18}$$

since we can neglect the off-diagonal elements which are higher order of the coupling expansion. Equation (17) represents the Kadanoff–Baym equations for electric dipole fields in the two-energy-level approximation in 2 + 1 dimensions. Similarly, the Kadanoff–Baym equation for photon fields in Equation (14) is written by,

$$i D\_0^{-1} - i D^{-1} - i \Pi = 0,\tag{19}$$

with *i*Π ≡ −*δ*Γ<sup>2</sup> *δD* . Equation (15) is given by,

$$
\partial^{\nu} F\_{\nu i} = f\_{i\nu} \tag{20}
$$

with,

$$\delta f\_1(\mathbf{x}) \quad = \quad -\varepsilon d\_\varepsilon \frac{\partial}{\partial \mathbf{x}^0} \sum\_{\mathbf{a} = -1, 1} \left( \Delta\_{0\mathbf{a}}(\mathbf{x}, \mathbf{x}) + \Delta\_{\mathbf{a}0}(\mathbf{x}, \mathbf{x}) + \bar{\Psi}(\mathbf{x})\bar{\Psi}\_\mathbf{a}^\*(\mathbf{x}) + \bar{\Psi}\_\mathbf{a}(\mathbf{x})\bar{\Psi}\_\mathbf{0}^\*(\mathbf{x}) \right), \tag{21}$$

$$J\_2(\mathbf{x}) \quad = \left. -\varepsilon d\_\varepsilon \frac{\partial}{\partial \mathbf{x}^0} \sum\_{a=-1,1} \left( -ia(\Delta\_{0a}(\mathbf{x}, \mathbf{x}) - \Delta\_{00}(\mathbf{x}, \mathbf{x}) + \bar{\Psi}\_0(\mathbf{x})\bar{\Psi}\_a^\*(\mathbf{x}) - \bar{\Psi}\_a(\mathbf{x})\bar{\Psi}\_0^\*(\mathbf{x})) \right) \right. \tag{22}$$

Equation (20) represents the Klein–Gordon equations for spatial dimensions *i* = 1 and 2. Equation (16) is written by,

$$\left(i\frac{\partial}{\partial x^0} + \frac{\nabla\_i^2}{2m}\right)\Psi\_0 + \sum\_{a=-1,1} \epsilon d\_\epsilon (E\_1 + iaE\_2)\Psi\_a = 0,\tag{23}$$

$$\left(i\frac{\partial}{\partial x^0} + \frac{\nabla\_i^2}{2m} - \frac{1}{2I}\right)\psi\_a + ed\_c(E\_1 - iaE\_2)\psi\_0 = \,\_0\tag{24}$$

and their complex conjugates. They are Schrödinger-like equations for coherent dipole fields. Equations (23) and (24) and their complex conjugates give the following probability conservation,

$$\frac{\partial}{\partial x^0} \left( \bar{\psi}\_0^\* \bar{\psi}\_0 + \sum\_{a=-1,1} \bar{\psi}\_a^\* \bar{\psi}\_a \right) + \frac{1}{2mi} \nabla\_i \left( \bar{\psi}\_0^\* \nabla\_i \bar{\psi}\_0 - \bar{\psi}\_0 \nabla\_i \bar{\psi}\_0^\* + \sum\_{a=-1,1} \left( \bar{\psi}\_a^\* \nabla\_i \bar{\psi}\_a - \bar{\psi}\_a \nabla\_i \bar{\psi}\_a^\* \right) \right) = 0. \tag{25}$$

We shall define *J*0(*x*) as,

$$J\_0(\mathbf{x}) = -ed\_\varepsilon \frac{\partial}{\partial \mathbf{x}^1} \sum\_{\mathbf{a} = -1, 1} \left( \Delta\_{0\mathbf{a}}(\mathbf{x}, \mathbf{x}) + \Delta\_{\mathbf{a}0}(\mathbf{x}, \mathbf{x}) + \bar{\psi}\_0(\mathbf{x}) \bar{\psi}\_\mathbf{a}^\*(\mathbf{x}) + \bar{\psi}\_\mathbf{a}(\mathbf{x}) \bar{\psi}\_0^\*(\mathbf{x}) \right)$$

$$-ed\_\varepsilon \frac{\partial}{\partial \mathbf{x}^2} \left( -ia(\Delta\_{0\mathbf{a}}(\mathbf{x}, \mathbf{x}) - \Delta\_{\mathbf{a}0}(\mathbf{x}, \mathbf{x}) + \bar{\psi}\_0(\mathbf{x}) \bar{\psi}\_\mathbf{a}^\*(\mathbf{x}) - \bar{\psi}\_\mathbf{a}(\mathbf{x}) \bar{\psi}\_0^\*(\mathbf{x})) \right). \tag{26}$$

Then since we can use *∂*0 *J*0 − ∇*iJi* = 0 with *i* = 1, 2,

$$
\delta \partial\_0 l\_0 = \nabla\_{\bar{i}} l\_{\bar{i}} = -\partial^{\bar{i}} \partial^{\nu} F\_{\nu \bar{i}} = \partial^{\mu} \partial^{\nu} F\_{\nu \mu} - \partial^{i} \partial^{\nu} F\_{\nu \bar{i}} = \partial^{0} \partial^{\nu} F\_{\nu 0 \nu}
$$

$$
\text{or, } \partial^{\nu} F\_{\nu 0} = f\_{0 \nu} \tag{27}
$$

where the time dependent term in the time integral might be interpreted as an initial charge, but it is set to be zero. This equation represents the Poisson equation for scalar potential *A*<sup>0</sup> given by ∇2*A*<sup>0</sup> = ∇ · *μ* with the vector of dipole moments −*μ* on the right-hand side in Equation (26). (Since the Fourier transformed *A* ˜ <sup>0</sup>(**q**) is written by *A*˜ <sup>0</sup>(**q**) ∝ (*qiμ*˜*i*)/**q**<sup>2</sup> with *μi* = *μ*˜*i<sup>δ</sup>*(**r**), the electric field *Ej* = −∇*jA*<sup>0</sup>(**r**) is proportional to **q** *ei***q**·**<sup>r</sup>** *<sup>q</sup>jqiμ*˜*i* **q**2 . If we can also apply the analysis in this section to the case in 3 + 1 dimensions, we find *Ej* ∝ *<sup>∂</sup>j∂i μ* ˜ *i r* . Then we obtain dipole–dipole interaction potential −*μ* ˜ *jEj* ∼ *<sup>μ</sup>*˜*jμ*˜*j r*3 − <sup>3</sup>(*riμ*˜*i*)(*rjμ*˜*j*) *r*5 in 3 + 1 dimensions.)

#### **3. Kinetic Entropy Current in the Kadanoff–Baym Equations and the H-Theorem**

In this section, we derive a kinetic entropy current from the Kadanoff–Baym equations with first order approximation of the gradient expansion and show the H-theorem for the leading-order approximations in the coupling expansion based on [56–58]. The analysis in this section is similar to that in open systems (the central region connected to the left and the right region) [59]. Since (−1, 1) and (1, −<sup>1</sup>) components in *i***Δ**−<sup>1</sup> 0 (*<sup>x</sup>*, *y*) in Equation (9) are zero, the same procedures to rewrite the Kadanoff–Baym equations as those in open systems [59–63] can be adopted. We set *t*0 → <sup>−</sup>∞.

First, we shall write the Kadanoff–Baym equations in Equation (17) for each components. By multiplying the matrix **Δ** from the right in Equation (17) and taking the (0, 0) component, we can write it as,

$$\mathrm{i}\left(\Delta\_{0,00}^{-1} - \Sigma\_{00}\right)\Delta\_{00} + \sum\_{a=-1,1} \mathrm{ed}\_{\mathfrak{c}}(E\_1 + iaE\_2)\Delta\_{a0} = \mathrm{i}\delta\_{\mathcal{C}}(\mathbf{x} - \mathbf{y})\_{\prime} \tag{28}$$

where the (0, 0) component of the matrix **Δ**−<sup>1</sup> 0 represents *i*Δ−<sup>1</sup> 0,00(*<sup>x</sup>*, *y*) = *i ∂∂x*0 + ∇2*i* 2*m δ*C (*x* − *y*). By taking (*<sup>α</sup>*, 0) component, we can write it as,

$$\mathrm{ch}(\Delta\_{0,aa}^{-1} - \Sigma\_{aa})\Delta\_{a0} + \mathrm{ed}\_{\mathfrak{c}}(E\_1 - i\mathfrak{a}E\_2)\Delta\_{00} = 0.\tag{29}$$

It is convenient to introduce the Green's functions <sup>Δ</sup>*g*,*αα* as,

$$i\Delta\_{\mathfrak{J},aa}^{-1} = i\Delta\_{0,aa}^{-1} - i\Sigma\_{\mathfrak{a}\mathfrak{a}}.\tag{30}$$

Then by using Equations (29) and (30), we can write Δ*α*<sup>0</sup> as,

$$
\Delta\_{\rm td}(\mathbf{x}, \mathbf{y}) = -\frac{\varepsilon d\_{\varepsilon}}{i} \int\_{\mathcal{C}} dw \Delta\_{\rm g,aa}(\mathbf{x}, w)(E\_1(w) - iaE\_2(w))\Delta\_{00}(w, \mathbf{y}).\tag{31}
$$

Equation (31) means the propagation from *y* to *x* with zero angular momentum, change of angular momentum at *w* and the propagation from *w* to *x* with angular momentum *α* = ±1. By using Equation (31), we can rewrite Equation (28) as,

$$\begin{aligned} i \int\_{\mathcal{C}} dw (\Lambda\_{0,0}^{-1}(\mathbf{x}, w) - \Sigma\_{00}(\mathbf{x}, w)) \Lambda\_{00}(w, y) \\ i + i \sum\_{\mathbf{z} = -1, 1} (cd\_{\mathbf{z}})^2 \int\_{\mathcal{C}} dw (E\_1(\mathbf{x}) + iaE\_2(\mathbf{x})) \Lambda\_{\mathbf{g}, \mathbf{z}a}(\mathbf{x}, w) (E\_1(w) - iaE\_2(w)) \Lambda\_{00}(w, y) = i \delta\_{\mathcal{C}}(\mathbf{x} - y). \end{aligned} \tag{32}$$

The second term on the left-hand side in Equation (32) represents the propagation from *y* to *w* with zero angular momentum, the change of the angular momentum to *α* = ±1 at *w* due to the coherent electric fields, the propagation from *w* to *x* and the change of the angular momentum from *α* = ±1 to zero due to the coherent electric fields. In a similar way to *φ*4 theory in open systems [59], we can derive,

$$\begin{aligned} \mathrm{i} \int\_{\mathcal{C}} dw \Delta \alpha (\mathbf{x}, w) (\Delta\_{0, 00}^{-1} (w, y) - \Sigma \alpha (w, y)) \\ + \mathrm{i} \sum\_{w = -1, 1} (ed\_{\mathcal{C}})^2 \int\_{\mathcal{C}} dw \Delta \alpha (\mathbf{x}, w) (E\_1(w) + iaE\_2(w)) \Delta\_{\text{g.a.}} (w, y) (E\_1(y) - iaE\_2(y)) = \mathrm{i} \delta\_{\mathcal{C}} (\mathbf{x} - y), \end{aligned} \tag{33}$$

where we have used,

$$
\Delta \alpha\_{\rm tr}(\mathbf{x}, \mathbf{y}) = -\frac{1}{i} \int\_{\mathcal{C}} dw \Delta \alpha\_{\rm 0}(\mathbf{x}, w) (ed\_{\rm t}) (E\_1(w) + iaE\_2(w)) \Delta\_{\rm g.a.t}(w, \mathbf{y}). \tag{34}
$$

The (*<sup>α</sup>*, *α*) components of the Kadanoff–Baym equations are written by,

$$\begin{aligned} \mathrm{i} \int\_{\mathcal{C}} dw \left( \Delta\_{0,\mathrm{ad}}^{-1}(\mathbf{x}, w) - \Sigma\_{\mathrm{ad}}(\mathbf{x}, w) \right) \Delta\_{\mathrm{ad}}(w, y) \\ \mathrm{i} + \mathrm{i}(\mathrm{ad}\_{\mathfrak{C}})^2 \int\_{\mathcal{C}} dw (E\_1(\mathbf{x}) - i\mathrm{a}E\_2(\mathbf{x})) \Delta\_{\mathrm{00}}(\mathbf{x}, w) (E\_1(w) + i\mathrm{a}E\_2(w)) \Delta\_{\mathrm{\mathcal{S},\mathrm{ad}}}(w, y) &= i\delta\_{\mathcal{C}}(\mathbf{x} - y), \end{aligned} \tag{35}$$

and,

$$\begin{aligned} i \int\_{\mathcal{C}} dw \Delta\_{\text{nat}}(\mathbf{x}, w) \left( \Delta\_{0, \text{aa}}^{-1}(w, y) - \Sigma\_{\text{nat}}(w, y) \right) \\ + i (\text{ed}\_{\mathbb{C}})^2 \int\_{\mathcal{C}} dw \Delta\_{\text{g.a.}}(\mathbf{x}, w) (\mathcal{E}\_1(w) - i a \mathcal{E}\_2(w)) \Delta\_{00}(w, x) (\mathcal{E}\_1(\mathbf{x}) + i a \mathcal{E}\_2(\mathbf{x})) &= i \delta\_{\mathbb{C}}(\mathbf{x} - y), \end{aligned} \tag{36}$$

where we have used Equations (31) and (34).

Next, we shall perform the Fourier transformation ( *d*(*x* − *y*)*eip*·(*<sup>x</sup>*−*y*)) with the relative coordinate *x* − *y* of the (0, 0) and (*<sup>α</sup>*, *α*) components of the Kadanoff–Baym equations. We use the 2 × 2 matrix notation in the closed-time path with *a*, *b*, *c*, *d* = 1, 2. Equations (32) and (33) are transformed as,

$$\mathrm{id}\left(\Delta\_{0,00}^{-1}(p) - \Sigma\_{00}(X,p)\sigma\_z + \sum\_a \mathrm{l}\!\!\!L\_{aa}(X,p)\sigma\_z\right)^{ac} \circ \Delta\_{00}^{cb}(X,p) = \mathrm{i}\sigma\_z^{ab},\tag{37}$$

$$i\Lambda\_{00}^{ac}(X,p)\circ\left(\Lambda\_{0,00}^{-1}(p)-\sigma\_z\Sigma\_{00}(X,p)+\sigma\_z\sum\_a \mathcal{U}\_{aa}(X,p)\right)^{cb}=i\sigma\_z^{ab},\tag{38}$$

*Entropy* **2019**, *21*, 1066

where *X* = *x*+*y* 2 , *σz* = diag(1, <sup>−</sup><sup>1</sup>),

$$i\Delta\_{0,00}^{-1}(p) = p^0 - \frac{\mathbf{p}^2}{2m'} \tag{39}$$

and the *<sup>U</sup>αα*(*<sup>X</sup>*, *p*) is the Fourier transformation,

$$\begin{split} \mathrm{d}\mathrm{Lan}(X,p) &= \, ^1\mathrm{ed}\, ^2\int d(\mathbf{x}-\mathbf{y}) \epsilon^{ip\cdot(\mathbf{x}-\mathbf{y})} (E\_1(\mathbf{x}) + \mathrm{ia}E\_2(\mathbf{x})) \Delta\_{\mathrm{g},\mathrm{ad}}(\mathbf{x},\mathbf{y}) (E\_1(\mathbf{y}) - \mathrm{ia}E\_2(\mathbf{y})) \\ &= \, ^1\mathrm{ed}\, ^2\mathrm{E}(X)^2 \Delta\_{\mathrm{g},\mathrm{ad}}(X,p+\mathrm{a}\partial\_{\mathbf{y}}^{\prime}) + \left(\frac{\partial^2}{\partial X^2}\right) \, , \end{split} \tag{40}$$

with the definition of *ζ* and |**E**|,

$$E\_1(\mathbf{x}) + i\mathfrak{a}E\_2(\mathbf{x}) = |\mathbf{E}(\mathbf{x})|e^{i\mathbf{a}\tilde{\gamma}(\mathbf{x})},\tag{41}$$

and,

$$\left(\mathcal{U}\_{aa}(X,p)\sigma\_z\right)^{ac} = \mathcal{U}\_{aa}^{ad}(X,p)\sigma\_z^{dc} \,,\tag{42}$$

The ◦ is expanded by the derivative of *X* [64–67] as,

$$H(X,p)\circ I(X,p) = H(X,p)I(X,p) + \frac{i}{2}\left\{H,I\right\} + \left(\frac{\partial^2}{\partial X^2}\right),\tag{43}$$

with the definition of the Poisson bracket,

$$\{H, I\} \equiv \frac{\partial}{\partial p^{\mu}} \frac{\partial}{\partial X\_{\mu}} - \frac{\partial}{\partial X^{\mu}} \frac{\partial}{\partial p\_{\mu}}.\tag{44}$$

We find that the *Uαα* represents the change of momenta of dipoles as shown in Figure 2a.

**Figure 2.** Diagrams of (**a**) *<sup>U</sup>αα*(*<sup>X</sup>*, *p*) and (**b**) self-energy <sup>Σ</sup>00(*<sup>X</sup>*, *p*).

In a similar way to [59], in the 0th and the first order in the gradient expansion in Equations (37) and (38), we can derive the following retarded Green's function,

$$\Delta\_{00,R}(X,p) = \frac{-1}{p^0 - \frac{\mathbf{p}^2}{2m} - \Sigma\_{00,R} + \sum\_{a=-1,1} \mathcal{U}\_{aa,R}},\tag{45}$$

with the retarded parts (the subscript '*R*') Δ00,*R* = *i*(Δ1100 − <sup>Δ</sup>1200), Σ00,*R* = *i*(Σ1100 − Σ1200) and *Uαα*,*<sup>R</sup>* = *<sup>i</sup>*(*U*11*αα* − *<sup>U</sup>*12*αα*). By taking the imaginary part of the retarded Green's function <sup>Δ</sup>00,*<sup>R</sup>*(*<sup>X</sup>*, *p*), we can derive the spectral function *ρ*00 = *i*(Δ2100 − Δ1200) = <sup>2</sup>*i*ImΔ00,*<sup>R</sup>*(*<sup>X</sup>*, *p*) which represents the information of dispersion relations. Similarly, the (*<sup>α</sup>*, *α*) components of the Kadanoff–Baym equations are written as,

$$i\left(\Delta\_{0,\text{aa}}^{-1}(p) - \Sigma\_{\text{aa}}(X,p)\sigma\_z\right) \circ \Delta\_{\text{aa}}(X,p) + iV\_{\text{aa}}(X,p)\sigma\_z \circ \Delta\_{\text{\textdegree{g}},\text{aa}}(X,p) = i\sigma\_z. \tag{46}$$

*Entropy* **2019**, *21*, 1066

and,

$$i\Delta\_{\text{flat}}(X,p)\diamond \left(\Delta\_{0,\text{aa}}^{-1}(p) - \sigma\_{\text{z}}\Sigma\_{\text{aa}}(X,p)\right) + i\Delta\_{\text{g.out}}(X,p)\diamond \sigma\_{\text{z}}V\_{\text{aa}}(X,p) = i\sigma\_{\text{z}}\tag{47}$$

where,

$$i\Delta\_{0,aa}^{-1}(p) = p^0 - \frac{\mathbf{p}^2}{2m} - \frac{1}{2I},\tag{48}$$

and,

$$V\_{\rm an}(X, p) = (\epsilon d\_{\epsilon})^2 \int d(\mathbf{x} - y) \epsilon^{ip \cdot (\mathbf{x} - y)} (E\_1(\mathbf{x}) - iaE\_2(\mathbf{x})) \Delta\_{00}(\mathbf{x}, y) (E\_1(y) + iaE\_2(y))$$

$$= (\epsilon d\_{\epsilon})^2 \mathbf{E}(X)^2 \Delta\_{00}(X, p - a\partial \zeta) + \left(\frac{\partial^2}{\partial X^2}\right) \,. \tag{49}$$

We can also write for <sup>Δ</sup>*cbg*,*αα*(*<sup>X</sup>*, *p*) as,

$$\mathrm{id}\left(\Delta\_{0,aa}^{-1}(p) - \Sigma\_{aa}(X\_\prime p)\sigma\_z\right)^{ac} \circ \Delta\_{\mathrm{g,aa}}^{cb}(X\_\prime p) \quad = \quad i\sigma\_z^{ab},\tag{50}$$

$$\Delta\_{\rm g.aa}^{\rm ac}(X, p) \diamond i \left(\Delta\_{0,aa}^{-1}(p) - \sigma\_z \Sigma\_{\rm aa}(X, p)\right)^{cb} = \quad \rm i\sigma\_z^{ab}.\tag{51}$$

In the 0th and the first order in the gradient expansion in Equations (46) and (47), we can derive,

$$
\Delta\_{\text{na},R} = \Delta\_{\text{\\$},\text{na},R} + \Delta\_{\text{\\$},\text{na},R} V\_{\text{an},R} \Delta\_{\text{\\$},\text{na},R} \tag{52}
$$

with Δ*αα*,*<sup>R</sup>* = *<sup>i</sup>*(Δ11*αα* − <sup>Δ</sup>12*αα*) and *Vαα*,*<sup>R</sup>* = *i*(*V*<sup>11</sup> *αα* − *V*<sup>12</sup> *αα* ). Here we have used the solution in the 0th and the first order in the gradient expansion in Equations (50) and (51) given by,

$$\Delta\_{\mathbb{S},aa,R} = \frac{-1}{p^0 - \frac{p^2}{2m} - \frac{1}{27} - \Sigma\_{aa,R}},$$

with Σ*αα*,*<sup>R</sup>* = *<sup>i</sup>*(Σ11*αα* − <sup>Σ</sup>12*αα*). The derivation is the same as [59]. The imaginary part of the retarded Green's function <sup>Δ</sup>*αα*,*<sup>R</sup>*(*<sup>X</sup>*, *p*) multiplied by 2*i* represents the spectral function *ραα* = *<sup>i</sup>*(Δ21*αα* − <sup>Δ</sup>12*αα*) = <sup>2</sup>*i*ImΔ*αα*,*<sup>R</sup>*(*<sup>X</sup>*, *p*) which represents the information of dispersion relations. In addition, the Kadanoff–Baym equations for photons in Equation (19) are written by,

$$\mathrm{id}\left(D\_{0,ij}^{-1}(k) - \Pi\_{ij}(X,k)\sigma\_z\right)^{ac} \circ D\_{j\mathrm{l}}^{cb}(X,k) \quad = \quad i\delta\_{\mathrm{il}}\sigma\_z^{ab},\tag{54}$$

$$i D\_{ij}^{a\varepsilon}(X,k) \diamond \left( D\_{0,jl}^{-1}(k) - \sigma\_{\overline{z}} \Pi\_{jl}(X,k) \right)^{\varepsilon b} = \quad i \delta\_{il} \sigma\_z^{ab} \tag{55}$$

with,

$$i D\_{0,ij}^{-1}(k) = k^2 \delta\_{ij}.\tag{56}$$

Next we shall derive the self-energy in the leading-order (LO) of the coupling expansion in Equation (6). The (*a*, *b*)=(1, 2) and (2, 1) component of *i* Γ22 are given by,

$$\begin{split} \Lambda^{\frac{\Gamma\_{110}}{2}} &= -\frac{1}{2} (\boldsymbol{e} \boldsymbol{l}\_{\varepsilon})^{2} \int \boldsymbol{d} \boldsymbol{u} \boldsymbol{w} \sum\_{a=-1,1} \left( \Lambda^{21}\_{\rm{at}} (\boldsymbol{w}, \boldsymbol{u}) \Lambda^{12}\_{00} (\boldsymbol{u}, \boldsymbol{w}) (1, -\boldsymbol{u} \boldsymbol{i}) \right) \partial^{0}\_{\boldsymbol{u}} \partial^{0}\_{\boldsymbol{w}} \left( D^{12}\_{\parallel l} (\boldsymbol{u}, \boldsymbol{w}) + D^{21}\_{l\parallel} (\boldsymbol{w}, \boldsymbol{u}) \right) (1, \boldsymbol{u} \boldsymbol{i})^{\dagger}\_{l} \\ &+ \Lambda^{12}\_{\rm{at}} (\boldsymbol{w}, \boldsymbol{u}) \Delta^{21}\_{00} (\boldsymbol{u}, \boldsymbol{w}) (1, -\boldsymbol{u} \boldsymbol{i}) \, \partial^{0}\_{\boldsymbol{u}} \partial^{0}\_{\boldsymbol{w}} \left( D^{21}\_{\parallel l} (\boldsymbol{u}, \boldsymbol{w}) + D^{12}\_{l\parallel} (\boldsymbol{w}, \boldsymbol{u}) \right) (1, \boldsymbol{u} \boldsymbol{i})^{\dagger}\_{l}, \end{split} \tag{57}$$

where *t* represents the transposition. It is convenient to rewrite,

$$D\_{ij}^{ab}(k) \quad = \left(\delta\_{ij} - \frac{k\_i k\_j}{\mathbf{k}^2}\right) D\_T^{ab}(k) + \frac{k\_i k\_j}{\mathbf{k}^2} D\_L^{ab}(k), \tag{58}$$

$$
\Pi\_{ij}^{ab}(k) \quad = \begin{pmatrix} \delta\_{ij} - \frac{k\_i k\_j}{\mathbf{k}^2} \end{pmatrix} \Pi\_T^{ab}(k) + \frac{k\_i k\_j}{\mathbf{k}^2} \Pi\_L^{ab}(k), \tag{59}
$$

where *T* and *L* represent the transverse and the longitudinal part, respectively. The LO self-energy *<sup>i</sup>*Π21*ji* (*y*, *x*) = − *<sup>δ</sup>*Γ2,LO *<sup>δ</sup>D*12*ij*(*<sup>x</sup>*,*y*) is,

$$i\Pi\_{jl}^{21}(y,x) = -i(cd\_{\boldsymbol{\ell}})^2 \sum\_{a=-1,1} \left(\partial\_x^0 \partial\_y^0 \left(\Delta\_{aa}^{21}(y,x)\Delta\_{00}^{12}(x,y)\right)(1,-ai)\_l(1,ai)\_j^t\right)$$

$$+\partial\_x^0 \partial\_y^0 \left(\Delta\_{00}^{21}(y,x)\Delta\_{aa}^{12}(x,y)\right)(1,-ai)\_j(1,ai)\_l^t\Big).\tag{60}$$

By Fourier-transforming with the relative coordinate *x* − *y* and multiplying *δij* − *kikj* **k**<sup>2</sup> or *kikj* **k**<sup>2</sup> , we arrive at,

$$\begin{split} \Pi\_{T}^{21}(\mathbf{X},k) &= -(\mathrm{c}\mathbf{d}\_{\varepsilon})^{2} \left(\mathbf{k}^{0}\right)^{2} \int\_{\mathbb{P}\_{\mathbf{z}\mathbf{a}}} \sum\_{\mathbf{a}=-1,1} \left(\Delta\_{\mathbf{a}\mathbf{a}}^{21}(\mathbf{X},k+p)\Lambda\_{00}^{12}(\mathbf{X},p) + \Delta\_{00}^{21}(\mathbf{X},k+p)\Lambda\_{\mathbf{a}\mathbf{a}}^{12}(\mathbf{X},p)\right) \\ &+ \left(\frac{\partial^{2}}{\partial X^{2}}\right), \end{split} \tag{61}$$

$$\begin{array}{rcl} \Pi\_L^{21}(X,k) & = & \Pi\_T^{21}(X,k), \\ \end{array} \tag{62}$$

with *p* = *dd*+<sup>1</sup> *p* (<sup>2</sup>*π*)*<sup>d</sup>*+<sup>1</sup> . The second equation is due to the spatial dimension *d* = 2. Similarly, we arrive at,

$$\begin{split} \|\Pi\|\_{\mathrm{T}}^{2}(\mathcal{X},k) &= -(\mathrm{c}\mathrm{d}\mathrm{c})^{2} \left(\mathrm{k}^{0}\right)^{2} \int\_{\mathbb{P}\_{\mathrm{ad}}} \sum\_{a=-1,1} \left(\Delta\_{\mathrm{ad}}^{12}(\mathcal{X},k+p)\Delta\_{00}^{21}(\mathcal{X},p) + \Delta\_{00}^{12}(\mathcal{X},k+p)\Delta\_{\mathrm{ad}}^{21}(\mathcal{X},p)\right) \\ &+ \left(\frac{\partial^{2}}{\partial\mathcal{X}^{2}}\right), \\ \|\Delta\_{\mathrm{ad}}^{-21}(\mathcal{X},k)\|\_{\mathrm{L}}^{2} &= -\Delta\_{\mathrm{ad}}^{-21}(\mathcal{X},k). \end{split} \tag{63}$$

$$
\Pi\_L^{12}(X,k) \quad = \quad \Pi\_T^{12}(X,k). \tag{64}
$$

The Fourier transformation of the LO self-energy *<sup>i</sup>*Σ1200(*<sup>x</sup>*, *y*) = −12 *<sup>δ</sup>*Γ2,LO *<sup>δ</sup>*Δ2100(*y*,*<sup>x</sup>*) is,

$$\Sigma\_{00}^{12}(X,p) = -(cd\_\ell)^2 \int\_{k\_{\rm id}} \sum\_{a=-1,1} \left(k^0\right)^2 \Delta\_{\rm tot}^{12}(X,p-k) \left[D\_T^{12}(X,k) + D\_L^{12}(X,k)\right] + \left(\frac{\partial^2}{\partial X^2}\right). \tag{65}$$

Similarly,

$$\Sigma\_{00}^{21}(X,p) = -(\epsilon d\_{\epsilon})^2 \int\_{k} \sum\_{a=-1,1} \left(k^0\right)^2 \Delta\_{\text{tot}}^{21}(X,p-k) \left[D\_T^{21}(X,k) + D\_L^{21}(X,k)\right] + \left(\frac{\partial^2}{\partial X^2}\right). \tag{66}$$

This self-energy is shown in Figure 2b. Similarly we can derive,

$$
\Sigma\_{\rm aa}^{12}(X, p) = -(\epsilon d\_{\varepsilon})^2 \int\_{k} \left(k^0\right)^2 \Lambda\_{00}^{12}(X, p - k) \left[D\_T^{12}(X, k) + D\_L^{12}(X, k)\right] + \left(\frac{\bar{\partial}^2}{\partial X^2}\right), \tag{67}
$$

and,

$$
\Delta\_{\rm rad}^{21}(X, p) = -(\alpha d\_{\epsilon})^2 \int\_{k} \left(k^0\right)^2 \Delta\_{00}^{21}(X, p - k) \left[D\_T^{21}(X, k) + D\_L^{21}(X, k)\right] + \left(\frac{\partial^2}{\partial X^2}\right). \tag{68}
$$

Finally we derive a kinetic entropy current in the first order approximation in the gradient expansion and show the H-theorem in the LO approximation in the coupling expansion. By taking a difference of Equations (32) and (33), we arrive at,

$$i\left\{p^0 - \frac{\mathbf{p}^2}{2m'}\Delta\_{00}^{ab}\right\} = i\left[\left(\Sigma\_{00} - \sum\_a \mathcal{U}\_{aa}\right)\sigma\_z \diamond \Delta\_{00}\right]^{ab} - i\left[\Delta\_{00} \diamond \sigma\_z \left(\Sigma\_{00} - \sum\_a \mathcal{U}\_{aa}\right)\right]^{ab}.\tag{69}$$

We use the Kadanoff–Baym Ansatz Δ1200 = *ρ*00*i f*00, Δ2100 = *ρ*00*i* (*f*00 + <sup>1</sup>), Σ1200 = <sup>Σ</sup>00,*ρ i γ*00, Σ2100 = <sup>Σ</sup>00,*ρ i* (*<sup>γ</sup>*00 + <sup>1</sup>), *<sup>U</sup>*12*αα* = *<sup>U</sup>αα*,*<sup>ρ</sup> i γU*,*αα* and *<sup>U</sup>*21*αα* = *<sup>U</sup>αα*,*<sup>ρ</sup> i* (*<sup>γ</sup>U*,*αα* + 1) with *ρ*00 = *i*(Δ2100 − Δ1200) = <sup>2</sup>*i*ImΔ00,*R*, <sup>Σ</sup>00,*ρ* = *i*(Σ2100 − Σ1200) = <sup>2</sup>*i*ImΣ00,*R* and *<sup>U</sup>αα*,*<sup>ρ</sup>* = *<sup>i</sup>*(*U*21*αα* − *<sup>U</sup>*12*αα*) = <sup>2</sup>*i*Im*Uαα*,*<sup>R</sup>* where we just rewrite the (1, 2) and the (2, 1) components with the spectral parts *ρ*00, <sup>Σ</sup>00,*ρ* and *<sup>U</sup>αα*,*<sup>ρ</sup>* and distribution functions *f*00, *γ*00 and *γU*,*αα*. The distribution functions *f*00, *γ*00 and *γU*,*αα* approach the Bose–Einstein distributions near equilibrium states. In the first order approximation in the gradient expansion in Equation (69) for (*a*, *b*)=(1, 2) and (2, <sup>1</sup>), we can derive,

$$f\_{00} = \gamma\_{00} + O\left(\frac{\partial}{\partial X}\right), \quad \text{and} \quad f\_{00} = \gamma\_{\mathcal{U}, \text{at}} + O\left(\frac{\partial}{\partial X}\right). \tag{70}$$

(Rewrite (*a*, *b*)=(1, 2) and (2, 1) components in Equation (69), then we can show the collision terms Δ2100Σ1200 − Δ1200Σ2100 ∝ *f*00 − *γ*00 = *O ∂∂X* and *f*00 − *γU*,*αα* = *O ∂∂X* .) We shall multiply ln *i*Δ1200 *ρ*00 in (*a*, *b*)=(1, 2) component in Equation (69) and ln *i*Δ2100 *ρ*00 in (2, 1) component in Equation (69), take the difference of them and integrate with *p*. By the use of Equation (70), we arrive at,

$$\begin{split} \partial\_{\boldsymbol{p}}s^{\boldsymbol{\Pi}}\_{\text{matter,00}} &= \ -\int\_{\boldsymbol{p}} \left( \Sigma^{21}\_{00}(\boldsymbol{X},\boldsymbol{p}) \Delta^{12}\_{00}(\boldsymbol{X},\boldsymbol{p}) - \Sigma^{12}\_{00}(\boldsymbol{X},\boldsymbol{p}) \Delta^{21}\_{00}(\boldsymbol{X},\boldsymbol{p}) \right) \ln \frac{\Delta^{12}\_{00}(\boldsymbol{X},\boldsymbol{p})}{\Delta^{21}\_{00}(\boldsymbol{X},\boldsymbol{p})} \\ &+ \sum\_{\mathbf{a}} \int\_{\boldsymbol{p}} \left( \mathcal{U}^{21}\_{\mathbf{a}\mathbf{a}}(\boldsymbol{X},\boldsymbol{p}) \Delta^{12}\_{00}(\boldsymbol{X},\boldsymbol{p}) - \mathcal{U}^{12}\_{\mathbf{a}\mathbf{a}}(\boldsymbol{X},\boldsymbol{p}) \Delta^{21}\_{00}(\boldsymbol{X},\boldsymbol{p}) \right) \ln \frac{\Delta^{12}\_{00}(\boldsymbol{X},\boldsymbol{p})}{\Delta^{21}\_{00}(\boldsymbol{X},\boldsymbol{p})}, \end{split} \tag{71}$$

with the definition of entropy current *<sup>s</sup>μ*matter,00 for (0, 0) component,

$$s\_{\text{matter},00}^{\mu} \equiv \int\_{p} \left[ \left( \delta\_{0}^{\mu} + \frac{\delta\_{i}^{\mu} \mathbf{p}^{i}}{m} - \frac{\partial \text{Re} \left( \Sigma\_{00,R} - \sum\_{\alpha} \mathcal{U}\_{aa,R} \right)}{\partial p\_{\mu}} \right) \frac{\rho\_{00}}{i} \right]$$

$$+ \frac{\partial \text{Re} \Delta\_{00,R}}{\partial p\_{\mu}} \frac{\Sigma\_{00,\rho} - \sum\_{\alpha} \mathcal{U}\_{aa,\rho}}{i} \Big] \sigma[f\_{00}],\tag{72}$$

$$\text{where} \qquad (i\_{\text{max}} \text{ is the } i\_{\text{max}} \text{ -- the } i\_{\text{max}}).\tag{73}$$

$$\sigma[f\_{00}] \equiv \ (1 + f\_{00})\ln(1 + f\_{00}) - f\_{00}\ln f\_{00}.\tag{73}$$

We can derive the Boltzmann entropy **p** [(1 + *n*)ln(<sup>1</sup> + *n*) − *n* ln *n*] with the number density *<sup>n</sup>*(*<sup>X</sup>*, **p**) in the quasi-particle limit Im*Uαα*,*<sup>R</sup>* = ImΣ00,*R* → 0 in the same way as in [58]. Similarly, we can derive a kinetic entropy current for (*αα*) components. >From Equations (46) and (47), we can derive

$$\begin{split} \mathrm{i}\left\{\mathrm{p}^{0} - \frac{\textbf{p}^{2}}{2m} - \frac{1}{2\mathsf{I}}, \Delta\_{\mathsf{aa}}^{ab}\right\} &= \mathrm{i}\left[\Sigma\_{\mathsf{aa}}\sigma\_{z} \diamond \Delta\_{\mathsf{aa}} - \Delta\_{\mathsf{aa}} \circ \sigma\_{z}\Sigma\_{\mathsf{aa}}\right]^{ab} \\ &- \mathrm{i}\left[V\_{\mathsf{aa}}\sigma\_{z} \diamond \Delta\_{\mathsf{\mathsf{g}},\mathsf{aa}} - \Delta\_{\mathsf{\mathsf{g}},\mathsf{aa}} \circ \sigma\_{z}V\_{\mathsf{aa}}\right]^{ab}. \end{split} \tag{74}$$

We use the Kadanoff–Baym Ansatz <sup>Δ</sup>12*αα* = *ρααi fαα*, <sup>Δ</sup>21*αα* = *ρααi* (*fαα* + <sup>1</sup>), <sup>Δ</sup>12*g*,*αα* = <sup>Δ</sup>*g*,*αα*,*<sup>ρ</sup> i γg*,*αα*, <sup>Δ</sup>21*g*,*αα* = <sup>Δ</sup>*g*,*αα*,*<sup>ρ</sup> i* (*<sup>γ</sup>g*,*αα* + <sup>1</sup>), <sup>Σ</sup>12*αα* = <sup>Σ</sup>*αα*,*<sup>ρ</sup> i γαα*, <sup>Σ</sup>21*αα* = <sup>Σ</sup>*αα*,*<sup>ρ</sup> i* (*γαα* + <sup>1</sup>), *V*<sup>12</sup> *αα* = *<sup>V</sup>αα*,*<sup>ρ</sup> i γV*,*αα* and *V*<sup>21</sup> *αα* = *<sup>V</sup>αα*,*<sup>ρ</sup> i* (*<sup>γ</sup>V*,*αα* + 1) with *ραα* = *<sup>i</sup>*(Δ21*αα* − <sup>Δ</sup>12*αα*) = <sup>2</sup>*i*ImΔ*αα*,*R*, <sup>Σ</sup>*αα*,*<sup>ρ</sup>* = *<sup>i</sup>*(Σ21*αα* − <sup>Σ</sup>12*αα*) = <sup>2</sup>*i*ImΣ*αα*,*<sup>R</sup>* and *<sup>V</sup>αα*,*<sup>ρ</sup>* = *i*(*V*<sup>21</sup> *αα* − *V*<sup>12</sup> *αα* ) = <sup>2</sup>*i*Im*Vαα*,*R*. In Equation (74), we can show,

$$f\_{\mu\alpha} \sim \gamma\_{\mu\alpha\nu} \quad \gamma\_{\mathcal{K}\mu\alpha} \sim \gamma\_{\mathcal{V}\mathcal{A}\mu\nu} \tag{75}$$

for distribution functions *fαα*, *γαα* and *γV*,*αα* by writing the (*a*, *b*)=(1, 2) and (2, 1) components in the Kadanoff–Baym equations (74). We can also show,

$$
\gamma\_{\text{ann}} \sim \gamma\_{\text{g,aux}} \tag{76}
$$

from Equations (50) and (51). We shall multiply ln *<sup>i</sup>*Δ12*αα ραα* in (*a*, *b*)=(1, 2) component in Equation (74) and ln *<sup>i</sup>*Δ21*αα ραα* in (2, 1) component in Equation (74), take the difference of them and integrate with *p*. By using Equations (75) and (76), we arrive at,

$$\begin{split} \partial\_{\boldsymbol{p}}s^{\mu}\_{\text{matter,aa}} &= \ -\int\_{\boldsymbol{p}} \left( \Sigma^{21}\_{\text{aa}}(\mathcal{X},\boldsymbol{p}) \Delta^{12}\_{\text{aa}}(\mathcal{X},\boldsymbol{p}) - \Sigma^{12}\_{\text{aa}}(\mathcal{X},\boldsymbol{p}) \Delta^{21}\_{\text{aa}}(\mathcal{X},\boldsymbol{p}) \right) \ln \frac{\Delta^{12}\_{\text{aa}}(\mathcal{X},\boldsymbol{p})}{\Delta^{21}\_{\text{aa}}(\mathcal{X},\boldsymbol{p})} \\ &+ \int\_{\boldsymbol{p}} \left( V^{21}\_{\text{aa}}(\mathcal{X},\boldsymbol{p}) \Delta^{12}\_{\text{ga},\text{aa}}(\mathcal{X},\boldsymbol{p}) - V^{12}\_{\text{aa}}(\mathcal{X},\boldsymbol{p}) \Delta^{21}\_{\text{ga},\text{aa}}(\mathcal{X},\boldsymbol{p}) \right) \ln \frac{\Delta^{12}\_{\text{aa}}(\mathcal{X},\boldsymbol{p})}{\Delta^{21}\_{\text{aa}}(\mathcal{X},\boldsymbol{p})}, \end{split} \tag{77}$$

with the definitions of entropy current *<sup>s</sup>μ*matter,*αα* for (*αα*) components,

$$s\_{\text{matter},\mu\alpha}^{\mu} \equiv \int\_{p} \left[ \left( \delta\_{0}^{\mu} + \frac{\delta\_{i}^{\mu} \mathbf{p}^{i}}{m} - \frac{\partial \text{Re}\Sigma\_{\text{ax},R}}{\partial p\_{\mu}} \right) \frac{\rho\_{\text{ax}}}{i} + \frac{\partial \text{Re}\Delta\_{\text{ax},R}}{\partial p\_{\mu}} \frac{\Sigma\_{\text{ax},\rho}}{i} \right. $$

$$+ \frac{\partial \text{Re} V\_{\text{ax},R}}{\partial p\_{\mu}} \frac{\Delta\_{\text{\text{\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\textdegree\text\text$$

In this derivation, we have used the same way as that in open systems in [59]. We can also derive the following equations for the Kadanoff–Baym equations for photons with the Kadanoff–Baym Ansatz *<sup>D</sup>*21*T* = *ρTi* (1 + *fT*), *<sup>D</sup>*12*T* = *ρTi fT*, *<sup>D</sup>*21*L* = *ρLi* (1 + *fL*) and *<sup>D</sup>*12*L* = *ρLi fL* with distribution functions *fT* and *fL* and spectral functions *ρT* and *ρ<sup>L</sup>*,

$$\begin{split} \partial\_{\boldsymbol{\mu}}s^{\mu}\_{\text{photon}} &= -\frac{1}{2} \int\_{\boldsymbol{k}} \left[ \Pi^{21}\_{\text{T}}(\mathbf{X},\boldsymbol{k}) D^{12}\_{\text{T}}(\mathbf{X},\boldsymbol{k}) - \Pi^{12}\_{\text{T}}(\mathbf{X},\boldsymbol{k}) D^{21}\_{\text{T}}(\mathbf{X},\boldsymbol{k}) \right] \ln \frac{D^{12}\_{\text{T}}(\mathbf{X},\boldsymbol{k})}{D^{21}\_{\text{T}}(\mathbf{X},\boldsymbol{k})} \\ &- \frac{1}{2} \int\_{\boldsymbol{k}} \left[ \Pi^{21}\_{\text{L}}(\mathbf{X},\boldsymbol{k}) D^{12}\_{\text{L}}(\mathbf{X},\boldsymbol{k}) - \Pi^{12}\_{\text{L}}(\mathbf{X},\boldsymbol{k}) D^{21}\_{\text{L}}(\mathbf{X},\boldsymbol{k}) \right] \ln \frac{D^{12}\_{\text{L}}(\mathbf{X},\boldsymbol{k})}{D^{21}\_{\text{L}}(\mathbf{X},\boldsymbol{k})}, \end{split} \tag{79}$$

with the entropy current for photons,

$$\begin{split} s\_{\text{photon}}^{\mu} & \equiv \quad \int\_{k} \left[ \left( k^{\mu} - \frac{1}{2} \frac{\partial \text{Re}\Pi\_{T,R}}{\partial k\_{\mu}} \right) \frac{D\_{T,\rho}}{i} + \frac{1}{2} \frac{\partial \text{Re}D\_{T,R}}{\partial k\_{\mu}} \frac{\Pi\_{T,\rho}}{i} \right] \sigma[f\_{T}] \\ & \quad + \int\_{k} \left[ \left( k^{\mu} - \frac{1}{2} \frac{\partial \text{Re}\Pi\_{L,R}}{\partial k\_{\mu}} \right) \frac{D\_{L,\rho}}{i} + \frac{1}{2} \frac{\partial \text{Re}D\_{L,R}}{\partial k\_{\mu}} \frac{\Pi\_{L,\rho}}{i} \right] \sigma[f\_{L}]. \end{split} \tag{80}$$

As a result, the total entropy current *sμ* = *<sup>s</sup>μ*matter,00 + ∑*α <sup>s</sup>μ*matter,*αα* + *sμ*photon satisfies,

*∂μs<sup>μ</sup>* = (*ede*)<sup>2</sup> *p*,*k k*<sup>0</sup><sup>2</sup> ∑*α* Δ21*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ1200(*p*)*D*21*T* (*k*) − <sup>Δ</sup>12*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ2100(*p*)*D*12*T* (*k*) × ln <sup>Δ</sup>21*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ1200(*p*)*D*21*T* (*k*) <sup>Δ</sup>12*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ2100(*p*)*D*12*T* (*k*) +(*ede*)<sup>2</sup> *p*,*k k*<sup>0</sup><sup>2</sup> ∑*α* Δ21*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ1200(*p*)*D*21*L* (*k*) − <sup>Δ</sup>12*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ2100(*p*)*D*12*L* (*k*) × ln <sup>Δ</sup>21*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ1200(*p*)*D*21*L* (*k*) <sup>Δ</sup>12*αα*(*<sup>p</sup>* − *<sup>k</sup>*)Δ2100(*p*)*D*12*L* (*k*) +(*ede*)<sup>2</sup>(**E**(*X*))<sup>2</sup> ∑ *α p* Δ21*g*,*αα*(*<sup>p</sup>* + *α∂ζ*)Δ1200(*p*) − <sup>Δ</sup>12*g*,*αα*(*<sup>p</sup>* + *α∂ζ*)Δ<sup>2100</sup> × ln <sup>Δ</sup>21*g*,*αα*(*<sup>p</sup>* + *α∂ζ*)Δ1200(*p*) <sup>Δ</sup>12*g*,*αα*(*<sup>p</sup>* + *α∂ζ*)Δ2100(*p*) ≥ 0, (81)

where we have used the inequality (*x* − *y*)ln *xy* ≥ 0 for real variables *x* and *y* with *x* > 0 and *y* > 0. The equality is satisfied in *f*00 = *fαα* = *fT* = *fL* = 1 *ep*0/*<sup>T</sup>*−<sup>1</sup> . Here we have used <sup>Δ</sup>21*αα* <sup>Δ</sup>12*αα* ∼ <sup>Δ</sup>21*g*,*αα* <sup>Δ</sup>12*g*,*αα* with *γg*,*αα* ∼ *fαα* in first order in the gradient expansion. We have shown the H-theorem in the LO approximation in the coupling expansion and in the first order approximation in the gradient expansion. There is no violation in the second law in thermodynamics in the dynamics.

#### **4. Time Evolution Equations in Spatially Homogeneous Systems and Conserved Energy**

In this section, we write time evolution equations in spatially homogeneous systems and show a concrete form of the conserved energy density.

It is convenient to introduce the statistical functions *F*00 = <sup>Δ</sup>2100+<sup>Δ</sup>1200 2 , *Fαα* = <sup>Δ</sup>21*αα*+Δ12*αα* 2 , *FT* = *<sup>D</sup>*21*T* <sup>+</sup>*D*12*T* 2 , *FL* = *<sup>D</sup>*21*L* <sup>+</sup>*D*12*L* 2 , which represent the information of how many particles are occupied in (*p*0, **p**) (particle distributions) and statistical parts, *Uαα*,*<sup>F</sup>* = *<sup>U</sup>*21*αα*+*U*12*αα* 2 , *Vαα*,*<sup>F</sup>* = *V*<sup>21</sup> *αα*+*V*<sup>12</sup> *αα* 2 , <sup>Δ</sup>*g*,*αα*,*<sup>F</sup>* = <sup>Δ</sup>21*g*,*αα*+Δ12*g*,*αα* 2 , Σ00,*F* = <sup>Σ</sup>2100+<sup>Σ</sup>1200 2 , Σ*αα*,*<sup>F</sup>* = <sup>Σ</sup>21*αα*+Σ12*αα* 2 , Π*T*,*<sup>F</sup>* = <sup>Π</sup>21*T* <sup>+</sup>Π12*T* 2 and Π*L*,*<sup>F</sup>* = <sup>Π</sup>21*L* <sup>+</sup>Π12*L* 2 . The variables of these functions are (*X*0, *p*0, **p**) with the center-of-mass coordinate *X*<sup>0</sup> = *x*0+*y*<sup>0</sup> 2 and *p* given by the Fourier transformation with the relative coordinate *x* − *y* in variables (*<sup>x</sup>*, *y*) in Green's functions and self-energy in Section 2. The statistical functions and parts are real at any time when we start with real statistical functions at initial time. The spectral functions are given by taking the difference of (2, 1) and (1, 2) components multiplied by *i*, namely *ρ*00 = *i*(Δ2100 − <sup>Δ</sup>1200). They represent the information of which states can be occupied by particles in (*p*0, **p**) (dispersion relations). The spectral parts in self-energy are given by taking the difference of (2, 1) and (1, 2) components multiplied by *i* (and written by the subscript *ρ*), namely <sup>Δ</sup>*g*,*αα*,*<sup>ρ</sup>* = *<sup>i</sup>*(Δ21*g*,*αα* − <sup>Δ</sup>12*g*,*αα*), <sup>Σ</sup>00,*ρ* = *i*(Σ2100 − Σ1200) and so on. The spectral functions and parts are pure imaginary at any time when we start with pure imaginary spectral functions at initial time. We can use the real statistical parts labeled by the subscripts *F* and the pure imaginary spectral parts labeled by the subscript *ρ* in self-energy in the time evolution. We use the subscript '*R*', '*F*' and '*ρ*' to represent the retarded, statistical and spectral parts in self-energy, respectively.

The Kadanoff–Baym equation for the statistical and spectral functions are given by,

$$\begin{aligned} \left\{ p^0 - \frac{\mathbf{p}^2}{2m} - \text{Re}\Sigma\_{00,\mathbb{R}} + \sum\_{a=-1,1} \text{Re}\mathcal{U}\_{aa,\mathbb{R},\prime} F\_{00} \right\} &+ \left\{ \text{Re}\Delta\_{00,\mathbb{R}}\Sigma\_{00,\mathbb{F}} - \sum\_a \mathcal{U}\_{aa,\mathbb{F}} \right\} \\ &= \frac{1}{i} \left( F\_{00}\Sigma\_{00,\mathbb{\rho}} - \rho\_{00}\Sigma\_{00,\mathbb{F}} \right) - \frac{1}{i} \sum\_a \left( F\_{00}\mathcal{U}\_{aa,\mathbb{\rho}} - \rho\_{00}\mathcal{U}\_{aa,\mathbb{F}} \right), \end{aligned} \tag{82}$$

$$\left\{ p^0 - \frac{\mathbf{p}^2}{2m} - \mathrm{Re}\Sigma\_{00,\mathbb{R}} + \sum\_{a=-1,1} \mathrm{Re}\mathcal{U}\_{aa,\mathbb{R},\rho} \rho\_{00} \right\} + \left\{ \mathrm{Re}\Delta\_{00,\mathbb{R},\mathbb{Z}}\Sigma\_{00,\rho} - \sum\_a \mathcal{U}\_{aa,\rho} \right\} = 0,\tag{83}$$

$$\begin{split} \left\{ p^{0} - \frac{\mathbf{p}^{2}}{2m} - \frac{1}{2I} - \text{Re}\Sigma\_{\text{aa},\text{R},\text{F}}\Sigma\_{\text{aa}} \right\} &+ \left\{ \text{Re}\Lambda\_{\text{aa},\text{R},\text{F}}\Sigma\_{\text{aa},\text{F}} \right\} + \left\{ \text{Re}V\_{\text{aa},\text{R},\text{A}\_{\text{g},\text{a},\text{F}}} \right\} - \left\{ \text{Re}\Lambda\_{\text{g},\text{a},\text{R},\text{F}}V\_{\text{aa},\text{F}} \right\} \\ &= \frac{1}{\bar{l}} \left( \mathcal{F}\_{\text{aa}}\Sigma\_{\text{aa},\text{g}} - \rho\_{\text{aa}}\Sigma\_{\text{aa},\text{F}} \right) - \frac{1}{\bar{l}} \left( \Delta\_{\text{g},\text{a},\text{F}}V\_{\text{aa},\text{g}} - \Delta\_{\text{g},\text{a},\text{g}}V\_{\text{aa},\text{F}} \right), \end{split} \tag{84}$$

$$\begin{aligned} \left\{ p^0 - \frac{\mathbf{p}^2}{2m} - \frac{1}{2I} - \text{Re}\Sigma\_{aa,R}, \rho\_{aa} \right\} &+ \left\{ \text{Re}\Delta\_{aa,R}, \Sigma\_{aa,\rho} \right\} \\ + \left\{ \text{Re}V\_{aa,R}, \Delta\_{\mathfrak{F},aa,\rho} \right\} - \left\{ \text{Re}\Delta\_{\mathfrak{F},aa,R}, V\_{aa,\rho} \right\} &= 0, \end{aligned} \tag{85}$$

$$\left\{p^{0} - \frac{\mathbf{p}^{2}}{2m} - \frac{1}{2I} - \text{Re}\Sigma\_{aa,R'}\Delta\_{\text{g.a.},F}\right\} + \left\{\text{Re}\Delta\_{\text{g.a.},R,\prime}\Sigma\_{aa,F}\right\}$$

$$= \frac{1}{i}\left(\Delta\_{\text{g.a.},F}\Sigma\_{aa,\rho} - \Delta\_{\text{g.a.},\rho}\Sigma\_{aa,F}\right),\tag{86}$$

$$\left\{p^0 - \frac{\mathbf{p}^2}{2m} - \frac{1}{2I} - \text{Re}\Sigma\_{\text{aa},R\prime}\Delta\_{\text{\textdegree},\text{aa},\rho}\right\} + \left\{\text{Re}\Delta\_{\text{\textdegree},\text{aa},R\prime}\Sigma\_{\text{\textdegree},\text{\textdegree},\rho}\right\} = 0,\tag{87}$$

$$\begin{Bmatrix} p^2 - \text{Re}\Pi\_{\mathbb{R},T} F\_T \\ \text{~} \end{Bmatrix} + \begin{Bmatrix} \text{Re} D\_{\mathbb{R},T} \Pi\_{\mathbb{F},T} \\ \text{~} \end{Bmatrix} = \begin{array}{c} \frac{1}{i} \left( F\_T \Pi\_{\mathbb{R},T} - \rho\_T \Pi\_{\mathbb{F},T} \right) \\ \text{'} \end{array} \tag{88}$$

$$\left\{p^2 - \text{Re}\Pi\_{\text{R},T}\rho\_T\right\} + \left\{\text{Re}D\_{\text{R},T}\Pi\_{\rho,T}\right\} \quad = \quad 0,\tag{89}$$

and longitudinal parts given by changing the label *T* to *L* in Equations (88) and (89). We can use Equation (69) in the previous section to derive Equations (82) and (83), for example.

We can write,

$$\mathcal{U}\_{\rm an,F}(X,p) = (\varepsilon d\_{\varepsilon})^2 \mathbf{E}(X)^2 \Delta\_{\mathfrak{g},\rm an,F}(p+a\mathfrak{d}\zeta), \quad \mathcal{U}\_{\rm an,\rho}(X,p) = (\varepsilon d\_{\varepsilon})^2 \mathbf{E}(X)^2 \Delta\_{\mathfrak{g},\rm an,\rho}(p+a\mathfrak{d}\zeta), \tag{90}$$

$$V\_{\rm an,F}(X,p) = (\varepsilon d\_\varepsilon)^2 \mathbf{E}(X)^2 F\_{00}(p - a\partial\_\circ^\circ), \qquad V\_{\rm an,p}(X,p) = (\varepsilon d\_\varepsilon)^2 \mathbf{E}(X)^2 \rho\_{00}(p - a\partial\_\circ^\circ). \tag{91}$$

In case we start with initial condition *<sup>E</sup>*2(*X*<sup>0</sup> = 0) = 0, *<sup>∂</sup>*0*E*2(*X*<sup>0</sup> = 0) = 0 and symmetric Green's functions for *α* → −*α* in spatially homogeneous systems, we can use *∂ζ* = 0 in the above equations at any time. We can write the self-energy as,

$$\Sigma\_{00,\mathbb{F}}(p) \quad = \quad - (\epsilon d\_{\ell})^2 \sum\_{a=-1,1} \int\_{\mathbf{k}} \left(\mathbf{k}^0\right)^2 \left[F\_{\text{ad}}(p-k)(F\_T(\mathbf{k}) + F\_L(\mathbf{k})) + \frac{1}{4} \frac{\rho\_{\text{ad}}(p-k)}{i} \frac{\rho\_T(\mathbf{k}) + \rho\_L(\mathbf{k})}{i}\right], \quad \text{(92)}$$

$$\Sigma \alpha\_{\vartheta}(p) = -(\epsilon d\_{\ell})^2 \sum\_{a=-1,1} \int\_{\mathbf{k}} \left( \mathbf{k}^0 \right)^2 \left[ F\_{\mathbf{k}\mathbf{k}}(p-k)(\rho \mathbf{r}(k) + \rho \boldsymbol{\upmu}(k)) + \rho\_{\mathbf{k}\mathbf{k}}(p-k)(\mathbf{F}\mathbf{r}(k) + \mathbf{F}\boldsymbol{\upmu}(k)) \right],\tag{93}$$

$$\Sigma\_{\text{aux},\mathbb{F}}(p) = -(\varepsilon d\_{\varepsilon})^2 \int\_{k} \left(k^0\right)^2 \left[F\_{00}(p-k)(F\_T(k) + F\_L(k)) + \frac{1}{4} \frac{\rho\_{00}(p-k)}{i} \frac{\rho\_T(k) + \rho\_L(k)}{i}\right],\tag{94}$$

$$\Sigma\_{\text{an},\rho}(p) \quad = -(\epsilon d\_{\text{\textell}})^2 \int\_{\mathcal{k}} \left(k^0\right)^2 \left[\text{Fo}(p-k)(\rho\_T(k) + \rho\_L(k)) + \rho\_{\text{00}}(p-k)(F\_T(k) + F\_L(k))\right],\tag{95}$$

$$\Pi\_{T,\mathcal{F}}(k) = \Pi\_{l,\mathcal{F}}(k) \ = \ - (cd\_{\epsilon})^2 \left(k^0\right)^2 \sum\_{a=-1,1} \int\_p \left[F\_{aa}(k+p)F\_{00}(p) - \frac{1}{4} \frac{\rho\_{aa}(k+p)}{i} \frac{\rho\_{00}(p)}{i}\right]$$

$$+ F\_{00}(k+p)F\_{aa}(p) - \frac{1}{4} \frac{\rho\_{00}(k+p)}{i} \frac{\rho\_{aa}(p)}{i}\Big|\_{}.\tag{96}$$

$$\Pi\_{T, \theta}(k) = \Pi\_{L, \varphi}(k) \ = \ -(ed\_{\epsilon})^2 \left(k^0\right)^2 \sum\_{a=-1,1} \int\_p \left[\rho\_{\rm{ax}}(k+p)F\_{00}(p) - F\_{\rm{ax}}(k+p)\rho\_{00}(p)\right] \tag{97}$$

$$+ \rho \mathbf{u}(k+p)F\_{\rm{ax}}(p) - F\_{\rm{0}}(k+p)\rho\_{\rm{ax}}(p)\right],\tag{98}$$

where we have omitted the label of the center-of-mass cordinate *X* in Green's functions and self-energy. We find that the <sup>Π</sup>*T*,*<sup>F</sup>*(*k*) = <sup>Π</sup>*L*,*<sup>F</sup>*(*k*) are symmetric (Π*T*,*<sup>F</sup>*(−*k*) = <sup>Π</sup>*T*,*<sup>F</sup>*(*k*)) under *k* → −*k* and that <sup>Π</sup>*T*,*<sup>ρ</sup>* = <sup>Π</sup>*L*,*<sup>ρ</sup>* are anti-symmetric (Π*T*,*<sup>ρ</sup>*(−*k*) = <sup>−</sup>Π*T*,*<sup>ρ</sup>*(*k*)) under *k* → <sup>−</sup>*k*, for any Green's functions for dipole fields. When we prepare initial conditions with symmetric *FT*,*<sup>L</sup>* and anti-symmetric *ρ<sup>T</sup>*,*<sup>L</sup>* for photons, we can derive symmetric *FT*,*<sup>L</sup>* and anti-symmetric *ρ<sup>T</sup>*,*<sup>L</sup>* at any time. In addition, since Π(*k*)'s are proportional to (*k*<sup>0</sup>)2, there is no mass gap for incoherent photons for the leading-order self-energy in the coupling expansion. The velocity of gapless modes of incoherent photons will decrease when we increase the density of dipoles in this theory.

Finally, we show the energy density *E*tot. In the spatially homogeneous system in the 2 + 1 dimensions, we can derive *∂E*tot *∂X*<sup>0</sup> = 0 with the energy density given by,

$$\begin{split} \text{E}\_{\text{tot}} & \equiv \quad \frac{1}{2I} \sum\_{a=-1,1} \overline{\psi}\_{a}^{\*} \overline{\psi}\_{a} + \frac{1}{2} \left( \partial\_{0} A\_{i} \right)^{2} + \int\_{p} p^{0} \left( F\_{00} + \sum\_{a=-1,1} F\_{aa} \right) + \frac{1}{2} \int\_{p} \left( p^{0} \right)^{2} \left( F\_{T} + F\_{L} \right) \\ & \quad + 2 \langle ed\_{\epsilon} \rangle^{2} \text{E}^{2} \sum\_{a=-1,1} \int\_{p} \left( F\_{00} (p) \text{Re} \Delta\_{\mathbb{S},aa,R} (p + a \partial\_{\epsilon}^{\prime}) + \text{Re} \Delta\_{\mathbb{t}0,R} (p) \Delta\_{\mathbb{S},aa,F} (p + a \partial\_{\epsilon}^{\prime}) \right) \\ & \quad - \int\_{p} \left( \text{Re} \Sigma\_{00,R} F\_{00} + \text{Re} \Lambda\_{00,R} \Sigma\_{00,F} \right) - \sum\_{a=-1,1} \int\_{p} \left( \text{Re} \Sigma\_{aa,R} F\_{\text{an}} + \text{Re} \Delta\_{\text{aa},R} \Sigma\_{aa,F} \right) \\ & \quad - \frac{1}{2} \int\_{p} \left( \text{Re} \Pi\_{R,T} F\_{T} + \text{Re} D\_{R,T} \Pi\_{F,T} + \text{Re} \Pi\_{R,L} F\_{L} + \text{Re} D\_{R,L} \Pi\_{F,L} \right), \end{split} \tag{98}$$

where we have used the KB equations in this section, the Klein–Gordon Equation (20) and the Schödinger-like Equations (23) and (24) in Section 2. The first term represents the contribution of nonzero angular momenta for coherent dipole fields. The second term represents the contribution by electric fields *Ei* = *∂*0*Ai*. The third and the fourth terms represent the contribution by quantum fluctuations for dipoles and photons, respectively. When the temperature is nonzero *T* = 0 at equilibrium states and the spectral width in the spectral functions is small enough, statistical functions which are proportional to the Bose–Einstein distributions 1 *ep*0/*<sup>T</sup>*−<sup>1</sup> give temperature-dependent terms *mT*<sup>2</sup> for dipole fields and ∝ *T*<sup>3</sup> for photon fields in 2 + 1 dimensions. The fifth term represents the potential energy in processes in Figure 2a. The sixth, seventh and eighth terms represent the potential energy in processes in Figure 2b. The coefficients in the sixth and seventh terms are not 13 but 1. While the factor 1 might look like a contradiction with the preceding research in [68,69] which sugges<sup>t</sup> that the factor 13 appears in the interaction with 3-point-vertex, the factor 1 appears due to time derivative (*∂*<sup>0</sup>)<sup>2</sup> in self-energy for dipole fields and photon fields.

#### **5. Dynamics of Coherent Fields**

In this section, we show that our Lagrangian describes the super-radiance phenomena in time evolution equations of coherent fields. We shall assume that all the coherent fields are independent of *x*1 (dependent on *x*0 and *x*2). We also assume the symmetry for *α* = −1 and *α* = 1, namely *ψ*¯(∗) 1 = *ψ*¯(∗) −1 , Δ01 = Δ0−1, and Δ10 = Δ−10. We set initial conditions *E*2 = 0 and *∂*0*E*2 = 0 at *x*0 = 0.

We define *Z* ≡ 2|*ψ*¯1| 2 − | *ψ*¯0| 2. It is possible to derive the following equations from time evolution Equations (20), (23) and (24) with their complex conjugates for background coherent fields in Section 2.

$$\begin{array}{rcl}\partial\_0 Z &=& \mathbf{i}4cd\_c E\_1 \left(\ddot{\psi}\_1^\* \ddot{\psi}\_0 - \ddot{\psi}\_0^\* \ddot{\psi}\_1\right), \end{array} \tag{99}$$

$$\partial\_0 \left( \bar{\psi}\_1^\* \bar{\psi}\_0 \right) \quad = \quad \frac{i}{2I} \bar{\psi}\_1^\* \psi\_0 + i \epsilon d\_\epsilon E\_1 Z \tag{100}$$

$$\left[\left(\partial\_{0}\right)^{2} - \left(\partial\_{2}\right)^{2}\right]E\_{1} = -2\varepsilon d\_{\varepsilon}(\partial\_{0})^{2}\left[\bar{\psi}\_{1}^{\*}\bar{\psi}\_{0} + \bar{\psi}\_{0}^{\*}\bar{\psi}\_{1} + \Delta\_{01}(\mathbf{x},\mathbf{x}) + \Delta\_{10}(\mathbf{x},\mathbf{x})\right].\tag{101}$$

We have used moderately varying spatial dependence |∇<sup>2</sup> *i ψ*¯−1,0,1/*m*||*∂*0*ψ*¯−1,0,1|. We derive aspects of the super-radiance and the Higgs mechanism in the above three equations.
