Simple views on surface stress and surface energy concepts

Classical theory of elasticity treats solids as continuum mediums. Let us consider a piece of matter (of volume V limited by a surface S) in an infinite body. This volume is submitted to two forces fields: the forces ${\bf f}^{{\rm{surf}}}$ exerted by the remaining body on the surface boundary of the piece of matter and the external forces ${\bf f}^{{\rm{ext}}}$ originating from an external field such as, for instance, gravity. At equilibrium, these forces are written as $f_i^{{\rm{surf}}} = \sum\nolimits_j {\sigma _{ij} n_j }$ and $f_i^{{\rm{ext}}} = - \sum\nolimits_j {\frac{{\partial \sigma _{ij} }}{{\partial x_j }}}$, where σ_ij are the components of the bulk stress tensor and n is a normal unit vector directed towards the exterior of the surface. The virtual work against an infinitesimal displacement field δu can thus be written $\delta W^{{\rm{bulk}}} = \oint_S {\sum\nolimits_i {f_i^{{\rm{surf}}} \delta u_i \,{\rm{d}}S} } + \int_V {f_i^{{\rm{ext}}} \delta u_i \,} {\rm{d}}V$. Using the above-mentioned equilibrium expressions of f_i^surf and f_i^ext and the divergence theorem to transform the surface integral to an integral on the whole volume, there is [2, 17]

$$\begin{equation} \delta W^{{\rm{bulk}}} = \int {\sum\limits_{i,j} {\sigma _{ij} } } \delta \varepsilon _{ij} \,{\rm{d}}V\quad (i,{\rm{ }}j = 1,{\rm{ }}2,3). \end{equation} \tag{ 1 }$$

This elastic energy is an extensive quantity. As described by Gibbs [18], any extensive quantity may present an excess at the interface between two materials. When one of the materials is in vacuum, this excess quantity is a pure surface-excess and is defined as the difference between the energy of the whole system (bulk + surface) minus the energy of a piece of bulk having the same volume (without surface). In the case under study we can thus define the surface elastic energy as

$$\begin{equation} \delta W^{{\rm{surf}}} = \int {\sum\limits_{ij} {\sigma _{ij} \left( z \right)\delta \varepsilon _{ij} \left( z \right)\,{\rm{d}}V - \sum\limits_{ij} {\sigma _{ij}^{\rm{B}} } \delta \varepsilon _{ij}^{\rm{B}} V} } , \end{equation} \tag{ 2 }$$

where σ_ij (z) and ε_ij (z) are the true stress and strain profiles for the complete system (bulk + surface), whereas σ_ij^B and ε_ij^B are the homogeneous stress and strain tensors in the homogeneous bulk. When considering mechanical equilibrium and non-gliding conditions this equation can be simplified [2]:

$$\begin{equation} \delta W^{{\rm{surf}}} = A\sum\limits_{\alpha \beta } {s_{\alpha \beta } \delta \varepsilon _{\alpha \beta } } \quad (\alpha ,\beta = 1,2), \end{equation} \tag{ 3 }$$

where

$$\begin{equation} s_{\alpha \beta } = \frac{1}{A}\left[ {\int {\sigma _{_{\alpha \beta } } (z)\,{\rm{d}}V - \sigma _{\alpha \beta }^{\rm{B}} V} } \right] \end{equation} \tag{ 4 }$$

defines the quantity s_{α β} as an excess quantity of the bulk stress component tangential to the surface [1, 2]. This quantity is thus called surface stress³. Surface stress components may be positive (tensile) or negative (compressive). Obviously for an isotropic surface s_{α β}  =  sδ_{α β} where δ_{α β} is the Kronecker symbol.

For illustrative purposes, figure 1 shows the calculated energy versus normal deformation ε_zz (figure 1(a)) and versus in-plane deformation ε_xx (figure 1(b)). Three calculations have been performed [2]: for a bulk material, for an un-relaxed slab of n layers then for an elastically relaxed slab. The calculation has been performed using the Wien2k code [20] (ab initio calculations by Wien code for a cubic material). For the slab, two calculations have been performed, the first one is for an un-relaxed slab the second one for a relaxed slab.

Zoom In Zoom Out Reset image size

Figure 1 shows that (i) for such small deformations Hooke's law is valid (quadratic shape of the bulk energy), (ii) for an unrelaxed slab, the surface stress shifts the parabola of the bulk energy and (iii) for an elastically relaxed slab, the bulk parabola is shifted for the in-plane deformation but remains centred on zero for the normal deformation. This last point (iii) simply illustrates that, at mechanical equilibrium, surface stress components s_iz must be zero so that surface stress can be considered as a degenerated two-dimensional second rank tensor (in the referential of the flat surface)

$$\begin{equation} \left[ s \right] = \left( {\begin{array}{@{}cc@{}} {s_{11} } & {s_{12} } \\ {s_{21} } & {s_{22} } \\\end{array}} \right) \cdot \end{equation} \tag{ 5 }$$

For fluids, first and second principles lead to the classical expression of the differential of the internal energy change

$$\begin{equation} {\rm{d}}U = T\,{\rm{d}}S + \sum\limits_i {\mu _i \,{\rm{d}}N_i } - P\,{\rm{d}}V, \end{equation} \tag{ 6 }$$

where S, V, N_i, T, P and μ_i are the usual thermodynamics quantities entropy, volume, number of atoms, temperature, pressure and chemical potential of species i.

The internal energy U(S,N_i ,V) is a homogeneous function of first order so that equation (6) can be integrated by means of the Euler procedure giving thus

$$\begin{equation} U = TS + \sum\limits_i {\mu _i N_i } - PV \cdot \end{equation} \tag{ 7 }$$

For solid systems containing surfaces, the analogue of equations (6) and (7) has to be found [2, 21–24].

Let us consider, at first, the differential of the superficial energy change due to a change of heat, a change of the number of surface particles and a surface area change dA. It is tempting to write

$$\begin{equation} {\rm{d}}U^{{\rm{surf}}} = T\,{\rm{d}}S^{{\rm{surf}}} + \sum\limits_i {\mu _i \,{\rm{d}}N_i^{{\rm{surf}}} } + g\,{\rm{d}}A. \end{equation} \tag{ 8 }$$

However, since the days of Gibbs [18] it is known that the surface area can be changed by two different processes: creation at constant strain or deformation at constant number of surface atoms. The isothermal work to create a surface at constant strain reads dW^cre = γ dA (where the quantity γ is connected to the surface energy density f^surf by f^surf = γ - N_i^surf μ_i /A ), whereas the isothermal work necessary to deform a surface at a constant number of surface atoms obtained in equation (3) can be written ${\rm{d}}W^{{\rm{def}}} = \sum\nolimits_{\alpha ,\beta } {s_{\alpha \beta } } A_0 \,{\rm{d}}\varepsilon _{\alpha \beta }$ where A₀ is the undeformed surface area and dε_{α β} the bulk strain extrapolated to the surface).

Expression (8) is thus difficult to use since g dA has to be identified to γ dA^def or $A_0 \sum\nolimits_{\alpha \beta } {s_{\alpha \beta } } \,{\rm{d}}\varepsilon _{\alpha \beta }$ according to the thermodynamic transformation that is considered. Obviously, one can imagine a process in which there is surface creation and surface deformation so that from a general viewpoint g in (8) is defined by

$$\begin{equation} g\,{\rm{d}}A = \gamma \,{\rm{d}}A^{{\rm{cr}}} + A_0 \sum\limits_{\alpha \beta } {s_{\alpha \beta } \,{\rm{d}}\varepsilon _{\alpha \beta } } \cdot \end{equation} \tag{ 9 }$$

Obviously g dA has no physical meaning and cannot be measured, whereas γ and s_{α β} can be measured by specific experiments.

Equations (8) and (9) could lead us to think that U^surf is a function of the variables S^surf, N_i^surf, A^cre and A^def. However, U^surf (S^surf ,N_i^surf ,A^cre ,A^def ) is not a state function since it depends on the thermodynamic process (creation or deformation). The true state function remains U^surf (S^surf ,N_i^surf ,A). The differential form

$$\begin{eqnarray*}{\rm{d}}U^{{\rm{surf}}} & =& \left. {{{\partial U}/{\partial S^{{\rm{surf}}} }}} \right|_{N_i^{{\rm{surf}}} ,A} \,{\rm{d}}S^{{\rm{surf}}} + \left. {{{\partial U}/{\partial N_i^{{\rm{surf}}} }}} \right|_{S^{{\rm{surf}}} ,A} \\ &&\times{\rm{d}}N_i^{{\rm{surf}}} + \left. {{{\partial U}/A}} \right|_{N_i^{{\rm{surf}}} ,S^{{\rm{surf}}} } \,{\rm{d}}A\end{eqnarray*}$$

cannot be directly compared to the physical variation

$$\begin{eqnarray*}&&{\rm{d}}U^{{\rm{surf}}} = T\,{\rm{d}}S^{{\rm{surf}}} + \mu _i \,{\rm{d}}N_i^{{\rm{surf}}} + \gamma \,{\rm{d}}A^{{\rm{cr}}} + A_0 \sum\limits_{\alpha \beta } {s_{\alpha \beta } \,} {\rm{d}}\varepsilon _{\alpha \beta } . \end{eqnarray*}$$

Let us define U^surf as the surface excess (in the Gibbs meaning) of the bulk internal energy U

$$\begin{equation} U^{{\rm{surf}}} = TS^{{\rm{surf}}} + \mu _i N_i^{{\rm{surf}}} + \gamma A, \end{equation} \tag{ 10 }$$

where we used the well known fact that γ A is the excess of the grand potential - PV.

It is now possible to differentiate equation (10) with dA  = dA^def + dA^cr. Using (8) and (9) gives the Gibbs Duhem equation valid at the surface [2, 21]:

$$\begin{equation} A\,{\rm{d}}\gamma = - S^{{\rm{surf}}} \,{\rm{d}}T + N_i \,{\rm{d}}\mu _i + A_0 \sum\limits_{\alpha \beta } {\left( {s_{\alpha \beta } - \gamma \delta _{\alpha \beta } } \right)\,} {\rm{d}}\varepsilon _{\alpha \beta } , \end{equation} \tag{ 11 }$$

where δ_{α β} is the Kronecker symbol. From equation (11) can be extracted the so-called Shuttleworth relation [25] (here written in Euler coordinates)

$$\begin{equation} s_{\alpha \beta } = \gamma \delta _{\alpha \beta } + \left. {{{\partial \gamma }/{\partial \varepsilon _{\alpha \beta } }}} \right|_{T,\mu _i } . \end{equation} \tag{ 12 }$$

For incompressible liquids surface energy and surface stress are equal and are then usually merged under the single term of 'surface tension'.

For completeness, notice that one can write equation (10) under the form U^surf =  TS^surf + μ_i N_i^surf + γ (ε )A(ε ) where we have explicitly written the strain dependence of γ and A. Then equation (10) can be expanded up to the first order in strain so that using $\gamma (\varepsilon ) = \gamma _0 + \varepsilon \left. {{{\partial \gamma }/{\partial \varepsilon }}} \right|_{\varepsilon = 0}$ and the surface stress definition s  = (γ (ε )A(ε ) - γ₀ A₀ )/ε A₀ (here assumed to be isotropic s_{α β}  =  sδ_{α β} ) one obtains U^surf ≈ TS^surf + μ_i N_i^surf + γ₀ A^cr +  sA₀ ε in agreement with Hermann structure of thermodynamics (for a discussion see [26–29]). 

Let us underline the difference of nature between surface stress and surface energy.

The surface energy γ can be associated with the energy needed to break bonds. In reality there is an electronic and a structural surface relaxation that we will neglect in this section so that we will suppose that the surface energy can be calculated, at least at zero kelvin, by counting the number of broken bonds and thus by adding interatomic potentials [30]. The surface energy thus is a scalar that scales as an energy per unit area.

The surface stress, which is a tensor, is associated with the work against surface deformation, it is the result of forces acting at the material surface and thus calculated as a sum of the first derivatives of the interatomic potentials (figures 2 and 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As an illustration, let us consider a solid whose total surface energy can be written as a sum of pair potentials. In this case the surface energy γ can be expressed as $\gamma = \frac{1}{{2a^2 }}\sum\nolimits_j {m_{\!j} } \Phi _{\!j} (a_j )$ the surface energy expressed from broken bonds where m_j is the number of broken jth neighbours and Φ _j (a_j ) the interaction energy of atoms at a distance a_j. In the following we will write a_j =  K_j a where a is an atomic distance. For the sake of simplicity we only consider an isotropic surface stress so that the Shuttleworth relation reads s  = γ + ∂ γ /∂ ε with ε = dA/A . We thus obtain $s = \frac{1}{{2a}}\sum\nolimits_j {K_j } \Phi '(a_j )$ where Φ ' is the first derivative of Φ,⁴. The main difference between γ (a function of interatomic bonds) and s (a function of the derivatives of interatomic bonds) is illustrated in figure 3.

In other words, if the surface is divided into parts by a curved boundary, where n is a unit vector normal to the boundary in the surface plane, the surface stress tensor gives the force $f_i = \sum\nolimits_j {s_{ij} } n_j$ exerted across the boundary curve (see figure 2).

From a physical point of view, the surface energy plays a role in the processes where a surface is created, whereas surface stress plays a role each time a surface is deformed. It is illustrated in table 1 where we report, in the simplified case of an isotropic sphere of radius R, the excess of chemical potential Δ μ and the excess of pressure Δ P induced in the bulk by its surrounding surface [28, 32].

Table 1. Basic thermodynamics equations valid for isotropic solid and liquid spheres of radius R illustrating the deep difference between surface energy and surface stress.

	Excess of pressure	Excess of chemical potential
Solid	$\sigma _{33} = \frac{{2s}}{R}$	$\Delta \mu = \frac{{2\gamma }}{R}$
Liquid	$\Delta P = \frac{{2\gamma }}{R}$	$\Delta P = \frac{{2\gamma }}{R}$

For a solid, the normal component of the bulk stress σ₃₃ is due to the surface stress that exerts a force on the underlying volume phase, whereas the chemical potential change, connected to a variation of particle number, is due to the surface energy. For a liquid, that cannot be stretched at constant number of surface atoms, one cannot distinguish s from γ and Δ μ = Δ P  = 2γ /R .

The deep difference of nature between surface stress and surface energy can also be illustrated by comparing surface energy and surface stress anisotropy. Indeed, a surface is characterized by its surface energy γ (a scalar) and by the two principal components of its surface stress (a tensor). In figure 4 we show surface energy and surface stress values calculated for various surface orientations [2]. The γ values are distributed on a single branch, whereas surface stress values are distributed on two branches corresponding to the two principal components of the surface stress. Moreover, it can be seen that local minima in the surface energy correspond to maxima in the surface stress. It is quite normal since close to low index orientations there is an energy cost to form steps that increase the surface energy, whereas the steps allow relaxation and thus reduce the surface stress.

Zoom In Zoom Out Reset image size

Vicinal surfaces can be described as stairs-like surfaces, where monatomic steps separate two microscopically flat terraces characterized by their own surface stress [2–4, 33]. Since the atoms belonging to the step edges have a different number of nearest neighbours than the atoms in the underlying bulk, steps give rise to a lattice distortion in the underlying bulk. In the simplest scheme, this elastic field may be described replacing the steps by a distribution of point forces applied on a planar surface [2, 34, 35]. In this case, mechanical equilibrium conditions naturally lead to distinguishing two cases: (i) if the two neighbouring terraces separated by one step have the same surface stress, the step can be described as an elastic dipole whose normal component A_⊥^dip is simply the surface stress value multiplied by the step height. (ii) If there is a surface stress discontinuity at the step, the step has to be described as an elastic monopole whose tangential component A_∥ ^mono is proportional to the surface stress difference in between both terraces.

In figure 5(a) is reported the case of a Si(111) surface whose (111) terraces are characterized by an isotropic surface stress. From an elastic viewpoint, steps can be thus modelled by rows of elastic dipoles distributed along the step edge.

Zoom In Zoom Out Reset image size

In figure 5(b) we show the case of a Si(001) surface misoriented in the [110] direction. Due to the diamond structure of silicon, one terrace reconstructs with dimers parallel to the $[\bar 110]$ direction, whereas the other terrace reconstructs with dimers parallel to the [110] direction. The surface stress of the (001) terraces is thus a second rank tensor which reads $({{s_{xx}}\atop{0}} {{0}\atop{s_{yy}}})$ for one terrace and $({{s_{yy}}\atop{0}} {{0}\atop{s_{xx}}})$ for the other (when written in the [110], $[\bar 110]$ surface axis). As a consequence, a vicinal surface misoriented in the [110] direction is characterized by a surface stress difference $\pm | {s_{xx} - s_{yy} } |$ in the direction normal to the step. It is then characterized by an elastic monopole.

A facetted surface can be described by its mesoscopic facets (black and green facets in figure 6). Each facet can be characterized by its own surface stress. The facetted surface can thus be described as an assembly of L-apart elastic monopoles located at the facet boundaries [36, 37]. The surface stress discontinuities at the facet boundaries can be modelled by rows of elastic monopoles perpendicular to the discontinuities. At equilibrium the elastic relaxation (energy gain) induced by these surface forces balances the energy cost (β) due to the domain boundaries. It results that the period of the facetted structure is fixed by this equilibrium. In the simplest case depicted in figure 6 this period reads [36, 37]

$$\begin{equation} \lambda = \frac{{\pi a_0 e}}{{\sin (\pi p)}}{\rm exp}\,\left( {\frac{{\pi E\beta }}{{2(1 - \nu ^2 )\left| {{\bf s}_1 - {\bf s}_2 } \right|^2 }}} \right), \end{equation} \tag{ 13 }$$

where p  =  N/M is the relative coverage of one phase with respect to the other, E and ν the Young modulus and the Poisson ratio of the solid, a₀ a cut-off introduced to avoid the local divergences due to the use of point forces. The surface stress components which appear in equation (13) are the surface stress components perpendicular to the edge which separates the α and β facets.

Zoom In Zoom Out Reset image size

The same facetted surface may also be described as an assembly of L-apart bunches of monoatomic steps characterized by elastic dipoles located at the steps edge (on the green facet). In this case the period reads [15]

$$\begin{equation} \lambda = \frac{{\pi a}}{{\sin \left( {\pi p} \right)}}{\rm exp}\,\left( {\frac{{\beta a^2 }}{A}} \right), \end{equation} \tag{ 14 }$$

where $A = 2\frac{{1 - \nu ^2 }}{{\pi E}}( {( {A_\parallel ^{{\rm{dip}}} })^2 + \left( {A_\parallel ^2 } \right)} )$ with A_⊥^dip and A_∥ ^dip the normal and tangential components of the elastic dipole.

The nanoscale and microscale models are equivalent if a₀ e  =  a and if $( {A_\parallel ^{{\rm{dip}}} })^2 + \left( {A_\parallel ^2 } \right) = ( {{\bf s}_1 - {\bf s}_2 } )^2$, introducing the step height h, so that $a = {h/{\tan \theta }}$ (where θ is the angle of the vicinal facet) and A_⊥^dip =  hs₁ where s₁ =  s_xx is the surface stress component perpendicular to the step. For weak values of θ one obtains [15]

$$\begin{equation} s_2 = s_1 - \frac{{\left| \theta \right|}}{h}A_\parallel ^{{\rm{dip}}} . \end{equation} \tag{ 15 }$$

This expression is analogous to the one found by Salanon for stressed solids [38] where the steps are described by the sum of rows of dipoles and monopoles (both perpendicular to the step) and the surface stress expression is developed up to second order in θ.

Equation (7) means that since the presence of steps leads to surface stress relaxation, the surface stress is maximum for a low index surface and thus decreases with $\left| \theta \right|$. In contrast, the energy cost to create surface steps implies that the surface energy increases with $\left| \theta \right|$. In other words, local minima (cusps) of the surface energy plot (gamma-plot) correspond to local maxima (anticusps) of the surface stress plots.

Finally note that the elastic description of the surface stress discontinuity at the step has been widely checked from grazing incidence x-ray diffraction (GIXD) [39, 40] used to measure the elastic displacement field created by steps.

In the classical Marchenko approach used to calculate the bulk deformation induced by steps, the steps are replaced by elastic forces only characterized by components normal and parallel to the average surface [36, 37]. A better description of elastic displacements and energy is obtained when taking into account the crystalline anisotropy and describing steps as 'buried dipoles', i.e. dipoles with a lever arm not necessary parallel to the surface [41].

In figure 7 is reported the surface free energy anisotropy obtained by Wulf inversion of the experimental equilibrium shape (T = 1373 K) [14, 42]. This method does not allow obtaining the absolute values of the free surface energies but only the relative ones. However, the absolute surface energy of the Si(111) surface at 1373 K has been obtained independently [43].

Zoom In Zoom Out Reset image size

The surface stress anisotropy of silicon is also reported in figure 7. To obtain the surface stress anisotropy plot, electromigration has been used to facet Si surfaces that belong to the equilibrium shape [15]. In this work [15] the authors have assumed that, though electromigration was the driving force for faceting, the period of faceting is still fixed by elasticity.

In this case equation (13) has to be generalized to take into account the new geometry [2] (see figure 8)

$$\begin{eqnarray*}&&\lambda = \frac{{2\pi c}}{{\sin \pi \delta }}{\rm exp}\,\left( {1 + \frac{{\pi E\beta }}{{2f^2 \left( {1 - \nu ^2 } \right)}}} \right) \end{eqnarray*}$$

with

$$\begin{equation} f^2 = \left( {s_\tau ^\alpha - s_\tau ^\beta } \right)^2 + 4s_\tau ^\alpha s_\tau ^\beta \sin ^2 \left( {{\textstyle{{\theta _1 + \theta _2 } \over 2}}} \right), \end{equation} \tag{ 16 }$$

θ₁ and θ₂ the angles the facets form with the original orientation, $\delta = \frac{{{\rm{tg}}\,\theta _1 }}{{{\rm{tg}}\,\theta _1 + {\rm{tg}}\,\theta _2 }}$ a geometrical factor and c an atomic unit. The surface stress components s_τ ⁱ which appear in (16) are the surface stress components perpendicular to the edge τ which separates the facets α and β.

Zoom In Zoom Out Reset image size

Measuring the period of faceting thus gives access to the surface stress difference between the formed facets. In our experiments we have measured the period of faceting for various initial surface orientations. Considering a whole set of vicinal surfaces, it is thus formally possible to obtain surface stress differences between many facets and thus the whole surface stress anisotropy. More precisely in [13], the surface stress differences have been calculated by using equation (13) with a common edge energy β (since all the facets we have studied have the same common zone axis) and by calculating the bulk elastic properties (E and ν) for each crystallographic orientation. The so-obtained data lead to a set of equations whose mathematical solution should give access to the surface stress polar dependence. However, only a few of the many numerical solutions have a physical meaning. Taking into account what precedes, these solutions are the positive solutions that lead to maxima of surface stress at minima of surface energy. We find a single solution [13]. We compare in figure 7 the surface energy-plot to the surface stress. We also report their polar-plot in the (110) plane. Note that the surface stress anisotropy is much more important than the surface energy anisotropy in agreement with what we calculated for various metals in [2]. This behaviour is quite normal since it is well known that surface stress is much more sensitive to surface relaxation than surface energy.

As above-mentioned, the amplitude of the elastic forces exerted by the steps on their underlying substrate can be obtained from grazing incidence x-ray diffraction (GIXD) measurements [39]. Many results have been published for metals [39–41]. For silicon, to the best or our knowledge, only the Si(7710) surface, which is a vicinal of the Si(111) 10° misoriented towards the [1–12] direction, has been studied [16]. This surface exhibits triple steps that induce an elastic field in their underlying substrate. This elastic field has been recently studied by GIXD. The strain field induced by the triple steps is well described by an elastic model in which the steps are modelled by parallel rows of buried elastic dipoles. The best fit of the experimental results is reached on the basis of the (7710) surface structure model proposed by Teys et al [44].

The value of the torque component P_T of the dipole is known to be directly related to the surface stress by sh where s is the isotropic surface stress of the nominal surface and h the step height. Since the Si(7710) surface is formed by triple steps one should expect P_T = 3sh where h is the height of a monoatomic step. The experimental value extracted from GIXD experiments is P_T = 1.5 nN [15]. This value is not three times the value expected for a monoatomic step: P_T^mono =  sh  ≈  0.9 θ₁ nN [15]. It is quite normal since a triple step on a (7710) surface is not a simple stack of three monoatomic steps but appears to be quite defective [15, 44].