May 24, 2018 - Schulz et al.43 of 0.34 Ã
. LEED-I(V ) and nc-AFM data for the moirÃ© ...... Center (Aalto NMC) facilities and was supported by the. Ac...

0 downloads 4 Views 1MB Size

Fabian Schulz† and Peter Liljeroth

arXiv:1805.08880v2 [cond-mat.mtrl-sci] 24 May 2018

Department of Applied Physics, Aalto University School of Science, P.O.Box 15100, FI-00076 Aalto, Finland (Dated: May 25, 2018) There has been enormous interest in weak, van der Waals-type interactions due to their fundamental relevance in the field of two-dimensional materials and the so-called van der Waals heterostructures. Tackling this problem using computer simulation is extremely challenging due to the non-trivial, non-local nature of these interactions. We benchmark different treatments of London dispersion forces within the density functional theory (DFT) framework, using a layer of h-BN adsorbed on Ir(111) as a prototypical weakly-bound interface. We calibrate our calculations with non-contact atomic force microscopy (nc-AFM) measurements as well as previous experiments. In addition, we provide results for graphene adsorbed on the same substrate, gr/Ir(111). Our results show strong variations in the calculated atomic geometry, originating from the approximative character of treatments of the the weak interactions. While some approximations reproduce the experimental structure, this is rather based on “a posteriori” comparison with the “target” results. Thus, we are forced to conclude that van der Waals treatment in DFT is currently at an empirical stage and achieving true predictive power calls for new approaches.

The tremendous success of density functional theory (DFT) together with the continuously growing computing power have laid out the path for the rapidly expanding field of computational materials science. Rather than being employed only in a complementary manner to support the interpretation of experimental data, DFT is now widely used to screen large numbers of structures and compounds,1–7 with the goal to guide experimentalists in their search and synthesis of new functional materials. In addition, DFT for computational materials science is increasingly combined with machine learning methods,8–12 which further increases automation and decreases human oversight. This growing trust in and reliance on DFT calculations suggests wide-spread predictive power. However, at the current state of DFT, such an assumption should not be made a priori but checked and validated for each class of materials and their properties separately. For example, DFT-based materials discovery has been extensively used for two-dimensional materials4,7 , which have recently sparked intense theoretical and experimental interest. These layered materials can be combined in the so-called van der Waals (vdW) heterostructures,13–17 which enables to tune their structural and electronic behaviour. Because of the absence of covalent bonds between the layers in such heterostructures, vdW interactions are crucial in determining their properties, in contrast to normal solids.18–21 For example, band hybridization due to interlayer hopping can lead to emergence of superconductivity or strongly-correlated insulating behaviour in certain vdW stacks.22,23 Another field of research where the interplay of vdW interactions with wave function hybridization is of paramount importance is the adsorption of lowdimensional objects on metal surfaces. Here, the balance

between these two interactions governs whether an adsorbate is physisorbed or chemisorbed, and thus also determines adsorption geometry and interfacial electronic properties.24–26 The importance of the vdW interactions - or more precisely the London dispersion forces - in such systems has been recognized for over a decade. However, their computational modelling with the DFT method remains challenging, because these interactions are purely nonlocal in nature and much weaker than the interactions between chemically bound atoms. Nevertheless, in recent years there has been considerable progress in the treatment of vdW interactions within DFT (in the following referred to as vdW DFT methods). This ranges from semi-empirical approaches with an explicit interaction term depending on the ionic coordinates and usually no dependency on the electronic structure, to density functionals that contain a non-local dependency of the correlation energy on the density.27–29 Briefly, in the former case, the system is forced into a geometry that does not constitute a minimum of the total energy expression of the electronic structure. In the latter, the exact form of the vdW energy is known only in the asymptotic limit, which does not correspond to the realistic case where the tails of the electronic densities of the constituents already start to overlap. Therefore, approximations have to be made, and the final expression is not unique, similarly to the collection of generalised gradient approximations (GGAs). Consequently, these approximations and their combinations with different exchange and correlation functionals need to be tested and validated. Many benchmark studies of vdW DFT focus on small molecular systems, where reference data can be obtained from highly accurate quantum chemistry calculations

2 (see Refs. 27–29 and references therein). Assessing vdW DFT methods for larger, more complex systems is challenging, because their size is prohibitive for quantum chemistry methods and thus reliable experimental data for comparison is required. In addition, the large system size usually allows only for a small set of vdW DFT methods to be tested. Examples of such studies include the calculation of the structure of sheet silicates30 or the adsorption of benzene31 or graphene32,33 on different metal surfaces for cases where the lateral extension of the unit cell is small. Here, we apply DFT to two experimentally wellcharacterized systems of two-dimensional layers adsorbed on a metal surface, monolayer hexagonal boron nitride (h-BN) on Ir(111) and graphene (gr) on Ir(111), to benchmark the different approximations to the exchangecorrelation and vdW terms. These two systems are very challenging for vdW DFT for two reasons: (i) Due to the lattice mismatch between h-BN and gr and the substrate, moir´e patterns are formed where the stacking between overlayer and metal surface varies continuously. This leads to very large unit cells, which makes the calculations computationally very expensive. (ii) In addition, the alternating atomic registry with the substrate results in different adsorption strengths along the moir´e unit cell and consequently, a corrugation of the h-BN and gr overlayers. The final geometry of the moir´e superstructure then crucially depends on the subtle interplay between vdW interactions and wave function overlap leading to charge transfer and chemical bonding. This second point makes these two systems ideal candidates to compare and benchmark vdW DFT, because the modelling of these effects depends on the chosen exchange-correlation functional and the treatment of the vdW interactions therein. We test a large number of different vdW DFT methods, and for some of them, we also investigate the influence of varying parametrization, e.g. such as different starting geometries or number of substrate layers. The calculated geometries are compared to experimental results, in particular the adsorption height and corrugation of the overlayers, which represent very sensitive measures for the predictive power of vdW DFT. For the experimental reference values, we draw upon previously published data as well as a new set of noncontact atomic force microscopy experiments conducted specifically for this study.

RESULTS

Previous studies of h-BN/Ir(111) and gr/Ir(111) Before we can assess the different vdW DFT methods, we need to establish a reference to compare the DFT results with. In the following, we provide a brief overview of experimental studies on the adsorption height and moir´e corrugation of h-BN/Ir(111) and gr/Ir(111), as well as some of the previous DFT results (see Table I). Determining the geometry of such systems is challeng-

ing not only from a computational point of view; accessing the adsorption configuration experimentally in a reliable manner is difficult as well. Typical approaches include sample averaging techniques such as dynamic low-energy electron diffraction (LEED-I(V )) and x-ray standing wave measurements (XSW) and local probes such as scanning tunneling microscopy (STM) and noncontact atomic force microscopy (nc-AFM). These techniques have been extensively demonstrated and compared with DFT calculations for gr/Ir(111).34–40 Here, XSW measurements by Busse et al.35 yielded a mean adsorption height of 3.38 ˚ A, in good agreement with DFT calculations presentend in the same publication. For the moir´e corrugation, XSW yielded 0.4 to 1.0 ˚ A, depending on the coverage (where the largest value corresponds to a full gr monolayer).38 While the lowest value of 0.4 ˚ A is in agreement with LEED-I(V ) and ncAFM results (to be discussed below) and some of the DFT data35–37,39 , the remaining values are at odds with most other results, a discrepancy which was attributed to stress in the graphene layer.38 The analysis of XSW data is also not straight-forward, as it requires assumptions about the height distribution within the gr (or h-BN) layer and the quality of the prepared surface. While the local probes allow in principle direct access to the moir´e corrugation, STM topography is always a convolution of structural and electronic sample properties (for gr/Ir(111) and h-BN/Ir(111), this effect can lead even to an inverted apparent moir´e corrugation41–43 ). In contrast to STM, nc-AFM measurements are expected to yield values approaching the topographic corrugation.36,41,44,45 The most sophisticated nc-AFM experiments on gr/Ir(111) were carried out by H¨am¨al¨ainen et al.36 , who used a carbon monoxidefunctionalized tip to probe the repulsive force regime above the sample and measured a moir´e corrugation 0.47 ˚ A. In the same publication, a LEED-I(V ) study was reported, which yielded a corrugation of 0.43 ˚ A and a mean adsorption height of the gr layer of 3.39 ˚ A.36 The corrugation is in very good agreement with the nc-AFM value, and there is overall good agreement with many of the DFT results. In the case of h-BN/Ir(111), XSW measurements carried out by zum Hagen et al.40 yielded a moir´e corrugation of 1.55 ˚ A, in agreement with DFT calculations presented in the same publication of 1.50 ˚ A. It is also close to the DFT value reported by Liu et al.39 of 1.40 ˚ A, but significantly larger than the corrugation obtained by Schulz et al.43 of 0.34 ˚ A. LEED-I(V ) and nc-AFM data for the moir´e adsorption geometry of h-BN/Ir(111) are lacking thus far. nc-AFM measurements on h-BN/Ir(111). To add a local probe measurement of the h-BN/Ir(111) moir´e corrugation, we have carried out low-temperature nc-AFM experiments, where tip-sample forces are measured as the frequency shift ∆f of a deliberately oscillating cantilever.46,47 In order to minimize the influence of local variations in chemical reactivity over the h-BN

3 TABLE I. Structural parameters from previous DFT calculations and experiments on h-BN and gr on Ir(111). Method h-BN vdW-DF2-rB86 revPBE+D3 PBE+D2 XSW nc-AFM gr

vdW-DF vdW-DF2-rB86 PBE+D2 XSW LEED-IV nc-AFM

corrugation ∆BN /∆CC 1.50 0.338 1.4 1.55 1.65

distance z h−BN -z Ir1 3.24 3.187 3.8 2.20, 3.72a N/A

Reference 40 43 39 40 present work

0.35 0.36 0.2 0.4 - 1.0b 0.43 0.47

3.41 3.43 4.2 3.38 3.39 N/A

35 40 39 35 and 38 36 36

a

Average adsorption height for the strongly and weakly interacting species, respectively. b Coverage-dependent, where 1.0 ˚ A was found for a full monolayer of gr.

FIG. 1. (a) Constant-∆f nc-AFM image of h-BN/Ir(111) acquired with a CO-functionalized tip. Setpoint: -12.0 Hz. Scale bar is 1 nm. (b) Color-coded ∆f (z) curves recorded at the positions marked in panel (a). Inset: Second-order polynomial fits (red) around the ∆f minima.

moir´e48,49 , we have chemically passivated the tip apex by controlled functionalization with a carbon monoxide (CO) molecule.36,44,50,51 The details of the nc-AFM experiments and the sample preparation can be found in the Methods section. Figure 1a shows an nc-AFM image of h-BN/Ir(111) recorded in the constant-∆f mode with a CO-passivated tip, yielding the expected moir´e superstructure composed of depressions arranged in an hexagonal lattice. In addition, the weakly adsorbed regions of the h-BN show atomic contrast, revealing the honeycomb lattice of the B and N atoms. There is no atomic contrast in the moir´e depressions due to long-range vdW interactions between the tip and Ir(111) substrate, which affect the tip-sample distance feedback.41 Figure 1b shows the measured frequency shift as a function of the tip-sample distance [∆f (z)] at a moir´e depression and in the surrounding region. It can be seen that the ∆f minimum above the depression is not only located at a smaller tip-sample distance but also yields a more negative ∆f value. In regions where the adsorption height is smaller, the contribution of long-range vdW interaction between tip and

iridium substrate is increased, shifting the ∆f (z) curve to more negative frequency shift values.41 Thus, when imaging on the attractive branch of the ∆f (z) curve, even in the constant-∆f mode with a chemically inert tip apex, the tip-h-BN distance will be larger over the moir´e depressions than over the surrounding regions. Consequently, ∆f setpoints sufficient to achieve atomic resolution on the latter region might not yield atomic contrast on the former. At very small tip-sample distances, however, the frequency shift is governed by repulsive interactions and it steeply increases as the distance is further reduced; thus the influence of long-range vdW interactions becomes less pronounced. Indeed, it was suggested that the tip-sample distance corresponding to the minimum in ∆f (z) curves could be used to measure relative differences in adsorption heights.52 Following this reasoning, the inset in Figure 1b shows second-order polynomial fits to the minima, from which we extract their vertical difference as 1.41 ˚ A. Note that this value is significantly larger than the corrugation observed in the constant ∆f image in Figure 1a. We have carried out analogous measurements with two other CO tips on two other regions of the sample, which yielded differences of 1.70 and 1.85 ˚ A, respectively. Taking the average of the three values, we get a moir´e corrugation of (1.65±0.23) ˚ A, in agreement with the value obtained by zum Hagen et al. from XSW (1.55 ˚ A)40 , as well as with recent DFT results.39,40,53 DFT: (1 × 1)-h-BN/Ir(111). We now turn to the DFT calculations and how well they can reproduce the experimental geometries. We will devote the major part to h-BN/Ir(111), because for this system there is a larger spread in the reported DFT results, and thus some controversy. First, we explore the importance of the adsorption site, and the convergence of the results in the moir´e structure by first considering the artificially commensurate (1 × 1)-h-BN/Ir(111) structure. Despite the strain caused in the h-BN layer, this procedure should give information on the differences in the layer height, local buckling and energies at different lateral positions of the overlayer.54 This will naturally be tested later with the calculations of the full moir´e structure. We also show results obtained with the local density approximation (LDA) even if it does not include vdW interactions or they are not approximated with an additional term, but we want to include the LDA because it is still used in calculations involving graphene. Figure 2 contains the binding energy Eb and the height zN of the N atom above the topmost Ir(111) substrate layer. Here the experimental bulk lattice constant aexp of iridium was used (further details are given in the Supplementary Information (SI)). The immediate observation is that there is not only a large variation in the adsorption strengths of the hBN layer on various lateral adsorption configurations, but also with different treatments of the exchange and correlation (XC). Values of Eb range from almost −0.1 to −0.9 eV. zN varies similarly: Depending on the choice of the XC at the preferred lateral adsorption sites, where N takes the on-top site above the substrate atom, we ob-

4

FIG. 2. Binding energy Eb (top) and height zN of N above the top-most substrate layer (bottom) of (1 × 1) commensurate layer, at the experimental lattice constant of Ir(111) from QE-DFT calculations; thus h-BN is stretched; lateral adsorption sites are abbreviated as “f” = fcc, “h” = hcp and “o” = on-top with respect to the fcc(111) termination of the surface.

tain values from 2.25 to 3.75 ˚ A. The variations within Eb and zN are correlated, so that the approximations to XC that yield the lowest average Eb and/or smallest variations among the different lateral arrangements also place the h-BN layer furthest away from the substrate; the weakest adsorption is found with the “original” vdW-DF and vdW-DF2 density functionals. These are also the approximations that yield the largest lattice constants (cf SI). These two approximations lead in general to too weak binding among the atoms within the system. These shortcomings have been addressed in subsequent approximations by “fitting” the lattice constants and later more sophisticated quantities, e.g. the adsorption energies of molecules on surfaces. Considering the variations in Eb and zN with a given treatment of XC across the six different high-symmetry adsorption sites, we see that irrespective of the approximation, the placement of N atom on the on-top site yields the strongest binding, with the Bf No slightly preferred over Bh No in energy. The other sites are energetically practically degenerate, even if the Bo Nf and Bo Nh tend

to reside somewhat closer to the substrate than the arrangements where both atomic species are located at the hollow sites. Thus, the on-top site of either species yields closer adsorption geometry. Independent of the treatment of XC and the lateral adsorption registry, the B atom is closer to the substrate than the N atom. Given that in the moir´e structures the B and N atoms explore the lateral potential energy surface, not only the high-symmetry sites, in a more or less continuous manner, we can use the trends in Eb and zN to “predict” the expected trends in the moir´e structures: This suggests small corrugation and large adsorption height with the “original” vdW-DF and vdW-DF2 approximations to XC, and the strongest binding and thus the smallest adsorption height with the PBE+D2 and revPBE+D2 approximations; largest corrugations would then be expected with PBE+TS, vdW-DF with C09, Cx, optB86b exchange functionals, and with PBE+rVV10 approximation. The results from our two set of calculations with Quantum ESPRESSO (QE) and the PBE+TS treatment

5 (“PBE+TS”, “PBE+TS 120 Ry”) are independent of the cut-off energy used, indicating a good convergence already at the lower value. Also increasing the number of layers from four to seven, tested with the “vdW-DFoptB88” treatment of the XC, leads to the same results. Scaling the lattice constant of the substrate by 12/11 in “vdW-DF-optB88-comm” (compared to the the equilibrium lattice constant of the bulk), so that the h-BN is close to its equilibrium, leads to a very different structure, with almost no preference of the lateral registry. Reducing the number of k points to a grid of 8×8 points leads to small changes in the relative adsorption energies and corrugation. DFT: 12-on-11 h-BN/Ir(111) moir´ e superstructure. We model the moir´e structure as commensurate, such that 12×12 cells of h-BN match 11×11 cells of Ir(111); we use the notation 12-on-11 in this case. We focus here on the most central results from the DFT calculation, the height of the B and N atoms above the substrate and the corrugation within the h-BN layer, given in Figure 3. With some treatments of the XC, we include two results, originating from the two different procedures of atomic relaxation described in the Methods section. We also include the experimental XSW determination40 of the minimal and maximal height of the h-BN layer with dashed horizontal lines, and the corrugation exctracted from the nc-AFM experiments described above; these act as reference values for the DFT results. We note that to be more consistent with the experimental XSW analysis, we should use distances from the outer-most layer coordinates of the substrate that have been extrapolated from the bulk coordinates rather than the actual relaxed surface coordinates. We do, however, refer to the latter in trying to minimise the possible confusion of different values, and since the difference is small: For example with the CP2k code and the vdW-DF2-rB86 treatment of the XC approximation, the average distance between the two outer-most substrate layers turns out to be 2.20 ˚ A, whereas in bulk the distance is 2.22 ˚ A. Also the magnitude of the corrugation is independent of the reference point of the individual minimum and maximum. We begin to decipher the results with two general observations that characterize our qualitative DFT results: (i) There is a wide range of heights depending on the different treatments of the XC; the h-BN layer is furthest away from the substrate with the same approximations to XC as in the (1 × 1) structure in Figure 2, namely vdW-DF and vdW-DF2. In addition, BEEF-DF2 leads to a large average height. (ii) With some treatments of the XC, we find two structures where the forces vanish. Thus there is one stable and one metastable structure. These differ in the magnitude of the corrugation in the h-BN layer. The large variance of the adsorption heights indicates that the approximations to the exchange-correlation term are not reliable per se, but a “calibration” with good experimental data is necessary. Afterwards, the transferability of the chosen approximation needs to be tested.

That there are two structures which are (meta-)stable in the calculations is an interesting result. This raises the question if both structures are indeed realistic. The approximations that yield the two structures are PBE+D2, PBE+D3 and BEEF-DF2. Some approximations, in particular vdW-DF-rB86, vdW-DF2-C09 and vdW-DF2rB86 only yield the structure with large corrugation after ionic relaxation. These functionals also result in the best agreement with the experimental corrugation and minimal and maximal heights as measured with nc-AFM and XSW. This suggests that they have the highest “accuracy” in the present system and obtaining two different structures is probably unrealistic. The value of the binding energy per BN unit, defined as Eb = − {E(h-BN/Ir) − [E(h-BN) + E(Ir)]} /#nBN with the total energies E from the relaxed calculations of the adsorbate system and the free constituents, is 0.164, 0.106 and 0.145 eV with the DFT+D3, vdW-DF and vdW-DF2-rB86 treatments of the XC effects, respectively. These values are intermediate of the results in the (1 × 1) cell above. The values are somewhat smaller than the 0.174 eV obtained in recent calculations.40 DFT: 13-on-12 h-BN/Ir(111) moir´ e superstructure. Even if the closest commensurate periodicity of the aligned moir´e pattern has now been established to be 12-on-11 (Ref. 40 and 43), we studied also the larger 13-on-12 moir´e pattern to compare with an earlier DFT calculation.43 The results are also illustrated in Figure 3; the numerical values are also listed in the SI. The resulting geometry is practically the same in both 12-on-11 and 13-on-12 structures with a given treatment of the XC term. Thus, it is clear that the origin of the different corrugations in different DFT calculations is the XC term, and not the moir´e periodicity as was suggested earlier40 (with a possible contribution from the different numerical parameters in the calculations). DFT: 12-on-11 h-BN/Ir(111) electronic structure. In order to investigate the influence of different treatments of the XC to the electronic structure, we evaluated the difference in the electronic density due to the adsorption of the h-BN onto the substrate; this density is then averaged along the surface: With the electron density nh−BN/Ir (~r) of the full system, we subtract the individual electron densities nIr (~r) and nh−BN (~r), with the respective atomic coordinates of the full, adsorbed case, Z n h io ∆n(z) = nh−BN/Ir (~r) − nIr (~r) + nh−BN (~r) dx dy . x,y

∆n(z) evaluated with this formula from calculations with PBE+D3 and vdW-DF2-rB86 treatments of the XC are shown in Figure 4. Also PBE+D3 electron density at the atomic positions of vdW-DF2-rB86 is shown (“PBE+D3(@vdW-DF2-rB86)”). The differences between the treatments are relatively small, even considering that in the case of PBE+D3 the electronic structure

6

FIG. 3. Height of all B and N atoms in 12-on-11 structure above the average height of the top-most layer of the substrate: Top panel, experimental, bottom panel DFT-derived lattice constant; the XSW result by zum Hagen et al. from Ref. 40 is marked with the lateral dashed lines, minimum and maximum. Blue horizontal lines mark the minimum and maximum values for a 13-on-12 structure using the revPBE+D3 approximation from Ref. 43, and the blue vertical bar on the right indicates the corrugation obtained from the analysis of the AFM data in this publication. With some treatments of the XC we have performed two calculations with different starting geometries – see Methods; in this case the second structure is drawn with green markers. The grey-shadowed regions refer to calculations in the 13-on-12 cell, the blue to different calculations with the same approximation.

corresponds to the one from PBE, but with atomic positions forced away from the DFT minimum with the D3 term. A notable difference is the larger enhancement of the electron density in between the ad-layer and the substrate with vdW-DF2-rB86, in particular when the same atomic coordinates have been used. DFT: 10-on-9 gr/Ir(111) moir´ e superstructure. We have also checked the structure of graphene on Ir(111)

with a variety XC functionals and treatments of vdW interactions. Here, we focus on periodicity 10-on-9 (calculations on 1 × 1 structure can be found in the SI). The height of the carbon atoms above the outer-most layer of Ir(111) are shown in Figure 5 together with the results from the experimental derivation of the atomic heights from the LEED-IV, XSW and nc-AFM experiments.36,40 Again several treatments of the XC term were employed,

7

FIG. 4. Difference in electron density ∆n(z) averaged parallel to the surface; PBE+D3 and vdW-DF2-rB86 treatments of the XC were used, in PBE+D3(@vdW-DF2-rB86) the atomic positions were taken from the latter. The vertical positions of the atoms are shown with blue and black circles, respectively. Positive and negative values indicate enhancement and dehancement of electron electron density upon the adsorption of h-BN on the surface

some at the experimental, some at the equilibrium DFT value of the lattice constant. The variation in the heights is less pronounced than in the case of h-BN/Ir(111), due to the overall smaller corrugation of the graphene layer. Still the same approximations lead to too weak binding, and thus too large adsorption heights (including vdW-DF, vdW-DF2 and BEEF-DF2), therefore binding too little in both type of investigated moir´e structures. As the lateral unit cell is here smaller than in hBN/Ir(111), the importance of the proper k point sampling in the Brillouin zone is enhanced; also the calculations with more k points are less expensive. We used QE in convergence studies, by increasing the cut-off energy and k point sampling to 2 × 2. The former does not lead to any noticeable difference in the height of the carbon atoms, whereas increasing the k point sampling leads to a quantitative difference, where some carbon atoms are adsorbed closer to the surface than with only Γ point. Most of the treatments of the exchange-correlation term yield reasonably accurate adsorption heights.

DISCUSSION

With the present results the structure of h-BN on Ir(111) becomes clearer, resolving most of the controversies in the past literature: Two different moir´e structures, one with smaller, one with larger corrugation, can be found with some approximations to the exchange and correlation effects in the DFT calculations. The details of the calculations, whether the consistently optimised DFT-XC lattice constant or 12-on-11 or 13-on-12 lateral cell is used, do not affect the results much. The energy difference ∆E between these two kind

of structures indicates how delicate the geometry is: Sometimes the larger corrugation is preferred – axc PBE+D2 (∆E = -0.13 eV), axc -PBE+D2/13-on-12 (∆E = -0.26 eV), sometimes the smaller corrugation axc PBE+D3 (∆E = +0.13 eV), aexpr -PBE+D3 (∆E = +0.13 eV), aexpr -PBE+D2 (∆E = +0.18 eV). Where the vdW-functionals yield two structures, they are energetically almost degenerate, with the ∆E difference being -40 meV or lower in favour of the larger corrugation. That the ∆E is always so small, at most ≈ 0.2 eV in magnitude in such a large moir´e cell, indicates that the energy balance is indeed very sensitive, and we cannot exclude that even slightest numerical issues, or a finite but low temperature, might turn the balance toward the other kind of corrugation. Overall, there are in both approaches to model the vdW forces in DFT, semi-empirical and nonlocal, treatments that yield good or very good agreement with the experiment and some that miserably fail. That there is such a large scatter in the structures is somewhat worrying. In particular the original vdW-DF and vdWDF2 yield too weak interaction between the h-BN or gr and the substrate, respectively. While this underbinding transfers from one system to the other for certain vdW DFT methods, there is no general trend in the transferability. This is highlighted by the results for vdWDF-optB88, which show good agreement for gr but completely underestimate the corrugation for h-BN. Similarly, while none of the semi-empirical methods performs particularly well for h-BN/Ir(111), PBE+D3 in combination with the DFT lattice constant yields very good agreement in the case of gr/Ir(111). Results closer to the experiments are nowadays obtained mainly by “tuning” the approximation used in the exchange term, even if the vdW interactions originate purely from the correlation, due to the unpredictable error cancellation between the exchange and correlation. This kind of “fitting” to obtain correct structures is an indication of the challenge that the community currently faces while looking for an efficient yet most accurate and “reliable” approximation to use. Overall it is thus not sufficient to state that the calculations were performed using a treatment of the vdW effects, as this is no guarantee of the accuracy of the results. In conclusion, we have performed DFT calculations on h-BN and graphene adsorbed on Ir(111) using different approximations to incorporate the vdW interactions. Overall we find a large variation in the corrugation and distance of the h-BN layer from the substrate with different treatments of the vdW interactions. This supports the “common knowledge” that the choice of the treatment is more “ad hoc” than “predicted”, and the DFT community is in the process of searching for the “best” approach, much like the GGAs have been “fitted” over the years – and like there, the approach to be “preferred” can depend on the kind of system under study. Future progress requires both quantitative experiments and developments in vdW DFT methodology.

8

FIG. 5. Height of C atoms above the average top-most layer of the substrate in the 10-on-9 structure of gr/Ir(111): Top panel, experimental, bottom panel DFT-derived lattice constant; the LEED-I(V ) result by H¨ am¨ al¨ ainen et al. from Ref. 36 is marked with the dashed lines, minimum and maximum. The blue bar on the right indicates the corrugation obtained from the analysis of the AFM data in Ref. 36 and the green ones the extremal values from XSW data in Ref. 38.

METHODS

nc-AFM measurements. Monolayer h-BN on Ir(111) was grown by low-pressure high-temperature chemical vapour deposition under ultra-high vacuum (UHV) conditions (base pressure 10−10 mbar) as described in Ref. 43. nc-AFM measurements were carried out in a Createc LT-STM/AFM housed within the same UHV system. The microscope was equipped with a qPlus tuning fork sensor55 and operated at a temperature of 5 K. The qPlus sensor had a resonance frequency f0 of ∼30.68 kHz, a quality factor Q of ∼98k and a stiffness k of ∼1.8 kN/m. In order to minimize attractive shortrange interactions between the probe tip and the h-BN surface, the tip apex was passivated by deliberate pickup of a carbon monoxide molecule (CO) from a Cu(111)

surface50,51 prior to all measurements. After successful CO pick-up, the Cu(111) sample was exchanged for the hBN/Ir(111) sample.44 All subsequent nc-AFM measurements were acquired in the frequency modulation mode46 using an oscillation amplitude of 50 pm and at a sample bias voltage of 0 V. DFT calculations. We performed total energy calculations using density functional theory (DFT)56 within the Kohn-Sham formalism.57 We have used two codes for the DFT calculation, CP2k (http://www. CP2k.org/) and Quantum ESPRESSO (QE) (http://www. Quantum-ESPRESSO.org/); the choice of two separate codes gave us the possibility to include more approximations in the total energy functional, and verify the results obtained with the two different kinds of numerical implementations. If not otherwise mentioned, the code

9 used was CP2k. The details of the calculation are given in Supporting Material, where also the different approximations, acronyms, and basic bulk and surface properties are listed. In general, we include the vdW interactions to the total energy either in a semi-empirical manner, ie an additional term in the total energy that includes or does not the electron density, or by employing a density functional in the exchange and correlation (XC) term. The XC/vdW treatments employed are listed in Table II.

TABLE II. DFT approaches employed in the present work XC LDA PBE+D2 PBE+D3 PBE+TS revPBE+D2 revPBE+D3 vdW-DF vdW-DF-C09 vdW-DF-Cx vdW-DF-optB86b vdW-DF-optB88 vdW-DF-rB86 vdW-DF2 vdW-DF2-C09 vdW-DF2-rB86 BEEF-DF2 rVV10 PBE+rVV10

Exchange Slater PBE PBE PBE revPBE revPBE revPBE C09 Cx optB86b optB88 rB86 rPW86 C09 rB86 BEEF rPW86 PBE

Correlation PZ PBE PBE PBE PBE PBE LDA+DF1 LDA+DF1 LDA+DF1 LDA+DF1 LDA+DF1 LDA+DF1 LDA+DF2 LDA+DF2 LDA+DF2 scPBE+DF2 PBE+rVV10 PBE+rVV10

vdW References 58 D2 59 and 60 D3 59 and 61 TS 59 and 62 D2 59, 60, and D3 59, 61, and 63 and 64 64 and 65 64 and 66 64 and 67 64 and 68 64 and 69 70 and 71 65 and 71 69 and 71 59, 71, and 59 and 73 59, 70, and

In addition, we tested different lattice constants and investigated how they influence the final geometry: (i) the experimental lattice constant aexp = 3.840 ˚ A,74 (ii) the DFT-optimised value of bulk Ir axc or (iii) the same lattice constant as used in Ref. 54 of aIr = 3.801 ˚ A, obtained by optimising the lattice constant of h-BN, ah−BN = 2.688 ˚ A. We started the optimisation of the structure with axc from a flat layer of h-BN or gr at either 2.0 or 2.7 ˚ A above the outer-most layer of the truncated, four-layer slab of Ir(111), and kept the two lowest layers of Ir fixed at the bulk positions; to start the calculation with aexp we rescaled isotropically the previously obtained geometry to the experimental lattice constant. In some calculations we initially kept the N atom fixed at the distance 2.0 ˚ A, and when the structure was otherwise converged, we relieved this constraint. Some of these calculations led to a structure with a larger corrugation and a smaller height of the pores above the substrate than the calculation where also this N atom was allowed to fully relax from the beginning; thus two different geometries were obtained. Some calculations yielded the same geometry independent of the constraint. In calculations of h-BN we mainly used the moir´e periodicity of 12-on-11, but did some calculations also with 13-on-12 so as to be able to compare with the previous results43 and the dependency of the calculated properties on the periodicity.

∗

†

1

Corresponding Author. [email protected]; http://www.iki.fi/˜apsi/ Present address: IBM Research - Zurich Research Laboratory, S¨ aumerstrasse 4, CH-8803 R¨ uschlikon, Switzerland A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fischer, T. Mueller, K. A. Persson, and G. Ceder, Comp. Mater. Sci. 50, 2295 (2011).

2

3

4

K. Yang, W. Setyawan, S. Wang, M. B. Nardelli, and S. Curtarolo, Nat. Mater. 11, 614 (2012). A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, APL Mater. 1, 011002 (2013). S. Leb`egue, T. Bj¨ orkman, M. Klintenberg, R. M. Nieminen, and O. Eriksson, Phys. Rev. X 3, 031002 (2013).

63 63

72 73

10 5

6

7

8 9

10

11

12

13 14

15

16

17

18

19 20

21

22

23

24

25

26

27 28

29

30

31

32

A. Jain, Y. Shin, and K. A. Persson, Nat. Rev. Mater. 1, 15004 (2016). B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017). N. Mounet, M. Gibertini, P. Schwaller, D. Campi, A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Cepellotti, G. Pizzi, and N. Marzari, Nat. Nano. 13, 246 (2018). J. Behler, J. Chem. Phys. 145, 170901 (2016). M. Todorovi´c, M. U. Gutmann, J. Corander, and P. Rinke, arXiv:1708.09274 . T. L. Jacobsen, M. S. Jørgensen, and B. Hammer, Phys. Rev. Lett. 120, 026102 (2018). V. L. Deringer, C. J. Pickard, and G. Cs´ anyi, Phys. Rev. Lett. 120, 156001 (2018). T. Yamashita, N. Sato, H. Kino, T. Miyake, K. Tsuda, and T. Oguchi, Phys. Rev. Mater. 2, 013803 (2018). A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013). C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, A. V. Kretinin, J. Park, L. A. Ponomarenko, M. I. Katsnelson, Y. N. Gornostyrev, K. Watanabe, T. Taniguchi, C. Casiraghi, H. J. Gao, A. K. Geim, and K. S. Novoselov, Nat. Phys. 10, 451 (2014). K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. Castro Neto, Science 353, aac9439 (2016). M. Yankowitz, K. Watanabe, T. Taniguchi, P. San-Jose, and B. J. LeRoy, Nat. Commun. 7, 13168 (2016). G. Argentero, A. Mittelberger, M. Reza Ahmadpour Monazam, Y. Cao, T. J. Pennycook, C. Mangler, C. Kramberger, J. Kotakoski, A. K. Geim, and J. C. Meyer, Nano Lett. 17, 1409 (2017). J. J. Rehr, E. Zaremba, and W. Kohn, Phys. Rev. B 12, 2062 (1975). P. Pyykk¨ o, Angew. Chem. Intern. Ed. 43, 4412 (2004). I.-C. Lin, A. P. Seitsonen, M. Coutinho-Neto, I. Tavernelli, and U. Roethlisberger, J. Phys. Chem. B 113, 1127 (2009). J. Schmidt, J. VandeVondele, I.-F. W. Kuo, D. Sebastiani, J. I. Siepmann, J. Hutter, and C. J. Mundy, J. Phys. Chem. B 113, 11959 (2009). Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018). Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. JarilloHerrero, Nature 556, 80 (2018). G. Ertl, Reactions at Solid Surfaces (John Wiley & Sons, Inc., Hoboken, 2009). A. Zangwill, Physics at Surfaces (Cambridge University Press, Cambridge, 1988). S. Braun, W. R. Salaneck, and M. Fahlmann, Adv. Mater. 21, 1450 (2009). S. Grimme, WIREs-Comp. Mol. Sci. 1, 211 (2011). J. Klimeˇs and A. Michaelides, J. Chem. Phys. 137, 120901 (2012). K. Berland, V. R. Cooper, K. Lee, E. Schroeder, T. Thonhauser, P. Hyldgaard, and B. I. Lundqvist, Rep. Prog. Phys. 78, 066501 (2015). D. Tunega, T. Buˇcko, and A. Zaoui, J. Chem. Phys. 137, 114105 (2012). H. Yildirim, T. Greber, and A. Kara, J. Phys. Chem. C 117, 20572 (2013). F. Mittendorfer, A. Garhofer, J. Redinger, J. Klimes,

33 34

35

36

37

38

39

40

41

42

43

44

45

46

47 48

49

50

51

52

53

54

55 56 57 58

J. Harl, and G. Kresse, Phys. Rev. B 84, 201401 (2011). I. Hamada and M. Otani, Phys. Rev. B 82, 153412 (2010). W. Moritz, B. Wang, M. L. Bocquet, T. Brugger, T. Greber, J. Wintterlin, and S. G¨ unther, Phys. Rev. Lett. 104, 136102 (2010). C. Busse, P. Lazi´c, R. Djemour, J. Coraux, T. Gerber, N. Atodiresei, V. Caciuc, R. Brako, A. T. N’Diaye, S. Bl¨ ugel, J. Zegenhagen, and T. Michely, Phys. Rev. Lett. 107, 036101 (2011). S. K. H¨ am¨ al¨ ainen, M. P. Boneschanscher, P. H. Jacobse, I. Swart, K. Pussi, W. Moritz, J. Lahtinen, P. Liljeroth, and J. Sainio, Phys. Rev. B 88, 201406 (2013). E. N. Voloshina, E. Fertitta, A. Garhofer, F. Mittendorfer, M. Fonin, A. Thissen, and Y. S. Dedkov, Sci. Rep. 3, 1072 (2013). S. Runte, P. Lazi´c, C. Vo-Van, J. Coraux, J. Zegenhagen, and C. Busse, Phys. Rev. B 89, 155427 (2014). M. Liu, Y. Li, P. Chen, J. Sun, D. Ma, Q. Li, T. Gao, Y. Gao, Z. Cheng, X. Qiu, Y. Fang, Y. Zhang, and Z. Liu, Nano Lett. 14, 6342 (2014). F. H. F. zum Hagen, D. M. Zimmermann, C. C. Silva, C. Schlueter, N. Atodiresei, W. Jolie, A. J. Mart´ınezGalera, D. Dombrowski, U. A. Schr¨ oder, M. Will, P. Lazi´c, V. Caciuc, S. Bl¨ ugel, T.-L. Lee, T. Michely, and C. Busse, ACS Nano 10, 11012 (2016). Z. Sun, S. K. H¨ am¨ al¨ ainen, J. Sainio, J. Lahtinen, D. Vanmaekelbergh, and P. Liljeroth, Phys. Rev. B 83, 081415 (2011). A. T. N’Diaye, J. Coraux, T. N. Plasa, C. Busse, and T. Michely, New J. Phys. 10, 043033 (2008). F. Schulz, R. Drost, S. K. H¨ am¨ al¨ ainen, T. Demonchaux, A. P. Seitsonen, and P. Liljeroth, Phys.Rev. B 89, 235429 (2014). M. P. Boneschanscher, J. van der Lit, Z. Sun, I. Swart, P. Liljeroth, and D. Vanmaekelbergh, ACS Nano 6, 10216 (2012). F. Schulz, S. K. H¨ am¨ al¨ ainen, and P. Liljeroth, “Noncontact atomic force microscopy: Volume 3,” (Springer International Publishing, 2015) Chap. Atomic-scale contrast formation in AFM images on molecular systems, pp. 173– 194. T. R. Albrecht, P. Gr¨ utter, D. Horne, and D. Rugar, J. Appl. Phys. 69, 668 (1991). F. J. Giessibl, Rev. Mod. Phys. 75, 949 (2003). F. Schulz, R. Drost, S. K. H¨ am¨ al¨ ainen, and P. Liljeroth, ACS Nano 7, 11121 (2013). F. Schulz, M. Ij¨ as, R. Drost, S. K. H¨ am¨ al¨ ainen, A. Harju, A. P. Seitsonen, and P. Liljeroth, Nat. Phys. 11, 229 (2015). L. Bartels, G. Meyer, and K.-H. Rieder, Appl. Phys. Lett. 71, 213 (1997). L. Gross, F. Mohn, N. Moll, P. Liljeroth, and G. Meyer, Science 325, 1110 (2009). B. Schuler, W. Liu, A. Tkatchenko, N. Moll, G. Meyer, A. Mistry, D. Fox, and L. Gross, Phys. Rev. Lett. 111, 106103 (2013). F. Schulz, J. Ritala, O. Krejˇc´ı, A. P. Seitsonen, A. S. Foster, and P. Liljeroth, in review. J. G. D´ıaz, Y. Ding, R. Koitz, A. P. Seitsonen, M. Iannuzzi, and J. Hutter, Theor. Chem. Acc. 132, 1350 (2013). F. J. Giessibl, Appl. Phys. Lett. 73, 3956 (1998). P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).

11 59

60 61

62

63 64

65 66

67

68

69 70

71

72

73

74

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). S. Grimme, J. Comp. Chem. 25, 1463 (2004). S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010). A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998). M. Dion, H. Rydberg, E. Schr¨ oder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004). V. R. Cooper, Phys. Rev. B 81, 161104 (2010). K. Berland and P. Hyldgaard, Phys. Rev. B 89, 035412 (2014). J. Klimeˇs, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011). J. Klimeˇs, D. R. Bowler, and A. Michaelides, J. Phys: Cond. Matt. 22, 022201 (2010). I. Hamada, Phys. Rev. B 89, 121103 (2014). E. D. Murray, K. Lee, and D. C. Langreth, J. Chem. Theo. Comp. 5, 2754 (2009). K. Lee, D. E. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, Phys. Rev. B 82, 081101 (2010). J. Wellendorff, K. T. Lundgaard, A. Mogelhoj, V. Petzold, D. D. Landis, J. K. Nørskov, T. Bligard, and K. W. Jacobsen, Phys. Rev. B 85, 235149 (2012). R. Sabatini, T. Gorni, and S. de Gironcoli, Phys. Rev. B 87, 041108 (2013). N. W. Ashcroft and N. D. Mermin, Solid State Physics

(HRW International Editions, 1976).

ACKNOWLEDGEMENTS

This research made use of the Aalto Nanomicroscopy Center (Aalto NMC) facilities and was supported by the Academy of Finland (projects no. 311012 and 314882). APS acknowledges the computational resources at CSC, Espoo, project 2000606, and Centro Svizzero di Calcolo Scientifico (CSCS), Lugano, project uzh11.

AUTHOR CONTRIBUTIONS

The project was initiated by APS, who also carried out all the DFT calculations. FS carried out the nc-AFM experiments. APS and FS wrote the manuscript, with additional input from PL. All authors contributed to the discussion and interpretation of the results.

ADDITIONAL INFORMATION