# A new method for directly computing reduced density matrices

Christian Käding<EMAIL_ADDRESS>Technische Universität Wien,
Atominstitut, Stadionallee 2, 1020 Vienna, Austria  Mario Pitschmann
<EMAIL_ADDRESS>Technische Universität Wien, Atominstitut,
Stadionallee 2, 1020 Vienna, Austria

###### Abstract

We demonstrate the power of a first principle-based and practicable method
that allows for the perturbative computation of reduced density matrix
elements of an open quantum system without making use of any master equations.
The approach is based on techniques from non-equilibrium quantum field theory
like thermo field dynamics, the Schwinger-Keldsyh formalism, and the Feynman-
Vernon influence functional. It does not require the Markov approximation and
is essentially a Lehmann-Szymanzik-Zimmermann-like reduction. In order to
illustrate this method, we consider a real scalar field as an open quantum
system interacting with an environment comprising another real scalar field.
We give a general formula that allows for the perturbative computation of
density matrix elements for any number of particles in a momentum basis.
Finally, we consider a simple toy model and use this formula to obtain
expressions for some of the system’s reduced density matrix elements.

open quantum dynamics, density matrix, non-equilibrium quantum field theory

##  I Introduction

Most realistic quantum systems can be regarded as open, which means that they
are surrounded by a number of (uncontrollable) environments whose effects on
the system cannot entirely be neglected [1]. Interactions between a system and
its environments can lead to phenomena like momentum or energy diffusion [2],
which may induce open quantum dynamical effects like phase shifts or
decoherence [3] if the system is in a superposed state. Applications of open
quantum systems and their associated phenomena can be found in a variety of
areas in non-relativistc quantum physics, see e.g. Ref. [4, 5, 6, 7, 8, 9,
10], as well as in quantum field theory [2, 11, 12, 13, 14, 15, 16, 17, 18,
19] and related areas, as, for example, Early Universe cosmology [20, 21, 22,
23, 24, 25, 26, 27, 28], black holes [29, 30, 31, 32, 33, 34, 35], or heavy-
ion physics [36, 37, 38, 39, 40, 41, 42, 43]. In addition, in recent years,
there have been investigations of gravitationally-induced decoherence [44, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54] and decoherence due to time dilation [55,
56], as well as studies of the induction of open quantum dynamics due to
interactions with environments comprising light scalar fields originating from
modified gravity theories [16, 57, 58].
In the theory of open quantum systems we usually deal with a system and its
surrounding environment, which together form the total, closed system. As is
the case in the present article, the open system could, for example, be a real
quantum scalar field $\phi$ which is interacting with an environment in form
of another real scalar field $\chi$. The total system consisting of $\phi$ and
$\chi$ may be described via a density operator $\hat{\rho}_{\phi\chi}(t)$
fulfilling the quantum Liouville equation [1]. Tracing out the environmental
degrees of freedom allows us to write down a quantum master equation for the
elements of the reduced density matrix $\rho_{\phi}(t)$ in a basis of our
choice (see Ref. [16] for a proper, first-principle way of doing this). While
the evolution of the closed system is strictly unitary, $\rho_{\phi}(t)$
usually displays also non-unitary behavior leading to typical open quantum
dynamical effects like phase shifts, or changes of coherence in form of re- or
decoherence. Imposing the Born approximation, which assumes a weak coupling
between system and environment, and/or the Markov approximation, which
basically ignores memory effects in the environment due to its interaction
with the system, might significantly simplify a master equation [1].
Nevertheless, actually solving a master equation often poses an intricate or
even analytically impossible task. Therefore, especially if we are interested
in making actual physical predictions for experiments, it might prove to be
more useful to have a way of directly computing reduced density matrix
elements without even having to write down a quantum master equation. That
such a way actually exists, will be demonstrated by example in the present
article. For this, we use the Lehmann-Szymanzik-Zimmermann- or LSZ-like
reduction [59] that was first presented in Ref. [16], but also discussed in
Refs. [57, 58]. There, the authors used methods from non-equilibrium quantum
field theory like thermo field dynamics (TFD) [60, 61, 62], see also Ref.
[63], the Schwinger-Keldysh formalism [64, 65] and the Feynman-Vernon
influence functional [66] in order to derive a master equation for a reduced
density matrix in a single-particle momentum subspace. Using the well-known
fact that the Schrödinger picture quantum Liouville equation may be expressed
within TFD as a Schrödinger-like equation, which has well-understood solutions
at the operator level, we apply the LSZ-like reduction for perturbatively
finding explicit expressions for elements of the reduced density matrix
$\rho_{\phi}$ in a momentum basis for any number of superposed or non-
superposed particles in terms of the Feynman-Vernon influence functional.
Clearly, for this, we must impose the Born approximation. However, our
approach allows for the description of non-Markovian open quantum systems.
In principle, our method presented here also allows for the treatment of
systems at finite temperatures by replacing the usual propagators by thermal
ones (see e.g. Ref. [67]), as was also done in Ref. [16]. Though, for
simplicity and since we only want to demonstrate the general technique, we
will not consider thermal effects in what follows. In addition, if we were
particularly interested in discrete qudit systems, using the Schrödinger-like
form of the quantum Liouville equation in TFD in combination with the results
from Ref. [68] might allow us to obtain exact, non-perturbative solutions in
some situations.
The present article is structured as follows: in Sec. II we first introduce
the necessary theoretical concepts and then use those in order to derive a
general formula for the superposed single-particle element of the reduced
density matrix $\rho_{\phi}(t)$ in a momentum basis. From this result we can
then easily extrapolate a more general formula, which is valid for any number
of particles. Next, in Sec. III, we consider a simple, concrete example for
the system field $\phi$ and the environmental field $\chi$. More specifically,
we apply the derived formula in order to find explicit expressions for the
associated vacuum, single-particle and superposed single-particle matrix
elements of $\rho_{\phi}(t)$. Finally, in Sec. IV, we draw our conclusions.

##  II Derivation

Here we are going to derive an expression that allows for the direct
computation of reduced density matrix elements in terms of the Feynman-Vernon
influence functional. Certainly, the methods presented in what follows can be
applied to a variety of systems and environments. However, for simplicity and
with applications to field theoretical systems in mind, we just consider a
real scalar field $\phi$ as the system and another real scalar field $\chi$ as
its surrounding environment. In order to better motivate the system-
environment split, we may assume that $\phi$ has a constant mass $M$ much
larger than $\chi$’s constant mass $m$, such that $m/M\ll 1$. This situation
is, for example, comparable to the one described in Ref. [16], where the
system scalar field was used as a proxy for a heavy atom in atom
interferometry that is interacting with an environment comprising fluctuations
of a very light scalar field.
In Secs. II.1-II.3 we will first present the necessary mathematical and
conceptual preliminaries before finally combining them for the actual
derivation in Sec. II.4.

###  II.1 Density matrices in Fock space

Since we are working with a field theoretical system and want to be as general
as possible, we must consider density matrix elements in Fock space, that
allow for any number of particles. In order to closely follow the treatment in
Ref. [16], we choose to also work in a momentum basis and expand density
operators accordingly. Consequently, the most general expression for a density
operator expanded in this basis for any occupation number in Fock space is
given by

$\displaystyle\hat{\rho}(t)$ $\displaystyle=$
$\displaystyle\sum\limits_{i,j=0}^{\infty}\int\left(\prod\limits_{a=1}^{i}d\Pi_{\mathbf{k}^{(a)}}\right)\left(\prod\limits_{b=1}^{j}d\Pi_{\mathbf{l}^{(b)}}\right)\rho_{i,j}(\mathbf{k}^{(1)},...,\mathbf{k}^{(i)};\mathbf{l}^{(1)},...,\mathbf{l}^{(j)};t)\ket{\mathbf{k}^{(1)},...,\mathbf{k}^{(i)}}\bra{\mathbf{l}^{(1)},...,\mathbf{l}^{(j)}}\,\,\,,\,\,\,\,\,\,\,\,\,\,\,$
(1)

where the cases $i=0$ or $j=0$ correspond to the static vacuum state $\ket{0}$
or $\bra{0}$, respectively, each density matrix element is given by

$\displaystyle\rho_{i,j}(\mathbf{k}^{(1)},...,\mathbf{k}^{(i)};\mathbf{l}^{(1)},...,\mathbf{l}^{(j)};t)$
$\displaystyle=$
$\displaystyle\bra{\mathbf{k}^{(1)},...,\mathbf{k}^{(i)}}\hat{\rho}(t)\ket{\mathbf{l}^{(1)},...,\mathbf{l}^{(j)}}\,\,\,,$
(2)

and we use

$\displaystyle\int d\Pi_{\mathbf{k}}$ $\displaystyle:=$
$\displaystyle\int_{\mathbf{k}}\frac{1}{2E_{\mathbf{k}}}$ (3)

with

$\displaystyle\int_{\mathbf{k}}$ $\displaystyle:=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\,\,\,.$ (4)

We note that, as usual, the following has to hold:

$\displaystyle\rho_{i,j}(\mathbf{k}^{(1)},...,\mathbf{k}^{(i)};\mathbf{l}^{(1)},...,\mathbf{l}^{(j)};t)$
$\displaystyle=$
$\displaystyle\rho_{j,i}^{\ast}(\mathbf{l}^{(1)},...,\mathbf{l}^{(j)};\mathbf{k}^{(1)},...,\mathbf{k}^{(i)};t)\,\,\,.$
(5)

Physically these density matrix elements can be interpreted as follows:
$\rho_{0,0}$ describes a $0$-particle or vacuum state, $\rho_{i,0}$ and
$\rho_{0,j}$ stand for $i$\- and $j$-particle states that are not in a
superposition, while $\rho_{i,j}$ represents states with $i$ particles out of
which $j$ are in superposition, if $i\geq j$, or states with $j$ particles out
of which $i$ are in superposition if $i\leq j$.
As will become more apparent at a later point, it will be useful for us to
initially describe the density operator in the Schrödinger picture. In this
case, its time evolution is given by the quantum Liouville equation [1] (we
use $\hbar\equiv 1$ throughout the entire article)

$\displaystyle\frac{\partial}{\partial t}\hat{\rho}_{S}(t)$ $\displaystyle=$
$\displaystyle-\mathrm{i}[\hat{H}_{S}(t),\hat{\rho}_{S}(t)]\,\,\,,$ (6)

which is solved by

$\displaystyle\hat{\rho}_{S}(t)$ $\displaystyle=$ $\displaystyle
e^{-\mathrm{i}\hat{H}_{S}t}\hat{\rho}(0)e^{\mathrm{i}\hat{H}_{S}t}\,\,\,$ (7)

if the Hamiltonian is time-independent, or by

$\displaystyle\hat{\rho}_{S}(t)$ $\displaystyle=$
$\displaystyle(\mathrm{T}e^{-\mathrm{i}\int^{t}_{0}d\tau\hat{H}_{S}(\tau)})\hat{\rho}(0)(\mathrm{T}e^{\mathrm{i}\int^{t}_{0}d\tau\hat{H}_{S}(\tau)})\,\,\,$
(8)

if the Hamiltonian is time-dependent due to an external source, where the
index $S$ labels objects in the Schrödinger picture and $\mathrm{T}$ stands
for time-ordering. Note that Eqn. (8) is more general than the solution in
Eqn. (7) and recovers the latter if $\hat{H}_{S}(t)$ is constant.

###  II.2 The Feynman-Vernon influence functional

When describing an open system in field theory, the Feynman-Vernon influence
functional [66] based on the Schwinger-Keldysh closed-time-path formalism [64,
65] is often the tool of choice, see e.g. Ref. [2]. It 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 Fig. 1 between an initial
$t_{\text{initial}}=0$ and a final time $t_{\text{final}}=t$. Physically the
Feynman-Vernon influence functional captures the evolution of the open system
under influence of the environment, while not describing the environmental
degrees of freedom directly.

Figure 1: Schematic depiction of the closed time path for a density matrix
$\rho$ evolving from an initial time $0$ to a final time $t$ and backwards

Since the Feynman-Vernon influence functional will be an essential object in
the equation that we are later going to derive, following the elaborations in
and using the notations from Ref. [16], we will now sketch how it comes into
play when tracing out the environmental degrees of freedom of the total
density operator. For this, we begin with the usual assumption that system and
environment can be separated at the initial time, such that

$\displaystyle\hat{\rho}_{\phi\chi}(0)$ $\displaystyle=$
$\displaystyle\hat{\rho}_{\phi}(0)\otimes\hat{\rho}_{\chi}(0)\,\,\,.$ (9)

Since we are only interested in the evolution of the reduced density operator
$\hat{\rho}_{\phi}(t)$, we trace out the environmental degrees of freedom:

$\displaystyle\hat{\rho}_{\phi}(t)$ $\displaystyle=$
$\displaystyle\mathrm{Tr}_{\chi}\hat{\rho}_{\phi\chi}(t)\,\,\,.$ (10)

Projecting this expression into a field basis living on the closed-time-path,
we can express the reduced density functional as

$\displaystyle\rho_{\phi}[\phi^{\pm}_{t};t]$ $\displaystyle:=$
$\displaystyle\langle\phi^{+}_{t}|\hat{\rho}_{\phi}(t)|\phi^{-}_{t}\rangle$
(11) $\displaystyle=$ $\displaystyle\int
d\chi^{\pm}_{t}\delta(\chi^{+}_{t}-\chi^{-}_{t})\rho_{\phi\chi}[\phi^{\pm}_{t},\chi^{\pm}_{t};t]\,\,\,,$

where $\pm$ indicates a dependence on both $+$\- and $-$-type operators,
subscript $t$ labels the time slice on which the field eigenstates have to be
taken, and

$\displaystyle\rho_{\phi\chi}[\phi^{\pm}_{t},\chi^{\pm}_{t};t]$
$\displaystyle:=$
$\displaystyle\langle\phi^{+}_{t},\chi^{+}_{t}|\hat{\rho}_{\phi\chi}(t)|\phi^{-}_{t},\chi^{-}_{t}\rangle\,\,\,.$
(12)

Furthermore, taking into account Eq. (9), the reduced density functional at
time $t$ can be expressed as

$\displaystyle\rho_{\phi}[\phi^{\pm}_{t};t]$ $\displaystyle=$
$\displaystyle\int
d\phi_{0}^{\pm}\mathcal{I}[\phi^{\pm}_{t},\phi^{\pm}_{0};t,0]\rho_{\phi}[\phi^{\pm}_{0};0]\,\,\,,$
(13)

where the time translation is given in terms of the so-called influence
functional (IF) propagator, see e.g. Ref.[2],

$\displaystyle\mathcal{I}[\phi^{\pm}_{t},\phi^{\pm}_{0};t,0]$ $\displaystyle=$
$\displaystyle\int^{\phi^{\pm}_{t}}_{\phi^{\pm}_{0}}\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\text{eff}}[\phi;t]}\,\,\,.$
(14)

The latter contains two path integrals over an effective action

$\displaystyle\widehat{S}_{\text{eff}}[\phi;t]$ $\displaystyle=$
$\displaystyle\widehat{S}_{\phi}[\phi;t]+\widehat{S}_{\phi,\text{int}}[\phi;t]+\widehat{S}_{\text{IF}}[\phi;t]\,\,\,,$
(15)

where $\widehat{}$ indicates functionals that depend on both field variables
$\phi^{+}$ and $\phi^{-}$. For example, the free action and the self-
interaction of $\phi$ are given by

$\displaystyle\widehat{S}_{\phi}[\phi;t]$ $\displaystyle:=$ $\displaystyle
S_{\phi}[\phi^{+};t]-S_{\phi}[\phi^{-};t]\,\,\,,$
$\displaystyle\widehat{S}_{\phi,\text{int}}[\phi;t]$ $\displaystyle:=$
$\displaystyle
S_{\phi,\text{int}}[\phi^{+};t]-S_{\phi,\text{int}}[\phi^{-};t]\,\,\,.$ (16)

However, the influence action $\widehat{S}_{\text{IF}}[\phi;t]$ has a more
involved definition by the Feynman-Vernon influence functional

$\displaystyle\widehat{\mathcal{F}}[\phi;t]$ $\displaystyle=$
$\displaystyle\exp\left\\{\mathrm{i}\widehat{S}_{\text{IF}}[\phi;t]\right\\}$
(17) $\displaystyle=$ $\displaystyle\int
d\chi^{\pm}_{t}d\chi^{\pm}_{0}\delta(\chi^{+}_{t}-\chi^{-}_{t})\rho_{\chi}[\chi^{\pm}_{0};0]$
$\displaystyle\times\int^{\chi^{\pm}_{t}}_{\chi^{\pm}_{0}}\mathcal{D}\chi^{\pm}\exp\left\\{\mathrm{i}\Big{(}\widehat{S}_{\chi}[\chi;t]+\widehat{S}_{\chi,\text{int}}[\chi;t]+\widehat{S}_{\text{int}}[\phi,\chi;t]\Big{)}\right\\}\,\,\,,$

where $\widehat{S}_{\chi}[\chi;t]$ is the free action,
$\widehat{S}_{\chi,\text{int}}[\chi;t]$ the self-interaction action for $\chi$
and $\widehat{S}_{\text{int}}[\phi,\chi;t]$ an action describing the
interaction between system and environment. Since we are working within a
finite time interval, we define all actions used here only on
$\Omega_{t}:=[0,t]\times\mathbb{R}^{3}$, which means they all have the
structure

$\displaystyle S[t]$ $\displaystyle=$
$\displaystyle\int_{x\in\Omega_{t}}\mathcal{L}[x]$ (18)

with

$\displaystyle\int_{x}$ $\displaystyle:=$ $\displaystyle\int d^{4}x\,\,\,.$
(19)

In addition, note that we are working with two types of functional integrals:
one, denoted $d$, considers fields over all of $\mathbb{R}^{3}$, but only at
one particular time slice, while the other, denoted $\mathcal{D}$, is over the
whole $\Omega_{t}$. See Eq. (17) for an example in which both types appear
together.
Finally, we define an expectation value with respect to $\chi$ as (cf. Refs.
[2, 16])

$\displaystyle\langle A[\chi^{a}]\rangle_{\chi}$ $\displaystyle:=$
$\displaystyle\int
d\chi^{\pm}_{t}d\chi^{\pm}_{0}\delta(\chi_{t}^{+}-\chi_{t}^{-})\rho_{\chi}[\chi^{\pm}_{0};0]\int^{\chi^{\pm}_{t}}_{\chi^{\pm}_{0}}\mathcal{D}\chi^{\pm}A[\chi^{a}]\exp\left\\{\mathrm{i}\widehat{S}_{\chi}[\chi;t]\right\\}\,\,\,,$
(20)

and find for the Feynman-Vernon influence functional:

$\displaystyle\widehat{\mathcal{F}}[\phi;t]$ $\displaystyle=$
$\displaystyle\left\langle\exp\big{\\{}\mathrm{i}\big{(}\widehat{S}_{\chi,\text{int}}[\chi;t]+\widehat{S}_{\text{int}}[\phi,\chi;t]\big{)}\big{\\}}\right\rangle_{\chi}\,\,\,.$
(21)

###  II.3 Thermo field dynamics

Next, we look at thermo field dynamics (TFD) [60, 61, 62] (see also Ref. [63])
as the final ingredient we need. Again we are following the description laid
out in Ref. [16]. To some extent, we can understand TFD as an algebraic
formulation of the Schwinger-Keldysh formalism, which itself works with
operators on a positive $(+)$ or negative $(-)$ complex time axis, in a
doubled Hilbert space

$\displaystyle\widehat{\mathcal{H}}$ $\displaystyle:=$
$\displaystyle\mathcal{H}^{+}\otimes\mathcal{H}^{-}\,\,\,,$ (22)

where $\mathcal{H}^{\pm}$ are the Hilbert spaces corresponding to the
$\pm$-branches of the closed time path, respectively. Operators living on the
closed time path can be expressed within TFD as

$\displaystyle\hat{\mathcal{O}}^{+}$ $\displaystyle=$
$\displaystyle\hat{\mathcal{O}}\otimes\hat{\mathbb{I}}\,\,\,,$
$\displaystyle\hat{\mathcal{O}}^{-}$ $\displaystyle=$
$\displaystyle\hat{\mathbb{I}}\otimes\hat{\mathcal{O}}^{\mathcal{T}}\,\,\,,$
(23)

where $\hat{\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

$\displaystyle|0\rangle\\!\rangle$ $\displaystyle:=$
$\displaystyle\ket{0}\otimes\ket{0}$ (24)

via creation operators

$\displaystyle\hat{a}^{+\dagger}_{\mathbf{k}}|0\rangle\\!\rangle\,=\,\ket{\mathbf{k}}\otimes\ket{0}\,=:\,|\mathbf{k}_{+}\rangle\\!\rangle\,\,\,,\,\,\,\,\,\,\,\,\,\hat{a}^{-\dagger}_{\mathbf{k}}|0\rangle\\!\rangle\,=\,\ket{0}\otimes\ket{\mathbf{k}}\,=:\,|\mathbf{k}_{-}\rangle\\!\rangle\,\,\,.$
(25)

Consequently, the corresponding annihilators act like

$\displaystyle\hat{a}^{\pm}_{\mathbf{k}}|\mathbf{p}_{+},\mathbf{p}_{-}\rangle\\!\rangle\,=\,(2\pi)^{3}2E^{\phi}_{\mathbf{k}}\delta^{(3)}(\mathbf{p}-\mathbf{k})|\mathbf{p}_{\mp}\rangle\\!\rangle\,\,\,.$
(26)

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

$\displaystyle|1\rangle\\!\rangle$ $\displaystyle:=$
$\displaystyle|0\rangle\\!\rangle+\int
d\Pi_{\mathbf{k}}|\mathbf{k}_{+},\mathbf{k}_{-}\rangle\\!\rangle+\frac{1}{2!}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k^{\prime}}}|\mathbf{k}_{+},\mathbf{k^{\prime}}_{+},\mathbf{k}_{-},\mathbf{k^{\prime}}_{-}\rangle\\!\rangle+...\,\,\,,$
(27)

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

$\displaystyle\braket{\hat{\mathcal{O}}(t)}$ $\displaystyle=$
$\displaystyle\text{Tr}\hat{\mathcal{O}}(t)\hat{\rho}(t)\,=\,\langle\\!\langle
1|\hat{\mathcal{O}}^{+}(t)\hat{\rho}^{+}(t)|1\rangle\\!\rangle\,\,\,.$ (28)

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

$\displaystyle\frac{\partial}{\partial
t}\hat{\rho}_{S}^{+}(t)|1\rangle\\!\rangle_{S}$ $\displaystyle=$
$\displaystyle-\mathrm{i}\widehat{H}_{S}(t)\hat{\rho}_{S}^{+}(t)|1\rangle\\!\rangle_{S}\,\,\,,$
(29)

where
$\widehat{H}_{S}(t):=\hat{H}_{S}(t)\otimes\mathbb{I}-\mathbb{I}\otimes\hat{H}_{S}(t)$.
This can generally be solved by

$\displaystyle\hat{\rho}_{S}^{+}(t)|1\rangle\\!\rangle_{S}$ $\displaystyle=$
$\displaystyle\text{T}\exp\left\\{{-\mathrm{i}\int\limits_{0}^{t}\widehat{H}_{S}(\tau)d\tau}\right\\}\hat{\rho}^{+}(0)|1\rangle\\!\rangle_{S}\,\,\,.$
(30)

Note that at time $0$ the different pictures coincide and we therefore dropped
the label $S$.
Analogously, we write down the Schrödinger-like form of the quantum Liouville
equation for the reduced density operator that can be obtained by tracing out
the $\chi$ degrees of freedom from the total density operator
$\hat{\rho}_{\phi\chi}$:

$\displaystyle\frac{\partial}{\partial
t}\hat{\rho}_{\phi,S}^{+}(t)|1\rangle\\!\rangle_{S}$ $\displaystyle=$
$\displaystyle-\mathrm{i}\widehat{H}_{\phi,S}(t)\hat{\rho}_{\phi,S}^{+}(t)|1\rangle\\!\rangle_{S}\,\,\,,$
(31)

where the effective Hamiltonian is given by the usual unitary evolution due to
the Hamiltonian $\hat{H}_{\phi,S}$ and a non-unitarian
$\widehat{H}_{\text{IF},S}$:

$\displaystyle\widehat{H}_{\phi,S}(t)$ $\displaystyle=$
$\displaystyle\hat{H}_{\phi,S}(t)\otimes\mathbb{I}-\mathbb{I}\otimes\hat{H}_{\phi,S}(t)+\widehat{H}_{\text{IF},S}(t)\,\,\,.$
(32)

We recall that, generally, Hamiltonian and action are related via
$H(t)=-\frac{\partial}{\partial t}S(t)$. Furthermore, if we consider that the
effective action in Eq. (15) describes the full evolution of the reduced
density matrix elements in the $\phi$-basis, the interaction picture
$\widehat{H}_{\phi,I}$ must be its corresponding Hamiltonian, such that
$\hat{H}_{\phi,I}(t):=\hat{H}_{0,I}(t)+\hat{H}_{\text{int},I}(t)$ represents
$S_{\phi}(t)+S_{\phi,\text{int}}(t)$ and $\widehat{H}_{\text{IF},I}(t)$ the
influence action given in Eq. (17). This will later be of greater importance
when we write down a path integral expression for the reduced density matrix
elements. For now, however, we just remind ourselves that the free Hamiltonian
$\hat{H}_{0}$ is the same in Schrödinger and interaction picture, and
consequently drop the subscript.
Since this will later be useful for us, we will now transform Eq. (31) into
the interaction picture. For this, we remind ourselves that an operator
transforms like

$\displaystyle\hat{\mathcal{O}}_{S}$ $\displaystyle=$ $\displaystyle
e^{-\mathrm{i}\hat{H}_{0}t}\hat{\mathcal{O}}_{I}(t)e^{\mathrm{i}\hat{H}_{0}t}$
(33)

from the Schrödinger into the interaction picture, and conveniently define a
unitary evolution operator

$\displaystyle\widehat{U}$ $\displaystyle:=$ $\displaystyle
e^{-\mathrm{i}\hat{H}_{0}t}\otimes e^{\mathrm{i}\hat{H}_{0}t}\,\,\,,$ (34)

such that operators acting on the TFD doubled Hilbert space transform like

$\displaystyle\widehat{\mathcal{O}}_{S}$ $\displaystyle:=$
$\displaystyle\widehat{U}\widehat{\mathcal{O}}_{I}(t)\widehat{U}^{\dagger}\,\,\,.$
(35)

The state introduced in Eq. (27) transforms like a base vector:

$\displaystyle\widehat{U}^{\dagger}|1\rangle\\!\rangle_{S}$ $\displaystyle=$
$\displaystyle|1(t)\rangle\\!\rangle_{I}\,\,\,.$ (36)

However, it can straightforwardly be seen that $\widehat{U}$ and
$\widehat{U}^{\dagger}$ act like unit operators on $|1\rangle\\!\rangle_{S}$
(and consequently on $|1(t)\rangle\\!\rangle$ as well):

$\displaystyle\widehat{U}^{\dagger}|1\rangle\\!\rangle_{S}$ $\displaystyle=$
$\displaystyle\ket{0}\otimes\ket{0}+\int
d\Pi_{\mathbf{k}}e^{iE_{\mathbf{k}}t}\ket{\mathbf{k}}\otimes
e^{-iE_{\mathbf{k}}t}\ket{\mathbf{k}}+\frac{1}{2!}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k^{\prime}}}e^{i(E_{\mathbf{k}}+E_{\mathbf{k^{\prime}}})t}\ket{\mathbf{k}}\otimes
e^{-i(E_{\mathbf{k}}+E_{\mathbf{k^{\prime}}})t}\ket{\mathbf{k}}+...$ (37)
$\displaystyle=$ $\displaystyle|1\rangle\\!\rangle_{S}\,\,\,,$

which implies that the state is time- and picture-independent, and will
therefore simply be denoted by $|1\rangle\\!\rangle$ in what follows.
Using all this, Eq. (31) becomes in the interaction picture:

$\displaystyle\frac{\partial}{\partial
t}\left(\widehat{U}\hat{\rho}_{\phi,I}^{+}(t)\widehat{U}^{\dagger}\right)|1\rangle\\!\rangle$
$\displaystyle=$
$\displaystyle-\mathrm{i}\widehat{U}\widehat{H}_{\phi,I}(t)\hat{\rho}_{\phi,I}^{+}(t)|1\rangle\\!\rangle\,\,\,.$
(38)

Evaluating the partial derivative on the left-hand side of Eq. (38) leaves us
with

$\displaystyle\left[(\partial_{t}\widehat{U})\hat{\rho}_{\phi,I}^{+}(t)+\widehat{U}\partial_{t}\hat{\rho}_{\phi,I}^{+}(t)+\widehat{U}\hat{\rho}_{\phi,I}^{+}(t)\partial_{t}\widehat{U}^{\dagger}\right]|1\rangle\\!\rangle$
$\displaystyle=$
$\displaystyle\widehat{U}\left[-\mathrm{i}\widehat{H}_{0}\hat{\rho}_{\phi,I}^{+}(t)+\partial_{t}\hat{\rho}_{\phi,I}^{+}(t)+\mathrm{i}\hat{\rho}_{\phi,I}^{+}(t)\widehat{H}_{0}\right]|1\rangle\\!\rangle\,\,\,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,$
(39)

where
$\widehat{H}_{0}:=\hat{H}_{0}\otimes\mathbb{I}-\mathbb{I}\otimes\hat{H}_{0}$.
Acting with $\widehat{H}_{0}$ on $|1\rangle\\!\rangle$ gives nil, as can
straightforwardly be seen. Consequently, we find for Eq. (38):

$\displaystyle\left[-\mathrm{i}\widehat{H}_{0}\hat{\rho}_{\phi,I}^{+}(t)+\partial_{t}\hat{\rho}_{\phi,I}^{+}(t)\right]|1\rangle\\!\rangle$
$\displaystyle=$
$\displaystyle-\mathrm{i}\widehat{H}_{\phi,I}(t)\hat{\rho}_{\phi,I}^{+}(t)|1\rangle\\!\rangle\,\,\,.$
(40)

Since $\widehat{H}_{\phi,I}(t)$ also contains $\widehat{H}_{0}$, we finally
arrive at

$\displaystyle\partial_{t}\hat{\rho}_{\phi,I}^{+}(t)|1\rangle\\!\rangle$
$\displaystyle=$
$\displaystyle-\mathrm{i}[\widehat{H}_{\text{int},I}(\tau)+\widehat{H}_{\text{IF},I}(\tau)]\hat{\rho}_{\phi,I}^{+}(t)|1\rangle\\!\rangle\,\,\,,$
(41)

where
$\widehat{H}_{\text{int},I}(t):=\hat{H}_{\text{int},I}(t)\otimes\mathbb{I}-\mathbb{I}\otimes\hat{H}_{\text{int},I}(t)$,
and which can be solved by

$\displaystyle\hat{\rho}_{\phi,I}^{+}(t)|1\rangle\\!\rangle$ $\displaystyle=$
$\displaystyle\text{T}\exp\left\\{{-\mathrm{i}\int\limits_{0}^{t}[\widehat{H}_{\text{int},I}(\tau)+\widehat{H}_{\text{IF},I}(\tau)]d\tau}\right\\}\hat{\rho}_{\phi}^{+}(0)|1\rangle\\!\rangle\,\,\,.$
(42)

###  II.4 Reduced density matrix elements

Finally, we have all we need for deriving an equation that allows us to
directly compute elements of the reduced density matrix $\rho_{\phi}(t)$ in a
momentum basis for any occupation number in Fock space. For notational
convenience, we will from now on drop the index $\phi$ when talking about the
reduced density matrix and its elements. In addition, we will from now on only
work in the interaction picture and therefore not use any corresponding labels
for operators and states anymore.
At first, we will derive an expression for the density matrix element
$\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$, which represents a single
particle in superposition in momentum space. From the resulting equation and
its derivation it will be straightforward for us to infer for expressions for
the other density matrix elements.
As in Ref. [16], we start with

$\displaystyle\langle\mathbf{p};t|\hat{\rho}(t)|\mathbf{p}^{\prime};t\rangle$
$\displaystyle=$
$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)\,\,\,.$ (43)

Note that, as was pointed out in Refs. [69, 70], the matrix element
$\rho_{1,1}(\mathbf{p},\mathbf{p}\,^{\prime};t)$ is picture-independent. In
TFD, Eq. (43) can be expressed as

$\displaystyle\text{Tr}|\mathbf{p}^{\prime};t\rangle\langle\mathbf{p};t|\hat{\rho}(t)$
$\displaystyle=$ $\displaystyle\langle\langle
1|(|\mathbf{p}^{\prime};t\rangle\langle\mathbf{p};t|\otimes\hat{\mathbb{I}})\hat{\rho}^{+}(t)|1\rangle\rangle\,\,\,.$
(44)

Substituting Eq. (42), we are left with

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$ $\displaystyle=$
$\displaystyle\langle\langle\mathbf{p}_{+},\mathbf{p}^{\prime}_{-};t|\text{T}\exp\left\\{{-\mathrm{i}\int\limits_{0}^{t}[\widehat{H}_{\text{int}}(\tau)+\widehat{H}_{\text{IF}}(\tau)]d\tau}\right\\}\hat{\rho}^{+}(0)|1\rangle\rangle\,\,\,.$
(45)

Next, we expand the density operator in Eq. (45) as in Eq. (1). However, we
assume that the exact number of system particles and their superposedness at
the initial time are practically well-understood, such that we only have to
consider the $\rho_{i,j}(0)$ elements for one particular choice of $i$ and
$j$. In our case, we choose $i=j=1$ and assume all other initial elements to
be nil. Though, if there was a reason to assume that the number of system
particles is not certain or it is not clear whether a superposition was
created, then the derivation could be easily modified. So, for the sake of
readability, we do not consider the most general case here, but include it
later when we extrapolate the most general expression. In our current case we
find

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$ $\displaystyle=$
$\displaystyle\langle\langle\mathbf{p}_{+},\mathbf{p}^{\prime}_{-};t|\text{T}\exp\left\\{{-\mathrm{i}\int\limits_{0}^{t}[\widehat{H}_{\text{int}}(\tau)+\widehat{H}_{\text{IF}}(\tau)]d\tau}\right\\}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)|\mathbf{k}_{+},\mathbf{k}^{\prime}_{-}\rangle\rangle\,\,\,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,$
(46)

which, using Eqs. (25) and (26) becomes

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$ $\displaystyle=$
$\displaystyle\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
(47) $\displaystyle\times\langle\langle
0|\text{T}\hat{a}^{+}_{\mathbf{p}}(t)\hat{a}^{-}_{\mathbf{p}^{\prime}}(t)\exp\left\\{{-\mathrm{i}\int\limits_{0}^{t}[\widehat{H}_{\text{int}}(\tau)+\widehat{H}_{\text{IF}}(\tau)]d\tau}\right\\}\hat{a}^{+\dagger}_{\mathbf{k}}(0)\hat{a}^{-\dagger}_{\mathbf{k}^{\prime}}(0)|0\rangle\rangle\,\,\,.\,\,\,\,\,\,\,$

After replacing the creation and annihilation operators by (cf. Ref. [16])

$\displaystyle\hat{a}^{+}_{\mathbf{p}}(t)$ $\displaystyle=$
$\displaystyle+\mathrm{i}\int_{\mathbf{x}}e^{-\mathrm{i}\mathbf{p}\cdot\mathbf{x}}\partial_{t,E^{\phi}_{\mathbf{p}}}\hat{\phi}^{+}(t,\mathbf{x})\,\,\,,\,\,\,\,\,\,\,\,\,\hat{a}^{+\dagger}_{\mathbf{p}}(t)\,=\,-\mathrm{i}\int_{\mathbf{x}}e^{+\mathrm{i}\mathbf{p}\cdot\mathbf{x}}\partial_{t,E^{\phi}_{\mathbf{p}}}^{*}\hat{\phi}^{+}(t,\mathbf{x})\,\,\,,$
$\displaystyle\hat{a}^{-}_{\mathbf{p}}(t)$ $\displaystyle=$
$\displaystyle-\mathrm{i}\int_{\mathbf{x}}e^{+\mathrm{i}\mathbf{p}\cdot\mathbf{x}}\partial_{t,E^{\phi}_{\mathbf{p}}}^{*}\hat{\phi}^{-}(t,\mathbf{x})\,\,\,,\,\,\,\,\,\,\,\,\,\hat{a}^{-{\dagger}}_{\mathbf{p}}(t)\,=\,+\mathrm{i}\int_{\mathbf{x}}e^{-\mathrm{i}\mathbf{p}\cdot\mathbf{x}}\partial_{t,E^{\phi}_{\mathbf{p}}}\hat{\phi}^{-}(t,\mathbf{x})\,\,\,,$
(48)

where
$\partial_{t,E^{\phi}_{\mathbf{p}}}:=\overset{\rightarrow}{\partial}_{t}-\mathrm{i}E^{\phi}_{\mathbf{p}}$,
we obtain

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$ $\displaystyle=$
$\displaystyle\lim_{\begin{subarray}{c}x^{0(\prime)}\,\to\,t^{+}\\\
y^{0(\prime)}\,\to\,0^{-}\end{subarray}}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
(49)
$\displaystyle\times\int_{\mathbf{x}\mathbf{x}^{\prime}\mathbf{y}\mathbf{y}^{\prime}}e^{-\mathrm{i}(\mathbf{p}\cdot\mathbf{x}-\mathbf{p}^{\prime}\cdot\mathbf{x}^{\prime})+\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{x^{0\prime},E^{\phi}_{\mathbf{p}^{\prime}}}^{*}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}$
$\displaystyle\times\langle\langle 0|{\rm
T}[\hat{\phi}^{+}(x)\hat{\phi}^{-}(x^{\prime})\exp\left\\{{-\mathrm{i}\int\limits_{0}^{t}[\widehat{H}_{\text{int}}(\tau)+\widehat{H}_{\text{IF}}(\tau)]d\tau}\right\\}\hat{\phi}^{+}(y)\hat{\phi}^{-}(y^{\prime})]|0\rangle\rangle\,\,\,.\,\,\,\,\,$

Here we introduced limits in which $x^{0(\prime)}$ approach $t$ from above and
$y^{0(\prime)}$ approach $0$ from below in order to recover the correct time
ordering.
Translating the resulting expression in Eq. (49) into the path integral
formalism, gives us

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$ $\displaystyle=$
$\displaystyle\lim_{\begin{subarray}{c}x^{0(\prime)}\,\to\,t^{+}\\\
y^{0(\prime)}\,\to\,0^{-}\end{subarray}}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
(50)
$\displaystyle\times\int_{\mathbf{x}\mathbf{x}^{\prime}\mathbf{y}\mathbf{y}^{\prime}}e^{-\mathrm{i}(\mathbf{p}\cdot\mathbf{x}-\mathbf{p}^{\prime}\cdot\mathbf{x}^{\prime})+\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{x^{0\prime},E^{\phi}_{\mathbf{p}^{\prime}}}^{*}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}$
$\displaystyle\times\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\phi^{+}(x)\phi^{-}(x^{\prime})$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\exp\left\\{{\mathrm{i}\int\limits_{0}^{t}\frac{d}{d\tau}[\widehat{S}_{\phi,\text{int}}[\phi;\tau]+\widehat{S}_{\text{IF}}[\phi;\tau]]d\tau}\right\\}\phi^{+}(y)\phi^{-}(y^{\prime})\,\,\,.$

For obtaining the total time derivatives of the actions, we assumed that the
actions are not depending on the partial time derivatives of $\phi^{\pm}$. In
this case, the corresponding conjugate momenta $\pi^{\pm}$ vanish and
consequently leave

$\displaystyle\frac{d}{d\tau}\widehat{S}_{\text{IF}}(\tau)$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\tau}\widehat{S}_{\text{IF}}(\tau)\,\,=\,\,-\widehat{H}_{\text{IF}}(\tau)\,\,\,,$
(51)

and the same for $\widehat{S}_{\phi,\text{int}}$. Finally, we perform the time
integration in the last line of Eq. (50), substitute Eq. (17), and use that,
by construction, all actions vanish at the initial time, which also means that
the initial influence functional equals $1$. This leads us to

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$ $\displaystyle=$
$\displaystyle\lim_{\begin{subarray}{c}x^{0(\prime)}\,\to\,t^{+}\\\
y^{0(\prime)}\,\to\,0^{-}\end{subarray}}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
(52)
$\displaystyle\times\int_{\mathbf{x}\mathbf{x}^{\prime}\mathbf{y}\mathbf{y}^{\prime}}e^{-\mathrm{i}(\mathbf{p}\cdot\mathbf{x}-\mathbf{p}^{\prime}\cdot\mathbf{x}^{\prime})+\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{x^{0\prime},E^{\phi}_{\mathbf{p}^{\prime}}}^{*}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}$
$\displaystyle\times\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\phi^{+}(x)\phi^{-}(x^{\prime})\exp\left\\{\mathrm{i}\widehat{S}_{\phi,\text{int}}[\phi;t]\right\\}\widehat{\mathcal{F}}[\phi;t]\phi^{+}(y)\phi^{-}(y^{\prime})\,\,\,.$

What we just found is an equation that enables us to directly compute the
momentum basis superposed single-particle elements of the reduced density
matrix $\rho_{\phi}(t)$ without having to explicitly solve any master
equation. Using the Born approximation would now allow us to perturbatively
evaluate Eq. (52) for a particular system-environment model in the same way as
was done in Ref. [16] for a quantum master equation. Doing so will be subject
of Sec. III.
Though, before concluding the present discussion, we will extrapolate the
formula for a general density matrix element with an arbitrary choice of non-
vanishing initial density matrix elements. That this result must be correct
can easily be seen by retracing the steps we took in order to derive Eq. (52).
We find:

$\displaystyle\rho_{g,h}(\mathbf{k}^{(1)},...,\mathbf{k}^{(g)};\mathbf{l}^{(1)},...,\mathbf{l}^{(h)};t)=$
$\displaystyle\sum\limits_{i,j=0}^{\infty}\lim_{\begin{subarray}{c}x_{(1)}^{0},...,x_{(g)}^{0},x_{(1)}^{0\prime},...,x_{(h)}^{0\prime}\,\to\,t^{+}\\\
y_{(1)}^{0},...,y_{(i)}^{0},y_{(1)}^{0\prime},...,y_{(j)}^{0\prime}\,\to\,0^{-}\end{subarray}}\int\left(\prod\limits_{a=1}^{i}d\Pi_{\mathbf{r}^{(a)}}\right)\left(\prod\limits_{b=1}^{j}d\Pi_{\mathbf{s}^{(b)}}\right)\rho_{i,j}(\mathbf{r}^{(1)},...,\mathbf{r}^{(i)};\mathbf{s}^{(1)},...,\mathbf{s}^{(j)};0)$
$\displaystyle\times\int_{\mathbf{x}_{(1)}...\mathbf{x}_{(g)}\mathbf{x}_{(1)}^{\prime}...\mathbf{x}_{(h)}^{\prime}\mathbf{y}_{(1)}...\mathbf{y}_{(i)}\mathbf{y}_{(1)}^{\prime}...\mathbf{y}_{(j)}^{\prime}}\exp\left\\{-\mathrm{i}\Bigg{(}\sum\limits^{g}_{a=1}\mathbf{k}^{(a)}\mathbf{x}_{(a)}-\sum\limits^{h}_{a=1}\mathbf{l}^{(a)}\mathbf{x}_{(a)}^{\prime}\Bigg{)}+\mathrm{i}\Bigg{(}\sum\limits^{i}_{a=1}\mathbf{r}^{(a)}\mathbf{y}_{(a)}-\sum\limits^{j}_{a=1}\mathbf{s}^{(a)}\mathbf{y}_{(a)}^{\prime}\Bigg{)}\right\\}$
$\displaystyle\times\left(\prod\limits_{a=1}^{g}\partial_{x^{0}_{(a)},E^{\phi}_{\mathbf{k}^{(a)}}}\right)\left(\prod\limits_{b=1}^{h}\partial_{x^{0\prime}_{(b)},E^{\phi}_{\mathbf{l}^{(b)}}}^{*}\right)\left(\prod\limits_{c=1}^{i}\partial_{y^{0}_{(c)},E^{\phi}_{\mathbf{r}^{(c)}}}^{*}\right)\left(\prod\limits_{d=1}^{j}\partial_{y^{0\prime}_{(d)},E^{\phi}_{\mathbf{s}^{(d)}}}\right)$
$\displaystyle\times\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\phi^{+}_{x_{(1)}}...\phi^{+}_{x_{(g)}}\phi^{-}_{x^{\prime}_{(1)}}...\phi^{-}_{x^{\prime}_{(h)}}e^{\mathrm{i}\widehat{S}_{\phi,\text{int}}[\phi;t]}\widehat{\mathcal{F}}[\phi;t]\phi^{+}_{y_{(1)}}...\phi^{+}_{y_{(i)}}\phi^{-}_{y^{\prime}_{(1)}}...\phi^{-}_{y^{\prime}_{(j)}}\,\,\,.\,\,\,$
(53)

Clearly, by taking $g=h=1$ and letting all intitial elements apart from the
ones for $i=j=1$ vanish, we recover Eq. (52).
When evaluating Eqs. (52) or (II.4) it is important to understand that, at
least for the system field $\phi$, only contractions of two $+$\- or two
$-$-fields are permitted, while for computing the influence functional via Eq.
(21) all types of contractions of $\chi$ are allowed. That the former is true
for single-particle elements as in Eq. (52) was already pointed out in Ref.
[16]. There the authors showed that the $2\times 2$ matrix propagator used in
their calculation is diagonal in the single-particle subspace at zero
temperature. This is also the case in the calculation presented here. However,
it can be seen that the matrix propagator is not only of diagonal form in the
single-particle but also for all other $n$-particle subspaces. For this, we
consider (schematically) the $g+h+i+j$-point function that was used in the
derivation of Eq. (II.4) as the generalisation of the four-point function in
Eq. (49)

$\displaystyle\langle\langle 0|{\rm
T}[\hat{\phi}^{+}_{x_{(1)}}...\hat{\phi}^{+}_{x_{(g)}}\hat{\phi}^{-}_{x^{\prime}_{(1)}}...\hat{\phi}^{-}_{x^{\prime}_{(h)}}\exp\left\\{...\right\\}\hat{\phi}^{+}_{y_{(1)}}...\hat{\phi}^{+}_{y_{(i)}}\hat{\phi}^{-}_{y^{\prime}_{(1)}}...\hat{\phi}^{-}_{y^{\prime}_{(j)}}]|0\rangle\rangle\,\,\,,$
(54)

which, after using Wick’s theorem [71], reduces to a product of two-point
functions. Taking into account the definition of the TFD vacuum state
$|0\rangle\\!\rangle$, it can be seen that the two-point functions that give
the off-diagonal elements of the matrix propagator must vanish:

$\displaystyle\langle\\!\langle
0|\hat{\phi}^{+(-)}_{x}\hat{\phi}^{-(+)}_{y}|0\rangle\\!\rangle$
$\displaystyle=$ $\displaystyle 0\,\,\,.$ (55)

The diagonal elements, however, give the usual Feynman and Dyson propagators

$\displaystyle\langle\\!\langle
0|\mathrm{T}[\hat{\phi}^{+}_{x}\hat{\phi}^{+}_{y}]|0\rangle\\!\rangle$
$\displaystyle=$ $\displaystyle
D^{++}_{xy}\,=\,D^{\mathrm{F}}_{xy}\,=\,-\mathrm{i}\int_{k}\frac{e^{\mathrm{i}k\cdot(x-y)}}{k^{2}+M^{2}-i\epsilon}\,\,\,,$
(56) $\displaystyle\langle\\!\langle
0|\mathrm{T}[\hat{\phi}^{-}_{x}\hat{\phi}^{-}_{y}]|0\rangle\\!\rangle$
$\displaystyle=$ $\displaystyle
D^{--}_{xy}\,=\,D^{\mathrm{D}}_{xy}\,=\,+\mathrm{i}\int_{k}\frac{e^{\mathrm{i}k\cdot(x-y)}}{k^{2}+M^{2}+i\epsilon}\,\,\,.$
(57)

##  III A simple toy model

Here we are going to apply the results from Sec. II, i.e. Eqs. (52) and
(II.4), to a toy model for system $\phi$ and environment $\chi$. For this, we
choose a model that is as simple as possible, but still allows for a distinct
system-environment split due to an enormous difference in both field’s
constant masses, and gives rise to at least some interesting open quantum
dynamics.
In Sec. III.1 we will introduce the model and compute the resulting Feynman-
Vernon influence functional perturbatively up to second order in the coupling
constant. Next, in Sec. III.2, we will consider the evolution of the matrices
$\rho_{0,0}$, $\rho_{1,0}$ and $\rho_{1,1}$ each for the case that their own
respective initial density matrix is the only one that is non-vanishing. In
addition, we will also investigate how the matrices $\rho_{0,0}$ and
$\rho_{1,0}$ evolve under the assumption that only the initial $\rho_{1,1}$ is
non-vanishing, which corresponds to the system being prepared in a superposed
single-particle state, but then decaying into the vacuum or loosing its
superposedness due to decoherence.
After explicitly evaluating Eq. (II.4) for each case, we will repeatedly
encounter terms with time-dependent divergences. Dealing with such peculiar
terms was briefly discussed in Ref. [16]. However, since we only intend to
demonstrate how to use the formulas we derived in this article and not make
actual experimental predictions for the model we use as an example, we will
leave the divergent terms in the final results. An actual renormalization of
such time-dependent divergences deserves a much more extensive discussion,
which is beyond the scope of the present article, but will be the topic for a
future one [72].

###  III.1 The model

As an example we choose a two-scalar field model with no self-interactions,
and the following actions:

$\displaystyle S_{\phi}[\phi]$ $\displaystyle=$
$\displaystyle\int_{x}\left[-\frac{1}{2}(\partial\phi)^{2}-\frac{1}{2}M^{2}\phi^{2}\right]\,\,\,,$
(58) $\displaystyle S_{\chi}[\chi]$ $\displaystyle=$
$\displaystyle\int_{x}\left[-\frac{1}{2}(\partial\chi)^{2}-\frac{1}{2}m^{2}\chi^{2}\right]\,\,\,,$
(59) $\displaystyle S_{\text{int}}[\phi,\chi]$ $\displaystyle=$
$\displaystyle\int_{x\in\Omega_{t}}\left[-\frac{\alpha}{2}m\chi^{2}\phi\right]\,\,\,,$
(60)

where $m/M,\alpha\ll 1$ and $\Omega_{t}:=[0,t]\times\mathbb{R}^{3}$. Even
though it is rather simple, this model already includes some interesting
features that we would like to illustrate. Taking into account self-
interactions and more coupling terms is therefore not much more illuminating,
but more computationally challenging. A master equation for a density matrix
element derived with the same LSZ-like reduction method used here but for a
model with a self-interaction for $\chi$ and different coupling terms can be
found in Ref. [16].
Using Eq. (21), we expand the Feynman-Vernon influence functional up to second
order in $\alpha$:

$\displaystyle\widehat{\mathcal{F}}[\phi;t]$ $\displaystyle=$
$\displaystyle\langle\exp\left\\{\mathrm{i}\widehat{S}_{\text{int}}[\phi,\chi;t]\right\\}\rangle_{\chi}$
(61) $\displaystyle=$ $\displaystyle 1+\mathrm{i}\sum\limits_{a=\pm}a\langle
S_{\text{int}}[\phi^{a},\chi^{a};t]\rangle_{\chi}-\frac{1}{2}\sum\limits_{a,b=\pm}ab\langle
S_{\text{int}}[\phi^{a},\chi^{a};t]S_{\text{int}}[\phi^{b},\chi^{b};t]\rangle_{\chi}+\mathcal{O}(\alpha^{3})\,\,\,.$

In contrast to the system fields, where only Feynman and Dyson propagators are
allowed, evaluating the expectation value with respect to $\chi$ includes all
types of contractions and therefore also leads to the appearance of
negative/positive frequency Wightman propagators, such that we have to use the
following:

$\displaystyle\contraction{}{\chi}{{}^{+}_{x}}{\chi}\chi^{+}_{x}\chi^{+}_{y}$
$\displaystyle=$ $\displaystyle\langle
T\chi_{x}\chi_{y}\rangle\,=\,\Delta^{++}_{xy}\,=\,\Delta^{\rm
F}_{xy}\,=\,-\mathrm{i}\int_{k}\frac{e^{ik\cdot(x-y)}}{k^{2}+m^{2}-\mathrm{i}\epsilon}\,\,\,,$
(62)
$\displaystyle\contraction{}{\chi}{{}^{-}_{x}}{\chi}\chi^{-}_{x}\chi^{-}_{y}$
$\displaystyle=$
$\displaystyle\langle\tilde{T}\chi_{x}\chi_{y}\rangle\,=\,\Delta^{--}_{xy}\,=\,\Delta^{\rm
D}_{xy}\,=\,+\mathrm{i}\int_{k}\frac{e^{ik\cdot(x-y)}}{k^{2}+m^{2}+\mathrm{i}\epsilon}\,\,\,,$
(63)
$\displaystyle\contraction{}{\chi}{{}^{+}_{x}}{\chi}\chi^{+}_{x}\chi^{-}_{y}$
$\displaystyle=$
$\displaystyle\langle\chi_{y}\chi_{x}\rangle\,=\,\Delta^{+-}_{xy}\,=\,\Delta^{<}_{xy}\,=\,\int_{k}e^{ik\cdot(x-y)}2\pi\Theta(-k^{0})\delta(k^{2}+m^{2})\,\,\,,$
(64)
$\displaystyle\contraction{}{\chi}{{}^{-}_{x}}{\chi}\chi^{-}_{x}\chi^{+}_{y}$
$\displaystyle=$
$\displaystyle\langle\chi_{x}\chi_{y}\rangle\,=\,\Delta^{-+}_{xy}\,=\,\Delta^{>}_{xy}~{}=~{}\Delta^{<}_{yx}\,=\,(\Delta^{<})^{*}_{xy}\,\,\,.$
(65)

These propagators fulfil the greatest time equation for any positive integer
$n$:

$\displaystyle\sum_{a,b=\pm}ab(\Delta_{xy}^{ab})^{n}$ $\displaystyle=$
$\displaystyle 0\,\,\,.$ (66)

Using Eqs. (62)-(65), we can separately compute each term in Eq. (61) and find

$\displaystyle\langle S_{\text{int}}[\phi^{a},\chi^{a};t]\rangle_{\chi}$
$\displaystyle=$
$\displaystyle-\frac{\alpha}{2}m\int_{x}\phi_{x}^{a}\Delta^{\rm
F}_{xx}\,\,\,,$ (67) $\displaystyle\langle
S_{\text{int}}[\phi^{a},\chi^{a};t]S_{\text{int}}[\phi^{b},\chi^{b};t]\rangle_{\chi}$
$\displaystyle=$
$\displaystyle\frac{\alpha^{2}}{4}m^{2}\int_{xy}\phi_{x}^{a}\phi_{y}^{b}\left[\Delta^{\rm
F}_{xx}\Delta^{\rm F}_{yy}+2(\Delta^{ab}_{xy})^{2}\right]\,\,\,.$ (68)

Finally, substituting Eqs. (67) and (68) into Eq. (61), we obtain for the
Feynman-Vernon influence functional:

$\displaystyle\widehat{\mathcal{F}}[\phi;t]$ $\displaystyle=$ $\displaystyle
1-\mathrm{i}\frac{\alpha}{2}m\sum\limits_{a=\pm}a\int_{x}\phi_{x}^{a}\Delta^{\rm
F}_{xx}-\frac{\alpha^{2}}{8}m^{2}\sum\limits_{a,b=\pm}ab\int_{xy}\phi_{x}^{a}\phi_{y}^{b}\left[\Delta^{\rm
F}_{xx}\Delta^{\rm
F}_{yy}+2(\Delta^{ab}_{xy})^{2}\right]+\mathcal{O}(\alpha^{3})\,\,\,.\,\,\,\,\,\,\,\,\,\,\,\,$
(69)

###  III.2 Density matrix elements

Now we are going compute the elements of the density matrices $\rho_{0,0}$,
$\rho_{1,0}$ and $\rho_{1,1}$ under different initial conditions. For this, we
will use Eqs. (II.4) and (69).

####  III.2.1 Vacuum state $\rho_{0,0}$

We begin with the system-vacuum state under the assumption that only itself
was non-vanishing at the initial time. In this case, Eq. (II.4) becomes

$\displaystyle\rho_{0,0}(t)$ $\displaystyle=$
$\displaystyle\rho_{0,0}(0)\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\widehat{\mathcal{F}}[\phi;t]\,\,\,.$
(70)

Substituting Eq. (69) and contracting the system fields, leads us to

$\displaystyle\rho_{0,0}(t)$ $\displaystyle\approx$
$\displaystyle\rho_{0,0}(0)\bigg{\\{}1-\frac{\alpha^{2}}{8}m^{2}\sum\limits_{a=\text{F},\text{D}}\int_{zz^{\prime}}\left[\Delta^{a}_{zz}\Delta^{a}_{z^{\prime}z^{\prime}}+2(\Delta^{a}_{zz^{\prime}})^{2}\right]D^{a}_{zz^{\prime}}\bigg{\\}}$
(71) $\displaystyle=:$ $\displaystyle\rho_{0,0}(0)[1+\aleph(t)]\,\,\,,$

where in the last line we summarized all terms beyond zeroth order as
$\aleph(t)$. As expected, the considered corrections to $\rho_{0,0}(0)$ are
disconnected bubble diagrams, see Fig. 2. It is obvious that this must also
hold for all terms beyond second order in $\alpha$ since the initial vacuum
density matrix in Eq. (71) can never connect to the Feynman-Vernon influence
functional, leaving the latter’s fields to contract only with themselves.
Since we will also consider other cases, in which such disconnected diagrams
will appear, we introduced the short-hand notation $\aleph(t)$, which can
explicitly be evaluated as

$\displaystyle\aleph(t)$ $\displaystyle\approx$
$\displaystyle-\frac{\alpha^{2}m^{2}}{2}\mathrm{vol}\left\\{\frac{\sin^{2}(Mt/2)}{M^{3}}\Delta^{\rm
F}_{zz}\Delta^{\rm
F}_{z^{\prime}z^{\prime}}+\frac{1}{2}\int_{\mathbf{k}\mathbf{q}}\frac{\sin^{2}[(E^{\phi}_{\mathbf{k}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{k}+\mathbf{q}})t/2]}{E^{\phi}_{\mathbf{k}}E^{\chi}_{\mathbf{q}}E^{\chi}_{\mathbf{k}+\mathbf{q}}(E^{\phi}_{\mathbf{k}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{k}+\mathbf{q}})^{2}}\right\\}\,\,\,,$
(72)

where the factor $\mathrm{vol}$ denotes the volume of our spatial domain of
integration. The first term on the right-hand side of Eq. (72) corresponds to
the tadpole diagram in Fig. 2(a), while the second term corresponds to
$\phi$-particle creation represented in Fig. 2(b), where we have to keep in
mind that the closed $\chi$-loops only appear due to us tracing out all
environmental degrees of freedom, which disallows external legs of $\chi$.
Particle annihilation does not appear here since, in the present case, it is a
phenomenon that arises at fourth order in $\alpha$. For annihilation to take
place, first, a particle has to be created within the time interval $[0,t)$
before it can decay into the $\chi$-vacuum before or at the final time.

(a)

(b)

Figure 2:  Diagrammatic representation of the terms contributing to
$\aleph(t)$, as given in Eq. (72); solid/dashed lines represent
$\phi$/$\chi$-propagators.

Next, we consider the case in which there existed a superposed single particle
at the initial time represented by the matrix $\rho_{1,1}(0)$, and see how
this modifies the vacuum density matrix at the final time $t$. For this, we
write down the respective version of Eq. (II.4) and substitute Eq. (69):

$\displaystyle\rho_{0,0}(t)$ $\displaystyle=$
$\displaystyle\lim_{y^{0(\prime)}\,\to\,0^{-}}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
(73)
$\displaystyle\phantom{}\times\int_{\mathbf{y}\mathbf{y}^{\prime}}e^{\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\widehat{\mathcal{F}}[\phi;t]\phi^{+}(y)\phi^{-}(y^{\prime})$
$\displaystyle\approx$
$\displaystyle\frac{\alpha^{2}}{4}m^{2}\lim_{y^{0(\prime)}\,\to\,0^{-}}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
$\displaystyle\times\int_{\mathbf{y}\mathbf{y}^{\prime}}e^{\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}\int_{zz^{\prime}}\big{[}\Delta^{\mathrm{F}}_{zz}\Delta^{\mathrm{F}}_{z^{\prime}z^{\prime}}+2(\Delta^{+-}_{zz^{\prime}})^{2}\big{]}D^{\mathrm{F}}_{zy}D^{\mathrm{D}}_{z^{\prime}y^{\prime}}\,\,\,.$

In contrast to the previous case in Eq. (71), the fields coming from the
Feynman-Vernon influence functional can now also contract with other fields,
namely those coming from the initial density matrix. Therefore, we find:

$\displaystyle\rho_{0,0}(t)$ $\displaystyle\approx$ $\displaystyle\frac{\alpha
m^{2}}{8}\Bigg{\\{}\frac{\sin^{2}(Mt/2)}{M^{4}}\Delta^{\mathrm{F}}_{zz}\Delta^{\mathrm{F}}_{z^{\prime}z^{\prime}}\rho_{1,1}(\mathbf{0},\mathbf{0};0)$
(74) $\displaystyle\phantom{\approx\frac{\alpha
m^{2}}{8}}\left.+\frac{1}{2}\int_{\mathbf{k}\mathbf{q}}\frac{\sin^{2}[(E^{\phi}_{\mathbf{k}}-E^{\chi}_{\mathbf{q}}-E^{\chi}_{\mathbf{k}-\mathbf{q}})t/2]}{(E^{\phi}_{\mathbf{k}})^{2}E^{\chi}_{\mathbf{q}}E^{\chi}_{\mathbf{k}-\mathbf{q}}(E^{\phi}_{\mathbf{k}}-E^{\chi}_{\mathbf{q}}-E^{\chi}_{\mathbf{k}-\mathbf{q}})^{2}}\rho_{1,1}(\mathbf{k},\mathbf{k};0)\right\\}\,\,\,.$

These terms are diagrammatically represented in Fig. 3. The first one on the
right-hand side of Eq. (74) corresponds to the two tadpoles in Fig. 3(a) and
the second one describes particle decay of the two superposed states into the
traced-out $\chi$-fields as depicted in Fig. 3(b). We see that the right-hand
side of Eq. (74), including the time-dependent divergences, vanishes for
$t\equiv 0$, which is in concordance with the initial conditions.

(a)

(b)

Figure 3:  Diagrammatic representation of the terms contributing to
$\rho_{0,0}(t)$ if only $\rho_{1,1}(0)$ was initially non-vanishing, as given
in Eq. (74); a crossed box depicts an insertion of an initial density matrix
element, here of the matrix $\rho_{1,1}(0)$.

####  III.2.2 Single-particle state $\rho_{1,0}$

We now move away from the vacuum state and have a look at the evolution of a
single particle that is not in a momentum superposition. At first, we assume
that the non-superposed single-particle state was the only one present at the
initial time. Under this assumption, we find for Eq. (II.4) after substituting
Eq. (69):

$\displaystyle\rho_{1,0}(\mathbf{p};t)$ $\displaystyle=$
$\displaystyle\lim_{\begin{subarray}{c}x^{0}\,\to\,t^{+}\\\
y^{0}\,\to\,0^{-}\end{subarray}}\int
d\Pi_{\mathbf{k}}\rho_{1,0}(\mathbf{k};0)\int_{\mathbf{x}\mathbf{y}}e^{-\mathrm{i}\mathbf{p}\cdot\mathbf{x}+\mathrm{i}\mathbf{k}\cdot\mathbf{y}}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\phi^{+}(x)\widehat{\mathcal{F}}[\phi;t]\phi^{+}(y)$
(75) $\displaystyle\approx$
$\displaystyle\lim_{\begin{subarray}{c}x^{0}\,\to\,t^{+}\\\
y^{0}\,\to\,0^{-}\end{subarray}}\int
d\Pi_{\mathbf{k}}\rho_{1,0}(\mathbf{k};0)\int_{\mathbf{x}\mathbf{y}}e^{-\mathrm{i}\mathbf{p}\cdot\mathbf{x}+\mathrm{i}\mathbf{k}\cdot\mathbf{y}}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}$
$\displaystyle\times\bigg{\\{}[1+\aleph(t)]D^{\rm
F}_{xy}-\frac{\alpha^{2}}{4}m^{2}\int_{zz^{\prime}}\Big{(}\Delta^{\mathrm{F}}_{zz}\Delta^{\mathrm{F}}_{z^{\prime}z^{\prime}}+2(\Delta^{\mathrm{F}}_{zz^{\prime}})^{2}\Big{)}D^{\mathrm{F}}_{xz}D^{\mathrm{F}}_{z^{\prime}y}\bigg{\\}}\,\,\,.$

Notice how the term $\aleph(t)$, that we already computed in Eq. (72), appears
here as a dressing of the system propagator. Evaluating Eq. (75) leads us to

$\displaystyle\rho_{1,0}(\mathbf{p};t)$ $\displaystyle\approx$ $\displaystyle
e^{-\mathrm{i}E^{\phi}_{\mathbf{p}}t}\rho_{1,0}(\mathbf{p};0)\Bigg{\\{}1+\aleph(t)+\mathrm{i}\frac{\alpha^{2}m^{2}}{16E^{\phi}_{\mathbf{p}}}\sum\limits_{s=\pm}\int_{\mathbf{q}}\frac{1}{E^{\chi}_{\mathbf{q}}E^{\chi}_{\mathbf{p}-\mathbf{q}}(sE^{\phi}_{\mathbf{p}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{p}-\mathbf{q}})}$
(76)
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\Bigg{[}t+\frac{\mathrm{i}}{sE^{\phi}_{\mathbf{p}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{p}-\mathbf{q}}}\Bigg{(}1-e^{-\mathrm{i}(sE^{\phi}_{\mathbf{p}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{p}-\mathbf{q}})t}\Bigg{)}\Bigg{]}\Bigg{\\}}$
$\displaystyle+\frac{\alpha^{2}(2\pi)^{3}}{8}\frac{m^{2}}{M^{3}}\delta^{(3)}(\mathbf{p})\Delta^{\mathrm{F}}_{zz}\Delta^{\mathrm{F}}_{z^{\prime}z^{\prime}}\left(1-e^{-\mathrm{i}Mt}\right)^{2}\rho_{1,0}(\mathbf{0};0)\,\,\,.$

In the first two lines of Eq. (76) we observe corrections to the unitary
evolution of $\rho_{1,0}$ in form of disconnected diagrams corresponding to
$\aleph(t)$, see Fig. 2, and a loop correction to the propagator as depicted
in Fig. 4(a). However, the third line is only depending on the initial
$\rho_{1,0}$ for a particle at nil momentum (in the chosen, fixed reference
frame), and can be interpreted as a $\phi$-particle decaying within the time
interval $[0,t)$ and another one with momentum $\mathbf{p}$ being produced
before or at the final time $t$. This process is represented diagrammatically
in Fig. 4(b).

(a)

(b)

Figure 4:  Diagrammatic representation of the terms contributing to
$\rho_{1,0}(t)$ if only $\rho_{1,0}(0)$ was initially non-vanishing, as given
in Eq. (76); a crossed box depicts an insertion of an element of the initial
density matrix $\rho_{1,0}(0)$.

Next, we assume that the single particle was initially in a momentum
superposition. In this case, we have

$\displaystyle\rho_{1,0}(\mathbf{p};t)$ $\displaystyle=$
$\displaystyle\lim_{\begin{subarray}{c}x^{0}\,\to\,t^{+}\\\
y^{0(\prime)}\,\to\,0^{-}\end{subarray}}\int d\Pi_{\mathbf{k}}\int
d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)\int_{\mathbf{x}\mathbf{y}\mathbf{y}^{\prime}}e^{-\mathrm{i}\mathbf{p}\cdot\mathbf{x}+\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}$
(77)
$\displaystyle\times\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\phi^{+}(x)\widehat{\mathcal{F}}[\phi;t]\phi^{+}(y)\phi^{-}(y^{\prime})$
$\displaystyle\approx$
$\displaystyle\mathrm{i}\frac{\alpha}{2}m\lim_{\begin{subarray}{c}x^{0}\,\to\,t^{+}\\\
y^{0(\prime)}\,\to\,0^{-}\end{subarray}}\int d\Pi_{\mathbf{k}}\int
d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
$\displaystyle\times\int_{\mathbf{x}\mathbf{y}\mathbf{y}^{\prime}}e^{-\mathrm{i}\mathbf{p}\cdot\mathbf{x}+\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}\int_{z}D^{\rm
F}_{xy}D^{\rm D}_{zy^{\prime}}\Delta^{\rm F}_{zz}\,\,\,.$

Note that here, in contrast to the previously considered cases, no terms of
zeroth or second order in $\alpha$ appear, but only the one of first order.
Evaluating Eq. (77) leaves us with a single term:

$\displaystyle\rho_{1,0}(\mathbf{p};t)$ $\displaystyle\approx$
$\displaystyle\mathrm{i}\frac{\alpha m}{4M^{2}}\Delta^{\rm
F}_{zz}\left(1-e^{\mathrm{i}Mt}\right)e^{-\mathrm{i}E^{\phi}_{\mathbf{p}}t}\rho_{1,1}(\mathbf{p},\mathbf{0};0)\,\,\,,$
(78)

which is diagrammatically represented in Fig. 5 and corresponds to decoherence
due to the decay of one of the two superposed momentum states into
$\chi$-fields.

(a)  Figure 5:  Diagrammatic representation of the term contributing to
$\rho_{1,0}(t)$ if only $\rho_{1,1}(0)$ was initially non-vanishing, as given
in Eq. (78)

####  III.2.3 Superposed single-particle state $\rho_{1,1}$

Finally, we look at the evolution of a single particle in a momentum
superposition under the assumption that it already and only it existed at the
initial time. In this case, we can use Eq. (52), substitute Eq. (69), and find

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}\,^{\prime};t)$ $\displaystyle=$
$\displaystyle\lim_{\begin{subarray}{c}x^{0(\prime)}\,\to\,t^{+}\\\
y^{0(\prime)}\,\to\,0^{-}\end{subarray}}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
(79)
$\displaystyle\times\int_{\mathbf{x}\mathbf{x}^{\prime}\mathbf{y}\mathbf{y}\,^{\prime}}e^{-\mathrm{i}(\mathbf{p}\cdot\mathbf{x}-\mathbf{p}\,^{\prime}\cdot\mathbf{x}^{\prime})+\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}\,^{\prime})}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{x^{0\prime},E^{\phi}_{\mathbf{p}\,^{\prime}}}^{*}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}$
$\displaystyle\times\int\mathcal{D}\phi^{\pm}e^{\mathrm{i}\widehat{S}_{\phi}[\phi]}\phi^{+}(x)\phi^{-}(x^{\prime})\widehat{\mathcal{F}}[\phi;t]\phi^{+}(y)\phi^{-}(y^{\prime})$
$\displaystyle\approx$
$\displaystyle\lim_{\begin{subarray}{c}x^{0(\prime)}\,\to\,t^{+}\\\
y^{0(\prime)}\,\to\,0^{-}\end{subarray}}\int
d\Pi_{\mathbf{k}}d\Pi_{\mathbf{k}^{\prime}}\rho_{1,1}(\mathbf{k},\mathbf{k}^{\prime};0)$
$\displaystyle\times\int_{\mathbf{x}\mathbf{x}^{\prime}\mathbf{y}\mathbf{y}^{\prime}}e^{-\mathrm{i}(\mathbf{p}\cdot\mathbf{x}-\mathbf{p}^{\prime}\cdot\mathbf{x}^{\prime})+\mathrm{i}(\mathbf{k}\cdot\mathbf{y}-\mathbf{k}^{\prime}\cdot\mathbf{y}^{\prime})}\partial_{x^{0},E^{\phi}_{\mathbf{p}}}\partial_{x^{0\prime},E^{\phi}_{\mathbf{p}^{\prime}}}^{*}\partial_{y^{0},E^{\phi}_{\mathbf{k}}}^{*}\partial_{y^{0\prime},E^{\phi}_{\mathbf{k}^{\prime}}}$
$\displaystyle\times\bigg{\\{}[1+\aleph(t)]D^{\rm F}_{xy}D^{\rm
D}_{x^{\prime}y^{\prime}}-\frac{\alpha^{2}}{4}m^{2}\int_{zz^{\prime}}\bigg{[}\Big{(}\Delta^{\rm
F}_{zz}\Delta^{\rm F}_{z^{\prime}z^{\prime}}+2(\Delta^{\rm
D}_{zz^{\prime}})^{2}\Big{)}D^{\rm F}_{xy}D^{\rm D}_{x^{\prime}z}D^{\rm
D}_{z^{\prime}y^{\prime}}$ $\displaystyle+\Big{(}\Delta^{\rm
F}_{zz}\Delta^{\rm F}_{z^{\prime}z^{\prime}}+2(\Delta^{\rm
F}_{zz^{\prime}})^{2}\Big{)}D^{\rm D}_{x^{\prime}y^{\prime}}D^{\rm
F}_{xz^{\prime}}D^{\rm F}_{zy}\bigg{]}\bigg{\\}}\,\,\,.$

Note the reappearance of $\aleph(t)$, which is now dressing the two
$\phi$-propagators, each of which is corresponding to either one of the two
superposed states. Evaluating Eq. (79), we find terms that show a structural
resemblance to those in Eq. (76):

$\displaystyle\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};t)$
$\displaystyle\approx$ $\displaystyle
e^{-\mathrm{i}(E^{\phi}_{\mathbf{p}}-E^{\phi}_{\mathbf{p}^{\prime}})t}\rho_{1,1}(\mathbf{p},\mathbf{p}^{\prime};0)$
(80)
$\displaystyle\times\Bigg{\\{}1+\aleph(t)+\mathrm{i}\frac{\alpha^{2}m^{2}}{16}\sum\limits_{s=\pm}\int_{\mathbf{q}}\frac{1}{E^{\chi}_{\mathbf{q}}}\Bigg{[}\frac{1}{E^{\phi}_{\mathbf{p}}E^{\chi}_{\mathbf{p}-\mathbf{q}}(sE^{\phi}_{\mathbf{p}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{p}-\mathbf{q}})}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\Bigg{(}t+\frac{\mathrm{i}}{sE^{\phi}_{\mathbf{p}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{p}-\mathbf{q}}}\Bigg{(}1-e^{-\mathrm{i}(sE^{\phi}_{\mathbf{p}}+E^{\chi}_{\mathbf{q}}+E^{\chi}_{\mathbf{p}-\mathbf{q}})t}\Bigg{)}\Bigg{)}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-(\mathbf{p}\longleftrightarrow\mathbf{p^{\prime}})^{\ast}\Bigg{]}\Bigg{\\}}$
$\displaystyle+\frac{\alpha^{2}m^{2}(2\pi)^{3}}{8M^{3}}\Delta^{\rm
F}_{zz}\Delta^{\rm
F}_{z^{\prime}z^{\prime}}\left[\delta^{(3)}(\mathbf{p}^{\prime})\left(1-e^{\mathrm{i}Mt}\right)^{2}e^{-\mathrm{i}E^{\phi}_{\mathbf{p}}t}\rho_{1,1}(\mathbf{p},\mathbf{0};0)+(\mathbf{p}\longleftrightarrow\mathbf{p^{\prime}})^{\ast}\right]\,\,\,.$

Again we find corrections to the unitary evolution, in this case, of
$\rho_{1,1}$. Those are given in the second to fourth lines of Eq. (80) and
comprise disconnected diagrams, as depicted in Fig. 2, and loop corrections to
the two propagators corresponding to the diagrams in Figs. 6(a) and (b). In
the last line of Eq. (80) we only have terms proportional to initial density
matrix elements with one of the two superposed states representing a
motionless particle (in the chosen, fixed reference frame), and which
correspond to the decay and creation of particles, as shown in Figs. 6(c) and
(d). The processes described by Eq. (80) result in coherence changes and/or
phase shifts.

(a)

(b)

(c)

(d)

Figure 6:  Diagrammatic representation of the terms contributing to
$\rho_{1,1}(t)$ if only $\rho_{1,1}(0)$ was initially non-vanishing, as given
in Eq. (80)

##  IV Conclusions and outlook

The theory of open quantum systems has various applications in many areas of
physics. It often makes use of density matrices as its tool of choice for
describing open quantum systems and their environments. Tracing out the
environmental degrees of freedom leads to reduced density matrices, whose time
evolution can be investigated with master equations. However, in many cases,
analytically solving a master equation poses an intricate or even impossible
task especially if Markovianity is not assumed.
In the present article we used a previously presented path-integral–based
formalism in order to derive, from first principles, a formula that allows us
to directly compute reduced density matrices without having to solve a master
equation. We only used the Born approximation, but allowed for non-Markovian
open quantum systems. In addition, the derived formula can be applied to any
number of superposed or non-superposed particles in Fock space. We chose to
present the formula for density matrix elements in a momentum basis, and for a
real scalar field open system interacting with another real scalar field as
its environment. However, it is also possible to extend this formalism to
other bases like position, and other field species.
First, we presented the description for density matrices in Fock space and
introduced the necessary techniques from non-equilibrium quantum field theory,
namely the Feynman-Vernon influence functional, which is based on the
Schwinger-Keldysh formalism, and TFD. Next, we stated the well-known fact
that, in the Schrödinger picture TFD, the quantum Liouville equation can be
expressed in a Schrödinger-like form, which in turn can be solved for a
general Hamiltonian. We then applied the LSZ-like reduction formalism that was
developed in Ref. [16] in order to derive the equation that provides a new way
of directly computing reduced density matrix elements. Finally, we considered
a simple example for an open system and its environment, and applied the
formula for computing a few selected density matrix elements for different
initial conditions. In this way, we demonstrated the practicability of the
presented formalism, and observed open quantum dynamical effects like
decoherence. However, we kept the time-dependent divergences, whose curing by
renormalization will be discussed in a future article.
Even though we presented the formula in a field theoretical setting, it could
well be adapted to non-relativistic quantum physics and find a variety of
applications there. In its present form, the formalism discussed here can be
applied for providing new insights into phenomena like gravitational or scalar
field-induced decoherence, as well as for entirely new investigations in the
field theory of open quantum systems.

###### Acknowledgements.

We thank P. Millington and J. Minář for useful comments and helpful
discussions. This work was supported by the Austrian Science Fund (FWF): P
34240-N.

## References

  * [1] H.-P. Breuer and F. Petruccione, _The Theory of Open Quantum Systems_. Oxford University Press, Oxford, 2002.
  * [2] E. A. Calzetta and B.-L. Hu, _Nonequilibrium Quantum Field Theory_. Cambridge University Press, Cambridge UK, 2008.
  * [3] M. Schlosshauer, _Decoherence and the Quantum-To-Classical Transition_. Springer-Verlag Berlin Heidelberg, 2007.
  * [4] H. Carmichael, _An Open Systems Approach to Quantum Optics: Lectures Presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_. Springer Berlin Heidelberg, 1993.
  * [5] C. Gardiner and P. Zoller, _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_ , Springer Series in Synergetics. Springer, 2004.
  * [6] D. Walls and G. Milburn, _Quantum Optics_. Springer Berlin Heidelberg, 2008\.
  * [7] L. Aolita, F. de Melo and L. Davidovich, _Open-system dynamics of entanglement:a key issues review_ , _Reports on Progress in Physics_ 78 (2015) 042001.
  * [8] J. Goold, M. Huber, A. Riera, L. d. Rio and P. Skrzypczyk, _The role of quantum information in thermodynamics—a topical review_ , _Journal of Physics A: Mathematical and Theoretical_ 49 (2016) 143001.
  * [9] A. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert et al., _Positive Tensor Network Approach for Simulating Open Quantum Many-Body Systems_ , _Physical Review Letters_ 116 (2016) .
  * [10] J. Huber, P. Kirton, S. Rotter and P. Rabl, _Emergence of PT-symmetry breaking in open quantum systems_ , _SciPost Physics_ 9 (2020) .
  * [11] J. F. Koksma, T. Prokopec and M. G. Schmidt, _Decoherence in an interacting quantum field theory: The vacuum case_ , _Phys. Rev. D_ 81 (2010) 065030.
  * [12] J. F. Koksma, T. Prokopec and M. G. Schmidt, _Decoherence in an interacting quantum field theory: Thermal case_ , _Phys. Rev. D_ 83 (2011) 085011.
  * [13] L. M. Sieberer, M. Buchhold and S. Diehl, _Keldysh field theory for driven open quantum systems_ , _Reports on Progress in Physics_ 79 (2016) 096001.
  * [14] J. Marino and S. Diehl, _Quantum dynamical field theory for nonequilibrium phase transitions in driven open systems_ , _Phys. Rev. B_ 94 (2016) 085150.
  * [15] A. Baidya, C. Jana, R. Loganayagam and A. Rudra, _Renormalization in open quantum field theory. Part I. Scalar field theory_ , _JHEP_ 11 (2017) 204 [1704.08335].
  * [16] C. Burrage, C. Käding, P. Millington and J. Minář, _Open quantum dynamics induced by light scalar fields_ , _Phys. Rev. D_ 100 (2019) 076003 [1812.08760].
  * [17] S. Nagy and J. Polonyi, _Renormalizing Open Quantum Field Theories_ , _Universe_ 8 (2022) 127 [2012.13811].
  * [18] C. Jana, _Aspects of open quantum field theory_ , Ph.D. thesis, Tata Inst., 2021.
  * [19] H. C. Fogedby, _Field theoretical approach to open quantum systems and the Lindblad equation_ , 2202.05203.
  * [20] F. Lombardo and F. D. Mazzitelli, _Coarse graining and decoherence in quantum field theory_ , _Phys. Rev. D_ 53 (1996) 2001.
  * [21] F. C. Lombardo and D. L. Nacir, _Decoherence during inflation: The generation of classical inhomogeneities_ , _Phys. Rev. D_ 72 (2005) 063506.
  * [22] F. C. Lombardo, _Influence functional approach to decoherence during inflation_ , _Braz. J. Phys._ 35 (2005) 391 [gr-qc/0412069].
  * [23] D. Boyanovsky, _Effective field theory during inflation: Reduced density matrix and its quantum master equation_ , _Phys. Rev._ D92 (2015) 023527 [1506.07395].
  * [24] D. Boyanovsky, _Effective field theory during inflation. II. Stochastic dynamics and power spectrum suppression_ , _Phys. Rev._ D93 (2016) 043501 [1511.06649].
  * [25] D. Boyanovsky, _Fermionic influence on inflationary fluctuations_ , _Phys. Rev._ D93 (2016) 083507 [1602.05609].
  * [26] D. Boyanovsky, _Imprint of entanglement entropy in the power spectrum of inflationary fluctuations_ , _Phys. Rev._ D98 (2018) 023515 [1804.07967].
  * [27] C. P. Burgess, R. Holman, G. Tasinato and M. Williams, _EFT beyond the horizon: stochastic inflation and how primordial quantum fluctuations go classical_ , _Journal of High Energy Physics_ 2015 (2015) 90.
  * [28] T. J. Hollowood and J. I. McDonald, _Decoherence, discord, and the quantum master equation for cosmological perturbations_ , _Phys. Rev. D_ 95 (2017) 103521.
  * [29] H. W. Yu, J. Zhang, H.-w. Yu and J.-l. Zhang, _Understanding Hawking radiation in the framework of open quantum systems_ , _Phys. Rev. D_ 77 (2008) 024031 [0806.3602].
  * [30] F. C. Lombardo and G. J. Turiaci, _Dynamics of an Acoustic Black Hole as an Open Quantum System_ , _Phys. Rev. D_ 87 (2013) 084028 [1208.0198].
  * [31] C. Jana, R. Loganayagam and M. Rangamani, _Open quantum systems and Schwinger-Keldysh holograms_ , _JHEP_ 07 (2020) 242 [2004.02888].
  * [32] K. Agarwal and N. Bao, _Toy model for decoherence in the black hole information problem_ , _Phys. Rev. D_ 102 (2020) 086017.
  * [33] G. Kaplanek and C. P. Burgess, _Qubits on the Horizon: Decoherence and Thermalization near Black Holes_ , _JHEP_ 01 (2021) 098 [2007.05984].
  * [34] C. P. Burgess, R. Holman and G. Kaplanek, _Quantum Hotspots: Mean Fields, Open EFTs, Nonlocality and Decoherence Near Black Holes_ , 2106.10804.
  * [35] G. Kaplanek, C. P. Burgess and R. Holman, _Qubit heating near a hotspot_ , _JHEP_ 08 (2021) 132 [2106.10803].
  * [36] N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, _Quarkonium suppression in heavy-ion collisions: an open quantum system approach_ , _Phys. Rev._ D96 (2017) 034021 [1612.07248].
  * [37] N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, _Heavy quarkonium suppression in a fireball_ , _Phys. Rev._ D97 (2018) 074009 [1711.04515].
  * [38] Y. Akamatsu, _Quarkonium in quark–gluon plasma: Open quantum system approaches re-examined_ , _Prog. Part. Nucl. Phys._ 123 (2022) 103932 [2009.10559].
  * [39] W. A. De Jong, M. Metcalf, J. Mulligan, M. Płoskoń, F. Ringer and X. Yao, _Quantum simulation of open quantum systems in heavy-ion collisions_ , _Phys. Rev. D_ 104 (2021) 051501 [2010.03571].
  * [40] X. Yao, _Open quantum systems for quarkonia_ , _Int. J. Mod. Phys. A_ 36 (2021) 2130010 [2102.01736].
  * [41] N. Brambilla, M. A. Escobedo, M. Strickland, A. Vairo, P. Vander Griend and J. H. Weber, _Bottomonium production in heavy-ion collisions using quantum trajectories: Differential observables and momentum anisotropy_ , _Phys. Rev. D_ 104 (2021) 094049 [2107.06222].
  * [42] P. V. Griend, _Bottomonium observables in an open quantum system using the quantum trajectories method_ , _EPJ Web Conf._ 258 (2022) 05005 [2111.13520].
  * [43] X. Yao, _Quarkonium Suppression in the Open Quantum System Approach_ , in _19th International Conference on Hadron Spectroscopy and Structure_ , 1, 2022, 2201.07702.
  * [44] M. P. Blencowe, _Effective Field Theory Approach to Gravitationally Induced Decoherence_ , _Phys. Rev. Lett._ 111 (2013) 021302 [1211.4751].
  * [45] C. Anastopoulos and B. L. Hu, _A Master Equation for Gravitational Decoherence: Probing the Textures of Spacetime_ , _Class. Quant. Grav._ 30 (2013) 165007 [1305.5231].
  * [46] T. Oniga and C. H. T. Wang, _Quantum gravitational decoherence of light and matter_ , _Phys. Rev. D_ 93 (2016) 044027 [1511.06678].
  * [47] J. Minář, P. Sekatski and N. Sangouard, _Bounding quantum-gravity-inspired decoherence using atom interferometry_ , _Phys. Rev. A_ 94 (2016) 062111.
  * [48] J. Minář, P. Sekatski, R. Stevenson and N. Sangouard, _Testing unconventional decoherence models with atoms in optical lattices_ , 2016.
  * [49] A. Bassi, A. Großardt and H. Ulbricht, _Gravitational Decoherence_ , _Class. Quant. Grav._ 34 (2017) 193002 [1706.05677].
  * [50] L. Asprea, G. Gasbarri and A. Bassi, _Gravitational decoherence: A general nonrelativistic model_ , _Phys. Rev. D_ 103 (2021) 104041 [1905.01121].
  * [51] L. Asprea, A. Bassi, H. Ulbricht and G. Gasbarri, _Gravitational decoherence and the possibility of its interferometric detection_ , _Phys. Rev. Lett._ 126 (2021) 200403.
  * [52] L. Asprea and G. Gasbarri, _Gravitational decoherence: A nonrelativistic spin 1/2 fermionic model_ , _Phys. Rev. D_ 104 (2021) 024043.
  * [53] M. Lagouvardos and C. Anastopoulos, _Gravitational decoherence of photons_ , _Classical and Quantum Gravity_ 38 (2021) 115012.
  * [54] C. Anastopoulos and B.-L. Hu, _Gravitational Decoherence: A Thematic Overview_ , 2111.02462.
  * [55] I. Pikovski, M. Zych, F. Costa and Č. Brukner, _Universal decoherence due to gravitational time dilation_ , _Nature Phys._ 11 (2015) 668 [1311.1095].
  * [56] I. Pikovski, M. Zych, F. Costa and Č. Brukner, _Time dilation in quantum systems and decoherence_ , _New Journal of Physics_ 19 (2017) 025011.
  * [57] C. Burrage, C. Käding, P. Millington and J. Minář, _Influence functionals, decoherence and conformally coupled scalars_ , _J. Phys. Conf. Ser._ 1275 (2019) 012041 [1902.09607].
  * [58] C. Käding, _Astro- and Quantum Physical Tests of Screened Scalar Fields_ , Ph.D. thesis, University of Nottingham, Nottingham NG7 2RD, UK, 10, 2019\.  1910.05738.
  * [59] H. Lehmann, K. Symanzik and W. Zimmermann, _On the formulation of quantized field theories_ , _Nuovo Cim._ 1 (1955) 205.
  * [60] Y. Takahasi and H. Umezawa, _Thermo field dynamics_ , _Collect. Phenom._ 2 (1975) 55.
  * [61] T. Arimitsu and H. Umezawa, _A General Formulation of Nonequilibrium Thermo Field Dynamics_ , _Prog. Theor. Phys._ 74 (1985) 429.
  * [62] T. Arimitsu and H. Umezawa, _Non-Equilibrium Thermo Field Dynamics_ , _Prog. Theor. Phys._ 77 (1987) 32.
  * [63] F. C. Khanna, A. P. C. Malbouisson, J. M. C. Malbouisson and A. E. Santana, _Thermal Quantum Field Theory: Algebraic Aspects and Applications_. World Scientific, Singapore, 2009.
  * [64] J. S. Schwinger, _Brownian Motion of a Quantum Oscillator_ , _J. Math. Phys._ 2 (1961) 407.
  * [65] L. V. Keldysh, _Diagram technique for nonequilibrium processes_ , _Zh. Eksp. Teor. Fiz._ 47 (1964) 1515.
  * [66] R. P. Feynman and F. L. Vernon, _The theory of a general quantum system interacting with a linear dissipative system_ , _Annals of physics_ 24 (1963) 118.
  * [67] M. Le Bellac, _Thermal field theory_. Cambridge University Press, 1996.
  * [68] E. R. Loubenets and C. Käding, _Specifying the Unitary Evolution of a Qudit for a General Nonstationary Hamiltonian via the Generalized Gell-Mann Representation_ , _Entropy_ 22 (2020) .
  * [69] P. Millington and A. Pilaftsis, _Perturbative nonequilibrium thermal field theory_ , _Phys. Rev._ D88 (2013) 085009 [1211.3152].
  * [70] P. Millington and A. Pilaftsis, _Perturbative Non-Equilibrium Thermal Field Theory to all Orders in Gradient Expansion_ , _Phys. Lett._ B724 (2013) 56 [1304.7249].
  * [71] G. C. Wick, _The evaluation of the collision matrix_ , _Phys. Rev._ 80 (1950) 268.
  * [72] C. Käding, P. Millington and M. Pitschmann, _to be published_ , 22xx.xxxxx.