A Molecular Dynamics Simulation Study of In- and Cross-Plane Thermal Conductivity of Bilayer Graphene

1. Introduction

Thermal management at the nanoscale devices is a significant challenge due to the rapid increase in power densities [1,2,3]. For state-of-the-art electronics, having the appropriate thermal conductivity of the material, whether it’s high or low, is crucial for thermal management. Development of the next generations of high-power light-emitting diodes (LEDs) and integrated circuits (ICs) demands efficient heat dissipation. This can be achieved using materials with high thermal conductivity to maintain their operational performance and long-term reliability. On the other hand, materials with reduced thermal conductivity are essential for thermoelectric devices, thermal insulators, and even phononic computing devices to ensure a high figure of merit in such applications [4,5].

In recent years, novel materials, especially low-dimensional materials like quasi-one-dimensional nanowires [6,7], nanotubes [8], and two-dimensional (2D) materials [9,10], have attracted significant attention all over the world. Many materials exhibit unique properties in their low-dimensional state that are not observed in their three-dimensional structures. Graphene, the first 2D material discovered via exfoliation in 2004 by Novoselov [11], has exceptional properties in stability, conductivity, and flexibility [12]. Graphene is a single layer of carbon atoms (one atom thick), densely packed in a 2D honeycomb lattice structure which characterized by sp² hybridization. The strong sp² bonding, combined with the lightweight nature of the carbon atoms in the microscopic structure, imparts to graphene exceptional physiochemical characteristics. These include a high Young’s modulus (∼1.0 TPa), a large specific surface area (2630 m²/g), and an optical transmittance of 97.7% [13,14,15]. Compared to other mentioned unique features of graphene, its novel thermal properties have gained the most interest. The pioneering experimental measurements, utilizing optothermal Raman method, revealed that single-layer graphene (SLG) could potentially exhibit a thermal conductivity ranging from 2000 to 5000 W/m.K at room temperature [16]. These values surpass the bulk graphite limit of K=2000 W/m.K for basal planes at room temperature, exceeding even that of diamond. The remarkable thermal conductivity of graphene arises from the strongest covalent sp² bonds in nature. Experimental research has shown that the thermal transport properties of SLG depend on various factors, such as the geometry size [16,17,18,19], axial strain [20], and functionalization [21].

In addition to SLG, bilayer graphene (BLG) and multilayer graphene sheets also exhibit unique behaviors and possess properties that differ from those of SLG. For example, the room temperature thermal conductivity of few-layer graphene varies from 1300 W/m.K to 2800 W/m.K [22]. Multilayer graphene often exhibits notable thermal conductivity anisotropy, wherein the conductivity perpendicular to the layers (cross-plane) is typically lower than the conductivity along the interfaces (in-plane). Given this anisotropic behavior of graphene, it is necessary to conduct separate investigations into both the in-plane thermal conductivity (k_ip) and cross-plane thermal conductivity (k_cp) of BLG and multilayer graphene. While there have been numerous experimental studies on in-plane thermal transport, research on the cross-plane heat transfer of graphene is relatively limited. In a recent study, Sood et al. [23] investigated cross-plane heat transfer in superlattices assembled from single layers of graphene and Mo${\mathrm{S}}_2$ employing a combinatorial experimental approach.

The novel thermal properties of SLG, BLG, and multilayer graphene have been extensively studied by both experimental methods and numerical simulations. While the mentioned experimental research has yielded direct measurements of these properties of graphene, the substantial costs associated with experimentation and the limited ability of nanoscale temperature probing have presented considerable challenges for conducting more extensive investigations. Numerical methods are promising tools to complement experiments and enhance our comprehension of the distinctive thermal properties exhibited by graphene. Molecular dynamics (MD) simulation is one of the main numerical methods for studying the properties of 2D materials, such as graphene. In recent years, there has been a dramatic increase in MD simulations to explore the unique thermal properties of SLG [24,25,26,27,28,29]. References [30] and [31] provide an extensive review of the latest developments in MD simulation concerning the remarkable thermal properties of graphene.

While there have been numerous studies dedicated to MD simulations of SLG, the attention given to MD simulations of BLG or multilayer graphene is relatively limited. Wei et al. [32] employed the non-equilibrium molecular dynamics method (NEMD) to calculate the k_ip of single- and few-layered graphene nanoribbons. Their findings indicated a reduction in the k_ip of graphene as the number of layers increased. Cao et al. [33] investigated the k_ip of multilayer graphene using NEMD simulations. Their results showed that as the size of multilayer graphene decreased, its k_ip also decreased. Guo et al. [34] determined the k_ip of interlayer-bonded BLG, tri-layer graphene, and pyrolytic graphite using NEMD. They discovered that interlayer sp³ bonding reduces thermal conductivity, with reductions of up to 80%. Additionally, they observed that tensile strain further decreases the thermal conductivity of interlayer-bonded graphene by 50 %. Zhan et al. [35] calculated the k_ip of a BLG sheet approximately 64.8 W/m.K employing the reverse non-equilibrium molecular dynamics (RNEMD) method. They also observed an 80% reduction in thermal conductivity with a 1.39% increase in bridge defect. Si et al. [36] used NEMD method to study the k_ip of both SLG and multilayer graphene sheets, utilizing Tersoff, REBO, and AIREBO potentials. They concluded that for SLG, the optimized Tersoff potential produced the highest thermal conductivity, closely aligning with experimental results. Regarding multilayer graphene, they concluded that the Tersoff potential provides the best simulation of its thermal conductivity. By using the MD method with REBO and Leonard Jones potentials, Li et al. [37] calculated the k_ip of misoriented BLG. They explored the effect of interlayer misorientation angle and temperature using NEMD. The study’s conclusion highlighted that thermal conductivities exhibit an approximately linear decrease with the increasing lattice constant of the commensurate misoriented BLG unit cell. Nie et al. [38] examined the k_ip and k_cp of twisted multilayer graphene sheets utilizing NEMD. The results revealed that the highest thermal conductivity could be achieved when no twist angles were present, for both in-plane and cross-plane conduction. Zhang et al. [39] studied the k_ip of BLG by the NEMD, employing different in-layer and interlayer potentials. The Leonard Jones and Kolmogorov-Crespi potentials described the interlayer interactions, while the in-layer interaction potentials, namely Tersoff, REBO, and AIREBO potentials, were compared. Their results showed significant differences in both the magnitude and temperature-related trends of thermal conductivity when different potentials were applied in the simulations.

Reviewing the aforementioned research reveals that, in comparison to the research on in-plane thermal transport, there has been limited focus on investigating the cross-plane conduction at the interface of graphene-graphene. This signifies an existing gap in our understanding that requires resolution. To the best of the authors’ knowledge, a comprehensive investigation into the in-plane and cross-plane conductivity of BLG has been rarely documented in published studies. Therefore, this work aims to explore the thermal conductivity of BLG through NEMD, considering both k_ip and k_cp. Moreover, the effects of several parameters, such as sample dimensions, edge shapes, temperature, and strain, on thermal conductivity are investigated. This study is expected to yield valuable insights for the design and application of BLG in nanoscale devices, particularly where the manipulation of thermal conductivity, whether to suppress it or enhance it, in different directions is essential.

2. Simulation procedures

In this research, we utilized NEMD simulations to calculate the thermal conductivity of BLG. These simulations were executed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [40]. Atom visualization was accomplished through the use of the Open Visualization Tool (OVITO) [41]. Moreover, the lattice structure for the input script was provided using Visual Molecular Dynamics (VMD) [42].

The schematic diagram of BLG is shown in Figure 1a. This figure illustrates the hexagonal structure of zigzag BLG, depicting the relative position of the two layers. Figure 1b illustrates the simulation model for the NEMD simulations of the k_ip. The simulation domain comprises three regions: the main simulation region, reservoir regions, and fixed walls. Two fixed walls made up of fixed atoms are positioned at the model’s ends, in order to prevent atom escape from the system and limit relative sliding motion between graphene layers. Next to these fixed walls, cold and hot reservoirs (heat sink and heat source) are established to create a temperature gradient. Once the simulation structure reaches an equilibrium state, the atom velocities in the reservoir regions are artificially modified to provide a constant heat flux into the main simulation region, situated between the heat source and heat sink. During each time step, the particle velocities within the heat source and heat sink are uniformly scaled, ensuring that a uniform quantity of kinetic energy is consistently added to or deducted from both reservoirs. This process provides a constant heat flux within the simulation domain, establishing a temperature gradient aligned with the direction of heat flow. From the calculated temperature gradient and heat flow, the thermal conductivity of the graphene can be then determined from the Fourier law. To assess the localized temperatures, the model is segmented into 25 equal divisions along the X-axis direction.

To simulate the k_cp of the BLG with the NEMD method, the initial configuration of a 6-layer structure needs to be established. The top and bottom layers are set as the fixed layers, and then the heat sink and heat source are positioned as shown in Figure 1c.

The optimized Tersoff potential was employed to describe in-layer interactions [43]. In addition, the van der Waals interaction between distinct layers was modeled using the Lennard-Jones potential, with parameters obtained from Girifalco et al. [44]. The details of the interatomic potentials are provided in Appendix A. These potentials have been proved the most appropriate potentials for modeling the thermal conductivity of multi-layer graphene [38]. Incorporating the Lennard-Jones term into the optimized Tersoff potential to account for interlayer atomic interactions, the optimized Tersoff potential predicts a broader frequency range and a higher value for the high-frequency peak. Furthermore, the graphene’s temperature was determined by calculating the average kinetic energy of every atom within the layer, as indicated in Equation (1):

$$\mathrm{T}=\frac{2}{3\mathrm{N}{\mathrm{k}}_{\mathrm{b}}}\sum _{\mathrm{i}}\frac{{\mathrm{p}}_{\mathrm{i}}^2}{2{\mathrm{m}}_{\mathrm{i}}}$$

(1)

where ${\mathrm{k}}_{\mathrm{b}}$ represents the Boltzmann constant, p is the momentum of atoms, m is the mass of atoms, and N represents the number of atoms.

While simulating k_ip and k_cp, the simulation box was first relaxed in an isothermal-isobaric (NPT) ensemble for 100 ps, which proved sufficient to achieve equilibrium. The heat source and heat sink were controlled by Langevin thermostats to maintain constant temperatures. Then, an additional 1 ns was conducted within a microcanonical (NVE) ensemble to ensure the attainment of a desirable linear temperature profile for k_ip calculations. For simulating k_cp, a longer simulation time of 2 ns in the NVE ensemble was necessary to reach a good equilibrium of the temperature profile. This extended duration was required due to the lower thermal conductivity in the cross-plane direction compared to the in-plane direction. Periodic boundary conditions were applied along the X, Y, and Z axis directions.

A typical temperature profile of the BLG for the in-plane conduction calculations is depicted in Figure 2a. The temperature gradient was determined using the linear region of this profile. Notably, a temperature jump is observed at the two boundaries. This temperature jump near the heat baths can be attributed to the presence of localized edge modes of phonons [45,46].

The k_ip was computed using Fourier law as specified in Equation (2):

$$\mathrm{q}=-{\mathrm{k}}_{ip}\mathrm{A}\frac{\partial \mathrm{T}}{\partial \mathrm{x}}$$

(2)

where q is the heat flux obtained from the NEMD simulation, A is the cross-section area, ∂T/∂x is the temperature gradient, and k_ip is the in-plane thermal conductivity coefficient.

Figure 2b displays a representative temperature profile for cross-plane calculations in BLG. To determine the cross-plane thermal conductance, Equation (3) was used:

$$\mathrm{G}=\frac{\mathrm{q}}{\mathrm{A}∆\mathrm{T}}$$

(3)

where G represents the cross-plane thermal conductance, which is the inverse of the cross-plane thermal resistance (${\mathrm{R}}_{\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}})$:

$${\mathrm{R}}_{\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}}=\frac{1}{\mathrm{G}}$$

(4)

3. Results and discussion

3.1. In-plane thermal conductivity

In order to check the accuracy of the NEMD simulation results for the k_ip of BLG, we conducted simulations following the dimensions and operating conditions outlined in the study by Nie et al. [38]. The investigated bilayer zigzag graphene had a size of 22×10 nm² with a 0.35 nm interlayer distance. The comparison of results is presented in Table 1 for different temperatures. We performed error analysis using four different heat fluxes across the simulation region. As indicated in Table 1, the k_ip calculated in this study exhibits only slight deviations from the findings of Nie et al. [38], with a maximum error of 3.46%.

In addition, comparing the range of k_ip obtained for other sizes of BLG, as reported in the upcoming sections, with the values from prior studies, indicates that our NEMD simulation results are consistent with the findings of previous research [36,39].

In our study, we investigated BLG sheets with dimensions of 10×10 nm² and an interlayer distance of 0.335 nm for the calculation of k_ip. To quantify the impacts of various factors on the k_ip of BLG, we varied several key parameters such as length, width, temperature, strain, and edge shapes.

Figure 3 illustrates the influence of sample length on the k_ip of BLG with armchair and zigzag edges at 300 K. The thermal conductivities of 10 nm × 10 nm zigzag and armchair BLG are estimated to be 203.6 W/m.K and 124.3 W/m.K, respectively. Contrary to bulk materials, where thermal conductivity is typically size-independent, the results in Figure 3 show that the thermal conductivity of BLG increases with the length of the BLG sheets. The average mean free path of the primary heat carriers, namely acoustic phonons, is in the range of 100-600 nm in graphene [47]. In the present study, the simulation domain is not sufficiently extensive to enable dominant phonons to undergo umklapp scattering. Previous reports indicate that phonons in graphene can travel without scattering, exhibiting ballistic transport [48]. As a result, it is the heat source and the heat sink that ultimately limit the phonon mean free path, leading to the observed increase in thermal conductivity with the sample length. This behavior is in agreement with previously reported experimental findings, where the k_ip of SLG was observed to increase with the length of the graphene sheet [19]. The length-dependent thermal conductivity of 2D materials is not limited to cases where the sample size is comparable to the average phonon mean free path; it’s observed even when sample lengths significantly exceed the average phonon mean free path [19]. However, for sample lengths exceeding the average phonon mean free path, the phenomenon can be attributed to the 2D nature of phonons in graphene.

Furthermore, Figure 3 shows that the thermal conductivity of the BLG with a zigzag edge is 30-40% higher than that of the armchair edge. Previous studies using various methods, such as MD [24,49,50], Green’s function method [51], or solving the Boltzmann transport equation [52], have shown that the thermal conductivity of graphene with zigzag edges is greater than that of armchair graphene, with reported differences ranging from 15% to 50%. Similar trends have been observed in other 2D nanostructures as well [53,54]. The dependence of thermal conductivity on edge shape in finite-sized graphene sheets can be attributed to several factors, including different phonon scattering rates at the armchair and zigzag edges [49,50,52], edge roughness scattering [52], and phonon localization at edges [24,50]. However, Wang et al. [24] argued that edge roughness scattering should be excluded as the cause for the large thermal conductivity discrepancy in graphene with different edge shapes. They supported their argument with a cross-sectional decomposition of the steady-state heat flux in graphene nanoribbons, revealing a significant suppression of thermal transport at the edges, particularly in armchair configurations. Their phonon spectra analyses and observations of nonuniform heat flux distribution suggested that the strong edge localization of phonons in regions near and at the edges, particularly in armchair edges, is the primary factor responsible for the observed thermal conductivity discrepancy between armchair and zigzag edges [24].

Figure 4 illustrates the effect of BLG width on the kip for three different BLG lengths. It is evident from this figure that the thermal conductivity of BLG sheets with zigzag edges is higher at very small widths, particularly for the smallest BLG lengths (10 nm). This phenomenon could be linked to the limited number of phonons present in the system at these small sizes. This scarcity results in fewer available phonon-phonon combinations that satisfy the momentum and energy conservation rules for scattering, leading to reduced phonon-phonon scattering. The higher thermal conductivity at very small widths aligns with findings from other research [37], where the thermal conductivity of BLG becomes width-independent for widths greater than 6 nm at 300 K.

In addition, Figure 4 reveals that the width of armchair BLG sheets does not exhibit any noticeable effect on their thermal conductivity. These findings are ascribed to the influence of periodic boundary conditions in the simulation, consistent with prior research studies [55].

Figure 5 presents the effect of temperature variation on the thermal conductivity of zigzag and armchair BLG, covering the temperature range from 300 to 600 K for three different BLG lengths. As it has been proved that MD simulations are not well-suited for accurately predicting the thermal conductivity of graphene with temperature below 300 K [38,56], hence we focus on temperatures above this threshold. Based on this figure, the thermal conductivity of BLG decreases with increasing temperature, showing a slightly different rate for armchair and zigzag configurations. For a 10 nm× 10 nm sample, the thermal conductivity of the zigzag configuration drops from 203.6 W/m.K to 170.8 W/m.K as the temperature rises from 300 K to 600 K, while for the armchair configuration, it decreases from 124.3 W/m.K to 111.1 W/m.K. The behavior of thermal conductivity for BLG follows a pattern similar to that of bulk dielectric material, decreasing as the temperature rises beyond room temperature. This temperature dependence at high temperatures is attributed to umklapp phonon-phonon scattering. Experimental studies have also addressed the effect of temperature on the thermal conductivity of BLG, with results confirming that BLG’s thermal conductivity decreases as the temperature rises from 300 to 600 K [57,58].

Moreover, as observed in Figure 5, zigzag configurations consistently exhibit higher thermal conductivity than the armchair configurations at all temperatures, as discussed previously.

Prior research has demonstrated the substantial influence of strain on modulating the thermal conductivity of various 2D materials. The application of strain presents a method to precisely tune the thermal conductivity of materials, as evidenced by these investigations [20,34,59]. Figure 6 illustrates the impact of uniaxial tensile strain on the thermal conductivity of armchair and zigzag BLG sheets. The range of strain considered in this study spans from 0 to 12.5%. This figure demonstrates that increasing strain can result in a decrease in thermal conductivity for both armchair and zigzag configurations. This trend is in good agreement with previous NEMD simulations [59] and experimental studies [20] conducted on SLG confirming the sensitivity of thermal conductivity to tensile strain. The decrease in thermal conductivity with tensile strains can be ascribed to the increase in lattice anharmonicity and the reduction of the stiffness tensor.

Furthermore, it’s noteworthy that the thermal conductivity of zigzag BLG demonstrates a higher sensitivity to tensile strain compared to armchair BLG. For a tensile strain of 12.5%, the thermal conductivity reduction is 51% in 10 nm, 52% in 14 nm, and 45% in 18 nm compared to the strain-free values for the zigzag BLG. On the other hand, the thermal conductivity reduction is 29% in 10 nm, 32% in 14 nm, and 30% in 18 nm armchair BLG sheets. This behavior is in agreement with the findings in reference [59], where the thermal conductivity of zigzag graphene nanoribbons exhibited greater sensitivity to the tensile strain compared to armchair graphene nanoribbons. The different sensitivity observed can be attributed to the difference in bond orientation along the direction of loading. Indeed, the effect of tensile strain on the thermal conductivity of graphene with zigzag and armchair edges has been the subject of debate due to the presence of contradictory findings. In another MD simulation, it was observed that the thermal conductivity of armchair graphene nanoribbons exhibited greater sensitivity to tensile strain compared to that of zigzag graphene nanoribbons [60]. This heightened sensitivity was attributed to the greater stress and larger deformation of carbon-carbon bonds along the tensile direction in armchair graphene nanoribbons than zigzag graphene nanoribbons under the same level of strain [60]. Many of these discrepancies can be attributed to computational details, such as the specific interatomic potentials used in simulations and sample size. The true behavior will remain uncertain until relevant experimental results become available.

3.2. Cross-plane thermal conductivity

The characteristics of cross-plane thermal conduction differ significantly from those of in-plane conduction in certain materials. For example, in the case of graphite, the k_cp has been shown to be two to three orders of magnitude lower than its k_ip [61]. The thermal properties of graphene are significantly influenced by those of graphite, given that graphene can retain the anisotropic nature of the graphite crystal. Previous measurements indicate that the cross-plane thermal properties of few-layered graphene are also considerably smaller than its in-plane thermal properties [23,62].

The effects of temperature on the cross-plane thermal conductance of BLG are illustrated in Figure 7. The graph shows that the cross-plane thermal conductance of BLG at 300 K is 162.3 MW/m².K. To the best of our knowledge, there have been no studies conducted on the thermal conductance of the graphene/graphene interface within BLG. Wei et al. [63] studied the cross-plane thermal resistance in multilayer graphene structures using NEMD simulations, reporting a value of 4.1×10^-9 m².K/W for the cross-plane thermal resistance of 6-layer graphene at 300 K, as shown in Figure 4 of Ref.[63], which corresponds to a cross-plane thermal conductance of 244 MW/m².K. In another research, Ni et al. [64] calculated the cross-plane thermal resistance of 5-layer graphene as 4.83×10^-9 m².K/W using equilibrium molecular dynamics simulations, corresponding to a cross-plane thermal conductance of 207 MW/m².K. In addition, Ding et al. [65] explored interfacial thermal conductance in graphene/MoS₂ heterostructures using NEMD simulations. For 7-layer graphene, the interfacial thermal conductance at the graphene/graphene interface was calculated as 212.6 MW/m².K, as shown in Table 2 of Ref. [65]. Considering the strong layer number dependence of cross-plane thermal resistance/conductance in multilayer graphene structures and the documented increase in cross-plane thermal conductance with layer number [63], it is reasonable that our obtained cross-plane thermal conductance for BLG is lower than that of 5-, 6-, and 7-layer graphene. It’s worth highlighting that when few-layer graphene samples on a substrate are studied, the cross-plane thermal resistance does not significantly change with the number of graphene layers [23,62], indicating that the thermal resistance between graphene and its substrate has a more notable effect than the resistance between individual graphene sheets.

Our simulation results reveal that the corresponding k_cp of BLG at 300 K is 0.054 W/m·K, which closely aligns with both the reported values of 0.05 W/m·K for BLG (Figure 5 of Ref. [66]) and 0.07 W/m·K for 6-layer graphene [63]. These findings highlight that the k_cp is significantly lower than the in-plane results by approximately four orders of magnitude. This anisotropic behavior arises from the high k_ip of graphene, which is ascribed to the strong covalent sp² bonding between adjacent carbon atoms. This type of bonding is known to be one of the strongest in nature, marginally stronger than the sp³ bonds in diamond. In contrast, the adjacent graphene layers in multilayer graphene are linked by weak van der Waals interactions.

Furthermore, Figure 7 demonstrates that the cross-plane thermal conductance of BLG decreases with rising temperatures in the range of 300 to 600 K. This trend suggests that temperature impacts the cross-plane thermal conductance similarly to the k_ip, which monotonically decreases with increasing temperature. This observation aligns qualitatively with prior research on the temperature effect on cross-plane heat transfer in few-layered graphene [63,66,67]. The underlying physics behind this phenomenon lies in the dominance of phonon scattering in graphene heat transfer. Consequently, heat transfer should decrease as the temperature rises, independent of the direction of heat transport. In addition, it’s evident from Figure 7 that the temperature-induced reduction is not particularly pronounced: the thermal conductance at 600 K is approximately 16% of that at 300 K, consistent with the behavior reported for 5-layer graphene, where cross-plane resistance was found to be almost unaffected by temperature [64].

4. Conclusion

This study investigated the in-plane and cross-plane thermal conductivity of BLG sheets using NEMD. We employed the optimized Tersoff potential to describe in-layer interactions, while the van der Waals interaction between different layers was modeled using the Lennard-Jones potential. The obtained results for both k_ip and k_cp of BLG sheets exhibited strong agreement with published data. Furthermore, we conducted a parametric study and thoroughly analyzed the resulting outcomes.

Our results showed that the k_ip of 10 nm × 10 nm zigzag and armchair BLG sheets at room temperature are around 203.6 W/m.K and 124.3 W/m.K, respectively. The k_ip of BLG increased with the length of the sheets: For 10 nm × 30 nm sheets, the values for the zigzag and armchair configurations increased to 1214.15 W/m.K and 773.15 W/m.K, respectively. The influence of sheet width, due to the applied periodic boundary conditions, was less evident on BLG’s thermal conductivity. In addition, the k_ip of BLG decreased with increasing temperature, with a slightly divergent rate observed between armchair and zigzag configurations. A 10 nm × 10 nm zigzag BLG sample, for instance, saw a reduction from 203.6 W/m.K to 170.8 W/m.K as the temperature increased from 300 K to 600 K. Conversely, the armchair configuration’s thermal conductivity decreased from 124.3 W/m.K to 111.1 W/m.K over the same temperature range. Furthermore, applying uniaxial tensile strain on the BLG sheet led to a reduction in thermal conductivity for both armchair and zigzag configurations, with zigzag BLG displaying a higher sensitivity compared to armchair BLG. A tensile strain of 12.5% on a 10 nm× 10 nm sample resulted in a thermal conductivity reduction of 51% for zigzag BLG and 29% for armchair BLG when compared to strain-free values.

Our findings also demonstrated that the cross-plane thermal conductance of BLG at 300 K is 162.3 MW/m².K with the corresponding k_cp of 0.054 W/m·K. These findings highlighted that the k_cp is significantly lower than the in-plane results by about four orders of magnitude. As temperatures increased in the range of 300 to 600 K, the cross-plane thermal conductance exhibited a reduction. However, the temperature-induced reduction was not highly pronounced and the thermal conductance at 600 K was approximately 16% of that at 300 K.

The high anisotropic heat conduction properties observed in BLG sheets hold promise for applications in the thermal management of nanoscale devices that require controlled regulation of thermal conductivity in different directions.