Interface engineering of graphene nanosheet reinforced ZrB$_2$ composites by tuning surface contacts

TL;DR

This study uses density functional theory to analyze and compare the structure and mechanical behavior of heterophase interfaces in graphene-reinforced ZrB$_2$ composites, highlighting how surface chemistry influences interface properties.

Contribution

It systematically investigates the atomic-scale interface structures and interactions, revealing how surface chemistry and contact types affect mechanical response and bonding in graphene/ZrB$_2$ composites.

Findings

Zr-C-Zr interfaces exhibit strong chemical bonding.

B-C-B interfaces show weak π-π stacking interactions.

Surface chemistry controls interface structure and mechanical response.

Abstract

The mechanical properties of heterophase interfaces are critically important for the behaviour of graphene-reinforced composites. In this work, the structure, adhesion, cleavage and sliding of heterophase interfaces, formed between a ZrB$_2$ matrix and graphene nanosheets, are systematically investigated by density functional theory, and compared to available experimental data. We demonstrate that the surface chemistry of the ZrB$_2$ matrix material largely shapes the interface structures (of either Zr-C-Zr or B-C-B type) and the nature of the interfacial interaction. The Zr-C-Zr interfaces present strong chemical bonding and their response to mechanical stress is significantly influenced by graphene corrugation. In contrast B-C-B interfaces, interacting through the relatively weak $\pi$-$\pi$ stacking, show attributes similar to 2D materials heterostructures. Our theoretical results provide insights into the interface bonding mechanisms in graphene/ceramic composites, and emphasize the prospect for their design via interface engineering enabled by surface contacts.

Tables

Table 1. Table 1: Mechanical parameters of ZrB 2 matrix, graphene and ZrB 2 /graphene interface. W sep subscript 𝑊 sep W_{\mathrm{sep}} : work of separation, σ c subscript 𝜎 c \sigma_{\mathrm{c}} : cleavage strength, δ c subscript 𝛿 c \delta_{\mathrm{c}} : critical separation distance, δ f subscript 𝛿 f \delta_{\mathrm{f}} final separation distance, γ us subscript 𝛾 us \gamma_{\mathrm{us}} : unstable stacking fault energy, γ sf subscript 𝛾 sf \gamma_{\mathrm{sf}} : stacking fault energy, τ c subscript 𝜏 c \tau_{\mathrm{c}} : shear strength.

Mode I loading (opening)		$W_{\mathrm{sep}}$ (J/m2)	${\sigma }_c$ (GPa)	${\delta }_c$ (Å)	${\delta }_f$ (Å)
ZrB2	[10$\overline{1}$0]B-B	10.29	82.16	0.43	2.50
	[10$\overline{1}$0]Zr-B	9.97	56.06	0.62	3.56
	[0001]Zr-B	9.38	49.54	0.66	3.79
Graphene	[0001]C-C	0.36	2.27	0.55	3.17
Interfaces	IZr-C	1.79	15.95	0.39	2.24
	IIIZr-C	1.73	14.67	0.41	2.35
	IIZr-C	1.33	11.34	0.41	2.35
	IB-C	1.20	8.96	0.46	2.68
	IIB-C	0.85	4.02	0.73	4.23
Mode II loading (sliding)			${\gamma }_{\mathrm{sf}}$ (J/m2)	${\gamma }_{\mathrm{us}}$ (J/m2)	${\tau }_c$ (GPa)
ZrB2			-	3.47	49.62
			3.03	3.47	24.45
Graphene	armchair		-	0.09	1.00
	zigzag		-	0.04	0.48
Interfaces	IZr-C		-	0.57	4.55
	IZr-C		0.05	0.57	8.37
	IB-C		0.04	0.10	1.87
	IB-C		0.04	0.15	3.35

Equations

$$E_{\mathrm{ad}}=\frac{E_{\mathrm{tot}}(a_i)-E_{\mathrm{top}}(a_0)-E_{\mathrm{bot}}(a_0)-E_{\mathrm{Gr}}(l_{\mathrm{C}-\mathrm{C}})}{2A}.$$

$$E_{\mathrm{ad}}=-E_{\mathrm{b}}+U_i.$$

$$U_i(a_i)=\frac{[E_{\mathrm{top}}(a_i)-E_{\mathrm{top}}(a_0)]+[E_{\mathrm{bot}}(a_i)-E_{\mathrm{bot}}(a_0)]+[E_{\mathrm{Gr}}(a_i)-E_{\mathrm{Gr}}(l_{\mathrm{C}-\mathrm{C}})]}{2A},$$

$$-E_{\mathrm{b}}(d_i;a_i)=E_{\mathrm{ad}}(d_i;a_i)-U_i(a_i)=\frac{E_{\mathrm{tot}}(d_i;a_i)-E_{\mathrm{top}}(a_i)-E_{\mathrm{bot}}(a_i)-E_{\mathrm{Gr}}(a_i)}{2A}.$$

$$E_{\mathrm{c}}(z)=\frac{E_{\mathrm{tot}}(d_0+z)-E_{\mathrm{tot}}(d_0)}{A}.$$

$$E_{\mathrm{c}}(z)=W_{\mathrm{sep}}\times [1-ex{p}^{-a^{\prime }z}{]}^2,$$

$$\sigma (z)=20a^{\prime }{W}_{\mathrm{sep}}\times [ex{p}^{-a^{\prime }z}-e^{-2a^{\prime }z}],$$

$${\delta }_{\mathrm{f}}=\frac{2\Delta W_{\mathrm{sep}}}{{\sigma }_{\mathrm{c}}}.$$

$$\gamma =(E_{\mathrm{tot}}^{\mathrm{sh}}-E_{\mathrm{tot}}^0)/A,$$

$${\tau }_x=\frac{d\gamma }{dx},$$

Full text

\SectionNumbersOn\altaffiliation

School of Physics and CRANN, Trinity College, Dublin 2, Dublin, Ireland

\altaffiliationSchool of Physics and CRANN, Trinity College, Dublin 2, Dublin, Ireland

{tocentry}

1 Introduction

Within the last ten years, the use of graphene as nanofiller in ceramic matrix composites (CMCs), the so called GCMC materials, has attracted plenty of research interest. They find application in various industry sectors such as aerospace, automotive, energy & power, micro-electronics and pharmaceutical. 1, 2, 3, 4 In addition to the excellent mechanical (a tensile strength of 130 GPa and a Young’s modulus of 1 TPa), electronic and thermal properties, the extremely high specific surface area (2630 m2g*-1*) of graphene provides great capacity for functionalizing and bonding to the surrounding ceramic matrices. 5, 6 For instance, the fracture toughness parameter, $K_{IC}$, can be improved by as much as 235 % for only a 1.5 vol% addition of graphene in a Si3N4 matrix 4. Toughness improvement is also found for zirconium diboride (ZrB2) 7, 8, 9 silicon carbide 10, tantalum carbide 11 and alumina 12. At the same time, the addition of graphene can also suppress the growth of unwanted oxide layers and refine the ceramic grains 7, 13. Last but not least, the GCMCs developed with hierarchical architectures can simultaneously improve the mechanical and functional properties 14, 8, 15.

Among various ceramic materials that can be benefited from graphene-based nanofillers, ZrB2 classified as ultra-high temperature ceramic (UHTC), is one of the most promising structural ceramics for aerospace propulsion systems 16, 17. It exhibits unique combination of high melting point (T${}_{m}\sim$ 3246 $\degree$C), chemical inertness, effective wear and environment resistance. However, the relatively weak fracture toughness and the drop of flexural strength and oxidation resistance at high temperatures 18 are awaiting further improvement. Adding continuous fibers (for enhancing fractural toughness and flexural strength) 19 and nano-particles (such as SiC for improving oxidation resistance) 20 can partially overcome these issues. Very recently, the incorporation of graphene into ZrB2 matrix 7, 8, 9 found great prospect of property enhancement via interfacial impacts. Although it was reported that the interfacial shear strength can be enhanced by $\sim$ 236% and the tensile strength by $\sim$ 96% via coating graphene materials. 21 In most cases, the interfacial characteristics is still vague since the characterization techniques of interfaces are generally in an early development stage 22, 23, 24. The difficulty here comes from the nano-size and morphology of interfaces, the presence and variation of defects along interfaces, the sophisticated interface alignments during tests, as well as the complexity associated to data deconvolution and deviations from physical models.

To crack the nuts related with interfacial mechanics, the atomistic simulation method offers a valid alternative. The technique of molecular dynamics (MD) has been applied to examine NbC/Nb25, ZrB2/ZrC26 and ZrB2(0001)/graphene interfaces 27. For instance, the bonding energies of the ZrB2(0001)/graphene interface have been predicted by using a universal potential function (a purely diagonal harmonic force field). The first-principles calculations were adopted to investigate new mechanical and chemistry phenomena related with interfaces and interphase. 28, 29, 30 The large enhancement of strength, ductivlity and resilience of nano-layered h-BN/silcates was demonstrated by Shahsavarito using the horizontally stacked nanolaminate model.

Here, we exhibited that the interfacial strength of graphene-reinforced ZrB2 nanocompoiste can be largely engineered for more than one order of magnitude by tuning the contact surfaces. This origins in the variation of interfacial bonding mechanism, covalent bonded Zr-C-Zr interfaces or B-C-B interfaces with weak $\pi$-$\pi$ stacking. Also, the corrugation of graphene can further modify the deformation behavior of Zr-C-Zr interfaces, which is different for the processes of interfacial opening and sliding. In comparison, B-C-B interfaces are not that sensitive to the rippling of graphene. We highlighted that the enhancement of interfacial properties of graphene/ZrB2 nanocomposite (a typical example of GCMC) is viable by tuning the chemical environment (leading to a rich variety of Zr-, B- and mix-terminated surfaces) and the interfacial strains (resulting in various extent of graphene ruga). This kind of GCMC materials when properly tailored can be a multifunctional nanocomposites with superior characteristics such as mechanical and electrical properties, thermal and radiation tolerance.

2 Methodology

In this work, we build twelve interface models presenting the tri-layer structures ZrB2/Gr /ZrB2 (Gr = graphene), as schematically shown in Fig. 1(a). Since the direct description of misfit dislocations is beyond the current calculation capability of DFT (supercells with too many atoms are needed), we adopt here the commensurate interface model, in which the ZrB2 slab and the graphene monolayer are constrained to have a common lateral lattice parameter. In brief, three types of interface supercells (I, II and III) are constructed by exposing graphene to the three most stable surfaces of the ZrB2 matrix. These are respectively the Zr- and B-terminated (0001) and the Zr-terminated (10$\bar{1}$0) surfaces. The two ends of such hybrid structures are separated by a vacuum region of 16 Å, in order to prevent the fictitious interaction between the periodic replicas. The two surfaces facing the vacuum regions are all Zr terminated so to reduce the possible effects arising from surface dipoles and ruffling. The structural details of these interface models will be elaborated in section 3.1.

First-principles calculations are performed within the DFT framework using the plane-wave basis projector augmented wave method 31 as implemented in the VASP code 32. The generalised gradient approximation (GGA) parameterised by Perdew, Burke and Ernzerhof (PBE) 33 provides the exchange-correlation energy and potential. In addition, damped van der Waals (vdW) corrections (DFT-D2) 34 are included to account for dispersion interactions. The reliability of the PBE+D2 method in describing transition metal di-borides and graphite has been established before 35, 36.

The Brillouin zone of our interface models are sampled by using the Monkhorst-Pack k-point method, with the following k-meshes, 16$\times$16$\times$1, 14$\times$14$\times$1 and 12$\times$7$\times$1, respectively for supercells I, II and III. The plane-wave kinetic energy cutoff is set to 500 eV. These convergence parameters have all been tested to ensure an energy convergence of 1 meV/atom.

3 Results

3.1 Interface configurations

The ZrB2/Gr/ZrB2 hybrid structures have been constructed with two possible orientations, namely (0001)ZrB2//(0001)Gr and (10$\bar{1}$0)ZrB2//(0001)Gr. When constructing the supercells, we have considered both Zr- and B-terminated ZrB2 facing the graphene layer for the (0001)ZrB2 surface, while only Zr-termination is investigated for (10$\bar{1}$0)ZrB2. As such the interface supercells bind the (0001)Zr, (0001)B and (10$\bar{1}$0)Zr surfaces of ZrB2 to graphene (the subscript here indicates the chemical termination of a given surface). These surfaces have been identified as having high thermodynamic stability in previous DFT calculations 35 and have been frequently found in experiments 37, 38. They are the most likely ones to be exposed to the bonding with graphene nanosheets. In addition, the surfaces of (11$\bar{2}$2)Zr+B with mixed termination of Zr and B and (11$\bar{2}$3)Zr show a relatively high stability when the chemical environment is properly adjusted 35.

Graphene and ZrB2 slabs are joined to form a interface model by using the coincidence lattice method 39. In a nutshell this consists in rotating and straining graphene and ZrB2 supercells so to obtain one supercell with common lattice vectors and little lattice mismatch that, as the same time, with a limited number of atoms. Such exercise has returned us three optimum supercells, where the ZrB2 slabs are as follows: I) $\sqrt{3}\times\sqrt{3}$ (0001); II) $2\times 2$ (0001) and III) $2\times 3$ (10$\bar{1}$0). These define three different interface models, labelled as I, II and III [see figure 1 panels (b), (c) and (d)], with a lattice mismatch of $\Delta a=10$%, $\Delta a=2$% and $\Delta a=3$% ($\Delta b=6$%), respectively ($a$ and $b$ are the in-plane lattice parameters). A second characteristic defining the interfaces is the stacking order at the ZrB2 ends. We denote as AA the situation where the two ZrB2 blocks neighbouring graphene are symmetric, so that Zr atoms of the upper slab is on top of the same kind of atoms on the lower one [see Figs. 1(b), 1(c), 1(d), 1(h) and 1(l)]. In contrast, the stacking is called AB type when the Zr atoms of top ZrB2 slab are in a bridge position with respect to those at the bottom [see 1(e), 1(f), 1(g), 1(j), 1(k) and 1(m)]. Finally AC stacking denotes the situation where the atoms in top ZrB2 slab are in a hollow position [see 1(c)]. The thicknesses of the various interfaces are set between 17 and 26 atomic layers (see the thickness test in the SI of Ref. 35). As already mentioned the two ends of such hybrid structures are separated by a vacuum region of 16 Å, in order to prevent the fictitious interaction between the periodic replicas. It is known that the Zr-terminated surfaces have low surface strains 35. Therefore, the two ZrB2 surfaces facing the vacuum regions are all Zr terminated so to reduce the possible effects arising from surface dipoles and ruffling.

Taking all this into consideration our notation used to describe the hybrid structures is based on 1) the supercell types, 2) the atomic species facing at the interface, and 3) the stacking sequence. For example, II${}_{\mathrm{Zr-C}}^{\mathrm{AA}}$, II${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$ and II${}_{\mathrm{Zr-C}}^{\mathrm{AC}}$ describe the interface models having type II supercell, Zr-C-Zr facing species and the stacking sequences AA, AB and AC, respectively. The top view of their interface regions are illustrated in the panels (c), (f) and (i) of Fig. 1, which contains the same plots for all the interfaces investigated. To be more specific, the atomic structures of our interface models have the following characteristics:

• •

I${}_{\mathrm{Zr-C}}^{\mathrm{AA}}$ and I${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$ [86 atoms with chemical composition (Zr3B6)8C8Zr6] have Zr-C-Zr bonding across both interfaces. They have the same in-plane arrangement as the Ni(111)/graphene interface shown in Fig. 1(b) of Ref. 40. The supercells are made of two fragments of the $\sqrt[]{3}$ $\times$ $\sqrt[]{3}$ (0001) ZrB2 slab, sandwiching a graphene monolayer rotated by 30° with respect to the ZrB2 borophene plane [see Figs. 1(b) and 1(e)];

• •

I${}_{\mathrm{B-C}}^{\mathrm{AA}}$ and I${}_{\mathrm{B-C}}^{\mathrm{AB}}$ [89 atoms with composition (Zr3B6)9C8], displayed in Figs. 1(h) and 1(k), are similar to I${}_{\mathrm{Zr-C}}^{\mathrm{AA}}$ and I${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$, except that the interfaces are B-C-B type;

• •

II${}_{\mathrm{Zr-C}}^{\mathrm{AA}}$, II${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$ and II${}_{\mathrm{Zr-C}}^{\mathrm{AC}}$ [118 atoms with composition (Zr4B8)8C14Zr8], shown in Figs. 1(c), 1(f) and 1(i), are built from two blocks of the 2 $\times$ 2 (0001) ZrB2 slab and one graphene monolayer. Here, the graphene is rotated by 19.1° with respect to the ZrB2 borophene plane. They have the Zr-C-Zr [0001]${}_{\mathrm{ZrB_{2}}}$ stack similarly to I${}_{\mathrm{Zr-C}}^{\mathrm{AA}}$ and I${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$. However, the Zr atoms are misaligned with the C atoms in graphene, resulting in a smaller lattice mismatch (2% against 10%);

• •

II${}_{\mathrm{B-C}}^{\mathrm{AA}}$ and II${}_{\mathrm{B-C}}^{\mathrm{AB}}$ [110 atoms with composition (Zr4B8)8C14] are similar to II${}_{\mathrm{Zr-C}}^{\mathrm{AA}}$ and II${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$, except that the interface structure is B-C-B type [see Figs. 1(l) and 1(m)];

• •

III${}_{\mathrm{Zr-C}}^{\mathrm{AA}}$, III${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$ and III${}_{\mathrm{Zr-C}}^{\mathrm{AB^{\prime}}}$ [148 atoms with composition (Zr6B12)6C28Zr12], shown in Figs. 1(d), 1(g) and 1(j), are constructed from the 2 $\times$ 3 (10$\bar{1}$0) surface slab and one graphene monolayer. They have the Zr-C-Zr [10$\bar{1}$0]${}_{\mathrm{ZrB_{2}}}$ termination, and the Zr atoms are misaligned with respect to the C atoms in graphene.

Finally, the ZrB2/graphene interlayer distance, $d_{i}$, and the common lattice parameters ($a_{i}$ and $b_{i}$) are optimized by searching for the lowest energy points of the $E(a_{i},b_{i},d_{i})$ potential energy surface (see section 3.2.1). The so-calculated parameters are tabulated in Table S1 of the supplementary information (SI), together with the in-plane strains. At such optimized lattice parameters the ionic positions are fully relaxed into their ground state by using the quasi-Newton algorithm to relieve the residue stresses.

3.2 Interface Energetics

3.2.1 Interfacial adhesion energy ($E_{ad}$)

The adhesion energy, $E_{\mathrm{ad}}$, is defined as the energy per unit area released when forming a multi-layered structure from the isolated surface slabs, namely

[TABLE]

Here, Etot(ai), Etop(a0), Ebot(a0) and EGr($l_{\mathrm{C-C}}$) are, respectively, the total energy of the hybrid structure, that of the top and bottom ZrB2 slab and of the graphene monolayer. Note that the reference ZrB2 slabs have an in-plane lattice parameter, $a_{0}$, corresponding to that of their strain-free surface configuration. In contrast the reference graphene monolayer has a C-C bond length of $l_{\mathrm{C-C}}=1.42$ Å. As a consequence of this choice $E_{\mathrm{ad}}$ for the various hybrid structures is computed relatively to the same reference states, namely the strain-free surface slabs. An alternative choice is to take the bulk structures as reference, a choice that will include the surface formation energy into the definition of $E_{\mathrm{ad}}$. Finally, $A$ is the interface area and the pre-factor 2 takes into account the fact that our hybrid structures have two interfaces.

The calculated $E_{\mathrm{ad}}(a_{i};d_{i})$ curves are presented in Fig. 2(a) for all the three kinds of interface supercells I, II and III, while the interlayer distances, $d_{i}$, are kept fixed at 2.5 Å. A parabolic behaviour is observed in all cases, similarly to what recently reported for the interface between graphene and Ti2C MXene 41. The in-plane lattice parameters of the various interfaces are thus determined from the minima of Fig. 2(a). Then the optimal interlayer distances, $d_{0}$, are computed by looking at the minima of the $E_{\mathrm{ad}}(d_{i};a_{i})$ curves taken at the optimized $a_{i}$, see Fig. 2(b).

Several comments can be made by looking at Fig. 2(b). Firstly, we note that the overall $E_{\mathrm{ad}}$ curves move to a lower energy as we go across the series IB-C, IIB-C, IIIZr-C, IZr-C and IIZr-C. Their thermal stability therefore has to be ranked in the reverse order. Secondly, it is quite clear that all the Zr-C interfaces have a deeper $E_{\mathrm{ad}}$ potential well than those with B-C bonding, suggesting that the Zr-C interfaces are energetically more favorable than the B-C ones. Thirdly, we find that the exact stacking order has little effect on the $E_{\mathrm{ad}}$ curves, in particular on their energy minimum, indicating that the local bonding environment plays only a minor role in the interface stability. Finally, one has to note that all the $E_{\mathrm{ad}}$ curves have a long distance tail that asymptotically converges to a positive $E_{\mathrm{ad}}$ value. In particular we have all type I curves converging to $U_{\mathrm{I}}=1.288$ J/m2, the type II to $U_{\mathrm{II}}=0.261$ J/m2 and the type III to $U_{\mathrm{III}}=0.838$ J/m2. Such asymptotic values, $U_{i}$, correspond to the misfit strain energies introduced by imposing a common in-plane lattice parameter. Thus the $U_{i}$’s are proportional to the misfit strains associated to the various interface models. In fact the strain energies are ranked in the order $U_{\mathrm{I}}>U_{\mathrm{III}}>U_{\mathrm{II}}$, which is the same order of the in-plane strains.

Although it is too computationally expensive to include misfit dislocations in our DFT calculations, because of the large supercells required, the misfit strain can be effectively released by considering configurations where graphene presents vertical corrugation, an intrinsic feature of graphene flakes 42. After full structural relaxation, the Zr-C interfaces with AB and AC stacking orders exhibit a more pronounced graphene corrugation than that corresponding to the AA stacking. In contrast, for B-C interfaces more pronounced graphene rippling is present for the AA order. The underlying mechanisms leading to these structural differences will be analyzed later when discussing the electronic structure of the interfaces.

3.2.2 Interfacial binding energy

As discussed before our interface supercells are constructed using the commensurate model so that the heterogeneous layers have common in-plane lattice parameters, a fact that introduces misfit strain. As a consequence, contributions to $E_{\mathrm{ad}}$ originating from the misfit strain energy, $U_{i}$, add to those coming from the formation of the chemical bond at the interface, $E_{\mathrm{b}}$. In order to decouple the two contributions, we assume that $E_{\mathrm{ad}}$ can be written as

[TABLE]

The misfit energy $U_{i}$ can be approximated by the following expression

[TABLE]

where $E_{\mathrm{top}}(a_{i})$, $E_{\mathrm{bot}}(a_{i})$ and $E_{\mathrm{Gr}}(a_{i})$ are, respectively, the total energy of the top and bottom ZrB2 slab and that of the graphene monolayer. The other terms are the same as mentioned in Eq. (1). Hence one has

[TABLE]

The binding energy curves $E_{\mathrm{b}}(d_{i};a_{i})$ are then presented in Fig. 3. In comparison of $E_{\mathrm{ad}}(d_{i})$, they all asymptotically converge to zero. Interestingly all the $E_{\mathrm{b}}(d_{i})$ curves seem to cluster into three main groups: (i) the IZr-C interface has the deepest $E_{\mathrm{b}}$ well ($\sim 2.6$ J/m2), indicating strong interfacial interaction; (ii) the interfaces IIZr-C and IIIZr-C have an intermediate $E_{\mathrm{b}}$ minimum at around 1.7 J/m2, which is 35% lower than that of IZr-C, suggesting a moderate interfacial interaction; (iii) all the interfaces of B-C type show a very shallow $E_{\mathrm{b}}$ profile, with a minimum at around 0.5 J/m2. This latter value is close to that calculated for graphite (0.4 J/m2).

It is worthy to mention that $E_{\mathrm{ad}}$ and $E_{\mathrm{b}}$ for the B-C interfaces are about one order of magnitude smaller than those of the Zr-C interfaces. Such differences are related to the interfacial bonding mechanism. Our guess is that graphene/borophene layers in the B-C interface are coupled by weak physical adsorption, while the Zr-C interfaces are bonded by strong chemical interaction, similar to that at play in metal/graphene contacts 40. These hypotheses will be investigated further in the electronic structures section 4.1.

In concluding this section we remark our main finding, namely that the interface adhesion energy, $E_{\mathrm{ad}}$, is strongly affected by the chemical species facing each other across the interface. This, however, also includes the misfit strain energy, $U_{i}$. When $U_{i}$ is subtracted from $E_{\mathrm{ad}}$ we are left with the interface binding energy, $E_{\mathrm{b}}$, which ranks the stability of our interface models as I${}_{\mathrm{Zr-C}}>\mathrm{II/III}_{\mathrm{Zr-C}}>\mathrm{I/II}_{\mathrm{B-C}}$. It is worthy to mention that, besides the thermodynamically stable Zr-C interfaces, also the B-C ones may be relevant for composites, since B-terminated ZrB2 43 can form in B-rich growth conditions.

3.3 Interfacial cleavage

The interfacial fracture behavior is investigated by computing the traction-displacement curves, which are obtained from the derivative of the total energy changes upon displacing two materials adjacent to the interface 44. Calculations are performed for two different displacement modes, namely opening (mode I, the two sides of the interface are displaced orthogonally to the interface plane) and sliding (mode II, the two sides of the interface are displaced along the interface plane). In this section we will present results for the loading mode I, for which we have adopted the same cleavage model of Lazar and Podloucky 45 for bulk materials, while the displacement mode II (sliding) will be discussed in section 3.4. In the fracture mode I the tri-layer structure is separated by introducing an initial displacement of length $z$ between the top ZrB2 slab and the graphene layer, so that a pre-existing crack is introduced at one side of the interface. Thus, the inter-layer distances of the graphene layer with the bottom and top ZrB2 slabs are, respectively $d_{0}$ and $d_{0}+z$, where $d_{0}$ is the equilibrium interlayer distance calculated before. The corresponding energy change per unit area defines the cleavage energy, $E_{\mathrm{c}}$, which writes

[TABLE]

In Eq. (5) $E_{\mathrm{tot}}(d_{0})$ and $E_{\mathrm{tot}}(d_{0}+z)$ are the total energies of the pristine hybrid structure and of the cleaved one, respectively.

We have then fitted the $E_{\mathrm{c}}(z)$ curves with a Morse function

[TABLE]

where $W_{\mathrm{sep}}$ is the work of separation. Here the parameter $a^{\prime}$ determines the width of the $E_{\mathrm{c}}(z)$ curve, while $2W_{\mathrm{sep}}\times(a^{\prime})^{2}$ controls its curvature at $z=0$. The traction curve, $\sigma(z)$, is calculated as the first derivative of $E_{\mathrm{c}}(z)$ with respect to $z$

[TABLE]

where $\sigma(z)$ is in GPa. Then, the interfacial cleavage strength, $\sigma_{\mathrm{c}}$, and the critical crack length, $\delta_{\mathrm{c}}$, are defined as the values of $\sigma$ and $z$ at the maximum of the $\sigma(z)$ curve. The general behaviour of $\sigma(z)$ is rather simple. As the pre-opening crack grows ($z$ gets larger), $E_{\mathrm{c}}$ continuously increases until it reaches its maximum value, $W_{\mathrm{sep}}$. Thereafter a crack between free (non-interacting) surfaces is formed. The final separation, $\delta_{\mathrm{f}}$, thus can be written as

[TABLE]

The typical mode of cleavage used in the calculation is presented in Fig. 4(a), our calculated cleavage energies are shown in Fig. 4(b) for cleaving across the Zr-C and B-C interfaces, and in Fig. 4(c) for cleaving Zr-B and B-B atomic layers in the ZrB2 matrix, while the corresponding traction-separation curves are displayed in Fig. 5. Note that three kinds of structural relaxation strategies are considered when calculating the relevant cleavage quantities, and these are explained in the SI. We start our analysis with the brittle cleavage due to a sharp fracture surface and then we move our attention to the effects of structural relaxation and of the graphene layer corrugation.

Among all the cleavage modes, cleaving B-B bonds along [10$\bar{1}$0] in bulk ZrB2 has the highest $W_{\mathrm{sep}}$ value of 10.29 J/m2, while cleaving between Zr and B both along [10$\bar{1}$0] and [0001] has rather similar $W_{\mathrm{sep}}$ ($\sim$9.50 J/m2). In comparison, cleaving across the interfaces with graphene requires works of separation one order of magnitude smaller (values in the range 0.70-1.79 J/m2). This confirms that the ZrB2/graphene interfaces as the weak parts can deflect cracks during the fracture of composites. Such feature allows the energy to be released at the ZrB2/graphene interface, so that the structural integrity of the ZrB2/graphene composite can be preserved.

In more detail, both $W_{\mathrm{sep}}$ and $\sigma_{\mathrm{c}}$ systematically decrease when going through the interfaces IZr-C, IIIZr-C, IIZr-C, IB-C and IIB-C, namely they follow the ranking obtained from the binding energy curves of Fig. 3. This correspondence was also observed in the past when studying interfaces in lamellar TiAl alloys 46. Since the crack is initiated by breaking the interfacial Zr-C or B-C bonds, it appears that the bonding mechanism across the interfaces is responsible for the interface strength. In other words, the values of $\sigma_{\mathrm{c}}$, $W_{\mathrm{sep}}$ and $E_{b}$ all reflect the ease of interfacial debonding. Although in real composites additional features, such as defects and impurities, may modify the strength of the interfacial interaction and then affect the way of crack propagation, the ranking calculated here already provides a clear map of the mechanical stability of the various components of a graphene/ZrB2 composite.

Several important parameters, such as the work of separation, $W_{\mathrm{sep}}$, the traction strength, $\sigma_{\mathrm{c}}=\mathrm{max}[\sigma(z)]$, the critical separation, $\delta_{\mathrm{c}}$, and the final separation, $\delta_{f}$, can all be extracted from these curves and are summarized in Table 1. These quantities, together with the shear parameters presented in Table 1, can then be used as inputs in continuum simulation models of cohesive zone to provide a complete understanding of the interface de-bonding of GCMC materials. It is worth mentioning that the $\sigma_{\mathrm{c}}$ varies from 15.95 GPa to 4.02 GPa when going from IZr-C to IIB-C. This provides a chance of tuning the mechanical behavior of interfaces in GCMCs over one order of magnitude.

We have then investigated the effects of the structural relaxation on the energetics of the fracture (a detail discussion is provided in the SI). The atomic positions within a region 4.5 Å vertical to the fracture surfaces are relaxed in order to minimize the internal stresses and the total energy. The results are labelled as IIIZr-C(r) in Figs. 4(b) and 5, while the final geometry is shown in the inset of Fig. 5, where it is evident that relaxation produces a corrugation of the graphene layer. It is experimentally known that graphene sheets may exhibit large corrugations when adsorbed on metal surfaces 42. This is observed here, since the buckling of graphene can efficiently relieve the misfit strains across the interface and provides an energy reduction channel alternative to misfit dislocations. As illustrated in Fig. 4(b) with the open circles, the corrugation of graphene can lower down $W_{\mathrm{sep}}$ by $\sim$ 0.4 J/m2. At the same time it makes the critical cleavage stress going from 15 GPa to 10 GPa. This demonstrates that, in general, graphene buckling has the effect of weakening the interface adhesion. The effect originates from the fact that the interface, namely the surface of graphene in direct contact with the ZrB2 surface, is partially detached so that the contact area is reduced by the corrugation. Note that the rippling period of graphene is constrained by the in-plane dimensions of the supercell, and that the III${}_{\mathrm{Zr-C}}^{\mathrm{AB}}$ supercell (148 atoms) is the largest studied here. One then has to expect that for planar structures with a larger in-plane cell, and consequently smaller lattice misfit and internal stress, the effect of graphene buckling will be in general less pronounced.

3.4 Interface sliding

Interfacial sliding processes are studied in order to extract the traction-separation curves under the loading mode II 44. This is modelled by displacing the top ZrB2 slab along a direction parallel to the interface plane, while monitoring the total energy as a function of the sliding vectors. The sliding profile, $\gamma$, can be defined as the change in total energy with respect to the energy of the undistorted structure as a function of the sliding vector, namely

[TABLE]

where $E^{0}_{\mathrm{tot}}$ and $E^{\mathrm{sh}}_{\mathrm{tot}}$ are the total energies of the undistorted and of the distorted structure, respectively.

The rigid energy landscape is derived by monitoring the energy of the distorted structure without performing any structural relaxation. In addition, we have also calculated the effects of full structural relaxation by using the nudged elastic band (NEB) method 47, 48. In practice we allow atomic relaxation both perpendicular to the gliding plane and in-plane and track the minimum energy path (MEP) while keeping the Burgers vectors fixed. This fully relaxed calculations have been performed only for the interface models IZr-C (86 atoms in the supercell) and IB-C (89 atoms). As for the interface types II and III (118 and 148 atoms, respectively), only the rigid sliding profiles have been studied because of the heavy computational costs associated to the DFT-based NEB method. The shear stress along a given direction $x$ is then calculated as the slope of the energy profile along that direction, namely

[TABLE]

where the maximum $\tau$ is defined as the interfacial shear strength, \tau_{\mathrm{c}}^{x}=\mathrm{max}\big{\{}\tau_{x}\big{\}}.

Five slip systems, namely basal $\langle$a$\rangle$: (0001)\langle$$\bar{1}2$\bar{1}$0$\rangle$, basal $\langle$b$\rangle$: (0001)$\langle$10$\bar{1}$0$\rangle$, prismatic $\langle$a$\rangle$: (10$\bar{1}$0)\langle$$\bar{1}2$\bar{1}$0$\rangle$, prismatic $\langle$c$\rangle$: (10$\bar{1}$0)$\langle$0001$\rangle$ and prismatic $\langle$a+c$\rangle$: (10$\bar{1}$0)$\langle$11$\bar{2}$3$\rangle$, are studied. These all concern interfaces of type I and III with [0001]${}_{\mathrm{ZrB_{2}}}$ and [10$\bar{1}$0]${}_{\mathrm{ZrB_{2}}}$ orientation. The equilibrium geometry of the Zr-C interfaces has an AB stacking sequence, while that of the B-C structure is AA type. The top panels of Fig. 6 displays the energy profiles for rigid sliding along the Zr-C interfaces. In general the energy barriers are observed to be relatively small in magnitude ($<0.15$ J/m2). In addition we find that the sliding curves present a number of maxima and minima, marked in the figure by arrows. These correspond to different stacking configurations of the supercell (also indicated in the figure), which are connected by the sliding process. As the different stacking orders have a rather similar binding energy (see Section 3.2.2), we expect the sliding energy profile to be relatively shallow, as confirmed by our calculations.

Then, the traction curves for mode II are derived as numerical derivative of the corresponding sliding energy profiles and they are presented in the bottom panels of Fig. 6. For the basal sliding, namely the slips along $\langle\bar{1}2\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ and $\langle 10\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ of interface IZr-C, the shear strengths are calculated as 1.47 GPa and 0.74 GPa, respectively. In contrast, those along $\langle\bar{1}2\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$, $\langle 0001\rangle_{\mathrm{ZrB_{2}}}$ and $\langle 11\bar{2}3\rangle_{\mathrm{ZrB_{2}}}$ for interface IIIZr-C are, respectively, 1.15 GPa, 1.13 GPa and 0.66 GPa. Therefore, the slip system {10$\bar{1}$0}$\langle 11\bar{2}3\rangle_{\mathrm{ZrB_{2}}}$ is prone to be activated first. The easy activation of the {10$\bar{1}$0}$\langle 11\bar{2}3\rangle_{\mathrm{ZrB_{2}}}$ slip system was previously reported for bulk ZrB2 at over 700$\degree$C 49. The interfacial sliding along $\langle\bar{1}2\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ of IZr-C shows comparable shear strengths with those of $\{0001\}_{\mathrm{ZrB_{2}}}$ and $\{11\bar{2}0\}_{\mathrm{ZrB_{2}}}$ of IIIZr-C, although the latter two exhibit slightly lower values.

In order to trace the the sliding energy profile along the minimum energy path, we perform full relaxation calculations using the NEB method. This also allows us to extract some key shear parameters. The energy maximum, namely the unstable stacking fault energy, $\gamma_{\mathrm{us}}$, governs the dislocation nucleation at sites of stress concentrations such as at the crack tips. The metastable points (local energy minima) correspond to stable stacking faults with their energies, $\gamma_{\mathrm{sf}}$, determining the dislocation core dissociation, the Peierls stress, the dislocation energy, and the primary slip planes. These shear parameters ($\gamma_{\mathrm{us}}$ and $\gamma_{\mathrm{sf}}$) are summarized in Table 1.

The relaxed sliding energy profiles are then shown in Figs. 7(a) and 7(b) for the interfaces IZr-C and IB-C, respectively. Let us look at the IZr-C case [panel (a)] first. We find that $\gamma_{\mathrm{us}}$ for both the basal $\langle\bar{1}2\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ and $\langle 10\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ shear is 0.57 J/m2, a value that is almost four times larger than those obtained for rigid shearing. This increase is due to the atomic rearrangement, especially the out-of-plane corrugation of the graphene layer, which effectively obstructs the interfacial sliding. In order to understand such dramatic increase we note that the total energy reduction due to atomic relaxation differs depending on the precise stacking order. In the case of Zr-C interface, the AB stacking corresponds to the initial position of its sliding energy landscape (see Fig. 6), which shows a deeper energy well after atomic relaxation since some of the stress is released. In contrast, the AA configuration is at a peak of the $\gamma$ curve and it is relatively stress-free. As a consequence atomic relaxation has little effect on its energetics. Thus, the final result of the relaxation is that of increasing the energy barrier, $\gamma_{\mathrm{us}}=E_{\mathrm{AA}}-E_{\mathrm{AB}}$. This suggests that the rôle of corrugated graphene is twofold. On the one hand, it reduces the interface adhesion strength during interface debonding, on the other hand it increases the interfacial friction during the sliding process.

As for the IB-C interface, we note that the $\gamma_{\mathrm{us}}$ values (0.10-0.15 J/m2) are similar to those calculated without performing atomic relaxation. Now $\gamma_{\mathrm{us}}$ is four times smaller than that of IZr-C. This is the result of the shallow $\gamma$ surface of the IB-C interface, consistent with the ranking given before for the adhesive and binding energies. Since the B-C interfaces are bonded by weak physical adsorption (see later discussion), also their low energy states have the symmetric AA stacking, rather than the non-symmetric AB or AC ones. Graphene corrugation has little effect on the interfacial sliding of B-C interfaces, when compared to the case of Zr-C interfaces. We then conclude that the B-C interfaces are much more favourable for sliding that the more adhesive Zr-C ones.

Finally, we compare the ideal shear strength of Zr-C and B-C interfaces with the ZrB2 matrix. The calculated $\tau_{\mathrm{c}}$ for IZr-C are 4.55 GPa and 8.37 GPa, when sliding respectively along $\vec{b}_{1}=\frac{1}{3}\langle\bar{1}2\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ and $\vec{b}_{2}=\langle 10\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$. The same quantities are reduced to 1.87 GPa and 3.35 GPa for the corresponding IB-C interface and they are increased to 49.62 GPa and 24.45 GPa for bulk ZrB2. Thus the interfacial shear strength of ZrB2/graphene interfaces is at least one order of magnitude lower than those of the ZrB2 matrix and the same observation is valid for the cleavage strength. This suggests that the $\gamma$-surfaces are much less corrugated at a heterophase interface, which indicates relatively ease of interfacial sliding. When examining specific interfaces, it is clear that IB-C is more prone to host deformation pathways, than IZr-C.

4 Discussion

4.1 Interfacial bonding mechanism

One of the main result from the previous sections is that the Zr-C and B-C interfaces present very different strengths of adhesion and binding, as well as the traction curves under the two loading modes investigated. This suggests that different bonding mechanisms are at play in these two interfaces. Here, we analyse the interfacial interaction by looking at the projected density of states (PDOSs) and the projected charge density, as shown in Figs. 8(a) through 8(c).

As expected, the partially filled Zr $d$ bands dominate the PDOS for energies starting at 4 eV below the Fermi level, $E_{\mathrm{F}}$, and they represent the main contribution to the Fermi surface for Zr-C interfaces [see Fig. 8(a) and 8(b)]. In this case a distinctive feature is the appearance of a PDOS reduction at about 1 eV below $E_{\mathrm{F}}$, which separates a region ($E>E_{\mathrm{F}}-1$ eV) with little contribution from either B or C, from another region ($E<E_{\mathrm{F}}-1$ eV), where the PDOS from the C 2$p_{z}$ orbital is significant. This is a clear signature of the strong hybridization between the C 2$p_{z}$ and the Zr 4$d$ orbitals, and the consequent formation of Zr-C covalent bonds. Such attribution is consistent with the strong affinity between Zr and C and the existence of covalent ZrC1-x compounds. Covalent bonding across interfaces has also been suggested for the interfaces of Ti2C/graphene 41, ZrC/SiC 50 and metal graphene contacts 40.

Across the Zr-C interfaces the electron-acceptor B ions (with valence configuration $s^{2}p^{1}$) are replaced by C ($s^{2}p^{2}$). This effectively shifts the Fermi level upwards in energy, and in fact the PDOS pseudogap moves from $E_{\mathrm{F}}$ in bulk ZrB2 to 0.8 eV below $E_{\mathrm{F}}$ at the Zr-C interface. The effective charges, as calculated from Bader analysis, of Zr and B in bulk ZrB2 are +1.54$e$ and -0.77$e$ ($e$ is the electron charge) 51. These become +1.23$e$ (Zr) and -0.35$e$ (C) at the Zr-C interface, suggesting that the Zr-C bond has a mixture of ionic and covalent nature.

The situation for the B-C interfaces is rather different, as one can easily conclude by looking at Fig. 8(c). In this case the most prominent feature of the DOS is a peak just below the Fermi level, with significant projections originating from the B and C 2$p_{z}$ orbitals. This suggests that the bonding has a $\pi$-$\pi$ stacking nature, a fact that is confirmed by the projected charge density shown in the inset of Fig. 8(c). Such $\pi$-$\pi$ interaction is commonly present in aromatic compounds 52.

As expected from elementary chemistry the B-C $\pi$-$\pi$ interaction is weaker than the covalent bond between Zr and C, as shown by our previous calculations of the mechanical properties of the various interfaces. However, such $\pi$-$\pi$ interaction is still stronger than the van der Waals coupling between graphene layers in graphite and, in general, among monolayers of van der Waals layered compounds. This is demonstrated in Fig. 9(a), where we compare the binding energy curves, $E_{\mathrm{b}}(d_{i})$, of the B-C interfaces with that of graphite. The binding energies are calculated to be 0.488 J/m2 and 0.663 J/m2 for IB-C and IIB-C, respectively. In contrast we compute $E_{\mathrm{b}}(d_{0})=0.374$ J/m2 for a 6-layer graphene nanosheet. In addition the equilibrium interlayer distance at the B-C interfaces is between 2.5 Å and 3.0 Å, while that in graphite is 3.23 Å. The $\pi$-$\pi$ stacking interaction here is made stronger than that in graphite by two features of the electronic structure:

1. the electrostatic interaction between the negatively charged C and the positively charged B atoms ;
2. the offset stacking of the borophene and graphene layers contributing to reduce the repulsive interaction. Since the van der Waals interaction among the graphene sheets is believed to be the reason behind its agglomeration in composites 53, the stronger adhesion between boraphene and graphene can be exploited as a tool for breaking the agglomeration and promoting a more uniform dispersion of graphene nanofillers in ceramic matrices. This can be beneficial to the improvements of the overall structural stability and functionality of the nanocomposites.

Interestingly, the $\pi$-$\pi$ bonding mechanism across the B-C interfaces gives to these interfaces mechanical properties analogous to those of 2D materials. This is illustrated in Fig. 9(b), where we compare the $\gamma$ sliding curves of the two IB-C interfaces with those calculated for graphene. We find that the $\gamma_{\mathrm{us}}$ values for graphene are 0.04 J/m2 and 0.09 J/m2, respectively for shearing along the zigzag and armchair edges. The slipping along the basal $\langle\bar{1}2\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ and $\langle 10\bar{1}0\rangle_{\mathrm{ZrB_{2}}}$ directions returns the $\gamma_{\mathrm{us}}$ values of 0.10 J/m2 and 0.15 J/m2, respectively, which are comparable to those of graphene.

4.2 Interfacial mechanics

In the previous sections we have observed a large difference between the Zr-C and the B-C interface with respect to their adhesive and mechanical properties. In particular the Zr-C interfaces display stronger adhesion, cleavage and shear strength than those of B-C interfaces. This can be well explained by the difference in interfacial bonding mechanism, which is covalency-dominated in Zr-C and $\pi$-$\pi$-stacking type in B-C interfaces. The stronger Zr-C interfaces thus present high resistance to both cracks and interfacial sliding, while the weak B-C interfaces are much more easy to allow interfacial debonding and shear.

The ZrB2 matrices, have a rich variety of stable surfaces, either Zr-, B- or mix-terminated surfaces with different orientations. They can be somehow engineered by acting on the ZrB2 growth conditions. Note the existence of the Zr/B-terminated (0001) surface has been confirmed by both experiments 43, 37, 38 and theory 54, 35. The (0001) surfaces with Zr/B-terminations emphasized here are the two extreme cases, which exhibit completely different mechanical behavior. In general, Zr-C interfaces can be prepared by stabilising the Zr-terminated surfaces in Zr-rich conditions. While a B-rich environment will favour the formation of B-C interfaces. Therefore, surface chemistry can be used as an efficient tool for tailoring the interfacial mechanics by tuning the surface contacts.

It is important to note that the toughening effect of graphene fillers is correlated with the interface debonding mechanism. Since in composites the interfacial debonding will favour the dissipation of energy via crack bridging, crack deflection and fillers pulling-out, the predominance of a mechanism over the others will affect the crack propagation behavior along various interfaces, so that different ZrB2 growth conditions may result in tuning possibility of the mechanical properties of GCMCs. This observation also provides hints on why there is a wide variation of the fracture-toughness-parameter values in experiments on graphene-reinforced composites 55.

Last but not least, monolithic ZrB2 and graphene nanoplatelets are used as starting materials during spark plasma sintering 7 of ZrB2/graphene composites. The nanoplatelets are made of short stacks of ribbon-shaped graphene sheets. These are functionalised with groups like ethers, carboxyls or hydroxyls, which may further modify the interfacial bonding 56. Therefore, a scrupulously designed surface treatment is vital for the accurate control of the interfacial interaction, and hence of the interface mechanics and the failure mode. Our theoretical results suggest that the preferred interface mechanics can be engineered by acting on the interfacial interactions through the tuning of the surface chemistry of ZrB2 and graphene.

5 Conclusion

In this work we have studied the heterophase interfaces of graphene-reinforced ZrB2 composites based on DFT calculations. A number of atomistic models for the various possible interfaces has been constructed by using the most thermodynamically stable surfaces of the ZrB2 matrix. Then, we have systematically investigated their interface adhesion, mechanical behaviour and bonding mechanism. We have demonstrated that two kinds of interfaces, namely Zr-C-Zr and B-C-B, offer a wide spectrum of mechanical properties due to their dissimilar interaction between the constituent materials. In particular Zr-terminated surfaces bind to graphene in a covalent way, while the interaction with the borophene planes in ZrB2 has a $\pi$-$\pi$-stacking nature. By tuning the surface chemistry of the ZrB2 matrix one can prepare composites that expose graphene predominantly to a specific ZrB2 termination (either Zr or B or mixed), so it is possible to go from a regime of weak graphene/matrix interaction to one where the interaction is relatively strong. This provides an important design scheme in the synthesis of ceramic composites. Furthermore, the fact that the borophene/graphene interface is more adhesive than the graphene/graphene one may offer an opportunity to tackle flakes agglomeration, a issue encountered when processing most graphene-reinforced composites. Finally, we have analysed the rôle of graphene rippling over the mechanical properties of the interfaces. Interestingly we have found that rippling drastically increases the friction for sliding graphene over ZrB2 in the case of Zr termination, while it has no significant effect for that with B termination.

{acknowledgement}

This work is supported by the European Union’s Horizon 2020 “Research and innovation programme” under the grant agreement No. 685594 (C3HARME). Computational resources have been provided by the Irish Center for High-End Computing (ICHEC) and the Trinity Centre for High Performance Computing (TCHPC). The authors also thank Rui Dong for the clattice code used to merge surface slabs.

{suppinfo}

The Supporting Information is available free of charge on the ACS Publications website for: Geometry details of interface models, Structural relaxation effects on interfacial cleavage.