Thermal Conductance along Hexagonal Boron Nitride and Graphene Grain Boundaries

: We carried out molecular dynamics simulations at various temperatures to predict the thermal conductivity and the thermal conductance of graphene and hexagonal boron-nitride (h-BN) thin ﬁlms. Therefore, several models with six different grain boundary conﬁgurations ranging from 33–140 nm in length were generated. We compared our predicted thermal conductivity of pristine graphene and h-BN with previously conducted experimental data and obtained good agreement. Finally, we computed the thermal conductance of graphene and h-BN sheets for six different grain boundary conﬁgurations, ﬁve sheet lengths ranging from 33 to 140 nm and three temperatures (i.e., 300 K, 500 K and 700 K). The results show that the thermal conductance remains nearly constant with varying length and temperature for each grain boundary.


Introduction
Due to its exceptional material properties graphene is a promising material for numerous industrial applications ranging from nanoelectronics to the aerospace industry. As the strongest man-made material, graphene possesses a tensile strength of 130 ± 10 GPa, an elastic modulus of 1 ± 0.1 TPa, and also a high thermal conductivity of 4100 ± 500 W/m·K [1][2][3][4][5][6].
The discovery of graphene has motivated scientists to study other 2D materials such as hexagonal boron-nitride (h-BN), which also exhibits noteworthy mechanical, thermal, and electrical properties but in contrast to graphene, it possesses a natural band gap. The physical structure of graphene; i.e., a thin-layer of carbon atoms forming a honey-comb lattice structure-enables it to be used in applications such as hetero-structured photo detectors, strain sensors, and even composites combined with aluminum to improve heat-dissipation properties [7][8][9][10]. Possessing a similar physical structure as that of Graphene, h-BN shows very good mechanical, thermal insulating, and dielectric properties, in addition to chemical stability [11]. By using different combinations of stiffness, ultimate tensile strength and strain, the mechanical properties of h-BN can be tuned and used in advanced applications such as strain-engineered nano-ribbons, bio-compatible nano-devices and nano-electronics [12]. The heterogeneous integration of both materials, due to their similar physical structures, results in the increase of the overall tensile strength, thermal, and electrical conductivities, which is useful for developing semiconductor-insulators [13][14][15][16][17][18][19].
Thermal chemical vapor deposition (CVD) is the most common technique used to produce thin-film sheets of graphene and hexagonal boron-nitride. During the fabrication process, polycrystalline growth leads to grain boundaries and defects along the boundaries of the lattice surface. These defects are of different shapes; the most common type being pentagon-heptagon shaped. Less common shapes include square-octagon pairs. The pentagon-heptagon shaped pairs are comprised of homonuclear boron-boron or nitrogen-nitrogen pairs rather than heteronuclear boron-nitrogen pairs, which constitute the rest of the hexagonal lattice structure in the entire thin-film sheet. These defects act as dislocations and cause stress concentrations and phonon scattering of the grains reducing the tensile strength as well as the thermal and electrical conductivity of the materials, thereby affecting its structural stability during loading conditions [13][14][15][16][17][18][19][20].
Computational methods such as molecular dynamics (MD) have been employed to study the effect of grain boundaries. Defect patterns are modelled for symmetric and non-symmetric grain boundaries that occur at different angles depending on the armchair and zigzag orientation of the defect-patterns along the lattice structure. The number of grain boundaries increases with increasing misorientation angle along the boundaries of the polycrystalline lattice structure. This rise in the number of grain boundaries leads to an enhanced thermal and electrical resistance, which might affect their efficiency. Hence, a basic understanding of the grain boundaries and their impact on thermal, mechanical, and electrical properties is essential for potential industrial applications [21,22].
In this study, we carried out MD simulations in order to investigate the thermal properties of h-BN and graphene sheets containing ultra-fine grains up to 1 nm. The thermal properties were computed for different sheet dimensions in order to fully understand the effect of the grain boundaries on the surface of the lattice sheets and their influence on the heat flow across the surface of the materials.

Modeling Using Molecular Dynamics
In this study, we used our own C++ script to construct the required grain boundaries while the MD simulations were done with the LAMMPS package (Large-scale Atomic/Molecular Massively Parallel Simulator) [23]. Due to the combination of the planar geometry and its high thermo-mechanical and oxidation-resistant properties, hexagonal boron-nitride can be employed for several other industrial applications including high temperature membrane filters, as fillers in reinforced polymeric materials, and high oxidation-resistant coatings for metals such as nickel. Moreover, graphene may be employed as transparent conducting layer for high temperature thin film device applications where temperatures can reach up to 1100 • C [24][25][26]. Hence, the influence of various sheet lengths ranging from 33 to 140 nm and temperatures of 300, 500 and 700 Kelvin on the thermal conductivity of pristine models were investigated. Studying the thermal behavior of the materials at higher temperatures helped us to better understand the effect of phonon scattering. In the following sections, the effect of the aforementioned parameters on thermal conductance for six different grain boundaries will be assessed.
In past numerical studies, it was observed that the optimized Tersoff potential, proposed by Lindsay and Broido [27], accurately predicts the thermal conductivity of pristine h-BN and Graphene [13][14][15][16][17][18][19][20]. Hence, this study employs the Tersoff potential for modeling the atomic interaction, which takes the form [28][29][30][31]: f c being a cutoff function, r ij is the distance between atom i and j, b ij accounts for the bond order and the local environment, f R and f A are the repulsive and attractive pair potentials, respectively. Figure 1 shows the pristine molecular dynamic (MD) model for extracting the thermal conductivity of single-layer h-BN and Graphene films with the length (denoted as L) of 33 nm. Similar to the work in [20], periodic boundary conditions were employed along the width of the h-BN and graphene sheets in their planar directions for all the developed MD models. Hence, the results obtained implicate an infinite system of grain boundaries.  Figure 2 shows the MD model developed for computing the thermal conductance for both materials, which is composed of two pristine sheets (red and blue sheets) that are connected by six types of grain boundary (as shown in Figure 3). We then divided the sample into 20 slabs excluding the fixed zones at both ends. In Figure 2, L L and L R are the length of the left and right sheets, respectively.  To study the thermal conductance between two separate h-BN and Graphene sheets in Figure 2, we constructed six grain boundaries shown in Figure 3. The 6 plots represent the defect patterns with varying gap sizes filled with hexagonal boron-nitrogen or graphene-graphene bonds. The rate of phonon scattering along the grain boundaries has a direct correlation with the increase in defect concentration. Hence, from several possible combinations of heptagon-pentagon shaped defect pairs, we have chosen the grain boundaries in which the defect concentration ranges from nil to a few dislocations located in-between the defect pairs. Furthermore, our choice of grain boundary selection allows us to study the change in thermal resistivity as the distance between the defect pairs gradually changes. In the following section, the procedure to construct the grain boundaries is described: The grain boundaries as depicted in Figure 3 are constructed based on the number of hexagonal lattices presented between each heptagon-pentagon (5-7 pair) pair along the line of defect.
As explained in [32], the complete polycrystalline h-BN lattice sheets are a combination of two individual pristine h-BN lattice sheets joined together. The common side of the sheet forms the molecular grain boundaries based on the angle at which one sheet or both sheets are tilted.
The misorientation angle θ is the sum of both angles of the left and right domains and can be computed by θ = θ L + θ R , where θ L and θ R are the corresponding angles in the left or right domain of the polycrystalline h-BN sheets in their respective crystallographic orientations. Each domain is characterized by the translation vectors, n L , m L , n R and m R . These vectors represent the left and right domains of the lattice structure and are represented as (n L ,m L ) | (n R ,m R ) [22,32].
According to [33], the translation vectors are responsible for the periodic arrangement of the defects along the grain boundaries. The following expression explains the technique used to calculate the angle of misorientation by using the corresponding translation vectors.
The misorientation angle of the armchair oriented GB's can be determined from the misorientation angle of the zigzag oriented GB's, through the relation: From the values of the translational vectors representing the left and right domains and using them in Equation (2), the misorientation angles were calculated and the symmetric grain boundaries (Figure 3a-d) could be constructed. GB.1-4 represent the symmetric GB's.
Since the angles calculated for the symmetric grain boundaries are oriented in a zigzag manner, they must be converted to their armchair orientation using Equation (3).
The non-symmetric grain boundaries (Figure 3e,f) were also constructed using Equation (2), but were not converted into their armchair orientation due to the irregularity in the positioning of the heptagon-pentagon pairs along the line of defect. From their translation vectors, the misorientation angles for the non-symmetric grain boundaries were calculated. GB.5 and 6 represent the non-symmetric GB's. During the construction of the chosen non-symmetric grain boundaries, the formation of adatoms (additional atoms in the heptagon ring) was unavoidable and hence, a single carbon or boron atom had to be implemented in the molecular structure to overcome this issue. Edge-dislocations were particularly observed during the construction of non-symmetric grain boundaries and by reorganizing the atomic structure it was possible to successfully construct them at the desired misorientation angle. Strong out-of-plane corrugations may result from the presence of dislocations on the structure. This has been further investigated in [34]. By choosing the above-mentioned misorientation angles the gap between neighboring defects were smaller, which facilitates our calculations. More details can be found in [22].
In this study, we employed the non-equilibrium molecular dynamics (NEMD) method [35]. In order to avoid atomic discretion, the atoms are fixed at the beginning and at the end of the model. (Depicted as the black atoms in Figures 1 and 2). Then, we divided the atomic sheet into several strips and imposed heat conduction to atomic sheets. The temperature at each strip may be determined as: where, T i denotes the temperature of ith strip, N i is the number of atoms in ith strip, k B is the Boltzmann's constant, m j is atomic mass of atom jth and P j indicates the atomic momentum of atom j, respectively. h-BN and Graphene sheets were firstly relaxed at 300 K using the Nosé-Hoover thermostat (NVT) method [36] for 50 ps. Then, the first and last strips were subjected to a temperature difference of 20 K using the NVT method, while the other strips were governed at a desired temperature and simulated at constant energy level using the NVT method and NVE ensemble, respectively. Hot and cold sources were assigned between the first and twentieth slabs, respectively (Figures 1 and 2). After 200 ps the system reached (heat transfer) a steady-state condition including a transient temperature difference. After this stage the heat conduction was imposed to the models for 4800 ps, while the energy of the atoms at the hot strip increases slightly and simultaneously, the energy of the atoms at the cold strip decreases slightly. This way, the thermal properties can be computed using the average of temperature at each strip.

Obtained Results for Thermal Conductivity
To obtain more accurate results from the NEMD simulations, the total energy of the system should remain constant. In other words, the quantity of the increased and decreased energy in the hot and cold sources must be the same. Figure 4 explains the possible fluctuation in energy of the atoms computed at the hot-and cold-strips for both models as well as the pristine sheets. The slope value for both hot-strip and cold-strips are nearly identical by means of adding and removing energy from the hot and cold reservoirs constantly, the energy of the atoms is constant at the conducting strips. Figure 4 depicts the evolution of the energy added to the hot reservoir and removed from the cold reservoir per unit length of the h-BN and graphene sheets. The curves are following a linear pattern, which demonstrates that the total energy of the system remains constant during the analysis that in turn is caused by applying a constant heat flux given by where, q is the added or removed energy, A is the cross sectional area of the sheet and dt is the time increment. Furthermore, the temperature gradient can be calculated based on the slope of the dotted red line in Figure 5. It is worth noting that, for the sake of better accuracy of dT/dx, the results of the average temperature of the first and last three layers (as shown in Figure 5) were ignored.  Based on the calculated heat fluxes (q x ) and the temperature gradient ( dT dx ), the thermal conductivity (k) of h-BN, and Graphene sheets was calculated using the 1D form of Fourier law as follows: where q x is the heat flux, A cross is cross sectional area of the sheet, and dT dx is the temperature gradient. In Figure 5, it is assumed that the temperature remains constant at both ends of the lattice sheet. Hence, in a pristine lattice sheet, the temperature tends to drop gradually as the heat flows from the hotter end to the colder end. This assumption is confirmed by the correlation values obtained from the MD simulations, showing a linear flow of heat across the pristine lattice sheets.
We carried out MD simulations for models of lengths ranging from 33 to 140 nm. Figure 6 shows the inverse of computed thermal conductivity for h-BN and graphene models as a function of inverse length. The thermal conductivity pristine models with infinite length was calculated by fitting a linear relation. It can be clearly seen in Figure 6 that the thermal conductivity increases with increasing length of the sheets. In contrast, Figure 6 demonstrates that the thermal conductivity of the sheets decreases at higher temperature due to the increase in phonon-phonon scattering. As the temperature increases, the phonon lifetime drops and this leads to the convergence of the thermal conductivity at lower correlation times [13].  Table 1 represents the computed thermal conductivity of h-BN and Graphene models for 3 temperatures. It should be noted that the thermal conductivity of infinitely long sheets can be computed by extrapolating the results for finite lengths. Hence, the term effective phonon mean free path (Λ eff ) can be explained as: 1 where, L is the length of the pristine sheet. With attention to proportionality of thermal conductivity to Λ eff , the term thermal conductivity for an infinite sheet can be computed by extrapolating the values to 1/L = 0 [37][38][39]. Consequently, the effective thermal conductivity can be determined by the one-dimensional form of the Fourier law: where, L denotes the length of the cross-section, q is the heat flux of the system, and ∆T is the temperature difference along the laminate. The predicted thermal conductivity for a single-layer h-BN is 606.6 W/m· K, which is in good agreement with the predictions of 600 W/m·K in [13], and the measured experimental value of 390 W/m·K for in-plane thermal conductivity of bulk h-BN [14]. We obtained a thermal conductivity of single-layer Graphene of 2941.2 W/m·K, which agrees well with the reported experimental data ranging from 3000 to 5000 W/m· K at room temperature (300 K) [40][41][42]. Our predictions also agree well with other predictions in the literature; see e.g., the contribution in [38] reporting a value for the thermal conductivity of graphene around 2650 W/m·K. In addition, Hahn et al. predicted the thermal conductivity of 262 W/m·K for a single-crystalline graphene with the sheet length of 200 nm using approach-to-equilibrium molecular dynamics (AEMD) simulation [43].

Obtained Results for Thermal Conductance
The thermal conductance of h-BN and Graphene sheets along the six modeled grain boundaries depicted in Figure 3 was calculated based on the variation of the sheet length and temperature. The thermal conductance for each model takes the form: where, q x is the heat flux, A cross is the cross sectional area of the sheet and ∆T is the temperature difference along the interface of two sheets. The heat flux can be calculated using the slope of the energy curves as shown in Figure 7. As explained in the previous section, the slope of two curves must be close to each other since the total amount of energy must be constant during the analysis.  (5) can be calculated based on the slope of each line. The slope of two drawn linear curves must be the same because the energy added in the hot source is equal to the energy removed in cold source during the simulations. Figure 8 illustrates the temperature gradient for 20-layer at room temperature. For the sake of conciseness, only the temperature gradient for a sheet with a length of 33 nm is presented. The temperature difference (∆T) at the grain boundaries can be computed based on the temperature-position graph. A sudden drop in the temperature is seen along the transition regions of the two domains, which is denoted by the gap in the plotted area. The gap indicates the presence of the grain boundaries and confirms the non-linearity in the heat flow across a polycrystalline sheet. After the heat is transferred across the line of defects, the flow of heat is linear again. This is confirmed by their respective slopes, which shows a linear correlation. Kelvin. This temperature difference should be used in Equation (5) for the calculation of the thermal conductance. Figure 9 exhibits the computed thermal conductance for h-BN models for six grain boundaries at different temperatures and lengths. It can be seen in Figure 9 that that the term of thermal contact conductance for h-BN material is almost independent of temperature and length [13]. As can be seen from Figure 9, the thermal conductance for h-BN for different sheet lengths, temperatures, and grain boundaries varies between 11.1 to 30.3 GW/m 2 ·K, which agrees well with the reported values of 11 ± 1 GW/m 2 ·K in [13]. The insignificant influence of the temperature on the thermal conductance at grain boundaries can be explained as the heat transfer along the grain boundaries is dominated by the phonon-defect scattering. As expected, Figure 9 demonstrates that the variation of the sheet length does not influence the thermal conductance of h-BN films significantly. The maximum variation of thermal conductance for different sheet lengths with identical grain boundary and temperature is 16%. Figure 9. Predicted the term of thermal conductance of h-BN material for six grain boundaries, sheet lengths (330 A, 500 A, 670 A, 1000 A and 1400 A) and three temperatures (300 K, 500 K and 700 K). Figure 10 depicts the thermal conductance of graphene different grain boundaries in terms of different sheet lengths and temperatures. It can be seen in Figure 10 that the predicted thermal conductance of Graphene varies between 34.3 to 68.1 GW/m 2 ·K. The experimentally measured thermal conductance of graphene was reported around 133 GW/m 2 ·K [4], while the predictions vary from 15 to 45 GW/m 2 ·K [44] and from 0.14 to 0.58 GW/m 2 ·K [43] using the NEMD and AEMD approaches, respectively. In the case of identical grain boundary and temperature, the maximum variation thermal conductance of graphene material for five different sheet lengths is around 8%. For the sake of better understanding, the variation of thermal conductance for the six studied grain boundaries in terms of temperature is illustrated in Figure 11. In Figure 11, the term of thermal conductance was calculated using the same methodology for computing the thermal conductivity (cf. Equation (7)). In Figure 11a,b, the thermal conductance for six employed grain boundaries varies between 14.3 to 31.5 GW/m 2 ·K and 38.5 to 69.2 GW/m 2 ·K for h-BN and Graphene, respectively.  Figure 11 demonstrates that the thermal conductance of h-BN and graphene is independent of the temperature and sheet length. However, small discrepancies in the predicted thermal conductivities for h-BN films can be observed in Figure 11a.

Summary and Conclusions
In this study, we performed numerical simulations to study the influence of temperature and sheet length on thermal properties (i.e., thermal conductivity and thermal conductance) of two materials; i.e., hexagonal Boron Nitride (h-BN) and graphene. Six grain boundaries with different defect concentrations were used to investigate the thermal conductance along hBN-hBN and Graphene-Graphene interfaces. The non-equilibrium molecular dynamics method (NEMD) was employed to simulate and analyze the thermal properties. To verify our results, we first calculated the thermal conductivity of pristine h-BN and Graphene models and compared these values with the previous results from the literatures. The computed thermal conductivities of h-BN and Graphene sheets were close to the measured thermal conductivities of pristine h-BN and Graphene. The results show that by increasing the sheet lengths, the thermal conductivity in h-BN and Graphene models increases and by increasing the temperature, the thermal conductivity decreases. On the other hand, the thermal conductance of the material is independent of the sheet length and temperature.

Conflicts of Interest:
The authors declare that there is no conflict of interest regarding the publication of this paper.