Point source in a phononic grating: stop bands give rise to phonon-focusing caustics

We use locally-excited gigahertz surface phonon wavepackets in microscopic line structures of different pitches to reveal profound anisotropy in the radiation pattern of a point source in a grating. Time-domain data obtained by an ultrafast optical imaging technique and by numerical simulations are Fourier transformed to obtain frequency-filtered real-space acoustic field patterns and k-space phononic band structure. The numerically-obtained k-space images are processed to reveal an intriguing double-horn structure in the lowest-order group-velocity surface, which explains the observed non-propagation sectors bounded by caustics, noted at frequencies above the bottom of the first stop band. We account for these phonon-focusing effects, analogous to collimation effects previously observed in two- and three-dimensional lattices, with a simple analytical model of the band structure based on a plane wave expansion. As the frequency is increased, a transition to dominant waveguiding effects along the lines is also documented.


Introduction
Collimation-or non-diffractive propagation-of waves inside periodic structures has been investigated for two-dimensional (2D) and three-dimensional (3D) lattices for both electromagnetic [1][2][3] and acoustic [4][5][6] waves because of the interesting possibilities for waveguiding. Analogous effects in 2D lattices have been observed for water waves diffracting over a bottom of periodically-drilled holes, [7] and are predicted to occur for matter waves in Bose-Einstein condensates [8]. Collimation occurs when points on constant-frequency contours in k-space (wave-vector space) exhibit zero curvature. Collimation has not been investigated for in-plane propagation in one-dimensional (1D) or grating structures, the only studies having been done on beam steering in optical gratings [9,10].
In contrast to optics, in acoustics it is relatively easy to produce point (i.e. sub-wavelengthsized) sources. This led to a series of investigations into the strong anisotropy in the radiation pattern from point sources in natural crystals, [11] a phenomenon known as phonon focusing that is closely analogous to the above-mentioned collimation effects. The most striking feature is the existence of caustics, i.e. directions in which the angular energy density approaches infinity at large distances from the source. These caustics correspond to the same zero-curvature points of the constant-frequency contours that yield the collimation effects. Considering the occurrence of phonon focusing in natural crystals, the existence of this phenomenon in phononic crystals (PCs) may not seem particularly surprising. In the long-wavelength limit, the acoustic properties of a PC are described by an effective elastic constant tensor governed by the symmetry of the structure, and long-wave propagation effects should therefore be similar to those encountered in natural crystals; this has indeed been observed in synthetic anisotropic materials such as fiber-reinforced composites [12]. However, at wavelengths comparable to the structure period, acoustic properties of periodic structures become profoundly different from those of natural materials-the most pronounced effect being the formation of stop bands. Narita et al [13] pointed out that anisotropic propagation from a point source is made more conspicuous by stop band formation; however, the nature of the connection between phonon focusing and the acoustic band structure has not been elucidated. Does the stop band formation result in specific anisotropic effects that are not observed in phonon focusing in natural anisotropic materials? Here, we address this question by putting a point source inside a 1D periodic structure fabricated on a solid surface and studying its radiation pattern by both experiment and simulation. We are concerned here with point-source radiation anisotropy inside the structure, leaving outside the scope of this paper any phenomena involving refraction at interfaces [14]. A second goal of this paper is to develop a calculation technique capable of both analyzing band structure and simulating time-domain data. Gratings provide a generic test case for investigating scattering in periodic structures. Although much work has been done on acoustic in-plane propagation in gratings, [15][16][17][18][19][20][21][22][23][24][25][26][27] stop bands above the first two or three have not been investigated in detail, and neither has their dependence on the periodicity of the grating for a given grating thickness.
Using point-generated gigahertz surface-phonon wavepackets, we investigate in detail two different gratings composed of microscopic alternating strips of copper and silica on crystalline silicon substrates. By temporal Fourier transforms we obtain frequency-filtered images in real space and demonstrate phonon-focusing effects, revealing non-propagation sectors bounded by caustics (forming an X shape). By spatiotemporal Fourier transforms we analyze the acoustic band structure. We then derive the topology of the lowest-order sheet of the group-velocity surface, important for understanding the speed at which phonon pulses propagate in a given direction. The basic physics of these effects is accounted for by a simple model based on a plane wave expansion. At higher frequencies we document a transition to waveguiding effects for acoustic propagation along the grating lines. Many of these results, in particular the role of stop bands in the formation of caustics in the radiation pattern, will be applicable to electromagnetic and other types of waves.
In section 2.1 we describe the experimental setup and samples, and in section 2.2 we explain the simulations. Then in section 3 we present the results for the 4 µm period sample: in section 3.1 in real space, including frequency-filtered data, and in sections 3.2 and 3.3 in k-space. In section 4.1 we present the analytical model, and compare it with the results for the 4 µm grating structure. We then compare the shape of the group-velocity surface derived using the numerical results and the analytical model; in section 4.2 we use this model to interpret the phonon focusing. We then describe in section 5 results for a sample with 10 µm period and discuss the effect of the different spatial period. We conclude in section 6. The appendix describes the Fourier analysis of a Bloch wave in a periodic phononic structure.

Experimental setup and samples
Surface phonon imaging is done with an optical pump and probe technique [28]. In brief, light pulses of wavelength 400 nm, duration ∼200 fs and pulse energy 0.2 nJ are focused to the center of a Cu strip to a small spot (∼2 µm in diameter), non-destructively launching surface phonons with displacements in the 10 pm range. The detection is performed with two scanned 800 nm probe beams focused to a similar spot size using a common-path interferometer configuration.
The 76.3 MHz repetition rate of the laser corresponds to a period of 13.1 ns. We typically record 40 images, equally spaced in time within this period, which represent out-of-plane surface velocity maps (denoted in this paper by f (r, t), where r is the in-plane position vector). Frequency-filtered images obtained by temporal Fourier transforms are therefore obtained in steps of 76.3 MHz. The broadband surface-phonon pulses contain frequencies up to 1.3 GHz, showing a maximum amplitude between 200 and 700 MHz and a dominant acoustic wavelength ∼ 10 µm. The range of wavelengths significantly excited is 2-20 µm.
We investigate 1D PCs with periods a = 4 and 10 µm in the form of surface gratings. The former sample and raw data set is identical to that used by [19]. We present more results from this data set, in particular in real space, and provide an in-depth analysis. The latter sample, taken from the same wafer as the first sample and very close to it, was measured using the same apparatus. Alternating polycrystalline copper and amorphous silicon oxide lines of thickness 800 nm and width 2 or 5 µm are deposited perpendicular to the [011] direction on silicon (100) substrates (see figure 1(a) for a cross section for a = 4 µm). A layer of tantalum of thickness 25 nm serves as a diffusion barrier for the copper. On top of this specimen a 30 nm polycrystalline gold layer is sputtered to achieve a uniform optical reflectivity. The thin gold layer and the thin tantalum diffusion barrier do not significantly affect the surfacewave dispersion [29,30]. (We verified this by simulations with and without these components.) The imaged area is 150 µm × 150 µm for the a = 4 µm sample and 100 µm × 100 µm for the a = 10 µm sample. The PC region of the sample is sufficiently large for surface-wave reflections from the boundaries to be negligible.

Time-domain simulations
Simulation methods for PCs in the frequency domain using a single unit cell and periodic boundary conditions are common [23,31]. Although less computationally efficient, modeling the full structure in the time-domain is better suited for comparison with time-resolved measurements (see, e.g., [32] and [33]), as in the present paper. We use a commercial timedomain finite element modeling (FEM) package PZFlex (Weidlinger A, Inc.). This exploits 3D elements, each consisting of eight nodes arranged on an orthogonal grid, and uses an explicit time-integration method. We make use of the material and geometrical symmetry of the 1D PCs: in the simulations only one quarter of the sample needs to be modeled by the use of symmetrical boundary conditions. Details of the simulation geometry for a = 4 µm are shown in figures 1(b) and (c).
The thermoelastic laser excitation is represented by a simplified elastic-dipole model [34][35][36] that makes use of a radially-directed horizontal force with a spatial distribution given by r e −r 2 /d 2 applied over a circular region with a radius of 2 µm, where r is the radius and excitation used. The applied initial force is radially distributed over a 4 µmdiameter circle in the center of the sample. In the simulation we include the gold layer and the grating, but we neglect the tantalum layer (see text). The lateral dimension of the substrate L shown is taken as 127 and 122.5 µm for a = 4 and 10 µm, respectively. and d = 1 µm. The temporal variation of the excitation is a steplike function (a quarter period of a sinusoid) with a 1 ns rise time. (The precise choice of this function and rise time are ad hoc, but suffice to yield a frequency spectrum in the simulation similar to that observed in experiment.) To reduce the effect of any acoustic reflections from the boundaries, the sample is continued outside the PC area for a further 30 µm as a silicon substrate with copper and gold coatings. At the closing vertical and bottom planes absorbing boundaries are applied. Larger regions of the sample than imaged in experiment, and sufficiently long simulation times (39.3-78.6 ns) were chosen for the simulation (see figure 1(c) and caption) in order to improve the wave-number and frequency resolutions, respectively. This yielded corresponding resolutions of ∼2 and 3-6 times greater than in experiment. (The multiple reflection of surface acoustic waves within the sample ensured that the acoustic field did not significantly die away and thereby effectively reduce the frequency resolution associated with the longer simulation times.) We also checked the effect of multiple-pulse excitation [37] at the laser repetition period (13.1 ns) to see what effect this had on Fourier-transformed data; the results within a window of 13.1 ns (including the effect of five previous pulses) agreed very well with those making use of a single pulse, and so we have retained the latter type of simulation for comparison with experiment in this paper. The accuracy of the simulation, making use of ∼10 8 nodes, was checked with different discretizations. 7 Since only displacements were calculated in the simulations, we evaluated temporal Fourier transforms of the surface velocity by multiplying the temporal Fourier transforms of the displacement by the factor i ω, where ω is the angular frequency. (When implementing the spatial Fourier transform for the numerical and experimental data in this paper, the imaginary part of the temporal Fourier transform is used 8 .) The material properties were taken from the literature 9 . Polycrystalline Cu, amorphous SiO and polycrystalline Au are assumed to be isotropic; silicon was modeled as a cubic material. (The Ta layer was not included, for simplicity.) In the simulations intrinsic material ultrasonic attenuation is neglected; at the frequencies concerned this is a very good approximation for the distances of propagation (∼100 µm) in question and for the materials concerned [38][39][40].

Analysis in real space
Figure 2(a) shows a simulated image for out-of-plane surface displacement (denoted by u(r, t)) at approximately 19.6 ns after excitation of the a = 4 µm grating. The effect of the cubic Si substrate on the surface-phonon propagation is evident in the rounded-square shape of the wave fronts far from the source [41]. One can also detect the expected twofold symmetry nearer to the source caused by the scattering at the line pattern.
Vertical sections of the z component of the simulated displacement field in the x and y directions are shown in figures 2(b) and (c), with excitation on the left hand side. Waves with a higher phase velocity and longer wavelength than the grating period, located at the front of the wave packets, exhibit a larger penetration depth, whereas waves with a lower phase velocity and shorter wavelength than the period, located in the middle and latter parts of the wave packets, exhibit a smaller penetration depth, and are more influenced by grating scattering 7 The PC is reproduced on an orthogonal grid with elements of dimensions 0.095 × 0.095 × 0.10 µm 3 (δx × δy × δz) at the top; the gold layer is modelled with one element over its thickness resulting in the smallest element dimension δz = 0.03 µm. The vertical dimension of the elements is gradually increased to δz ≈ 0.38 µm at the bottom of the substrate. Numerical simulations of wave propagation are influenced by numerical dispersion [54], dependent on the element size and time step. This particularly affects the higher frequencies in broadband signals, where it leads in general to a lower wave velocity. Si and Cu have the lowest wave velocities among the materials considered in this paper, and so a Cu layer on a Si substrate suffices to check the simulation accuracy. Two test simulations were carried out with average top-element dimensions 0.13 and 0.19 µm. The dispersion curves in the [001] direction were evaluated up to 1.5 GHz and compared to a semi-analytical solution [47]. Differences in the wave velocity compared to those predicted for the present simulations (that make use of an average topelement dimension of 0.097 µm) remained below 1 and 5% for the finer and coarser meshes, respectively. The temporal discretization of the model (5-10 ps) is chosen for numerical stability, and is derived from the smallest element dimension and the maximum wave velocity [33,55]. Since this high temporal resolution is unnecessary for the acoustic frequencies (<1.5 GHz) investigated, we only save data for every 15th time step (corresponding to a ∼100 ps effective time resolution). The simulations were carried out on the University of Strathclyde High Performance Cluster. 8 Either the real or imaginary parts of the temporal Fourier transform can be used, or combinations of the two, but both give similar information in the present study. 9 For the deposited materials, density ρ and elastic constants c 11    in the x direction. The wave packet in figure 2(b) (corresponding to the y direction) contains waves that also have ∼ 4 µm or shorter. However, in the x-direction the wave packet only contains significant amplitude for waves with > 4 µm. In both the x and y directions the scattering of waves into the bulk can be observed.
The spatial distributions of the magnitude of the temporal Fourier transforms |F(r, ω)| of the simulated z-directed surface velocity, denoted for simplicity by A, along the x and y axes are shown in figures 2(d) and (e) (axes along which plane waves do not scatter to other directions). The Fourier magnitude spectrum for the y direction shows a maximum between 600 and 700 MHz, with frequencies up to ∼1.5 GHz determined by the duration of the simulated excitation (see figure 1(a)). The x-directed and y-directed waves have a similar magnitude below 500 MHz. Above this frequency the x-directed waves are highly attenuated in certain frequency ranges, indicating the presence of several stop bands: the first large stop band occurs between 500 and 600 MHz, a second at ∼900 MHz, and a third at ∼1.25 GHz. These arise from interactions between different wave modes (see figure 4 and the next subsection). Figure 2(f) shows sections of the acoustic spectra at a distance of 7 µm from the excitation for both the x and y directions.     (5) and (6): maps of the magnitude of the spatiotemporal Fourier transform |F(k, ω)| of the z-directed surface velocity, corresponding to constant-frequency contours, for the simulation and experiment, respectively. Row (7): results of an analytical model for the Rayleigh-wave branch based on a plane wave expansion. Row (8): interpretation with Bloch harmonics (BH) in the extended zone scheme. Vertical dashed lines: boundaries of the first BZ at ±π/a. The data in row (6) are represented on slightly different scales compared to [19]. The interpretation for 687 MHz here is revised compared to [19].
To visualize the spatial distributions of the temporal Fourier transform it is convenient to plot frequency-filtered images; denoting the magnitude and phase of the temporal Fourier transform of the z-directed surface velocity by A and ψ respectively, the real part of this Fourier transform, A cos ψ, represents a snapshot of the instantaneous wave field [5]. In other words A(ω, r)cos ψ(ω, r) is the acoustic field pattern at a particular frequency.
Images of the simulated spatial distributions in the x-y and x-z planes are presented in rows (1) and (2)  Immediately obvious in the x-y plane images for frequencies above 382 MHz are characteristic X-shaped phonon-focusing patterns, seen more clearly in the simulations. As the frequency increases, the angle between the two straight-line components making up the X shape decreases. The detailed origin of this effect is described in section 4, and is shown for an analytical model in row (3) of figure 3. Similar effects should be observable in the analogous system of oblique propagation in a superlattice, [13,42] although this has not been experimentally confirmed.
The frequencies chosen for figure 3 lie around the range of the first large stop band at 500-600 MHz (see figure 2(f)). At 382 MHz (figure 3(a)), the orthotropic nature of the grating leads to elliptical wave fronts with the major axis along the y direction.
The better frequency resolution of the simulation (∼25 MHz) allows extra x-y plane images (at 483 and 508 MHz) to be obtained. These images and the one at 611 MHz are influenced by the proximity or coincidence of their frequencies with the first large stop band, and show a reduced amplitude for propagation in the ±x directions. (See both the x-y and x-z images in Rows 1, 2 and 4.) At 611 MHz, radiation to the bulk is evident in the x-z plane image in Row 2.
The numerical and experimental x-y plane images show approximate agreement over the whole frequency range shown, but at the higher frequencies some strong scattering is apparent in experiment. At 687 MHz both the simulated and experimental images hint at an independent mode along the ±y directions with a different phase compared to the remaining field. We interpret this to be caused by a guided mode in the central Cu line, which behaves as a waveguide for wavelengths commensurate with the dimensions of the line [43]. This phenomenon will be demonstrated more clearly for the a = 10 µm grating in section 5. Because of the lower wave velocities in the Cu line compared to the Si substrate or SiO lines, critical angles for total internal reflection exist at every boundary of the Cu line, and hence waveguiding can occur, as previously demonstrated in a similar structure on millimeter length scales [17] and also in single Cu lines [44]. This waveguiding is not expected to be perfect owing to the evanescent coupling between adjacent Cu lines.

Dispersion relations for the k x and k y directions
The complete dispersion relation of the surface phonons is a set of surfaces in 3D (k, ω) space, where k = (k x , k y ) is the 2D wave vector and ω is the angular frequency. In the appendix we show in detail how a temporal series of images for the wave field f (r, t) allows direct access to the complete dispersion relation using a combination of 2D spatial and 1D temporal Fourier transforms. The analysis is general, applying to wave fields other than acoustic fields, and is used to interpret our results in k-space. In brief we expand the wave field in the periodic structure in Bloch-wave form. For a given frequency and for a given mode, the field is a sum of terms called Bloch harmonics (also termed spatial harmonics), with wave vectors differing by a reciprocal lattice vector. These Bloch harmonics are important in understanding the results obtained by applying Fourier analysis to the time-domain data.
In the case of coated isotropic surfaces, when the shear-wave velocity in the coating is lower than that in the substrate, two branches exist with acoustic polarization in the sagittal plane, known as Rayleigh-like waves and faster Sezawa waves [30,45]. (Similar considerations apply to propagation along symmetry directions in coated crystals.) Plots of the magnitude of the spatiotemporal Fourier transform |F(k, ω)| of the z-directed surface velocity as a function of k x and k y for our experimental and numerical results, shown in figures 4(a)-(d), show the effect of these two types of waves. (The Sezawa waves were referred to as leaky waves in [19], but the naming in the present paper is more appropriate.) The horizontal axes are expressed in units of k B , where k B = π/a = 0.79 µm −1 is the position of the first Brillouin-zone (BZ) edge. Rayleigh-like and Sezawa-like portions of the dispersion are evident [19,21,25]. The details show up more clearly in the simulations, which were obtained for a simulation time of 39.3 ns. The cross sections of |F(k, ω)|, representing the experimental dispersion curves (figures 4(a) and (b)), show modes up to ∼1.3 GHz. 10 Ignoring Bloch harmonics, lossless propagation is limited to an area on the f -k x or f -k y plots between two thresholds [46]. The upper threshold (upper dashed lines in (figures 4(c) and (d)) corresponds to the vertically-polarized shear-wave velocity in the Si substrate parallel to the surface in the [011] or [101] directions (v Si T ≈ 5850 m s −1 ) and the lower threshold (lower dashed lines in (figure 4(c) and (d)) corresponds to the Rayleigh wave velocity in the slowest material Cu (v Cu R ≈ 2160 ms −1 ). The curvature of the branches of the dispersion relation is fairly pronounced owing to the presence of the Cu-SiO coating [30,45]. To understand these curvatures, the dispersion relations of both Rayleigh and Sezawa waves were calculated with a semi-analytical method [47] for a Cu-SiO layer using effective material properties based on orthotropic symmetry [21,48], and also for a Cu layer, both on Si substrates (see figures 4(c) and (d), respectively). The 30 nm Au layer and the orientation of the Si substrate ([011] direction) were taken into account in   [21,25]. This results in each sheet of the dispersion surface in 11 The density of the layer was taken as the average of those of Cu and SiO, whereas the components of the anisotropic elastic constant tensor were calculated according to the method of [48], yielding six independent values. Strictly speaking, this effective medium model can only be relied upon when all dimensions of the structure are much smaller than the wavelength. 12 These models do not predict the observed faster y-directed velocity at frequencies below ∼500 MHz.
(k x , k y , ω) space possessing regions with both Rayleigh and Sezawa character depending on the values of k x and k y [25].
The evident second and third RW-SW avoided crossings at ∼900 MHz and ∼1.25 GHz, respectively, are the result of the interaction between the first and second Bloch harmonics, respectively, of these modes. The gaps shown shaded in figure 4(h) correspond to best estimates of the position of the first three RW-SW stop bands. The first such band agrees in position with results obtained elsewhere on an uncoated but otherwise identical sample [21,25]. The decreasing slope of the curves with increasing frequency results in a flattening of the avoidedcrossing angles and also explains the decrease in the RW-SW gap widths with frequency.
There are also avoided crossings of the same wave modes (Rayleigh-Rayleigh or Sezawa-Sezawa-see the crossing of lines with the same weight on the right hand side of figure 4(g)) below the first RW-SW stop band and between the first and second such stop bands at ∼700 and ∼800 MHz, and also above. (Because the corresponding gaps are small [25], we have not represented them in figure 4(g).) These avoided crossings appear at the boundaries of the BZs [21,45]. Separate stop bands for Rayleigh and Sezawa waves are formed at the corresponding frequencies (see figure 4(h)), but our frequency resolution is not sufficient to estimate their width. The first such stop band is the first Rayleigh-Rayleigh stop band [25], and occurs near 450 MHz, below the first RW-SW stop band.
The curves for |F(k, ω)| for the k y direction (figures 4(b) and (d)) show both Rayleigh and Sezawa branches below ∼800 MHz. Above this frequency the effect of Bloch harmonics from the x periodicity produces branches starting at thresholds of ∼850 MHz and ∼1.2 GHz, as identified in figure 4(f).  (5)) and experiment (row (6)-taken for the most part from [19]), as well as their interpretation with Bloch harmonics (row (8)). (Results from an analytical model, shown in row (7), are discussed in section 4.) The boundaries of the first BZ are indicated by the vertical dashed lines. The topology of the first three sheets of the surface-wave dispersion relation for our grating structure is discussed in detail elsewhere [25].

Constant-frequency planes
At 382 MHz (figure 3(h)), the Rayleigh mode forms an approximately elliptical ring. The faster Sezawa mode forms a rounded-square shape inside, it being more affected by the Si substrate owing to its larger penetration depth [30]. At this frequency the constant-frequency plane intersects with two sheets in (k x , k y , ω) space. Bloch harmonics of the Rayleigh ring are also visible, shifted in the ±k x directions by the reciprocal lattice vector G 0 .
At a frequency below but close to 458 MHz (figure 3(i)) these rings intersect, resulting in the previously-mentioned Rayleigh-Rayleigh stop band. The much larger RW-SW stop band is visible at 534 MHz in both simulation and experiment in figure 3(j), strongly affecting propagation in and near the x direction. The overall shape of the outer parts of constantfrequency contours, that arise from the Rayleigh branch, becomes more and more flattened through the effect of the periodicity as the frequency increases. This flattening affects the phonon-focusing in real space, in particular the X-shaped pattern. These effects are discussed in detail in section 4.
The RW-SW stop band ends around 611 MHz (figure 3(m)), and above this frequency the propagation of acoustic waves in the x direction again becomes possible. However, the amplitudes remain small in this direction compared to the y direction, as is evident in figures 2 and 3.
At 687 MHz (figure 3(n)), the band structure is complicated by the intersection with a third sheet in (k x , k y , ω) space [25], the full details of which remain to be elucidated. Although the RW-SW stop band is closed at 687 MHz, the Rayleigh portions of the constant-frequency contour remain discontinuous. The flattened parts parallel to the k x directions-that we shall later term 'the outer parts' of the Rayleigh branch-imply propagation parallel to the lines, as observed in real space. The curved parts, together with the corresponding Bloch harmonics, form loops at the boundaries of the first BZ, and are responsible for the remaining acoustic propagation over all angles. A similar, but more complicated splitting is evident for the Sezawa portions. Further research into the detailed shape of the dispersion sheets would require a higher frequency resolution than available in the present experiments or simulations.

Group-velocity surfaces
To understand the surface-wave propagation and focusing behavior produced by a localized impulsive acoustic source in real space, it is necessary to first derive the dependence of group velocity v g on direction [49]. This is described by the acoustic ray surfaces, referred to here as group-velocity surfaces, that can be calculated from the dispersion relation by a gradient operation v g = ∇ k (ω(k)). The group-velocity vectors point in the direction normal to the constant-frequency contours. Figure 5(b) shows, for a selection of frequencies, the angular dependence of the group velocity calculated from interpolated fits to the numerically simulated constant-frequency contours for the outer parts of the Rayleigh branch, shown in figure 5(a). At 382 MHz, where the constant-frequency contour is slightly elliptical, we obtain an elliptical group-velocity contour. However, at higher frequencies, the group-velocity contour breaks up into two round pockets that decrease in size as the frequency increases. Related pockets were reported by Tanaka et al [42] and by Narita et al [13] in theoretical investigations on oblique propagation in superlattices.
In order to understand the general topology of the low-order sheets of the dispersion relation and of the corresponding group-velocity surfaces, it is useful to develop a simple analytical model for the 1D phononic structure based on a plane wave expansion [9,50]. For the sake of simplicity we consider only one mode, the Rayleigh wave, which we describe by a scalar wave equation; in 2D, the acoustic disturbance f (here the z-directed surface velocity) in a structure with a spatially-varying sound velocity in the x direction can be written in the form of Mathieu's equation [51] ∂ where v 0 is the sound velocity, and an exp(−iωt) temporal dependence has been assumed. The quantity V (x), that takes the role of an effective potential, represents spatial variations in the sound velocity. Consider the case in which V (x) is a perturbing periodic sinusoidal potential in the form [51] V where α 1/v 2 0 and G 0 is the magnitude of the reciprocal lattice vector G 0 = G 0 i = 2πi/a. We shall seek an approximate solution based on the inclusion of two Bloch harmonics (with k x >0):  where a 0 and a 1 are constant coefficients of the same order of magnitude. Substituting equation (3) into equation (1), and using equation (2), one may obtain two simultaneous equations for a 0 and a 1 to first order in α. These have solutions provided that the secular determinant vanishes, i.e. In order to adapt this model to our data we first choose v 0 = v 0 (ω) to fit the intrinsic curvature of the numerical Rayleigh dispersion in the k x direction shown in figure 4(e). 13 Then the value of α = 3.4 × 10 −8 s 2 m −2 was chosen to give reasonable correspondence of the outer part of the Rayleigh branch with experiment; with this value, we find, for example, α 0 v 2 0 ≈ 0.5 for a typical frequency of f = 500 MHz. This implies that we should not expect very accurate agreement with the numerical approach.
The resulting predictions for some constant-frequency contours and the corresponding angular dependence of the group velocity are shown in figures 5(c) and (d) for the outer part of the Rayleigh branch. We have made use of symmetry to extend the predictions from −k B to +k B . One can see that, in spite of the simplicity of the model and the relatively large value of α 0 v 2 0 , the agreement with the constant-frequency contours is quite good, in particular at the higher frequencies. This can also be seen by comparing rows (5) and (7) in figure 3. There are, however, some differences. For example, the two Rayleigh Bloch harmonics first appear at a higher frequency (in fact at ∼530 MHz) in the analytical model. The agreement with the group velocity (see (b) and (d) in figure 5) is also reasonable. In particular, the two pockets formed in the group-velocity plot and their frequency dependence are reproduced fairly well. Because we chose the analytical model to fit to the outer part of the Rayleigh branch, the predictions for the Rayleigh-Rayleigh stop band are not fitted to experiment. 14 The analytical model allows a wider range of frequencies to be plotted for the group velocity than in the case of the numerical simulation, because of the difficulty of accurately extracting interpolated constant-frequency contours from the latter. Figure 6 shows the shape of the group-velocity surface in (v gx , v gy , ω,) space (where v gx and v gy are the x and y components of the group velocity) corresponding to the outer parts of the Rayleigh-wave constant-frequency contours. The surface transforms into one with two horns as the frequency is increased. This can be explained as follows: as the frequency goes up, the relevant constant-frequency contours approach straight lines (along the x direction). This means that the group velocity is nearly the same for all wave vectors, hence the horns. To estimate this single group velocity value, one should consider the group velocity of the lowest waveguide mode in an individual Cu line, which, in turn, approaches the Rayleigh velocity in Cu (v Cu R = 2160 m s −1 ). A second sheet of the group-velocity surface corresponding to the second sheet of the dispersion surface [25] is predicted by the analytical model. (See the parts of the Rayleigh branch corresponding to the closed loops at 534-687 MHz in row (7) of figure 3.) This groupvelocity surface also has a two-horned shape at frequencies higher than that, ∼800 MHz, at which the closed loops first touch, with a predicted cutoff for this sheet below ∼530 MHz (corresponding to the frequency at the top of the Rayleigh-Rayleigh stopband). Narita et al have reported related group-velocity contours that arise in superlattices [13,42]. We shall see in the next section that the first sheet is responsible for the X-shaped focusing pattern in real space. 13 We used a polynomial cubic fit in the form k x (ω) The present approach has the advantage of simplicity. 14 The predicted extent for the Rayleigh-Rayleigh stop band is from ∼430 to 530 MHz. In reality the stop band is much smaller, and lies near ∼450 MHz.

Phonon focusing
The effect of the phononic band structure on the directionality of the surface-phonon field excited by a point source in real space is closely related to the angular variation of the group velocity. This surface-phonon focusing effect is characterized by an enhancement factor A p (φ) defined by A p = 1/|dφ/dθ |, where θ is the polar angle defining the direction of k and φ is that defining v g . A p (φ) characterizes how the surface waves are focused in real space, and is calculated numerically.
The result of this procedure for the analytical model (for the first group-velocity sheet) is shown in row (3) of figure 3. The X-shaped focusing pattern and its dependence on frequency are reproduced well. The angular variation A p (φ) predicts very sharp features as a function of φ. The directions of strong focusing correspond to inflection points (points of zero curvature) in k-space, as shown in figures 5(a) and (b). The observed frequency dependence is a result of these inflection points moving to positions of lower |k x |-and hence to a higher angle |φ| (for k x >0)-as frequency increases. These phonon caustics first appear at the bottom of the first Rayleigh stop band when the dispersion relation first touches the first BZ edge, and continue to exist above the top of this stop band because of the continuation at these frequencies of the lowest sheet of the surface-wave dispersion relation [25]. Because of the finite wavelength, the focusing is smoothed out to some extent in our images. Infinitely far from the source the angular distribution should exhibit true caustics. (The effect of the second sheet of the group-velocity surface is smaller than that of the first, and below ∼800 MHz shows focusing primarily along the x-axis, partly precluded by the RW-SW stop band, so we have not presented its focusing effects here.) We conclude this section by pointing out some significant differences between phonon focusing in gratings and traditional phonon focusing in anisotropic crystals: [11] 1. In anisotropic solids, surface phonon caustics arise when the anisotropy is strong enough: there is a threshold for their onset [49]. In gratings they arise whenever there is a periodic perturbation, no matter how weak. 2. The caustics for gratings arise above some frequency threshold defined by the bottom of a stop band. 3. The caustics for gratings bound non-propagation sectors for frequencies in the associated stop band. 4. Group-velocity contours form closed loops rather than cuspidal structures typical for phonon focusing in non-dispersive anisotropic media. 5. In anisotropic solids, caustics are only observed for certain kinds of waves. They do not, for example, occur for longitudinal acoustic waves. In gratings the phenomenon is general and should be observable for waves of any kind.

Results and discussion for the 10 µm-period sample
It has previously been found by simulation that the stop bands are sensitive to the grating thickness in these flat-grating structures [52]. Here we investigate the effect of a different periodicity. Experiments and simulations were carried out with a sample of period a = 10 µm, corresponding to Si and SiO line widths of 5 µm. The longer spatial period lowers both the frequency of the relevant mode interactions and the scattering potential, allowing the study of further features of the phonon scattering in the periodic structure. The experimental conditions were the same as for the a = 4 µm case. In particular, the optical pump spot was again centered on a Cu line. Unless otherwise stated the simulation was conducted with a simulation time of 39.3 ns. In the xdirection, however, the frequency spectrum is uniformly weakened over a wide frequency range, as opposed to the appearance of clear stop bands; most frequency components are attenuated within a few periods (∼20-30 µm) in the a = 10 µm grating. Figure 7(f) shows sections of these spectra at a distance of 17.5 µm from the excitation. The reduction in amplitude is very clear for x-directed propagation. At 229 MHz, the Rayleigh mode again forms an approximately elliptical ring, whereas at 280 and 305 MHz, at frequencies sufficiently above the bottom of Rayleigh-Rayleigh stop band, we again find in the simulations the characteristic X-shaped focusing pattern. At 280 MHz the wave fronts in the x-y plane are also visibly affected by the RW-SW stop band, which is responsible for regions of reduced amplitude around the ±x directions far from the source. At 305 MHz this stop band is no longer present, but the reduced-amplitude effect remains. At this frequency the acoustic field within angles ∼ ± 30 • to the x axis in the x-y plane show a  different phase compared to the remaining directions, leading to well defined directions of low amplitude. This effect may arise because of the interference between wave fields dominated by the waveguiding effects in the y direction and those related to scattering from the periodic structure. Similar phenomena occur at higher frequencies. At 611 MHz, radiation to the bulk is evident in the x-z plane image, very similar to the case for the a = 4 µm grating at this frequency. Over the whole frequency range the expected decreasing wave penetration in the z direction with increasing frequency is evident in the x-z plane images. Waveguiding effects in the central Cu line for this a = 10 µm sample are much more apparent than for the a = 4 µm sample. At 458, 611 and 763 MHz a guided mode along the ±y directions in the center Cu line becomes clearly visible. In particular, at 763 MHz the pattern of this guided mode shows a more complex lateral acoustic field pattern associated with multiple guided modes. This is also evident in the animations.

Analysis in real space
The shorter wavelengths at 611 and 763 MHz result in an interesting effect for propagation in the ±x directions. Since the wavelength is shorter than the 5 µm line width, the line pattern becomes visible in the acoustic field as an alternation between the different wavelengths, with the shorter wavelengths in the Cu lines and the longer wavelengths in the SiO lines.

Dispersion relations for the k x and k y directions
The thickness of the a = 10 µm grating is identical to that of the a = 4 µm grating, and so the basic modes underlying the features of the dispersion relation not associated with the periodicity remain the same. The previously identified Rayleigh and Sezawa portions of the dispersion relations are visible in k x and k y directions in the cross sections of the dispersion relation, as shown in figures 9(a)-(d) together with the previous fits superimposed in figures 9(c) and (d). The reciprocal lattice vector has magnitude |G 0 | = 2k B = 2π/a = 0.628 µm −1 . The boundary of the first BZ, k B = 2π/a = 0.314 µm −1 , is consequently much smaller than for the a = 4 µm grating (see figure 9(e)).
The curves for |F(k, ω)| for the k x direction (figures 9(a) and (c)) can be interpreted with the help of Bloch harmonics, some of which are shown in figure 9(e). These are clearly visible in the simulation. As with the a = 4 µm grating, Bloch harmonics are also evident in the experiment, in particular at the higher frequencies, but again the resolution is better in the simulation. The first RW-SW stop band occurs at a lower frequency than in the previous case, here at ∼250 MHz where the dispersion curves are almost straight. This stop band is narrower than in the case of a = 4 µm because of the lower effective scattering potential for this case. Mixing of the Rayleigh and Sezawa modes is again present, as shown in more detail in figure 9(g). As with the case of figures 4(g) and (h), we used a longer time for the simulation (78.6 ns) to improve the frequency resolution. The maximum in the frequency spectrum for x-directed propagation is ∼500-600 MHz, as reproduced in figure 9(h) at a specific point on the sample 17.5 µm from the source. The effects of the first three RW-SW stop bands on the frequency spectrum at this point are relatively small because of the strong scattering.
The first three RW-SW stop bands are also marked in figure 9(h). Because of the limited frequency resolution their position is not as certain as for the case of a = 4 µm, and their marked positions and widths are just best estimates. In the case of the second RW-SW stop band, the Rayleigh branch does not show a decrease in intensity, possibly owing to the narrow width of the stop band, and we have inferred its presence from the gap in the Sezawa branch (see left hand side of figure 9(g)).    Further avoided crossings as a result of the interaction between the first and second Bloch harmonics as well as those resulting from the same wave modes (Rayleigh-Rayleigh or Sezawa-Sezawa) are visible. The different avoided crossings and the stop bands lie much closer in frequency and wave number compared to the a = 4 µm case. The Rayleigh and Sezawa Bloch harmonics are visible in the 1st BZ, similar to the case for a = 4 µm. The first and second Rayleigh-wave Bloch harmonics are especially strong.
The curves for |F(k, ω)| for the k y direction (figures 9(b) and (d)) show both Rayleigh and Sezawa branches below ∼400 MHz. Above this frequency the effect of Bloch harmonics from the x periodicity produces branches starting at thresholds of ∼400 and ∼700 MHz, as identified in figure 9(f).  (4)) and experiment (row (5)), as well as their interpretation with Bloch harmonics (row (6)). The boundaries of the first BZ are indicated by the vertical dashed lines. At 229 MHz both the Rayleigh and Sezawa branches form concentric rings, the latter just reaching the boundary of the first BZ; the effect of the anisotropic Si substrate is less obvious than for the case of a = 4 µm at 458 MHz owing to the lower frequency here. Bloch harmonics are also visible, particularly in the simulations. The previously observed mixing of the Rayleigh and Sezawa modes is again present in the region of the stop band at 280 MHz (a frequency accessed in the simulation only). As mentioned above, this stop band closes at 305 MHz. The simulations and experiment agree well at 229 and at 305 MHz. In the range 458-763 MHz several Bloch harmonics are visible. This can be understood with the help of the interpretation in row (6) of figures 8(g)-(l). Experiment and simulation also agree well in this range.

Constant-frequency planes
Over this frequency range of 458-763 MHz, horizontally-oriented, flat, outer contours in the constant-frequency planes responsible for the waveguiding are evident both in the simulation and in experiment. As previously noted in the context of the real-space data, the wider Cu line (5 µm) here results in these effects being stronger, and leads to the appearance of higher-order guided modes. These modes are visible at 763 MHz in the constant-frequency planes as extra horizontal lines.

Conclusions
In conclusion, we have investigated the anisotropy of the radiation pattern of gigahertz surface phonons excited at point sources in two different 1D PC structures composed of thin, flat Cu-SiO gratings of different pitches by a time-domain imaging method and by numerical simulations. The acoustic field patterns were first mapped at individual frequencies by temporal Fourier transforms, revealing in 2D (real) space the phonon-focusing effects of the band structure as well as waveguiding effects in the Cu lines.
In addition, the dispersion relations were obtained by spatiotemporal Fourier transforms in ω-k section form in the x and y directions as well as in 2D k-space. Many features of the surface-wave band structure could be extracted. Both experiments and simulations, showing broad agreement, exhibit Bloch harmonics and branches composed of two different wave modes: Rayleigh-like waves and Sezawa waves. The higher frequency-and wave-number resolution of the simulations reveal many features of the surface-wave band structure up to 1.5 GHz, in particular the width and position of the first three phononic stop bands arising from the mixing of Rayleigh and Sezawa modes in the 4 µm-period structure.
The numerical simulations for the 4 µm-period structure were processed to extract the angular dependence of the group velocity corresponding to the outer parts of the Rayleigh-wave constant-frequency contours. This revealed pockets in (v gx , v gy ) space symmetrically positioned on the v gy axis. The most important features of the band structure in the region of the first stop band could be reproduced using a simple analytical model involving a plane wave expansion. This model confirmed that the lowest-order group-velocity surface (v gx , v gy , ω) exhibits two horns. We showed using this analytical model that this leads to non-propagating sectors in real space bounded by phonon caustics.
More work remains to be done to better elucidate the topology of the higher-order dispersion sheets and group-velocity surfaces, and the form of the series of mini-gaps formed by avoided crossings in the dispersion relation of the same polarization (Rayleigh-Rayleigh or Sezawa-Sezawa). This will require a higher frequency resolution than in the present study.
The simple analytical model presented here is independent of the form of the 1D periodic structure. We have only considered surface acoustic waves on flat gratings, but the simple analytic model presented, depending only on the assumption of waves in a weak 1D periodic potential as in equations (1) and (2), is general, so that the conclusions concerning phonon focusing and the overall form of the group-velocity surfaces should apply also to acoustic waves in other forms of 1D structure such as corrugated surfaces or superlattices. Our approach could also be extendable to 2D and 3D phononic structures, where similar concepts of caustics arising from stop bands should apply.
Finally, our findings should hold for other forms of waves: caustics arising from the presence of stop bands in periodic materials should be observable, for example, for electromagnetic waves, surface plasmons, water waves and matter waves.
Having established the form of f (r, t), we proceed to take the spatiotemporal Fourier transform: When acoustic attenuation is negligible over the imaged region, as in the present case, ω j (k ) is real, that is ω j (k ) = ω * j (k ). Equation (A.7) can then be further simplified by multiplying both sides by exp(iω t)/(2π ) and integrating with respect to t. Replacing ω by ω for convenience in the final result, we obtain F(k, ω) = 1 2 j A j (k )C j (k)δ(ω − ω j (k )) + A * j (−k )C * j (−k)δ(ω + ω j (−k )) . (A.8) The delta functions in this equation show that the Fourier amplitude F(k, ω) is only non-zero when the combination (k , ω) satisfies the relation ω = ω j (k ) or ω = −ω j (−k ) for the branch j, and so the dispersion relation can be determined from F(k, ω). Moreover, the complex amplitude of F(k, ω) is proportional to the amplitude of the relevant Bloch harmonic C j (k). In other words, for (k , ω) satisfying the reduced-zone dispersion relation, that is ω = ω j (k ), the ratios between values of F(k + G, ω) for different G exactly match those between the values of C j (k + G). The amplitude (and phase) of the Bloch harmonics derived from F(k, ω) are therefore directly related to the coefficients C j in the Bloch expansion of equation (A.1), and so yield valuable information about the periodic scattering potential. The second term in equation (A.8) ensures the following relation holds: which is expected for the Fourier transform of the real function u(r, t) (see equations (A.4) and (A.5)). This implies that half of (k, ω) space contains all the information on F(k, ω). If ω j (k ) is chosen to be positive, the first term in the sum of equation (A.8) corresponds to ω > 0, whereas the second term corresponds to ω < 0. Because of the dual nature of the information contained in F(k, ω), it suffices to consider only the first term: For a lossless medium, one can also derive extra conditions resulting from time-reversal symmetry: by applying the time-reversal operation t → −t to equation (A.1), we get f k, j (r, −t) = G C j (k + G) exp{i((k + G) · r + ω j (k)t)}, (A.11) which should also be an eigenmode of the lossless wave equation (which obeys time-reversal symmetry). This can be compared with the eigenmode for −k: where substitution of −G instead of G has been carried out without affecting the sum (over all positive and negative G). The complex conjugate of this equation is The exponential terms exp(i(k + G) · r) in equations (A.11) and (A.13) are the same, so without losing generality one may choose f k, j (r, −t) = f * −k, j (r, t), which implies that ω j (k) = ω j (−k) and C j (k) = C * j (−k) in the lossless case.