Electrical and thermal transport in coplanar polycrystalline graphene-hBN heterostructures

We present a theoretical study of electronic and thermal transport in polycrystalline heterostructures combining graphene (G) and hexagonal boron nitride (hBN) grains of varying size and distribution. By increasing the hBN grain density from a few percents to $100\%$, the system evolves from a good conductor to an insulator, with the mobility dropping by orders of magnitude and the sheet resistance reaching the M$\Omega$ regime. The Seebeck coefficient is suppressed above $40\%$ mixing, while the thermal conductivity of polycrystalline hBN is found to be on the order of $30-120\,{\rm W}{\rm m}^{-1}{\rm K}^{-1}$. These results, agreeing with available experimental data, provide guidelines for tuning G-hBN properties in the context of two-dimensional materials engineering. In particular, while we proved that both electrical and thermal properties are largely affected by morphological features (like e.g. by the grain size and composition), we find in all cases that nm-sized polycrystalline G-hBN heterostructures are not good thermoelectric materials.

Introduction. Owing to a small lattice mismatch (2%), graphene and hexagonal boron nitride can be assembled in coplanar twodimensional heterostructures (1 ) . Such atomic sheets, covering a wide range of compositions, result in new materials with properties complementary to those of graphene and hBN, such as tunable bandgap optoelectronic materials (2 ) . Graphene is well appreciated for its high electrical (3 ) and thermal conductivities (4 ) , whereas hBN is an electrical insulator with to date an unmeasured thermal conductivity (5 ,6 ) . Large-scale coplanar G-hBN heterostructures have been successfully fabricated using chemical vapor deposition (CVD), enabling the possible control of periodic arrangements of domains whose sizes range from tens of nanometers to millimeters (7 -10 ) . Their charge transport properties can be, however, quite surprising, such as the presence of a metal-insulator transition (11 -13 ) and anomalous transport phenomena, that is not fully understood. (14 ) Additionally, fast CVD growth results in polycrystalline materials with grains of varying sizes and morphologies, and the electronic and thermal properties of these materials are limited by the presence of grain boundaries (GBs) (15 -19 ) .
In polycrystalline graphene, GBs are characterized by Van Hove singularities near the Dirac point (20 -22 ) , whereas in hBN the GBs reduce the bandgap and introduce gap states generated by the presence of B-B or N-N bonds (23 ) . The interface between G and hBN is also expected to give rise to local boundary states, especially at low energies (24 ,25 ) . GBs are also usually accompanied by local structural deformation, which enhances phonon scattering and thus lowers thermal conduction. The thermal properties of polycrystalline graphene have been theoretically calculated using molecular dynamics simulations as a function of average grain size (17 ,18 ,26 ,27 ) , in fair agreement with experimental results (4 ) .
Recently, a sample of CVD-grown graphene was gradually converted into hBN, and it was observed that chemical substitutions are initiated around structural defects. This process of conversion demonstrated a fine tunability between highly conductive graphene and insulating hBN (13 ) . To date however, the electronic and thermal properties of CVD-grown hybridized G-hBN heterostructures are poorly understood, and their potential use in energy harvesting, optoelectronic, or nanoelectronic applications remains unclear.
Here we use quantum transport and molecular dynamics (MD) simulations to calculate the electronic and thermal properties of polycrystalline G-hBN heterostructures with varying grain size and distribution. The electronic mobility and sheet resistance are studied as a function of the density of hBN grains, which ranges from a few percent to full coverage. The contribution of GB interface states to the transport properties is also illustrated and quantified. By performing a complete calculation of thermal and electrical transport, we estimate the thermoelectric conversion ratio and find that it remains far too low to be useful for energy harvesting applications.
Generation of samples. Polycrystalline G-hBN heterostructures with uniform average grain size were generated using a Voronoi algo-rithm, resulting in large square periodic samples containing up to 3 million atoms (17 ,28 ) . The algorithm starts with a random selection of nucleation centers within a square cell of predefined dimension, which dictates the average grain size as L grain = L 2 /n grains , where L is the sample length and n grains is the number of grains. Next we set a random crystal orientation for each nucleation site and we use a Voronoi method to construct the grains. The atoms along the GBs with separation below 0.1 nm are removed, and an MD annealing process is used to construct the GBs, setting all the atoms as carbon. We use the LAMMPS simulation package (29 ) , the second-generation reactive empirical bond order potential (30 ) , and a small time increment of 0.1 fs. The annealing starts with a 3-ps equilibration at room temperature using the Nosé-Hoover thermostat, continues with a heating up to 3000 K for 12 ps and keeping this temperature for 3 ps, and ends with a cooling back down to room temperature for 10 ps. Finally, based the concentration of hBN, we assign which grains are graphene and which ones are hBN ( Figure 1).  Electronic properties. We describe the electronic properties of the G-hBN heterostruc-tures with a tight-binding Hamiltonian where ε i (r i ) is the on-site potential of each atom and t i,j is the hopping between nearest neighbors. In systems containing 1D interfaces between two different 2D materials, the electronic properties are sensitive to the interface termination, and thus care must be taken when describing the GBs between graphene and hBN grains. For example, zigzag BN nanoribbons are polar, presenting bound charge of opposite signs at the B and N edges. In hybrid systems, mobile electrons from the graphene will tend to screen the excess interfacial charge, which changes the potential profile across the GB. Therefore, we modify the on-site term of the Hamiltonian to include a position-dependent electrostatic potential, which can be derived from the screened Poisson equation considering point charges midway between the C-B or C-N interfacial bonds. The on-site term of the TB Hamiltonian can then be written as (31 ) where ε i (r i ) denotes the on-site energy for an atom of type i (either carbon, boron, or nitrogen) at position r i , ε i0 is the on-site energy of atoms far from the GBs, is the position of the excess charge at the C-B (C-N) interface, λ i is the decay length of the interface potential, and the sum is done for all n q charges within a radius of 1 nm. The onsite potential and nearestneighbor hopping parameters have been derived from a Wannierization of DFT calculations and are given in Table 1 (see Supplementary Information for more details). Finally, because the GBs contain non-hexagonal rings, B-B or N-N bonds will be present. For these bonds we set t BN as the hopping parameter, while the on-site energy is taken as 1.1ε i0 .
We calculate the electronic density of states (DOS) using the Lanczos recursion method with an energy resolution of η = kT = 26 meV (T = 300 K). Figure 2(a) shows the DOS with increasing hBN grain density in steps of 20%, for an average grain size of 40 nm. The gap is seen to progressively widen with increasing hBN concentration, but with a faster decay on the electron side of the spectrum. This electron-hole asymmetry stems from the GB states, which generate more resonances on the electron side. This can be seen more clearly for 100% hBN, where the formation of boundary states, with energy lying inside the gap, is illustrated by the local density of states projected over all the GB sites (LDOS GB ; Figure 2(b)). The energy resonances at -1.2 and 2 eV has been observed experimentally, which can be associated to homoelemental bonds in the GB. (23 ) Besides, we observe other peaks at E = 0.0 and 0.76 eV, both are found for polycrystalline graphene and hBN, which suggest specific fingerprints of the structural morphology of grain boundaries. These states are mainly localized at the GBs, as visualized in the inset of Figure 2(b), with stronger energy resonances on the electron side of the spectrum (see additional LDOS GB projected around a G-hBN interface in Supplementary Information). The presence of such states could be at the origin of the finite electrical conductivity computed for polycrystalline hBN (see below).
We next evaluate the electronic transport properties using a real-space order-N wave packet propagation method (32 ,33 ) . The core of this method is to calculate the time-dependent diffusion coefficient as where ∆X 2 is the mean-square displacement of the wave packet and ρ(E) = Tr[δ(E −Ĥ)] is the DOS. We eval-uate the trace using the Lanczos recursion and the same parameters as the DOS. We calculate the energy-dependent semiclassical conductivity as σ(E) = e 2 ρ(E)D(E), whereD(E) is the value of the diffusion coefficient when the mean displacement has reached six times the average grain size (see Supplementary Information).
In Figure 3(a) we report σ(E), where a drop of more than two orders of magnitude is observed near the charge neutrality point with increasing hBN concentration. To further clarify the impact of the density of hBN grains, we fix the carrier concentration to n = 0.3 × 10 12 cm −2 , which is a typical value for graphene on SiO 2 (34 ) , and evaluate the charge mobility µ = σ(n)/n, shown in Figure 3(b). The sheet resistance R is shown in the inset of Figure 3(b), where one can see that the maximum value for 100% hBN is about 1 MΩ. Experimentally, a sheet resistance of a few GΩ has been measured, (13 ) ; this value is consistent with our calculations as a consecuence of the grain size scaling of the electrical properties. (19 ) Additionally, we estimate the GB-resistivity, ρ GB , using an ohmic scaling analysis (16 ,19 ) , where R and R 0 are the sheet resistances of the polycrystalline sample and the individual grains, respectively. The estimated resistivity for the G-G interface is 0.12 and for hBN-hBN is 5.93 kΩ · µm (see Supplementary Information).
To complement the information about the electronic properties, we evaluate the Seebeck coefficient where G is the sheet conductance and f is the Fermi distribution. As shown in Figure 3(c), the Seebeck coefficient of the polycrys-talline samples is reduced compared to pristine graphene (35 ) , but is insensitive to hBN concentrations below 40%. However, beyond 40% the thermoelectric capability is strongly suppressed. Thermal properties. In order to evaluate the thermal conductivity as a function of the grain size, we construct a finite element (FE) model in the ABAQUS package with 4000 grains constructed as Voronoi cells (right panel Figure 4(a)). Using six representative pentagon-heptagon GB structures, we extract the GB thermal conductance for G-G, G-hBN and hBN-hBN interfaces by performing a nonequilibrium molecular dynamics (NEMD) calculation with LAMMPS (see Supplementary Information for details); which are introduced as contact conductances between interfaces. In the FE model, we include two highly conductive strips at the two ends of the structure (17 ,28 ) and fix the ingoing (outgoing) heat flux on the left (right) side, h f . Then, we evaluate the steadystate temperature profile along the sample and use the ∆T between the strips to evaluate the effective thermal conductivity of the sample as where L is the sample length. We calculate the themal conductivity for 16 grain sizes between 1-1000 nm while changing the concentration of hBN (Figure 4(b)). The scaling of κ shows that the impact of the GBs on thermal transport becomes negligible for grain sizes above 100 nm, which suggest that heat carriers with mean free path longer than 100 nm bring low contribution to κ. Figure 4(c) displays the thermal conductivity as a function of the hBN grain density where we observe that, for small average grain size, the minimum of thermal conductivity occurs near 70% hBN, similar to prior estimates (36 ) . This minimum can be rationalized by the fact that the thermal conductance for the G-hBN interface is lower than that of the hBN-hBN and G-G interfaces. For larger grain sizes, where the GBs no longer dominate the thermal transport, we observe a monotonic scaling of κ with hBN grain density, as the thermal conductivity of pristine hBN is lower than that of pristine graphene. In order to validate the above FE analysis we perform an independent investigation based on MD simulations. The goal is to provide evidence that the FE analysis, although missing most of the atomic-scale details, nevertheless provides the correct gross features on thermal transport across hBN and graphene GBs. The fully atomistic study of the thermal conductivity employs an approach-to-equilibrium molecular dynamics (AEMD) method following the same approach as in Ref. (18 ) using the Tersoff BNC potential (37 ) (see Supplementary Information). While NEMD provides direct access to the temperature drop across the GB, which is the relevant quantity needed to calculate the interface thermal resistance (and, therefore, the GB conductance), the AEMD ap- proach is better suited to calculate the effective κ in a large system, since it requires a comparatively smaller computational effort (38 ) . We observe a quantitative difference between the approach described above and AEMD, which is reflected in the extracted value of the thermal conductance of the hBN-hBN interface, C hBN−hBN = 5.27 GW/m 2 K. From the data reported in Ref. (18 ) , we also estimate the thermal conductance of the G-G interface to be C G−G = 12.66 GW/m 2 K. We attribute these lower values to the structure of the GBs investigated; the GBs in the AEMD calculations tend to be disordered and meandering, as shown in Figure 1, while the GBs used in the NEMD method were mirror symmetric and perfectly periodic arrays of pentagon-heptagon pairs. The inset of Figure 4(b) shows the thermal conductivity of polycrystalline graphene and hBN using the FE method and the GB conductances extracted from the AEMD method. The smaller values of GB conductivity manifest themselves in a lower overall thermal conductivity, but the main trend holds, and a grain size of 100 nm still appears to be the crossover where thermal transport is no longer dominated by the GBs.
To summarize, we have presented an electrical and thermal characterization of coplanar G-hBN heterostructures. The tight-binding model includes a refined description of the G-hBN interfaces, and is used to describe the electrical properties of polycrystalline structures with varying percentages of graphene and hBN. Our results reproduce the transition from graphene to insulating hBN, with an electrical conductivity change of more than two orders of magnitude and a strong suppression of the Seebeck coefficient. Additionally, the thermal conductivity of these polycrystalline structures has been investigated using a combination of atomistic MD simulations and a FE evaluation of the heat equation. We find that for small-grain structures, the thermal conductivity is minimized for a hBN grain density of 70%. From our study, we can evaluate the upper value of the thermoelectric figure of merit, ZT = σS 2 T /κ. For example, in the case of 40 nm average grain size and 20% hBN, ZT ∼ 1 × 10 −4 for a carrier concentration n = 5 × 10 12 cm −2 , which is quite small. Even for energies near the edge of the gap, where the Seebeck coefficient should be maximized, the value of ZT only reaches ∼ 1 × 10 −2 . dad de Cataluna and the Severo Ochoa Program (MINECO SEV-2013-0295). M.P. and L.C. acknowledge Spanish MINECO (FIS2015-64886-C5-3-P) and Generalitat de Catalunya (2014SGR301). B.M. and T.R. greatly acknowledge the financial support by European Research Council for COMBAT project (Grant number 615132).

Supporting Information Available
The following files are available free of charge. Details of tight-binding model, details on the numerical evaluation of electrical and thermal conductivity, and scaling analysis to estimate electrical and thermal GB resistivity.