Scattering of acoustic waves by a nonlinear resonant bubbly screen

Abstract Our study concerns the propagation of acoustic waves through a thin screen made of a periodic arrangement of air bubbles in water. The bubbles are oscillators of the Minnaert type whose dynamics is modified by the containment. This nonlinear dynamics is obtained in the time domain using asymptotic analysis and a homogenization technique involving three scales, those being the scale of a bubble, that of the array and eventually that of the wavelength. The resulting effective model is set in the water (the screen has disappeared) and it encapsulates the effect of the screen in a jump of the normal acoustic velocity. The jump is linked to the continuous version of the bubble radius which satisfies an equation of the Rayleigh–Plesset type. This allows us to highlight two important effects. Firstly, a bubble within the array has a much larger radiative damping than an isolated bubble. Secondly it perceives a pressure which differs from the acoustic pressure imposed by the source due to bubble–bubble interactions; it results in a term of mass correction deduced from the Green's function for a Laplace problem which accounts for the bubble arrangement. Our findings are exemplified by numerical experiments of the scattering of a short pulse in the linear and nonlinear regimes.


Introduction
In 1985, Caflisch and co-workers analysed the wave propagation in bubbly liquids (Caflisch et al. 1985a,b) offering a rigorous mathematical framework to the former analysis developed by Van Wijngaarden (1968). Using asymptotic analysis, they derived an effective wave equation involving a continuous version of the bubble radius satisfying In some practical situations, the region of the bubbles is reduced to one or a few layers (figure 1). The most famous examples are the anechoic tiles, codenamed Alberich, developed during the Second World War by the German Navy; these tiles are the building blocks of a rubber net containing bubbles used to block the signals of sonars (Gaunaurd 1977;Hladky-Hennion & Decarpigny 1991). Other examples are the bubble nets used to protect underwater structures from damage by underwater explosions (Domenico 1982) or those created by some of the marine mammals to catch fish (Leighton 2004). Recently, these bubbly screens with subwavelength thickness have been revisited in the context of metamaterials. Leroy et al. (2009) realized controlled experiments to study the transmission of ultrasound by a single layer of bubbles in water; Bretagne, Tourin & Leroy (2011) realized a bubble raft close to an interface with air or sandwiched in air and analysed the effects of such interfaces; Leroy et al. (2015) demonstrate the ability of the metascreens to produce superabsorption and Lanoy et al. (2018) obtain perfect absorption by critical coupling. Eventually an original metasurface has been recently proposed by Schnitzer, Brandão & Yariv (2019) with cylindrical bubbles trapped in the grooves of a microstuctured hydrophobic wall which is shown to be capable of supporting guided waves analogues to spoof plasmon waves.
From a theoretical point of view, the first study on the interaction of acoustic waves with an array of bubbles dates from 1966 (Weston 1966). Although Weston concludes that 'a plane array behaves as a plane screen of gas -there is no resonance at all', his study contains almost all the elements of the model of Leroy et al. (2009). In both cases, the analysis is based on multiple scattering theory in the linear regime but Weston consider an over-simplified scattering function of a single bubble. In contrast, the model of Leroy et al. (2009) shows that bubbles within a screen experience a radiative damping much larger than an isolated bubble does; also that the resonant frequency of the screen is higher than that of an isolated bubble due to bubble-bubble interactions. A different approach has been proposed by Ng & Ting (1986) and Miksis & Ting (1989) who use the effective wave equation derived in Caflisch et al. (1985a,b) for a bubbly liquid. By considering the limiting case of a thin layer, they derive transmission conditions through an effective screen in the nonlinear regime and in the time domain. This model recovers the large damping of the screen but misses the shift of the resonant frequency. Eventually, Ammari et al. (2017a) provides almost exact analytical solution of the scattering coefficients in the linear harmonic regime within a rigorous mathematical framework when the screen is placed in the vicinity of a Dirichlet wall; their work applies mutatis mutandis to an isolated screen.
In parallel to these works, one approach has become popular to describe the collective dynamics of a cluster of bubbles. It consists of enriching the Rayleigh-Plesset (RP) equation with a term encapsulating the effect of bubble-bubble interaction. Starting with Harkin, Kaper & Nadim (2001) and Doinikov (2001), who consider the interaction of two bubbles, the model has been further extended heuristically to a system of N bubbles. With neglect of the surface tension and liquid viscosity, this heuristic RP equation reads where ρ is the mass density of the liquid, R i is the radius of one bubble within the cluster and R eq its value at equilibrium, p eq the pressure at equilibrium and p a the forcing acoustic pressure (γ stands for the adiabatic index). In this equation, the sum in the right-hand side term is assumed to account for the effect of the (N − 1) bubbles with radii R j at distances h ij of the ith bubble, see e.g. Bremond et al. (2006), Guédra, Cornu & Inserra (2017). One obvious drawback of (1.1) is that 1/h ij diverges for large N while R 2 jṘ j is bounded. The reason is that this sum is the contribution of the bubbles j to the pressure seen by the bubble i; as such, it should contain the term e ikh ij with k the wavenumber, as in Weston (1966) and Leroy et al. (2009). Hence, (1.1) holds for a cluster extension much smaller than the typical wavelength.
In this study, we adapt the asymptotic analysis of Miksis & Ting (1989) to the case of a thin bubbly screen. The calculations are performed in the time domain and preserving the nonlinear dynamics of the bubble oscillations. To do so, we focus on the limit of sparse arrays with spacing much larger than the bubble radius, which allows for the resonances to take place; in the other limit, a dense array would behave basically as a wall. As the typical wavelength remains larger than the spacing, three scales can be defined. In § 2, we discuss the meanings of the problems obtained at these three scales and their matchings; the technical calculations have been collected in appendix A. The whole analysis results in effective transmission conditions with a jump of the normal acoustic velocity dictated by a modified RP equation, see forthcoming (2.5) (figure 2). In § 3 results of the model are discussed in the linear and nonlinear regimes for a short incident pulse. In particular, it is shown that for one-dimensional propagation, a different (although equivalent) formulation of the effective problem provides a simple physical interpretation of the mechanisms of interaction, see forthcoming (3.2); our modified RP x r FIGURE 2. The effective problem: an equivalent screen at x = 0 across which jump conditions (2.5) apply, with R(r, t) being a continuous version of the bubble radius.
equation is similar although not equivalent to (1.1) and it generalizes to the nonlinear regime the result of Leroy et al. (2009). The properties of the model in terms of energy conservation are discussed in § 4. Eventually, concluding remarks and perspectives are given in § 5. The numerical results reported throughout the paper use the following parameters: ρ = 10 3 kg m −3 , c = 1500 m s −1 , p eq = ρ g c 2 g γ = 0.225 atm., γ = 1.4, R eq = 10 μm, with ρ g, and c g, the mass density and speed of sound in the gas and in the liquid at equilibrium. Different spacings h from 0.1 to 2 mm will be considered. The separation of the three scales is verified with a wavelength λ M 1 cm at the Minnaert resonance (1.3)

The actual problem and the effective screen problem
The effective screen problem aims to simplify the actual problem of the interaction of acoustic waves with an array of gas bubbles. In the actual problem, the array is located at x = 0 and it has the periodicity h in both directions (figure 1). We assume a weak compressibility of the liquid, an irrotational flow and the effects of the viscosity and surface tension are disregarded. Let ρ, p and u be the density, the acoustic pressure and the acoustic velocity which are functions of the time t and space x = xe x + r (with r · e x = 0).
In the liquid and in the gas, (ρ, p, u) satisfy the Euler equations where p = ( p tot − p eq ) with p tot the total pressure. At the equilibrium p tot = p eq , u = 0 and the mass densities in the liquid and in the gas are ρ ,g . For a weakly compressible liquid and a perfect gas under adiabatic transformations (see e.g. Bergamasco & Fuster (2017) for a discussion on the validity of the adiabatic regime), the equations of state are written in the form in the gas, with c 2 g = γ p eq ρ g .
For an acoustic wavelength that is large compared to the bubble radii, each bubble with index i is set in forced radial oscillations R i (t). Thus, at the interface gas/liquid, the pressure and the velocity are continuous and where the dot means the time derivative.

Formulation of the effective screen problem
The effective problem is set in the liquid only; because of the weak compressibility of the liquid and because of the low Mach number, the propagation is linear, hence ( p, u) satisfy the Euler equations The bubbly screen has disappeared, its effects being now encapsulated in a jump of the normal velocity, specifically where R(r, t) is a continuous version of the bubble radius. We have defined [ the jumps of the acoustic pressure and normal velocity across the effective screen (figure 2). In (2.5), the RP equation contains two important contributions. Firstly, the forcing term is the acoustic pressure p(0, r, t) at the equivalent interface; from (2.4a-c) and the form of [[u x ]] in (2.5), p(0, r, t) has a contribution ρ c (R 2Ṙ /h 2 ) corresponding to the radiative damping of the screen. It is termed superradiative damping in Leroy et al. (2015) since it is much greater than the radiative damping of an isolated bubble considered in, for example, Keller & Miksis (1980). Next, the term δ (ρ /h) ∂ t (R 2Ṙ ) is attributable to the bubble-bubble interaction, with δ a constant depending on the arrangement of the bubbles within the array; it will be discussed further in the following section. For the time being, we can notice that the RP equation in (2.5) differs from (1.1); it applies to an infinite number of bubbles since δ is finite with a sign a priori free while δ is negative by construction in (1.1).

2.2.
Results of the asymptotic analysis at the three scales As previously said, the effective screen problem is obtained using asymptotic analysis combined with homogenization. The underlying assumptions are (i) a low frequency regime, (ii) a dilute screen and (iii) a low Mach number, specifically where ω is the typical frequency imposed by the source. These scalings produce a separation of three scales, the scale of a single bubble, that of the screen spacing and that of the typical wavelength λ = 2π(c /ω); each scale is associated with a problem much simpler than the complete actual one and which can be solved at the dominant order and then iteratively at higher orders in ε. At each order, the analysis aims to put those three problems together using matching conditions. Our analysis has been conducted at the dominant, order 0 in ε, see appendix A.2 and at order 1 in ε, see appendix A.3; the results at the two first orders are then gathered to form a unique problem. Below, we summarize the results at the three scales. The microscopic problem results from a zoom in on a single bubble of centre r n . It corresponds to the classical problem of an isolated bubble submitted to a pressure at infinity p ∞ (r n , t), primarily considered by Rayleigh (1917) in the incompressible case, see also Plesset & Prosperetti (1977). Due to the low Mach number, the pressure and mass density are uniform within the bubble, with ∂ t ρ g + ρ g div u = 0. In the liquid, the flow satisfies the incompressible nonlinear Euler equation ρ (∂ t u + u · ∇u) = −∇p. We get at the microscopic scale: A classical form of the Rayleigh-Plesset equation is obtained of the form ρ (RR + 3 2Ṙ 2 ) = p n − p ∞ , where p n (t) is the pressure in the liquid at the interface |x − r n | = R, and in the absence of surface tension p n (t) = p(t). However, at this scale p ∞ (t), being the pressure at infinity seen by a single bubble, is unknown. In figure 3, we report an illustration of the patterns of the radial velocity u · e r and of the projection u · e x for a bubble at r n = 0 and t = 0. We have considered a plane wave at incidence 45 • in the linear regime p inc (x, t) = Δp exp(ik · x − iωt) with Δp = 0.1p eq and ω = 0.7ω M ; the array spacing is h = 50R eq . The velocity fields (2.7) in the bubble and in the liquid are easily obtained owing to the resolution of the linearized version of (2.5), see appendix B. The order of magnitude of the velocity of 0.1 m s −1 is given byṘ ∼ ωΔp/ρ R eq (ω 2 M − ω 2 ). The mesoscopic problem corresponds to the scale of the array spacing, intermediate between the scale of the bubbles and that of the wavelength (figure 4). At this scale, the . Solution of the microscopic problem -a single bubble with pressure p ∞ at infinity (a). Radial velocity (b) and velocity along the x-axis (b) in m s −1 for an incident plane wave at oblique incidence in the linear regime with amplitude Δp = 0.1p eq , see main text.
bubbles are reduced to points at r n . The flow in the liquid is still incompressible but now it is linear. Specifically, for x ∈ Y n the periodic cell containing r n , we have at the mesoscopic scale (in the liquid): and G(x) ∼ |x|→∞ 2π/h|x| (i.e. without constant at infinity), see also (A 34). For a square array, the Green's function can be calculated explicitly using the decomposition (2.10) with = √ n 2 + m 2 and (n, m) ∈ Z 2 \(0, 0). The Green's function is singular at x = 0 since S(x) ∼ x→0 −(h/|x|) + δ. We find δ = 3.9 in excellent agreement with Weston (1966) and Leroy et al. (2009). In both references the authors evaluate the pressure seen by a single bubble within the array using multiple scattering theory in the linear harmonic regime. They introduce a cutoff distance b as the effective inter-bubble distance; in Weston (1966), b = h log γ with log γ = 0.58 the Euler constant and in Leroy et al. (2009) the cutoff is linked to δ through δ = 2π(b/h) which allows us to conclude on the agreement. Interestingly, both authors justify this value by making reference to a problem similar to the Green's function problem (2.9). However, this problem is set for the Helmholtz problem while ours is set for the Laplace problem since k → 0 at this intermediate scale. An illustration of the mesoscopic solution is shown in figure 4. From (2.8), the order of magnitude of the velocity and it diverges at x = r n due to the singularity of the Green's function (see appendix B). -1 FIGURE 4. Solution of the mesoscopic problem -the bubble is reduced to a singularity at r n within a cell Y n infinite along e x and periodic in the two other directions (a). Velocity (u − u ∞ ) · e x along the x-axis from (2.8) resulting from the microscopic solution in figure 3. The profile along x is taken at a distance 0.15 h from the centreline; the profile at the centreline is singular at x = 0 according to (2.9) (b).
Micro Meso FIGURE 5. Micro-to-meso matching -the solution at the mesoscopic scale is the limit of that at the microscopic scale when 2R eq /h vanishes. Two-dimensional fields of the u x velocity at the microscopic scale (a) and at the mesoscopic scale (b); the dotted white lines show the bubble boundary; velocity profiles along the centreline (c).
The singularity of the mesoscopic solution results from the micro-to-meso matching. Indeed, the mesoscopic solution in (2.8) when x/h → 0 has to coincide with the microscopic solution (2.7) when |x − r n |/R → ∞. Hence, from (2.7), the singular behaviour of the pressure has to be of the form p ∼ p ∞ (r n , t) + ρ (∂ t (R 2Ṙ )/|x − r n |), and that of u is given by ρ ∂ t u = −∇p; this is what we recover with (2.8). This matching is illustrated in figure 5. Reporting the fields of both solutions for a small but non-zero R eq /h shows that the mesoscopic solution coincides with the microscopic one in the liquid that it extends up to 0. At this stage p ∞ (r n , t) is still unknown.
The compressibility of the liquid appears in the macroscopic problem. At this wavelength scale, the discrete set of points r n is seen as a continuum hence p ∞ (r n , t) can be replaced by p ∞ (r, t), the same for u ∞ . Now, we aim to know the meanings of p ∞ and u ∞ and this is done by considering the limits x/h → ±∞ of the mesoscopic solution which provide the limits x/λ → 0 ± of the macroscopic solution. Owing to the behaviour of G in (2.8), the macroscopic solution has to satisfy The meaning of the relations for u is unambiguous: FIGURE 6. Meso-to-macro matchings -the solution at the macroscopic scale for x/λ → 0 is the limit of that of the mesoscopic scale for large x/h; the solution at the mesoscopic scale accounts for the jump in the velocity (2.5). Profiles and fields of the solutions (a); panel (b) shows a zoom near the screen. The incident wave is a plane wave at incidence 45 • and ω = 0.7ω M .
as announced in (2.5). Owing to this result, and with ρ ∂ t u = −∇p which applies for ( p, u) and ( p ∞ , u ∞ ), the relations for p appear as the Taylor expansions of the pressure p being continuous at x = 0 with a discontinuous first derivative with respect to x. Specifically, the relations on p in (2.11) can be written as This tells us that [[p]] = 0, as announced in (2.5), and provides the meaning of p ∞ . The pressure p ∞ (r, t) = p(0, r, t) − ρ (∂ t (R 2Ṙ )/h)δ seen by a single bubble has a contribution due to the arrangement of the array and which has been captured by the Green's function at the mesoscopic scale; again, as the Green's function is that for the Laplace problem, δ is independent of the frequency. It is now sufficient to replace p ∞ in the RP equation that we recover from the microscopic solution (2.7) to get the final RP equation in (2.5).
The meso-to-macro matching is illustrated in figures 6 (for ω = 0.7ω M ) and 7 (for ω = ω M ). We have represented the macroscopic solution in reflection and transmission for |x|/λ > h/λ and the mesoscopic ones for |x|/λ < h/λ; this choice of the x-range is somehow arbitrary, it simply means the close vicinity of the array. It is visible in both cases that the mesoscopic velocities are basically reduced to constant values u ∞ (r n , t) ± (2πR 2Ṙ /h 2 )e x from (2.8); also visible is the fact that these values coincide with the macroscopic velocities u(x → 0 ± , r ∈ Y n , t).

Analysis of the effective model
In this section we shall analyse the effective model (2.5) thanks to numerical simulations. The one-dimensional compressible Euler equations in the liquid are solved using the finite element method presented in Fuster, Dopazo & Hauke (2011) and the FIGURE 7. Same representation as in figure 6 for ω = ω M resulting in almost perfect reflection.
transmission condition is accounted for by using dual nodes at the location of the bubbly screen, each node representing the two side solution u x (0 ± , t). Next, R is obtained by integrating the RP equation in (2.5) with a classical Runge-Kutta method. Following Doc et al. (2016), we consider an initial Gaussian pulse of the form and we use σ = λ M /10. Hence, with a duration T/10, and T = λ M /c ∼ 6 μs, the pulse acts almost as a delta Dirac function and excites a wide frequency band around ω M . Results are reported in figure 8 against h/R eq in the linear and nonlinear regimes. The two time scales, the duration of pulse T/10 and the free oscillation period T, are visible in the response R(t) of the bubbly screen and in the resulting acoustic pressure p(0, t). Next, two limiting cases are observed. For large inter-bubble distances h/R eq , each bubble behaves as if it were on its own. The radiative damping is small, which allows for long time oscillations but as a counterpart, the screen has a weak interaction with the liquid (p(0, t) ∼ p inc (0, t)). In the opposite limit of small inter-bubble distances, the radiative damping is large, resulting in over-damped oscillations and strong interaction with the liquid (p(0, t) significantly departs from p inc (0, t)). This is consistent with the intuitive idea that dense arrays of bubbles act basically as perfect shields while sparse ones are transparent. From figure 8, the intermediate regime is for h/R eq equal to a few dozen.
The figure 9 show snapshots of the acoustic pulse p(x, t) in the liquid in this intermediate regime (h/R eq = 20) for increasing Δp. While the propagation in the liquid is linear, increasing Δp affects the response of the nonlinear oscillator (the screen). It is visible that the signal emitted by the screen, with nonlinear shape, is transmitted to the liquid, in reflection and also, but less visible, in transmission. What is also visible is the fact that increasing the nonlinearities weakens the interaction of the screen with the liquid.

Analysis in a one-dimensional problem
We shall now simplify the screen model (2.5) to get physical insights into the interaction mechanisms. Restricting ourselves to one-dimensional problem allows for the emergence of a modified RP equation which contains explicitly the radiative damping term and which is independent of the propagation in the liquid; this later will appear as a simple byproduct 3.1.1. Interpretation in terms of pressure radiated by the screen One-dimensional propagation involves acoustic pressure of the form p( The same function f is involved for x ∈ R ± since in our problem (2.5) the pressure at x = 0 is continuous; whence f (t) = p(0, t) − p inc (0, t). Next, the propagation in the liquid is governed by the Euler equations (2.4a-c) which Eventually the jump in the velocity in (2.5) provides f (t) = 2πρ c (R 2 /h 2 )Ṙ. It follows that the problem (2.4a-c) and (2.5) can be written as for x ∈ R ± (figure 10). This formulation provides a different although equivalent interpretation of the mechanism of interaction. In a first step, the incident pulse excites the screen; p inc (0, t) is the forcing term in the RP equation, (3.2), and it is worth noting that the superradiative damping is now explicitly accounted for by the term in R 2Ṙ . In a second step, the screen radiates pressure fluctuations in the liquid, symmetrically on the right and on the left. These fluctuations p rad (t ∓ x/c ), initially generated through a nonlinear mechanism at x = 0, propagate linearly in the liquid; they contain the two time scales, that of the forcing term p inc and that of the natural oscillation (under-or over-damped). For our short incident pulse p inc the radiated pressure lasts much longer than the pulse as sketched in figure 10. Incidentally, and from a practical point of view, the resolution of (3.2) is simpler than that of (2.4a-c) and (2.5). Indeed, the RP equation is set in time only, hence it can be solved once and for all, and afterwards the solution in the liquid is explicitly known. In figure 10(b), the two numerical solutions of (2.4a-c) and (2.5) and of (3.2) are shown to coincide.

One-dimensional problem in the linear regime
In the linear regime, with R = R eq + r, r R eq , (3.2) simplifies further to = 10 x σ 2), the screen radiates pressure wave trains symmetrically on the right and on the left after the incident signal has passed through it. The solution is the sum of the two contributions p = p inc + p rad from (3.2). (b) Resulting pressure profiles p in the linear and nonlinear regimes, blue lines from (3.2) and dotted red line from the resolution of (2.5).
The analysis can be conducted frequency by frequency by simple Fourier transform. We use p inc (x, t) = p inc (ω) exp(−i(ωt − kx)) dω and p rad (t ± x/c ) = p rad (ω) exp(−i(ωt ± kx)) dω, with k = ω/c the wavenumber. The reflection R and transmission T coefficients read and [ω 2 M − ω 2 (1 − δR eq /h + iKR eq )]r = −p inc /ρ R eq ,p rad = −iKρ R 2 eq ω 2r with K = 2π/kh 2 , hence with the radiative damping contained in KR eq and a mass correction contained in δ(R eq /h). (Obviously, the same result for (R, T ) is found by linearizing (2.5) with u x = −i/(ρ ω)∂ x p from (2.4a-c).) As previously said, the damping of the array is much larger than that of an isolated bubble; this latter is given by kR eq and we have k/K = (kh) 2 1. The mass correction is the product of δ with R eq /h. The constant δ depends on the arrangement of the bubbles within the array and it is independent of R eq and h. Hence, for the same arrangement, the mass correction vanishes for sparse arrays R eq /h → 0. From (3.5a,b), the terms of bubble-bubble interaction and of radiative damping vanish for large h resulting in R 0 (the screen is transparent for the acoustic waves). In the opposite limit of dense arrays R eq /h ∼ 1, hence KR eq 1, resulting in R −1; the screen becomes equivalent to a wall of air with vanishing acoustic pressure. In the intermediate regime, R has a maximum at a frequency ω M / 1 − δR eq /h shifted above the Minnaert frequency as observed in experiments (Leroy et al. 2009. It is worth noting that nonlinearities also produce a shift of the resonance frequency. This is illustrated in figure 11 where we report the Fourier transforms of the reflected and transmitted pressure fields for increasing Δp/p eq ∈ (1, 50) in our Gaussian pulse (3.1). By analogy with the linear regime, the nonlinear resonance frequency is identified at the maximum in the reflected signal. Expectedly, the frequency shift in the linear regime FIGURE 11. Fourier transforms of the reflectedp(x < 0, ω) and transmittedp(x > 0, ω) pressures for increasing nonlinearities Δp/p eq (h/R eq = 40).
(Δp/p eq = 1) for h/R eq = 40 is small and we recover (3.4). In the nonlinear regime, the screen imparts to the liquid acoustic pressure fluctuations with complex frequency content. It is in particular visible that (i) for moderate nonlinearities (up to Δp/p eq ∼ 10) the shift of the resonance to higher frequency is significantly increased from ω M to approximately 3ω M in the reported case, (ii) strong nonlinearities produce the appearance of harmonics and multiple local maxima and minima in the reflection spectrum; in this case, the fundamental resonance is shifted back to ω M .

Conservation of the energy in the screen model
In the absence of viscous loss, the sum of the energy in the liquid and of that in the bubbles is conserved in time. This has to be true in the effective model (2.4a-c) and (2.5) too. To check that this is indeed the case, we consider the equation of energy conservation associated with the Euler equations in the liquid for any bounded domain Ω; E is the acoustic energy and Φ = ∂Ω p u · n dS the flux of the Poynting vector. However, since x = 0 is a surface of discontinuity, Φ has a contribution 2) It is easy to see that this flux is the time derivative of an effective energy Φ scr = (d/dt)E scr with where e p = 0 at equilibrium (note that, to derive (4.3), we used notably that 2 )), and Here, M g is the (constant) mass of a bubble and M rad , sometimes termed the 'radiation mass', is the effective, or dynamic, mass obtained in the linear regime in (3.3). Expectedly, the effective energy is similar to the classical energy of a single bubble in an incompressible liquid, with the work necessary to change the bubble radius from R eq to R. Next, e c is the kinetic energy invested in a (incompressible) liquid during bubble oscillations; it takes the form e c = 1 2 ∞ R u 2 (4πρ r 2 dr) for an isolated bubble. For the array, it can be written where we recover the cutoff distance h/δ introduced by Weston (1966) and Leroy et al. (2009). Recent works discuss the conservation of energy for a compressible liquid, see e.g. Devaud et al. (2008) and Wang (2016), and the effects of confinement (Leighton 2011). As in these references, we find that the energy which is conserved is not the 'local energy' E scr of the bubbles; the liquid and the bubbles can exchange energy hence, in the absence of viscous losses, only (E scr + E) is conserved. We report in figure 12(a) an example of the time variations of E scr (t) along with those of E R (reflected, i.e. computed in x < 0) and E T (transmitted, i.e. computed in x > 0) in the linear regime. In (4.1) E = E R + E T is the energy in the liquid and the energy conservation applies to (E + E scr ). In agreement with the snapshots in the figure 8(a), the screen takes a significant part of the acoustic energy during the transit of the incident pulse, up to 40 % of the total energy in the reported case. It then releases this energy over a time which scales with T. In the present case, the regime is over-damped and the energy release varies as ∼ exp(−5.6(t/T)), in agreement with (3.3) (δR eq /h 0.2, KR eq 2.4).
FIGURE 12. Normalized energies. (a) Example of time variations of E R and E T (energies in the liquid, computed for x < 0 and x > 0, with E = E R + E T in (4.1)), and E scr from (4.3); grey dashed line shows (E + E scr ) 1 (h/R eq = 20, Δp/p eq = 10 −3 ); (b) E R and E T against h/R eq in the linear regime (Δp/p eq = 10 −3 ); (c) E R and E T against Δp/p eq (h/R eq = 20).
Panels (b,c) further quantify the influence of the array spacing and of the nonlinearities. We have considered a sufficiently large time so that almost all the energy is in the liquid. As previously said, going from sparse to dense array increasing h/R eq produces a transition from a perfect shield (total reflection) to a transparent screen (total transmission). Next, incident signals with high nonlinearities are less sensitive to the screen; for Δp/p eq larger than approximately 50, the reflected energy E R is negligible.

Concluding remarks
We have derived an effective model which aims to reproduce the scattering properties of a bubbly screen with resonances of the Minnaert type. The model enables to reduce the effect of the screen to a jump of the acoustic velocity in the liquid coupled to a modified RP equation. The RP equation contains the super-radiative damping of the array being much larger than the classical damping of an isolated bubble. Next, a more significant result is the derivation of the term of bubble-bubble interaction in the RP equation dictating the nonlinear dynamics of the array. This contribution was absent from the analysis of Caflisch et al. (1985a,b) for a bubbly liquid and its extension to a bubbly screen (Ng & Ting 1986;Miksis & Ting 1989); indeed its derivation requires the analysis to be conducted one order further than the dominant one considered in these references. Our result generalizes the findings of Leroy et al. (2009) to the nonlinear regime and to the time domain. It was shown that this contribution was able to describe the shift in the resonance frequency to the higher frequency observed experimentally (Leroy et al. 2009Lanoy et al. 2018); we have stressed that this tendency can be significantly accentuated in the weakly nonlinear regime and inverted with the occurrence of harmonics.
Some extensions of the present work are easy, some others are more tricky. We have neglected the effects of surface tension and viscosity; as in Schnitzer et al. (2019), this is partially justified by the large, dominant, radiative damping of the array. However, there are no particular difficulties in including them, as has been done in Caflisch et al. (1985a,b), and there is little doubt that this will affect the Rayleigh-Plesset equation only. This may become necessary, for instance to explain the perfect absorption by critical coupling (a balance between viscous loss and radiative damping) studied by Lanoy et al. (2018). It is also straightforward to account for different arrangements within the array which would affect the mesoscopic solution by modifying the cross-section of the unit cell for the Green's function in (2.9). Next, accounting for different shape of the bubbles would affect the microscopic solution, namely the response of a single bubble e.g. a disk with monopolar resonance (Calvo et al. 2015). Following Ammari et al. (2017aAmmari et al. ( , 2018Ammari et al. ( , 2019, it should be possible to extend our present work to any geometry. The case of arrays of cylindrical bubbles involved in anechoic tiles, see e.g. Hladky-Hennion & Decarpigny (1991) and Sharma et al. (2017), with a long or infinite length in one direction is slightly more involved but of interest; it requires us to use a different scaling adapted to the two-dimensional problem, as done in Schnitzer et al. (2019). For this application, accounting in addition for a viscoelastic matrix rather than for a liquid may be quite involved depending on the considered viscoelastic model. Next, a different extension concerns the case of bubbly liquids. Pursuing the analysis of Caflisch et al. (1985a,b), and using the mathematical framework of Ammari & Zhang (2017), Ammari et al. (2017b), would provide the bubble-bubble interaction in these media by use of a Green's function that is periodic in the three directions; this would allow us to test the validity of (1.1) and in particular to inspect under which circumstances the dynamic mass M rad in (4.4a,b) is higher or lower than the actual bubble mass (δ positive or negative). Eventually, we have assumed that the transformations undergone by the gas are adiabatic and the result holds for isothermal transformations by simply setting γ = 1. A much more challenging extension is to account for the thermal exchanges between the gas and the liquid since the whole asymptotic analysis has to be reconsidered including the complete thermodynamics.

Declaration of interests
The authors report no conflict of interest.

Appendix A. Asymptotic analysis
In this appendix, we shall use non-dimensional variables t → ωt and x → (ω/c )x (with x = (x, r)), where ω is the typical frequency imposed by the source. Denoting by U a the amplitude of the acoustic velocity, we have, from (2.1a-c) and (2.2), u ∼ U a , p ∼ ρ c U a , (ρ − ρ ) ∼ (ρ /c )U a . With Ma = U a /c the Mach number, we shall use dimensionless quantities Note that the bubble radius is not assumed to scale with the Mach number. This choice prevents a simple linear response where the solution would be proportional to Ma. In the low frequency regime, the maximum wavelength imposed by the source is much larger than h and R eq ; this is accounted for with ε = (ω/c )h 1. Next, we impose the

Microscopic scale Mesoscopic scale
Macroscopic scale . The separation of the scales -at the microscopic scale, the problem is that of a single bubble with a radius of unity in a liquid with x μ ∈ R 3 . At the mesoscopic scale, the bubble is reduced to a singularity within a periodic cell Y with x m = x m e x + r m , x m ∈ R and r m ∈ Y , Y = (− 1 2 , 1 2 ) 2 . At the final, macroscopic scale the whole bubbly screen is reduced to an effective screen at x = 0 across which jump conditions apply (with x = xe x + r, x ∈ R 2 ).
following scalings: Doing so, we assume that the bubbly screen is dilute since R eq /h = O(ε). Such separation of scales was already used in Caflisch et al. (1985b) for bubbly liquid, and in this case, the separation was even more demanding with R eq /h = O(ε 2 ). The Minnaert resonance is dictated by the contrast in the mass densities only, hence, we set ρ g /ρ = O(ε 4 ) and c g /c = O(1). The scalings are such that the resonance takes place in the low frequency regime, specifically ω M h/c = O(ε). Eventually, the scaling for the Mach number produces a linear propagation of waves at large scale in the liquid while keeping nonlinear responses of the bubbles, see Lombard, Barrière & Leroy (2015). The asymptotic analysis requires the use of the rescaled dimensionless coordinates x at the macroscopic scale of the wavelength, x m at the mesoscopic scale of the array spacing and x μ at the microscopic scale of the bubble radius (figure 13), with A.1. Non-dimensional form of the problem The problem (2.1a-c) and (2.2) now reads as in the gaz: in the liquid: with initial conditions (i.c.) and boundary conditions at the bubble interface (b.c.) b.c. p, u · n are continuous and u · n = 1 ε 2 r eq mṘ .
A.2. Derivation of effective transmission conditions at the dominant order The separation of the scales requires the asymptotic analysis to be performed at three scales, the microscopic scale being that of a single bubble, the mesoscopic scale of the array where the bubbles are reduced to points and eventually the macroscopic scale of the waves, the scale at which the whole screen is reduced to an interface x = 0 (figure 2).

A.2.1. Order 0 -resolution at the microscopic scale: the scale of the bubble
The microscopic scale is the scale corresponding to a zoom on a single bubble at a given time t. At this scale, the spherical coordinate x μ = (1/r eq R(r, t))x/ε 2 (see (A 3)) where, following Caflisch et al. (1985b), the discrete version of the bubble radius R i (see e.g. (2.3)) has been replaced by a continuous counterpart R(r, t). Doing so, the bubble has a unitary radius and the nearby bubbles have been sent to infinity. We consider the following expansions of the solution in the gas: (A 8a,b) Accordingly, denoting ∇ μ = ∂/∂ x μ and ∇ r = ∂/∂r, the differential operators read Resolution inside the bubble -using the expansions (A 7) in (A 4) along with (A 9), the problem at the dominant order gives uniform pressure and mass density with ∇ μp and initial and boundary conditions, from (A 6), of the form i.c. at t = 0,p 0 where e r = x μ /|x μ |. Now, we shall expressp 0 μ (r, t) andρ 0 μ (r, t) as a function of R 0 (r, t). By integrating the balance of mass over the gas bubble (first equation in (A 10a-c)) we get μ · e r ds = 0. This differential equation onρ 0 μ is solved using in (A 11) (i) the boundary condition onû −2 μ and (ii) the initial conditionsρ 0 μ = 0, R 0 = 1 at t = 0. Also, making use of the equation of state in (A 10a-c), we obtain the pressurep 0 μ . Eventually, the velocityû −2 μ is a byproduct ofρ 0 μ in (A 10a-c). The result iŝ given in dimensional form in (2.7). Resolution in the liquid -in the liquid, we repeat the exercise, using the expansions (A 8) in (A 5) and (A 6) along with (A 9). We get at the dominant order with initial and boundary conditions i.c. at t = 0, p 0 For the time being, the continuity of the pressure at |x μ | = 1 is not needed, but we have introduced the (unknown) pressure field P 0 (r, t); its meaning will appear later on. We have also assumed that the bubble remains spherical. From (A 13), the flow is incompressible and irrotational. Hence it is associated with a velocity potential φ and, the bubble remaining spherical, φ = (A/|x μ |) + B. Next, the boundary condition in (A 14) for u −2 μ = ∇ μ φ provides A = −r eqṘ 0 /m. Using u −2 μ in (A 13), the pressure p 0 μ is found to satisfy a Bernoulli equation ∇ μ ( p 0 μ + (r 2 eq /m) F(R 0 ,Ṙ 0 ,R 0 , |x μ |)) = 0, that we integrate using p 0 μ → |x μ |→+∞ P 0 (r, t) from (A 14). Doing so, the fields (u −2 μ , p 0 μ ) are found to be functions of R 0 and of their We have used u −2 μ in (A 15a,b) and we have assumed that the terms u n μ , when increasing n > −2, decrease faster than u −2 μ ∝ 1/|x μ | 2 , which is consistent with non-diverging mass flux through a spherical cap. It follows that the main results provided by the mesoscopic problem are: the pressureP 0 (r, t) introduced in (A 14) is the pressure p 0 m at |x m | = 0 (at the mesoscopic scale, the bubbles are reduced to points), the velocity u 0 m is singular at the origin, and it encapsulates the effect of the oscillating bubble point (ifṘ 0 = 0, the bubble is not seen and u 0 m is regular). In dimensional form, these results contribute to the solution given in (2.8) at the dominant order; they will be complemented with the solution found at the following order.
The mesoscopic problem -using (A 17) in (A 5) along with (A 18a,b), we get that ∇ m p 0 m = 0. At this scale again, the pressure p 0 m is uniform as announced in (A 19a,b) resulting in (A 20a,b). It follows that the problem at the dominant order is set on and u 0 m is singular at x m = 0 from (A 20a,b). We have three unknown macroscopic fields: P 0 that we reach back from the microscopic scale (see (A 14)) and U 0± in (A 21). Note that we have assumed that U 0± do not depend on r m ; this will be justified in the forthcoming section.
(A 22a-c) At this scale, the bubbly screen is reduced to an interface but the boundary conditions when x → 0 ± are still missing. These are the conditions we are looking for and they are provided by matching the solutions of the mesoscopic and macroscopic problems. The intermediate region where both solutions are valid corresponds to x ∼ ± √ ε hence x m ∼ ±1/ √ ε. Hence, at the dominant order, the matchings read p 0 (0 ± , r, t) = p 0 m (r, t) = P 0 (r, t), u 0 (0 ± , r, t) = lim with the boundary conditions (A 23a,b) and the initial conditions p 0 (x, 0) = 0, u 0 (x, 0) = 0. From (A 23a,b), the pressure p 0 is continuous at x = 0. This is not the case for the normal velocity, and to show this result, we have to come back to the problem satisfied by u 0 m . The equation of incompressibility, div m u 0 m = 0 in (A 21), is integrated in the domain Y minus a ball B r of radius r (in order to avoid the singularity of u 0 m ). Using the boundary conditions in (A 21) and the singular behaviour of u 0 m in (A 20a,b), we get 0 = Y/B r div m u 0 m dx m from which u 0 x (0 + , r, t) − u 0 x (0 − , r, t) = ∂B r (r 3 eq /m)((R 0 ) 2Ṙ0 /|x m | 2 ) dx m . The integral over ∂B r gives a finite but non-zero contribution (independent of the radius of B r ), and eventually we obtain a non-trivial transmission condition Conclusion of the dominant order -at this stage the dominant order provides a model in which the bubbly screen is reduced to an effective screen across which the normal velocity experiences a jump dictated by R 0 , and R 0 satisfies (A 16) with P 0 (r, t) = p 0 (0, r, t). The result is identical to that obtained in Miksis & Ting (1989). We shall see that conducting the analysis at the following order makes the contribution of the bubble-bubble interaction appear.

A.3. A result at the following order: the bubble-bubble interaction
The analysis at order 1 is conducted following exactly the same steps as at order 0. For most of the results, we simply recover the same prediction in terms of Taylor expansions of the fields with respect to the small bubble radius, see e.g. (A 12a-c) and forthcoming (A 28a-c); this will be specified whenever it happens. The single but significant difference concerns the pressure seen by a single bubble at infinity. It is P 0 in (A 15a,b) and we have established that P 0 (r, t) = p 0 (0, r, t) the macroscopic pressure. It will beP 1 in (A 30) and we shall establish thatP 1 (r, t) = p 1 (0, r, t)+ a 'δ-contribution' due to the presence of the other bubbles, with δ defined by (A 34) and (A 36). This result is given in (A 49); it results in the modified Rayleigh-Plesset equation (A 52) whose dimensional form is (2.5).
Incidentally, solving at order 1 allows us to obtain the explicit form of the solution at the mesoscopic scale (A 33), which coincides with (2.8) in dimensional form.