INTRODUCTION

Various modalities have been used for diagnostic imaging, such as clinical radiography, computed tomography (CT), magnetic resonance imaging (MRI) and ultrasounds (US). Being safe, noninvasive and relatively cheap, US imaging techniques have been much improved by the introduction of ultrasound micron size contrast agents (UCAs). The need of image enhancers was essential because first, human blood within an organ has poor scattering properties and low signal amplitude relative to human tissues, that generate strong echoes, and second, the old-school Doppler technique [START_REF] Campbell | New doppler technique for assessing uteroplacental blood flow[END_REF] could not satisfy anymore the demands, especially in some more complex and confined geometries.

The presence of air bubbles inside an injected hand-agitated saline solution during an echocardiography was the first ultrasonic image enhancement technique that was proposed [START_REF] Gramiak | Echocardiography of the aortic root[END_REF]). These bubbles are known as the first generation of UCAs.

As they dissolve rapidly in the liquid, a second generation of UCAs was developed, that are made of air bubbles encapsulated by a thin shell: galactose as in Echovist R (1991) or albumin (a human protein) as in Albunex R (1995) or galactose and palmitic acid as in Levovist R (1995). Finally, the third generation of UCAs includes microbubbles with higher life-time, air being simply replaced by a gas with higher molecular weight, responsible of decreased solubility: SF 6 as in Sonovue R (2001), C 3 F 8 as in Definity R (2001) or C 4 F 10 as in Sonazoid R (2007). All of these gases are encapsulated by phospholipids. The resulting shells are known as soft-shell UCAs, while the ones made with polymers are known as hard-shell UCAs.

UCAs react to high amplitude pulses (1 MPa) of short duration (of the order of the µs). Even in presence of such agents, axial and lateral resolutions of ultrasonic devices used in clinical applications are limited by diffraction, such that the resolution is fixed by the typical wavelength, which lies between 100 µm and 1 mm in practice. This limit has been strikingly overcome recently: by analysing the transient signal re-emitted by UCAs, in-vivo resolution has been decreased to about 10 µm, at a detection frequency high enough to also allow velocity measurement in blood flow by image correlation [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF]. Such a super-resolved technique can also be implemented through a photoacoustic device, where the UCAs are excited by light rather than by sound [START_REF] Vilov | Overcoming the acoustic diffraction limit in photoacoustic imaging by the localization of flowing absorbers[END_REF]. The scattered ultrasound signal of UCAs has also been recently used to discriminate between two networks topologies, with application to cancerous tumor detection [START_REF] Mohanty | Ultrasound multiple scattering with microbubbles can differentiate between tumor and healthy tissue in vivo[END_REF].

These recent advances are based on a complex interplay between hardware development and post-processing to extract the relevant information from the acquired signal. The response of a shelled bubble is strongly dependent on its size and on the shell material properties. While commercial UCAs are quite polydisperse in size, narrowing the size distribution of UCAs appears then as a way to match better the relatively narrow frequency bandwidth of ultrasonic devices with that of the UCAs, thus leading to better sensitivity of the whole detection process. Recent works go in that direction, that make use of shell material of various types, such as polymers [START_REF] Liu | Fabrication of uniform sized polylactone microcapsules by premix membrane emulsification for ultrasound imaging[END_REF][START_REF] Song | Controllable formation of monodisperse polymer microbubbles as ultrasound contrast agents[END_REF], phospholipids (Gong et al., 2014; Helfield et al., 2014; Lum et al., 2016; Parrales et al., 2014; Segers et al., 2016 , 2020; van Rooij et al., 2015) -sometimes forming more than 2 layers [START_REF] Shafi | Probing phospholipid microbubbles by atomic force microscopy to quantify bubble mechanics and nanostructural shell properties[END_REF], silica [START_REF] Hu | Biocompatiable hollow silica microspheres as novel ultrasound contrast agents for in vivo imaging[END_REF] or proteins [START_REF] Wang | The fabrication of protein microbubbles with diverse gas core and the novel exploration on the role of interface introduction in protein crystallization[END_REF]. This calls for models of bubble oscillations that are able to describe a wide variety of shell materials. While previous ones are focused on incompressible and isotropic material, we present here a model that includes compressibility and the possibility for spherical UCAs to present different material properties in the radial and orthoradial directions ("transverse isotropic" material), a feature that would naturally occur for layered shells like lipidic ones.

Such anisotropy has been shown to greatly influence the buckling process of shells [START_REF] Munglani | Collapse of orthotropic spherical shells[END_REF][START_REF] Pitois | On the collapse pressure of armored bubbles and drops[END_REF][START_REF] Quemeneur | Gel phase vesicles buckle into specific shapes[END_REF].

II. PREVIOUS MODELS

Since the early work of [START_REF] Besant | A Treatise on Hydrostatics and Hydrodynamics[END_REF] who was concerned by the time needed to fill up the empty space of a collapsed bubble and the pressure generated at any point in an incompressible liquid, forced vibrations of bubbles have attracted attention for decades.

Giving a simpler derivation of Besant's results, Lord Rayleigh (1917) generalized the case to a cavity with nonzero pressure, i.e. to a gas-filled bubble. The surface tension and the viscosity of the surrounding fluid were taken into account (see [START_REF] Plesset | Bubble dynamics and cavitation[END_REF] for a review), leading to the famous Rayleigh-Plesset equation.

In order to take into account the shell encapsulating the microbubble, a semi-empirical model was developed [START_REF] De | Higher harmonics of vibrating gas-filled microspheres. part one: simulations[END_REF][START_REF] De Jong | Absorption and scatter of encapsulated gas filled microspheres: theoretical considerations and some measurements[END_REF] by way of the introduction of two ad-hoc quantities S p and S f that account for the effective elastic and dissipative properties of the interface. Assuming a zero-thickness shell, which is motivated by the proximity between the shell thickness and the molecular scale, other models have introduced rheological constants that are explicitly related to the expected properties of the shell material.

The first approach by [START_REF] Chatterjee | A newtonian rheological model for the interface of microbubble contrast agents[END_REF] was followed by the models of [START_REF] Sarkar | Characterization of ultrasound contrast microbubbles using in vitro experiments and viscous and viscoelastic interface models for encapsulation[END_REF] and [START_REF] Marmottant | A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture[END_REF]. In the latter, a non-linear model is proposed, presenting the elasticity of the shell as an effective surface tension. Its linearized form is equivalent to the de Jong model.

In [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF], a finite thickness shell was considered. It was assumed to be made of a homogeneous, incompressible and isotropic material, that was described by a Kelvin-Voigt model. This model was linearized relatively to the thickness to radius ratio in [START_REF] Hoff | Oscillations of polymeric microbubbles: Effect of the encapsulating shell[END_REF], giving rise to the Church-Hoff model. In [START_REF] Morgan | Experimental and theoretical evaluation of microbubble behavior: Effect of transmitted phase and bubble size[END_REF], thin-shell UCAs were described by a constant thickness model using bulk elasticity and viscosity.

Following [START_REF] Marmottant | A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture[END_REF], other nonlinear models have been proposed, with a more complex rheological behavior like strain softening and strain hardening [START_REF] Paul | Material characterization of the encapsulation of an ultrasound contrast microbubble and its subharmonic response: Strain-softening interfacial elasticity model[END_REF][START_REF] Tsiglifis | Nonlinear radial oscillations of encapsulated microbubbles subject to ultrasound: The effect of membrane constitutive law[END_REF], or shear thinning [START_REF] Doinikov | Modeling of nonlinear viscous stress in encapsulating shells of lipid-coated contrast agent microbubbles[END_REF][START_REF] Li | Modeling complicated rheological behaviors in encapsulating shells of lipid-coated microbubbles accounting for nonlinear changes of both shell viscosity and elasticity[END_REF].

III. CONFRONTATION WITH EXPERIMENTS

Vibration experiments on UCAs should a priori allow to determine the rheological constants of the material, through the chosen model among those cited above, as long as they are not too many. The final goal is usually to choose the best fitting couple of one elastic and one viscous parameter to describe the observed damped signal. This couple is unique for the model selected, for instance, (S p , S f ) in de [START_REF] De Jong | Absorption and scatter of encapsulated gas filled microspheres: theoretical considerations and some measurements[END_REF] model, (G S , µ S ) in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] model and so on. Note that using finite thickness shell models requires to make assumptions, or additional measurements, to determine the value of the shell thick-ness. In all models, additional assumptions are generally made in order not to consider the inner gas pressure as an unknown to be determined.

Several techniques can be used to determine the shell oscillations. In [START_REF] Gorce | Influence of bubble size distribution on the echogenicity of ultrasound contrast agents: A study of sonovue[END_REF], a batch of encapsulated microbubbles are insonated at frequencies up to 30 MHz, and the viscoelastic parameters are deduced by measuring the attenuation expression. The spectroscopy approach relies on using a high speed camera to directly measure the radial displacement of the UCAs, which is fitted with the theoretical one (van der [START_REF] Van Der Meer | Microbubble spectroscopy of ultrasound contrast agents[END_REF]. Light scattering methods were also developed [START_REF] Li | Modeling complicated rheological behaviors in encapsulating shells of lipid-coated microbubbles accounting for nonlinear changes of both shell viscosity and elasticity[END_REF][START_REF] Tu | Estimating the shell parameters of sonovue R microbubbles using light scattering[END_REF], where the scattering cross section is related to the resonance frequency containing the viscoelastic properties using the Mie scattering theory. A photoacoustic measurement technique was developed in [START_REF] Lum | Single microbubble measurements of lipid monolayer viscoelastic properties for small-amplitude oscillations[END_REF]. Readers can refer to [START_REF] Helfield | A review of phospholipid encapsulated ultrasound contrast agent microbubble physics[END_REF][START_REF] Van Rooij | Non-linear response and viscoelastic properties of lipid-coated microbubbles: DSPC versus DPPC[END_REF] for recent reviews on linear models theory and experimental measuring methods.

All the existing linear models are virtually the same, with 2D moduli that can be expressed explicitly in terms of 3D moduli and thickness. Tables I-IV summarize some experimental estimations of shell properties using the De Jong, Marmottant, Church-Hoff and Sarkar models respectively for different UCAs and using various techniques.

Such experiments may also be used to validate the model that is used to describe the results. This requires to determine by another means the rheological properties of the shell material. Such validations are scarce in the literature, and yield only accurate order of magnitude so far. The atomic force microscopy (AFM) is a direct approach used to estimate UCAs' properties. However, depending on the model used to extract elastic constants from the force-displacement curve of an AFM, very different values can be found (Abou-Saleh et al., 2013;[START_REF] Buchner Santos | Nanomechanical properties of phospholipid microbubbles[END_REF][START_REF] Lytra | Modeling atomic force microscopy and shell mechanical properties estimation of coated microbubbles[END_REF][START_REF] Shafi | Probing phospholipid microbubbles by atomic force microscopy to quantify bubble mechanics and nanostructural shell properties[END_REF]. This makes the validation of spherical oscillation models a tricky task so far. As an example, in Buchner [START_REF] Buchner Santos | Nanomechanical properties of phospholipid microbubbles[END_REF] and [START_REF] Lytra | Modeling atomic force microscopy and shell mechanical properties estimation of coated microbubbles[END_REF], values between 8 and 38 MPa are found for the Young modulus E of a Definity R UCA probed by an AFM. For an incompressible material the 2D compression modulus χ 2D is E d 0 , where d 0 is the shell thickness, estimated to be around 5 nm for such UCAs. This leads to 0.04 < χ 2D < 0.2 N/m, which is not in agreement with the values around 1 N/m found with the de Jong (Table I) or Marmottant (Table II) model. Note however, that static values of the shell may differ considerably from dynamic values measured in the MHz range.

In addition, experimental determinations have lead to unexpected dependencies of the viscoelastic parameters on shell radius, as also shown in Tables I This dependence on radius of the material properties was not substantiated by physical arguments, suggesting that extra modeling was required.

So far, the models have not considered the possible compressibility or anisotropy of the material constituting the shell. The purpose of the present article is to include these effects in the model of bubble oscillations and to quantify their influence on the linearized oscillation properties, i.e. the eigenfrequency and the damping coefficient.

IV. MODEL

We consider an encapsulated gas bubble immersed in an incompressible fluid with a density ρ f and a shear viscosity µ f . The effect of the liquid compressibility could be further included as described in the work of [START_REF] Prosperetti | The equation of bubble dynamics in a compressible liquid[END_REF]. The bubble shell is modeled as a visco-elastic solid of initial thickness d e . Furthermore, it was shown according to thin shell theory that a shell made of an homogeneous material with Poisson ratio ν may sustain a

(∆V /V ) b = 3(1 -ν) 1 + ν d 0 R 20 , (1)
before it buckles [START_REF] Hutchinson | Imperfection sensitivity of externally pressurized spherical shells[END_REF][START_REF] Quilliet | Numerical deflation of beach balls with various poisson's ratios: from sphere to bowl's shape[END_REF], where d 0 and R 20 are the shell thickness and external radius at rest, respectively. The first fraction is of order 1, except for exotic values of ν close to -1. Even for shells happening to be thicker than the commercially available ones, (∆V /V ) b is hence reasonably expected not to exceed 1/10. This point, plus recent experimental results having suggested that pressure-volume relationships obtained within the framework of thin shell theory apply also for thick shells [START_REF] Coupier | Let's deflate that beach ball[END_REF], indicates that we may safely consider, here and in the following, that linear elasticity framework is sufficient to describe the spherical behaviour of a wide range of UCAs in the unbuckled regime.

A. Quasi-static approximation

In the absence of body forces, the equation of motion in the solid [START_REF] Landau | Theory of Elasticity[END_REF] reads

ρ S ∂ 2 u ∂t 2 -∇ • σ = 0, (2)
where ρ S is the initial density of the solid, u the displacement field, and σ is the Cauchy stress tensor calculated on the actual configuration.

If we consider only elastic contributions to the stress, the dimensional analysis of (2)

shows that if the parameter ε = ω 2 0 d 2 e ρ s /E, that compares the orders of magnitude of the first and second term in (2), is small, then acceleration can be neglected (see e.g. [START_REF] Langtangen | Scaling of differential equations[END_REF]). Here E is a typical elastic constant of the material, and ω 0 is the (unknown) shell pulsation. Physically, √ ε is the ratio of the typical time scale τ 0 = d e ρ s /E needed for an elastic wave to travel across the shell thickness d e over the time scale ω -1 0 of the motion of the boundary. In general, E is not smaller than 100 MPa for a polymeric material where d e ∼ 100 nm, but for lipid shells of thickness of order 5 nm which are made of the type of anisotropic material that we treat later on in this paper, orders of magnitude as low as 100 kPa were proposed for an effective isotropic Young modulus [START_REF] Shafi | Probing phospholipid microbubbles by atomic force microscopy to quantify bubble mechanics and nanostructural shell properties[END_REF].

Hence, with ρ S ∼ 1000 kg/m 3 , τ 0 is expected to be smaller than 5 × 10 -10 s. This implies that, with ω 0 usually measured or found according to previous models lower than 10 MHz, ε is lower than 10 -4 .

The acceleration term can therefore be neglected for actual UCAs and will be so in the rest of this paper. The resolution of Eq.( 2) under this assumption will serve to determine the boundary conditions for the stress in the fluid, in order to determine its acceleration.

A problem similar to ours has been widely studied recently, that of a bubble oscillating in a liquid confined by a visco-elastic solid (Doinikov et al., 2018;Doinikov and Marmottant, 2018;[START_REF] Vincent | On the statics and dynamics of fully confined bubbles[END_REF]Wang, 2017). A simplifying hypothesis, that is used in [START_REF] Vincent | On the statics and dynamics of fully confined bubbles[END_REF] and Wang (2017) is to consider that the surrounding solid is not accelerated by the pressure waves. Here, we have shown that this hypothesis holds for our problem, due in particular to the thinness of the shells.

Note that the resonance frequency ω 0 is the unknown of this problem, so the validity of the hypothesis has to be checked a posteriori.

B. Stress-strain relation in the solid

We consider the shell as being made of a transverse isotropic material, i.e. whose properties in the orthoradial plane do not depend on the direction considered but can be different from that in the radial direction. The elastic properties of such a material are characterized by five independent elastic constants. The stress-strain relationship can be written as [START_REF] Lubarda | On the elastic moduli and compliances of transversely isotropic and orthotropic materials[END_REF]:

σ el ij = λ kk δ ij + 2µ ij + 2(µ 0 -µ) (δ i 0 i i 0 j + δ i 0 j i 0 i ) + α ( i 0 i 0 δ ij + δ i 0 i δ i 0 j kk ) + β δ i 0 i δ i 0 j i 0 i 0 , (3)
where is the strain tensor, λ is the first Lamé coefficient, µ the shear modulus in the plane of isotropy, µ 0 the out-of-plane shear modulus, and α and β two other coefficients. The direction i 0 points the axis of transverse isotropy. For an isotropic material, α = β = 0 and µ 0 = µ. For radial displacements, the elastic Cauchy stress tensor has only diagonal components given by

           σ el rr = (λ + 4µ 0 -2µ + 2α + β) rr + 2(λ + α) θθ σ el θθ = σ el φφ = (λ + α) rr + 2(λ + µ) θθ , (4)
with rr = ∂u/∂r and θθ = φφ = u/r, where u = u(r, t) is the Eulerian radial displacement in the shell.

The viscoelastic properties of the material are described by the generalized Kelvin-Voigt model [START_REF] Thompson | On the elasticity and viscosity of metals[END_REF][START_REF] Voigt | Ueber innere reibung fester körper, insbesondere der metalle[END_REF] where the complete strain tensor reads σ = σ el + σ visc , where σ visc is the viscous stress. For a transverse anisotropic material, integrating a thermodynamical consistent model [START_REF] Dalenbring | An explicit formulation of a three-dimensional material damping model with transverse isotropy[END_REF] based on the augmented Hooke's law (AHL) [START_REF] Dovstam | Augmented hooke's law in frequency domain. a three dimensional, material damping formulation[END_REF] in this fluid-structure interaction problem requires finite elements implementation. Another approach may be to consider viscosity effect for only some components of the stress tensor [START_REF] Lubarda | Viscoelastic response of anisotropic biological membranes. part II: Constitutive models[END_REF].

We will therefore consider two cases in this paper, both going one step further compared to the model by [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] that considers an isotropic and incompressible material:

1. A visco-elastic isotropic material, that can be compressible,

2. An anisotropic, purely elastic, material which is transversely isotropic and compressible.

In the case of an isotropic linear material, the elastic stress reads

           σ el rr = (K + 4 3 G ) rr + 2(K -2 3 G ) θθ σ el θθ = σ el φφ = (K -2 3 G ) rr + 2(K + 1 3 G ) θθ , (5)
where we have introduced the shear modulus G = µ and the bulk modulus

K = λ + 2 3 µ.
Both are a-priori functions of the oscillation frequency, which would call for the resolution of a self-consistency equation when the oscillation frequency will be eventually found as a function of, in particular, these elastic constants. We introduce the Kelvin-Voigt viscous stress σ visc whose expression is similar to that of the elastic stress:

           σ visc rr = (µ K + 4 3 µ G ) ˙ rr + 2(µ K -2 3 µ G ) ˙ θθ σ visc θθ = σ visc φφ = (µ K -2 3 µ G ) ˙ rr + 2(µ K + 1 3 µ G ) ˙ θθ (6)
The viscosities µ K and µ G describe the friction losses due to volume changes and shear, respectively. Little is known, in general, about the values of the loss moduli and, in particular, the "viscous Poisson ratio" whose definition may vary depending on the authors [START_REF] Lakes | On poisson's ratio in linearly viscoelastic solids[END_REF]. Its determination generally requires to perform two independent experiments aiming at determining, e.g., a shear loss modulus G and a traction loss modulus E (see, e.g. [START_REF] Guillot | Complete elastic characterization of viscoelastic materials by dynamic measurements of the complex bulk and young's moduli as a function of temperature and hydrostatic pressure[END_REF]). From a modelling perspective, one approach consists in following [START_REF] Lemaitre | Mechanics of solid materials[END_REF] where it is assumed , with no explicit justification that the ratio µ K /µ G is equal to K /G which amounts to say that the viscous Poisson ratio that would characterize a ratio of strain rates is equal to the elastic Poisson ratio that characterises the ratio of strains [START_REF] Linn | Derivation of a viscoelastic constitutive model of kelvin-voigt type for cosserat rods[END_REF][START_REF] Von Ende | On the thermodynamically consistent fractional wave equation for viscoelastic solids[END_REF]. Without this assumption, and considering an AHL model as in [START_REF] Tschoegl | Poisson's ration in linear viscoelasticity -a critical review[END_REF], [START_REF] Pritz | Relation of bulk to shear loss factor of solid viscoelastic materials[END_REF] have proposed bounds for the potential values of the loss moduli for materials with a positive Poisson ratio and a low enough shear loss factor. They show that 2/3 < K /G < 1 which, for a sinusoidal signal of given pulsation w 0 , amounts to the tight inequalities 2/3 < µ K /µ G < 1. We discuss these two assumptions in the discussion (Sec.

V J), but one should keep in mind that the difficulties in characterizing accurately two dissipation constants in viscoelastic materials, whose properties are often frequency dependent, must lead to consider the aforementioned relationships as pure hypotheses as for now.

V. ISOTROPIC COMPRESSIBLE SHELL

A. Deformation in the solid

The Eulerian radial displacement u within the shell is defined on the actual configuration as the variation from an unstrained position holding no stress within the shell:

u(r, t) = r -r e , ( 7
)
where r is the actual position of a material particle located at r e in the reference configuration.

The radial displacement u(r, t) is then calculated by solving Eq. ( 2) in the quasi-static approximation:

[∇ • (σ el + σ visc )] r = 0. ( 8
)
The ratio between the viscous and the elastic terms in the above equation is given by the ratio between the loss and storage moduli. Previous experimental studies on existing UCAs

show that the ratio between the viscosity and the storage modulus is of order 10 -8 -10 -9 s (see values in tables I-IV) therefore ω 0 τ S is often small, which we will take as a hypothesis in the following.

For an isotropic solid, from [∇ • σ] r = ∂σrr ∂r + 2 r (σ rr -σ θθ ) and using Eqs. ( 5) and ( 6), we are lead to solve:

∂ 2 ∂r 2 + 2 r ∂ ∂r - 2 r 2 (u + τ S u) = 0, (9)
with:

τ S = µ M M , M = K + 4 3 G , µ M = µ K + 4 3 µ G . ( 10
)
The solutions of Eq. ( 9) can be written:

u(r, t) = a(t)r + b(t) r 2 + A(r)e -t/τ S , (11)
where the term in A(r) characterizes the internal relaxation within the shell. Note that since ω 0 τ S is small, this term will marginally contribute to the overall response of the shell, and we shall therefore place ourselves in the conditions where it is zero.

The two variables a(t) and b(t) depend on the long time t τ S associated with the variations of the boundary conditions. We first express them as functions of R 1 (t) and R 2 (t), respectively the internal and external radii of the shell, which are our variables of R 1e and R 2e are the values of the radii in the unstrained case, and R 10 and R 20 their values at equilibrium in the fluid, which may differ from R 1e and R 2e , notably because of hydrostatic pressure. We find:

a(t) = R 2 2 (t)[R 2 (t) -R 2e ] -R 2 1 (t)[R 1 (t) -R 1e ] R 3 2 (t) -R 3 1 (t) , (12)
and:

b(t) = R 2 1 (t)R 2 2 (t) R 3 2 (t) -R 3 1 (t) × {R 2 (t)[R 1 (t) -R 1e ] -R 1 (t)[R 2 (t) -R 2e ]}. ( 13
)
Note that a and b are of order 1 relatively to the displacements at the boundaries, in agreement with the linear elastic theory used here to characterize the deformation tensor.

B. Velocity in the solid

The velocity field in the shell v s is the material derivative of the Eulerian displacement u(r, t):

v s (r, t) = Du Dt = ∂u ∂t + ∇u • v s (r, t). ( 14
)
For small deformations, |∇u| 1, the radial component of the velocity v s thus can be approximated to

v s ≈ ∂u ∂t = ȧr + ḃ r 2 , ( 15
)
where ȧ and ḃ are the time derivative of the variables a and b. Direct calculation of ȧ and ḃ, shown in Appendix, leads to expressions which violate the kinematic boundary conditions, i.e. v s (r = R 1 ) = U 1 and v s (r = R 2 ) = U 2 , where we define U 1 = Ṙ1 and U 2 = Ṙ2 .

However, the deviations from the kinematic boundary conditions remain of order |R i -R ie |/R ie , consistently with the assumption of small deformation and linear elastic behavior.

Hence, we can restrict the velocity to its leading order expression, where R 1e coincides with R 1 and R 2e with R 2 . This leads to an expression for v s that can otherwise be obtained directly from Eq. ( 15) by applying continuity condition at R 1 and R 2 :

v s (r, t) ≈ a v r + b v r 2 , ( 16
)
where

a v = R 2 2 (t)U 2 -R 2 1 (t)U 1 R 3 2 (t) -R 3 1 (t) , (17)
and

b v = [R 2 (t)U 1 -R 1 (t)U 2 ]R 2 1 (t)R 2 2 (t) R 3 2 (t) -R 3 1 (t) . ( 18
)
The above calculated displacement and velocity generalize the ones found in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] where an incompressible solid material is considered. Such materials are characterised by a traceless deformation tensor:

rr + θθ + φφ = 3a = 0, (19)
leading to the following displacement in the solid:

u inc (r, t) = R 2 1 (t)[R 1 (t) -R 1e ] r 2 , ( 20
)
where we have also used the relation a = 0 to reformulate the expression of b. The velocity v s given by Eq. ( 16) then becomes:

v inc s (r, t) = R 2 1 (t)U 1 r 2 . ( 21
)
Eqs. ( 20) and ( 21) are identical to the ones found in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] where the solid velocity was calculated directly from the law conservation of the mass for an incompressible fluid ∇ • v s = 0, while the displacement was deduced from volume conservation that reads, in the small deformation limit, ∇ • u = 0. Note that the two approaches are cross-consistent only in the small deformation framework: then in this case ∇ • u = 0, and moreover the 1/r 2 behaviour of v s is recovered only if Eq. ( 14) is approximated to Eq. ( 15).

Note finally that the displacement can also be defined on the reference configuration i.e.

using Lagrangian formalism without significant difference (Altenbach et al., 2008).

C. Equations of motion in the liquid

The conservation of mass for an incompressible fluid in a spherical coordinate system gives:

1 r 2 ∂ ∂r (r 2 v) = 0, ( 22
)
where v = (v f (r), 0, 0) is the radial Eulerian velocity vector in the fluid. For r = R 2 ,

v f (r = R 2 ) ≡ U 2 .
The velocity profile of the fluid is then:

v f (r) = U 2 R 2 2 r 2 . ( 23
)
The Navier-Stokes equation for an incompressible fluid and irrotational flow writes (Landau and Lifschitz, 1987):

ρ f ∂v f ∂t + v f ∂v f ∂r = - ∂P ∂r .
(24)

Integration of Eq. ( 24) between R 2 and +∞, using Eq. ( 23), leads to:

ρ f R 2 U2 + 3 2 U 2 2 = P f | r=R 2 -P ∞ , (25)
where P f | r=R 2 is the pressure in the fluid near the shell boundary, and P ∞ is the sum of the applied acoustic pressure P ac (t) and the ambient pressure P 0 .

In addition, conservation of radial momentum at the external surfaces of the shell imposes:

-P G (t) = (σ el rr + σ visc rr )| r=R 1 - 2γ 1 R 1 , (26) and:
(

σ el rr + σ visc rr )| r=R 2 = -P f | r=R 2 + σ f rr | r=R 2 - 2γ 2 R 2 , ( 27
)
where γ 1 and γ 2 are the surface tensions respectively at the internal and external boundaries of the shell, and P G (t) is the pressure of the gas inside the bubble. We assume the gas to obey a polytropic law, such that P G (t) = P g 0 (R 10 /R 1 ) 3κ , where P G 0 is the equilibrium gas pressure and κ is the polytropic exponent of the gas. The radial component σ f rr of the viscous stress equals:

σ f rr = 2µ f ∂v f ∂r = -4µ f U 2 R 2 2 (t) r 3 . ( 28
)
The normal stresses in the shell are obtained from Eqs. ( 5), ( 6), ( 11) and ( 15), noting that as done for the velocity, the strain rate ˙ is approximated in the linear elasticity limit 20 to: ˙ = ∂ /∂t, therefore the relation between ˙ and v s is similar to that between and u.

We have then:

σ el rr = 3K a -4G b r 3 , (29)
σ visc rr = 3µ K a v -4µ G b v r 3 . (30)
Inserting ( 29) and ( 30) in the first boundary condition ( 26) leads to a first equation for R 1 and R 2 :

-P G + 2 γ 1 R 1 = 3K a -4G b R 3 1 + 3µ K a v -4µ G b v R 3 1 . ( 31
)
Replacing a, b, a v and b v by their values ( 12), ( 13), ( 17) and ( 18) in the above equation, one eventually gets:

-P G + 2 γ 1 R 1 = 3K R 2 2 (R 2 -R 2e ) -R 2 1 (R 1 -R 1e ) R 3 2 -R 3 1 -4G [R 2 (R 1 -R 1e ) -R 1 (R 2 -R 2e )]R 2 2 (R 3 2 -R 3 1 )R 1 + 3µ K R 2 2 U 2 -R 2 1 U 1 R 3 2 -R 3 1 -4µ G (R 2 U 1 -R 1 U 2 ) R 2 2 (R 3 2 -R 3 1 )R 1 . (32)
We use the second boundary condition ( 27) to get rid of the unknown fluid pressure in Eq. ( 25), such that:

ρ f R 2 U2 + 3 2 U 2 2 = - 2γ 2 R 2 -P ∞ -4µ f U 2 R 2 -3K a + 4G b R 3 2 -3µ K a v + 4µ G b v R 3 2 . ( 33
)
This equation can be rewritten in a form that resembles a Rayleigh-Plesset equations by replacing the term R 2 -R 2e in a and b thanks to Eq. ( 32) :

ρ f R 2 U2 + 3 2 U 2 2 = -P 0 -P ac (t) -2 γ 2 R 2 -4µ f U 2 R 2 + P G -2 γ 1 R 1 1 - 4G 3K + 4G R 3 2 -R 3 1 R 3 2 -4G 3K 3K + 4G R 3 2 -R 3 1 R 3 2 R 1 -R 1e R 1 + 4 U 2 R 2 µ G 1 - 4G 3K + 4G -µ K 3G 3K + 4G -4 U 1 R 1 µ G 1 - 4G 3K + 4G -µ K 3G 3K + 4G R 3 1 R 3 2 . ( 34
)
In this expression, it is interesting to observe that the elastic contribution of the internal gas is modulated by the intrinsic elastic properties of the shell. This feature will disappear in the incompressible limit. Eq. ( 32) and Eq. (34) constitute a system of differential equations for the two unknowns R 1 and R 2 . For incompressible materials, the Rayleigh-Plesset equation is sufficient as R 1 and R 2 are simply linked through the incompressibility condition.

D. Unstrained Vs Initial Radii

As mentioned before, the unstrained radii may be different from the initial radii: R ie = R i0 . The radius R ie is defined by the unstrained state of the shell before it is plunged into the liquid, after what stresses within the shell take place, due to the surface tension at the interfaces, and the internal and external pressures.

Taking Eqs. ( 32) and (34) at equilibrium (U 1 = U 2 = 0, P ∞ = P 0 ), one can extract the displacements R i0 -R ie . They can be written R 1e = R 10 (1 + Z 1 ) and R 2e = R 20 (1 + Z 2 ), with:

Z 1 = 1 3K R 3 20 VS P 0 + 2γ 2 R 20 - R 3 10 VS P G 0 - 2γ 1 R 10 1 4G R 3 20 VS P 0 -P G 0 + 2γ 2 R 20 + 2γ 1 R 10 , (35)
Z 2 = 1 3K R 3 20 VS P 0 + 2γ 2 R 20 - R 3 10 VS P G 0 - 2γ 1 R 10 1 4G R 3 10 VS P 0 -P G 0 + 2γ 2 R 20 + 2γ 1 R 10 , (36)
where VS = R 3 20 -R 3 10 .

These formulations highlight the effect of compressibility, which is the same for the two radii.

If the shell is incompressible (K P G 0 , P 0 , γ i /R 10 ), one has:

Z inc i = P 0 -P G 0 + 2γ 1 R 10 + 2γ 2 R 20 R 3 20 + R 3 10 -R 3 i0 4G VS , (37)
which is identical to the expression found in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] when P G 0 = P 0 , which was hypothesized in that paper.

In [START_REF] Doinikov | Spatio-temporal dynamics of an encapsulated gas bubble in an ultrasound field[END_REF], where incompressible shells are also considered, the authors find the same relation as Eq. ( 37), which is the first order of their Eq. ( 33).

However in a second step they go further in the calculation using deformation profiles that are valid in the compressible case and find expressions (Eqs. ( 40) and (41) in their paper) that contradicts our findings, and the ones in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] and in [START_REF] Sarkar | Characterization of ultrasound contrast microbubbles using in vitro experiments and viscous and viscoelastic interface models for encapsulation[END_REF] in that they find the counter-intuitive result that surface tension tends to increase the equilibrium radius. Here, we are satisfied with the observation that an increase of surface tension leads to a shrinkage of the shell. This altogether suggests that care must be taken not to mix expressions from the compressible case with expressions from the incompressible case.

It is worth emphasizing that in the incompressible case, the ratio of the volume in the unstressed configuration to that after the shell is plunged in the fluid, namely [R 3 20 (1 +

Z inc 2 ) 3 -R 3 10 (1 + Z inc 1 ) 3 ]/ R 3 20 -R 3 10
, is equal to 1 in this model or in the other models [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF][START_REF] Doinikov | Spatio-temporal dynamics of an encapsulated gas bubble in an ultrasound field[END_REF][START_REF] Khismatullin | Radial oscillations of encapsulated microbubbles in viscoelastic liquids[END_REF] only to first order in Z inc i . This corresponds to the domain of validity of the linear elasticity framework.

In the general case, one must therefore restrict the obtained expressions to the first order in Z i for consistency.

E. Linear analysis

Assuming a small-amplitude oscillation, linear equations for the R i can be obtained using the following relations:

R 1 (t) = R 10 [1 + x(t)], |x(t)| 1; R 2 (t) = R 20 [1 + y(t)], |y(t)| 1; U 1 (t) = R 10 ẋ; U 2 (t) = R 20 ẏ; (38)
To the first order in x, y, Z 1 and Z 2 , Eq. ( 32) becomes, after using Eqs ( 35) and ( 36):

-3κP G 0 - 2γ 1 R 10 + 4G R 3 20 + 3K R 3 10 R 3 20 -R 3 10 x + (4G + 3K )R 3 20 R 3 20 -R 3 10 y - 4µ G R 3 20 + 3µ K R 3 10 R 3 20 -R 3 10 ẋ + (4µ G + 3µ K )R 3 20 R 3 20 -R 3 10 ẏ = 0. ( 39
)
Dividing this equation by K and taking the limit K → ∞, one gets x = yR 3 20 /R 3 10 , which is the relationship obtained for an incompressible material as in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF]. Eq. ( 39)

is therefore a generalization of this relationship for the case of a viscoelastic compressible material.

Eq. ( 39) together with the linearized Rayleigh-Plesset-like equation obtained from Eq.

(34) constitutes the following linear system:

M Ẍ + B Ẋ + KX = F (t), (40)
where:

X =       x y       , F (t) =       -P ac (t) 0       , M =       0 ρ f R 2 20 0 0       B =       b 11 b 12 -4µ G R 3 20 -3µ K R 3 10 R 3 20 -R 3 10 (4µ G +3µ K )R 3 20 R 3 20 -R 3 10       (41) with b 11 = 12 K µ G -G µ K R 3 10 R 3 20 3K + 4G , b 12 = 4 µ f + 3 K µ G -G µ K 3K + 4G , and
K =       k 11 -2γ 2 R 20 k 21 (4G +3K )R 3 20 R 3 20 -R 3 10       (42)
with

k 11 = 3κP G 0 - 2γ 1 R 10 1 - 4G 3K + 4G R 3 20 -R 3 10 R 3 20 + 12G K 3K + 4G R 3 20 -R 3 10 R 3 20 , k 21 = -3κP G 0 + 2γ 1 R 10 - 4G R 3 20 + 3K R 3 10 R 3 20 -R 3 10
.

The free oscillations of the shells (P ac = 0) are described by non-trivial harmonic solutions of the above system X = X 0 e λt , where λ = -δ + iω, that are obtained by setting det(λ 2 M + λB + K) = 0. This leads to a polynomial equation of order 3 for λ, which can be solved

analytically (yet leading to very long expressions) or numerically. This equation reads

c 1 λ 3 + c 2 λ 2 + c 3 λ + c 4 = 0, (43)
where

c 1 = -b 21 m 12 , c 2 = det[B] -k 21 m 12 , (44)
c 3 = b 11 k 22 -b 12 k 21 -b 21 k 12 + b 22 k 11 , c 4 = det[K].
For sake of comparison with the literature, and making use of the observation that damping coefficient is usually small, we present the leading order approximation and the first-order correction with respect to this damping coefficient in the following. Note that this is a sec-ond and independent approximation, based on the usual values of dissipation factors, that is added to that of small deformation.

F. Leading order approximation

For B = 0, Eq. ( 43) becomes det[K] -k 21 m 12 λ 2 = 0, therefore λ = iω 0 where the undamped resonance frequency ω 0 is given by:

ω 2 0 = 3K + 4G ρ f R 2 20 4G 3K 3K + 4G R 3 20 -R 3 10 R 3 20 + 3κP G0 - 2γ 1 R 10 1 - 4G 3K + 4G R 3 20 -R 3 10 R 3 20 × 3κP G0 - 2γ 1 R 10 R 3 20 -R 3 10 R 3 20 + 4G + 3K R 3 10 R 3 20 -1 - 2γ 2 ρ f R 3 20 . (45)
This constitutes the central result of this paper. The last term in the above expression is the classical contribution of the surface tension of the outer surface, which acts against an effective mass of fluid whose scale is given by the shell size. By contrast, the contribution of the shell elasticity and of the elastic forces acting on the inner side of the shell (the gas pressure and the surface tension) are strongly coupled. As discussed later on, this coupling disappears in the incompressibility limit. As in the Rayleigh-Plesset expression for a free bubble -

ω RP 0 = 1 ρ f R 2 20 3κP G0 -2γ 1 R 20 -2γ 2 R 20 1/2
, that is recovered here with Eq. (45) taken in the limit of vanishing shell volume (R 10 → R 20 ) -adding surface tension makes the shell pulsation decrease, at fixed P G0 . In practice, P G0 is not known nor measurable and it would be preferable to express the pulsation as a function of the external pressure P 0 . While this is easily done for a free bubble, leading to an increase of pulsation with surface tensions, this is more complex in the present situation: P G0 and P 0 also couple through the elastic stress within the shell, which depends on the reference configuration (R 1e , R 2e ), which is not known in general.

In this context, measuring oscillation frequency cannot be sufficient to determine the elastic constants of the shell material. Even if surface tensions are assumed to be zero, and considering that the external radius is known, we are left with four unknowns which are the two elastic constants, the internal pressure and the internal radius. This is one more than in Church model and two more than in zero-thickness shell models. Even in these simpler model, and in all cases, one needs to know more on the fabrication process of the shell to know their stress-free state or to make additional assumptions. In [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF], it is for instance assumed that permeability of the shell under study allows to consider that P G0 = P 0 , which may be true for thin lipid shells, but not for thicker shells, as pointed out in [START_REF] Doinikov | Spatio-temporal dynamics of an encapsulated gas bubble in an ultrasound field[END_REF].

For an incompressible shell, the undamped natural frequency becomes:

ω inc 0 = (ρ S R 2 10 α inc ) -1/2 3κP G0 - 2γ 1 R 10 - 2γ 2 R 20 R 3 10 R 3 20 +4G R 3 20 -R 3 10 R 3 20 1/2 , with α inc = ρ f ρ S R 10 R 20 . ( 46
)
This differs from the expression proposed in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF]:

ω Ch 0 = (ρ S R 2 10 α Ch ) -1/2 3κP G0 - 2γ 1 R 10 - 2γ 2 R 20 R 3 10 R 3 20 +4G R 3 20 -R 3 10 R 3 20 1 + Z Ch 1 1 + 3R 3 10 R 3 20 1/2 , ( 47
)
with

α Ch = ρ f ρ S R 10 R 20 + 1 - R 10 R 20 . ( 48
)
A first difference lies in the effective mass characterized by the coefficient α, since we neglected the inertia of the shell. Note that it introduces a correction on ω 2 0 of order d 0 /R 20 , where d 0 = R 20 -R 10 , that is of at most a few percent for actual UCAs.

The other difference lies in the presence of a Z Ch 1 term in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF]. This is due to a subtle inconsistency in the linearizing process: as discussed in Sec. V D, Z Ch 1 must be considered as a small parameter in order to keep the validity of linear elasticity framework. It characterizes the difference between the unstrained state and the equilibrium state, the same way as x and y in Eq. 38 characterize the difference between the actual and the equilibrium state. Terms like xZ 1 should therefore not be included in the linearized equation, contrary

to what is done in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] between his Eqs. ( 12) and ( 17). Replacing Z Ch 1 by its value in ω Ch 0 one gets:

ω Ch 0 = (ρ S R 2 10 α Ch ) -1/2 3κP G0 + 2γ 1 R 10 3R 3 10 R 3 20 + 2γ 2 R 20 1 + 2R 3 10 R 3 20 + 4G R 3 20 -R 3 10 R 3 20 1/2 . ( 49
)
One can see that the contributions of the surface tension are incorrectly estimated with this contested expression by Church, as this expression does not converge to the Rayleigh-Plesset pulsation ω RP 0 in the vanishing volume limit.

G. First order approximation

If Eq. ( 43) is expanded to the first orders in b ij and δ, one gets that ω = ω 0 and which can be reformulated as:

δ = - c 3 + b 21 m 12 ω 2 0 2k 21 m 12 , (50)
δ = 1 2m 12 b 12 -b 11 k 22 k 21 + k 11 k 2 21 (b 21 k 22 -b 22 k 21 ) . ( 51
)
As seen in Eq. ( 40 For an incompressible shell, the damping ratio δ in Eq. ( 50) simply becomes:

δ inc = 2 (R 3 20 -R 3 10 )µ G + R 3 10 µ f ρ S R 2 10 R 3 20 α inc . ( 52
)
In [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] it reads:

δ Ch = 2 (R 3 20 -R 3 10 )µ G + R 3 10 µ f ρ S R 2 10 R 3 20 α Ch . ( 53
)
As for the pulsation a difference of a few percents remains, which is related to the absence of shell mass in our model.

H. Discussion: effect of compressibility on ω 0

We discuss in this section to which extent the frequency is modified when the material is compressible. We first consider a reference configuration, denoted R, which is considered in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF]: d 0 = 15 nm, P G 0 = 101.3 kPa, ρ f = 1000 kg/m 3 , ρ S = 1100 kg/m 3 , µ f = 0.001 Pa.s, G = 88.8 MPa, γ 1 = 0.04 N/m, γ 2 = 0.005 N/m and κ = 7/5. For such a shell whose external radius lies in the range 1 -10 µm, we find that 0.99 < ω inc 0 /ω Ch 0 < 1, which indicates that while our model has led us to neglect the inertia of the shell, this assumption will modify the final result by a negligible amount. Note that in this example, since γ i /R i0 G , the inaccuracy that we exhibited in the [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] model has no quantitative consequence. In the following, we consider ω inc 0 as the reference value for discussion.

We now discuss the effect of compressibility together with an evaluation on the impact of the contribution of gas compressibility. For most commercial shells, G is actually 10 to 1000 times the ambient pressure (see e.g. Table III). Since the G contribution is weighted by d/R (which is roughly the ratio of the shell material volume over the volume of gas), both contributions are likely to contribute with comparable weight.

As the contribution of the external surface tension is purely additive, for simplicity we set γ 2 = 0 and consider several values of P = P G 0 -2γ 1 3κR 10 , that characterizes the contribution of the inner gas to UCA oscillations. In this case, regarding space variables, R 20 × ω 0 depends only on d 0 /R 20 .

In Fig. 1, the ratio ω 0 /ω inc 0 is calculated for different values of G and K , that are set relatively to P .

The influence of compressibility is significant: it reduces the resonance frequency, all the more that the relative thickness d 0 /R 20 increases, the Poisson ratio decreases, and the shear modulus increases. In particular, when the bubble radius decreases at fixed thickness, this effect of compressibility will become relatively more important. Compressibility thus and2 3 G , and the corresponding Poisson

P G 0 -2γ 1 3κR 10 . We varied K as 100 G , 10 G , 3 G , 2 G , G
ratios ν = (3K -2G )/(6K + 2G ) = (3µ K -2µ G )/(6µ K + 2µ G )
are shown on each curve. Note that curves do not vary by more than 1% for G ≥ 5 P .

introduces a dependency of the frequency on the shell radius that is more complex than in the incompressible case, where ω 0 ∝ 1/R 20 in the thin shell limit. and2 3 G , and the corresponding Poisson ratios ν = (3K -2G )/(6K + 2G ) are shown on each curve.Note that curves do not vary by more than 1% for G ≥ 5 P . In all cases, µ K = 0.7µ G , following [START_REF] Pritz | Relation of bulk to shear loss factor of solid viscoelastic materials[END_REF].

P G 0 -2γ 1 3κR 10 . We varied K as 100 G , 10 G , 3 G , 2 G ,G
For large values of G , Eq. ( 45) yields, after setting

K = 2(1+ν) 3(1-2ν) G ≡ f (ν)G , 471 ω 2 0,G P = 1 ρ f R 2 20 4G R 3 20 -R 3 10 R 3 20 3f (ν) 4 + 3f (ν) R 3 10 R 3 20 . (54)
Then,

ω 0,G P ω inc 0,G P 2 = 3f (ν)(1 -d 0 R 20 ) 3 4 + 3f (ν)(1 -d 0 R 20 ) 3 ≡ g(ν, d 0 , R 20 ).
(55)

In the thin shell limit, compressibility leads to a decrease of the pulsation squared by a factor f (ν)/[f (ν) + 4/3].

Eq. ( 54) can also be interpreted from the following practical viewpoint: if one measures a shell pulsation and deduces from this measurement a value G 0 for the shell, assuming incompressibility, the same measurement can also be obtained with a shell of shear modulus

G and Poisson ratio ν obeying G 0 = G g(ν, d 0 , R 20 ).
The consequences are two-fold: as g(ν, d 0 , R 20 ) is significantly smaller than 1 as soon as ν < 1/2, the existence of unforeseen compressibility will lead to an underestimation of the shear modulus. For instance, for a shell of estimated thickness 15 nm and external radius 2 µm, if ν happens to be 0.4 instead of 0.5, g(ν, d 0 , R 20 ) = 0.77, which means that the shear modulus will be underestimated by 23%. This value reaches 28% if d 0 = 200 nm.

Second, as g is an increasing function of R 20 , using a model for incompressible material can lead to an artificial increase of the (apparent) shear modulus with the radius, a feature regularly pointed out in the literature.

These compressibility effects are more pronounced for thick shells, and we are not aware of oscillation measurements in the literature based on thick shells like polymeric shells. In addition, a more quantitative analysis of the impact of compressibility on the radius dependency of the frequency, by comparison with other suggestions like non-linear effects, requires to know more about the inner pressure inside the considered shells, which depends on their manufacturing process and also potentially on the allotted time for pressure equal-ization through transmembrane diffusion. This point becomes even more evident in the zero-thickness shell limit that is discussed in the following.

I. Discussion: from finite thickness to zero-thickness shell

For vanishing thickness, considering the corresponding limit in our finite thickness model or in that of [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF] leads to models that can be compared to zero-thickness models.

In particular, in [START_REF] Hoff | Oscillations of polymeric microbubbles: Effect of the encapsulating shell[END_REF], the vanishing thickness limit of Church model is considered and the resulting frequency is shown to be similar with that obtained in de [START_REF] De Jong | Absorption and scatter of encapsulated gas filled microspheres: theoretical considerations and some measurements[END_REF] or in the linearized version of [START_REF] Marmottant | A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture[END_REF]:

ω 0-thickness 0 = (ρ f R 2 0 ) -1/2 3κP G0 + 4 χ 0 R 0 , 1/2 , ( 56
)
where R 0 is the shell radius and χ 0 has the dimension of a surface tension and includes in-plane elasticity effects as well as surface tension effects on both sides of the interface [START_REF] De Jong | Absorption and scatter of encapsulated gas filled microspheres: theoretical considerations and some measurements[END_REF][START_REF] Hoff | Oscillations of polymeric microbubbles: Effect of the encapsulating shell[END_REF][START_REF] Marmottant | A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture[END_REF][START_REF] Sarkar | Characterization of ultrasound contrast microbubbles using in vitro experiments and viscous and viscoelastic interface models for encapsulation[END_REF][START_REF] Van Der Meer | Microbubble spectroscopy of ultrasound contrast agents[END_REF]. In [START_REF] Hoff | Oscillations of polymeric microbubbles: Effect of the encapsulating shell[END_REF] when surface tension effects are neglected, χ 0 is shown to be equal to χ 2D = 3G d 0 , the in-plane surface contraction modulus.

We examine here the small thickness limit of our model. We consider only the incompressibility limit, which is already an interesting source for discussion and allows direct comparison with the actual zero-thickness models.

Keeping only the 0 th and 1 st orders in d 0 /R 20 in Eq. ( 45), we find the following expansion:

ω inc 0 = (ρ f R 2 20 ) -1/2 3κP G0 - 2γ 1 R 20 - 2γ 2 R 20 + 12G + 9κP G0 - 8γ 1 R 20 d 0 R 20 + o d 0 R 20 1/2 . ( 57
)
By comparison with Eq. ( 56), this introduces a correction that implies that pressure and inner surface tension have a more complex space-dependency than that proposed in the Church-Hoff model, where the 1 st order in d 0 /R 20 was neglected in the inertial term.

We attempt to discuss the implication of our modeling regarding the interpretation of experimental data. Authors generally consider a given experiment for a set of shells of different sizes, which they either watch (measuring thus the radius oscillation [START_REF] Chetty | Highspeed optical observations and simulation results of sonovue microbubbles at low-pressure insonation[END_REF][START_REF] Doinikov | Modeling of nonlinear viscous stress in encapsulating shells of lipid-coated contrast agent microbubbles[END_REF][START_REF] Li | Modeling complicated rheological behaviors in encapsulating shells of lipid-coated microbubbles accounting for nonlinear changes of both shell viscosity and elasticity[END_REF][START_REF] Tu | Estimating the shell parameters of sonovue R microbubbles using light scattering[END_REF][START_REF] Van Der Meer | Microbubble spectroscopy of ultrasound contrast agents[END_REF] or listen (measuring thus the acoustic transmission [START_REF] Parrales | Acoustic characterization of monodisperse lipid-coated microbubbles: Relationship between size and shell viscoelastic properties[END_REF]). The obtained curves are then fitted according to the chosen model, which results in the determination of the corresponding elastic modulus for each shell radius. It is then generally observed that this constant increases with the radius, which highlights the limit of the chosen model.

Other parameters are generally considered as known but they are not always given by the authors. In particular, the inner pressure P G0 is sometimes set to atmospheric pressure with not much justification [START_REF] Doinikov | Modeling of nonlinear viscous stress in encapsulating shells of lipid-coated contrast agent microbubbles[END_REF][START_REF] Li | Modeling complicated rheological behaviors in encapsulating shells of lipid-coated microbubbles accounting for nonlinear changes of both shell viscosity and elasticity[END_REF][START_REF] Tu | Estimating the shell parameters of sonovue R microbubbles using light scattering[END_REF][START_REF] Van Der Meer | Microbubble spectroscopy of ultrasound contrast agents[END_REF] but some authors do not specify their choice [START_REF] Parrales | Acoustic characterization of monodisperse lipid-coated microbubbles: Relationship between size and shell viscoelastic properties[END_REF]. On the other hand, the descriptions of fabrication processes of home-made microshells often mention initial gaz pressure larger than 1 bar [START_REF] Parrales | Acoustic characterization of monodisperse lipid-coated microbubbles: Relationship between size and shell viscoelastic properties[END_REF][START_REF] Segers | Stability of monodisperse phospholipid-coated microbubbles formed by flow-focusing at high produc-tion rates[END_REF], which questions the hypothesis of atmospheric pressure inside the shells. Though diffusion may favor this hypothesis, such a phenomenon will also induce stresses inside the shell reaching its new equilibrium, resulting in uncertainties on the exact state around which the oscillations take place.

Finally, it is generally observed in all papers that while the radius varies by a factor 2 to 3, the corresponding elastic modulus varies by a factor 3 to 4. In [START_REF] Parrales | Acoustic characterization of monodisperse lipid-coated microbubbles: Relationship between size and shell viscoelastic properties[END_REF] this is the case but contrary to most other papers where only the values of the elastic constants are given, the measured frequencies are also mentioned. We therefore use these raw data to make the following comments. In Fig. 4, the pulsations found in the experiments are plotted as a function of shell radius. Those shells are lipidic shells, therefore the small-thickness limit holds. The fit of their data by the usual zero-thickness law (Eq. 56), assuming P G0 = 1 bar and κ = 1.4, is not that good, which illustrates the conclusions of the authors who, considering each radius separately, showed that the elastic modulus must be an increasing function of the radius. We note however that the fit yields χ 0 = 0.21 N/m, which is smaller than all the values reported by the author for the different shell radii, which questions the (implicit) choice of inner pressure or of polytropic constant they made. Interestingly, letting P G0 free leads to a better fit, with P G0 = 1.6 bar. This shows the importance of the knowledge of the inner pressure or, equivalently, of the polytropic coefficient that depends on the chosen gas and on the details of the thermodynamics process, as discussed in [START_REF] Parrales | Acoustic characterization of monodisperse lipid-coated microbubbles: Relationship between size and shell viscoelastic properties[END_REF].

In the expression for the zero-thickness limit that we established (Eq. 57), we show that the contribution of pressure is more complex and that it is important to decouple, in the elastic contribution of the interface, bulk effects from surface tension effects: they do not sum up in a simple χ 0 parameter. Using this expression we find an even better fit for Full black line: fit with Eq. ( 56) with fixed inner pressure P G0 = 1 bar and χ 0 as a free parameter.

Red dotted line: fit with the same equation but the pressure is also a free parameter. Blue dashed line: fit with Eq. ( 57) with also the surface tension being a free parameter. Using Eq. ( 57) allows to recover the full spatial dependency of the data, with a 1/R 0 and a 1/R 3/2 0 contribution.

the data of [START_REF] Parrales | Acoustic characterization of monodisperse lipid-coated microbubbles: Relationship between size and shell viscoelastic properties[END_REF], illustrating thus the complex interplay between all the parameters of these models. Note that we do not claim here that the parameters we find are those that actually characterize the considered shell. Our discussion simply highlights the need for a good knowledge of a maximum of parameters, if one wishes to extract one unknown parameter from the sole measurement of oscillation frequencies.

J. Discussion: Effect of compressibility on the damping Reminder: our model assumes that τ S = µ M /M is much smaller than ω -1 0 . In the following examples, we checked that τ S ω 0 is always lower than 0.01. We set here the fluid viscosity µ f = 0.001 Pa.s and the shear viscosity µ G = 0.002 Pa.s.

As for the discussion on pulsation, we set γ 2 = 0 and consider several values of P = P G 0 -2γ 1 3κR 10 . In Fig. 2 we show the ratio of the damping constants δ/δ inc under the hypothesis that µ K varies with µ G the same way K varies with G that is, the viscous and elastic

Poisson ratio are equal [START_REF] Lemaitre | Mechanics of solid materials[END_REF]. In Fig. 3 µ K is chosen to be equal to 0.7µ G following [START_REF] Pritz | Relation of bulk to shear loss factor of solid viscoelastic materials[END_REF], where it is shown that 2/3 < µ K /µ G < 1 for thermodynamic consistency.

Compressibility has the effect to make the damping constant decrease. As for the elastic constant determined through the frequency, this may lead to an underestimation of the shear viscosity if an incompressible model is used. For large values of G and K = G f (ν) compared to P , as for the frequencies, the damping depends only on elastic properties through the Poisson ratio, as can be seen through Eqs. ( 41) and (42) taken in the limit G , K P . In practice, as seen in Figs. 2 and 3, this limit is reached as soon as G > 5 P , which is generally the case for actual commercial UCAs. Interestingly, the choice of the model for the viscous Poisson ratio has little impact on the final results: for high values of K , µ K is not expected to contribute much for both models as its contribution vanishes in the incompressibility limit (see Eq. ( 52)), and for values of K comparable G , µ K is close to µ G in both models. If the [START_REF] Pritz | Relation of bulk to shear loss factor of solid viscoelastic materials[END_REF] model is assumed though, the coupling between elastic and viscous terms is such that the damping is not a monotonous function of the Poisson ratio ν.

We now examine the effect of anisotropy in the properties of a purely elastic material. We reformulate Eq. ( 4) using elastic constants corresponding to standard deformations [START_REF] Itskov | Elastic constants and their admissible values for incompressible and slightly compressible anisotropic materials[END_REF][START_REF] Lempriere | Poisson's ratio in orthotropic materials[END_REF]:

             σ el rr = (1-ν )E r 1-ν -2 E E r ν 2 θr rr + 2ν θr E 1-ν -2 E E r ν 2 θr θθ σ el θθ = σ el φφ = ν θr E 1-ν -2 E E r ν 2 θr rr + E 1-ν -2 E E r ν 2 θr θθ , (58)
where E r is the Young modulus for traction in the radial direction while E is the Young modulus in the orthoradial plane. ν is the Poisson ratio in this same plane, and ν θr the Poisson ratio governing deformations in the orthoradial plane when there is a radial load 1 .

Thermodynamical consistency imposes [START_REF] Lempriere | Poisson's ratio in orthotropic materials[END_REF]:

-1 ≤ ν ≤ 1 -E r /E ≤ ν θr ≤ E r /E ν ≤ 1 -2ν 2 θr E E r . ( 59
)
For an isotropic material of Poisson ratio ν, these inequalities reduce to -1 ≤ ν ≤ 1/2.

The case ν = 1/2 corresponds to incompressible material as considered in [START_REF] Church | The effects of an elastic solid surface layer on the radial pulsations of gas bubbles[END_REF].

A. Displacement within the shell

Following the same steps as in Sec. V A, the displacement now obeys the following equation:

d 2 u dr 2 + 2 r du dr - 2ku r 2 = 0, (60)
with:

k = E (1 -ν θr ) E r (1 -ν ) , (61)
which is the equivalent of Eq. ( 9) for this purely elastic case.

The solutions of Eq. ( 60) have the form

u T r (r) = a T r r β + + b T r r β -, (62)
with

β ± = 1 2 (-1 ± √ 1 + 8k).
Note that by virtue of Eqs. ( 59), it can be shown that k ≥ -1/8 whatever the material properties and the exponents β ± are real. The isotropic case corresponds to k = 1 then β -= -2 and β + = 1. The variables a T and b T are related to the boundary conditions thanks to Eq. ( 7):

a T r = (R 2 -R 2e )R β - 10 -(R 1 -R 1e )R β - 20 R β - 10 R β + 20 -R β + 10 R β - 20 , (63) and b
T r = (R 1 -R 1e )R β + 20 -(R 2 -R 2e )R β + 10 R β - 10 R β + 20 -R β + 10 R β - 20 . (64)
The Rayleigh-Plesset-like equation can be derived following the same steps as in Sec. V C.

For the sake of simplicity we calculate directly the resonance frequency in the following section.

B. Resonance frequency

Following the same steps as in Sec. V E, one gets the following system:

M Ẍ + KX = F (t), (65)
where

X =       x y       , F (t) =       -P ac (t) 0       , M =       0 ρ f R 2 20 0 0       , and
K T r =       k T r 11 -2γ 2 R 20 k T r 21 k T r 22      
, where

k T r 11 = 3κP G 0 - 2γ 1 R 10 1 -2 (R β + -1 20 -R β + -1 10
)(β + -1) -1 R β - 10 A + -(R β --1 20 -R β --1 10
)(β --1) -1 R β + 10 A - E r (1 -ν )(β + -β -) + 2(R β + -1 20 -R β + -1 10
)(β + -1) -1 R β - 10 A + [E r (1 -ν )β -+ 2ν θr E ] 1 -ν -2 E E r ν 2 θr E r (1 -ν )(β + -β -) - 2(R β --1 20 -R β --1 10
)(β --1) -1 R β + 10 A -[E r (1 -ν )β + + 2ν θr E ] 1 -ν -2 E E r ν 2 θr E r (1 -ν )(β + -β -) , (66 with:
A + = β + (1 -ν )E r -[1 -(2 -β + )ν θr ]E , (67)
A -= β -(1 -ν )E r -[1 -(2 -β -)ν θr ]E , (68)
and:

k T r 21 = -3κP G 0 + 2γ 1 R 10 + R β + 20 R β - 10 [E r (1 -ν )β -+ 2ν θr E ] 1 -ν -2 E E r ν 2 θr (R β - 10 R β + 20 -R β + 10 R β - 20 ) - R β - 20 R β + 10 [E r (1 -ν )β + + 2ν θr E ] 1 -ν -2 E E r ν 2 θr (R β - 10 R β + 20 -R β + 10 R β - 20 ) , ( 69
)
k T r 22 = R 20 R β -+β + -1 10 E r (1 -ν )(β + -β -) 1 -ν -2 E E r ν 2 θr (R β - 10 R β + 20 -R β + 10 R β - 20 ) . ( 70
)
Then, the undamped resonance frequency is: In what follows, the couple (E = 2G (1 + ν),ν) will be used to describe the elastic properties of an isotropic solid for the sake of comparison with the elastic properties of a transversely isotropic material.

ω T r 0 = det[K T r ] m 12 k T
We first discuss what is the impact of anisotropy keeping the material incompressible. In such a situation, it is shown in [START_REF] Itskov | Elastic constants and their admissible values for incompressible and slightly compressible anisotropic materials[END_REF] that ν θr = 1/2 while ν = 1 -

E 2E r
.

Thermodynamics constraints (Eq. ( 59)) then impose E r > E /4.

We remark that Eq. ( 61) becomes k = 1 that is, the deformation function is the same as in the isotropic case. Second, the terms in the K T r matrix implying the shell elastic constants are all proportional to E r and do not depend on E . We conclude that incompressible shells oscillate exactly the same way whatever the value of their in-plane Young modulus that is, they oscillate like isotropic incompressible shells of Young modulus E r

For anisotropic compressible material, in general, hydrostatic stress does not necessarily induces a uniform dilatation in the three directions. It is interesting for comparison with the isotropic case to consider the situation where this is true. In such a situation of isotropic volumetric response, a bulk modulus can be defined as a material constant [START_REF] Itskov | Elastic constants and their admissible values for incompressible and slightly compressible anisotropic materials[END_REF]. It is equal to κ = E r /[3(1 -2ν θr )]. In addition, it can be shown that ν is given by 1 -

E 2 1 E r + 1 3κ
; the material properties are thus described, for the radial motion considered here, by 3 independent variables (e.g. κ, E r and E ) instead of 4 in the general case and 2 in the incompressible case. In this case also, the pulsation is that of the isotropic material of moduli K = κ and E = E r , though the deformation inside the shell is not the same: in Eq. ( 61), k is independent from E (but not necessarily equal to 1) and because ν appears only under the pattern 1 -ν in K T r , it can be easily seen that the contribution of E vanishes.

Finally, in Fig. 5, we consider a general (arbitrary) case, based on the test case the configuration R as it refers to a lipidic shell, which we may expect to exhibit transverse anisotropic properties. We have fixed ν = ν θr = 0.35 and varied E as a function of E r , within the bounds allowed by thermodynamics. Here, E also influences the frequency, which increases as E increases.

In all cases, as for the isotropic case, these results show that for a given measure of pulsation frequency, several sets of elastic parameters can yield the same result. More complex dependency with the external radius is also expected.

We have developed a finite thickness shell model for the oscillations of an encapsulated bubble whose material can be compressible and/or present different elastic properties in the radial and orthoradial directions. The main hypothesis is that we have neglected the mass of the shell, leading to infinite velocity for wave propagation in the material, in order to simplify the equations. The next step would be to consider the complete problem of wave propagation, as done for instance in (Doinikov et al., 2018;Doinikov and Marmottant, 2018) for a bubble oscillating in a liquid confined by a visco-elastic solid.

We have found exact expressions for the free pulsation of an encapsulated bubble, that could be used to interpret more accurately experimental characterization of UCAs. Our results suggest that neglecting compressibility will lead to underestimation of the shear modulus, and that adding some compressibility in the model may explain the apparent growth of the elastic moduli with the shell radius.

Due to the growing interest in the development of new generation UCAs, made of various material and built with well defined radii, we expect that several opportunities to test our model will emerge in the near future. Our predictions can also be used to build more complete theories accounting for the response of the shells to external signals.

  -III. van der Meer et al. (2007) observed a dependence of the shell viscosity on the initial bubble radius using Marmottant model for BR-14 R . Chetty et al. (2008) measured an increase of the shear modulus G with the radius using Church-Hoff model for Sonovue R . Tu et al. (2009) and Li et al. (2013) measured an increase of the elasticity and the viscosity parameters of the shell with the shell radius, using the linearized Marmottant model for Sonovue R . Identical observation where made in Doinikov et al. (2009) (lipid encapsulated bubbles with De Jong model), Helfield and Goertz (2013) (Definity R with Marmottant model), and Parrales et al. (2014) (home-made monodisperse encapsulated microbubbles with the linearized Marmottant model).

  FIG. 1. Ratio of the undamped resonance frequencies ω 0 /ω inc 0 as a function of d 0 /R 20 , in the

  ), the first term b 12 in the above expression represents the damping directly affecting the motion of the external radius of the shell, through the fluid viscosity and a contribution of the shell viscosity. The second term stems from the damping affecting the motion of the internal radius, which is weighted by the elastic contribution k 22 /k 21 . The third term stems from the coupling between dissipation and elastic deformation inside the shell.

  FIG.2. Ratio of the damping ratios δ/δ inc as a function of d 0 /R 20 , in the absence of external

  FIG. 3. Ratio of the damping ratios δ/δ inc as a function of d 0 /R 20 , in the absence of external

  FIG. 4. Dots: experimental pulsations found in Parrales et al. (2014) as a function of shell radius.

  : effect of anisotropy on ω 0

FIG. 5 .

 5 FIG. 5. Ratio of the undamped resonance frequencies in the compressible case ω T r 0 /ω iso 0 . The

TABLE I .

 I Shell properties estimations using de Jong model[START_REF] De Jong | Absorption and scatter of encapsulated gas filled microspheres: theoretical considerations and some measurements[END_REF]. f and p are the characteristic frequencies and amplitudes of the acoustic waves used to excite the UCAs. S p and S f are the elastic and viscous ad-hoc parameters that are introduced in the model. The intervals for the viscoelastic parameters correspond to cases where dependency on the radius was reported.

	UCA	R 20 (µm) f (MHz)	p(kPa) S p (N/m) S f (10 -6 N/m.s) Method	Reference
	SonoVue R 0.6 -4.5 1 -10	< 10	0.35 -2.61	0.46 -3.42	Attenuation	(Gorce et al., 2000)
	Albunex R 2.5 -6 0.7 -12.5 not known	8	4	Attenuation (de Jong and Hoff, 1993)
	Definity R 0.5 -2.5 12 -28	25	1.71	0.015	Attenuation (Goertz et al., 2007)
		1 -3	7 -15	25	1.64	0.15	Attenuation	(Faez et al., 2011)

TABLE II .

 II Shell properties estimations using Marmottant model[START_REF] Marmottant | A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture[END_REF]. f and p are the characteristic frequencies and amplitudes of the acoustic waves used to excite the UCAs.

	Analysis of experiments through the model allow to determine the 2D compression modulus χ 2D
	and the surface dilatationnal viscosity κ S . In the linearized version of the model, they are related
	to the constants introduced by de Jong through S p = 2χ 2D and S f = 12πκ S . The intervals for the
	viscoelastic parameters correspond to cases where dependency on the radius was reported.
	UCA	R 20 (µm) f (MHz) p(kPa) χ 2D (N/m) κ S (10 -8 N/m.s)	Method	Reference
	SonoVue R	0.975	2.9	130	1	1.5	Spectroscopy (Marmottant et al., 20
		0.8 -3.25	2.5	150 0.024 -0.87 a	0.1 -3	Light scattering	(Tu et al., 2009)
		0.75 -3.25	2.5	150	0.39 -0.55	0.05 -2	Light scattering	(Tu et al., 2011)
		0.8 -3.25	2.5	150	0.4 -0.55	0.1 -3	Light scattering	(Li et al., 2013)
	BR14 R	1.9	1.5 -2.5 < 40	0.54 a	2.3	Spectroscopy (van der Meer et al., 20
	Definity R	0.72 -1.4	1	308	0.5 -0.97	0.01 -0.9	Light scattering	(Tu et al., 2011)
		1.4 -2.8 4 -13.5 6 -25 0.5 -2.5	0.02 -0.6	Spectroscopy (Helfield and Goertz, 20
	Home-made 2.9 -6.3 0.5 -4	-	0.28 -0.85	3 -6	Attenuation	(Parrales et al., 2014
	lipid shell
	a Linearized version of the model

TABLE III .

 III Shell properties estimations using Church-Hoff model[START_REF] Hoff | Oscillations of polymeric microbubbles: Effect of the encapsulating shell[END_REF]. The thickness d 0 is an estimation which is made in each paper. f and p are the characteristic frequencies and

	amplitudes of the acoustic waves used to excite the UCAs. The intervals for the viscoelastic
	parameters correspond to cases where dependency on the radius was reported.
	UCA	R 20 (µm) d 0 (nm) f (MHz) p(kPa) G (MPa) µ G (Pa.s)	Method	Reference
	SonoVue R	1.78	4	2.5	150	20	0.6	Light scattering	(Tu et al., 2009)
		3 -5.5	2.5 6.8 -7.3	40	1.9 -105	1	Microscopy	(Chetty et al., 2008)
	Sonazoid R	1.6	4	2 -6 300 -800	52	0.99	Attenuation	(Sarkar et al., 2005)
	Optison R	1.5	5 -10 3.6 -4.3	100	20.7	1.7	Attenuation (Chatterjee and Sarkar, 20

TABLE IV .

 IV Shell properties estimations using Sarkar model[START_REF] Sarkar | Characterization of ultrasound contrast microbubbles using in vitro experiments and viscous and viscoelastic interface models for encapsulation[END_REF]. f and p are the characteristic frequencies and amplitudes of the acoustic waves used to excite the UCAs. E S and κ S are the surface dilatational elasticity and viscosity respectively introduced in the model.

		UCA	R 20 (µm) f (MHz) p(kPa) E S (N/m) κ S (10 -8 N/m.s) Method	Reference
		Sonazoid R	3.2	2 -6 200 -600	0.51	1	Attenuation (Sarkar et al., 2005)
		Home-made 0.7 -1.5 2.5 -3 100 -150	0.02	0.85	Attenuation (Paul et al., 2013)
		PLA shell
	148	maximum relative loss of volume

ACKNOWLEDGMENTS

The authors appreciated fruitful discussions with Philippe Marmottant and Salvatore Federico.

APPENDIX: DEFORMATION VELOCITIES

We present here the expressions obtained for ȧ and ḃ by directly derivating a and b, see

and:

1 Note that ν rθ is used by some authors instead of ν θr [START_REF] Lempriere | Poisson's ratio in orthotropic materials[END_REF]