Adsorption of dysprosium on the graphite (0001) surface: Nucleation and growth at 300 K

The interaction of magnetic materials with graphitic surfaces, from single-layer graphene to few-layer graphene to bulk graphite, holds possibilities for designing and controlling magnetic properties at the nanoscale. For graphene, in particular, it has been proposed that the carbon sheet may play a valuable role in new spintronics devices, or as a support or modifier of magnetic materials.^1–5 In most cases, it is necessary to achieve a good contact between the magnetic material and the carbonaceous substrate, in the form of a thin, essentially two-dimensional (2D) metal layer.¹ However, this has proven difficult since the adsorption energy of a metal on a graphitic surface is often much less than its cohesive energy, a factor which promotes three-dimensional (3D) growth of metals.⁶ 

In terms of magnetic properties, rare earths such as Dy attract special interest because of their high magnetic moment. This paper is a study of Dy islands that form when Dy is deposited on bulk graphite at room temperature. It sheds light on the fundamental energetics and mechanisms involved in adsorption, diffusion, nucleation, and growth, which may contribute to an ability to predict and manipulate Dy nanostructures, and the conditions under which they form.

There are two prior sets of studies of Dy’s interaction with related substrates. In one set, the substrate was graphene grown on SiC(0001).^7–9 There, it was shown that deposition of Dy at 660 K produced compact three-dimensional (3D) islands, with (mainly) triangular shapes indicative of fcc atomic structure, rather than the bulk hcp structure.⁷ Smaller, more irregular islands were formed at 300 K and were very stable against coarsening when annealed at 600 K.^8,9 In the second set, the substrate was highly oriented pyrolytic graphite (HOPG).¹⁰ A 20 nm film of Dy was deposited, corresponding to about 70 atomic layers if distributed uniformly. Annealing to 1200 K produced a carbide-like surface species, consistent with the high propensity for rare earths to form carbides.¹¹ 

The experiments were performed in an ultrahigh vacuum (UHV) chamber with base pressure 5 × 10⁻¹¹ mbar. The graphite samples were HOPG of grade ZYA in experiments at low coverage, i.e., below 0.8 monolayer (ML) of Dy, or ZYB in experiments at higher coverage. The clean graphite surface was prepared by tape-cleavage in air, followed by transfer into the UHV chamber where all subsequent experiments took place. Samples were heated in the manipulator to either 800 K for 60 min (low coverage experiments) or 1000 K for 20 min (high coverage experiments) in UHV to remove contaminants and then transferred to the scanning tunneling microscopy (STM) stage for subsequent Dy deposition and STM imaging. Dy was deposited via physical vapor deposition from a Mantis QUAD-EV-C Mini e-beam evaporator, using a Mo crucible lined with pyrolytic boron nitride. For heating, the crucible was biased at +2 kV with respect to an electron filament mounted parallel to and near the top of the evaporation target, with total power typically 25 W.

In a study of Cu/graphite, we showed that this evaporator, under similar conditions, produced a significant fraction of metal ions that damaged the graphite surface and influenced nucleation. We also showed that this could be circumvented by heating the evaporator, then shutting off the high voltage and/or filament current during deposition, with only a small drop in flux over times up to 10 s.¹² For all experiments reported herein, we employed this protocol. The deposition time was constant at 10 s; coverage was varied by adjusting power and hence flux.

In STM, tunneling parameters were typically in the range −1.0 V to −1.2 V tip bias and 0.15 nA–0.25 nA tunneling current for low coverage experiments, and 0.70 V–1.3 V and 0.13 nA–0.77 nA for higher coverage experiments. The STM tip was electrochemically etched W. All STM images were planed for data analysis. Other details of STM experiments and data analysis were the same as reported elsewhere.¹² 

First-principles energy calculations with non-local van der Waals correction were performed based on density functional theory (DFT) using VASP.^13,14 The exchange and correlation energy functional adapted the opt-B88 scheme developed by Klimeš et al.^15–17 This functional has been verified to describe accurately the energy and other properties of graphite and metals.^15,17 The electron-ion interaction was described by the projector augmented wave method.¹⁸ The energy cutoff for the plane wave basis set used in the calculation was 600 eV. The calculated fundamental properties of graphite agreed well with the corresponding experimental values.¹⁹ The calculated lattice constants of hexagonally close-packed Dy were 0.357 nm for a and 0.558 nm for c, respectively, in good agreement with experimental data of 0.359 nm for a and 0.565 nm for c.²⁰ The calculated cohesive energy was 3.58 eV/atom, somewhat larger than the experimental value, 3.04 eV/atom.²⁰ 

The graphite substrate was modeled by a slab with a 6 × 6 unit cell in the xy plane and four layers along the (0001) direction, plus enough vacuum (1.57 nm) to avoid interaction between the slab and its images under the periodic boundary condition. Spin polarization and dipole correction were considered in all calculations. A Γ-centered k-point grid of 15 × 15 × 1 was used for Brillouin zone sampling to ensure energetic convergence. During geometric optimization, the bottom three layers were fixed at their bulk positions, while carbon atoms in the top layer and metal adatom were relaxed fully with a force tolerance of 0.1 eV/nm.

Table I shows calculated values of the adsorption energy of a Dy atom at different sites on the graphite surface. The adsorption energy at the favored site is 1.90 eV. This is higher than the value of 1.47 eV calculated for Dy on graphene,²¹ which reinforces a trend noted elsewhere: For a given metal, its calculated adsorption energy on graphite exceeds that on graphene.¹⁹ The adsorption energy is about 2/3 the cohesive energy of bulk Dy, 3.04 eV. The adsorption energy falls in the range of chemisorption and implies an essentially infinite surface lifetime for a Dy atom at 300 K. By all measures, the interaction of Dy with graphite is reasonably strong.

The Dy atom is most stable at the hollow (H) site, in the center of the hexagon of carbon atoms. (See inset in Table I for identification of adsorption sites.) This is the same as the site predicted for Dy and other rare earths on free-standing graphene.⁹ The diffusion barrier, E_d, can be approximated as the difference in adsorption energies along a path that links minimum-energy sites. In this case that path leads from the H site, over T_α and back to H. The energy difference between the H and T_α sites is 0.048 eV.

Fig. 1 shows representative STM images of Dy islands on graphite terraces, for coverages spanning 0.12–2.5 ML. The islands consistently have a flat base with an irregular footprint. The flat base is decorated with smaller upper features, many of which are found near the island edges. The base is 0.87 ± 0.02 nm high, based on a sample size, N, of 36 islands. This height is invariant with island size. To first order, it can be compared with the 0.282 nm spacing between close-packed planes in bulk Dy, which yields a thickness of 3 atomic layers. Thus, the base can be described as quasi two-dimensional (2D), which is not unreasonable in light of the relatively high adsorption energy found in Sec. III.

In many cases a honeycomb pattern with amplitude 0.04 ± 0.01 nm (N = 144 islands) can be observed on the base, as shown in Fig. 2. This is a moiré pattern between the Dy film and the graphite support, which indicates that the Dy film is atomically well-ordered. The modulation periodicity is 1.23 nm. In moiré patterns such as this, a rotation between the unit cell vectors of the overlayer and the substrate, along with possible distortion of the overlayer, brings the two lattices into coincidence. Unfortunately, atomic-scale resolution of the graphite support could not be obtained in the presence of Dy islands, so the rotational angle is unknown. However, the existence of a coincidence lattice provides further evidence of significant interaction between Dy and the carbon substrate.

At high Dy coverage, Dy islands become very close but often retain their individual identity, as illustrated in Fig. 3. For example, the oval in Fig. 3(a) encompasses 3 islands separated by thin gaps, indicating a barrier to coalescence. We attribute this to the different epitaxial orientations of the Dy base, associated with the moiré pattern. Only when orientations happen to match can adjacent islands merge easily.

The bases often have taller dot-like features around their edges. The height of these dots corresponds to 1 Dy atomic layer above the base. Examples are visible in Figs. 1 and 4. Larger, layer-like features can also cover part of the islands’ interiors. These are 1.13 ± 0.02 and 1.42 ± 0.02 nm above the graphite. Figure 4 provides examples of these features, including line profiles that illustrate heights. The fact that each layer is higher by 0.28 nm indicates that single Dy layers are populated successively, above the base. With increasing Dy coverage, it becomes more probable to find these partially filled layers. The layers can exist near the center of an island, or emanate from an edge. When more than one layer exists atop the base, layers are stacked in a wedding-cake morphology such that each level below the top one is partially exposed.

For comparison, Fig. 5 shows the result when Dy is deposited at 800 K. Now the Dy islands are much taller and more three-dimensional; the three-layer mesa-like base is no longer evident. This indicates that the quasi-2D shape at 300 K is metastable. Indeed, the equilibrated Wulff shape of crystalline islands should be relatively tall compared to their width for a metal weakly bound to the substrate. The footprint should also be geometric, e.g., hexagonal for hcp crystal structure. Presumably, upward diffusion beyond the first three layers becomes facile at the higher temperature, allowing islands to achieve a more equilibrated shape during growth, than that which can be attained at 300 K.

The island density, N_isl, is shown as a function of Dy coverage in Fig. 6. The density increases from a value of (2.4 ± 0.3) × 10⁻⁴ nm⁻² at 0.12-0.15 ML and saturates at 6.0 × 10⁻⁴ nm⁻² at 1.2 ML.

The island size distribution is shown in Fig. 7 for two coverages, 0.13 and 1.2 ML. In both cases the distribution is monomodal.

These observations are analyzed in Sec. V.

First we assess the island nucleation mechanism by comparing the observed value of N_isl and its coverage dependence for deposition at 300 K with theories for homogeneous nucleation. It is useful to note that the area of the surface unit cell for graphite is Ω = 0.052 nm². Multiple experiments for deposition of 0.12-0.15 ML provide quite consistent estimates for N_isl (see Fig. 6). Averaging these data yields N_isl = (2.4 ± 0.3) × 10⁻⁴ nm⁻² (or 1.3 × 10⁻⁵ per unit cell or per adsorption site) for coverage θ = 0.13 ± 0.1 ML corresponding to a deposition flux of F = 0.013 ML/s.

To analyze this behavior, we use the estimate for the terrace diffusion barrier of E_d ≈ 48 meV, obtained from DFT analysis described in Sec. III. Assuming an attempt frequency of ν = 10^12.5 s⁻¹ yields a rate of h = νexp[ − E_d/(k_BT)] for adatom hopping associated with terrace diffusion. A key parameter for nucleation is h/F. Using the experimental values of T = 300 K and F = 0.013 ML/s, one has h = 10^11.7 s⁻¹ and h/F = 10^13.6 s⁻¹. For homogeneous nucleation and growth of islands with critical size i (where islands of size >i atoms are stable), the island density satisfies N_isl ∼ θ^(1−χ)/2(h/F)^χ where χ = i/(i + 2) in the regime of low coverage.²² 

For this system, it is reasonable to anticipate that homogeneous nucleation and growth of islands is at least close to irreversible, corresponding to critical size i = 1. Theoretical analysis indicates that the temperature for transition to reversible island formation is proportional to E_rev = (2/3)E_d + E_NN, where E_NN > 0 is the strength of the nearest-neighbor interaction.²² More precisely, to ensure a transition temperature above 300 K for a standard deposition flux of 0.005 ML/s and attempt frequency of ν = 5 × 10¹² s⁻¹, one requires that E_rev ≈ 0.55 eV or above. Thus, for Dy on HOPG with E_d ≈ 0.05 eV, one requires that E_NN ≈ 0.52 eV or above. A crude estimate of E_NN as 1/6 of the bulk cohesive energy yields a value of 0.51 eV, and it is plausible that the true E_NN for a pair of Dy on HOPG is stronger.

For i = 1, one has that N_isl ∼ (h/F)^−1/3 for fixed coverage, and simulations indicate that N_isl ≈ 2.18 × 10⁻⁵ per site for θ = 0.13 ML.^22,23 Thus, the experimental value of N_isl, 1.26 × 10⁻⁵ per site, is close to but slightly reduced from classic i = 1 prediction (reduced by a factor of r = 0.58).

The lower experimental value could be achieved by choosing an attempt frequency which is higher by a factor of 5.1. Such a prefactor, just above 10¹³ s⁻¹, cannot be excluded, but one generally expects low prefactors for low diffusion barriers (as in this system). Alternatively, the lower experimental N_isl could be naturally explained by deviations from classic i = 1 behavior due to either significant mobility of adsorbed Dy dimers (and perhaps other small Dy clusters) or due to the onset of reversibility in island formation. To obtain a sense of the impact of dimer mobility or reversibility in reducing N_isl, we note the following. If monomers and dimers have comparable mobility, then N_isl ∼ (h/F)^{−2/5 22,24} which implies a reduction in N_isl from the classic i = 1 value by a factor of r = (h/F)^−1/15 = 0.12 for the above h/F = 10^13.6, i.e., the reduction is too strong. For facile dissociation of dimers but stable trimers corresponding to i = 2, one has that N_isl ∼ (h/F)^−1/2, reducing N_isl for the classic i = 1 value by a factor of r = (h/F)^−1/6 = 0.005, which is much too strong. In conclusion, the observed reduction by r = 0.58 relative to classic i = 1 could be induced by some dimer mobility²⁴ (still significantly lower than adatom mobility) or by slight reversibility in island formation.

Another way of accounting for these effects is by introducing a scaling exponent χ = 1/3 + δχ which is slightly above the value of 1/3 for i = 1, so that (h/F)^−δχ = r = 0.58. This implies that δχ = 0.018 so χ = 0.351 and (1 − χ)/2 = 0.325. As a result, N_isl = c(h/F)^−0.351θ^0.325. Choosing c to recover the experimental N_isl at 0.13 ML and recognizing that F = θ/10 in experiment, one obtains N_isl ≈ 5.01 × 10⁻⁴ θ^0.676 per site, which should apply in the lower coverage regime where the footprint of the islands covers a small fraction of the substrate area (i.e., the “point island” regime²²). This curve is shown in Fig. 6 up to a coverage of 0.4 ML, beyond which one expects that the point island model will significantly overestimate N_isl because it underestimates the capture probability at existing islands. (We emphasize that for the Dy islands with a 3-layer base, the areal coverage or fraction of the substrate covered is only ∼θ/3 in contrast to θ for 2D islands. Thus, nucleation persists to higher coverages in the 3-layer case.) There is a good agreement between the model and the data in this regime.

Next, we analyze experimental observations for the island size distribution, N_s, which gives the density of islands of s atoms, so that $Nisl=∑s>1Ns$ and $θ=∑ssNs$⁠. If s_av ≈ θ/N_isl denotes the average island size (measured in atoms), then the island size distribution is naturally written in scaled form, N_s ≈ (N_isl/s_av) f(s/s_av, θ), where the scaling or shape function f(x) satisfies ∫_x>0 dxf(x) = 1 and ∫_x>0 dxxf(x) = 1.²² For homogeneous nucleation at lower coverages around 0.1 ML, f(x) is monomodal with a peak height around x = 1 increasing from about 0.75 for i = 1 to 1.1, 1.3, and 1.5 for i = 2, 3, and 6, respectively (and also f(x) narrows with increasing i). The peak also increases to about 1.1 with high dimer mobility.²⁴ There is a significant population of smaller islands reflected in f(0) ≈ 0.3 for i = 1, with f(0) quickly decreasing with increasing i (e.g., to 0.15 for i = 2). The experimental f(x) shown in Fig. 7(a) for 0.13 ML has a peak height around 1.0 consistent with slight deviations from classic i = 1 behavior. However, there is a significant depletion in the population of small islands relative to classic i = 1 behavior. This feature together with an increase in peak height and shift to smaller x-values has been seen previously and attributed to post-deposition diffusion and coalescence of small clusters.²⁵ Fig. 7(b) shows the experimental f(x) for a much higher coverage of 1.2 ML. The peak height is similar to 0.13 ML, but the peak and the weight of the distribution have shifted to higher (scaled) island sizes. This feature might be expected, as inspection of STM images reveals that island growth is impeded by impingement upon other islands (without coalescence) in this regime. As noted in Sec. IV, we attribute this to epitaxial mismatch between different islands. This would tend to induce a sharper cutoff on the right side of the distribution.

Despite the complexity of multilayer growth in this system (described in more detail below), simplified modeling can produce significant insights, e.g., into the extent of interlayer transport. To provide some context, such models have been traditionally developed for systems with equivalent layers where growth in higher layers proceeds before growth in lower layers is complete. If θ_j denotes the coverage of layer j, then P_j = θ_j − θ_j+1 denotes the exposed fraction of layer j. In the extreme case where interlayer transport is absent, one finds a Poisson distribution for P_j and film roughness W satisfies W² = d²θ, where d is the interlayer spacing.²² This modeling has been refined to account for downward interlayer transport by specifying that a fraction, α, of atoms deposited in each layer reach the next lower layer (so α = 0 recovers the above case).^26–28 Then, for rough growth with small α, one can show that W² ≈ (1 − 2α) θd².^27,28 Although not well-recognized, this result must be modified for larger α [see  Appendix A].

A special feature of the Dy/graphite system is that atoms deposited on the substrate form bases consisting of three atomic layers (rather than single layer islands). This implies facile transport from the substrate upward to the second and third layers. However, on top of these bases, growth occurs in single layers. Thus, it is convenient to refer to the fraction, P_j, of exposed level, j (rather than layer), where j = 0 is the substrate, j = 1 is the top of the 3-layer base, j = 2 is the top of the single layer island on the base, etc. Experimental values for these populations or fractions of exposed levels for coverage θ = 1.2 ML are given in Table II.

To elucidate this behavior, we develop a theory in the spirit of the more traditional one mentioned above, for multilayer growth with equivalent layers. However, here we naturally specify that all atoms deposited on the graphite substrate are incorporated immediately into 3-layer bases. A fraction α^∗ of atoms deposited on top of bases hop down to become incorporated in the base and the remainder stay on top of the base forming single-layer islands. Of atoms deposited in higher layers, a fraction α hop down to be incorporated in the next lower layer, and the remainder stay in the layer where they were deposited, forming or joining single-layer islands in that layer. See the schematic in Fig. 8. With this prescription, the rate equations describing the evolution of the level populations are

Simple modification of these equations is possible if one wishes to exclude population of level 4 (i.e., to enforce P₄ = 0) consistent with experiment.

It is natural to first consider the extreme case of no interlayer transport (other than upward transport of atoms deposited on the substrate to form 3-layer bases). In this case α^∗ = α = 0, and the above equations partially decouple in the sense that P_j is coupled only to P_k<j so they can be recursively solved exactly to obtain

As shown in Table II, the resultant values of P_j agree roughly with experimental values, but the fact that P₀ is higher in the model than in experiment (0.670 vs. 0.631, respectively) is a clear indication that net downward 3 transport with α^∗ > 0 is operative in the Dy/graphite system.

To refine the above simplest model, we introduce interlayer transport, initially setting α^∗ = α. We also enforce P₄ = 0. Numerical integration of the associated rate equations indicates that good agreement with experiment comes from choosing α = 0.65, with the resultant values of P_j shown in Table II. A slightly better fit is achieved by allowing independent values α^∗ = 0.69 and α = 0.85. In either case, this analysis indicates that there is significant net downward transport in the Dy/graphite system, which impacts the film height distribution.

The existence of small single-layer features at the edges of islands was noted in Sec. IV. Their presence indicates that nucleation is somewhat more likely near islands’ edges than in islands’ interiors. This could divert some Dy atoms from reaching and crossing over the edge of the base, leading to α^∗ < α as suggested by the modeling.

One major conclusion from this work is that Dy island nucleation is homogeneous. This is based upon the measured value of N_isl at 0.13 ML, in comparison with the value predicted from a point-island growth model using the diffusion barrier from DFT and i = 1. The agreement is quite good, with the experimental value being slightly lower than the predicted value. Elsewhere, we have argued that homogeneous nucleation of metals on graphite terraces may be less common than previously thought, with heterogeneous nucleation playing a larger role than expected due to the susceptibility of graphite substrates to inadvertent ion beam damage.¹² We have argued that if the measured island density is significantly higher than the predicted value, then heterogeneous nucleation should be considered. If the opposite is true, as in this case, then homogeneous nucleation is viable with some adjustments. Here, we have shown that the degree of adjustment is too extreme if i = 2 or if dimer diffusion is fully operative, though a small contribution from one or both of these mechanisms could explain the discrepancy.

The island size distributions are also consistent with homogeneous nucleation. Their monomodal shapes are quite different from the monotonically decreasing shape found in our previous study of Cu on graphite, where Cu ions in the metal atom flux introduced defects in the graphite surface during deposition. Precautions are taken in the present work to prevent such damage, consistent with the radically different distributions observed.

As noted in Sec. I, Hupalo et al. studied a related system, Dy deposited on monolayer graphene supported on SiC, at 300 K.⁸ They reported N_isl = 2 × 10⁻² nm⁻² at 0.4 ML, which is two orders of magnitude higher than our value of N_isl at 0.13 ML. While a somewhat higher value of E_d = 0.125 eV was calculated for that system, this and other differences—in experimental flux and coverage—are insufficient (by far) to explain the difference in N_isl. Instead, we propose that the explanation lies in the influence of the number of carbon layers, n, on nucleation, growth, and stability of metal nanoclusters. There is now broad evidence for this dependence, with N_isl generally decreasing as n increases.⁶ This has been demonstrated in studies of Au, Ag, and Pd nanoclusters on n-layer graphene sheets supported on SiO₂/Si substrates. This trend, decreasing N_isl with increasing n, describes exactly the trend in going from Dy/graphene (n = 1) in the prior work, to Dy/graphite (n → ∞) in this study. One explanation involves charge transfer from the metal to the underlying carbon and consequent Coulomb and dipole repulsion terms, which have significant effect at small n due inefficient screening.²⁹ 

Finally, we comment on the morphology of the Dy islands. Experiment indicates that growth is quasi-2D. The islands exhibit a 3-layer base, and single atomic layers grow on top of the base. DFT shows that the interaction between Dy and the graphite is reasonably strong, providing a rationale for the quasi-2D base. Analysis of the level populations reveals significant downward interlayer transport, which facilitates growth of the base. On the other hand, the quasi-2D structure is metastable, based on comparison with the taller, more compact 3D islands that form at 800 K; the most natural explanation is that upward transport (beyond the 3-layer base) is kinetically limited at 300 K.

It is possible that the base is limited to 3 atomic layers because of strain between Dy and graphite lattices. It is well-established that strain influences growth in heteroepitaxial systems. For the growth mode known as Stranski-Krastanov (SK), the film wets the substrate (growth is smooth) up to a critical thickness, beyond which growth is 3D.^30,31 In SK systems strain is often manifest in a moiré pattern at or below the critical thickness.^32–35 Strictly speaking, SK is a thermodynamic (equilibrium) picture of growth, but the same factors that lead to SK growth can influence structures formed under kinetically limited conditions.

We have shown that experimental measurements of the value of Dy island density, and coverage-dependence of Dy island density at low coverage, are well described by analytic theory for homogeneous nucleation, using energetics provided by DFT. The Dy islands have a quasi-2D 3-layer base which exhibits atomic order coincident with the graphite substrate. This natural tendency toward quasi-2D growth may be useful in situations where good contact between a magnetic metal and a graphitic substrate is desired.

The experimental and DFT components of this work (Secs. III and IV) were conducted by the following authors: E.J.K., H.L., A.L.-R., M.W., Y.Z., C.-Z.W., M.C.T., and P.A.T. Experimental and DFT effort was supported by the US Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. Research was performed in part at the Ames Laboratory, which is operated by Iowa State University under Contract No. DE-AC02-07CH11358. DFT was performed, in part, with a grant of computer time at the National Energy Research Scientific Computing Centre (NERSC). NRSEC is a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. H.L.’s DFT was supported in part by the National Science Foundation of China (NSFC) Grant No. 11575230, and Y.Z.’s participation in the experimental effort was supported by the China Scholarship Council. The modeling and analyses described in Sec. V were performed by J.W.E., with support from NSF Grant No. CHE-1111500.

For equivalent layers with deposition at rate F so θ = Ft and a fraction α of deposited atoms reaching the next lower layer, one has that

Separate equations are needed for j = 0 and j = 1. Neglecting these different equations yields the simple result for W² = (1 − 2α) θd² quoted in the text. However, reliable results for smoother growth with larger α must account for the modified equations, noting that clearly this expression for W² cannot apply for α > 1/2.

Data used to generate this manuscript are available via the DOI: http://dx.doi.org/10.17039/ameslab.dmse.2016.DS2/1240605.