Exact dynamics and thermalization of open quantum systems coupled to reservoir through particle exchanges

In the realistic world, physical systems are inevitably coupled to environments, which makes the theory of open quantum systems vastly used in many fields of physics, chemistry, biology, and engineering. The systematic study of open quantum systems has aroused the researchers’ interest since the 1960s Feynman and Vernon (1963); Schwinger (1961); Zwanzig (1960); Nakajima (1958); Breuer et al. (2002); Weiss (2012); Gardiner and Zoller (2004), and becomes more and more important for the prosperously developing field of quantum information processing Gardiner and Zoller (2004); Wiseman and Milburn (2009); Nielsen and Chuang (2002), quantum transport theory Haug and Jauho (2008); Yang and Zhang (2016), and rapidly improving time-resolved measurement technologies Wickenhauser et al. (2005); Kaldun et al. (2016). One of the most crucial problems in dealing with open quantum systems is how to determine explicitly the evolution of open quantum system states, through which all the information about the system dynamics can be obtained. However, due to the contamination of the huge environment, the system dynamics is non-unitary and always involves complicated fluctuation and dissipation. As a consequence, up to date, most open quantum systems can only be dealt with perturbative methods, such as Born-Markov approximation or cutting-off in Nakajima-Zwanzig operator projective technique Nakajima (1958); Zwanzig (1960); Breuer et al. (2002). Only a few classes of open quantum systems, e.g., the harmonic oscillator in quantum Brownian motion Caldeira and Leggett (1983); Hu et al. (1992) and the open systems interacting with particle-exchange coupling  Zhang et al. (2012), can an exact master equation be obtained, let alone the exact state evolution.

In a series of papers Tu and Zhang (2008); Jin et al. (2010); Lei and Zhang (2012); Zhang et al. (2012), we have obtained the general exact master equation for the open quantum systems whose interaction only involves particle exchanges. This large class of open systems is an extension of the famous Fano-Anderson model Anderson (1961); Fano (1961) and characterizes many physical phenomena in different systems, such as Fano resonance in atomic and condensed matter physics Miroshnichenko et al. (2010), Anderson localization in many-body systems Anderson (1958), photon-atom bound states in photonic crystal Yablonovitch (1987); John (1987); Kofman et al. (1994), quantum transport in quantum dot junctions Haug and Jauho (2008), impurity defects in solids Mahan (2013), etc. Although it has been studied for decades, the exact general master equation has not been derived until very recently Tu and Zhang (2008); Jin et al. (2010); Lei and Zhang (2012); Zhang et al. (2012). With the progress, various perspectives of the open systems can be studied, e.g., memory effects Yang et al. (2013), entanglement dynamics Tan et al. (2011); Lin et al. (2016); An and Zhang (2007), decoherence Tu and Zhang (2008); Liu et al. (2016); Yang and Wu (2014), the fluctuation-dissipation theorem Zhang et al. (2012), exact transient quantum transport Yang and Zhang (2018); Lo et al. (2015); Yang and Zhang (2016), etc.

In this paper, one main result we obtained is the solution of the exact master equation, i.e., the exact form of the reduced density matrix at an arbitrary later time. Note that the time-dependence of the reduced density matrix includes complete information about the system dynamics, it offers more perspectives in studying the memory effects Li et al. (2018); Breuer et al. (2016); Rivas et al. (2014), entanglement dynamics and dynamical phase transition Heyl (2018) of such class of open systems. Using the result of the reduced density matrix, we also study the general thermalization in this paper, which is another important and hot topic for both experimentalists Trotzky et al. (2012); Gring et al. (2012) and theorists Deutsch (1991); Kosloff (2013); Ali et al. (2018); Xiong et al. (2015); Srednicki (1994); Linden et al. (2009); Short and Farrelly (2012); Reimann (2008); Rigol (2009); Polkovnikov et al. (2011); Cazalilla and Rigol (2010); Hsiang et al. (2018) in recent years. In equilibrium statistic mechanics, the thermal distribution is based on the assumption of equal probability of all permissible microstates Callen (1998). However, the foundation of thermalization has always been a tough problem and a long-term goal of physicists Polkovnikov et al. (2011). In the past decade, by taking advantage of the eigenstate thermalization hypothesis, researchers study the thermalization of closed quantum systems and make interesting contributions in understanding the physics beneath it Deutsch (1991); Kosloff (2013); Ali et al. (2018); Srednicki (1994); Linden et al. (2009); Short and Farrelly (2012); Reimann (2008); Rigol (2009); Polkovnikov et al. (2011); Cazalilla and Rigol (2010). In contrast, in this paper, we discuss the thermalization of open systems from the dynamical perspective. With the exact evolution of the open quantum systems, such an aim is achieved.

The rest of the paper is organized as follows. In Sec. II, we introduce the system we concerned about and briefly review the previous results related to our present work. In Sec. III, we shall derive the exact solution of the reduced density matrix. In Sec. IV, we explain the physical consequences of the solution and discuss the general thermalization of the open systems. A brief summary will be given in Sec. V.

The system we are interested are characterized by the Hamiltonian $H_{{\mathrm{S}}}=\sum_{i,j=1}^{d}\epsilon_{ij}a^{\dagger}_{i}a_{j}$, and interacts with bosonic (fermionic) environment described by $H_{\mathrm{E}}=\sum_{k}\epsilon_{k}b^{\dagger}_{k}b_{k}$, via the interaction $H_{\mathrm{int}}=\sum_{j,k}(V_{jk}a^{\dagger}_{j}b_{k}+\mathrm{H.c.})$, where $a_{j}^{\dagger}$ ($a_{j}$) is the creation (annihilation) operator of the $j$th level in the system; $b_{k}^{\dagger}$ ($b_{k}$) is the creation (annihilation) operator of the $k$-mode in the environment; $\epsilon_{ij}$ and $\epsilon_{k}$ characterize the energy levels of the system and the environment, and $V_{jk}$ is the interaction strength between them. If the creation and annihilation operators satisfy the commutation (anti-commutation) relations, the corresponding systems is bosonic (fermionic).

To make the notations more compact, we denote that $\bm{a}^{\dagger}=(\!\begin{array}[]{ccc}a^{\dagger}_{1}&\cdots&a^{\dagger}_{d}\end{array}\!)$, $\bm{a}=(\!\begin{array}[]{ccc}a_{1}&\cdots&a_{d}\end{array}\!)^{{\mathrm{T}}}$, and the matrix constituted of ${\epsilon}_{ij}$ as $\bm{\epsilon}_{\mathrm{S}}$. The spectrum of the environment and the system-environment interaction can be specified with the spectral density function $\bm{J}(\epsilon)$, whose elements are defined through $J_{ij}(\omega)=2\pi\sum_{k}V_{ki}V^{*}_{kj}\delta(\epsilon-\epsilon_{k})$. We assume that initially the system and the environment is decoupled Leggett et al. (1987), i.e., $\rho_{{\mathrm{tot}}}(0)=\rho(0)\otimes\rho_{\mathrm{E}}(0)$, where $\rho(0)$ is an arbitrary physical state of the system and $\rho_{\mathrm{E}}(0)=\otimes_{k}e^{-\beta(\epsilon_{k}-\mu)b_{k}^{\dagger}b_{k}}/{\cal Z}$ is the thermal state of the environment with inverse temperature $\beta$ and chemical potential $\mu$. Here, ${\cal Z}=\prod_{k}|1\mp e^{-\beta[\epsilon_{k}-\mu]}|^{\mp 1}$ denotes the grand partition function of the environment with the upper and lower sign corresponding to the bosonic and fermionic case respectively. Hereafter, we also use this convention.

The reduced density matrix at latter time $t$ is connected to the initial state through the propagating function ${\cal{J}}(\bm{\eta}^{\ast},\bm{\eta},t;\bm{\zeta}^{\prime\ast},\bm{\zeta},0)$ in the coherent state representation, i.e.,

where $\bm{\zeta}^{\ast}=(\!\begin{array}[]{ccc}\zeta^{\ast}_{1}&\cdots&\zeta^{\ast}_{d}\end{array}\!)$ and $\bm{\zeta}=(\!\begin{array}[]{ccc}\zeta_{1}&\cdots&\zeta_{d}\end{array}\!)^{{\mathrm{T}}}$, with the components being complex (Grassmannian) for bosons (fermions); $\bra{\bm{\zeta}^{*}}=\bra{0}e^{\bm{\zeta}^{*}\bm{a}}$ ($\ket{\bm{\zeta}}=e^{\bm{a}^{\dagger}\bm{\zeta}}\ket{0}$) is the unnormalized coherent state with $\bra{0}$ ($\ket{0}$) standing for the vacuum state Zhang et al. (1990); ${\mathrm{d}}\mu(\bm{\zeta})$ is the integral measure, and ${\mathrm{d}}\mu(\bm{\zeta})=\prod_{j=1}^{d}({{\mathrm{d}}^{2}\zeta_{j}}/{\pi})e^{-\zeta^{*}_{j}\zeta_{j}}$ for bosons and ${\mathrm{d}}\mu(\bm{\zeta})=\prod_{j=1}^{d}{{\mathrm{d}}\zeta^{*}_{j}{\mathrm{d}}\zeta_{j}}e^{-\zeta^{*}_{j}\zeta_{j}}$ for fermions. With coherent state path-integral approach, the propagating function can be derived as Jin et al. (2010); Lei and Zhang (2012)

where ${\widetilde{\bm{A}}}=1-\bm{u}^{{\dagger}}\frac{1}{1\pm\bm{v}}\bm{u}$, $\widetilde{\bm{u}}=\frac{1}{1\pm\bm{v}}\bm{u}$, with the matrices $\bm{u}$ and $\bm{v}$ being respectively the abbreviations of the nonequilibrium Green functions $\bm{u}(t)$ and $\bm{v}(t,t)$. More specifically, $\bm{u}(t)$ is the spectral Green function with elements defined as $u_{ij}(t)=\langle[a_{i}(t),a_{j}^{\dagger}(0)]_{\mp}\rangle$ Kadanoff (2018), and $\bm{v}(t,t)$ is related to the system’s particle number originated from the environment, reading Zhang et al. (2012)

Here, $\widetilde{\bm{g}}(t_{1},t_{2})=\int\frac{{\mathrm{d}}\epsilon}{2\pi}f(\epsilon)\bm{J}(\epsilon)e^{-i\epsilon(t_{1}-t_{2})}$ denotes a system-bath correction, with $f(\epsilon)=\frac{1}{e^{\beta(\epsilon-\mu)}\mp 1}$ standing for the initial particle distribution of the environment.

Taking advantage of Eqs. (1)-(2), the master equation has been derived Tu and Zhang (2008); Jin et al. (2010); Lei and Zhang (2012); Zhang et al. (2012), reading

In this master equation, $\widetilde{H}_{\mathrm{S}}(t)=\sum_{ij}\widetilde{\epsilon}_{ij}(t)a^{\dagger}_{i}a_{j}$ is the renormalized system Hamiltonian; $\gamma_{ij}(t)$ and $\widetilde{\gamma}_{ij}(t)$ characterize the dissipation and fluctuation rate induced by the system’s coupling with the environment, respectively. The coefficients $\widetilde{\epsilon}_{ij}(t)$, $\gamma_{ij}(t)$ and $\widetilde{\gamma}_{ij}(t)$ are determined only by the non-equilibrium Green functions $\bm{u}(t)$ and $\bm{v}(t,t)$ Zhang et al. (2012).

Furthermore, the Green function $\bm{u}(t)$ satisfies the Dyson equation and the solution can be expressed as Zhang et al. (2012)

The quantities $\epsilon_{l}$, $\bm{Z}_{l}$ and $\bm{D}(\epsilon)$ are all connected with the Green function in the energy domain

where $\bm{\Sigma}(z)=\int\frac{{\mathrm{d}}\epsilon}{2\pi}\frac{\bm{J}(\epsilon)}{z-\epsilon}$ is the self-energy correction. $\epsilon_{l}$ is the energy of the $l$th localized mode and is determined by the singularities of $\bm{U}(z)$, i.e., it satisfies the equation $\det(\epsilon_{l}-\bm{\epsilon}_{\mathrm{S}}-\bm{\Sigma}(\epsilon_{l}))=0$. $\bm{Z}_{l}$ is the Hermitian matrix characterizing the amplitude of the $l$th localized mode, reading $\bm{Z}_{l}=\frac{1}{2\pi\mathrm{i}}\oint_{C_{l}}\!{\mathrm{d}}z\,\bm{U}(z)$ , where $C_{l}$ is the positively-oriented curve in the neighbourhood of $z=\epsilon_{l}$. $\bm{D}(\epsilon)$ is the spectrum of the system broadened by the environment and reads $\bm{D}(\epsilon)=\bm{U}(\epsilon+\mathrm{i}0^{+})\bm{J}(\epsilon)\bm{U}(\epsilon-\mathrm{i}0^{+})$ ¹¹1The $\bm{Z}_{l}$’s and $\bm{D}(\epsilon)$ are the discrete and continuous components of the modified spectrum of the system $\bm{D}_{\mathrm{S}}(\epsilon)$, respectively. They are connected through the relation $\bm{D}_{\mathrm{S}}(\epsilon)=2\pi\sum_{l}\bm{Z}_{l}\delta(\epsilon-\epsilon_{l})+\bm{D}(\epsilon)$, where $\bm{D}_{\mathrm{S}}(\epsilon)$ is defined through $\bm{D}_{\mathrm{S}}(\epsilon):=-2\imaginary[\bm{U}(\epsilon+\mathrm{i}0^{+})]$ (See Ref. Zhang et al. (2012))..

With the initial system state $\rho(0)$, and the Green functions $\bm{u}(t)$ and $\bm{v}(t,t)$, the reduced density matrix $\rho(t)$ at arbitrary instant can be determined from Eq. (II) in principle. However, because the coefficients are all time-dependent, the problem is too complicated to solve directly. Luckily, one can solve it by taking advantage of the propagating function in Eqs. (1)-(2), as is shown in this section. First, in Sec. III.1, we shall introduce some conventions and derive $\rho(0)$’s elements in the coherent state representation. In Sec. III.2, we obtain $\rho(t)$’s elements at arbitrary later time. In Sec. III.3, the exact form of $\rho(t)$ will be given.

In this section, we shall discuss the physics beneath the solution. The physical interpretation of the solution will be given in Sec. IV.1, and the thermalization problem will be discussed in Sec. IV.2.

In this paper, we have investigated a general solution of open quantum systems interacting with the environment through particle exchanges. The exact evolution of the reduced density matrix is given in terms of the nonequilibrium Green functions. We explained the physical consequences of the solution. With the exact density matrix, we study the thermalization process. We obtain the result of equilibrium statistical mechanics from the dynamical perspective and go beyond it. That is, when there are no localized modes and the system-environment coupling is very weak, the final state would be as expected from the conventional statistical mechanics; for no localized modes but relatively strong coupling regime, the steady state would be thermal-like, which departures from the prediction of conventional statistical mechanics; when there are localized modes, the system keeps the memory of the initial state and can not be thermalized.

With the explicit expression of the reduced density matrix, one can obtain the complete information about the system dynamics, which is quite important for the rapidly developing quantum thermodynamics and quantum information, because their central physical quantity, entropy, is directly related to the state. It is also noteworthy that the model studied in our work involves non-Markovian nature. With the explicit form of the density matrix evolution, one can study the memory effects from more perspectives, e.g., quantum coherence, entanglement, and dynamical phase transition. Although we only consider the single-reservoir case, our result can be directly extended to the multi-reservoir case by just extending the corresponding expressions of the nonequilibrium Green functions $\bm{u}(t)$ and $\bm{v}(t,t)$ to multi-reservoirs (multi-leads in nano/quantum devices). Therefore, it is also easy to apply to quantum transport theory.