A new method for directly computing reduced density matrices

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.

In the theory of open quantum systems we usually deal with a system and its surrounding envi-discrete qudit systems, using the Schrödinger-like form of the quantum Liouville equation in TFD in combination with the results from Ref. [71] 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 ρ ϕ (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 ϕ and the environmental field χ. More specifically, we apply the derived formula in order to find explicit expressions for the associated vacuum and single-particle matrix elements of ρ ϕ (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 ϕ as the system and another real scalar field χ as its surrounding environment. In order to better motivate the system-environment split, we may assume that ϕ has a constant mass M much larger than χ's constant mass m, such that m/M ≪ 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 A-II C we will first present the necessary mathematical and conceptual preliminaries before finally combining them for the actual derivation in Sec. II D.

A. 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 bŷ where the cases i = 0 or j = 0 correspond to the static vacuum state |0⟩ or ⟨0|, respectively, each density matrix element is given by and we use Later we will also make use of We note that, as usual, the following has to hold: Physically these density matrix elements can be interpreted as follows: ρ 0;0 describes a 0-particle or vacuum state, ρ i;0 or ρ 0;j are correlations between the vacuum and i-or j-particle states, while ρ i;j represents correlations between i-and j-particle states.
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 which is solved byρ if the Hamiltonian is time-independent, or bŷ if the Hamiltonian is time-dependent due to an external source, where the index S labels objects in the Schrödinger picture, and T andT stand for time-ordering and anti-time-ordering, respectively.
Note that Eqn. (9) is more general than the solution in Eqn. (8) and recovers the latter ifĤ S (t) is constant.
When describing an open system in field theory, the Feynman-Vernon influence functional [69] based on the Schwinger-Keldysh closed-time-path formalism [67,68] 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 initial = 0 and a final time t final = t. Physically the 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 Since we are only interested in the evolution of the reduced density operatorρ ϕ (t), we trace out the environmental degrees of freedom:ρ Projecting this expression into a field basis living on the closed-time-path, we can express the reduced density functional as where ± indicates a dependence on both +-and −-type operators, subscript t labels the time slice on which the field eigenstates have to be taken, and Furthermore, taking into account Eq. (10), the reduced density functional at time t can be expressed where the time translation is given in terms of the so-called influence functional (IF) propagator, see e.g. Ref. [2], The latter contains two path integrals over an effective action where indicates functionals that depend on both field variables ϕ + and ϕ − . For example, the free action and the self-interaction of ϕ are given by However, the influence action S IF [ϕ; t] has a more involved definition by the Feynman-Vernon influence functional where S χ [χ; t] is the free action, S χ,int [χ; t] the self-interaction action for χ and S int [ϕ, χ; 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 Ω t := [0, t] × R 3 , which means they all have the structure In addition, note that we are working with two types of functional integrals: one, denoted d, considers fields over all of R 3 , but only at one particular time slice, while the other, denoted D, is over the whole Ω t . See Eq. (18) for an example in which both types appear together.
Finally, we define an expectation value with respect to χ as (cf. Refs. [2,16]) and find for the Feynman-Vernon influence functional: C. Thermo field dynamics Next, we look at thermo field dynamics (TFD) [63][64][65] (see also Ref. [66]) 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 where H ± are the Hilbert spaces corresponding to the ±-branches of the closed time path, respectively. Operators living on the closed time path can be expressed within TFD aŝ whereÎ is the unit operator and T means time reversal. States in a momentum basis in TFD can be reached from the doubled vacuum state via creation operatorŝ Consequently, the corresponding annihilators act likê In addition, we can construct a special state corresponding to the unit operator, see Ref. [64], which allows us to express the expectation value of an operator as Having all this, we can now use the fact that Eq. (7) can be re-written in a Schrödinger-like form: where . This can generally be solved bŷ Note that at time 0 the different pictures coincide and we therefore dropped the label S.
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 χ degrees of freedom from the total density operatorρ ϕχ : where the effective Hamiltonian is given by the usual unitary evolution of ϕ due to the Hamiltonianŝ H 0,S (t) +Ĥ int,S (t) and a non-unitarian H IF,S : We recall that, generally, Hamiltonian and action are related via H(t) = − ∂ ∂t S(t). Furthermore, if we consider that the effective action in Eq. (16) describes the full evolution of the reduced density matrix elements in the ϕ-basis, the interaction picture H ϕ,I must be its corresponding Hamiltonian, such thatĤ 0,I (t) +Ĥ int,I (t) represent S ϕ (t) + S ϕ,int (t) and H IF,I (t) the influence action given in Eq. (18). 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Ĥ 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. (32) into the interaction picture.
For this, we remind ourselves that an operator transforms likê from the Schrödinger into the interaction picture, and conveniently define a unitary evolution such that operators acting on the TFD doubled Hilbert space transform like The state introduced in Eq. (28) transforms like a base vector: However, it can straightforwardly be seen that U and U † act like unit operators on |1⟩⟩ S (and consequently on |1(t)⟩⟩ as well): which implies that the state is time-and picture-independent, and will therefore simply be denoted by |1⟩⟩ in what follows.
Using all this, Eq. (32) becomes in the interaction picture: Evaluating the partial derivative on the left-hand side of Eq. (39) leaves us with where H 0 :=Ĥ 0 ⊗Î −Î ⊗Ĥ 0 . Acting with H 0 on |1⟩⟩ gives nil, as can straightforwardly be seen.

D. 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 ρ ϕ (t) in a momentum basis for any occupation number in Fock space.
For notational convenience, we will from now on drop the index ϕ 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 ρ 1;1 (p; p ′ ; t), which represents a single particle in momentum space. From the resulting equation and its derivation it will be straightforward for us to infer expressions for the other density matrix elements.
Note that, as was pointed out in Refs. [72,73], the matrix element ρ 1;1 (p; p ′ ; t) is pictureindependent. In TFD, Eq. (44) can be expressed as Substituting Eq. (43), we are left with Next, we expand the density operator in Eq. (46) as in Eq. (1). However, we assume that the exact number of system particles and their correlations at the initial time are practically well-understood, such that we only have to consider the ρ 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, 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 which, using Eqs. (26) and (27) becomes After replacing the creation and annihilation operators by (cf. Ref. [16]) Here we introduced limits in which x 0(′) approach t from above and y 0(′) approach 0 from below in order to recover the correct time ordering, and started to use the short-hand notation ϕ x := ϕ(x) etc.
Translating the resulting expression in Eq. (50) into the path integral formalism, gives us where we used Using the definition in Eq. (18), we are finally led to What we just found is an equation that enables us to directly compute the momentum basis singleparticle elements of the reduced density matrix ρ ϕ (t) without having to explicitly solve any master equation. Using the assumption of a weak coupling would now allow us to perturbatively evaluate Eq. (53) 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. (53). We find: vanish, we recover Eq. (53).
When evaluating Eqs. (53) or (54) it is important to understand that, at least for the system field ϕ, only contractions of two +-or two −-fields are permitted, while for computing the influence functional via Eq. (22) all types of contractions of χ are allowed. That the former is true for singleparticle elements as in Eq. (53) was already pointed out in Ref. [16]. There the authors showed that the 2 × 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. (54) as the generalisation of the four-point function in Eq. (50) where exp {...} represents an exponentiation of a combination of ϕ-field operators. After using Wick's theorem [74], Eq. (55) reduces to a product of two-point functions. Taking into account the definition of the TFD vacuum state |0⟩⟩, it can be seen that the two-point functions that give the off-diagonal elements of the matrix propagator must vanish: The diagonal elements, however, give the usual Feynman and Dyson propagators

III. A SIMPLE TOY MODEL
Here we are going to apply the results from Sec. II, i.e. Eqs. (53) and (54), to a toy model for system ϕ and environment χ. 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.
After explicitly evaluating Eq. (54) for each case, we will repeatedly encounter terms with timedependent 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 [75].

A. The model
As an example we choose a two-scalar field model with no self-interactions, and the following actions: where m/M, α ≪ 1 and Ω t := [0, t] × R 3 . Even though it is rather simple, this model already includes some interesting features that we would like to illustrate. Taking into account selfinteractions 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 χ and different coupling terms can be found in Ref. [16].
Using Eq. (22), we expand the Feynman-Vernon influence functional up to second order in α: In contrast to the system fields, where only Feynman and Dyson propagators are allowed, evaluating the expectation value with respect to χ includes all types of contractions and therefore also leads to the appearance of negative/positive frequency Wightman propagators. Note that here we consider zero temperature and the trivial χ-vacuum, such that we have to use in what follows: These propagators fulfil the greatest time equation for any positive integer n: Using Eqs. (63)- (66), we can separately compute each term in Eq. (62) and find Finally, substituting Eqs. (68) and (69) into Eq. (62), we obtain for the Feynman-Vernon influence functional:

B. Density matrix elements
Now we are going to compute the elements of the density matrices ρ 0;0 , ρ 1;0 and ρ 1;1 under different initial conditions. For this, we will use Eqs. (54) and (70). How to evaluate these equations in detail was briefly outlined in Appendix D of Ref. [61].

Vacuum state ρ 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. (54) becomes Substituting Eq. (70) and contracting the system fields, leads us to where in the last line we summarized all terms beyond zeroth order as ℵ(t). As expected, the considered corrections to ρ 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 α since the initial vacuum density matrix in Eq. (72) 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 ℵ(t), which can explicitly be evaluated as The first term on the right-hand side of Eq. (73) corresponds to the tadpole diagram in Fig. 2(a), while the second term corresponds to the vacuum bubble represented in Fig. 2(b).
In contrast to the previous case in Eq. (72), 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: These terms are diagrammatically represented in Fig. 3. The first one on the right-hand side of Eq. (75) corresponds to the two tadpoles in Fig. 3(a) and the second one describes the bubble diagram given in Fig. 3(b). We see that the right-hand side of Eq. (75), including the time-dependent divergences, vanishes for t ≡ 0, which is in concordance with the initial conditions. 2. Correlation between single particle and vacuum ρ 1;0 We now move away from the vacuum state and have a look at the evolution of the correlation between a single particle and the vacuum. For this, we assume ρ 1;0 (0) to be non-vanishing. However, if such an initial correlation between a single particle and the vacuum existed, then we are required to also take into account non-vanishing ρ 0;0 (0) and ρ 1;1 (0). Furthermore, Eq. (6) tells us that also ρ 0;1 (0) ̸ = 0. Consequently, we find for Eq. (54) after substituting Eq. (70): system propagator. Evaluating Eq. (76) leads us to In the first line of Eq. (78) we see the creation of a single ϕ-particle with momentum p out of the vacuum, see Fig. 4(a). This process is possible since the χ-degrees of freedom were traced out and are therefore not visible even though they are actually responsible for this creation via the process χχ → ϕ. In the next two lines of Eq. (78) we observe corrections to the unitary evolution of ρ 1;0 in form of disconnected diagrams corresponding to ℵ(t), see Fig. 2, and a loop correction to the propagator as depicted in Fig. 4(b). The first term in the fourth line of Eq. (78) is only dependent on the initial ρ 1;0 for a particle at zero momentum in the chosen, fixed reference frame, and can be interpreted as a ϕ-particle decaying within the time interval [0, t) and another one with momentum p being produced before or at the final time t. This process is represented diagrammatically in between a unitarily evolving single particle and a particle decaying into the vacuum (via the process ϕ → χχ), and is represented in Fig. 4(f).
3. Single-particle state ρ 1,1 Finally, we look at the evolution of a single particle in momentum space under the assumption that it already and only it existed at the initial time. In this case, we can use Eq. (53), substitute Eq. (70), and find which becomes we find: Again we find corrections to the unitary evolution, including de-and recoherence, in this case, of

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 assumed a weak coupling, but allowed for non-Markovian open quantum systems. In addition, the derived formula can be applied to any number of 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 wellknown fact that, in 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.