O Van Der Sluis
email: olaf.van.der.sluis@philips.com

Y Y Hsu

H M Timmermans

M Gonzalez

P M Hoefnagels

J P M Hoefnagels

Stretching-induced interconnect delamination in stretchable electronic circuits

Keywords: numbers: 6835, 4655, 4650

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Typical examples of stretchable electronics are electronic textiles and biomedical applications, as discussed in e.g. [START_REF] Li | Compliant thin film patterns of stiff materials as platforms for stretchable electronics[END_REF][START_REF] Khang | A stretchable form of single-crystal silicon for high-performance electronics on rubber substrates[END_REF][START_REF] Song | Buckling of a stiff thin film on a compliant substrate in large deformation[END_REF][START_REF] Gonzalez | Design of metal interconnects for stretchable electronic circuits[END_REF] and the recently published review paper by Rogers et al [START_REF] Rogers | Materials and mechanics for stretchable electronics[END_REF] and references therein. In order to provide conformability to the body and a high comfort during wearing, mechanical stretchability is a prerequisite for these electronics applications. In general, stretchable electronics are composed of small rigid semiconductor islands interconnected by thin metal conductor lines, see Figure 1a. The stretchability of the thin interconnects can be realized by depositing serpentine-shaped interconnects [START_REF] Li | Compliant thin film patterns of stiff materials as platforms for stretchable electronics[END_REF] on compliant flat substrates using planar patterning technologies, of which an example is shown in Figure 1b [START_REF] Hsu | In situ observations on deformation behavior and stretching-induced failure of fine pitch stretchable interconnect[END_REF].

The structural reliability of the electronic circuits is one of the major concerns for stretchable electronics applications: (i) cohesive fracture of the interconnects directly leads to failure of the product due to the loss of its electrical functionality. This failure has been addressed by Li et al [START_REF] Li | Compliant thin film patterns of stiff materials as platforms for stretchable electronics[END_REF] and Gonzalez et al [START_REF] Gonzalez | Design of metal interconnects for stretchable electronic circuits[END_REF]. Here, it is conclusively shown that a serpentine shape of the metal conductor permits large deformation of the substrate while at the same time small strains are present in the metal. Hsu et al [START_REF] Hsu | The effect of pitch on deformation behavior and the stretching-induced failure of a polymer-encapsulated stretchable circuit[END_REF] point out that cohesive interconnect fracture is accelerated by interface failure between interconnect and substrate during stretching of fully encapsulated interconnects; (ii) adhesive fracture between the metal lines and the stretchable substrate may lead to shorts, while also the delaminated region may be a source of cohesive fracture (necking of the metal might occur, as discussed in [START_REF] Li | Ductility of thin metal films on polymer substrates modulated by interfacial adhesion[END_REF][START_REF] Lu | Metal films on polymer substrates stretched beyond 50%[END_REF]) and thus electrical failure of the stretchable device.

Recently, horseshoe and zigzag interconnect patterns have been proposed as metal conductor designs [START_REF] Gonzalez | Design of metal interconnects for stretchable electronic circuits[END_REF][START_REF] Hsu | In situ observations on deformation behavior and stretching-induced failure of fine pitch stretchable interconnect[END_REF] involving in-plane patterned copper conductors onto poly(dimethylsiloxane) (PDMS) rubber flat substrates. The horseshoe shape offers high stretchability (strains up to 100% have been achieved without metal failure [START_REF] Hsu | The effect of pitch on deformation behavior and the stretching-induced failure of a polymer-encapsulated stretchable circuit[END_REF]) combined with a rather low I/O density. In contrast, the zigzag shape provides moderate stretchability (strains up to 60%) while a high I/O density can be achieved [START_REF] Hsu | In situ observations on deformation behavior and stretching-induced failure of fine pitch stretchable interconnect[END_REF].

In order to understand the mechanisms leading to cohesive failure of the conducting lines, insight into the key factors that promote delamination is required [START_REF] Lu | Delamination of stiff islands patterned on stretchable substrates[END_REF]. For this purpose, 3D numerical models are developed to describe the delamination process during stretching. The mechanical behaviour of the PDMS rubber substrate is described by a hyperelastic model while for the copper conducting lines a nonlinear elasto-plastic model is employed. The constitutive parameters have been established by performing mechanical characterization experiments. In the finite element model, the interface behaviour is described by cohesive zone elements (see e.g. [START_REF] Ortiz | Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis[END_REF]). The cohesive zone properties have been characterized by means of a 90 ∘ peel test during which actual rubber fibrillation and subsequent fracture is observed. Accordingly, the cohesive zone element as formulated by Van den Bosch et al [START_REF] Van Den Bosch | A cohesive zone model with a large displacement formulation accounting for interfacial fibrilation[END_REF] is employed which is capable of dealing with large deformations occurring at the interface.

Although the adhesive toughness Γ is rather straightforward to determine from peel test experiments [START_REF] Freund | Thin film materials; Stress, defect formation and surface evolution[END_REF], the adhesive strength max is less obvious. In literature, the use of an extended Kalman filter [START_REF] Corigliano | Parameter identification of a time-dependent elastic-damage interface model for the simulation of debonding in composites[END_REF] and variational or inverse analyses [START_REF] Ferracin | On the determination of the cohesive zone properties of an adhesive layer from the analysis of the wedge-peel test[END_REF][START_REF] Zhao | Determination of interface properties between micron-thick metal film and ceramic substrate using peel test[END_REF][START_REF] Gustafson | The influence of adhesive constitutive parameters in cohesive zone finite element models of adhesively bonded joints[END_REF] have been proposed to extract cohesive zone parameters. In [START_REF] Geißler | Peel process simulation of sealed polymeric film computational modelling of experimental results[END_REF], the interface strength is determined by the ratio of peel force and cross-sectional area of the peel arm. Liechti and Wu [START_REF] Liechti | Mixed-mode, time-dependent rubber/metal debonding[END_REF] select the value of the maximum traction by matching experimental and numerical results of the crack-tip opening displacement fields during the initial stages of crack propagation. In this paper, a physically based approach is pursued wherein use is made of the local deformation of the compliant substrate and the copper film at the delamination front. This is achieved by combining finite element results with experimentally measured local deformations, of which first results have been reported in [START_REF] Van Der Sluis | Analysis of the three-dimensional delamination behavior of stretchable electronics applications[END_REF].

The thus obtained material and cohesive zone properties are then applied in the numerical models to calculate the delamination behaviour of three-dimensional horseshoe and zigzag shaped interconnects. To assess the obtained accuracy of the numerical results, the evolution of the interconnect shape during stretching is compared with in-situ SEM and optical measurements.

The paper is organized as follows: in Section 2, the mechanical properties of the metal lines and the PDMS substrate are determined from experiments. Furthermore, cohesive zone properties are established by means of 90 ∘ peel testing. In Section 3, the stretching induced delamination in the horseshoe and zigzag shape interconnect samples are predicted by means of 3D finite element simulations thereby employing the experimentally determined material and cohesive zone properties. The numerically obtained deformation behaviour of the sample designs is compared with experimental in-situ measurements. The paper ends with conclusions and recommendations in Section 4.

Characterization

Peel test experiments

Several testing methods are available to extract thin film interface properties and are well-documented in literature, see e.g. [START_REF] Freund | Thin film materials; Stress, defect formation and surface evolution[END_REF][START_REF] Volinksy | Interfacial toughness measurements for thin films on substrates[END_REF][START_REF] Kolluri | Insitu characterization of interface delamination by a new miniature mixed mode bending setup[END_REF]. In this paper, a 90 ∘ peel test is used to establish adhesion properties. In this specific test, a thin film is peeled-off a substrate at a fixed angle of 90 ∘ from the substrate. The force is measured as a function of the clamp displacement and during peeling the force reaches a more or less constant value, called the peel force. In general, this force is dependent on film thickness, peeling angle, materials and interface properties [START_REF] Thouless | A parametric study of the peel test[END_REF]. At the steady-state delamination process, the following relation holds:

Γ 0 + Γ = (1 -cos ), (1)
where Γ 0 the energy dissipation due to interface fracture (also called the intrinsic interface toughness) and Γ the energy dissipation due to any dissipative mechanism except interface delamination, such as plasticity [START_REF] Zhao | Determination of interface properties between micron-thick metal film and ceramic substrate using peel test[END_REF]. Moreover, is the peel force, the width of the sample, and the peeling angle. When assuming that the film and substrate deform elastically, the substrate stiffness is high compared to the film stiffness, and only interface fracture contributes to dissipation, (1) reduces to

Γ 0 = (1 -cos ). (2)
A more detailed treatise of film peeling can be found in [START_REF] Kim | Elastoplastic analysis of the peel test[END_REF][START_REF] Freund | Thin film materials; Stress, defect formation and surface evolution[END_REF][START_REF] Thouless | A parametric study of the peel test[END_REF].

The peel test samples studied in this paper consist of a substrate material poly(dimethylsiloxane) (PDMS) (Sylgard 186 Ⓡ , Dow Corning) with a copper (electrodeposited, grade TW-YE, Circuit Foil) film. The rubber is cast in a mould onto the copper film. After casting the samples are cured at room temperature. The rubber has a thickness of approximately 1.0 mm while the thickness of the copper film is 18 ± 0.5 m, and consists of an untreated shiny side with low roughness and a treated side with high roughness values ranging from 6 to 8 m, used for improved rubbercopper adhesion; in addition an adhesion promoter is added. The samples are 84 mm in length and 18 mm in width. Prior to testing, the rubber side of the sample is glued on a metal plate that is mounted in the test set-up. For initial clamping of the sample, a part of the copper is peeled-off manually after mounting. To achieve a constant peeling angle, the device is mounted under an angle of 45 ∘ with respect to the tensile device, while the load cell is positioned perpendicular to the sample, see Figure 2a.

Even though several of the aforementioned assumptions underlying the relation between peel force and adhesion energy (2) do not hold for the studied samples, this relation is used to calculate a first estimate for the adhesion energy or work of separation Γ according to Γ = / , when assuming that the peeling angle remains approximately 90 ∘ during testing. From the resulting force-displacement curves, discussed later on in the paper and depicted in Figure 5b, the average work of separation is 1343 ± 51 Jm -2 . Extraction of the interface strength value from interface experiments is however not as straightforward. In literature, several approaches have been proposed to determine the interface strength value max , and are mentioned in Section 1 [START_REF] Corigliano | Parameter identification of a time-dependent elastic-damage interface model for the simulation of debonding in composites[END_REF][START_REF] Ferracin | On the determination of the cohesive zone properties of an adhesive layer from the analysis of the wedge-peel test[END_REF][START_REF] Zhao | Determination of interface properties between micron-thick metal film and ceramic substrate using peel test[END_REF][START_REF] Gustafson | The influence of adhesive constitutive parameters in cohesive zone finite element models of adhesively bonded joints[END_REF][START_REF] Geißler | Peel process simulation of sealed polymeric film computational modelling of experimental results[END_REF][START_REF] Liechti | Mixed-mode, time-dependent rubber/metal debonding[END_REF]. Alternatively, in this paper the local deformation of the substrate and the film is used to extract the interface strength [START_REF] Van Der Sluis | Analysis of the three-dimensional delamination behavior of stretchable electronics applications[END_REF]. Three local geometry parameters are defined: the width of the lifted rubber , the height of the lifted rubber ℎ, and the radius of the copper foil , see Figure 2c. The width is defined as the distance between the points where the lifted rubber crosses the original rubber upper plane. The height is defined as the distance from the rubber upper plane to the top of the rubber peak. For this purpose, a camera is added to the experimental set-up, as shown in Figure 2a. The camera is fixed with respect to the metal plate on which the sample is glued, and is positioned perpendicular to the side of the sample. An example of an image recorded with the camera system of the local geometry near the delamination front during steady-state peeling is shown in Figure 2c. The averaged values over 15 measurements of different samples are: ℎ = 1.31 ± 0.02 mm, = 2.1 ± 0.3 mm, = 1.43 ± 0.05 mm. The relatively large range in the width value can be attributed to uncertainty of the position of the rubber upper plane in the images.

During peeling of the copper film from the rubber substrate, rubber fibrillation is observed at the copper-rubber interface, which is illustrated by the ESEM picture in Figure 3a. In addition, actual fracture of the rubber fibrils occurs during peel testing as evidenced by surface analysis on the copper film that reveals remaining traces of rubber after peeling, see Figure 3b. This rubber fibrillation and subsequent fracture is hypothesized to be the main explanation for the high work of separation value mentioned above. It is remarked that actual fracture toughness testing of the PDMS rubber reveals values of 16 ± 2.3 kJm -2 , obtained by following the test methods proposed in [START_REF] Rivlin | Rupture of rubber. I. Characteristic energy for tearing[END_REF] and [START_REF] Kim | Single specimen test method for determining fracture energy ( ) of highly deformable materials[END_REF].

Determination of bulk mechanical properties

The mechanical behaviour of the silicon elastomeric rubber substrate will be described by a hyperelastic material model, commonly used to describe the nonlinear elastic behaviour of rubber materials at large strains. In hyperelasticity theory, the relation between the stress and strain measures is derived from the strain energy density function ( 1 , 2 , 3 ) with the principal stretch ratios. It can be shown that these models exhibit real elastic behaviour meaning that no energy is dissipated. In a finite element context, no remaining stresses are present after a closed deformation path [START_REF] Weber | Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids[END_REF]. A simple neo-Hookean hyperelastic constitutive model is selected to describe the mechanical behaviour of the PDMS material as more sophisticated models did not reveal any significant improvements in their ability to match the experimental data. The strain energy equation for this model reads

= ( 1 -3), (3)
in which = 0.5 with = is the rubber elastic (shear) modulus, with the chain density, Boltzmann's constant, the absolute temperature, and 1 is the first stretch invariant according to 1 = 2 1 + 2 2 + 2 3 . In the neo-Hookean approach, the material is assumed to be isotropic and to undergo sufficiently small deformation not to reach the finite extensibility of the rubber polymer chains. In that case, the probability distribution of chain extension is Gaussian [START_REF] Ward | An Introduction to the Mechanical Properties of Solid Polymers[END_REF]. Uniaxial tensile experiments and planar extension tests are performed under room conditions (23 ± 2 ∘ C and 50 ± 5% Relative Humidity) on a Zwick 1474 tensile testing device. To check the applicability of the neo-Hookean model, hysteresis loops are performed by loading/unloading cycles. By fitting the Neo-Hookean model on the measured uniaxial and planar extension results, a value of = 0.165 [MPa] is established, see Figure 4a. The obtained accuracy is acceptable.

The mechanical behaviour of the copper foils is described by a hypo-elasto-plastic material model including large deformations, and relies on the linear relation between the objective Jaumann rate of the Cauchy stress tensor

deformation rate tensor ∘ = ( -), (4)
in which the Hookean stiffness tensor = ( + ) + , with , Lamé's constants, and the additive decomposition of the strain rate tensor is applied, = + [START_REF] Nemat-Nasser | On finite deformation elasto-plasticity[END_REF]. The occurrence of plastic yielding depends on a yield function, denoted by ( , ), where is a history parameter, in this case the total equivalent plastic strain. As is common in metal plasticity models, the deviatoric von Mises yield function is applied, expressed as

= -( ), (5)
where ( ) defines the current yield stress and is determined from uniaxial tensile experiments. For more details on finite strain elasto-plasticity, the reader is referred to [START_REF] Simo | Computational Inelasticity[END_REF]. The experiments on the copper foils are performed under room conditions using a Kammrath & Weiss 10 kN micro-loading stage equipped with a 100 N load cell. The displacements are measured on the specimen surface by digital image correlation (DIC). Subsequently, the yield stress function ( ) is obtained from the measurements by converting the measured force to the Cauchy stress and the displacement to logarithmic strain values. The Cauchy stress is calculated by assuming incompressibility which is allowed due to the fact that the elastic part is very small for metals: true = ( / 0 )( / 0 ), with the measured force, 0 the initial cross section of the specimen, and and 0 the actual length and initial length of the specimen, respectively. The logarithmic strain is determined from ln = ln( / 0 ). The experimental and numerical stress-strain curves are depicted in Figure 4b. To prevent negative stiffness values during the numerical simulations, the thus determined yield curve was smoothened by removing the noise from the curve. It is remarked that a simplified isotropic model is assumed even though a columnar grain structure over the thickness of the foil was observed.

Determination of cohesive zone properties

In order to establish the cohesive zone properties Γ and max from the peel test, results of a finite element peel test model will be matched with experimental results. The 5a. To reduce the calculation time, the length of the peel sample is taken to be 15 mm instead of the specimen length of 84 mm, which is still sufficiently large to ensure reaching the steady-state plateau value . At the bottom of the substrate all degrees of freedom are suppressed, which corresponds to the experiment where the sample bottom is glued to the test setup. The peeling of the copper film is simulated by an initially delaminated region of 5 mm and prescribing a displacement in -direction at the upper left corner point of the copper foil: = 13 mm.

As already mentioned, large interface separations accompanied by fibrillation of the PDMS occur during the peel test. Therefore, the delamination process is properly captured by the large displacement cohesive zone elements as formulated by Van den Bosch et al [START_REF] Van Den Bosch | A cohesive zone model with a large displacement formulation accounting for interfacial fibrilation[END_REF]. In this formulation, first Piola-Kirchhoff tractions are employed. Consequently, the required cohesive zone element length corresponds to the initial element length, which can be determined unambiguously. The opening and traction vectors are no longer decomposed with respect to a local basis but resolved globally and thus, no distinction is made between normal and tangential direction. A single exponential relation between the traction and separation is formulated according to

= Γ ( ) exp ( - ) , (6)
in which the traction vector ⃗ = ⃗ and the opening displacement vector ⃗ = ⃗ , where ⃗ is the unit vector along the line between corresponding, opposite points of the interface, Γ the work of separation, and is the separation value at which the maximum traction

max is reached, according to = Γ max exp(1) . ( 7
)
Irreversibility is taken into account by assuming linear elastic unloading to the origin, i.e. elasticity-based damage. In [START_REF] Van Den Bosch | A cohesive zone model with a large displacement formulation accounting for interfacial fibrilation[END_REF], it is noted that in the case of interfacial fibrillation the presence of a traction vector parallel to the opening displacement vector and a single work-of-separation are valid assumptions, however, for other failure mechanisms this might not always be the case. In the peel test model illustrated in Figure 5a, these cohesive zone elements are located at the interface of the copper film and the rubber substrate, except at the position of the manually delaminated initial imperfection in the peel experiment. The geometry is discretized using a sufficient amount of linear quadrilateral finite elements that account for geometrical and material non-linearities, such that a unique, converged solution is obtained.

To determine the cohesive zone parameters, the work of separation Γ is varied between 1200 and 1400 Jm -2 , based on the preliminary value mentioned earlier in this section. The interface strength max is varied between 1.0 and 4.0 MPa. To illustrate the deformation behaviour of the PDMS-copper samples obtained from the numerical peel test model, Figure 6 shows typical deformation stages during loading.

Upon using Γ as the work of separation in the cohesive zone formulation, consequently, all dissipation processes are 'lumped' onto the cohesive zone elements except for the plastic dissipation in the copper film which is captured by the continuum plasticity model from ( 4) and ( 5), and is denoted as Γ . Hence, this lumping approach requires the redefinition of the energy balance (1)

Γ = Γ 0 + Γ -Γ . (8)
By matching the numerical and experimental steady-state peel forces values, the critical work of separation value is obtained as Γ = 1330 Jm -2 with a resulting numerical peel force of = 24.3 N (Figure 5b). Moreover, from the three local geometry parameters, max is obtained as max = 3.0 MPa which results in ℎ = 1.2 mm, = 1.4 mm, and = 1.5 mm. Interestingly, it is found that the effect of the value of Γ on the local geometry is negligible, that is, in the selected range of parameters. The dependency of the peel test results on deformation rate appears to be insignificant. With the thus determined parameters, the critical opening can directly be obtained from ( 7):

= 163 m. In [START_REF] Hoefnagels | Copperrubber interface delamination in stretchable electronics[END_REF], the relation between and observed fibril lengths is discussed. Although the geometry parameters are in overall agreement with the experimental values, the following factors have a pronounced impact on these values: (i) inaccuracies in the experimental values due to the fact that a reference line for the original rubber upper plane has to be constructed to obtain the values (see Figure 2c); (ii) the neo-Hookean model does not fit the experiments exactly (see Figure 4a). In fact, it is known that the validity of the neo-Hookean model is limited to strain values up to 40%; (iii) the elasto-plastic model of the copper is isotropic while a columnar microstructure is present over the thickness of the foil; (iv) all dissipative mechanisms are lumped entirely onto the CZ elements, as formulated in [START_REF] Li | Ductility of thin metal films on polymer substrates modulated by interfacial adhesion[END_REF], thereby reducing the CZ parameters to model parameters instead of intrinsic interface parameters. Although this macroscopic approach is very well suited for the analysis of the interface delamination behaviour at comparable mode angles, as can be observed from the results presented in the next section, detailed micro-mechanical analyses combined with extensive peel testing under different loading angles are required to be able to predict delamination for the whole range of mode angles.

Application to two characteristic three-dimensional structures

In this section, the deformation behaviour of two typical interconnect structures will be evaluated: the zigzag structure, as proposed in [START_REF] Hsu | In situ observations on deformation behavior and stretching-induced failure of fine pitch stretchable interconnect[END_REF], and the horseshoe shape [START_REF] Gonzalez | Design of metal interconnects for stretchable electronic circuits[END_REF][START_REF] Hsu | The effect of pitch on deformation behavior and the stretching-induced failure of a polymer-encapsulated stretchable circuit[END_REF]. As argued previously, the horseshoe shape offers high stretchabilities (strains up to 100% have been achieved without metal failure) combined with a rather low I/O density. In contrast, the zigzag shape provides moderate stretchability (strains up to 60%) while a high I/O density can be achieved [START_REF] Hsu | In situ observations on deformation behavior and stretching-induced failure of fine pitch stretchable interconnect[END_REF].

The zigzag patterned interconnect

To predict the stretching induced delamination behaviour of the zigzag interconnect structure depicted in Figure 7a, numerical simulations are performed using the obtained material and cohesive zone properties from the preceding section. As boundary conditions, in accordance with the experiments reported in [START_REF] Hsu | In situ observations on deformation behavior and stretching-induced failure of fine pitch stretchable interconnect[END_REF], a uniaxial strain of 60% is prescribed along the -axis. The geometry of the sample is defined by the following values: 0 = 60 ∘ , l Cu = 2 mm, w Cu = 0.1 mm, W S = 10 mm, and L S = 8 mm. Furthermore, the thickness of the interconnect is 18 m, the thickness of the PDMS is 0.5 mm, while its inner and outer radii are 0.1 mm. As can be seen in the picture, the zigzag sample contains 4 zigzag interconnects. Of course, the geometry is discretized using a sufficient amount of elements such that mesh convergence is ensured: the model contains 75200 elements of which a detail is given in Figure 7b.

In order to assess the results of the numerical model, implicitly determined by the accuracy of the obtained model parameters, the evolving geometry of the interconnects during sample loading is used. In Figure 8, the simulated deformation and the corresponding in-situ optical microscope images are depicted at different deformation stages. Three geometry parameters have been defined according to Figure 8a: is the angle of the zigzags with respect to the -axis, the length and the height of the zigzag structure. However, optical inaccuracies prevent absolute comparison of these values: indeed, the measured initial value of is 1.91 mm from the optical images while according to the specifications, this value is 2.00 mm (which is the actual value in the numerical model). Alternatively, the relative change of the geometry parameters between subsequent deformation stages is used at 25%, 50% and 60% strain, hereby implicitly assuming that the mentioned optical inaccuracies remain equal during the experiment. In Table 1, the results are given. From these results, it is found that the numerical delamination accurately captures the deformation of the interconnects during stretching as a maximum deviation of 6% is obtained.

Remarkably, the numerically predicted values are larger than the experimental values, except during the strain change from 50% to 60%. Here, the numerical values are smaller than the experimental values. During the elongation from 0 to 50%, the interconnects are deforming mainly in opening mode. During this stage, the zigzags are deforming in a more or less uniform way. It is observed that at around 60% deformation, the interconnect starts to twist in order to accommodate for the increasing deformation. This can be seen in Figure 9a: here, the local deformation and rotation below the crest is clearly visible. This increased deformation explains the higher values of the geometry parameters compared to the numerical results, in which this localized deformation appears not to be present. This can be attributed to the assumptions of the applied elasto-plastic model for the copper, as mentioned earlier.

Fibrillation of the PDMS can be observed in Figure 9a. In the numerical model this process is described by the cohesive zone elements in an averaged sense. To illustrate this, the calculated region of the deforming cohesive zone elements is identified in black in Figure 9b. It appears that the damaged or fibrillation regions are in good agreement from which it can be concluded that the lumping approach results in a meaningful, albeit averaged, description of the fibrillation process.

The horseshoe patterned interconnect

Next, the stretching induced delamination behaviour of the horseshoe interconnect structure depicted in Figure 10a is considered. The material and interface properties and boundary conditions are similar to the zigzag structure, except that in this case, the value of the prescribed uniaxial strain is 100% along the -axis. This is due to the fact that this structure is able to withstand larger strains before final failure occurs. The geometry of the sample is as follows: R Cu = 0.75 mm, 0 = 120 ∘ , w Cu = 0.1 mm, W S = 10 mm, and L S = 10.4 mm. The thickness of the copper is 18 m while the thickness of the PDMS is 0.5 mm. As can be seen in Figure 10a, the horseshoe sample contains 4 entire horseshoe interconnects, while mesh convergence is achieved with 82000 elements of which a detail is depicted in Figure 10b. In Figure 11, several in-situ deformation stages obtained in the SEM are depicted. At 30% strain, only local PDMS deformation occurs around the copper layer, indicated by the arrow. Continuing to 50% strain, the PDMS stretches globally which can be concluded from the change in the surface texture of the sputtered carbon layer of a few nanometers thick. In addition, local fibrillation occurs which is however hardly visible in the picture. At 100% stretch, large-scale fibrillation takes place around the copper lines. Notice the apparent S-shaped delamination regions around the interconnects. Application of different strain rates, ranging from 0.1%sec -1 to 10.0%sec -1 , did not result in any noticeable difference in deformation of the horseshoe structures.

In order to validate the numerical model, the evolving geometry of the interconnects during sample loading is used in this case as well. However, for this geometry, it's not straightforward to define parameters that describe the interconnect geometry in a unique way as this shape evolves in a non-circular way. Although in [START_REF] Hsu | The effect of pitch on deformation behavior and the stretching-induced failure of a polymer-encapsulated stretchable circuit[END_REF], the horseshoe width at two different height values is considered for this purpose, in this paper a discrete geometrical description is generated from the SEM pictures by digital image processing and used instead. Consequently, coordinates of points located on the outer and inner boundaries of the horseshoes are extracted from the pictures given in Figure 11. Of course, the coordinates of these corresponding points from the numerical analyses are readily available from the simulation results. The thus extracted interconnect shapes are depicted in Figure 12. From the fact that the initial experimental and numerical shapes are almost identical (but not shown in the picture), it is concluded that this method is indeed sufficiently accurate even though extraction of the relevant coordinates from the SEM images is a manual procedure which inherently introduces a certain degree of inaccuracy. From the pictures, the stretching in horizontal direction and lateral contraction in vertical direction can be recognized. It can be concluded that a good agreement is obtained between the numerically calculated and experimentally obtained horseshoe shapes, as function of increasing stretch ratio. Noticeably, the shapes start to deviate from 80% straining. After careful inspection of the test setup, it appears that slippage at the clamps occurs which clearly explains the observed 'delay' in deformation of the horseshoes. There are several other reasons that introduce inaccuracies, which are in fact mentioned in Section 2.3, but it is believed that the occurring slippage is indeed the main reason for the observed deviation.

From Figure 11c, PDMS fibrillation is apparent. It is interesting to verify this observed fibrillation area with the numerically determined damaged zone, similar to the zigzag structure in Figure 9. In Figure 13a, the calculated damaged area is illustrated in black. From Figure 11c, it can be concluded that the damage location and shape resemble in a qualitative way. On the other hand, the experimental damage zone appears more pronounced at the sides and not present at the bottom, which is different from the numerical simulations. This can again be explained by the followed lumping approach. Another explanation is perhaps the difference in mode angles between the peel test and the actual application. The dependence of interface properties on mode angle is not considered in the current approach.

-2.7 -1.8 -0.9 0 0.9 1.8 2.7 0 0.5 Finally, to illustrate the dependency of the damaged zone on max , Figure 13b,c show the shape of the damaged zone for values of 1 MPa and 4 MPa. The lowest value results in the largest damage area while the highest strength value severely limits the damage zone in the PDMS. This confirms that the fitted value of 3 MPa yields the most realistic description of the fibrillation zone in an averaged sense.

Discussion

In order to predict the stretching induced interconnect delamination that might occur in stretchable electronic circuits, three-dimensional finite element models were developed.

Mechanical properties of the PDMS substrate and Cu film were determined from tensile testing from which hyperelastic and hypo-elasto-plastic material models were established, respectively. The cohesive zone parameters, used to model the transient interface delamination, were characterized by analyzing results of the 90 ∘ peel test at two distinct scales: the global force-displacement curve was used to identify the work of separation Γ , while local deformation measures at the delamination front were used to extract the strength value max . Consequently, all dissipation mechanisms occurring during peeling of the copper are lumped onto the cohesive zone elements, including experimentally observed rubber fibrillation and subsequent fracture, but excluding the continuum plasticity, as formulated in the expression for the energy balance [START_REF] Li | Ductility of thin metal films on polymer substrates modulated by interfacial adhesion[END_REF].

The thus obtained mechanical and cohesive zone properties were used to predict the stretching induced delamination of two characteristic interconnect structures: the zigzag and horseshoe designs. To validate the obtained results from the delamination models, the evolving interconnect shapes during stretching were used. For both designs, a good agreement was obtained between the experimentally measured and the numerically obtained shapes. Furthermore, the experimentally observed region in which fibrillation occurred during loading of the structure agreed well with the numerically calculated process or damaged zone in an averaged sense due to the lumping approach. The dependency of this region with respect to the value of the interface strength was illustrated.

In [START_REF] Hoefnagels | Copperrubber interface delamination in stretchable electronics[END_REF] a significant mismatch is reported between the lumped critical separation value, , and actually measured fibril lengths through in-situ microscopic measurements which is attributed to the lumping approach. Naturally, the current approach is not a-priori suited for other loading and testing conditions due to the implicit phenomenological nature. In order to identify the parameters leading to rubber fibrillation, detailed multi-scale analyses are required. Only in this case, the cohesive zone parameters correspond to the intrinsic interface properties. For instance, Lin et al [START_REF] Lin | Modeling the failure of an adhesive layer in a peel test[END_REF] study the failure behaviour of an adhesive layer from which it is suggested that cavitation is a prerequisite to the fibrillation process. It appears that cavitation is governed by surface roughness of the opposite material (in our case, the copper film) and large hydrostatic tension in the highly constrained adhesive layer. A possible approach to consider the effect of surface roughness on adhesion properties is given in the works of Yao and Qu [START_REF] Yao | Interfacial versus cohesive failure on polymer-metal interfaces in electronic packaging-effects of interface roughness[END_REF] and Noijen et al [START_REF] Noijen | Numerical prediction of failure paths at a roughened metal/polymer interface[END_REF]. Supplementary to these deterministic approaches, the stochastic nature of the geometry of the roughness profiles, being of major importance on the occurring interface phenomena, could be taken into account by adopting a perturbation method within a stochastic finite element framework, as formulated in [START_REF] Gutierrez | Stochastic finite element methods[END_REF].

Finally, the local rubber geometry at the delamination front shown in Figure 2b suggests that delamination does not occur purely at mode 1 (i.e. opening mode). In fact, a relatively large mode 2 contribution (i.e. shear mode) is present which indicates that the interface is loaded in mixed mode. To account for this, a mode angle dependent work of separation Γ ( ), with the mode angle, may have to be taken into account for more accurate prediction of delamination phenomena in arbitrary 3D stretchable electronics designs, which can be characterized by either varying the peel angle in the current peel test or by using a dedicated mixed-mode setup, as for instance proposed in [START_REF] Kolluri | Insitu characterization of interface delamination by a new miniature mixed mode bending setup[END_REF] and selected references therein.

Figure 2 .

 2 Figure 2. (a,b) The peel test set-up with camera to study the (c) local geometry at the delamination front with the geometry parameters

Figure 3 .

 3 Figure 3. ESEM pictures of (a) fibrillation at the copper-rubber interface during peel testing; (b) remaining PDMS (black areas) at the copper surface (light areas) after peel testing

∘

  and the elastic part of the

Figure 4 .

 4 Figure 4. Measured and fitted mechanical response of (a) PDMS and (b) copper

Figure 5 .

 5 Figure 5. (a) Illustration of the peel test model (not to scale); (b) experimental and numerical force-displacement curves

Figure 6 .

 6 Figure 6. Illustration of subsequent deformation stages during peeling at a vertical displacement of (a) 2.5 mm; (b) 5 mm; and (c) 10 mm for Γ = 1330 Jm -2 and max = 3.0 MPa. The rubber deformation at the delamination front is apparent (for visualization purposes the cohesive zone elements and the element edges are made invisible).

Figure 7 .

 7 Figure 7. Illustration of the patterned zigzag interconnect structure on top of the PDMS substrate: (a) entire sample; (b) mesh details of a part of the numerical model

Figure 8 .

 8 Figure 8. Deformation of the zigzag structure at a loading of 25%, 50% and 60%: (a,b,c): numerical results in which the mesh is removed for visualisation purposes; (d,e,f): micrographs from experiments. The loading is in horizontal direction.

Figure 9 .Figure 10 .

 910 Figure 9. (a) SEM micrograph after rotating and tilting the sample at 60% stretching in which the fibrillation is clearly visible; (b) the simulated result in which the black area depicts the damaged region as determined by the cohesive zone elements (the mesh is removed from the picture).

Figure 11 .

 11 Figure 11. In-situ SEM images of horseshoe interconnects while stretching to (a) 30%, (b) 50%, (c) 100%

Figure 12 .Figure 13 .

 1213 Figure 12. Deformation of the horseshoe structure at a loading of (a) 30%, (b) 50%, (c) 80%, and (d) 100%. The × symbols correspond to the interconnect outline from the SEM results while the continuous lines are results from the numerical simulations. The loading is in horizontal direction.

Table 1 .

 1 The

	0 → 25% strain 25 → 50% strain 50 → 60% strain
	exp	num	exp	num	exp	num
	20%	23%	30%	33%	31%	25%
	25%	28%	18%	20%	8 %	6 %
	10%	12%	19%	21%	18%	15%

relative changes of the geometrical parameters, defined in Figure

8a

, between different loading stages from experimental measurements ('exp') and numerical simulations ('num')

John Wiley & Sons, Ltd, Chichester, 2004.

Acknowledgments

The authors would like to thank Edwin van der Zanden and Marc van Maris (Eindhoven University of Technology), Sergei Shulepov and Ruud van Mullekom (Philips Applied Technologies) and Fabrice Axisa and Jan Vanfleteren (IMEC). We thank the European Commission for partial funding of this work under project Stella (IST-028026).


Valentin Leroy

Anatoliy Strybulevych

John H Page

Martin G Scanlon

Influence of positional correlations on the propagation of waves in a complex medium with polydisperse resonant scatterers

Keywords: numbers: 43, 35, +d, 43, 20, +g

We present experimental results on a model system for studying wave propagation in a complex medium exhibiting low frequency resonances. These experiments enable us to investigate a fundamental question that is relevant for many materials, such as metamaterials, where low-frequency scattering resonances strongly influence the effective medium properties. This question concerns the effect of correlations in the positions of the scatterers on the coupling between their resonances, and hence on wave transport through the medium. To examine this question experimentally, we measure the effective medium wave number of acoustic waves in a sample made of bubbles embedded in an elastic matrix over a frequency range that includes the resonance frequency of the bubbles. The effective medium is highly dispersive, showing peaks in the attenuation and the phase velocity as functions of the frequency, which cannot be accurately described using the Independent Scattering Approximation (ISA). This discrepancy may be explained by the effects of the positional correlations of the scatterers, which we show to be dependent on the size of the scatterers. We propose a self-consistent approach for taking this "polydisperse correlation" into account and show that our model better describes the experimental results than the ISA.

I. INTRODUCTION

Propagation of waves in complex media continues to be a very active subject of research. Of particular interest are complex media containing scatterers with resonances at low frequencies. Indeed, when the wavelength is large compared to the typical size of the inhomogeneities, one can use an effective medium approach, i.e., consider that the wave propagates as it would do in a homogeneous medium with an effective wavevector whose value is related to the composition and the structure of the material. In most cases, the effective medium has properties that are close to the average of that of its components. But when there are resonant scatterers, their contribution to the wave field can be strong enough to significantly alter the propagation. In some cases, the effective media show intriguing properties, such as a negative refractive index, which are not encountered in nature; such materials are called metamaterials [1][2][3][4].

One of the key issues for wave propagation in strongly scattering materials is the effect of coupling between the scatterers near resonance, an effect that makes the traditional independent scattering approximation (ISA) inadequate. Indeed, when the resonances are strong and/or the concentration of scatterers high, positional corre- * Electronic address: valentin.leroy@univ-paris-diderot.fr lations and multiple scattering loops may have a nonnegligible contribution. To gain a better understanding of the underlying physics, experiments on a wellcontrolled system with strong resonances are needed. One such system can be found in acoustics, where bubbly liquids are often regarded as model systems for the study of wave propagation with low frequency resonances. Indeed, air bubbles in a liquid have a particularly low resonance, known as the Minnaert resonance, with an angular frequency given by

ω M = v 1 3ρ 1 /ρ 0 R ,
where R is the bubble radius, v 1 is the velocity of longitudinal waves in air, and ρ 1 and ρ 0 are the mass densities of the air and the liquid, respectively. As the ratios of the densities and the velocities are small, it is easy to see that this resonance is at very low frequency:

ω M R v 0 = 3ρ 1 /ρ 0 v 1 v 0 ≪ 1.
Here, v 0 is the longitudinal sound velocity in the liquid. This equation shows that, at resonance, the wavelength is large compared to the bubble size. This is important because, then, at the resonance frequency, not only is the response of the bubble strong, but several bubbles (and potentially many) can be driven in phase, thus yielding a collective constructive contribution to the pressure field. This explains why a minute quantity of gas bubbles can dramatically change the effective acoustic properties of a liquid.

If the bubble is no longer in a liquid, but in an elastic medium with shear velocity u 0 , its resonance frequency becomes [5,6]:

ω M = 3v 2 1 ρ 1 /ρ 0 + 4u 2 0 R , (1)
so the same condition of small radius with respect to wavelength persists as long as the elastic medium is a soft solid in which u 0 ≪ v 0 . As a result, elastic bubbly media are also good candidates as model systems for investigating wave propagation in the presence of low frequency resonances. In addition, they offer two practical advantages compared to bubbly liquids: (1) bubbles do not move, which means that a precise knowledge of their positions can be obtained; (2) the frequency of the resonance can be shifted by tuning the shear velocity in the elastic matrix. Bubbly soft media with negligible shear velocity have already been studied [7][8][9], and the ISA was found to reliably predict the experimental results, at least for concentrations of bubbles up to 1%. The present article focuses on a bubbly polydimethylsiloxane (PDMS) sample, in which the shear modulus is expected to significantly change the bubble's response to the incident ultrasonic field and the propagation of waves in the medium. Our aim is to see if the ISA is still relevant in this case. Beyond the fundamental aspect of the study, and the potential applications for design of metamaterials, it should be noted that a better understanding of the acoustic properties of bubbly elastic media is interesting per se, since many practical applications, in the medical or the food science fields for instance, involve the presence of bubbles in a soft elastic matrix. The present paper is organized as follows. The next section gives a brief overview of the ISA, with a focus on its application to acoustic propagation in bubbly media. Then we describe how our bubbly samples were made and characterized. Section IV reports the experimental measurements of velocity and attenuation in one typical sample, and compares the experimental results with the ISA. Two main limitations of the model are also discussed: the assumption of sample homogeneity on length scales larger than the scatterer size (e.g., gradients in the concentration of the scatterers) and the neglect of correlations in the positions of the scatterers. Methods for overcoming these limitations are proposed.

II. THE INDEPENDENT SCATTERING APPROXIMATION

Based on a diagrammatic representation of scattering events [START_REF] Sheng | Introduction to Wave Scattering, Localization and Mesoscopic Phenomena[END_REF], the ISA predicts the following effective wavenumber k for a medium with n scatterers per unit volume defined by their scattering function f s :

k 2 = k 2 0 + 4πnf s , (2)
where k 0 is the wavenumber in the pure medium (matrix). The ISA is actually equivalent to the model developed, with another approach, by Foldy in 1945 [START_REF] Foldy | [END_REF].

In the following, we will use equally "ISA" or "Foldy's model" when refering to equation (2). In eq. ( 2), the effect of the scatterers on propagation in the pure medium is described by a term proportional to nf s , i.e., a simple addition of the individual scattering events, which are assumed to be independent. Models have been proposed for incorporating either the correlations between the scatterings [START_REF] Keller | Proc. Symp. Appl. Math[END_REF] or the loops, i.e., the events in which the same scatterer or spot inside the sample is visited more than once by the wave [START_REF] Henyey | [END_REF]14]. We discuss only the ISA in this section, since it is the simplest model and it gives reasonable predictions of wave transport in many cases. We will introduce later, in section IV C and the appendix, extensions to the ISA for dealing with the positional correlations of the scatterers.

Equation ( 2) considers a monodisperse assembly of scatterers. In practical applications, one is often dealing with polydisperse media, in which case the ISA predicts that

k 2 = k 2 0 + 4πn(R)dRf s (R), (3)
where n(R)dR is the number of bubbles per unit volume whose radius is between R and R + dR, and f s (R) is the scattering function for a scatterer of radius R.

Let us apply the ISA to a bubbly medium. When the scatterer is a bubble, its isotropic scattering function is given by

f s (ω, R) = R (ω M /ω) 2 -1 + iδ , ( 4
)
where δ is the damping constant due to thermal, viscous and radiative losses [15] and ω M the resonance angular frequency given by Eq. ( 1), which can be written

ω 2 M = 3κP 1 + 4µ ′ ρ 0 R 2 , (5)
with κ the polytropic index for the transformations undergone by the gas [15] (κ = 1.4 for air), P 1 the pressure of the gas in the bubbles, ρ 0 the mass density of the matrix, and µ ′ the real part of the shear modulus of the matrix. Surface tension effects are negligible for the bubble size considered in our experiments. Mode conversion is not considered in this model since it has been shown to be negligible for this system on account of the large difference between the shear modulus of the matrix and the longitudinal one (u 0 ≪ v 0 ) [16]. Hence the only significant contribution of the nonzero shear modulus is its effect on the resonance frequency of the bubble, which is increased by the stiffness of the matrix. Another approximation we make is to reduce the scattering function to its isotropic component, neglecting the dipolar and higher contributions. This approximation is generally good when the wavelength is large compared to the size of the scatterers, and it will be shown to be accurate for describing our experiments.

In figure 1, we give two examples of the effective attenuation and velocity as functions of the frequency, as calculated by Foldy's model. For the sake of generality, we consider log-normal distributions of bubble size:

n(R) = n 0 √ 2πǫR exp - (lnR/R 0 ) 2 2ǫ 2 , (6)
where R 0 is the median radius, ǫ the width of the distribution, and n 0 the bubble concentration. Note that for distributions that are not too broad (ǫ < 0.2), log-normal distributions are close to normal distributions: then R 0 and ǫ can be viewed as the mean radius and coefficient of variation (i.e., standard deviation divided by the mean), respectively. The latter quantity, which we will refer to as the polydispersity factor, measures the range of bubble sizes in the distribution. The bubble concentration n 0 is related to the gas volume fraction Φ by

n 0 = 3Φ 4πR 3 0 exp(9ǫ 2 /2) . ( 7
)
The bubble size distribution in Fig. 1 is described by R 0 = 150 µm, ǫ = 25 % and Φ = 2 %. Two host media are considered: water (dashed lines), and an elastic matrix with µ ′ = 0.7 MPa, µ ′′ = 0.4 MPa (solid lines).

For both media, very dispersive behavior is predicted: the attenuation and velocity show peaks, which occur at higher frequencies in the elastic medium. At low frequencies (ω < ω M ), sound propagation in the bubbly medium is weakly attenuated, and its very low phase velocity is accurately predicted by Wood's relation [START_REF] Wood | A Textbook of sound[END_REF], modified to account for the elastic effects. For 5×10 -4 < Φ < 10 -1 , this relation can be approximated (within 5 %) as

v lf ≃ κP 1 + 4µ ′ /3 ρ 0 Φ , (8)
where P 1 is the pressure of the gas in the bubbles and ρ 0 the density of the matrix. Eq. ( 8) gives, in our example, v lf = 122 and 248 m/s in water and the elastic matrix, respectively.

Near the Minnaert frequencies of the bubbles, the acoustic waves are highly attenuated by the bubbly medium: the attenuation reaches a peak α max at frequency f max , corresponding to the resonance frequency for a monodisperse distribution. In water, the peak is very sharp and associated with the bubble resonances (f max = 22 kHz as expected by Eq. ( 5) for a 150-µmradius bubble). Eqs. ( 2) and (4) give the physical origin of this sharp peak: at resonance, the scattering function f s is large and imaginary, which leads to a large imaginary part of the effective wave number k. This sharp resonant effect is followed by a broader regime (roughly 60-400 kHz in our example) in which the attenuation is still large. This is explained by the negative response of

f max v hf v lf a max a hf
FIG. 1: (Color online) Predictions of the ISA for the attenuation α (top) and phase velocity v (bottom) as functions of frequency for bubbles in water (dashed lines) and in an elastic matrix (solid lines) with µ ′ = 0.7 MPa, µ ′′ = 0.4 MPa; the bubble concentration is 2 % and the mean radius is taken to be 150 µm with a 25 % polydispersity. The characteristic features (fmax, αmax, α hf , v lf , and v hf ) are indicated.

the bubbles. Indeed, when driven at a frequency higher than its resonance frequency, a bubble has a scattering function which reduces approximately to f s = -R when δ ≪ 1(see Eq. ( 4)). Even though the bubble is not resonating, its contribution is far from negligible. Indeed, one can rewrite Eq. (2) as

k 2 = k 2 0 1 - Rλ 2 πd 3 , ( 9
)
where d is the typical distance between bubbles. As long as the wavelength is large compared to this distance (λ/d ≫ 1), the square of the effective wavenumber is negative: waves are evanescent. In our example of 2 % bubbly water with 150 µm-radius bubbles, the typical distance between bubbles is 0.9 mm, which means that Rλ 2 /(πd 3 ) > 1 for f < 380 kHz.

In our example of an elastic medium, the high viscosity (µ ′′ of 0.4 MPa corresponds, at 0.1 MHz, to a viscosity 640 times larger than that of water) tends to damp the oscillations of the bubbles at resonance, so the attenuation peak is lower and less sharp (see figure 1). As a paradoxical consequence, the more viscous the medium in which bubbles are embeded, the lower the maximum of attenuation of sound. Note that the sharp peak of attenuation also dissapears when the bubble size distribution is too broad. Interestingly, even if the resonant effects are too weak to give rise to significant resonant attenuation, they still lead to an evanescent regime. In this case, the position of the attenuation peak does not correspond to the resonance of the average size bubbles, but rather to a frequency that is just above the resonance frequency of the small bubbles of the distribution, as it corresponds to all the bubbles being driven at a frequency higher than their resonance frequencies. For instance, in our example f max = 90 kHz for the bubbly elastic medium, which corresponds to a radius of 100 µm according to Eq. (5).

When the resonances are smoothed out, an approximate formula can be established for the maximum of the attenuation:

α max = 2 √ 3Φ R 0 exp(2ǫ 2 ) , (10)
which has the benefit of depending only on the bubble concentration and size distribution (not on either the damping constant of the bubbles or the shear modulus of the matrix). At higher frequencies, the evanescent regime disappears: the attenuation and velocity reach constant values, which can be approximated by:

α hf = 3Φ R 0 exp(5ǫ 2 /2) , (11)
v hf = v 0 . (12)
Interestingly, if Foldy's model applies, one can determine the 5 parameters Φ, R 0 , ǫ, µ ′ , and v 0 from the measurement of the 5 quantities v lf , f max , α max , α hf , and v hf in a bubbly medium.

III. MATERIALS AND METHODS

A. Sample preparation

Samples were prepared by mixing 100 g of RTV 615 (GE Silicones, 90 g monomers and 10 g hardener) in a 600 mL beaker with a Sorvall Omni Mixer (Ivan Sorvall, Inc., Norwalk, CT; nominal operating speed 16,000 rpm) for 1 minute to incorporate bubbles. Then, after waiting for an interval of 10-20 minutes for the biggest bubbles to rise and disappear, the bubbly liquid was poured into a cell (thickness 2.6 mm, diameter 11 cm), which was rotated at 5 rpm to prevent creaming (bubbles rising to the top of the container). At room temperature, a solid sample was obtained in 1 day.

B. Sample characterization

Testing of Foldy's model requires that some of the parameters of the sample are known: the bubble size distribution (n(R)), the matrix rheology (µ) and the velocity of sound in the matrix (v 0 ).

a. Optical characterization was made by taking pictures at different positions in the sample. The depth of field was larger than the thickness of the sample, which means that the total volume of the image was known (16 × 11 × 2.6 mm 3 ). A typical image is shown in figure 2; the contrast was good enough for an automatic size analysis to be performed (ImageJ). The radius measurement was possible with a one pixel accuracy, i.e., 5 µm. A total of 839 bubbles was analyzed. To avoid biased measurements of big bubbles, overlaid bubbles were excluded from the size analysis. However, the total number of bubbles was determined (1098) for a correct estimation of the bubble volume fraction, assuming the same size distribution for the overlaid bubbles. The radius distribution was centred around 165 ± 5 µm with a polydispersity of 21 ± 1 % and a volume fraction of 2 ± 0.2 % (see histogram in figure 3).

b. Rheological characterization was necessary to measure the shear modulus over the frequency range used for the ultrasonic experiments (30-800 kHz). We used a shear wave reflection technique (see [START_REF] Longin | [END_REF] or [19] for details) to measure the complex shear modulus µ = µ ′ +iµ ′′ from 300 to 500 kHz. For lower frequencies (around 30 kHz) an estimation of µ was possible by measuring the acoustic response of a single bubble (see [START_REF] Strybulevych | Proceedings of Symposium on Ultrasonic Electronics[END_REF] for a description of the technique). Over the 30-500 kHz frequency range, our measurements are consistent with the following behaviour of the real (µ ′ ) and imaginary (µ ′′ ) parts of the shear modulus:

µ ′ = 0.6 + 0.7 × f, (13a) µ ′′ = 0.2 + 1.8 × f, (13b)
where the moduli are in MPa and the frequency f in MHz. Equations ( 13) only give the order of magnitude of the shear modulus for a given RTV615 sample. In- deed, the precise protocol used for preparation, such as the temperature and time of the hardening process for instance, may have a non-negligible impact on the modulus [START_REF] Schneider | [END_REF]. As a consequence, deviations from the predictions of Eqs. [START_REF] Henyey | [END_REF] can be observed for a given sample.

The stationary limit (f = 0 MHz) of our estimate of µ ′ is in good agreement with the value reported in the literature [START_REF] Schneider | [END_REF]. c. Ultrasonic characterization was also performed by standard techniques [8] to measure the velocity and attenuation of sound in the pure PDMS. The longitudinal velocity was found to be constant over the frequency range investigated (30 kHz -5 MHz): v 0 = 1.02 mm/µs, in accordance with previously reported values [22]. For the attenuation, the following law was found: α 0 = 0.023f 1.54 mm -1 (with frequency f in MHz). d. X-ray characterization. To investigate the spatial distribution of the bubbles in one of the samples, x-ray tomography images were acquired, using an Xradia Mi-croXCT 3D x-ray transmission microscope (18 µm/px). After the ultrasonic measurements, the sample was cut in pieces small enough to fit in the apparatus (7.5 × 10 × 2.6 mm 3 ), and 21 pieces were imaged, for a total of 3700 bubbles. A software program was written to enable the reconstruction of the 3D structure of each piece: for every bubble i, its radius R i and position [x i , y i , z i ] were determined with 1 pixel accuracy. The size distribution was determined and successfully compared to the results of the optical characterization: we found a mean radius of 160 ± 18 µm with a polydispersity of 25 ± 1 % and a volume fraction of 2.2 ± 0.3 % (see figure 3 andtable I). Moreover, the 3D data gave insight into the homogeneity of the sample, as well as on the correlations in the bubble positions, as described in sections IV B and IV C.

Φ (%) R (µm) √ R 2 -R 2 R (%) optics 2.0 ± 0.2 165 ± 5 21 ± 1 x-ray 2.2 ± 0.3 160 ± 18 25 ± 1 C. Ultrasonic measurements
The acoustic properties of the samples were measured with the following set-up. In a large tank (60 × 60 × 120 cm 3 ) filled with reverse osmosis water, a piezoelectric transducer generated a pulse that propagated through water, traversed the sample and was detected by a hydrophone. Because the attenuation in the sample was large, and because the divergence of the beam was not negligible (especially at low frequencies), the use of a screen was essential to reduce spurious signals. The aperture of the screen (6 cm) was larger than the wavelength of the pulse (3.75 cm at 40 kHz) to limit diffraction effects, but smaller than the diameter of the sample (11 cm).

Gaussian pulses, with central frequencies ranging from 40 to 600 kHz, were generated by an Arbitrary Wave Generator, and three different transducers, having central frequencies of 100, 250 and 500 kHz, were used to cover the range of frequencies. The pulses were recorded, with a hydrophone, in two different cases: when the sample was mounted on the screen (s 1 (t)), and when the sample was absent (s 2 (t)). The signals were averaged over 100 acquisitions when the attenuation was low (for reference measurements, for example), and up to 5000 acquisitions for highly attenuated signals. Signals s 1 (t) and s 2 (t) were then truncated to eliminate spurious echoes, and Fourier transformed into S 1 (ω) and S 2 (ω) respectively. Then, T (ω), the ratio of the transmission with and without the sample in the path of the acoustic beam at a given angular frequency, was calculated.

For an incoming plane monochromatic wave exp(ikxiωt), T (ω) is given by

T (ω) = 4ZZ w (Z w + Z) 2 × e i(k-kw)d 1 -Z-Zw Z+Zw e ikd 2 , (14)
where k and k w are the wave numbers in the sample and in the water respectively, Z and Z w the acoustic impedances, and d the thickness of the sample. Because the impedance Z depends on k (Z = ρω/k, with ρ the mass density of the sample), equation ( 14) is not directly invertible to extract k: an iteration method was used [8].

From k, the phase velocity v = ω/Re(k) and the attenuation α = 2Im(k) can be calculated. Error-bars on the measurements were evaluated by taking into account the following sources of uncertainty: thickness of the sample, noise in the Fourier transforms, truncation of the time pulse.

An important issue in measuring the acoustic properties of a medium is the statistical relevance of results for a single sample. Theories predict the average properties, i.e., what one can measure when averaging over different samples having the same general parameters. To obtain insight on this issue, we measured the transmitted signal not only in the centre of the acoustic beam, but also on a 7 × 7 grid with a 5 mm step. From this scanning, a "typically averaged signal" could be obtained (with correction for the geometrical differences), and this was compared to the centre pulse.

IV. EXPERIMENTAL RESULTS AND DISCUSSION

Several bubbly PDMS samples were tested, with different void fractions (0.5-3.5 %) and polydispersities (20-50 %), but with the same average radius (around 150 µm). In this article, we focus on representative results obtained with one selected sample, which was characterized by x-ray tomography. Figure 4 shows the measured attenuation and phase velocity in the sample as functions of the frequency. The general aspect of the curves follows what one would expect for a bubbly medium (see figure 1): both the attenuation and the velocity exhibit peaks in frequency with very large values at their maxima (3.5 mm -1 and 4.5 mm/µs, respectively), the velocity is low (0.26 mm/µs) at low frequencies, and close to the value in the pure PDMS at higher frequencies. Note that two cases are considered: when only one pulse is analysed (open circles) and when the average is made over the 49 positions of the hydrophone (solid circles). Except for the peak in velocity, both analyses give similar results. It indicates that substantial information can be obtained from a single pulse measurement, which is interesting for practical purposes. Note that measurements are difficult when the velocity is high because small time shifts have to be extracted from the phase shift which is dominated by the complex impedance rather than the transit time.

The five typical quantities introduced in section II can be measured and compared to equations (5), and ( 8)-( 12), calculated with the parameters found in section III B. As shown in table II, the agreement is rather satisfying: f max and α hf are within 10% and the other quantities within 20%. Note that for v lf and v hf the measured values may not correspond exactly to the asymptotic values because of the limited range of frequencies in the measurements. a Calculated with Eq. ( 5) by taking the radius of the small bubbles of the distribution (i.e., 100 µm, see figure 3). This is appropriate since for polydisperse distributions of sizes, the peak in attenuation corresponds to the resonances of the smaller bubbles, as noted in section II.

Parameters for the model were taken as measured in section III B, the histogram measured in the x-ray experiment being chosen for n(R). In the low and high limits of the frequency range, the comparison is good. However, the attenuation and velocity peaks are not correctly predicted: their frequencies and amplitudes are both smaller than the experimental data. Note that this discrepancy was observed for all the samples we investigated. The comparison is slightly better when the dipole terms in the bubbles' responses are considered (solid lines) but the discrepancy in the peaks' positions and magnitude is still large. Higher order scattering has negligible contributions. Solid circles correspond to experimental data (same as in Fig. 4). Comparison with Foldy's model is shown for monopolar scattering of the bubbles (dotted lines) and monopolar + dipolar scattering (solid lines). Note that there are no adjustable parameters in this comparison of the theory with experiment. Results of fitting the model to the data by allowing either the size distribution (dashed lines) or the shear modulus (triangles) to vary are also shown. In the inset, the size distribution determined by the x-ray experiment (solid line) is compared to the one inferred from the fit to the ultrasonic data (dashed lines).

To check whether the discrepancy could be explained by the uncertainty in the measured size distribution, we determined the n(R) function that was needed to achieve a better fit. As shown in the inset of Fig. 5, this leads to a totally unrealistic distribution, with smaller bubbles (mean radius 120 µm), a narrower distribution (10 % polydispersity), and a smaller void fraction (1.1 %). One can understand why such a size distribution is obtained by looking at equation (2). If the peaks are at too low frequencies and with too low amplitudes, it means that the nf s term is not high enough compared to k 2 0 , i.e., that the bubble contribution is too weak. A way of fixing this is to consider smaller bubbles, in order to shift the resonance to higher frequencies, and to increase the number of bubbles per unit volume to magnify their contribution. This explains the difference between the fitted and the actual size distribution in figure 5 inset. An important practical application is that if one uses Foldy's model for fitting the experimental data, an incorrect bubble size distribution will be obtained. Note that the same phenomenon was observed for bubbles in bread dough, with the same trend of Foldy's model predicting bubbles smaller than those determined from x-ray measurements [START_REF] Leroy | Bubbles in Food 2: Novelty, Health and Luxury[END_REF].

Another possible source of error was the shear modulus, whose measurement was not very precise. For each frequency, we determined the complex shear modulus that gave the best fit to the experimental data. As shown by the triangles in Fig. 5, a good fit was possible for the attenuation. For the velocity, the agreement was good for lower frequencies, but poor at higher frequencies. Moreover, as shown in figure 6, the fitting leads to unrealistic values for µ ′ and µ ′′ . Again, one can qualitatively understand these values: the real part of the fitted shear modulus is high in order to shift the peaks to higher frequencies, and the imaginary part is low (even zero) so that the oscillations of the bubbles are not damped by viscosity, and the scattering function f s is large.

These results suggest that Foldy's model is inadequate for describing wave propagation in the PDMS sample. In the following we examine two possible explanations, and propose a model to account for these discrepancies. The lines correspond to the results for experimental µ ′ (solid line) and µ ′′ (dashed line), as given in equations [START_REF] Henyey | [END_REF]. Open triangles show the fitted values for µ ′ , solid ones for µ ′′ .

B. Homogeneity of the sample

An important hypothesis made in all the effective medium models consists in assuming that the sample is homogeneous, i.e., any part of the sample can be considered as equivalent to the others. The x-ray inspection of the sample reveals a homogeneous structure in the transverse directions, but a non-homogeneous one along the thickness. Probably because of the upward buoyancy force exerted on the bubbles during the filling of the cell, the concentration of bubbles is varying, as shown in Fig. 7: the bottom part of the sample is free of bubbles, whereas the top part has a higher concentration. The bubble size distribution is also slightly different from one layer to the other.

We investigated the influence of this heterogeneity of the sample on the acoustic propagation. Each layer i (1 < i < 50) of thickness d = 52 µm was considered as an effective medium following Foldy's model with a concentration n ′ /d and a bubble size distribution as measured by the x-ray tomography. Following Brekhovskikh [START_REF] Brekhovskikh | Waves in layered media[END_REF], we computed the transmission through this multi-layer system, and compared it to the experimental transmission. Figure 8 shows the result: accounting for the inhomogeneity of the bubble concentration with the 50 layer model does not change the prediction much, implying that the heterogeneity of the bubble concentration plays a negligible role in acoustic wave propagation through the medium.

C. Positional correlations

Another approximation in Foldy's model is the assumption that the positions of the scatterers are uncorrelated. Recent experiments in a 2D system of steel rods in water have shown that Keller's approach results in a very good correction to the model when correlations needed to be included [START_REF] Derode | [END_REF]. In Keller's approach, the effective wave vector is related to the correlation function by:

k 2 = k 2 0 + 4πnf s -(4πnf s ) 2 (1 -g(r))
sin kr k e ik0r dr (15) where the correlations are taken into account by a radial function g(r) such that the local concentration of scat- 5, the same discrepancy between Foldy's prediction for a homogenous sample (solid lines) and the experimental data (circles) is obtained. The prediction for a layered medium (dashed lines) does not bring better agreement.

terers at a distance r from a scatterer is n(r) = n 0 g(r), where n 0 is the average concentration. For point-like scatterers, g = 1 and Foldy's approximation is recovered. For hard spheres, g = 0 for r < 2R and 1 at larger distances. In general, g can take more complicated forms, on the correlation mechanisms involved in the medium.

From the x-ray tomography data, we were able to estimate the function g for our sample. The general procedure for measuring such a function involves the following steps. First, one considers a sphere of radius r N around one bubble and counts the number of bubbles in each spherical layer between r i and r i+1 . Then, averaging over different central bubbles gives a statistical estimate of g(r). However, a limit of this technique is that the bubbles chosen as centres of the spheres cannot be at less than r N from an edge of the sample. Thus, if one wants to investigate long-range correlations by taking high values of r N , the statistics are poor because the number of bubbles that are far enough from the borders is small. For our thin sample, this limitation was severe. To circumvent this problem, we used 1/8 of each sphere, with an orientation chosen so that its volume did not cross any boundaries of the sample. Figure 9 shows the result of these measurements up to a maximum distance of 0.9 mm. It appears that the g function is very similar to the predictions of the hard-sphere approximation (solid line). Note that the oscillations in the g function that are predicted by the frequently used Perkus-Yevick approximation [START_REF] Derode | [END_REF] are not expected to be significant here because the concentration of scatterers is low. Numerical simulations with random positions of scatterers in a box confirm that a simple hard-spheres law is a good approximation for a 2 % volume fraction.

The predictions of Keller's model for the attenuation and velocity in the sample are shown in Fig. 10 for three different g functions. When the hard sphere approximation is used (i.e., the solid curve of Fig. 9 is taken), one obtains the dashed lines, which only give a slightly better agreement than the ISA. From the measured correlations, we also consider the high limit (i.e., the experimental points plus one standard deviation) and the low limit (i.e., the experimental points minus one standard deviation) for the g function. They give different corrections: while the former gives a worse agreement, the latter improves it. One can understand this difference qualitatively by re-writing equation (15) as

k 2 = k 2 0 1 + 4πnf s k 2 0 1 -4πnf s (1 -g)rdr , (16)
where the reasonable approximations kr ≪ 1 and k 0 r ≪ 1 have been made. It is important to note that, in the frequency range over which the discrepancy is high, the bubbles are responding in phase opposition to the pressure field, i.e. Re(f s ) < 0. So, according to equation ( 16), the correlations increase the effect of the bubbles on k if g < 1, whereas they decrease it if g > 1. Note that a simple criterion can be deduced from Eq. ( 16) to estimate 15)) is plotted for the hard-sphere correlation (dashed lines), for the high estimation of g measured from the x-ray data (dotted lines), and for the low estimation of g (thick solid lines).

the importance of the corrections due to correlation. If g is roughly described as a step function going from 0 to 1 at distance r c , the term with the integral in Eq. ( 16), which we will denote as β, reduces to

β = 3 2 Φ f s R r c R 2 . (17)
Thus, the magnitude of the correction depends on the ratio of r c to the radius of the bubbles R. Long-range correlations are thus expected to have a stronger effect. It appears that Keller's approach, with a reasonable estimate of the g function, introduces a correction that goes in the right direction for better agreement. But the agreement is still unsatisfactory. Moreover, Keller's expression applies to a monodisperse medium, whereas our bubbly PDMS is a polydisperse sample. So far, we have assumed that the polydispersity could be taken into account by changing nf s into n(R)dRf s (R) in the equations. However, as we show in appendix A, a selfconsistent approach provides a more general expression (see equation (A7)), in which the correlations depend on the size of the scatterer: g(r) becomes g(r, R). To exam- ine this polydispersity effect, we considered 3 different populations of "central" bubbles: 0-140 µm, 140-180 µm, and 180-600 µm. For each population, we were able to average over approximately 1000 bubbles, ensuring good statistics. Figure 11 shows the g functions obtained for these 3 populations. For the small bubbles, there is a clear maximum of the g function. This is an indication of clusters. A close inspection of the picture in figure 2 indeed gives the impression that the small bubbles are clustered. For the average size bubbles, the g function is similar to the hard-sphere approximation. For the big bubbles, we find evidence of depletion: there is a lower probability of finding a bubble close to a big one. Again, note that the photo in Fig. 2 shows isolated big bubbles.

Figure 12 shows the predictions of our self-consistent approach (SCA) when 3 different possibilities are considered for g(r, R): hard sphere approximation (solid lines in Fig. 11), and the measured values of g including the low and the high limits (bottom and top of the error bars in Fig. 11). It appears that the effects observed in Fig. 10 are amplified when the polydisperse correlation is taken into account. Interestingly, even for the hard sphere approximation, slightly better agreement is found (see table III). This can be explained by the following argument: large bubbles scatter more than small bubbles (see equation ( 4)), and their average correlation distance r c is also larger.

When the range of measured values of g is considered, the agreement is better or worse, depending on whether the high or the low limit of the possible g function is considered. It follows that, given the precision of our correlation measurements, we cannot conclude whether 12: (Color online) The effect of polydisperse correlations as predicted by the self-consistent approach (SCA). Three g(R, r) functions are considered (see figure 11): hardsphere correlation (dashed lines), high and low estimates of g measured from the x-ray data (dotted and thick lines, respectively).

TABLE III:

A simple criterion for estimating how well the predictions of the different models agree with the experimental data: the relative mean squared differences between theoretical and experimental effective wavenumbers (average of the imaginary and real parts) over the 30-300 kHz frequency range. Different g functions are considered: from the hardsphere approximation (hs), and from the measured correlation (see Fig. 11 or not the polydisperse correlations are sufficient to fully explain the discrepancy. Nonetheless, we find that plausible values for g are able to give good agreement. A more precise determination of the correlations would be necessary for a definitive conclusion to be drawn.

V. CONCLUSIONS

We have investigated wave propagation through a complex medium in which low-frequency scattering resonances lead to strong dispersion, with large peaks in the velocity and attenuation. By choosing to work with a relatively simple acoustic system -bubbles in the elastic medium PDMS -we have shown how these signatures of strong resonant scattering are influenced by the coupling between the scatterers. This leads to the breakdown of the independent scatterer approximation (ISA), which is commonly used to interpret experimental data. A crucial step in investigating the failure of the ISA in this system has been our ability to characterize the physical properties of bubbly PDMS very carefully, including independent measurements of the elastic properties of the matrix as well as the concentration and size distribution of the scatterers. Remarkably, we find that this failure of the ISA occurs even though the concentration of scatterers is low, only 2% in our case. It is worth noting one important consequence: if the ISA is used to estimate the size distribution of the scatterers, significantly smaller sizes are found than by direct imaging methods.

In this paper, we have proposed that this discrepancy can be explained by accounting for the role of correlations in the scatterers' positions, which we were able to probe directly using x-ray tomography. Following Keller's approach, we have shown that these correlations brought a non-negligible correction to the ISA, and that the agreement with the experimental data was improved, although not perfect. We have proposed an alternative approach, based on a self-consistent argument, with which we were able to take into account the effect of polydispersity on the correlations: bubbles of different radii are not correlated with their neighbours in the same way. Again, x-ray tomography allowed us to investigate this polydisperse correlation, which was found to be different from the monodisperse estimate: small bubbles are more likely to be close to other bubbles, whereas big bubbles are often isolated. We have shown that, within the possible polydisperse correlations that are consistent with our xray measurements, the SCA was able to provide a satisfactory explanation for the observed shift to higher frequencies in the acoustic resonances in the sample. It should also be noted that our model could be refined to account for more complex situations. In particular, when the correlations are long-range in a thin sample, one can expect the boundaries to play a role in the correlation, an effect we did not incorporate. Also, we looked at correlation effects or inhomogeneity effects, but we have not considered the case in which both effects are coupled.

An interesting question is the following: why are correlations important in the bubbly medium studied here, whereas they have not been detected before for bubbles in water [26] or bubbles in a yield stress fluid [8]? According to equation [START_REF] Wood | A Textbook of sound[END_REF], the magnitude of the correction due to the positional correlations is proportional to the volume fraction of scatterers Φ, which was generally lower in pre-vious studies (5 × 10 -4 to 10 -2 at most), so that β would be expected to be low. However, it is interesting to note that some previous studies did report deviations between Foldy's prediction and the experimental data that might be interpreted as a manifestation of correlation effects. In the historical data by Silberman [27], for instance, the measured attenuation was much larger than Foldy's estimate. Feuillade proposed that the discrepancy could be resolved by introducing an arbitrary fitting parameter [28] that one can interpret as a correlation length (see appendix). Wilson et al. also observed deviations from Foldy's estimate, as their attenuation peak was at a lower frequency than what was expected from the size analysis (see Fig. 10 in [26]). They attributed this discrepancy to the uncertainty in their measurement of the bubble size distribution, but one can imagine that correlation effects were responsible for this shifting of the attenuation peak. Indeed, large depletion zones around bubbles have been observed in bubbly liquids [29], which would lead to positional correlations. In our sample, the exact mechanism leading to the positional correlations has not been identified unambiguously. Ostwald ripening might induce the disappearance of small bubbles close to bigger bubbles, hence leading to the isolated big bubbles we observed. Another possible scenario could involve the rotation of the cell when the sample was prepared. If one considers that only the buoyancy and the Stokes drag force are relevant (inertial forces were negligible in our set-up), a bubble is expected to move along a circle whose radius is inversely proportional to the rotation frequency, and proportional to the square of the bubble radius. The movement of the bubbles may induce hydrodynamic interactions leading to positional correlations. This hypothesis is supported by the observation of strong clustering when the cell was rotated at a lower speed (0.7 rpm), i.e., when the radii of the circles followed by the bubbles were larger.

to the size of the scattering inclusions, which means we can consider that the waves scattered by the inclusions are spherical: when scatterer i is excited by a wave with amplitude p i , it generates a wave p i f i s exp(ik 0 r)/r at distance r, where f i s is its scattering function and k 0 the wave number in the pure medium. The total field experienced by inclusion i, due to the other scatterers, is given by the self-consistent relation

p i = j =i p j f j s e ik0rij r ij , (A1)
where r ij is the distance between scatterers i and j. In the case of a finite number of scatterers, Eq. (A1) can be solved and all the p i exactly computed. For an infinite number of inclusions, however, taking into account the full multiple scattering is a difficult task and the p i are not easily accessible. In particular, the 1/r range of the interaction means that no pair (i, j) can be neglected a priori in Eq. (A1).

Let us further simplify the problem by assuming that all the scatterers are identical (f i s = f s ) and look for modes that take the form of plane waves propagating in the x direction: p i = P exp(ikx i ), where k is the effective wave number we want to determine, and P the amplitude of the mode. Then, Eq. (A1) reduces, for scatterer i = 0 arbitrarily chosen as the central one located at r = 0, to

1 = j =0 e ikxj f s e ik0rj r j = f s j =0 e i(krj cosθj+k0rj ) r j , (A2)
where r j is the distance to scatterer j, and θ j the angle between the direction of propagation and the position of scatterer j (see Fig. 13).

One can approximate the discrete equation by a continuous equation, taking into account the positional correlations with the pair correlation function g(r): Invoking a small imaginary part of k 0 to preserve the convergence of the integral, one can then calculate the effective wave number: with experimental data for bubbly liquids (see equation (30) in [30]). Feuillade's model has been subject to criticisms [START_REF] Henyey | [END_REF]31], in particular because, for unclear reasons, it assumes that the interactions between the scatterers are confined within a finite range. It might be that this finite range is related to correlations. Note that the approach we develop here is similar to Feuillade's. The differences are that (1) we abandon the assumption that all the scatterers are in phase, which is true only locally, and (2) we include the positional correlations. For a polydisperse assembly of scatterers, the same self-consistent scheme can be used. In this case, the total field p i depends a priori on the radius of the bubble, p i = P (R i )exp(ikx i ), meaning that equation (A3a) must be generalized into P (R) = f s (R) n(R ′ )dR ′ P (R ′ ) × g(r, R, R ′ ) e i(krcosθ+k0r) r d 3 r, (A5)

k 2 = k 2 0 + 4πnf s 1 + 4πnf s (1 -g(r))
where g(r, R, R ′ ) gives the probability of finding a bubble with radius R ′ at a distance r from a central bubble with radius R. In the general case, equation (A5) cannot be directly simplified into a dispersion relation. However, if we assume that the correlation does not depend on R ′ (g(r, R, R ′ ) = g(r, R)), a further integration over R yields 1 = n(R)f s (R)dR g(r, R) e i(krcosθ+k0r) r d 3 r, (A6) from which we obtain the dispersion relation

k 2 = k 2 0 + 4πn(R)f s (R)dR 1 + A , (A7)
with A = 4πn(R)f s (R)dR (1 -g(r, R)) sin kr k e ik0r dr.

FIG. 2 :

 2 FIG. 2: (Color online) Example of an image taken for the bubble size analysis (optical characterization).

3 FIG. 3 :

 33 FIG.3:(Color online) Histogram of the bubble radius distribution for the sample in figure2, obtained from the image analysis. The superposed curve (circles plus connecting lines) corresponds to the result of the x-ray tomography analysis.

Figure 5

 5 Figure 5 offers a further step in the comparison process. Predictions of the ISA model are plotted as dotted lines.

FIG. 4 :

 4 FIG.4: (Color online) Attenuation α and phase velocity v in the bubbly PDMS sample, measured from a single pulse (open circles) and from an average pulse (solid circles). The velocity in the pure PDMS is shown by the solid line. The attenuation in pure PDMS is, on average, 2 orders of magnitude less than in the bubbly sample (too small to be visible on the graph).

FIG. 5 :

 5 FIG.5:(Color online) Close-up of the attenuation and velocity peaks. Solid circles correspond to experimental data (same as in Fig.4). Comparison with Foldy's model is shown for monopolar scattering of the bubbles (dotted lines) and monopolar + dipolar scattering (solid lines). Note that there are no adjustable parameters in this comparison of the theory with experiment. Results of fitting the model to the data by allowing either the size distribution (dashed lines) or the shear modulus (triangles) to vary are also shown. In the inset, the size distribution determined by the x-ray experiment (solid line) is compared to the one inferred from the fit to the ultrasonic data (dashed lines).

FIG. 6 :

 6 FIG.6: (Color online) Shear moduli as functions of frequency. The lines correspond to the results for experimental µ ′ (solid line) and µ ′′ (dashed line), as given in equations[START_REF] Henyey | [END_REF]. Open triangles show the fitted values for µ ′ , solid ones for µ ′′ .

FIG. 7 :

 7 FIG. 7: (Color online) The layered structure of the sample is revealed by the x-ray tomography. Along the thickness h of the sample, we report the number of bubbles per unit area n ′ (left) and the radii distribution histogram (right) for each layer.

FIG. 8 :

 8 FIG.8: (Color online) Amplitude (top) and phase (bottom) of the acoustic transmission through the sample. Because the data are the same as in figure5, the same discrepancy between Foldy's prediction for a homogenous sample (solid lines) and the experimental data (circles) is obtained. The prediction for a layered medium (dashed lines) does not bring better agreement.

FIG. 9 :

 9 FIG.9: (Color online) Function g as measured from the xray tomography data (circles). Error-bars correspond to one standard deviation. The solid line corresponds to the hardsphere approximation for the bubble size distribution measured in the sample.

FIG. 10 :

 10 FIG. 10: (Color online) Investigation of the effect of correlations. The experimental attenuation and velocity (circles) and the ISA prediction (solid thin lines) are shown again, as in figure 4. Keller's approach (Eq. (15)) is plotted for the hard-sphere correlation (dashed lines), for the high estimation of g measured from the x-ray data (dotted lines), and for the low estimation of g (thick solid lines).

FIG. 11 :

 11 FIG.11:(Color online) Function g(r, R) as measured from the x-ray tomography data (circles) on three populations of bubble radii. Error-bars correspond to one standard deviation. The solid line corresponds to the hard-sphere approximation for the bubble size distribution measured in the sample.

  FIG.12:(Color online) The effect of polydisperse correlations as predicted by the self-consistent approach (SCA). Three g(R, r) functions are considered (see figure11): hardsphere correlation (dashed lines), high and low estimates of g measured from the x-ray data (dotted and thick lines, respectively).

1 =

 1 f s d 3 r × ng(r) e i(krcosθ+k0r) ) sin(kr)e ik0r dr.(A3b)

FIG. 13 :

 13 FIG.13:(Color online) An infinite 3D medium with spherical scatterers is considered. The position of scatterer 0 is arbitrarily chosen as the centre of the coordinates and we look for modes of the system consisting of plane waves propagating in the x direction. The dispersion relation is obtained by calculating the total field experienced by scatterer 0 due to the other scatterers.

TABLE I :

 I Parameters of the bubble size distribution, extracted from the optical and x-ray measurements.

	volume fraction average radius polydispersity

TABLE II :

 II Five characteristic quantities measured experimentally (first line), and predicted by the model (second line).

	fmax αmax	α hf	v lf	v hf
	(kHz) (mm -1 ) (mm -1 ) (m/s) (m/s)
	exp. 100	3.5	0.40	260 1180
	mod. 90 a	2.8	0.35	224 1020

  ) with the lower (xp-) middle (xp) and upper (xp+) estimates.

	Model ISA	Keller	SCA
	g	/	hs xp -xp xp + hs xp -xp xp +
	(%) 14.8 11.3 7.3 12.4 17.7 10.1 4.3 10.8 19.6

  Note that equation (A4) reduces to Foldy's equation if the denominator is approximated by 1, and to Keller's equation if it is expanded to second order in nf s . Interestingly, it has some similarities with the equation proposed by Feuillade, which was found to give excellent agreement

	sin kr k e ik0r dr	. (A4)

Acknowledgments

Support from the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. The authors wish to thank A. Derode, M. Devaud, T. Hocquet, and J.-C. Bacri for numerous discussions, as well as R. Jayaraman for his help in enabling the x-ray experiments to be carried out.

Note that, as in the ISA, polydispersity modifies the equation by changing nf s into n(R)dRf s (R). But polydispersity also brings another non-trivial modification: the correlation can depend on the radius of the cen-tral bubble considered, which makes different correlations possible for different bubble sizes, as in our bubbly sample.

Appendix A: Self-consistent approach for multiple scattering of waves

We propose a self-consistent approach for calculating the dispersion relation of a random multiply scattering medium, taking into account the positional correlations of the scatterers. Let us consider an infinite medium with n inclusions per unit volume. We limit ourselves to low frequencies such that the wavelength is large compared