Nonequilibrium processes from generalized Langevin equations: Realistic nanoscale systems connected to two thermal baths

We extend the generalized Langevin equation (GLE) method [L. Stella, C. D. Lorenz, and L. Kantorovich, Phys. Rev. B 89, 134303 (2014)] to model a central classical region connected to two realistic thermal baths at two different temperatures. In such nonequilibrium conditions a heat flow is established, via the central system, in between the two baths. The GLE-2B (GLE two baths) scheme permits us to have a realistic description of both the dissipative central system and its surrounding baths. Following the original GLE approach, the extended Langevin dynamics scheme is modified to take into account two sets of auxiliary degrees of freedom corresponding to the mapping of the vibrational properties of each bath. These auxiliary variables are then used to solve the non-Markovian dissipative dynamics of the central region. The resulting algorithm is used to study a model of a short Al nanowire connected to two baths. The results of the simulations using the GLE-2B approach are compared to the results of other simulations that were carried out using standard thermostatting approaches (based on Markovian Langevin and Nosé-Hoover thermostats). We concentrate on the steady-state regime and study the establishment of a local temperature profile within the system. The conditions for obtaining a flat profile or a temperature gradient are examined in detail, in agreement with earlier studies. The results show that the GLE-2B approach is able to treat, within a single scheme, two widely different thermal transport regimes, i.e., ballistic systems, with no temperature gradient, and diffusive systems with a temperature gradient.

Figure 1

Schematic representation of the system. It includes the finite size central system (in blue) where the GLE dynamics is performed, and the two bath $\nu =1,2$ regions at temperature $T_{\nu }$. Because the forces $f_{b_{\nu }}$ and the quantities $g_{i\alpha ,b_{\nu }}$ are of finite range (not necessarily short ranged), one can perform the mapping of ${\mathrm{\Pi }}_{b_{\nu }{b}_{\nu }^{\prime }}\left(\omega \right)$ on a finite region of space: the bath reduced region (in pink), one for each bath.

Figure 2

Model system for a short cylindric Al wire (central system containing 59 atoms) connected to two left $\left(L\right)$ and right $\left(R\right)$ baths. For the GLE calculations, the $L$ and $R$ bath reduced region contains 68 atoms each. For the calculations of the dynamical matrix and the mapping of the polarization matrix, larger baths (each containing 203 atoms) were considered [53]. The embedded atom method (EAM) is used for the interatomic effective potential. For convenience, we now call the baths $\nu =L,R$ instead of $\nu =1,2$.

Figure 3

Schematic representation of systems under consideration for the GLE calculations. The central system consists of 7 layers of Al (labeled A to G) and 3 extra layers embedded in the left bath (labeled L1, L2, and L3) and 3 extra layers embedded in the right bath (labeled R1, R2, and R3). The layers L3, L1, B, D, F, R1, and R3 (L2, A, C, E, G, R2) each contain 5 (4) Al atoms. The bath reduced regions are represented by the brown and blue regions for the left and right baths, respectively. Each of these regions contains 68 atoms. We recall that, during the GLE simulations, the positions of the bath atoms are fixed. The left and right baths are at their own temperature $T_L$ and $T_R$, respectively. For the other thermostatting approaches, the bath reduced regions are described by their own stochastic dynamics [Langevin-Gauss (LG) or Nosé-Hoover (NH) thermostatting] at their own temperatures $T_L$ and $T_R$. The central system (L1–L3, A–G, R1–R3) evolves according to a Newtonian (NVE) dynamics. The remaining outer atoms are kept at fixed positions to ensure the structural stability of the system in the LG and NH thermostatting calculations.

Figure 4

Histograms of the velocity distribution for different layers of atoms for the system shown in Fig. 3. Top left panel: Layers L3, L2, L1 of the central system. Top right panel: Layers R1, R2, R3 of the central system. Bottom panel: Layers A to G. The GLE calculations are performed for $T_L=T_R=200$ K. All distributions (histograms) fit perfectly the equilibrium Maxwell distribution (solid red line) as expected.

Figure 5

Histograms of the velocity distribution for the different layers of atoms in the system shown in Fig. 3. The GLE calculations are performed for two different bath temperatures $T_L=200$ K and $T_R=125$ K [76]. The central layer D is frozen, i.e., atomic positions fixed, and acts as a perfect thermal barrier. The distributions (histograms) for the atomic layers on both sides of the frozen layer D, namely layers (L3, L2, L1, A, B, C) and layers (E, F, G, R1, R2, R3), fit perfectly the equilibrium Maxwell distributions obtained for the two temperatures $T_L$ (lines in the top left and bottom panels) and $T_R$ (lines in the top right and bottom panels).

Figure 6

Histograms of the velocity distribution for different layers of atoms for the system shown in Fig. 3. The GLE calculations are performed for proper nonequilibrium conditions with the bath temperatures being $T_L=150$ K and $T_R=50$ K. The line shape of the different distributions (histograms) is similar to that of the equilibrium Maxwell-Boltzmann distribution. All distributions fit well onto the Maxwell distribution obtained from the average temperature $T_{\mathrm{av}}=(T_L+T_R)/2=100$ K.

Figure 7

Temperature profiles across the central system for different bath temperatures. The labels Lbath (Rbath) represent the temperatures at the left (right) bath. The other labels indicate different layers in the central system. For different sets of $T_{L,R}$, there is no temperature gradient in the system; i.e., it behaves like a perfectly integrable system. There is an inherent uncertainty in the evaluation of the local temperature from a fit to the Maxwell distribution. All temperatures have an error bar of $\pm 5$ K. As a visual guide, this corresponds to the size of the symbols. Note that we also have performed calculations swapping $T_L\leftrightarrow T_R$ and obtained the same flat temperature profiles.

Figure 8

Temperature profiles across the system for different bath temperatures $T_L=150,T_R=50$. Comparison between the GLE-2B approach and the other LG and NH thermostatting approaches. While the GLE-2B calculations provide a uniform temperature profile inside the system, the LG and NH thermostats show either the building up of a temperature gradient in the central system or a flatter temperature profile depending on the value used for the damping parameter. The LG calculations performed with the damping parameter ${\tau }_{\mathrm{damp}}=1/\gamma =7$ ps show a temperature gradient, while almost no gradient is obtained for ${\tau }_{\mathrm{damp}}=100$ ps. The NH calculations show a strong temperature gradient for ${\tau }_{\mathrm{damp}}^{\mathrm{NH}}=0.7$ ps, and a flatter temperature gradient for ${\tau }_{\mathrm{damp}}^{\mathrm{NH}}=20$ ps. There is still a gradient for NH thermostats (with ${\tau }_{\mathrm{damp}}^{\mathrm{NH}}=0.7$ ps) even for smaller $\mathrm{\Delta }T$ with $T_L=125$ and $T_R=75$ (see orange curve).

Figure 9

Temperature profiles across the system from the GLE-2B calculations under different $T_{L,R}$ and different atomic masses. It is possible to break the integrability of the system by considering larger temperatures (red curves) and/or by introducing configurational disorder (blue curve). The disorder is introduced by changing the mass $\left(\pm 20\%\right)$ of some atoms picked up randomly in the central system. For higher temperatures, parts of the system are “driven” beyond the harmonic limit. In such cases, the system is not entirely ballistic and presents a finite thermal conductivity, leading to the buildup of a temperature gradient.

Figure 10

Temperature profiles across the system from the GLE-2B calculations. A simplified form of anharmonic effects in the baths is introduced by modifying the values of the fitted parameters ${\tau }_{k_{\nu }}$. Such anharmonic effects lead to a small temperature gradient on the side of the hotter bath $(L$ side). The building up of the temperature gradient is however not as strong as for the LG/NH thermostatting calculations.

Figure 11

Schematic description of the GLE-2B and the LG/NH calculations. Upper panel: The GLE dissipative dynamics Eq. (7) is performed only on the atoms of the central region. The $L$ and $R$ baths (displaced farther to the left and right for clarity) enter into the GLE calculations via the sets of fitting parameters $\{c_{b_{\nu }}^{\left(k_{\nu }\right)},{\tau }_{k_{\nu }},{\omega }_{k_{\nu }}\}$ and the $g_{i\alpha ,b_{\nu }}\left(\left\{r_{i\alpha }\right\}\right)$ quantities. In such calculations, the positions of the bath atoms are fixed at their equilibrium positions shown in Fig. 3. Lower panel: The atoms in the central region follow Newton's EOM and the atoms in the $L$ and $R$ baths follow a dissipative LG or NH dynamics. The atoms with a black cross are kept fixed to ensure the overall stability of the system. The two approaches, GLE-2B and LG/NH, represent two different kinds of stochastic processes.

Figure 12

(a) Evolution of the total kinetic energy $E_{\mathrm{kin}}^{\mathrm{TOT}}$ and of the total potential energy $E_{\mathrm{pot}}^{\mathrm{TOT}}$ of the system shown in Fig. 2 versus time $(\mathrm{\Delta }t=2$ fs) [115]. The calculations are performed with the bath temperatures $T_L=150,T_R=50$. The horizontal dashed blue lines correspond to the averages of the energy over the time range 100 to 300 ps. The convergence is fairly well obtained after $\sim 80$ ps (40 000 time steps). (b) Evolution of the total kinetic energy $E_{\mathrm{kin}}^{\mathrm{TOT}}$ of the central system only (59 atoms from layers L3 to R3). For the LG calculations, the convergence towards the stationary value occurs after roughly 80 ps for the damping time ${\tau }_{\mathrm{damp}}=7$ ps. The stationary state is reached after more time steps, as expected, for the larger damping time ${\tau }_{\mathrm{damp}}=100$ ps (see inset where it is apparent that the stationary value is reached after 500 ps or 250 000 time steps).

Figure 13

Histograms of the total kinetic energy of the central system containing 59 atoms for the bath temperatures $T_L=150,T_R=50$. The averaged kinetic energy $\langle E_{\mathrm{kin}}^{\mathrm{TOT}}\rangle \sim 0.8$ eV is the energy reference.

Figure 14

Temperature profiles across the system consisting only of the 7 layers labeled A to G in Figs. 3 and 11. The different bath temperatures are $T_L=200,T_R=100$. While the GLE-2B calculations provide a uniform temperature profile inside the system, the LG thermostats show either the building up of a temperature gradient in the central system $({\tau }_{\mathrm{damp}}=7$ ps) or a flat temperature profile $({\tau }_{\mathrm{damp}}=100$ ps).

Figure 15

Model for a one-dimensional Al chain (central system containing 11 atoms labeled A to K) connected to the left $\left(L\right)$ and right $\left(R\right)$ baths. For the GLE calculations, the $L$ and $R$ bath reduced region contains 30 atoms each. For the calculations of the dynamical matrix and the mapping of the polarization matrix, larger baths (containing each 95 atoms) were considered. The mapping of the polarization matrix elements is done using 78 aDOFs with a set of 78 parameters $\{c_{b_{\nu }}^{\left(k_{\nu }\right)},{\tau }_{k_{\nu }},{\omega }_{k_{\nu }}\}$ for each bath $\nu =L,R$.

Figure 16

Temperature profiles across the central 1D chain for different bath temperatures $T_{L,R}$. The labels Lbath (Rbath) represent the temperatures at the left (right) bath. The other labels A to K indicate different atoms in the central chain. For “low” temperatures $T_{L,R}$, the system is harmonic and one gets a flat temperature profile as expected for integrable systems. For “higher” temperatures, the motion of some of the atoms samples the anharmonic part of the interatomic potential, and therefore the system is nonintegrable and a temperature gradient is built up.

Figure 17

Temperature profiles across the central 1D chain for different defects present in the chain. The labels Lbath (Rbath) represent the temperatures at the left (right) bath. The other labels A to K indicate different atoms in the central chain. The mass of atom H in Fig. 15 is increased by 10%, 20%, 30%, or 40%. The case of two defects (atoms C and H with their mass increased by 40%) is also shown. For the temperatures considered, the chain is in the harmonic regime. However the presence of defects, and their associated more localized vibration modes, leads to the build up of a temperature gradient across the chain.