Acoustic interaction between 3D-fabricated cubic bubbles.

Spherical bubbles are notoriously difficult to hold in specific arrangements in water and tend to dissolve over time. Here, using stereolithographic printing, we built an assembly of millimetric cubic frames overcoming these limitations. Indeed, each of these open frames holds an air bubble when immersed into water, resulting in bubbles that are stable for a long time and are still able to oscillate acoustically. Several bubbles can be placed in any wanted spatial arrangement, thanks to the fabrication process. We show that bubbles are coupled acoustically when disposed along lines, planes or in 3D arrangements, and that their collective resonance frequency is shifted to much lower values, especially for 3D arrangements where bubbles have a higher number of close neighbours. Considering that these cubic bubbles behave acoustically as spherical bubbles of the same volume, we develop a theory allowing one to predict the acoustical emission of any arbitrary group of bubbles, in agreement with experimental results.


Introduction
Acoustic metamaterials present extraordinary properties such as negative index of refraction, or enhanced absorption, giving the hope that it would be possible to create an invisibility cloak around an object cancel acoustic echoes from it. Typically, such metamaterials contain sub-wavelength resonators that give them unique properties in terms of effective density or compressibility. 1,2 Gas bubbles are good candidates for these sub-wavelength resonators, because of their remarkable resonance, explained by the much greater compressibility of the gas they contain compared to the surrounding liquid. A small amount of bubbles have a huge acoustic effect on sound propagation, 3 which can be experienced in every day life with the hot chocolate effect, 4 where tiny bubbles dramatically change the frequency of sound when the mug is lightly struck. Bubbles can be arranged in static configurations within a solid material, embedded in an aqueous gel 5 or silicone elastomer 6,7 with specific positions and sizes. Free gas bubbles in water can be trapped under a net, but with random positions. 8 Bubbles then give super absorption properties to these materials. 5,9 Here we would like to introduce metamaterials with gas bubbles in water in order to have as little viscous friction as possible while still being able to precisely control their sizes and positions. This would enable interesting underwater applications for providing surfaces that could be totally absorbant or acoustically transparent to ultrasonic waves. On a more fundamental level, it is of importance to characterize the resonance of groups of interacting bubbles.
Thanks to the recent evolution of 3D-printing techniques, we have shown in a previous work 10 that it is possible to (i) maintain the position of a gas bubble in water, and (ii) stabilise its size, using a single cubic frame where a bubble is trapped.
The present manuscript proposes to incorporate cubic bubbles as building blocks to be arranged into a metamaterial. Our purpose is to find out how bubbles couple acoustically to each other, and then to investigate what are the fundamental laws dictating the resonance as a function of the spatial arrangement of bubbles.

Methods
We built open cubic frames (external size 3 mm) supported by a loose scaffold (see Fig. 1). This structure was fabricated using a stereolithographic (SLA) 3D-printer (Titan from Kudo3D, ''Hard and Tough'' type resin, 50 mm resolution). The faces of the cube have square openings 1.94 AE 0.05 mm in height and 1.75 AE 0.05 mm in width (slightly smaller than the 2 mm in design). The actual dimensions of the cubic frames were determined by taking photographs of several structures using a macroscope (Leica Z6). The frames were silanized (30 minutes of vapourphase deposition of trichloro(perfluoro-octyle)silane) in order to make them hydrophobic.
Upon immersion in a water tank, bubbles spontaneously stay inside the cubic frames, because openings are small enough for capillary forces to be stronger than gravity and prevent water from entering. On the other hand, the spaces in the scaffold are large enough to let water invade the structure. In order to minimise the number of parasitic bubbles forming outside the cubes, the number of supporting frames was minimised as much as possible in the limit of two supports per cube. In particular, these frames were put in place so that no closed volumes are delimited, giving the ''X'' shapes that can be seen in Fig. 1 for complex networks.
Each bubble has flat interfaces located on the same plane than the external faces of the cube. The gas volume V g of these bubbles is equal to 17 mm 3 (using measured dimensions). As previously shown 10 they can, with a good approximation, be considered acoustically equivalent to a spherical bubble of radius R eq = 1.6 mm, with the same volume. Because the interfaces are not curved, the Laplace overpressure inside these bubbles is null and therefore bubbles are much more stable over time. We could study complex arrangements of bubbles, like the one shown in Fig. 1, for a long time before dissolution, which can occur after a day.
Experiments were performed in a tank (29 Â 29 Â 50 cm 3 ) made of PMMA filled with tap water and a small amount of bleach to prevent the development of microorganisms. Experiments were performed a few days after filling the tank in order for the dissolved gas to reach equilibrium. The sample, denoted (1) in Fig. 1, stood on a mesh cage placed on a steel plate, denoted (5) in Fig. 1, in order for the sample to be situated at the middle of the tank.
The acoustical response of such arrangements was measured by sending acoustic waves with an underwater speaker, broadcasting repeated frequency sweeps. Frequency sweeps ranging from 0.1 kHz to 5 kHz and lasting 1 second were generated by an arbitrary waveform generator (Handyscope HS5, TiePie) at a sampling frequency of 100 kHz. After being amplified (amplifier 7600M, Krohn Hite Corporation) they were sent to a waterproof loudspeaker (FR 13 WP, Visaton) denoted (2) in Fig. 1. The sound was measured by a hydrophone (8103, Brüel & Kjaer) denoted (3) in Fig. 1 and amplified (Nexus Conditioning Amplifier 2692, Brüel & Kjaer). This hydrophone can either be attached to a 3D moving stage, denoted (4) in Fig. 1, made from scavenged parts of a 3D printer (Prusa i3 -eMotionTech), allowing us to measure acoustic signals at different positions in the tank, or fixed at a given position for the whole experiment.
A signal V mic was recorded by placing a hydrophone in the vicinity of the structure. An additional recording V 0 mic was also performed at the same position, in the absence of the structure. Following Leroy et al., 11 we compute the relative bubble contribution using the normalised spectrum where Ṽ mic and Ṽ 0 mic are the Fourier transforms of the signals acquired and f is the frequency.

Acoustical response
The emission spectrum of a single bubble is shown in Fig. 2. It features a marked resonant behaviour with a maximum of the amplitude of the normalised emitted pressure |A| and a shift of its phase with respect to the external exciting field from 0 to p, meaning a shift from cos(j(A)) = 1 to cos(j(A)) = À1. Experimentally, we will define the resonant frequency by the condition cos(j(A)) = 0 at resonance. Here it is equal to 2050 Hz. This value is very close to the Minnaert's prediction 12 for a spherical bubble of equivalent gas  volume: f Minnaert = A/R eq = 2040 Hz with A = 3.24 m s À1 , the Minnaert constant and a radius R eq = 1.6 mm for a spherical bubble having the same gas volume as for the cubic bubble. The resonance frequency slowly increases with time (1.7% per hour), which we interpret as a slow dissolution of the bubble. These bubbles have a good quality factor, around Q = 20, meaning that the damping ratio d = 1/Q is around 5 Â 10 À2 .
In order to understand the coupling between bubbles, we started with two identical bubbles. One cubic bubble was placed at the sample area of the tank and the other was fixed to a 3D moving stage. Starting from a position of the stage where the bubbles are in contact, the distance d between the centers was increased by 0.5 mm steps. The hydrophone was placed at the same height as the bubbles, approximately 1 cm away from the fixed bubble and outside the path of the moving bubble. Fig. 3 (dots) shows the experimental values of the resonance frequency of such a system with the evolution of distance d: the resonance frequency decreases when bubbles are approaching each other, as is the case for spherical bubbles. 13 To further study the interactions between bubbles, we performed experiments with various numbers of bubbles arranged along a 1D line. We varied the number of bubbles from 2 to 8 bubbles, spaced by d = 8 mm. Fig. 4 shows that the experimental resonance frequency diminishes with the number of bubbles.
2D networks of identical bubbles were also studied. Bubbles were arranged in a matrix configuration, the closest neighbours spaced by 10 mm. Similar to lines of bubbles, it was found that the resonant frequency decreases with the number of bubbles present in the system (see Fig. 4). The frequency is more reduced with an increasing number of bubbles compared to the linear case because of the larger number of neighbours in 2D arrangements, which is highlighted in Fig. 4.
In 3D arrangements the number of neighbours increases dramatically. As every bubble will interact with more bubbles one can expect an even greater reduction in the resonance frequency compared to in 1D/2D geometries. In order to specifically pinpoint the influence of the dimensions, we have designed arrangements with a central bubble, surrounded by two neighbours in a line (1D line called 1 + 2), four neighbours in a plane (2D cross 1 + 4), and six neighbours in three directions (3D cross 1 + 6), see the photographs in Fig. 5a-c.
The table in Fig. 5 shows experimental observations up to a 3D arrangement of 3 Â 3 Â 3 bubbles. This last system has a resonant frequency of around 1 kHz, which represents a huge   reduction of the resonant frequency, which is 2 times smaller than the resonant frequency of a single bubble.

Modelling
We now attempt to model those observations.

Bubbles as spherical pulsators
We have previously shown 10 that the acoustic resonance frequency of a cubic bubble is very close to that of a spherical bubble of the same volume. In addition, the emitted field is spatially very close to that of a monopolar source, provided the distance d of the observation to the cube centre is larger than the cube size. For these reasons we model bubbles as spherical, and the sound emitted by each bubble is described as the one emitted by a spherical bubble whose radius evolves as R n (t) = R n0 + a n exp(iot) where R n0 is the bubble radius at the rest of the bubble n and a n is the amplitude of pulsation (|a n | { R n0 ), in response to a driving acoustic pressure P ac = P a exp(iot). The mechanical equations for the fluid are written in Appendix A. For simplicity we have included the full derivation in the appendix and we give only the end results here. We found that each bubble, when alone, behaves as an oscillator with a resonant pulsation frequency o n and damping ratio d n (given by eqn (25) and (26)). Note that the resonant frequency is well approximated by the Minnaert equation f n ¼ o n =2p ¼ A=R n for millimetric bubbles in water, with a = 3.24 m s À1 . Damping is the consequence of three phenomena: thermal effects, viscosity and radiation.

Interactions between pulsators
Each bubble receives the pressure emitted by neighbouring bubbles in addition to the incoming driving field. The model predicts that bubble oscillators when oscillating with an amplitude a n are coupled by the following set of equations: see Appendix A for a full derivation. The coupling term is the second term on the left hand side. In this equation the distance between the centres of bubbles n and m is noted as d nm . The wavevector is k = o/c, with c being the speed of sound and neglecting here fluid viscosity. The acoustic pressure received by the bubbles is assumed to be uniform and equal to P a , and r 0 is the water density. Inverting this system of equations provides the value of the amplitudes a n . The pressure emitted by all these bubbles is p bubbles ðr; tÞ ¼ Àr 0 c 2 k 2 expðiotÞ X N n¼1 R n0 2 a n exp Àikd rn ½ d rn (4) with d rn = |r À r n |, the distance between the microphone of position vector r and the nth bubble of position vector r n .

Numerical predictions of the acoustic response
Here we assume all bubbles to have the same rest radius R n = R 0 , and we will adjust this parameter to describe the experiments. The total acoustic pressure amplitude field received by a microphone is p mic = P ac + p bubbles , while it is p 0 mic = P ac in the absence of bubbles. We can thus model the measured relative spectrum as: This prediction always showed a peak in amplitude and a change of phase from 0 to p. As in the experiments, we define the numerical resonant frequency f num 0 as the first frequency for which the phase of the response crosses 0, which is the condition cos(j(A num )) = 0 at resonance.

Comparison with experiments
One bubble. The model predicts a spectral response A num ( f ) for single bubbles, in agreement with one bubble experiments when choosing the parameter R 0 = 1.57 mm (Fig. 2) close to the expected value from the gas volume.
Two bubbles. We correctly predict the resonance of a couple of bubbles (Fig. 3), as a function of distance between the bubble centers, with a slightly larger R 0 = 1.65 mm.
In this special case, it is possible to give a simple analytical prediction (see details in Appendix B, eqn (30)) and it gives a good agreement with experimental and numerical data as it is shown in Fig. 3 that is classical in the litterature. 14 Arrangements in lines or matrices. For 1D lines of bubbles and 2D networks (from 2 Â 2 to 6 Â 6) the numerical predictions are in good agreement with experiments, still taking R 0 = 1.6 mm, see Fig. 4.
Arrangements in volumes. We also find a good agreement for 3D grids, see the table in Fig. 5. Note that because some of these systems have a lot of symmetries, their resonant frequency can be expressed analytically, assuming two groups of bubbles oscillating with the same amplitude: the central one, and the peripheral ones, see Appendix B for formulas.

Conclusion
In conclusion, we have shown that the interactions which take place between the bubbles downshift the global resonant frequency of the system and can be predicted using a model of spherical bubbles. This interaction has tremendous effects in 3D arrangements, even when bubbles are parted by distances several times their own sizes (here d/R 0 = 5). We found that the radius of the spherical bubbles used in the simulations varied between 1.57 mm and 1.65 mm depending on the experiments, which we believe is due to imperfections in the fabrication process that slightly changes the inner gas volume, or to a different attachment of the contact line during immersion.
Perspectives will be the global study of larger scale metamaterials, with different inter-bubble distances, in order to understand the transmission and absorption properties. In addition further work is needed to detect higher order modes where not all bubbles oscillate in phase.

Conflicts of interest
There are no conflict of interest to declare.

A Theory for an assembly of bubbles
We have previously shown 10 that the acoustic resonance frequency of a cubic bubble is very close to that of a spherical bubble of the same volume. In addition, the emitted field is spatially very close to that of a monopolar source, provided the distance d of the observation to the cube centre is larger than the cube size. For these reasons, we will model bubbles as spherical, which should be valid if they are a sufficient distance apart.
We consider a cluster of spherical bubbles located arbitrarily in space. We introduce a global Cartesian coordinate system. The position vector of an arbitrary space point is denoted by r and has coordinates (x, y, z). The position vector of the centre of the nth bubble is denoted by r n and has coordinates (x n , y n , z n ).
The time-varying radius of the nth bubble is represented as R n (t) = R n0 + a n exp(iot) (6) where R n0 is the bubble radius at rest and a n is the amplitude of the bubble pulsation. We assume that |a n | { R n0 and solve the problem in the linear approximation. The liquid around the bubbles is assumed to be viscous and compressible. In the linear approximation, the liquid has a velocity v and a perturbed density r, which obey the following equations: 15 where r 0 is the equilibrium liquid density, p is the perturbed liquid pressure, and Z and z are the shear viscosity and the bulk viscosity, respectively, while c is the speed of sound. These equations are the linearized versions of the compressible Navier-Stokes equation, the continuity equation and the equation of state of the liquid. The liquid motion is assumed to be irrotational with v = rf, where f is the velocity potential. With a time dependence proportional to exp(iot), eqn (7)-(9) then provide the Helmoltz equation where k is given by The velocity potential around the nth bubble is a spherically symmetrical solution of eqn (10), f n ðr; tÞ ¼ A n jr À r n j exp iot À ikjr À r n j ð Þ : The velocity is then v n ðr; tÞ ¼ À A n 1 þ ikjr À r n j ð Þr À r n ð Þ jr À r n j 3 Â exp iot À ikjr À r n j ð Þ (13) and the pressure produced by the bubble is p n ðr; tÞ ¼ À ir 0 c 2 k 2 A n ojr À r n j exp iot À ikjr À r n j ð Þ The influence of bubbles on the pulsation of each other is taken into account with an accuracy up to leading terms with respect to interbubble distances and compressibility effects. Within the framework of this accuracy, from the boundary condition for the liquid velocity at the surface of the nth bubble, one obtains To find a n , we apply the boundary condition for the normal stress at the surface of the nth bubble, which is given by s m ð Þ jrÀrnj¼Rn (16) where P gn is the equilibrium pressure of the gas inside the nth bubble, s is the surface tension, P 0 is the equilibrium pressure in the liquid, P ac = P a exp(iot) is the driving acoustic pressure, N is the number of bubbles and s m (r,t) is the normal stress produced by the mth bubble in the liquid. Note that in eqn (16), s m (r,t) is taken at the surface of the nth bubble. In order to take into account deviations from the adiabatic law, the exponent k n is defined as 16 k n = g(a n + ib n ) (17) where g is the specific heat ratio of the gas. The quantities a n and b n (which are introduced to describe the phase shift between variations of the gas pressure and the bubble volume, caused by heat losses) are the real and imaginary part of k n and are calculated by where w n ¼ 3ðg À 1Þ Â X n sinh X n þ sin X n ð Þ À 2 cosh X n À cos X n ð Þ X n 2 cosh X n À cos X n ð Þ þ 3ðg À 1ÞX n sinh X n À sin X n ð Þ (20) , D gn = K/r n c pg is the thermal diffusivity of the gas inside the nth bubble, K is the thermal gas conductivity, c pg is the specific gas heat at constant pressure, r n = r A P gn /P A is the equilibrium gas density inside the nth bubble and r A is the gas density at the atmospheric pressure P A . The complex exponent k n in eqn (16) makes it possible to take into account of thermal effects whose contribution to the damping of bubble oscillations is known to dominate the viscous and radiation contributions over a wide range of bubble radii -from a few microns to several hundred microns -unless the driving frequency is considerably above the bubble monopolar resonance frequency. [16][17][18] The normal stress produced by the nth bubble is calculated by where v n = |v n |. Substitution of eqn (13)-(15) into eqn (21) yields s n ðr; tÞ ¼ À ioR n0 2 a n ir 0 o jr À r n j þ 4ikZ jr À r n j 2 þ 4Z jr À r n j 3 Â exp iot À ikjr À r n j ð Þ Substituting eqn (6) and (22) into eqn (16), one obtains which is eqn (3) in the main text, where we have introduced the resonance frequency of bubbles (when isolated) ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3ga n P gn r 0 À 2s r 0 P n0 s (25) their damping constant and the inter-bubble distances d nm = |r n À r m |. It should be noted that we assume the bubble radii to be much smaller than the acoustic wavelength, kR n { 1. Therefore, terms of this order are neglected in eqn (24), whereas the distances between the bubbles, d nm , can be comparable to and even greater than the acoustic wavelength. Thus, terms of the order kd nm are kept in eqn (24). We also assume the bubble radii to be small compared to the distances between the bubbles, R n0 { d nm . Therefore, terms of a higher order than R n0 /d nm are omitted in eqn (24).
The system of eqn (24) is a system of N algebraic equations in the unknown a n . Its solutions give the values of a n . The total pressure produced by all bubbles at the point r is calculated by pðx; y; z; tÞ ¼ Àr 0 c 2 k 2 expðiotÞ P N n¼1 R n0 2 a n exp Àikd rn ½ d rn (27) with d rn = |r À r n |, which is eqn (4) in the main text.
B Analytical predictions for the resonance frequency of 1 + N bubbles For the specific case of a bubble surrounded by N bubbles at the same distance d, it is possible to find analytical predictions for the resonance frequency, still assuming spherical bubbles. We have tested several configurations. B.1 Two bubbles: 1 + 1. Note that for the case of a couple of bubbles we can derive an analytical formula for the resonance frequency. Eqn (3) gives, in a matrix form: where we have introduced o n = o 0 and d n = d. The solution is: This suggests a natural resonance at B.2 Three aligned bubbles: 1 + 2. Now, for simplicity we neglect damping factors d, and we have assumed the distance d to be extremely small compared to the wavelength kd { 1 and e Àikd C 1.
Owing to the symmetry of the network of three aligned bubbles, we can assume that the amplitudes of the first and last bubbles are equal (a 1 = a 3 ), while the amplitude of vibration of the central bubble (a 2 ) might be different. Eqn (24) gives The eigenvalues of the matrix on the left-hand side are the frequency o ana 0 = o À and corresponds to the mode where all bubbles oscillate in phase, while o + is a higher frequency mode that is not excited here.