Density Matrix Formalism for Interacting Quantum Fields

Abstract

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.

1. Introduction

Density operators and matrices are powerful tools for describing the quantum state of a physical system since they not only allow for the description of pure states, but also of mixed ones. The latter is a feature not accessible with descriptions by wave functions [1]. Consequently, besides their uses in non-relativistic quantum mechanics, density matrices also find ample applications in (open) quantum field theory [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] and related areas, including cosmology [24,25,26,27,28,29,30,31,32,33,34,35,36,37], black holes [38,39,40,41,42] or heavy-ion physics [43,44,45,46,47,48].

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 theoretical 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 $\varphi$ and $\phi$. 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: Section 2 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 $\varphi$ and $\phi$. A more general formula for any number of scalar field species will also be extrapolated. Subsequently, in Section 3, a few examples are investigated using the derived formula before the article is concluded in Section 4.

2. 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 $\varphi$ and $\phi$, but generalizing the resulting equation to more scalar fields is straightforward.

Initially, we will introduce all required mathematics and concepts in Section 2.1, Section 2.2 and Section 2.3. 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 Section 2.4.

2.1. 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:

$$\begin{array}{ccc}\hfill \widehat{\rho }\left(t\right)& =& \sum _{I,J=0}^{\infty }\frac{1}{I!J!}\int \left(\prod _{A=1}^{I}d{\mathsf{\Pi }}_{K^A}\right)\left(\prod _{B=1}^{J}d{\mathsf{\Pi }}_{L^B}\right){\rho }_{I;J}(K^I|L^J\left|t\right)|K^I\rangle \langle L^J|,\hfill \end{array}$$

(1)

where we defined multi-indices

$$\begin{array}{ccc}\hfill I& :=& (i_{\alpha },i_{\beta },\cdots ),J:=(j_{\alpha },j_{\beta },\cdots )\hfill \end{array}$$

(2)

with $\alpha$, $\beta$, … representing different field or particle species, such that

$$\begin{array}{ccc}\hfill I!& :=& i_{\alpha }!i_{\beta }!\cdots .\hfill \end{array}$$

(3)

Notice that the cases $i_{\alpha }=0$ or $j_{\alpha }=0$ correspond to the static vacuum state $|0\rangle$ or $\langle 0|$ for the respective species $\alpha$. 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

$$\begin{array}{ccc}\hfill K^I& :=& {\mathbf{k}}_{\alpha }^{\left(1\right)},\cdots ,{\mathbf{k}}_{\alpha }^{\left(i_{\alpha }\right)};{\mathbf{k}}_{\beta }^{\left(1\right)},\cdots ,{\mathbf{k}}_{\beta }^{\left(i_{\beta }\right)};\cdots ,L^J:={\mathbf{l}}_{\alpha }^{\left(1\right)},\cdots ,{\mathbf{l}}_{\alpha }^{\left(j_{\alpha }\right)};{\mathbf{l}}_{\beta }^{\left(1\right)},\cdots ,{\mathbf{l}}_{\beta }^{\left(j_{\beta }\right)};\cdots \hfill \end{array}$$

(4)

for the 3-momenta, and

$$\begin{array}{ccc}\hfill \prod _{A=1}^{I}d{\mathsf{\Pi }}_{K^A}& :=& \left(\prod _{a=1}^{i_{\alpha }}d{\mathsf{\Pi }}_{{\mathbf{k}}_{\alpha }^{\left(a\right)}}\right)\left(\prod _{a=1}^{i_{\beta }}d{\mathsf{\Pi }}_{{\mathbf{k}}_{\beta }^{\left(a\right)}}\right)\cdots ,\hfill \end{array}$$

(5)

where

$$\begin{array}{ccc}\hfill \int d{\mathsf{\Pi }}_{{\mathbf{k}}_{\alpha }}& :=& {\int }_{{\mathbf{k}}_{\alpha }}\frac{1}{2E_{{\mathbf{k}}_{\alpha }}^{\alpha }}\hfill \end{array}$$

(6)

with1

$$\begin{array}{ccc}\hfill {\int }_{\mathbf{k}}& :=& \int \frac{d^{3}k}{{\left(2\pi \right)}^3}.\hfill \end{array}$$

(7)

The density matrix elements that appear in Equation (1) are given by

$$\begin{array}{ccc}\hfill {\rho }_{I;J}(K^I|L^J\left|t\right)& =& \langle K^I|\widehat{\rho }\left(t\right)|L^J\rangle ,\hfill \end{array}$$

(8)

and, as usual, have to fulfill:

$$\begin{array}{ccc}\hfill {\rho }_{I;J}(K^I|L^J\left|t\right)& =& {\rho }_{J;I}^{*}(L^J|K^I\left|t\right).\hfill \end{array}$$

(9)

In addition, as was pointed out in Refs. [54,55], density matrix elements are picture-independent. We physically interpret them in the following way: ${\rho }_{0;0}$ describes a 0-particle or vacuum state, ${\rho }_{(i_{\alpha },\cdots );0}$ or ${\rho }_{0;(j_{\alpha },\cdots )}$ stand for correlations between $i_{\alpha }$- or $j_{\alpha }$-$\alpha$-particles with the vacuum, while ${\rho }_{(i_{\alpha },\cdots );(j_{\alpha },\cdots )}$ represent correlations between $i_{\alpha }$- and $j_{\alpha }$-$\alpha$-particle states.

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 equation2 [49].

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}{\widehat{\rho }}_S\left(t\right)& =& -\mathrm{i}[{\widehat{H}}_S\left(t\right),{\widehat{\rho }}_S\left(t\right)],\hfill \end{array}$$

(10)

which is generally solved by

$$\begin{array}{ccc}\hfill {\widehat{\rho }}_S\left(t\right)& =& \left(\mathrm{T}{e}^{-\mathrm{i}{\int }_0^{t}d\tau {\widehat{H}}_S\left(\tau \right)}\right)\widehat{\rho }\left(0\right)\left(\tilde{\mathrm{T}}{e}^{\mathrm{i}{\int }_0^{t}d\tau {\widehat{H}}_S\left(\tau \right)}\right)\hfill \end{array}$$

(11)

if the Hamiltonian is time-dependent due to an external source.

2.2. The Schwinger-Keldysh Formalism

When computing expectation values for operators at finite times, for example the density matrix elements in Equation (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 in–in formalism, is used. A good introduction to this topic can, for example, be found in Ref. [5]. The Schwinger–Keldysh formalism is essentially relying on doubling the degrees of freedom, where the two copies are distinguished by labels +/−, and letting those evolve on the positive/negative branch of the closed time path depicted in Figure 1 between an initial $t_{\mathrm{initial}}=0$ and a final time $t_{\mathrm{final}}=t$ [21]. While for scattering amplitudes in- and out-states are used, which are taken to be asymptotic, the path-integral representation of the trace of an operator at a finite time requires us to consider two copies of field-states and therefore naturally leads to the closed time path. This formalism is often applied in the context of field theoretical descriptions of open quantum systems, 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 $\varphi$ and $\phi$ interacting with each other, then, within the the Schwinger–Keldysh formalism, we can define the density functional3

$$\begin{array}{ccc}\hfill \rho [{\varphi }_t^{\pm };{\phi }_t^{\pm }|t]& :=& \langle {\varphi }_t^{+},{\phi }_t^{+}\left|\widehat{\rho }\left(t\right)\right|{\varphi }_t^{-},{\phi }_t^{-}\rangle ,\hfill \end{array}$$

(12)

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

$$\begin{array}{ccc}\hfill \rho [{\varphi }_t^{\pm };{\phi }_t^{\pm }|t]& =& \int d{\varphi }_0^{\pm }d{\phi }_0^{\pm }\mathcal{P}[{\varphi }_t^{\pm },{\varphi }_0^{\pm };{\phi }_t^{\pm },{\phi }_0^{\pm }|t,0]\rho [{\varphi }_0^{\pm };{\phi }_0^{\pm }|0]\hfill \end{array}$$

(13)

is given in terms of the functional propagator4,

$$\begin{array}{ccc}\hfill \mathcal{P}[{\varphi }_t^{\pm },{\varphi }_0^{\pm };{\phi }_t^{\pm },{\phi }_0^{\pm }|t,0]& =& {\int }_{{\varphi }_0^{\pm }}^{{\varphi }_t^{\pm }}\mathcal{D}{\varphi }^{\pm }{\int }_{{\phi }_0^{\pm }}^{{\phi }_t^{\pm }}\mathcal{D}{\phi }^{\pm }{e}^{\mathrm{i}\widehat{S}[\varphi ;\phi |t]}\hfill \end{array}$$

(14)

with the action

$$\begin{array}{ccc}\hfill \widehat{S}[\varphi ;\phi |t]& =& {\widehat{S}}_{\varphi }\left[\varphi \right|t]+{\widehat{S}}_{\phi }\left[\phi \right|t]+{\widehat{S}}_{\varphi ,\mathrm{int}}\left[\varphi \right|t]+{\widehat{S}}_{\phi ,\mathrm{int}}\left[\phi \right|t]+{\widehat{S}}_{\mathrm{int}}[\varphi ;\phi |t].\hfill \end{array}$$

(15)

Here, $\widehat{}$ indicates functionals that depend on both field variables ${\varphi }^{+}$ and ${\varphi }^{-}$, such that

$$\begin{array}{ccc}\hfill \widehat{S}[\varphi ;\phi |t]& :=& S[{\varphi }^{+};{\phi }^{+}|t]-S[{\varphi }^{-};{\phi }^{-}|t].\hfill \end{array}$$

(16)

The action ${\widehat{S}}_{\alpha }\left[\alpha \right|t]$ is for a free $\alpha$-field, ${\widehat{S}}_{\alpha ,\mathrm{int}}\left[\alpha \right|t]$ is the corresponding self-interaction action, and ${\widehat{S}}_{\mathrm{int}}[\varphi ;\phi |t]$ is the action describing the interaction between the fields $\varphi$ and $\phi$. Since we are only interested in finite times, we restrict all actions to ${\mathsf{\Omega }}_t:=[0,t]\times {\mathbb{R}}^3$, such that, schematically, we use

$$\begin{array}{ccc}\hfill S\left[t\right]& =& {\int }_{x\in {\mathsf{\Omega }}_t}\mathcal{L}\left[x\right]\hfill \end{array}$$

(17)

with

$$\begin{array}{ccc}\hfill {\int }_{x\in {\mathsf{\Omega }}_t}& :=& {\int }_{{\mathsf{\Omega }}_t}{d}^{4}x.\hfill \end{array}$$

(18)

Furthermore, notice the notational distinction between functional integrals over all of ${\mathbb{R}}^3$ on a particular time slice, denoted by d as in Equation (13), and over the whole ${\mathsf{\Omega }}_t$, denoted by $\mathcal{D}$ as in Equation (14).

2.3. Thermo-Field Dynamics

The final concept we need to discuss is thermo-field dynamics (TFD) [58,59,60] (see also Ref. [61]). 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

$$\begin{array}{ccc}\hfill \widehat{\mathcal{H}}& :=& {\mathcal{H}}^{+}\otimes {\mathcal{H}}^{-},\hfill \end{array}$$

(19)

where ${\mathcal{H}}^{\pm }$ 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 as

$$\begin{array}{ccc}\hfill {\widehat{\mathcal{O}}}^{+}& =& \widehat{\mathcal{O}}\otimes \widehat{\mathbb{I}},\hfill \\ \hfill {\widehat{\mathcal{O}}}^{-}& =& \widehat{\mathbb{I}}\otimes {\widehat{\mathcal{O}}}^{\mathcal{T}},\hfill \end{array}$$

(20)

where $\widehat{\mathbb{I}}$ is the unit operator and $\mathcal{T}$ means time reversal. States in a momentum basis in TFD can be reached from the doubled vacuum state

$$\begin{array}{ccc}\hfill |0\rangle \rangle & :=& |0\rangle \otimes |0\rangle \hfill \end{array}$$

(21)

via creation operators

$$\begin{array}{c}\hfill {\widehat{a}}_{\mathbf{k}}^{+†}|0\rangle \rangle =|\mathbf{k}\rangle \otimes |0\rangle =:|{\mathbf{k}}_{+}\rangle \rangle ,{\widehat{a}}_{\mathbf{k}}^{-†}|0\rangle \rangle =|0\rangle \otimes |\mathbf{k}\rangle =:|{\mathbf{k}}_{-}\rangle \rangle .\end{array}$$

(22)

Consequently, the corresponding annihilators act like

$$\begin{array}{c}\hfill {\widehat{a}}_{\mathbf{k}}^{\pm }|{\mathbf{p}}_{+},{\mathbf{p}}_{-}\rangle \rangle ={\left(2\pi \right)}^{3}2E_{\mathbf{k}}^{\varphi }{\delta }^{\left(3\right)}(\mathbf{p}-\mathbf{k})|{\mathbf{p}}_{\mp }\rangle \rangle .\end{array}$$

(23)

In addition, we can construct a special state corresponding to the unit operator, see Ref. [59],

$$\begin{array}{ccc}\hfill |1\rangle \rangle & :=& |0\rangle \rangle +\int d{\mathsf{\Pi }}_{\mathbf{k}}|{\mathbf{k}}_{+},{\mathbf{k}}_{-}\rangle \rangle +\frac{1}{2!}\int d{\mathsf{\Pi }}_{\mathbf{k}}d{\mathsf{\Pi }}_{{\mathbf{k}}^{\prime }}|{\mathbf{k}}_{+},{\mathbf{k}}_{+}^{\prime },{\mathbf{k}}_{-},{\mathbf{k}}_{-}^{\prime }\rangle \rangle +\cdots ,\hfill \end{array}$$

(24)

which allows us to express the expectation value of an operator as

$$\begin{array}{ccc}\hfill \langle \widehat{\mathcal{O}}\left(t\right)\rangle & =& \mathrm{Tr}\widehat{\mathcal{O}}\left(t\right)\widehat{\rho }\left(t\right)=\langle \langle 1|{\widehat{\mathcal{O}}}^{+}\left(t\right){\widehat{\rho }}^{+}\left(t\right)|1\rangle \rangle .\hfill \end{array}$$

(25)

Having all this, we can now use the fact that Equation (10) can be re-written in a Schrödinger-like form:

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}{\widehat{\rho }}_S^{+}{\left(t\right)|1\rangle \rangle }_S& =& -\mathrm{i}{\widehat{H}}_S\left(t\right){\widehat{\rho }}_S^{+}{\left(t\right)|1\rangle \rangle }_S,\hfill \end{array}$$

(26)

where ${\widehat{H}}_S\left(t\right):={\widehat{H}}_S\left(t\right)\otimes \mathbb{I}-\mathbb{I}\otimes {\widehat{H}}_S\left(t\right)$. This can generally be solved by

$$\begin{array}{ccc}\hfill {\widehat{\rho }}_S^{+}{\left(t\right)|1\rangle \rangle }_S& =& \mathrm{T}exp\left\{-\mathrm{i}\underset{0}{\overset{t}{\int }}{\widehat{H}}_S\left(\tau \right)d\tau \right\}{\widehat{\rho }}^{+}{\left(0\right)|1\rangle \rangle }_S.\hfill \end{array}$$

(27)

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\left(t\right)=-\frac{\partial }{\partial t}S\left(t\right)$. Furthermore, if we consider that the action in Equation (15) describes the full evolution of the density matrix elements in the field-basis, the interaction picture ${\widehat{H}}_I$ must be its corresponding Hamiltonian, such that ${\widehat{H}}_I\left(t\right):={\widehat{H}}_{0,I}\left(t\right)+{\widehat{H}}_{\varphi ,\mathrm{int},I}\left(t\right)+{\widehat{H}}_{\phi ,\mathrm{int},I}\left(t\right)+{\widehat{H}}_{\mathrm{int},I}\left(t\right)$ with ${\widehat{H}}_{0,I}\left(t\right):={\widehat{H}}_{\varphi ,0,I}\left(t\right)+{\widehat{H}}_{\phi ,0,I}\left(t\right)$ representing the action $S[\varphi ;\phi |t]$ as appearing in Equation (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 ${\widehat{H}}_0$ is the same in Schrödinger and interaction picture, and consequently drop the subscript.

Translating Equation (26) into the interaction picture, and using the fact that the state $|1\rangle \rangle$ is actually time- and picture-independent [21], therefore leaves us with

$$\begin{array}{ccc}\hfill {\partial }_t{\widehat{\rho }}_I^{+}\left(t\right)|1\rangle \rangle & =& -\mathrm{i}[{\widehat{H}}_{\varphi ,\mathrm{int},I}\left(t\right)+{\widehat{H}}_{\phi ,\mathrm{int},I}\left(t\right)+{\widehat{H}}_{\mathrm{int},I}\left(t\right)]{\widehat{\rho }}_I^{+}\left(t\right)|1\rangle \rangle ,\hfill \end{array}$$

(28)

which can be solved by

$$\begin{array}{ccc}\hfill {\widehat{\rho }}_I^{+}\left(t\right)|1\rangle \rangle & =& \mathrm{T}exp\left\{-\mathrm{i}\underset{0}{\overset{t}{\int }}[{\widehat{H}}_{\varphi ,\mathrm{int},I}\left(\tau \right)+{\widehat{H}}_{\phi ,\mathrm{int},I}\left(\tau \right)+{\widehat{H}}_{\mathrm{int},I}\left(\tau \right)]d\tau \right\}{\widehat{\rho }}^{+}\left(0\right)|1\rangle \rangle .\hfill \end{array}$$

(29)

2.4. 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 $\varphi$ and $\phi$ 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 ${\rho }_{1,1;1,1}(\mathbf{p};\mathbf{k}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime }\left|t\right)$, which represents the state of a single $\varphi$-particle (momenta $\mathbf{p}$ and ${\mathbf{p}}^{\prime }$) and a single $\phi$-particle (momenta $\mathbf{k}$ and ${\mathbf{k}}^{\prime }$). 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

$$\begin{array}{ccc}\hfill \langle \mathbf{p};\mathbf{k};t\left|\widehat{\rho }\left(t\right)\right|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime };t\rangle & =& {\rho }_{1,1;1,1}(\mathbf{p};\mathbf{k}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime }\left|t\right).\hfill \end{array}$$

(30)

Using Equation (25), Equation (30) can be rewritten as

$$\begin{array}{ccc}\hfill \mathrm{Tr}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime };t\rangle \langle \mathbf{p};\mathbf{k};t|\widehat{\rho }\left(t\right)& =& \langle \langle 1|\left(\right|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime };t\rangle \langle \mathbf{p};\mathbf{k};t|\otimes \widehat{\mathbb{I}}){\widehat{\rho }}^{+}\left(t\right)|1\rangle \rangle \hfill \end{array}$$

(31)

in TFD language. Next, we substitute Equation (29) and find

$$\begin{array}{ccc}\hfill {\rho }_{1,1;1,1}(\mathbf{p};\mathbf{k}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime }\left|t\right)& =& \langle \langle {\mathbf{p}}_{+},{\mathbf{p}}_{-}^{\prime };{\mathbf{k}}_{+},{\mathbf{k}}_{-}^{\prime };t|\hfill \\ & & \times \mathrm{T}exp\left\{-\mathrm{i}\underset{0}{\overset{t}{\int }}[{\widehat{H}}_{\varphi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\phi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\mathrm{int}}\left(\tau \right)]d\tau \right\}{\widehat{\rho }}^{+}\left(0\right)|1\rangle \rangle .\hfill \end{array}$$

(32)

In order to proceed, we need to expand the density operator on the right-hand side of Equation (32) as in Equation (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 $\varphi$- and one $\phi$-particle. Consequently, we only have to consider the initial density matrix element for the case $i_{\varphi }=i_{\phi }=j_{\varphi }=j_{\phi }=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

$$\begin{array}{ccc}\hfill {\rho }_{1,1;1,1}(\mathbf{p};\mathbf{k}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime }\left|t\right)& =& \langle \langle {\mathbf{p}}_{+},{\mathbf{p}}_{-}^{\prime };{\mathbf{k}}_{+},{\mathbf{k}}_{-}^{\prime };t|\mathrm{T}exp\left\{-\mathrm{i}\underset{0}{\overset{t}{\int }}[{\widehat{H}}_{\varphi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\phi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\mathrm{int}}\left(\tau \right)]d\tau \right\}\hfill \\ & & \times \int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}d{\mathsf{\Pi }}_{\mathbf{l}}d{\mathsf{\Pi }}_{{\mathbf{l}}^{\prime }}{\rho }_{1,1;1,1}(\mathbf{q};\mathbf{l}|{\mathbf{q}}^{\prime };{\mathbf{l}}^{\prime }\left|0\right)|{\mathbf{q}}_{+},{\mathbf{q}}_{-}^{\prime };{\mathbf{l}}_{+},{\mathbf{l}}_{-}^{\prime }\rangle \rangle ,\hfill \end{array}$$

(33)

which, using Equations (22) and (23) becomes

$$\begin{array}{ccc}\hfill {\rho }_{1,1;1,1}(\mathbf{p};\mathbf{k}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime }\left|t\right)& =& \int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}d{\mathsf{\Pi }}_{\mathbf{l}}d{\mathsf{\Pi }}_{{\mathbf{l}}^{\prime }}{\rho }_{1,1;1,1}(\mathbf{q};\mathbf{l}|{\mathbf{q}}^{\prime };{\mathbf{l}}^{\prime }\left|0\right)\langle \langle 0|\mathrm{T}{\widehat{a}}_{\mathbf{p}}^{+}\left(t\right){\widehat{a}}_{{\mathbf{p}}^{\prime }}^{-}\left(t\right){\widehat{b}}_{\mathbf{k}}^{+}\left(t\right){\widehat{b}}_{{\mathbf{k}}^{\prime }}^{-}\left(t\right)\hfill \\ & & \times exp\left\{-\mathrm{i}\underset{0}{\overset{t}{\int }}[{\widehat{H}}_{\varphi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\phi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\mathrm{int}}\left(\tau \right)]d\tau \right\}\hfill \\ & & \times {\widehat{a}}_{\mathbf{q}}^{+†}\left(0\right){\widehat{a}}_{{\mathbf{q}}^{\prime }}^{-†}\left(0\right){\widehat{b}}_{\mathbf{l}}^{+†}\left(0\right){\widehat{b}}_{{\mathbf{l}}^{\prime }}^{-†}\left(0\right)|0\rangle \rangle .\hfill \end{array}$$

(34)

Here we denoted creators/annihilators for $\varphi$ with ${\widehat{a}}^{†}$/$\widehat{a}$ and those for $\phi$ with ${\widehat{b}}^{†}$/$\widehat{b}$. Next, we replace them by (cf. Ref. [14])

$$\begin{array}{ccc}\hfill {\widehat{a}}_{\mathbf{p}}^{+}\left(t\right)& =& +\mathrm{i}{\int }_{\mathbf{x}}{e}^{-\mathrm{i}\mathbf{p}·\mathbf{x}}{\partial }_{t,E_{\mathbf{p}}^{\varphi }}{\widehat{\varphi }}^{+}(t,\mathbf{x}),{\widehat{a}}_{\mathbf{p}}^{+†}\left(t\right)=-\mathrm{i}{\int }_{\mathbf{x}}{e}^{+\mathrm{i}\mathbf{p}·\mathbf{x}}{\partial }_{t,E_{\mathbf{p}}^{\varphi }}^{*}{\widehat{\varphi }}^{+}(t,\mathbf{x}),\hfill \\ \hfill {\widehat{a}}_{\mathbf{p}}^{-}\left(t\right)& =& -\mathrm{i}{\int }_{\mathbf{x}}{e}^{+\mathrm{i}\mathbf{p}·\mathbf{x}}{\partial }_{t,E_{\mathbf{p}}^{\varphi }}^{*}{\widehat{\varphi }}^{-}(t,\mathbf{x}),{\widehat{a}}_{\mathbf{p}}^{-†}\left(t\right)=+\mathrm{i}{\int }_{\mathbf{x}}{e}^{-\mathrm{i}\mathbf{p}·\mathbf{x}}{\partial }_{t,E_{\mathbf{p}}^{\varphi }}{\widehat{\varphi }}^{-}(t,\mathbf{x}),\hfill \end{array}$$

(35)

where ${\partial }_{t,E_{\mathbf{p}}^{\varphi }}:={\overrightarrow{\partial }}_t-\mathrm{i}{E}_{\mathbf{p}}^{\varphi }$, and there exist corresponding expressions for ${\widehat{b}}^{†}$/$\widehat{b}$. This leaves us with

$$\begin{array}{ccc}\hfill {\rho }_{1,1;1,1}(\mathbf{p};\mathbf{k}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime }\left|t\right)& =& \underset{\begin{array}{c}{x}_{\varphi ,\phi }^{0\left({}^{\prime }\right)}\to t^{+}\\ y_{\varphi ,\phi }^{0\left({}^{\prime }\right)}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}d{\mathsf{\Pi }}_{\mathbf{l}}d{\mathsf{\Pi }}_{{\mathbf{l}}^{\prime }}{\rho }_{1,1;1,1}(\mathbf{q};\mathbf{l}|{\mathbf{q}}^{\prime };{\mathbf{l}}^{\prime }\left|0\right)\left(\prod _{a=\varphi ,\phi }{\int }_{{\mathbf{x}}_a{\mathbf{x}}_a^{\prime }{\mathbf{y}}_a{\mathbf{y}}_a^{\prime }}\right)\hfill \\ & & \times e^{-\mathrm{i}(\mathbf{p}·{\mathbf{x}}_{\varphi }-{\mathbf{p}}^{\prime }·{\mathbf{x}}_{\varphi }^{\prime })+\mathrm{i}(\mathbf{q}·{\mathbf{y}}_{\varphi }-{\mathbf{q}}^{\prime }·{\mathbf{y}}_{\varphi }^{\prime })}{e}^{-\mathrm{i}(\mathbf{k}·{\mathbf{x}}_{\phi }-{\mathbf{k}}^{\prime }·{\mathbf{x}}_{\phi }^{\prime })+\mathrm{i}(\mathbf{l}·{\mathbf{y}}_{\phi }-{\mathbf{l}}^{\prime }·{\mathbf{y}}_{\phi }^{\prime })}\hfill \\ & & \times {\partial }_{x_{\varphi }^0,E_{\mathbf{p}}^{\varphi }}{\partial }_{x_{\varphi }^{0\prime },E_{{\mathbf{p}}^{\prime }}^{\varphi }}^{*}{\partial }_{y_{\varphi }^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y_{\varphi }^{0\prime },E_{{\mathbf{q}}^{\prime }}^{\varphi }}{\partial }_{x_{\phi }^0,E_{\mathbf{k}}^{\phi }}{\partial }_{x_{\phi }^{0\prime },E_{{\mathbf{k}}^{\prime }}^{\phi }}^{*}{\partial }_{y_{\phi }^0,E_{\mathbf{l}}^{\phi }}^{*}{\partial }_{y_{\phi }^{0\prime },E_{{\mathbf{l}}^{\prime }}^{\phi }}\hfill \\ & & \times \langle \langle 0\left|\mathrm{T}\right[{\widehat{\varphi }}^{+}\left(x_{\varphi }\right){\widehat{\varphi }}^{-}\left(x_{\varphi }^{\prime }\right){\widehat{\phi }}^{+}\left(x_{\phi }\right){\widehat{\phi }}^{-}\left(x_{\phi }^{\prime }\right)\hfill \\ & & \times exp\left\{-\mathrm{i}\underset{0}{\overset{t}{\int }}[{\widehat{H}}_{\varphi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\phi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\mathrm{int}}\left(\tau \right)]d\tau \right\}\hfill \\ & & \times {\widehat{\varphi }}^{+}\left(y_{\varphi }\right){\widehat{\varphi }}^{-}\left(y_{\varphi }^{\prime }\right){\widehat{\phi }}^{+}\left(y_{\phi }\right){\widehat{\phi }}^{-}\left(y_{\phi }^{\prime }\right)\left]\right|0\rangle \rangle .\hfill \end{array}$$

(36)

As in Ref. [21], we introduced limits in which $x_{\varphi ,\phi }^{0\left({}^{\prime }\right)}$ approach t from above and $y_{\varphi ,\phi }^{0\left({}^{\prime }\right)}$ approach 0 from below in order to recover the correct time ordering.

We translate the expression in Equation (36) into the path integral formalism, which yields

$$\begin{array}{ccc}\hfill {\rho }_{1,1;1,1}(\mathbf{p};\mathbf{k}|{\mathbf{p}}^{\prime };{\mathbf{k}}^{\prime }\left|t\right)& =& \underset{\begin{array}{c}{x}_{\varphi ,\phi }^{0\left({}^{\prime }\right)}\to t^{+}\\ y_{\varphi ,\phi }^{0\left({}^{\prime }\right)}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}d{\mathsf{\Pi }}_{\mathbf{l}}d{\mathsf{\Pi }}_{{\mathbf{l}}^{\prime }}{\rho }_{1,1;1,1}(\mathbf{q};\mathbf{l}|{\mathbf{q}}^{\prime };{\mathbf{l}}^{\prime }\left|0\right)\left(\prod _{a=\varphi ,\phi }{\int }_{{\mathbf{x}}_a{\mathbf{x}}_a^{\prime }{\mathbf{y}}_a{\mathbf{y}}_a^{\prime }}\right)\hfill \\ & & \times e^{-\mathrm{i}(\mathbf{p}·{\mathbf{x}}_{\varphi }-{\mathbf{p}}^{\prime }·{\mathbf{x}}_{\varphi }^{\prime })+\mathrm{i}(\mathbf{q}·{\mathbf{y}}_{\varphi }-{\mathbf{q}}^{\prime }·{\mathbf{y}}_{\varphi }^{\prime })}{e}^{-\mathrm{i}(\mathbf{k}·{\mathbf{x}}_{\phi }-{\mathbf{k}}^{\prime }·{\mathbf{x}}_{\phi }^{\prime })+\mathrm{i}(\mathbf{l}·{\mathbf{y}}_{\phi }-{\mathbf{l}}^{\prime }·{\mathbf{y}}_{\phi }^{\prime })}\hfill \\ & & \times {\partial }_{x_{\varphi }^0,E_{\mathbf{p}}^{\varphi }}{\partial }_{x_{\varphi }^{0\prime },E_{{\mathbf{p}}^{\prime }}^{\varphi }}^{*}{\partial }_{y_{\varphi }^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y_{\varphi }^{0\prime },E_{{\mathbf{q}}^{\prime }}^{\varphi }}{\partial }_{x_{\phi }^0,E_{\mathbf{k}}^{\phi }}{\partial }_{x_{\phi }^{0\prime },E_{{\mathbf{k}}^{\prime }}^{\phi }}^{*}{\partial }_{y_{\phi }^0,E_{\mathbf{l}}^{\phi }}^{*}{\partial }_{y_{\phi }^{0\prime },E_{{\mathbf{l}}^{\prime }}^{\phi }}\hfill \\ & & \times \int \mathcal{D}{\varphi }^{\pm }\mathcal{D}{\phi }^{\pm }{e}^{\mathrm{i}{\widehat{S}}_0[\varphi ;\phi |t]}{\varphi }^{+}\left(x_{\varphi }\right){\varphi }^{-}\left(x_{\varphi }^{\prime }\right){\phi }^{+}\left(x_{\phi }\right){\phi }^{-}\left(x_{\phi }^{\prime }\right)\hfill \\ & & \times exp\left\{\mathrm{i}[{\widehat{S}}_{\varphi ,\mathrm{int}}[\varphi ;t]+{\widehat{S}}_{\phi ,\mathrm{int}}[\phi ;t]+{\widehat{S}}_{\mathrm{int}}[\varphi ;\phi ;t]]\right\}\hfill \\ & & \times {\varphi }^{+}\left(y_{\varphi }\right){\varphi }^{-}\left(y_{\varphi }^{\prime }\right){\phi }^{+}\left(y_{\phi }\right){\phi }^{-}\left(y_{\phi }^{\prime }\right),\hfill \end{array}$$

(37)

where we defined

$$\begin{array}{ccc}\hfill {\widehat{S}}_0[\varphi ;\phi |t]& :=& {\widehat{S}}_{\varphi }\left[\varphi \right]+{\widehat{S}}_{\phi }\left[\phi \right],\hfill \end{array}$$

(38)

and used

$$\begin{array}{ccc}\hfill {\widehat{H}}_{\varphi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\phi ,\mathrm{int}}\left(\tau \right)+{\widehat{H}}_{\mathrm{int}}\left(\tau \right)& =& \widehat{H}\left(\tau \right)-{\widehat{H}}_0\left(\tau \right)\hfill \\ & =& \frac{\partial }{\partial \tau }{\widehat{S}}_0\left(\tau \right)-\frac{\partial }{\partial \tau }\widehat{S}\left(\tau \right)\hfill \\ & =& \frac{\partial }{\partial \tau }{\widehat{S}}_0\left(\tau \right)+\sum _{a=\pm }a({\dot{\varphi }}^a{\pi }_{\varphi }^a+{\dot{\phi }}^a{\pi }_{\phi }^a)\hfill \\ & & -\left[\frac{\partial }{\partial \tau }\widehat{S}\left(\tau \right)+\sum _{a=\pm }a({\dot{\varphi }}^a{\pi }_{\varphi }^a+{\dot{\phi }}^a{\pi }_{\phi }^a)\right]\hfill \\ & =& \frac{d}{d\tau }{\widehat{S}}_0\left(\tau \right)-\frac{d}{d\tau }\widehat{S}\left(\tau \right)\hfill \\ & =& -\frac{d}{d\tau }[{\widehat{S}}_{\varphi ,\mathrm{int}}\left(\tau \right)+{\widehat{S}}_{\phi ,\mathrm{int}}\left(\tau \right)+{\widehat{S}}_{\mathrm{int}}\left(\tau \right)].\hfill \end{array}$$

(39)

The Equation (37) we just found allows us to directly compute the density matrix element that describes the interaction between two single particles $\varphi$ and $\phi$. Taking Equation (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 Equation (37) accordingly. We find:

$$\begin{array}{ccc}\hfill {\rho }_{G;H}(K^G|L^H\left|t\right)& =& \sum _{I,J=0}^{\infty }\frac{{\mathrm{i}}^{G+J}{(-\mathrm{i})}^{H+I}}{I!J!}\underset{\begin{array}{c}{X}^{0,G},X^{0\prime ,H}\to t^{+}\\ Y^{0,I},Y^{0\prime ,J}\to 0^{-}\end{array}}{lim}\int \left(\prod _{A=1}^{I}d{\mathsf{\Pi }}_{R^A}\right)\left(\prod _{B=1}^{J}d{\mathsf{\Pi }}_{P^B}\right){\rho }_{I;J}(R^I|P^J\left|0\right)\hfill \\ & & \times {\int }_{X^G{X}^{\prime H}{Y}^I{Y}^{\prime J}}exp\left\{-\mathrm{i}\left(K^G{X}^G-L^H{X}^{\prime H}\right)+\mathrm{i}\left(R^I{Y}^I-P^J{Y}^{\prime J}\right)\right\}\hfill \\ & & \times \left(\prod _{A=1}^G{\partial }_{X^{0,A},E_{K^A}}\right)\left(\prod _{B=1}^H{\partial }_{X^{0\prime ,B},E_{L^B}}^{*}\right)\left(\prod _{C=1}^I{\partial }_{Y^{0,C},E_{R^C}}^{*}\right)\left(\prod _{D=1}^J{\partial }_{Y^{0\prime ,D},E_{P^D}}\right)\hfill \\ & & \times \int \mathcal{D}{\mathsf{\Phi }}^{\pm }{e}^{\mathrm{i}{\widehat{S}}_{\mathsf{\Phi }}\left[\mathsf{\Phi }\right]}{\mathsf{\Phi }}_{X^G}^{+}{\mathsf{\Phi }}_{X^{\prime H}}^{-}exp\left\{\mathrm{i}[{\widehat{S}}_{\mathsf{\Phi },\mathrm{int}}[\mathsf{\Phi };t]+{\widehat{S}}_{\mathrm{int}}[\mathsf{\Phi };t]]\right\}{\mathsf{\Phi }}_{Y^I}^{+}{\mathsf{\Phi }}_{Y^{\prime J}}^{-}.\hfill \end{array}$$

(40)

In Equation (40), we used some short-hand notations from Section 2.1 and introduced a number new ones. More precisely, we are using

$$\begin{array}{ccc}\hfill {\mathrm{i}}^G& :=& {\mathrm{i}}^{g_{\alpha }+g_{\beta }+\cdots },\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill X^{0,G}& :=& x_{\alpha ,\left(1\right)}^0,\cdots ,x_{\alpha ,\left(g_{\alpha }\right)}^0;x_{\beta ,\left(1\right)}^0,\cdots ,x_{\beta ,\left(g_{\beta }\right)}^0;\cdots ,X^G:={\mathbf{x}}_{\alpha ,\left(1\right)},\cdots ,{\mathbf{x}}_{\alpha ,\left(g_{\alpha }\right)};{\mathbf{x}}_{\beta ,\left(1\right)},\cdots ,{\mathbf{x}}_{\beta ,\left(g_{\beta }\right)};\cdots ,\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill K^G{X}^G& :=& {\mathbf{k}}_{\alpha }^{\left(1\right)}{\mathbf{x}}_{\alpha ,\left(1\right)}+\cdots +{\mathbf{k}}_{\alpha }^{\left(g_{\alpha }\right)}{\mathbf{x}}_{\alpha ,\left(g_{\alpha }\right)}+{\mathbf{k}}_{\beta }^{\left(1\right)}{\mathbf{x}}_{\beta ,\left(1\right)}+\cdots +{\mathbf{k}}_{\beta }^{\left(g_{\beta }\right)}{\mathbf{x}}_{\beta ,\left(g_{\beta }\right)}+\cdots ,\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill \prod _{A=1}^G{\partial }_{X^{0,A},E_{K^A}}& :=& {\partial }_{x_{\alpha ,\left(1\right)}^0,E_{{\mathbf{k}}_{\alpha }^{\left(1\right)}}^{\alpha }}\cdots {\partial }_{x_{\alpha ,\left(g_{\alpha }\right)}^0,E_{{\mathbf{k}}_{\alpha }^{\left(g_{\alpha }\right)}}^{\alpha }}{\partial }_{x_{\beta ,\left(1\right)}^0,E_{{\mathbf{k}}_{\beta }^{\left(1\right)}}^{\beta }}\cdots {\partial }_{x_{\beta ,\left(g_{\beta }\right)}^0,E_{{\mathbf{k}}_{\beta }^{\left(g_{\beta }\right)}}^{\beta }}\cdots ,\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill \mathcal{D}{\mathsf{\Phi }}^{\pm }& :=& \mathcal{D}{\alpha }^{\pm }\mathcal{D}{\beta }^{\pm }\cdots ,\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill {\widehat{S}}_{\mathsf{\Phi },\mathrm{int}}[\mathsf{\Phi };t]& :=& {\widehat{S}}_{\alpha ,\mathrm{int}}[\alpha ;t]+{\widehat{S}}_{\beta ,\mathrm{int}}[\beta ;t]+\cdots ,{\widehat{S}}_{\mathrm{int}}[\mathsf{\Phi };t]={\widehat{S}}_{\mathrm{int}}[\alpha ;\beta ;\cdots ;t]+\cdots \hfill \end{array}$$

(46)

and

$$\begin{array}{ccc}\hfill {\mathsf{\Phi }}_{X^G}^{+}& :=& {\alpha }_{x_{\alpha ,\left(1\right)}}^{+}\cdots {\alpha }_{x_{\alpha ,\left(g_{\alpha }\right)}}^{+}{\beta }_{x_{\beta ,\left(1\right)}}^{+}\cdots {\beta }_{x_{\beta ,\left(g_{\beta }\right)}}^{+}\cdots \hfill \end{array}$$

(47)

with

$$\begin{array}{ccc}\hfill {\alpha }_x& :=& \alpha \left(x\right).\hfill \end{array}$$

(48)

Taking $G=H=(1,1)$ and assuming all initial density matrix elements except for the $I=J=(1,1)$ one to be vanishing, we recover Equation (37) from Equation (40).

It is important to remember that, as was already pointed out in Ref. [21], only the diagonal elements of the $2\times 2$-matrix propagator are permitted when evaluating Equation (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. Figure 1) are possible, such that we only have to work with Feynman and Dyson propagators

$$\begin{array}{ccc}\hfill \langle \langle 0|\mathrm{T}\left[{\widehat{\varphi }}_x^{+}{\widehat{\varphi }}_y^{+}\right]|0\rangle \rangle & =& D_{xy}^{\mathrm{F}}=-\mathrm{i}{\int }_k\frac{e^{\mathrm{i}k·(x-y)}}{k^2+M^2-iϵ},\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill \langle \langle 0|\mathrm{T}\left[{\widehat{\varphi }}_x^{-}{\widehat{\varphi }}_y^{-}\right]|0\rangle \rangle & =& D_{xy}^{\mathrm{D}}=+\mathrm{i}{\int }_k\frac{e^{\mathrm{i}k·(x-y)}}{k^2+M^2+iϵ},\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill \langle \langle 0|\mathrm{T}\left[{\widehat{\phi }}_x^{+}{\widehat{\phi }}_y^{+}\right]|0\rangle \rangle & =& {\Delta }_{xy}^{\mathrm{F}}=-\mathrm{i}{\int }_k\frac{e^{\mathrm{i}k·(x-y)}}{k^2+m^2-iϵ},\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill \langle \langle 0|\mathrm{T}\left[{\widehat{\phi }}_x^{-}{\widehat{\phi }}_y^{-}\right]|0\rangle \rangle & =& {\Delta }_{xy}^{\mathrm{D}}=+\mathrm{i}{\int }_k\frac{e^{\mathrm{i}k·(x-y)}}{k^2+m^2+iϵ},\hfill \end{array}$$

(52)

while the other two-point functions vanish:

$$\begin{array}{ccc}\hfill \langle \langle 0|{\widehat{\varphi }}_x^{+(-)}{\widehat{\varphi }}_y^{-(+)}|0\rangle \rangle & =& \langle \langle 0|{\widehat{\phi }}_x^{+(-)}{\widehat{\phi }}_y^{-(+)}|0\rangle \rangle =0.\hfill \end{array}$$

(53)

3. Example

In this section, we will apply the formulas in Equations (37) and (40) to a toy model example with two real scalar field species $\varphi$ and $\chi$. More precisely, we will consider solely a single particle $\varphi$ at the initial time 0, which is described by the only non-vanishing density matrix ${\rho }_{1,0;1,0}\left(0\right)$. We will assume that due to an interaction potential $V\sim {\phi }^2\varphi$ the state of two $\phi$-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 ${\rho }_{0,2;0,2}\left(t\right)$, ${\rho }_{1,0;0,2}\left(t\right)$ and ${\rho }_{1,0;1,0}\left(t\right)$. Each of these matrices has its own physical interpretation: the first one corresponds to a $1\to 2$-particle decay event, the second describes the correlation between a single $\varphi$-particle state and a two $\phi$-particles state, while the third stands for the single $\varphi$-particle state continuing to exist, but with a modified probability. The last of these events, i.e., the change of ${\rho }_{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

$$\begin{array}{ccc}\hfill S_{\varphi }\left[\varphi \right]& =& {\int }_x\left[-\frac{1}{2}{(\partial \varphi )}^2-\frac{1}{2}{M}^2{\varphi }^2\right],\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill S_{\phi }\left[\phi \right]& =& {\int }_x\left[-\frac{1}{2}{(\partial \phi )}^2-\frac{1}{2}{m}^2{\phi }^2\right]\hfill \end{array}$$

(55)

for the scalar fields $\varphi$ and $\phi$. As interaction action we choose to use

$$\begin{array}{ccc}\hfill S_{\mathrm{int}}[\varphi ;\phi ]& =& {\int }_{x\in {\mathsf{\Omega }}_t}\left[-\frac{\alpha }{2}\mathcal{M}{\phi }^2\varphi \right],\hfill \end{array}$$

(56)

where $\alpha \ll 1$ and $\mathcal{M}$ are some sufficiently small mass scales. For simplicity, we do not consider self-interactions. 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 Equation (56). We decide to work up to second order in $\alpha$, such that we find:

$$\begin{array}{ccc}\hfill exp\left\{\mathrm{i}{\widehat{S}}_{\mathrm{int}}[\varphi ;\phi ;t]\right\}& =& 1-\mathrm{i}\frac{\alpha }{2}\mathcal{M}\sum _{a=\pm }a{\int }_z{\left({\phi }_z^a\right)}^2{\varphi }_z^a\hfill \\ & & -\frac{{\alpha }^2}{8}{\mathcal{M}}^2\sum _{a,b=\pm }ab{\int }_{z{z}^{\prime }}{\left({\phi }_z^a\right)}^2{\left({\phi }_{z^{\prime }}^b\right)}^2{\varphi }_z^a{\varphi }_{z^{\prime }}^b+\mathcal{O}\left({\alpha }^3\right),\hfill \end{array}$$

(57)

where in this and all following equations z and $z^{\prime }$ are integrated over the domain ${\mathsf{\Omega }}_t$. We begin our investigations with the decay of a single $\varphi$-particle into two $\phi$-particles, and compute the density matrix elements for this process. Obviously, the only non-vanishing density matrix at time 0 must be ${\rho }_{1,0;1,0}$, while we want to find the matrix ${\rho }_{0,2;0,2}$ at time t. Following Equation (40), we must evaluate

$$\begin{array}{ccc}\hfill {\rho }_{0,2;0,2}(;\mathbf{p},\mathbf{k}|;{\mathbf{p}}^{\prime },{\mathbf{k}}^{\prime }\left|t\right)& =& \underset{\begin{array}{c}{x}_{\left(1\right),\left(2\right)}^{0{(}^{\prime })}\to t^{+}\\ y^{0{(}^{\prime })}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}{\rho }_{1,0;1,0}(\mathbf{q};|{\mathbf{q}}^{\prime };\left|0\right)\hfill \\ & & \times {\int }_{{\mathbf{x}}_{\left(1\right)}{\mathbf{x}}_{\left(1\right)}^{\prime }{\mathbf{x}}_{\left(2\right)}{\mathbf{x}}_{\left(2\right)}^{\prime }\mathbf{y}{\mathbf{y}}^{\prime }}{e}^{-\mathrm{i}(\mathbf{p}·{\mathbf{x}}_{\left(1\right)}+\mathbf{k}·{\mathbf{x}}_{\left(2\right)}-{\mathbf{p}}^{\prime }·{\mathbf{x}}_{\left(1\right)}^{\prime }-{\mathbf{k}}^{\prime }·{\mathbf{x}}_{\left(2\right)}^{\prime })+\mathrm{i}(\mathbf{q}·\mathbf{y}-{\mathbf{q}}^{\prime }·{\mathbf{y}}^{\prime })}\hfill \\ & & \times {\partial }_{x_{\left(1\right)}^0,E_{\mathbf{p}}^{\phi }}{\partial }_{x_{\left(1\right)}^{0^{\prime }},E_{{\mathbf{p}}^{\prime }}^{\phi }}^{*}{\partial }_{x_{\left(2\right)}^0,E_{\mathbf{k}}^{\phi }}{\partial }_{x_{\left(2\right)}^{0^{\prime }},E_{{\mathbf{k}}^{\prime }}^{\phi }}^{*}{\partial }_{y^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y^{0^{\prime }},E_{{\mathbf{q}}^{\prime }}^{\varphi }}\hfill \\ & & \times \int \mathcal{D}{\varphi }^{\pm }\mathcal{D}{\phi }^{\pm }{e}^{\mathrm{i}{\widehat{S}}_{\varphi }\left[\varphi \right]+\mathrm{i}{\widehat{S}}_{\phi }\left[\phi \right]}{\phi }_{x_{\left(1\right)}}^{+}{\phi }_{x_{\left(1\right)}^{\prime }}^{-}{\phi }_{x_{\left(2\right)}}^{+}{\phi }_{x_{\left(2\right)}^{\prime }}^{-}\hfill \\ & & \times exp\left\{\mathrm{i}{\widehat{S}}_{\mathrm{int}}[\varphi ;\phi ;t]\right\}{\varphi }_y^{+}{\varphi }_{y^{\prime }}^{-}.\hfill \end{array}$$

(58)

Substituting Equation (57), we are left with

$$\begin{array}{ccc}\hfill {\rho }_{0,2;0,2}(;\mathbf{p},\mathbf{k}|;{\mathbf{p}}^{\prime },{\mathbf{k}}^{\prime }\left|t\right)& \approx & \frac{{\alpha }^2{\mathcal{M}}^2}{4}\underset{\begin{array}{c}{x}_{\left(1\right),\left(2\right)}^{0{(}^{\prime })}\to t^{+}\\ y^{0{(}^{\prime })}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}{\rho }_{1,0;1,0}(\mathbf{q};|{\mathbf{q}}^{\prime };\left|0\right)\hfill \\ & & \times {\int }_{{\mathbf{x}}_{\left(1\right)}{\mathbf{x}}_{\left(1\right)}^{\prime }{\mathbf{x}}_{\left(2\right)}{\mathbf{x}}_{\left(2\right)}^{\prime }\mathbf{y}{\mathbf{y}}^{\prime }}{e}^{-\mathrm{i}(\mathbf{p}·{\mathbf{x}}_{\left(1\right)}+\mathbf{k}·{\mathbf{x}}_{\left(2\right)}-{\mathbf{p}}^{\prime }·{\mathbf{x}}_{\left(1\right)}^{\prime }-{\mathbf{k}}^{\prime }·{\mathbf{x}}_{\left(2\right)}^{\prime })+\mathrm{i}(\mathbf{q}·\mathbf{y}-{\mathbf{q}}^{\prime }·{\mathbf{y}}^{\prime })}\hfill \\ & & \times {\partial }_{x_{\left(1\right)}^0,E_{\mathbf{p}}^{\phi }}{\partial }_{x_{\left(1\right)}^{0^{\prime }},E_{{\mathbf{p}}^{\prime }}^{\phi }}^{*}{\partial }_{x_{\left(2\right)}^0,E_{\mathbf{k}}^{\phi }}{\partial }_{x_{\left(2\right)}^{0^{\prime }},E_{{\mathbf{k}}^{\prime }}^{\phi }}^{*}{\partial }_{y^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y^{0^{\prime }},E_{{\mathbf{q}}^{\prime }}^{\varphi }}\hfill \\ & & \times {\int }_{z{z}^{\prime }}{D}_{zy}^{\mathrm{F}}{D}_{z^{\prime }{y}^{\prime }}^{\mathrm{D}}[{\Delta }_{x_{\left(1\right)}{x}_{\left(2\right)}}^{\mathrm{F}}{\Delta }_{zz}^{\mathrm{F}}{\Delta }_{x_{\left(1\right)}^{\prime }{x}_{\left(2\right)}^{\prime }}^{\mathrm{D}}{\Delta }_{z^{\prime }{z}^{\prime }}^{\mathrm{D}}+2{\Delta }_{x_{\left(1\right)}{x}_{\left(2\right)}}^{\mathrm{F}}{\Delta }_{zz}^{\mathrm{F}}{\Delta }_{x_{\left(1\right)}^{\prime }{z}^{\prime }}^{\mathrm{D}}{\Delta }_{x_{\left(2\right)}^{\prime }{z}^{\prime }}^{\mathrm{D}}\hfill \\ & & +2{\Delta }_{x_{\left(1\right)}z}^{\mathrm{F}}{\Delta }_{x_{\left(2\right)}z}^{\mathrm{F}}{\Delta }_{x_{\left(1\right)}^{\prime }{x}_{\left(2\right)}^{\prime }}^{\mathrm{D}}{\Delta }_{z^{\prime }{z}^{\prime }}^{\mathrm{D}}+4{\Delta }_{x_{\left(1\right)}z}^{\mathrm{F}}{\Delta }_{x_{\left(2\right)}z}^{\mathrm{F}}{\Delta }_{x_{\left(1\right)}^{\prime }{z}^{\prime }}^{\mathrm{D}}{\Delta }_{x_{\left(2\right)}^{\prime }{z}^{\prime }}^{\mathrm{D}}].\hfill \end{array}$$

(59)

Using the propagators given in this expression, we can easily depict the contributing physical processes in Figure 2a–d. We observe that one (Figure 2b,c) or two (Figure 2a) of the $\varphi$-propagators end in $\phi$-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 Figure 2d is finite and corresponds to the actual decay of a single $\varphi$-particle into two $\phi$-particles without the appearance of any unphysical terms.

Of course, Equation (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:

$$\begin{array}{ccc}\hfill {\rho }_{0,2;0,2}(;\mathbf{p},\mathbf{k}|;{\mathbf{p}}^{\prime },{\mathbf{k}}^{\prime }\left|t\right)& \approx & \frac{{\alpha }^2{\mathcal{M}}^2}{4}\{\frac{{\rho }_{1,0;1,0}(\mathbf{0};|\mathbf{0};\left|0\right)}{M^4}{\delta }^2\left(0\right){\Delta }_{zz}^{\mathrm{F}}{\Delta }_{z^{\prime }{z}^{\prime }}^{\mathrm{D}}\hfill \\ & & \times {\left(2\pi \right)}^6{\delta }^{\left(3\right)}(\mathbf{p}+\mathbf{k}){\delta }^{\left(3\right)}({\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime }){sin}^2(Mt/2)\hfill \\ & & -\mathrm{i}\frac{{\rho }_{1,0;1,0}(\mathbf{0};|{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime };\left|0\right)}{2M^2{E}_{{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime }}^{\varphi }(E_{{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime }}^{\varphi }-E_{{\mathbf{p}}^{\prime }}^{\phi }-E_{{\mathbf{k}}^{\prime }}^{\phi })}\delta \left(0\right){\Delta }_{zz}^{\mathrm{F}}{\left(2\pi \right)}^3{\delta }^{\left(3\right)}(\mathbf{p}+\mathbf{k})\hfill \\ & & \times \left[1-e^{-\mathrm{i}Mt}\right]\left[e^{\mathrm{i}(E_{{\mathbf{p}}^{\prime }}^{\phi }+E_{{\mathbf{k}}^{\prime }}^{\phi })t}-e^{\mathrm{i}{E}_{{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime }}^{\varphi }t}\right]\hfill \\ & & +\mathrm{i}\frac{{\rho }_{1,0;1,0}(\mathbf{p}+\mathbf{k};|\mathbf{0};\left|0\right)}{2M^2{E}_{\mathbf{p}+\mathbf{k}}^{\varphi }(E_{\mathbf{p}+\mathbf{k}}^{\varphi }-E_{\mathbf{p}}^{\phi }-E_{\mathbf{k}}^{\phi })}\delta \left(0\right){\Delta }_{z^{\prime }{z}^{\prime }}^{\mathrm{D}}{\left(2\pi \right)}^3{\delta }^{\left(3\right)}({\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime })\hfill \\ & & \times \left[1-e^{\mathrm{i}Mt}\right]\left[e^{-\mathrm{i}(E_{\mathbf{p}}^{\phi }+E_{\mathbf{k}}^{\phi })t}-e^{-\mathrm{i}{E}_{\mathbf{p}+\mathbf{k}}^{\varphi }t}\right]\hfill \\ & & +\frac{{\rho }_{1,0;1,0}(\mathbf{p}+\mathbf{k};|{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime };\left|0\right)}{E_{\mathbf{p}+\mathbf{k}}^{\varphi }{E}_{{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime }}^{\varphi }(E_{\mathbf{p}+\mathbf{k}}^{\varphi }-E_{\mathbf{p}}^{\phi }-E_{\mathbf{k}}^{\phi })(E_{{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime }}^{\varphi }-E_{{\mathbf{p}}^{\prime }}^{\phi }-E_{{\mathbf{k}}^{\prime }}^{\phi })}\hfill \\ & & \times \left[e^{-\mathrm{i}(E_{\mathbf{p}}^{\phi }+E_{\mathbf{k}}^{\phi })t}-e^{-\mathrm{i}{E}_{\mathbf{p}+\mathbf{k}}^{\varphi }t}\right]\left[e^{\mathrm{i}(E_{{\mathbf{p}}^{\prime }}^{\phi }+E_{{\mathbf{k}}^{\prime }}^{\phi })t}-e^{\mathrm{i}{E}_{{\mathbf{p}}^{\prime }+{\mathbf{k}}^{\prime }}^{\varphi }t}\right]\}.\hfill \end{array}$$

(60)

We observe that, as expected, the terms in the fourth, sixth and eighth lines of Equation (60) contain differences in the unitary evolutions of a single $\varphi$-particle and two $\phi$-particles.

Next, we compute one of the two density matrices that describe correlations between a single $\varphi$-particle and two $\phi$-particles. The other one can easily be obtained by making use of the property given in Equation (9). Starting with Equation (40), we find:

$$\begin{array}{ccc}\hfill {\rho }_{1,0;0,2}(\mathbf{p};|;{\mathbf{k}}^{\prime },{\mathbf{r}}^{\prime }\left|t\right)& =& -\mathrm{i}\underset{\begin{array}{c}{x}_{\varphi }^0,x_{\phi ,\left(1\right)}^{0^{\prime }},x_{\phi ,\left(2\right)}^{0^{\prime }}\to t^{+}\\ y^{0\left({}^{\prime }\right)}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}{\rho }_{1,0;1,0}(\mathbf{q};|{\mathbf{q}}^{\prime };\left|0\right)\hfill \\ & & \times {\int }_{{\mathbf{x}}_{\varphi }{\mathbf{x}}_{\phi ,\left(1\right)}^{\prime }{\mathbf{x}}_{\phi ,\left(2\right)}^{\prime }\mathbf{y}{\mathbf{y}}^{\prime }}{e}^{-\mathrm{i}\mathbf{p}·{\mathbf{x}}_{\varphi }+\mathrm{i}(\mathbf{q}·\mathbf{y}-{\mathbf{q}}^{\prime }·{\mathbf{y}}^{\prime })}{e}^{\mathrm{i}({\mathbf{k}}^{\prime }·{\mathbf{x}}_{\phi ,\left(1\right)}^{\prime }+{\mathbf{r}}^{\prime }·{\mathbf{x}}_{\phi ,\left(2\right)}^{\prime })}\hfill \\ & & \times {\partial }_{x_{\varphi }^0,E_{\mathbf{p}}^{\varphi }}{\partial }_{y^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y^{0\prime },E_{{\mathbf{q}}^{\prime }}^{\varphi }}{\partial }_{x_{\phi ,\left(1\right)}^{0\prime },E_{{\mathbf{k}}^{\prime }}^{\phi }}^{*}{\partial }_{x_{\phi ,\left(2\right)}^{0\prime },E_{{\mathbf{r}}^{\prime }}^{\phi }}^{*}\hfill \\ & & \times \int \mathcal{D}{\varphi }^{\pm }\mathcal{D}{\phi }^{\pm }{e}^{\mathrm{i}{\widehat{S}}_{\varphi }\left[\varphi \right]+\mathrm{i}{\widehat{S}}_{\phi }\left[\phi \right]}{\varphi }_{x_{\varphi }}^{+}{\phi }_{x_{\phi ,\left(1\right)}^{\prime }}^{-}{\phi }_{x_{\phi ,\left(2\right)}^{\prime }}^{-}\hfill \\ & & \times exp\left\{\mathrm{i}{\widehat{S}}_{\mathrm{int}}[\varphi ;\phi ;t]\right\}{\varphi }_y^{+}{\varphi }_{y^{\prime }}^{-}.\hfill \end{array}$$

(61)

Evaluating the path integrals leads us to

$$\begin{array}{ccc}\hfill {\rho }_{1,0;0,2}(\mathbf{p};|;{\mathbf{k}}^{\prime },{\mathbf{r}}^{\prime }\left|t\right)& \approx & \frac{\alpha \mathcal{M}}{2}\underset{\begin{array}{c}{x}_{\varphi }^0,x_{\phi ,\left(1\right)}^{0^{\prime }},x_{\phi ,\left(2\right)}^{0^{\prime }}\to t^{+}\\ y^{0\left({}^{\prime }\right)}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}{\rho }_{1,0;1,0}(\mathbf{q};|{\mathbf{q}}^{\prime };\left|0\right)\hfill \\ & & \times {\int }_{{\mathbf{x}}_{\varphi }{\mathbf{x}}_{\phi ,\left(1\right)}^{\prime }{\mathbf{x}}_{\phi ,\left(2\right)}^{\prime }\mathbf{y}{\mathbf{y}}^{\prime }}{e}^{-\mathrm{i}\mathbf{p}·{\mathbf{x}}_{\varphi }+\mathrm{i}(\mathbf{q}·\mathbf{y}-{\mathbf{q}}^{\prime }·{\mathbf{y}}^{\prime })}{e}^{\mathrm{i}({\mathbf{k}}^{\prime }·{\mathbf{x}}_{\phi ,\left(1\right)}^{\prime }+{\mathbf{r}}^{\prime }·{\mathbf{x}}_{\phi ,\left(2\right)}^{\prime })}\hfill \\ & & \times {\partial }_{x_{\varphi }^0,E_{\mathbf{p}}^{\varphi }}{\partial }_{y^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y^{0\prime },E_{{\mathbf{q}}^{\prime }}^{\varphi }}{\partial }_{x_{\phi ,\left(1\right)}^{0\prime },E_{{\mathbf{k}}^{\prime }}^{\phi }}^{*}{\partial }_{x_{\phi ,\left(2\right)}^{0\prime },E_{{\mathbf{r}}^{\prime }}^{\phi }}^{*}\hfill \\ & & \times {\int }_z{D}_{x_{\varphi }y}^{\mathrm{F}}{D}_{z{y}^{\prime }}^{\mathrm{D}}\left[{\Delta }_{x_{\phi ,\left(1\right)}^{\prime }{x}_{\phi ,\left(2\right)}^{\prime }}^{\mathrm{D}}{\Delta }_{zz}^{\mathrm{D}}+2{\Delta }_{x_{\phi ,\left(1\right)}^{\prime }z}^{\mathrm{D}}{\Delta }_{x_{\phi ,\left(2\right)}^{\prime }z}^{\mathrm{D}}\right],\hfill \end{array}$$

(62)

which can be used to depict the processes as in Figure 3a,b. As in the previous case, we observe one of the $\varphi$-propagators ending in a $\phi$-tadpole, while an equal-time $\phi$-propagator is appearing, in Figure 3a. In contrast, Figure 3b contains only physical terms without any divergences.

Continuing the computation, we find:

$$\begin{array}{ccc}\hfill {\rho }_{1,0;0,2}(\mathbf{p};|;{\mathbf{k}}^{\prime },{\mathbf{r}}^{\prime }\left|t\right)& \approx & \frac{\alpha \mathcal{M}{e}^{-\mathrm{i}{E}_{\mathbf{p}}^{\varphi }t}}{2}\{-\mathrm{i}\delta \left(0\right){\Delta }_{zz}^{\mathrm{D}}{\left(2\pi \right)}^3{\delta }^{\left(3\right)}({\mathbf{k}}^{\prime }+{\mathbf{r}}^{\prime })\frac{{\rho }_{1,0;1,0}(\mathbf{p};|\mathbf{0};\left|0\right)}{2M^2}\left[1-e^{\mathrm{i}Mt}\right]\hfill \\ & & +\frac{{\rho }_{1,0;1,0}(\mathbf{p};|{\mathbf{k}}^{\prime }+{\mathbf{r}}^{\prime };\left|0\right)}{E_{{\mathbf{k}}^{\prime }+{\mathbf{r}}^{\prime }}^{\varphi }(E_{{\mathbf{k}}^{\prime }+{\mathbf{r}}^{\prime }}^{\varphi }-E_{{\mathbf{k}}^{\prime }}^{\phi }-E_{{\mathbf{r}}^{\prime }}^{\phi })}\left[e^{\mathrm{i}{E}_{{\mathbf{k}}^{\prime }+{\mathbf{r}}^{\prime }}^{\varphi }t}-e^{\mathrm{i}(E_{{\mathbf{k}}^{\prime }}^{\phi }+E_{{\mathbf{r}}^{\prime }}^{\phi })t}\right]\},\hfill \end{array}$$

(63)

where, in the last line, we again find the difference in the unitary evolutions of a single $\varphi$-particle and two $\phi$-particles.

Finally, we study the change of ${\rho }_{1,0;1,0}$ over time. We use Equation (40) to find:

$$\begin{array}{ccc}\hfill {\rho }_{1,0;1,0}(\mathbf{p};|{\mathbf{p}}^{\prime };\left|t\right)& =& \underset{\begin{array}{c}{x}^{0\left({}^{\prime }\right)}\to t^{+}\\ y^{0\left({}^{\prime }\right)}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}{\rho }_{1,0;1,0}(\mathbf{q};|{\mathbf{q}}^{\prime };\left|0\right){\int }_{\mathbf{x}{\mathbf{x}}^{\prime }\mathbf{y}{\mathbf{y}}^{\prime }}{e}^{-\mathrm{i}(\mathbf{p}·\mathbf{x}-{\mathbf{p}}^{\prime }·{\mathbf{x}}^{\prime })+\mathrm{i}(\mathbf{q}·\mathbf{y}-{\mathbf{q}}^{\prime }·{\mathbf{y}}^{\prime })}\hfill \\ & & \times {\partial }_{x^0,E_{\mathbf{p}}^{\varphi }}{\partial }_{x^{0^{\prime }},E_{{\mathbf{p}}^{\prime }}^{\varphi }}^{*}{\partial }_{y^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y^{0\prime },E_{{\mathbf{q}}^{\prime }}^{\varphi }}\int \mathcal{D}{\varphi }^{\pm }\mathcal{D}{\phi }^{\pm }{e}^{\mathrm{i}{\widehat{S}}_{\varphi }\left[\varphi \right]+\mathrm{i}{\widehat{S}}_{\phi }\left[\phi \right]}{\varphi }_x^{+}{\varphi }_{x^{\prime }}^{-}\hfill \\ & & \times exp\left\{\mathrm{i}{\widehat{S}}_{\mathrm{int}}[\varphi ;\phi ;t]\right\}{\varphi }_y^{+}{\varphi }_{y^{\prime }}^{-},\hfill \end{array}$$

(64)

which becomes

$$\begin{array}{ccc}\hfill {\rho }_{1,0;1,0}(\mathbf{p};|{\mathbf{p}}^{\prime };\left|t\right)& \approx & \underset{\begin{array}{c}{x}^{0\left({}^{\prime }\right)}\to t^{+}\\ y^{0\left({}^{\prime }\right)}\to 0^{-}\end{array}}{lim}\int d{\mathsf{\Pi }}_{\mathbf{q}}d{\mathsf{\Pi }}_{{\mathbf{q}}^{\prime }}{\rho }_{1,0;1,0}(\mathbf{q};|{\mathbf{q}}^{\prime };\left|0\right){\int }_{\mathbf{x}{\mathbf{x}}^{\prime }\mathbf{y}{\mathbf{y}}^{\prime }}{e}^{-\mathrm{i}(\mathbf{p}·\mathbf{x}-{\mathbf{p}}^{\prime }·{\mathbf{x}}^{\prime })+\mathrm{i}(\mathbf{q}·\mathbf{y}-{\mathbf{q}}^{\prime }·{\mathbf{y}}^{\prime })}\hfill \\ & & \times {\partial }_{x^0,E_{\mathbf{p}}^{\varphi }}{\partial }_{x^{0^{\prime }},E_{{\mathbf{p}}^{\prime }}^{\varphi }}^{*}{\partial }_{y^0,E_{\mathbf{q}}^{\varphi }}^{*}{\partial }_{y^{0\prime },E_{{\mathbf{q}}^{\prime }}^{\varphi }}\{[1+\aleph \left(t\right)]D_{xy}^{\mathrm{F}}{D}_{x^{\prime }{y}^{\prime }}^{\mathrm{D}}\hfill \\ & & -\frac{{\alpha }^2}{4}{\mathcal{M}}^2{\int }_{z{z}^{\prime }}[\left({\Delta }_{zz}^{\mathrm{F}}{\Delta }_{z^{\prime }{z}^{\prime }}^{\mathrm{F}}+2{\left({\Delta }_{z{z}^{\prime }}^{\mathrm{D}}\right)}^2\right)D_{xy}^{\mathrm{F}}{D}_{x^{\prime }z}^{\mathrm{D}}{D}_{z^{\prime }{y}^{\prime }}^{\mathrm{D}}\hfill \\ & & +{(x,y⟷x^{\prime },y^{\prime })}^{*}\left]\right\}\hfill \end{array}$$

(65)

after evaluating the path integrals. Here, we introduced

$$\begin{array}{ccc}\hfill \aleph \left(t\right)& :=& -\frac{{\alpha }^2}{8}{\mathcal{M}}^2\sum _{a=\mathrm{F},\mathrm{D}}{\int }_{z{z}^{\prime }}\left[{\Delta }_{zz}^a{\Delta }_{z^{\prime }{z}^{\prime }}^a+2{\left({\Delta }_{z{z}^{\prime }}^a\right)}^2\right]D_{z{z}^{\prime }}^a\hfill \\ & \approx & -\frac{{\alpha }^2{\mathcal{M}}^2}{2}{\left(2\pi \right)}^3{\delta }^{\left(3\right)}\left(\mathbf{0}\right)\{\frac{{sin}^2(Mt/2)}{M^3}{\Delta }_{zz}^{\mathrm{F}}{\Delta }_{z^{\prime }{z}^{\prime }}^{\mathrm{F}}\hfill \\ & & +\frac{1}{2}{\int }_{\mathbf{q}\mathbf{l}}\frac{{sin}^2[(E_{\mathbf{q}}^{\varphi }+E_{\mathbf{l}}^{\phi }+E_{\mathbf{q}+\mathbf{l}}^{\phi })t/2]}{E_{\mathbf{q}}^{\varphi }{E}_{\mathbf{l}}^{\phi }{E}_{\mathbf{q}+\mathbf{l}}^{\phi }{(E_{\mathbf{q}}^{\varphi }+E_{\mathbf{l}}^{\phi }+E_{\mathbf{q}+\mathbf{l}}^{\phi })}^2}\},\hfill \end{array}$$

(66)

which only contains disconnected loop diagrams. Interestingly, at this stage, Equation (65) is exactly the same as the equation obtained for the time evolution of the single-particle reduced density matrix elements in Ref. [21]. Therefore, we point to the results and figures presented in this reference. The reasons for the formalism presented here giving the same results as the one in Ref. [21] lie in us only working up to second order in $\alpha$ and in our choice of the interaction and self-interaction potentials. Even though the formalism in Ref. [21] generally permits contractions for some +-fields with −-fields (the ones representing the traced-out environmental degrees of freedom), the usage of the interaction potential $V\sim {\phi }^2\varphi$ (Ref. [21] names the $\phi$-field $\chi$ instead) does not allow for them to appear in this particular computation to order ${\alpha }^2$. It is important to note that considering higher orders and/or including self-interactions would lead, in addition to Feynman and Dyson propagators, to the appearance of Wightman propagators for the $\chi$-fields in the computation in 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.

4. 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,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 Section 2.1, the Schwinger–Keldysh formalism, Section 2.2, and thermo-field dynamics, Section 2.3. Those were the ingredients required for the actual derivation, which we concluded in Section 2.4. Finally, in Section 3, 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 of low-energy but high-precision phenomenology as, e.g., in quantum optics experiments.