Elasto-buoyant heavy spheres: a unique way to test non-linear elasticity

Extra-large deformations in ultra-soft elastic materials are ubiquitous, yet systematic studies and methods to understand the mechanics of such huge strains are lacking. Here we investigate this complex problem systematically with a simple experiment: by introducing a heavy bead of radius $a$ in an incompressible ultra-soft elastic medium. We find a scaling law for the penetration depth ($\delta$) of the bead inside the softest gels as $\delta \sim a^{3/2}$. While this result is inconsistent with an ideal neo-Hookean model of elastic deformation, according to which the displacement fields must diverge, it is vindicated by an original asymptotic analytic model developed in this article. This model demonstrates that the observed relationship is precisely at the demarcating boundary of what would be required for the field variables to either diverge or converge. This correspondence between a unique mathematical prediction and the experimental observation ushers in new insights into the behavior of the deformations of strongly non-linear materials.

of the containers are grafted with a thin layer (∼5 nm) of polydimethyl siloxane chains so that the gel solution contacted the walls at 90 o to ensure that the surface of the cured gel is flat. All the experiments are performed after two hours of gelation. For the estimation of the shear modulus of the gel, we use a linear elastic model in order to ensure consistency of quantification for all the gels. The shear modulus is determined from the resonant mode of vibration of a gel slab confined between two parallel glass slides [22]. Steel spheres (density 7.8 g/cc, diameters 1 -10 mm) are gently placed one by one on the gel surface and its side-view image is captured by a camera. Dissipative processes within the gel dampen any oscillations, and the spheres sink until they become stagnant in the polyacrylamide gel, the whole process taking only fractions of a second. The depth of submersion, δ, is measured from the upper surface of the gel till the base of sphere, that denoting the net downward displacement due to the inclusion of the spherical particle by the surface (Fig.1). The cells are large enough to avoid any finite boundary effects, including the side walls and the bottom. The measurement of depth for each sphere is made in the central region of the container. After each measurement, the sphere is gently removed from the gel using a magnet, held slightly away from the free surface. We wait for a few minutes between each measurement that ensures there is no memory of the position of the previous sphere inside the gel. If the bead is too small and the gel is stiff [7,23], the surface bends slightly under the weight of the bead and δ a ( Fig.1-a). By increasing the bead radius or decreasing the elastic modulus, the particle submerges itself to a considerable depth inside the gel. The surface of the gel wraps around the particle and closes to create a line singularity connecting the particle to the free surface of the gel ( Fig.1-b). Strings of tiny air bubbles appear in this thin channel, which soon coalesce and escape through it while the channel further closes due to the auto-wetting forces of the gel's surface. If the surface of the gel is premarked with ink spots, it is easy to visualize that the surface of the gel becomes appreciably stretched while the sphere sinks through the gel while being still connected to the free surface via a thin channel. It is also possible to release ink inside the gel in the form of thin vertical lines with the help of a fine needle, which bend toward the sphere in a dramatic way when the sphere is released inside the gel demonstrating a substantial amount of tensile strain developed in the gel network due to the inclusion of the bead. These basic experiments were reported in a previous article [14], but without a detailed analysis. Here we report a detailed set of experiments, in which the shear modulus of the gel was varied from 13 Pa to about 3000 Pa.
Prior to subjecting these experimental results to a comparison with a theoretical analysis, we needed to verify how meaningful it is to consider only the effect of elasticity by ignoring the surface tension of the gel in predicting the depth of submersion of the bead and how reversible is the deformation of the gel. First question is partly philosophical that rests upon the distinction between surface free energy γ and surface tension. The latter differs from the former by surface stress dγ/de, e being the surface strain. As the major constituent of these amorphous gels is water, we expect that the surface stress is negligible. In fact, several recent studies that measured surface tensions of various amorphous soft polymers strongly suggest that their surface tensions are practically same as their surface free energies [24]. Thus we need to figure out only if the surface free energies of the gels play any role. Part of the answer can be obtained by comparing them with the equivalent spring constant of the sample. The latter can be obtained by slightly raising the height of the bead by an electromagnet and releasing it so that the bead undergoes an under-damped oscillation and reaches the neutral position. For three gels, in which the beads were completely submerged, the equivalent spring constants were estimated from the frequency of oscillation to be 0.2 N/m, 5.6 N/m and 13 N/m, which increase systematically with their shear moduli (13 Pa, 140 Pa and 360 Pa). Fig.2-a shows a typical profile of such an oscillation. Comparing these spring constants with the surface tension (0.07 N/m) of water, we conclude that the contribution of surface tension can be safely neglected for all the gels used in this study except, perhaps, for the lowest modulus gel for which the spring constant is three times that of the gel's surface tension. However, when a bead is completely submerged in the gel, any variation of the height of the sphere does not alter the area and the excess energy of the free surface of the gel. We thus believe that the surface tension spring does not play a significant role in determining the depth of submersion of the sphere as long as it is completely engulfed by the gel. Further support to this viewpoint, i.e. the dominant role of the elasticity over surface tension was gathered from the experiment described below that also exemplified the reversibility of the softest gel (µ ∼ 13 Pa) employed in our experiments.
By varying only the system temperature, the elastic modulus of polyacrylamide gels changes, which increases with temperature. A thin layer of Paraffin oil is poured over the gel to avoid its surface from drying. The temperature of the gel is monitored by placing a thermometer inside an identical sample of gel in a similar sized container, placed inside the oven. After the gel is heated to 70 o C, a 5 mm diameter steel sphere is released into it through the layer of paraffin oil. The depth of the sphere is measured at this temperature while it is in the oven. As the gel is gradually cooled, the sphere sinks deeper inside it. We wait for an hour between acquiring data for the depth of sphere at each temperature to allow the gel to equilibrate reasonably well. After the gel is cooled to about 5 o C, it is heated again that decreases the depth of the embedded steel sphere (Heating Cycle, Fig.2-b). The depth of submersion of the sphere plotted as a function of the temperature for both the cooling and the heating cycles ( Fig.2-b) shows that there is a little hysteresis in this system, in that the difference in the depths of the sphere for a given temperature is within 5%. We conclude that the deformations of the gel generated by a bead are predominantly reversible and controlled by elastic forces. The penetration depth, δ, is plotted versus the radius a of the steel spheres in log-scales in upper inset of Fig.3, for ten different gels (µ : 13 Pa -2930 Pa) and different radii. The upper white region of this inset shows the data points for beads that were completely below the surface of the gel. When we non-dimensionalize the depths as well as the radii of each bead in each gel by the material lengthscale δ 0 defined as δ 0 = µ ∆ρg (Fig.3), where µ is the shear modulus of the gel, ∆ρ is the effective density of the buoyant spheres, and g is the gravitational acceleration, we see that all the data cluster around a mean master curve with two distinctive asymptotic limits. For each sphere-gel system, ∆ρ was precisely estimated from the image on the basis of the Archimedes principle of buoyancy analogous to the scenario for floatation on a liquid [25]. The normalized data can be divided more or less into two regimes: one with the non-dimensional radii a/δ 0 < 1 and other one with a/δ 0 > 1 with an intermediate transition regime. We fit the non-dimensional depths as a function of the non-dimensional radii for the regime a/δ 0 < 1 with the power law function (δ/δ 0 ) = k(a/δ 0 ) α with adjustable parameters, α and k and find α = 1.96 ± 0.06 and k = 1.09 ± 0.1. The error bars are obtained from a 95% confidence limit analysis. We conclude that for a/δ 0 < 1, the depths as a function of the radii follows δ ∼ a 2 within the error limits. On careful examination of the experimental points corresponding to large deformations, i.e. a > δ 0 , we find that the data corresponding to the two softest gels (µ : 13 Pa and 25 Pa) are shifted from the rest of the data due to their multiplicative factor (k) being significantly larger than the rest of the data for the other gels. This indicates that the lengthscale δ 0 , which is defined with the elastic linear properties of the sample, is therefore not sufficient to describe the whole data accurately, and non linear effects have to be taken into account. We conclude that all the data in this regime (a/δ 0 > 1) cannot be investigated together. It is more appropriate to investigate the data for each gel composition separately i.e, for a given non-linear material stress-strain relationship. By fitting the data for the two softest gels (µ : 13 Pa and 25 Pa) where all the beads are completely engulfed, with the power law function, we find that α 13 P a = 1.42 ± 0.05 (see lower inset of Fig.3), k 13 P a = 2.52 ± 0.3 and α 25 P a = 1.52 ± 0.09, k 25 P a = 1.63 ± 0.3. The exponents for the fits for the gels (a/δ 0 ∼ 1) in the intermediate regime lie between α ∼ 1.5 and α ∼ 2. Thus, from the experimental observations, we infer that in the limit of a significantly greater than δ 0 , the general trend is close to δ ∼ a p , where p is in the range of 1.4-1.5, within the error limits. Thus, we conclude that the power law observed (δ ∼ a 2 ) in the gels of higher shear moduli is closer to the regime already studied before what one would expect with an analogy with the elastic Stokes equation in the Hookean limit. What is astounding is the observation of δ ∼ a 1.5 in the case of the gels where δ > 2a for each tested bead radius. In order to interpret these novel observations in the ultra-soft gels, we develop a new model to tackle such extra-large deformations in the following part of the paper.

II. ASYMPTOTIC ANALYTIC MODEL FOR LARGE ELASTIC DEFORMATIONS
From the outset, a motivating picture of the problem can be gleaned from the comparison of the potential energy of the bead in the gravitational field and the energy of the elastic deformation of the gel. In the limit of Hookean deformation, an analogy with the Stokes equation suggests that the elastic energy is of the form µδ 2 a, which is to be compared with the gravitational potential energy ∆ρga 3 δ. This leads to the scaling: δ ∼ a 2 . However, the experimentally observed scaling δ ∼ a 3/2 implies that the elastic energy in extremely large deformation must scale as δ 3 . The detailed analysis based on a model, presented in this paper for the first time, shows that the above scaling is non-trivial. Moreover, the scaling result that follows is independent of the constitutive laws of an elastic material.

A. Gravity energy of engulfed spheres
Consider deformations that preserve the volume, as is the case with the elastic gels used in the experiments. Since the bead is supposed to be totally engulfed with δ a, the surface of the deformed gel is fairly flat and horizontal. A supplementary (virtual) downshift δ produces an opposite rise of the same volume of gel. Therefore the gravitational energy variation of the system consisting of the sphere plus the gel is 4 3 πa 3 ∆ρgδ . We conclude that the gravitational energy shift is E gr 4 3 πa 3 ∆ρgδ in the limit we are considering.
B. Outline of calculating the elastic energy in the limit δ a The axis of symmetry being the vertical axis, the two coordinates changed by the deformation are the radius in the horizontal plane, r, and the vertical coordinate, z. The deformation maps the undisturbed state with coordinates (r, z) to a disturbed state (R(r, z), Z(r, z)), or, in radial coordinates from coordinates (r; θ) to (R(r, θ); Z(r, θ)) withr radius and θ polar angle of the (r; z) plane (see Fig.4-(d)). The elastic energy is a function of the strain tensor (also called right Cauchy-Green deformation  (r, z).r is the distance from the initial contact point of the bead. In the deformed state, the point that was at (r, z) is located at R(r, z), Z(r, z).
tensor), which gives the square of local change in distances due to deformation: C = F T F , F being the deformation gradient tensor [26]. In the absence of any preferred direction, the elastic energy of an incompressible solid may depend only on two scalars (invariant under global rotation) that can be made out of the strain tensor: the trace (I 1 ) and the sum of the square of its components (I 2 ). In cylindrical coordinates with an azimuthal invariance, these invariants read I 1 = R 2 ,r + R 2 r 2 + R 2 ,z + Z 2 ,r + Z 2 ,z and I 2 = (R 2 ,r + Z 2 ,r ) 2 + (R 2 ,z + Z 2 ,z ) 2 + 2(R ,r R ,z + Z ,z Z ,r ) 2 + R 4 r 4 , where indices preceded by a comma denote partial derivatives. For the neo-Hookean model, which is commonly used to describe soft gels, W = µ 2 I 1 . However, beyond a certain deformation, the neo-Hookean model cannot be a fair representation of rubber-like materials whose elasticity originates from unfolding of polymer chains: once these chains reach their full extension, the energy cost for a supplementary unfolding diverges so that a further deformation is accompanied by a diverging additional elastic energy, in contrast with the ideal neo-Hookean law (see Fig. 6). In what follows, we do not restrict ourselves to the neo-Hookean case.
The deformation field of the gel can be described with two characteristic lengths, the downshift of the sphere δ, and the radius of the sphere a. Here, we are dealing with the limit of large downshifts, that is δ far larger than a. Below we demonstrate that in this asymptotic case, the elastic energy of the system does not depend on a in the limit a δ. As a first step we consider a normal point-load applied on the free elastic surface at (r, z) = (0, 0) (Fig 4-(b)) and we find a scaling law for the elastic energy per unit volume at a given distancer from the loading point. The condition for the elastic energy to be convergent is established. Then as a second step we consider a hard sphere with a finite radius, producing a displacement δ at the contact point (r, z) = (0, 0) equal to the displacement at (r, z) = (0, 0) induced by the previously considered point-like load. We demonstrate that for δ a, the elastic energy associated with the displacement field generated by the sphere or the point-like load are equal: These are independent of a and proportional to δ 3 .

C. Elastic energy for a point-like load
Incompressibility of the elastic medium is imposed by writing that the determinant D of the first derivatives of R(x, y, z) is equal to one. For axisymmetric deformations, D = R r (R ,r Z ,z − R ,z Z ,r ). The total energy of the system {gel + sphere} reads : where the Lagrange multiplier q(r, z) imposes the incompressibility condition D = 1 [27]. For the point-like heavy sphere we first consider (see Fig. 4-(b)), the weight the bead can be reduced to a point force located at (r, z) = (0, 0). Eq. 1 simplifies into: where δ 3 (r, z) is the 3D Dirac distribution. The Euler-Lagrange conditions of minimization of the energy (Eq. 2) read [5]: Let δ be the displacement of the gel at (r; z) = (0; 0). δ is the unique relevant length for the displacement field. For points located at a distance to the loadr much smaller than δ, δ is no more relevant at these short length-scales. Therefore no length-scales are expected to occur in the scaling laws and one assumes a power law ofr for the displacements forr δ: Z = δ +r β f 1 (θ) and R =r γ f 2 (θ). Since the gel is vertically stretched and horizontally squeezed in the vicinity of the bead one assumes that γ > β. This hypothesis will be checked afterwards (see Eq. 16 below).
We assume that the strain energy density function is W ∼ I α1 1 I α2 2 (the exponents α 1 and α 2 are constant) in the range of strains undergone by the elastic solid in the vicinity of the application point, i.e. forr far smaller than δ. This choice for W is crucial to obtain the scaling law, but it has no effect on the final results, as argued in section II E: it does not limit the general nature of the theory.
The scaling laws for the first and the second invariants of the Cauchy deformation tensor are: These scaling laws for I 1 and I 2 are used to obtain the scaling laws for the various terms of the Cauchy-Poisson equations (Eqs. 3 and 4). We first deal with Eq. 3: In the same way one finds: From Eq. 3 and Eqs. 8, 9 and 10 one obtainsr , and the scaling law for the Lagrange multiplier: The last Cauchy-Poisson equation (Eq. 4) is now used to get an expression of β as a function of α 1 and α 2 . It writes: Denoting as σ zr = ∂(W−qD) ∂Z,r and σ zz = ∂(W−qD) ∂Z,z , the two quantities present on the left-hand side of Eq.12, one finds that this left-hand side is the divergence (expressed in cylindrical coordinates) of the z-components of the stress tensor σ [26]. We first calculate the scaling laws for σ zr and σ zz using Eqs. 5, 6, 11: The second term of the right-hand side of Eq. 14 scales asr (β−1)(2α1+4α2−3)r1−β . It is negligible whenr tends to zero when compared to the first term of Eq. 14 provided that the exponent of the second term is larger than the exponent of the first one, i.e.: which is formally equivalent to β < 1. Since it is assumed from the beginning that β < 1, we conclude that σ zr ∼ r (β−1)(2α1+4α2−1) and in the same way, σ zz ∼r (β−1)(2α1+4α2−1) : σ zr and σ zz follow the same power law whenr tends to zero. The right-hand side of Eq. 12 being a delta-like charge density, we conclude from Gauss's theorem that the scaling law for σ zr and σ zz is also σ zr ∼ σ zz ∼ 1 r 2 . Therefore,r (β−1)(2α1+4α2−1) ∼r −2 , and: β is an increasing function of (α 1 + 2α 2 ). It is negative for α 1 + 2α 2 < 3/2, and it is always smaller than 1. The displacement at the application point is finite if the exponents are positive, i.e. if α 1 + 2α 2 > 3/2. Otherwise the material cannot withstand a point-load. In the following we assume that α 1 + 2α 2 > 3/2, an ansatz that will be legitimized later.
Since δ is the unique length-scale of the deformation, the coordinates in the deformed configuration can be written as: R = δf 0 r δ and Z = δg 0 r δ , where f 0 and g 0 are two numerical functions depending only on the constitutive law of the elastic medium. The invariants I 1 and I 2 , and thus the elastic energy density are dimensionless functions depending only on r/δ. The elastic energy due to the deformation, (subscripts 0 are for the point-load problem) is therefore proportional to δ 3 with a coefficient C 0 , which has the dimension of shear modulus and depends on the mechanical properties of the elastic solid, as: Note that the convergence of the integral in Eq.17 is ensured since the work done by the applied force is finite for α 1 + 2α 2 > 3 2 , yielding the validity condition of Eq. 18.

D. Elastic energy for a finite sphere
We consider now a heavy bead of radius a, and denote δ its vertical downshift, which is supposed to be far larger than a (Fig.  4-(c)).
• The unique relevant length-scale forr a being δ, the displacement field reduces in this range ofr to the displacement field for a point-load with the same penetration depth δ: Z = δg 0 r δ , where g 0 has been defined in section II C.
• The power law for Z − δ with the exponent given by Eq. 16 applies in the intermediate range a r δ. Combining it with the previous expression of Z (valid forr a) yields: • Close to the bead (r δ) the unique relevant length-scale is a. Z − δ can be expressed as where g 1 is a numerical function and the constant K (with the dimension of a length) depends on a and also on the far field deformation, i.e. on δ. Comparing Eq. 19 with the general expressions for Z in the ranger δ, one obtains K = δ a δ β .
A similar expression for R can be obtained with the exponent (3 − β)/2, leading to negligible contributions for the elastic energy.
Forr δ, the first invariant I 1 scales as Z 2 ,z and Z 2 ,r , and the second invariant I 2 scales as Z 4 ,z and Z 4 ,r . The strain energy density function W ∼ I α1 1 I α2 2 can thus be written forr δ (using Eq. 20 with the expression of K) as W ∼ a δ β−3 g 2 ( r a ). In the ranger a, W W 0 , W 0 being the elastic energy density of the point-load problem (Fig. 5).
Moreover, in the ranger δ, one finds from the scaling laws of section II C: One concludes that, for anyr, W − W 0 ∼ a δ β−3 g 3 ( r a ), with g 3 a numerical function (depending neither on δ nor on a) whose limit as its argument approaches ∞ equals 0. Thus, we obtain: Since one assumes again that α 1 + 2α 2 > 3 2 , β > 0 and the difference of the elastic energy with the bead of radius a, E el , to the elastic energy with the point-load, E el 0 , is negligible with respect to δ 3 . One concludes from Eq. 18 that the elastic energy E el is proportional to δ 3 . Note that if α 1 + 2α 2 were smaller than 3/2 the near-field part of the elastic energy would not be negligible anymore compared to the far field contribution.

E. The strain energy density function
In sections II C and II D the strain energy density function has been assumed to scale as W ∼ I α1 1 I α2 2 for the strains encountered atr δ. It has been then demonstrated that the elastic energy corresponding to the downshift of a sphere of finite radius a over the distance δ is proportional to δ 3 in the limit δ a if α 1 + 2α 2 > 3/2.
It is worth considering the case of an elastic material for which the strain energy density function increases softer than I α1 1 I α2 2 with α 1 + 2α 2 = 3/2. Forr δ we demonstrated that Z − δ = δ a δ β g 1 r a and then Z ,z ∼ δ a 1−β , showing that the strain is arbitrarily large in the limit δ a , if β < 1, i.e. if α 1 + 2α 2 < 3/2. However, due to the finite maximum stretch of the polymer chains constituting the material, the maximum stretching of any real elastic rubber-like material is bounded. To this maximum stretching corresponds a divergence of the strain energy function. This divergence results in a steeper and steeper increase of W with I 1 , I 2 , or both, which is associated with increasing values of the exponents α 1 and/or α 2 [26].
In order to illustrate this last point, let us consider a particular strain energy density function, for instance the Gent hyper-elastic model [21,28]. The strain energy density function of this model has a singularity when the first invariant I 1 reaches a limiting value. It is plotted in Fig. 6: deviations from the initial neo-Hookean behaviour (α 1 = 1 and α 2 = 0) yield an increasingly stiffer and stiffer strain energy density function. At any value I 10 of I 1 one can define the local exponent α 10 = I1 W dW dI1 so that at the vicinity of I 10 the strain energy density function is W ∼ I α10 1 , locally. Beyond a certain finite value I * 10 of I 10 the exponent is larger that 3/2, ensuring that for any I 1 larger than I * 10 the condition α 1 + 2α 2 > 3/2 is fulfilled. Note that values of I 1 larger than I * 10 are automatically reached, otherwise the strain would diverge aroundr ∼ 0, leading to arbitrarily high values of I 1 (as explained above).
What is illustrated using the example of the Gent model is general, and can be applied to any elastic constitutive law involving I 1 and I 2 of a real elastic material: starting from low to moderate strains for which the neo-Hookean model is expected to apply far from the bead, energy density functions stiffer than I α1 1 I α2 2 with α 1 + 2α 2 > 3 2 are necessarily encountered next to the bead for any real elastic material.
A material for which the failure limit is reached before the strain energy density is stiffer than I α1 1 I α2 2 with α 1 + 2α 2 > 3/2 would not be able to sustain the heavy sphere, and thus would be drilled by the sphere. These cases can be definitively excluded for the experiments we are dealing with since no fracture, plasticity nor creep have been observed during these experiments. One concludes that for any elastic material, the elastic energy of the gel is proportional to δ 3 within the limit δ a. The proportionality constant is C 0 (see Eq. 18) depends on the whole elastic constitutive law of the elastic material.

F. Equilibrium condition
The energy is minimum in the equilibrium state. The δ-derivatives of E gr and E el are therefore equal giving the scaling law valid for δ a, for a given elastic material: Note that the beads has been assumed to be spherical although this is not a crucial point in this theory: if it is compact with all typical lengths of the same order a the previous scaling arguments apply.

III. DISCUSSION AND CONCLUSION
The theoretical analysis presented above provides an understanding sui generis of the experimentally observed scaling law.
The elastic material being chosen (so that the proportionality constant C 0 is fixed), this law predicts the scaling behaviour of the downshift of spheres with different radii and/or different densities. These predictions are in quantitative agreement with the observations: one observes experimentally a scaling law with an exponent close to 3/2 for the radius, and the prefactor is not only related to the linear elastic properties (µ), but also to non-linear features of the elastic material. These non-linear features being a priori distinct from one gel to another one, this explains why plotting δ/δ 0 as function of a/δ 0 for different gel compositions yields to different curves. The prefactor in the scaling law depends not only on the shear modulus of the solid, but also on the elastic behaviour at large deformations, in agreement with the experiments ( [15]). This provides a way to assess some characteristics of elastic materials under large strains, as strain stiffening properties.
The range of applications of the derived scaling goes far beyond the description of the elasto-buoyancy phenomenon. For instance, it can be directly checked that the recent indentation experiments performed on various types of compliant gels, by Fakhouri et al [15], are in quantitative agreement with our prediction, however their very naive model based on the neo-Hookean equation, misses the driving non-linearities.
In this paper, we have shown that the exponent 3/2 for the radius is independent of the strain-stress relation provided that the increase in the elastic energy density (W) with the strain is stiff enough. Generic behaviors for elastic materials undergoing large and complex deformations can therefore be identified, going beyond scaling arguments that is blind to the crucial effect of strain stiffening. This work sheds new light on the mechanics of extremely large elastic deformations that are crucial in various emerging techniques that even encompass an important procedure such as the computer-assisted surgery involving human organs [29][30][31], which are indeed non-linear materials and undergo large deformations. In all these important fields, an ideal neo-Hookean model is used that, according to our current work, is not valid sensu stricto. The danger is that one may miss the divergences of the solution that is inherent in the ideal neo-Hookean model because the full solution is not used. The work presented here is a panacea to the possible pitfalls that one may embark upon in studying extreme deformations of nonlinear materials, thereby also offering possible benchmarks for numerical simulations and even opportunities to formulate new simulations methods.