Shape, orientation and magnitude of the curl quantum flux, the coherence and the statistical correlations in energy transport at nonequilibrium steady state

We consider the molecular vibrations (i.e., C=O stretching) described by two quantum-mechanically coupled oscillators with different frequency ${\omega }_{1}$ and ${\omega }_{2}$ which are immersed into the solvent environment at the interface. The free Hamiltonian for the system and solvent environment in terms of the displacements of oscillations reads

$$\begin{eqnarray}{H}_{{ssf}} & = & \displaystyle \sum _{j=1}^{2}(\displaystyle \frac{{p}_{j}^{2}}{2{m}_{j}}+\displaystyle \frac{1}{2}{m}_{j}{\omega }_{j}^{2}{x}_{j}^{2})+\kappa (\sqrt{{m}_{1}{m}_{2}}{\omega }_{1}{\omega }_{2}{x}_{1}{x}_{2}+\displaystyle \frac{1}{\sqrt{{m}_{1}{m}_{2}}}{p}_{1}{p}_{2})\\ & & +\displaystyle \sum _{{\bf{f}},\sigma }(\displaystyle \frac{{\tilde{p}}_{{\bf{f}}\sigma }^{2}}{2{\tilde{m}}_{{\bf{f}}}}+\displaystyle \frac{1}{2}{\tilde{m}}_{{\bf{f}}}{\omega }_{{\bf{f}}\sigma }^{2}{\tilde{x}}_{{\bf{f}}\sigma }^{2}),\end{eqnarray} \tag{ 1 }$$

where the low-energy fluctuation solvent environment is treated as a set of harmonic oscillators, neglecting the anharmonic effect hereafter. Actually this effect becomes important under some conditions which goes beyond the scope of this article and we would address this issue in the forthcoming studies. m_j and ${\tilde{m}}_{{\bf{k}}}$ are the effective masses of the molecular-vibrations and environmental vibration mode ${\bf{k}},$ respectively; p_j and x_j are the canonical momentum and coordinate (describing the displacement of the oscillation), respectively. κ is dimensionless that characterizes the coupling strength between two oscillators. In equation (1) we only pick up the minimal coupling between oscillation modes, in order to elucidate our model. It is in general replaced by the dipole–dipole interaction in many chemical systems. The interaction between the system and solvent environment is in similar form as that between the vibrations

$$\begin{eqnarray}&&{H}_{\mathrm{int}}^{{ss}}=\displaystyle \sum _{i=1}^{2}\displaystyle \sum _{{\bf{f}},\sigma }{\lambda }_{{\bf{f}}\sigma }^{{ss}}(\sqrt{{m}_{i}{\tilde{m}}_{{\bf{f}}}}{\omega }_{i}{\omega }_{{\bf{f}}\sigma }{x}_{i}{\tilde{x}}_{{\bf{f}}\sigma }+\displaystyle \frac{1}{\sqrt{{m}_{i}{\tilde{m}}_{{\bf{f}}}}}{p}_{i}{\tilde{p}}_{{\bf{f}}\sigma }).\end{eqnarray} \tag{ 2 }$$

Additionally one oscillator interacts with a heat source characterizing the energy pump from the chemical reactions, and the other one interacts with a cooler environment which harvests the energy transfered by molecular vibrations, as shown in figure 1. The quantum-mechanically effective couplings of system to these reservoirs can be realized by exchanging the energy quanta. Generally the interaction between molecule stretchings is mediated via the dipole–dipole coupling, therefore in the formalism of quantum field theory, the free Hamiltonian of system+reservoirs+environment is

$$\begin{eqnarray}&&{H}_{0}={\bar{\varepsilon }}_{1}{a}_{1}^{\dagger }{a}_{1}+{\bar{\varepsilon }}_{2}{a}_{2}^{\dagger }{a}_{2}+\bar{{\rm{\Delta }}}({a}_{1}^{\dagger }{a}_{2}+{a}_{2}^{\dagger }{a}_{1})+\displaystyle \sum _{\nu =1}^{3}\;\displaystyle \sum _{{\bf{k}},\sigma }{\hslash }{\omega }_{{\bf{k}}\sigma }{b}_{{\bf{k}}\sigma }^{(\nu ),\dagger }{b}_{{\bf{k}}\sigma }^{(\nu )}\end{eqnarray} \tag{ 3 }$$

and the interactions are

$$\begin{eqnarray}&&{H}_{\mathrm{int}}=\displaystyle \sum _{i=1}^{2}\;\displaystyle \sum _{{\bf{k}},p\;}\;{g}_{{\bf{k}}p}({a}_{i}^{\dagger }{b}_{{\bf{k}}p}^{(i)}+{a}_{i}{b}_{{\bf{k}}p}^{(i),\dagger })+\displaystyle \sum _{{\bf{q}}\;,\sigma }\;{g}_{{\bf{q}}\sigma }({c}^{\dagger }{b}_{{\bf{q}}\sigma }^{(3)}+c\ {b}_{{\bf{q}}\sigma }^{(3),\dagger }),\end{eqnarray} \tag{ 4 }$$

where $c\equiv {a}_{1}+{a}_{2}$ and ${b}_{{\bf{k}}\sigma }^{(1)},$ ${b}_{{\bf{k}}\sigma }^{(2)}$ and ${b}_{{\bf{k}}\sigma }^{(3)}$ are the bosonic annihilation operators for heat source, cool reservoir and solvent environment, respectively. σ and p denote the polarizations of bosons in the reservoirs. The rotating-wave approximation [39, 40] has been appiled to the vibration-bath interactions, owing to the dominant contribution by real absorption and emission of quanta in long time limit. The effect of VP interaction can be approximately described by the renormalized frequency gap $\delta \bar{\varepsilon }={\bar{\varepsilon }}_{1}-{\bar{\varepsilon }}_{2}$ and coupling strength $\bar{{\rm{\Delta }}}$ of the molecular vibrations [41]. This in particular, indicates that the tuning between frequencies $\delta \bar{\varepsilon }\ll \bar{{\rm{\Delta }}}$ (strong VP) makes the wave packet of vibrations extended while the large detuning between frequencies $\delta \bar{\varepsilon }\gg \bar{{\rm{\Delta }}}$ (weak VP) makes the wave packet localized, as elaborately illustrated in Anderson localization mechanism [42–44], in which the disorder (both diagonal and off-diagonal) is imaged to be connected with the presence of impurities, vacancies and dislocations in an ideal crystal lattice, or the random distributions of atoms and molecules [45].

Zoom In Zoom Out Reset image size

The dynamics of the system is governed by the reduced quantum master equation (RQME), which is obtained by tracing out the degree of freedoms of the baths. As we pointed out above and also discussed in previous papers [41, 46, 47], the strong interactions between the system and some discrete vibrational modes (i.e., the hydrogen bond) owing to the quasi-resonance between frequencies, leads to the comparable time scales between the vibrational modes and system, which subsequently acquires us to include the dynamics of these discrete vibrational modes together with the system. In other words, these vibrational modes must be seperated from the reservoirs. Thereby the remaining modes consisting of the low-energy fluctuations can be reasonably treated as the baths, which are effectively in weak coupling to the systems due to the mismatch of frequencies between these continous modes and the system. On the basis of perturbation theory, the whole solution of the density operator can be written as ${\rho }_{\mathrm{SR}}={\rho }_{s}(t)\otimes {\rho }_{R}(0)+{\rho }_{\delta }(t)$ with the traceless term in higher orders of couplings between system and reservoirs. Because the time scale associated with the environmental correlations is much smaller than that of system over which the state varies appreciably, the RQME for the reduced density matrix of the systems in the interaction picture can be derived under the so-called Markoff approximation

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}{\rho }_{s}}{{\rm{d}}t} & = & \displaystyle \frac{1}{2{{\hslash }}^{2}}\{\displaystyle \sum _{\nu =1}^{2}\;\displaystyle \sum _{p=1}^{2}[{\gamma }_{p}^{{T}_{\nu },+}({a}_{p}{\rho }_{s}{a}_{\nu }^{\dagger }-{a}_{\nu }^{\dagger }{a}_{p}{\rho }_{s})+{\gamma }_{p}^{{T}_{\nu },-}({a}_{p}^{\dagger }{\rho }_{s}{a}_{\nu }-{a}_{\nu }{a}_{p}^{\dagger }{\rho }_{s})]\\ & & +\displaystyle \sum _{j=1}^{2}\;\displaystyle \sum _{p=1}^{2}[{\gamma }_{p}^{{T}_{3},+}({a}_{p}{\rho }_{s}{a}_{j}^{\dagger }-{a}_{j}^{\dagger }{a}_{p}{\rho }_{s})+{\gamma }_{p}^{{T}_{3},-}({a}_{p}^{\dagger }{\rho }_{s}{a}_{j}-{a}_{j}{a}_{p}^{\dagger }{\rho }_{s})]\}+{\rm{h.c.}},\end{eqnarray} \tag{ 5 }$$

where the reservoirs and solvent environment are in thermal equilibrium. The dissipation rates ${\gamma }_{...}$ are given in appendix A. ${n}_{\omega }^{T}={[{\rm{exp}}({\hslash }\omega /{k}_{{\rm{B}}}T)-1]}^{-1}$ is the Bose occupation on frequency ω at temperature T. The mixture angle reads

$$\begin{eqnarray}&&{\rm{cos}}2\theta =\displaystyle \frac{{\bar{\varepsilon }}_{2}-{\bar{\varepsilon }}_{1}}{\sqrt{{({\bar{\varepsilon }}_{1}-{\bar{\varepsilon }}_{2})}^{2}+4{\bar{{\rm{\Delta }}}}^{2}}};\quad {\rm{sin}}2\theta =-\displaystyle \frac{2\bar{{\rm{\Delta }}}}{\sqrt{{({\bar{\varepsilon }}_{1}-{\bar{\varepsilon }}_{2})}^{2}+4{\bar{{\rm{\Delta }}}}^{2}}}.\end{eqnarray} \tag{ 6 }$$

Conventionally the QME in equation (5) was solved in Liouville space [48, 49], by writing the density matrix as a supervector. This strategy however, seems to be unrealizable due to the infinite dimension of Fock space for bosons. Here we will solve the QME in the coherent representation, which was first developed by Glauber [50]. It is alternatively named as P-representation, according to the terminology in quantum optics [39, 40, 51]. As is known, the eigenstate of the annihilation operators is $| {\alpha }_{1},{\alpha }_{2}\rangle$ where the eigenequation ${a}_{j}| {\alpha }_{1},{\alpha }_{2}\rangle ={\alpha }_{j}| {\alpha }_{1},{\alpha }_{2}\rangle$ is satisfied. In terms of these components, the density matrix can be expanded into the following form [39]

$$\begin{eqnarray}&&{\rho }_{s}(t)=\int P({\alpha }_{\beta },{\alpha }_{\beta }^{*},t)| {\alpha }_{1},{\alpha }_{2}\rangle \langle {\alpha }_{1},{\alpha }_{2}| {{\rm{d}}}^{2}{\alpha }_{1}{{\rm{d}}}^{2}{\alpha }_{2},\end{eqnarray} \tag{ 7 }$$

where $P({\alpha }_{\beta },{\alpha }_{\beta }^{*},t)$ is called the quasi-probability, due to the overcompleteness of the coherent basis and the non-positive definition of $P({\alpha }_{\beta },{\alpha }_{\beta }^{*},t).$ But this quasi-probability is always non-negative everywhere if the quantum system has classical analog [52]. Using equation (7) we can project the QME into coherent space and then obtain the following dynamical equation for $P({\alpha }_{\beta },{\alpha }_{\beta }^{*},t)$

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}P({\alpha }_{\beta },{\alpha }_{\beta }^{*}) & = & \gamma [2(\displaystyle \frac{\partial }{\partial {\alpha }_{1}}{\alpha }_{1}+\displaystyle \frac{\partial }{\partial {\alpha }_{2}}{\alpha }_{2})+\displaystyle \frac{\partial }{\partial {\alpha }_{1}}{\alpha }_{2}+\displaystyle \frac{\partial }{\partial {\alpha }_{2}}{\alpha }_{1}+{\rm{c.c.}}]P({\alpha }_{\beta },{\alpha }_{\beta }^{*})\\ & & +\gamma [2{Y}_{1}^{1}\displaystyle \frac{{\partial }^{2}}{\partial {\alpha }_{1}^{*}\partial {\alpha }_{1}}+2{Y}_{2}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {\alpha }_{2}^{*}\partial {\alpha }_{2}}\\ & & +{Y}_{12}^{21}(\displaystyle \frac{{\partial }^{2}}{\partial {\alpha }_{1}^{*}\partial {\alpha }_{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {\alpha }_{1}\partial {\alpha }_{2}^{*}})]P({\alpha }_{\beta },{\alpha }_{\beta }^{*}),\end{eqnarray} \tag{ 8 }$$

where the coefficients ${Y}_{...}$ are given in appendix A. Equation (8) is of the same formalism as the classical Fokker–Planck equation [53], apart from the complex variables. Moreover, this equation describes the Ornstein–Uhlenbeck process in the phase space, so that it is exactly solvable. Now we focus on the steady-state case with $t\to \infty ,$ and therefore the steady-state solution is of Guassian type

$$\begin{eqnarray}&&{P}_{{ss}}({\alpha }_{\beta },{\alpha }_{\beta }^{*})=\displaystyle \frac{1}{Z}{{\rm{e}}}^{-[a{| {\alpha }_{1}| }^{2}+b{| {\alpha }_{2}| }^{2}+2c{\rm{Re}}({\alpha }_{1}^{*}{\alpha }_{2})]}\end{eqnarray} \tag{ 9 }$$

with the a, b and c being

$$\begin{eqnarray}a & = & \displaystyle \frac{4({Y}_{1}^{1}+7{Y}_{2}^{2}-2{Y}_{12}^{21})}{{({Y}_{1}^{1}+{Y}_{2}^{2})}^{2}+4[3{Y}_{1}^{1}{Y}_{2}^{2}-{({Y}_{12}^{21})}^{2}]},b=\displaystyle \frac{4(7{Y}_{1}^{1}+{Y}_{2}^{2}-2{Y}_{12}^{21})}{{({Y}_{1}^{1}+{Y}_{2}^{2})}^{2}+4[3{Y}_{1}^{1}{Y}_{2}^{2}-{({Y}_{12}^{21})}^{2}]}\\ c & = & \displaystyle \frac{8({Y}_{1}^{1}+{Y}_{2}^{2}-2{Y}_{12}^{21})}{{({Y}_{1}^{1}+{Y}_{2}^{2})}^{2}+4[3{Y}_{1}^{1}{Y}_{2}^{2}-{({Y}_{12}^{21})}^{2}]},Z=\displaystyle \frac{{\pi }^{2}}{12}\{{({Y}_{1}^{1}+{Y}_{2}^{2})}^{2}+4[3{Y}_{1}^{1}{Y}_{2}^{2}-{({Y}_{12}^{21})}^{2}]\}.\end{eqnarray} \tag{ 10 }$$

In the forthcoming sections, we will discuss the nonequilibrium behaviors and heat transport based on this steady-state solution of quasi-probability in equations (9) and (10).

The quantitative description of nonequilibriumness, especially the far-from-equilibrium case in quantum systems was devoid, although it is known that macroscopical current would be observed if the system deviates from equilibrium. Recently the curl flux in discrete space was developed to describe the nonequilibrium behaviors in quantum systems [37, 38]. Here we will alternatively develope a curl quantum flux in the continous space, in analogous with the classical case [35]. We are able to write the probabilistic evolution of diffusion equation (8) covering the whole coherent space into the form of ${\partial }_{t}P+{{\rm{\nabla }}}_{\alpha }\cdot {\bf{J}}=0$ which gives the steady-state case ${{\rm{\nabla }}}_{\alpha }\cdot {\bf{J}}=0.$ For the nonequilibrium systems at steady state in general, this does not necessarily mean that the flux ${\bf{J}}$ has to vanish, due to the detailed-balance-breaking. Instead, the divergence-free nature implies that the flux is a rotational curl field in coherent space. To further elucidate the nonequilibrium nature of the curl quantum flux, we will reduce our discussion into the space spanned by (${x}_{1},{x}_{2}$) where ${x}_{j}={\rm{Re}}[{\alpha }_{j}].$ On the other hand, the issue of even and odd variables then does not arise, because of the time-reversibility of x_j.

Assuming ${\alpha }_{j}={x}_{j}+{\rm{i}}{p}_{j}$ where ${x}_{j}=\sqrt{2{m}_{j}{\omega }_{j}/{\hslash }}\langle {\alpha }_{1},{\alpha }_{2}| {\hat{x}}_{j}| {\alpha }_{1},{\alpha }_{2}\rangle$ are the dimensionless means of the displacements of the oscillators in coherent space, the dynamical equation of probability in (${x}_{1},{x}_{2}$) space becomes

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}P({x}_{1},{x}_{2}) & = & \gamma (2\displaystyle \frac{\partial }{\partial {x}_{1}}{x}_{1}+2\displaystyle \frac{\partial }{\partial {x}_{2}}{x}_{2}+\displaystyle \frac{\partial }{\partial {x}_{1}}{x}_{2}+\displaystyle \frac{\partial }{\partial {x}_{2}}{x}_{1})P({x}_{1},{x}_{2})\\ & & +\displaystyle \frac{\gamma }{2}({Y}_{1}^{1}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{1}^{2}}+{Y}_{2}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{2}^{2}}+{Y}_{12}^{21}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{1}\partial {x}_{2}})P({x}_{1},{x}_{2})\end{eqnarray} \tag{ 11 }$$

the steady-state solution to which reads

$$\begin{eqnarray}&&{P}_{{ss}}({x}_{1},{x}_{2})=\displaystyle \frac{\sqrt{{ab}-{c}^{2}}}{\pi }{{\rm{e}}}^{-({{ax}}_{1}^{2}+{{bx}}_{2}^{2}+2{{cx}}_{1}{x}_{2})}\end{eqnarray} \tag{ 12 }$$

Thus the curl quantum flux is of the form

$$\begin{eqnarray}{{\mathcal{J}}}_{{x}_{1}} & = & \displaystyle \frac{\gamma ({Y}_{1}^{1}-{Y}_{2}^{2}){P}_{{ss}}}{{({Y}_{1}^{1}+{Y}_{2}^{2})}^{2}+4[3{Y}_{1}^{1}{Y}_{2}^{2}-{({Y}_{12}^{21})}^{2}]}\\ & & \times [2({Y}_{1}^{1}+{Y}_{2}^{2}-2{Y}_{12}^{21}){x}_{1}+(7{Y}_{1}^{1}+{Y}_{2}^{2}-2{Y}_{12}^{21}){x}_{2}]\\ {{\mathcal{J}}}_{{x}_{2}} & = & -\displaystyle \frac{\gamma ({Y}_{1}^{1}-{Y}_{2}^{2}){P}_{{ss}}}{{({Y}_{1}^{1}+{Y}_{2}^{2})}^{2}+4[3{Y}_{1}^{1}{Y}_{2}^{2}-{({Y}_{12}^{21})}^{2}]}\\ & & \times [({Y}_{1}^{1}+7{Y}_{2}^{2}-2{Y}_{12}^{21}){x}_{1}+2({Y}_{1}^{1}+{Y}_{2}^{2}-2{Y}_{12}^{21}){x}_{2}].\end{eqnarray} \tag{ 13 }$$

The behaviors of curl quantum flux are illustrated in figure 3, which provides a description of nonequilibriumness on microscopic level. It is straightforward to prove that the stream lines of the curl flux in 2D space in our case form a set of ellipses, with the major axis in the vicinity of the anti-diagonal line. Therefore the polarization of the curl flux can be quantified in terms of geometric language, namely the eccentricity $\bar{e}$ and rotation angle $\beta /2$

$$\begin{eqnarray}&&\bar{e}=\sqrt{\displaystyle \frac{2\sqrt{{(a-b)}^{2}+4{c}^{2}}}{a+b+\sqrt{{(a-b)}^{2}+4{c}^{2}}}},\quad {\rm{sin}}\beta =-\displaystyle \frac{2c}{\sqrt{{(a-b)}^{2}+4{c}^{2}}}.\end{eqnarray} \tag{ 14 }$$

In particular, based on equations (12) and (13) we know that the shape and orientation of curl flux is governed by the eccentricity $\bar{e}$ and rotation angle $\beta /2$ where the polarization is in slender-cigar shape along the vicinity of the line $| {x}_{1}| =| {x}_{2}|$ as $\bar{e}$ increases and $\beta /2$ aprroaches $\pi /4.$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 3(a) and (b) show the effect of VP coupling on the curl flux with the eccentricities $\bar{e}=0.978$ and $\bar{e}=0.957,$ respectively. As we can see, the strong-VP-bond-contributed delocalization of the excited vibrational modes causes the flux to be more polarized in the vicinity of anti-diagonal than the weak-VP-bond-contributed

localization does, which as explored later, means that the correlation between the molecular vibrations is much stronger as quasi-particles become delocalized, rather than the localization of quasi-particles. The quality of heat transport will be promoted as the correlation between vibrations becomes strong, as shown in the following discussion. Figures 3(c) and (d) show that the thermal fluctuations in heat source which in some sense dictates the effective thermal voltage, can definitely strengthen the molecular-vibration correlation, which will considerably raise the heat transport, as will be shown in figure 4(c). Figures 3(e) and (f) illustrate the effect of solvent environment on the curl flux, which shows that the correlation between molecular vibrations is unavoidably suppressed by the thermal fluctuations induced by solvent environment.

Zoom In Zoom Out Reset image size

Due to the vector feature of the curl flux, here we will use the Wilson loop to quantify the value of the curl flux, based on the integral of curl flux along a specific closed path

$$\begin{eqnarray}&&W=\displaystyle \frac{1}{L}{\displaystyle \int }_{{\rm{\Sigma }}}{J}_{{x}_{1}}{\rm{d}}{x}_{1}+{J}_{{x}_{2}}{\rm{d}}{x}_{2},\end{eqnarray} \tag{ 15 }$$

where L represents the length of closed path Σ. Currently we choose the closed path as one of the stream lines of curl flux in 2D space, to maximize the quantity W in equation (15). After some manipulations the Wilson loop of curl flux reads

$$\begin{eqnarray}&&W=\displaystyle \frac{\gamma | {Y}_{1}^{1}-{Y}_{2}^{2}| (a+b)}{16E(\displaystyle \frac{\pi }{2}| {\bar{e}}^{2})\sqrt{e}}\sqrt{\displaystyle \frac{{ab}-{c}^{2}}{a+b+\sqrt{{(a-b)}^{2}+4{c}^{2}}}}.\end{eqnarray} \tag{ 16 }$$

Notice that e = 2.71828... is the Euler's number and $E(\phi | {k}^{2})$ is the elliptic integral of the second kind. As is shown in figure 2(g), the detailed-balance is more broken as the two molecules approaches each other, governed by the increase of $\frac{\bar{{\rm{\Delta }}}}{\delta \bar{\varepsilon }}.$ This can be understood by an extreme case that the two molecular vibrations will equilibrate individually with the environments as they becomes infinitely distanced with no correlation between each other ($\bar{{\rm{\Delta }}}=0$).

The correlations between the molecular vibrations are measured by the correlation functions. In particular, we mainly focus on the four-point correlation functions here, which corresponds to the density–density correlations

$$\begin{eqnarray}&&{C}_{4}=\displaystyle \frac{\langle {a}_{1}^{\dagger }{a}_{1}{a}_{2}^{\dagger }{a}_{2}\rangle }{\langle {a}_{1}^{\dagger }{a}_{1}\rangle \langle {a}_{2}^{\dagger }{a}_{2}\rangle }-1=\displaystyle \frac{{\bar{e}}^{4}{{\rm{sin}}}^{2}\beta }{4(1-{\bar{e}}^{2})+{\bar{e}}^{4}{{\rm{sin}}}^{2}\beta },\end{eqnarray} \tag{ 17 }$$

where $\bar{e}$ shares the same definition of the eccentricity as before, and $\beta /2$ is the angle between the major axis of the ellipse and x-axis. Equation (17) follows the definition of density–density correlation function in the site basis, as given in [39]. ${C}_{4}=0$ indicates that the occupations on the two vibrational modes are uncorrelated. Equation (17) in fact uncovers the connection between the microscopic nonequilibriumness (flux) and the macroscopic observables in a geometric manner. Hence it is evident to say that the correlations between molecular vibrations are charaterized by the polarization of the curl flux in coherent space, where the slender-cigar type of polarization of flux along the anti-diagonal means the strong vibrational correlations while the cake type of polarization with tiny inhomogeneity or polarization along the axis means the weak vibrational correlations. Possibly the microscopic curl quantum flux can be explored in experiments, by measuring these geometric parameters through the measurement of density–density correlation function, which has been recently probed in degenerate quantum gas [54, 55].

The large figure in figure 2(e) show that (i) the thermal fluctuations in the heat source do strengthen the vibration correlations and (ii) the detuning between the frequencies of molecular vibrations distroys the vibration correlations, as long as the heat transport is on track (this will be examined in detail in the section of heat transport). These properties are also microscopically reflected by curl flux, as shown in figures 2(a), (b) and the paths (dashed lines) on the landscape of density–density correlation C₄ in the small figure in figure 2(e), in spite of the reduction of the contribution by rotation angle. Besides the geometry of the curl flux, the promotion of the magnitude of curl flux (quantified by Wilson loop) is also strongly correlated to the increase of correlation function, by improving the coupling strength between vibrational modes, as shown in figure 2(h).

As will show later in the section of coherence effect, the site-basis coherence has no contribution to the heat transport in the secular approximation since it is decoupled from the population dynamics. As approaching this regime, one can easily demonstrate that ${\mathrm{lim}}_{\epsilon \to 0\;}\tilde{c}=0$ in equation (26), which subsequently gives $\beta \to 0.$ Hence ${\mathrm{lim}}_{\epsilon \to 0}{C}_{4}=0,$ which is physically reasonable owing to the polarization of flux at the moment orientatied in the vicinity of x₁-axis. As the entanglement between coherence and population dynamics adiabatically increases, the macriscopic correlation C₄ is considerably promoted by coherence from the microscopic curl flux as it becomes significantly polarized with the angle approaching $-\frac{\pi }{4},$ as illustrated in figures 2(c), (d) and (f). Therefore we can conclude that the coherence generates and considerably improves the correlation between the molecular vibrations.

We first introduce the heat-current operators J₁, J₂ for the heat currents pumping into and flowing out from the system, respectively. By ignoring the back influence of system to reservoirs, the current operators are defined as ${J}_{\nu }=\frac{1}{{\rm{i}}{\hslash }}[{H}_{s},{H}_{\mathrm{int}}^{(\nu )}],$ and furthermore in our system

$$\begin{eqnarray}{J}_{1} & = & \displaystyle \frac{1}{{\rm{i}}{\hslash }}\displaystyle \sum _{{\bf{k}},\sigma \;}{g}_{{\bf{k}}\sigma }[({\bar{\varepsilon }}_{1}{a}_{1}^{\dagger }+\bar{{\rm{\Delta }}}{a}_{2}^{\dagger }){b}_{{\bf{k}}\sigma }^{(1)}-{\rm{h.c.}}],\\ {J}_{2} & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \sum _{{\bf{k}},\sigma \;}{g}_{{\bf{k}}\sigma }[({\bar{\varepsilon }}_{2}{a}_{2}^{\dagger }+\bar{{\rm{\Delta }}}{a}_{1}^{\dagger }){b}_{{\bf{k}}\sigma }^{(2)}-{\rm{h.c.}}].\end{eqnarray} \tag{ 18 }$$

To calculate the heat current up to the 2nd order of coupling strength, we need to carry out the 1st order correction to the density matrix so that ${\rho }_{s}(t)={\rho }_{s}(0)+\frac{{\rm{i}}}{{\hslash }}{\displaystyle \int }_{0}^{t}[{\rho }_{s}(t),{\tilde{H}}_{\mathrm{int}}(\tau )]{\rm{d}}\tau .$ Substituting this into the product ${\rho }_{s}(t){\tilde{J}}_{1}(t)$ and after a lengthy derivation shown in SM the subsequent heat current pumping into the molecules is

$$\begin{eqnarray}{\langle {J}_{1}\rangle }_{{ss}} & = & \underset{t\to \infty }{\mathrm{lim}}{\rm{Tr}}[{\rho }_{s}(t){\tilde{J}}_{1}(t)]\\ & = & (-\gamma )[2{\bar{\varepsilon }}_{1}\langle {a}_{1}^{\dagger }{a}_{1}\rangle +\bar{{\rm{\Delta }}}\langle {a}_{1}^{\dagger }{a}_{2}+{a}_{2}^{\dagger }{a}_{1}\rangle -2({E}_{1}{n}_{{\nu }_{1}}^{{T}_{1}}{{\rm{cos}}}^{2}\theta +{E}_{2}{n}_{{\nu }_{2}}^{{T}_{1}}{{\rm{sin}}}^{2}\theta )]\\ & = & -\displaystyle \frac{\gamma }{6}\{[{E}_{1}(7{{\rm{cos}}}^{2}\theta -2{\rm{sin}}\theta {\rm{cos}}\theta )+{E}_{2}(7{{\rm{sin}}}^{2}\theta +2{\rm{sin}}\theta {\rm{cos}}\theta )]{Y}_{1}^{1}\\ & & +[{E}_{1}({{\rm{cos}}}^{2}\theta -2{\rm{sin}}{\rm{cos}}\theta )+{E}_{2}({{\rm{sin}}}^{2}\theta +2{\rm{sin}}{\rm{cos}}\theta )]({Y}_{2}^{2}-2{Y}_{12}^{21})\\ & & -12({E}_{1}{n}_{{\nu }_{1}}^{{T}_{1}}{{\rm{cos}}}^{2}\theta +{E}_{2}{n}_{{\nu }_{2}}^{{T}_{2}}{{\rm{sin}}}^{2}\theta )\},\end{eqnarray} \tag{ 19 }$$

where ${\hslash }{\nu }_{1}=\frac{1}{2}[{\bar{\varepsilon }}_{1}+{\bar{\varepsilon }}_{2}-\sqrt{{({\bar{\varepsilon }}_{1}-{\bar{\varepsilon }}_{2})}^{2}+4{\bar{{\rm{\Delta }}}}^{2}}]$ and ${\hslash }{\nu }_{2}=\frac{1}{2}[{\bar{\varepsilon }}_{1}+{\bar{\varepsilon }}_{2}+\sqrt{{({\bar{\varepsilon }}_{1}-{\bar{\varepsilon }}_{2})}^{2}+4{\bar{{\rm{\Delta }}}}^{2}}]$ are the eigenenergies of the coupled oscillators. Moreover the coherence contribution to heat currents ${J}_{1},{J}_{2}$ is governed by the term $\langle {a}_{1}^{\dagger }{a}_{2}+{a}_{2}^{\dagger }{a}_{1}\rangle$ which will vanishes in secular approximation as shown later. $\langle {a}_{1}^{\dagger }{a}_{2}\rangle ={\displaystyle \sum }_{{n}_{1},{n}_{2}}\sqrt{{n}_{1}{n}_{2}}\ \langle {n}_{1}-1,{n}_{2}| {\rho }_{s}| {n}_{1},{n}_{2}-1\rangle .$ The similar manner for calculating the heat current tranfered by the molecular vibrations and we can simply replace the ${\tilde{H}}_{\mathrm{int}}^{(1)}$ by ${\tilde{H}}_{\mathrm{int}}^{(2)}$ in J₁, to reach the following result

$$\begin{eqnarray}{\langle {J}_{2}\rangle }_{{ss}} & = & \displaystyle \frac{\gamma }{6}\{[{E}_{1}({{\rm{sin}}}^{2}\theta -2{\rm{sin}}\theta {\rm{cos}}\theta )+{E}_{2}({{\rm{cos}}}^{2}\theta +2{\rm{sin}}\theta {\rm{cos}}\theta )]({Y}_{1}^{1}-2{Y}_{12}^{21})\\ & & +[{E}_{1}(7{{\rm{sin}}}^{2}\theta -2{\rm{sin}}\theta {\rm{cos}}\theta )+{E}_{2}(7{{\rm{cos}}}^{2}\theta +2{\rm{sin}}\theta {\rm{cos}}\theta )]{Y}_{2}^{2}\\ & & -12({E}_{1}{n}_{{\nu }_{1}}^{{T}_{2}}{{\rm{sin}}}^{2}\theta +{E}_{2}{n}_{{\nu }_{2}}^{{T}_{2}}{{\rm{cos}}}^{2}\theta )\}\end{eqnarray} \tag{ 20 }$$

As a QHE, the stationary working efficiency of the vibratioanl energy transport is naturally defined as $\eta =\frac{{\langle {J}_{2}\rangle }_{{ss}}}{{\langle {J}_{1}\rangle }_{{ss}}}.$ Figure 4(a) shows the effect of external pumping on the heat flowing into the system according to equation (19), which demonstrates the improvement of energy pumping into system by the activity of external heat source. Since the excitation energy in our system is around 1 eV, the considerable excitation with respect to this energy scale needs the effective temperature of the environment to be around 5000 K. On the other hand, the analog is the light-harvesting complex in which the temperature of radiations is around 5800 K. Figures 4(c) and (e) illustrate the heat flow (quantifying the vibrational energy transfer) and working efficiency η under the influence of external energy pumping by heat source, according to equations (19) and (20). As seen first, the thermal fluctuation and pumping of the external heat source causes a considerable improvement of the energy transfered by the molecular vibrations and the working efficiency reflected in figure 4(e) as well. More importantly, it is also shown in figure 4(c) that the energy transport relates to a critical value of ($\delta \bar{\varepsilon },{T}_{1}$), under which the vibrational energy transport is suspended. The critical points are determined by acquiring ${\langle {J}_{2}\rangle }_{{ss}}=0$ which will be shown in Supplementary Information (SI). In figure 4(c) the critical values of T₁ are 3372, 2514 and 2257 K with respect to $\delta \bar{\varepsilon }=0.3,0.1,0.01$ eV, respectively. Above the critical point, the promotion of energy current flowing through molecules and the efficiency by thermal fluctuation of the heat source and also the frequency detuning between molecular vibrations can be critically demonstrated by the polarization of the curl flux, illustreated in figures 3(a)–(d), respectively, since the correlation between the molecular vibrations is enhanced which will be reached later. Therefore we can evidently claim that the curl quantum flux in equations (12) and (13) on microscopic level significantly correlates to and characterizes the vibrational energy transport on macroscopic level. We will come back to this issue when discussing the correlation functions later.

To explore the statistical distribution of heat current, one dose not only necessarily calculate the mean heat current as what we did above, but also need to uncover the higher order properties, i.e., the correlations between the currents. By replacing the summation over different modes in reservoir by the integration

$$\begin{eqnarray}&&\displaystyle \sum _{{\bf{k}},\sigma \;}{g}_{{\bf{k}}\sigma }^{2}\langle {b}_{{\bf{k}}\sigma }^{(2),\dagger }{b}_{{\bf{k}}\sigma }^{(2)}\rangle {{\rm{e}}}^{{\rm{i}}(\nu -{\omega }_{{\bf{k}}\sigma }){\rm{\Delta }}t}\\ &&\quad \longrightarrow \;\gamma {\rm{sin}}(\nu {\rm{\Delta }}t){\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{\rm{sin}}({\omega }_{{\bf{k}}\sigma }{\rm{\Delta }}t)}{{{\rm{e}}}^{{\hslash }{\omega }_{{\bf{k}}\sigma }/{kT}}-1}{\rm{d}}{\omega }_{{\bf{k}}\sigma }=\displaystyle \frac{\pi \gamma }{2}{\rm{sin}}(\nu {\rm{\Delta }}t)({\rm{coth}}\displaystyle \frac{\pi {\rm{\Delta }}t}{\beta {\hslash }}-\displaystyle \frac{\beta {\hslash }}{\pi {\rm{\Delta }}t})\end{eqnarray} \tag{ 21 }$$

and after a lengthy calculation we arrive at the current–current correlation function at steady state

$$\begin{eqnarray}&&{\rm{Re}}[{\langle {\tilde{J}}_{2}(t){\tilde{J}}_{2}(t^{\prime} )\rangle }_{{ss}}]\\ &&\quad =\displaystyle \frac{\pi \gamma }{{{\hslash }}^{2}}\{({E}_{1}^{2}\langle {\eta }_{1}^{\dagger }{\eta }_{1}\rangle {{\rm{sin}}}^{2}\theta +{E}_{2}^{2}\langle {\eta }_{2}^{\dagger }{\eta }_{2}\rangle {{\rm{cos}}}^{2}\theta +{E}_{1}{E}_{2}\langle {\eta }_{1}^{\dagger }{\eta }_{2}+{\eta }_{2}^{\dagger }{\eta }_{1}\rangle {\rm{sin}}\theta {\rm{cos}}\theta )\delta ({\rm{\Delta }}t)\\ &&\qquad +\displaystyle \frac{1}{2}[{E}_{1}^{2}{\rm{sin}}({\nu }_{1}{\rm{\Delta }}t)\langle 2{\eta }_{1}^{\dagger }{\eta }_{1}+1\rangle {{\rm{sin}}}^{2}\theta +{E}_{2}^{2}{\rm{sin}}({\nu }_{2}{\rm{\Delta }}t)\langle 2{\eta }_{2}^{\dagger }{\eta }_{2}+1\rangle {{\rm{cos}}}^{2}\theta \\ &&\qquad +{E}_{1}{E}_{2}({\rm{sin}}({\nu }_{1}{\rm{\Delta }}t)+{\rm{sin}}({\nu }_{2}{\rm{\Delta }}t))\langle {\eta }_{1}^{\dagger }{\eta }_{2}+{\eta }_{2}^{\dagger }{\eta }_{1}\rangle {\rm{sin}}\theta {\rm{cos}}\theta ]({\rm{coth}}\displaystyle \frac{\pi {\rm{\Delta }}t}{{\beta }_{1}{\hslash }}-\displaystyle \frac{{\beta }_{1}{\hslash }}{\pi {\rm{\Delta }}t})\},\end{eqnarray} \tag{ 22 }$$

where ${\rm{\Delta }}t\equiv t-t^{\prime} ,{\beta }_{i}\equiv \frac{1}{{{kT}}_{i}}.$ Equation (22) roughly illustrates the beat oscillation will occur if the detuning between the frequencies of vibrations is suppressed. To demonstrate this, in figure 5 we show both the cases $\delta \bar{\varepsilon }\ll \bar{{\rm{\Delta }}}$ and $\delta \bar{\varepsilon }\gg \bar{{\rm{\Delta }}}$ in current–current correlation function (apart from the crusp peak described by the δ-function). The former one characterizing the strong vibrational coupling of molecular chain to the molecule stretching (i.e., C=O stretching) reveals the beat oscillation and subsequently a coherent regime of correlation, in which the correlation is strong. This is in contrast to the latter one characterizing the weak vibrational coupling of molecular chain to the molecule stretching, with an incoherent way of correlation in which the correlation is much weaker.

Zoom In Zoom Out Reset image size