Thermal rectification and negative differential thermal conductance in harmonic chains with nonlinear system-bath coupling

TL;DR

This paper demonstrates thermal rectification and negative differential thermal conductance in harmonic chains by incorporating nonlinear system-bath coupling, revealing how coupling strength influences heat transport and rectification direction.

Contribution

It introduces a generalized Caldeira-Leggett model with nonlinear system-bath coupling, extending previous linear models to explain thermal phenomena in harmonic chains.

Findings

Weak linear coupling enhances heat current via multiphonon processes.

Asymmetric nonlinear coupling induces thermal rectification.

Strong linear coupling suppresses heat current and reverses rectification direction.

Abstract

Thermal rectification and negative differential thermal conductance were realized in harmonic chains in this work. We used the generalized Caldeira-Leggett model to study the heat flow. In contrast to the most previous studies considering only the linear system-bath coupling, we considered the nonlinear system-bath coupling based on recent experiment [A. Eichler \emph{et al.}, Nat. Nanotech. \textbf{6}, 339 (2011)]. When the linear coupling constant is weak, the multiphonon processes induced by the nonlinear coupling allow more phonons transport across the system-bath interface and hence the heat current is enhanced. Consequently, thermal rectification and negative differential thermal conductance are achieved when the nonlinear couplings are asymmetric. However, when the linear coupling constant is strong, the umklapp processes dominate the multiphonon processes. Nonlinear coupling suppresses the heat current. Thermal rectification is also achieved. But the direction of rectification is reversed comparing to the results of weak linear coupling constant.

Equations

$$H_S=\sum _{l=1}^N\frac{p_l^2}{2m_l}+\sum _{l=0}^{N+1}\frac{(x_l-x_{l+1}{)}^2}{2},$$

$$H_B=\sum _{\alpha =1}^M\frac{P_{\alpha }^2}{2}+\sum _{\alpha ,\beta }\frac{1}{2}{K}_{\alpha \beta }{Q}_{\alpha }{Q}_{\beta },$$

$$H_I=-g(x_1)Q_L-f(x_N)Q_R.$$

$$m_l{\ddot{x}}_l=-(2x_l-x_{l-1}-x_{l+1})-{\gamma }_l(x_l){\dot{x}}_l+{\xi }_l(x_l).$$

$$\frac{\partial P}{\partial t}=-\sum _{l=1}^N\frac{\partial }{\partial x_l}(v_{l}P)-\sum _{l=1}^N\frac{\partial }{\partial v_l}\left[\frac{-(2x_l-x_{l-1}-x_{x+1})}{m_l}P\right]$$

$$J^{st}=\frac{1}{2}(J_L^{st}-J_R^{st})$$

$$\displaystyle J_{L(R)}^{st}=\frac{\gamma_{L(R)}k_{B}T_{L(R)}}{m_{1(N)}}\Big{(}k_{L(R)}^{2}+\mu_{L(R)}^{2}\big{\langle}x_{1(N)}^{2}\big{\rangle}\Big{)}$$

$$\displaystyle-\gamma_{L(R)}\Big{(}k_{L(R)}^{2}\big{\langle}v_{1(N)}^{2}\big{\rangle}+\mu_{L(R)}^{2}\big{\langle}x_{1(N)}^{2}v_{1(N)}^{2}\big{\rangle}\Big{)}.$$

$$x_l(t)$$

$$\hat{Z}$$

$$\hat{Y}$$

$${\hat{\Phi }}_{lm}$$

$${\hat{M}}_{lm}$$

$${\hat{\Gamma }}_{lm}$$

$${\hat{\xi }}_l$$

$$J_L^{st}$$

$${\hat{M}}^{-1}=\frac{1}{\pi }{\int }_{-\infty }^{\infty }d\omega \hat{Z}(\omega ){\omega }^2\hat{\Gamma }\hat{Z}(-\omega )$$

$$J_L^{st}={\kappa }_0(T_L-T_R)\left[1+\frac{{\mu }_L^2}{\pi k_L^2}({\gamma }_L{k}_L^2{k}_B{T}_L{I}_{11}+{\gamma }_R{k}_R^2{k}_B{T}_R{I}_{1N})\right].$$

$$J_R^{st}={\kappa }_0(T_R-T_L)\left[1+\frac{{\mu }_R^2}{\pi k_R^2}({\gamma }_R{k}_R^2{k}_B{T}_R{I}_{NN}+{\gamma }_L{k}_L^2{k}_B{T}_L{I}_{N1})\right].$$

$$(A)T_L=T_0+\Delta T/2,T_R=T_0-\Delta T/2.$$

$$J^{st}$$

$$(B)T_L=T_s,T_R=T_s-\Delta T.$$

$$J^{st}$$

$$(\Delta T{)}_m=\frac{T_s}{2}+\frac{2\pi +({\mu }_L^2{I}_{11}+\frac{k_L^2}{k_R^2}{\mu }_R^2{I}_{1N}){\gamma }_L{k}_B{T}_s}{2({\mu }_R^2{I}_{NN}+\frac{k_R^2}{k_L^2}{\mu }_L^2{I}_{1N}){\gamma }_R{k}_B}>\frac{T_s}{2}.$$

$$Q_{\alpha }$$

$$P_{\alpha }$$

$$\sum _{\beta =1}^M{K}_{\alpha \beta }{U}_{\beta s}={\omega }_s^2{U}_{\alpha s}$$

$$\sum _{\alpha =1}^M{U}_{\alpha s}{U}_{\alpha s^{\prime }}={\delta }_{s{s}^{\prime }}.$$

$$H_B=\sum _{s=1}^M(\frac{1}{2}{\tilde{P}}_s^2+\frac{1}{2}{\omega }_s^2{\tilde{Q}}_s^2).$$

$$H_I=-\sum _{s=1}^{M}g(x_1)U_{Ls}{\tilde{Q}}_s-\sum _{s^{\prime }=1}^{M}f(x_N)U_{R{s}^{\prime }}{\tilde{Q}}_{s^{\prime }}.$$

$${\hat{Z}}_{11}$$

$${\hat{Z}}_{1N}$$

$${\hat{Z}}_{NN}$$

$$A_{1,N}$$

$$A_{1,N-1}$$

$$A_{2,N}$$

$$A_{1,N-1}(\omega )A_{1,N-1}(-\omega )$$

Full text

I Introduction

In the past decade, phononics – a science and engineering of manipulating heat – has attracted intense interest from fundamental research as well as applied researchLi et al. (2012); Pop (2010); Marconnet et al. (2013); Balandin and Nika (2012); Maldovan (2013); Roberts and Walker (2011). There are two essential effects in phononics, thermal rectification and negative differential thermal conductance (sometimes referred to as negative differential thermal resistance). Thermal rectification allows heat current to flow preferably in one direction. The first nanoscale thermal rectifier was proposed theoretically in 2002 based on an one-dimensional (1D) nonlinear chainTerraneo et al. (2002). Since then, a variety of theoretical thermal rectifiers were proposed based on diverse nonlinear systemsLi et al. (2004, 2005); Hu et al. (2006); Segal and Nitzan (2005a, b); Segal (2006, 2008); Wu and Segal (2009); Eckmann and Mejía-Monasterio (2006); Casati et al. (2007). Inspired by the seminal experimental work demonstrating that heat current flow preferentially along the direction of decreasing mass density in asymmetrically mass-loaded nanotubesChang et al. (2006), a lot of studies on thermal rectification were performed in nonlinear mass graded systemsYang et al. (2007); Zeng and Wang (2008); Pereira (2010a, 2011); Wang et al. (2012a); Romero-Bastida and González-Alarcón (2014); Liu et al. (2014a) and asymmetric carbon based nanostructuresWu and Li (2007); Yang et al. (2008); Lee et al. (2012); Bui et al. (2012) including asymmetric graphene nanoribbonsHu et al. (2009a); Wang et al. (2012b, 2014); Yang et al. (2009); Zhang and Zhang (2011); Zhong et al. (2011a, 2012). All studies attribute thermal rectification to the intrinsic nonlinearity (anharmonicity) of the studied systems.

Negative differential thermal conductance refers to the effect that the heat current decreases as the applied temperature difference increases. This effect is the critical element to realize thermal transistorsLi et al. (2006), thermal logic gatesWang and Li (2007) and thermal memoryWang and Li (2008). Negative differential thermal conductance can be obtained in many nonlinear latticesLan and Li (2006); Li et al. (2006); Yang et al. (2007); He et al. (2009); Shao et al. (2009); Shao and Yang (2011); Zhong et al. (2009, 2011b); Hu and Chen (2013); Pereira (2010a); Ai et al. (2011); Ai and Hu (2011); He et al. (2010); Chan et al. (2014). Graphene nanoribbons are also the suitable platforms for practically realizing the negative differential thermal conductanceHu et al. (2011); Ai et al. (2012); Zhong et al. (2012). The intrinsic nonlinearities of the systems are the necessary conditions to achieve negative differential thermal conductance, although the interface resistance between two-segment nonlinear systemsLan and Li (2006); Li et al. (2006); He et al. (2009); Zhong et al. (2009, 2011b); Shao et al. (2009); Shao and Yang (2011); Hu and Chen (2013) or the boundary resistance between heat baths and nonlinear systemsAi and Hu (2011); He et al. (2010) are important.

However, the phonon mean-free path in graphene ($\sim 775$ nm near room temperatureBalandin (2011)) is much longer than the sizes of graphene nanoribbons. Therefore the intrinsic nonlinearity is insignificant and thus graphene nanoribbons can be regarded as harmonic systems. In harmonic systems, although exceptions existSegal and Nitzan (2005a, b); Segal (2006), negative differential thermal conductance cannot be achieved and thermal rectification can only be obtained in quantum regime by asymmetric coupling with an additional self-consistent heat bathMing et al. (2010); Zhang et al. (2010); Pereira (2010b); *Pereira2011b; *Pereira2011c; Bandyopadhyay and Segal (2011); Ouyang et al. (2010); Xie et al. (2012). However, the nanoscale self-consistent heat bath is hard to realizeBergfield et al. (2013). Therefore, harmonic systems did not receive much interest in researches on thermal rectification and negative differential thermal conductance.

In almost all theoretical models mentioned above, the system-bath couplings were supposed to be linear. This is because the energy dissipations (dampings) in previous studied systems were supposed to be linear in general. However, recent researches have revealed that the nanostructures with high aspect ratio such as nanotubes and graphene nanoribbons can easily be driven into nonlinear dissipation regimeEichler et al. (2011). As shown in Ref. [Eichler et al., 2011], nonlinear dissipation can be treated as a generalized Caldeira-Leggett model with nonlinear system-bath couplingBarik and Ray (2005); Zaitsev et al. (2012). Nonlinear system-bath coupling corresponds to the inelastic boundary phonon scattering. In low-dimensional systems, thermal boundary conductance (also referred to as interfacial thermal conductance) becomes increasingly importantMarconnet et al. (2013); Pop (2010); Cahill et al. (2003); *Cahill14. At high temperature, many experimentalStoner et al. (1992); *Stoner93; Stevens et al. (2005); Lyeo and Cahill (2006); Hopkins et al. (2007, 2008); Chen et al. (2009); Panzer et al. (2010); Oh et al. (2011); Zhang et al. (2012); Norris et al. (2012); Duda et al. (2013); Dechaumphai et al. (2014); Hohensee et al. (2015); Wang et al. (2015), computationalTwu and Ho (2003); Chen et al. (2004); Stevens et al. (2004, 2007); Hu et al. (2009b); Landry and McGaughey (2009); Ong and Pop (2010); Duda et al. (2011a); Chalopin et al. (2012); Pei et al. (2012); Khosravian et al. (2013); Sääskilahti et al. (2014); Liu et al. (2014b, c); Zhang et al. (2015), and empiricalDuda et al. (2010, 2011b); Hopkins (2009); Hopkins and Norris (2009); Hopkins et al. (2011) approaches have uncovered that the thermal boundary conductance at weakly bonded interface (or interface between highly dissimilar materials) exceeds the upper bound of elastic thermal conductance and nearly increases linearly with temperature. These results reveal that inelastic phonon scattering at interface contributes significantly to thermal boundary conductance. Therefore, nonlinear system-bath coupling is non-trivial for studying heat transport in low-dimensional systems.

In contrast to linear system-bath coupling, in thermal transport community, nonlinear system-bath coupling has received a little consideration in the previous works. In Refs. [Segal and Nitzan, 2005a, b; Segal, 2006], using small polaron transformation and based on master equation analysis, nonlinear system-bath coupling was considered and consequently the thermal rectification and negative differential thermal conductance were achieved in nonlinear two-level system and even in a harmonic molecular junction (a single harmonic oscillator). We noted that the nonlinear coupling is so strong that the linear coupling is omitted as shown in the Appendix of Ref. [Segal, 2006]. Therefore, the relative contributions of nonlinear coupling and linear coupling to thermal rectification and negative differential thermal conductance were not addressed.

In this work, we modeled the high-aspect-ratio nanostructure as an 1D harmonic chain. Two heat baths are coupled to it at the ends. The couplings are allowed to be nonlinear in addition to linear. Heat transport in the chain is studied at high temperature limit. When the linear system-bath couplings are weak, the numerical results reveal four effects of nonlinear system-bath coupling on heat current. Firstly, heat current is enhanced when nonlinear couplings are taken into account. When the nonlinear coupling constant is weak, heat current is proportional to the square of the nonlinear coupling constant. Secondly, heat current increases linearly with the average temperature of the two baths when the temperature difference is fixed. When the nonlinear coupling constant is weak, the slope of increasing is proportional to the square of the nonlinear coupling constant. Thirdly, negative differential thermal conductance can be obtained in any temperature region by properly choosing the coupling constants. Lastly, thermal rectification is also obtained when the chain asymmetrically couples to the two baths. When both linear couplings are weak, the higher heat current is obtained when the hot bath couples to the chain with the stronger nonlinear coupling. All numerical results are consistent with our analytical results by approximately solving the generalized Langevin equation. When the linear system-bath couplings are strong, there is no available analytical result. Heat current is calculated numerically. Comparing with the results of weak coupling, heat current is suppressed by the nonlinear couplings. And heat current decreases with the average temperature of the two baths when the temperature difference is fixed. The slope of decreasing is also dependent on the nonlinear coupling constant. Moreover, the direction of thermal rectification is reversed. The higher heat current is obtained when the hot bath couples to the chain with the weaker nonlinear coupling. However, negative differential thermal conductance is not achieved in strong linear coupling regime.

The rest of the paper is organized as follows. In Sec. II, the model and the analytical formulas and results are presented. The numerical results are presented in Sec. III. Finally, we draw the conclusions and discuss the potential experimental realization of the nonlinear system-bath coupling as well as the range of validity of our approximate analytical results in Sec. IV.

II Model and methods

II.1 Model

Harmonic chain contains $N$ particles is considered. These particles are connected by harmonic springs with equal spring constants which are chosen as equal to $1$. Then Hamiltonian of the chain is

[TABLE]

where $x_{l}$, $p_{l}$ and $m_{l}$ denote respectively the displacement of the $l$th particle from its equilibrium position, the momentum and the mass of the $l$th particle. The fixed boundary conditions are chosen as $x_{0}=x_{N+1}=0$. Two uncorrelated heat baths ($L$ and $R$) which are initially in thermal equilibrium at temperatures $T_{L}$ and $T_{R}$ are connected to the $1$st particle and the $N$th particle. Each bath is modeled by a collection of $M$ oscillators with harmonic interactions. The Hamiltonian of each bath is

[TABLE]

where $Q_{\alpha}$ and $P_{\alpha}$ are the displacement and the momentum of the $\alpha$th unit-mass oscillator of the bath, $K_{\alpha\beta}$ is the spring constant between the $\alpha$th and the $\beta$th oscillator of the bath.

The $L$th ($R$th) oscillator of the left (right) bath is connected to the $1$st ($N$th) particle of the chain. The system-bath coupling Hamiltonian is

[TABLE]

Where $g(x_{1})$ and $f(x_{N})$ are functions of $x_{1}$ and $x_{N}$ for describing the coupling strength. If $g(x_{1})$ (or $f(x_{N})$) is proportional to the higher exponent of $x_{1}$ (or $x_{N}$) than $1$, the coupling is nonlinear in the coordinate of chain but linear in the heat-bath coordinates. Then the bath coordinates can be integrated out and we can obtain the generalized Langevin equations with multiplicative noises for the coordinates of chain. It should be mentioned that the coupling Hamiltonian in here is equal to those in Refs. [Barik and Ray, 2005] and [Zaitsev et al., 2012] by just transforming the heat-bath coordinate into its normal-mode coordinate space (as show in Appendix A). However, Hamiltonian (3) exhibits the direct meaning of coupling between two particles.

In the Markovian limit, the following generalized Langevin equations of the chain can be obtained according to the standard procedureBarik and Ray (2005); Zaitsev et al. (2012); Dhar and Sriram Shastry (2003); Dhar and Roy (2006); Dhar (2008) as (see Appendix A)

[TABLE]

Where $\gamma_{l}(x_{l})=\gamma_{L}[g^{\prime}(x_{1})]^{2}\delta_{l,1}+\gamma_{R}[f^{\prime}(x_{N})]^{2}\delta_{l,N}$ denotes dissipation and $\xi_{l}(x_{l})=g^{\prime}(x_{1})\eta_{L}\delta_{l,1}+f^{\prime}(x_{N})\eta_{R}\delta_{l,N}$ is the noise term. The prime ($\prime$) indicates derivative with respect to the corresponding argument. At high temperature (classical limit), the fluctuation-dissipation relations $\langle\eta_{L}(t_{1})\eta_{L}(t_{2})\rangle_{\eta}=2k_{B}T_{L}\gamma_{L}\delta(t_{2}-t_{1})$ and $\langle\eta_{R}(t_{1})\eta_{R}(t_{2})\rangle_{\eta}=2k_{B}T_{R}\gamma_{R}\delta(t_{2}-t_{1})$ are satisfied. Where $\langle\cdots\rangle_{\eta}$ denotes an average over the noise.

To analytically study the heat current flowing in the chain, the Fokker-Planck equation corresponding to Eq. (4) is expressed asReichl (1998)

[TABLE]

Where $v_{l}=p_{l}/m_{l}=\dot{x}_{l}$ is the velocity of the $l$th particle of the chain. $P=P(\{x_{l}\},\{v_{l}\},t)$ is the phase-space probability density function. From the Fokker-Planck equation (5), the time derivative for the energy of the chain $\partial\langle H_{S}\rangle/\partial t$ can be calculated. Where $\langle\cdots\rangle$ implies an ensemble average over the whole phase space of the chain. Then the heat current in steady state with $\partial\langle H_{S}\rangle_{st}/\partial t=0$ can be defined via the continuity equation. For simplicity, in this work, we choose the coupling functions as polynomial in $x_{1}$ and $x_{N}$ only up to the quadratic terms as in Ref. [Zaitsev et al., 2012]: $g(x_{1})=k_{L}x_{1}+\mu_{L}x_{1}^{2}/2$ and $f(x_{N})=k_{R}x_{N}+\mu_{R}x_{N}^{2}/2$. Where $k_{L}$ and $k_{R}$ are linear coupling constants. $\mu_{L}$ and $\mu_{R}$ are nonlinear coupling constants. As a consequence, the steady-state heat current is defined as

[TABLE]

with

[TABLE]

Where the subscripts without brackets corresponding to the heat current flowing into the chain from the left bath and the subscripts in brackets corresponding to the heat current flowing into the chain from the right bath.

II.2 Perturbation approximation

For nonlinear system-bath coupling, there are no rigorous results about $J^{st}$ even for harmonic chain. However, when the couplings are linear, heat current flowing through harmonic chain can be obtained exactly as shown in Refs. [Dhar, 2001; Dhar and Sriram Shastry, 2003; Dhar and Roy, 2006; Dhar, 2008]. Therefore, approximate analytical results can be obtained by using a proper perturbation scheme when the nonlinear couplings are weak. Only considering the linear system-bath coupling by letting $\mu_{L}=\mu_{R}=0$ as the zeroth approximation, $x_{1(N)}$ and $v_{1(N)}$ can be calculated for harmonic chain via the equations of motion (4)Dhar (2001); Dhar and Sriram Shastry (2003); Dhar and Roy (2006); Dhar (2008). Then the heat current can be obtained by inserting these zeroth approximations of $x_{1(N)}$ and $v_{1(N)}$ into Eqs. (6) and (II.1).

In the linear coupling approximation with $\mu_{L}=\mu_{R}=0$, following Refs. [Dhar, 2001; Dhar and Sriram Shastry, 2003; Dhar and Roy, 2006; Dhar, 2008], the equations of motion (4) can be solved exactly by taking the Fourier transformation. The results can be expressed as

[TABLE]

Letters with hat symbol represent the matrices. The superscript $-1$ means the inversion of the corresponding matrix. The corresponding fluctuation-dissipation relations in Fourier space are $\langle\eta_{L(R)}(\omega_{1})\eta_{L(R)}(\omega_{2})\rangle_{\eta}=4\pi\gamma_{L(R)}k_{B}T_{L(R)}\delta(\omega_{1}+\omega_{2})$.

Substituting Eq. (II.2) into Eq. (II.1), the results can be expressed as

[TABLE]

Where

[TABLE]

is used to obtain the second equalityCasher and Lebowitz (1971); Ishii (1973). We should mention that Eq. (10) is satisfied only when $T_{L}=T_{R}$. Therefore, Eqs. (11), (12) and (14) are approximate results for $T_{L}\approx T_{R}$. Furthermore, the second term of the last equality equals to zero because $\omega\hat{Z}_{11}(\omega)\hat{Z}_{11}(-\omega)$ and $\omega\hat{Z}_{1N}(\omega)\hat{Z}_{1N}(-\omega)$ are odd functions of $\omega$. Then Eq. (9) can be simplified as

[TABLE]

Similarly, one can obtain

[TABLE]

Where $\kappa_{0}=\gamma_{L}k_{L}^{2}\gamma_{R}k_{R}^{2}k_{B}K_{1N}/\pi$, $K_{lm}=\int_{-\infty}^{\infty}d\omega\omega^{2}\hat{Z}_{lm}(\omega)\hat{Z}_{lm}(-\omega)$ and $I_{lm}=\int_{-\infty}^{\infty}d\omega\hat{Z}_{lm}(\omega)\hat{Z}_{lm}(-\omega)$. $\kappa_{0}(T_{L}-T_{R})$ is just the heat current obtained in Refs. [Dhar, 2001; Dhar and Sriram Shastry, 2003; Dhar and Roy, 2006; Dhar, 2008] for linear system-bath coupling. Which is valid for high temperature difference. Therefore, in spite of the approximation (10) is used, our results are supposed to valid also for high temperature difference.

As shown in Ref. [Segal, 2006], the temperature difference $\Delta T$ can be imposed in two different ways. Firstly, we set

[TABLE]

Then the heat current (6) can be expressed as

[TABLE]

When the system-bath couplings are linear with $\mu_{L}=\mu_{R}=0$, the heat current $J^{st}=\kappa_{0}\Delta T$ depends linearly on the temperature difference $\Delta T$ but is independent on $T_{0}$ as shown in Eq. (14). However, when the couplings are nonlinear (e.g., $\mu_{L}=\mu_{R}=\mu\neq 0$, $\gamma_{L}=\gamma_{R}=\gamma$ and $k_{L}=k_{R}=k$, then $J^{st}=\kappa_{0}\Delta T+\kappa_{0}\mu^{2}\gamma(I_{11}+I_{1N})k_{B}T_{0}\Delta T/\pi$), heat current is enhanced relative to $\kappa_{0}\Delta T$ and it increases linearly with $T_{0}$ because $I_{lm}\geq 0$ (as shown in Appendix B, $\hat{Z}_{lm}(-\omega)=\hat{Z}^{*}_{lm}(\omega)$). These results are consistent with the results in Refs. [Stoner et al., 1992; *Stoner93; Stevens et al., 2005; Lyeo and Cahill, 2006; Hopkins et al., 2007, 2008; Chen et al., 2009; Panzer et al., 2010; Oh et al., 2011; Zhang et al., 2012; Norris et al., 2012; Duda et al., 2013; Dechaumphai et al., 2014; Hohensee et al., 2015; Wang et al., 2015; Twu and Ho, 2003; Chen et al., 2004; Stevens et al., 2004, 2007; Hu et al., 2009b; Landry and McGaughey, 2009; Ong and Pop, 2010; Duda et al., 2011a; Chalopin et al., 2012; Pei et al., 2012; Khosravian et al., 2013; Sääskilahti et al., 2014; Liu et al., 2014b, c; Zhang et al., 2015; Duda et al., 2010, 2011b; Hopkins, 2009; Hopkins and Norris, 2009; Hopkins et al., 2011].

As one can expect, thermal rectification is absent when the couplings are linear with $\mu_{L}=\mu_{R}=0$. But the situation becomes very different when the system-bath couplings are nonlinear. Thermal rectification is achieved when the second term of Eq. (14) is nonzero. If the system-bath couplings are symmetric with $k_{L}=k_{R}$, $\gamma_{L}=\gamma_{R}$ and $\mu_{L}=\mu_{R}$, thermal rectification is absent. However, if the couplings are asymmetric, thermal rectification can be achieved. These approximate analytical results are verified by the following numerical results.

For the second way to impose the temperature difference, we set

[TABLE]

Where $\Delta T\leq T_{s}$ to ensure $T_{R}\geq 0$. The heat current (6) is thus expressed as

[TABLE]

It indicates that, when $\Delta T>0$, heat current first increases with $\Delta T$, and then decreases after reaching a maximum. From Eq. (16), the temperature difference corresponding to the maximum heat current is

[TABLE]

When letting $k_{R}/k_{L}\gg 1$, by choosing the suitable $\mu_{L}$ and $\mu_{R}$, one can expect that $(\Delta T)_{m}<T_{s}$ from Eq. (17) and thus $T_{R}>0$. Therefore, negative differential thermal conductance occurs with the onset temperature difference being $(\Delta T)_{m}$. One should note that $T_{R}$ has to less than $T_{s}/2$ to achieve the negative differential thermal conductance. However, there is no limitation on temperature region to realize negative differential thermal conductance. In contrast to the results of Ref. [Segal, 2006], if the system-bath couplings are symmetric, the temperature difference calculated from Eq. (17) is larger than $T_{s}$ and thus negative differential thermal conductance cannot be achieved.

III Numerical results

To obtain the heat current in harmonic chains, we use the implicit midpoint algorithmBurrage et al. (2007) to integrate the equations of motion (4). (The results have been compared with those obtained by using Mannella’s leapfrog algorithmBurrage et al. (2007) and the velocity Verlet algorithmAllen and Tildesley (1987). The differences are negligible.) Equilibration times ranged from $10^{8}$ - $10^{9}$ time steps of step size $0.05$ and steady-state averages were taken over another $10^{8}$ - $10^{9}$ time steps. (The results have been compared with those obtained by setting time step size as $0.005$. The differences are negligible.) The steady state is reached by checking whether the results of different equilibration times are equal and checking whether the local heat currents are constant along the chain. In all simulations, we study the equal-mass harmonic chains with $m_{l}=1$. Moreover, we set $k_{B}=1$ and $N=24$. Therefore, the cut-off frequency of the harmonic chain is $2$. Comparing with the cut-off frequency of the out-of plane acoustic (ZA) phonon polarization branches of graphene ($\sim 14$ THzBalandin (2011)), the real temperature $T_{real}$ is related to the dimensionless temperature $T$ through the relation $T_{real}\approx 336\,T$ (K). In this work, $T$ is chosen in the range from $1$ to $2$, thus the corresponding $T_{real}$ is in the range from $336$ K to $672$ K.

The local heat current at site $l$ is defined as $J_{l}=\langle(\dot{x}_{l}+\dot{x}_{l-1})f_{l,l-1}\rangle/2$, where $f_{l,l-1}$ is the force exerted by the $(l-1)$th particle on the $l$th particle and $\langle\cdots\rangle$ denotes a steady-state average. At steady state, $J_{l}$ is independent on the site position $l$. In our simulations, the heat current flows from the left bath to the right bath is defined as $J_{+}=\sum_{l=2}^{N}J_{l}/(N-1)$. Reversing the temperature difference, the heat current flows in the reverse direction is denoted as $J_{-}$.

III.1 Weak linear coupling constant

We set $\mu_{R}=0$ in Fig. 1 and the inset. Only the left system-bath coupling is nonlinear. The linear coupling constants are set as $k_{L}=0.1$ and $k_{R}=1$ in Fig. 1. As revealed in Ref. [Zhang et al., 2011], without the nonlinear coupling, phonons in the whole frequency domain can transport across the right system-bath interface with the transmission equals to one. However, only the low-frequency phonons can transport across the left system-bath interface. If the left system-bath coupling is nonlinear, high-frequency phonons can transport across the left interface via the multiphonon process and thus the heat current is enhanced just as depicted in Fig. 1 and the inset. The enhancement is consistent with our analytical results in Sec. II.2 and the results of weak coupling in Refs. [Landry and McGaughey, 2009; Chen et al., 2009; Hopkins and Norris, 2009; Duda et al., 2013; Norris et al., 2012; Duda et al., 2010; Hohensee et al., 2015; Hopkins et al., 2011; Sääskilahti et al., 2014; Stoner et al., 1992; *Stoner93; Stevens et al., 2004; Lyeo and Cahill, 2006; Stevens et al., 2005]. Moreover, $J_{+}$ and $J_{-}$ are proportional to the square of $\mu_{L}$ in Fig. 1 when $\mu_{L}$ is small as predicted in Eq. (14).

Thermal rectification is apparent in Fig. 1. Heat current flow preferably from the right bath to the left bath (i.e. $J_{-}>J_{+}$) when $\mu_{L}\neq 0$. This is because the right system-bath interface is transparent to phonons. If the right bath is hotter, more phonons can be excited in the chain to participate in the multiphonon processes. Therefore, thermal rectification with $J_{-}>J_{+}$ is obtained. With the increasing of nonlinear coupling constant $\mu_{L}$, the transmission of phonons across the left interface approach saturation. Consequently, the heat currents saturate as shown in Fig. 1. However, when the right coupling is weak with $k_{R}=0.1$, thermal rectification is not evident. $J_{+}$ is only little higher than $J_{-}$ when $\mu_{L}<0.6$ as shown in the inset of Fig. 1. This is because $k_{R}=0.1$ allows only low-frequency phonons transport across the right interface. Benefiting by the left nonlinear system-bath coupling, more low-frequency phonons can be excited in harmonic chain. Therefore, $J_{+}$ is little higher than $J_{-}$ when $\mu_{L}<0.6$. When $\mu_{L}>0.6$, the heat currents ($J_{+}$ and $J_{-}$) approach the saturated value. This saturated value at $\mu_{L}=1$ is marked as an open diamond point in Fig. 1 by corresponding the right interface in the inset to the left interface in Fig. 1. The agreement reveals the fact that the transmissions of low-frequency phonons across the left interface approach one when $\mu_{L}>0.6$ in the inset.

Heat current is further enhanced when the right coupling is nonlinear. Comparing with the inset of Fig. 1, we set $\mu_{R}=0.1$ in Fig. 2. Thermal rectification is apparent with $J_{+}/J_{-}\approx 1.336$ at $\mu_{L}=1$. In addition, the direction of thermal rectification is reversed at $\mu_{L}=\mu_{R}=0.1$. This reversing of thermal rectification indicates that heat current is higher when the hot bath is coupled to the chain with stronger nonlinear coupling constant. It is consistent with the results of nonlinear coupling in Ref. [Segal, 2006]. This is based on the aforementioned mechanism that the stronger the nonlinear coupling is, the more the phonons can be excited in the chain to participate in the multiphonon process. The saturated values of $J_{+}$ and $J_{-}$ at $\mu_{L}=1$ are marked as the open up-triangle and the open down-triangle in Fig. 1 by corresponding the right interface in Fig. 2 to the left interface in Fig. 1. We can find that the corresponding values are equal. This indicates that the left interface is transparent for all the phonons which can transport across the right interface when $\mu_{L}>0.6$. Based on this mechanism, thermal rectification will absent (i.e. $J_{+}=J_{-}$) when both $\mu_{L}$ and $\mu_{R}$ are higher than $0.6$. This is confirmed in the inset of Fig. 2 with $\mu_{R}=0.7$. Moreover, in the inset, $J_{-}>J_{+}$ when $\mu_{L}<0.6$ because $\mu_{R}>\mu_{L}$.

In the high temperature limit, the phonon population in heat bath increases linearly with temperature. As a consequence, the heat conductivity increases with temperature when the system-bath coupling is nonlinear. The heat conductivity is defined as $\lim\limits_{\Delta T\rightarrow 0}J_{st}/\Delta T$. When the system-bath couplings are symmetric with $\mu_{L}=\mu_{R}$ and $k_{L}=k_{R}$, heat current $J_{st}$ is proportional to the temperature difference $\Delta T$. Therefore, fixing $\Delta T$, heat current $J_{st}$ increases linearly with the average temperature $T_{0}$, which is predicted by Eq. (14). The numerical results in Fig 3 confirm this prediction. Moreover, $a_{1}\approx c_{1}$ and $a_{2}\approx c_{2}$ are obtained. Which indicates that the slopes of the linear fits are $\mu_{L(R)}$ dependent when $k_{L}=k_{R}$. The same $\mu_{L(R)}$ corresponds to the same slope. The higher the $\mu_{L(R)}$ is, the higher the slope is. This is consistent with Eq. (14) and Ref. [Duda et al., 2011a]. As predicted in Eq. (14), when the system-bath couplings are symmetric, the slope is proportional to $\mu_{L}^{2}$ and the intercept equals to $\kappa_{0}\Delta T$. The intercept $d_{1}$ is plotted in the inset of Fig. 1 as a solid down-triangle point. It approaches the predicted value $\kappa_{0}\Delta T$ for $k_{L}=k_{R}=0.1$. However, the intercept $d_{2}$ corresponding to $\mu_{L}=\mu_{R}=0.2$ is larger than $\kappa_{0}\Delta T$ for $k_{L}=k_{R}=0.1$. In addition, the ratios of the slopes are $a_{2}/a_{1}\approx 3.43$ and $c_{2}/c_{1}\approx 3.52$. They are smaller than the corresponding ratio of $\mu_{L}^{2}$, which is $(0.2/0.1)^{2}=4$. Besides the fitting errors and the numerical errors, we attribute these discrepancies to the fact that the approximation in Eq. (14) is crude and it is only valid for small $\mu_{L(R)}$.

When the temperature difference is imposed as model B (15), negative differential thermal conductance can be obtained as predicted in Eq. (17). According to Eq. (17), we set $\mu_{R}=0$, $k_{R}=1$, $k_{L}=0.05$ and thus $k_{R}/k_{L}=20$, the negative differential thermal conductance is obtained in Fig. 4. As predicted in Eq. (17), all the onset temperature differences are larger than $T_{L}/2=1$. In addition, the higher the $\mu_{L}$ is, the higher the onset temperature difference is. The appearance of negative differential thermal conductance can be attributed to the same mechanism as in Ref. [Hu et al., 2011]. With the increasing of $\Delta T$, the average temperature $(T_{L}+T_{R})/2$ decreases. When the heat conductivity decreases with the decreasing $(T_{L}+T_{R})/2$, the negative differential thermal conductance may be obtained. This mechanism is verified in Fig. 5. The heat current $J_{+}$ decreases with the decreasing $T_{0}$. The higher the $\mu_{L}$ is, the faster the decreasing of $J_{+}$ is, and thus the faster the decreasing of the corresponding $J_{+}$ in the negative differential thermal conductance region is as shown in Fig. 4.

III.2 strong linear coupling constant

With the increasing of the linear coupling constant, more high-frequency phonons can transport across the system-bath interfaceZhang et al. (2011). Therefore, the heat current increases with the linear coupling constant as shown in the inset of Fig. 6. In addition, when the linear coupling is weak, the heat current is proportional to the square of the linear coupling, which is predicted by $\kappa_{0}\Delta T$ in Eq. (14) and is consistent with the results in Refs. [Hu et al., 2009b; Shen et al., 2011; Rieder et al., 1967]. (One should note that $\gamma_{L}k_{L}^{2}$ and $\gamma_{R}k_{R}^{2}$ in this work correspond to the friction constant $\lambda$ in Ref. [Rieder et al., 1967].) Without the nonlinear coupling, thermal rectification is absent as shown in the inset. However, when the asymmetric nonlinear couplings are present, thermal rectification appears in Fig. 6. Moreover, the direction of thermal rectification reverses as $k_{L}(=k_{R})$ increasing. When the linear coupling is weak, the nonlinear system-bath coupling can enhance the heat current as aforementioned (if $k_{L}$ is weak, heat currents depicted as square symbols are higher than heat currents depicted as up-triangle symbols in Fig. 6). The higher nonlinear coupling enables more phonons transport across the interface, and thus the heat current is larger when the hot bath is coupled to the system with higher nonlinear coupling constant. With the increasing of $k_{L}$ and $k_{R}$, the linear coupling becomes strong. Consequently, more high-frequency phonons can transport across the interface. Therefore, the umklapp process becomes dominating and thus the nonlinear coupling suppresses the heat current (if $k_{L}$ is strong, heat currents depicted as square symbols are lower than heat currents depicted as up-triangle symbols in Fig. 6). The higher the nonlinear coupling is, the more the heat current is suppressed and thus $J_{+}<J_{-}$ when $k_{L}>0.4$. It should be mentioned that even without the nonlinear couplings, the heat current will also be suppressed when the linear couplings are strong enough as shown in Ref. [Rieder et al., 1967]. However, the mechanism of suppression is potentially the mismatching between frequencies of the bath and the system. Which is different from the mechanism of suppression by the nonlinear couplings.

The impacts of nonlinear system-bath coupling on heat current is shown in Fig. 7 when the linear couplings are strong ($k_{L}=k_{R}=1$). The heat currents $J_{+}$ and $J_{-}$ are suppressed by the nonlinear couplings and thus decrease with $\mu_{L}$. Higher nonlinear coupling suppresses more heat current. Hence, thermal rectification appears when $\mu_{L}\neq\mu_{R}$ and the direction of thermal rectification reverses at $\mu_{L}=\mu_{R}$ with the increasing of $\mu_{L}$.

At high temperature, the population of phonons increases linearly with temperature. Consequently, when the nonlinear coupling is present and the linear coupling constant is strong, there are more high frequency phonons participate the umklapp processes with the increasing temperature. Therefore, the heat conductivity decreases linearly with temperature. As described above, heat current $J_{+}$ will decrease linearly with the average temperature $T_{0}$. This is confirmed in Fig. 8. We obtained that $c_{1}\approx c_{2}$. Which means the slopes of the linear fits for $\mu_{L}=\mu_{R}=0.2$ are equal. Although $a_{1}\neq a_{2}$, the corresponding fitting lines (the long-dash line and the short-dash line) in Fig. 8 approach parallel. We attribute this discrepancy to the fitting errors and the numerical errors. Therefore, the slopes of the linear fits are also $\mu_{L(R)}$ dependent when $k_{L}=k_{R}$. The higher the $\mu_{L(R)}$ is, the higher the slope is. This is coincide with the results of weak linear coupling. Furthermore, the intercepts are equal, i.e., $b_{1}\approx d_{1}$ and $b_{2}\approx d_{2}$. According to the results of weak linear coupling, these intercepts approach the corresponding heat currents of harmonic chain with only the linear system-bath couplings. This is confirmed in Fig. 6. Where the cross symbol and the star symbol correspond to the intercepts $b_{1}$ and $b_{2}$ respectively.

Although the heat conductivity decreases with the temperature when the linear coupling is strong, we have not achieved the negative differential thermal conductance by fixing $T_{L}$ but increasing $T_{R}$ (or fixing $T_{R}$ but increasing $T_{L}$). We attribute the absence of negative differential thermal conductance to that the slope of decreasing (see Fig. 8) is much lower than the slope of increasing (see Fig. 3).

IV Conclusion and discussion

In summary, heat flow in harmonic chain with nonlinear system-bath coupling is studied based on the generalized Caldeira-Leggett model in this work. The obtained Langevin-like equations of motion are solved analytically and numerically. When the linear coupling constant is weak, the numerical results are consistent with the predictions of the approximate analytical results. The heat current is enhanced by the nonlinear system-bath coupling. This is attributed to the fact that the weak linear system-bath coupling allows only the low-frequency phonons to transport across the system-bath interface. When the nonlinear coupling is present, the high-frequency phonons can transport across the interface through the multiphonon processes. Hence the heat current is enhanced. The stronger nonlinear coupling enables more phonons to transport across the interface. Therefore, thermal rectification is obtained when the nonlinear couplings are asymmetric. When both linear couplings are weak, higher heat current is obtained when the hot bath is coupled to the chain with the stronger nonlinear coupling. Moreover, the populations of phonons increase linearly with temperature at high temperature. Therefore, the heat conductivity increases linearly with temperature when the nonlinear system-bath coupling is present. As predicted by the analytical results, by suitable choosing of coupling constants, the negative differential thermal conductance is achieved.

However, when the linear coupling constant is strong, high-frequency phonons can transport across the system-bath interface through linear coupling. The umklapp processes dominate the multiphonon processes when the nonlinear coupling is present. Hence the heat current is suppressed. The stronger nonlinear coupling suppresses more heat current. Therefore, in contrast to the results of weak linear coupling constant, the direction of thermal rectification is reversed, namely, higher heat current is obtained when the hot bath is coupled to the chain with the weaker nonlinear coupling. However, although the heat conductivity decreases with the temperature, the negative differential thermal conductance is not achieved in this work. The potential reason is attributed to the slow decreasing of heat conductivity with temperature.

As stated in Sec. II.2, the zeroth approximation is derived for weak nonlinear coupling constant. In addition, the numerical results indicate that the validity of the approximate analytical results is limited to the weak linear coupling constant. It is not valid for the strong linear coupling constant unless the nonlinear coupling constant equals to zero. In deriving the analytical results, Eq. (10) is used. It is satisfied only when $T_{L}=T_{R}$. However, when the nonlinear coupling constants equal to zero, the analytical results coincide with the reported results in Refs. [Dhar, 2001; Dhar and Sriram Shastry, 2003; Dhar and Roy, 2006; Dhar, 2008]. Therefore, we expect that the analytical results are valid for high temperature difference. This is consistent with the numerical results of weak linear coupling constant.

All the results obtained in this work are based on the nonlinear system-bath coupling. Nonlinear system-bath coupling is non-trivial. The experiments on thermal boundary conductance (interfacial thermal conductance)Stoner et al. (1992); *Stoner93; Stevens et al. (2005); Lyeo and Cahill (2006); Hopkins et al. (2007, 2008); Chen et al. (2009); Panzer et al. (2010); Oh et al. (2011); Zhang et al. (2012); Norris et al. (2012); Duda et al. (2013); Dechaumphai et al. (2014); Hohensee et al. (2015); Wang et al. (2015) reveal that the inelastic phonon scattering at interface between highly dissimilar materials is the dominant reason for the enhancement of heat current. Additionally, at the interface between similar materials, the suppression of heat current relative to the elastic thermal conductance is also observedHohensee et al. (2015); Stevens et al. (2005). Our results of strong linear coupling constant indicate that the nonlinear coupling is one potential reason for the suppression of heat current. Therefore, the nonlinear coupling between different materials is intrinsic. Especially for the nanostructures with high aspect ratio, the nonlinear dissipation is easy achieved. The nonlinear dissipation is significant for nanostructures under tensile stress, but is negligible for them with slackEichler et al. (2011). This is consistent with the results of thermal boundary conductance in Ref. [Sääskilahti et al., 2014]. In which, under tensile stress, the transmission of phonons through the linear coupling is suppressed but the transmission of inelastic energy is nearly unaffected. This can be understood as follows. The applied tensile stress weakens the linear coupling constant but almost does not impact the nonlinear coupling constant. Therefore, the relative strength of linear coupling and nonlinear coupling can be tuned by applying pressureHohensee et al. (2015); Sääskilahti et al. (2014).

We hope that our study motivates further research on thermal rectification and negative differential thermal conductance in nanostructures with nonlinear dissipation.

Acknowledgements.

We thank the referees for their constructive comments. Z. J. Ding is supported by the National Natural Science Foundation of China (No. 11274288), the National Basic Research Program of China (Nos. 2011CB932801 and 2012CB933702), Ministry of Education of China (No. 20123402110034) and “111” project (No. B07033). Some numerical calculations in this work were performed on the supercomputing system in the Supercomputing Center of University of Science and Technology of China.

Appendix A Transformation of Hamiltonian

According to Refs. [Dhar and Sriram Shastry, 2003; Dhar and Roy, 2006; Dhar, 2008], the Hamiltonian of each bath (2) can be transformed into the normal-mode form by a canonical transformation

[TABLE]

Where $U_{\alpha s}$ satisfies

[TABLE]

Hence, Hamiltonian (2) is transformed as

[TABLE]

Additionally, the coupling Hamiltonian (3) can be transformed into

[TABLE]

The transformed Hamiltonians (20) and (21) are coincide with them in Refs. [Zaitsev et al., 2012; Barik and Ray, 2005]. Therefore, the generalized Langevin equation (4) can be obtained according to Refs. [Zaitsev et al., 2012; Barik and Ray, 2005].

Appendix B Entries of matrix $\hat{Z}$

As shown in Eq. (II.2), matrix $\hat{Z}$ is just the inversion of matrix $\hat{Y}$. According to Ref. [Huang and McColl, 1997], the entries of matrix $\hat{Z}$ can be calculated as

[TABLE]

with

[TABLE]

Where $A_{l,m}$ and $D_{l,m}$ are defined as the determinants of the submatrices of $\hat{Y}$ and $\hat{\Phi}-\omega^{2}\hat{M}$ beginning with the $l$th row and column and ending with the $m$th row and column. Obviously, $A_{l,m}(-\omega)=A_{l,m}^{*}(\omega)$ is satisfied, where the star (∗) implies the complex conjugate. Therefore,

[TABLE]

with

[TABLE]

Hence, when $\gamma_{L}k_{L}^{2}=\gamma_{R}k_{R}^{2}$ and the harmonic chain is equal-mass, one can obtain $\hat{Z}_{11}=\hat{Z}_{NN}$ and then $\hat{I}_{11}=\hat{I}_{NN}$ in Eq. (14).