Heat transport along a chain of coupled quantum harmonic oscillators

TL;DR

This paper investigates heat transport in a chain of coupled quantum harmonic oscillators connected to heat reservoirs, deriving a closed-form expression for thermal conductance and revealing infinite thermal conductivity in the thermodynamic limit.

Contribution

It introduces a quantum approach based on a density operator evolution equation to analyze heat transport and derives a closed-form expression for conductance.

Findings

Thermal conductance is finite for small interactions.

Thermal conductivity diverges in the thermodynamic limit.

The approach uses covariances and stationary states.

Abstract

We study the heat transport properties of a chain of coupled quantum harmonic oscillators in contact at its ends with two heat reservoirs at distinct temperatures. Our approach is based on the use of an evolution equation for the density operator which is a canonical quantization of the classical Fokker-Planck-Kramers equation. We set up the evolution equation for the covariances and obtain the stationary covariances at the stationary states from which we determine the thermal conductance in closed form when the interparticle interaction is small. The conductance is finite in the thermodynamic limit implying an infinite thermal conductivity.

Equations

$$\mathcal{H}=\frac{1}{2m}\sum _{i=1}^L{p}_i^2+\frac{k_0}{2}\sum _{i=1}^L{q}_i^2-k_1\sum _{i=1}^{L-1}{q}_i{q}_{i+1},$$

$$i\hslash \frac{\partial \rho }{\partial t}=[\mathcal{H},\rho ]-[q_1,J_1]-[q_L,J_L],$$

$$J_i=-\frac{\gamma }{2}(\rho g_i+g_i^{†}\rho )-\frac{\gamma m}{i\hslash {\beta }_i}[q_i,\rho ],$$

$$g_i=-\frac{m}{i\hslash {\beta }_i}(e^{{\beta }_i\mathcal{H}}{q}_i{e}^{-{\beta }_i\mathcal{H}}-q_i).$$

$$g_i=p_i+\frac{{\beta }_i}{2!}[\mathcal{H},p_i]+\frac{{\beta }_i^2}{3!}[\mathcal{H},[\mathcal{H},p_i]]+\frac{{\beta }_i^3}{4!}[\mathcal{H},[\mathcal{H},[\mathcal{H},p_i]]]+\dots$$

$$g_i=\sum _{j=1}^L(a_{ij}{p}_j+i{b}_{ij}{q}_j),$$

$$i\hslash \frac{d}{dt}⟨f⟩=⟨[f,H]⟩-\mathrm{Tr}[f,q_1]J_1-\mathrm{Tr}[f,q_L]J_L,$$

$$\mathrm{Tr}[f,q_i]J_i=-\frac{\gamma }{2}⟨g_i[f,q_i]⟩-\frac{\gamma }{2}⟨[f,q_i]g_i^{†}⟩-\frac{\gamma m}{i\hslash {\beta }_i}⟨[[f,q_i],q_i]⟩.$$

$$\frac{d}{dt}⟨\mathcal{H}⟩={\Phi }_1+{\Phi }_L$$

$$\frac{d}{dt}{x}_{ij}=\frac{1}{m}(z_{ij}+{\tilde{z}}_{ij}),$$

$$\frac{d}{dt}{z}_{ij}=-\sum _{k=1}^L{x}_{ik}{G}_{kj}+\frac{1}{m}{y}_{ij}-\frac{\gamma }{2}\sum _{k=1}^L{a}_{jk}({\tilde{z}}_{ki}+z_{ik}),$$

$$\begin{gathered} \frac{d}{dt}{y}_{ij}=-\sum _{k=1}^L(G_{ik}{z}_{kj}+{\tilde{z}}_{ik}{G}_{kj})+2\gamma m{\delta }_{ij}\left(\frac{{\delta }_{j1}}{{\beta }_1}+\frac{{\delta }_{jL}}{{\beta }_L}\right) \\ -\gamma \sum _{k=1}^L(a_{jk}{y}_{ki}+a_{ik}{y}_{kj})+\frac{\gamma \hslash }{2}(b_{ji}+b_{ij}), \end{gathered}$$

$$\sum _{k=1}^L{x}_{ik}^{\mathrm{e}}{G}_{kj}=\frac{1}{m}{y}_{ij}^{\mathrm{e}},$$

$$\sum _{k=1}^L(A_{jk}{y}_{ki}^{\mathrm{e}}+A_{ik}{y}_{kj}^{\mathrm{e}})=\frac{\hslash }{2}(B_{ji}+B_{ij})+\frac{2m}{\beta }{\delta }_{ij}({\delta }_{j1}+{\delta }_{jL}),$$

$$\frac{d}{dt}{X}_{ij}=\frac{1}{m}(Z_{ij}+{\tilde{Z}}_{ij}),$$

$$\frac{d}{dt}{Z}_{ij}=-\sum _{k=1}^L{X}_{ik}{G}_{kj}+\frac{1}{m}{Y}_{ij}-\frac{\gamma }{2}\sum _{k=1}^L{A}_{jk}({\tilde{Z}}_{ki}+Z_{ik}),$$

$$\begin{gathered} \frac{d}{dt}{Y}_{ij}=-\sum _{k=1}^L(G_{ik}{Z}_{kj}+{\tilde{Z}}_{ik}{G}_{kj})+\gamma m{k}_B{\delta }_{ij}(-{\delta }_{j1}+{\delta }_{jL}) \\ -\gamma \sum _{k=1}^L(A_{jk}{Y}_{ki}+A_{ik}{Y}_{kj})-\gamma (C_{ji}+C_{ij}), \end{gathered}$$

$$C_{ij}\Delta T=\sum _{k=1}^L(a_{ik}-A_{ik})y_{kj}^{\mathrm{e}}-\frac{\hslash }{2}(b_{ij}-B_{ij}),$$

$$\frac{d}{dt}X=\frac{1}{m}(Z+Z^{†}),$$

$$\frac{d}{dt}Z=-XG+\frac{1}{m}Y-\gamma Z{A}^{†},$$

$$\frac{d}{dt}Y=-(GZ+Z^{†}G)+\gamma D-(Y{A}^{†}+AY)-\gamma (C^{†}+C),$$

$$Z=-Z^{†},$$

$$Y=m(XG+\gamma Z{A}^{†}),$$

$$Y=m(GX-\gamma AZ),$$

$$\gamma D-\gamma (Y{A}^{†}+AY)-\gamma (C+C^{†})=GZ-ZG.$$

$$Z_{ij}=\left\{\begin{array}[]{ll}\;\;\varphi_{j-i},&\;\;i<j,\\ \;\;0,&\;\;i=j,\\ -\varphi_{i-j},&\;\;i>j.\end{array}\right.$$

$$XG-GX=-\gamma (AZ+Z{A}^{†}).$$

$$X_{ij}=\left\{\begin{array}[]{ll}\;\;\psi_{i+j-1},&\;\;i+j<L+1,\\ \;\;0,&\;\;i+j=L+1,\\ -\psi_{2L+1-i-j},&\;\;i+j>L+1.\end{array}\right.$$

$$\frac{k_1}{\gamma }{\psi }_{\ell }=-\sum _{j=1}^L{\eta }_{\ell ,j-1}{A}_j{\phi }_{|\ell -j+1|},$$

$$Y_{ij}=\left\{\begin{array}[]{ll}\;\;\theta_{i+j-1},&\;\;i+j<L+1,\\ \;\;0,&\;\;i+j=L+1,\\ -\theta_{2L+1-i-j},&\;\;i+j>L+1,\end{array}\right.$$

Full text

Fifty years ago, Rieder, Lebowitz, and Lieb rieder1967 introduced and exactly solved a microscopic model for thermal conduction, that consisted of a chain of coupled classic harmonic oscillators with its ends in contact with heat reservoirs at distinct temperatures. Using this model they provided a rigorous proof of the well known result that the thermal conductance, the ratio between the heat current and the temperature difference, is finite regardless of the chain length peierls1955 . This result amounts to say that Fourier’s law does not hold because the conductivity, which is the product of the conductance and the chain length, becomes infinite when the length increases without bounds. The reason for the occurrence of a finite conductance is that the excitations in ordered systems with harmonic interactions travel balistically. To get the Fourier’s law, new ingredients should be added to the harmonic model in order to transform the ballistic into a diffusive motion. Such ingredients include anharmonic potentials lepri1997 ; aoki2001 ; eckmann2004 ; mai2007 ; bernardin2011 ; roy2012 , self-consistent reservoirs bolsterli1970 ; bonetto2004 ; pereira2004 ; delfini2006 , energy conserving noise dhar2011 ; landi2013 ; landi2014 , and others narayan2002 ; lepri2003 ; casati2003 ; basile2009 ; iacobucci2010 .

In the model studied by Rieder, Lebowitz, and Lieb rieder1967 , the oscillators were under the action of conservative forces except the first and the last which in addition were subject to dissipating-fluctuating forces representing the contact with heat reservoirs. They are composed by a dissipative force, proportional to the velocity, and a Gaussian white noise with zero mean and variance proportional to the temperature. The equations of motion are understood as Langevin equations, and the equation governing the time evolution of the probability density is a Fokker-Planck-Kramers (FPK) equation kampen1981 ; gardiner1983 ; risken1984 ; tome2015book .

Here, we study a quantum version of the model studied by Rieder, Lebowit, and Lieb rieder1967 . We have exactly calculated the thermal conductance in the regime of small interparticle interaction and reached a similar result that the conductance is finite regardless of the length of the chain. However, as should be expected the conductance is not independent of temperature, as is the case of the classical version. It vanishes in the limit of zero temperature and saturates at the classical value at high temperatures. Our approach is based on a quantum version of the FPK equation coming from a canonical quantization of the ordinary FPK equation, recently introduced oliveira2016 , and differs from other approaches regarding the treatment of quantum dissipation zurcher1990 ; saito2000 ; dhar2003 ; segal2003 ; dhar2006 ; roy2008 ; rivas2010 ; dhar2012 ; asadian2013 ; freitas2014 ; nicacio2015 . These approaches include the use of quantum Langevin equations zurcher1990 ; dhar2006 ; roy2008 , and the use of the rotating wave approximation and Lindblad master equation to describe the contact with heat reservoirs rivas2010 ; asadian2013 . These approaches as well as ours, when applied to the harmonic chain, predict a finite conductance regardless of the chain length zurcher1990 ; asadian2013 . It worth mentioning that the approach we use leads to a proper thermalization in the sense that the Gibbs equilibrium state is the stationary solution of the quantum FPK equation when the reservoirs have the same temperatures oliveira2016 .

The model we consider is a chain of $L$ particles of equal masses interacting through a harmonic potential. The quantum hamiltonian of the system is given by

[TABLE]

where $m$ is the mass of the particles, $k_{0}$ is the spring constant and $k_{1}$ is the interparticle interacting parameter. The position $q_{i}$ and momentum $p_{i}$ obey the usual commutation relation, $[q_{i},p_{j}]=i\hbar\delta_{ij}$. To describe the contact of the system with heat reservoirs we use the quantum Fokker-Planck-Kramers (FPK) equation oliveira2016 with the first particle of the chain in contact with a heat reservoir A at temperature $T_{A}$ and the last with a reservoir B at temperature $T_{B}$. The quantum FPK equation reads oliveira2016

[TABLE]

where $\rho$ is the density matrix, and $J_{1}$ and $J_{2}$ are given by

[TABLE]

where $\gamma$ is the dissipation parameter and $\beta_{i}=1/k_{B}T_{i}$, with $T_{1}=T_{A}$ and $T_{L}=T_{B}$. The operator $g_{i}$ is given by

[TABLE]

When the temperatures are the same, the Gibbs density $\rho_{0}=(1/Z)e^{-\beta{\cal H}}$ is the stationary equilibrium solution of the quantum FPK equation (2) because it makes each one of the $J_{i}$ to vanish and $[{\cal H},\rho_{0}]=0$.

The operator $g_{i}$ can be written in a form involving nested commutators

[TABLE]

Taking into account the bilinear form of ${\cal H}$ and that the coordinates are not coupled to the momenta, we notice that the odd terms in this expansion are linear combinations of the momenta only, and that the even terms are linear combinations of the coordinates only. From these properties, it follows that $g_{i}$ is a linear combination of the positions and momenta,

[TABLE]

where the coefficients $a_{ij}$ and $b_{ij}$ depend on the temperature $T_{i}$ and on the parameters of the Hamiltonian. In addition, using the fact that the odd terms are Hermitian and that the even terms are anti-Hermitian it follows that the coefficient of $p_{j}$ is real and the coefficient of $q_{j}$ is pure imaginary so that $a_{ij}$ and $b_{ij}$ are real.

Next, we wish to write down evolution equations for the averages of quantities of interest in our analysis. The evolution equation for a certain average $\langle f\rangle={\rm Tr}\{f\rho\}$ of an operator $f$ can be obtained from the quantum FPK equation (2) and it is given by

[TABLE]

[TABLE]

From this formula we get at once

[TABLE]

where $\Phi_{i}=(1/m){\rm Tr}p_{i}J_{i}$ is the energy flux from reservoir $i$ to the chain.

Using formula (7) we get the evolution equation for the correlations among the coordinates and momenta, $x_{ij}=\langle q_{i}q_{j}\rangle$, $y_{ij}=\langle p_{i}p_{j}\rangle$, and $z_{ij}=\langle q_{i}p_{j}\rangle$,

[TABLE]

[TABLE]

[TABLE]

where $\tilde{z}_{ij}=\langle p_{i}q_{j}\rangle$, and $G_{ij}$ are the elements of a tridiagonal matrix $G$, with $G_{ii}=k_{0}$, and $G_{i,i+1}=G_{i+1,i}=-k_{1}$. In these equations, we are setting the coefficients $a_{ij}$ and $b_{ij}$ to vanish unless $i=1$ or $i=L$. We wee that equations (10), (11), and (12) make up a closed set of equations for the correlations $x_{ij}$, $y_{ij}$, and $z_{ij}$ by recognizing that $\tilde{z}_{ji}=z_{ij}-i\hbar\delta_{ij}$.

In the stationary state, $\Phi_{L}=-\Phi_{1}$ and $\Phi=\Phi_{L}$ may thus be understood as the heat flux from reservoir B to reservoir A through the chain. Using equations (10) and (12), we can show that $z_{21}=z_{32}=\ldots=z_{L,L-1}$, that $z_{ji}=-z_{ij}$, and that $\Phi=(k_{1}/m)z_{i+1,i}$, a relation that will be used to determine the conductance.

To simplify the evolution equations, we will subtract the equilibrium solution, which we denote by $x_{ij}^{\rm e}$, $y_{ij}^{\rm e}$, and $z_{ij}^{\rm e}$. By equilibrium solution we mean the stationary solution of (10), (11), and (12) for the case when both temperatures of the reservoirs are the same and equal to $T=(T_{1}+T_{L})/2$. In equilibrium $z_{ij}^{\rm e}=\tilde{z}_{ij}^{\rm e}=0$ if $i\neq j$ and $z_{ii}^{\rm e}=-\tilde{z}_{ii}^{\rm e}=i\hbar/2$. The correlations $x_{ij}^{\rm e}$ and $y_{ij}^{\rm e}$ are given by

[TABLE]

[TABLE]

where $\beta=1/k_{B}T$, and $A_{ik}$ and $B_{ik}$ are the values of $a_{ij}$ and $b_{ij}$ that one obtains by replacing both $\beta_{1}$ and $\beta_{L}$ by $\beta$. The quantities $A_{ij}$ and $B_{ij}$ are nonzero only when $i=1$ or $i=L$ in which case they hold the property $A_{Lj}=A_{1,L+1-j}$ and $B_{Lj}=B_{1,L+1-j}$.

We define the deviations $X_{ij}\Delta T=x_{ij}-x_{ij}^{\rm e}$, $Y_{ij}\Delta T=y_{ij}-y_{ij}^{\rm e}$, $Z_{ij}\Delta T=z_{ij}-z_{ij}^{\rm e}$, and $\tilde{Z}_{ij}\Delta T=\tilde{z}_{ij}-\tilde{z}_{ij}^{\rm e}$ from the equilibrium solution, where $\Delta T=T_{L}-T_{1}$. We remark that $\tilde{Z}_{ij}=Z_{ji}$. We wish, in the following, to write down evolution equations for the variables $X_{ij}$, $Y_{ij}$, and $Z_{ij}$ for small values of $\Delta T$. In this regime the evolution equation reads

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

The quantities $C_{ij}$ are nonzero only when $i=1$ or $i=L$ in which case they hold the property $C_{Lj}=-C_{1,L+1-j}$.

These equations are written in matrix form as

[TABLE]

[TABLE]

[TABLE]

where the matrix $D$ has only two nonzero elements, which are $D_{11}=-mk_{B}$ and $D_{LL}=mk_{B}$. It should be noted that all the entries of $X$, $Y$, and $Z$ are real and that $X$ and $Y$ are symmetric matrices. All entries of matrices $G$, $A$ and $C$ are also real and $G$ is symmetric. The matrices $A$ and $C$ have nonzero values only on the first and the last row.

In the stationary state we are left with the equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the classical limit, $C$ vanishes and $A$ becomes a matrix whose only nonzero entries are $A_{11}=1$ and $A_{LL}=1$, and we recover the equations obtained by Rieder, Lebowitz, and Lieb rieder1967 .

Our next step is to seek the solution of equations (22), (23), (24), and (25). To this end we follow the reasoning put forward by Rieder, Lebowitz, and Lieb rieder1967 . We start by observing that the left-hand side of equation (25) is a bordered matrix, that is, a matrix whose nonvanishing entries are found only in the first and last rows and columns. Therefore, the right-hand side $GZ-ZG$ should also be bordered. Using the definition of $G$ and the relation $Z_{ji}=-Z_{ij}$, that comes from (22) and the property that $Z_{ij}$ is real, it follows that $Z$ is an antisymmetric Toeplitz matrix, that is, a matrix of the type

[TABLE]

From equations (23) and (24), we get the following relation between $X$ and $Z$,

[TABLE]

The right-hand side of this equation is again a bordered matrix and so is the matrix $XG-GX$. Since $X$ is required to be symmetric, a solution for $X$ is an antisymmetric Hankel matrix, that is, a matrix of the following form

[TABLE]

Replacing (28) into (27), we find $\psi_{i}$ in terms of $\varphi_{i}$

[TABLE]

$1\leq\ell<L$, where $A_{j}$ stands for $A_{1j}$ and $\eta_{\ell,j}=-1,0,1$ according to whether $j<\ell$, $j=\ell$, $j>\ell$, respectively.

From equations (23) and (24), it is straitforward to show that $Y$ is also an antisymmetric Hankel matrix, that is, a matrix of the following form

[TABLE]

and that $\theta_{\ell}$ is related to $\psi_{\ell}$ by

[TABLE]

$1\leq\ell<L$, where $\psi_{0}=0$ and $\psi_{L}=0$.

Since $\psi_{\ell}$ is related to $\varphi_{\ell}$ by expression (29), then equation (31) gives $\theta_{\ell}$ in terms of $\varphi_{\ell}$. To get a closed equation for $\varphi_{\ell}$ we use equation (25) to obtain another relation between $\varphi_{\ell}$ and $\theta_{\ell}$, namely

[TABLE]

$1\leq\ell<L$, where $C_{\ell}$ stands for $C_{1\ell}$. Therefore, equations (29), (31), and (32) constitute the desired closed equations for the variables $\varphi_{\ell}$. To solve them we need to know $A_{\ell}$ and $B_{\ell}$ because $C_{\ell}$ is related to these quantities by

[TABLE]

which follows from (18), where $B_{\ell}$ stands for $B_{1\ell}$. We recall that $A_{\ell}$ and $B_{\ell}$ depend on $T$ and are the values of $a_{1\ell}$ and $b_{1\ell}$ obtaining by setting $T_{1}$ equal to $T$.

Explicit solutions of equations (29), (31), and (32) are very cumbersome to find, but closed solutions can be found when the interparticle interaction is small. Thus, from now on we will confine ourselves to the case where the interparticle interacting parameter $k_{1}$ is small. To this end, we first notice that the quantities $A_{\ell}$ and $B_{\ell}$ are of the order $k_{1}^{\ell-1}$ and so is $C_{\ell}$, a result that follows from expansion (5). Thus, from equations (29), (31), and (32) it follows that $\varphi_{\ell}$, $\theta_{\ell}$, and $\psi_{\ell}$ are of order greater or equal $k_{1}$, except $\theta_{1}$ and $\psi_{1}$, which are

[TABLE]

Using equation (29), we see that $\varphi_{1}$ is of the order $k_{1}$ and given by

[TABLE]

The conductance $K$ is defined as the ratio $K=\Phi/\Delta T$ in the limit $\Delta T\to 0$. To determine $K$, we recall that, in the stationary state, the heat flux $\Phi=(k_{1}/m)z_{i+1,i}=(k_{1}/m)Z_{i+1,i}\Delta T$ so that the conductance is $K=(k_{1}/m)Z_{i+1,i}=-(k_{1}/m)\varphi_{1}$, leading us to following expression for the conductance

[TABLE]

In this formula the value of $C_{1}$ is

[TABLE]

Using formula (5), we obtain explicit expression for $A_{1}$ and $B_{1}$, which for $k_{1}=0$ reads

[TABLE]

where $\omega=\sqrt{k_{0}/m}$. When $k_{1}=0$, $y_{11}^{\rm e}=mE$ where

[TABLE]

A straightforward algebra leads us to the result

[TABLE]

In the classical limit, $E=k_{B}T$ and $A_{1}=1$ so that $K=K_{0}$. Thus $K_{0}$ is the classical conductance, a result obtained by Rieder, Lebowitz, and Lieb rieder1967 , in the regime of small $k_{1}$. An explicit form for $K$ is

[TABLE]

As regards the behavior with temperature, this result is qualitatively similar to the ones obtained by other approaches roy2008 ; asadian2013 ; nicacio2015 . At high temperatures all these approaches give the classical conductance obtained in referece rieder1967 , and vanish when $T\to 0$. However, the behavior at low temperature is distinct. Our results give the behavior $\beta^{3}e^{-2\beta\hbar\omega}$ for the conductance whereas the reference roy2008 , for instance, gives the behavior $\beta^{1/2}e^{-\beta\hbar\omega}$.

In conclusion, we have used a quantum FPK equation to cond the contact of a chain of coupled quantum harmonic oscillators with heat reservoirs at distinct temperatures. Starting from the quantum FPK equation, we have set up evolution equations for the covariances and solved them in the stationary regime to get the thermal conductance. An exact closed form for the conductance was obtained for small values of the interparticle interacting parameter. The conductance was found to be finite regardless of the chain length, implying an infinite conductivity and thus the absence of Fourier’s law. This is a consequence of the ballistic motion of phonons that occurs in a system with harmonic interactions where the phonons do not interact.

The quantum FPK equation we use holds two important properties with relevant consequences to our analysis. In equilibrium, that is, when the heat baths have the same temperatures, its stationary state is the Gibbs state. Second, it is a canonical quantization of the ordinary FPK equation oliveira2016 and, as a consequence it turns into this equation, in the classical limit. Thus, in the classical limit our approach becomes identical to that of Rieder, Lebowitz, and Lieb rieder1967 . We remark that our quantum approach differs from other approaches used to study the thermal transport such as the use of quantum Langevin equations with dissipation proportional do velocity or the use of Lindblad dissipators. Although the quantum FPK equation we used here has not been derived from a full quantum system that includes the reservoirs, we think that the canonical quantization might give support for the validity of the quantum FPK equation, although this procedure does not give an unambiguous prescription for obtaining a quantum version of a classical equation oliveira2016 . The present calculation of the conductance, on the other hand, might give indirectly the desired support if we uderstand that the conductance obtained here is a reasonable result.