Path Integral based Non-equilibrium Quantum Field Theory of Non-relativistic Pairs inside an Environment

We derive differential equations from path integral based non-equilibrium quantum ﬁeld theory, that cover the dynamics and spectrum of non-relativistic two-body ﬁelds for any environment. For concreteness of the two-body ﬁelds, we choose the full potential non-relativistic Quantum Electrodynamics Lagrangian in this work. After closing the correlation function hierarchy of these differential equations and performing consistency checks with previous literature under certain limits, we demonstrate the range of physics applications. This includes Cosmology such as Dark Matter in the primordial plasma, Quarkonia Physics inside a quark-gluon plasma, and Condensed and strongly Correlated Matter Physics such as Bose-Einstein condensation or Superconductivity. Since we always had to take limits or approximations of our equations in order to recover those known cases, our equations could contain new phenomena. In particular they are based on non-equilibrium Green’s function that can deal with non-hermite potentials as well as dynamical formation of different extreme phases. We propose a scheme for other Lagrangian based theories or higher N-body states such as molecules to derive analogous equations.

one-point functions in Section 5, needed to discuss extreme phenomena like the spontaneous symmetry breaking. In particular system transitions from local to global phases in the presence of critical electric and magnetic field values.
After all these deriving steps, in Section 6, we show to certain extend the power of some of our equations by explicitly demonstrating in a simple spin field model, that they can even cope with the dynamical description of phase transitions. In Section 7, we give a comprehensive consistency check of all our results (six coupled differential equations) by applying specific limits or neglecting contributions. For example, in the absolute vacuum and classic, static electric and magnetic field limit, where all the dynamical equations become trivial, quantum mechanics from the remaining spectral equation emerges. I.e. fine, and hyperfine corrections to two-body pair energy levels pop out directly from the theory in these limits and even higher order contributions. Furthermore, in the dilute limit of one of our dynamical equation and vacuum assumption of the spectral function, we prove that this dynamical equation reduces to the open quantum system treatment (QED analog) of Quarkonia inside the quark-gluon plasma. We also comment on Bose-Einstein condensation for general occupation number and also argue that our one-point function equations are consistent with literature on Superconductivity under some limits. The work is discussed in Section 8 before concluded in Section 9.

Potential non-relativistic effective field theory
In this work, we choose for simplicity Quantum Electrodynamics (P.A.M. Dirac 1928), i.e., a heavy Dirac fermion, χ, charged under the gauge group U(1): where L env [A] contains an arbitrary environment, such as a plasma with a temperature T . This plasma can, e.g., contain other ultra-relativistic fermions interacting with the gauge field A. The precise content of the environment is kept open and the form of our results are formulated independently. We are interested in the non-equilibrium dynamics of non-relativistic χ fields in this environment. To simplify Eq. (1) to the relevant part, we adopt potential non-relativistic effective field theory (pNREFT) [22][23][24] . This is a designated 'effective theory for non-relativistic two-body states'. The steps from the above relativistic theory to this potential non-relativistic theory are the following. First, one takes the non-relativistic expansion of χ fields (technically, 'integrating out the scale m χ '), thereby splitting them into particle and anti-particle contributions, and derives a 2 by 2 self-scattering term which sets in the static limit the effective potential between pairs. This leads to non-relativistic effective field theory 25,26 . Second, one projects this non-relativistic theory to the two-particle sub-space and performs a multipole expansion of the gauge field ('integrating out m χ v for continous and Bohr momenta for various discrete states'). Then, the relevant part of the Lagrangian density for the χ-χ particle and anti-particle pairs and its interaction through the gauge field at leading order in this multipole expansion are given by [22][23][24] : with the trace over spin indices, O being the two-body field operator of (χ) particle and anti-particle pair with reduced mass m and total mass M, depending on the center-of-mass coordinate x and relative distance r, and For O being the two-body field operator of a particle anti-particle pair as we consider by default, the effective static potential at the leading order in vacuum is given by V eff (r) = − α r , which makes our convention clear. The sign of this potential, as well as other higher order terms in the mass expansion as contained in [· · · ] terms in this h can be different for a particle particle, or an anti-particle anti-particle pair. For example, only the particle anti-particle pair would have an additional complex absorption coefficient, representing the annihilation (or decay) of the pair state. This term was already taken into account in Ref. 1 , and we can neglect it here as the leading order approximation. Note also that we did not write out explicitly these various other higher order contributions in the mass expansion, which can however be simply adapted.
The effective potential, V eff , at finite temperature can be complex in general 8 . It can be expressed as a general gauge field correlator, e.g., Ref. 1 . It is rotationally invariant, and therefore conserves angular momentum. It can however change the two-body momentum at such higher order, which can mix the scattering with the bound state once such interactions with the environment are efficient, described as an imaginary thermal width 8 . Notice that another reason why we have neglected other [...] terms is that they need to be expressed in terms of gauge field correlator to capture finite temperature effects, which is beyond the scope of this work.
Non-angular momentum conserving transitions among two-body fields are described via the electric dipole moment operator, r · gE, in the above Lagrangian. This interaction term describes two-body transitions through spin-conserving dipole electric moment emission and absorption, such as bound-state formation or dissociation. Spin and angular momentum two-body field transitions take place via the magnetic dipole moment operator, µ · gB, where µ is a three vector containing the orbital angular momentum operators (acting and depending on r only) and spin operators (Pauli matrices). This potential non-relativistic theory is valid for v rel ≪ 1, α = g 2 /(4π) ≪ 1 and T ≪ m for a plasma environment. Thus, it is a general non-relativistic theory of pairs. We leave the content of the environment as open in the course of our following derivation of the dynamics and spectrum of this theory.

Correlation function equations of motion
Let us start by defining a two-body field two-point correlation function, with four spatial arguments, on a general time contour: where T C ... ≡ Tr[ρT C (...)], ρ is the density matrix of the full system including also the environment, and the time contour with time ordering symbol, T C , can be any, e.g., the Keldysh-Schwinger contour 27, 28 (for comprehensive introduction into this non-equilibrium quantum field theory, see, e.g., Ref. 1 ). From the path integral (Feynman, Richard Phillips, 1948), we derive the equation of motions (EoMs) for this two-point function with the potential non-relativistic Lagrangian in Eq. (2), by using the measure invariance principle under infinitesimal field shifts: In this notation, the prime over h in the second line indicates that the relative position is r ′ . We shall make the spin indices implicit for notational simplicity from now in deriving steps, but include spin summation explicitly in some relevant main results. The three point functions with the electric field, as occurring on the right hand side, obey the EoMs: The latter two equations can be rewritten into integral form as Repeating the same for three point functions with the magnetic field, we obtain: Figure 1. Figure shows a self-energy diagram for the two-body field as a consequence of the leading order approximation of the four point correlators in Eq. (13) and Eq. (14). Black blobs stand for interacting electric, magnetic and two-body two-point correlator, respectively. The interaction with the environment is included in these. Gray blobs stand for the electric and magnetic dipole vertices, respectively.
Plugging these expression for the three point functions with magnetic and electric field into Eq. (5) and Eq. (6), leads to: These are the most general out-off equilibrium equations that we can write down in the general time formalism starting from the introduced pNR effective Lagrangian. Note that they apply for any environment. The implicit dependence on the environment is contained in the electric and magnetic fields, as well as in h. These equation represent the full path integral in terms of its field moments (correlation function). Thus, the trace and the state density ρ in the correlator definition should not only be read in the traditional quantum mechanical sense. Here, we use it as just a notation, representing the second moment in fields of Feynmans' path integral. In this sense, we have not used any quantum mechanical state or quantized operators living on Hilbert spaces. We work with classical fields (and Grassmann numbers or Grassmann valued fields), as present in the path integral.
Applying this to Eq. (13) and Eq. (14), gives: Here, we dressed the free propagators to account for some higher order terms that are neglected otherwise from the truncation. In relativistic QED, similar equations would occur. In this case, the dressing of the correlators allows to include the contribution of Compton scattering. This is why we believe on physical arguments, that free correlators, after truncating, should be dressed in order to contain some of the higher order neglected terms. It requires a rigorous mathematical proof, which is beyond the scope of this work. The truncation and dressing leads to the two-body field self-energy diagram as shown in Fig. 1. Note the dressed two-body propagator, that accounts, e.g., for several gauge field emissions or absorption.
In the next two subsection, we adopt, for concreteness, the Keldysh-Schwinger time contour and neglect memory effects for simplicity. These can be included by extending the Keldysh-Schwinger time contour, which is relevant in some Correlated Physics application like the Many Body Localization (MBL) theory, where due to memory effects the closed system without environment eventually never reaches an equilibrium state.
In the Keldysh-Schwinger formalism, each of these two equations turn into a 2 by 2 matrix, where the components are coupled and each component has spin indices, such that, when including spin indices explicitly, the overall equation is a coupled tensor equation in general. We will look at specific components (or combinations of components) of this matrix, and look into them individually. Concretely, we are interested in the statistical component G +− and the retarded equation G R = G ++ − G +− . While the statistical can be roughly speaking seen as a quantum Boltzmann equation, the imaginary part of the retarded correlator is related to the spectrum of the theory. Both are coupled in general as it may have been expected from the traditional Kadanoff-Bayam equations. We shall explicitly see that this somewhat similarly applies to our four space arguments two-time correlation functions considered in this work.

Dynamical equation of statistical component
In the following, we will assume that the system is sufficiently close to a homogeneous and isotropic state, such that we can simplify these derived equations. Let us introduce Wigner coordinates for the Center-of-Mass (CoM) coordinates only: In the following we shall add and subtract Eq. (16) from Eq. (15). To make the derivation steps simpler to follow, let us write down added (A) and subtracted (S) kinetic terms in terms of these introduced Wigner coordinates:

5/19
Using these, Wigner transforming, Fourier transforming, gradient expanding the correlators CoM coordinates, and subtracting the statistical component +− of Eq. (16) from Eq. (15), results in: where we introduced the two-body field self-energy: Σ E i j (x, y; r,r) ≡ g 2 G(x, y; r,r) T C E i (x)E j (y) for the electric component and Σ B i j (x, y; r,r) ≡ g 2 G(x, y; r,r) T C B i (x)B j (y) for the magnetic part. Notice that the potential on the l.h.s. as well as the Laplacian for relative coordinates have exactly canceled. We applied here the limit y → x and took a spin-unpolarized medium. This lead to the fact that all one-point function contributions of the electric and magnetic field canceled identically and that we can write the above collision term in terms of only one common self-energy. For a spin polarized system, the more general structure of the above collision term would be: and we can not choose the applied combination of spin indices that results to the above simplified collision term. The more general case may be relevant to discuss, e.g., the spin-flip temperature evolution at the Cosmic Dawn as discussed later. For the Dark Matter case, this assumption of an unpolarized medium is commonly adopted throughout the thermal history.
For cosmological applications the gradient expansion in the spatial part of the CoM coordinates is exact due to the cosmological principle (spatial homogeneity and isotropy). In this case, G is also independent of R and this spatial gradient on the l.h.s. vanishes. Note that instead for cosmological applications there is always a Hubble expansion term occuring on the l.h.s. For an estimate of the error introduced from the time component gradient expansion in such scenarios, see, e.g., Ref. 29 .
To proceed, let us assume that the environment can be described by a Kadanoff-Bayam Ansatz. Consequently, we can relate some of the electric field correlator components to the spectral function: Then, with these relations, the self-energy can be written as Plugging this into the dynamical equation of the statistical component, the equations are closed under the statistical and spectral correlation functions: To obtain the total two-body particle number, we need to take the limit r ′ → r, and integrate over P, relative position and summing over spins: n 2 ≡ P r G +− (T, R, P; r, r). This leads indeed to vanishing right hand side, meaning that the total two-body number n 2 is separately conserved under electric and magnetic dipole transitions in our equations. This is expected, since the pNR Lagrangian we started from has a global symmetry. If the thermal width is negligible, then this number reduces to: n 2 = ∑ B n B +scattering, meaning that one can distinguish continues from discrete states by separating the energy integration. For a large thermal width, one can however not do this decomposition, and Eq. (33) and Eq. (34) allows to account for such general quasi-particle excitation.
Let us come back to our statistical Eq. (32). We may take an Ansatz for this statistical two-body field two-point correlator and express it in terms of the spectral function as: The spectral function can be obtained from the imaginary part of the retarded components, implying equations are closed in G space (see, next Section).The dispersion relation entering the out-off-equilibrium phase-space distribution, f , with four spin d.o.f., is determined by the spectral function and can be non-trivial if a thermal width is present. While for single elementary particles, this would be exactly the Kadanoff-Bayam Ansatz, we have assumed here for two-body field correlators (composite object, with four spatial coordinate dependence), that the f does not depend on the relative distances r, r ′ .
Since this may be a crucial assumption, let us argue in more detail. In particular, under the vacuum limit of the retarded correlator in Eq. (40), this Ansatz leads to consistent non-equilibrium equations as shown later. At finite temperature (nonzero thermal width), and in kinetic and chemical equilibrium, this function would be an exact Bose-Einstein distribution, f = (e β (M+2P 0 ) − 1) −1 , which is set by the KMS condition. Therefore, also in this case, f does not depend on r, r ′ . Outof-chemical equilibrium but in kinetic and ionization equilibrium, this function would be a Bose-Einstein distribution with a chemical potential, µ, for the total pair number: f = (e β (M+2P 0 −µ) − 1) −1 . Also in this case, f does not depend on r, r ′ as fixed by the KMS relation with finite chemical potential. So the only remaining case, where this function could obtain a relative distance dependence, would be if the potential has strong complex contributions in out-off-ionization equilibrium situations. Currently, we can not exclude this possibility, but shall continue under the assumption that f is independent of relative differences. At least for temperatures up to T ≫ αm χ ≫ m D (where m D is the Debye mass), it still should be a good, if not an exact description. For m D αm χ , the thermal width becomes important and bound states melt 9 .
Inserting this Ansatz in Eq. (35) into the statistical Eq. (32), results in our final simplified form of the statistical function: Spin contraction causes a difference between the electric and magnetic collision term as expected. At the leading order in the G spectral function, it is proportional to spin delta, such that one recovers in the second and third line on the right hand side spin conserving bound-state formation and dissociation with electric dipole emission and absorption, respectively. This can however change at next-to-leading order, where, in addition to this outer E field emission/absorption, a second emission (or absorption) of a B field quanta can flip the spin. The spin flip is present already at the leading order in the magnetic collision term (fourth and fifth line).
This equation is linear in f if the two-body system is sufficiently dilute ( f ≪ 1, as for thermally produced Dark Matter and Quarkonia), so a development of an efficient numerical two-body phase-space density solution seems to be tractable (e.g., a good starting point would be 30 ). To discuss, e.g., Bose-Einstein condensation, however, we can not apply the dilute limit and the Eq. (37) is in general non-linear.
d 3 rd 3r r ir j Σ R,E i j (T, R, P; r,r)G +− (T, R, P;r, r) + G +− (T, R, P; r,r)Σ A,E i j (T, R, P;r, r) Spin indices with the magnetic dipole is here µ i G = (µ i ) oo G osso , hence, only the third Pauli matrix contributes.

Dynamical equation of spectral component
Subtracting the ++ minus +− correlator of Eq. (14) from Eq. (13), leads to the dynamical retarded equation: In the limit r ′ → r and free of interactions, it has plan wave solution.

Constraint equation of spectral component
Adding the ++ minus +− correlator of Eq. (14) and Eq. (13), leads to the constraint retarded equation for the two-body field two-point correlator: The left hand side is a kind of Schrödinger type of operator and the right hand side contains a delta function, characteristic for Green's functions, plus a contribution coming from electric and magnetic dipole interactions with the system, including the environment. Since the effective in-medium potential, V eff , enters directly in this retarded equation, let us be specific. For ultrarelativistic fermions in the plasma environment and dilute two-body sub-system, it is, for a certain temperature range much above the absolute value of the considered bound state binding energy given by 1,8,9 V eff (r) = −αm D − α r e −m D r − iαT φ (m D r).

8/19
The first term is the Salpeter correction, the second term is a Debye-mass screend Coulomb potential and the last term is purely complex, associated to static Coulomb scattering with ultra-relativistic bath particles from the environment with temperature T . The dimensionless function φ is zero for r = 0 and 1 for r → ∞. In between, the function is monotonically increasing. This term is complex since efficient elastic scattering with ultrarelativistic plasma particles mixes the energies of a two-body state with fixed quantum number (E.g., it conserves angular momentum quantum number due to rotational invariance). This uncertainty in the energy coming from outside appears here as a complex thermal width and is also known as Landau damping in the literature, since it originates from space-like exchange of dressed photon, whose Hard-Thermal-Loop self-energy has an imaginary part 31 .
In the limit r → ∞ the imaginary term exactly reduces to twice the thermal width of a single particle. This complex term gives rise to melting effects of two-body bound states, as observed in Quarkonia systems inside a quark-gluon plasma [10][11][12][13][14] .
The other terms on the right hand side in Eq. (40) originate from the ultra-soft electric and magnetic dipole contributions in the pNR Lagrangian and are new to us. They mix different two-body states with different quantum number with each other. This is why correlators, after angular momentum decomposition, will mix with each other through extremely efficient electric and magnetic dipole transitions inside the plasma encoded in that term. In the QCD case, with a color electric dipole operator, we expect that a similar term appears that mixes the singlet with the octet retarded correlators, while V eff can not introduce such mixing. V eff is elastic scattering, while the two new contributions are coming from dissociation and recombination like terms. Through these new terms, a specific bound state feels the presence of all the other states. In vacuum, and in particular at finite temperature.

One-point function equations
The two-body field two-point function equations as derived in the previous section depend on the one-point function of the electric and magnetic field. To close the system self-consistently, let us express these one-point functions in terms of two-body fields, by deriving the equation of motion similarly to the previous section from the path integral with the introduced pNREFT Lagrangian, obtaining: These expressions for the expectation of the electric and magnetic field are fully quantum. We have not used any classical description (Maxwell equations). As before, r and µ should be understood in this notation as usual coordinates. Consequently, one can take them out of the inner thermal expectation value, however, we would like to keep them inside because this makes spin summation more clear. The trace is over spin indices, and the dots stand for all possible other interactions with the environment or externally applied fields. Under the integral, this can be expressed as our G +− (x, x, r, r), which is noting but the pair particle number density once integrated over r. So the meaning of this equation is that the electric and magnetic field on the left hand side, is the sum over all electric and magnetic dipole moments present in the system. Let us also derive the EoM for the two-body field one-point function from the path integral with the introduced pNREFT Lagrangian, resulting in: The EoMs of the two-point function from the previous section do not depend on these one-point functions. However, the one-point functions may be used to characterize the state of a system. The goal now is to express the correlators on the right hand side in terms of two-point functions from the previous section, such that we get closed coupled equations. To do so, we first require O † O to be invariant under infinitesimal field shifts, O † → O † + ε † . This leads to the following EoM in integral form: and we neglected higher order contributions in the coupling and external fields. From this we get the needed correlators, by integrating over relative coordinates as follow: O ts (x, r)B(x) = − d 3r G tsab (x, x; r,r)μ bc g O ca (x,r) + ...    (Fig. 2b). The initial state is described by a plane wave state, where the absolute value is globally set to one, and the phase is locally rotating. At the origin, at x = 0, a small external magnetic field with a narrow Gaussian profile is applied short after the initial time. It is switched on and its value is increased slowly. The system reacts to this perturbation and creates, at some critical value of the applied external field, excess in the absolute value of the two-body field, de-located from the center. After this transition the external field is switched off at t = 150, and the system evolves under free dynamics. The phase shows chaotic behavior in the region where the applied field has passed the critical value.
These are the leading order expressions, that contribute as kind of mean field terms to the dynamics of the two-body one-point function in Eq. (44). In particular, without the dots they can be used to study closed systems. External fields in Eq. (44) can be included by super-positioning the self-induced electric and magnetic field with the external one, where the latter can probably be described purely classically. Finally, we would like at this stage to remind that the two-body one-point function depends on x and r. The two-body field may be written in terms of its absolute value and two local phases as O ∼ ρe −iφ 1 (t,x)−iφ 2 (t,r) . The two phases may be related from U(1) gauge invariance point of view, where it may depend on the considered type of pair (particle-antiparticle, particle-particle, particle-hole,...). The latter determines also the precise sign in the potential as well as other higher order terms in h, relevant to discuss the gauge invariance. While the φ 1 (t, x) phase is the one under which the two-body field transforms under so called ultra-soft gauge transformations, it is the phase φ 2 (t, r) which , e.g., contributes to vanishing B field in Eq.(43) if it turns spontaneously global. Such phenomena can occur in Superconductivity (Meissner effect) as discussed later. The phase φ 1 (t, x) is the only dynamical one remaining if only spin interactions (Pauli matrices in the magnetic moment) are considered. It could be that it becomes also spontaneously broken by critically applied B-fields.

Transition in a spin field toy model
In the previous two sections, we derived our set of six coupled differential equations. Let us now study the time evolution of the two-body one-point function for a simple case, to demonstrate that the equations may be used to explore the transition into different phases of the system (i.e., phase transitions). To make the problem tractable, it is assumed that the system has only spin non-conserving interactions, as contained in the magnetic dipole moment. In particular, focusing on spin non-conserving interactions only, allows us to get rid of the relative position dependence r. This is because, the only remaining operators are the Pauli matrices inside the magnetic dipole moment, which are space-time independent. Furthermore, we shall consider only interactions with a classical, external magnetic field. Consequently, we may write our two-body one-point function, in Eq. (44), as

10/19
To study the relevant issue, we work in 1 + 1 dimension, set µ to a number and there is only one O field with no spin indices. Note that this corresponds to a non-physical set up. Nevertheless, to solve this simple equation numerically and to finally see a transition when the external magnetic field passes a critical value, let us for demonstration choose M = 2 and µ = g = 1. Also, we choose boundary conditions as the free plane wave solution. We initialize the system with the absolute value of O(x) equals to one and the phase is local. The numerical solution of this initial value problem is shown in Figure 2. We control the value of the external magnetic field by hand.
In particular, the spatial dependence of the external field profile is chosen as a narrow gaussian profile with a spatial width of 10 and applied in the center (x = 0). The time dependence of this external field as a multiplicative factor of the spatial dependence has three stages: i) initially zero and switched on slowly ii) at about t ∼ 40 it stays constant with a maximum value of 5 at the center until t ∼ 150 and iii) afterwards switched off exponentially fast.
During i), the system does almost not recognize the externally applied field until it passes a critical value of 1/4. At this critical point, t ∼ 20, the absolute value somewhat collapses at a particular point that is not in the center, although the applied field is symmetric. Also the phase shows a transition into a chaotic state when the critical field value is passed. The chaotic phase region does not broaden during ii) and returns to certain extend into a local rotating object during stage iii).

Some Physics applications
In the previous sections, we derived our differential equations based on non-equilibrium quantum field theory and the pathintegral formalism. These in total six coupled differential equations are our general theory of non-relativistic pairs in an environment. They could be used to describe the dynamical evolution into extreme states of the system, as well as in the static case they can be used to show the existence of such states. Now, let us apply these formulas to some specific and realistic cases. At the same time, major parts of this section should be read as a comprehensive consistency check. Concretely, certain limits and simplifications are applied to our equations in order to recover known Physics.

Quantum Mechanics of a pair
Under the limit of absolute vacuum and time independent fields, the only remaining non-trivial equation is the constraint equation of the retarded function, Eq. (40). From this retarded equation one can infer the spectrum of the theory. Let us treat the E and B field corrections as perturbations. Then, the retarded Green's function is up to first order in perturbation theory where ψ fields are the solution of the two homogeneous retarded equations, where each is of the Schrödinger equation type, and here spin summation is implicit. The spectral function follows from G ρ = 2 Im iG R . At the second order in perturbation theory the new g 2 correction enters. This is true in vacuum, however, it could be that through the plasma the second order becomes leading order. This needs to be investigated. The fine splitting of the energy levels of the pairs emerges from the content of h in the retarded Eq. (40), while the hyper fine splitting is contained in the term proportional to the one-point function of B field once the classical limit of the B field expectation value is taken. This shows that under all the mentioned limits, that our equations are consistent with basic quantum mechanics of pairs. Eq. (39) is probably consistent with time-dependent quantum mechanics.
Last, we would like to emphasize that the retarded Eq. (40) and Eq. (39) may have some advantages compared to the traditional Schrödinger equation together with the Hilbert space of state description. Both, Schrödinger/Hilbert space description and our retarded, give identical Physics results for hermite potentials. The difference may be that the retarded equation can deal with the more general complex potentials, through the automatic normalization coming from the delta function on the right hand side. While in Schrödingers equation (homogeneous equation) normalization of states has to be done by hand, the retarded Green's function accounts for this through its inhomogeneity. This is needed to discuss non-hermite contributions arising from, e.g., the environment.
To give one example that potentials are in general non-hermitian even in this absolute vacuum limit considered in this section, we choose the top quark. The top quark is unstable and can decay. It would contribute to h with a complex part. And this non-hermitan part is so strong, that no top quark bound-state solutions exist as well known 32 . However, the spectrum is not empty at negative values of the binding energy for top quark pairs. It is continuous due to large width of the bound state. It is this 'vacuum melting' of the top quark bound state pair that may go beyond the traditional Schrödinger/Hilbert space description in this case. At least we do not understand yet how to treat non-hermite potentials in the traditional picture. So at the moment, it is fair to say that the Green's function approach as field moments of the path integral, seems to be slightly more practical.

Heavy Dark Matter in the early Universe plasma
In Ref. 1 , the total number density equation for any N-body state has been derived from non-relativistic effective field theory including pair annihilation: where the total particle number density is defined in terms of non-relativistic (without potential) particle two-component spinors as n η = T C η † (x)η(x) , and the four-point function on the right hand side consists of two particle and two anti-particle operators: . (σ v rel ) 0 is the tree-level annihilation cross section. Since non-relativistic effective field theory is not projected into any sub-space, this equation contains any N-body state. Under the assumption of ionization equilibrium, the non-equilibrium G ++−− ηξ was related to the spectral function via KMS relation with one single chemical potential for the total number of particles (including the ones in the bound states): Now, this assumption of ionization equilibrium is clearly too strong latest when temperature becomes comparable to the ground state binding energy. To go beyond this assumption, we can identify the product f eq G ρ ηξ with our spin singlet contribution of the out-of-equilibrium phase space density Ansatz from section 4.2 and the spin singlet spectral function of two-body fields evaluated at the origin (proportional to Sommerfeld-enhanced annihilation and bound-state decay). This replacement is justified, since the thermal trace can be seen as taking the trace in state space. Since it contains the two-body field contribution, we can set them approximately equal. Note that this is precisely what leads to the pNREFT action. It is the projection of the non-relativistic action into the two-particle sub-space. This completes the theory of Ref. 1

by including any electric and magnetic dipole transitions among two-body pairs.
Notice that this theory has also computational advantages. Instead of having an infinitely coupled system of Boltzmann equations for the particle and bound-state number densities, our formalism only contains three coupled equations. One for the total particle number density in Eq. (52), one for the phase-space distribution of two-body pairs, Eq. (37), and one for the solution of the spectral function from the imaginary part of the retarded Eq. (40). The usual, infinitely coupled, Boltzmann equation and our treatment are equivalent in the vacuum limit of the spectral function (on-shell dispersion relation only). However, our formalism allows to include all bound states and their transitions inside the environment including quasi-particle excitation and multi-mediator emission. Another important example is that it applies to higher excited states which can be melted for temperature around the ground state binding energy (but included, since the spectrum is not empty, rather continues, see e.g., Refs. [18][19][20][21]. Finally, let us remark that we dropped in our pNREFT action from the beginning the short distance annihilation contribution. Although it can also be included straightforwardly, we have treated it as a perturbation, where the leading order term occurs in the dynamical equation for the total number density evolution of particles, i.e., it is contained already in Eq. (52) from integrating out the scale m in NREFT, and we have just perturbed around this contribution by neglecting it in the pNREFT Lagrangian.

Heavy Quarkonia in a quark-gluon plasma
In the literature of Quarkonia inside the quark-gluon plasma (QGP) [33][34][35][36][37][38][39] and heavy Dark Matter (DM) in the early Universe 40,41 , the open quantum system treatment, for deriving the dynamical equations of pairs starting from von Neumann equation for the density matrix, has been established. Let us check, if our Eq. (37) is consistent with some of the previous open quantum system treatment in the expected limit of zero temperature and static limit of the spectral function of the two-body states.
To do so, let us take the static vacuum limit in G ρ . Then, we know that it can only be a function of the energy E = P 0 − p 2 /(2M). Substituting this in the K 0 and P 0 integrals gives for electric dipole transitions in the particle number evolution: with ∆E ≡ E − E ′ (which consistently turns into ∆E = mv 2 rel /2 + |E nlm | for bound-state formation). Now let us pick-up bound-state formation (BSF) by limiting the energy integration:

12/19
Applying this, giveṡ whereṅ 2 | BSF is the time derivative of the scattering particle number. This is because through the limitation of the energy integrals, we consider only positive energies leaving the phase space into the negative energy region, as well as negative energy states entering into the positive as the back reaction term. The ψ fields are under the assumed vacuum limit equivalent to zero temperature Schrödinger wave function field solutions. Changing variables, k → −q + p, we can reproduce Refs. 39,40 : To see this explicitly, notice that p is the integral over the CoM momentum and E can be recast into three dimensional integration over relative momentum (or velocity): +∞ 3 . The first term on the right hand side corresponds to bound-state formation via the emission of an electric field quanta while the last term is the dissociation. Note that through this observation, we can define from Eq. (37) the general dissociation rate from first principles which has physical interpretation. It is not divergent at the leading order in contrast to other previous definitions in finite temperature quarkonia effective field theoretical computations. It is also not the imaginary part of the self-energy, but the product of the imaginary part of the individual retarded functions inside the self-energy.
That our and the open quantum system treatment in this limit fall together was expected. In the latter, there are two main assumptions in addition to the quantum mechanical van Neumann equation: i) perfect vacuum Schrödinger wave function states and ii) factorization of the sub-system and environmental density matrix. Now we have shown here, that the limit of zero temperature in the two-body spectral functions is equivalent with taking both assumptions in the open quantum system treatment (vacuum wave functions plus density matrix factorization). This may be a good approximation for temperature below the Bohr momentum of the ground state, however, can be insufficient for excited states that are much shallower than the ground state and hence can be more probed by the plasma. As an additional remark, it can be seen from the reduction of energy integrals that not only BSF is contained in our Eq. (37), but also any level transitions, Bremsstrahlung (scattering-scattering state transitions) processes, and multi-mediator emission.
To evaluate Eq. (37) more generally by including environmental effects, one needs the two-point spectral functions of the electric and magnetic field correlator and the spectral function of the two-body correlator. All these can be obtained from the imaginary part of the corresponding retarded correlator. The spectral function of the electric field correlator at next-to-leading order, applicable for m D ∆E, was computed in Ref. 40 for QED and for any singlet-adjoint transition with SU(N) quanta in Ref. 41 . For QCD in the regime ∆E = 0 (diffusion aka kinetic equilibration), where Hard-Thermal-Loop resummation is needed, the results can be found in Ref. 42 . Note that in general both results, resummed and fixed NLO, are needed to capture the full energy range of inelastic ultrasoft transitions for any temperature of the medium.

Hydrogen recombination at CMB decoupling
The Hubble tension (H 0 tension) is mainly a mismatch between the todays expansion rate of the Universe, as extrapolated from the European Satelite Planck measurement 43 of the Cosmic Microwave Background (CMB), and indipendent other, local measurements as, e.g., done by Holicow group 44 . While various scrutenizations of our concordance ΛCDM model have been made over the past few years [45][46][47][48][49][50] , less scrutenization of the recombination codes such as RECFAST, CosmoSpec and Rico 51-54 was given. These codes may assume vacuum spectral functions, i.e., the spectral function of the electric, magnetic, or two-body field two-point function is approximated to be in vacuum with zero impact of the surrounding early Universe environment. For example, our one-point and two-point functions of the electric and magnetic field at finite temperature in the spectral function of the two-body field correlator are missing. This impact of the environment on the recombination physics through our equations should be investigated. It is unclear how the energy shifts of all bound state levels, induced from the environment in Eq. (40), influences the free electron number density evolution and hence the CMB evolution, and consequently the todays Hubble value prediction.

Spin-flip temperature evolution at the Cosmic Dawn
The EDGES experiment has pointed out a detection of an anomalous absorption in the 21 cm signal 55 , which is about twice as deep as expected from our concordance ΛCDM model. This lead to intense speculation of Physics beyond the standard model of particle physics and cosmology that could explain this discovered anomaly. Interestingly, there are various close future experiments that have all the potential to either confirm or disprove this discovery.
For hydrogen recombination in the matter-dominated epoch or this spin-flip temperature physics at the Cosmic Dawn, quasi-particle effects may be small, however, our formalism includes infrared and collinear safe scattering of relativistic particles (produced from super-nova explosion) with hydrogen, which may not be included in the current spin-temperature evolution description to the best of our knowledge. Only in thermal field theory, this process is finite 40 , while collinear divergent in ordinary quantum field theory at zero temperature. Also here, the one-point functions in Eq. (40) (order g terms, where the B term gives the hyperfine splitting in vacuum) including finite temperature effects, are currently not taken into account in the predictions. The full impact of the environment needs to be investigated.

Bose-Einstein condensation and Superconductivity
Our six coupled differential equations have features that support the argument that they can be used to study phase transitions dynamically, as well as to find the existence of such extreme phases in, e.g., certain stationary or eqilibrium limits. In particular, introducing one common chemical potential by KMS relation for the total number of two-body states, Eq. (37) has clearly a Bose-Einstein pole 56,57 . Bound states seem to prefer this region, since their (binding) energy has the same sign as the chemical potential for high occupation, which contributes more to the divergence at low temperature. It seems that the system could dynamically move into this state. A numerical solution is required to explicitly show the formation of the Bose-Einstein condensate.
Furthermore, our equations could contain superconductivity phenomena. First theoretical understanding of the superconductivity phenomena in electron many body systems goes back to the Ginzburg Landau theory (GL) 58,59 . They empirically constructed, based on a macroscopic wave function description, one differential equation, that covers many of the first superconductivity features. In summary, the outcome of their work was to phenomenologically understand superconductivity in terms of the magnetic field penetration depth (which goes back to Fritz London 1935, which explains the Meissner effect), as well as the healing length or coherence length of the pairs themselves. Based on the value of the ratio of these two length scales, superconductivity can be classified into type I and type II.
In the ideal case, the electric and magnetic field vanishes in the superconductor, where the latter is the Meissner effect. From Eq. (42) we see that the expectation of the electric field in the system vanishes if the two-body fields O and O † become r independent. This could spontaneously occur in the case of a transition into a superconducting phase, where the coherence length may be constantly large over large relative distances. Note that although the electric field vanishes in such a case, there can be still an electric current due to zero resistance. In case where O and O † become r and spin independent, the expectation of the B field, in Eq. (43), vanishes. In this ideal case, it would imply that the penetration depth is zero. In detail, the B field vanishes because i) the angular momentum operator inside µ is proportional to a derivative with respect to r. Consequently, if the phase φ (t, r) becomes spontaneously global (and the absolute value of O has zero angular momentum and is spin independent) then this contribution is zero. ii) If the system is spin independent, the remaining trace over spin operator (trace over Pauli matrices) vanishes. These two aspects may explain the Meissner effect for superconductivity on a quantum level with two-body field correlators.
Later pioneering works, by Matsubara 60 , Gor'kov 61 , and Baarden et al. 62 , the GL equation could be understood on a more fundamental level. In particular, Matsubara invented the imaginary time formalism (a special case of the general time contour considered in Section 2), so his equation from the original we may recover under some limits including in particular the assumption of kinetic equilibrium and finite chemical potential. Notice that the finite chemical potential in our retarded equation, as well as an explicit mass term can be simply included by shifting P 0 . Consistently, one needs to shift the P 0 integration in the statistical equation as well. One of the main differences of our work to the original by Gor'kov and Baarden et al. is that we do not assume a semi-classical description of the gauge fields. Note also that it seems that the original GL theory, as explained in Eq. (21.6.23) and Eq. (21.6.24) in Weinbergs book the 'Quantum theory of fields Vol.2.', may rely on the classical Euler-Lagrange equation of motion of the Lagrange density. There are of course more recent works like Ref. 63 , which may come close to our coupled Green's functions.
While still, we haven't shown numerically or by finding stationary points of our equations that they contain superconductivity, it is maybe fair to say that in all references or books that we know so far to our best, it is either the semi-classical description, or, missing equations (6 coupled differential equations in total), or missing terms (we have all terms up to g 2 in the correlation function hierarchy, and to all orders in g on the level of interacting two-body field two-point functions), or assumptions on equilibrium which makes them less general compared to this work. Also, we would like to remark at the end that in almost all references, the starting point is more close to NREFT. In this case four point functions are typically needed to be approximated by products of two-point functions, which neglects correlation. This correlation, is fully contained in pNREFT. The projection of the four spinor term, that leads to the four point function in NREFT (potential term), into the two-particle sub-space leads to a quadratic term in pNREFT, which is included in this work.
To repeat, unless our equations are numerically solved and we have a numerical proof, the author of this work does not claim that any of our equations capture the full out-of-equilibrium formation of a Bose-Einstein condensate nor the formation of superconducting phases. To proof that such states form dynamically, we need to solve the equations.

Discussion
We derived field moments of the potential non-relativistic path integral and closed the equations at the two-point function level to all orders in g. We have seen that this generalized eventually multiple physics disciplines, at least, the Dark Matter case and the standard open quantum system treatment of Quarkonia QED interactions inside a plasma.
One interesting aspect, that could need some further discussion (since it is also a recent field), is the co-existence of multiple correlated, condensed, or collapsed states of the same system or mixed systems. For instance, we have seen that Bose-Einstein condensation requires only high occupation number. This may not conflicting with the superconductivity or purely magnetic effects requirement that the two local phases become spontaneously broken and turn into global phases. It could eventually be that all 3 states could exist at the same time. This probably could depend on the order how the state is transformed. We could imagine that bringing the system first into a magnetic phase by applying an external magnetic field, then into superconductivity where the electric field expectation is zero (but non-zero electric current transport) with finite magnetic field (may be possible) and then increasing the occupation number by compressing or cooling the material to ultimately form, e.g., a Ferromagnetic Superconducting Bose-Einstein Condensate.
Currently, we can not make such concrete predictions in this work. The equations need to be solved or sufficiently approximated to proof i) the (co-)existence of such (multiple) extreme states and ii) all the possible formation direction that such states can indeed dynamically be reached. We believe that our work has to some extend contributed into such fascinating direction.

Summary and Conclusion
In summary, we started from the rather recent potential non-relativistic effective field theory [22][23][24] for the case of Quantum Electrodynamics, and derived, based on non-equilibrium quantum field theoretical techniques, evolution equations of the system. This combination, based on pNREFT and non-eq. QFT, is new to us. E.g., in contrast to the standard Kadanoff-Bayam equations which depend on two space arguments, our equations depend on four space arguments since they capture the dynamics of two-body fields as a consequence of the adopted potential non-relativistic effective field theory of pairs. In this sense, our in general 6 new coupled differential equations, representing up to g 2 contributions of the correlator hierarchy the full path integral, needed to undergo a detailed consistency check with many previous physics literature results under certain limits and approximations of our equations. As a result, it seems that too many low energy physics phenomena of QED can be reproduced under certain limits and approximations that this daring attempt is incorrect.
In particular, we started in section 7 by taking the absolute vacuum, static and semi-classical limit (gauge field) of the system. This corresponds also to the historical order the author started slowly to believe in such kind of equations. The outcome of taking such approximations in our equations, is the static quantum mechanics of a pair. Namely, we recovered from the spectral function fine and hyper fine corrections to the pairs energy levels, as well as higher order corrections.
Subsequently, we discussed the time dependent vacuum theory, generalized dark matter pair annihilation at finite temperature 1 for any out-off equilibrium state of the system in the Keldysh-Schwinger formalism, and recovered the open quantum system treatment in quarkonia literature for dilute particles inside an ultra-relativistic environment (such as a quark-gluon plasma) under the approximation of zero temperature in the spectral function. In the latter section, we also demonstrated that not only bound-state formation or dissociation is contained, but any ultra-soft transition, like Bremsstrahlung or any bound-bound transition. Including finite temperature corrections in the spectral function, most of our achieved theory is now new in both fields.
We also emphasized that with the insight of our work, the current free electron decoupling at CMB and spin-flip temperature evolution description at the Cosmic Dawn, deserves a detailed investigation of contributions from the environment (e.g., shifts in energy levels), or, in the latter also contributions from relativistic electrons produced from super-nova explosion to the spin-flip of hydrogen. This is possible now, since the collision terms in the evolution equations are well defined in the framework of thermal field theory where such processes do not suffer from collinear divergences as in ordinary vacuum QFT, see Ref. 40 .
Our consistency check chapter continued with condensed and correlated matter physics. We argued that our equations should be general enough to describe the transition of the system into the Bose-Einstein condensation. It was discussed that this