Shape Dynamics of Bouncing Droplets.

Oscillating shape motion of a freely falling and bouncing water droplet has long fascinated and inspired scientists. We propose dynamic non-linear equations for closed, two-dimensional surfaces in gravity and apply it to analyze shape dynamics of freely falling and bouncing drops. The analytic and numerical solutions qualitatively well explain why drops oscillate among prolate/oblate morphology and display a number of features consistent with experiments. In addition, numerical solutions for simplified equations indicate nonlinear effects of nonperiodic/asymmetric motion and the growing amplitude in the surface density oscillations and well agree to previous experimental data.

Introduction

Droplets are very peculiar systems mostly due to their shape dynamic properties[1]. They conserve the shape at the rest due to a surface tension, but are very deformable. Small perturbation of an equilibrium shape may induce large deformations in surface morphologies. One can trigger droplet inertial motion as well as morphological dynamics by inducing substrate surface oscillation[2]. Gradient of free substrate surface energy or angle difference between droplets leading and trailing edges provokes self-propulsion[3,4]. Droplets dynamical properties and their vibration modes may control their motion over substrate surface[5,6] and transiently defy gravity[7]. The contact with the substrate can be minimized by the development of super hydrophobic surfaces[8], allowing droplets to bounce on the surface like elastic balls[9].

Recent experiments demonstrate that an oscillating surface can launch wobbly water drop into the air at higher speed than it would launch a hard ball of the same mass[2]. Therefore, a synchronization between the internal vibration of a drop projectile and the frequency of the rising and falling surface can more than double kinetic energy of the droplets[2].

The behavior of bubbles differs from drops in gas-liquid systems. They display much broad range of shape deformations including turbulence[10], thickness variations[11-13], the Marangoni effect[14], draining[15], and ejection of droplets[16], rupture[17], self-adaptation[18], and chaotic behavior[19]. From this broad range of effects, numerical solutions of dynamic nonlinear equations for free thin fluid films explained thickness variations[20].

Despite of numerous experimental data[1-9], description of droplets shape dynamics remain analytically largely unsolved. Some dynamic effects in linear regimes have led to the classical wave equations[21,22]. Numerical simulations of Navier-Stokes proved to be effective tool in the path of describing shape motions of drops and bubbles[23-25]. Furthermore, numerical solution of one-dimensional Navier-Stokes equations showed that simulated shape dynamics agrees with experiment qualitatively[26]. In addition, combination of experiments and numerical simulations explained the shapes of singularities around the neck for a drop falling from a faucet[27].

In the letter, we apply closed, two-dimensional surface dynamics equations[28,29] and analytically address dynamic morphological patterns of droplets observed in experiments[1-9]. Our theoretical model explains why droplets can transiently defy gravity. Even though the results well agree with previous ones, the formalism we apply is new and different from previous theories. The purpose of the paper is to show the effectiveness of the differential variational surface (DVS) formalism in handling the shape dynamic problems and pinpoint existence of relatively simple dynamic solutions. We discuss shape dynamics of droplets, though due to generality of our dynamic solutions the analyses can be extended to wobbly dynamics of bubbles too.

Theory

To address shape dynamics properties of droplets, instead of using classically well-defined fluid dynamics (or we shall call it Navier-Stokes dynamics) we examine the surface dynamics of droplets. As far as the surface is closed and two-dimensional it has to follow to our generic equations of surface motions[29]. In this section, we briefly review fundamental principles of the equations. We base on the extension of differential geometry for moving surfaces. Because of the surface is moving, all parameters: base vectors (Vectors are designated by bold letters throughout the paper), the surface area S, the surface topology, the enclosed volume Ω are all time dependent functions.

Procedures for defining base vectors, metric tensor, and contravariant/covariant derivatives follow fundamentals of differential geometry. The base vectors are defined as partial derivative of the position vector  = ∂ = ∂/∂S, i = 1, 2, where S, are the surface coordinate for an arbitrary chosen point and the position vector.

The α = 1, 2, 3 component of the surface velocity is defined as V = ∂X/∂t (X is a space coordinate for the point on the surface). The surface velocity has the normal C and tangential V velocities, so thatwhere V = V1 + V2 implies Einstein summation convention throughout the text: repeated indexes in upper and down indexes is shortly written summation by the index.

Definitions of the surface velocity and the base vectors form fundamentals for defining Christoffel symbols leading to the definition of curvilinear derivatives and invariant time derivativewhere the symbol is defined for the moving surfaces as and B = ∇ is the curvature tensor. Definitions (1–3) along with the time derivatives of integrals form basics of calculus for moving surfaces and are necessary tools for derivation of dynamic equations[29].

As far as we have published exact derivation few times already[28,29], to avoid repetition but keep self-consistence, we provide generic equations for two-dimensional surface dynamics here and give only introduction to the derivation. The derivation follows from the minimal action principal of the Lagrangian:where u is a potential energy density of a field defined on the volume Ω and ρ is the surface mass density. Details about the surface mass density is given below in the results section. If there are n ‘micro’ fields lived on the volume, then due to additivity of potential energies

One may represent the potential energy density as a negative internal pressure applied on the surface, then and the Lagrangian can be rewritten as:

Trivial boundary condition dictates that in the absence of the shape dynamics the surface mass must be conserved, therefore the first boundary condition leads to:

(6) is continuity equation rewritten for curved surfaces, for derivations see refs[28,29].

According to the time derivative of the space integral theorem, a variation of the potential energy term is [28,29]. The normal term of the variation of the kinetic energy leads to and therefore

Here we used C = NV, where N is the α component of the surface normal, and applied the Gauss theorem to the left hand side of the integral equation, converting the surface integral to the space integral. The tangential term of the kinetic energy variation is and it must be equal to the term coming from the tangential gradient of the surface pressure, which can be modeled as (here p = pN)[28,29] and hence

The combination of equations (6–8) gives the set for closed, two-dimensional surface dynamics. The equation (6) displays how the surface mass density may change when the surface moves and (7 and 8) show the motion in the normal and the tangent directions. The equations are complete set for the surface dynamics, as far as have four ρ, C, V1, V2 unknowns and four differential equations. All information about how internal processes may effect on the surface dynamics is stored in the surface pressure term, which can be subject of modeling dependently on the nature of the problem. Because the Lagrangian (5) is invariant and the variation is taken by tensor calculus the equations are fully covariant.

Results

Before we start formal derivations, we should give a formulation of the problem. We assume following scenario: a freely falling droplet in a gravity is bouncing vertically from a super-hydrophobic substrate, where the substrate itself is freely rising and falling surface. Find a shape dynamics of the droplet if the tangential velocities compared to the normal velocity, interaction between the droplet and the substrate, internal friction and path instabilities are all negligibly smalls. Note that according to experiments internal friction has insignificant effect on shape dynamics of droplets[2]. To tackle the problem one may numerically solve the Navier-Stokes equations[23-25] that has already proven to be difficult (not to mention that the Navier-Stokes equations follow from simplifications of our surface dynamics equations)[28]. Instead, we use dynamic nonlinear equations for moving two-dimensional surfaces[28-30] (referred as DVS equations throughout the text) and solve it analytically, therefore producing very short solution.

For the surface modeling we assume following characterization: the droplet is continuum medium of molecules and has a surface boundary. The surface boundary is continuum medium too, but is made from fast diffusing layer. Formally, molecules from the droplet also diffuse with an environment, though the diffusion of surface molecules is quicker because they are in direct chemical contact with the environment. Technically, it means that the surface is transient layer phase separated from the internal droplet and the environment. Given this, we assign ρ to the surface mass density, while the droplet mass density is designated by ρ. Because the surface is continuum medium the surface mass density directly implies the thickness of the ‘diffusive’ layer. For simplicity, we assume that a friction between the surface and the droplet is negligibly small, so that the surface does not slip on the droplet. In other words, we assume single surface formalism, though generalization for many surfaces is not conceptually difficult and can explain buckling of the surface when the droplet evaporates. Though buckling is entirely different phenomena and we do not address it in this paper. Note that invoking many surfaces formalism would also provoke concept of internal friction directly linked to a visco-elastic effect. Here, we ignore viscosity of the droplet because for water, according to experiments, viscosity has low or no effect on shape dynamics of water drops[2]. Therefore limitations by single surface should be sufficient, though it might not be enough for description of more viscous drops.

Dynamic system: a surface

The formalism leading to the derivation of the differential variation of the surface (DVS) equations[28,29], for freely falling droplets, generalizes Eulerian representation of fluid dynamics and in contrast to Navier-Stokes, as it is demonstrated in this paper, is analytically solvable. The formalism is fully covariant and its analytic as well as numerical solutions qualitatively exactly reproduces the surface motion.

As far as the body falls freely in the gravity, we need to add gravitation to the equations of the two-dimensional surface motion and therefore, the DVS equations[28-30] should be rewritten accordingly. The Lagrangian of the motion reads:where ρ is the surface mass density, is the surface velocity, p is an internal pressure across the surface, ρgh is a hydrodynamic pressure applied by the gravity and ρ is the drop mass density. Note that the hydrodynamic pressure is added by negative sign because it is external one.

The variation of the Lagrangian (9) modifies the DVS equations so that the gravitation is taken into account. As far as DVS equations are already derived[28,29], instead of taking brute mathematical steps, we just mention how gravitational term ρgh modifies final equations that reads:where p = (p − ρgh)N, α = 1, 2, 3. As we have already mentioned above the first equation is a generalization of the conservation of mass, the second and the third equations display a motion in the normal and the tangent directions of the closed surface, and all satisfy conservation of mass and energy. is the curvilinear, invariant time derivative; ∇ is the curvilinear derivative; C is the surface normal velocity; V is the component of the surface tangent velocity; V stands for ambient component for the surface velocity; N is the component of the unit surface normal; B is the curvature tensor and is the mean curvature. The surface velocities are illustrated on Fig. 1. Repeated indexes indicate the Einstein summation convention. The equations (10) are covariant and are valid not only for capillary surfaces, but for moving surfaces of molecules too[28].

Figure 1

(a) Schematic drawing of the z = sinR/R surface and the illustration of the surface velocity with the normal and tangential components. The surface velocities C, V1, V2 for the point at the top are shown by arrows. The surface z = sinR/R is taken because it is reminiscent of the liquid substrate dynamics when the falling drop jumps. (b) Graphical illustration of a droplet sitting on a super hydrophobic substrate. Equilibrium shape of a droplet in a rest is a sphere.

Note that the equations (10) contradict Newton laws. According to Newton laws if there is no potential energy acting on the surface (no force), then the surface is in the rest or has constant surface velocities in the inertial reference frame. In contrast, what we find is that if one removes potential fields, i.e. sets ρgh and p naught in the equations, then one should expect that a solution to the equations must be constant normal C and constant tangent V velocities. It is easy to check that the constant surface velocity indeed satisfies the equations of motions (10) for freely floating droplets, but they are not the only solutions. Therefore, the DVS equations (10) predict that freely floating water drops will continue shape dynamics (if they were moving a priori) with retaining non-linearity even in the case when no potential field (no force) acts on them.

Simplifications

(10) remains numerically unsolved. We give simplifications in this subsection and solve simplified equations numerically and analytically. We assume that the tangent components of the surface velocity is infinitesimally small compared to normal velocity and hence the tangential gradient of the surface pressure can be modeled as negligibly small:

This approximation is acceptable because in the most cases the surface pressure acts in the normal direction so that the pressure gradient in the tangent directions is zero. As we have already shown in our previous works[28,29], the second equation significantly simplifies if the surface is homogeneous and can be described with time invariable surface tension σ, then

The equation (12) is commonly known as the dynamic fluid film equation for the surface normal motion and has been reported before[20,31]. Using (12) in (10), we end up with

The equation (13) is a solution for the second equation indicating the surface motion in the normal direction. We now assume that deformations along tangent directions compared to normal ones are negligibly small, then equations (10) with condition (11) simplify as:where (16) comes from (12) with assumption that tangent velocities are infinitesimally small. Note that, in this case, the (15) is the solution of (14, 16).

When the droplet touches the ground (the substrate surface), for very short interval of time, the gravity becomes compensated by a substrate surface reaction forces, so that, the drop instantaneously comes in equilibrium with the substrate and passes instantaneous equilibrium shape. The equilibrium shape satisfies conditions: C = 0, ∂V = 0 and ∂(p − ρgh)/∂t = 0, then solutions to (14, 15) are:

The equation (17) dictates that water molecules are homogeneously distributed on the surface, and (18) shows that the shape adopts constant mean curvature. Water droplets are closed and compact surfaces, therefore by the Alexandrov theorem[32] the shape with constant mean curvature must be a sphere. Therefore, instantaneous equilibrium shape of the droplet as well as the droplet sitting on a super hydrophobic substrate is indeed perfect sphere[1].

Infinitesimal linear models

Before we proceed further, note that (14–17) are consistent with existing infinitesimal models[21,22]. Indeed, let ρ0 be equilibrium surface density and assume that both C and are infinitesimally small, so that linearized conservation of mass (14) reads ∂ρ/∂t = 0 with the solution ρ = ρ0. This indicates that for small enough C, the density of the diffusive layer remains constant and infinitesimal models, in the limit of small oscillations and small mean curvature, are consistent with our framework.

Furthermore, if the mean curvature is the time dependent function, then according to calculus of moving surfaceswhere Δ = ∇∇ is the surface Laplacian[31]. For surfaces with vanishing mean curvature, BB becomes negative twice the Gaussian curvature K [20]. Therefore, by differentiating (16) linearized acceleration becomeswhere γ = σ/ρ0. In the limit of planar cut off, i.e for flat equilibrium configurations, the Gaussian curvature is zero K = 0, and infinitesimal deformations are governed by the wave equation:

With the sign convention, the equation (20) is exactly a case for linear models and is consistent with previous results[21,22].

Dynamic solutions

Now we proceed with extracting approximate analytic solutions for (14–16) equations and give numerical solution for (14, 16) equation set.

After the droplet spends some energy on the surface deformation, taking into account that it found instant equilibrium shape it jumps again, but since it has already found equilibrium, which is a sphere, it will start fluctuations around the equilibrium shape. Since, the shape oscillates near to equilibrium; we can suggest that according to (18) the term becomes infinitesimally small constant, therefore (15) can be modified as

Since cannot be zero and its combination with the surface and hydrodynamic pressure is quasi infinitesimally small constant, then the solution to (21) is

This is the continuity condition in the Navier-Stokes equations for incompressible fluids[33]. However, setting as infinitesimally small constant is an approximation and therefore, (22) is not an exact solution. (22) shows that zero divergence of the surface velocity is the near equilibrium solution of the DVS for water drops under the gravity. (22) can be trivially handled by the  = k/R3 function (where k is some constant and is the position vector), which corresponds to the sphere in equilibrium stationary case. Therefore, in the (0 ≤ θ ≤ 2π, 0 ≤ ϕ ≤ π) spherical coordinates, assuming designation  = (sinϕ sinθ, sinϕ cosθ, cosθ) the solution reads:where ω, ξ = x, y, z are some frequency like functions and for the sake of simplicity we consider them as constants. is the base vector of and R is the projection of the position vector in the ξ direction (for example  = (sinϕ sinθ, 0, 0),  = (0, sinϕ cosθ, 0),  = (0, 0, cosθ) and R = R + R + R). Trivially, (22, 23) leads to the solutions: for some R constants with initial condition  =  =  = R0 and R = R (t) time variable function, reading respectively:where is defined as the vector components of which are exponents of frequency times time and the spherical coordinate: . (Figure 2a–e) are a graphical representation of (24, 25) solutions. In addition to (24, 25), according to the vector calculus the analytic solution to (22) iswhere vector has the same dimension as angular velocity. Suggesting that should behave same way as , with conserving full generality, the solution (26) can be rewritten aswhere v0 is some constant wave propagation velocity. After taking the curl of (27) and using the curl of the curl (∇ × ∇) identity one ends up with a wave equation for the position vector:

Figure 2

Illustration of analytically solved shapes. From left to right: (a) sphere with unit radius  = 1 (initial condition of (24) solution), (b) prolate shape drop  = (sinϕ sinθ, 0.86sinϕ cosθ, 1.75cosθ) (solution (24) at some given time), (c) oblate shape  = (sinϕ sinθ, 1.72sinϕ cosθ, 0.75cosθ) (solution (24) at other moment), (d) non-prolate/oblate shape  = (esin, esin, ecos) (solution (25) at some given time), (e) summation of (a–d) shapes indicate development of continues singularities (left hand side of the droplet). Panel (f) displays plot of analytic solution and indicates how spacial coordinates (x, y, z) non-periodically, asymmetrically and with growing amplitudes fluctuate in time. The graph also indicates that (z) direction (black graph) diverges from (x, y) rapidly while (x) (blue graph) separates from (y) slowly, meaning: the drop will more rapidly develop singularities in (z) direction while will keep stability in (x, y) directions for relatively long time.

According to curl of curl identity ∇ × (∇ × ⋅) = ∇(∇⋅)−∇2⋅, applying this to (27) one obtains . Taking into account (22) ∇ = const, therefore first term vanishes and one ends up with (28).

Taking into account that as the time evolves all spatial directions become linearly independent R = ψ (ξ, t) (where ξ = x, y, z and ψ is some function), then (28) transforms as one dimensional wave equation . Applying boundary ψ(0, t) = 0, ψ(L, 0) = 0 and initial conditions ψ(ξ, 0) = f(ξ), ∂ψ/∂t(ξ, 0) = g(ξ) one obtains a solution

Taking into account that at the initial condition (t = 0) =  the water drop is a sphere with the radius R0, then the solution can be written as

The solution (29) precisely explains oblate-prolate oscillation of the drop observed in experiments. (Figure 2f) shows dependence of the position vector components (x, y, z) on the parametric time and is plotted based on the linear combination of the (24, 25, 29) solutions. The graph is plotted for L = 1, f(ξ) = ξ parameters on the ϕ = θ = π/4 spherical coordinate and is designated to show phenomenological dependence of the point spacial coordinates (x, y, z) on the parametric time. The solution clearly indicates non-periodic and asymmetric behavior in all directions. z direction diverges from x, y rapidly while x, y diverges slowly from each other. Consequently, the drop more rapidly develops singularities in one direction while keeps stability in others for relatively long time. After long enough time, there is no overlap in amplitudes resulting in drop division.

We have also solved (14, 16) equations numerically by the 4th generation Runge-Kutta method for 0.0001 time lapse and interval of 30 nominal parametric time. The numerical solution in (ρ, C) phase space (Fig. 3) displays nonlinear thickening with growing amplitudes of diffusive layer, once again indicates nonlinear effects of nonperiodic/asymmetric motion and the growing amplitude in the surface density oscillations, and qualitatively well agrees to experimental data[1-19].

Figure 3

The numerical solution in (ρ, C) phase space of (13, 15) equations reminiscent to the analytic solutions of the (14) equation. Black, blue and red lines are plotted for different amplitude values of C = fsin(πt): f = 0.5 black, f = 1 blue and f = 2 red. The initial condition for the interface velocity C is set to be zero C(0) = 0 and the initial value for the surface mass density for initially flat configurations is fixed ρ = 0.5. The mean curvature is set to be one corresponding to the sphere with unit radius. (a–c) Panels display simulations for three different tensions.

In addition, to indicate how different values of surface tensions may effect on the surface density oscillations, we have plotted (ρ, C) phase space for three different surface tensions (see Fig. 3). The initial value of the surface tension (Fig. 3a) σ = 0.073N/m is set to be the standard value for water drops. (Fig. 3b) is plotted for the surface tension of the inkjetprinted picoliter droplet σ = 0.063N/m and (Fig. 3c) displays oscillations for the microliter droplet with tension σ = 0.035N/m. The different values for the surface tensions is taken from the ref.[34]. (Fig. 3) shows that decrease in the surface tension decreases oscillation amplitudes and increases frequency.

Since (24, 25, 29) are particular solutions, then their linear combination is a general solution. Therefore, it predicts non-linear oscillations among spherical, oblate/prolate and non-oblate/prolate shapes, tendency towards increasing the radius, and inducing shape instabilities by forming singularities, and with combination of the numerical solution to (14, 16) indicates the growing amplitude in the oscillations of the surface mass density, meaning increasing the mass instabilities in the diffusive layer of the surface (Figs 2f and 3). Mass instabilities in the diffusive layer and linear combination of the (24, 25, 29) solutions ultimately lead to a development of singularities in the drop, which may induce drop division.

Relevance to experiments

To explain shape dynamics of axisymmetrically deformed drop Rayleigh[35] proposed expansion of absolute value of radius vector R in Legendre polynomials L so that where θ is the polar angle, a is the time dependent amplitude coefficients and n is so called the mode number. Later on Lamb[36] suggested that for a freely oscillating drop with small amplitude of oscillation the coefficients can be estimated aswhere ω = 2π/T is an oscillation eigenfrequency with decaying period T, ϕ is the phase angle, b is the eccentricity of the drop, which may arise from oblation of the droplets due to hydrodynamic effects, β is the dumping rate of the n mode. We refer to this model as Rayleigh’s exponentially decaying model. The model was successfully used to describe recent experiments[34,37-40].

The eigenfrequency ω depends on the liquid density ρ and the surface tension σ, therefore a = a(ρ) and hence in the first approximation . On other hand, in the first approximation a is proportional to the surface interface velocity C, therefore (ρ, C) phase space has to mimic the Rayleigh’s exponentially decaying model.

Indeed, for comparison to previous experimental works we mention that in our numerical calculations for interface velocity amplitudes we used f = asin(πt). However, instead of harmonic f if one uses Rayleigh’s decaying amplitude model, which can be expressed as f(t) = f(0)exp(−t/T) sin(2πt/T + ϕ) + b[35,36], then one should expect that (ρ, t) would mimic overall morphology of experimentally obtained decaying amplitude oscillation of droplets mean radius. Indeed, as the (Fig. 4) shows plotted graphs for different f: f(t) = 0.5exp(−t/5) sin(πt + π/4) red, f(t) = exp(−t/5) sin(πt/2 + π/4) green and f(t) = 2exp(−t/2) sin(πt/2 + π/4) blue curves indicate similarity of overall morphology of numerically calculated (Fig. 4a) and experimentally observed (Fig. 4b)[34] graphs. Note that in the (Fig. 4a) the surface mass density is modeled as ρ = ka(t)/R0, where k is a coupling constant and is set as unit in the numerical solution, a is the oscillation amplitude of mode n and R0 is the initial radius of the spherical droplet. The numerical solution for the (a(t)/R0, t) space were obtained by ode45 numerical solver. These calculations do imply that if one carefully picks and fits parameters in theoretical (24, 25, 29) and numerical calculations of (14, 16), then one will reproduce experimentally observed dependence of oscillation amplitude on the time. For clarifications that are more rigorous it is necessary to be made theoretical calculations, experimental observations and parameter fitting simultaneously, which is not possible now but is left out as a source for future collaborations for experimental suggestions and clarifications.

Figure 4

(a) The numerical simulation of the dependence of amplitudes on time. (a) Red, green and blue curves correspond to different f(0) interface velocity amplitudes and produce remarkable morphological similarities with (b) experimentally observed decaying oscillation amplitudes of water drops. The panel (c) provides photographs of micro sized droplets shape oscillations at different microseconds. Experimentally observed shapes well agree to our analytic and numerically solved shapes. (b,c) Experimental data is obtained from the ref.[34].

Discussion

Despite neglecting uniform translational motion of the drop (one may argue that for sufficiently small periods of time the inertial effects on drop shape motion is negligibly small), we might still speculate why oscillating surfaces can lunch water drop at higher speed, than hard sphere[2]. Due to the generality of the DVS equations, an oscillating surface undergoes the same shape motion at the interface of substrate/droplet as water drops. If frequencies of the substrate surface and water drop matched, within an order of magnitude, then a resonance effect takes place and the drop will be launched by higher speed. We should also note that according to (12–18) and (21) equations gravitational term ρgh has no effect on governing shape equations. Gravity has no role therefore surface motion would be the same if one would remove gravity. This explains why water drops can move against gravity[7] on oscillating super hydrophobic substrate.

To conclude we have proposed a system of nonlinear equations to describe the nonlinear features of the dynamics of two-dimensional surfaces, with large deformations and large variations in density, both spatial and temporal[28-30]. Analytic solutions to simplified version of two-dimensional surface equations in gravity describe freely falling water drop’s shape dynamics and precisely explain why the shape non-linearly oscillates among prolate/oblate and non-prolate/oblate forms and display a wide range of shape instabilities. Surprisingly, despite the fact that the proposed system purposefully disregards the tangential components of surface velocities, the solution qualitatively explains experimentally observed shape dynamics[1-9]. By this analytic solution, we have reproduced numerically solved Navier-Stokes equations results[24], but avoided intensive modeling for the surface pressure, which would be necessary for numerically solving Navier-Stokes. We have shown that the continuity equation is an approximate solution of the DVS equations, which, in contrast to Navier-Stokes, can be trivially handled in this particular case. In addition, we have shown that linear models leading to wave equations[21,22] are generally correct for near to equilibrium approximations.

The numerical solutions of (14, 16) combined with analytical ones (24, 25, 29) directly imply that: the surface mass density is non-linearly increasing while the shape adopts diverse forms. Impact of this statement in biology is that: if one starts formation of cell from a mixture of organic molecules in a droplet and lets evolution of the shape dynamics by applying some potential field (like hydrophobic-hydrophilic interactions for instance), then the droplet will ultimately develop some membrane like structure (because surface mass density is nonlinearly increasing) and will adopt rich divers shapes. Such active dynamics may explain why cells got membranes[41] and is in contradiction with generally thought idea that membranes were formed first.