An Atomistic Framework for Time-Dependent Thermal Transport - The

Subscriber access provided by University of South Dakota

C: Physical Processes in Nanomaterials and Nanostructures

An Atomistic Framework for Time-Dependent Thermal Transport Leonardo Medrano Sandonas, Alexander Croy, Rafael Gutierrez, and Gianaurelio Cuniberti J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b06598 • Publication Date (Web): 16 Aug 2018 Downloaded from http://pubs.acs.org on August 17, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

An Atomistic Framework for Time-Dependent Thermal Transport Leonardo Medrano Sandonas,†,‡ Alexander Croy,∗,† Rafael Gutierrez,† and Gianaurelio Cuniberti†,¶,§ †Institute for Materials Science and Max Bergmann Center of Biomaterials, TU Dresden, 01062 Dresden, Germany ‡Max Planck Institute for the Physics of Complex Systems, 01187 Dresden, Germany ¶Dresden Center for Computational Materials Science (DCMS), TU Dresden, 01062 Dresden, Germany §Center for Advancing Electronics Dresden, TU Dresden, 01062 Dresden, Germany E-mail: [email protected] Phone: +49 (0)351 4633 1419. Fax: +49 (0)351 4633 1422

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Phonons play a major role for the performance of nanoscale devices and, consequently, a detailed understanding of phonon dynamics is required. Using an auxiliarymode approach, which has successfully been applied for the case of electrons, we develop a new method to numerically describe time-dependent phonon transport. This method allows one to gain insight into the behavior of local vibrations in molecular junctions, which are driven by time-dependent temperature differences between thermal baths. Exemplarily, we apply the method to study the non-equilibrium dynamics of quantum heat transport in an 1D atomic chain as well as in realistic molecular junctions made of poly-acetylene and poly-ethylene chains, in which case the vibrational structure of the junction is described at the density-functional theory level. We calculate the transient energies and heat currents and compare the latter to the standard Landauer approach in thermal equilibrium. We show that the auxiliary-mode representation is a powerful and versatile tool to study time-dependent thermal transport in nanoscale systems.

INTRODUCTION The feasibility of nanoscale functional devices depends crucially on the ability to control the relevant microscopic degrees of freedom: electrons for nanoelectronics, 1,2 spins in spintronics, 3,4 and phonons in nanophononics. 5,6 While many problems belonging to the former two fields have already been in the focus of intensive experimental and theoretical efforts for a long time, the issue of tuning thermal transport by influencing vibrational degrees of freedom in nanophononic devices is a relatively new field. 5–10 So far, the majority of theoretical studies of quantum phonon transport in nanoscale systems deals with stationary situations within the harmonic approximation. Here, the Landauer approach, based on the calculation of transmission functions, can be used. 11–14 To include the anharmonic contributions one needs to resort to full nonequilibrium Green’s function (NEGF) schemes. 15–18 Both methodologies can eventually be combined with density functional-based modelling to achieve an 2


Page 2 of 21

Page 3 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


atomistic level of description. 18–24 In non-stationary situations, where time-dependent (TD) external parameters can strongly affect a nanoscale system, the dynamics of the vibrational system becomes crucial. For example, a time-varying temperature bias 25–27 or strong local heating by laser fields 28–30 can be used to exert additional control over thermal transport. In this context, novel nonequilibrium effects like molecular heat pumping, 31,32 cooling, 33 and rectification 34–36 have been proposed. The description of such phenomena often requires working directly in the time domain, which is very challenging from a numerical point of view. In this respect, considerable progress has been achieved in the description of TD transport of spins 37–39 and electrons. 40–53 However, less attention has been paid to a similar treatment of vibrational degrees of freedom. Thus, only recently new approaches have been developed to deal with TD thermal phenomena in nanojunctions. 54–56 Yet, the application of those approaches is mostly limited to simple model systems, so that the issue of combining the time-dependent approach with an atomistic material description remains open and very challenging. In this article, we present a novel atomistic method capable of studying time-dependent phonon transport in nanosystems by combining a TD-NEGF approach with DFT-based modelling. This computational method is based on the solution of the equation of motion of the phonon density matrix by developing an auxiliary-mode approach. This very efficient scheme has already been successfully used for various investigations of time-resolved electron transport. 41–46 As we show, the efficiency of the auxiliary-mode scheme can be used to study phonon transport in open nanoscale systems. As a proof-of-principle, we first consider an one-dimensional atomic chain and, in a second step, carbon-based molecular junctions to demonstrate the potential of the method to describe time-dependent transport in realistic nanoscale materials.

3



Page 4 of 21

MODEL AND COMPUTATIONAL METHOD The basic setup of our approach is shown in Figure 1. Two thermal baths consisting of non-interacting harmonic oscillators are connected to a generic scattering region, whose vibrational properties are assumed to be well represented by a purely harmonic Hamilton operator. The total system is then described by the following Hamiltonian

H = HC +

X 1 αk

2

p2αk

1 2 2 + ωαk uαk 2

+

X1 α,k

2

T uT · Vαk uαk + uαk Vαk ·u

.

(1)

The first term HC = (1/2)pT · p + (1/2)uT · Keff · u is the Hamiltonian of the central region, u is a column vector consisting of all the atomic displacement variables in the region, and p contains the corresponding momenta. Both vectors have length N , where N is the number of vibrational degrees of freedom in the central region. We have chosen renormalized √ displacements ui = mi xi , where mi is the mass associated with the ith vibrational degree of freedom and xi is the actual displacement having the dimension of length. The effective force-constant matrix Keff = K + Kct has dimension N × N and includes the force constant matrix of the central region, K, as well as a counter-term Kct , whose meaning is explained below. The index α ∈ {L, R} labels the left (L) and right (R) heat baths, while k denotes their corresponding vibrational mode with frequency ωαk . The second term of Equation (1) is the Hamiltonian of the heat bath. The last term represents the interaction between the central region and the baths, which is given by coupling vectors Vαk that are assumed to vanish before time t0 → −∞. Written in this form, the coupling leads to a renormalization of the bare dynamical matrix, which can be canceled by the previously mentioned counter P T 2 term Kct = α,k (Vαk · Vαk )/ωαk . Consequently, the coupling to the thermal baths will solely introduce dissipation rather than a shift in the vibrational spectrum of the system. 57 The energy of the central region, EC (t), can be expressed in terms of the phonon density matrix σ(t) = iG < (t, t), where G < (t, t0 ) is the lesser Green function (GF), as EC (t) = T (1/2)Tr Keff · Q · σ(t) . Here and in the following we adopt units with ~ = 1. Moreover, we 4


Page 5 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


have introduced the 2N × 2N dimensional auxiliary matrices (we use caligraphic symbols to denote matrices with dimension 2N × 2N troughout the manuscript) 





 I 0  I≡ , 0 I



 0 I  Q≡ , −I 0

  Keff ≡ 

 0 −Keff

I   . 0

If no external force is acting on the system, the total energy must be conserved and the ∂ heat current coming from both heat baths can be defined as J(t) = − EC (t). Therefore, ∂t the time evolution of the heat current is related to the lesser GF, which can be written as a 2N × 2N block matrix, 55,56 



0 T 0 T u(t )p (t)   u(t )u (t) G < (t, t0 ) = −i 

0 T  . p(t0 )uT (t) p(t )p (t)

(2)

Note that the GF includes displacement u and momentum p vectors on the same footing. By using the Hamiltonian given by Equation (1) and the GF relations given in the Supporting Information (see Equations (S1-S3)), we find that the density matrix satisfies the following equation of motion (EOM) X ∂ T σ(t) = Keff · σ(t) + σ(t) · Keff +i Πα · QT − h.c. . ∂t α

(3)

Here, we have defined the phonon current matrices Πα (t) as Z

t

Πα (t) =

dt0 (G > (t, t0 )Sα< (t0 , t) − G < (t, t0 )Sα> (t0 , t)) .

(4)

t0

Consequently, the heat current from reservoir α into the system can be conveniently expressed in terms of Πα as follows: i T Jα (t) = − Tr Keff · Q · Πα (t) · QT − h.c. . 2 5


(5)


Page 6 of 21

The current matrices contain memory effects through the self-energies Sα . Since the thermal baths are harmonic, the self-energies Sα (t, t0 ) can be explicitly calculated Sα≶ (t, t0 )

Z = −i 0

∞

dω ω 0 0 coth cos ω(t − t ) ± i sin ω(t − t ) Lα (ω) . π 2kB Tα

Here, Tα is the bath temperature and Lα (ω) is the bath spectral density of reservoir α, which characterizes the action of the enviroment on the open system. 57,58 Together, a complicated set of coupled integro-differential equations for realistic systems is obtained (Equations (3) and (4)). To proceed further, we use an auxiliary-mode representation of the self-energies in terms of exponential functions, which has previously been developed for electrons 41–46 and for vibrations 59–61 (details are provided in the Supporting Information). The bath spectral density is given by the specifc form of the coupling to the heat reservoirs and can be computed using a recursive GF scheme. 57–59 Following the auxiliary-mode approach, we consider a spectral density with a Drude regularization (i.e., adding a cut-off frequency ωc ) (0)

(0)

(0)

of the form Lα (ω) = (ωc2 ω)/(ω 2 + ωc2 )Lα . 57 Here, Lα ≡ diag(Λα , 0) contains the coupling (0)

strength matrix Λα of the central region to the leads (see Section 2 in the Supporting Information for details). This spectral density contains the wide-band limit as a special case for ωc → ∞. Moreover, using a linear combination of several Lorentzians, any spectral density can in principle be approximated. In this way, the self-energies can be written as Sα (t, t0 ) =

∂ ∂t

PNP

p=0

0

(0)

−bα,p (t−t ) a Lα , α,p e

where the coefficients a α,p and bα,p are related to the auxiliary-mode decomposition (see Equations (S10) and (S11) in the Supporting Information). The phonon current matrices Πα (t) can then be expressed as

Πα (t) =

NP X

(0) Φpα (t) + a∗,< , α,p Q · Lα

(6)

p=0

where NP is the number of poles used in the auxiliary-mode expansion. The functions Φpα (t) (see Equation (S13) in the Supporting Information) are also obtained through the solution 6


Page 7 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of their EOMs, ∂ p p Φ (t) = Apα · Φpα (t) − Bαp · KT · L(0) α + Q · Ωα (t) , ∂t α p where Apα = K−b∗α,p I and Bαp = iωc δp,0 σ(t)+a∗,< α,p Q. The functions Ωα (t) =

(7) P PNP α0

p0 =0

0

Ωpα0pα (t)

correlate the main features of both heat baths and their values are obtained by performing 0

the time derivative of Ωpα0pα (t), 0 ∂ p0 p (0) p0 p p0 p Ωα0 α (t) = Cαp0pα Lα0 · Q · KT · L(0) α − Dα0 α Ωα0 α (t) ∂t 0 (0) − ωc δp0 ,0 Lα0 · K · Φpα (t) − δp,0 Φpα0† (t) · KT · L(0) . α

0

(8)

0

pp > ∗,> ∗,< ∗ Cαp0pα = a< α0 ,p0 aα,p − aα0 ,p0 aα,p and Dα0 α = bα0 ,p0 + bα,p are real numbers. Thus, the original

integro-differential equation for the reduced density matrix is mapped onto a closed set of ordinary differential equations (Equations (3), (7), and (8)), which are solved in this study by a fourth-order Runge-Kutta method with a timestep of 0.1 fs.

RESULTS AND DISCUSSION One dimensional chain To illustrate this general method, we first consider, for the central region, an one-dimensional (1D) linear chain consisting of N atoms coupled by force constants with magnitude λ = 1.0 eV/µÅ2 . In this case the coupling matrix to the reservoirs has only one non-zero entry (0)

(0)

η which characterizes the coupling strength, i.e., (ΛL )N N = η and (ΛR )11 = η. After extensive benchmarking of the parameters related to the Drude spectral density (see Section 2 in the Supporting Information), we have chosen the coupling parameter η = λ/2 and ωc = 1000 THz for both baths. These parameters provide an accurate description of the phonon dynamics in this 1D model and the heat current values in the steady state are in

7



good agreement with those obtained by the Landauer approach, as shown below. In Figure 2, we show the total energy after thermal equilibration (T0 = TL = TR ) as a function of the number of poles used in the auxiliary mode expansion for the case of a dimer (N = 2). Although the heat current tends to zero for each NP (see the inset in Figure 2), the system energy in the steady state still grows as NP increases. It converges for a given NP value depending on the temperature, e.g. 16 and 8 poles are needed to reach convergence at 50 K and 300 K, respectively. This is related to the fact that only few poles are sufficient to describe high temperature effects in the auxiliary-mode approach. 59 The final energy values are close to those obtained for an ideal harmonic oscillator in thermal equilibrium, P EC = N i=1 ωi (n(ωi ) + 1/2), where n(ω) is the Bose-Einstein distribution function and ωi are the frequencies of the isolated central system. Thus, we expect to find EC = 8.83 × 10−2 eV and EC = 9.56 × 10−2 eV at 50 K and 300 K, respectively. The discrepancies observed in Figure 2 relate to the coupling to the baths, i.e., they vanish in the weak coupling limit η → 0, as it is shown in Figure S1 of the Supporting Information. To gain further insight, we perform an eigenvalue decomposition of the quantity Z(t) = T · Q · σ(t), which is related to the total energy of the central region after tracing out (1/2)Keff

the bath degrees of freedom. To this end, we solve the eigenvalue problem for the N × N matrix resulting from adding the diagonal blocks of Z(t). In Figure 3a, we present the time evolution of the two eigenvalues for the diatomic chain at T0 = 300 K. Both eigenvalues increase as a function of time and reach their time-independent limit once the system reaches equilibrium. Moreover, by analyzing the corresponding eigenvectors, one finds that the lower mode (E1 ) corresponds to the case when the atoms vibrate with the same amplitude and in-phase (acoustic-like mode). This mode corresponds to the zero frequency translational mode of the isolated dimer; since translational symmetry is now broken due to the coupling to the baths, the mode acquires a finite frequency. The other mode (E2 ) corresponds to the stretching mode of the dimer, where the atoms vibrate out of phase. Once the system has reached thermal equilibrium at time te , we apply a symmetric 8


Page 8 of 21

Page 9 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


temperature bias ∆T = TL − TR = 2ξT0 (ξ > 0) with T0 corresponding to the mean temperature at which the system was previously equilibrated. The temperatures for the left and right leads can be written as TL = (1 + ξ)T0 and TR = (1 − ξ)T0 , respectively. Here, we define t = t − te as the simulation time after the temperature bias is switched on. In the insets of Figure 3a, we show the time evolution of the eigenvalues associated to the acoustic and optical modes of the diatomic chain for two different values of the temperature bias, ∆T = 60 K and 180 K. As ∆T increases, the fluctuations during the transient get larger because the central region exchanges a larger amount of energy with the heat baths in order to stabilize the entire system. As a function of time, the energy E1 first decreases and then increases until reaching the steady state, while the optical mode E2 is strongly affected during the transient, and it displays larger fluctuations than E1 . The time evolution of the heat current in both leads for linear chains containing two and eight atoms is presented in Figure 3b. In both cases, the transients show similar features: the heat current grows, oscillates and then saturates. One can also see that the time for reaching the steady state is longer for larger systems. The heat current in the steady state for several chain sizes at 300 K and for different values of the temperature bias is shown in Figure 3c. Independently of the temperature bias, the heat current increases with the number of sites, but it saturates for N > 8. The current values in the steady state are in agreement to those calculated by using the standard Landauer approach (see dashed lines in Figure 3c).

Poly-acetylene and poly-ethylene dimers In the next step, we demonstrate that our approach can be efficiently combined with firstprinciples methods to address phonon dynamics in realistic nanoscale systems. As paradigmatic examples we focus on poly-acetylene (PA, 4 atoms) and poly-ethylene (PE, 6 atoms) dimers connected to two harmonic baths, as shown in Figure 4a,b. The main difference between them is the presence of double carbon bonds in the PA junction, while the PE junction only contains single carbon bonds. Structural optimizations were performed by us9



ing density functional theory as implemented in the Gaussian09 code. 62 B3LYP functionals were used together with the 6-31G basis set. Then, a finite difference method has been used to get the corresponding force constant matrices, which are the main input in our method. (0)

Thus, Λα takes force constant values corresponding to realistic bonds connecting the end of the central region to the reservoirs (see Section 3 in the Supporting Information). For both molecular junctions, a Drude cut-off frequency of wc = 100 THz was used. This value is around twice the maximum frequency corresponding to their vibrational spectrum. The number of poles in the auxiliary-mode expansion at 100 K, 300 K, and 500 K are 10, 8, and 4, correspondingly. The variation of the energy density D(E, t) during the thermal equilibration at T0 = 300 K for PA and PE dimers is shown in Figures 4a and b, respectively. This density has been √ 2 √ P calculated as D(E, t) = N i=1 (1/γ 2π) exp[− (E − Ei )/γ 2 ] with γ = 0.001 eV as the width of the Gaussian distribution and {Ei } as the set of eigenvalues of Z(t). For both systems, all eigenmodes of the Z(t) matrix display very low energy at the beginning of the transient. Then, the lowest lying modes gain energy and achieve a maximum once the system reaches equilibrium. However, the eigenvalues in PE need a longer time to converge comparing to PA. This difference arises from the different coupling strengths to the leads (0)

(related to the matrices Λα ). The coupling of the PA dimer takes place through double bonds which are stiffer than the coupling produced by single bonds (PE dimer). The values of the coupling matrix are given in Section 3 in the Supporting Information. As a consequence, the magnitude of the oscillations in the heat current during the transient after applying a temperature bias is different (larger for PA as shown in Figure 4c). This difference in covalently bonded configurations also results in a larger heat current for the PA dimer. This is better seen by comparing the heat current as a function of ξ for both molecular junctions at T0 = 300 K, see Figure 4d. Here, we have also found that the heat current values for PE are very close to those found for the 1D atomic chain model, which is expected because of the similarity in configuration. Hydrogen atoms do not play an important role in the 10


Page 10 of 21

Page 11 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


thermal transport of these molecular junctions. The heat current displays a nearly linear dependence with respect to ξ for different mean temperatures T0 as one can see in Figure 4d in accordance with the linear response regime.

CONCLUSIONS In summary, we have developed a novel atomistic method which combines a TD-NEGF formalism with first-principle based modeling to study phonon dynamics in nanoscale junctions. The main idea of the method is to solve the equation of motion of the phonon density matrix with the help of an efficient auxiliary-mode approach. We have first applied the method to study thermal transport properties in the transient regime of an 1D atomic chain. The results for the steady states are in agreement with those obtained by using the Landauer approach. By computing the force constant matrix and coupling matrices using density-functional theory, we have been able to describe the phonon dynamics of molecular junctions consisting of poly-acetylene and poly-ethylene dimers. Although our examples were based on the Drude regularization of the spectral density, more realistic scenarios can be easily investigated. Our computational scheme allows it to study a variety of topical questions related to nanophononics, such as heat pumping, on an atomistic basis.

Acknowledgement L.M.S. thanks to the International Max Planck Research School Dynamical processes in atoms, molecules and solids and the Deutscher Akademischer Austauschdienst (DAAD) for the financial support. This work has also been partly supported by the German Research Foundation (DFG) within the Cluster of Excellence "Center for Advancing Electronics Dresden". The authors thank Prof. C. Lewenkopf and Dr. R. Biele for helpful discussions. We acknowledge the Center for Information Services and High Performance Computing (ZIH) at TU Dresden for providing computational resources. 11



Supporting Information Available The following files are available free of charge. • supp-info.pdf: detailed derivation of the equations of motion and benchmarking of the spectral density parameters.

References (1) Akinwande, D.; Petrone, N.; Hone, J. Two-Dimensional Flexible Nanoelectronics. Nat. Commun. 2014, 5, 5678. (2) Lu, W.; Lieber, C. M. Nanoelectronics from the Bottom-Up. Nat. Mater. 2007, 6, 841. (3) Žutić, I.; Fabian, J.; Das Sarma, S. Spintronics: Fundamentals and Applications. Rev. Mod. Phys. 2004, 76, 323–410. (4) Bogani, L.; Wernsdorfer, W. Molecular Spintronics Using Single-Molecule Magnets. Nat. Mater. 2008, 7, 179. (5) Li, N.; Ren, J.; Wang, L.; Zhang, G.; Hänggi, P.; Li, B. Colloquium: Phononics: Manipulating Heat Flow with Electronic Analogs and Beyond. Rev. Mod. Phys. 2012, 84, 1045–1066. (6) Volz, S.; Ordonez-Miranda, J.; Shchepetov, A.; Prunnila, M.; Ahopelto, J.; Pezeril, T.; Vaudel, G.; Gusev, V.; Ruello, P.; Weig, E. M. et al. Nanophononics: State of the Art and Perspectives. Eur. Phys. J. B 2016, 89, 15. (7) Balandin, A. A.; Nika, D. L. Phononics in Low-Dimensional Materials. Mater. Today 2012, 15, 266 – 275. (8) Maldovan, M. Sound and Heat Revolutions in Phononics. Nature 2013, 503, 209–217.

12


Page 12 of 21

Page 13 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(9) Cui, L.; Jeong, W.; Hur, S.; Matt, M.; Klöckner, J. C.; Pauly, F.; Nielaba, P.; Cuevas, J. C.; Meyhofer, E.; Reddy, P. Quantized Thermal Transport in Single-Atom Junctions. Science 2017, 355, 1192–1195. (10) Mosso, N.; Drechsler, U.; Menges, F.; Nirmalraj, P.; Karg, S.; Riel, H.; Gotsmann, B. Heat Transport Through Atomic Contacts. Nat. Nanotech. 2017, 12, 430–433. (11) Ozpineci, A.; Ciraci, S. Quantum Effects of Thermal Conductance Through Atomic Chains. Phys. Rev. B 2001, 63, 125415. (12) Yamamoto, T.; Watanabe, K. Nonequilibrium Green’s Function Approach to Phonon Transport in Defective Carbon Nanotubes. Phys. Rev. Lett. 2006, 96, 255503. (13) Dhar, A. Heat Transport in Low-Dimensional Systems. Adv. Phys. 2008, 57, 457–537. (14) Zhang, W.; Fisher, T. S.; Mingo, N. The Atomistic Green’s Function Method: An Efficient Simulation Approach for Nanoscale Phonon Transport. Numer. Heat Tr. BFund 2007, 51, 333–349. (15) Wang, J.-S.; Agarwalla, B. K.; Li, H.; Thingna, J. Nonequilibrium Green’s Function Method for Quantum Thermal Transport. Front. Phys. 2014, 9, 673–697. (16) Mingo, N. Anharmonic Phonon Flow Through Molecular-Sized Junctions. Phys. Rev. B 2006, 74, 125402. (17) Galperin, M.; Nitzan, A.; Ratner, M. A. Heat Conduction in Molecular Transport Junctions. Phys. Rev. B 2007, 75, 155312. (18) Wang, J.-S.; Zeng, N.; Wang, J.; Gan, C. K. Nonequilibrium Green’s Function Method for Thermal Transport in Junctions. Phys. Rev. E 2007, 75, 061128. (19) Li, Q.; Duchemin, I.; Xiong, S.; Solomon, G. C.; Donadio, D. Mechanical Tuning of Thermal Transport in a Molecular Junction. J. Phys. Chem. C 2015, 119, 24636–24642. 13



(20) Bürkle, M.; Hellmuth, T. J.; Pauly, F.; Asai, Y. First-Principles Calculation of the Thermoelectric Figure of Merit for [2,2]Paracyclophane-Based Single-Molecule Junctions. Phys. Rev. B 2015, 91, 165419. (21) Medrano Sandonas, L.; Gutierrez, R.; Pecchia, A.; Seifert, G.; Cuniberti, G. Tuning Quantum Electron and Phonon Transport in Two-Dimensional Materials by Strain Engineering: A Green’s Function Based Study. Phys. Chem. Chem. Phys. 2017, 19, 1487–1495. (22) Medrano Sandonas, L.; Sevincli, H.; Gutierrez, R.; Cuniberti, G. First-Principle-Based Phonon Transport Properties of Nanoscale Graphene Grain Boundaries. Adv. Sci. 2018, 5, 1700365. (23) Klöckner, J. C.; Bürkle, M.; Cuevas, J. C.; Pauly, F. Length Dependence of the Thermal Conductance of Alkane-Based Single-Molecule Junctions: An ab-initio Study. Phys. Rev. B 2016, 94, 205425. (24) Sevinçli, H.; Sevik, C.; Çagin, T.; Cuniberti, G. A Bottom-Up Route To Enhance Thermoelectric Figures of Merit in Graphene Nanoribbons. Sci. Rep. 2013, 3, 1–5. (25) Kim, P.; Shi, L.; Majumdar, A.; McEuen, P. L. Thermal Transport Measurements of Individual Multiwalled Nanotubes. Phys. Rev. Lett. 2001, 87, 215502. (26) Park, W.; Romano, G.; Ahn, E. C.; Kodama, T.; Park, J.; Barako, M. T.; Sohn, J.; Kim, S. J.; Cho, J.; Marconnet, A. M. et al. Phonon Conduction in Silicon Nanobeam Labyrinths. Sci. Rep. 2017, 7, 6233. (27) Seol, J. H.; Jo, I.; Moore, A. L.; Lindsay, L.; Aitken, Z. H.; Pettes, M. T.; Li, X.; Yao, Z.; Huang, R.; Broido, D. et al. Two-Dimensional Phonon Transport in Supported Graphene. Science 2010, 328, 213–216.

14


Page 14 of 21

Page 15 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(28) Wagner, M. R.; Graczykowski, B.; Reparaz, J. S.; El Sachat, A.; Sledzinska, M.; Alzina, F.; Sotomayor Torres, C. M. Two-Dimensional Phononic Crystals: Disorder Matters. Nano Lett. 2016, 16, 5661–5668. (29) Anufriev, R.; Ramiere, A.; Maire, J.; Nomura, M. Heat Guiding and Focusing Using Ballistic Phonon Transport in Phononic Nanostructures. Nat. Commun. 2017, 8, 15505. (30) Maire, J.; Anufriev, R.; Yanagisawa, R.; Ramiere, A.; Volz, S.; Nomura, M. Heat Conduction Tuning by Wave Nature of Phonons. Sci. Adv. 2017, 3 . (31) Segal, D.; Nitzan, A. Molecular Heat Pump. Phys. Rev. E 2006, 73, 026109. (32) Ren, J.; Hänggi, P.; Li, B. Berry-Phase-Induced Heat Pumping and Its Impact on the Fluctuation Theorem. Phys. Rev. Lett. 2010, 104, 170601. (33) Galperin, M.; Saito, K.; Balatsky, A. V.; Nitzan, A. Cooling Mechanisms in Molecular Conduction Junctions. Phys. Rev. B 2009, 80, 115427. (34) Segal, D.; Nitzan, A. Heat Rectification in Molecular Junctions. J. Chem. Phys. 2005, 122, 194704. (35) Segal, D. Heat Flow in Nonlinear Molecular Junctions: Master Equation Analysis. Phys. Rev. B 2006, 73, 205415. (36) Díaz, E.; Gutierrez, R.; Cuniberti, G. Heat Transport and Thermal Rectification in Molecular Junctions: A Minimal Model Approach. Phys. Rev. B 2011, 84, 144302. (37) Wei, D.; Obstbaum, M.; Ribow, M.; Back, C. H.; Woltersdorf, G. Spin Hall Voltages from A.C. and D.C. Spin Currents. Nat. Commun. 2014, 5, 3768. (38) Weiler, M.; Shaw, J. M.; Nembach, H. T.; Silva, T. J. Phase-Sensitive Detection of Spin Pumping via the AC Inverse Spin Hall Effect. Phys. Rev. Lett. 2014, 113, 157204.

15



(39) Bocklage, L. Coherent THz Transient Spin Currents by Spin Pumping. Phys. Rev. Lett. 2017, 118, 257202. (40) Gaury, B.; Weston, J.; Santin, M.; Houzet, M.; Groth, C.; Waintal, X. Numerical Simulations of Time-Resolved Quantum Electronics. Phys. Rep. 2014, 534, 1 – 37. (41) Croy, A.; Saalmann, U. Propagation Scheme for Nonequilibrium Dynamics of Electron Transport in Nanoscale Devices. Phys. Rev. B 2009, 80, 245311. (42) Popescu, B.; Woiczikowski, P. B.; Elstner, M.; Kleinekathöfer, U. Time-Dependent View of Sequential Transport through Molecules with Rapidly Fluctuating Bridges. Phys. Rev. Lett. 2012, 109, 176802. (43) Popescu, B. S.; Croy, A. Efficient Auxiliary-Mode Approach for Time-Dependent Nanoelectronics. New J. Phys. 2016, 18, 093044. (44) Popescu, B. S.; Croy, A. Emergence of Bloch Oscillations in One-Dimensional Systems. Phys. Rev. B 2017, 95, 235433. (45) Leitherer, S.; Jäger, C. M.; Krause, A.; Halik, M.; Clark, T.; Thoss, M. Simulation of Charge Transport in Organic Semiconductors: A Time-Dependent Multiscale Method Based on Nonequilibrium Green’s Functions. Phys. Rev. Mater. 2017, 1, 064601. (46) Lehmann, T.; Croy, A.; Gutierrez, R.; Cuniberti, G. Time-Dependent Framework for Energy and Charge Currents in Nanoscale Systems. Chem. Phys. 2018, https://doi.org/10.1016/j.chemphys.2018.01.011. (47) Xie, H.; Jiang, F.; Tian, H.; Zheng, X.; Kwok, Y.; Chen, S.; Yam, C.; Yan, Y.; Chen, G. Time-Dependent Quantum Transport: An Efficient Method Based on Liouville-VonNeumann Equation for Single-Electron Density Matrix. J. Chem. Phys. 2012, 137, 044113.

16


Page 16 of 21

Page 17 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(48) Morzan, U. N.; Ramirez, F. F.; Lebrero, M. C. G.; Scherlis, D. A. Electron Transport in Real Time from First-Principles. J. Chem. Phys. 2017, 146, 044110. (49) Kurth, S.; Stefanucci, G.; Almbladh, C.-O.; Rubio, A.; Gross, E. K. U. Time-Dependent Quantum Transport: A Practical Scheme Using Density Functional Theory. Phys. Rev. B 2005, 72, 035308. (50) Zhang, L.; Xing, Y.; Wang, J. First-Principles Investigation of Transient Dynamics of Molecular Devices. Phys. Rev. B 2012, 86, 155438. (51) Zheng, X.; Wang, F.; Yam, C. Y.; Mo, Y.; Chen, G. Time-Dependent DensityFunctional Theory for Open Systems. Phys. Rev. B 2007, 75, 195127. (52) Yam, C.; Zheng, X.; Chen, G.; Wang, Y.; Frauenheim, T.; Niehaus, T. A. TimeDependent Versus Static Quantum Transport Simulations Beyond Linear Response. Phys. Rev. B 2011, 83, 245448. (53) Oppenländer, C.; Korff, B.; Niehaus, T. A. Higher Harmonics and AC Transport from Time Dependent Density Functional Theory. J. Comput. Electron. 2013, 12, 420–427. (54) Biele, R.; D’Agosta, R.; Rubio, A. Time-Dependent Thermal Transport Theory. Phys. Rev. Lett. 2015, 115, 056801. (55) Sena-Junior, M. I.; Lima, L. R. F.; Lewenkopf, C. H. Phononic Heat Transport in Nanomechanical Structures: Steady-State and Pumping. J. Phys. A-Math. Theor. 2017, 50, 435202. (56) Tuovinen, R.; Säkkinen, N.; Karlsson, D.; Stefanucci, G.; van Leeuwen, R. Phononic Heat Transport in the Transient Regime: An Analytic Solution. Phys. Rev. B 2016, 93, 214301. (57) Ulrich, W. Quantum Dissipative Systems (Second Edition); Series in Modern Condensed Matter Physics; World Scientific Publishing Company, 1999. 17



(58) de Vega, I.; Alonso, D. Dynamics of Non-Markovian Open Quantum Systems. Rev. Mod. Phys. 2017, 89, 015001. (59) Ritschel, G.; Eisfeld, A. Analytic Representations of Bath Correlation Functions for Ohmic and Superohmic Spectral Densities Using Simple Poles. J. Chem. Phys. 2014, 141, 094101. (60) Kleinekathöfer, U. Non-Markovian Theories Based on a Decomposition of the Spectral Density. J. Chem. Phys. 2004, 121, 2505–2514. (61) Meier, C.; Tannor, D. J. Non-Markovian Evolution of the Density Operator in the Presence of Strong Laser Fields. J. Chem. Phys. 1999, 111, 3365. (62) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian−09 Revision E.01. 2009.

18


Page 18 of 21

Page 19 of 21

Right lead

Left lead TL

VL

TR

VR

Molecule K

Figure 1: Schematic representation of the target molecular junctions in the present study. A molecular system is connected to two harmonic thermal baths, which drive the heat flow in the molecule. 0.12

T0 = 50 K T0 = 100 K

0.10 J [eV/ps]

EC [eV]

T0 = 300 K

0.08

0

NP = 2 NP = 4 NP = 16

-0.8 -1.6 0.001 0.01

4

0.1

t [ps] 12 16

8

1

20

NP

E2

0.8

6.4

E1

0.4

-2

t [ps]

2.0

4.6

t [ps]

0.001

0.01

t [ps]

(a)

0.1

Right lead

0.0 Left lead

-0.4

4.4 0.01

0.0

0.2

1

1

J [eV/ps]

1

J [eV/ps]

0.01

E [x10 eV]

-2

4.0

6.5

ξ = 0.1 ξ = 0.3 ξ = 0.5

-2

6.0

E [X10 eV]

Figure 2: Energy of a diatomic chain vs number of poles NP of the auxiliary-mode approach at different temperatures. The arrows indicates the convergence point for each temperature. Inset: Dynamics of the heat current for the diatomic chain at T0 = 50 K.

E [x10 eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.1

N=2 N=8

-0.8 0.001

0.01

0.1

1

t [ps]

(b)

0.0

4

8

N

12

16

(c)

Figure 3: (a) Time evolution of the eigenvalues Ei for a diatomic chain at T0 = 300 K. Inset: Their corresponding transient states after applying symmetric temperature bias around T0 (TL = (1 + ξ)T0 and TR = (1 − ξ)T0 ), ξ = 0.1 (∆T = 60 K, solid lines) and ξ = 0.3 (180 K, dashed lines). (b) Dynamics of the heat current from the left and right leads for atomic chains with 2 and 8 atoms with ξ = 0.1. (c) Variation of the steady heat current after increasing the number of atoms in the atomic chain for different ξ. The mean temperature is 300 K. The dashed lines represent the values obtained by using Landauer approach. 19



Poly-acetylene (PA)

Poly-ethylene (PE)

(a)

(b) 1.2

0.8 Right lead

J [eV/ps]

0.4

J [eV/ps]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 21

0.0 -0.4 PA PE

-0.8 0.001

0.01

0.1

Left lead

1

0.8

0.4

0.0 0

10

t [ps] (c)

PA, T0 = 100 K PA, T0 = 300 K PA, T0 = 500 K PE, T0 = 300 K

0.2

0.4

0.6

ξ

(d)

Figure 4: Energy density plot for molecular junctions made of (a) poly-acetylene (PA) and (b) poly-ethylene (PE) dimers. (c) Dynamics of the heat current for both leads for PA and PE dimers after applying a temperature bias of ∆T = 60 K at T0 = 300 K. (d) Variation of the steady-state heat current as a function of ξ for the PA dimer at different T0 . We also show the results for the corresponding PE dimer at T0 = 300 K (dashed lines).

20


Page 21 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

21


An Atomistic Framework for Time-Dependent Thermal Transport - The

Recommend Documents