Nonlinear self-collimated sound beams in sonic crystals

We report the propagation of high-intensity sound beams in a sonic crystal, under self-collimation or reduced-divergence conditions. The medium is a fluid with elastic quadratic nonlinearity, where the dominating nonlinear effect is harmonic generation. The conditions for the efficient generation of narrow, non-diverging beam of second harmonic are discussed. Numerical simulations are in agreement with the analytical predictions made, based on the linear dispersion characteristics in modulated media and the nonlinear interaction in a quadratic medium under phase matching conditions.


I. INTRODUCTION
The beams of different kind of waves, such as electromagnetic or sonic, experience diffraction and broaden as they propagate in a homogeneous medium. This spreading of energy in space diminishes the wave amplitude of the beam in the axis along propagation, unless the spreading is balanced by some focusing mechanism. Still, it is possible to prevent this fundamental wave propagation property and create non-diverging beams. Among the beam patterns without divergence the most popular are the linear Bessel beams [1], and the nonlinear solitonic or self-trapped beams [2]. Recently, another method for creating linear non-diverging beam was proposed, for waves propagating in a periodic medium. Such beams have been named self-collimated beams, and were first proposed for light beams in photonic crystals [3], and later extended to other type of waves. The phenomenon of self-collimation has attracted much attention, as a technique to propagate optical, acoustical or even matter wave beams at long distances without a sensible loss of amplitude. Self-collimation is highly sensitive to the frequency: the divergence of beams can be reduced or even suppressed only at particular frequencies, those presenting particular dispersion characteristics, namely flat regions in the isofrequency contours [3]. The size of the self-collimated beam is also limited by the extension of such flat region in angular space. Self-collimation of low amplitude (linear), monochromatic acoustic waves has been demonstrated in 2D [4] and 3D [5] sonic crystals. More recently, the simultaneous self-collimation of two beams of different frequencies was also demonstrated experimentally [6]. These results show that the conditions for self-collimation can be achieved also for nonmonochromatic beams; in particular the case of the superposition of beams of one frequency and its second harmonic was considered in [6]. The latter results are valid in linear regime; actually in [6] both frequency components were present in the input beam, and the corresponding beams propagated in the crystal without nonlinear interaction between them.
In the linear case, the propagation of light and sound beams obey similar equations, and similar propagation characteristics are expected. The similarities between photonic and sonic crystals are well established [7], and have motivated many studies, where analogous effects in both systems have been investigated. The analogy, however, breaks for high amplitude waves, where nonlinear effects appear. For example, second and higher harmonic generation processes may be essentially different in optics and acoustics. One reason is the absence of intrinsic dispersion for acoustic waves propagating in homogeneous media. Nonlinear acoustical waves in nondispersive media as homogeneous fluids, eventually generate shock waves, which are not observed in optics. Also, the type and strength of nonlinearity may be different. While most common optical nonlinearities are cubic (kerr-type), in fluids and homogeneous solids, quadratic nonlinearity is dominant in acoustics. Even, nontraditional acoustic nonlinearities (power-law, hysteretic,...) are typical of some complex or microstructured acoustic media. In this sense, nonlinear effects of acoustic waves in periodic media, and in particular the self-collimation problem considered here are not a direct extension of the same effects in the optical case. Furthermore, the propagation of nonlinear acoustic beams in sonic crystals has never been addressed before.
The basic effect in nonlinear acoustics is harmonic generation [8]. It is known that efficient harmonic generation is only possible under fulfilment of phase matching conditions. For acoustic waves in fluids, this condition is rather natural, being always fulfilled for all harmonics due to the absence of dispersion, however in optics this requires special materials and special phase matching techniques [9].
Acoustic harmonic generation has been studied in a variety of highly-dispersive nonlinear media, as bubbly liquids, or acoustic waveguides [8], and weakly dispersive media as elastic plates [10,11], nonlinear porouselastic media [12] and in granular media [13]. It has been proven as a useful effect in different applications, as material characterization [14,15], ultrasound imaging and echography, [16], biological tissue characterization [17] and other medical ultrasound applications.
The purpose of this paper is to study nonlinear propagation of high intensity sound beams in periodic media, and in particular to demonstrate the formation of nonlinear self-collimated acoustic beams, and discuss the conditions under which this process occurs with maximal efficiency. A sonic crystal is designed, by using an iterative method, to fulfil the three conditions for optimal energy transfer between harmonics: flatness of isofrequency contour for each harmonic, phase matching and large overlap between distributions of the interacting Bloch modes. The predictions are checked by FDTD simulations of the nonlinear problem, that demonstrate the efficient generation of fundamental and second harmonic acoustic narrow beams.

II. NONLINEAR SOUND BEAM PROPAGATION MODEL
Several models can be used to describe nonlinear sound wave propagation though a fluid medium, with different levels of accuracy. An accurate description, when thermal and viscous effects are negligible, follows from the conservation laws of mass and momentum, can be written, respectively, in a Eulerian form [8]: where v is the particle velocity vector, p is the acoustic pressure, ρ is the total density field that can be expressed as ρ = ρ 0 +ρ ′ , where ρ 0 the ambient fluid density and ρ ′ is the acoustic density. The system is closed by the equation of state of the fluid, that under our assumptions is a pressure-density relation, p = p(ρ). A commonly used Top: Band structure (left) of a square lattice of rigid cylinders with r = 0.11a, where a is the lattice constant, immersed in water (right). Red and blue lines mark the bands (2nd and 8th) for which simultaneous self-collimation for both fundamental and second harmonic searched, along Γ − X direction. Bottom: Isofrequency contours for the 2nd (left) and 8th bands (right).Points denote the wavevectors for both waves, lying on flat segments respectively expression is obtained after Taylor expansion, keeping nonlinear terms up to second order. Then where B/A is the nonlinear parameter of the medium (which is known for most of materials, see e.g. [18]) and c 0 the sound speed in the medium. Note that Eqs. (1) and (2) also contain nonlinearities related to (1) mass and momentum advection , or (2) geometrical nonlinearities. However, for the second harmonic generation they are of minor importance compared with the nonlinear terms in the equation of state, Eq. (3).

III. SELF-COLLIMATION OF INTENSE ACOUSTIC BEAMS
We consider a narrow, intense acoustic beam incident on a 2D sonic crystal made of cylindrical scatterers of radius r embedded in a fluid, arranged in a square-lattice with lattice constant a. The corresponding filling fractio is f = π(r/a) 2 . The beam width is roughly 6 lattice periods. For the sake of simplicity the scatterers are considered perfectly rigid (the sound field is totally reflected from the wall of scatterer). Assuming water as a host fluid, the material parameters are ρ 0 = 1000 kg/m 3 , c 0 = 1490 m/s. The special conditions required for a sound beam to propagate without diffraction are presented in this section. The problem of self-collimation has been already discussed for linear, monofrequency [4,19] and bifrequency [6] beams. Since a nonlinear beam is composed by a fundamental frequency component and its high frequency harmonics, self-collimation of the nonlinear beam requires self-collimation of its constituent frequency components. We remind that in self-collimation regime the sonic beam does not spread diffractively because Bloch wave vectors lying on the flat segment of the spatial dispersion curve have equal longitudinal components and thus do not dephase mutually in propagation. In general, this flatness of the dispersion curve appears at a particular frequency, but as shown in [6] it can be also obtained for a wave and its second harmonic regarded they propagate in different propagation bands.
To illustrate this case, we show in Fig. 1 the dispersion diagram of the sonic crystal for small-amplitude excitations, obtained using the Plane Wave Expansion (PWE) method on a linearized version of Eqs. (1)-(3). The conventional form of the band diagram is represented on the trajectory along the principal directions of the crystal, Γ − X − M − Γ, which are the boundary of the irreducible Brillouin zone (BZ). Figure 1 (top) shows the dispersion diagram for Γ − X direction. The red line in Fig. 1 denotes the 2nd propagation band. The fundamental (driving) field lies on this band. Its frequency Ω is chosen such that the corresponding isofrequency contour contains flat regions. The blue line in Fig. 1 denotes the 8th propagation band. The second harmonic frequency, 2Ω, lies in this band for the particular crystal parameters considered. As shown in [6], for a given crystal, it is possible to choose the fundamental frequency such that the isofrequency contours for both frequencies present flat regions. In our particular crystal, this happens when Ω = 0.125 [these are dimensionless frequencies, related to physical frequencies ω as ω = Ω(2πc 0 /a)]. Similarly, a normalized Bloch wavevector is defined as K = k x (a/π) Such condition is necessary to achieve self-collimated propagation of the nonlinear beam. In order to obtain an efficient generation of the second harmonic, together with simultaneous self-collimation for both waves two additional geometric conditions must be fulfilled, related to the wavenumber and spatial shape of the interacting shows that for the self-collimated second harmonic, two solutions (modes A and B, with distinct K) are found, phase and non-phase matched, respectively beams. These conditions have been discussed for photonic crystals [20,21].

A. Phase matching conditions
It is well known from nonlinear optics [9] that a proper phase relationship between the fundamental and second harmonic waves must be satisfied for an efficient nonlinear frequency conversion along the propagation direction. In a dispersive medium, the wavenumbers of first harmonics do not combine to result precisely in wavevector of second harmonics, and a phase mismatch ∆K = 2K(Ω) − K(2Ω) occurs. As a consequence, the second harmonic field is limited in amplitude: it does not grow linearly but oscillates in propagation, with a characteristic period given by the coherence length l c = π ∆k = a ∆K [9]. The conversion efficiency into second harmonics generally is smaller in optics than in acoustics, because of the inherent material dispersion for light waves (absent for sound waves in fluids), that causes the fundamental and second harmonic waves to travel along the crystal with different phase velocities. Thus, the presence of the scatterers is the only important source of dispersion in the acoustic case.
Phase matching corresponds to ∆K = 0. Figure 2 shows that it can be actually achieved for the pair of frequencies where self-collimation occurs, as follows from the previous analysis. There, we represent the dispersion branches involved in self-collimation, along Γ − X direction, as in Fig. 1. Fundamental and second harmonic modes correspond to the crossings of the dotted horizontal lines with the corresponding dispersion branches. For a given fundamental frequency Ω, in order to check the fulfillment of the phase matching condition, the curve corresponding to the "double" of the 2nd band dispersion curve, 2Ω(2K), has been represented in Fig. 2 as a red dashed line. Phase matching is satisfied at the intersection between this curve with the corresponding curve at the 8th band. This corresponds to the mode labeled A in Fig. 2, which is phase-matched with the fundamental mode.
Note that due to the concavity of the 8th band, at frequency 2Ω a second mode labelled with B in Figure. 2 can be also excited. This solution presents a large phase mismatch with the fundamental mode, and its contribution to the second harmonic field is negligible.
The simultaneous fulfillment of both conditions is obtained by an iterative procedure, which implies a redesign of the crystal parameters. The procedure is as follows: we start from a pair of frequencies (Ω, 2Ω) showing self-collimation for a given crystal parameters. Around this doublet, we seek the closest pair of frequencies Ω ′ = Ω + δΩ and 2Ω ′ showing phase-matching. Then the isofrequency curves (Fig. 1) are again calculated in order to evaluate the deviation of flatness in the isofrequency contours. The sonic crystal parameters are then modified, e.g. by a slight variation of the filling fraction, in order to tune the dispersion relations to get again selfcollimation conditions. The process is repeated again until both conditions (flatness and phase matching) are simultaneously satisfied. We note that, despite this is out of the scope of this paper, optimization techniques as genetic algorithms can be applied here in order to find an optimal structure.

B. Nonlinear coupling of Bloch modes
Efficient energy transfer between harmonics requires also a strong mode coupling, which depends on the spatial overlapping between the two interacting waves. For plane waves in a homogeneous medium a perfect spatial nonlinear coupling between first and second harmonic is assured, since mode overlapping is maximal. The propagation eigenmodes in a periodic medium are Bloch waves, whose amplitudes are spatially modulated and does not necessarily overlap. If two modes do not overlap in space, the energy transfer is less efficient even if they are phase matched. The amount of energy transfer can be estimated evaluating the spatial overlap between the envelopes of the corresponding Bloch modes. Let B 1 and B 2 be the spatial envelopes of the Bloch modes of the fundamental and and second harmonic waves, respectively. The nonlinear coupling coefficient is calculated as the cross-correlation between the functions B 2 1 and B 2 normalized in such a way that unity would correspond to the perfect matching of the modes [20]. We define the coupling coefficient as where the upper integral is calculated in the nonlinear medium from one unit cell, while the lower integrals are taken over the entire unit cell. To calculate B 1 and B 2 we solve the eigenvalue problem for the pressure field by means of the PWE method, which converts the differential equation to an infinite matrix eigenvalue problem that can be truncated and solved numerically. For that, we follow the same procedure as in [19] however the problem is solved inversely, i.e. for a given frequency, the corresponding wave vectors, satisfying the phase matching condition, are obtained. Then B 1 and B 2 are obtained as the eigenvectors corresponding to fundamental and second harmonic frequencies, respectively. In Fig. 3 we plot the spatial distributions of B 1 and B 2 , respectively, for the selected final design where a coupling coefficient of κ = 0.85 is obtained. This value is of the same order as the coupling in homogeneous media, κ = 1, and therefore sufficiently large for an efficient harmonic generation.

C. Numerical simulation
A full-wave nonlinear simulation was performed, using the FDTD method, to validate the efficiency of the second harmonic generation in the proposed structure. The crystal parameters, obtained after the iterative procedure describe above, are as in Figs. (1) and (2). The source is a plane piston with a width of 6a, located near the crystal, and radiating a harmonic wave with normalized frequency Ω = 0.125. In order to minimize numerical dispersion a computational grid with N λ = 45 elements per wavelength was used, and a Courant-Friedrich-Levy number of S = 0.95. In Fig. 4 we present the numerically obtained spatial distributions of the fundamental and its nonlinearly generated second harmonic. As predicted, both beams are nearly collimated. For comparison, the beam spatial distribution calculated for an homogeneous material (removing the crystal) are represented in Fig. 4 (e-f), where the diffractive broadening of the beams is visible. Transversal intensity distributions are shown in Fig. 4 (a-b) for a distance 80a, where beam widths are compared with the reference beams in the homogeneous medium, broadened by diffraction. The pressure amplitudes along the beam axis are shown in Fig. 5 for each harmonic. Here, the analytic so- lution for a nonlinear plane wave propagating in a homogeneous (nondispersive) medium is plotted for reference (dashed lines) [18], given by p n /p 0 = 2J n (nσ)/nσ, where J n is the Bessel funcion of order n, σ = (εωp 0 /ρ 0 c 3 0 )x is coordinate normalized to the shock formation distance, p 0 is the pressure at the source and β = 1 + B/2A the nonlinearity parameter. Such analytical solution is valid in the pre-shock region σ < 1. The growth rate of the self-collimated second harmonic beam propagating in the crystal matches well the growth rate of a plane wave in a homogeneous medium in such preshock region, which is a consequence of the weak divergence of the beam and the high degree of phase matching. Also, the second harmonic field can reach even higher amplitudes than those corresponding to nondispersive media (where harmonics decay beyond the shock formation distance). The latter effect can be understood in terms of phase mismatch of higher harmonics: in homogeneous media all harmonics are phase matched while in the crystal only the second harmonic is phase matched. In Fig. (5) the third harmonic is also plotted, where its small contribution to the beam is evident. The phase mismatch in third and higher harmonics decrease the energy flow into these components. Finally, the amplitude in Fig. (5) decays because non-perfect conditions for self-collimations, that makes the beam to start diverging after a long distance, or nonperfect phase matching, which results in a beating period with long coherence length.

IV. CONCLUSIONS AND REMARKS
We have demonstrated the possibility of efficient second harmonic generation of sound in a sonic crystal, by means of the formation of narrow, weakly diverging nonlinear acoustic beams. Three conditions must be simultaneously present for a efficient second harmonic generation, which are: 1) simultaneous self-collimation, 2) phase mathing and 3) high spatial coupling of interact-ing harmonics. The use of simultaneous self-collimation regime limits the diffraction of both harmonic beams, maintaining the amplitude at the axis and therefore the nonlinear interaction. Under ideal conditions (no divergence and losses), the decrease of the first harmonic beam is mainly attributed to the energy transfer to second and higher harmonics. The sonic crystal parameters can be chosen to fulfill phase matching with the second harmonic, maximizing second harmonic generation due to synchronous cumulative interaction. Finally, the spatial coupling (overlapping) between interacting modes is also analyzed by calculating a nonlinear coupling coefficient. It is shown that its value (κ = 0.85 for the case studied) is not far from the ideal case, revealing a strong spatial overlap between both Bloch modes that leads to high energy transfer.
The study show that linear dispersion characteristics (band structures, isofrequency contours) can be used to predict the behaviour of nonlinear beams propagating in periodic media. This opens the possibility of extending the study of nonlinear sound beam propagation in sonic crystals to other cases of interest. For example, crystals with higher filling factors present full frequency bandgaps, that may be used to filter out the propagation of selected higher harmonics. In this sense, sonic crystals can be a way to control the spectrum of intense acoustic waves, using the strong dispersion properties introduced by the periodicity.