{'fraction_non_alphanumeric': 0.029407701265123255, 'fraction_numerical': 0.011999921053150966, 'mean_word_length': 4.50978686385385, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Post-buckling Dynamics of Spherical Shells

We explore the intrinsic dynamics of spherical shells immersed in a fluid in the vicinity of their buckled state, through experiments and 3D axisymmetric simulations. The results are supported by a theoretical model that accurately describes the buckled shell as a two-variable-only oscillator. We quantify the effective"softening"of shells above the buckling threshold, as observed in recent experiments on interactions between encapsulated microbubbles and acoustic waves. The main dissipation mechanism in the neighboring fluid is also evidenced.


Introduction
Buckling of elastic structures has recently emerged as a powerful mechanism to trigger fast motion at small scale. This includes fast reorientation of microswimmers [1,2], thrust generation in fluids [3,4] or, through solid friction [5], valves actuation for flow control [6,7] or fast actuation of optical lenses [8]. Smart design of materials in order to obtain the desired buckling behavior has now become an intense field of research [9][10][11][12]. From a modeling perspective, stable configurations of structures prone to buckling have been widely explored, as well as the stress or strain thresholds that have to be overcome to switch between states. The existence of (at least) two stable states separated by energy barriers allows for the design of robust devices that can maintain their state within a certain range of external perturbation without any external energy input. While the picture in terms of equilibrium states is now quite clear, at least for simple geometries (rods, half-spheres, closed spheres), full control of soft structures by external signals requires to know more about their dynamics. Recent papers have shed light on the complexity of the first-stage dynamics, close to the buckling threshold, where the response time of the material depends on classical dissipative mechanisms coupled to intrinsic slowing down observed in such a critical phenomenon [13][14][15]. The goal of the present paper is to explore the second-stage dynamics, around the buckled state, where the geometry is often more complex than that of the initial state.
Because of their ubiquitousness in nature but also of their simplicity in terms of fabrication and of modeling, spherical closed shells enclosing a compressible fluid are particularly in the spotlight [3,[16][17][18][19][20][21][22][23][24][25][26][27][28][29][30]. Existing studies are essentially focused on understanding the scenario of the buckling instability that occurs beyond a certain threshold of compression or deflation, and on characterizing the stability branches [22,24,25,29]. More recently, shells made of non isotropic material have also attracted some attention [31][32][33]. In terms of dynamics, the reaction of shells to a steep increase of pressure have been recently studied [15,34] while an experimental study has highlighted the role of dissipation within the shell while reaching the stable buckled state [29].
Hollow microshells have been used for decades as ultrasound contrast agents (UCAs) and their resonance frequencies in the spherical configuration have been widely studied [35]. Noteworthy, even in such a simple configuration, the existing models lack to describe accurately all experimental observations [36]. Depending on the applied acoustic field, UCAs may also buckle. In [37] the current state of a suspension of UCAs is controlled by a low frequency acoustic field while the propagation velocity of a high frequency acoustic signal is measured. The authors observe a decrease of this sound speed while the ambient pressure is increased above a certain threshold, in marked contrast with the standard results in a simple fluid. As in other preceding works [38,39], this is interpreted as a "softening" of the shell due to its buckling. This interpretation is consolidated by the existence of a hysteretic loop as the ambient pressure is varied, which is also a signature of the buckling-unbuckling transitions. A very different study reaches the same conclusion: in [40], primary Bjerkness forces on hollow micrometric shells are measured ; a strong rise of this force above a given amplitude of the applied acoustic field is interpreted again as a signature of the sudden "softening" of the shell. This interpretation is supported by independent measurements of the buckling pressure by AFM.
In all above mentioned studies, the data are not quantitatively fitted by a model. Indeed, to our knowledge, the sole model accounting for shell response in the buckled state is that of Marmottant et al. [41], which has been refined in [42]. These models assume that in the buckled configuration, the elastic response of the whole shell is simply that of the encapsulated gas while that due to shell material has disappeared, as if the shell was broken.
In the present work, we show that this rough approach is not valid. While our results confirm the effective softening due to buckling, we highlight a more complex interplay between gas response and shell material response, reaching the unexpected result that softening is even more pronounced than that obtained through neglecting shell material response.

Statement of the problem (a) Description of relevant parameters
We consider spherical elastic shells made of an isotropic, incompressible elastic material of Young's modulus E and initial thickness d. We have performed experiments on home-made shells of centimetric size (see Appendix A), which allowed to check the consistency of our numerical   Table 1. Some orders of magnitude for two typical objects of interest: ultrasound contrast agents (UCAs) [35,36] and polymeric macroscopic shells [3,5]. Initial inner pressure in UCAs is often assumed to be close to atmospheric pressure while it may be varied on purpose in macroscopic shells. The ranges of parameters covered by the present experiments and simulations are also indicated. See text for the definition of the dimensionless parameters (hatted symbols).
simulations where the different relevant parameters could be varied on a wider range. In these simulations, we consider isolated zero-thickness elastic shells of initial radius R 0 whose elastic constants (compression and curvature modulus) are functions of E and d (see Appendix B). These shells are immersed in a Newtonian fluid in which the Navier-Stokes equation is solved.
They are filled with a gas at pressure P which is assumed to be instantaneously set by the shell volume according to an adiabatic process, a reasonable hypothesis considering the high velocities encountered in this problem: P V κ = P 0 V κ 0 , where V is the shell volume, V 0 its initial value, P 0 the initial pressure, and κ is the polytropic coefficient. While the fabrication process of shells makes it difficult to obtain a measurable initial pressure other than the atmospheric one, the simulations have allowed to vary it so as to explore the relative contributions of gas compressibility and shell elasticity on the overall response. The thin shell limit that is considered here, though it may appear as a strong simplification, is indeed a relevant model even to describe thick shells, as evoked in [29] and confirmed in the following. The ambient pressure is initially equal to P 0 and is suddenly increased to a constant value P ext , which is large enough to trigger buckling.
The large range of parameters we consider here will allow us to establish post-buckling dynamics for a wide range of objects and scales, from the thin colloidal shells that are used, in particular, as UCAs [41,[43][44][45] or photoacoustics contrastagents [46], to macroscopic shells in elastomer which are among the favorite building blocks in soft robotics [3,5,9].
To describe this whole range, we consider the dimensionless problem obtained by considering the shell radius R 0 as the lengthscale (its midplane radius for a real shell), Ed/R 0 as the pressure and elastic modulus scale, and ρ f as the density scale. The time scale is then ρ f R 3 0 /(Ed). This scaling is that of the undamped period of an oscillating shell when the contribution of inner pressure is neglected.
The problem now depends on four parameters: the reduced thickness of the shelld = d/R 0 , the initial pressureP 0 = P 0 /(Ed/R 0 ) in the shell when it is in its stress-free spherical configuration, the applied pressureP ext and the dimensionless viscosity of the fluidη = η f / Edρ f R 0 that will characterize the damping in the system. This last number is the equivalent of an Ohnesorge number where surface tension has been replaced by the 2D elastic modulus Ed. We detail in Table  1 the typical values of the three parametersd,P 0 andη (characterizing the initial state) one may find when considering microshells and macroshells used within the current research context, as well as the range covered by our experiments and simulations.
Both in numerical simulations [22,24,26] or in experiments [29], it is now well established that at equilibrium, the pressure difference P ext − P quasi-plateaus to a constant ∆P pl as a function of equilibrium volumes Ve. Therefore, varying the fourth parameterP ext strictly amounts to varying the equilibrium pressurePe inside the shell, which we already set here by varying the initial pressure. In our experiments and simulations,P ext is typically chosen such that the pressure difference is right above the buckling threshold ∆P b = 4d/3 [22,26]. Increasing too much the external pressure, or starting with very low pressure inside the shell, leads to a full collapse of the shell, with the two opposite poles being in contact. This pertains to a new kind of physics taking into account solid friction and requires additional development in the numerical method. We will avoid these extreme situations in the present work. Oscillations of the width W and height H of the convex envelop of a shell with external radius R ext,0 = 25 mm, thickness d = 5 mm, Young modulus E = 0.5 MPa, initial internal pressure P 0 = 1 bar, external pressure Pext = 1.77 bar, immersed in glycerol (experiments in the reference configuration R). Black lines correspond to the fit by the sum of two damped sinusoidal functions, which are obviously necessary to recover the full curve for the width W . The two frequencies are of the order 30 and 90 Hz, and are much lower than that typical of spherical configuration, which is Hz (see e.g. [41]). Here R 0 is the midplane radius R ext,0 − d/2. Note that the sole gas contribution leads to a Minnaert frequency 1 Hz which is also larger than the observed frequencies. These observations call for a finer modeling of the oscillation mechanisms around the buckled state.

(b) A 2-frequency oscillator
In the experiments as in the simulations, we consider the dynamics of a shell which is brought to the onset of buckling by setting the external pressure P ext to a value determined by a preliminary study. The external pressure is then kept fixed at this value and the dynamics of deformation of the shell is recorded. The sudden collapse of a shell at the buckling transition is followed by oscillations around the new equilibrium position. We show in Fig. 1 the oscillations of two geometrical characteristics of a shell in configuration R, corresponding to one of our experiments. A noticeable effect is the apparition of a second frequency in the oscillation pattern ; moreover, in agreement with the widely admitted softening of buckled shells, these two frequencies are much lower than that typical of the spherical state.
The apparition of a second frequency invalidates the previous models where the oscillations come from the sole contribution of gas compressibility. Note that we cannot expect more complex behavior, like the apparition of a second frequency, to emerge from a bubble yet having a nonspherical shape: it has been shown recently [47], in agreement with [48], that the generalized Minnaert model, where the radius of the spherical bubble is replaced by the effective radius extracted from the volume, is robust against geometry changes.
In the experiments, buoyancy issues have required to restrict the displacement of the shell, which was attached to a fixed support through a suction pad located at the pole opposite to the buckling pole. In order to build up a general view of the oscillation mechanism devoid of any suspicion of strong influence of the boundary conditions, we turn to the numerical simulations 5 rspa.royalsocietypublishing.org Proc R Soc A 0000000  which allow to consider free shells and a larger range of parameters. We shall come back in section 5(d) to the impact of having a part of the shell that is fixed.

Post-buckling oscillations of a free shell
We illustrate in Fig. 2 the buckling process as obtained from the simulations. From the many available data characterizing the shell shape, we focus on its volume V and surface S and plot in Fig. 3 their oscillations towards their equilibrium values Ve and Se. The drop in volume is all the more pronounced as the initial pressure is low. By contrast, the shell surface increases along the buckling process, and reaches a value close to the initial value of the sphere, which hardly depends on the initial pressure. This validates the vision of the buckling process as a mechanism that makes the shell switch from a spherically compressed state to a more favorable state where the in-plane compression energy is released, the cost to pay being the curvature energy which is localized in the rim [19].
Right after the buckling transition, the signals are, as in the experiments, well fitted by the sum of two sinusoidal functions damped by the same decreasing exponential, as shown in Fig.  3(b). For each simulation, we then obtain three characteristic parameters of the shell dynamics: two frequencies ω − and ω + , with ω − < ω + , and a single relaxation time τ . We define ω ± as the intrinsic frequency, obtained from the measured frequency Ω ± through ω 2 ± = Ω 2 ± + 1/τ 2 . These  three parameters are similar for both V and S, as well as for other characteristic parameters such as the width or the height of the shell (as discussed with Fig. 9 further in the text). They all depend, a priori, on the 4 control parameters of the problem. In all the simulations we considered, the contribution of the highest frequency to the signal is between 50% to 80% smaller than that of the lowest frequency. We shall then focus first on this latter. In Fig. 4 we show how this frequency ω − varies with some control parameters, starting from the reference configuration R. In Fig. 4(a) the initial pressureP 0 is kept fixed and the external pressure is set right above the threshold for buckling, and the dependency with space parameters is explored. Decreasing the thickness d of a shell of given radius leads to the intuitive result that the frequency decreases, as a result of the decrease of the shell elasticity. For a given reduced thicknessd, Fig. 4(a) shows that the reduced frequency does not vary significantly for small enough values of the reduced viscosityη. When damping becomes more important, the frequency decreases, illustrating the threshold at which viscous stresses start to influence the dynamics of the shell beyond a mere linear damping. Less intuitive is the variation of the frequency with the equilibrium pressure inside the shell, shown in Fig. 4(b). Starting from shell in configuration R, we varied the initial pressureP 0 from 0.36 to 1.8. The resulting equilibrium pressurePe increases conjointly with this initial pressure. The result is that in the meantime, the frequency decreases. One can question this observation by considering, in a first approach, that the contribution of the gas to the frequency would scale like the Minnaert frequency ω M of a free spherical bubble: where we have simply replaced the usual radius of the original Minnaert expression by a generalized radius based on the volume, as in [47,48]. Since the volume at equilibrium Ve is an increasing function of internal pressure Pe, it is not clear a priori that the Minnaert frequency increases with pressure. In the spherical case, it can however be calculated that V 3/2 e does not increase as quickly as Pe. In our case of buckled configuration, data from the simulations show also that such a Minnaert contribution increases with the equilibrium (or the initial) pressure. Therefore the origin of the decrease of ω − with increasing internal pressure has to be found elsewhere.
In order to separate the direct effect of pressure on compressibility from its indirect effect through its influence on the shell volume, it is insightful to consider the following configuration: from the reference configuration R in its buckled equilibrium state, we have varied the inner and external pressure by the same amount δP , which does not modify the shape of the shell. Then, using this state as the reference state for inner pressure calculation (i.e.PV κ = (Pe + δP )V κ e ), we apply a small perturbation to the shell and measure the induced oscillations, thus characterizing the pressure contribution to the elastic response, at fixed geometry. As shown in Fig. 4(b), while this contribution is now increasing with the pressure, this increase of the frequency squared with the inner pressure is sublinear, in contrast with the spherical case. Note that the procedure amounts to varyP ext together withP 0 such that the final shape is the same in all simulations. We checked, for a zero offset δP , that these oscillations with small amplitudes are similar to that obtained right after the shell has buckled, which initially implies larger deformations. This can be seen by the proximity between the two data points forPe 1.5 on Fig. 4 In order to decipher these behaviors, we introduce a reduced model which is then used to fit the dataset obtained from the simulations.

Model
We describe the dynamics of the system through Lagrangian mechanics which is an energy approach that allows to introduce the generalized coordinates that are the more relevant to the problem (see e.g. [49] for such an approach in the case of spherical oscillations).
In this frame, the dynamics is given by: The Lagrangian function L is defined as L = T − U with T and U the kinetic and potential energies of the system, respectively. The Rayleigh dissipative function D allows to introduce non conservative forces in the formalism, providing they originate from viscous friction and that pressure drag is negligible. We shall see later that, due to the high buckling velocities, this assumption is fragile. The q i andq i are the generalized coordinates and velocities that characterize the state of the system. The observation that the dynamics is well described by two frequencies calls for a modeling of the system as a pair of coupled oscillators. We then choose to describe the shell through two parameters, that naturally appear as relevant parameters in this problem: while the volume V is intimately related to the gas behavior and is also the usual control or output parameter in shell buckling studies, the surface S appears in the stretching contribution of the elastic energy. While in the spherical case both parameters are interdependent, in the buckled state that introduces at least one more degree of freedom they can vary independently, as depicted in Fig. 5.
The modeling consists in writing how the elastic energy of the shell material and the potential energy of the inner gas depend on these two variables V and S, as well as the kinetic energy and the dissipation due to viscous friction.

(a) Potential energy
The potential energy can be written U = Ug + U ext + U S . The internal energy Ug of the gas reads, assuming an adiabatic behavior: The term U ext is the opposite of the work done by the external pressure on the shell, therefore Finally, U S is the elastic energy in the shell. For a given value of external pressure, equilibrium conditions ∂L ∂V = 0 and ∂L ∂S = 0 lead to the following set of equations that determine the equilibrium configuration (Ve, Se): As recalled before, at equilibrium the pressure difference P ext − P quasi-plateaus to a constant ∆P pl as a function of equilibrium volumes Ve. The plateau pressure ∆P pl depends on shell elastic properties and thickness over radius ratio. While there is to our knowledge no available theoretical expression for it, an ad-hoc expression was proposed in [29] based on numerical simulations.
Buckling is the mechanism by which the compressed spherical shell switches into an energetically more favorable configuration. Since bending energy is located in the rim, one can consider in a first approximation that at equilibrium, the surface in the buckled state is not compressed, as shown in Fig. 3(a) where the surface area at equilibrium is only 5% lower than the initial one. In the scheme of Fig. 5 (a), it corresponds to an inversion of the right spherical cap towards the left pole. This implies that we can consider in first approximation that Se = 4πR 2 0 , independently from Ve. As a consequence of this and of the existence of the plateau pressure This assumption has been validated afterwards by letting the above partial derivative be a free parameter in the fitting procedure of the data, that eventually lead to very small values compared to the other elastic constants of the problem. Around an equilibrium configuration, the elastic energy of the shell finally reads . These parameters depend, a priori, on the equilibrium state.
As discussed above, we expect the variations of elastic energy with S to be that of the stretching energy in the spherical case, since the curvature energy is localized in the rim. The elastic energy for a sphere of radius R can be directly obtained from the general expressions that are recalled in  is a convenient simplified model that has proven its efficiency to describe buckling dynamics of shallow arch [52] or elastic hemisphere (the "jumping popper" toy) [14], in particular when the system is subjected to oscillating external forces. In this paper, we do not go further in the analogy which will be developed in subsequent papers, but we simply make the remark that in the 3D shape obtained by axial rotation along the horizontal axis of the closed Mises truss depicted in (b) and (c), with the two lateral masses always located at fixed distance 2R 0 , V variations can be made while keeping S constant, as in (b), while in (c) surface variations are obtained while keeping V constant. This is obtained by a similar horizontal motion of both poles, keeping their interdistance constant. This results in an hydrodynamic dipole p at long range. Its intensity is related to the volume changes on each side of the midsurface S = S/2 indicated in grey. The assumption of fixed distance 2R 0 between the two flanks is supported by the numerical simulations where we saw that most of the deformation occurs along the axis of symmetry (Fig. 2(h)).
the Appendix B (Eqs. A 6 and A 7): We have introduced the area bulk modulus K A = dE and the bending stiffness K B = d 3 E/9. From d 2 U stretch /dS 2 | Se = Ed/Se in the spherical case, we therefore postulate that A SS = Ed/Se. Similar reasoning is not possible for A V S , as V and S are interdependent in the spherical case. We simply assume that the scaling A V S = Ed/Ve k V S , where k V S is a dimensionless constant, will capture the dependency with space variables. Since there are two relevant space variables that are R 0 (or Se) and Ve, other scalings could be possible. Complementary tests with the general scaling )k V S while proceeding with the fitting procedure of the data according to the model have indeed shown that β 1 is the best choice.

(b) Kinetic energy
We now turn to the determination of the kinetic energy T . The determination of the fluid flow around the oscillating shell is a complex problem, for which we only consider the long distance solution, obtained by adding the contributions of the constant surface and of the constant volume motions of the shell.
Far from the shell (i.e. neglecting the shape details), shell volume variations lead to a flow induced by a point source, imposing a flow rate q =V . By incompressiblity of the fluid, one obtains a first (radial) contribution v 1 to the fluid velocity:  Determining the flow due to surface variations is more complex. We wish to establish its dependency with the two dynamical variables V and S. To that end, the geometry of the shell As sketched in Fig. 5(c), this simplification of the shell geometry leads to assimilate the deformation at constant volume to the synchronized motion of the two poles, thus considered as a pair of source/sink of flow rates q and −q . This gives rise to an hydrodynamic dipole producing the flow field where θ is the angle between up that points in the pole to pole direction and the radial unit vector ur, and u θ is the associated unit vector (see Fig. 5(c)). The dipolar moment is p = q , where is the pole to pole distance. The dipolar approximation is well supported by the PIV observations made in [3], though it has been shown there that more complex patterns can take place due to shear waves near the flanks. The flow rate q is related to the deformation of the shell. As we have symmetrized the problem in a first approximation, by considering that the source and the sink have the same absolute intensity, we relate the flow rate q to the motion of the midsurface (of area S = S/2) going through the midpoint between the two poles, depicted in grey on Fig.  5(c). As the poles move with velocity vp, the rate of change of volume on each side of the midsurface is q /2 = πR 2 0 vp/3. A bit of geometry gives the relation between the pole velocity vp and the variationṠ of the midsurface, eventually leading to Finally, as V = πR 2 0 /3, p can be written under the form Eventually, the kinetic energy is obtained by integrating 1 2 ρ f (v 1 + v 2 ) 2 on the whole fluid volume Ω. By lack of knowledge on the effective shape, we integrate from R 0 to ∞ and compensate our successive approximations by introducing numerical prefactors α 1 and α 2 to the monopolar and dipolar terms, that we expect to be of order unity. They constitute the cost to pay in this simplified -then tractable -model.

(c) Dissipative function
The dissipative function D reads D = η f Ω (u ij ) 2 dV , where u ij is the symmetric strainrate tensor in the fluid. We implicitly assume here that dissipation originates from viscous friction and that pressure drag is negligible ; this is the framework in which the Rayleigh dissipative function has been introduced and validated so far. From the above expression for the fluid velocity, D can be straightforwardly calculated and we find by integrating from R 0 to ∞ where the β i have been introduced, as α i for the kinetic energy, to account for the lack of precision in the integration domain.
Together with the equilibrium conditions that set Pe and Ve, these equations allow to determine the motion of the shell as a function of the initial pressure P 0 inside the shell, the external pressure P ext , the initial shell radius R 0 , its thickness d, Young's modulus E and the fluid density and viscosity ρ f and η f . The impact of the geometry differs for each term considered in the equations, where R 0 and Ve appear with different powers. Switching to the dimensionless parameters, Eqs. 4.13 and 4.14 become: In most cases of interest (see Table 1),η is small. This means that the fluid does not influence the shell dynamics, which would simplify the design and optimization of such a system, as well as its modeling: we can considerη as a small parameter and look for solutions of Eqs Taking the leading order inδ andη (which is the 0th order) of the real part of the above equation we obtain two eigenpulsationŝ (4.20) This constitutes the central result of this study. The leading order of the imaginary part of Eq. 4.19 gives the two relaxation timesτ ± =δ −1 ± associated with the two above pulsations. A priori, these two times can be different. Nevertheless, we directly make use of the observation that all obtained oscillation signals in the simulations are very well fitted by the sum of two sinusoidal functions damped by the same exponential function. Consequently, there is only one damping time. It can be shown thatτ + −τ − is proportional to 5α 2 β 1 − 18α 1 β 2 , which is therefore taken to be 0. Injecting the obtained value of β 2 in the expression forτ ≡τ ± , we obtain

Discussion (a) Oscillation frequencies
We discuss the agreement between the model and the simulations. According to these, for the maximal value exploredη = 1.3, the dynamics is fully damped. Belowη = 0.13, few oscillations are seen and we are only able to determine the frequencyω − that contributes the most to the signal. Belowη = 0.013 both frequencies can be accurately determined. This range corresponds also to the small damping limit that lead to Eq. 4.20. This equation suggests thatV 2 eω 2 ± plotted as a function ofVePe follows a mastercurve. In Fig. 6, we follow this idea and find that indeed the data collapse reasonably on a single curve considering the large range of control parameters considered. The set of data forV 2 eω 2 − andV 2 eω 2 + are well fitted by the proposed equation (Eq. 4.20), the free parameters being α 1 , α 2 , and k V S .
Fitting the two frequencies together leads to α 1 = 3.6 ± 3.5%, α 2 = 1.1 ± 4.1% , and k V S = 0.23 ± 6.6%. These are acceptable ranges: α i larger than 1 was expected since it is intended to be a correction of the fact we haven not included the fluid inside the sphere of radius R 0 in the calculation of the kinetic energy.
Interestingly, ω 2 − vanishes for a finite pressure, which sets a stability limit of the buckled state, which is given byV This provides the minimal pressure needed to maintain the buckled shape. The physical meaning of the two frequencies is particularly clear in the largePe limit. From Eq. 4.17, we see that ω + and ω − can then be interpreted as the contributions of the (now decoupled) volume and surface oscillations, respectively. From Eq. 4.20, we obtain   [41] where it is assumed that the pulsation is the Minnaert pulsation for the shell of radius R 0 at ambient pressure: Though this expression ressembles a Minnaert contribution, as in the model of Marmottant et al. [41], it indeed differs from it by a factor R 0 /(α 1 Re) and will be smaller in most cases, as illustrated in the next section.
The smallest pulsation ω − becomes, in the largePe limit: Contrary to the spherical case, where the remaining contribution at high pressure is that of the gas pressure, we end up here with the sole contribution of the shell material compressibility. The high pressure makes the whole shell incompressible, but the buckled geometry makes it possible for the surface to evolve independently, with its own associated stiffness A SS . The equivalent picture is that of two springs in series, one of them becoming infinitely stiff: the dynamics is then given by the smooth one. This contrasts with the case of of the spherical geometry, where the springs are in parallel (see Eq. A 13).
The above considerations help explaining the behavior shown in Fig. 4(b): asP 0 increases,Ve increases and Eq. 5.3 shows that the pulsation is eventually a decreasing function of the initial -or of the equilibrium -pressure. For fixed volumeVe and varying inner pressure, Eq. 5.3 tells us that the frequency converges towards a finite limit, thus explaining why it does not increase linearly with the pressure (as in the spherical case), as shown in Fig. 4(b).

(b) "Softening" of buckled shells
Our results demonstrate that the oscillation frequencies of a shell drop after it has buckled. This so-called "softening" is illustrated in Fig. 7 where we compare the obtained frequencies with that of the spherical case before buckling and that proposed in Ref. [41] for the buckled case. To do so, we consider increasing values of the external pressure P ext starting at P 0 such that the shell is first in its spherical configuration until the pressure difference reaches the buckling pressure. In the buckled state, we consider either the Marmottant et al. model [41] or ours. The equilibrium configurations are calculated for each value of P ext , with the internal pressure given by the adiabatic condition (Eq. 4.16) and the resulting pressure difference being equilibrated by where Cx is the drag coefficient [53] (that depends on Re) and Re the Reynolds number.
The full line indicates this quantity, which matches perfectly the measured damping time for the oscillations without any adjusted parameter.
the elastic restoring force, which is given by Eq. 4.15 in the buckled case and by −dU S,sph /dV in the spherical case, where U S,sph is given by Eq. 4.7.
We show in particular that the "softening" is characterized by a drop of the oscillation frequencies by factors 5 (for ω − ) and more than 2 (for ω + ) for parameters that are typical of the usual commercial UCAs, and is much more pronounced than that predicted in previous models, over a large range of parameters.

(c) Damping
The scaling of the damping time τ = 2 α1 β1 ρ f η f R 2 0 predicted by Eq. 4.21 is the only one that can be built when considering only viscous effects in a fluid of density ρ f , viscosity η f in the vicinity of an object of size R 0 . Indeed, the assumption of viscous friction only was needed to calculate a dissipative function. It turns out that, due to the high Reynolds numbers that characterize the flow, this scaling does not hold in most of the considered range of parameters. We consider the Reynolds number based on the maximal flow velocity v pole , that is located near the buckling pole, which is then defined as Re = v pole R 0 ρ f /η f . As shown in Fig. 8, the dimensionless relaxation time η f τ /(ρ f R 2 0 ) indeed depends on the Reynolds number, while it would be constant if the dissipation was only of viscous origin (see Eq. 4.21). The dependency with the Reynolds number is indeed pretty close (with no adjusted parameter) to that of the relaxation time towards constant velocity for a settling shell of radius R 0 , for which the drag coefficient Cx is known [53]. This leads to the following general expression for the relaxation time: The scaling of Eq. 4.21, that was valid only in the low Reynolds limit, is preserved, but a dependency of the prefactor with the Reynolds number must be introduced, in order to account for the hydrodynamics regime. In particular, in the heart of our space parameter of interest, the dissipation regime is thus that of the Allen regime [54], where a combination of viscous and pressure effects leads to a drag coefficient that roughly scales like R −0.6 e , in the 1 < Re < 10 3 range.
As for the purely viscous case, considering the data forη = 0.13, where Re < 0.1, Eq. 4.21 gives α 1 /β 1 = 0.17, quite close to the case of spherical oscillations where it would be equal to 0.125 according to [49].
We eventually make the remark that the model predicts that, for not too highη -which will be the case in most situations -the frequencies scale with 1/R 0 , for a given reduced thickness d/R 0 . This implies that the typical maximal velocity at the pole v pole , that scales like ωR 0 , is scale independent. In particular, on a large range of inner pressures, we can use the large pressure limit of Eq. 5.3 as an estimate for the leading frequency ω − , leading to With typical values E d R0 = 1 bar, ρ f = 1000 Kg/m 3 andVe = 1/2 × 4π/3, we find vp ∼ 15 m·s −1 . This implies that, even at small scales, the Reynolds number will be quite high, in the range 1 − 10 3 discussed above, that was obtained in the simulations.

(d) Fixed shells: the effect of boundary conditions
In practical applications, a shell may be linked to other components. For our experiments on six different shells (see Appendix A), buoyancy issues have required to restrict the displacement of the shell. For the three shells of external radii 7.5 mm, 25 mm and 75 mm, but identicald we checked the scalings ω ∝ 1/R. We have also simulated the same six shells in the same attachment conditions, assuming the part of the shell connected to the suction pad is fixed. In spite of this simplification, we find good agreement between the simulations and the experiments within 15% maximal variation for the measured frequencies, with no adjustable parameters.
In addition, these simulations show that the oscillation frequency is smaller by 25 to 50 % than that of the free shell case. This is in qualitative agreement with the decrease of oscillation frequency of an initially free beam when it is attached at one end [19]. Yet, this fact is not that intuitive since the attachment may lead to less fluid motion in the vicinity of the attachment patch, therefore less accelerated mass. Noteworthy, the additional local constraint also modifies the relative weights of the two sinusoidal responses, depending on the considered quantity ( Fig.  9). While both frequencies are well present whatever the boundary condition when considering the volume oscillations, the picture is different when considering the oscillation of the height H of the shell or of its width W , as defined in Fig. 1. In particular, for the height, a single frequency becomes predominant in the attached case, as seen in Fig. 1 for the experiments. As a signature of robustness of the modeling, the same pair of frequencies is sufficient to describe the oscillations whatever the considered quantity. A more systematic study of the dependency of each frequency on the boundary condition, and of its relative weight is the full signal, would require to define and parametrize the attachment area precisely. We simply conclude here that the two-oscillator modeling is a robust description that remains valid for shells with fixed parts.

Conclusion
In spite of the potentially huge number of deformation modes of a buckled shell, its dynamics can be described by a pair of coupled oscillators. The two resulting modes decouple in the high inner pressure limit and can be identified as the surface oscillation mode and the volume oscillation mode. The former corresponds to the lowest frequency and contributes more to the overall signal than the high frequency mode. Both frequencies are much lower than that in the spherical case. This confirms and quantifies the apparent softening of buckled shells observed in the literature.  By contrast with previous models [41,42], we show that this softening is mainly due to the absence of contribution of the inner gas.
We have also shown that over a large range of parameters, the dissipation mechanism is similar to that of a sphere in translation in a fluid, whose drag coefficient is Reynolds-dependent.
The established model can be used to anticipate the behavior of shells under more complex actuations than a fixed applied pressure. Using the model of Marmottant [41], it has recently been shown that buckling, by introducing an additional non-linearity, may trigger complex dynamics characterized by subharmonic behavior or even chaos [55]. More work is also needed to identify what controls the relative weight of each mode, and how this weight will be influenced by the boundary conditions, like local bonds to another component. We plan to dig in that direction in a near future. The model could also be further enriched by incorporating the contribution of the bending energy, which may help improving the collapse of data of Fig. 6 since it would introduce another dependency with the shell thickness. This contribution is already included in the description of equilibrium through the expression of the plateau pressure (Eq. 4.15) but it is still poorly understood. Another open question is that of the viscous dissipation within the shell material.

A. Experimental set-up
As in Ref. [3], the spherical shells were realized by molding two half hollow spheres made of commercial elastomer resin "Dragon Skin ©  glued together using the same material. Flat disks of diameter of order 20% the shell diameter and thickness around 4% the shell thickness were added at the pole of one of the two external moulds so as to create a weak point on one hemisphere, where buckling will systematically occur.
In order to explore the effect of size and of thickness over radius ratio, we have considered 6 different shells whose external radii R ext,0 lie between 7.5 and 75 mm while their dimensionless thicknesses d/(R ext,0 − d/2) are between 0.08 and 0.3. We cover one order of magnitude in size and the relative thickness was varied by a factor almost 4. Reaching broader ranges poses technical issues in the manufacturing process.
The shell under study was attached at the level of its pole opposite to the weak one to a fixed support and immersed in a 60 × 60 × 60 cm tank in anodized aluminum with polycarbonate polymer windows that could bear pressure differences of + 2 bar. The tank was filled with glycerol. Pressure variations were obtained with a pressure controller (OB1 by Elveflow) connected to the thin layer of air that was left in the tank on top of the liquid. Shell deformations were recorded using a fast camera. Pictures of an immersed grid of known characteristics allows to characterize the deformation of pictures due to window deformation by the overpressure. As in Fig. 1, the obtained images of the convex envelops of the shells are well contrasted and their contours could be directly extracted and analyzed using home-made routines written in Python. The problem is numerically solved with finite elements making usage of the finite element toolbox AMDiS [56]. An axisymmetric arbitrary Lagrangian-Eulerian (ALE) method is employed according to [57], where the membrane points move with the fluid velocity v. This movement is harmonically extended to the internal part of the grid using one of the mesh smoothing approaches presented in [57].

B. Numerical method
The shell is modeled as an infinitely thin axisymmetric elastic surface, immersed in a Newtonian viscous fluid of viscosity η f and density ρ f and filled with air. The hydrodynamics in the surrounding fluid Ω are governed by the Navier-Stokes equations. The pressure due to the internal air is implemented as boundary conditions on the shell surface Γ . Assuming axisymmetric flow conditions reduces the problem to a 2D meridian x − r− domain, describing half of the 3D domain's cross section. The simulation domain is shown in Fig. 10. The problem in axisymmetric formulation reads [58,59] in Ω (A 1) where v = (vx, vr) T and P f are the fluid velocity and fluid pressure, respectively. The axisymmetric terms read The elastic shell reacts to in-plane (stretching) and out-of-plane (bending) deformations, quantified by the stretching energy E stretch and the bending energy E bend . The stretching energy in an axisymmetric setting can be described in terms of the two principal stretches λ 1 and λ 2 , which provide information about relative changes of surface lengths in lateral and circumferential direction, respectively. It is with the arc length parameter s of the interface curve Γ and the distance of the shell r S to the symmetry axis. The subscript 0 refers to the quantities in initial state. The shell energies then read with area bulk modulus K A , area shear modulus K S , bending stiffness K B , mean curvature H and spontaneous curvature H 0 , which is the mean curvature in the initial state. These three quantities can all be calculated using Young's modulus E, the Poisson ratio ν (which will be set to 0.5), and the (initial) thickness d of the shell: .
The force exerted by the shell is composed of the sum of the first variations of these energies and of the force due to the pressure of the internal air. Note that our thin shell model is a simplification that is made for technical purpose, but the elastic model describes a shell of any thickness. Since the strains are small compared to the deformations of a thin shell, the in-plane force is linearized in λ i (see [60], Appendix). In-plane and out-of-plane forces read [57] where Kg is the Gaussian curvature. For the contribution of internal pressure, we assume homogeneous air pressure P inside and use adiabatic gas theory to relate this pressure to the inner shell volume V by P = P 0 (V 0 /V ) κ , where P 0 and V 0 denote the respective initial values, and κ = 1.4 by assumption of an adiabatic process. Accordingly, the stress exerted by the air is reduced to P n = P 0 (V 0 /V ) κ n whereupon the boundary conditions become  The numerical scheme in time step m + 1 can be summarized as follows. An implicit Euler method is employed for the Navier-Stokes equations. In the following, time steps indices are denoted by superscripts. The scheme reads: (i) Calculate λ 1 , λ 2 , H, Kg and n from the position of boundary grid points along Γ (see [57] for details). This allows computation of δE stretch (iv) Calculate the grid velocity v grid by using one of the mesh smoothing algorithms described in [57] to harmonically extend the interface movement into Ω and move all internal grid points of Ω with velocity v grid .
A function P ext (t) imposes variable pressure at the boundary ∂Ω\Γ of the computational domain, excluding the shell (see Eq. A 12). The values for P ext start for t = 0 at P 0 . Afterwards, P ext is increased slowly until the maximum desired value Pmax is reached. The slow increase of pressure circumvents the occurrence of shell oscillations in the spherical compression stage before buckling. In order to analyze the post buckling oscillations, the maximum pressure value is kept constant afterwards.
For the simulation results shown in Fig. 4, we first compute the equilibrium state after buckling as explained above. Then, we impose an offset δP to both, the internal pressure P and the external pressure P ext . The updated internal pressure acts as the new "initial" pressure P 0 for the rest of the simulation. Subsequently, a perturbation to the shell is applied by adding another prescribed pressure P perturb to the new P ext . If the offset dP is zero, we could confirm that these oscillations, whose initial amplitude is smaller, are similar to those obtained right after the shell has buckled, which initially implies larger deformation.
A weak spot is imposed into the formula when calculating the elastic moduli K A , K S , and K B by decreasing the thickness d locally around the center of the weak spot, which is the membrane point touching the symmetry axis (r = 0) on the right. This weak spot provides the perturbation needed to trigger buckling, otherwise our code is stable enough to let the shell stay in its spherical configuration, though unstable. The oscillation frequencies are the main quantitative output of the simulations. They directly depend on the total accelerated mass and because of that, the size of the simulation box may bias the results. To estimate this, we consider the spherical oscillations of the shell under low pressure difference, such that it has not buckled.
To that purpose, we considered simulations of the shell in configuration R but considered different initial pressure inside and outside the shellP 0 going from 0.18 to 1.8. Then we applied suddenly an external pressureP ext =P 0 + δP to the shell. For not too high values of δP , the shell shrinks isotropically with damped oscillations towards its new configuration, characterized by an equilibrium volumeVe.
The radial oscillations of a coated shell have been documented by several authors, for 2D shells (thin shells) [41,61] and for real shell of finite thickness [36,49,62,63]. For 2D shells, the expressions for the oscillation frequency that are proposed in the literature are generally obtained considering oscillations around the stress-free configuration of volume V 0 . Here, we are interested in oscillations around any state, which results in a modified expression for the pulsation, which we derive hereafter. In addition, models for thin shells generally consider only the stretching energy, which scales with d, and discard the curvature energy, which scales with d 3 . Here, in Eq. 4.7, we included both contributions in order to account for the oscillations of the mid-plane of shells of any thickness. This elastic energy in the spherical configuration given in Eq. 4.7 can be rewritten in terms of V and V 0 , and we can follow the Lagrangian formalism that has lead to Eq. 4.13. The pressure term is identical and the acceleration term can be calculated exactly by integrating the kinetic energy arising from the monopolar term from the equilibrium radius Re to ∞.
The pulsation in the absence of damping then reads In the absence of bending energy, and for oscillations around the stress-free state (Ve = V 0 ), we recover the usual expression [41,61] ρ f R 2 0 ω 2 sph = 3κP 0 + 4Ka/R 0 . In Fig. 11, we plot the obtained pulsations in the simulation and compare them with the expected values obtained from Eq. A 13. It is shown that the dependency with the shell parameters is recovered but the theoretical expression for ω sph has to be multiplied by a prefactor 1.07 to obtain an agreement between numerical data and theory.
This means that while the physics is well described by the simulations, they lead to an overestimation of the pulsation by 7%. This overestimation factor does not depend on the chosen configuration so we associate it with the finite size of the simulation box Ω, which we checked by decreasing its size. The considered size in the paper is the result of a compromise between accuracy of the simulations and computation time.
Data Accessibility. All data available from the authors upon reasonable request.
Authors' Contributions. MM participated in the design of the study, developed the simulation code, carried out the numerical simulations, and critically revised the manuscript; AD performed the experiments and the image analysis and critically revised the manuscript; CQ participated in the design of the study and critically revised the manuscript; SA developed the numerical model, participated in the design of the study and critically revised the manuscript. GC designed and coordinated the study, performed the data analysis, developed the theoretical model and drafted the manuscript. All authors gave final approval for publication and agree to be held accountable for the work performed therein.
Competing Interests. We declare we have no competing interests.