Dynamics of Small Solid Particles on Substrates of Arbitrary Topography
Quan Zhao Wei Jiang Yan Wang David J. Srolovitz Tiezheng Qian Weizhu Bao




1]Quan Zhao
[1]School of Mathematical Sciences, University of Science and Technology of China,  Hefei, China


2]Wei Jiang8
[2]School of Mathematics and Statistics,  Wuhan University,  Wuhan, China
jiangwei1007@whu.edu.cn


[8]Corresponding author.

3]Yan Wang
[3] School of Mathematics and Statistics, and Key Lab NAA–MOE Central China Normal University, Wuhan, China


4]David J. Srolovitz
[4]Department of Mechanical Engineering, The University of Hong Kong,  Hong Kong SAR, China


5]Tiezheng Qian
[5]Department of Mathematics, Hong Kong University of Science and Technology, Hong Kong, China


6]Weizhu Bao
[6]Department of Mathematics, National University of Singapore, Singapore



§ ABSTRACT 
We study the dynamics of a small solid particle arising from the dewetting of a thin film on a curved substrate driven by capillarity, where mass transport is controlled by surface diffusion.
We consider the case when the size of the deposited particle is much smaller than the local radius of curvature of  the substrate surface.
The application of the Onsager variational principle leads to a reduced-order model for the dynamic behaviour of  particles on arbitrarily curved substrates.
We demonstrate that  particles move toward region of the substrate surface with lower mean curvature with a determined  velocity.
In particular, the velocity is proportional to the substrate curvature gradient and inversely proportional to the size of the particle, with a coefficient that depends on material properties that include the surface energy, surface diffusivity, density, and Young's (wetting) angle.
The reduced model is  validated by comparing with numerical results for the full, sharp-interface model in both two  and three dimensions.




Surface diffusion, Onsager principle, Solid-state dewetting, Substrate curvature gradient, Wasserstein distance






§ INTRODUCTION


The large and increasing diversity of thin film-based technological applications has led to growing interest in how thin films dewet from or islands forms on substrates (e.g., see <cit.>). Initially continuous thin films are often unstable and  dewet or agglomerate to form isolated particles/island in order to minimize the total interfacial energies.
These observed morphology changes are mainly due to capillary effects, which most commonly occurs via diffusional mass transport along surface <cit.>.
In the case of isotropic surface energy, the normal velocity of the evolving interface is proportional to the surface Laplacian of the local mean curvature <cit.>.
The dynamics of the contact line, where the thin film/vapor interface meets the substrate, is an additional, important  kinetic feature in the evolution of the morphology.
In particular, at the contact line, the force or line tension balance along the substrate implies an equilibrium contact angle (i.e., Young's law <cit.>).


A growing body  of experimental and theoretical efforts has focused on the solid state dewetting mechanisms (see e.g., <cit.>).
In general, the film surface morphology evolution is influenced by many parameters including  film thickness,   film microstructure, surface tension anisotropy, relative surface and interface energies, stress in the film, and the elastic constants of the film and substrate.
More recently, templated dewetting has drawn significant attention; this refers to patterning the shape of the substrate or the film to produce the desired dewetted island microstructures.
Indeed, experiments have shown how topographically patterned substrates can be employed to produce particles of near-uniform size in patterned arrays <cit.>. Moreover, recent experiments on deposited thin films of platinum on a sinusoidally modulated alumina substrate  have provided the clear evidence about the migration of particles from convex to concave sites of the substrates <cit.>. While most theoretical and simulation studies focus on  flat substrates, relatively little attention has been devoted to dewetting on not-flat substrates (despite their growing interest in experimental studies).
In <cit.>, a continuum simulation approach was used to study dewetting on a sinusoidal substrate leading to a period array of islands,  and some related research works can also be found in <cit.>. A mathematical understanding of island dynamics on topographically patterned substrates would be valuable  to guide the  precise control of the dewetting process to produce desired self-assembled islands through templating <cit.>.

Recent years have seen a great deal of experimental and theoretical activity focused on the dynamics of liquid droplets on curved substrates, see e.g., <cit.>.
In such a case, the driving force for evolution is related to the local substrate curvature gradient which drives rapid droplet motion.
Chen and Xu <cit.> derived a reduced-order model for this type of dynamical system and demonstrated quantitative agreement with experiment.
Their  model was based on the Onsager variational principle <cit.>, which provides a  useful approximate framework for describing irreversible thermodynamic processes.
This approach is based on the minimization of the Rayleighian, which has contributions from  the free energy change rate and dissipation function which are described in terms of a small number of suitable state variables.
This  leads to a set of ordinary differential equations for these state variables and forms a reduced-order model (see below Section <ref> for a brief description of this  variational principle).


In this paper, we apply the Onsager principle to provide a reduced-order model for the motion  of solid particles on curved substrates, where the dynamics is controlled via surface diffusion (unlike fluid droplets evolving under viscous momentum transport).
For ease of analysis, we  assume that the interface energies are isotropic, the elastic effects of the thin films are negligible and no chemical reactions or phase transformations occur.
We focus on the case of a small particle on a topographically-patterned substrate, in which the length scale of the particle is much smaller than the radius of curvature of any point on the substrate surface.
We show how the total free energy of the system can be related to the substrate curvature  <cit.> and show how the dissipation function is  uniquely determined by the surface normal velocity in two dimensions (2D).
In three dimensions (3D), we develop an alternative approach to compute the dissipation function by connecting it with the Wasserstein metric <cit.>, which leads to a constrained minimization problem.
We demonstrate that in both 2D and 3D that the new derived model is both simple and quantitative  via several examples.

The remainder of the paper is organized as follows.
In Section <ref>, we introduce the full sharp-interface model for the dewetting of solid thin films on a curved substrate and include a brief review of the Onsager variational principle.
Next, in Section <ref>, we apply this principal to derive a reduced-order model for the dynamics of a 2D particle and validate the model by the numerical results from solving the full model.
In Section <ref>, we generalize this approach to three dimensions.  In Section <ref>, we provide some generalizations of the proposed approach to the cases
of chemically inhomogeneous flat substrates and anisotropic solid particles.
Finally, we draw some conclusions in Section <ref>.






§ PREMINARY




In this section, we first introduce the full sharp-interface model for the dewetting system and then give a short review of the Onsager principle.


 §.§ The full sharp-interface model




Consier a solid thin film deposited on a curved substrate, as shown in Fig. <ref>.
The film/vapor interface is represented by an open hypersurface (t) in d=2 or 3 dimensions,  d∈{2,3}.
We assume  a parameterization of the film surface (t) on the reference domain 𝒪⊂^d-1 (1 dimension lower than the space)  given by

    X⃗(ρ⃗, t): 𝒪×[0,T]↦^d.

The velocity of the interface (t)  is

    V⃗(X⃗(ρ⃗,t), t) = ∂_tX⃗(ρ⃗,t)    X⃗∈(t).

Let n⃗ and ℋ be the unit normal and mean curvature of the interface (t), respectively.
Assuming isotropic diffusion and surface tension, the interface dynamics can be expressed as <cit.>


    V_n    =  -Ω_0∇_s·j⃗,
    j⃗   = -D_sν/k_b T∇_sμ,
    μ   = Ω_0γ_0ℋ,


where V_n=V⃗·n⃗ is the velocity of the  interface in the direction of n⃗, j⃗ is the flux of the surface atoms, μ is the chemical potential of a film atom on the surface and ∇_s is the surface gradient operator.
The physical parameters/constants are the volume per atom of the film material Ω_0 , D_s is the surface diffusivity, k_b T is the thermal energy,
ν is the number of diffusing atoms per unit area (in the direction normal to the surface flux vector), and γ_0 is the isotropic surface energy density.
Combining terms,  (<ref>) can be rewritten as

    V_n = Bγ_0 Δ_sℋ,

where B=D_sνΩ_0^2/k_b T is a material constant and Δ_s is the Laplace-Beltrami (Laplacian) operator.

At the contact line  Γ(t) where the film/varpor interface meets the substrate, we impose the following boundary conditions:



  (i)attachment condition

    V⃗·n⃗_w = 0;


  (ii)contact angle condition

    n⃗·n⃗_w + cosθ_i=0;


  (iii)zero-flux condition

    n⃗_c·j⃗ = 0.



Here, n⃗_w is the substrate unit normal (positive pointing towards the substrate interior) and n⃗_c is the conormal vector of Γ(t), as shown in Fig. <ref>.
Moreover, θ_i is the equilibrium, isotropic Young angle, which satisfies

    cosθ_i=(γ__-γ__)/ γ_0,

with γ__ and γ__ representing the varpor/substrate and film/substrate surface energy densities, respectively. Condition (ii) can be interpreted as the contact angle condition; this  leads to the Young angle  θ_i between n⃗ and -n⃗_w at the contact line.
We note that under some conditions, the dynamic contact angle may differ from the Young's  contact angle condition (e.g., see <cit.>).

The total free energy of the dynamic system is given by

    W(t) = γ_0|(t)|-γ_0cosθ_i A_ sub(Γ(t)),

where |(t)| is the surface area of (t), and A_  sub(Γ(t)) represents the substrate surface area covered by the film/island (i.e., enclosed by Γ(t)).
The free energy of the evolving film/substrate system satisfies

    W(t) = -k_b T/D_s ν∫_(t)|j⃗|^2  S≤ 0.

We also assume that the volume of the film material is conserved (we do not consider phase transformations or strains); i.e., (Ω(t)) = (Ω(0)) for all time t≥0.








 §.§ Onsager variational principle

The Onsager variational principle was first formulated <cit.> based on the reciprocal symmetry in a linear irreversible thermodynamic process.
This fundamental principle provides a general framework to describe non-equilibrium kinetics in cases where linear response  is applicable and has found  wide applications in fluid dynamics <cit.>,  soft matter physics <cit.> and solid-state deweting <cit.>. We first provide a short review of this principle.

Consider an isothermal system described by a set of time-dependent state variables

    β(t) = (β_1(t),β_2(t),…, β_n(t)),

and  β̇(t)=(β̇_1(t),β̇_2(t),…,β̇_n(t)) are the rates of change of these state variables (the raised dot “·” denotes a time derivative).
We further introduce W(β) as the total free energy of the system, then
the rates {β̇_i} are  determined by miniminzing the Rayleighian <cit.>

    ℛ(β̇,β)=Ẇ(β, β̇)+Φ(β̇,β̇),

where Ẇ(β, β̇)=∑_i=1^n(∂ W / ∂β_i)β̇_̇i̇ is the rate of change of the total free energy W
and Φ(β̇,β̇) is the dissipation function.


In the linear response regime, the dissipation function is a quadratic function of the rates {β̇_i}

    Φ(β̇,β̇)=1/2∑_i=1^n∑_j=1^nλ_ij(β)β̇_iβ̇_j,

where the friction coefficients {λ_ij} form a positive definite, symmetric  matrix.
Minimizing the Rayleighian  (<ref>) with respect to  rates {β̇_i}  yields the kinetic equations

    -∂ W/∂β_i=∑_j=1^nλ_ijβ̇_j,     i=1,2,…,n,

which precisely gives the force balance between the reversible force -∂ W/ ∂β_i and the dissipative force -∂Φ / ∂β̇_i. Multiplying (<ref>) by β_i and summing (and  recalling (<ref>))  yields

    Φ(β̇,β̇)=-1/2Ẇ(β,β̇).

This means that the dissipation function is half the rate of the free energy dissipation.
Physically, the variational principle  for isothermal systems can  also be derived from the maximization of the Onsager-Machlup action that is used for more general non-isothermal systems <cit.>.


In the present work, we apply the Onsager variational principle to describe the dynamics
of a small particle migrating on a curved substrate.
We  focus on  cases in which the full model, introduced in <ref>, can be approximately described by a finite set of suitable state variables.
Application of the Onsager principle  enables us to obtain a reduced model for a continuous dissipative system which is governed by a set of ordinary differential equations for a few state variables.
Note that the full model satisfies the energy dissipation law in (<ref>) and  can be obtained as well by applying the Onsager principle in the linear response regime,
where the constitutive equation (<ref>) is derived for the flux j⃗.




§ DYNAMICS IN TWO DIMENSIONS


We now apply the Onsager variational principle in <ref> to derive a reduced-order model for the dynamics of a particle on a substrate in 2D.
The reduced model is then numerically validated by comparisons  with the full model in <ref>.



 §.§ A reduced-order model


We assume that the particle/island  is much smaller than the radius of curvature  of the substrate.
The separation of these two length scales leads to two distinct time scales:
one for the island to establish the circular shape and the other for the island to migrate along the substrate.
In particular, the latter is much slower than the former.
Therefore, it is reasonable to assume that the particle will remain circular during its (relatively) slow motion along the substrate.

As shown in Fig. <ref>(a), we assume that the film/vapor interface of the particle is given by section of a  (small) circle of radius r(t) and the substrate is locally approximated by a (large) circle of radius R(P), where P(t) represents the intersection point of the substrate and the straight line that connects the centers of the two circles (at all times R≫ r).
This gives rise to a parameterization of the interface profile X⃗(θ, t)=(x(θ, t), y(θ, t))^T as

    {[    x(θ,t)=P(t)+r(t)sinθ,; y(θ,t)=r(t)(cosθ-cosθ̂), ].  θ∈[-θ̂,θ̂],

where θ̂=θ_i+α(P) with  rsinθ̂ = Rsinα.
The conserved particle area is

    A_0= r^2ζ(θ̂) - R^2ζ(α),

where ζ(α) = α -cosαsinα.
For 0 ≤α≪ 1, this yields the following  identities


    √(A_0)/R   := √(sin^2α/sin^2θ̂ ζ(θ̂) - ζ(α)) ,
    r/√(A_0)   :=√(sin^2α/sin^2α ζ(θ̂) - sin^2θ̂ ζ(α)).



Taylor expanding (<ref>) about α = 0 gives


    √(A_0)/R   = √(ζ(θ_i))/sinθ_iα + O(α^2), 
    r/√(A_0)   = 1/√(ζ(θ_i)) + O(α).



Recalling (<ref>), the total free energy of the approximate 2D  system can be written as

    W = 2γ_0(rθ̂ - Rcosθ_i α).

Using (<ref>) in (<ref>)  yields

    W(P)=2γ_0√(A_0) (√(ζ(θ_i)) + sin^2θ_i/3 √(ζ(θ_i))α + O(α^2)).

Taking the time derivative of the total free energy and using (<ref>), we obtain

    Ẇ(P, Ṗ)    =  2γ_0 √(A_0) sin^2θ_i/3 √(ζ(θ_i))α^'(P)Ṗ +O(α^2)
       ≈2γ_0 A_0 sin^3θ_i /3 ζ(θ_i)κ^'(P) Ṗ,

where we introduced the substrate curvature (in 2D) κ(P) = R(P)^-1  and   primes  denote  derivatives with respect to P.

To this point, we  restricted ourselves to the case when the curvature of the substrate is positive at the point of interest, implying that R(P(t)) > 0 and α≥ 0, see Fig. <ref>(a).
When the substrate curvature is negative, R(P(t)) < 0;  equations (<ref>), (<ref>), (<ref>) and (<ref>) hold as well except that  0<-α≪ 1, as shown in Fig. <ref>(b).

Recalling (<ref>), we write the dissipation function  (<ref>) as

    Φ=1/2k_b T/D_s ν∫_(t)|j⃗|^2  S =1/2k_b T/D_s ν∫_(t)|J|^2  S,

where J(θ)=j⃗·τ is the magnitude of the flux and τ is the unit interface tangent.
Using the kinematic equation (<ref>), we find

    ∂_s J(θ) = -1/Ω_0V_n(θ)     J(±θ̂)=0,

where θ∈[-θ̂,θ̂] and V_n(θ) is

    V_n(θ) = Ṗsinθ + ṙ(1-cosθ_icosθ),

on recalling (<ref>) and n⃗ = (sinθ, cosθ)^T.
Using (<ref>) and integrating (<ref>) then yields

    J(θ)     = -1/Ω_0∫_-θ̂^θ r V_n(θ) θ
       = -1/Ω_0Ṗ r(cosθ_i - cosθ)+ O(α),

where we note ṙ = O(α) because of (<ref>).
Inserting (<ref>) into (<ref>) and using (<ref>),  we find

    Φ(Ṗ)    = 1/2k_b T/D_s ν∫_-θ̂^θ̂|J^2| rθ
       = C_2^0(θ_i) B^-1 r^3Ṗ^2 +O(α^2)
       ≈ C_2^0(θ_i) B^-1√(A_0^3/ζ^3(θ_i))Ṗ^2,

where B=D_sνΩ_0^2/k_B T and 
    C_2^0(θ_i)=1/2(θ_i + 2θ_icos^2θ_i - 3sinθ_icosθ_i).

Note that for the dynamics in 2D, the flux J(θ) is completely determined from
Ṗ and V_n(θ), and hence the coefficient of Φ(Ṗ)∝Ṗ^2 is
readily obtained from (<ref>)
(this is not the case in 3D).



Using (<ref>) and (<ref>) in (<ref>), and applying the Onsager principle by minimizing the Rayleighian ℛ with respect to Ṗ, we thus obtain the following ODE for the state variable P (see also Equation (16) in <cit.>):

    d P/ d t = - Bγ_0 C_2(θ_i)/√(A_0) κ ' (P),
    
        P(0) = P_0,

where

    C_2(θ_i) =2sin^3(θ_i) √(θ_i- cosθ_isinθ_i)/3(θ_i+2θ_icos^2θ_i - 3sinθ_icosθ_i).


Equation (<ref>) prescribes the velocity of the particle on the substrate  when |α|≪ 1.
It shows that the velocity is proportional to the substrate curvature gradient and is inversely proportional to the length scale of the particle, i.e., √(A_0).
Moreover, the coefficient C_2(θ_i) is a monotone decreasing function for θ_i∈[0,π].
This implies that the velocity tends to  zero as θ_i  approaches π. Physically, a larger contact angle θ_i implies a weaker coupling between the particle and the substrate (with the contact length or area approaching zero for θ_i→π), and hence the migration velocity Ṗ becomes less responsive to substrate curvature gradients.




 §.§ Numerical  validation


We choose the length scale L_0, time scale L_0^4/Bγ_0 and use the quantities with hats (·̂) to denote dimensionless physical quantities.
The full model (<ref>) can then be rewritten in  dimensionless form  as

    V̂_n = Δ̂_sℋ̂,

with boundary conditions at the contact line:



  (i)attachment condition

    V̂⃗̂·n⃗_w = 0;


  (ii)contact angle condition

    n⃗·n⃗_w + cosθ_i=0;


  (iii)zero-flux condition

    n⃗_c·∇̂_sℋ̂ = 0.



The reduced model, in  dimensionless form, is

    dP̂/ dt̂ = -  C_2(θ_i)/√(Â_0) κ̂ ' (P̂),
    P̂(0) = P̂_0.

We then conduct numerical comparisons between the two models.
We employ a parametric finite element method to solve the full sharp-interface model (see <cit.>).
The system of ODEs for the reduced model are solved via the forward Euler method.


We first perform a full model simulation of the migration of a small particle (Â_0=0.4) on a substrate, which is modeled by a curve satisfying κ̂^'(s)=-0.01, where  s is the arc length parameter.
Initially, we place a square particle at a position with κ̂=0.05 and then capture the dynamics of the particle  to a position with κ̂ = -0.05.
As shown in Fig. <ref>, we see observe that the particle quickly adopts a circular shape and maintains it while gradually migrates towards a position of lower curvature.

We next conduct a series of simulations for a particle on substrates with different curvature gradients, where the other parameters are fixed (Â_0=1, θ_i=π/3).
The numerical results obtained from the full model are presented in Fig. <ref>(a).
We observe that the position of the particle varies approximately linearly with time and the speed of the particle is proportional to the substrate curvature gradient.
In particular, the speed of the particle is very similar to the analytical results from the reduced model, as shown in Fig. <ref>(b).

To further validate the reduced model, we next study the dependence of the particle velocity on A_0 and θ_i, and the curvature gradient of the substrate is fixed as κ̂^' = -0.01. It can be readily seen that the velocity of the particle is inversely proportional to √(A_0) (Fig. <ref>(c)) and proportional to C_2(θ_i) (Fig. <ref>(d)). In particular, the results show excellent quantitative agreement with the reduced model.








Finally, we investigate the dynamics of a particle on a general sinusoidal substrate described by y = 4sin(x/4) with  Â_0=1 and θ_i=π/3.
The numerical results from the full sharp-interface model are compared with those from solving the system of ODEs in the reduced model, (<ref>).
The agreement between the full and reduced model results are in excellent quantitative agreement as observed in Fig. <ref>.
This not only validates our reduced-order model in <ref>, but also
verifies the accuracy of the numerical results from the full sharp-interface model. Moreover, we note that these theoretical results are also consistent with those reported in <cit.>.





§ DYNAMICS IN THREE DIMENSIONS




 §.§ The reduced model



Inspired by the work in <ref>, we now consider the extension of the reduced model from 2D to 3D.
Given the position P⃗ =(P_x, P_y, 0)^T on the substrate surface, we  locally approximate the substrate at P⃗  as a sphere of radius R(P⃗)=1/κ(P⃗), where κ(P⃗)=1/2(κ_1 + κ_2) is the mean curvature of the substrate with κ_1 and κ_2 being the two principal curvatures. We parameterize the interface (t) by X⃗:=X⃗(θ,ϕ,t)

    {[ x(θ,φ,t)=P_x(t)+r(t)sinθcosϕ,; y(θ,φ,t)=P_y(t)+r(t)sinθsinϕ,;   z(θ,φ, t)=r(t)(cosθ-cosθ̂), ].

for θ∈[0,θ̂] and ϕ∈[0,2π], where θ̂= θ_i + α with α satisfying Rsinα = rsinθ̂.
The volume of the particle is the volume difference of two spherical caps

    V_0 =π/3(r^3 η(θ̂)-R^3 η(α)),

with η(α) = (1-cosα)^2(2+cosα).
The total free energy  of the 3D system is

    W = 2πγ_0[r^2 (1-cosθ̂) - R^2 cosθ_i (1-cosα)].

Using (<ref>) in (<ref>) leads to

    W(P⃗) =    2γ_0√(9 V_0^2 π/η^2(θ_i)) (1-cosθ_i-1/2cosθ_isin^2θ_i)
       +3γ_0 V_0/2(1+cosθ_i)^2/(2+cosθ_i)1/R(P⃗) + O(α^2).




Using the fact that the substrate curvature κ(P⃗) = 1/R(P⃗), the time derivative of the energy becomes

    Ẇ(P⃗, Ṗ⃗̇) ≈3γ_0 V_0/2(1+cosθ_i)^2/(2+cosθ_i)∇_Γκ·Ṗ⃗̇,

where ∇_Γ represents the  curvature gradient along the substrate surface.


In contrast to the 2D case,  the flux contribution on the spherical cap in  3D cannot be uniquely determined by Eq. (<ref>) although the normal velocity of the interface is given.
Nevertheless, the dissipation function can be connected with and interpreted as the Wasserstein distance in the framework of minimum dissipation <cit.>.
This leads to a constrained minimization problem:


   min_j⃗Φ(Ṗ⃗̇)=1/2k_b T/D_s ν∫_(t)|j⃗|^2 S       

     ∇_s·j⃗ = -Ṗ⃗̇·n⃗/Ω_0          (t),

     j⃗·n⃗_c=0          Γ(t).


in which the dissipation function Φ(Ṗ⃗̇) is obtained through minimization
with respect to j⃗ subject to the constraint imposed by the continuity equation.
This system satisfies  rotational invariance, meaning that

    Φ(Ṗ⃗̇)∝ |Ṗ⃗̇|^2.

We then introduce the dimensionless flux

    j̃⃗̃=j⃗/r Ω^-1 |Ṗ⃗̇|,

to obtain

    Φ(Ṗ⃗̇)   = 1/2k_b T/D_s ν∫_(t)r^2Ω^-2|Ṗ⃗̇|^2 |j̃⃗̃|^2  S
       =1/2k_b T/D_sνΩ_0^2 r^4 |Ṗ⃗̇|^2 ∫_(t)|j̃⃗̃|^2 
       =k_b T/D_sνΩ_0^2 (3 V_0/π η(θ_i))^4/3 |Ṗ⃗̇|^2 m(θ_i),

where  (t) is  S(t) scaled by the dimension r and m(θ_i)=1/2∫_(t)|j̃⃗̃|^2 is a dimensionless function of the Young angle θ_i.
In practice, we compute m(θ_i) via the  minimization problem (<ref>) as described in <ref>.

Combining Eqs. (<ref>) and (<ref>) with Eq. (<ref>) and applying the Onsager principle (minimizing the Rayleighian ℛ with respect to Ṗ⃗̇), we  obtain our reduced-order model for the motion of a particle in 3D.
The particle velocity is

    dP⃗/ d t = - Bγ_0 C_3(θ_i)/√(V_0) ∇_Γκ(P⃗),
    P⃗(0) = P⃗_0,

where

    C_3(θ_i) = π/4(1-cos^2θ_i)^2/m(θ_i)√(π η(θ_i)/3).

Note that the application of the Onsager variational principle to the dynamics in 3D
consists of two steps:
(i) Obtaining the dissipation function Φ(Ṗ⃗̇)  by minimizing the rate of dissipation
with respect to j⃗ subject to the constraint imposed by the continuity equation
for a prescribed Ṗ⃗̇;
(ii) Determination of the migration velocity Ṗ⃗̇ by minimizing the Rayleighian ℛ with respect to Ṗ⃗̇.


The particle moves in a direction opposite the curvature gradient, meaning that the trajectory of the particle is purely determined by the substrate topography.
While the other physical parameters (e.g.,B, γ_0, V_0 and θ_i) only affect how fast the particle moves in this particular trajectory.



 §.§ Numerical validation

To numerically confirm the reduced model for particle motion in 3D, Eq. (<ref>), we  compare these and the full model, Eq. (<ref>).
Similar to the 2D case, we normalize all lengths by L_0 and times by L_0^4/Bγ_0 such that the reduced model can be written in the dimensionless form as

    dP̂⃗̂/ dt̂ = - C_3(θ_i)/√(V̂_0) ∇̂_Γκ̂(P̂⃗̂),
    P̂⃗̂(0) = P̂⃗̂_0.

We solve for the particle trajectory via the forward Euler method.
For the full model, as described in <ref>, we employ parametric finite element approximations - see  <cit.>.
As shown in Fig. <ref>(a), we consider an egg-carton shape  substrate surface

    z(x,y) = 1/5[sin^2(x/2)+sin^2(y/2)].

We choose θ_i=π/2 and place a small volume particle  (V̂_0=20^-3)  at the point (x_0, y_0, z(x_0,y_0)).
In our computation, we start with two different initial positions: (x_0,y_0)=(π+1/10, π+1/2) and (x_0,y_0)=(π+3/10, π+1/10).
The trajectory of the particle in the reduced and full models are compared in Fig. <ref>(b) and Fig. <ref>(c).
Note that the two models show  excellent agreement over the entire trajectories,  confirming our reduced model for the dynamics of  the particle in 3D.

To further assess the reduced model, we next consider the more general cases by varying the volume of the particle and Young angle θ_i.
We compare the time history of the x-position of the particle (the initial particle position was (x_0,y_0)=(π+1/10, π+1/2)) from the two models by fitting m(θ_i) in (<ref>).
The results are shown in Figs. <ref>(d)-(f).
Again,  we observe excellent consistency between the full and reduced models. 








In our final numerical test, we focus on the verification of  the dissipation function (<ref>) in the reduced model.
m(θ_i) in (<ref>) can be computed in a minimization framework, as discussed in <ref>.
For the full model,  we introduce an analogue of m(θ_i) as

    m̂(θ_i)= 1/2 |P̂⃗̂/t̂|^2(π η(θ_i)/3 V̂_0)^4/3∫_(t)|∇̂_sĤ|^2 Ŝ.

If the dissipation law for the migrating particle is adequately approximated by (<ref>), then m̂(θ_i) should be similar to the constant m(θ_i) in time.
We test this for several different Young angles (V_0=L_0^3/80^3).
The time history of m̂(θ) as well as the constant m(θ_i) are shown in Fig. <ref>.
Indeed, we observe a small oscillation of m̂(θ_i) about a constant which is  slightly larger than m(θ_i).
These discrepancies may be associated with the errors in the asymptotic approximations in the reduced model or  numerical  errors in the computation of the full sharp-interface model, see <ref>.


Finally, we assess the quantitative agreement between m(θ_i) and m̂(θ_i) as functions of θ_i in Fig. <ref>.
Again, the agreement is good  and consistent across two orders of magnitude from 10^-2 to 1. Furthermore, m(θ_i), obtained from the constrained minimization, is everywhere very close to m̂(θ_i) from the full model computation, as expected since m(θ_i) is a theoretical value in the minimum dissipation framework, while m̂(θ_i) is  a numerical value measured in simulation.  The relative error  may be attributed to the following:
(i) In the full model computation, the particle is not small enough compared to
the radius of  curvature  of the substrate, implying that the equilibration of the particle shape is not fast enough compared with the migration along the substrate as assumed in the reduced model;
(ii) The numerical results produced in the full model computation are not sufficiently accurate;
we  observe that when the parameter θ_i is near π (e.g., 3π/4),
a large oscillation in m̂(θ_i) occurs which may make the full model numerical computation (see <ref>) unreliable.

This further suggests the utility of the reduced-order model  (<ref>) for such particle migration studies.








§ GENERALIZATIONS




 §.§ For chemically inhomogeneous substrates


We employed our framework of the Onsager principle to the case of deposited thin films on a flat and chemically inhomogeneous substrate. This will drive the particle to transport over the substrate in a similar behaviour of the droplet, see <cit.>. We assume that the material parameter θ_i of the substrate is a function of the position, and is varying pretty slowly. Thus there is a fast time scale for the thin film to form a circular shape and a slow time scale for the migration of the particle.

In the 2D case, using (<ref>) and (<ref>) and letting R→+∞, it is not difficult to obtain that

    W(P) = 2γ_0 √(A_0 ζ(θ_i)),

of which taking time derivative gives us

    Ẇ(P,Ṗ) = γ_0 √(A_0) ζ^'(θ_i) θ_i^'(P) Ṗ/√(ζ(θ_i)).

We note that the dissipation function for the interface will stay unchanged and can be computed by (<ref>) as well. Now combining (<ref>) and (<ref>) and  applying the Onsager principle yields the following ODE for P,

    d P/ d t = - 1/2B γ_0/A_0 ζ(θ_i) ζ^'(θ_i)/C_2 d^0(θ_i) θ_i^'(P),
    
        P(0) = P_0.

This implies that the velocity of the particle is proportional to the gradient of the Young's angle θ_i. Moreover, the parameter θ_i will also determine how fast the particle move over the substrate.

For the 3D particle on the flat and chemically inhomogeneous substrates, on recalling (<ref>), we have

    W(P⃗) = 2γ_0√(9 V_0^2 π) C_0(θ_i),

where

    C_0(θ_i) = [η(θ_i)]^-2/3(1-cosθ_i-1/2cosθ_isin^2θ_i).

Taking time derivative of W(P⃗) then leads to

    Ẇ(P⃗, Ṗ⃗̇) = 2γ_0 √(9 V_0^2 π) C_0^'(θ_i)∇_Γθ_i(P⃗)·Ṗ⃗̇,

where ∇_Γ is again the gradient along the substrate surface. On recalling the dissipation function (<ref>) in 3D and using the Onsager principle, we obtain the following ODE system for the particle velocity

    dP⃗/ d t = - Bγ_0/√(V_0^2) (π^5η^4(θ_i)/9 m^3(θ_i))^1/3C_0^'(θ_i) ∇_Γθ_i(P⃗),
    P⃗(0) = P⃗_0,

which implies that the particle moves in a direction that is aligned with the gradient of the parameter θ_i.




 §.§ Discussions on anisotropic surface energies

We next consider the case when the surface energy of the thin film is anisotropic and modeled by a convex anisotropy function γ(n⃗). Again, we assume that the size of the particle is much smaller than the curvature radius of the substrate surface so that there exists a short time scale for the particle to form an anisotropic quasi-static shape and a long time scale for the particle to migrate along the substrate.

On the short time scale, we assume that the particle forms a shape that is no longer circular or spherical but depends on γ(n⃗) and also the local substrate topology, i.e., the unit normal n⃗_w to the substrate surface.
Starting from a flat substrate surface with a constant n⃗_w, we note that the zeroth-order particle shape 
𝒮_γ can be constructed via the Winterbottom construction <cit.> or the well-known Cahn-Hoffman vector formulation <cit.>.
As the substrate surface becomes gently curved with a nonzero ∇_Γn⃗_w, 
the particle shape and the free energy W(P⃗) will be slightly modified by ∇_Γn⃗_w. 
To the first order, we consider the linear dependence of W(P⃗) on ∇_Γn⃗_w by noting that 
the orientation of 𝒮_γ relative to the principal axes of ∇_Γn⃗_w is involved. 
This is totally different from the isotropic case, where the particle shape possesses a rotational symmetry so that the linear dependence of W(P⃗) on ∇_Γn⃗_w is
through Tr[∇_Γn⃗_w], i.e., the mean curvature of the substrate surface, see (<ref>).

To compute the dissipation function Φ(Ṗ⃗̇),  we can employ a technique similar to the isotropic case in the framework of constrained minimum dissipation. This leads to the minimization problem (<ref>), 
where the anisotropic effects are manifested in the integration over the anisotropic interface 𝒮_γ 
and the constraint which involves the unit normal n⃗ of 𝒮_γ.

In summary, our variational reduced-order modeling approach can also be employed to deal with solid thin films with anisotropic surface energies on topologically patterned substrates, although 
the particle migration will exhibit much more complicated trajectories.




§ CONCLUSIONS


We studied the dynamics of a small particle on  curved substrates driven by surface/interface energies and controlled by surface diffusion  with a moving contact line.
We observed two distinct time scales for the dynamical behaviour of the particle.
Small particles evolves towards their equilibrium shape (circular or spherical section in 2D or
 3D) on a fast time scale.
The small particle move  along the substrate in order to adapt itself to the substrate topography on a slow time scale.
We derived a reduced-order model to describe this particle migration using the Onsager variational principle.
The particle moves in a direction in which the curvature of the substrate decays most quickly.
The particle velocity is proportional to a material constant B=D_sνΩ_0^2/k_b T, the surface energy density γ_0, and inversely proportional to the size of the particle.
The reduced model was numerically confirmed by comparing with  numerical results obtained by numerical solution of the full dynamical model.

The main overall conclusion is that the trajectory of  small particle on a substrate is solely determined by the substrate topography, while only the rate of motion is controlled by material properties. 
This may provide some insights into the understanding of the templated dewetting on  patterned substrate.









§ ACKNOWLEDGEMENTS

This work was partially supported by the National Natural Science Foundation of China Nos. 12271414 and 11871384 (W.J.) and 
No. 12371395 (Y.W.), the NSF Division of Materials Research through Award 1507013 (D.J.S.), 
Hong Kong RGC grants CRF No. C1006-20WF and GRF No. 16306121 (T.Q.), 
and the Ministry of Education of Singapore under its AcRF Tier 2 
funding MOE-T2EP20122-0002 (A-8000962-00-00) (W.B.).





§ COMPUTATIONAL METHOD FOR THE FULL MODEL (<REF>)

The computation of the full model (<ref>) is based on the following formulation  <cit.>:


    X̂⃗̂·n⃗ = Δ̂_sĤ,
    Ĥn⃗ = -Δ̂_sX̂⃗̂.


A weak formulation is introduced for (<ref>),  where the contact angle condition (<ref>) and the zero-flux condition (<ref>) are implemented via the variational formulation (see <cit.> for details).
To  ensure satisfaction of the attachment condition (<ref>), the velocity of the contact line is forced to be tangential to the substrate surface.
We  employ a piecewise linear element method in space and a backward Euler discretization in time to obtain a parametric approximation, where the interface surface is approximated as a polyhedron.
This  discretization implies that  contact points may not exactly lie on the curved substrate at the following time step.
Therefore, an orthogonal projection of these points onto the substrate surface is required; this may incur some additional numerical errors beyond that associated with the discretization of the geometric equation.



§ THE DISSIPATION FUNCTION

We  introduce a dimensionless flux on a spherical section of the surface

    j̃⃗̃ = (θ,ϕ) e⃗_θ +(θ,ϕ)e⃗_ϕ,

where e⃗_θ and e⃗_ϕ represent the unit vector in the axial and azimuthal directions, respectively.

To calculate the dissipation function in 3D, we consider the constrained minimization problem (<ref>) which is rewritten in the spherical coordinate and in terms of j̃⃗̃ as


    min_j̃⃗̃1/2k_b T r^4/D_s ν Ω_0^2∫_0^2π∫_0^θ_i(^2 + ^2)sinθθϕ,
        ∂( sinθ)/∂θ +∂/∂ϕ = -(Ṗ⃗̇/|Ṗ⃗̇|·n⃗^ϕ)sin^2θ,
         (θ_i,ϕ) = 0,


for θ∈[0,θ_i] and ϕ∈[0,2π], where n⃗^ϕ = (cosϕ, sinϕ, 0)^T.

We then discretize the problem in (θ,ϕ)∈[0,θ_i]×[0,2π].
We introduce the following notations:


    h_θ = θ_i/N,   θ_j = j h_θ     0≤ j≤ N, 
         h_ϕ = 2π/M,  ϕ_k = k h_ϕ     0≤ k≤ M,


and _jk≈(θ_j,ϕ_k), _jk≈(θ_j,ϕ_k) and  rewrite (<ref>) in the discrete form.
In particular, the objective function, up to a multiplicative constant, is approximated as

    m(θ_i)   =1/2∫_0^2π∫_0^θ_i(^2+^2) sinθθϕ
       ≈1/2∑_j=1^N∑_k=1^M(_jk + _jk)2πθ_i/NM,

where

    _jk=1/4∑_l_1=0,1∑_l_2=0,1_(j-l_1)(k-l_2)^2sinθ_j-l_1,

and _jk is defined similarly.

We  introduce the finite difference discretization operators as follows:


    δ_θ_jk   =_j+1,k-_j-1,k/2h_θ,  δ_ϕ_jk=_j,k+1-_j,k-1/2h_ϕ,
    δ_θ^+_jk   =_j+1,k-_j,k/h_θ,  δ_θ^-_jk=_j,k-_j-1,k/h_θ,


and observe that  _0,k and _0,k make no contributions to the objective function (because sin0=0).
The constraints are discretized naturally as


  * for j=1, 1≤ k≤ M,


    δ_θ^+(sinθ)_1k+δ_ϕ_1k=-[Ṗ⃗̇/|Ṗ⃗̇|·n⃗^ϕ_1k]sin^2θ_1k,


  * for 1< j < N and 1≤ k≤ M,

    δ_θ(sinθ)_jk+δ_ϕ_jk=-[Ṗ⃗̇/|Ṗ⃗̇|·n⃗^ϕ_jk]sin^2θ_jk,


  * for j=N, 1≤ k≤ M,

    _N,k=0.




This is the discrete minimization problem for the objective function (<ref>) with the constraints (<ref>); it is a quadratic minimisation problem with linear constraints for {_jk} and {_jk} for j=1,…, N and k=1,…, M.
Thus it can be solved directly as a linear system corresponding to the Karush-Kuhn-Tucker (KKT) condition.








§ REFERENCES


10
url<#>1urlprefixURL href#1#2#2 #1#1

Thompson12
C. V. Thompson, Solid-state dewetting of thin films, Annu. Rev. Mater. Res. 42
  (2012) 399–434.

Leroy16
F. Leroy, F. Cheynis, Y. Almadori, S. Curiotto, M. Trautmann, J. Barbé,
  P. Müller, How to control solid state dewetting: A short review, Surf.
  Sci. Rep. 71 (2) (2016) 391–409.

Naffouti17
M. Naffouti, R. Backofen, M. Salvalaglio, T. Bottein, M. Lodari, A. Voigt,
  T. David, A. Benkouider, I. Fraj, L. Favre, A. Ronda, I. Berbezier,
  D. Grosso, M. Abbarchi, M. Bollani, Complex dewetting scenarios of ultrathin
  silicon films for large-scale nanoarchitectures, Sci. Adv. 3 (11) (2017)
  eaao1472.

Srolovitz86a
D. J. Srolovitz, S. A. Safran, Capillary instabilities in thin films: I.
  energetics, J. Appl. Phys. 60 (1) (1986) 247–254.

Mullins57
W. W. Mullins, Theory of thermal grooving, J. Appl. Phys. 28 (3) (1957)
  333–339.

Young1805
T. Young, An essay on the cohesion of fluids, Philos. Trans. R. Soc. London 95
  (1805) 65–87.

Ye11b
J. Ye, C. V. Thompson, Templated solid-state dewetting to controllably produce
  complex patterns, Adv. Mater. 23 (13) (2011) 1567–1571.

Amram12
D. Amram, L. Klinger, E. Rabkin, Anisotropic hole growth during solid-state
  dewetting of single-crystal Au–Fe thin films, Acta Mater. 60 (6-7) (2012)
  3047–3056.

Jiang12
W. Jiang, W. Bao, C. V. Thompson, D. J. Srolovitz, Phase field approach for
  simulating solid-state dewetting problems, Acta Mater. 60 (15) (2012)
  5578–5592.

Jiang16
W. Jiang, Y. Wang, Q. Zhao, D. J. Srolovitz, W. Bao, Solid-state dewetting and
  island morphologies in strongly anisotropic materials, Scripta Mater. 115
  (2016) 123–127.

Zucker16
R. V. Zucker, G. H. Kim, J. Ye, W. C. Carter, C. V. Thompson, The mechanism of
  corner instabilities in single-crystal thin films during dewetting, J. Appl.
  Phys. 119 (12) (2016) 125306.

Jiang19xi
W. Jiang, Q. Zhao, Sharp-interface approach for simulating solid-state
  dewetting in two dimensions: a Cahn-Hoffman ξ-vector formulation,
  Physica D 390 (2019) 69–83.

Jiang2018curved
W. Jiang, Y. Wang, D. J. Srolovitz, W. Bao, Solid-state dewetting on curved
  substrates, Phys. Rev. Mater. 2 (11) (2018) 113401.

JiangZB20
W. Jiang, Q. Zhao, W. Bao, Sharp-interface model for simulating solid-state
  dewetting in three dimensions, SIAM J. Appl. Math. 80 (4) (2020) 1654–1677.

Boc22stress
F. Boccardo, F. Rovaris, A. Tripathi, F. Montalenti, O. Pierre-Louis,
  Stress-induced acceleration and ordering in solid-state dewetting, Phys. Rev.
  Lett. 128 (2) (2022) 026101.

Garcke23diffuse
H. Garcke, P. Knopf, R. Nürnberg, Q. Zhao, A diffuse-interface approach for
  solid-state dewetting with anisotropic surface energies, J. Nonlinear Sci.
  33 (2) (2023) 34.

Giermann05
A. L. Giermann, C. V. Thompson, Solid-state dewetting for ordered arrays of
  crystallographically oriented metal particles, Appl. Phys. Lett. 86 (12)
  (2005) 121903.

Cheng06templated
J. Y. Cheng, C. A. Ross, H. I. Smith, E. L. Thomas, Templated self-assembly of
  block copolymers: top-down helps bottom-up, Adv. Mater. 18 (19) (2006)
  2505–2521.

Wang11
D. Wang, R. Ji, P. Schaaf, Formation of precise 2D Au particle arrays via
  thermally induced dewetting on pre-patterned substrates, Beilstein J.
  Nanotechnol. 2 (1) (2011) 318–326.

Wang2013solid
D. Wang, P. Schaaf, Solid-state dewetting for fabrication of metallic
  nanoparticles and influences of nanostructured substrates and dealloying,
  Physica Status Solidi (a) 210 (8) (2013) 1544–1551.

Lu16nanostructure
L.-X. Lu, Y.-M. Wang, B. M. Srinivasan, M. Asbahi, J. K. Yang, Y.-W. Zhang,
  Nanostructure formation by controlled dewetting on patterned substrates: A
  combined theoretical, modeling and experimental study, Sci. Rep. 6 (1) (2016)
  32398.

Ruffino17experimental
F. Ruffino, Experimental analysis on the molten-phase dewetting characteristics
  of AuPd alloy films on topographically-structured substrates, Metals 7 (9)
  (2017) 327.

Ahn80
T.-M. Ahn, J. K. Tien, P. Wynblatt, Coarsening kinetics of platinum particles
  on curved oxide substrates, J. Catal. 66 (2) (1980) 335–346.

Klinger12
L. Klinger, E. Rabkin, Capillary-driven motion of nanoparticles attached to
  curved rigid substrates, Acta Mater. 60 (17) (2012) 6065–6075.

lv14substrate
C. Lv, C. Chen, Y.-C. Chuang, F.-G. Tseng, Y. Yin, F. Grey, Q. Zheng, Substrate
  curvature gradient drives rapid droplet motion, Phys. Rev. Lett. 113 (2)
  (2014) 026101.

Galatola18s
P. Galatola, Spontaneous capillary propulsion of liquid droplets on substrates
  with nonuniform curvature, Phys. Rev. Fluids 3 (10) (2018) 103601.

Mccarthy19
J. McCarthy, D. Vella, A. A. Castrejón-Pita, Dynamics of droplets on cones:
  self-propulsion due to curvature gradients, Soft Matter 15 (48) (2019)
  9997–10004.

Chen21self
Y. Chen, X. Xu, Self-propulsion dynamics of small droplets on general surfaces
  with curvature gradient, Phys. Fluids 33 (8) (2021) 082107.

Sykes22droplet
T. C. Sykes, B. D. Fudge, M. A. Quetzeri-Santiago, J. R. Castrejón-Pita,
  A. A. Castrejón-Pita, Droplet splashing on curved substrates, J. Colloid
  Interface Sci. 615 (2022) 227–235.

Onsager31a
L. Onsager, Reciprocal relations in irreversible processes. I., Phys. Rev.
  37 (4) (1931) 405.

Onsager31b
L. Onsager, Reciprocal relations in irreversible processes. II., Phys. Rev.
  38 (12) (1931) 2265.

Reina15entropy
C. Reina, J. Zimmer, Entropy production and the geometry of dissipative
  evolution equations, Phys. Rev. E 92 (5) (2015) 052117.

Van23thermodynamic
T. Van Vu, K. Saito, Thermodynamic unification of optimal transport:
  Thermodynamic uncertainty relation, minimum dissipation, and thermodynamic
  speed limits, Phys. Rev. X 13 (1) (2023) 011013.

Cahn94
J. W. Cahn, J. E. Taylor, Surface motion by surface diffusion, Acta Metall.
  Mater. 42 (4) (1994) 1045–1063.

Karim22
A. Mohammad Karim, A review of physics of moving contact line dynamics models
  and its applications in interfacial science, J. Appl. Phys. 132 (8) (2022)
  080701.

Qian06
T. Qian, X.-P. Wang, P. Sheng, A variational approach to moving contact line
  hydrodynamics, J. Fluids Mech. 564 (2006) 333–360.

Qian17
X. Xu, T. Qian, Hydrodynamic boundary conditions derived from Onsager's
  variational principle, Procedia IUTAM 20 (2017) 144–151.

Xu16
X. Xu, Y. Di, M. Doi, Variational method for liquids moving on a substrate,
  Phys. Fluids 28 (8) (2016) 087101.

Di18
Y. Di, X. Xu, J. Zhou, M. Doi, Thin film dynamics in coating problems using
  Onsager principle, Chin. Phys. B 27 (2) (2018) 024501.

Man16
X. Man, M. Doi, Ring to mountain transition in deposition pattern of drying
  droplets, Phys. Rev. Lett. 116 (6) (2016) 066101.

Zhang22effective
Z. Zhang, X. Xu, Effective boundary conditions for dynamic contact angle
  hysteresis on chemically inhomogeneous surfaces, J. Fluid Mech. 935 (2022)
  A34.

Doi11
M. Doi, Onsager's variational principle in soft matter, J. Phys. Condens.
  Matter 23 (28) (2011) 284118.

Doi13book
M. Doi, Soft matter physics, Oxford University Press, 2013.

Doi15
M. Doi, Onsager principle as a tool for approximation, Chin. Phys. B
  24 (020505) (2015) 1674–1056.

JZOnsager
W. Jiang, Q. Zhao, T. Qian, D. J. Srolovitz, W. Bao, Application of onsager's
  variational principle to the dynamics of a solid toroidal island on a
  substrate, Acta Mater. 163 (2019) 154–160.

Suo97
Z. Suo, Motions of microscopic surfaces, Adv. Appl. Mech. 33 (1997) 193–294.

Barrett20
J. W. Barrett, H. Garcke, R. Nürnberg, Parametric finite element
  approximations of curvature driven interface evolutions, Handb. Numer. Anal.
  (Andrea Bonito and Ricardo H. Nochetto, eds.) 21 (2020) 275–423.

BGNZ23
W. Bao, H. Garcke, R. Nürnberg, Q. Zhao, A structure-preserving finite
  element approximation of surface diffusion for curve networks and surface
  clusters, Numer. Methods Partial Diff. Equ. 39 (1) (2023) 759–794.

Malinowski20advances
R. Malinowski, I. P. Parkin, G. Volpe, Advances towards programmable droplet
  transport on solid surfaces and its applications, Chem. Soc. Rev. 49 (22)
  (2020) 7879–7892.

Winterbottom67
W. Winterbottom, Equilibrium shape of a small particle in contact with a
  foreign substrate, Acta Metall. 15 (2) (1967) 303–310.

Hoffman72
D. W. Hoffman, J. W. Cahn, A vector thermodynamics for anisotropic surfaces: I.
  fundamentals and application to plane surface junctions, Surface Science 31
  (1972) 368–388.




