Nonequilibrium Nonlinear Open Quantum Systems I. Functional Perturbative Analysis of a Weakly Anharmonic Oscillator

We introduce a functional perturbative method for treating weakly nonlinear systems such as a chain of coupled quantum oscillators with own baths. We demonstrate using this method to obtain the correlation functions of a quantum anharmonic oscillator interacting with a heat bath. These results are useful for studying the nonequilibrium physical processes of nonlinear quantum systems such as heat transfer or electron transport.


Introduction
Background Notes. These notes containing original research results were written in August 2014. We have reorganized the raw materials with minor revisions in wording for improved readability. For this Paper I we examine the nonequilibrium dynamics of a quantum Brownian oscillator [1,2,3,4,5,6] in a nonlinear potential. We use the functional perturbative method originated in [7] to study the open quantum system of a weakly anharmonic oscillator linearly coupled to a thermal field bath. A companion Paper II [8] studies two weakly nonlinearly coupled quantum oscillators each with its own bath. The correlation functions obtained in these notes are useful for the study of heat transfer or charge transport in these systems.

Nonequilibrium Dynamics of Nonlinear Open Quantum Systems
While linear analysis has occupied the attention of physicists and mathematicians for the longest time in theoretical science we know very well that large scale structures of Nature are more aptly and accurately depicted by nonlinear systems. Likewise, while the equilibrium condition is often assumed and most studied in statistical mechanics and thermodynamics, it cannot capture the real time dynamics of a many-body system. Nonequilibrium statistical physics is needed. Studies of the nonequilibrium dynamics of nonlinear systems in the last half century have spawned new fields, such as driven diffusive and dissipative systems. This is further enriched when noise and stochastic processes are incorporated -the interplay of nonlinearity and stochasticity reveals many interesting phenomena such as those pertaining to forms and growth, and created new fields, such as stochastic thermodynamics, etc. Studies of the nonequilibrium dynamics of quantum systems pose a new level of challenge, yet bring abundant rewards.
In fact, such studies rapidly became a necessity, as quantum devices are becoming smaller, the functional numbers (e.g., of photon) are becoming fewer, and new resources in quantum coherence and entanglement enabled innovations with wide ranging applications such as quantum information processing.
Studies of the nonequilibrium dynamics of quantum systems has its own long history -quantum Brownian motion (QBM) is probably the best known paradigm in the description of quantum dissipative and stochastic processes. For linear systems our understanding of QBM has progressed from Markovian processes a la the Lindblad-GKS master equations [9,10] in the 60s-70s, the Caldeira-Leggett master equation [3] in the 80s to nonMarkovian processes a la the Hu-Paz-Zhang master equation [6] in the 90s. These foundational work have invigorated new disciplines like open quantum systems [11,12] and enriched new fields such as quantum thermodynamics [13,14]. They also serve as an effective platform to construct theories of nonlinear open quantum systems, the goal of this series of papers.
There are many lines of inquiry into the fundamental issues of quantum and statistical physics of nonlinear quantum open systems. We mention two examples here: 1) Quantum decoherence [15] 2) Quantum chaos [16]. Using the decoherent history formalism [17], following what was done earlier [18] for quasi-classical domains with linear Brownian motion models, Brun [19] derived the quasiclassical equations of motion for nonlinear Brownian systems. The case of linear interactions was treated exactly, and nonlinear interactions are also compared, using classical and quantum perturbation theory. On the issue of quantum chaos and quantum-classical transition, the systematic work of Habib's group is noteworthy [20]. For the most recent developments of quantum chaos involving OTOC, see e.g., the work of Ueda's group [21] (where earlier references of AdS-CFT relevance can be found).
This work. In this work (Paper I) we treat the nonequilibrium dynamics of a quantum Brownian oscillator in a nonlinear potential. Shy of tackling the full nonlinearity we appeal to perturbation theory for the study of a weakly anharmonic oscillator linearly coupled to a thermal field bath. In an earlier paper (Paper 0) [22] we apply perturbation theory to the quantum Langevin equation. We examine the late time behavior of this open system and calculate the power delivered to the system to determine whether it approaches an equilibrium state. From there we extract a fluctuation-dissipation relation. Here we use the functional perturbative method developed in [7]. With a driven quantum oscillator linearly coupled to an environment as the background solution, taking the functional derivatives of physical observables with respect to the source brings down the expectation values of these variables. In [7] this method was applied to the cases where the coupling of the system with the environment is (weakly) nonlinear. Here we treat the case where the system itself is (weakly) nonlinear. A companion paper (Paper II, also based on results obtained in August 2014) studies two weakly nonlinearly coupled quantum oscillators each with its own bath. The correlation functions obtained are useful for the study of heat transfer or charge transport in these systems.

Canonical versus Functional Approaches
To see the power and convenience of the functional approach (despite its seemingly intimidating appearance) we briefly comment on the ambiguities in the operator ordering issue in the canonical approach, then elucidate the functional approach.
Operator ordering issue. The ambiguity in operator ordering often occurs in the context of higher moments of the canonical variables of a quantum system, linear or nonlinear. In the case of multiplyinteracting systems, the subsystems may start off as independent, commutating components, but as the interaction sets in, the operator of one subsystem becomes dependent on the operators of the other subsystems. This introduces ambiguity in the arrangement of the operators associated with different subsystems because they are no longer commutating. Calculations based on the canonical operator approach for these configurations may be undermined by this ambiguity. In contrast, the path integral formulation 1 seems to give a more definite answer; it chooses the time ordering for the in-out formalism 1 Here we only deal with systems that are at most quadratic in the canonical momentum and are not constrained, and with the case that the Hamiltonian/Lagrangian operator itself has a definite operator ordering in terms of the canonical and the path ordering for the in-in or closed-time-path (CTP) formalism [23]. The latter can in addition introduce various operator orderings in sync with the path-ordering. In this sense, the path integral formulation seems to provide in some cases an easier way to carry out quantum calculations 2 .
Functional perturbative approach. The correction to the expectation values of quantum observables due to the nonlinear potential in principle can be found order by order by taking the functional derivatives of the corresponding quantities of the linear system supplemented with external sources. We can obtain the expectation values of the latter system once we have the density operator 3 of such a system and then compute the trace of the product of the observable operator and the density operator. However, this procedure, which involves performing integrations, can be quite tedious in evaluating the expectation values of the various elements of the covariance matrix. A similar but more efficient approach is to calculate only the generating functional of the linear system, driven by external sources. The expectation values of the covariance matrix elements or the two-point functions of the canonical variables can then be found by taking the functional derivatives. We will focus on this approach but still offer the explicit calculations based on the density matrix in the Appendices for comparison.

Driven Oscillators: Nonequilibrium Evolution
To illustrate how the functional techniques are applied to nonlinear open quantum systems it is convenient to first deal with the case of a linear oscillator, driven by an external source j, bilinearly coupled to a quantum scalar field. Taking the functional derivatives with respect to this external source one can derive quantities associated with the linear oscillator, and in perturbative orders, the nonlinear oscillator. Thus we begin with a description of the nonequilibrium evolution of a linear oscillator driven by an external source j.
The action for a harmonic oscillator of mass m and natural frequency ω coupled to a massless Klein-Gordon field φ in 1 + 3 dimensional Minkowski spacetime takes on the form where in the dipole approximation the evolution of the oscillator's internal degree of freedom (idf) χ is decoupled from its external trajectory z(s). Because the whole system is Gaussian we can solve the problem exactly with no need to assume small coupling strength e between the oscillator and the field.
We allow the detector to move along a prescribed trajectory z(s), and thus this action includes that which describes the Unruh effect [24] whence the linear oscillator is referred to as an Unruh-DeWitt detector. Note the difference between the oscillator's idf χ and external degree of freedom (edf) z(s).
For our purpose here, χ is also driven by an external source j.
We now derive the reduced density matrix or the in-in generating functional for this driven linear open quantum system. Perturbative treatment of driven nonlinear open quantum systems will be treated in the next section.

Driven Harmonic Oscillator
The reduced density operator of a free harmonic oscillator linearly coupled with a scalar field bath has been derived before [4,25,26]. The corresponding reduced density operator for the driven harmonic oscillator can be formally obtained in the same manner, and ρ χ (χ a , χ ′ a , 0) is the density matrix of the oscillator at the initial time t a . The two kernel functions G H,β (s, s ′ ) are respectively the retarded Green's function and the Hadamard function of the scalar field evaluated at spacetime points (s, z) and (s ′ , z). They are defined by if the initial state of the scalar field is a thermal state ρ β (t a ) at the temperature β −1 . The function θ(τ ) is the unit-step function.

Coarse-Grained Effective Action
Now introduce the center-of-mass coordinate r and the relative coordinate q the reduced density operator becomes with the coarse-grained effective action S CG [27] given by This effective action describes all the influences from the environmental field in the form of convolution integrals in terms of the kernel functions.
Since the integrand of the path integrals in (2.5) is Gaussian, the path integration will give a result proportional the exponential of the classical value of the coarse-grained effective action where S CG , evaluated in Appendix A, is given by The functional forms of r(s) and q(s) can be found in (A.14) and (A.15).
If we take the trace of the final reduced density operator, it gives an expression which is the in-in version of the generating functional (different from the more familiar in-out version): The normalization constant N determined by Z[0, 0] = 1 takes on the value . (2.10) The function D 2 (t), given in Appendix A shown later, is related to the retarded Green's function of the linear damped oscillator.

'In-In' Generating Functional Z for Driven Linear Systems
It proves useful to have an explicit form of the generating functional Z[j q , j r ] when the initial state of the oscillator has a Gaussian form. Starting from (2.9), we have If, for example, the initial state of the oscillator is a wavepacket of width σ ρ χ (q a , r a , 0) = 1 πσ 2 then carrying out the Gaussian integrals for the generating functional Z[j q , j r ] gives Details of these derivations are given in Appendix B. Note that when j q = 0 = j r , we indeed have If we set e = 0, the generating functional looks exactly the same as an influence functional of some system when the environment is modeled by oscillators and linearly coupled with the system variable j.
In this case, the kernel D 2 (s − s ′ )/m plays the role of the retarded Green's function of the environment, corresponds to the Hadamard function of the environment.

Quantum Statistical Averages of Operators in Linear Systems
Since we already have the explicit expression for the reduced density matrix elements ρ χ (q b , r b ; t) of the driven oscillator at any time t, we may straightforwardly find the quantum (expectation value) statistical average Ô j of the operator associated with a physical variableÔ at that time. For example, the quantum statistical average of χ b = χ(t) is formally given by (2.14) We keep a subscript j to remind ourselves of the dependence on the external source. This allows us to find the quantum statistical averages of physical variables associated with the nonlinear oscillator by performing functional differentiations.
In Appendix C, we compute the quantum statistical averages of the elements of the covariance matrix, χ 2 b j and χ b p b j , for a driven linear oscillator by evaluating the trace such as (2.14). In the limit j → 0, we recover the quantum expectation values of the covariance matrix elements for the free linear oscillator coupled to the quantum scalar field. In fact, the former can also be found by taking the functional derivative of the generating functional (2.12) with respect to various combinations of j r and j q . We will take this approach, and derive in next section the two-point functions of a free harmonic oscillator coupled to a scalar field.
We have argued that the expectation values of the elements of the system covariance matrix can be found with the help of the reduced density operator. This can be quite tedious and repetitive if the reduced density operator is very complicated, and can pose a real challenge for the perturbative calculations. Fortunately these tasks can be accomplished in a more economical way by taking the functional derivatives of the generating functional Z (2.12) with respect to the external sources.

Covariance matrix element χ 2 (τ )
Here we compute χ 2 (τ ) , the second moment of the displacement of the free harmonic oscillator coupled to the field, and it involves a few variants of taking the second-order functional derivatives with respect to j ± .
It is useful to first introduce a couple of shorthand notations In the case of δ 2 Z/δj 2 We find that in this case only does not vanish in the limit j → 0. We may carry out similar functional differentiations of Z[j; t) with respect to j + or j − up to the second order .

(2.23)
They all give the same result in the limit j → 0 because the time-ordering is irrelevant, and it is identical to the result (C.8) obtained by the reduce density matrix according to (2.14), evaluated in Appendix C.
Thus we see it seems much easier to obtain the covariance matrix elements like χ 2 b by taking the functional derivatives of the generating functional. We can also calculate, say, p 2 b in a similar way.
Before showing that, we first discuss how to derive the various real-time two-point functions of χ.

Real-time Two-point Functions of the system variable χ
According to the 'in-in' formalism [1] in terms of the closed-time path integrals [23], we may introduce path-ordered two-point functions by where C +/− represents the forward/backward time branch and T , T + denote time-ordering and antitime-ordering.
Let us first find the Feynman propagator by evaluating the second derivatives of the generating functional with respect to j + at two different times, that is, with 0 < τ, τ ′ < t, where · · · represents terms that will vanish in the j → 0 limit. Compare (2.26) with the decomposition of the Feynman propagator of the free Klein Gordon field φ, we can make these identifications It is consistent with the fact that D 2 (τ ) is a form of retarded function according to its definition in terms of the Laplace transformation.
We may likewise compute the Pauli-Jordan function where · · · are terms that will vanish in the limit j → 0. Since for the free Klein-Gordon field the Pauli-Jordan function can be decomposed into the Hadamard function, the retarded and the advanced Green's functions by, where the overdot denotes taking the derivative of a function with respective to its argument. The coincident limit of (2.33) gives This is the well-known result for the momentum dispersion of the free harmonic oscillator coupled to a scalar field.

Equal-time commutation relation
As a consistency check, we examine whether the uncertainty relation is satisfied at all times in the sense that We first find the expectation values of the commutator of χ at two different times. From the Pauli-Jordan function of χ, we have The lefthand side of (2.35) can then be written as Thus it shows that the equal-time commutation relation is satisfied at all times when the oscillator interacts with the surrounding scalar field and undergoes nonequilibrium time evolution.
In this section we show that when the linear oscillator is coupled to a scalar quantum field, by taking suitable functional derivatives of the 'in-in' generating functional of the oscillator driven by the external sources, we can find the time evolution of the covariance matrix elements, and various real-time two-point functions associated with the free linear oscillator. It can be less painstaking and more efficient than the corresponding calculations based on the reduced density matrix of the free linear oscillator when it couples with the field. These results can be further extended to the case when the nonlinear oscillator.

Nonequilibrium Anharmonic Oscillator
In the previous section we have calculated the expectation values associated with a free harmonic oscillator in contact with a thermal bath, modeled by a quantum scalar field initially in a thermal state.
Here we wish to treat the nonequilibrium evolution of a nonlinear oscillator based on the perturbative functional approach. We also would like to calculate the corresponding expectation values for the anharmonic oscillator. Suppose the action is given by Here the nonlinear potential is represented by V (χ) and contains a small parameter. Since the final where the density matrix element ρ χ (q b , r b , t) of the free oscillator has been given by (2.7) and Thus the expectation value of an operator O consisting only the system variables is given by to ensure proper normalization. Here, and p r , p q are momenta conjugated to r and q, respectively. The subscript 0 on an expectation value O 0 refers to the results without the anharmonic potential.
For example, consider the anharmonic potential where the self-coupling strength λ is assumed to be small. This configuration inherently can induce dynamical instability; however, interaction with the scalar field environment renders it a damped nonlinear oscillator driven by the quantum fluctuations of the field. If the nonlinear potential is weak, in the sense that after a proper rescaling by ω we have λ ≪ γ < ω, and the displacement χ is sufficiently small, then the oscillator's dynamics can remain stable within a certain range where the higher-order contributions may remain small even in the long time limit.
The functional derivatives associated with the nonlinear potential is explicitly given by (3.9) Thus the leading order nonlinear contribution to the expectation value from the nonlinear potential is (3.10) Next we apply this to find the corrections to the in-in generating functional.

Generating Functional Z V for Driven Nonlinear Systems
The full generating functional Z V in the presence of nonlinear action can be expanded by where Z 1 is the leading order correction of Z due to the nonlinear potential V . Thus Z 1 is given by (3.12) First, we show from (2.15) that and they all vanish in the limit j = 0, so we find Z 1 = 0 and then that is, there is no first-order correction of Z V from the nonlinear potential.

Quantum Statistical Averages of Operators in Nonlinear Systems
The quantum statistical average of operator O consisting of the system variables, as noted earlier, is given by In contrast to the generating functional, these expectation values may still have non-zero first order nonlinear contributions. In particular, we are interested in the expectation values of the covariance matrix elements because they can be used to construct quantities we are interested in for the Gaussian system, or the nonGaussian system in terms of the perturbative treatment from the Gaussian system.
Here we will use as an example the first-order correction of χ b due to the nonlinear potential.
It is straightforward to show that for 0 < τ < t. Therefore we have a nonvanishing first-order nonlinear contribution to χ b due to the anharmonic potential (3.8), At first sight this may seem unusual because from perturbative quantum field theory there should be no first-order nonlinear contribution to the expectation value of the position operator of the anharmonic oscillator. However, the presence of the environment makes a difference.
It may be easier to see the difference from an alternative approach. The Langevin equation of this anharmonic system is where the damping constant γ = e 2 /8πm is associated with the dissipative dynamics of the system from the backreaction of the environment, while the stochastic noise is Gaussian and its second moment is related to the expectation value of the anti-commutator of the environmental field 20) with ξ(t) = 0 in the dipole approximation and x = (t, x). The double brackets · · · represent taking the ensemble average of some quantity with respect to the probability distribution of the noise.
Let us solve the Langevin equation perturbatively by assuming the solutions take the form x(t) = x (0) (t) + λ x (1) (t) + · · · . Then we have The solution to (3.21) is given by Plugging it in the righthand side of (3.22) we get Thus the corresponding mean value is given by assuming there is no correlations between the conjugated variables of the system, nor cross correlations between the system and the environment at the initial time. This is exactly what we have obtained in Similarly we can show that p b also have nonvanishing first-order correction from the potential The expressions and their derivations are collected in Appendix D. They are useful for exploring many interesting physical problems, such as the approach to equilibrium and the existence of a fluctuationdissipation relation at late times in weakly anharmomic oscillator open quantum systems [22], or energy transport and the existence of nonequilibrium steady state in weakly nonlinear quantum systems with multiple baths, using the results of the linear cases explored before (e.g., [6,25]). These problems are investigated in the present series of papers.

Appendix A. Derivation of the Classical Values of the Coarse-Grained Effective Action S CG
To carry out the path integral in (2.5), we will need the classical trajectories of r(s) and q(s), which can be obtained by taking the variations of the coarse-grained effective action with respect to q(s) and r(s), They will be used to find the classical value S CG of the coarse-grained effective action (2.6), where r(s), q(s) satisfy the equations of motion. From (A.1), we see that in general r(s) is a complex function with real and imaginary parts, r(s) = r R (s) + i r I (s), each of which satisfies Since the equations of motion for q(s) and r R (s) depend on the sources j q and j r , their general solutions will also have contributions from the sources, Hereafter, we will suppress the subscript R in r R (s) unless confusion arises.
The fundamental solutions D i (s) are the homogeneous solutions to (A.4), satisfying As for q(s), it can be related to r(t − s), but since it has a source term in the equation of motion (A.2), this correspondence is less transparent. Let us assume that q(s) is given by Then we want to show that it satisfies (A.2). Substituting (A.10) into the lefthand side of (A.2) we have with τ = t − s and τ ′ = t − s ′ . We note that t > τ > τ ′ > 0 by construction, and the fundamental solution D 2 (s) is zero when s < 0, so the last term in (A.11) can be re-written to Therefore we have (A.11) given by It can be shown that the second term in (A.12) vanishes, and (A.2) follows. Therefore (A.10) is indeed the general solution to (A.2). For future reference we can also write (A.10) as Next we would like to express (A.9) and (A.13) in terms of the endpoints q b , r b and q a , r a , so that it will facilitate calculations of the evolutionary operators used to evaluate the density matrix (2.2).
After some algebraic manipulation, we have r(s) = α(s) r a + β(s) r b + J r (s) , (A.14) The evolutionary propagator J(r b , q b , t; r a , q a , 0) is then given by with r(s) and q(s) given by (A.14) and (A.15).

Appendix B. Evaluation of the 'In-in' Generating Functional
It is instructive to give an explicit calculation of the generating functional Z[j q , j r ] when the initial state of the oscillator assumes a Gaussian form. Starting from (2.9), we have where we have introduced the following shorthand notationṡ to make the generating functional Z[j q , j r ] more manageable. We also assume that the initial density operator of the system is given by In fact, it can be any Gaussian form more general than (B.2).
We first perform the integral over q b , and then over r b , so (B.1) becomes Carrying out the integral of q a is straightforward and leads to Before proceeding, we first re-write some of the expressions in the exponential, where we have defined due to the fact that D 2 (s ′ − s) = 0 for s ′ < s. Combining these results together gives The normalization is chosen such that when j q = 0 = j r , we have Z[0, 0]=1.

Appendix C. Quantum Statistical Averages in the Linear Model in Sec. 2.2
Here we present detailed and explicit calculations of the quantum statistical averages of the covariance matrix elements for the driven, linear oscillator model, based on the reduced density matrix of χ.
Evaluating multiple Gaussian integrations is straightforward, even a bit tedious, but they are handy for future references and extensions. Among these calculations, the evaluation of χ b p b j is particularly important because it offers a self-consistent check on the equal-time commutation relations.
Here we compute χ b j with j collectively representing {j r , j q }. Thus χ b j is the quantum expectation value of χ when the external source j is present. The procedure is very similar to the calculations of the generating functional. We start with where we have usedq a to denoteq The relevant notations and the initial state of the oscillator have been introduced in Sec. 2.1.2. Thus the expectation value of χ b becomes where we have used Therefore we arrive at In the case that the external sources j q , j r vanish, we recover χ b = 0 for the case of a free harmonic oscillator bilinearly coupled to a massless scalar field initially in a thermal state.
Here we calculate χ 2 b j , From the results in the previous section, we have χ 2 b given by In the absence of external sources, this reverts to the well-known result so we may conclude p b j = m χ b j . This is odd because one expects the equal sign should hold when the oscillator is driven by an external classical source. The resolution lies in the observation that if j describes a physical, external source, then at the end of the calculation, we need to set j q to zero, that is, j + = j − = j exactly like what we do for χ ± . Doing so we obtain the expected result that (1) given by We break the whole expression into smaller components.
Appendix D.1.1. (δ/δj r ) 3 : We then have because all terms that are proportional to J q or Ξ will definitely vanish in the limit j q = 0 = j r . This result can also be verified with the help of the Langevin equation (3.19). Since χ(t) = χ (0) (t) + λ χ (1) (t) + · · · , the second moment of χ(t) with the leading order nonlinear contribution from the anharmonic potential is given by Now for χ (0) (t) χ (1) (t) , we have For the initial density matrix of the Gaussian form we have chosen, it implies that there is no correlation between χ(0) andχ(0), and that there is no odd moments for the initial position or the initial momentum. In addition, since the noise ξ is Gaussian by construction, its three-point function vanishes too.
Therefore we conclude there is no nonlinear contribution from the anharmonic potential to the second moment of χ(t) to order O(λ). This is consistent with (D.3).

Appendix D.2. p b
For p b , the leading order nonlinear contribution is given by  We now evaluate the leading order nonlinear contribution term by term.
Since we have δ δj q (τ ) in the limit j → 0, we arrive at j=0 , (D.15) and thus Thus again we obtain (C.15) by the functional method.