On the Stability of Surface Growth: The Effect of a Compliant Surrounding Medium

In a recent paper Abeyaratne et al. (J. Mech. Phys. Solids 167:104958, 2022) concerning the stability of surface growth of a pre-stressed elastic half-space with surface tension, it was shown that steady growth is never stable, at least not for all wave numbers of the perturbations, when the growing surface is traction-free. On the other hand, steady growth was found to be always stable when growth occurred on a flat frictionless rigid support and the stretch parallel to the growing surface was compressive. The present study is motivated by these somewhat unexpected and contrasting results. In this paper the stability of a pre-compressed neo-Hookean elastic half-space undergoing surface growth under plane strain conditions is studied. The medium outside the growing body resists growth by applying a pressure on the growing surface. At each increment of growth, the incremental change in pressure is assumed to be proportional to the incremental change in normal displacement of the growing surface. It is shown that surface tension stabilizes a homogeneous growth process against small wavelength perturbations while the compliance of the surrounding medium stabilizes it against large wavelength perturbations. Specifically, there is a critical value of stretch, λcr∈(0,1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda _{\mathrm{cr}} \in (0,1)$\end{document}, such that growth is linearly stable against infinitesimal perturbations of arbitrary wavelength provided the stretch parallel to the growing surface exceeds λcr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda _{\mathrm{cr}}$\end{document}. This stability threshold, λcr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda _{\mathrm{cr}}$\end{document}, is a function of the non-dimensional parameter σκ/G2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sigma \kappa /G^{2}$\end{document}, which is the ratio between two length-scales σ/G\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sigma /G$\end{document} and G/κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$G/\kappa $\end{document}, where G\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$G$\end{document} is the shear modulus of the elastic body, σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sigma $\end{document} is the surface tension, and κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\kappa $\end{document} is the stiffness of the surrounding compliant medium. It is shown that (a)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(a)$\end{document}λcr→1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda _{\mathrm{cr}} \to 1$\end{document} as κ→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\kappa \to 0$\end{document} and (b)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(b)$\end{document}λcr→0+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda _{\mathrm{cr}} \to 0^{+}$\end{document} as κ→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\kappa \to \infty $\end{document}, thus recovering the results in Abeyaratne et al. (J. Mech. Phys. Solids 167:104958, 2022) pertaining to the respective limiting cases where growth occurs (a)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(a)$\end{document} on a traction-free surface and (b)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(b)$\end{document} on a frictionless rigid support. The results are also generalized to include extensional stretches.


Surface Growth Is Ubiquitous
Surface growth -the addition of new material to the surface of a solid body -occurs in a variety of contexts, perhaps the most familiar being the solidification of water on the surface of an ice cube below the freezing temperature.Other examples include 3D-printing, chemical vapor deposition, the growth of hard tissue like bone, and the polymerization of actin networks in the cytoskeleton of a biological cell.

Compliant Resistance to Surface Growth
Actin polymerization plays a central role in cell biology; see, for example, the review articles [5], [26], [32] and the references therein.In a series of intriguing experiments, Parekh et al. [30], Chaudhuri et al. [12] and Bieling et al. [7], polymerized an actin network between an AFM (atomic force microscope) cantilever and a fixed surface.One end of the cantilever was functionalized with an actin nucleating agent and this led to polymerization at that end.Once the growing actin network reached the fixed surface, it pushed against the cantilever during further growth, with the elasticity of the cantilever providing a compliant resistance to growth.Similarly, an actin network growing inside a cell pushes against the cell membrane which provides a compliant resistance to growth.

Instability of Surface Growth
Surface instability is well-known in materials science where growing surfaces can become unstable due to a coupling between growth and stress, e.g., Mullins and Serkeka [27], Asaro and Tiller [6] and Grinfeld [19].
Symmetry-breaking instabilities are frequently encountered in biology during surface growth.An example of this was seen in the striking experiments of Noireaux et al. [29], Cameron et al. [11] and others in which they chemically treated the surface of a micronsize rigid bead immersed in a solution of free actin monomers.The actin monomers were attracted to the surface of the bead where they polymerized and attached onto the actin network growing around the bead.Each newly formed layer of solid at the bead surface pushed out the previously formed solid which induced stress in the network; for modeling of this spherically symmetric growth process see Noireaux et al. [29], Dafalias et al. [14], [15] and Tomassetti et al. [35].When the growing actin shell reached a certain critical thickness it lost spherical symmetry and developed a "comet tail"; see van der Gucht et al. [20]; Prost et al. [31]; and John et al. [23].

Continuum Mechanical Modeling of Surface Growth
The addition of new material points to the surface of a growing body leads to a timedependent reference configuration, with the evolution of the reference configuration being intimately tied to the accretive growth of the body.Modeling this is challenging.Inspired by the seminal paper of Skalak et al. [33], some progress was made by Tomassetti et al. [35] when the geometry is simple, e.g., one dimensional, planar or cylindrically/spherically symmetric.An alternative Eulerian approach, which seeks to by-pass consideration of the reference configuration, has been developed by Naghibzadeh et al. [28].
The growth process itself is characterized by a kinetic law relating the propagation speed of the boundary in the reference configuration (which is a measure of the rate at which new material is added to the body) and a power conjugate thermodynamic driving force.

A Recent Study
In a recent paper [4], we examined the stability of steady surface growth of a pre-stressed elastic half-space with surface tension.While several mechanical boundary conditions at the growing surface were considered, the two of interest to us here are (a) the case where the growing surface is traction-free, and (b) the case where growth occurs on a flat frictionless rigid support.These cases were motivated by the boundary conditions at (a) the outer and (b) the inner surfaces of the actin shell growing around a bead mentioned above.By constructing steady solutions and examining their linearized stability, we showed that the steady solution in the former case, (a), is never stable, at least not for all wave numbers of the perturbations, and that it is always stable in the latter case (b) when the stretch parallel to the growth surface is compressive.The present study was motivated by these somewhat unexpected and contrasting results.
In the absence of growth, the problem studied in [4] is the classical Biot instability problem for a half-space, albeit with somewhat more general boundary conditions.In the most familiar case, where the boundary of the half-space is traction-free and surface tension is ignored, Biot [8] showed that the body became unstable through a surface instability when the stretch parallel to the boundary reached the value λ Biot (≈ 0.5437) given by the unique real, positive root of the polynomial equation λ 6 + λ 4 + 3λ 2 − 1 = 0.In such a surface instability mode, the infinitesimal deformation superposed on the homogeneous one involves periodic oscillations parallel to the boundary and exponential decay away from it.There is a substantial literature on various questions and generalizations of this problem, e.g., the stability of the deformation1 [13]; generalization beyond the neo-Hookean model to arbitrary incompressible isotropic materials [16]; the rich phenomena seen beyond plane strain [34]; generalization to hydrogels [24]; loading by electrostatic stress [22]; and so on.

The Present Study and Its Key Results
In the present study we again consider the growth of a pre-stressed neo-Hookean half-space, but now allow the medium outside the body to be compliant and to afford some resistance to growth.Specifically, the medium on the outside is assumed to apply a pressure on the boundary of the growing body such that, at each increment of growth, the incremental change in pressure is proportional to the incremental change in normal displacement of the growing surface.This was motivated in part by the examples described previously in Sect.1.2, and by the simple fact that our previous study [4] corresponds to the limiting cases of a compliant environment when its stiffness tends (a) to zero and (b) to infinity respectively.
We first construct a time-dependent spatially homogeneous growth process conforming to the boundary-initial value problem, that for brevity, we refer to as a "homogeneous solution".During such a growth process, the deformation is spatially homogeneous and timeindependent, but the body undergoes an unsteady time-dependent evolution.This is because when, for example, the driving force for growth is positive at the initial instant, the body starts to grow through accretion at its boundary.As it grows, the boundary moves outwards and so the pressure applied on it by the surrounding medium increases.As the pressure increases, the driving force decreases, and so does the rate of growth.Eventually, as t → ∞, the driving force vanishes, growth stalls, and the body reaches equilibrium.
Next, we study the linearized stability of a homogeneous solution and show that there is a critical stretch λ cr ∈ (0, 1) such that growth is stable when the compressive stretch parallel to the boundary of the body is greater than this value.There are two length-scales in the problem, σ/G and G/κ where G is the shear modulus of the neo-Hookean elastic body, σ the surface tension and κ the stiffness of the surrounding medium.We show that the stability threshold λ cr is a function of the nondimensional parameter σ κ/G 2 that is the ratio of these two length-scales.Moreover, we find that λ cr → 1 as κ → 0 and λ cr → 0 + as κ → ∞, thus recovering the aforementioned results in [4] as special cases.
In Appendix C we generalize these results to include extensional stretches.It is shown there that in general there are two critical values of stretch, λ ± cr , which demarcate the range of stability.The stretch λ − cr ∈ (0, 1) is compressive (and coincides with the stretch λ cr referred to above) while λ + cr ∈ (1, 1/λ Biot ) is extensional.The homogeneous solution is stable against all perturbations when the stretch parallel to the free surface lies between these two critical values of stretch.The stretch λ + cr → 1 when κ → 0 and λ + cr → 1/λ Biot when κ → ∞.

Organization of the Rest of the Paper
The basic equations pertaining to the problem of interest are presented in the preliminary Sect. 2. Next, in Sect.3, we specialize the field equations and boundary conditions to a semi-infinite neo-Hookean body whose boundary is a flat plane.We then construct a homogeneous solution to the problem which, as already mentioned, involves a homogeneous time-independent deformation but a non-steady time-dependent evolution of stress, driving force and boundary location.Thereafter, in Sect.4, we perturb the homogeneous solution, keeping in mind that as part of perturbing the solution, one must also perturb the domain on which it is defined, i.e., the reference configuration.This is related to the fact that the problem at hand is a free boundary-initial value problem and the (time-dependent) domain of the solution is one of the unknowns.An inhomogeneous time-dependent process is now considered and linearized about the homogeneous solution.This leads to the linearized problem.Finally in Sect. 5 we analyze the linearized problem with a focus on the question of whether perturbations grow or decay with time.In Appendix A we establish some formulae used elsewhere and in Appendix B we derive the expression for the driving force in the presence of both the boundary pressure and surface tension.The results of this paper are generalized to include extensional stretches in Appendix C.

Preliminaries
We will be concerned with plane strain throughout and will not need to refer to the out-ofplane components of vector and tensor fields.
In the reference configuration, the body occupies a time-dependent (two-dimensional) region R R (t) that is bounded by the curve S R (t).Because we are allowing for surface growth, i.e., the addition of material particles to the body at S R (t), the region R R (t) and its boundary S R (t) will be time-dependent in general.
At each instant t during an evolution process, the particle located at ξ in the reference configuration is mapped to the position y(ξ , t) in the current configuration.In the current configuration the body occupies the region R(t) = y(R R (t), t) that is bounded by the curve S(t) = y(S R (t), t).The regions occupied by the body in the reference and current configurations, and their boundaries, are shown schematically in Fig. 1.We denote the unit outward normal vectors on S R and S by n R and n respectively.
The motion of the boundary S R (t) in the reference configuration is solely due to surface growth, and when we characterize the kinetics of surface growth below, we will make use Fig. 1 Schematic depiction of the time-dependent regions R R (t) and R(t) occupied by a body in the reference and current configurations and their boundaries S R (t) and S(t).The velocities of points on the moving curves S R (t) and S(t) are V R and V, and the unit outward normal vectors are n R and n.The body is surrounded by a compliant medium of the propagation speed V R • n R of S R .On the other hand, the motion of the boundary S(t) in the current configuration is due to both growth and deformation.As the body grows and the curve S(t) moves outwards, it pushes against the surrounding medium which we assume applies a pressure on the growing body.In order to characterize this mechanical boundary condition we will need the propagation speed V • n of S(t).When speaking of velocities, it is important to distinguish between the velocities of points that move with these curves and the velocity v of the material particle that is instantaneously located on the boundary.The relation between them can be derived by differentiating y * (t) = y(ξ * (t), t) with respect to time where ξ * (t) and y * (t) are corresponding points on the propagating surfaces S R (t) and S(t).This leads to where v = ∂y(ξ , t)/∂t is particle velocity and is the deformation gradient tensor.In the absence of growth V R = o, and therefore V = v as usual.
The dissipation rate associated with surface growth can be shown to be f V R where is the thermodynamic driving force and Here μ > 0 is the difference between the chemical energy (per unit reference volume) of a free material particle unattached to the growing body and that of a material particle bound to the body; S is the Piola stress and W (F) is the strain energy density.A derivation of (3) is sketched in Appendix B, generalizing the calculation of Tomassetti et al. [35] to include both the pressure loading and surface tension.
In general irreversible processes, driving forces and their conjugate fluxes are identified from the dissipation (entropy-production) rate, and the kinetics of such processes are then taken to be characterized by relations between these fluxes and driving forces; see for example Chap.14 of Callen [10], Chap.14 of Kestin [25], and Abeyaratne and Knowles [1], [2], [3].Following this approach, we assume the kinetics of surface growth to be described by a relation between the normal growth speed V R and the conjugate driving force f : Note that this is essentially a kinetic law for the mass flux ρ R V R associated with growth where ρ R is the referential mass density.The kinetic response function V(f ) characterizes growth and the dissipation inequality f V R ≥ 0 requires that V(f )f ≥ 0. We shall assume V(f ) to be monotonically increasing with V(0) = 0: Thus, the larger the driving force, the faster the growth.When f > 0 the speed V R > 0 and so the boundary S R moves outwards and material is added to the body (accretion); when f < 0 the speed V R < 0 and so S R moves inwards and material is removed from the body (ablation).There is no growth at zero driving force.We now turn to the mechanical boundary condition.In the experiments described in [30], [12] and [7], the authors grew an actin network between an AFM cantilever and a fixed surface where the growth was resisted by the elasticity of the cantilever.Likewise, when an actin network grows inside a cell during cell locomotion, the growth is resisted by the compliant cell wall.Motivated by these examples, suppose that the growing body here is surrounded by a compliant medium that applies a pressure on the boundary S(t) as it moves outwards.The increment of pressure is assumed to be linearly related to the increment of normal displacement of the boundary, and so, in the absence of surface tension, we would take the boundary condition on S(t) to be Tn = −pn where ṗ = κV • n.
Here T is the Cauchy stress tensor and the constant parameter κ > 0 is the stiffness of the surrounding medium.In order to account for surface tension, let T + n and T − n denote the limiting values of Cauchy traction on either side of S(t), with plus denoting the outside of the growing body.Then where σ > 0 denotes the surface tension on the boundary, t is the unit tangent vector in the direction of increasing arc length s, and κ is the curvature of S(t); see for example [9].In the setting we have in mind here, T + n = −pn on the outside and T − = T is the limiting stress from within the solid body.Combining this with the two preceding equations yields the following mechanical boundary condition: The surface tension is assumed to obey the constitutive relation where λ = 1/|F −1 t| is the stretch along the boundary and W(λ) is the surface energy per unit reference area (per unit reference length in plane strain).We remark that even though σ and W do not appear explicitly in (3), the model does account for surface tension; see Appendix B. Finally, the Piola and Cauchy stress tensors satisfy the equilibrium equations, having neglected body forces.They also obey the constitutive relations for an incompressible elastic solid.Here q is the reactive pressure and incompressibility requires det Given the strain energy function W (F), the surface energy function W(λ), the kinetic response function V(f ), the stiffness κ and suitable initial conditions, the problem at hand requires us to determine the region R R (t), the motion y(ξ , t) and the stress T(y, t) conforming to the field equations ( 2), ( 9), ( 10), (11), and boundary conditions ( 4), (7).
In the rest of this paper we will limit attention to an incompressible neo-Hookean solid.In this case the strain energy function is where the constant parameter G > 0 is the infinitesimal shear modulus, and the constitutive relations (10) for the Piola and Cauchy stresses specialize to Moreover, we limit attention to the case of constant surface tension σ where W(λ) = σ λ.

Spatially Homogeneous Time-Dependent Growth Process
We now specialize the problem described in Sect. 2 to the case where the body is semiinfinite, the material neo-Hookean, the boundary S R (t) a straight line, and the deformation spatially homogeneous and time-independent.It is worth pointing out at the outset that, even though the deformation is timeindependent, the body undergoes an unsteady time-dependent process in general.This is because, as mentioned earlier, when for example the driving force for growth is positive at the initial instant, the body starts to grow through accretion on the boundary.As it grows, S(t) moves outwards and so the pressure applied on it by the surrounding medium increases.As the pressure increases, the driving force decreases, and so does the rate of growth.Eventually, as t → ∞, the driving force vanishes, growth stalls, and the body reaches equilibrium.Note that the stress is time-dependent during this process through the reactive pressure q(t) though the deformation is not.For simplicity, we will refer to this spatially homogeneous growth process as the "homogeneous solution".
In order to distinguish between the various quantities in this section and their inhomogeneous counterparts in Sect.4, mathematical precision requires that we use different symbols for them.However, in order to avoid a vast amount of notation, we will avoid doing this as much as possible, and hope that the context makes clear as to what we are referring to.There will be a handful of exceptions: in the homogeneous solution in this section, we shall append a zero to the reactive pressure q 0 , the pressure p 0 applied on the boundary by the surrounding medium, and the driving force f 0 .All components of vectors and tensors will be taken with respect to a fixed orthonormal basis {e 1 , e 2 } so that in particular, (ξ 1 , ξ 2 ) are the coordinates of a generic particle in the reference configuration.The region occupied by the body in the reference configuration at time t is the half-plane so that its boundary S R (t) is the horizontal straight line ξ 2 = Z(t) as shown in Fig. 2(a).The unit outward normal vector is e 2 and the propagation speed of a point on S R (t) is Ż e 2 .Note that Z(t) is one of the unknowns to be determined.The body undergoes a homogeneous deformation that maps2 (ξ 1 , ξ 2 ) → (x 1 , x 2 ): the stretches λ 1 and λ 2 being positive and constant.We consider λ 1 to be given.Therefore by (20) below, λ 2 is also known and constant.In the current configuration, the body occupies the half-plane where The boundary S(t) of the current configuration is the straight line x 2 = X(t) as shown in Fig. 2(b).The unit outward normal is e 2 and the propagation velocity of S(t) is Ẋ e 2 .At the initial instant, we assume The deformation gradient tensor associated with (15) has components and incompressibility, det F = 1, requires The components of the Piola and Cauchy stress tensors are found using ( 13) and ( 19) to be where, as mentioned above, q 0 (t) is the reactive pressure.
The boundary condition (7) 2 in the present context reads ṗ0 (t) = κ Ẋ(t).Integrating this with respect to t , using the initial condition X(0) = 0, and assuming the boundary to be traction-free at the initial instant, leads to p 0 (t) = κX(t) (22) as expected.The boundary condition (7) 1 with (21) 2 gives where the surface tension plays no role since the curvature of the straight boundary is zero and we are limiting attention to the case where σ is constant.Keep in mind that q 0 is the reactive pressure entering through the constitutive relation while p 0 is the pressure on the boundary of the body.From (3), ( 12), ( 19), (21) 1 and ( 23) we find the driving force to be and in view of the initial condition (18), Note from ( 24) that as the boundary moves outwards, i.e., as Z increases, the driving force decreases and vice versa.The kinetic relation ( 4) can be written using V R = Ż and (24) as the first order ordinary differential equation Figure 3 shows the f 0 , ḟ0 -phase plane for (26).The trajectory passes through the origin and has the monotonicity shown in the figure because of the properties (5) of the kinetic response function V.It is seen from the phase plane that in the case f 0 (0) > 0 (resp.f 0 (0) < 0) the driving force decreases (resp.increases) monotonically from its initial value to the value 0. Since V(0) = 0, growth stalls when f 0 (t) → 0. The phase plane shows that all motions are attracted to the equilibrium point (0, 0), and so, within the context of the present section, are stable.
Fig. 3 Schematic phase plane for the first-order ordinary differential equation (26).The monotonicity of the curve shown follows from the properties of V in (5); the curvature may be different to what is shown but that has no effect on the results In fact, one can integrate (26) to obtain3 This, together with ( 24) and ( 25), provides an algebraic equation that gives t as a function of Z.It can be readily verified that the right-hand side of ( 27) is a monotonically increasing function of Z and so this algebraic relation is uniquely invertible to give Z(t).Moreover, if the integral in (27) tends to infinity as the lower limit tends to zero, implying that it takes infinite time for f 0 (t) to approach zero.For example in the special case of linear kinetics, V(f ) = f/b where b > 0 is a constant parameter, ( 27) leads to the explicit expressions where we used ( 24) and ( 25) in getting to the second equation.Thus, as t → ∞, we see that Ż → 0 and so, eventually, the boundary stops moving and growth stalls.In summary, once Z(t) is determined from (27), the region R R (t) in the reference configuration is known from (14).Moreover, X(t) is then given by ( 17), p 0 (t) by ( 22) and q 0 (t) by (23).All other quantities can then be calculated, both their time evolution and their limiting values as t → ∞.The limiting solution is in equilibrium.

Perturbed Configuration. Linearized Problem
In order to investigate the stability of the time-dependent spatially homogeneous growth process, i.e., the homogeneous solution, found in Sect.3, we now consider a perturbation of that solution.This involves perturbing both the deformation and the reference configuration.We start by linearizing the various equations about the homogeneous solution.The linearized field equations will therefore hold on the region x 2 < X(t), −∞ < x 1 < ∞, see (16), and the boundary conditions will be enforced on the straight line x 2 = X(t).We will use the symbol .
= in an equation to imply that quadratic and smaller terms have been neglected.The solution of the linearized problem will be examined in Sect. 5.
Suppose that in the perturbed reference configuration, the body occupies the region as depicted in the left-hand figure in Fig. 4. Its boundary is the curve The perturbation ζ of the boundary is assumed to be small.If we introduce the function , then the boundary is the zero level set of g and therefore the unit outward normal vector to S R (t) is where in the last step we have introduced The propagation speed We now perturb the deformation studied in Sect. 3 so that at each instant t , (ξ 1 , ξ 2 ) → (y 1 , y 2 ) according to where The displacement components u 1 , u 2 and their derivatives are assumed to be small.The boundary S(t) of the body in the current configuration can be determined by setting x 2 = X(t) + λ 2 ζ(x 1 , t) in (32) and linearizing.This leads to the following parametric characterization of S(t): where x 1 ∈ (−∞, ∞) is the parameter and X(t) = λ 2 Z(t).The unit outward normal vector on S(t) is where here, and henceforth, we use the notation for any function h(x 1 , x 2 , t).We also need to calculate the curvature, κ, of the curve S(t) for which we use the formula 4 Substituting ( 33) into (35) and linearizing leads to

Field Equations
The deformation gradient tensor F = grad ξ y associated with the perturbed deformation (32) has components Incompressibility, det[F ] = 1, requires λ 1 λ 2 = 1 and Since the Piola stress S = GF − qF −T can be written as where we have expressed the reactive pressure as q = q 0 + q.
for any function ψ(x 1 , x 2 , t).The remaining field equation (48) 2 can now be written as On substituting (50) into the boundary conditions (49) we get The problem has now been reduced to finding ψ(x 1 , x 2 , t) and ζ(x 1 , t).Keeping in mind that the coefficient b(t) in (52) 3 is time-dependent, we seek a solution in the separable form where Substituting (53) 1 into (51) yields the ordinary differential equation whose solution that is bounded for y → −∞ is here, α and β are constants to be determined, and with no loss of generality we have assumed k > 0. Substituting (53) into the boundary conditions (52) leads to Equations ( 57) and (58) both imply that C(t)/Q(t) is time-independent.Thus, let where γ is a constant to be determined.Substituting (60) into (57)-(59) gives In writing (63) (and (59) above) we have explicitly reminded ourselves that in general, b is time-dependent: b = b(t) = 1/V (f 0 (t)); see (45) and (26).Equation (63) tells us that b(t) Ċ(t)/C(t) is time-independent, and so we can write where the constant ν is to be determined.Substituting (64) into (63) yields The 4 unknown parameters α, β, γ and ν are to be determined from ( 56), ( 61), ( 62) and (65).
The time evolution of the perturbation is governed by (64).In the special case where the homogeneous solution is in equilibrium, the associated driving force vanishes, i.e., f 0 (t) = 0 for all time, and so b(t) is constant.Then (64) yields (66) (to within a multiplicative constant).Therefore perturbations decay with time provided ν < 0. When the homogeneous solution is not in equilibrium, equation (64) leads to and therefore, (again to within a multiplicative constant), Recall from the discussion surrounding the phase plane in Figure (3) that V(f 0 (t)) approaches zero as t → ∞, and so again, C(t) (and therefore perturbations) decay when ν < 0.
Our goal is to study the sign of ν where stability demands ν < 0 for all perturbations.In principle, ν can be a function of, at most, λ 1 , G, κ, σ and k.However, observe that the stiffness κ and the surface tension σ do not appear in (61), and they enter as a pair, κ + σ k 2 , in (62) and (65).Thus ν can in fact be expressed as a function of λ 1 , Gk and κ + σ k 2 .The Pi theorem of dimensional analysis can now be used to infer that Once the expression (56) for h(y) is substituted into (61), ( 62) and (65), we have a system of 3 homogeneous linear algebraic equations for α, β and γ that we can write as: Necessary and sufficient for this system to have a nontrivial solution {α, β, γ } is that the determinant of the coefficient matrix vanish.This leads to an algebraic equation that is linear in N .Solving for N gives We now examine the sign of ν on the λ 1 , K-plane where, because our interest is in compression, we restrict attention to 0 < λ 1 ≤ 1.(In Appendix C we extend the results to the case of extensional stretches.)First, observe that the numerator of (68) vanishes on the curve where This corresponds to the thick solid curve in Fig. 5 which is monotonic, with K(λ 1 ) → ∞ as λ 1 → 0 + .Second, one finds that ν given by (68) is negative in the shaded region above this curve: There is another region of this plane on which ν < 0 corresponding to the left-hand side of the dashed curve in Fig. 5.The dashed curve corresponds to the vanishing of the denominator of (68) and so ν is unbounded on it.This region involves compression that is more severe than that associated with the gray region.Thus we shall not consider it from here on.Accordingly, consider points (λ 1 , K) in the interior of the gray region of the λ 1 , K-plane shown in Fig. 5. Since ν < 0 in this region, perturbations associated with this region decay with time.In general, the homogeneous solution corresponding to some value of λ 1 ∈ (0, 1), is stable for sufficiently large values of K, specifically for K > K(λ 1 ).Observe from the figure that the particular steady solution corresponding to λ 1 = 1 is stable for all K > 0.
It is important to emphasize that K defined in (67) 3 is not a simple scaling of the wave number.It is illuminating to write K as Fig. 6 The homogeneous solution corresponding to a value of λ 1 ∈ (λ cr , 1) is stable for all wave numbers k.The critical value of the stretch λ cr at instability is given by (75) where = σ κ/G 2 in which one term involves the wave number k nondimensionalized with the length-scale σ/G, and the other term involves the reciprocal of the wave number nondimensionalized with the length-scale G/κ.Consider the special case where the compliant ambient medium is absent: κ = 0. Then K ∝ k and so we conclude from the previous paragraph that the homogeneous solution is stable for sufficiently large wave numbers k.On the other hand, in the special case without surface tension, σ = 0, we have K ∝ 1/k.Thus in this case we conclude that the homogeneous solution is stable for sufficiently small wave numbers k.Thus we see that κ and σ work in complementary ways.
In order to study the sign of ν for various wave numbers k, observe that K is strictly positive and that its smallest value (as a function of k) is where the non-dimensional parameter is the ratio of the two length-scales G/κ and σ/G.The wave number at which K takes its minimum value is k = k cr := √ κ/σ .Figure 6 shows Fig. 5 again but now with the line K = K min added to it.Let λ cr be the value of stretch at which the bold solid curve K = K(λ 1 ) intersects the horizontal line K = K min .This critical value of stretch is given by the unique root in (0, 1) of the equation Collecting the preceding results, we know from (71) 2 and (72) that the homogeneous solution associated with a stretch λ 1 is stable against a perturbation with wave number k if κ/(Gk) + σ k/G > K(λ 1 ).Therefore it is stable against perturbations of all wave numbers provided the smallest value of κ/(Gk) + σ k/G over all k > 0 exceeds K(λ 1 ), i.e., provided K min > K(λ 1 ).From (73) and (75) this requires K(λ cr ) > K(λ 1 ) which, by the monotonicity of K(λ), is equivalent to λ cr < λ 1 ≤ 1; see Fig. 6.We thus conclude that ν < 0 for all wave numbers k provided λ cr < λ 1 ≤ 1 and therefore that the threshold for instability is λ 1 = λ cr .The corresponding critical wave number is k cr = √ κ/σ introduced below (74).Observe that as = σ κ/G 2 increases, (73) tells us that K min increases, and so from Fig. 6 we conclude that λ cr decreases.In fact, λ cr → 1 as → 0 and λ cr → 0 + as → ∞.
Stated differently, suppose the ambient environment offers no resistance to growth.Then, the stiffness κ → 0 and so → 0 according to (74) 2 and K min → 0 by (73).Then we see from Fig. 6 that λ cr → 1.Thus the only homogeneous solution that is stable against all perturbation is the one with λ 1 = 1.On the other hand if the ambient environment is rigid, then κ → ∞, → ∞ and K min → ∞ and so λ cr → 0 + .Thus all homogeneous solutions are stable in this case.The results of these two special cases coincide with what was found in [4].
Finally, we remark that λ cr need not exceed λ Biot .By substituting λ cr = λ Biot into (75) 1 we find that λ cr = λ Biot when ≈ 1.4196.For greater than this value λ cr < λ Biot and vice versa.In our two-dimensional setting, the surface energy per unit reference area, W(λ), is related to the surface energy per unit current area, w(λ), by W = λw.Therefore the term between the first pair of parentheses in the second integral can be written as W (λ) − σ .In view of the constitutive equation σ = W (λ) for the surface tension, see (8), the second integral in (B.6) drops out and we can write the dissipation rate as (B.7) The speed V R with which the referential boundary propagates is the volumetric flux of new material being added to the body.The factor multiplying it, i.e., is therefore the conjugate driving force.Observe that the explicit effect of the surface tension has dropped out, but it still does affect the driving force implicitly through the stress, see (40) and (46).
Fig. 7 The homogeneous solution corresponding to a value of λ 1 ∈ (λ − cr , λ + cr ) is stable for all wave numbers k.The critical values of the stretch λ ± cr at instability are given by (C.3) where = σ κ/G 2 separated from the shaded region by a region of instability; see the corresponding discussion pertaining to Fig. 5.
The minimum value of K is again given by (73) and the horizontal straight line K = K min on the λ 1 , K-plane intersects the bold curve at two points as shown in Fig. 7.The corresponding values of stretch λ ± cr are given by where was defined in (74).The discussion below (75) goes through again and we conclude that the homogeneous solution is stable against perturbations of all wave lengths provided λ − cr < λ 1 < λ + cr .The critical wave number, k cr = √ κ/σ , is the same for both stretches λ − cr and λ + cr .
Observe that when the stiffness κ → 0 we have → 0 and so K min → 0 whence λ − cr → 1 and λ + cr → 1; see Fig. 7.This agrees with the result in [4] concerning the case of a tractionfree boundary.At the other extreme when κ → ∞ one has → ∞ and thus K min → ∞ and therefore λ − cr → 0 − and λ + cr → 1/λ Biot .This agrees with the result in [4] pertaining to the case of growth on a smooth rigid surface.article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 2
Fig. 2 Homogeneous solution: (a) The region R R (t) associated with the reference configuration at time t , and (b) the corresponding region R(t) associated with the current configuration.Note that the corresponding boundaries S R (t) and S(t) have propagation velocities Że 2 and Ẋe 2 respectively