Thermal transport in a 2D stressed nanostructure with mass gradient

Thermal transport in a 2D stressed nanostructure with mass gradient

R. Barreto,[1, 2]   M. F. Carusela,[1, 2]   A. Mancardo Viotti,[1]   A. G. Monastra[1, 2] E-mail:


    • [1] Instituto de Ciencias, Universidad Nacional de Gral. Sarmiento, J. M. Gutiérrez 1150, 1613 Los Polvorines, Argentina.
    • [2] Consejo Nacional de Investigaciones Científicas y Técnicas, Argentina.

I. Introduction and motivation

Efficient energy consumption is one of the biggest challenges for modern societies. Moreover, miniaturization of electronic devices together with their increasing computing capabilities make heat production per surface/volume unit tend to increase constantly. New technologies are necessary for efficient heat transport. At nanometric scales and for low-dimensional systems, Fourier´s law is not fullfilled and new phenomena arise. Among them, thermal rectification would make possible to build thermal diodes, being low dimensional systems the best candidates (e.g., atomic chains, graphene) [3, 4, 5]. In particular, graphene nanoribbons (GNR) are promising candidates in nanoelectronics. However, at the chip-level integrated circuits, the power density highly increased, making the thermal managment vital to ensure a stable operation of any practical graphene-based device [6, 7]

Thermal conductance can be modified by defects, impurities, shapes, geometries, mechanical strains, asymmetries, etc. It is known that it is necessary to have an asymmetry on the system to achieve heat rectification [8]. Thus, systems with mass gradients due to dopping concentrations, deposition of heavy molecules or a variable width are candidates to present this phenomena [3, 9, 10, 11].

On the other hand, thermal conductance of 2D systems, as GNRs, is remarkably affected by tensile strain. Moreover, mechanical tension is relatively easy to control experimentally at nanoscale, being a good candidate to tune the phononic heat transport in a system [12, 13, 14, 15].

From this motivation, we present a simple model for a 2D system to better understand heat transport and the possibility of thermal rectification in a layered-device with variable widths subject to a mechanical longitudinal tension.

Ii. The model

The system is composed of a rectangular array of particles. They are identified by an index . We consider a linear mass gradient along the direction, with , being and the masses of the left and right columns respectively (see Fig. 1). The particles interact isotropically to nearest neighbors by potentials that only depend on their relative distance , as:

Figure 1: Schematic of the system. Particles bordered by dash-dot lines are coupled to Langevin thermal baths.

The particles on the left and right borders also interact by the same potential with two substracts that can be thought as a left and right columns of fixed particles. Therefore, there are interactions or bonds in the -direction, and bonds in the -direction that contribute to the total potential


that only depends on the positions.

The natural equilibrium position of the particles is a rectangular array with lattice constants . However, if the distance between the fixed rows is bigger than , there will be a tension in the system along the -direction that we parametrize by a change in the lattice constant . We are specially interested on the effects of tension on the thermal properties of the system because it is an external parameter that can be easily controlled. The bottom and top rows do not interact with any external substrate or reservoir, therefore there could not be any tension in the -direction, being for all cases. The equilibrium position of the particles is . We characterize the motion of the particles by the displacement with respect to their equilibrium position . Moreover, the particles on the left and right columns are coupled to two thermal reservoirs respectively. We consider a Langevin interaction by a viscous term proportional to velocity, and a decorrelated random force acting on the particles in contact with the reservoirs. The equation of motion for each particle is


where for or , and zero otherwise. The random forces have the correlations , where is and for and , the temperatures of the left and right reservoirs, respectively, or zero otherwise. Indices and run over and directions.

For a given realization of the random forces, equations of motion are integrated from some initial condition. After some time, position and velocity of the particles attain a stationary regime where their statistical behavior is constant. Averaging a given realization over time, we are mainly interested on the following quantities: i) Mean quadratic velocities: they define a kinetic site temperature , even if the system is not in thermodynamic equilibrium; ii) Heat current: for two interacting particles and , energy flows at a rate given by , where is the force that does on . The energy current could be different for each bond, being the only conserved quantity the total current through any transversal section of the system; iii) Time correlations: for any component of position or velocity of particle over time , we are interested on its autocorrelation . This function contains very useful information about characteristic time scales of the system, how long time averages should be done to be significant, and to estimate the statistical errors of averaged quantities; iv) Spatial correlations : the motion of particle will be in general correlated with the motion of particle . It is expected that neighbor particles will have stronger correlation than distant particles. A correlation lenght arises, depending on parameters of the system and thermal baths.

Iii. Fokker–Planck equations - Harmonic approximation

Expanding interatomic potential (II.) for two neighboring particles, it has the form . The term on is fully quadratic on the components of and . The problem arises with the second term proportional to ,

where, for simplification, is the vector joining the equilibrium positions of the particles, and the relative displacement between them. We see that this term is not linear on the components of or


due to the square root. In the last equality, we have used , , and . Taking into account that the proposed model is valid for small displacements, we can make an expansion considering , or equivalently . Collecting terms of the same order in , we arrive to


We see explicitly in this expansion that the potential energy is not quadratic on the coordinates. Nevertheless, for very low temperatures, where displacements are small with respect to the lattice constants, we can neglect the cubic and higher order terms. Making this approximation for all bonds in and direction, after collecting terms, we arrive to a quadratic approximation for the total potential energy of the system


where , , and for or 1 otherwhise. In this harmonic approximation, the directions and are completely decoupled. The dependence on the tension comes through the transversal spring constant .

For every particle, we have two degrees of freedom, so the system has total number of degrees of freedom. If we call this coordinates as , we can rewrite the equations of motion as


where are the momenta, and

are elements of a force matrix. This is a set of stochastic linear equations, the so called Fokker–Planck (FP) equations, that can be exactly integrated for a given realization of the random forces . We explain in more detail the general solution following [16].

By defining -dimensional vectors and we can further simplify the notation


The matrix has the structure


with the force matrix, and the diagonal matrices and . This matrix can be diagonalized, obtaining complex eigenvalues that come in complex conjugate pairs. Their real part is always positive. The eigenvectors are also complex, and arranging them as columns, we obtain the diagonalizing unitary matrix , which fulfills , where is a diagonal matrix with the eigenvalues as elements. Transforming and defining , and , the matrix Eq. (III.) transforms to

Now the FP equations are decoupled, although each component of is a linear combination of all components .

For a given realization of the random forces, and some initial condition, each equation can be integrated, obtaining

Taking the solutions for each , and using that , we finally obtain the solutions for positions and momenta . For sufficiently long times, the term proportional to the initial condition will vanish, provided that

Now we are interested in the statistical behavior of the system, particularly in the stationary regime. With these solutions, one can compute the mean values doing averages over the ensemble of random forces, which are necessary to estimate time and spatial correlations, currents per bond and site temperatures.

Taking into account the correlations of the random forces that depend mostly on bath temperatures, and after some calculations, the time and spatial correlations of all dynamical variables can be computed. First, we define the diagonal matrix , where corresponds to the temperature for the moment component which are coupled to a thermal bath, or zero otherwise. Then, transforming this matrix as , and defining a new matrix


the general correlation function is


We plot on Figs. 2 and 3 the time correlation functions computed by FP equations. We observe that particles coupled to a heat bath, or near to it, decorrelate in shorter times, as expected due to the random forces. Particles in the middle of the system have longer correlation times. The frequencies on these functions are related to normal modes weakly coupled to the heat baths. We observe that for times of the order of 200, most of these functions decay significantly, except for the velocity in direction. Anyway, the rapid and non periodic fluctuations make it difficult to predict on average the behavior for long times. From these time correlation functions, we conclude that the system can attain a stationary regime at times of the order of 500, and measurements of dynamical variables separated by this time scale can be considered independent.

Figure 2: Normalized time autocorrelation function for (up) and (down) positions of some selected particles. For all figures , , spring constant , and natural length . , , , . See final discussion for units.
Figure 3: Normalized time autocorrelation function for (up) and (down) velocities of some selected particles, same values as Fig. 2.

Iv. Numerical results

The full potential FP equations are difficult to solve. The only way to compute correlations, temperature profiles and heat currents in the stationary state is to perform a direct numerical integration of the equations of motion, as in molecular dynamics (MD) computations, using a Runge-Kutta stochastic integrator. At very low temperatures, where non-linear terms are negligible, FP results coincide with MD within the statistical error. At intermediate temperatures, we first compare the time correlation functions computed by both methods and we observe only minor differences. In general, the functions decay faster for MD, meaning that non-linear terms induce more decorrelation. Nevertheless, we can take the correlation times from FP results as a good upper estimation for all temperature regimes.

Spatial correlation functions for the displacements are shown in Fig. 4. The position of a given particle is correlated with other particles on the same row, while it is almost decorralated with particles in a different row because the coupling is only through higher order terms in the potential (in FP calculations it is strictly zero). On the other hand, the position of a given particle is strongly correlated with all other particles. Therefore, the correlation length is of the order of the system size for the studied parameters. The spatial correlation of velocities in both directions strongly decays even for two neighboring particles.

Figure 4: Normalized spatial correlation function for (left) and (right) position of the particle. , , , .

Averaging the kinetic energy of each particle, we compute temperature profiles, seeing an example in Fig. 5. There is an expected temperature decay in direction from the hot to the cold bath, although it is not uniform. There are also some fluctuations in the transversal direction, even at very low temperatures. This is produced by the asymmetry of particles on the top and bottom rows, which are coupled to 3 neighbors, compared to particles on the bulk rows with 4 neighbors.

Figure 5: Kinetic temperature profile. , , , .

The heat currents are computed for every bond. For bonds in direction, average currents are of the same order of statistical errors, and we can induce that they are almost zero, which is the result given by FP equations. On direction, the average current is the same for every bond in a given row. On the other hand, the currents on each row are slightly different, depending on strain and other parameters. We finally compute an average heat current per row through the system summing up the currents of all rows , and dividing by . We show on Fig. 6 this quantity as a function of the temperature bias , for a fixed value of an average bath temperature . We show both results for positive and negative bias, which produce positive and negative heat currents. In both cases, the behavior is linear on (allowing us to later define a conductance). However, slopes are different depending the sign of . This comes from the mass gradient that stablishes an asymmetry on the system, and by the non-linear part of the potential. This effect is not observed in the harmonic approximation used to solve FP equations. The asymmetry in the heat transmition is called thermal rectification, and it could be a very useful property.

Figure 6: Total current per row as a function of . , , .

To study the heat current and thermal rectification as a function of strain , we plot in Fig. 7 the conductance per row , for both signs of the temperature bias. We observe an increasing current as strain is increased, although the thermal rectification seems to decrease as both curves approach.

Figure 7: Average conductance per row as a function of strain. , .

As a funtion of system width , the conductance quickly achieves a constant value, being compatible with Fourier law, where heat current is proportional to transversal area. In Fig. 8, we plot the conductivity and conductance per row as a function of system length . FP calculations give a decreasing conductance although it seems to achieve a constant value, breaking the classical Fourier law. Indeed, the conductivity diverges as it is expected for ordered harmonic crystals. For MD computations, we observe a smaller conductance with respect to FP, and it decreases faster with system size. We observe that thermal rectification increases with length, at least for this small system regime, which can be a very useful effect to build a thermal diode. Whether or not these tendencies continue in the thermodynamic limit () cannot be inferred from these simulations; and answering this question is beyond the scope of the present work.

Figure 8: Average conductivity (top) and conductance (bottom) per chain as a function of system length. Full (blue) line corresponds to Fokker–Planck calculation. , , .

V. Conclusion and discussion

We studied thermal properties of a 2D atomic model with mass gradient under tension. Although it is a simple theoretical model, we can make some relations with graphene, a real 2D atomic system, at least for some orders of magnitude. For thermal and mechanical properties, graphene models use the effective classical Tersoff-Brenner potential. The interaction along carbon-carbon direction can be approximated by a quadratic potential with constant = 652 nN/nm and equilibrium distance =0.1421 nm [14]. Taking into account the carbon mass , we obtain time and temperature scales = 5.5 fs , and K.

Autocorrelation functions in our model would suggest a time scale from 3 to 6 ps for loss of information, at least for systems of length up to 30 atoms. Characteristic frequencies can be related to those of normal modes weakly coupled to thermal baths.

Room temperature of around 300 K corresponds to dimensionless parameter . From some experimental works [15], it is known that suspended graphene could be stable up to 2600 K. In our model, this temperature corresponds to , where the mean particle displacement becomes of the order of lattice constant. In this case, quadratic approximation for the atomic interaction and nearest neighbor approximation would not be valid anymore. A typical thermal current per atomic chain, as in Fig. 6, would correspond to W. We observed an increasing conductance with uniaxial strain, contrary to MD computations for AGNR (armchair) graphene[13]. However, as in our model, an increase of conductance was obtained for ZGNR (zigzag) graphene nanoribbons at small strain and room temperature, due to an increased of the phonon velocity of some modes in the low and high frequencies for small strains.

We observed thermal rectification in a model with an interatomic potential quadratic on the distance, although as a two dimensional model, the potential is non-linear in the coordinates. Also the mass gradient, which makes the system asymmetric, is essential for the thermal rectification. The effect could be even stronger incorporating cubic and quartic terms in the interatomic potentials (Fermi–Pasta–Ulam models), directional terms (strongly present in carbon-carbon interaction), and flexural modes.


  • [3] N A Roberts, D G Walker, A review of thermal rectification observations and models in solid materials, Int. J. Thermal Sci. 50, 648 (2011).
  • [4] A Dhar, Heat transport in low-dimensional systems, Adv. in Phys. 57, 457 (2008).
  • [5] N Li, J Ren, L Wang, G Zhang, P Hänggi, B. Li, Colloquium: Phononics: Manipulating heat flow with electronic analogs and beyond, Rev. Mod. Phys. 84, 1045 (2012).
  • [6] Z Li, H Qian, J Wu, B-L Gu, W Duan, Role of symmetry in the transport properties of graphene nanoribbons under bias, Phys. Rev. Lett. 100, 206802 (2008).
  • [7] E Pop, S Sinha, K E Goodson, Heat generation and transport in nanometer-scale transistors, Proc. IEEE 94, 1587 (2006).
  • [8] S Das, A Dhar, O Narayan, Heat conduction in the Fermi–Pasta–Ulam chain, J. Stat. Phys. 154, 204 (2014).
  • [9] J Hu, X Ruan, Y P Chen, Thermal conductivity and thermal rectification in graphene nanoribbons: A molecular dynamics study, Nano Lett. 9, 2730 (2009).
  • [10] E Pereira, Graded anharmonic crystals as genuine thermal diodes: Analytical description of rectification and negative differential thermal resistance, Phys. Rev. E 82, 040101(R) (2010).
  • [11] C W Chang, D Okawa, A Majumdar, A Zettl, Solid-state thermal rectifier, Science 314, 1121 (2006).
  • [12] N Wei, L Xu, H-Q Wang, J-C Zheng, Strain engineering of thermal conductivity in graphene sheets and nanoribbons: a demonstration of magic flexibility, Nanotechnology 22, 105705 (2011).
  • [13] K G S H Gunawardana, K Mullen, J Hu, Y P Chen, X Ruan, Tunable thermal transport and thermal rectification in strained graphene nanoribbons, Phys. Rev. B 85, 245417 (2012).
  • [14] G I Giannopoulos, I A Liosatos, A K Moukanidis, Parametric study of elastic mechanical properties of graphene nanoribbons by a new structural mechanics approach, Physica E 44, 124 (2011).
  • [15] K Kim et al., High-temperature stability of suspended single-layer graphene, Phys. Status Solidi 4, 302 (2010).
  • [16] H Risken, The Fokker–Planck equation, Springer, New York (1989).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description