# Collimated beam formation in 3D acoustic sonic crystals

A L Vanel, M Dubois, Clément Tronche, S Fu, Y-T Wang, G Dupont, A D Rakić, K Bertling, R Abdeddaim, S Enoch, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

A L Vanel, M Dubois, Clément Tronche, S Fu, Y-T Wang, et al., Collimated beam formation in 3D acoustic sonic crystals. New Journal of Physics, 2024, 26 (7), pp.073021. 10.1088/1367-2630/ad5c94. hal-04912298

HAL Id: hal-04912298

https://amu.hal.science/hal-04912298v1

Submitted on 8 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Collimated beam formation in 3D acoustic sonic crystals

To cite this article: A L Vanel et al 2024 *New J. Phys.* **26** 073021

View the [article online](#) for updates and enhancements.

You may also like

- The mean velocity of nylon spheres transported in a horizontal water pipe
  A M Binnie and T-C Kuo

- Manipulation and welding of metal spheres above 10 m using needle-like probe
  Takeshi Konno, Mikihiko Kobayashi, Mitsuru Egashira et al.

- [Foreword](#)
# New Journal of Physics

The open access journal at the forefront of physics

Deutsche Physikalische Gesellschaft DPG

IOP Institute of Physics

Published in partnership
with: Deutsche Physikalische
Gesellschaft and the Institute
of Physics

<figure><img src="image_3.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

21 March 2024

REVISED

14 June 2024

ACCEPTED FOR PUBLICATION

27 June 2024

PUBLISHED

12 July 2024

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_4.png" /></figure>

## PAPER

### Collimated beam formation in 3D acoustic sonic crystals

A L Vanel<fnref n="1" />,*, M Dubois<fnref n="2" />,<fnref n="3" /> ID, C Tronche<fnref n="4" /> ID, S Fu<fnref n="5" />, Y-T Wang<fnref n="6" />,<fnref n="7" />,<fnref n="13" /> ID, G Dupont<fnref n="8" />, A D Rakić<fnref n="9" /> ID, K Bertling<fnref n="9" /> ID, R Abdeddaim<fnref n="2" />, S Enoch<fnref n="2" /> ID, R V Craster<fnref n="6" />,<fnref n="10" /> ID, G Li<fnref n="11" /> ID, S Guenneau<fnref n="10" />,<fnref n="12" /> ID and J Perchoux<fnref n="4" /> ID

**E-mail:** alice.lucie.vanel@cern.ch

**Keywords:** phononic crystal, collimation, acoustic, endoscope

Supplementary material for this article is available [online](#)

### Abstract

We demonstrate strongly collimated beam formation, at audible frequencies, in a three-dimensional acoustic phononic crystal where the wavelength is commensurate with the crystal elements; the crystal is a seemingly simple rectangular cuboid constructed from closely-spaced spheres, and yet demonstrates rich wave phenomena acting as a canonical three-dimensional metamaterial. We employ theory, numerical simulation and experiments to design and interpret this collimated beam phenomenon and use a crystal consisting of a finite rectangular cuboid array of 4 × 10 × 10 polymer spheres 1.38 cm in diameter in air, arranged in a primitive cubic cell with the centre-to-centre spacing of the spheres, i.e. the pitch, as 1.5 cm. Collimation effects are observed in the time domain for chirps with central frequencies at 14.2 kHz and 18 kHz, and we deployed a laser feedback interferometer or Self-Mixing Interferometer – a recently proposed technique to observe complex acoustic fields—that enables experimental visualisation of the pressure field both within the crystal and outside of the crystal. Numerical exploration using a higher-order multi-scale finite element method designed for the rapid and detailed simulation of 3D wave physics further confirms these collimation effects and cross-validates with the experiments. Interpretation follows using High Frequency Homogenization and Bloch analysis whereby the different origin of the collimation at these two frequencies is revealed by markedly different isofrequency surfaces of the sonic crystal.

## 1. Introduction

A challenge for acoustics, in particular for phononic crystals, is to generate exotic wave phenomena at low-frequencies, i.e. at audible frequencies, for sound propagation in the human hearing range; the ability to achieve precise control would enable the design of novel devices. We address this challenge motivated from the two-dimensional analysis of simple closely-spaced two-dimensional crystal arrays which operate as subwavelength tuneable metamaterials [1] enabling the cut-off frequency for the acoustic branch to be lowered to pull the entire acoustic branch subwavelength; this enables exotic wave phenomena at low frequencies. Two-dimensional models are a convenient abstraction, but this behaviour needs to carry

<footnote n="1">¹ CERN, 1211 Geneva 23, Switzerland</footnote>

<footnote n="2">² Aix-Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France</footnote>

<footnote n="3">³ Multiwave Imaging SAS, Marseille, France</footnote>

<footnote n="4">⁴ LAAS-CNRS, Université de Toulouse, CNRS, INP, Toulouse, France</footnote>

<footnote n="5">⁵ Eastern Institute for Advanced Study, NingBo, People's Republic of China</footnote>

<footnote n="6">⁶ Department of Mathematics, Imperial College London, London SW7 2AZ, United Kingdom</footnote>

<footnote n="7">⁷ Photonics Initiative, Advanced Science Research Center, City University of New York, New York, NY 10031, United States of America</footnote>

<footnote n="8">⁸ Aix-Marseille Univ, CNRS, Centrale Marseille, IRPHE, Marseille, France</footnote>

<footnote n="9">⁹ School of Electrical Engineering and Computer Science, The University of Queensland, Brisbane 4072, Australia</footnote>

<footnote n="10">¹⁰ UMI 2004 Abraham de Moivre-CNRS, Imperial College London, London SW7 2AZ, United Kingdom</footnote>

<footnote n="11">¹¹ Department of Mathematics, The University of Hong Kong, Pokfulam Road, Hong Kong, Hong Kong Special Administrative Region of China, People's Republic of China</footnote>

<footnote n="12">¹² The Blackett Laboratory, Physics Department, Imperial College London, London SW7 2AZ, United Kingdom</footnote>

<footnote n="13">¹³ Department of Photonics, National Sun Yat-sen University, Kaohsiung 804201, Taiwan</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
through to full three-dimensions and be verified both numerically and experimentally before moving into device design; in this paper we investigate collimated beam formation for arrays of closely-packed three dimensional spheres in acoustics.

In two-dimensions a closely packed array of cylinders creates a network of coupled Helmholtz-like resonators with large voids coupled to each other by narrow gaps and the asymptotic approach in [1] generates a discrete mass-spring network analogue. Moving this to three dimensions is non-trivial: Numerical simulations of meaningfully large finite arrays are hampered by the need to accurately resolve the extremely thin air-filled regions between the spheres leading to numerical issues that are overcome here. Experiments are also challenging as, until recently, it was not possible to image the sound field within a crystal, and lastly, the asymptotic matching methods [1] employed are no longer valid since there are regions where the spheres nearly touch, but equally other regions where the acoustic medium is not as constrained as in two-dimensions.

We comprehensively investigate a crystal constructed by closely packed spheres using the following approaches: We use a specially designed numerical scheme to overcome the computational issues [2]. An experimental optical interferometry approach [3, 4], developed by some of the team, to enable the imaging, and we develop the physical modelling and intuition of these closely spaced media to explain the observed phenomena. Combining these approaches leads to a coherent design capability that we use to create collimated beams. Consequently, closely-packing a simple phononic crystal makes metamaterial behaviours appear, which means each element within the crystal is subwavelength and yet critically affects the macroscale behaviour.

There is a vast literature on acoustic metamaterials [5–7], including analysis of their 3D effective properties [8–10] and the fabrication of three dimensional mechanical [11] and acoustic [12] metamaterials with complex geometry. There is, however, a richness of wave phenomena in simple three-dimensional crystals that has not yet been fully explored, and exploited; examples of the behaviour possible include acoustic sonic crystals designed to achieve a flat lens effect in waterborne acoustics [13] and airborne acoustics [14]. For the collimated beam formation in acoustics we draw upon the extensive literature, primarily in two dimensions and in electromagnetism, using iso-frequency contours to predict extreme dynamic anisotropy; polygonal contours leading to self-collimation [15–19]. Self-collimation is a powerful concept, readily observed in 3D for discrete mass-spring systems [20], but apart from some notable examples in electromagnetism, i.e. [21, 22], its use in three dimensional continuum models has been limited. The richness of behaviours possible through dynamic anisotropy is also illustrated in [23] for electromagnetic crystals notably cubic arrays of split ring resonators and also [24–26] for dynamic anisotropy in 2D dielectric photonic crystals. We demonstrate that these approaches can be taken into acoustics and, for sources exterior to a crystal, a lattice medium generates wide or narrow acoustic beams. By doing so we highlight the versatility and tunability that is achievable using a basic phononic structure, without requiring significant complexity, and bring numerical, experimental and theoretical analysis to bear to provide a comprehensive analysis of the collimation effects created.

The plan of the article is as follows: section 2 presents the experimental setup and measurement based on the Self-Mixing Interferometry (SMI) sensing methodology [3, 27]. Section 3 is devoted to the higher-order Multi-scale Finite Element Method (HMsFEM) [2, 28] used to numerically solve the scattering problem. Results obtained are compared to experiments and standard Finite Element Method (FEM) references. Section 4 shows that High Frequency Homogenization theory (HFH) allow us to propose physical interpretations of the collimation effect at different frequencies. Unlike most of the studies that plot dispersion relations around the edges of the Brillouin zone [29, 30], we analyse the Fermi surfaces to get a fuller picture of the wave phenomena responsible for the observed self-collimation effects at 14.2 and 18 kHz. Finally some concluding remarks are drawn together in section 5.

## 2. Experiments

The finite crystal used here is a rectangular cuboid array of 400 (4×10×10) 3D-printed polymer spheres 1.38 cm in diameter with a centre-to-centre spacing $a = 1.5$ cm and is shown in figure 1. A commercial loudspeaker (Visaton SC5) is set 16 mm from the crystal—which corresponds roughly to the wavelength of the acoustic wave—and facing the thinner facet of the crystal (see figure 1). It is driven with square pulses of sinusoidal signals with frequencies centered at 14 kHz and 18 kHz. The pulses have a duration of 12 periods of the sinus and a spacing in time of 21 ms thus avoiding echoes appearing in the measurements while producing a relatively large frequency band signal. To visualize the acoustic propagation in a non-intrusive manner we use a methodology based on the observation of the refractive index change in the air induced by pressure changes [14]. This acousto-optic effect is captured through a compact interferometric sensor designed for this measurement with the principle first described in [3]. It consists of a commercial laser
<figure><img src="image_5.png" /><figcaption>Figure 1. Experimental setup: the crystal and the associated acoustic source are mounted on a rigid structure built with metallic rods. The structure is moved along X2 and X3 axis using motorized translation stages while the sensor and the distant reflector remain in a fixed position. (a) Top-view. (b) Side view. (c) Photograph of the setup.</figcaption></figure>

diode (Thorlabs L1550P5DFB) with single-mode emissions at $\lambda = 1550$ nm that points at a reflector so that a part of the emitted light re-enters the laser's cavity producing interferences between the inner optical wave and the back-reflected one. This interferometric approach is known as SMI [27, 31] and it is sensitive to any changes in the time of flight in the external cavity that realizes the laser and the reflector. In the context of acoustic sensing, this change is due to the acousto-optic effect that links the refractive index of air to the pressure among other parameters [32]. With pressure levels encountered in the acoustic domain, the relationship change in refractive index with pressure can be considered as perfectly linear. As expressed in [3] the laser emitted power variations $P(t)$ and the refractive index ones $\delta n(t)$ along the laser propagation axis $z$ are linked as

$$P(t) = P_0 \cos \left( \frac{4\pi\nu}{c} \int_{0}^{L_{\text{ext}}} \delta n(t, z) dz + \Psi \right), \quad (1)$$

where $P_0$ is the average emitted power, $\nu$ is the laser frequency, $L_{\text{ext}}$ is the length of the external cavity and $\Psi$ is a constant phase term. A part of the light emitted by the laser diode is collected using the monitoring photodiode that is included in the laser package and the photogenerated current is then amplified and converted thanks to a custom-made amplifier to a voltage variation that is acquired using a National Instrument acquisition card (USB-6361). This terminal voltage is then an image of the integration of the pressure field along the laser axis in function of time. Images of the acoustic wave propagation come from the sensor (laser + reflector) which is translated with a fixed step of 1.5 mm on both $x$ and $y$ axes; each position of the sensor represents a pixel of the image shown in figures 2(c) and 3(c). The experiment is repeated for each laser position in order to represent the integrated pressure field in two-dimensional maps.

For the temporal results shown in figures 2(a), (b) and 3(a), (b) the experimental data is temporally downsampled by a fourfold factor leading to a time resolution of 4 µs. It is then filtered spectrally with a Gaussian window to obtain the acoustic response with chosen central frequency and bandwidth. The Gaussian window presents a full width at half maximum of 5.25 kHz and is centered in 14.5 kHz and 17.7 kHz as shown in figures 2(a) and 3(a).

# 3. Simulation using a higher-order multi-scale finite element method

## 3.1. Mathematical formulation

For the modelling we use the (non-dimensionalised) Helmholtz equation that holds for time-harmonic acoustic waves:

$$\frac{\partial^2 u}{\partial x_1^2} + \frac{\partial^2 u}{\partial x_2^2} + \frac{\partial^2 u}{\partial x_3^2} + \Omega^2 u = 0, \quad (2)$$

for $u(\mathbf{x}) = u(x_1, x_2, x_3)$ with $r = \sqrt{x_1^2 + x_2^2 + x_3^2} > s$ (i.e. everywhere outside of the spheres in the absence of a source) with $\Omega = \omega a/c$, with $a$ the cubic array pitch (in unit of meters) and where $\omega$ is the angular frequency (in unit of rad s⁻¹) and $c = 340$ m s⁻¹. Consistent with the acoustically sound-hard boundary conditions on
<figure><img src="image_6.png" /><figcaption>Figure 2. Experimental observation at 14 kHz showing endoscopic effect: (a) spectrum of the transmitted signal for an excitation centered at 14 kHz; the dotted line denotes the Gaussian window used. (b) Transmitted time trace in solid red with the signal from the input side shown in solid blue. (c) Snapshot of the amplitude of integrated pressure at time t = 1.72 ms (see movie at [Multimedia 1]).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3. Experimental observation at 18 kHz of the endoscopic effect: (a) Spectrum of the transmitted signal for an excitation centered at 18 kHz; the dotted line denotes the Gaussian window used. (b) Transmitted time trace in solid red with the signal from the input side is shown in solid blue. (c) Snapshot of the amplitude of the integrated pressure at time t = 0.88 ms (see movie at [Multimedia 2]).</figcaption></figure>

each sphere, we take Neumann boundary conditions, the normal derivative of the function being zero, to hold on the surface of each sphere

$$ \frac{\partial u}{\partial \mathbf{n}} |_{r=s} = 0, \qquad (3) $$

where **n** denotes the outwards pointing vector to the spheres. Section 4 also presents the Bloch eigenvalue problem (2)–(3) posed in a periodic cell with Floquet-Bloch conditions. Here, we use a standard finite element scheme implemented in the Comsol Multiphysics package, since the computational domain is just a cubic cell containing a single sphere.
## 3.2. Method

We use the HMsFEM [2, 28] to rapidly solve the three-dimensional Helmholtz equation for the finite array excited by a source. The underlying methodology is to solve the Helmholtz equation in a relatively coarse grid and exploit carefully constructed multi-scale basis functions that capture local heterogeneity information and wavenumber information; [2, 28] compare and contrast HMsFEM with classical polynomial-based FEM highlighting the advantages in terms of speed and accuracy. A major advantage of HMsFEM is that it significantly reduces the degrees of freedom of the FEM defined on the fine scale mesh and thereby saves both a huge amount of CPU time and memory. The construction of the multiscale basis in a coarse element, K, for two types of multiscale basis associated with the boundary and interior interpolation points of the polynomial basis functions is the key building block required for HMsFEM. The first type is obtained by solving

$$ \Delta\phi_i + \Omega^2\phi_i = 0, \quad \text{in } K, \qquad (4) $$

$$ \phi_i = \Phi_i, \quad \text{on } \partial K, \qquad (5) $$

in which $\Phi_i$ (a Dirichlet datum) is the polynomial shape function corresponding to the $i^{th}$ interpolation point that lies on the boundary $\partial K$ of K. For the second type of multi-scale basis, one needs to solve

$$ \Delta\phi_j + \Omega^2\phi_j = \Phi_j, \quad \text{in } K, \qquad (6) $$

$$ \phi_j = 0, \quad \text{on } \partial K, \qquad (7) $$

where $\Phi_j$ (now a volumetric forcing) is the interior polynomial shape function associated with the $j^{th}$ interpolation point that lies inside K.

The power of the HMsFEM approach is that we replace each standard polynomial basis function in the classical FEM with a multiscale basis function constructed by solving a Helmholtz equation dependent upon the local problem. Since the basis functions are constructed in a local domain, the computational cost is affordable. Importantly, these multiscale bases are defined in a very coarse mesh, therefore, the degrees of freedom of the resulting linear system with HMsFEM is greatly reduced compared to a classical FEM. This enables us to solve the three-dimensional sphere array, with close spacing, both accurately and rapidly. In our numerical simulations, the size of each coarse element (1.5 × 1.5 × 1.5 cm) is set to be exactly the same as the cubic cell that contains only one polymer sphere and the computational domain consists of 20 × 10 × 20 coarse elements. Therefore the size of the simulation domain is 30 × 15 × 30 cm. Here, for simplicity, we use cubic elements to approximate the spheres and the length of the fine-scale cubic element is 1/60th the length of the coarse-scale element. Perfectly matched layers [33] with a thickness of one coarse element are implemented to absorb outgoing waves on all sides of the computational domain. The dimension of the final linear system to be solved is merely 119 629 after utilizing the static condensation algorithm [34] and the CPU time for solving this linear system is just under a minute using Matlab. The accuracy for all numerical tests are within a 5% relative error, with the reference solution calculated by FEM on an extremely fine mesh as the reference solution. Hence, we can avoid mesh independent study. We emphasize that it is prohibitive to simulate using FEMs directly since they have to be solved over the fine-scale mesh that resolves all the interfaces. Consequently, this results in a very large linear system. In this case, even if we use just first-order FEM subject to the fine-mesh to perform the numerical simulations, one needs to solve linear systems whose degrees of freedom add up to over 10 millions. Furthermore, this large linear system usually exhibits a very large condition number [35], making its resolution extremely expensive.

## 3.3. Results

The results of both the physical experiments and the HMsFEM numerical simulations for 14.2 and 18 kHz are shown in figures 4 and 5, respectively; the left columns correspond to the experiments, whereas the right ones show the numerics. The upper rows correspond to the top view of the crystal ($x_1x_2$ plane), the lower row to the side view ($x_2x_3$ plane): there is remarkable agreement between the numerical simulations and the experiments. Both frequencies demonstrate distinctive collimation effects, at 14.2 kHz we see a plane wave exit the crystal that is roughly the width of the crystal, whereas at 18 kHz we observe a focused beam exit the crystal. The fidelity of the field maps shown in figures 4 and 5 also confirm a negligible perturbation from the thin rods supporting the spheres in the experimental setup that are not included in the simulation model.

# 4. Periodic acoustic wave problem

We now move to investigate the band spectrum of the sonic crystal since these allow us to look at iso-frequency contours and surfaces and hence gain physical understanding and interpretation of the collimation shown in figures 4 and 5.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Comparison between experiments (left panels) and numerical simulations (right panels) at 14.2 kHz from top (upper panels) and side (lower panels) view.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Comparison between experiments (left panels) and numerical simulations (right panels) at 18 kHz from top (upper panels) and side (lower panels) view.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 6. Bloch dispersion curves around the edges of the irreducible Brillouin zone (shown inset) for a cubic array of rigid spheres of diameter 1.38 cm with periodicity 1.5 cm. The global minimum of the second (optical) band at point X corresponds to the frequency of the source in figures 2 and 4 which is plotted as the dash-dotted horizontal blue line at frequency f = 14.2 kHz. The frequency of the source in figures 3 and 5 is drawn as a dash-dotted horizontal purple line at frequency f = 18 kHz. The red dashed line show HFH approximations around the X point from equation (10). The yellow dashed lines represent free-space propagation.</figcaption></figure>

For acoustic waves propagating through an infinite, perfect, triply periodic medium, one can invoke the Floquet-Bloch theorem [29, 30] and consider the cubic cell, in non-dimensional coordinates, −1 < x₁, x₂, x₃ < 1 with quasi-periodic conditions applied to the faces. The quasi-periodic Bloch boundary conditions are given as

$$u(\mathbf{x} + \mathbf{d}_i) = e^{i\kappa \cdot \mathbf{d}_i} u(\mathbf{x}), \qquad (8)$$

where dᵢ is the lattice vector along xᵢ, i = 1, 2, 3. For instance d₁ = (2, 0, 0). Equation (8) also involves the Bloch wave-vector κ = (κ₁, κ₂, κ₃) characterizing the phase-shift as one moves from one cell to the next.

This Bloch problem is solved numerically and the dispersion relations that link the frequency ω and Bloch wavenumber |κ| are deduced; as is conventional in solid state physics [29] only a limited range of wavenumbers need be considered in order to detect band-gaps, namely the wavenumbers along the right-angled tetrahedron ΓXMR shown in the irreducible Brillouin zone (IBZ) in figure 6; there are some exceptions to this limitation such as operators on graphs [36]. The vertex coordinates are Γ = (0, 0, 0), X = (π/a, 0, 0), M = (π/a, π/a, 0) and R = (π/a, π/a, π/a).

The dispersion curves for a triply periodic array of rigid spheres, shown in figure 6, illustrate several interesting features: a partial stop-band where wave propagation is disallowed along ΓX and XM between 13.2 and 14.2 kHz, coalescing bands at high-symmetry points M, X and R and regions of flat dispersion curves where the group velocity is zero and features of slow sound occur [37]. One can see that the group velocity is vanishing at the second band near the high-symmetry point X in figure 6; the frequency 14.2 kHz indicated by the dash-dotted blue line should thus be associated with an interesting wave phenomenon. Excitation at, or very close, to the frequency predicted would lead to oscillations that resemble a standing wave as it seems to be the case in figure 4; crucially this standing wave has directionality that can be identified from the wavenumber description in the Brillouin zone shown in figure 6.

To further confirm the physical mechanism leading to the collimation effect experimentally and numerically observed at 14.2 kHz, we propose to replace the sonic crystal by an effective medium in the high frequency regime. The property of this effective medium should depend both on the frequency and the Bloch wavenumber, as the phenomenon occurs well beyond the quasi-static limit.

## 4.1. Interpretation of collimated beam forming at 14.2 kHz with high frequency homogenization

The collimation effect observed at 14.2 kHz is related to strong anisotropy as shown by the isofrequency surface of figure 7 which is wrapped around X point with strong ellipticity in κ₁ and κ₃ directions. The slight mismatch of frequency to 14.2 kHz for the infinite periodic crystal is attributed the finite size of the crystal used in the experiment and simulations. This leads to a plane wave like behavior in the crystal and results in the wave leaving the crystal with a flat phase profile. This is, in effect, a collimation experiment where the
<figure><img src="image_11.png" /></figure>

compact source is transformed into a plane wave similar to the Luneburg lens based on index gradient media [38–47]. The collimation effect is related to the excitation of a single Bloch mode of the crystal leading to a plane wave like propagation in the crystal as illustrated in figure 4. This is further confirmed by two experimental observations: (i) we observe a clear single frequency contribution at the exit of the crystal in figure 2(a) and (ii) the group velocity is significantly reduced as observed in figure 2(b) with the pulse delay. The movie of the experiment [Multimedia 1] is instructive in showing the combined effect of the low group velocity and the resultant delay allowing the wavefield to spread across the crystal.

In figure 6, we approximate the dispersion curves computed via Finite Elements (Comsol Multiphysics) with asymptotic curves originated in the High-Frequency Homogenization theory (HFH) developed in [48]. The multiple scales used for HFH are a short-scale $\xi_i = x_i/l$ and a long-scale $X_i = x_i/L$ for $i = 1, 2, 3$ where $l$ and $L$ represent respectively the characteristic small scale (half length of a cell) and the long scale. A small parameter is formed as $\epsilon := l/L$ and expansions of the frequency $\Omega^2 = \Omega_0^2 + \epsilon\Omega_1^2 + ...$ and the solution $u = u_0 + \epsilon u_1 + ...$ are substituted in (2). Here $\Omega_0$ is a standing wave frequency given at the high-symmetry points X, $\Gamma$, M and R of the irreducible Brillouin zone and a perturbation scheme can be developed about the band edges. The leading-order term for $u$ is $u_0(\xi, X) = f_0(X)U_0(\xi)$ where the function $f_0$ representing an envelope of the solution $u_0$ is obtained by an equation only on the long scale: The theory is detailed in [48]. Changing back to the original $x_i$ coordinates the effective medium equation reads,

$$T_{ij} \frac{\partial^2 f_0}{\partial x_i \partial x_j} - \frac{(\Omega^2 - \Omega_0^2)}{l^2} f_0 = 0. \qquad (9)$$

$T_{ij}$'s are integrated quantities of the leading order and first order short scale solutions with $T_{ij} = 0$ for $i \neq j$ in the present illustrations. As in [48], assuming Bloch waves, the asymptotic dispersion relation for $\Omega$ reads,

$$\Omega \sim \Omega_0 - \frac{T_{ij}}{2\Omega_0} \kappa_i \kappa_j, \qquad (10)$$

where $\kappa_i = K_i - d_i$ and $d_i = 0, -\pi/2, \pi/2$ depending on the band edge in the Brillouin zone about which the asymptotic expansion originates; these asymptotics give the red dashed curves in figure 6 and the red surface in figure 7. For the case of multiple modes originating from the same point, as is the case for the second and third bands at point M and the first and second bands at point R, equation (10) is no longer valid and one obtains two coupled equations for $f_0^{(i)}$, $i = 1, 2$. This homogenization theory is not limited to long-waves relative to the microstructure, one apparent failing is that the asymptotics appear to be only valid near the band edge frequencies but further refinements are possible, using foldings of the Brillouin zone, that extend the theory to provide complete coverage of the dispersion curves and provide accuracy at all frequencies, see e.g. [49] and further developments in [50, 51].
<figure><img src="image_12.png" /><figcaption>**Figure 8.** The isofrequency contours for the second band at 18 kHz. Left: zoomed-in 3D representation shows concentric near-squares centred around $\Gamma$ and flat perpendicular to the $\Gamma - X$ direction. Right: 2D representation of the isofrequency contour in the plane $\kappa_3 = 0$ with the free-space isofrequency contour at 18 kHz in yellow dashed line. The large cube shown on the left panel appears as smaller squares on the right panel.</figcaption></figure>

Figure 7 shows an ellipsoid flattened in the $\Gamma - X$ direction. The black lines are obtained from direct computation with Comsol Multiphysics. This isofrequency surface is superimposed by a red surface obtained from HFH. We point out that this demonstrates the sharpness of HFH approximation not only along the IBZ boundary as in figure 5, but also within the IBZ. The effective medium at 14.2 kHz is described by the matrix coefficients $T_{11} = T_{33} = 31.4$ and $T_{22} = 244.4$, which confirms the large anisotropy along the $x_2$ axis in agreement with the observed guiding effect in that direction with the sonic crystal in figure 4. Interestingly, the 2D acoustic magnifying hyperlens of Zhang's group consisting of brass fins in air embedded on a brass substrate [52] depicts elliptical isofrequency contours similar to those in figure 7 at a working frequency of 6.6 kHz. When viewed from the $\Gamma$ point, elliptical contours within the dashed circle in figure 7 appear locally hyperbolic. One indeed observes a collimated beam on the exit of the hyperlens in [52], analogous to that on the exit of the sonic crystal in figure 4. In electromagnetism [22] demonstrate self-collimation for a source inside a finite crystal in the microwave regime where the isofrequency contours are also hyperbolic and here we see a similar effect, but now for acoustic waves.

## 4.2. Interpretation of collimated beam forming at 18 kHz with Fermi surfaces and slowness curves

We directly use the isofrequency curves for this case as they give a self-explanatory insight of the physics. Indeed, by looking at the isofrequency surfaces in figure 8 we see that they form a slightly deformed cube. The concentric near-squares have almost flat edges in the $\Gamma - X$ direction and notably the isofrequency contour in the plane $\kappa_3 = 0$ surrounds the $\Gamma$ point in figure 8 leading to anti colinear phase and group velocity (see also movie at [Multimedia 2]). In the crystal this is reminiscent of an endoscope or self-guiding effect where the narrow support of the source is conserved across the crystal and then leaves the crystal as a focused beam.

Unlike the 14.2 kHz observation of collimation, at 18 kHz we experimentally observe the excitation of several Bloch modes in the crystal as expected from figure 3(a). The field is dominated by the mode giving these cubic equifrequency surfaces that lead to the self-guiding effect in the crystal. We note the similarity of the isofrequency contour at 18 kHz with that in [16] that leads to a self-guiding effect in 2D photonics whereas we propose a 3D experimental realization of the effect in acoustics. We further refer to [21] for validation in the microwave regime of self-collimation and to [53] for a 3D focusing effect via all-angle-negative refraction.

## 4.3. Quantitative measurement of the collimation effect

The acoustic source is a 48 mm diameter emitter corresponding to two times the acoustic wavelength at 14 kHz and 2.5 times at 18 kHz. It is located at a 20 mm distance from the crystal input side. This distance is 0.8 times the wavelength at 14 kHz and one wavelength at 18 kHz. Thus, the source is considered to be within the near field region for each of the frequency studied in the article. The source diameter is larger than the wavelength, leading to a non-isotropic radiation pattern from the source. We present additional data processing to quantify the collimation effect of the crystal. Figure 9 shows the radiation pattern of the source
<figure><img src="image_13.png" /></figure>

both in free space and in presence of the phononic crystal. At 14.2 kHz the source shows a beam width of 83° at half maximum, this value is reduced to 32° in presence of the crystal. Similar data processing at 18 kHz shows a reduction of beam width from 59° in free space to 28° in presence of the crystal. This confirms that the crystal achieves a significant collimation effect, by reducing the beam width by a factor 2.6 and 2 at 14 kHz and 18 kHz respectively.

## 5. Concluding remarks

We have conclusively shown strongly collimated beam effects through a three-dimensional acoustic phononic crystal consisting of a finite rectangular cuboid array of 4 × 10 × 10 polymer spheres 1.38 cm in diameter in air, arranged in a primitive cubic cell with pitch 1.5 cm for chirps with central frequencies at 14.2 kHz and 18 kHz. Our study is making use of three-dimensional higher-order multi-scale finite element simulations that enable us to capture the essence of the wave phenomena for our scattering problem in the time domain. These simulations show excellent agreement with time-domain experiments based on an in-house laser feedback interferometer setup. The prediction of the two frequencies where these effects are achieved requires analysis of band spectra of a three-dimensional periodic structure, that reveals the prominent role played by the Floquet-Bloch waves in the infinite periodic counterpart of the finite rectangular cuboid array of rigid spheres. Visualization of four-dimensional Fermi surfaces in the infinite crystal is somewhat delicate and thus we further analyse effective properties of the periodic structure through the lens of a high-frequency homogenization method. Our analysis shows one needs to be extra-careful when just considering edges of the three-dimensional Brillouin zone instead of its volume. While Fermi surfaces (here 4D dispersion surfaces) contain within them all the information required to completely describe the stop band structure of a 3D periodic structure, dispersion curves such as in figure 6 can obscure important features such as the remarkable isofrequency contours in figures 7–8, and most notably at 18 kHz wherein the phase and group velocity are anti-linear. This a frequency regime where the dispersion curves in figure 6 are difficult to decipher. Moreover, it is not obvious to conclude from figure 6 that there should be highly directive wave propagation within the crystal at 14.2 kHz, as evidenced by figure 4. On the other hand, the combination of HFH that reveals strong effective anisotropy of the crystal through inspection of the $T_{ii}$ coefficients deduced from equation (9), and of the isofrequency contours that display a remarkably flat isosurface in figure 7, sheds light on the waveguiding effect through the phononic crystal at 14.2 kHz. The sonic crystal can be viewed as a waveguide counterpart to a uni-directive metamaterial antenna with a vanishing effective refractive index such as in [54]. It is well known that curved wave fronts emanating from a source in free space are transformed into planar ones in a zero index slab [55, 56]. In the present case, the combination of low group velocity in the $X_1$ direction and large dynamic anisotropy of the sonic crystal at 14.2 kHz makes it possible to design an effective highly anisotropic low-index waveguide. Quite unexpectedly for such a simple phononic crystal consisting of a finite rectangular cuboid array of rigid spheres, it turns out that isofrequency contours can design for practically implementable collimated beams at two seemingly non remarkable frequencies, at least according to the dispersion diagram. It is only through further inspection of isofrequency surfaces and contours, in conjunction with high frequency homogenization, that the true nature of the wave propagation within the crystal is unveiled at these frequencies. Our study thus emphasizes the richness of physical phenomena via anomalous dispersion in 3D phononic crystals. We hope our study will foster experimental efforts in manipulation of sound in miscellaneous phononic crystals whose dynamic effective properties of practical importance for wavefront shaping might have been overlooked.
Data availability statement

The data cannot be made publicly available upon publication because no suitable repository exists for
hosting data in this field of study. The data that support the findings of this study are available upon
reasonable request from the authors.

Acknowledgments

K B and A D R acknowledge financial support provided by the Australian Research Council (DP210103342).
G L acknowledges support from Newton International Fellowships Alumni following-on funding awarded by
The Royal Society and Early Career Scheme (Project Number: 27301921), RGC, Hong Kong. M D, R A and S
E acknowledge financial support from the Excellence Initiative of Aix-Marseille University—A*MIDEX, a
French ‘Investissements d’Avenir’ program under the Multiwave chair of Medical Imaging and Chair
Professor Junior.

ORCID iDs

A L Vanel https://orcid.org/0000-0002-3603-1330

M Dubois https://orcid.org/0000-0002-9195-7363

C Tronche https://orcid.org/0000-0002-1387-4781

Y-T Wang https://orcid.org/0000-0003-3283-4424

A D Rakić https://orcid.org/0000-0002-4615-2240

K Bertling https://orcid.org/0000-0003-1087-2315

S Enoch https://orcid.org/0000-0003-0335-726X

R V Craster https://orcid.org/0000-0001-9799-9639

G Li https://orcid.org/0000-0003-0615-2681

S Guenneau https://orcid.org/0000-0002-5924-622X

J Perchoux https://orcid.org/0000-0002-5097-452X

References

[1] Vanel A L, Schnitzer O and Craster R V 2017 *Europhys. Lett.* **119** 64002

[2] Gao K, Fu S and Chung E T 2018 *Bull. Seismol. Soc. Am.* **108** 966–82

[3] Bertling K, Perchoux J, Taimre T, Malkin R, Robert D, Rakić A D and Bosch T 2014 *Opt. Express* **22** 30346–56

[4] Bertling K, Veidt M, Perchoux J and Rakić A D 2023 *Opt. Express* **31** 32761–71

[5] Craster R V and Guenneau S 2012 *Acoustic Metamaterials* (Springer)

[6] Kadic M, Bückmann T, Schittny R and Wegener M 2013 *Rep. Prog. Phys.* **76** 126501

[7] Cummer S A, Christensen J and Alù A 2016 *Nat. Rev. Mater.* **1** 1–13

[8] Dupont G, Movchan A, Enoch S and Guenneau S 2019 *Front. Mater.* **6** 50

[9] Liu L B and Albert D G 2006 *J. Acoust. Soc. Am.* **119** 2073–83

[10] Yves S and Alù A 2021 *J. Acoust. Soc. Am.* **150** 2040–5

[11] Bückmann T, Stenger N, Kadic M, Kaschke J, Frölich A, Kennerknecht T, Eberl C, Thiel M and Wegener M 2012 *Adv. Mater.* **24** 2710–4

[12] Liu Y, Xu W, Chen M, Yang T, Wang K, Huang X, Jiang H and Wang Y 2020 *Mater. Des.* **188** 108470

[13] Yang S, Page J H, Liu Z, Cowan M L, Chan C T and Sheng P 2004 *Phys. Rev. Lett.* **93** 024301

[14] Dubois M *et al* 2019 *Phys. Rev. B* **99** 100301

[15] Yu X and Fan S 2003 *Appl. Phys. Lett.* **83** 3251–3

[16] Chigrin D N, Enoch S, Torres C M S and Tayeb G 2003 *Opt. Express* **11** 1203–11

[17] Prather D W, Shi S, Murakowski J, Schneider G J, Sharkawy A, Chen C, Miao B and Martin R 2007 *J. Phys. D: Appl. Phys.* **40** 2635–51

[18] Bucay J, Roussel E, Vasseur J O, Deymier P A, Hladky-Hennion A C, Pennec Y, Muralidharan K, Djafari-Rouhani B and Dubus B 2009 *Phys. Rev. B* **79** 214305

[19] He Z, Heng Y, Peng S, Ding Y, Ke M and Liu Z 2009 *J. Appl. Phys.* **105** 116105

[20] Vanel A, Craster R, Colquitt D and Makwana M 2016 *Wave Motion* **67** 15–31

[21] Prather D W, Shi S, Venkataraman S, Lu Z, Murakowski J A and Schneider G J 2005 Self-collimation in 3D photonic crystals *Proc. SPIE* **5733** 84–93

[22] Lu Z, Shi S, Murakowski J A, Schneider G J, Schuetz C A and Prather D W 2006 *Phys. Rev. Lett.* **96** 173902

[23] Belov P A and Simovski C R 2005 *Phys. Rev. E* **72** 026615

[24] Gralak B, Enoch S and Tayeb G 2000 *J. Opt. Soc. Am.* **A** **17** 1012–20

[25] Luo C, Johnson S G, Joannopoulos J D and Pendry J B 2002 *Phys. Rev. B* **65** 201104

[26] Ceresoli L, Abdeddaim R, Antonakakis T, Maling B, Chmiaa M, Sabouroux P, Tayeb G, Enoch S, Craster R V and Guenneau S 2015 *Phys. Rev. B* **92** 174307

[27] Taimre T, Nikolić M, Bertling K, Lim Y L, Bosch T and Rakić A D 2015 *Adv. Opt. Photon.* **7** 570–631

[28] Fu S, Gao K, Gibson R L and Chung E T 2019 *Computat. Geosci.* **23** 997–1010

[29] Brillouin L 1953 *Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices* 2nd edn (Dover)

[30] Kittel C 1996 *Introduction to Solid State Physics* 7th edn (Wiley)

[31] Donati S 2011 *Laser Photon. Rev.* **2011** 1–25
[32] Ciddor P E 1996 *Appl. Opt.* **35** 1566–73

[33] Berenger J P 1994 *J. Comput. Phys.* **114** 185–200

[34] Wilson E L 1974 *Int. J. Numer. Methods Eng.* **8** 198–203

[35] Graham I, Spence E and Zou J 2020 *SIAM J. Numer. Anal.* **58** 2515–43

[36] Harrison J M, Kuchment P, Sobolev A and Winn B 2007 *J. Phys. A: Math. Theor.* **40** 7597–618

[37] Figotin A and Vitebskiy I 2006 *Waves Random Complex Media* **16** 293–392

[38] Boyles C 1969 *J. Acoust. Soc. Am.* **45** 356–64

[39] Climente A, Torrent D and Sánchez-Dehesa J 2010 *Appl. Phys. Lett.* **97** 104103

[40] Di Falco A, Kehr S C and Leonhardt U 2011 *Opt. Express* **19** 5156–62

[41] Chang T, Dupont G, Enoch S and Guenneau S 2012 *New J. Phys.* **14** 035011

[42] Romero-García V, Cebrecos A, Picó R, Sánchez-Morcillo V J, Garcia-Raffi L and Sánchez-Pérez J V 2013 *Appl. Phys. Lett.* **103** 264106

[43] Zhao Y Y, Zhang Y L, Zheng M L, Dong X Z, Duan X M and Zhao Z S 2016 *J. Opt.* **19** 015605

[44] Dong H Y, Cheng Q, Song G Y, Tang W X, Wang J and Cui T J 2017 *Appl. Phys. Express* **10** 087202

[45] Xie Y, Fu Y, Jia Z, Li J, Shen C, Xu Y, Chen H and Cummer S A 2018 *Sci. Rep.* **8** 16188

[46] Zhao L, Bi C, Huang H, Liu Q and Tian Z 2023 *Mech. Syst. Signal Process.* **199** 110468

[47] Zhao L, Lai C and Yu M 2020 *Mech. Syst. Signal Process.* **144** 106868

[48] Craster R V, Kaplunov J and Pichugin A V 2010 *Proc. R. Soc. A* **466** 2341–62

[49] Antonakakis T, Craster R and Guenneau S 2013 *New J. Phys.* **15** 103014

[50] Guzina B B, Meng S and Oudghiri-Idrissi O 2019 *Proc. R. Soc. A* **475** 20180547

[51] Assier R C, Touboul M, Lombard B and Bellis C 2020 *Proc. R. Soc. A* **476** 20200402

[52] Li J, Fok L, Yin X, Bartal G and Zhang X 2009 *Nat. Mater.* **8** 931–4

[53] Lu Z, Shi S, Schuetz C A, Murakowski J A and Prather D W 2005 *Opt. Express* **13** 5592–9

[54] Enoch S, Tayeb G, Sabouroux P, Guérin N and Vincent P 2002 *Phys. Rev. Lett.* **89** 213902

[55] Ziolkowski R W 2004 *Phys. Rev. E* **70** 046608

[56] Liberal I and Engheta N 2017 *Nat. Photon.* **11** 149–58