Thermal and transport properties of pristine single-layer hexagonal boron nitride: a first principles investigation

Molecular dynamics is used in combination with density functional theory to determine the thermal transport properties of the single-layer hexagonal boron nitride (SL h-BN) from ab initio calculations. Within this approach, the possible anisotropy in the thermal conductivity of SL h-BN was studied. For samples with finite length (of the order of 20 nm), we find a significant dependence of the conductivity on the transport direction. We make a direct comparison of the results obtained for 2D layers and for nanoribbons with similar size, and show that, as a consequence of edge scattering, the ribbon geometry induces a significant decrease in the conductivity, and produces a strong change in the anisotropy. For the zigzag and armchair transport directions, the dependence of the thermal conductivity on the system length was also obtained, as well as its value in the 2D-bulk limit case. A very small anisotropy was found for the limit of long samples, in contrast with the finite length ones. This is explained analyzing the dependence of the average square group velocities on the transport direction and the phonon frequency.


Introduction
The novel properties of carbon-based honeycomb structures [1,2,3,4,5] have stimulated strong interest in isomorphic two-dimensional (2D) materials such as hexagonal boron nitride (h-BN). The strong covalent bond between B and N atoms by sp2 hybridization [6] and its honeycomb atomic structure result in some of its physical properties being similar to those of graphene, like e.g.: strong mechanical properties, high thermal stability and superior thermal conductivity [7,8,9,10,11]. The fact that h-BN is an insulator, makes it an almost ideal insulating and dielectric layer for graphene based electronics. As a matter of fact, the feasibility of graphene/h-BN devices with improved electrical properties in comparison to the graphene on amorphous SiO 2 substrate counterpart has been demonstrated [12]. Moreover, the high thermal conductivity of h-BN makes it perfect for the thermal management in nanodevices where efficient heat dissipation is a key factor [13].
Understanding the lattice thermal transport properties of single layer (SL) h-BN is crucial for the design of novel nanodevices as well as for improving our fundamental insight into the behavior of phonon transport in 2D layered structures [14,15,16]. Recently, thermal conductivities in few-layer h-BN samples were measured obtaining: 360 W/mK for the 11-layer h-BN sample [17], 227-280 W/mK for the 9-layer sample [10], 250 W/mK for 5-layers [17] and 484 W/mK for 2-layers [18]. In contrast, the highest recorded thermal conductivity of bulk h-BN is around 400 W/mK [19] at room temperature. Despite sharing similar lattice structure with graphene, the obtained thermal conductivities of h-BN are significantly lower than for their carbon counterparts.
From the theoretical point of view, the thermal conductivity of h-BN has been extensively studied, but most of the studies were focused on nanoribbon structures. Thermal conductivity in two-dimensional layers was computed using model potentials from the Einstein relations, obtaining 460 W/mK [20], Green Kubo formalism 400 W/mK [21] and solving the linearized Boltzmann transport equation (BTE) 780-810 W/mK [22], where the Tersoff potential was adopted. Other approaches such as thigh-binding models have also been used in combination with molecular dynamics to study polycrystalline h-BN for different grain sizes reporting values close to 600 W/mK for the infinite limit [23]. More recently, Capellotti et al. have solved the linearized BTE equation using parameters obtained from Density Functional Perturbation Theory (DFPT) for SL h-BN reporting 1050 W/mK as the value of the bulk thermal conductivity [24]. This result is the most accurate value from first principles calculations to date since it is based on the exact solution of the BTE equation where phonon frequencies, phonon lifetimes, and scattering rates have been computed fully ab initio by DFPT. Concerning the thermal properties in nanorribons, the values reported so far are significantly lower than the 2D infinite monolayer counterpart. Moreover, strong dependence as a function of the width and the length of the ribbon has also been observed [25,26]. Although it was reported that the thermal conductivities for zigzag and armchair oriented nanoribbons are different [27,21,25] suggesting a possible anisotropy, the only work that has studied systematically the thermal conductivity as a function of the orientation of the ribbon was done by Y. -C. Chen et al. [26] by Non Equilibrium Molecular Dynamics (NEMD) using Tersoff-type potentials. They have reported anisotropy due to the boundary scattering in the free edges of the ribbon in a similar case to the graphene nanostructure [28].
In this work, the Approach to Equilibrium Molecular Dynamics (AEMD) [29,30,31] method was used in combination to density functional theory (DFT) as implemented in Siesta [32] to elucidate from first-principles the possible anisotropy in SL h-BN. The AEMD approach was used since it has demonstrated that its computational cost is reduced compared to other MD techniques [29,30] and, therefore, it optimally matches the high computational demands posed by DFT calculations. Siesta shows a comparatively high numerical efficiency due to the description of the electronic wavefunctions by means of strictly localized basis sets. By combining AEMD with Siesta we have reached an unprecedented computational efficiency which the possibility to determine accurately the thermal conductivity in h-BN from ab-initio MD.

Theory and computational setup
The evaluation of the thermal conductivity was done using the Approach to Equilibrium Molecular Dynamics (AEMD) methodology as described in [29,30]. This approach evaluates the thermal conductivity from the temperature transient regime using the exact solution of the heat transport equation for an initial step-like temperature profile. The system is initially decomposed in two different regions (subsystems 1 and 2) characterized by average temperatures T 1 and T 2 , respectively. The average temperature at given time t is defined as T = 1 L L 0 T (z, t)dz being L the subsystem length. From the solution of the heat equation for the above initial condition, the time evolution of the average temperature difference is expected to have the form being α n = 2πn/L z , L z the length of the total system and n an integer number. From the thermal diffusivityκ, the thermal conductivity κ can be obtained as where V is the total system volume and C v is the heat capacity. The quasi harmonic approximation [33] was used to evaluate the heat capacity C v needed for the calculation of the thermal conductivity (see eq. 3). C v at given temperature T can be directly calculated from the phonon density of states (DOS) g(ν) using the following expression [34] with the weighing factor W (x) = x 2 e x /(e x − 1) 2 . h is the Planck constant, k B is the Boltzmann constant, ν is the frequency of the phonons and r and N are the number of degrees of freedom and the number of atoms in the unit cell, respectively. From the heat capacity, the Debye temperature Θ D can be also obtained by fitting C QHA v (T ) to the expression for the heat capacity under the Debye approximation for 2D materials [35] This magnitude is a measurement of the temperature above which basically all modes are excited and below which some modes are instead frozen out. All the simulations were performed using the ab-initio Siesta code [32,36] taking advantage of its molecular dynamics package in combination to DFT for the forces determination (see Supporting Information for the code implementation [37]). Within this combination, the time dependent positions and velocities of the atoms are described by fully ab-initio interatomic forces in contrast to force fields methods based on empirical potential approaches.
The DFT-Siesta calculations were performed using norm-conserving Troullier-Martins [38] pseudopotentials with nonlinear core corrections within the local density approximation (LDA) with a Ceperley-Alder [39] exchange-correlation potential as parametrized by Perdew-Zunger [40]. A non-optimized Siesta single-ζ basis set [41] was used to describe the valence electrons, and a mesh cutoff of 800 Ry for the realspace grid [32]. The h-BN unit cell and the atomic positions were relaxed using conjugate gradients until the forces acting on each atom were smaller than 0.04 eV/Å. The first Brillouin zone (FBZ) was sampled by a 8×8×1 k-points mesh. The convergence criteria for the electron density matrix was chosen as 10 −4 between two consecutive steps. In order to properly describe a single layer h-BN, a large spacing of 10Å of vacuum was added in order to prevent interaction between periodic images. The optimized lattice constant and the B-N bond distance were a=2.511Å and d=1.450Å, respectively, in good agreement to the values reported on [42,43,44,45]. Phonon calculations were done using finite differences: we obtain the dynamical matrix from the atomic forces in a large enough supercell, caused by atomic displacements of the atoms in the unit cell of 0.04 Bohr from their equilibrium positions.
To perform the AEMD simulations, the simulation box was divided in two equal subsystems containing the same number of atoms and dimensions. Starting from a random Maxwell-Boltzmann velocity distribution, the initial step-like temperature profile was created during the equilibration run. At this stage, the subsystems were coupled to external Berendsen thermostats which rescale appropriately the atomic velocities of each subsystem in order to reach the target temperatures T 1 = 400 K and T 2 = 200 K, respectively. The average temperatures were obtained from the kinetic energy of the atoms in each subsystem. The number of equilibration steps was 500 and the equations of motion were integrated by the Velocity Verlet algorithm using 2 fs time step. Previous tests provided evidence that this number of steps is enough to reach the target temperatures and the calculated thermal conductivity is independent of the number of steps in this equilibration stage. Once the initial temperature profile was created with ∆T (0) = 200 K, the system was aged by a microcanonical run for which the time dependence of the temperature difference (eq. 1) was monitored. Each simulation was repeated 5 times starting from different initial distributions of the atomic velocities and, therefore, all calculated thermal conductivity values have been obtained as a configurational average. Finally, the number of fitting exponentials was set to n = 20 in eq. 1 and the thickness in the perpendicular direction is chosen as the experimental one of the corresponding three-dimensional (3D) material 3.306Å as is usually done in these simulations [24,26]. In order to assess any anisotropy in thermal transport, our simulation boxes are build from a 2D single layer, by defining rectangular samples of material with a given orientation angle σ (measured with reference to the zigzag direction), and with lengths L z and L y in the directions parallel and perpendicular to transport, respectively, as indicated in Fig. 1. The number of atoms in the direction of transport was made as large as possible, while the length in the perpendicular direction was reduced to a minimum. Periodic boundary conditions are used in the MD runs in which we simulate the 2D single layer. 2D samples were used to determine the thermal conductivity along the zigzag and armchair directions as a function of the system length. The width of the zigzag samples were L y = 0.869 nm while 10 nm ≤ L z ≤ 100 nm. The corresponding number of atoms in the systems varied from 320 to 3200. Similar numbers were obtained for the armchair case but the width of the samples was L y = 1.004 nm. Eight different orientation were studied from the zigzag (ZZ) configuration (σ = 0 • ), to the armchair (AC) configuration (σ = 30 • ), with system sizes around L z ≈ 20 nm and L y ≈ 0.86 nm (see Supporting Information section for further system details [37]).
Concerning the simulation of finite-width structures, the zigzag and the armchair BN nanoribbons (ZZBNNR and ACBNNR) were created starting from the previous systems including extra layer of vacuum (10Å) in the y-direction. Since our calculations are based on the DFT approach, the free dangling bonds have to be correctly passivated.
Thus, the edges of the ribbons were passivated including H atoms in order to maintain the appropriate atomic coordination. The H atomic mass was increased to a comparable value as the N and B to perform the integration of the equations of motion using the previous defined time step.  . This figure also shows a direct comparison between the present phonon bands obtained by finite differences from Siesta using a single-ζ basis set (solid line) and those (symbols) calculated within density functional perturbation theory (DFPT) using plane-wave methods as implemented in the PWSCF code [46] taken from ref. [42]. A rather good agreement is obtained specially for the flexural acoustic (ZA), longitudinal acoustic (LA) and transverse acoustic (TA) modes which are the principal phonon modes that contribute to the thermal transport. This agreement stands for the accuracy and reliability of the present computational setup. Similar to other 2D materials [47,48], the LA and the TA branches are linear near the Γ point, while the ZA branch is quadratic. Although the frequencies of the longitudinal and transverse optical branches (LO and TO) are slightly overestimated, these phonons do not contribute significantly to the thermal transport [49]. These small discrepancies, due to the use of a minimal basis set (single-ζ, selected to minimize the very high computational demands of the DFT-AEMD simulations) are expected to have a minor effect on the computed thermal transport properties, while the computational cost is greatly reduced.   3 shows the calculated C QHA v and the fitting to the Debye model (eq. 5). This model is used to fit the C v in the low temperature range, and the calculated Debye temperature Θ D lies between 1150 K and 1690 K, depending on the upper temperature limit of the fit. The discrepancies in the fit at low temperature arise from the assumption that, in the Debye model, the 3 acoustic bands are well described by a linear dispersion relation with a common group velocity v s . This assumption is not true for the ZA phonon band which has a parabolic dispersion and dominates the C v behaviour at low temperatures. The obtained Debye temperature is higher than the reported one for bulk h-BN which is around 1500 K [34]. On the other hand, the Debye temperature Θ D can also be calculated from the average sound velocity v s as [50]

Phonon and thermal properties
where is the reduced Planck constant, N is the number of atoms in the unit cell, S is the area of the unit cell and the avearge maximum sound velocity is given by which accounts the group velocities of the longitudinal (v l ), transverse (v t ) and out of plane (v z ) acoustic modes which are taken from Fig. 2. The calculated Debye temperature using eq. 6 was Θ D = 1030 K which is close to the lowest value calculated directly from the Debye model (1150 K) reflecting that the Debye model is only valid in the low temperature limit. This result can be compared to the Debye temperature of other 2D materials such as graphene, silicene and germanene (2539 K, 680 K and 352 K, respectively) [50,51]. Higher Debye temperatures mean larger phonon velocities and increased acoustic phonon frequencies which suppress phonon-phonon scattering by decreasing phonon populations [52] reflecting higher thermal conductivity. In all 2D materials, the thermal conductivity decreases monotonically with decreasing Θ D . Thus, being the conductivity of graphene in the range from 1000 W/mK to 8000 W/mK [53] and the conductivity of silicene κ ≈ 28 W/mK [54,55,56] and in view of the value of Θ D for SL h-BN, its thermal conductivity is expected to lie in the range between graphene and silicene.

Anisotropy of the thermal conductivity in samples of finite length
We now focus on the possible anisotropy of the thermal conductivity. We first center our study on 2D layers of finite length, which may represent the case of h-BN layers in nanoscale devices. In order to perform a feasible systematic study along several transport orientations, we build samples with different orientations and similar sizes, and compute the thermal conductivity of each of them. Although, due to geometrical constraints, it is not possible to build simulation cells with the exact same lengths for the different orientations, in all cases we have chosen cell sizes close to L y ≈ 1 nm and L z ≈ 20 nm (see the Supporting Information for further system details). Fig. 4 shows the calculated thermal conductivity for different orientations given by 0 • ≤ σ ≤ 30 • , i.e. ranging from zigzag to armchair. We find that there is a sizable dependence of the conductivity with the transport direction. The maximum thermal conductivity is found for the armchair direction, whereas the minimum one corresponds to the zigzag and the 15 • chiral angles. Moreover, there exists a local maximum around 7 • . In order to ensure that the obtained anisotropy corresponds to intrinsic properties of the SL h-BN and that the slight differences in the system sizes do not affect the observed trend, the two extremal chiral angles (zigzag and armchair) were simulated using the shortest and largest length (L z ) of the all orientated samples. This is reflected in two colored regions in Fig. 4 corresponding to the zigzag (red area) and the armchair (green area). Since both regions do not overlap, we can draw the conclusion that thermal anisotropy is present in finite h-BN samples for the sizes considered (around 20 nm). In particular, we obtain a ratio of the conductivities in the armchair and zigzag directions of κ AC /κ ZZ ≈ 1.6. It is interesting to contrast the results for the finite size 2D samples described above with those of h-BN nanoribbons (BNNRs) where edge effects (edges are parallel to the transport direction) are expected to affect the heat flux. Y. -C. Chen et al. [26] reported on the anisotropy of thermal conductivity of BNNRs along various transport directions by Non-Equilibrium Molecular Dynamics (NEMD) using a Tersoff-type potential. They observed a maximum value of thermal conductivity for the zigzag σ = 0 • , two local maxima in the armchair σ = 30 • and in the σ = 19 • nanoribbons, and minima of conductivity for the 13 • and 23 • chiral angles. Their thermal conductivity behavior was explained solving the linearized Boltzman transport equation (BTE) under the relaxation time approximation including an extra scattering term caused from the free lateral edges (parallel to transport direction). This scattering term depends on the edge specularity which is a parameter related to the edge roughness of the ribbons. From their results, the ratio between the armchair and the zigzag directions is κ N R AC /κ N R ZZ ≈ 0.68. This is further supported by the results from other authors, who have reported that, for the nanoribbon case, the most conductive transport direction is the zigzag one obtaining values around κ N R AC /κ N R ZZ ≈ 0.6 − 0.7 [27,21,25]. Our results for the anisotropy of 2D samples with finite length are not affected by any boundary scattering, as in the case of the nanoribbons, since 2D samples are simulated in periodic boundary conditions. Therefore, we conclude that Fig. 4 provides evidence of an intrinsic anisotropy in thermal transport for finite samples, not related to any boundary scattering physics. In order to further confirm this conclusion and to perform a direct comparison with the results for BNNRs from the literature, we have calculated the thermal conductivity of zigzag and armchair BN nanoribbons using our AEMD-Siesta methodology. From the previous systems for the study of the thermal transport along the zigzag and armchair transport directions, the nanoribbons were created adding extra layer of vacuum in the L y direction and including H atoms to pasivate the free edge dangling bonds (see Supporting Information for system details). The calculated thermal conductivity was κ N R ZZ = 11.47 ± 0.91 W/m K for the zigzag and κ N R AC = 8.63 ± 0.85 W/m K for the armchair nanoribbons, respectively. These values for the nanoribbon structures are lower than the corresponding 2D counterpart, implying the existence of extra scattering channels not present in the bulk samples, associated with the ribbon edges. Moreover, for the nanoribbon case, the zigzag direction has larger thermal conductivity than the armchair. The ratio is κ N R AC /κ N R ZZ = 0.75 ± 0.10 in good agreement to all the previous reported results and far from the value obtained for the fully periodic case. We conclude that in the nanoribbon case, the role played by boundaries is twofold, namely: (i) they reduce the thermal conductivity with respect to the bulklike samples (because of boundary scattering); (ii) they change completely the anisotropy behaviour in terms of the κ AC /κ ZZ ratio.

Bulklike thermal conductivity
We are now interested in the bulk-like behaviour of the thermal conductivity for the SL h-BN samples, i.e., the limit for large sample sizes. To obtain it, AEMD simulations must be repeated for samples with increasing length L z , and extrapolate to the limit of infinite L z . We use a second order Taylor expansion for the extrapolation [30,31]: where κ ∞ is the inverse of the thermal conductivity at the bulk limit L z → ∞, and A and B are the coefficients of the linear and quadratic terms, respectively. 1/κ ∞ , A and B are obtained from the fitting of the computed values of κ for different lengths L z to eq. 8. Fig. 5 shows the corresponding 1/κ vs. 1/L z plots and the fitted curves to the simulated data for the zigzag and the armchair orientations. The bulk thermal conductivity values are κ ZZ ∞ = 1003 ± 150 W/mK for the zigzag and κ AC ∞ = 1160 ± 254 W/mK for the armchair orientations. The errors of the fitted κ ∞ reflect both the dispersion around the functional form assumed for the fitting curve, and the error bars associated for each point. These values of the bulk thermal conductivity agree well with the most accurate one of 1050 W/mK obtained by Capellotti et al. [24] solving the linearized BTE equation using DFPT data for h-BN. This value is very close to the average over both directions (κ ZZ ∞ + κ AC ∞ )/2 obtained in this work. In comparison to other two-dimensional materials with the same honeycomb atomic structure such as silicene or graphene, the thermal conductivity of SL h-BN is one order of magnitude larger than that of silicene (28 W/mK) [54,55,56] and much lower than that of graphene (1000-8000 W/mK) [53], which is the most conductive material. It is interesting to note that the thermal conductivity of the SL h-BN is larger than the bilayer (∼ 484 W/mK) [18] and bulk (400 W/mK) [19] counterparts.
From the extrapolated value of 1/κ ∞ for both zigzag and armchair directions, we observe a noticeable but small anisotropy in the bulk limit However, the difference is much smaller than that obtained for finite samples, and for the nanoribbon structures. We, therefore, conclude that, for large area SL h-BN, thermal conduction is nearly isotropic, whereas small samples do present a marked anisotropy. In order to understand this finding, we have analyzed the average square group velocities of the acoustic branches along the high symmetry directions. From the solution of the BTE under the relaxation time approximation, the thermal conductivity for phonons with a given frequency ω can be written as [57] κ ω,α = 1 where N is the number of points of the k-sampling of the Brillouin zone, V is the volume of the unit cell, f 0 is the Bose-Einstein distribution function, τ λ is the relaxation time per phonon mode λ and v λ,α is the group velocity along the α direction. Since κ ω depends on the square of velocity v 2 , a detailed study of v 2 along the armchair (Γ − M k-point bandpath) and zigzag (Γ − K k-point bandpath) directions can reveal the origin of the anisotropy. The average square group velocity in one direction α was calculated with the following expression [57] v 2 The δ function was approximated with a Gaussian one using adaptive broadening parameter depending on the mode group velocity [58]. Fig. 6 shows the average square group velocity along the high symmetry lines calculated from the phonon dispersion curve presented before (see Fig. 2). Only the acoustic phonon branches (ZA, LA and TA) were used for this calculation since, in 2D materials, the thermal transport is practically mediated by acoustic phonons. It is clearly shown that the main differences are in the range of high phonon frequencies whereas for low frequencies, the square group velocity is quite similar. We remark that in 2D materials, the low frequency acoustic phonons play the dominant role in the thermal conductivity. Thus, when the length of the system increases, more low frequencies (long wavelength) phonons are excited and contribute to the thermal conduction [59]. In view of the above and taking into account the calculated v 2 behavior, for shorter systems the conductivity is mediated by phonons with high frequencies which have different group velocities, leading to anisotropic transport behavior along the zigzag and the armchair directions. However, when the length of the system increases, phonons with lower frequencies are excited which have comparable square group velocities along the two directions resulting in a similar thermal conductivity. This qualitative argument is supported by the fact that the value of v 2 in the armchair direction is greater than in the zigzag one in the high frequency range, leading to a higher thermal conductivity in that direction for short samples, as discussed above.

Conclusions
The thermal transport properties of single-layer hexagonal boron nitride were studied using AEMD in combination to DFT. The thermal conductivity along the zigzag and armchair directions as a function of the sample length was presented determining that the value of the bulk κ is around 1050 W/m K for both directions in agreement to recent publications. Moreover, in the bulk limit case, the results point out that the anisotropy in the thermal conductivity is small. A systematic study of the thermal conductivity along several orientations ranging from 0 • to 30 • rotation angle was done for finite samples. From the studied samples (20 nm length), the maximum conductivity is obtained for the armchair direction (30 • ) as well as a local maximum appears at 7 • whereas the zigzag (0 • ) one shows lower conductivity. The observed length dependent anisotropy in the SL h-BN is explained with the average square group velocities as follows: for shorter systems the phonons that contribute to the thermal conductivity are the ones with larger frequencies being v 2 quite different for the zigzag and armchair directions leading to anisotropy in the thermal conductivity. However, as the length of the system increases, more phonons with lower frequencies contribute to κ whose have comparable square group velocities resulting in a similar thermal conductivity.

Code implementation
Each AEMD simulation can be divided in two different parts: first, the temperature step-like profile has to be created. For this purpose, the total system is divided in two subsystems and coupled to external baths to reach the desired temperatures (canonical run). Once they are equilibrated to the initial temperature as it is shown in Fig. S. 1, a microcanonical run is done and the average temperature of each subsystem is monitored in the transient state.
In order to run the AEMD simulations, a modified Siesta version was used. Siesta is a code to perform efficient electronic structure calculations and ab inito molecular dynamics simulations within the DFT framework. It is based on strictly localized basis sets to describe the electron density. A very important feature of the code is that its accuracy and cost can be tuned in a wide range, from quick exploratory calculations to highly accurate simulations matching the quality of other approaches, such as planewave and all-electron methods.
Siesta has several options for MD simulations applying external constrains: temperature or/and pressure. For our purposes, the options to perform the canonical run at given target temperature are: Nosé MD with temperature controlled by a Nosé thermostat, or Anneal MD based on Berendsen thermostat. All the parameters associated to each thermostat can be easily tuned. Finally, for the microcanonical run, the velocity Verlet algorithm is also implemented in Siesta.
Basically, two main changes were applied to Siesta in the MD subroutines: (1) split the total system in two subsystems and (2), define two target temperatures (for the canonical run) and perform the kinetic energy summation for all atoms that belong to each subsystem (i.e. all the loop atoms inside these subroutines were split in two (Upper figure) System decomposition and preparation at given temperatures to create the step-like temperature profile (red and blue areas are the hottest and the coldest subsystems). A side view of the single layer BN is also reported, showing its thermal rippling. In all the simulations, the subsystems contain equal number of atoms. (Lower figure) Step-like temperature profile along the length of the system (L z ). Note that this structure is repeated in the space using periodic boundary conditions (PBC). different loops accounting for the atoms that compose the subsystems). Since the code is written in a modular form, the force calculation core (DFT self-consistent loop and energy derivatives) was not changed. Only the annealing, the Verlet and the initial maxwell Boltzmann velocity generator subroutines were changed to include the two subsystems. Concerning the annealing subroutine based on the Berendsen thermostat, the velocities are usually rescaled according to the scale factor where T is the target temperature, δt is the integration time, τ T is the rise time of the thermostat and T 0 is the instantaneous temperature. However, the pure velocity reescaling scheme is obtained if δt = τ T . In our simulations, the rise time constant was set equal to the time integration step reaching to the target temperature in only one MD step.

System definitions
The different systems under study were created starting from the relaxed unit cell. This unit cell was replicated several times along the y-z direction in order to create the BN supercell. An extra layer of vacuum was added in the direction normal to the plane in order to avoid interactions between consecutive h-BN layers. Since Siesta uses a localized set of orbitals which are non-zero up to a determined radius cut-off, the distance between layer images was set to twice the maximum cut-off radius of the orbitals.
Concerning the study of the thermal conductivity as a function of the length of the sample, the width of the zigzag systems were L y = 0.869 nm and L z ranges from ≈ 10 nm to ≈ 100 nm. The number of atoms in the systems varied from 320 to 3200. Similar numbers were obtained for the armchair case but the width of the sample was L y = 1.004 nm.
For the anisotropy study in finite samples, the system sizes were around L z ≈ 20 nm and L y = 0.86 nm for length and width, respectively. The sizes as well as the number of atoms of each system are presented in Table S   The two studied nanoribbons were created from the previous structures used in the anisotropy study (see the scheme shown in Fig. S. 2). Basically, a layer of vacuum was added in the y-direction (following the same prescription as before to avoid interaction between consecutive replicas). Moreover, H atoms were included in the lateral edges of the ribbon to passivate the free dangling bonds. Finally, the atomic mass of the H atoms was increased to a comparable value of the B and N mass in order to use the same integration time step to solve the equation of motions.