Exact Master Equation for Quantum Brownian Motion with Generalization to Momentum-Dependent System-Environment Couplings

In this paper, we generalize the quantum Brownian motion to include momentum-dependent system-environment couplings. The conventional QBM model corresponds to the spacial case $W_k = V_k$. The generalized QBM is more complicated but the generalization is necessary. This is because the particle transition and the pair production between the system and the environment represent two very different physical processes, and usually cannot have the same coupling strengths. Thus, the conventional QBM model, which is well-defined at classical level, is hardly realized in real quantum physical world. We discuss the physical realizations of the generalized QBM in different physical systems, and derive its exact master equation for both the initial decoupled states and initial correlated states. The Hu-Paz-Zhang master equation of the conventional QBM model is reproduced as a special case. We find that the renormalized Brownian particle Hamiltonian after traced out all the environmental states induced naturally a momentum-dependent potential, which also shows the necessity of including the momentum-dependent coupling in the QBM Hamiltonian. In the Hu-Paz-Zhang master equation, such a renormalized potential is misplaced so that the correct renormalization Hamiltonian has not been found. With the exact master equation for both the initial decoupled and and initial correlated states, the issues about the initial jolt which is a long-stand problem in the Hu-Paz-Zhang master equation is also re-examined. We find that the so-called"initial jolt", which has been thought to be an artificial effect due to the use of the initial decoupled system-environment states, has nothing do to with the initial decoupled state. The new exact master equation for the generalized QBM also has the potential applications to photonics quantum computing.

In this paper, we generalize the quantum Brownian motion to include momentum-dependent system-environment couplings. The resulting Hamiltonian is given by . The conventional QBM model corresponds to the spacial case W k = V k . The generalized QBM is more complicated but the generalization is necessary. This is because the particle transition and the pair production between the system and the environment represent two very different physical processes, and usually cannot have the same coupling strengths. Thus, the conventional QBM model, which is well-defined at classical level, is hardly realized in real quantum physical world. We discuss the physical realizations of the generalized QBM in different physical systems, and derive its exact master equation for both the initial decoupled states and initial correlated states. The Hu-Paz-Zhang master equation of the conventional QBM model is reproduced as a special case. We find that the renormalized Brownian particle Hamiltonian after traced out all the environmental states induced naturally a momentumdependent potential, which also shows the necessity of including the momentum-dependent coupling in the QBM Hamiltonian. In the Hu-Paz-Zhang master equation, such a renormalized potential is misplaced so that the correct renormalization Hamiltonian has not been found. With the exact master equation for both the initial decoupled and and initial correlated states, the issues about the initial jolt which is a long-stand problem in the Hu-Paz-Zhang master equation is also re-examined. We find that the so-called "initial jolt", which has been thought to be an artificial effect due to the use of the initial decoupled system-environment states, has nothing do to with the initial decoupled state. The new exact master equation for the generalized QBM also has the potential applications to photonics quantum computing.

I. INTRODUCTION
As it is well-known, any realistic system inevitably interacts with its environment. For nano-scale quantum devices or more general mesoscopic and quantum systems, such interactions are usually not negligible, and thus these systems must be treated as open quantum systems. The essential physics of open quantum systems is the entanglement generation between the system and the environment so that the system cannot be retained in pure quantum states. This makes the system lose its quantum coherence, known as quantum decoherence which is the key obstacle in the development of quantum technology. Up to date, theoretical description of open quantum systems is still a big challenge in both fundamental research and practical applications. The major difficulty for solving the dynamics of open quantum systems is how to find the physically reasonable systemenvironment coupling which is usually unknown in priori, and how to derive unambiguously the equation of motion, called the master equation, for the dynamical evolution of open quantum systems in terms of the reduced density matrix which is not always computable.
Historically, a prototype of open quantum systems is the quantum Brownian motion (QBM) [1,2]. Quantum mechanically, the Brownian particle is modeled as a harmonic oscillator that linearly couples to the environment which is made of a continuous distribution of infinite * Electronic address: wzhang@mail.ncku.edu.tw number of harmonic oscillators, as originally proposed in the seminal works of Feynman and Vernon [3] and also Caldeira and Laggett [4,5]. The system-environment coupling in QBM was assumed as a strict linear coupling k C k xq k , where x and q k are the position coordinates of the system and environment oscillators, respectively, and C k is the coupling amplitude. The exact master equation of the QBM for such a system-environment coupling has been derived during 1980's ∼ 1990's with different approaches [4,[6][7][8]. Correspondingly, the dissipative dynamics of QBM has been extensively investigated in a rather broad physical research areas [1,2]. In particular, the classical dissipative dynamics, the fluctuationdissipation theorem and the classical Langenvin equation have been deduced from this quantum Brownian motion in the weak-coupling and high temperature limits [1][2][3][4][5].
Microscopically, the above classically well-defined linear coupling k C k xq k depicts the Brownian particle dissipation and fluctuations through energy exchange between the system and the environment. However, probably one does not realize that this system-environment coupling is physically not self-consistent at quantum level. Specifically, the Hamiltonian of a harmonic oscillator can be written as H s = ω s a † a, where ω s is the quanta of the harmonic oscillator energy of the system, and a † and a are the corresponding creation and annihilation operators obeying commutation relation [a, a † ] = 1. Similarly, the Hamiltonian of the environment made by infinite number of harmonic oscillators can be expressed as H e = k ω k b † k b k , where b † k , b k are the creation and annihilation operators of the quantum oscillator mode k with frequency ω k . Consequently, the linear system-environment coupling Hamiltonian can be expressed as where M , m k are the masses of the system oscillator and the environmental oscillator mode k. It shows that the first two terms, proportional to a † b k + b † k a, depict the energy exchange (energy transition) between the system and the environment. While the last two terms, proportional to a † b † k + b k a, describe the processes of generating and annihilating separately two quanta of energy out of the blue, from nothing. As one knows, the energy quanta of harmonic oscillator can be realized with photons or phonons, and physically one can generate or annihilate two photons (phonons) through non-linear crystals or atomic virtual states, but these physical processes are very different from the single photon transition processes of the first two terms in linear quantum optics. In other words, in reality, the coupling strength V k for the single particle transition processes is usually not the same as that of the two particle generation or annihilation.
Analogy to the above linear coupling between the system and the environment for the QBM, a similar situation also occurs in the canonical quantization of lightmatter interaction. In quantum optics, under the dipole momentum approximation, the interesting light-matter interaction is mainly determined by −d · E(x, t), where d = qr is the atomic dipole moment and E(x, t) is the radiative electric field. Treating approximately the atom as a two-level system, the above light-matter interaction can be written quantum mechanically as −d · E(x, t) = where σ x = σ − + σ + is the Pauli matrix, and b † k , b k are the creation and annihilation operators of photons for the radiative field mode k. In this light-matter interaction, the first two terms proportional to σ + b k +σ − b † k depict the electron transitions in atoms from the low-energy state to the high-energy state by absorbing a photon, and from the high-energy state back to the low-energy state by emitting a photon. These are the fundamental physical processes for light-matter interaction. While, there are other two terms, proportional to σ + b † k + σ − b k , which describe the electron transitions from the low-energy state to the high-energy state by emitting a photon, and from the high-energy state to the low-energy state by absorbing a photon. Actually, the last two processes barely occur due to the law of energy conservation, see Fig. 1, also see more discussions given in the next section.
In the literature, these physically infeasible processes in the linear coupling of the QBM and in the lightmatter interaction of quantum optics are often called the counter-rotating terms. They can usually be eliminated under the so-called rotating wave approximation. In the high-temperature and weak-coupling limit, it can be shown that the contributions from these counter-rotating terms are negligible. However, in the strong coupling and low temperature regime, the counter-rotating terms could make a significant effect. In recent years, many effects have been devoted to search the counter-rotating effects theoretically and also experimentally [9][10][11][12][13][14][15]. Here we should point out that these counter-rotating terms either are forbidden by the energy conservation law for every individual microscopic physical processes, or can be modified by the comprehensive physical description of the system-environment couplings. Specifically, for an example, for the radiative electromagnetic field, E = −∂A/∂t and B = ∇ × A where the vector potential A and the electric field E are the generalized coordinate and momentum in quantum electrodynamics (QED). Besides the atomic dipole moment interaction with the electric field, there is also the atomic magnetic dipole moment interaction with the magnetic field. Both interactions come from the position-momentum dependent atom-light interaction − q m A · P e , where P e is the electron momentum [16]. For the QBM, besides the linear coupling between the oscillator coordinates of the system and the environment, there are also the momentum-dependent system-environment couplings. When these contributions are elaborated, the inconsistent processes in systemenvironment couplings of open quantum systems can be diminished.
In this paper, we shall provide first a detailed analysis of the system-environment couplings for open quantum systems in general, and then generalize the QBM by including the momentum-dependent coupling to get rid of the possible ambiguity of system-environment couplings at quantum level. We shall also discuss practical realizations of such a generalization. In Sec. III, we will derive the exact master equation of the generalized QBM for both the initial decoupled and initial correlated states. We shall also explore the relation of dissipation and fluctuation dynamics in open quantum systems with the decay and diffusions of the Brownian particle. In Sec. IV, we will show how the previous Hu-Paz-Zhang master equation for the QBM without including momentum-dependent coupling [4,[6][7][8] can be reproduced as a very special case. We will also re-examine some long-standing debated issues about QBM, such as the origin of the initial jolt phenomena and the initial decoupled state problem one often questioned for the study of the QBM with the exact master equation. Finally, discussions and perspectives are given in Sec. V. In Appendix A, the nonequilibrium Green functions for open quantum systems we previously introduced are generalized to open systems containing pairing couplings through the Heisenberg equation of motion. The generalized nonequilibrium fluctuation-dissipation theorem in the time-domain is also derived there.

II. A GENERAL ANALYSIS OF MOMENTUM-DEPENDENT SYSTEM-ENVIRONMENT COUPLINGS
In the early studies of quantum dissipative dynamics, system-environment couplings in open systems are mainly modeled based on classical mechanics. The prototype example is given in the pioneering work by Feynman and Vernon for developing their influence functional theory for QBM [3]. In the exploration of the dissipative and noise effects of other interaction systems on the system of interest, Feynman and Vernon proposed a general test system x interacting to linear systems {q k } (an oscillator bath) with the Lagrangian Without loss of generality, one usually takes a linear coupling U (x) = x. The oscillator bath is an environment consisting of a collection of harmonic oscillators with a continuous frequency distribution. Such an environment can be easily realized as radiative electromagnetic fields or vibrations of lattices in matter. The Lagrangian formulation in terms of the dynamical variables of positions and velocities is setup in order to conventionally use the Feynman's path integrals for the trace over all the environmental states [17]. In such circumstances, the path integrals over environmental degrees of freedom can be done exactly. Thus the environment effect on the test system is given by an influence functional of coordinates of the system only. It is equivalent to external semiclassical dissipative and fluctuating (random) forces acting on the test system. As Feynman and Vernon pointed out [3], if a general non-linear interaction system q k is weakly coupled to a test system x with a more general interaction potential U (x)U ({q k }) which is small, then the effect on the test system is that of a sum of oscillators whose frequencies correspond to the possible transitions of the interaction system. Thus, to the extent that second order perturbation theory yields sufficient accuracy, the effect of an interaction system is that of a linear system [3]. In 1980's, Caldeira and Leggett made further investigations to QBM and quantum tunneling of dissipative systems, in which similar linear couplings are made in order to derive a dissipative (damping) equations of motion from quantum mechanics [4,5]. Specifically, to consistently produce the semiclassical damped equation of motion, Mẍ + η(x)ẋ + dV /dx = F ext (t) where F ext (t) is an external force, Caldeira and Leggett also proposed the environment as a continuous frequency distribution of a set of harmonic oscillators interacting linearly with the system, When the system is a tunneling system, Eq. (2) is often called the Caldeira-Leggett model in the literature. Except for the last term which is added as a counterterm for the stability of potential renormalization, Eq. (2) is equivalent to Eq. (1). For x-independent dissipation η(x) = η, the semiclassical damped equation of motion can be realized with the strict linear coupling between the system and the environment F k (x) = C k x. Other more complicated coupling mechanisms are analyzed extensively in the Appendix C of [5] but these analysis are indeed superfluous from the microscopical level of interacting systems, as Caldeira and Leggett pointed out in their own work [5].
To have a clear physical picture on the microscopical level of interacting systems in the Feynman-Vernon model or the Caldeira-Leggett model, it will be more convenient to reformulate the Hamiltonian of Eq. (2) in the second quantization framework: Here, ε i is the spectrum of the system Hamiltonian H s = The operators a † i and a i are respectively the creation and annihilation operators that create the system state |ψ i from the vacuum state |0 : a † i |0 = |ψ i and a i |0 = 0. These creation and annihilation operators obey the commutation or anti-commutation relation, [a i , a † j ] ∓ = a i a † j ∓ a † j a i = δ ij , depending on the system being a bosonic or a fermionic system. The creation and annihilation operators b † k , b k of the harmonic bath can also be defined by The system-environment coupling V ijk = 2m k ω k ψ i |F k (x)|ψ j describes the transition of the system from the state |ψ j to the state |ψ i due to the linear coupling to the environmental oscillating mode ω k . Now one can see a physically troublesome process depicted by the linear system-environment coupling in the Feynman-Vernon model and the Caldeira-Leggett model. That is, the system emits or absorbs a "photon" with the same energy results in the same transition of the system from the state |ψ j to the state |ψ i . This is quantum mechanically infeasible, as shown in Fig. 1. This seems to indicate that Eqs. (1) and (2) are at least incomplete at quantum level if it is not unphyscial.
In fact, a general Hamiltonian of the system coupled to the oscillator-bath environment should contain both the position-dependent and the momentum-dependent system-environment couplings, as Caldeira and Leggett pointed out in their original paper [5], where Φ(p, x) is a counter-term for the renormalization of system Hamiltonian arisen from the systemenvironment couplings. Equation (4) contains various position and momentum dependent system-environment couplings. The second quantization form of the above A microscopic picture of the momentumindependent couplings between the system and the environment in Eq. (3). (a) The process of a particle decay from the high energy state |ψj to the low energy state |ψi by emitting an energy quanta described by the term V ijk a † i ajb † k , and its Hermitian conjugation process; (b) The process of the particle decay from the high energy state |ψj to the low energy state |ψi by absorbing an energy quanta described by the term V ijk a † i ajb k , and its Hermitian conjugation process, which are physically infeasible.

Hamiltonian is given by
Similarly, ε i is the spectrum of the system Hamiltonian Now it shows that the system-environment coupling (the last two terms in Eq. (5)) contains the physical process of the system transiting from the state |ψ j to the state |ψ i by emitting a photon, and the inverse process of transiting the state |ψ i back to the state |ψ j by absorbing a photon with the same energy. In other words, the original ambiguity in the Feynman-Vernon model and in the Caldeira-Leggett model can be diminished when a momentum-dependent system-reservoir coupling co-exists. The above analysis indicates that although Eqs. (1) and (2) are well-defined in classical mechanics, their quantum mechanical version is somehow inconsistent and contains some physically infeasible processes. But the inconsistency can be eliminated when a momentumdependent system-environment coupling is included, as it is shown above. In fact, for atoms interacting with radiative electromagnetic field, the momentum-dependent coupling is naturally manifested. Explicitly, the nonrelativistic atom-photon interaction Hamiltonian can be rigorously written as [16] where m and e are the electron mass and charge, P i = −i ∇ i is the i-th electron momentum operator in the atom. The vector field A(r) is the canonical coordinate of the electromagnetic field, and is quantized by with the Coulomb gauge k·ê λ = 0. The last term V coul = V e-N coul + V e-e coul in Eq. (7) contains all the electron-nuclei and electron-electron Coulomb interactions in the atom arisen from the charge-induced scalar potentials. The atom-light interaction in Eq. (7) is given by the term − e m i A(r i )·P i , when the Coulomb gauge ∇·A(r) = 0 is used. This is a typical momentum-dependent coupling between the system and the environment if the atom is considered as the system and the radiative field is treated as the environment.
Taking the second quantization to the electron energy eigenfunctions ψ α (r) of the atom under the interaction with nuclei, we have 1 2m P 2 +V e-N coul (r) ψ α (r) = E α ψ α (r), where the nuclei is considered to be static. Then the atom-light interaction Hamiltonian of Eq. (7) can be rewritten as Here we ignore the term e 2 2m i A 2 (r i ) because it only gives shift to the radiative frequencies for all the radiative field modes. The operators a † α creates the single particle state ψ α (r) from the vacuum state. The transition ampli- e −ik·r ψ * α (r)ê λ · P ψ β (r)d 3 r describes the photon-emitting transition from the atomic state ψ β to the atomic state ψ α under the constraints of energy, angular momentum and parity conservations (the transition selection rules). The transition amplitude V * αβkλ describes the photon-absorbing transition from the atomic state ψ α to the atomic state ψ β . In quantum optics, one usually takes the dipole moment approximation (e ik·r 1) and treats the atom as a two-level system, then the momentum-dependent atom-light coupling − e m i A(r i )·P i is approximated to an electric dipole in an electric field −d · E which would induce the counterrotating terms as we discussed in the introduction. These counter-rotating terms arisen from the dipole approximation is indeed artificial. A rigorous proof of no existence of counter-rotating terms can be found from the fundamental quantum theory of light-matter interaction, i.e., the quantum electrodynamics within the framework of quantum field theory [18].
On the other hand, even for the Feynman-Vernon model Eq. (1) and the Caldeira-Leggett model of Eq. (2), one cannot solve exactly the corresponding nonequilibrium decoherence dynamics in general. Only when one takes the system also as a simple harmonic oscillator, H s = 1 2M p 2 + 1 2 M Ω 2 x 2 , and focuses on the strict linear coupling between the system and the environment, i.e. U (x) = x in Eq. (1) and F k (x) = C k x in Eq. (2), can an exact master equation for the reduced density matrix of the system be derived [6][7][8]. In such circumstances, the simplified Feynman-Vernon model or Caldeira-Leggett model still contains infeasible physical processes at quantum level arisen from the canonical quantization. Explicitly, let us write the Hamiltonian of Eq. (1) or Eq. (2) for the system being a harmonic oscillator linearly coupled to the oscillator bath in the particle number representation, where ω 2 s = Ω 2 + Ω 2 c and Ω 2 c = k C 2 k /M m k ω 2 k is arisen from the counter-term to cancel the possible divergence of the renormalized frequency-shift [5], and As we have pointed out in the introduction, physically, the terms proportional to a † b † k + b k a correspond to the processes of generating and annihilating two energy quanta separately from nothing. Consider a practical realization of quantum harmonic oscillator by photons, these processes correspond to two-photon creation and annihilation in nonlinear optics. These twophoton processes could be physically realized, for example, by spontaneous parametric down-conversion which is very different from the other two terms proportional to a † b k + b † k a in Eq. (10). The later is a linear particle process corresponding to the quantum energy tunneling (energy quanta exchange) between the system and the environment. Consequently, the system-environment coupling strengths for the linear particle exchange and the nonlinear particle pairing production/annihilation must be different in the reality, and therefore Eq. (10) is physically incomplete at quantum level.
However, if the momentum-dependent couplings are included, the general strict linear coupling Hamiltonian can have the form which can resolve the above problem. In other words, the quantum mechanically feasible QBM Hamiltonian that includes momentum-dependent couplings should be where W k = V k in general. Obviously, the conventional QBM Eq. (10) in the literature is a special case of Eq. (12) with W k = V k . However, Eq. (12) can be easily realized in many different quantum systems. A typical example is the cavity coupled to the environment with an external squeezed driving, Taking a Bogoliubov transformation a = uc + vc † , a † = u * c † + v * c with |u| 2 − |v| 2 = 1, it is easy to show that this Hamiltonian becomes Eq. (12) with This shows explicitly how the system-environment coupling strengths for one-photon processes and two-photon processes become different in the reality. Thus, the generalization of the QBM with momentum-dependent coupling between the system and the environment is physically more acceptable.
Very recently, investigations of quantum photonics has made a great progress for solving boson sampling problem [19][20][21][22]. Quantum photonics is made of many squeezers and interferometers (beam splitters and phase shifters), which can be described by the following Hamiltonian including various photon loss, where index α labels different photon loss environments. After diagonalizing the Hamiltonian of the squeezers and interferometers (the first two terms in the above Hamiltonian) through a Bogoliubov transformation, the total Hamiltonian can be reduced to This is an extension of the generalized QBM Eq. (12) to multimodes that is practically used in the integrated quantum photonics for the current investigation of photonics quantum computing [23][24][25][26][27].
Analogically, in the current investigation of topological phases of matter, many quantum devices are made by superconductor-semiconductor, superconductorinsulator and superconductor-metal hybrid systems [28][29][30][31]. They are coupled to leads and thereby become open quantum systems. The corresponding Hamiltonian, after a Bogoliubov transformation, can be equivalently written as where index α labels different leads. This Hamiltonian has the same structure as the one for photonic quantum processor of Eq. (15), except that here it describes electron dynamics rather than photon dynamics in open quantum systems. Physically, Eq. (16) can be considered as a fermionic realization of the generalized QBM. The general exact master equation of Eq. (16) has been derived recently by us for initially correlated systemenvironment states [32]. Applications to the study of decoherence dynamics on exotic states, such as anyons and Majorana quasiparticles for topological quantum computer, have also been explored [33][34][35][36].
The above analysis indicates that there has been in the literature too much emphasis on quantization (i.e. general methods of obtaining quantum mechanics from classical methods) as opposed to the converse problem of the classical limit of quantum mechanics. This is unfortunate because the latter is an important question for various areas of modern physics while the former is a chimera [37]. From a more fundamental point of view, since the system plus the environment is considered as a closed system, the quantum physical process should ensure the compatibility between the unitarity and the conservation of energy-momentum transfer in every single microscopic quantum process. In the following, we shall focus on the more general and quantum mechanically self-consistent model with Eq. (12) as a generalization of the conventional QBM. We will derive its exact master equation for both initially decoupled and initially entangled systemenvironment states. For its applications in quantum photonics, such as the integrated quantum processor made of many squeezers and interferometers with hundred more photonic modes for photonics quantum computing that is described by Eq. (15), we will leave it for further investigation.

III. THE EXACT MASTER EQUATION OF THE GENERALIZED QBM
For an arbitrary initial state ρ tot (t 0 ), the time evolution of the total system (system plus environment) is determined by the Liouville-von Neumann equation of the quantum mechanics [38], Because the system and the environment together form a closed system, the Liouville-von Neumann equation is the same as the Schrödinger equation of quantum mechanics for the evolution of pure quantum states. But the Liouville-von Neumann equation is more general because it is also valid for statistically mixed states, such as Eqs. (22) and (23), where the Schrödinger equation is not applicable. The time evolution of the system is determined by the reduced density matrix which is defined as the partial trace of the total density matrix over all the environment states, Applying the partial trace over the environment states to Eq. (17) with the total Hamiltonian of Eq. (12), we obtain formally the exact master equation where and A[ρ s (t)] is a collective operator defined by Here L ± [ρ s (t)] are positive/negative supercurrent operators that describe particles flowing into/out of the system. It contains the usual current through particle tunnelings between the system and environments, and also the current via particle pair production and annihilation. The latter is similar to the Andreev reflections through the interface between superconductors and normal metals. Particle dissipation and fluctuations arisen from the system-environment back-actions can be naturally manifested after completed the partial trace over all the environment states in Eq. (21). Now the problem left is how to carry out explicitly the partial trace over all the environment states in Eq. (21). For the initial total density matrix ρ tot (t 0 ) being a partition-free thermal state (containing the initial system-environment correlation), because the total Hamiltonian is a bilinear function of particle creation and annihilation operators of the system and environment, the total density matrix and the reduced density matrix always evolve in Gaussian states after a dynamical quench. In this case, it is rather easy to trace over all the environment states with Gaussian integrals in the coherent state representation. Such configuration has been presented in our previous derivation of the exact master equation for dissipative topological systems [32]. Below we will focus on cases with initial decoupled state. In the end of the section, we will make a connection to the master equation with the above initial correlated state. Thus, we now assume that the system and the environment are initially decoupled, and the environment is initially in a thermal state [3,4], Such initial state can be prepared by tuning coupling strength between the system and the environment to zero at initial time. Although the initial total state is a system-environment decoupled state, the system initial state ρ s (t 0 ) can be arbitrary (non-Gaussian states). Therefore, the partial trace over the environment states in Eq. (21) is rather difficult to carry out. To complete this partial trace, we shall use the coherent state path integral method [39]. In the coherent state representation, the matrix element of the collective operator Eq. (21) can be expressed as In Eq. (24), |z ≡ e a † z |0 is the unnormalized coherent state. Correspondingly, the integral measure over the complex space is given by dµ(z * , z) ≡ dz * dz 2πi e −|z| 2 . The A-operator associated propagating function J A (z * t , z t , t; z 0 , z * 0 , t 0 ) is defined in a similar way as the propagating function for the reduced density matrix in the coherent state representation [40], The propagating function J (z * t , z t , t; z 0 , z * 0 , t 0 ) fully describes the time evolution of the reduced density matrix of Eq. (18). Similarly, J A (z * t , z t , t; z 0 , z * 0 , t 0 ) fully determine the evolution of the collective operator A[ρ s (t)] of Eq. (21).
Utilizing the coherent state path integrals, we have and similarly, Here, S s [z * , z] is the action of the system, given by The first two terms in the action are geometric terms arisen from the boundary (the end points of the coherent state path integral). This is because the paths in the coherent state path integrals are only fixed on one side: z(t 0 ) = z 0 , z * (t) = z * t and z * (t 0 ) = z * 0 , z (t) = z t , see Eqs. (24) and (25). While (z(t), z * (t)) and (z (t), z * (t)) are not the complex conjugate pairs of the above boundary conditions, and they must be determined by solving the path integral explicitly [41].
On the other hand, F[z * , z * , z, z ] in Eq. (26) is the influence functional arisen from the partial trace over all the environmental states [3]. It is obtained by rigorously integrating out all the environment degrees of freedom through the closed-loop path integrals [42], where we have introduced the new variables z q ≡ z − z , z q * ≡ z * − z * . The integral kernels G(τ, τ ) and G(τ, τ ) are the two-time system-environment correlation functions, Here we have also introduced the notations . The phase factors come from the particle evolution in the environment before it transits into the system through the tunneling or is created from the pairing coupling in Eq. (12). It shows that only the time correlation function G(τ, τ ) relates to the initial environment state through the initial particle occupations and particle squeezed parameters in the environment, n k (t . For the initial environment thermal state Eq. (23), the particle squeezed parameters s k (t 0 ) become zero but here we formally keep it so that it can also apply to other initial environment states. These time correlations depict back-action processes between the system and the environment. While, the operator A-associated influence functional, after taking the partial trace over the environment states, can be reduced as where Z = [ 1 0 0 −1 ] is the z-component Pauli matrix. Furthermore, the path integrals in the propagating function of Eq. (26) can be exactly carries out because of quadratic form of the effective action (including the influence functional). The result is Here z q (t), z q * (t) and z q (t 0 ), z q * (t 0 ) are determined by the stationary path, d dτ with t 0 ≤ τ ≤ t. This is a generalization of our previous work [42] to the open systems with paring correlations. To solve the above equation, we introduce the transformation with different boundary conditions of the initial and final ends in the coherent state path integrals, subjected to the boundary conditions U(t 0 , t 0 ) = I and V(t 0 , t) = 0. It is interesting to find that U(τ, t 0 ) and V(τ, t) are the generalization of the nonequilibrium Green functions for open quantum systems we previously introduced [40,[42][43][44] to the case involving pairing processes.
In terms of the usual definition of nonequilibrium Green functions, they can be expressed as They fully describe the transient dissipation and fluctuation dynamics of open quantum system incorporating various back-actions from the environment, and obey the integro-differential equations of motion (35) which can also be derived from the Heisenberg equation of motion, see the detailed derivation given in Appendix A. The solution of Eq. (35b) for the nonequilibrium noise-induced correlation Green function V(τ, t) can be obtained easily, (37) which is indeed the generalization of the Keldysh's correlation Green function with pairing couplings [45]. This solution is the generalized nonequilibrium fluctuationdissipation theorem [44], also see the derivation given in Appendix A.
With equation of motion Eq. (33), the A-associated propagating functions become ×J (z * t , z t , t; z 0 , z * 0 , t 0 ), (38) and using Eq. (34), the additional term can be fully determined in terms of the fixed end points of the path integrals. From the above result, we obtain the current superoperators Eq. (20) in the coherent state representation After some rearrangement, one can rewrite the exact master equation as the standard form as we did before where H s (t) = ω s (t)a † a + 1 2 ω s (t)a †2 + 1 2 ω * s (t)a 2 . Other coefficients are given by It shows that the coefficients ω s (t) and ω s (t) are the renormalized frequency and renormalized pairing strength modified by the environment. By diagonalizing the renormalized Hamiltonian with a Bogoliubov transformation, the renormalized eigen-frequency is ω r (t) = |ω s (t)| 2 − |ω s (t)| 2 . The term proportional to γ(t) represents a non-unitary evolution that describes the dissipation dynamics induced by the environment. Notice that the definition of γ(t) here is twice of that in our previous paper [32] because we rewritten the superoperator multiplied by a factor 1/2 so that all the superoperators have formally the same standard Lindblad form [46]. As a result, γ(t) corresponds precisely to the decay coefficient of the system. Also, there is a term proportional to γ (t) in our previous paper [32] which vanishes here for the single Brownian particle system. The other two terms proportional to γ(t) and γ(t) in the master equation describe fluctuation (diffusion) dynamics. They depend on the initial state of the environment. Note that the γ(t)associated fluctuations are induced by pair production and annihilation between the system and the environment.
To understand further the physical picture of the renormalized energy, the dissipation and fluctuation dynamics induced by the environment, let us calculate the mean values and their deviations (equivalent to the quadrature covariance in terms of the position and momentum) of the system variables from the exact master equation. From the exact master equation (41), it is easy to find that d dt where a(t) = Tr s [aρ s (t)] and a † (t) = Tr s [a † ρ s (t)].
The above equation shows that the Brownian particle undergoes a damping oscillation with the renormalized frequencies ω s (t) and ω s (t) and the decay coefficient γ(t) induced from the environment. The mean value a(t) is mixed with its complex conjugation which is related to the squeezed term 1 2 ω s (t)a †2 in the renormalized Hamiltonian, due to the pairing coupling between the system and the environment.
On the other hand, the mean-value deviations obey the following equations d dt where ∆n(t) = Tr s [δa † δaρ s (t)], ∆s(t) = Tr s [δaδaρ s (t)], ∆h(t) = Tr s [δaδa † ρ s (t)] = 1 + ∆n(t), and δa = a − a . The first term and its Hermitian conjugate in the right hand side of the above equation describe a damping oscillation of the quadrature variables, which follows the same dynamics as given by Eq. (43). The last term in Eq. (44) is the noise sources coming from the particle transitions and pair productions between the system associated with the initial state of the environment. Thus, γ(t) and γ(t) can be interpreted as the diffusion rates. We also find that the exact master equation (41) derived with the initial decoupled state has indeed the same form as the exact master equation for the initial correlated density matrix ρ tot (t 0 ) of Eq. (22) that we have derived in the previous work (i.e., Eq. (12) in Ref. [32]). In Eq. (22), H tot is the total Hamiltonian of the system and the environment plus their interaction, as given by Eq. (12). Thus, the initial state cannot be decoupled into a direct product of the system state with the environment state. Of course, to let the system and the environment evolve into a nonequilibrium evolution with such an initial thermal state, one needs to quench the system. A simple quench can be made by replacing the system frequency ω s with time-dependent one ω s (t), which is easy to be realized in experiments. All the formulas derived in this section remain the same forms with such a quenching modulation. As one can find from [32], with the initial correlated state Eq. (22), the only modification in the exact master equation is the fluctuation coefficients γ(t) and γ(t), which are still determined by the nonequilibrium Green function V(t, t) but V(t, t) of Eq. (40) is modified by Here G (es) (τ, τ ) is given by with n k (t 0 ) = a † (t 0 )b k (t 0 ) and s k (t 0 ) = a(t 0 )b k (t 0 ) being initial correlations between the system and the environment. The factor 2 comes from the integral t t0 dτ δ(τ − t 0 ) = 1/2. For the special case of W k = 0, Eq. (12) is reduced to the Fano Hamiltonian that describes a discrete state coupling to a continuous spectral bath. In this cases, the pairing-production or pairing-annihilation related coefficients vanish: γ(t) = 0 and ω s (t) = 0. Correspondingly, the exact master equation (41) is reduced to the one given in our earlier works [32,42,[49][50][51][52][53]. For another special case of W k = V k which corresponds to conventional QBM model [4,7], we will discuss the corresponding reduction in details in the next section.

IV. REPRODUCTION OF THE CONVENTIONAL QBM AND THE RESOLUTION OF THE INITIAL JOLT PROBLEM
In this section, we will show that the Hu-Paz-Zhang master equation of the QBM model [7] is a special case of our master equation Eq. (41) for the generalized QBM. We will also show that the long-standing initial jolt problem in the conventional QBM is artificial but it has nothing to do with the initial decoupled state.
Without including momentum-dependent couplings, the strength of the coupling and the pairing are the same W k = V k , as we shown in Eq. (10). This condition leads to the following relations between the renormalized frequency, dissipation (decay) and fluctuation (diffusion) coefficients Through these relations, the Hu-Paz-Zhang master equation of the conventional QBM can be reproduced easily from our general exact master equation Eq. (41) as a special limit. It can be expressed alternatively as where H R (t) is the complete renormalized system Hamiltonian, All the coefficients in the above master equation are deduced from the coefficients of the general master equation (41) under the condition W k = V k . The results are given by Notice that in the original Hu-Paz-Zhang master equation, the renormalized system Hamiltonian is incomplete, because not all commutator terms (associated with the unitary evolution) have been separated from the anticommutator terms (non-unitary evolution) in the master equation. This is how our complete renormalized system Hamiltonian is defined, as given by Eq. (49). It shows that the renormalized system Hamiltonian contains a potential of the position-momentum coupling. The reason that this momentum-dependent potential is misplaced into the dissipation (decay) term in the Hu-Paz-Zhang equation is because at the special case W k = V k , we have Im[ω s (t)] = γ(t)/2, see Eq. (47), so that this renormalization-induced potential has the same coefficient as the dissipation term in the master equation. On the other hand, this renormalization-induced momentum-dependent potential also indicates that including the momentum-dependent couplings between the system and the environment is a natural consequence, because even if only the position-position coupling is considered as in the conventional QBM, the renormalization not only changes the frequency of the system, but also induces a position-momentum coupling potential in the system Hamiltonian. As a consequence of this positionmomentum coupling potential, the physical frequency of the Brownian particle is ω 2 p (t) = ω 2 s + δω 2 s (t) − Γ 2 (t), which is consistent with the renormalized eigen-frequency ω r (t) = |ω s (t)| 2 − |ω s (t)| 2 in our general formalism, see the discussion after Eq. (42).
The decay and diffusion of the QBM can be manifested clearlier from the equation of motion for the quadrature covariance variables, which can be found from the master equation (48) directly: where ∆xp(t) = 1 2 {∆x(t), ∆p(t)} . Thus, the coefficients Γ(t)f (t) and Γ(t)h(t) can be interpreted as the diffusion rates of the quadrature covariance variables ∆xp(t) and ∆p 2 (t) , respectively. The missing diffusion rate for ∆x 2 (t) is due to the lack of positionmomentum coupling in the conventional QBM Hamiltonian.
Based on the above reproduction of Hu-Paz-Zhang master equation from our generalized QBM, we would also like to address a long-standing debate issue in the conventional QBM, called "initial jolt" [7,54,55]. The initial jolt is associated with a peak at the beginning of the diffusion coefficient Γ(t)h(t) at low temperature limit in the Hu-Paz-Zhang master equation, which is sensitive to the cut-off frequency Λ for the Ohmic-type spectral density. The inverse of this cut-off frequency, 1/Λ, is the typical time scale of the environment. The amplitude of the initial jolt is proportional to Λ, and it diverges when one takes a constant damping coefficient for the conventional QBM model. This corresponds to an Ohmic bath with an infinite cut-off frequency limit Λ → ∞. Such divergence causes covariances (∆p, ∆x) of the Brownian particle to increase immediately at the very beginning of the system evolution. In fact, all of Ohmic, sub-Ohmic and supra-Ohmic spectral baths are suffered from the same divergence [7,54,55].
In the literature [7,55], it has been argued that such initial jolt phenomenon is artificial and most likely comes from the assumption of the initial decoupled systemenvironment state one used in deriving the master equation. Here, we should prove that it is artificial but is not because of the use of the initial decoupled state. We find that not only the diffusion coefficients in the generalized QBM, the decay rate also has similar behavior. While, the decay rate is independent of the initial systemenvironment state. This can be seen from the combination of Eqs. (35), (40) and (42). It shows that the decay rate is determined purely by the same dissipation dynamics described by the retarded nonequilibrium Green function Eq. (35a) which is independent of the initial system-environment state. This indicates that the problem of "initial jolt" has nothing to do with the initial decoupled state that one often questioned [7,55]. Below we will show how does the serious divergence of the socalled "initial jolt" occurs and why it is not related to the choice of initial system-environment states.
Plug in the Ohmic spectral density with an infinite cut-off into our exact master equation, the serious divergence occurs in the decay coefficient γ(t) and the fluctuation (diffusion) coefficients γ(t), γ(t). To simplify the formulation, we take J V (ω) ≡ 2π k |V k | 2 δ(ω − ω k ) = πγ0 2Λ ωe −ω/Λ and ω s = 2γ 0 Λ/π. This setup corresponds to Ohmic spectral density with an exponential cut-off and zero renormalized frequency in the conventional QBM model. The Fourier transformation of the Ohmic spectral density is given by which has a pick in the time scale 1/Λ. This function is related to the integral kernel in the time-convolution equation for the generalized nonequilibrium Green function (35), which depicts back-action processes between the system and the environment. The strength and the time scale of memory effect of the system is directly related to the cut-off frequency Λ which is the key to making the decay and diffusion coefficients diverge.
To show how this divergence is purely associated with the cut-off frequency introduced in the Ohmic spectral densities, we evaluate the decay and diffusion coefficients γ(t) and γ(t) by varying the ratio of transition coupling strength and pairing strength |W k /V k | = α. Using Heisenberg equation of motion, one can easily find the equation of motion for the average occupation and the squeezed parameter in arbitrary Brownian particle state, d dt where n k (t) = a † (t)b k (t) and s k (t) = a(t)b k (t) are correlations between the system and the environment. The first term represents the free evolution of the system, and the last term represents the generalized information current flowing into the system. The later can be expressed explicitly as after solving the environmental equation of motion. Here is the quadrature correlation matrix of the average occupation and the squeezed parameter. The first part in the inegration is the diffusion current due to the initial distribution of the environment. The second part is the dissipation current due to the initial state of the system, and the backflow from the system. Comparing Eq. (44) with (54), the diffusion and decay coefficients γ(t), γ(t) and γ(t) can be also expressed as At the beginning of the evolution, due to the initial condition V(t 0 , t 0 ) = 0 and U(t 0 , t 0 ) = I, we only need to consider the first part in the integration, and the offdiagonal terms of U(τ, t) can be ignored. Then γ(t) and γ(t) in the time scale 1/Λ can be approximated as Here we only discuss the low temperature case k B T / Λ, so the terms proporsional tog V (t, τ ) ≡ k |V k | 2 n k (t 0 )e −iω k (t−τ ) and g W (t, τ ) ≡ k |W k | 2 n k (t 0 )e −iω k (t−τ ) have also been dropped.
We show numerically that the estimated results of the diffusion and decay coefficients γ(t) and γ(t) calculated from Eq. (56) are almost the same as the exact solutions of Eq. (42) for all the values of α: 0 ≤ α ≤ 1 at low temperature with a high cut-off frequency, see Fig. 2. When α = 1, the amplitudes of g V (t, τ ) and g W (t, τ ) become the same, then the divergences in the decay coefficient is cancelled, as analytically shown by Eq. (56b). Figure 2(a) clearly show that for α = 1, the decay coefficient γ(t) suddenly increases only a little bit within the time scale 1/Λ but then it approaches to γ 0 instead of decreasing. In other words, no "initial jolt" occurs in the dissipation γ(t) for the conventional QBM because of the accident cancellation of the divergences. For any other α values (which is more realistic in quantum optics), the divergence cannot be fully cancelled so that the decay coefficient γ(t) also shows the "initial jolt", as given in Fig. 2(a). Because the decay coefficient γ(t) is independent of the initial system-environment state. We conclude that the "initial jolt" has nothing to do with the initial decoupled state one used in deriving the exact master equation.
With the finding of this accident cancellation of the divergences in the decay coefficient γ(t) due to the special choice of α = 1 in the QBM, the origin of the initial jolt becomes clear. It comes purely from the Ohmictype spectral densities with the infinite cut-off frequency Λ → ∞. For a fundamental renormalizable theory, one requires that any renormalized quantity must be independent of regularization-scheme, so that the high energy cut-off scale can be taken to the infinite limit in the calculations of energy. However, both the conventional and the generalized QBM are build on phenomenological model Hamiltonians, no such renormalizability is required in priori. Furthermore, the cut-off frequency in the Ohmictype spectral densities is introduced to count the effective transition processes between the system and the environment. Quantum mechanically, the most active environmental modes ω k that effectively couple to the system are dominated by these not far away from the thermal frequency k B T / so that physical transition between the system and the environment can occur (with large probabilities). Physically, the spectral density is a summation of the probabilities over the physical processes between the system and the environment. For the Ohmic-type spectral densities, taking a very large cut-off frequency implies that the environmental modes with higher energy will have the larger probability amplitude for coupling to the system mode. This is obviously very unphysical. The physical reliable cut-off frequency should be estimated by the condition that: the mean frequency of a given spectral density should be the same order of the system characterized frequency, if the system-environment coupling cannot be derived from a more fundamental theory or no experimental data can be used to fix it. Under such a physical requirement, the problem of "initial jolt" will never occur.

V. DISCUSSIONS AND PERSPECTIVES
In this paper, we provide a detailed analysis on systemenvironmental couplings which are usually unknown in priori but are the key to the determination of dissipation and fluctuation dynamics in open quantum systems. In particular, we show that the quantum dissipative dynamics described by Feynman-Vernon model and Caldeira-Leggett model based on classical mechanics involve some physically inconsistent processes at quantum level. These inconsistent processes include the transitions from low-energy states to high-energy states by emitting an amount of energy or create particles from nothing, which should be forbidden in quantum mechanics. Within the framework of quantum mechanics, each individual microscopic particle process should be constrained by the conservations of energy-momentum and angular momentum, etc. We further show that one can get rid of these inconsistent processes by including momentumdependent system-environment couplings. We generalize the QBM to include such momentum-dependent couplings, which are indeed easier to be realized in realistic quantum world. We derive the exact master equation of such generalized QBM for both the initially decoupled and initially correlated system-environment states.
With the generalized QBM and its exact master equation for both the initially decoupled and the initially correlated state, we reproduce the Hu-Paz-Zhang master equation as a special case. We then find that in the Hu-Paz-Zhang equation for the conventional QBM with the position-dependent coupling only, the renormalized Brownian particle Hamiltonian actually contains a renormalization-induced momentum-dependent potential. This momentum-dependent potential is misplaced into the dissipation term in the Hu-Paz-Zhang equation, so that the correct renormalized Brownian particle Hamiltonian was not found before. The environment induced momentum-dependent potential in the Hu-Paz-Zhang equation also indicates that including the momentum-dependent system-environment coupling is a natural consequence of the QBM.
Furthermore, we re-examine the long-standing problem of the initial jolt in the conventional QBM and Hu-Paz-Zhang master equation. The initial jolt is related to the divergence arising from the Ohmic-type spectral density with an infinite cut-off frequency. It has been thought that the initial jolt is a result of using the initial decoupled state which is physically unacceptable, but we find that this is indeed a misunderstanding. The misunderstanding comes from the accident cancellation of the divergences in the decay coefficient in the Hu-Paz-Zhang equation, due to the same coupling strengths for the single particle transition and the pairing processes. With the generalized QBM including the momentumdependent coupling, we show that the initial jolt exists in both the decay and diffusion coefficients, while the decay coefficient is independent of the initial systemenvironment state. As a conclusion, the so-called "initial jolt" has nothing do to with the initial decoupled state. It is an artificial effect when one takes an extremely large cut-off frequency in comparison with the thermal energy of the environment, which is unphysical as we explained in the end of the last section.
With the generalized QBM, we have resolved these interesting and important problems for quantum dissipative dynamics, as we discussed in this work. The new exact master equation for the generalized QBM also has the potential applications to photonics quantum computing. With an extension of the generalized QBM to multimodes, one can apply it to the integrated quantum processor made of many squeezers and interferometers with hundred more photonic modes, that currently attracted a great attention in the development of quantum technology. We will leave this problem for further investigation.
The nonequilibrium Green functions we introduced for open quantum systems [40,[42][43][44] are generalized to open quantum systems involving pairing couplings in this work, and are denoted by U(τ, t 0 ) and V(τ, t) [32]. They are defined as These equations of motion determine all the dissipation and fluctuation dynamics of the Brownian particle. Formally, solving Eq. (A3b), we obtain This formal solution allows us to eliminate exactly and completely all the environmental degrees of freedom (equivalent to the trace over all the environment states). Substituting this solution into Eq. (A3a), we have d dt where is Eq. (30a). One can clearly see that Eq. (A5) is the generalized quantum Langevin equation, in which the third term in the left hand side is the dissipation (decay) term due to the system-environment coupling, and the term in the right hand side is the fundamental noise force generated by the environment through the systemenvironment coupling. It is obvious that the average of the noise force is zero if the environment is initially in a thermal state. Because Eq.(A5) is a linear differential equation, its general solution has the form Substituting this general solution into Eq. (A5), one can find that the retarded Green function U(t, t 0 ) describes the dissipation dynamics of the Brownian particle and obeys the Dyson equation of motion, d dτ U(τ, t 0 )+iω s ZU(τ, t 0 )+ The last term in Eq. (A7) describes the environmentinduced fluctuation (noise) dynamics that obeys the following equation of motion, d dt subjected to the initial condition F (t0) F † (t0) = 0. From this initial condition, the general solution of the equation for the noise dynamics is given explicitly, .

(A11)
The nonequilibrium correlation functions are defined by a † (t)a(τ ) a(t)a(τ ) a † (t)a † (τ ) a(t)a † (τ ) = U(τ, t 0 ) a † (t 0 )a(t 0 ) a(t 0 )a(t 0 ) a † (t 0 )a † (t 0 ) a(t 0 )a † (t 0 ) U † (t, t 0 ) The last term is the noise-induced correlation Green function V(τ, t), which has the solution This gives the same solution Eq. (37) obtained in our exact master equation derivation. The initialstate-dependent system-environment correlation function G(t, t ) is given by which is Eq. (30b). The solution Eq. (A13) is also the generalization of the Keldysh's correlation function in the nonequilibrium technique [45]. It is indeed the generalized nonequilibrium fluctuation-dissipation theorem, as a consequence of the unitarity principle of quantum mechanics for the whole system (the system plus the environment) [44]. The usual equilibrium fluctuationdissipation theorem can be easily obtained from this solution as a steady-state limit.