Shapes of Fe nanocrystals encapsulated at the graphite surface

We describe and analyze in detail the shapes of Fe islands encapsulated under the top graphene layers in graphite. Shapes are interrogated using scanning tunneling microscopy. The main outputs of the shape analysis are the slope of the graphene membrane around the perimeter of the island, and the aspect ratio of the central metal cluster. Modeling primarily uses a continuum elasticity (CE) model. As input to the CE model, we use density functional theory to calculate the surface energy of Fe, and the adhesion energies between Fe and graphene or graphite. We use the shaft-loaded blister test (SLBT) model to provide independent stretching and bending strain energies in the graphene membrane. We also introduce a model for the elastic strain in which stretching and bending are treated simultaneously. Measured side slopes agree very well with the CE model, both qualitatively and quantitatively. The fit is optimal for a graphene membrane consisting of 2–3 graphene monolayers, in agreement with experiment. Analysis of contributions to total energy shows that the side slope depends only on the properties of graphene/graphite. This reflects delamination of the graphene membrane from the underlying graphite, caused by upward pressure from the growing metal cluster. This insight leads us to evaluate the delamination geometry in the context of two related, classic models that give analytic results for the slope of a delaminated membrane. One of these, the point-loaded circular blister test model, reasonably predicts the delamination geometry at the edge of an Fe island. The aspect ratio also agrees well with the CE model in the limit of large island size, but not for small islands. Previously, we had speculated that this discrepancy was due to lack of coupling between bending and stretching in the SLBT model, but the new modeling shows that this explanation is not viable.


Introduction
There are numerous situations where it is important to understand or utilize the interaction between a metal, and a two-dimensional (2D) material (or by extension a heterostructure of 2D materials, or a three-dimensional (3D) van der Waals material). One example is electrical connections or heat sinks for device applications, where metal architectures with stable and high-area contacts to the 2D material are needed [1,2]. Another example is tuning the electronic properties of the 2D material, where dopants and intercalants can modify the Fermi level via charge transfer [3,4]. A third example is the area of magnetism, where there are many possible applications [5]. A particularly exciting one is the creation of spring magnets using a 2D material such as graphene as a sharp and stable interface between atomically-thin films of magnetic metals [6]. In almost all of these applications, it is desirable to maximize the contact area between the metal and the 2D material, i.e. to achieve 'flat' growth of the We analyze not only the side slope of the sides of the islands, but also their aspect ratio and a related dimensional ratio, which we term the lateral ratio, as a function of island height. All results for Fe are compared with those of Cu. The similarities and differences provide insights into the factors controlling shapes of these encapsulated islands.
In the following, section 2 provides some details of the experimental and computational methods. Section 3 presents results from experiment, from DFT, and from CE. In applying CE, we compare input from SLBT and a related model, which allows us to assess the bending contribution. Section 4 discusses these results. Auxiliary information is available in the supplemental information (SI) is available online at stacks.iop.org/NJP/22/ 023016/mmedia, as noted in the text. Also, because the paper contains many variables, we define them in table 1 for easy reference.

Experimental methods
We have described the experimental methods fully elsewhere [12]. In brief, all of the experiments were conducted in an UHV chamber. Encapsulated Fe clusters were grown beneath the surface of commercially available graphite (HOPG, ZYA grade). Necessary conditions for Fe encapsulation include (i) activating the graphite surface with atomic-scale defects via Ar + ion bombardment, and (ii) depositing Fe on the activated graphite surface that is held at elevated temperature (T dep ). Our group has demonstrated that these conditions are effective for encapsulating a variety of metals, including Cu [10], Ru [11], and Dy [9]. Fe encapsulation is operational in a narrow T dep window of 875-900 K, with 900 K being the optimal temperature where the extent of encapsulation is high with minimal bare Fe on top of graphite [12]. We have extensively characterized the encapsulated Fe clusters via scanning tunneling microscopy (STM) and X-ray photoelectron spectroscopy. From STM images we obtain key dimensions of encapsulated Fe clusters, including island height (h), island top width (d), and annular width (a). See figure 1. Specifically, h is determined from topographic STM images, while d and a are measured from derivative STM images. Derivative images offer vivid contrast of the islands and thus enable accurate lateral measurements. Additional experimental details are described in the SI.  [17] for the Fe-graphite and Fe-graphene systems. The projector-augmented-wave pseudopotentials [18] generated and released in 2013 by the VASP group were used for the electron-core interactions. The Γ-centered k mesh depended on the supercell size and will be listed for each specific system. The convergence tolerance for the force on each relaxed atom was set to be 0.1 eV nm −1 . For a surface system, which was modeled as a periodic slab with a specified lateral supercell size, the thickness of vacuum space between two adjacent slabs was not less than 1.5 nm. We also considered spin-polarization effects and dipole corrections unless noted otherwise. During the energy minimization with structure optimization for a Fe-C surface system using the optB88-vdW functional [19], we found that a large number of steps were needed for the relaxation of electronic degrees of freedom as well as all-band simultaneous update of orbitals.
Our goal was to obtain values for two distinct quantities for use in CE modeling: Surface energies per unit area (γ) of pure Fe, and adhesion energies per unit area (β FeG ) for mixed Fe-C systems. (These mixed systems are denoted FeGt for Fe-graphite and FeGn for Fe-graphene, or FeG for both together.) For γ of pure Fe, we considered both the Perdew-Burke-Ernzerhof (PBE) generalized-gradient approximation (GGA) [20] and the optB88-vdW functional, where the exchange functional is optimized for the correlation part [19], to approximately account for dispersion interactions. Based on comparison with an experimental value, we will conclude in section 3.2 that PBE is more appropriate than optB88-vdW for purpose of calculating γ of Fe. On the other hand, β FeG is a parameter of mixed Fe-C systems, where the functional must do well at describing both a metal and a van der Waals material. From our previous DFT calculations [16], the optB88-vdW functional can reproduce very well the experimental lattice constants, cohesive energy, and exfoliation energy of graphite (the AB-stacked hexagonal structure [21]) as well as experimental lattice constant of the GML. Furthermore, we have performed benchmark calculations for the bulk properties (lattice constants and cohesive energies) of α-Fe (bcc phase), γ-Fe (fcc phase), and ε-Fe (hcp phase), obtaining good agreement with experiment using optB88-vdW. The calculated bulk cohesive energies indicate a hierarchy of stabilities, wherein bcc-Fe is most stable and fcc-Fe is least stable. The bcc, hcp, and fcc structures of Fe are ferromagnetic, nonmagnetic, and antiferromagnetic, respectively. This result is consistent with Herper et al's DFT calculations [22]. The ability of the optB88-vdW functional to simultaneously give good results for bulk Fe, and for graphite and graphene, justifies its choice for obtaining β FeG in section 3.2.

Computational techniques: CE model
We adopt the same approach to CE modeling as in a previous paper [13]. In brief, we approximate the shape of the Fe cluster as a cylinder, with its top and bottom circular faces (corresponding to hcp(0001) facets) contacting a graphene membrane and a rigid graphite substrate, respectively. The graphene membrane undergoes elastic stretching and bending deformations to accommodate the cluster, causing a strain energy U e to be induced. An analytic expression for U e (treating the bending and stretching components independently) is obtained from the cylindrical SLBT model, which provides an excellent analog for the situation described here [14]. In section 3.3, we will first consider the case where U e has only a stretching component, derived from the SLBT model (equation (5)). We will then consider the effect of bending separately (equation (6)). Finally, we will present a different model for U e in which stretching and bending are incorporated simultaneously (equation (7)) but the real island shape is not approximated as closely as in the SLBT model.
It should be noted, here, that all of the expressions for U e treat the graphene with a Young's modulus, Y. In this approach, the in-plane stiffness of the membrane, E 2D , is related to the three-dimensional (3D) Young's modulus Y by where t is total membrane thickness and ν is Poisson's ratio. The flexural rigidity, D, is related to the Young's modulus Y by Many studies [23][24][25][26] experimentally characterizing Y obtain it by measuring E 2D , then back-calculating Y using equation (1). One group [27] experimentally measured D of multi-layered graphene (  L 8 c ), and found that it obeys equation (2), when Y is obtained as above. However, other groups have computed D of single-layer graphene from either DFT [28][29][30] or an empirical potential approach [31][32][33] and found it to be an order of magnitude lower than that predicted by equation (2). In this case, graphene cannot be considered as an isotropic material with a single Y value, and E 2D and D must be treated as separate properties that cannot be characterized by equations (1) and (2). In light of the above findings, it is quite possible that multilayered graphene thinner than eight monolayers ( < L 8 c ) behaves intermediate between the two above limits. In the absence of a well-defined experimentally-measured flexural rigidity in this range, we choose, in this work, to treat the graphene as a transverse isotropic material. As shown in equation (3), the total energy of the system, Π, is defined as U e plus a set of energy terms that represent the interfacial and surface (IS) components of Π: Here, U Fe >0 is the total surface energy of the clean Fe cluster; U GnGt >0 denotes the total energy cost due to reduced adhesion at the graphene-graphite interface; and U FeG >0 is the total adhesion energy of the Fegraphene and Fe-graphite interfaces combined. (Each of these terms involves γ or β as defined in section 3.2.) The expressions for individual U terms, and relevant materials parameters, are given in table 2.
The equilibrium shape is obtained by minimizing Π for fixed Fe cluster volume V. We then characterize cluster geometry for various sizes, repeating this analysis over a range of V (although we note that cluster shape is size-independent if one ignores the energy cost of membrane bending). However, rather than using V as a variable in the presentation of results, we use island height h, because this allows a more accurate and direct comparison with experiment. Calculations are carried out using Mathematica © . Details for the SLBT models are given in [13] and its associated SI. Details for the mixed bending and stretching model are given in the SI of the present paper.
For the SLBT model, one can in principle consider two variants which differ in the distribution of strain in the graphene membrane: Either the graphene in contact with the top face of the Fe cylinder is allowed to stretch (free), and strain is distributed over the entire island; or the graphene membrane is restricted from stretching (clamped) on top of the cylinder due to strong interaction, and strain is confined to the annulus [14]. Elsewhere, we have shown that there is minimal difference in the results from the two models for Cu clusters [13], and we find that to be true also for Fe clusters. Therefore, we focus on the free SLBT in this text. Results for Fe clusters in the clamped SLBT model are given in the SI. Figure 1 shows representative STM images and profiles of three Fe islands. It also includes a schematic that defines the key dimensions of an iron cluster encapsulated by a graphene layer: h is the island height; d is the diameter of the island top; a is the width of the sloping perimeter (annulus) of the graphene membrane; and t is the thickness of the graphene membrane, given by = t L t . C GML Here t GML is the thickness of a single GML in nm, and L C is membrane thickness in units of GMLs. In the schematic, the sides of the cluster are drawn vertical to be consistent with the cylindrical shape of the SLBT model. In reality, we believe that the sides are composed of various facets, and are only vertical on average.

Results from experiment
It should be noted that the profiles in figure 1 are plotted using different scales on the y-and x-axes, to make the representations compact. If scales were equal, it would become obvious that these islands are very wide and low. This is reflected in their aspect ratios (d/h), which range from about 5-40. Dimensions are shown in figure 2, for a dataset of 235 Fe islands. This is a significant expansion of the 140island dataset presented in an earlier report [13] of h versus a and h/a versus h for Fe. Island height h shows a strong linear correlation with annular width a (figure 2(a)). By extension, the ratio h/a (slope of annulus) is constant (figure 2(b)), at 0.27±0.04 for 235 Fe islands. This corresponds to an angle of 11.1°with respect to the surface plane. By contrast, the island top diameter d does not vary linearly with h (figure 2(c)). However, a different trend is revealed when the ratio d/h is plotted against h: d/h decays from large values at small h to an asymptotic plateau at large h (figure 2(d)). There is large scatter at small h. A similar trend is observed for the lateral ratio d/a, versus h (figure 2(f)). (We show and discuss d/a for completeness, even though it can be derived from the other two ratios.) The dimensional trends shown in figure 2 for encapsulated Fe clusters are very similar to those reported previously for encapsulated Cu clusters [13]. This reinforces a conclusion drawn elsewhere: encapsulation of Fe is very similar to that of Cu in general. However, there are two noteworthy differences. First, the heights of the Fe islands fall in a more limited range than those of Cu. Specifically, for Fe, h max =9 nm, whereas for Cu, h max =43 nm, so the Fe islands do not grow as tall. Second, for Cu, some islands have round tops rather than flat tops [10] whereas for Fe, only flat tops are observed. For Cu, the largest islands tend to have round tops, so these two differences may be related, e.g. perhaps taller growth requires adoption of the round top shape. At present the reasons for these differences between Fe and Cu are not understood and warrant further investigation.
Elsewhere, analysis of STM images of carbon atoms in the annulus region has shown that the graphene membrane consists of at least 2 GMLs, i.e. L C 2. Also, analysis of the footprints of the metal clusters has shown Table 2. Expressions for U-terms, and corresponding values of input parameters used in CE modeling. Here a, d, and h are dimensions of the islands as defined in figure 1. γ and β are surface and adhesion energies per unit area. Y and ν are Young's modulus and Poisson ratio, respectively. t GML is the thickness of a single graphene monolayer (GML) in nm, and L C is the thickness of the graphene membrane in units of GMLs. Subscripts FeTp and FeSd signify the top and side of the cylindrical Fe cluster, respectively. The components of the total strain energy, U i , are derived for models described and used in section 3.3. Additional subscripts on U e stand for stretching (s), bending (b), combined stretching+bending (s+b), and free (fr).  In (9), that the majority of clusters (62%) have hcp (hexagonally-close-packed) structure, with their basal plane parallel to the graphite interface. We therefore emphasize the hcp-Fe(0001) surface in the DFT calculations that follow.

Results from DFT
3.2.1. Surface energies, γ, of hcp-Fe(0001), bcc-Fe(100), and bcc-Fe(110) For reasons described in section 3.1, it is most relevant to consider the surface energy of hcp-Fe(0001). However, we include some results for bcc-Fe since bcc-Fe is the most stable phase under ambient conditions [36,37].
The surface energy g L of a slab (or an unsupported film) with two equivalent surfaces on both sides and with thickness L is calculated as where E L is the total energy of the slab in a supercell, N L is the total number of atoms in the slab, A is the area of the bottom or top face of the slab, and s bulk is the energy per atom in the bulk crystal. The slab thickness L is in units of atomic layers or monolayers. We note that to maximize accuracy in extracting surface energies for  ¥ L , we use our 'adjusting+observing' method [38][39][40] in which behavior of g L for a range of L are considered (see SI). This is different from the approach used by other groups, where the oscillatory behavior of g L with increasing L is not generally considered [41][42][43]. For a specific metal M, calculations over a range of the slab thickness L M are useful to eliminate quantum size effects, where the electronic properties of a metal slab oscillate as a function of the slab thickness [40,44,45]. Such oscillations are especially notable for hcp Fe(0001) slabs (see SI). Table 3 summarizes surface energies calculated in this work. Details of the calculations are given in the SI. An early experiment assessing the Fe surface energy gave a value of 2.417 J m −2 (near the melting point) [46], though the experiment did not yield values for specific surface planes. It appears that the surface energies from our PBE calculations are closer to the experimental value, and thus we will use the PBE values for surface energies of Fe. This conclusion is consistent with an analysis of the effect of van der Waals corrections on calculations of metallic surface properties [48]. Recently, Tran et al [43] generated an extensive database for various metals from DFT calculations using PBE GGA. Relevant values are also listed in table 3, showing good agreement between our PBE calculations and theirs. In addition, from table 3, the surface energy of bcc Fe(110) is lower than that of bcc Fe(100). This order is consistent with the result from Tran et al's database [43] but is opposite to early full charge density results from Vitos et al [47]. The value ( ) g =

Fe 0001
2.664 J m −2 for hcp Fe(0001) is used as input to the CE model, as described in section 3.3 and table 2.

FeGt and FeGn adhesion energies, β
The adhesion energy between two slabs (s1 and s2) bonded at an interface is a thermodynamic quantity defined as the energy required (per unit area) to separate the slabs and create two free surfaces. It can be calculated as [16] where E , s1 E , s2 and E s1s2 are the energies of slab s1, slab s2, and the s1s2 system, respectively; A is the s1s2 interface area.
For Fe(0001)-graphite (FeGt) or Fe(0001)-graphene (FeGn) systems, we use a1 1 supercell (in units of the graphitic lattice constant). For this configuration strain is very small, only 0.244% from optB88-vdW, which corresponds to a negligible tensile strain for the Fe film along the lateral direction. Our recent calculations for Cu(111)-graphite [16] have shown that the adhesion energy is not sensitive to the mismatch strain at the interface even when the strain is significantly large, and therefore we believe that the1 1 supercell calculations for the FeGt systems will be reliable. The k mesh is taken to be´51 51 1, and the energy cutoff is 600 eV. During relaxation for energy minimization, the bottommost GML of the graphite slab is always fixed. To cancel errors, we always use the same supercell and k mesh when separately calculating E , s1 E , s2 and E s1s2 for a specific hcp-Fe(0001) film thickness L Fe . For FeGt systems, the Fe(0001) slab s1 has variable thickness L Fe and the graphite substrate s2 always has L C =6. Also, we always choose the energetically most favorable adsorption site match between the Fe slab and the graphene/graphite slab for a given L Fe , as discussed in the SI.  (111) interacting with the graphene membrane and graphite, respectively [16]. This is reasonable, since Fe is generally found to have a stronger (more favorable) interaction with carbon than does Cu [49][50][51][52]. The β and γ values for Cu and Fe are compared directly in the SI. Table 3. Surface energies (in unit of J m −2 ) of bcc Fe(100), bcc Fe(110), and hcp Fe(0001) surfaces from different theories. For comparison, the experimental value sampling over various surface orientations is g =2.417 J m −2 near the melting point of Fe [46]. In the above calculations for adhesion energies, we do not consider spin polarization, because the hcp Fe(0001) slab is very weakly magnetic relative to bcc. (In fact, the magnetism per atom goes to zero with increasing L Fe , consistent with the nonmagnetic nature of bulk hcp Fe mentioned in section 2.2.) Furthermore, equation (5) shows that calculating β involves taking the difference between the energies of an hcp-Fe(0001) slab that is part of a FeGt or FeGn interface, and the energy of a free hcp-Fe(0001) slab, at given L Fe . We expect that the small magnetism of the hcp-Fe(0001) slab in both configurations, at given L Fe , will be nearly the same and will thus cancel out. This expectation has been verified for select values of L Fe .

3.3.
Results from the CE model 3.3.1. Results using the SLBT model for U e with stretching only With values for surface and adhesion energies from DFT in hand, the CE model can be implemented as described in section 2.3. Here we first consider a situation where the strain energy U e is due entirely to stretching, with negligible bending; this is described by equation (7) of table 2. For this model, the equilibrium shape is independent of island size as will be discussed in section 4.
The solid lines in figure 4 show CE model results in comparison to the experimental data. It can be seen that for a fixed value of L C , the CE model predicts constant values of slope h/a, aspect ratio d/h, and lateral ratio d/a, over most of the range of experimental heights, h. This agrees with the experimental trend for h/a, and it also agrees with the asymptotic trends for d/h and d/a in the limit of large h. (We take h 7 nm as representative of the asymptote; this cutoff is shown by the vertical red line in figures 4(b), (c).) Not only trends, but also absolute values of the ratios predicted by the CE model agree well with the experimental data. If the thickness of the top graphene membrane, L C , is treated as an independent variable, the range 1<L C <4 encompasses almost all the experimental values of h/a, and best agreement with the experimental average h/a is obtained with L C =2-3 in the free SLBT model, or with L C =2 in the clamped model (shown in the SI). This is the only significant difference between free and clamped SLBT. These values of L C agree with the experimental constraint, L C 2. As shown in figures 4(b), (c), for aspect ratio d/h or lateral ratio d/a, there is negligible variation in the CE results for different values of L C , compared with the scatter in the experimental data at large h. For both these ratios, CE values agree favorably with experimental data at large h.

Results using the SLBT model for U e with bending only
Previously, we speculated that the large values of d/h and d/a at small h (for Cu) were due to the contribution of bending to U e [13]. Qualitatively, bending should play a more important role for smaller islands, and this is exactly where the experimental data for d/h and d/a deviate from the CE results with SLBT-derived stretching only. While a mixed bending and stretching SLBT model can be derived for the SLBT model, it is analytically complex, and its application requires the use of numerical methods that are not straightforward to implement. The pure-bending and pure-stretching models are limiting cases of the mixed model. As such, if they do not fit the experimental data in the domains of island sizes in which they are valid (short islands for pure bending, tall islands for pure stretching), neither will the mixed model throughout the full domain of all island sizes. Based upon this rationale, we first consider the effect of bending as an independent contribution in the SLBT model (equation (8) in table 2). For the sake of clarity, we consider only L C =2. Conclusions would be qualitatively unaffected for any L C in the range 1L C 4, as shown in the SI. The long-dashed lines in figure 4 show the results of SLBT at small h, with bending strain only. With decreasing h, this model predicts that the ratio h/a drops steeply, d/h rises steeply, and d/a approaches a finite ordinate intersection. These trends can be rationalized as follows. As h decreases and bending becomes important, the graphene membrane stiffens and starts to act more like a solid 'plate'. This has two effects: (i) atop the island, the membrane exerts a greater downward force, increasing d/h; and (ii) at the outer edge of the overhanging annulus, the membrane has a greater tendency to pull up off the graphite, expanding the annulus width a and decreasing h/a. These trends are illustrated qualitatively in figure 5. For d/a, on the other hand, both d and a increase with decreasing h, and their ratio does not necessarily have a vanishing or asymptotic behavior, as h/a and d/h do, respectively. Based upon the result in figure 4(c), it appears that they increase almost proportionately.
There is little, if any, agreement between these results and the experimental data. In experiment, the ratio h/a does not dive downward at small h, though one could argue that clusters are not observed with sufficiently small h for this trend to be visible. The experimental ratio d/h does rise steeply with decreasing h, as predicted, but at much larger values of h than predicted. Quantitative agreement is thus lacking, even though the trend is reproduced. And for d/a, the experimental ratio rises steeply at small h, in complete disagreement with the predicted trend, both quantitatively and qualitatively. The SLBT model with bending strain only is thus inconsistent with the experimental data at small h.
3.3.3. Results using a model with coupled stretching and bending (s±b) to give U e In order to treat both stretching and bending simultaneously, we have formulated a model in which a general shape is assumed for the graphene membrane. Details are given in the SI. This stands in contrast to the SLBT Figure 5. Qualitative schematic showing the change in island dimensions as bending strain increases. The views are drawn such that a is constant. The cluster height h decreases from A to C, meaning that h/a also decreases. Thicker arrows at the edges of the metal cluster indicate stronger force. model, where the starting point was the governing differential equations, and the shape of the graphene membrane was derived. The SLBT model provides a superior analog to the real physical island, but as noted in section 2.2, it is not solvable simultaneously for bending and stretching. Results of CE modeling, using U e from the s+b model, are shown by the short-dashed lines in figure 4. One key observation is that for small h, these results fall very close to those of SLBT with bending only, both in terms of qualitative trends and quantitative values. This is reasonable, since bending should dominate at small h. The s+b model does not produce better agreement with the experimental data at small h. At large h, it agrees with the SLBT-stretching model, for d/h and d/a. For h/a, it shows an asymptote at large h (as observed) but the asymptotic value is in quantitative disagreement with the SLBT-stretching result. This is likely due to the fact that the s+b model requires the profile of the graphene membrane to be assumed, while the SLBT model derives the profile, ab initio.
We have considered the possibility that increasing D in the bending strain energy term might shift the upturn in d/h toward the right and toward better agreement with experiment. Increasing bending stiffness can be achieved by, for example, doubling the Young's modulus, Y, to 2.2 TPa in equation (8), or by adding a factor of two to the bending term in equation (9). Unfortunately, these have only a minor effect on d/h, and they enhance the downward dive in h/a with decreasing h. Neither of these changes improves agreement with experiment. We have also considered that perhaps smaller islands (smaller h) are associated with larger L C and increased bending stiffness. However, we find that increasing L C does not improve agreement with experiment significantly, as shown in the SI.
We conclude that the behavior of the graphene membrane is dominated by stretching, over the range of experimental observation. Bending strain cannot account for the deviations of d/h and d/a at small h.

Discussion
Insight into the above results is given by analyzing the energy terms in the CE model. Here, we return to the version of the CE model where U e is derived from the free SLBT with stretching only (equation (7)). Figure 6 shows the energy terms for a single choice of V, corresponding to 1.0×10 4 nm 3 . The choice of V does not matter, however, because results are independent of V, as will be discussed below. Figures 6(a), (b) show total energy ( ) P a h , .The minimum in Π defines the equilibrium state. Figures 6(c), (d) represent two 2D cuts across the 3D ( ) P a h , space, each passing through the minimum in P. One cut parallels the a-axis and the other parallels the h-axis. Each cut shows the variation in Π, as well as individual U IS and U e terms. The position of the minimum in Π is determined by those energy terms with highest curvature in the vicinity of the minimum. Hence inspection of figure 6(c) shows that U e and U GnGt determine a eq , while U Fe and U FeG have no effect since they are invariant with a. The former two terms are the elastic energy of the distorted graphene membrane and the adhesion energy between graphene and graphite, while the latter two are the surface and adhesion energies of the metal cluster. Therefore, the equilibrium value of a, a eq , depends solely on properties of graphene and graphite, and has nothing to do with the metal cluster except for the fact that the metal cluster displaces the graphene membrane to a height h. Thus, the ratio h/a is predicted to be independent of the nature of the metal or the lateral dimensions of the metal cluster, and to depend only on the nature of the membrane and substrate.
In accord with this prediction, we have previously presented a comparison of h/a for encapsulated Cu and Fe clusters on graphite, where we have shown that h/a is the same for the two systems within experimental error: 0.24±0.03 and 0.27±0.04, respectively [13]. If the slight difference between the average values is real, it can be attributed to different values of L C in the two systems, L C being slightly higher for Cu than for Fe. Those which provide the best fit to the data are L C =3-4 for Cu and L C =2-3 for Fe (see figure 4(a)).
The above perspective on h/a-that it is simply due to delamination of the top membrane by the metal cluster-leads us to consider the possibility that h/a may be predicted using a simpler framework than that of the full CE model. To this end we examine two established models for delamination of flexible membranes that yield analytic results for h/a. The first assumes a geometry in which a membrane originally lies on a flat substrate in the xy plane. At x=0, the height of the membrane is increased along z up to z=h and the membrane unbinds from the substrate up to a distance x=a from the origin. Essentially, in this model the shaft in the SLBT (which is equivalent to the Fe cluster in figure 1(g)) is replaced by a thin vertical plate parallel to the yz plane, exerting an upward line load on the membrane. This situation has been analyzed by Williams [53], and it can be called the rectangular strip peel test model. If h/a is substantially less than 1, which is true in our situation, h/a is given by where Y′ =Y/(1−ν 2 ). Using equation (12) and the values of variables listed in i.e. h/a is larger by a factor of 1.456 than the value predicted by equation (12). The assumption of a point source is expected to lead to an overestimation of h/a relative to the real geometry. This is due to non-linearity (i.e. apparent 'droop') in the delaminated graphene membrane, which is highly exaggerated when a point load is applied. This model yields h/a=0.287 or 0.259 for L C =2 or 3, respectively. Given the expected overestimation, L C =2 provides the best match to experiment, with an h/a that is only 6% higher than experiment and with a value of L C in agreement with the clamped SLBT model. In short, the point-loaded circular blister test model does a reasonable job of predicting delamination geometry at the edge of an Fe island.
Returning to the analysis of Π, we can similarly identify the energy terms that influence h eq . Here, figure 6(d) shows that all terms contribute, but the only term that increases strongly with increasing h near the minimum in Π is U e . Therefore, U e places the upper limit on island height, even though the value of U e is small at the minimum. This means that the resistance of the graphene membrane to strain strongly inhibits upward growth of the Fe clusters, leading to flat islands with high aspect ratio, d/h. The U Fe and U FeG terms also contribute, with the former having stronger curvature near the minimum. The involvement of these two terms leads to an expectation that the value of h eq predicted by the model will be metal-dependent. We can test this by comparing the limiting (asymptotic) value of d/h in the experimental data for Cu and Fe, where we find values of 7.3±2.8 and 5.3±1.2, respectively. Within experimental error these two values are equal, which is surprising. The slight difference between them, if real, is consistent with a higher value of surface energy γ for Fe than for Cu, driving the Fe clusters to be more compact. The similarity between the two asymptotic values is probably due to a partial cancellation between the increasing total adhesion energy (U FeG ) and the decreasing total surface energy (U Fe ) at the minimum of Π. Both adhesion energy β and surface energy γ are larger for Fe than for Cu, as shown in the SI, but it appears that the corresponding U-terms counteract each other to about the same degree for the two metals, leading to weak or negligible metal-dependence. The interpretation that strain energy in the membrane accounts for the high aspect ratios in the clusters, is borne out by comparison with equilibrium shapes of Fe clusters, which we have analyzed elsewhere. For an unsupported hcp-Fe cluster, d/h=0.70-0.81, depending upon which (0001) facets are chosen in defining h. For an hcp-Fe cluster supported on graphite, d/h=0.75-0.87. For an hcp-Fe cluster sandwiched between graphite and an unstrained graphene membrane, it is 0.96-1.1. Even the largest of these values is a factor of 5 lower than the CE result of d/h=5.3, or (equivalently) than the high-h asymptote in the experimental data. Hence the resistance to strain in the top graphene membrane flattens the Fe cluster by at least a factor of 5.
As we have discussed previously for Cu, CE modeling predicts that all three ratios h/a, d/h, and d/a are sizeindependent [13]. This is because the relative sizes and shapes of the energy graphs in figure 6 are invariant with volume. This, in turn, is rationalized by analysis of the relevant energetics, given by equations (3)-(7) in table 2, which show that each U-term scales as the square of linear dimension. The result is that the predicted profile of a shape-equilibrated encapsulated Fe island is size-invariant, just as the equilibrium crystal shape of a free (or supported) solid crystalline particle is size-invariant above the atomistic limit.
This prediction is borne out in the experimental data for large clusters, but not for small clusters, which show strong variation in d/h and d/a at small h. In experiment, similar trends in d/h and d/a are observed at small h also for encapsulated Cu clusters [13]. In this paper, we have shown that the contribution of bending strain at small h is not a viable explanation for the trends observed in d/h and d/a at small h. A different explanation must be found. We have previously considered that the continuum analysis may break down due to atomistic effects for small metal clusters, but we rejected this possibility because even the 'small' metal clusters in experiment are very large [13]. Two other explanations are currently being explored. One is that deformation of the graphite substrate plays a role. The other is that diffusion-mediated coalescence of small metal islands affects island shapes strongly at small h. There is evidence that such coalescence occurs. These two possibilities are being investigated in our group. For now, explaining the trends in d/h and d/a at small h remains an open challenge.

Conclusions
This paper gives a detailed presentation and analysis of the profiles of encapsulated Fe islands, as measured with STM. This complements a prior presentation and analysis of the shapes of encapsulated Cu islands. Comparison reveals many similarities between the two systems, which indicates that trends observed for the two metals have broad significance. The main trends are that the slope of the annulus h/a is constant, while the aspect ratio d/h of the central cluster falls sharply and then plateaus, as a function of increasing island height h. Quantitative agreement with CE theory is obtained for h/a over the entire range of h, and for d/h in the limit of large h, using values of L C =2-3 for Fe. Input for adhesion and surface energies is provided by DFT. Input for the elastic strain energy comes both from SLBT (which treats stretching and bending independently) and, in this paper, from a different model with coupled stretching+bending. Perhaps the most important conclusions, both for Fe and Cu, are that: (1) the side slope h/a is determined solely by delamination of the graphene membrane from the graphite substrate, and is independent of the metal cluster except in the sense that the cluster displaces the membrane upward by an amount h; (2) the metal clusters are extremely low and flat (they have high aspect ratio d/h), because they are squeezed by the graphene membrane; and (3) the profile of the islands is independent of size, in the limit of large islands.
An important result in this paper is that bending cannot account for the strong deviation of d/h (and d/a) from the CE prediction at low h, as had been suggested previously. In fact, the explanation remains an open challenge. Another interesting result is that the geometry of the annulus-the value of h/a-can be predicted well from an analytic expression derived from a simpler, classic model for debonding of a flexible membrane (the point-loaded circular blister test model), hence circumventing the full CE treatment.