Thermal transport in a 2D stressed nanostructure with mass gradient

Inspired by some recent molecular dynamics (MD) simulations and experiments on suspended graphene nanoribbons, we study a simplified model where the atoms are disposed in a rectangular lattice coupled by nearest neighbor interactions which are quadratic in the interatomic distance. The system has a mechanical strain, and the border atoms are coupled to Langevin thermal baths. Atom masses vary linearly in the longitudinal direction, modeling an isotope or doping distribution. This asymmetry and tension modify thermal properties. Although the atomic interaction is quadratic, the potential is anharmonic in the coordinates. By direct MD simulations and solving Fokker--Planck equations at low temperatures, we can better understand the role of anharmonicities in thermal rectification. We observe an increasing thermal current with an increasing applied mechanical tension. The temperatures and thermal currents vary along the transverse direction. This effect can be useful to establish which parts of the system are more sensitive to thermal damage. We also study thermal rectification as a function of strain and system size.


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, Fouriers 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) [1][2][3]. 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 [4,5] 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 [6]. Thus, systems with mass gradients due to dopping concentrations, deposition of heavy molecules or a variable width are candidates to present this phenomena [1,[7][8][9].
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 [10][11][12][13].
From this motivation, we present a simple model for a 2D system to better understand heat trans-port 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 N = N x × N y particles. They are identified by an index i = (i x , i y ). We consider a linear mass gradient along the x direction, with , being M L and M R 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 r = |R i − R j |, as: The particles on the left (i x = 1) and right (i x = N x ) borders also interact by the same potential with two substracts that can be thought as a left (i x = 0) and right (i x = N x + 1) columns of fixed particles. Therefore, there are (N x + 1)N y interactions or bonds in the x-direction, and N x (N y − 1) bonds in the y-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 (a x , a y ) = (l 0 , l 0 ). However, if the distance L x between the fixed rows is bigger than (N x + 1)l 0 , there will be a tension in the system along the xdirection that we parametrize by a change in the lattice constant a x > l 0 . 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 (i y = 1) and top (i y = N y ) rows do not interact with any external substrate or reservoir, therefore there could not be any tension in the y-direction, being a y = l 0 for all cases. The equilibrium position of the particles is R 0i = (a x i x , a y (i y − 1)). We characterize the motion of the particles by the displacement with respect to their equilibrium position r i = R i −R 0i . 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 γ i = 1 for i x = 1 or i x = N x , and zero otherwise. The random forces have the correlations where T i is T L and T R for i x = 1 and i x = N x , the temperatures of the left and right reservoirs, respectively, or zero otherwise. Indices µ and ν run over x and y 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 i and j, energy flows at a rate given by J i,j = F i,j · (v i + v j ) , where F i,j is the force that i does on j. 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 i over time Q i (t), 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 Q i (t)Q j (t) : the motion of particle i will be in general correlated with the motion of particle j. 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 (1) for two neighboring particles, it has the form v(r) = 1 2 k(r 2 − 2rl 0 + l 2 0 ). The term on r 2 is fully quadratic on the components of r i and r j . The problem arises with the second term proportional to r, where, for simplification, a is the vector joining the equilibrium positions of the particles, and ∆r the relative displacement between them. We see that this term is not linear on the components of r i or r j r = (a + ∆r) · (a + ∆r) due to the square root. In the last equality, we have used a = |a|,â = a/a, and δ = ∆r/a. Taking into account that the proposed model is valid for small displacements, we can make an expansion considering |∆r| a, or equivalently |δ| 1. 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 x and y direction, after collecting terms, we arrive to a quadratic approximation for the total potential energy of the system where , and α iy = 1/2 for i y = 1, N y or 1 otherwhise. In this harmonic approximation, the directions x and y are completely decoupled. The dependence on the tension comes through the transversal spring constant k ⊥ .
For every particle, we have two degrees of freedom, so the system has M = 2N total number of degrees of freedom. If we call this coordinates as {q n }, we can rewrite the equations of motion aṡ where p n = m nqn 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 {f n (t)}. We explain in more detail the general solution following [14].
By defining 2M -dimensional vectors X = (q 1 , . . . , q M ; p 1 , . . . , p M ) and F(t) = (0, . . . , 0; f 1 (t), . . . , f M (t)) we can further simplify the notationẊ The matrix A has the structure with K the force matrix, and the diagonal matrices M nm = δ nm m n and G nm = δ nm γ n /m n . This matrix A can be diagonalized, obtaining 2M complex eigenvalues λ i 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 U, which fulfills AU = UA , where A is a diagonal matrix with the eigenvalues as elements. Transforming and defining X = U −1 X, and F (t) = U −1 F(t), the matrix Eq. (9) transforms toẊ Now the FP equations are decoupled, although each component of F (t) is a linear combination of all components f n (t).
For a given realization of the random forces, and some initial condition, each equation can be integrated, obtaining Taking the solutions for each x i , and using that X(t) = UX (t), we finally obtain the solutions for positions q n and momenta p n . 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 D kl = 2γ k k B T k δ kl , where T k corresponds to the temperature for the moment component which are coupled to a thermal bath, or zero otherwise. Then, transforming this matrix as D = U −1 D (U −1 ) T , 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 y 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.

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 nonlinear 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 x 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 y 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.
Averaging the kinetic energy of each particle, we compute temperature profiles, seeing an example in Fig. 5. There is an expected temperature decay in x direction from the hot to the cold bath, although it is not uniform. There are also some fluctuations in the transversal y 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.
The heat currents are computed for every bond. For bonds in y direction, average currents are of the same order of statistical errors, and we can induce that they are almost zero, which is the result 070008-5 given by FP equations. On x direction, the average current is the same for every bond in a given i y 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 J, and dividing by N y . We show on Fig. 6 this quantity as a function of the temperature bias ∆T = T L − T R , for a fixed value of an average bath temperature T 0 = (T L + T R )/2. We show both results for positive and negative bias, which produce positive and negative heat currents. In both cases, the behavior is linear on ∆T (allowing us to later define a conductance). However, slopes are different depending the sign of ∆T . 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.
To study the heat current and thermal rectification as a function of strain a x , we plot in Fig. 7 the conductance per row C = J/(N y ∆T ), 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.
As a funtion of system width N y , 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 N x . 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 (N x → ∞) cannot be inferred from these simulations; and answering this question is beyond the scope of the

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 k r = 652 nN/nm and equilibrium distance a 0 =0.1421 nm [12]. Taking into account the carbon mass m 0 , we obtain time and temperature scales τ 0 = m 0 /k r = 5.5 fs , and T 0 = k r a 2 0 /k B ≈ 10 6 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 T ≈ 3 · 10 −4 . From some experimental works [13], it is known that suspended graphene could be stable up to 2600 K. In our model, this temperature corresponds to T ≈ 2.7 · 10 −3 , 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 2.5 · 10 −3 W. We observed an increasing conductance with uniaxial strain, contrary to MD computations for AGNR (armchair) graphene [11]. 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.