Nanoscale
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
COMMUNICATION
Cite this: DOI: 10.1039/c4nr07359b Received 12th December 2014, Accepted 16th March 2015 DOI: 10.1039/c4nr07359b
View Article Online View Journal
The unexpected non-monotonic inter-layer bonding dependence of the thermal conductivity of bilayered boron nitride Yufei Gao,†a Xiaoliang Zhang,†b Yuhang Jingc and Ming Hu*b,d
www.rsc.org/nanoscale
Hexagonal boron nitride (BN) and its bilayer form are very fascinating two-dimensional materials that have attracted tremendous interest recently. Their realistic applications in emerging nanoelectronics usually quest for manipulating the thermal transport properties in a precise manner. Using nonequilibrium molecular dynamics simulations, we herein studied the effect of inter-layer covalent bonding on the thermal conductivity of bilayered BN. We found that the in-plane thermal conductivity of bilayered BN, which can be largely tuned by introducing covalent bonding between the two BN layers, depends not only on the inter-layer bonding density, but also on the detailed topological configuration of the inter-layer bonds. For randomly distributed inter-layer bonding the thermal conductivity of bilayered BN decreases monotonically with inter-layer bonding density, the same behavior already found for bilayered graphene. However, for regularly arranged inter-layer bonding the thermal conductivity of bilayered BN surprisingly possesses a non-monotonic dependence on the inter-layer bonding density. This non-intuitive non-monotonic dependence is further explained by performing spectral energy density analysis, where the peak and valley values of the thermal conductivity are governed by different mechanisms. These results suggest the application of inter-layer covalent bonding in designing nanoscale devices with precisely tunable thermal conductivities.
Recently, there has been increased interest in the investigation on the physical and chemical properties of multilayered nanostructures, such as multilayered graphene and boron nitride (BN).1–3 A common feature of these materials is their layered
a School of Architecture & Civil Engineering, Shenyang University of Technology, Shenyang 110870, China b Institute of Mineral Engineering, Division of Materials Science and Engineering, Faculty of Georesources and Materials Engineering, Rheinisch-Westfaelische Technische Hochschule (RWTH Aachen University), 52064 Aachen, Germany. E-mail: [email protected] c Department of Astronautical Science and Mechanics, Harbin Institute of Technology, Harbin 150001, China d Aachen Institute for Advanced Study in Computational Engineering Science (AICES), RWTH Aachen University, 52062 Aachen, Germany † These authors contributed equally to this work.
This journal is © The Royal Society of Chemistry 2015
structures with high anisotropy, namely the atoms within the layers possess strong short-ranged covalent bonding and the interactions between different layers are described by weak long-ranged van der Waals forces.4 The van der Waals forces play an important role in making the layers maintain Bernal stacking, where the atoms in different layers do not eclipse each other. The previous research has demonstrated that the multilayered structures possess some fascinating properties compared to monolayered structures.5 Multilayered boron nitride is an analogue of graphite, in which the carbon atoms are replaced by boron and nitride atoms. The structure prefers to form a multilayered tubular structure or multilayered graphite-like flakes. Interestingly, the boron–nitride bonds in the in-plane possess a partially ionic character, which is quite different from that of graphene, whose bonds are completely covalent.6 Therefore, there are some distinct differences between these two structures. For example, the BN structures exhibit profound chemical and thermal stabilities.7,8 Moreover, multilayered BN structures are more promising for practical applications, such as optoelectronic technologies, tunnel devices, and field effect transistors.9–12 Nowadays, multilayered BN has been proposed as a potential material for thermal management applications,13–15 due to its high thermal conductivity and sensitivity to isotopic substitution. The multi-layered BN sheets and nanoribbons have even been foreseen as materials for designing microprocessors. Recently, the strong bonds between multi-walled carbon nanotubes and graphite layers can be created by irradiating the material with electron beam and femto-second laser excitation.16,17 Using similar ultraviolet laser irradiation, recently Komatsu et al.18 experimentally demonstrated that the structural relaxation of the sp2-bonded “metastable” phase of BN can be induced into more stable sp3-bonded phases at relatively low temperatures and below atmospheric pressure, which appears to be consistent with the recent pressure–temperature phase diagram of BN. Therefore, we speculate that the laser irradiation technique could also be used to create sp3bonds at the interface between bilayered BN. The inter-layer
Nanoscale
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Communication
bonds could cause strong cross-plane couplings between the layers and thus change the thermal transport properties of the materials significantly.19 Recently, the influence of inter-layer bonding in bilayered graphene on the thermal transport properties has been reported.19 The results showed that 5% bonding density monotonically leads to 75% decrease in the thermal conductivity, indicating that the thermal transport properties of bilayered graphene can be effectively controlled by modifying the bonding. Therefore, it is intuitive to anticipate that the emergence of the inter-layer bonds could also provide a controllable mechanism for tuning the physical properties of multilayered BN. In this paper, by performing nonequilibrium molecular dynamics simulation (NEMD) we studied the effect of interlayer covalent bonding on the thermal transport properties of bilayered BN. We demonstrate that the in-plane thermal conductivity of bilayered BN can be largely tuned by introducing covalent bonding between the two BN layers. However, the bonding density dependent thermal conductivity of bilayered BN exhibits unexpectedly different features, depending on the detailed topological configuration of the inter-layer bonds, which is counter-intuitive and has not been found in the literature so far. The underlying mechanism is provided by performing vibrational density of states and spectral energy density analysis. Our model system is armchair bilayered boron nitride sheets, consisting of m × n supercells, corresponding to a simulation box size of ∼2 × 540 nm2, with each supercell constructed by a bilayer honeycomb lattice and composed of 160 boron and nitride atoms. The structure of a supercell is shown in Fig. 1. The initial distance between the two BN layers is 0.33 nm. The inter-layer covalent bonding process is realized by the following steps: (1) the bonded atoms are chosen and moved closer to each other along the out-of-plane direction until the distance between the two bonded atoms reaches 0.145 nm (the equilibrium length of the B–N bond). (2) The period boundary condition is applied in the lateral (in-plane) directions. Initially, the bonded atoms are fixed and all other atoms in the system are relaxed with the NPT (constant temperature and constant pressure) ensemble for 100 000 steps. (3) The fix on the bonded atoms is removed and the entire system is relaxed with the NPT ensemble for additional 1 000 000 steps. Finally, the bilayered BN with inter-layer bonding is achieved. The interactions between boron and nitride atoms in the in-plane honeycomb structure and in the inter-layer bonding configuration are described by the Tersoff potential.20,21 The Tersoff potential with the optimized parameters has been proven suitable to investigate the properties of carbon, silicon, boron and nitride based nanomaterials.22–24 The interactions between boron and nitride atoms, that belong to different layers and are not covalently bonded, are presented with van der Waals forces. The standard LJ potential is used and the parameters can be found in ref. 25. It should be noted that the interfacial bonding created between BN layers should belong to the sp3-bonding configuration. In order to test whether the
Nanoscale
Nanoscale
Fig. 1 (a) Top and side views of bilayered BN without inter-layer covalent bonding. Color coding: light pink, boron; blue, nitride. From (b) to (d): a supercell of bilayered BN with inter-layer covalent bonding density 6.25%, 8.75%, 12.5% (left: randomly distributed bonding; right: regular bonding). The covalently bonded atoms are highlighted in green.
empirical Tersoff potential used is suitable for describing such sp3-bonding, we ran additional benchmark cases for bulk cubic BN (zinc blende structure) and wurtzite BN at room temperature, where boron and nitride atoms are explicitly bonded with sp3 bonding. The results of lattice constants and bulk modulus compared with experimental and ab initio calculations are shown in Table 1. It can be seen that the empirical
Table 1 Lattice constants and bulk modulus of cubic BN (zinc blende structure) and wurtzite BN calculated using the empirical Tersoff potential. The experimental and ab initio calculation results are also shown for comparison
Cubic BN Wurtzite BN
Tersoff potential Experiments or ab initio Tersoff potential Experiments or ab initio
Lattice constant, a (Å)
Lattice constant, c (Å)
Bulk modulus (GPa)
3.50 3.6226
— —
453.5 369,27 387.328
2.45 2.5529
4.11 4.2229
460.5 39730
This journal is © The Royal Society of Chemistry 2015
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Nanoscale
Communication
Tersoff potential, although originally optimized for planar sp2BN structures such as two-dimensional BN and BN nanotubes,20,21 can also simulate the bulk sp3-BN structures. This proves that the Tersoff potential can be effectively used for describing the interfacial sp3-like BN bonds in our model structures. All MD simulations herein are performed using the LAMMPS package31 with timestep 0.5 fs throughout. In our simulation, the velocity-Verlet algorithm is employed to integrate the equations of motion and the temperature is controlled by employing the Nosé–Hoover thermostat.32 First, we relaxed the entire system with the NPT ensemble for 50 ps, during which the pressure and temperature are kept at 1 atm and 300 K respectively. After that, we continued to relax the system with NVT (constant temperature and constant volume) and NVE (constant volume without a thermostat) ensembles, respectively. With the above processes, the system can reach the equilibrium state and a stable configuration can be obtained. Following equilibration, we used NEMD to calculate the thermal conductivity of bilayered BN. The constant heat flux was imposed by the Müller–Plathe method.33 During the NEMD simulation, the system is divided into slices along the longitudinal (heat flux) direction. The outmost two slices on the left and right side are set to be “cold” and “hot” regions, respectively. Then, the hottest atom in the “cold” region and the coldest atom in the “hot” region are selected and their energies are exchanged every 600–2000 steps, depending on the thermal conductivity of the system. The thermal conductivity of bilayered BN was calculated by Fourier’s law κ¼
JQ ; rT
ð1Þ
where JQ is the heat flux in the longitudinal (y) direction defined as JQ = (dQ/dt )/A, ∇T is the temperature gradient calculated by the linear fitting to the central part of the temperature profile of the system along the longitudinal direction, A is the cross-sectional area which is calculated as w·d, where w and d are the width and thickness of the bilayered BN, respectively. The thickness of monolayered BN is assumed to be 3.35 Å and the thickness of bilayered BN is assumed to be 6.7 Å. The total NEMD simulation time is 10 ns (2 × 107 steps). Normally, the temperature profile reaches steady state after 1–2 ns. Once the temperature profile becomes stable, we start to collect the heat flux and the temperature gradient to calculate the thermal conductivity every 100 ps. The reported thermal conductivity is averaged over all the data points collected during the steady state run. We have checked all the structures after the NPT and NVT relaxation runs and after imposing a constant heat flux during the NEMD simulation as well. We did not find a significant change in the structures, especially the interfacial covalent bonds created from the beginning remained until the end of the NEMD run. In addition, for selected cases we calculated the bond lengths for the sp3 hybridized (cross-plane) BN atoms and the regular sp2-bonded (in-plane) atoms to be about 1.491 Å and 1.444 Å, respectively. The in-plane BN bond
This journal is © The Royal Society of Chemistry 2015
length for bilayered BN is very close to that for monolayered BN calculated by the same classical potential (1.446 Å), both of which are in very good agreement with that obtained from ab initio calculations (1.446 Å).34,35 The trend that the bond length of sp3 bonding is slightly larger than that for sp2 bonding is also consistent with that found from ab initio calculation.28,35 Again, the results fully support that the optimized Tersoff potential is suitable to describe the pseudo-sp3 hybridized bonds and the honeycomb two-dimensional sp2 bonds as well. The relationship between the thermal conductivity of bilayered BN and the inter-layer bonding density is shown in Fig. 2. Here we consider two different configurations of interlayer bonding: regular and random distribution, with typical model structures shown in Fig. 1. The random model is the one used in the previous investigation on bilayered graphene.19 For comparison we also show the thermal conductivity of monolayered BN and bilayered BN without inter-layer covalent bonding, i.e. the inter-layer interaction is pure van der Waals forces. We first notice that generally there is a large decrease in the thermal conductivity of bilayered BN with the increase of bonding density, which is consistent with a previous study on bilayered graphene,19 where there is 75% decrease in thermal conductivity with 5% bonding density. The reduction is generally understood as the consequence of remarkable phonon coupling (scattering) between the two layers when strong covalent bonds are introduced at the interface. A more distinct feature in Fig. 2 is that in contrast to the monotonic dependence of the thermal conductivity of bilayered BN with randomly distributed interfacial bonding, which was also found for bilayered graphene,19 the thermal conductivity of bilayered BN with regular interfacial bonding exhibits a non-monotonic dependence on the bonding density.
Fig. 2 Dependence of the thermal conductivity of bilayered BN on the inter-layer covalent bonding density. The horizontal solid line represents the thermal conductivity of monolayered BN. The zero bonding density denotes bilayered BN with only van der Waals interactions between the layers.
Nanoscale
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Communication
Nanoscale
In general the thermal conductivity in the case of the regular bonding configuration (see the typical structure in the right panel of Fig. 1) is higher than that for random bonding, except for the two cases of the bonding density of 2.5% and 8.75%. For regular bonding the thermal conductivity has two peak values at the bonding density of 6.25% and 12.5%, where the thermal conductivity of the bilayered BN between regular and random bonding can differ 2–3 fold. These results indicate that when introducing covalent bonds at interfaces, the thermal conductivity of the bilayered BN not only depends on the overall bonding density, but also is very sensitive to the detailed topological arrangement of the interfacial bonding. Throughout our simulation we are aware that the different bonding configurations might lead to different thermal conductivities, especially for randomly distributed bonds. Therefore, in our paper we ran five different bonding configurations for each random bonding case in Fig. 1(b)–(d). All thermal conductivity results reported in Fig. 2 are the averaged value of these five deferent runs. In addition, we ran additional simulations with the bonding density of 6.25% (500 nm) and the five covalent bonds taking up two lines (in this case three covalent bonds form one line and the two other bonds form the other line). The thermal conductivity of this new configuration is 139.9 W mK−1. Compared with the results in Fig. 2, we can find that this value is very close to the average thermal conductivity of the bilayered BN with random interlayer bonding (130.9 W mK−1) but is far lower than that for regular interlayer bonding (275.7 W mK−1). Therefore, it should be noted that changing the saturated bond line in Fig. 1(b) into two unsaturated lines results in a significant drop in thermal conductivity, as the unsaturated bond lines result in lower phonon group velocities and also always have a very strong interaction between the lines and thus the phonon scattering is enhanced, the same mechanism as that for Fig. 1(d), as we will see later. We also considered the symmetric distribution of the interlayer bonding. We run additional cases of 4, 6, and 8 bonds that take up symmetrical positions but form two unsaturated lines, corresponding to the bonding density of 5%, 7.5%, and 10%, respectively. The thermal conductivity of bilayered BN with regular, symmetric, and random inter-layer bonding is compared in Table 2. From Table 2 we can see that for the bonding density of 5% and 7.5% the thermal conductivity of bilayered BN with sym-
Table 2 Comparison of the thermal conductivity of bilayered BN between different inter-layer bonding distributions: regular, symmetric, and random. The thermal conductivity is in the units of W mK−1
metric inter-layer bonding is lower than that for regular bonding. This trend is opposite for the case with the bonding density of 10%, which can be attributed to the less interference between the two bonding lines in the case of symmetric distribution. Nevertheless, the results clearly show that for all cases the thermal conductivity of bilayered BN with symmetric inter-layer bonding is higher than the corresponding randomly distributed cases. Since the symmetric distribution can also be regarded as regular distribution, once again the results in Table 2 confirm the trend in Fig. 2 that the thermal conductivity of bilayered BN strongly depends on the detailed topological configuration of the inter-layer bonds and the regular inter-layer bonding always yields higher thermal conductivity. The symmetric cases also confirm the nonmonotonic dependence in Fig. 2, when taking the 12.5% regular case shown in Fig. 1 as the symmetric distribution. Moreover, combining Fig. 1 and 2 we can see that the following three cases of inter-layer bonding arrangement can lead to a significant increase in thermal conductivity: (i) for low bonding density (2.5%–5%) the inter-layer bonds arrange only in one line but certainly cannot form one full line, i.e. only one unsaturated line can form at the interface; (ii) for higher bonding density (6.25%–12.5%) the inter-layer bonds form in one full (saturated) line along the heat flux direction and the rest try to fill up the other line; (iii) the inter-layer bonds form two symmetric but unsaturated lines. Other inter-layer bonding arrangements that do not fulfill any requirement of the above three cases has much lower thermal conductivity. It is well known that the thermal conductivity of many nanosystems calculated by the NEMD method exhibits size dependence. To address this issue, for selected cases with the bonding density of 8.75% and 12.5%, corresponding to the typical minima and maxima in Fig. 2, we ran additional simulations with the length doubled, i.e. 1078 nm. For each bonding density, we generated three different random bonding configurations and one regular bonding configuration, so the finally reported value for the random bonding configuration is the average of these three independent runs. The comparison of the calculated thermal conductivity between systems with lengths of 540 nm and 1078 nm is given in Table 3. We found that the thermal conductivity of bilayered BN only has a slight incremental change when the length is doubled. More importantly, the trend is exactly the same, i.e.
Table 3 Comparison of the thermal conductivity of bilayered BN with inter-layer bonding calculated by the NEMD method between systems with lengths of 540 nm and 1078 nm. The calculated thermal conductivity is in the units of W mK−1
Length (nm) 540
Distribution Bonding density
Regular
Symmetric
Random
5.0% 7.5% 10%
201.9 168.0 131.5
161.3 149.1 155.2
135.9 127.2 98.4
Nanoscale
1078
Bonding density
Random bonding
Regular bonding
Random bonding
Regular bonding
8.75% 12.5%
94.8 93.6
112.9 222.2
118.9 108.9
156.0 220.5
This journal is © The Royal Society of Chemistry 2015
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Nanoscale
for the bonding density of 8.75% the thermal conductivity is still insensitive to the detailed bonding configuration, while for the bonding density of 12.5% the thermal conductivity for the regular bonding configuration is at least twice as large as that for the random bonding configuration. Moreover, our spectral energy density (SED) analysis of the above two bonding density cases yields a similar trend to that found in NEMD simulation (see the results later). Since SED is based on equilibrium molecular dynamics (EMD) simulation, the size effect in SED analysis is minor. Therefore, we believe that the nonmonotonic dependence of the thermal conductivity on the bonding density revealed by NEMD simulation is size independent. We also run more NEMD simulations for the selected cases of bilayered BN with a fully random bonding density of 3.75%, 6.25%, and 8.75%. We created the fully randomly distributed interfacial covalent bonds by randomly selecting bonded atoms in the entire simulation domain, i.e. ∼2 × 540 nm2. The results show that, compared to the original construction method by duplicating the supercell containing irregularly aligned bonded atoms, the thermal conductivity of the bilayered BN with fully interfacial random bonding is only slightly lower. Nevertheless, the trend remains unchanged. Once again, it proves that the thermal transport of bilayered BN with random interfacial covalent bonding is not sensitive to the detailed bonding arrangement. The intuitive phenomenon of the reduction in the thermal conductivity of bonded bilayered BN in contrast to free bilayered BN can be understood in terms of enhanced phonon coupling (scattering) between the two bonded layers, as evidenced by the phonon density of state (PDOS) of BN atoms in our model structure. The PDOS was calculated by imposing a Fourier transform on the atomic velocity autocorrelation function.36 The PDOS values between the bonded and non-bonded atoms are compared in Fig. 3. Both PDOS values are similar for frequencies below 30 THz. However, there are distinct differences in the frequency between 30 and 50 THz. The high
Communication
frequency peak in the PDOS of the bonded atoms is approximately at 37 THz and for non-bonded atoms the peak is at 49 THz. That is to say, there is a significant redshift for the PDOS of the bonded atoms, which means the probability of phonon scattering increases, provided that the total number of phonon modes is exactly the same for both regular and random bonding structures. Therefore, the redshift of the high frequency phonons of the bonded atoms leads to a significant reduction of the thermal conductivity of the bonded structure. We also examined the effect of cross-linking on the contribution of the in-plane and out-of-plane phonons in bilayered BN. To this end, we quantified the relative contributions of the lattice vibrations in the x, y, and z directions to the total heat flux of the bilayered BN, corresponding to the in-plane longitudinal modes, in-plane transverse modes, and the out-ofplane flexural modes, respectively, by using the NEMD-based method we proposed recently37 J A!B;α ¼
1 XX F ijα ðυiα þ υjα Þ; 2A i[A j[B
ð2Þ
where α is the x, y, or z direction, A and B are the left and right side of a virtual interface located at the middle of the bilayered BN in the x (heat flux) direction, JA→B,α is the heat flux from A to B due to lattice vibrations in the α direction, A is the crosssectional area, Fijα is the component of the force acting on atom i due to atom j, and υiα is the α component of the velocity of atom i. Details of calculation in the form of the three-body Tersoff potential can be found in ref. 37. We performed the polarization analysis of three typical cases of bonding densities: 6.25%, 8.75%, and 12.5%. The relative contribution ( percentage) to the total heat flux from vibrations in x (longitudinal), y (transverse), and z (out-of-plane) directions for these three cases is compared in Table 4. It is clearly seen that for the cases of random interfacial bonding the relative contribution from in-plane and out-ofplane vibrations is nearly the same. For the case of the bonding density of 8.75%, the relative contribution percentage is also the same as those for random bonding. This is consistent with the results of their similar thermal conductivity. Comparing the results between regular and random bonding, we find that when the interfacial covalent bonds form one or two saturated lines, the interfacial bonds affect the in-plane phonon transport significantly, i.e. the contribution from the in-plane longitudinal modes is reduced while the contribution
Table 4 Comparison of the relative contribution ( percentage) to the total heat flux from vibrations in x (longitudinal), y (transverse), and z (out-of-plane) directions between different bonding densities
Fig. 3 Comparison of the phonon density of states between bonded atoms and non-bonded atoms in the structures of bilayered BN.
This journal is © The Royal Society of Chemistry 2015
6.25% 8.75% 12.5%
Regular interfacial bonding
Random interfacial bonding
x
y
z
x
y
z
45.4% 59.9% 44.4%
26.8% 27.7% 29.4%
27.8% 12.7% 26.2%
56.2% 56.5% 57.2%
28.6% 29.5% 28.7%
15.2% 14.0% 14.1%
Nanoscale
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Communication
Nanoscale
from the out-of-plane flexural modes is enhanced (highlighted in Table 4) and the contribution from the in-plane transverse modes remains the same. This is understandable considering that when the crosslinking covalent bonds are formed regularly at the interface, the atoms in different layers can vibrate in unison in the out-of-plane direction, and in this way the out-of-plane flexural modes have less phonon scattering, which results in higher thermal conductivity, since usually the flexural modes in graphene-like two-dimensional materials conduct heat more efficiently. We now explain the non-intuitive phenomenon that the regular bonding configuration results in substantially higher thermal conductivity than the random bonding configuration with exactly the same bonding density, e.g. for the bonding density of 12.5% in Fig. 2. To this end, we calculated the frequency dependent phonon lifetime by performing spectral energy density analysis. The phonon normal modes are obtained by38 rffiffiffiffiffiffi X mj ˙~ ~ e*j ð~ k; νÞexpð2πi~ k~ r l Þ: ð3Þ υjl ðtÞ~ Qð k; ν; tÞ ¼ N jl Then the spectral energy density is calculated by39 2 ð ˙~ Φð~ k; ν; f Þ ¼ Qð k; ν; tÞexpð2πiftÞdt ;
ð4Þ
and the phonon lifetime (relaxation time) is obtained by fitting the Lorentzian function.40 Finally, we investigated the mode dependent thermal conductivity of the bilayered BN, which can be derived from the phonon Boltzmann transport equation (BTE) under the relaxation time approximation41 XX κα ¼ cph υg;α 2 ð~ k; νÞτð~ k; νÞ; ð5Þ ~ k
ν
where ~ k is the wave vector in the first Brillouin zone, ν is the phonon branch, κα denotes the thermal conductivity in the αdirection, cph is the phonon volumetric specific heat calculated kB using cph ¼ , kB is the Boltzmann constant, V is the system V volume (a thickness of 3.35 Å and 6.7 Å is assumed for monolayered and bilayered BN, respectively), υg,α is the α-component of the phonon group velocity, and τ is the phonon lifetime. The frequency dependent phonon lifetimes of bilayered BN (inter-layer bonding density of 12.5%) with random and regular bonding configurations are compared in Fig. 4(a). We can see that the frequency dependent phonon lifetime of bilayered BN for both regular and random bonding roughly follows the power law of ω−α (α ∼ 0.71). The phonon lifetime in bilayered BN with the regular bonding configuration is only slightly higher than that with the random bonding configuration, which obviously cannot explain the three times higher thermal conductivity that occurs for the regular bonding configuration (Fig. 2). We also noticed that there is a significant increase in the lifetime of low frequency phonons (below 1 THz), but still this is not the main reason for the large thermal conductivity enhancement. Then, we further calcu-
Nanoscale
Fig. 4 Comparison of (a) phonon lifetime, (b) phonon group velocity, and (c) frequency dependent accumulated thermal conductivity of bonded bilayered BN between random and regular bonding configurations (inter-layer bonding density 12.5%). The dashed line is the power law fitting (ω−α) to the frequency dependent phonon lifetime.
lated the square of the phonon group velocity (υg,α2) in eqn (5) as a function of frequency, as shown in Fig. 4(b). Now it is clearly seen that the phonon group velocity for the case of regular inter-layer bonding is much higher than that for random bonding. According to eqn (5), finally we can plot the accumulated thermal conductivity of the bonded bilayered BN as a function of frequency, as shown in Fig. 4(c). Initially the accumulated thermal conductivity is almost the same for both random and regular bonding, when the phonon frequency is less than 2 THz, which means that the low frequency phonons are not sensitive to the detailed interfacial bonding configuration. This is understandable considering that the low frequency phonons usually have a considerably long wavelength, which is much larger than the characteristic distance between the bonded sites, thus they usually travel through the struc-
This journal is © The Royal Society of Chemistry 2015
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Nanoscale
tures without scattering. Above 2 THz, the thermal conductivity for the case of regular bonding increases dramatically until the phonon frequency reaches around 30 THz, while there is only a small increment in the thermal conductivity of bilayered BN with random inter-layer bonding. For even higher frequency phonon modes, both regular and random bonding show a minor increase in the overall thermal conductivity, suggesting that the very high frequency phonon modes do not contribute much to the thermal transport in bilayered BN. This is reasonable considering that the high frequency phonon modes are usually localized and propagate diffusively with very low phonon relaxation time [Fig. 4(a)]. Combining Fig. 4(a)–(c) we can conclude that the much higher thermal conductivity of the bilayered BN with regular bonding as compared with random bonding is mainly due to the large enhancement in the group velocity of the phonon modes in the medium frequency range (2 to 30 THz). We also ran SED simulations for the bonding density of 8.75% with regular and random bond patterns. The results are not reported here for brevity. First of all, the thermal conductivity of bilayered BN with random and regular bonding calculated by the SED method is 75.7 and 92.2 W mK−1, respectively, which are very close to the NEMD simulation results (94.8 and 112.9 W mK−1, respectively), as shown in Fig. 2. More importantly, the thermal conductivity trend predicted by SED calculation is in accord with that found from NEMD simulations, i.e. for the bonding density of 8.75%, the thermal conductivity of bilayered BN is insensitive to the bond pattern. This is because, as revealed by the SED calculations, for the bonding density of 8.75% there is no significant difference in both the phonon mode specific group velocity and lifetime between regular and random patterned interfacial bonds. Before closing, we explain the phenomenon that the thermal conductivity of bilayered BN with regular inter-layer bonding sometimes reaches a minimum at certain inter-layer bonding densities, i.e. the non-monotonic inter-layer bonding dependence is shown in Fig. 2. Taking bonding density of 8.75% as an example, we found that the thermal conductivity is only slightly higher than that for randomly distributed bonding and is significantly lower than that for the bonding density of 6.25% and 12.5%. By comparing the SED results of the phonon group velocity and phonon lifetime between the regular bonding density of 8.75% and 12.5%, we found that the significantly higher thermal conductivity for the regular bonding density of 12.5% than that for regular 8.75% stems from two factors: (1) the first factor, also the dominant factor, is that the phonon group velocities are largely enhanced, indicating that the formation of the full lines of the interfacial covalent bonds results in higher phonon group velocities. (2) “Interference” between bond lines: for the bonding density of 8.75% the inter-layer bonds only form one full line (we call it saturated line) along the heat flux direction and the rest do not fill up the other line (unsaturated line), as shown by the atomic structure in Fig. 1(c). In contrast, for the cases of 6.25% and 12.5% (regular bonding) all the inter-layer bonds either line up in one saturated line [6.25%, Fig. 1(b)] or form two
This journal is © The Royal Society of Chemistry 2015
Communication
saturated lines [12.5%, Fig. 1(d)]. The different bonding configurations between the saturated line and the unsaturated line for the bonding density of 8.75% will induce strong phonon scattering between the bonded atoms (reduced phonon lifetime) and therefore the heat conduction is significantly hindered. These results indicate that when controlling thermal transport of bilayered BN by inter-layer covalent bonds, one should not only consider the interfacial bonding density, but also carefully arrange the topological configuration of the inter-layer bonds. Finally, it is worth pointing out that the nonmonotonic behavior of the thermal conductivity on the density of regularly distributed interlayer bonds is not a particular property of boron nitride. We have also started the same calculation procedure as presented in this paper for bilayered graphene using the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential,42 considering the similarity of their honeycomb structures. Our preliminary results for two typical interlayer bonding densities of 6.25% and 12.5% show the same trend as that found for bilayered BN. Currently, this work is ongoing and the results will be submitted for publication separately. In summary, the effect of inter-layer covalent bonding on the thermal transport properties of bilayered BN is investigated by performing nonequilibrium molecular dynamics simulations. The in-plane thermal conductivity of bilayered BN can be largely tuned by introducing covalent bonding between the two BN layers. For randomly distributed inter-layer bonding the thermal conductivity of bilayered BN is reduced monotonically with the inter-layer bonding density. However, the thermal conductivity unexpectedly shows a non-monotonic dependence for regularly arranged inter-layer bonding. The mechanism underlying the counter-intuitive non-monotonic dependence is identified by performing spectral energy density analysis. We demonstrate that the peak values of the thermal conductivity of bonded bilayered BN with regular bonding are mainly due to the large enhancement in the group velocity of the phonon modes in the medium frequency range and the phonon lifetime augmentation only plays a minor role, while the valley values result from the low phonon group velocities and strong phonon scattering between the saturated and unsaturated bond lines. In short, when introducing chemical covalent bonding between the layers, the thermal conductivity of the bilayered BN not only depends on the inter-layer bonding density, but also is very sensitive to the topological configuration of how the inter-layer bonds are arranged. Our findings provide a guide to precisely modulate the thermal transport properties of bilayered BN with chemical engineering and they may be of use in tuning the electronic and optical properties for electronic, thermoelectric, photovoltaic, and opto-electronic applications.
Acknowledgements This research work was supported by the National Natural Science Foundation of China under grant no. 11272099 and
Nanoscale
View Article Online
Communication
11304059. Simulations were performed with computing resources granted by the Jülich Aachen Research Alliance-High Performance Computing (JARA-HPC) from RWTH Aachen University under project no. jara0104.
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
References 1 S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau and A. A. Balandin, Nat. Mater., 2010, 9, 555. 2 W. Q. Han, L. J. Wu, Y. M. Zhu, K. Watanabe and T. Taniguchi, Appl. Phys. Lett., 2008, 93, 223103. 3 D. Pacile, J. C. Meyer, C. O. Girit and A. Zettl, Appl. Phys. Lett., 2008, 92, 133107. 4 K. Gordiz, S. M. V. Allaei and F. Kowsary, Appl. Phys. Lett., 2011, 99, 251901. 5 D. Golberg, Y. Bando, Y. Huang, T. Terao, M. Mitome, C. C. Tang and C. Y. Zhi, ACS Nano, 2010, 4, 2979. 6 H. M. Ghassemi and R. S. Yassar, Appl. Mech. Rev., 2010, 63, 020804. 7 C. Y. Zhi, Y. Bando, C. C. Tang and D. Golberg, Mater. Sci. Eng., R, 2010, 70, 92. 8 J. S. Wang, C. H. Lee and Y. K. Yap, Nanoscale, 2010, 2, 2028. 9 K. Watanabe, T. Taniguchi and H. Kanda, Nat. Mater., 2004, 3, 404. 10 L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A. S. Mayorov, N. M. R. Peres, A. H. Castro Neto, J. Leist, A. K. Geim, L. A. Ponomarenko and K. S. Novoselov, Nano Lett., 2012, 12, 1707. 11 G. H. Lee, Y. J. Yu, C. Lee, C. Dean, K. L. Shepard, P. Kim and J. Hone, Appl. Phys. Lett., 2011, 99, 243114. 12 J. Beheshtian, A. Sadeghi, M. Neek-Amal, K. H. Michel and F. M. Peeters, Phys. Rev. B: Condens. Matter, 2012, 86, 195433. 13 T. Ouyang, Y. Chen, Y. Xie, K. Yang, Z. Bao and J. Zhong, Nanotechnology, 2010, 21, 245701. 14 D. A. Stewart, I. Savic and N. Mingo, Nano Lett., 2009, 9, 81. 15 C. W. Chang, A. M. Fennimore, A. Afanasiev, D. Okawa, T. Ikuno, H. Garcia, D. Li, A. Majumdar and A. Zettl, Phys. Rev. Lett., 2006, 97, 085901. 16 B. Peng, M. Locascio, P. Zapol, S. Y. Li, S. L. Mielke, G. C. Schatz and H. D. Espinosa, Nat. Nanotechnol., 2008, 3, 626. 17 J. Kanasaki, E. Inami, K. Tanimura, H. Ohnishi and K. Nasu, Phys. Rev. Lett., 2009, 102, 087402. 18 S. Komatsu, K. Kobayashi, Y. Sato, D. Hirano, T. Nakamura, T. Nagata, T. Chikyo, T. Watanabe, T. Takizawa, K. Nakamura and T. Hashimoto, J. Phys. Chem. C, 2010, 114, 13176.
Nanoscale
Nanoscale
19 A. Rajabpour and S. M. Vaez-Allaei, Appl. Phys. Lett., 2012, 101, 053115. 20 C. Sevik, A. Kinaci, J. B. Haskins and T. Cagin, Phys. Rev. B: Condens. Matter, 2011, 84, 085409. 21 A. Kinaci, J. B. Haskins, C. Sevik and T. Cagin, Phys. Rev. B: Condens. Matter, 2012, 86, 115410. 22 Y. F. Gao, Y. H. Jing, Q. Y. Meng, L. Zhang, J. Q. Liu and X. Qin, Phys. Status Solidi B, 2012, 249, 1728. 23 Y. H. Jing, M. Hu and L. C. Guo, J. Appl. Phys., 2013, 114, 153518. 24 M. Hu, X. L. Zhang and D. Poulikakos, Phys. Rev. B: Condens. Matter, 2013, 87, 195417. 25 B. Fakrach, A. Rahmani, H. Chadli, K. Sbai, M. Benhamou, M. Bentaleb, J.-L. Bantignie and J.-L. Sauvajol, J. Phys.: Condens. Matter, 2012, 24, 335304. 26 G. L. Doll, J. A. Sell, C. A. Taylor and R. Clarke, Phys. Rev. B: Condens. Matter, 1991, 43, 6816. 27 E. Knittle, R. M. Wentzcovitch, R. Jeanloz and M. L. Cohen, Nature, 1989, 337, 349. 28 B. Wen, J. Zhao, R. Melnik and Y. Tian, Phys. Chem. Chem. Phys., 2011, 13, 14565. 29 A. Nagakubo, H. Ogi, H. Sumiya, K. Kusakabe and M. Hirao, Appl. Phys. Lett., 2013, 102, 241909. 30 K. Kim, W. R. L. Lambrecht and B. Segall, Phys. Rev. B: Condens. Matter, 1996, 53, 16310. 31 S. Plimpton, J. Comput. Phys., 1995, 117, 1. 32 S. Nosé, J. Chem. Phys., 1984, 81, 511; W. G. Hoover, Phys. Rev. A, 1985, 31, 1695. 33 F. Müller-Plathe, J. Chem. Phys., 1997, 106, 6082. 34 Y. F. Zhukovskii, S. Piskunov, N. Pugno, B. Berzina, L. Trinkler and S. Bellucci, J. Phys. Chem. Solids, 2009, 70, 796. 35 S.-H. Lee and S.-H. Jhi, J. Phys.: Condens. Matter, 2015, 27, 075301. 36 Y. F. Gao, W. B. Bao, Q. Y. Meng, Y. H. Jing and X. X. Song, J. Non-Cryst. Solids, 2014, 387, 132. 37 X. Zhang, M. Hu and D. Tang, J. Appl. Phys., 2013, 113, 194307. 38 M. T. Dove, Introduction to Lattice Dynamics, Cambridge University Press, Cambridge, 1993. 39 J. M. Larkin, J. E. Turney, A. D. Massicotte, C. H. Amon and A. J. H. McGaughey, J. Comput. Theor. Nanosci., 2014, 11, 249. 40 J. A. Thomas, J. E. Turney, R. M. Iutzi, C. H. Amon and A. J. H. McGaughey, Phys. Rev. B: Condens. Matter, 2010, 81, 081411(R). 41 R. E. Peierls, Ann. Phys., 1929, 3, 1055. 42 S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472.
This journal is © The Royal Society of Chemistry 2015
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
COMMUNICATION
Cite this: DOI: 10.1039/c4nr07359b Received 12th December 2014, Accepted 16th March 2015 DOI: 10.1039/c4nr07359b
View Article Online View Journal
The unexpected non-monotonic inter-layer bonding dependence of the thermal conductivity of bilayered boron nitride Yufei Gao,†a Xiaoliang Zhang,†b Yuhang Jingc and Ming Hu*b,d
www.rsc.org/nanoscale
Hexagonal boron nitride (BN) and its bilayer form are very fascinating two-dimensional materials that have attracted tremendous interest recently. Their realistic applications in emerging nanoelectronics usually quest for manipulating the thermal transport properties in a precise manner. Using nonequilibrium molecular dynamics simulations, we herein studied the effect of inter-layer covalent bonding on the thermal conductivity of bilayered BN. We found that the in-plane thermal conductivity of bilayered BN, which can be largely tuned by introducing covalent bonding between the two BN layers, depends not only on the inter-layer bonding density, but also on the detailed topological configuration of the inter-layer bonds. For randomly distributed inter-layer bonding the thermal conductivity of bilayered BN decreases monotonically with inter-layer bonding density, the same behavior already found for bilayered graphene. However, for regularly arranged inter-layer bonding the thermal conductivity of bilayered BN surprisingly possesses a non-monotonic dependence on the inter-layer bonding density. This non-intuitive non-monotonic dependence is further explained by performing spectral energy density analysis, where the peak and valley values of the thermal conductivity are governed by different mechanisms. These results suggest the application of inter-layer covalent bonding in designing nanoscale devices with precisely tunable thermal conductivities.
Recently, there has been increased interest in the investigation on the physical and chemical properties of multilayered nanostructures, such as multilayered graphene and boron nitride (BN).1–3 A common feature of these materials is their layered
a School of Architecture & Civil Engineering, Shenyang University of Technology, Shenyang 110870, China b Institute of Mineral Engineering, Division of Materials Science and Engineering, Faculty of Georesources and Materials Engineering, Rheinisch-Westfaelische Technische Hochschule (RWTH Aachen University), 52064 Aachen, Germany. E-mail: [email protected] c Department of Astronautical Science and Mechanics, Harbin Institute of Technology, Harbin 150001, China d Aachen Institute for Advanced Study in Computational Engineering Science (AICES), RWTH Aachen University, 52062 Aachen, Germany † These authors contributed equally to this work.
This journal is © The Royal Society of Chemistry 2015
structures with high anisotropy, namely the atoms within the layers possess strong short-ranged covalent bonding and the interactions between different layers are described by weak long-ranged van der Waals forces.4 The van der Waals forces play an important role in making the layers maintain Bernal stacking, where the atoms in different layers do not eclipse each other. The previous research has demonstrated that the multilayered structures possess some fascinating properties compared to monolayered structures.5 Multilayered boron nitride is an analogue of graphite, in which the carbon atoms are replaced by boron and nitride atoms. The structure prefers to form a multilayered tubular structure or multilayered graphite-like flakes. Interestingly, the boron–nitride bonds in the in-plane possess a partially ionic character, which is quite different from that of graphene, whose bonds are completely covalent.6 Therefore, there are some distinct differences between these two structures. For example, the BN structures exhibit profound chemical and thermal stabilities.7,8 Moreover, multilayered BN structures are more promising for practical applications, such as optoelectronic technologies, tunnel devices, and field effect transistors.9–12 Nowadays, multilayered BN has been proposed as a potential material for thermal management applications,13–15 due to its high thermal conductivity and sensitivity to isotopic substitution. The multi-layered BN sheets and nanoribbons have even been foreseen as materials for designing microprocessors. Recently, the strong bonds between multi-walled carbon nanotubes and graphite layers can be created by irradiating the material with electron beam and femto-second laser excitation.16,17 Using similar ultraviolet laser irradiation, recently Komatsu et al.18 experimentally demonstrated that the structural relaxation of the sp2-bonded “metastable” phase of BN can be induced into more stable sp3-bonded phases at relatively low temperatures and below atmospheric pressure, which appears to be consistent with the recent pressure–temperature phase diagram of BN. Therefore, we speculate that the laser irradiation technique could also be used to create sp3bonds at the interface between bilayered BN. The inter-layer
Nanoscale
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Communication
bonds could cause strong cross-plane couplings between the layers and thus change the thermal transport properties of the materials significantly.19 Recently, the influence of inter-layer bonding in bilayered graphene on the thermal transport properties has been reported.19 The results showed that 5% bonding density monotonically leads to 75% decrease in the thermal conductivity, indicating that the thermal transport properties of bilayered graphene can be effectively controlled by modifying the bonding. Therefore, it is intuitive to anticipate that the emergence of the inter-layer bonds could also provide a controllable mechanism for tuning the physical properties of multilayered BN. In this paper, by performing nonequilibrium molecular dynamics simulation (NEMD) we studied the effect of interlayer covalent bonding on the thermal transport properties of bilayered BN. We demonstrate that the in-plane thermal conductivity of bilayered BN can be largely tuned by introducing covalent bonding between the two BN layers. However, the bonding density dependent thermal conductivity of bilayered BN exhibits unexpectedly different features, depending on the detailed topological configuration of the inter-layer bonds, which is counter-intuitive and has not been found in the literature so far. The underlying mechanism is provided by performing vibrational density of states and spectral energy density analysis. Our model system is armchair bilayered boron nitride sheets, consisting of m × n supercells, corresponding to a simulation box size of ∼2 × 540 nm2, with each supercell constructed by a bilayer honeycomb lattice and composed of 160 boron and nitride atoms. The structure of a supercell is shown in Fig. 1. The initial distance between the two BN layers is 0.33 nm. The inter-layer covalent bonding process is realized by the following steps: (1) the bonded atoms are chosen and moved closer to each other along the out-of-plane direction until the distance between the two bonded atoms reaches 0.145 nm (the equilibrium length of the B–N bond). (2) The period boundary condition is applied in the lateral (in-plane) directions. Initially, the bonded atoms are fixed and all other atoms in the system are relaxed with the NPT (constant temperature and constant pressure) ensemble for 100 000 steps. (3) The fix on the bonded atoms is removed and the entire system is relaxed with the NPT ensemble for additional 1 000 000 steps. Finally, the bilayered BN with inter-layer bonding is achieved. The interactions between boron and nitride atoms in the in-plane honeycomb structure and in the inter-layer bonding configuration are described by the Tersoff potential.20,21 The Tersoff potential with the optimized parameters has been proven suitable to investigate the properties of carbon, silicon, boron and nitride based nanomaterials.22–24 The interactions between boron and nitride atoms, that belong to different layers and are not covalently bonded, are presented with van der Waals forces. The standard LJ potential is used and the parameters can be found in ref. 25. It should be noted that the interfacial bonding created between BN layers should belong to the sp3-bonding configuration. In order to test whether the
Nanoscale
Nanoscale
Fig. 1 (a) Top and side views of bilayered BN without inter-layer covalent bonding. Color coding: light pink, boron; blue, nitride. From (b) to (d): a supercell of bilayered BN with inter-layer covalent bonding density 6.25%, 8.75%, 12.5% (left: randomly distributed bonding; right: regular bonding). The covalently bonded atoms are highlighted in green.
empirical Tersoff potential used is suitable for describing such sp3-bonding, we ran additional benchmark cases for bulk cubic BN (zinc blende structure) and wurtzite BN at room temperature, where boron and nitride atoms are explicitly bonded with sp3 bonding. The results of lattice constants and bulk modulus compared with experimental and ab initio calculations are shown in Table 1. It can be seen that the empirical
Table 1 Lattice constants and bulk modulus of cubic BN (zinc blende structure) and wurtzite BN calculated using the empirical Tersoff potential. The experimental and ab initio calculation results are also shown for comparison
Cubic BN Wurtzite BN
Tersoff potential Experiments or ab initio Tersoff potential Experiments or ab initio
Lattice constant, a (Å)
Lattice constant, c (Å)
Bulk modulus (GPa)
3.50 3.6226
— —
453.5 369,27 387.328
2.45 2.5529
4.11 4.2229
460.5 39730
This journal is © The Royal Society of Chemistry 2015
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Nanoscale
Communication
Tersoff potential, although originally optimized for planar sp2BN structures such as two-dimensional BN and BN nanotubes,20,21 can also simulate the bulk sp3-BN structures. This proves that the Tersoff potential can be effectively used for describing the interfacial sp3-like BN bonds in our model structures. All MD simulations herein are performed using the LAMMPS package31 with timestep 0.5 fs throughout. In our simulation, the velocity-Verlet algorithm is employed to integrate the equations of motion and the temperature is controlled by employing the Nosé–Hoover thermostat.32 First, we relaxed the entire system with the NPT ensemble for 50 ps, during which the pressure and temperature are kept at 1 atm and 300 K respectively. After that, we continued to relax the system with NVT (constant temperature and constant volume) and NVE (constant volume without a thermostat) ensembles, respectively. With the above processes, the system can reach the equilibrium state and a stable configuration can be obtained. Following equilibration, we used NEMD to calculate the thermal conductivity of bilayered BN. The constant heat flux was imposed by the Müller–Plathe method.33 During the NEMD simulation, the system is divided into slices along the longitudinal (heat flux) direction. The outmost two slices on the left and right side are set to be “cold” and “hot” regions, respectively. Then, the hottest atom in the “cold” region and the coldest atom in the “hot” region are selected and their energies are exchanged every 600–2000 steps, depending on the thermal conductivity of the system. The thermal conductivity of bilayered BN was calculated by Fourier’s law κ¼
JQ ; rT
ð1Þ
where JQ is the heat flux in the longitudinal (y) direction defined as JQ = (dQ/dt )/A, ∇T is the temperature gradient calculated by the linear fitting to the central part of the temperature profile of the system along the longitudinal direction, A is the cross-sectional area which is calculated as w·d, where w and d are the width and thickness of the bilayered BN, respectively. The thickness of monolayered BN is assumed to be 3.35 Å and the thickness of bilayered BN is assumed to be 6.7 Å. The total NEMD simulation time is 10 ns (2 × 107 steps). Normally, the temperature profile reaches steady state after 1–2 ns. Once the temperature profile becomes stable, we start to collect the heat flux and the temperature gradient to calculate the thermal conductivity every 100 ps. The reported thermal conductivity is averaged over all the data points collected during the steady state run. We have checked all the structures after the NPT and NVT relaxation runs and after imposing a constant heat flux during the NEMD simulation as well. We did not find a significant change in the structures, especially the interfacial covalent bonds created from the beginning remained until the end of the NEMD run. In addition, for selected cases we calculated the bond lengths for the sp3 hybridized (cross-plane) BN atoms and the regular sp2-bonded (in-plane) atoms to be about 1.491 Å and 1.444 Å, respectively. The in-plane BN bond
This journal is © The Royal Society of Chemistry 2015
length for bilayered BN is very close to that for monolayered BN calculated by the same classical potential (1.446 Å), both of which are in very good agreement with that obtained from ab initio calculations (1.446 Å).34,35 The trend that the bond length of sp3 bonding is slightly larger than that for sp2 bonding is also consistent with that found from ab initio calculation.28,35 Again, the results fully support that the optimized Tersoff potential is suitable to describe the pseudo-sp3 hybridized bonds and the honeycomb two-dimensional sp2 bonds as well. The relationship between the thermal conductivity of bilayered BN and the inter-layer bonding density is shown in Fig. 2. Here we consider two different configurations of interlayer bonding: regular and random distribution, with typical model structures shown in Fig. 1. The random model is the one used in the previous investigation on bilayered graphene.19 For comparison we also show the thermal conductivity of monolayered BN and bilayered BN without inter-layer covalent bonding, i.e. the inter-layer interaction is pure van der Waals forces. We first notice that generally there is a large decrease in the thermal conductivity of bilayered BN with the increase of bonding density, which is consistent with a previous study on bilayered graphene,19 where there is 75% decrease in thermal conductivity with 5% bonding density. The reduction is generally understood as the consequence of remarkable phonon coupling (scattering) between the two layers when strong covalent bonds are introduced at the interface. A more distinct feature in Fig. 2 is that in contrast to the monotonic dependence of the thermal conductivity of bilayered BN with randomly distributed interfacial bonding, which was also found for bilayered graphene,19 the thermal conductivity of bilayered BN with regular interfacial bonding exhibits a non-monotonic dependence on the bonding density.
Fig. 2 Dependence of the thermal conductivity of bilayered BN on the inter-layer covalent bonding density. The horizontal solid line represents the thermal conductivity of monolayered BN. The zero bonding density denotes bilayered BN with only van der Waals interactions between the layers.
Nanoscale
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Communication
Nanoscale
In general the thermal conductivity in the case of the regular bonding configuration (see the typical structure in the right panel of Fig. 1) is higher than that for random bonding, except for the two cases of the bonding density of 2.5% and 8.75%. For regular bonding the thermal conductivity has two peak values at the bonding density of 6.25% and 12.5%, where the thermal conductivity of the bilayered BN between regular and random bonding can differ 2–3 fold. These results indicate that when introducing covalent bonds at interfaces, the thermal conductivity of the bilayered BN not only depends on the overall bonding density, but also is very sensitive to the detailed topological arrangement of the interfacial bonding. Throughout our simulation we are aware that the different bonding configurations might lead to different thermal conductivities, especially for randomly distributed bonds. Therefore, in our paper we ran five different bonding configurations for each random bonding case in Fig. 1(b)–(d). All thermal conductivity results reported in Fig. 2 are the averaged value of these five deferent runs. In addition, we ran additional simulations with the bonding density of 6.25% (500 nm) and the five covalent bonds taking up two lines (in this case three covalent bonds form one line and the two other bonds form the other line). The thermal conductivity of this new configuration is 139.9 W mK−1. Compared with the results in Fig. 2, we can find that this value is very close to the average thermal conductivity of the bilayered BN with random interlayer bonding (130.9 W mK−1) but is far lower than that for regular interlayer bonding (275.7 W mK−1). Therefore, it should be noted that changing the saturated bond line in Fig. 1(b) into two unsaturated lines results in a significant drop in thermal conductivity, as the unsaturated bond lines result in lower phonon group velocities and also always have a very strong interaction between the lines and thus the phonon scattering is enhanced, the same mechanism as that for Fig. 1(d), as we will see later. We also considered the symmetric distribution of the interlayer bonding. We run additional cases of 4, 6, and 8 bonds that take up symmetrical positions but form two unsaturated lines, corresponding to the bonding density of 5%, 7.5%, and 10%, respectively. The thermal conductivity of bilayered BN with regular, symmetric, and random inter-layer bonding is compared in Table 2. From Table 2 we can see that for the bonding density of 5% and 7.5% the thermal conductivity of bilayered BN with sym-
Table 2 Comparison of the thermal conductivity of bilayered BN between different inter-layer bonding distributions: regular, symmetric, and random. The thermal conductivity is in the units of W mK−1
metric inter-layer bonding is lower than that for regular bonding. This trend is opposite for the case with the bonding density of 10%, which can be attributed to the less interference between the two bonding lines in the case of symmetric distribution. Nevertheless, the results clearly show that for all cases the thermal conductivity of bilayered BN with symmetric inter-layer bonding is higher than the corresponding randomly distributed cases. Since the symmetric distribution can also be regarded as regular distribution, once again the results in Table 2 confirm the trend in Fig. 2 that the thermal conductivity of bilayered BN strongly depends on the detailed topological configuration of the inter-layer bonds and the regular inter-layer bonding always yields higher thermal conductivity. The symmetric cases also confirm the nonmonotonic dependence in Fig. 2, when taking the 12.5% regular case shown in Fig. 1 as the symmetric distribution. Moreover, combining Fig. 1 and 2 we can see that the following three cases of inter-layer bonding arrangement can lead to a significant increase in thermal conductivity: (i) for low bonding density (2.5%–5%) the inter-layer bonds arrange only in one line but certainly cannot form one full line, i.e. only one unsaturated line can form at the interface; (ii) for higher bonding density (6.25%–12.5%) the inter-layer bonds form in one full (saturated) line along the heat flux direction and the rest try to fill up the other line; (iii) the inter-layer bonds form two symmetric but unsaturated lines. Other inter-layer bonding arrangements that do not fulfill any requirement of the above three cases has much lower thermal conductivity. It is well known that the thermal conductivity of many nanosystems calculated by the NEMD method exhibits size dependence. To address this issue, for selected cases with the bonding density of 8.75% and 12.5%, corresponding to the typical minima and maxima in Fig. 2, we ran additional simulations with the length doubled, i.e. 1078 nm. For each bonding density, we generated three different random bonding configurations and one regular bonding configuration, so the finally reported value for the random bonding configuration is the average of these three independent runs. The comparison of the calculated thermal conductivity between systems with lengths of 540 nm and 1078 nm is given in Table 3. We found that the thermal conductivity of bilayered BN only has a slight incremental change when the length is doubled. More importantly, the trend is exactly the same, i.e.
Table 3 Comparison of the thermal conductivity of bilayered BN with inter-layer bonding calculated by the NEMD method between systems with lengths of 540 nm and 1078 nm. The calculated thermal conductivity is in the units of W mK−1
Length (nm) 540
Distribution Bonding density
Regular
Symmetric
Random
5.0% 7.5% 10%
201.9 168.0 131.5
161.3 149.1 155.2
135.9 127.2 98.4
Nanoscale
1078
Bonding density
Random bonding
Regular bonding
Random bonding
Regular bonding
8.75% 12.5%
94.8 93.6
112.9 222.2
118.9 108.9
156.0 220.5
This journal is © The Royal Society of Chemistry 2015
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Nanoscale
for the bonding density of 8.75% the thermal conductivity is still insensitive to the detailed bonding configuration, while for the bonding density of 12.5% the thermal conductivity for the regular bonding configuration is at least twice as large as that for the random bonding configuration. Moreover, our spectral energy density (SED) analysis of the above two bonding density cases yields a similar trend to that found in NEMD simulation (see the results later). Since SED is based on equilibrium molecular dynamics (EMD) simulation, the size effect in SED analysis is minor. Therefore, we believe that the nonmonotonic dependence of the thermal conductivity on the bonding density revealed by NEMD simulation is size independent. We also run more NEMD simulations for the selected cases of bilayered BN with a fully random bonding density of 3.75%, 6.25%, and 8.75%. We created the fully randomly distributed interfacial covalent bonds by randomly selecting bonded atoms in the entire simulation domain, i.e. ∼2 × 540 nm2. The results show that, compared to the original construction method by duplicating the supercell containing irregularly aligned bonded atoms, the thermal conductivity of the bilayered BN with fully interfacial random bonding is only slightly lower. Nevertheless, the trend remains unchanged. Once again, it proves that the thermal transport of bilayered BN with random interfacial covalent bonding is not sensitive to the detailed bonding arrangement. The intuitive phenomenon of the reduction in the thermal conductivity of bonded bilayered BN in contrast to free bilayered BN can be understood in terms of enhanced phonon coupling (scattering) between the two bonded layers, as evidenced by the phonon density of state (PDOS) of BN atoms in our model structure. The PDOS was calculated by imposing a Fourier transform on the atomic velocity autocorrelation function.36 The PDOS values between the bonded and non-bonded atoms are compared in Fig. 3. Both PDOS values are similar for frequencies below 30 THz. However, there are distinct differences in the frequency between 30 and 50 THz. The high
Communication
frequency peak in the PDOS of the bonded atoms is approximately at 37 THz and for non-bonded atoms the peak is at 49 THz. That is to say, there is a significant redshift for the PDOS of the bonded atoms, which means the probability of phonon scattering increases, provided that the total number of phonon modes is exactly the same for both regular and random bonding structures. Therefore, the redshift of the high frequency phonons of the bonded atoms leads to a significant reduction of the thermal conductivity of the bonded structure. We also examined the effect of cross-linking on the contribution of the in-plane and out-of-plane phonons in bilayered BN. To this end, we quantified the relative contributions of the lattice vibrations in the x, y, and z directions to the total heat flux of the bilayered BN, corresponding to the in-plane longitudinal modes, in-plane transverse modes, and the out-ofplane flexural modes, respectively, by using the NEMD-based method we proposed recently37 J A!B;α ¼
1 XX F ijα ðυiα þ υjα Þ; 2A i[A j[B
ð2Þ
where α is the x, y, or z direction, A and B are the left and right side of a virtual interface located at the middle of the bilayered BN in the x (heat flux) direction, JA→B,α is the heat flux from A to B due to lattice vibrations in the α direction, A is the crosssectional area, Fijα is the component of the force acting on atom i due to atom j, and υiα is the α component of the velocity of atom i. Details of calculation in the form of the three-body Tersoff potential can be found in ref. 37. We performed the polarization analysis of three typical cases of bonding densities: 6.25%, 8.75%, and 12.5%. The relative contribution ( percentage) to the total heat flux from vibrations in x (longitudinal), y (transverse), and z (out-of-plane) directions for these three cases is compared in Table 4. It is clearly seen that for the cases of random interfacial bonding the relative contribution from in-plane and out-ofplane vibrations is nearly the same. For the case of the bonding density of 8.75%, the relative contribution percentage is also the same as those for random bonding. This is consistent with the results of their similar thermal conductivity. Comparing the results between regular and random bonding, we find that when the interfacial covalent bonds form one or two saturated lines, the interfacial bonds affect the in-plane phonon transport significantly, i.e. the contribution from the in-plane longitudinal modes is reduced while the contribution
Table 4 Comparison of the relative contribution ( percentage) to the total heat flux from vibrations in x (longitudinal), y (transverse), and z (out-of-plane) directions between different bonding densities
Fig. 3 Comparison of the phonon density of states between bonded atoms and non-bonded atoms in the structures of bilayered BN.
This journal is © The Royal Society of Chemistry 2015
6.25% 8.75% 12.5%
Regular interfacial bonding
Random interfacial bonding
x
y
z
x
y
z
45.4% 59.9% 44.4%
26.8% 27.7% 29.4%
27.8% 12.7% 26.2%
56.2% 56.5% 57.2%
28.6% 29.5% 28.7%
15.2% 14.0% 14.1%
Nanoscale
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Communication
Nanoscale
from the out-of-plane flexural modes is enhanced (highlighted in Table 4) and the contribution from the in-plane transverse modes remains the same. This is understandable considering that when the crosslinking covalent bonds are formed regularly at the interface, the atoms in different layers can vibrate in unison in the out-of-plane direction, and in this way the out-of-plane flexural modes have less phonon scattering, which results in higher thermal conductivity, since usually the flexural modes in graphene-like two-dimensional materials conduct heat more efficiently. We now explain the non-intuitive phenomenon that the regular bonding configuration results in substantially higher thermal conductivity than the random bonding configuration with exactly the same bonding density, e.g. for the bonding density of 12.5% in Fig. 2. To this end, we calculated the frequency dependent phonon lifetime by performing spectral energy density analysis. The phonon normal modes are obtained by38 rffiffiffiffiffiffi X mj ˙~ ~ e*j ð~ k; νÞexpð2πi~ k~ r l Þ: ð3Þ υjl ðtÞ~ Qð k; ν; tÞ ¼ N jl Then the spectral energy density is calculated by39 2 ð ˙~ Φð~ k; ν; f Þ ¼ Qð k; ν; tÞexpð2πiftÞdt ;
ð4Þ
and the phonon lifetime (relaxation time) is obtained by fitting the Lorentzian function.40 Finally, we investigated the mode dependent thermal conductivity of the bilayered BN, which can be derived from the phonon Boltzmann transport equation (BTE) under the relaxation time approximation41 XX κα ¼ cph υg;α 2 ð~ k; νÞτð~ k; νÞ; ð5Þ ~ k
ν
where ~ k is the wave vector in the first Brillouin zone, ν is the phonon branch, κα denotes the thermal conductivity in the αdirection, cph is the phonon volumetric specific heat calculated kB using cph ¼ , kB is the Boltzmann constant, V is the system V volume (a thickness of 3.35 Å and 6.7 Å is assumed for monolayered and bilayered BN, respectively), υg,α is the α-component of the phonon group velocity, and τ is the phonon lifetime. The frequency dependent phonon lifetimes of bilayered BN (inter-layer bonding density of 12.5%) with random and regular bonding configurations are compared in Fig. 4(a). We can see that the frequency dependent phonon lifetime of bilayered BN for both regular and random bonding roughly follows the power law of ω−α (α ∼ 0.71). The phonon lifetime in bilayered BN with the regular bonding configuration is only slightly higher than that with the random bonding configuration, which obviously cannot explain the three times higher thermal conductivity that occurs for the regular bonding configuration (Fig. 2). We also noticed that there is a significant increase in the lifetime of low frequency phonons (below 1 THz), but still this is not the main reason for the large thermal conductivity enhancement. Then, we further calcu-
Nanoscale
Fig. 4 Comparison of (a) phonon lifetime, (b) phonon group velocity, and (c) frequency dependent accumulated thermal conductivity of bonded bilayered BN between random and regular bonding configurations (inter-layer bonding density 12.5%). The dashed line is the power law fitting (ω−α) to the frequency dependent phonon lifetime.
lated the square of the phonon group velocity (υg,α2) in eqn (5) as a function of frequency, as shown in Fig. 4(b). Now it is clearly seen that the phonon group velocity for the case of regular inter-layer bonding is much higher than that for random bonding. According to eqn (5), finally we can plot the accumulated thermal conductivity of the bonded bilayered BN as a function of frequency, as shown in Fig. 4(c). Initially the accumulated thermal conductivity is almost the same for both random and regular bonding, when the phonon frequency is less than 2 THz, which means that the low frequency phonons are not sensitive to the detailed interfacial bonding configuration. This is understandable considering that the low frequency phonons usually have a considerably long wavelength, which is much larger than the characteristic distance between the bonded sites, thus they usually travel through the struc-
This journal is © The Royal Society of Chemistry 2015
View Article Online
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
Nanoscale
tures without scattering. Above 2 THz, the thermal conductivity for the case of regular bonding increases dramatically until the phonon frequency reaches around 30 THz, while there is only a small increment in the thermal conductivity of bilayered BN with random inter-layer bonding. For even higher frequency phonon modes, both regular and random bonding show a minor increase in the overall thermal conductivity, suggesting that the very high frequency phonon modes do not contribute much to the thermal transport in bilayered BN. This is reasonable considering that the high frequency phonon modes are usually localized and propagate diffusively with very low phonon relaxation time [Fig. 4(a)]. Combining Fig. 4(a)–(c) we can conclude that the much higher thermal conductivity of the bilayered BN with regular bonding as compared with random bonding is mainly due to the large enhancement in the group velocity of the phonon modes in the medium frequency range (2 to 30 THz). We also ran SED simulations for the bonding density of 8.75% with regular and random bond patterns. The results are not reported here for brevity. First of all, the thermal conductivity of bilayered BN with random and regular bonding calculated by the SED method is 75.7 and 92.2 W mK−1, respectively, which are very close to the NEMD simulation results (94.8 and 112.9 W mK−1, respectively), as shown in Fig. 2. More importantly, the thermal conductivity trend predicted by SED calculation is in accord with that found from NEMD simulations, i.e. for the bonding density of 8.75%, the thermal conductivity of bilayered BN is insensitive to the bond pattern. This is because, as revealed by the SED calculations, for the bonding density of 8.75% there is no significant difference in both the phonon mode specific group velocity and lifetime between regular and random patterned interfacial bonds. Before closing, we explain the phenomenon that the thermal conductivity of bilayered BN with regular inter-layer bonding sometimes reaches a minimum at certain inter-layer bonding densities, i.e. the non-monotonic inter-layer bonding dependence is shown in Fig. 2. Taking bonding density of 8.75% as an example, we found that the thermal conductivity is only slightly higher than that for randomly distributed bonding and is significantly lower than that for the bonding density of 6.25% and 12.5%. By comparing the SED results of the phonon group velocity and phonon lifetime between the regular bonding density of 8.75% and 12.5%, we found that the significantly higher thermal conductivity for the regular bonding density of 12.5% than that for regular 8.75% stems from two factors: (1) the first factor, also the dominant factor, is that the phonon group velocities are largely enhanced, indicating that the formation of the full lines of the interfacial covalent bonds results in higher phonon group velocities. (2) “Interference” between bond lines: for the bonding density of 8.75% the inter-layer bonds only form one full line (we call it saturated line) along the heat flux direction and the rest do not fill up the other line (unsaturated line), as shown by the atomic structure in Fig. 1(c). In contrast, for the cases of 6.25% and 12.5% (regular bonding) all the inter-layer bonds either line up in one saturated line [6.25%, Fig. 1(b)] or form two
This journal is © The Royal Society of Chemistry 2015
Communication
saturated lines [12.5%, Fig. 1(d)]. The different bonding configurations between the saturated line and the unsaturated line for the bonding density of 8.75% will induce strong phonon scattering between the bonded atoms (reduced phonon lifetime) and therefore the heat conduction is significantly hindered. These results indicate that when controlling thermal transport of bilayered BN by inter-layer covalent bonds, one should not only consider the interfacial bonding density, but also carefully arrange the topological configuration of the inter-layer bonds. Finally, it is worth pointing out that the nonmonotonic behavior of the thermal conductivity on the density of regularly distributed interlayer bonds is not a particular property of boron nitride. We have also started the same calculation procedure as presented in this paper for bilayered graphene using the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential,42 considering the similarity of their honeycomb structures. Our preliminary results for two typical interlayer bonding densities of 6.25% and 12.5% show the same trend as that found for bilayered BN. Currently, this work is ongoing and the results will be submitted for publication separately. In summary, the effect of inter-layer covalent bonding on the thermal transport properties of bilayered BN is investigated by performing nonequilibrium molecular dynamics simulations. The in-plane thermal conductivity of bilayered BN can be largely tuned by introducing covalent bonding between the two BN layers. For randomly distributed inter-layer bonding the thermal conductivity of bilayered BN is reduced monotonically with the inter-layer bonding density. However, the thermal conductivity unexpectedly shows a non-monotonic dependence for regularly arranged inter-layer bonding. The mechanism underlying the counter-intuitive non-monotonic dependence is identified by performing spectral energy density analysis. We demonstrate that the peak values of the thermal conductivity of bonded bilayered BN with regular bonding are mainly due to the large enhancement in the group velocity of the phonon modes in the medium frequency range and the phonon lifetime augmentation only plays a minor role, while the valley values result from the low phonon group velocities and strong phonon scattering between the saturated and unsaturated bond lines. In short, when introducing chemical covalent bonding between the layers, the thermal conductivity of the bilayered BN not only depends on the inter-layer bonding density, but also is very sensitive to the topological configuration of how the inter-layer bonds are arranged. Our findings provide a guide to precisely modulate the thermal transport properties of bilayered BN with chemical engineering and they may be of use in tuning the electronic and optical properties for electronic, thermoelectric, photovoltaic, and opto-electronic applications.
Acknowledgements This research work was supported by the National Natural Science Foundation of China under grant no. 11272099 and
Nanoscale
View Article Online
Communication
11304059. Simulations were performed with computing resources granted by the Jülich Aachen Research Alliance-High Performance Computing (JARA-HPC) from RWTH Aachen University under project no. jara0104.
Published on 26 March 2015. Downloaded by University of Sussex on 26/03/2015 18:10:03.
References 1 S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau and A. A. Balandin, Nat. Mater., 2010, 9, 555. 2 W. Q. Han, L. J. Wu, Y. M. Zhu, K. Watanabe and T. Taniguchi, Appl. Phys. Lett., 2008, 93, 223103. 3 D. Pacile, J. C. Meyer, C. O. Girit and A. Zettl, Appl. Phys. Lett., 2008, 92, 133107. 4 K. Gordiz, S. M. V. Allaei and F. Kowsary, Appl. Phys. Lett., 2011, 99, 251901. 5 D. Golberg, Y. Bando, Y. Huang, T. Terao, M. Mitome, C. C. Tang and C. Y. Zhi, ACS Nano, 2010, 4, 2979. 6 H. M. Ghassemi and R. S. Yassar, Appl. Mech. Rev., 2010, 63, 020804. 7 C. Y. Zhi, Y. Bando, C. C. Tang and D. Golberg, Mater. Sci. Eng., R, 2010, 70, 92. 8 J. S. Wang, C. H. Lee and Y. K. Yap, Nanoscale, 2010, 2, 2028. 9 K. Watanabe, T. Taniguchi and H. Kanda, Nat. Mater., 2004, 3, 404. 10 L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A. S. Mayorov, N. M. R. Peres, A. H. Castro Neto, J. Leist, A. K. Geim, L. A. Ponomarenko and K. S. Novoselov, Nano Lett., 2012, 12, 1707. 11 G. H. Lee, Y. J. Yu, C. Lee, C. Dean, K. L. Shepard, P. Kim and J. Hone, Appl. Phys. Lett., 2011, 99, 243114. 12 J. Beheshtian, A. Sadeghi, M. Neek-Amal, K. H. Michel and F. M. Peeters, Phys. Rev. B: Condens. Matter, 2012, 86, 195433. 13 T. Ouyang, Y. Chen, Y. Xie, K. Yang, Z. Bao and J. Zhong, Nanotechnology, 2010, 21, 245701. 14 D. A. Stewart, I. Savic and N. Mingo, Nano Lett., 2009, 9, 81. 15 C. W. Chang, A. M. Fennimore, A. Afanasiev, D. Okawa, T. Ikuno, H. Garcia, D. Li, A. Majumdar and A. Zettl, Phys. Rev. Lett., 2006, 97, 085901. 16 B. Peng, M. Locascio, P. Zapol, S. Y. Li, S. L. Mielke, G. C. Schatz and H. D. Espinosa, Nat. Nanotechnol., 2008, 3, 626. 17 J. Kanasaki, E. Inami, K. Tanimura, H. Ohnishi and K. Nasu, Phys. Rev. Lett., 2009, 102, 087402. 18 S. Komatsu, K. Kobayashi, Y. Sato, D. Hirano, T. Nakamura, T. Nagata, T. Chikyo, T. Watanabe, T. Takizawa, K. Nakamura and T. Hashimoto, J. Phys. Chem. C, 2010, 114, 13176.
Nanoscale
Nanoscale
19 A. Rajabpour and S. M. Vaez-Allaei, Appl. Phys. Lett., 2012, 101, 053115. 20 C. Sevik, A. Kinaci, J. B. Haskins and T. Cagin, Phys. Rev. B: Condens. Matter, 2011, 84, 085409. 21 A. Kinaci, J. B. Haskins, C. Sevik and T. Cagin, Phys. Rev. B: Condens. Matter, 2012, 86, 115410. 22 Y. F. Gao, Y. H. Jing, Q. Y. Meng, L. Zhang, J. Q. Liu and X. Qin, Phys. Status Solidi B, 2012, 249, 1728. 23 Y. H. Jing, M. Hu and L. C. Guo, J. Appl. Phys., 2013, 114, 153518. 24 M. Hu, X. L. Zhang and D. Poulikakos, Phys. Rev. B: Condens. Matter, 2013, 87, 195417. 25 B. Fakrach, A. Rahmani, H. Chadli, K. Sbai, M. Benhamou, M. Bentaleb, J.-L. Bantignie and J.-L. Sauvajol, J. Phys.: Condens. Matter, 2012, 24, 335304. 26 G. L. Doll, J. A. Sell, C. A. Taylor and R. Clarke, Phys. Rev. B: Condens. Matter, 1991, 43, 6816. 27 E. Knittle, R. M. Wentzcovitch, R. Jeanloz and M. L. Cohen, Nature, 1989, 337, 349. 28 B. Wen, J. Zhao, R. Melnik and Y. Tian, Phys. Chem. Chem. Phys., 2011, 13, 14565. 29 A. Nagakubo, H. Ogi, H. Sumiya, K. Kusakabe and M. Hirao, Appl. Phys. Lett., 2013, 102, 241909. 30 K. Kim, W. R. L. Lambrecht and B. Segall, Phys. Rev. B: Condens. Matter, 1996, 53, 16310. 31 S. Plimpton, J. Comput. Phys., 1995, 117, 1. 32 S. Nosé, J. Chem. Phys., 1984, 81, 511; W. G. Hoover, Phys. Rev. A, 1985, 31, 1695. 33 F. Müller-Plathe, J. Chem. Phys., 1997, 106, 6082. 34 Y. F. Zhukovskii, S. Piskunov, N. Pugno, B. Berzina, L. Trinkler and S. Bellucci, J. Phys. Chem. Solids, 2009, 70, 796. 35 S.-H. Lee and S.-H. Jhi, J. Phys.: Condens. Matter, 2015, 27, 075301. 36 Y. F. Gao, W. B. Bao, Q. Y. Meng, Y. H. Jing and X. X. Song, J. Non-Cryst. Solids, 2014, 387, 132. 37 X. Zhang, M. Hu and D. Tang, J. Appl. Phys., 2013, 113, 194307. 38 M. T. Dove, Introduction to Lattice Dynamics, Cambridge University Press, Cambridge, 1993. 39 J. M. Larkin, J. E. Turney, A. D. Massicotte, C. H. Amon and A. J. H. McGaughey, J. Comput. Theor. Nanosci., 2014, 11, 249. 40 J. A. Thomas, J. E. Turney, R. M. Iutzi, C. H. Amon and A. J. H. McGaughey, Phys. Rev. B: Condens. Matter, 2010, 81, 081411(R). 41 R. E. Peierls, Ann. Phys., 1929, 3, 1055. 42 S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472.
This journal is © The Royal Society of Chemistry 2015
