Mean Field Theory of Thermal Energy Transport in Molecular Junctions

TL;DR

This paper applies mean field theory and Ehrenfest trajectories to model thermal energy transport in molecular junctions, demonstrating robustness and accuracy in regimes challenging for other methods.

Contribution

It introduces a trajectory-based mean field approach for nonequilibrium thermal transport, validated against quantum benchmarks and applicable to complex regimes.

Findings

Ehrenfest approach yields accurate steady-state heat currents.

Method performs well in intermediate and strong coupling regimes.

Robustness surpasses many perturbative quantum transport methods.

Abstract

Mean field theory is applied to nonequilibrium thermal energy transport in a model molecular junction. An approximation to the total time-dependent heat current in the junction is constructed using an ensemble of Ehrenfest trajectories, and the average heat current in the steady state is obtained. The accuracy of this treatment is verified through benchmark comparisons with exact quantum mechanical results, and various approximate quantum transport theories, for the nonequilibrium spin-boson model. The performance of the multi-trajectory Ehrenfest approach is found to be quite robust, displaying good accuracy in intermediate cases that remain elusive to many perturbative approximations, and in the strong coupling limit where many methods break down. Thus, mean field theory and related trajectory-based approximate quantum dynamics methods emerge as a promising toolkit for the study of transport properties in nanoscale systems.

Equations

$$\hat{H}={\hat{H}}_s+{\hat{H}}_b+{\hat{H}}_{sb}.$$

$${\hat{H}}_s=\frac{ϵ}{2}{\hat{\sigma }}_z+\frac{\Delta }{2}{\hat{\sigma }}_x,$$

$$\hat{H}_{b}=\frac{1}{2}\sum_{\lambda=L,R}\sum_{k}\Big{[}\hat{P}_{\lambda,k}^{2}+\omega^{2}_{\lambda,k}\hat{Q}^{2}_{\lambda,k}\Big{]}.$$

$$\begin{array}{rl}{\hat{H}}_{sb}& =\sum _{\lambda =L,R}\sum _k{c}_{\lambda ,k}{\hat{Q}}_{\lambda ,k}{\hat{\sigma }}_z.\end{array}$$

$$g_{\lambda }(\omega )=\frac{\pi }{2}\sum _k\frac{c_{\lambda ,k}^2}{{\omega }_{\lambda ,k}}\delta (\omega -{\omega }_{\lambda ,k}),$$

$$g_{\lambda }(\omega )=\frac{\pi }{2}{\xi }_{\lambda }\omega \exp (-\omega /{\omega }_{c,\lambda }).$$

$$\frac{\partial }{\partial t}{\hat{\rho }}_W(X,t)=-i\mathcal{L}{\hat{\rho }}_W(X,t).$$

$${\hat{\rho }}_W(R,P)=\frac{1}{(2\pi \hslash {)}^{N_b}}\int dZ{e}^{iP\cdot Z}⟨R-\frac{Z}{2}|\hat{\rho }|R+\frac{Z}{2}⟩.$$

$$i\mathcal{L}\cdot =\frac{i}{\hslash }[{\hat{H}}_W,\cdot ]-\frac{1}{2}(\{{\hat{H}}_W,\cdot \}-\{\cdot ,{\hat{H}}_W\}),$$

$${\hat{\rho }}_W(X,t)={\hat{\rho }}_s(t){\rho }_{b,W}(X,t),$$

$$\hat{\rho}_{s}(t)=Tr_{b}\Big{(}\hat{\rho}(t)\Big{)}=\int dX\hat{\rho}_{W}(X,t),$$

$$\frac{d}{dt}\hat{\rho}_{s}(t)=-i\Big{[}\hat{H}_{s}+\hat{H}_{sb,W}(X(t)),\hat{\rho}_{s}(t)\Big{]}.$$

$$\frac{\partial R_{\alpha }}{\partial t}=\frac{\partial H_{b,W}^{Eff}}{\partial P_{\alpha }},\frac{\partial P_{\alpha }}{\partial t}=-\frac{\partial H_{b,W}^{Eff}}{\partial R_{\alpha }}.$$

$$H^{Eff}_{b,W}=H_{b,W}+Tr_{s}\Big{(}\hat{H}_{sb,W}(X,t)\hat{\rho}_{s}(X,t)\Big{)}.$$

$$⟨O(t)⟩=T{r}_s\int dX{\hat{O}}_W(X,t){\hat{\rho }}_W(X,0).$$

$$\langle H_{b,\lambda}(t)\rangle=Tr\Big{(}\hat{H}_{b,\lambda}(t)\hat{\rho}(0)\Big{)}.$$

$$\hat{\rho }(0)={\hat{\rho }}_{b,L}^{eq}\otimes {\hat{\rho }}_{b,R}^{eq}\otimes {\hat{P}}_s^{22},$$

$$\begin{split}(\rho_{b,\lambda}^{eq})_{W}&=\prod_{k}\frac{\tanh(\beta_{\lambda}\omega_{\lambda,k}/2)}{\pi}\\ &\ \times\exp\Bigg{[}-\frac{\tanh(\beta_{\lambda}\omega_{\lambda,k}/2)}{\omega_{\lambda,k}}\Big{[}P_{\lambda,k}^{2}+\omega_{\lambda,k}^{2}R_{\lambda,k}^{2}\Big{]}\Bigg{]},\end{split}$$

$$J_L$$

$$J_R$$

$$J=\lim_{t\rightarrow\infty}J(t)=\lim_{t\rightarrow\infty}\frac{1}{2}\Big{(}J_{L}(t)+J_{R}(t)\Big{)}.$$

Full text

I Introduction

With rapid recent developments in molecular-scale electronics and a plethora of emerging nanotechnologies Schwab et al. (2000); Wang et al. (2007a); Losego et al. (2012); Meier et al. (2014), theoretical modelling of nonequilibrium quantum charge and thermal energy transport processes is becoming of increasing interest and importance Montgomery, Todorov, and Sutton (2002); Chen, Zwolak, and Di Ventra (2003); Craven and Nitzan (2016); Carbogno, Ramprasad, and Scheffler (2017). From the perspective of molecular simulation, progress in this area requires moving toward the treatment of real-time dynamics in realistic systems. This motivates the development of accurate and highly efficient methods to treat open quantum systems that can account for the effects of anharmonicity and nonadiabaticity beyond the perturbative and Markovian limits.

In the cases where they can be converged, exact quantum mechanical solutions offer the benefit of unambiguous insight into the problem at hand. The current drawback is that range of systems for which exact solutions can be obtained remains rather narrow, and is essentially limited to generic systems with very few (less than twenty) degrees of freedom, or the class of quantum impurity models, which can have many degrees of freedom but a restricted functional form. Although the form of the impurity models may appear simple, they form the basis for simulations of a vast range of properties of condensed matter systemsGeorges et al. (1996); Kotliar et al. (2006). Specifically, here we consider the nonequilibrium spin-boson model, which is a quantum impurity model that has become a workhorse for investigating the thermodynamics of open quantum systemsEsposito, Harbola, and Mukamel (2009); Campisi et al. (2011); de Vega and Alonso (2017) and nonequilibrium thermal energy transferSegal and Nitzan (2005a); Segal and Agarwalla (2016). As such, a number of exact solvers have been developedVelizhanin, Wang, and Thoss (2008); Saito and Kato (2013); Chen, Cohen, and Reichman (2017); Makri (2017); Ridley et al. (2018) that allow access to ever an increasing range of parameter space within this model. Nevertheless the need for approximate techniques that can go beyond this class of models is also clear. As such, these systems also serve to provide a rich testing ground for assessing the accuracy of approximate methods, in addition to providing physical insight into more realistic systems

The application of approximate quantum dynamics techniques to real-time nonequilibrium thermal transport was pioneered by Segal and Nitzan, beginning with their work on molecular junctionsSegal, Nitzan, and Hänggi (2003); Segal and Nitzan (2005b, c). This opened up a broad interest in this class of nonequilibrium transport problems, which has now become an active area of development and application for both exact quantum dynamics solversVelizhanin, Wang, and Thoss (2008); Segal (2013a); Saito and Kato (2013); Kato and Tanimura (2015, 2016); Song and Shi (2017); Kilgour, Agarwalla, and Segal (2018) and approximate quantum transport theories including nonequilibrium Green’s functions (NEGF) approachesMingo (2006); Wang, Wang, and Zeng (2006); Wang et al. (2007b); Ruokola and Ojanen (2011); Velizhanin, Wang, and Thoss (2008, 2010), path integral techniques such as the noninteracting blip approximation (NIBA)Nicolin and Segal (2011); Segal (2013b) and other quantum master equationsWang, Ren, and Cao (2015); Liu et al. (2017); Agarwalla and Segal (2017), as well as an approximation to the quantum-classical Liouville equationLiu et al. (2018).

Here, mean field theory is employed to study thermal energy transport in a molecular junction, represented by the nonequilibrium spin-boson (NESB) model, in order to address the efficacy of quantum-classical trajectory-based dynamics methods for this problem. The remainder of this work is organized as follows: in Section II a brief overview of the NESB model is given, along with a description of the how mean field theory is employed to calculate time dependent heat currents. Simulation results for nonequilibrium steady state heat currents for various parameter regimes of the NESB model are reported in Section III, where a comparison with exact numerical benchmarks are reported where possible, as well with a range of other approximate quantum dynamics approaches. In Section IV an outlook to potential future work is offered, and some concluding remarks are given.

II Theory

II.1 Nonequilibrium Spin-Boson Model

The nonequilibrium spin - boson model model forms the basis of many recent studies of thermal energy transfer in nanoscale molecular junctions. Within this model, two vibrational states of the junction-bound molecule are considered, and each state is coupled to two phonon baths, labelled left (L) and right (R) respectively. The total Hamiltonian for the NESB system can be written in a standard system-bath form,

[TABLE]

where $\hat{H}_{s}$ describes the two-level subsystem,

[TABLE]

and $\epsilon$ and $\Delta$ are related to the coupling between the two vibrational states of the molecule, and their energy difference, respectively.

The Hamiltonian operator for the environment includes both the left and right baths, written using mass-weighted coordinate and momenta operators is

[TABLE]

The molecular subsystem is bilinearly coupled to each bath,

[TABLE]

This coupling is characterized by a spectral density function,

[TABLE]

and, as is common in the literatureSegal and Nitzan (2005a); Nicolin and Segal (2011); Velizhanin, Wang, and Thoss (2008, 2010), we choose the Ohmic form in this study.

[TABLE]

Here the Kondo parameter, $\xi_{\lambda}$, is a measure of the system-bath coupling strength at the $\lambda^{th}$ contact, and $\omega_{c,\lambda}$ is the inverse of the corresponding characteristic timescale of each bath. In this study the right and left bath spectral densities are chosen to be identical, i.e. $\xi_{L}=\xi_{R}=\xi$, and $\omega_{c,L}=\omega_{c,R}=\omega_{c}$, such that the baths only differ in terms of their temperature.

II.2 Mean Field Theory

The multi-trajectory Ehrenfest mean field theory (MFT) can be derived in a straightforward manner, via the quantum-classical Liouville (QCL) equationKapral and Ciccotti (1999); Grunwald, Kelly, and Kapral (2009). The QCL equation of motion for the density matrix is formally exact for an arbitrary quantum mechanical system that is bilinearly coupled to a harmonic environment, as is the case for the NESB problem, and many other quantum impurity models. Written in a compact form, the QCL equation is

[TABLE]

The QCLE describes the time evolution of $\hat{\rho}_{W}(X,t)$, the partial Wigner transform of the density operator taken over the coordinates of the $N_{b}$ bath degrees of freedom, which are represented by continuous phase space variables $X=(R,P)=(R_{1},R_{2},...,R_{N_{b}},P_{1},P_{2},...,P_{N_{b}})$. The partial Wigner transform of the density operator, $\hat{\rho}$, is defined as

[TABLE]

The QCL operator is defined as

[TABLE]

where $[\cdot,\cdot]$ is the commutator, and $\{\cdot,\cdot\}$ is the Poisson bracket in the phase space of the environmental variables. The Ehrenfest MFT equations of motion are obtained by requiring that the total density of the system can be written as an uncorrelated product of the system and bath reduced densities at all times,

[TABLE]

where the reduced density matrix of the system is

[TABLE]

and the bath density is $\rho_{b}(X,t)=Tr_{s}(\hat{\rho}_{W}(X,t))$. Requiring solutions to the QCL equation of this form yields the Ehrenfest MFT equations of motion for the subsystem:

[TABLE]

The evolution of the reduced Wigner density of the bath can be represented by an ensemble of multiple independent trajectories, $\rho_{b,W}(X,t)=\sum_{j}\delta(X_{j}-X(t))$, that evolve according to a set effective Hamilton’s equations of motion generated from the mean-field effective Hamiltonian,

[TABLE]

The effective Hamiltonian in the classical equations of motion is,

[TABLE]

The exact expression for the average value of any observable, $\langle O(t)\rangle$, can be written as

[TABLE]

The mean field limit of this expression simple corresponds to evaluating the integral by sampling initial conditions for an ensemble of independent trajectories from $\hat{\rho}_{W}(X,0)$, and then generating the time evolution for each trajectory by approximating $\hat{O}_{W}(X,t)$ by it’s mean-field counterpart.

II.3 Thermal Energy Transport: Observables of Interest

The transport of thermal energy in the NESB model can be monitored via time-dependence of the average energy of each bath,

[TABLE]

Initially, we imagine that the system is separable; each bath is assumed to be in a canonical equilibrium state with temperature $T_{\lambda}$, and the two-level system is initially in the excited state. This corresponds to the following density operator,

[TABLE]

where $\hat{P}_{s}^{\alpha\alpha^{\prime}}=|\alpha\rangle\langle\alpha^{\prime}|$ is a subsystem projector. The bath initial conditions are sampled from the Wigner transform of the initial canonical density operators, with inverse temperature $\beta_{\lambda}=(k_{B}T_{\lambda})^{-1}$,

[TABLE]

The nonequilibrium heat currents flowing through the junction are then defined as the time rate of change of the respective average bath energies,

[TABLE]

and the total nonequilibrium steady state (NESS) heat current is then defined as the long time limit of the total heat current through the junction,

[TABLE]

III Results and Discussion

The results of multi-trajectory Ehrenfest mean field theory simulations of the NESS thermal energy transfer rate in the NESB model, with $\epsilon=0$, are summarized. In all simulations, at $t=0$ the two uncoupled reservoirs are then brought into contact with the subsystem. After a short period of transient dynamics, thermal energy steadily flows through the junction along the direction of the thermal gradient, as shown in Fig. 2. The steady-state heat current is reached as the rate of the energy flux between the left and right baths reaches a constant value, which is depicted in panels (b - d) of Figure 2. Somewhat surprisingly, the statistical properties of the thermal energy transfer dynamics are different than those of the population dynamics of the junction-bound molecular system. For example, the MFT simulations are found to be quite well converged for molecular properties using $\sim 10^{4}$ trajectories, while the corresponding heat current calculation requires two to three orders of magnitude more trajectories for convergence.

The dependence of the NESS heat current on the vibrational splitting of the molecular system, $\Delta$, is depicted in Figure 3. In this case, the MFT results are compared with numerically exact multi-layer multi-configurational time-dependent Hartree (ML-MCTDH) benchmark calculationsVelizhanin, Wang, and Thoss (2008), as well as Redfield theorySegal and Nitzan (2005a), and a hybrid NEGF-Redfield approachVelizhanin, Wang, and Thoss (2010). All of these approaches capture the qualitative turnover behaviour of the current quite well in this case, with the MFT and NEGF approaches showing almost quantitative agreement across the full parameter range.

Next, in Fig. 4, the NESS heat current is shown as a function of the system-bath coupling strength, which is related to the vibrational reorganization energy of the molecular junction, $E_{R}=\frac{\omega_{c}\xi}{2}$. Both the Redfield and NEGF-based theories fail to capture the correct turnover behaviour for the heat current, while the MFT results nicely reproduce the exact trend. The MFT result consistently underestimates the total heat current through the intermediate coupling regimes, but the qualitative agreement with the ML-MCTDH benchmark from weak to strong coupling is impressive, and overall is a significant improvement over the other approaches.

A strongly nonadiabatic regime is depicted in Figure 5, with $\omega_{c}/\Delta=10$. In this regime the MFT approach is compared with the non-interacting blip approximation (NIBA) and it’s extension (xNIBA), which is expected to be more accurateNicolin and Segal (2011). Again we see very good agreement between the theories, across the entire parameter regime, with MFT and xNIBA showing a very similar shape, especially in the strong coupling regime.

In Fig. 6, the NESS heat current is investigated in another rather challenging scenario; a low-temperature system with a sizeable temperature gradient, and strong vibrational coupling, tuning between weak to strong system-bath coupling regimes. Overall, we see that the MFT approach emerges as a rather accurate approximate approach. Figure 7 shows the behaviour of the heat current through the junction as a function of the average temperature of the the two baths, at a fixed temperature bias. The turnover displayed here mirrors the trends shown in the other parameters, and is representative of the underlying resonant nature of the nonequilibrium transport process.

IV Conclusions and Outlook

Obviously not a panacea, this version of the mean field approach is expected to break down in some cases, for example at very low temperatures, or for very fast or strongly correlated baths. Furthermore, when $\epsilon$ becomes nonzero, MFT may also suffer from the breakdown of detailed balance, analogous to it’s performance for subsystem properties in the traditional, single bath, spin-boson model. In this respect, more accurate trajectory-based techniques based on mean field theoryMin, Agostini, and Gross (2015); Sato, Kelly, and Rubio (2018), or others stemming from the quantum-classical Liouville equationKapral and Ciccotti (1999), or the partially linearized path integral formalismDunkel, Bonella, and Coker (2008), may prove useful.

One technical item of note is that relatively large trajectory ensembles were needed to converge the average heat current, as compared to subsystem properties. Indeed, recent MFT studies of spontaneous emissionHoffmann et al. (2019) show very similar convergence behaviour for subsystem versus bath properties to that observed here. Potentially, other forms for the current operator could be considered, which may have different statistical convergence properties. For example, adopting a form which jointly depends on the junction and bath degrees of freedom may be useful in this respect. To further improve the efficiency, and potentially the accuracy, of this description the mean-field generalized quantum master equation approach, which has previously proven successful in treating the dynamics of systems coupled to harmonic Kelly, Brackbill, and Markland (2015); Montoya-Castillo and Reichman (2016); Kelly et al. (2016); Montoya-Castillo and Reichman (2017) or fully atomistic environmentsPfalzgraff, Kelly, and T.E. (2015), could be explored in this more general setting using an extensions of the Shi-Geva formalismShi and Geva (2003); Zhang, Ka, and Geva (2006); Mulvihill et al. (2019) for non-system operatorsCohen, Wilner, and Rabani (2013).

The multi-trajectory implementation of Ehrenfest mean field theory thus presents surprisingly accurate, and hence rather promising, tool for the description of nonequilibrium thermal energy transport in molecular junctions. This approach lays the foundation for future trajectory-based studies of nonequilibrium steady state transport, such as current fluctuations, higher order counting statistics, and thermal conductivity, in both impurity models for molecular junctions and more realistic systems.

Acknowledgements.

AK acknowledges financial support from the National Sciences and Engineering Research Council (NSERC) of Canada through the Discovery Grant program, and would like to thank Dvira Segal and Gabriel Hanna for helpful comments and insightful discussions.