Effect of Graphene on Nickel Surface Relaxation: Molecular Dynamics Simulation

Abstract

The effect of graphene (GR) on Ni surface relaxation and reconstruction in three different substrate orientations, {111}, {001}, and {011}, at two different temperatures, 300 K and 400 K, was studied using molecular dynamics simulation. The change in the interplanar distances of the substrate and redistribution of Ni and C atoms in a direction perpendicular to the surface was compared with the equilibrium state of GR and bulk Ni, in the absence of the counterpart. The surface reconstruction for the GR/Ni system was analyzed based on the calculated radial pair distribution functions of Ni and C atoms. The surface roughness was visualized using 2D atomic distribution maps. The introduction of GR on the Ni surface in any crystallographic orientation decreases the maximum modification of interplanar spacing compared to the bulk by less than 1%. For the studied substrate orientations and temperatures, it was found that the most densely packed {111} orientation of the Ni base provides minimal changes in the structural parameters of both counterparts at 400 K. Additionally, the system formed by GR deposition on Ni {111} at 400 K is characterized by the least roughness.

1. Introduction

Graphene (GR) is a favorable material for various applications due to the combination of its unique properties. In particular, high thermal (up to 5300 W/mK) [1,2] and electrical conductivity (~10⁶ S/m) [3], as well as high (97.7%) transmittance for a broad spectrum of solar radiation [4], makes GR application feasible for translucent and flexible electrodes that can be used in solar cell technologies [5]. Ideally, the solar cell electrode should demonstrate high transparency, low electrical resistance, chemical stability, and low surface roughness and provide all these at a reasonable cost [6]. Thus, intensive research to develop materials that meet these requirements is being carried out on carbon nanotubes [7], conductive polymers [8], metallic nanowires [9], metallic meshes [10], spintronic devices [11], etc. However, GR has the best chance of outperforming in such applications. It has recently been shown that due to the outstanding chemical stability of GR, its introduction into the perovskite solar cell can improve the functional characteristics and stability of the device [12]. In addition, due to the high conductivity of GR, its addition to the front electrodes of commercially used solar cells on the base of indium-tin-oxide [13] and fluorine-doped-tin-oxide [14] provides significant enhancement of the device performance. Numerous studies have recently been published on the incorporation of GR in solar cells. For instance, a semitransparent perovskite solar cell with a multilayer GR applied to the top electrode has been proposed [15]. GR monolayers can also be incorporated into flexible perovskite solar cells [16]. Furthermore, the combination of GR and transversely stacked carbon nanotubes may be used as the front and rear electrodes of a flexible perovskite solar cell [17].

However, GR synthesis or deposition on top of a metal substrate remains a challenging task, especially for some technology-relevant metals [18]. In this respect, GR surface morphology is mostly defined by the ratio of adsorption energy to the bulk cohesive energy of the substrate. When this ratio is lower, larger uniform GR islets can be formed on the metal surface [19]. Thus, a high-quality GR layer is generally easier to obtain on the Cu substrate due to the extremely low solubility of carbon in Cu (~7 × 10⁻⁶% at 1200 K) [20]. This limited solubility causes very low adsorption energy, which provides a relatively large GR island growth and ultimately allows a monolayer of GR to be produced and separated from the substrate [21]. However, a lattice mismatch of ~10% between GR and Cu leads to some GR surface deformation and buckling [22]. It has been reported [23] that the smaller lattice mismatch may be achieved by using the substrate of a Cu/Ni alloy for large-scale GR synthesis. Other substrate metals such as Au [24] and Ag [25] also demonstrate low adsorption and restrained mutual solubility with carbon, making them promising for single-layer GR production. Yet, these noble metals are of high cost, which limits their technological relevance. In turn, Ni provides strong interaction [26] and good lattice matching [27] with GR. However, GR deposited on top of Ni can lead to specific relaxation and structure modifications of the near-surface layers of the substrate [28]. This inevitably modifies the functional properties of both materials and ultimately the solar cell device, which incorporates them. The free surface of a metal is generally accepted as a 2D crystal structure defect where the structure of the near-surface regions can be significantly different compared to the bulk. Application of GR to the free surface of the metal can significantly modify its 2D defect properties, which can be classified as surface relaxation and surface reconstruction.

In the present paper, we report the results of a molecular dynamics simulation study of the surface structure modification of Ni induced by GR. Simulation was performed for the three most closely packed crystal orientations of a Ni substrate {111}, {001}, and {011} at two different temperatures that cover the operation temperature range of a solar cell: 300 K and 400 K. These crystal orientations were chosen due to the fact that for orientations with lower packing density, the reconstruction of surface layers is similar to one of the mentioned above closely packed orientations.

2. Simulation Procedure

The classical molecular dynamics approach was used in this study for the simulation of Ni surface modification by numerically solving Newton’s motion equations for a system of interacting particles. Forces between the particles and their potential energies were calculated using interatomic potentials for the Ni-C system [29]. The LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) software was used to conduct the simulation [30]. Equilibrium lattice parameters and coefficients of linear thermal expansion of Ni and GR were calculated to check the accuracy of applied interatomic potential. This simulation was conducted for a defect-free atomic cell, using periodic boundary conditions in all directions. Equilibrium values of the lattice constants were determined as parameters corresponding to zero stresses of the system, while thermal expansion coefficients were calculated taking into account the average distance between neighboring atoms, not considering changes in the lattice parameters. The selected potential is considered to be appropriate in this study, which is confirmed by less than 0.4% discrepancy between simulated and reference data shown in

Table 1. Structural parameters and coefficients of linear thermal expansion (CLTE) of Ni and GR calculated in this study using molecular dynamics simulation and corresponding reference data.

Material	Temperature, K	Lattice Constant, Å	CLTE, 10−6 1/K
Ni	300	3.5347 (this work) 3.523 [31]	12.75 (this work) 13.4 [32] 12.8 [33]
	400	3.5392
GR	300	1.4163 (this work) 1.42 [34]	8.5 (this work)
	400	1.4175

.

Simulation calibration and system stress analysis determine that the free surface introduces a stress source in Ni that gradually disappears at a distance of five lattice periods. Therefore, the simulation cell contains a GR monolayer, six lattice periods of an “active” Ni, and four lattice periods of a “frozen” Ni (Figure 1). The “active” Ni layer consists of atoms that can move during the simulation process, while the “frozen” Ni layer is composed of atoms that are strictly fixed in their equilibrium positions. The “frozen” area was introduced in order to mitigate the influence of the free surface at the bottom of the simulation cell. The boundaries of the cell are fixed along the Z-axis without periodic conditions, providing free surfaces at the bottom and top of the cell. Born–von Karman boundary conditions were applied along the X- and Y-axes. The size of this area was selected to be higher compared to the cut-off radius corresponding to the interatomic potential influence range. Generally, a simulation cell was chosen with the highest possible crystal structure matching between GR and Ni.

The simulation time step of 10⁻¹⁵ s was selected to make sure that force and atomic positions do not change too rapidly, which was done in order to negate system instability during the simulation process due to atoms overshooting in the direction of applied force.

3. Results

3.1. Surface Relaxation of Ni of Various Orientations with and without GR Layer

The equilibrium position of an individual atom in an ideal infinite crystal is determined by the acting forces from all other atoms in the crystal, resulting in a periodic structure. Laminating the surface with another material alters these forces, affecting the equilibrium positions of the remaining atoms. This changes the spacing and/or symmetry of the “active” layer of atoms compared to the bulk and forms a different surface structure. These modifications in the equilibrium positions of near-surface atoms can be categorized as either a relaxation or a reconstruction. Surface relaxation is characterized by the change in interplanar spacing Δd.

The simulation results of Ni relaxation in various crystal orientations at 300 K and 400 K are given in the form of relative interplanar spacing change compared to the equilibrium bulk values (Figure 2). The selected crystal orientations of Ni were evaluated by the packing density mismatch with GR and are presented in Figure 2 in the order from smallest to largest (top-to-bottom): {111}, {001}, {011}. It is important to note that interplanar spacing modification of less than 0.5% (Figure 2) could be considered as a temporary result of atom thermal fluctuations and were not taken into account for the surface relaxation analysis.

The {111} orientation provides the most densely packed surface plane in the face-centered cubic lattice. The packing density of this plane for Ni under normal conditions of room temperature and atmospheric pressure is 1.864 × 10⁻¹⁹ atoms/m² [35]. For ease of analysis, the packing density in all other orientations was correlated with respect to the {111} orientation. The corresponding values were 1, 0.866, and 0.612 for the {111}, {001}, and {011} planes, respectively.

Ni {111}. The calculated equilibrium values of the interplanar spacing for the Ni {111} orientation in the bulk are d_b³⁰⁰ = 2.041 ± 0.001 Å and d_b⁴⁰⁰ = 2.043 ± 0.001 Å at 300 K and 400 K, respectively. These are consistent with previously reported data [26]. The interplanar (d_ij) distances calculated using the molecular dynamics approach were compared with the d_b one: Δd_ij = (d_b − d_ij)/d_b. The i and j indexes in this notation correspond to the number of the atomic layers/planes counting from the surface; for example, d₁₂ is a spacing between the first and second atomic planes (Figure 3).

The largest absolute interplanar spacing modifications were observed between the first and second planes: d₁₂³⁰⁰ = 2.014 ± 0.001 Å and d₁₂⁴⁰⁰ = 2.016 ± 0.001 Å at 300 K and 400 K, respectively. For the first interplanar spacing (d₁₂), a modification of ~ −1.3% was observed for Ni without GR at both temperatures, as can be seen in Figure 2. The negative value of the spacing modification corresponds to the lattice compression. At 300 K, the relative relaxation Δd₁₂ changes from −1.29% to −0.13%, and Δd₂₃ changes from 0.38% to −0.59% in the presence of GR on top of Ni. At 400 K, the presence of GR almost completely eliminates the influence of the free surface, making interplanar spacing changes not exceeding 0.1%. As a consequence, it can be assumed that on Ni {111}, GR behaves like a missing Ni surface layer, negating overall surface relaxation.

Ni {001}. The calculated equilibrium values of interplanar spacing in the bulk in this orientation are d_b³⁰⁰ = 1.767 ± 0.001 Å and d_b⁴⁰⁰ = 1.769 ± 0.001 Å at 300 K and 400 K, correspondingly, which is also consistent with the previously reported data [36]. At 300 K, modification (Δd = −1.61%) is observed only for the first interplanar spacing from the Ni surface (Figure 3), giving an absolute calculated value of d₁₂³⁰⁰ = 1.7388 Å. An increase in the temperature up to 400 K leads to even more pronounced lattice compression (Δd = −1.99%). Compression of −2.7% was previously reported for this surface orientation, where it was experimentally measured using slow electron diffraction [26].

The presence of GR on top of Ni {001} leads to significant changes in the surface relaxation: the modification of d₁₂³⁰⁰ and d₁₂⁴⁰⁰ interplanar spacings decreases down to −0.78% and −0.85%, respectively.

Ni {011}. Compared to the {111} and {001} orientations, {011} has the lowest packing density, making bulk interplanar spacings of d_b³⁰⁰ = 1.250 ± 0.001 Å and d_b⁴⁰⁰ = 1.251 ± 0.001 Å at 300 K and 400 K, correspondingly. Significant modification of the interplanar distance between the first two atomic planes counting from the surface was registered for this Ni orientation as well. The corresponding calculated values are d₁₂³⁰⁰ = −5.4% and d₁₂⁴⁰⁰ = −4.4%. Significant compression (−4.8 ± 1.7%) between the first two atomic layers of Ni in this direction was also measured experimentally at room temperature and previously reported [37].

The addition of GR slightly changes the surface relaxation. The modification of the first interplanar spacing decreases down to −4.1%, again suggesting compression. In contrast, in the case of the second distance d₂₃, it increases to 1.48% and 2% at 300 K and 400 K, respectively, meaning expansion.

As already mentioned, the maximum value of interplanar spacing modification was registered for all orientations in question for the d₁₂, which is the first spacing from the surface. Relaxation has no impact on more than two first atomic layers of Ni. In addition, it was found that the maximum modification depends on the packing density of the plane, ranging from −1.3% for Ni {111} (the highest packing density) to −5.4% for Ni {011} (the lowest packing density) (Figure 4).

Additionally, it was found that the temperature increase from 300 K to 400 K does not significantly affect the surface relaxation.

Finally, it was shown that GR addition lowers the influence of the free surface regardless of the temperature: the maximum value of interplanar spacing modification decreases by approximately 1% as a result of GR presence. The Ni {111} plane laminated by GR reveals the smallest deviations from the bulk structure and has the smallest surface modification due to the best lattice matching, especially at 400 K.

3.2. Surface Reconstruction

Surface reconstruction is a process that changes the position of surface atoms from their equilibrium state, typical of the bulk, which is accompanied by the formation of a structure that differs from the bulk in periodicity and/or type of symmetry. The presence of GR on the metal surface can significantly modify reconstruction by altering both counterparts at the near-surface region. To analyze the result of such a change, the radial pair distribution function (RPDF), the rearrangement of atoms in a direction normal to the surface, and the 2D maps of atom reconstruction were plotted and evaluated.

Simulation results suggest that there is no surface reconstruction of Ni in the {001}, {011}, and {111} orientations before the introduction of GR. In this case, only relaxation processes were registered. However, Ni reconstruction of the upper first atomic plane and GR occurs in addition to the relaxation process, in the presence of GR. Details of these processes are given below.

3.2.1. GR Reconstruction

Structural changes of GR after its deposition on Ni in different orientations were compared to its standing-alone equilibrium state. Thus, GR structure simulation was first performed without Ni at 300 K and 400 K.

The RPDF, g(r), describes the probability of finding two atoms at distance r from each other. In this event, the first peak of the g(r) plot characterizes a short-range ordering, and its position corresponds to an average distance between the nearest atoms. The intensity of the second and other RPDF peaks allows for tracking the long-range ordering, while a change in these peaks’ positions represents structure deformation of the top plane in the considered orientation.

Results of the GR structure simulation suggest that distance r between C atoms is equal to r_(C-C) = 1.4179 ± 0.0005 Å (Figure 5a). This distance increases by 0.5–0.6% after GR deposition onto Ni in the {001}, {011}, and {111} orientations (Figure 5b–d). The second and the following RPDF peak analysis implies GR expansion of more than 0.2% in the case of the {111} orientation and less than 0.1% in the case of {001} and {011}.

Thus, the average distance between C atoms in the deposited GR changes is 3–5 times higher compared to the expansion of the GR layer itself. This effect can be explained by the increased displacement of C atoms in a direction perpendicular to the surface. The range of C atom positions in this direction is characterized by the ∆l_C parameter. Analysis of atom distribution in the direction perpendicular to the surface allowed us to determine the values: ∆l_C^{(Ni {111})} = 1.4 ± 0.05 Å, ∆l_C^{(Ni {001})} = 0.75 ± 0.05 Å, and ∆l_C^{(Ni {011})} = 0.9 ± 0.05 Å for Ni orientations {111}, {001}, and {011}, respectively. These values are 2–3 times larger compared to the GR in an equilibrium state. C atoms occupy some preferred positions with different distances to the Ni surface, which are schematically illustrated in Figure 3 as l₁₁ and l₁₂. This projection is confirmed by the character of the C peaks in Figure 6. For instance, the presence of the shoulders in Figure 6a,b and in Figure 6e,f and the doubling of peaks in Figure 6d are clear evidence for making such an account.

A temperature increase from 300 K to 400 K does not have a pronounced effect on the GR reconstruction.

Such inhomogeneous distribution of C atoms in a direction perpendicular to the surface causes a specific surface relief. The character of this relief is determined by the orientation of the Ni substrate used for GR application. This relief can be shown in 2D maps reflecting the atomic displacement perpendicular to the surface, where the displacement is color-coded (Figure 7).

If the crystal symmetry of the uppermost Ni layer and GR matches well, the 2D maps reveal clear periodicity, i.e., the frequency of the areas limited by the same height of atom positions (Figure 7a,b). The extent of this periodicity is about 20–25 atomic periods of Ni. If there is no symmetry match, the periodicity is less pronounced and reflected by the strips with a width of 4–5 Ni periods (Figure 7c–f).

The surface roughness of GR laminated systems can affect numerous functional properties [38,39]. That is why the following parameters were analyzed here: R_a—an arithmetic mean deviation of the C atom positions from the average height value; R_max—a maximum relief height, e.g., the difference between the maximum and minimum height of atom positions; S—an average distance between relief irregularities. It is worth noting that the integral surface roughness depends on both the height deviation and distance parameters of the surface relief; for instance, the lowest roughness corresponds to the minimal R_a and R_max values and maximal S. Therefore, the integral S/R_a ratio was used for comparative surface relief analysis; a higher S/R_a value corresponds to lower surface roughness.

The calculated surface roughness parameters of GR/Ni systems for different substrate orientations and temperatures are summarized in

Table 2. Surface roughness parameters of the GR/Ni systems depending on the Ni orientation and temperature.

Ni Plane	Temperature
	300 K	400 K
{111}	Ra = 0.7 ± 0.05 Å	Ra = 0.55 ± 0.05 Å
	Rmax = 1.95 ± 0.05 Å	Rmax = 1.75 ± 0.05 Å
	S = 146 ± 1.4 Å	S = 150 ± 1.4 Å
	S/Ra = 208	S/Ra = 272
{001}	Ra = 0.35±0.05 Å	Ra = 0.4 ± 0.05 Å
	Rmax =1.1 ± 0.05 Å	Rmax = 1.15 ± 0.05 Å
	S = 12 ± 1.4 Å	S = 13 ± 1.4 Å
	S/Ra = 34	S/Ra = 32
{011}	Ra = 0.45 ± 0.05 Å	Ra = 0.55 ± 0.05 Å
	Rmax = 1.2 ± 0.05 Å	Rmax =1.35 ± 0.05 Å
	S = 18 ± 1.4 Å	S = 19 ± 1.4 Å
	S/Ra = 40	S/Ra = 34

. The Ni {111} laminated with GR is characterized by maximum height difference (R_a and R_max parameters) and distance between relief irregularities (S parameter) compared to the other two orientations.

Figure 8 shows the dependence of the S/R_a parameter on the relative packing density of the Ni. At 300 K, the lowest roughness was found for the Ni {111} orientation, showing the maximum packing density of the surface layer. A temperature increase up to 400 K leads to the roughness decrease of the GR/Ni {111} system, giving the highest S/R_a ratio among all cases. This effect is not observed for the other two studied orientations.

3.2.2. Ni Reconstruction

Surface reconstruction of Ni is only observed in the presence of GR and only for the first atomic layer from the Ni surface. The RPDF provides summarized data on the structure of this Ni layer (Figure 9). In turn, structural changes in the normal direction to the surface can be tailed from Figure 6.

As follows from the analysis of the RPDF first peak position (Figure 9), an average distance between Ni atoms decreases by 0.5%, 0.6%, and 1% for the {111}, {001}, and {011} substrate orientations, respectively. According to other RPDF peak positions, the structure of the whole reconstructed Ni layer compacts not more than 0.5% after GR deposition onto the {111} plane. There are no relative changes for the other planes (registered changes are less than 0.1%). It should also be noted that the slight shift of Ni atoms toward the nearest C atoms shown with a dotted line in Figure 9 is taking place in the presence of GR.

The character of atomic distribution in the direction perpendicular to the surface (Figure 6) suggests that after lamination with GR, Ni atoms change their positions in the first atomic plane increasing the thickness of this layer by 79%, 62%, and 42% for the {111}, {001}, and {011} substrate orientations, respectively (Figure 10). Thus, an increase in the plane relative packing density leads to an increase in the width of the first layer from the surface. This can be explained by the fact that for densely packed crystalline planes, atoms have more freedom to displace from their regular positions in an out-of-plane direction (perpendicular to the surface) than in-plane. Surface relief formed at various substrate orientations is shown in Figure 6a,c,e.

A temperature increase from 300 K up to 400 K for the GR/Ni system leads to a lower variation of the thickness of the Ni upper layer and a smaller deviation of the average distance between its atoms (Figure 10). The change in the Ni upper layer thickness is lowered by ~30% for 400 K compared to 300 K for all studied substrate orientations. Modification of the average distance between Ni atoms is also lowered by 0.2–0.4% by the absolute value. This implies a higher stability of the GR/Ni system at a higher temperature. Most likely, such an improvement is a result of the thermal expansion effect, since GR has a lower coefficient of thermal expansion compared to Ni (Table 1). The Ni {111} orientation provides almost perfect matching between GR and metal lattice parameters at 400 K. Therefore, such conditions can be considered perfect for lamination with GR.

4. Discussion

A few words can be said about the implication of the simulation results presented here. It appears that the Ni surface in the {111} orientation provides the best conditions for GR deposition with the temperature of substrate close to 400 K. The structure of the substrate in this orientation is most close to the GR. As can be seen, the surface structure of Ni is restored after GR deposition in all studied cases; that is, the impact of the free surface of Ni decreases in the presence of GR. It suggests that for the GR/Ni system, the properties of Ni near-surface layers are closer to its bulk. Another implication of this simulation is the prediction of the surface roughness. It is well accepted that the oxidation of the metallic surface strongly depends on its morphology [40]. In turn, GR shows great potential for the corrosion resistance of metals (Me) [41]. Thus, it is very important to pay attention to the link between the structure of the GR/Me interface and corrosion resistance [12]. Various structure defects such as wrinkles, folds, and atomic defects lead to an interface roughness increase with the following insertion of oxygen and water molecules from the environment. Thus, GR deformation, warpage, and even the enfoldment of graphene sheets will affect corrosion properties. From this point of view, the considered integral parameter S/Ra allows us to analyze the roughness of metallic surfaces at different orientations and, thus, to predict their corrosion resistance. For instance, this parameter differs 8.5 times for the Ni {111} and {011} orientations at 400 K. However, calculated roughness values are quite small suggesting that GR could provide stable corrosion protection in all considered cases. Sufficient enhancement of the corrosion characteristics of the GR/Ni system compared to free Ni could be predicted as well. Of course, this simulation presents an ideal case of GR deposition on a Ni single crystal defect-free surface. Defects of both Ni (vacancies, grain boundaries, triple junctions, etc.) and GR (multilayer structure, GR sheets’ island size and their boundaries, impurities) can modify the outcome. It can be noted that in the case of polycrystalline Ni, the considered {111}, {011}, and {001} orientations are most probable for surface formation due to the highest packing density. As was shown by preliminary calculations, the structure of the planes with less packing density will transform into a cluster type, where each cluster approximately represents the three most packed planes. It means that for the simulation of a Ni polycrystalline structure, it is necessary to average the characteristics of the three different substrate orientations {111}, {011}, and {001}, taking into account some deviations in the areas when grain boundaries come to the surface.

5. Conclusions

Relaxation and reconstruction of the Ni surface in the {111}, {001}, and {011} orientations with and without GR lamination at 300 K and 400 K were analyzed using molecular dynamic simulation. The simulation results show that in the presence of GR, a higher packing density of substrate suppresses the structure modification in the plane with the surface but allows the change in direction perpendicular to the surface. It was found that the substrate in the {111} orientation provides the highest structural and thermal stability of the GR/Ni system. Relaxation becomes more pronounced when the orientation of the substrate changes according to the sequence {111} → {001} → {011}; as such, the maximum modification of the interplanar distance increases from 1% to 6%. Deposition of GR on the Ni substrate in any orientation lowers the maximum interplanar distance change compared to the bulk by less than ~1%. It was found that 400 K is the most suitable temperature for GR deposition on Ni since it provides minimum structure changes. At this temperature, the reconstruction of GR is taking place at all studied orientations of the Ni substrate, which is manifested in the expansion of GR layers up to 0.2% with an increase in the average interatomic distance of up to 0.6% compared to the equilibrium state of GR. The crystallographic orientation of the Ni substrate and temperature predetermine the surface roughness of the GR/Ni system. The system formed by GR deposition on Ni {111} at 400 K is characterized by the least roughness.