Strain distribution due to surface domains: a self-consistent approach with respect to surface elasticity

Elastically mediated interactions between surface domains are classically described in terms of point forces. Such point forces lead to local strain divergences that are usually avoided by introducing a poorly defined cut-off length. In this work, we develop a self-consistent approach in which the strain field induced by the surface domains is expressed as the solution of an integral equation that contains surface elastic constants, Sij. For surfaces with positive Sij the new approach avoids the introduction of a cut-off length. The classical and the new approaches are compared in case of 1-D periodic ribbons.


Introduction
The classical approach used to calculate the strain field that surface domains induce in their underlying substrate consists of modeling the surface by a distribution of point forces concentrated at the domain boundaries [1][2][3], the force amplitude being proportional to the difference of surface stress between the surface domains [3][4][5][6]. However, point forces induce local strain divergences, which are avoided by the introduction of an atomic cut-off length. Hu [7,8] stated that the concept of concentrated forces is only an approximation valid for infinite stiff substrates. Indeed if the substrate becomes deformed by the point forces acting at its surface, the substrate in turn deforms the surface and then leads to a new distribution of surface forces so that the surface forces have to be determined by a self-consistent analysis. In this paper, we show that when elastic surface properties are properly considered, the strain field induced by the surface domains may be expressed as the solution of a selfconsistent integro-differential equation.

Results and Discussion
Let us consider (see Figure 1a) a semi-infinite body whose surface contains two domains (two infinite ribbons) A and B characterized by their own surface stress s A and s B . The 1D domain boundary is located at x o = 0. Note that for the sake of simplicity only the surface stress components are taken to be different from zero (see Appendix I for the Voigt notation of tensors). In the classical approach [6][7][8] the strain field generated in the substrate is assumed to be generated by a line of point forces (with δ(x) being the Dirac function) and is given by: (1) where D xx (x/x',z) is the xx component of the Green tensor and where the component f x (x) = Δs 1 originates from the surfacestress difference at the boundary between the two surface domains. The Green tensor valid for a semi-infinite isotropic substrate can be found in many text books [1,2,9] so that the deformation at the surface ε 1 (x,z = 0) finally reads: (2) where with E subs and ν subs being the Young modulus and the Poisson coefficient of the substrate (supposed to be cubic). The strain at the surface (Equation 2) exhibits a local divergence at the boundary x = x 0 = 0. The elastic energy can thus be calculated after introduction of an atomic cut-off length to avoid this local divergence [6,10].
However, the concept of point forces is only an approximation. If the substrate is deformed by point forces acting at its surface, the substrate in turn deforms the surface and then leads to a new distribution of surface forces. In the following, we consider that, due to the elastic relaxation, the surface stress at equilibrium exhibits a Hooke's-law-like behavior along the surface [9,11,12]: with i = A, B according to whether x lies in region A or B. In Equation 3, is the surface stress far from the domain boundary (or in other words the surface stress before elastic relaxation) and the surface elastic constants properly defined in terms of excess quantities (see Appendix). The surface force distribution due to the surface stress variation (see Figure 1b) is obtained from force balance and reads f x (x,z = 0) = ds 1 /dx. By using the Green formalism again, we obtain at the surface, z = 0:

This equation replaces the classical result of Equation 2. Equation 4
is an integro-differential equation that has to be solved numerically. At mechanical equilibrium the absence of surface stress discontinuity at the domain boundary, combined to the constitutive Equation 3 leads to the following boundary condition . We also plot in Figure 2a the classical result calculated from Equation 2 (continuous red curve). It is clearly seen that the new expression avoids the local strain divergence that is now replaced by a local strain jump Δs 1 /S 11 at x 0 = 0.
Since the solutions of Equation 4 depend on the values of hS 11 and Δs 1 we report in Figure 2b the results obtained for different typical values of hS 11 and Δs 1 data obtained from [11]. More precisely, since the classical expression scales as 1/x, we plot ln ε versus x. As can be seen, in the limit of large x all solutions tends towards the classical one (common red asymptote in Figure 2b). Moreover we can clearly see that the classical approach is recovered in the limit S 11 → 0. The elemental solution of Equation 4 enables to describe more complex experimental configurations as the one that corresponds to the spontaneous formation of 1D periodic stripes by a foreign gas adsorbed on a surface (as for instance O/Cu(110) [13]). In the classical model each stripe (width 2d) is modeled by two lines of point forces one located at d and the other at −d with the opposite sign f x (x) = Δs 1 (δ(x − d) − δ(x + d))) so that for a set of periodic ribbons of the period L the elastic field is obtained by a simple superposition of the elemental solutions given in Equation 2. In the classical case it reads (6) whereas within the new approach the elastic field is solution of the integral equation: The results are shown in Figure 3 in which two cases are reported. In the first case d/L = 1/2, whereas in the second case d/L = 3/10. Again both solutions (classical and new approach) are quite similar since the only difference lies in the local divergences of the classical model (red curves in Figure 3) that are now replaced by local strain jumps. does not present stable solutions. It is quite normal since in this case, the surface is no more stable by itself but is only stabilized by its underlying layers (see Appendix I). From a physical point of view it means that, for mechanical reasons, we have to consider a "thick surface" or, in other terms, that the surface has to be modeled as a thin film the thickness a of which corresponds to the smaller substrate thickness necessary to stabilize the body (bulk + surface). It can be shown that this is equivalent to modify the integro-differential equation for S 11 < 0, by changing the kernel: (8) In Figure 4 we show the result obtained from numerical integration of Equation 8 for the test value hS 11 = −0.01. In this case a = |2hS 11 | is the minimum value necessary to stabilize Equation 8. Since s 1 is positive but S 11 is negative, there is a sign inversion of ε close to the boundary. For vanishing a this local oscillation propagates on the surface and is at the origin of the instabilities that do not allow to find stable solutions to Equation 4. However we cannot exclude that the total energy of materials with s 1 S 11 < 0 could be reduced by some local morphological modifications of their surface. In such a case, the Green tensor used for this calculation should be inadequate. In conclusion, the self-consistent approach expressed in terms of surface elastic constants is more satisfactory than the classical approach, particularly in the case of stable surfaces (characterized by positive surface elastic constants) for which there is no need to introduce a cut-off length. In case of unstable surfaces (negative surface elastic constants) a cut-off length is still necessary, its value is connected to the minimum substrate thickness necessary to stabilize the body (surface + underlying bulk). Even if the model only deals with 1D structures it can be generalized to other structures such as 2D circular domains. The so-obtained equations are less tractable but the main result remains the same (see Appendix II).

Appendix I: Surface elasticity
From a thermodynamic point of view all extensive quantities may present an excess at the interface between two media (for a review see [9]). For a system formed by a body facing vacuum the following excess quantities can be defined [9]: where (9) is the second order strain development of the energy of a body of volume V 0 limited by a surface of area A 0 and (10) is the second order development of a piece of body of same volume V 0 but without any surface. In these expressions are the bulk stress components and C ijkl the bulk elastic constant.
The so-defined surface quantities depend on a typical length scale at which surface effects are disentangled from bulk effects. Actually, in surface energy calculations, this length is unambiguoulsy determined by a Gibbs dividing surface construction [14]. Surface stress and surface elastic constants values can thus be calculated from strain derivatives of the welldefined surface energy quantity [11].
In contrast to surface energy density and bulk elastic constants, surface stress components and surface elastic constants do not need to be positive. [9,11]. This does not violate the thermodynamical stability condition since actually a surface can only exist when it is supported by a bulk material. Hence the stability of the solid is ensured only by the total energy (surface + volume).
Finally, in the body of the paper we use the Voigt notation so that the surface stress can be written as the components of a 3D vector s = (s xx ,s yy ,s xy ) = (s 1 ,s 2 ,s 6 ), while surface and bulk elastic constants are written as the components of 3D matrices S ij and C ij , respectively.

Appendix II: 2D circular domains
In case of a circular domain of radius R, the classical approach considers a force distribution f r (r) = Δs 0 δ(r−R) that generates a displacement field expressed in terms of complete elliptic integrals K(x) and E(x) as: (11) In the distributed force model, we use the stress-strain relations valid at the surface expressed in polar coordinates: (12) (13) again with the Voigt notation in polar coordinates A rr ≡A r , A θθ ≡A θ .
By using the classical mechanical equilibrium equation and strain-displacement relations expressed in polar coordinates we obtain the following force distribution (14) The displacement can thus be obtained from the self-consistent equation (which replaces Equation 11) The necessary boundary conditions, analog to Equation 5, must now be written for normal and tangential strains (16) (17) The integral equation for the displacement field, Equation 15, only needs the surface elastic constant S 11 , but the edge condition introduces the need of the other surface elastic constant S 12 . Qualitatively the result is similar to the one shown in Figure 2.