Density matrix formalism for interacting quantum fields

We provide a description of interacting quantum fields in terms of density matrices for any occupation numbers in Fock space in a momentum basis. As a simple example, we focus on a real scalar field interacting with another real scalar field, and present a practicable formalism for directly computing the density matrix elements of the combined scalar-scalar system. For deriving the main formula, we use techniques from non-equilibrium quantum field theory like thermo field dynamics and the Schwinger-Keldysh formalism. Our results allow for studies of particle creation/annihilation processes at finite times and other non-equilibrium processes including those found in the theory of open quantum systems.

The time evolution of a density operator is governed by the quantum Liouville equation, from which quantum master equations can be derived. Such master equations find heavy use in the theory of open quantum systems [49]. Especially in this context, it is common practice to trace or integrate out some degrees of freedom often referred to as environments. This process results in so-called reduced density matrices, which describe only the subsystems of the total system that were not traced out. In Ref. [14] the authors introduced a Lehmann-Szymanzik-Zimmermann-like reduction (LSZ-like reduction) [50] that provides a practicable and first principle-based formalism for deriving quantum master equations for reduced density matrices in a momentum basis in Fock space (also discussed in Refs. [51,52]). This was done with applications to quenched field theoret-ical open systems in mind. The formalism was then used in Ref. [21] in order to provide a way of directly computing reduced density matrix elements for such systems in terms of Feynman-Vernon influence functionals [53] without having to solve potentially intricate quantum master equations.
In the present article we will follow the path laid out by Ref. [21], but apply the formalism more generally to interacting quantum field theoretical systems without tracing out one of the degrees of freedom. Instead we will show how the formalism introduced in Ref. [14] can be used for computing total density matrix elements, which include information about all systems partaking in the described interactions, for any occupation in Fock space in a momentum basis. Having such a formalism can be a valuable tool, for example, when investigating decay or creation processes at finite times. Moreover, we see manifold direct applications at the high-precision frontier, i.e. in the area of low-energy but high-precision phenomenology as, e.g. in quantum optics experiments, which we will work on in the near future. The current article is a first step towards such applications.
As a simple working example, we will focus our discussion on the interactions of real scalar fields.
Extensions of the presented formalism to other types of fields are possible. For simplicity, for most of our discussion we will only consider the interaction between two scalar fields ϕ and φ. However, the formula that we will later present as a tool for directly computing their total density matrix elements will be extrapolated to a formula allowing to include any number of scalar field species. This article has the following structure: Sec. II deals with the required mathematical concepts, which are then applied for deriving a formula for directly computing total density matrix elements for the combined system of ϕ and φ. A more general formula for any number of scalar field species will also be extrapolated. Subsequently, in Sec. III, a few examples are investigated using the derived formula before the article is concluded in Sec. IV.

II. DERIVATION
In this section we will derive a formula for the direct computation of total density matrices in a momentum basis for any number of particles and species of interacting quantum fields. Since they constitute the simplest example, we will restrict our discussion to real scalar fields. However, similar formalisms can be developed for other types of fields as well. At first, we will only consider two interacting scalar fields ϕ and φ, but generalising the resulting equation to more scalar fields is straightforward.
Initially, we will introduce all required mathematics and concepts in Secs. II A-II C. For this, we will closely follow Ref. [21], which in turn is strongly based on the discussions in Ref. [14]. Finally, we will present the derivation of a formula for the direct computation of total density matrices in Sec. II D.

A. Density matrices in Fock space
Since we are interested in describing the interaction of quantum fields in terms of density matrices, we first need to introduce them in Fock space. For this, we expand the density operator in a momentum basis similarly to Ref. [21], but while taking into account multiple numbers of fields:ρ where we defined multi-indices with α, β,... representing different field or particle species, such that Notice that the cases i α = 0 or j α = 0 correspond to the static vacuum state |0⟩ or ⟨0| for the respective species α. If I = 0 or J = 0, then this describes the case of no particles being present, i.e. the total vacuum. Furthermore, we introduced the short-hand notation for the 3-momenta, and where 1 Later we will also make use of k := The density matrix elements that appear in Eq. (1) are given by and, as usual, have to fulfill: In addition, as was pointed out in Refs. [54,55], density matrix elements are picture-independent.
For later reference, it should also be noted that, in the Schrödinger picture (index S), the time evolution of a density operator is determined by the quantum Liouville equation 2 [49] ∂ ∂tρ which is generally solved bŷ if the Hamiltonian is time-dependent due to an external source.

B. The Schwinger-Keldysh formalism
When computing expectation values for operators at finite times, for example the density matrix elements in Eq. (8), the usual in-out formalism known from scattering amplitudes is not sufficient anymore. Instead, the Schwinger-Keldysh closed-time-path formalism [56,57], also known as inin formalism, is used. A good introduction to this topic can, for example, be found in Ref. [5]. e.g. in Ref. [14], in order to define the so-called Feynman-Vernon influence functional [53]. If we consider the example of two real scalar fields ϕ and φ interacting with each other, then, within the the Schwinger-Keldysh formalism, we can define the density functional 3 where ± indicates a dependence on both +-and −-type operators, and subscript t labels the time slice on which the field eigenstates have to be taken [21]. The time evolution of such a density functional is given in terms of the functional propagator 4 , with the action Here, indicates functionals that depend on both field variables ϕ + and ϕ − , such that The action S α [α|t] is for a free α-field, S α,int [α|t] is the corresponding self-interaction action, and S int [ϕ; φ|t] is the action describing the interaction between the fields ϕ and φ. Since we are only 3 The extension to more field species is straightforward. 4 Compare it to the influence functional propagator in the literature, e.g. in Ref. [5].
interested in finite times, we restrict all actions to Ω t := [0, t] × R 3 , such that, schematically, we use Furthermore, notice the notational distinction between functional integrals over all of R 3 on a particular time slice, denoted by d as in Eq. (13), and over the whole Ω t , denoted by D as in Eq. (14).
For this, we will directly quote from most of the corresponding elaborations in Ref. [21]. One way of understanding TFD is as an algebraic formulation of the Schwinger-Keldysh formalism, in the sense that it works with operators on a positive (+) or negative (−) complex time axis acting on a doubled Hilbert space where H ± are the Hilbert spaces corresponding to the ±-branches of the closed time path. Operators living on the closed time path can be expressed within TFD aŝ whereÎ is the unit operator and T means time reversal. States in a momentum basis in TFD can be reached from the doubled vacuum state via creation operatorŝ Consequently, the corresponding annihilators act likê which allows us to express the expectation value of an operator as Having all this, we can now use the fact that Eq. (10) can be re-written in a Schrödinger-like form: where . This can generally be solved bŷ Note that at time 0 the different pictures coincide and we therefore dropped the label S.
We recall that, generally, Hamiltonian and action are related via H(t) = − ∂ ∂t S(t). Furthermore, if we consider that the action in Eq. (15) describes the full evolution of the density matrix elements in the field-basis, the interaction picture H I must be its corresponding Hamiltonian, such that the action S[ϕ; φ|t] as appearing in Eq. (16). This will later be of greater importance when we write down a path integral expression for the density matrix elements. For now, however, we just remind ourselves that the free HamiltonianĤ 0 is the same in Schrödinger and interaction picture, and consequently drop the subscript.
Translating Eq. (26) into the interaction picture, and using the fact that the state |1⟩⟩ is actually time-and picture-independent [21], therefore leaves us with which can be solved bŷ

D. Density matrix elements
We can now use the concepts discussed in the previous sections in order to derive an equation for the direct computation of elements of the total density matrix for two interacting real scalar fields ϕ and φ in a momentum basis for any occupation number in Fock space. Doing so, we will only work in the interaction picture and therefore not use any corresponding labels for operators and states anymore.
Initially, we will derive a formula that allows for the computation of the density matrix element ρ 1,1;1,1 (p; k|p ′ ; k ′ |t), which represents the state of a single ϕ-particle (momenta p and p ′ ) and a single φ-particle (momenta k and k ′ ). Having found this expression, it will later be possible to extrapolate a more general formula for any number of particle species and respective occupation numbers in momentum space.
As in Refs. [14,21], our starting point is Using Eq. (25), Eq. (30) can be rewritten as in TFD language. Next, we substitute Eq. (29) and find In order to proceed, we need to expand the density operator on the right-hand side of Eq. (32) as in Eq. (1). However, in this way we would encounter the problem of potentially having to deal with an infinite number of initial density matrix elements. Luckily, we can make the in many experimental situations reasonable assumption that we know the initial setup sufficiently well, such that we can ignore all vanishing or suppressed occupations. In our particular example, we will assume that the total initial system already consists of only one ϕ-and one φ-particle. Consequently, we only have to consider the initial density matrix element for the case i ϕ = i φ = j ϕ = j φ = 1 and assume all others to be nil. We do this assumption for the sake of readability, but a more general result will later be presented, which could be obtained by modifying the present derivation accordingly.
Under the considered assumptions we are led to which, using Eqs. (22) and (23) becomes Here we denoted creators/annihilators for ϕ withâ † /â and those for φ withb † /b. Next, we replace them by (cf. Ref. [14]) where ∂ t,E ϕ p := → ∂ t − iE ϕ p , and there exist corresponding expressions forb † /b. This leaves us with As in Ref. [21], we introduced limits in which x 0(′) ϕ,φ approach t from above and y 0(′) ϕ,φ approach 0 from below in order to recover the correct time ordering.
We translate the expression in Eq. (36) into the path integral formalism, which yields where we defined and used The Eq. (37) we just found allows us to directly compute the density matrix element that describes the interaction between two single particles ϕ and φ. Taking Eq. (37), we extrapolate this formula for a general density matrix element with any number of particle species and occupations in Fock space. It can easily be seen that this result can be derived by extending the derivation for Eq. (37) accordingly. We find: In Eq. (40) we used some short-hand notations from Sec. II A and introduced a number new ones.
It is important to remember that, as was already pointed out in Ref. [21], only the diagonal elements of the 2 × 2-matrix propagator are permitted when evaluating Eqs. (37) or (40). This means, when applying Wick's theorem [62], only contractions of fields living on the same branch of the closed time path (cf. Fig. 1) are possible, such that we only have to work with Feynman and Dyson propagators while the other two-point functions vanish: In this section we will apply the formulas in Eqs. (37) and (40) to a toy model example with two real scalar field species ϕ and χ. More precisely, we will consider solely a single particle ϕ at the initial time 0, which is described by the only non-vanishing density matrix ρ 1,0;1,0 (0). We will assume that due to an interaction potential V ∼ φ 2 ϕ the state of two φ-particles can become excited within the interval between times 0 and t. This leads to other density matrices potentially being non-vanishing as well at time t. Of those we will compute the elements of the matrices ρ 0,2;0,2 (t), ρ 1,0;0,2 (t) and ρ 1,0;1,0 (t). Each of these matrices has its own physical interpretation: the first one corresponds to a 1 → 2-particle decay event, the second describes the correlation between a single ϕ-particle state and a two φ-particles state, while the third stands for the single ϕ-particle state continuing to exist, but with a modified probability. The last of these events, i.e. the change of ρ 1,0;1,0 over time, can be compared to the dynamics of an open quantum system.
It should be noted that when computing the density matrices we will sometimes encounter divergences which would usually require renormalization. However, we will only present the unrenormalized expressions here since, as was already mentioned in Ref. [21], the peculiar time-dependent divergences, which we will encounter, have only shortly [14] but not yet sufficiently been discussed in the literature. Dealing with such divergences is beyond the scope of the current article and will instead be the subject of a future work [63].
We will use the free actions for the scalar fields ϕ and φ. As interaction action we choose to use where α ≪ 1 and M is some sufficiently small mass scale. For simplicity, we do not consider selfinteractions. However, including them for more elaborate examples would of course be straightforward, but lead to longer and less readable results.
In every considered case we will have to use the perturbative expansion of the exponentiated interaction action from Eq. (56). We decide to work up to second order in α, such that we find: where in this and all following equations z and z ′ are integrated over the domain Ω t . We begin our investigations with the decay of a single ϕ-particle into two φ-particles, and compute the density matrix elements for this process. Obviously the only non-vanishing density matrix at time 0 must be ρ 1,0;1,0 , while we want to find the matrix ρ 0,2;0,2 at time t. Following Eq. (40), we must evaluate ρ 0,2;0,2 (; p, k|; p ′ , k ′ |t) = lim Substituting Eq. (57), we are left with Using the propagators given in this expression, we can easily depict the contributing physical processes in Figs. 2(a)-(d). We observe that one (Figs. 2(b) and (c)) or two ( Fig. 2(a)) of the ϕpropagators end in φ-tadpoles, while, in addition, equally many propagators at equal time but different positions are appearing. Those terms are of course divergent. However, the term diagrammatically represented in Fig. 2(d) is finite and corresponds to the actual decay of a single ϕ-particle into two φ-particles without the appearance of any unphysical terms.
Of course Eq. (59) can be further evaluated. A brief outline with more details on how to do this can be found in Appendix D of Ref. [52]. Following this procedure, we finally find: We observe that, as expected, the terms in the fourth, sixth and eighth lines of Eq. (60) contain differences in the unitary evolutions of a single ϕ-particle and two φ-particles.
Next, we compute one of the two density matrices which describe correlations between a single ϕ-particle and two φ-particles. The other one can easily be obtained by making use of the property given in Eq. (9). Starting with Eq. (40), we find: Evaluating the path integrals leads us to which can be used to depict the processes as in Figs. 3(a) and (b). As in the previous case, we observe one of the ϕ-propagators ending in a φ-tadpole, while an equal-time φ-propagator is appearing, in Fig. 3(a). In contrast, Fig. 3(b) contains only physical terms without any divergences.
Continuing the computation, we find: where, in the last line, we again find the difference in the unitary evolutions of a single ϕ-particle and two φ-particles.
Finally, we study the change of ρ 1,0;1,0 over time. We use Eq. (40) to find: (a) shows one of the ϕ-propagators ending in a φ-tadpole, while an equal-time φ-propagator is appearing, and another free ϕ-propagator. (b) contains one ϕ-propagator decaying into two φ-particles and a free ϕ-propagator.
which becomes Ref. [21]. A corresponding computation in this article, on the other hand, would lead to the appearance of Feynman and Dyson propagators only and consequently to different final results.

IV. CONCLUSIONS AND OUTLOOK
Density matrices are powerful tools in non-relativistic quantum mechanics, but also in quantum field theory. However, analytically solving the quantum master equations that describe the time evolution of density matrices can be extremely complicated or even hopeless.
In this article we used the ideas originally introduced in Refs. [14] and [21] for developing a practicable and first principle-based density matrix formalism for interacting quantum fields. It has the advantage of circumventing quantum master equations by providing a relatively simple formula for directly computing density matrices in a momentum basis for any number of field species and occupations in Fock space.
As a working example, we chose the setup of two different real scalar fields interacting with each other, but also provided an equation that allows for the inclusion of arbitrarily more field species.
Extending the formalism to other types of fields, e.g. complex scalar fields, vector fields etc., is rather straightforward. Before deriving the relevant formulas, we discussed density matrices in Fock space, see Sec. II A, the Schwinger-Keldysh formalism, Sec. II B, and thermo field dynamics, Sec. II C. Those were the ingredients required for the actual derivation, which we concluded in Sec. II D. Finally, in Sec. III, we applied the developed formalism to a selected example.
While the formalism in Ref. [21] is specifically designed for open quantum systems and therefore includes the tracing-out of environmental degrees of freedom, the formalism presented here captures all degrees of freedom equally without loss of any information. This is very useful, for example, for studies of creation and annihilation processes at finite times, including particle decays. Of course the formalism presented in this article also naturally leads to similar time-dependent divergences as in Ref. [21], which will be discussed in more detail in a future work [63]. However, while in Ref. [21] every computation led to the emergence of such divergences due to the process of tracing out the environmental degrees of freedom, in this article we encountered only a few of those at the considered order in perturbation theory since we allowed for external legs of all involved scalar fields to appear.
All in all, the present work is an important first step towards a description of physically relevant interacting quantum field theory systems at finite times. In particular, we see imminent applications of the presented formalism in a variety of areas at the high-precision frontier, i.e. in the area