Simulations of acoustic waves bandgaps in a surface of silicon with a periodic hole structure in a thin nickel film

Properties of elastic waves are particularly interesting in periodic media, characterized by periodic spatial distribution of density and elasticity coefficients. In such periodic media, a number of phenomena analogous to the behaviour of matter waves (electrons) in a periodic potential, described in terms of the band model, take place. Elastic waves propagate in the form of Bloch waves so that the vector of deformation is periodic in the reciprocal lattice,

$u⃗k⃗=u⃗k⃗+G⃗$⁠, so all inequivalent excitations are found in the first Brillouin zone. At the borders of this zone, a Bloch wave forms a standing wave and the group velocity of the wave is zero. Moreover, at the borders of the first Brillouin zone the energy bandgaps appear, i.e. there are the ranges of frequencies to which no wave vector corresponds in the dispersion relation, so no waves of such frequencies propagate in the medium. This property has been well-known in the band theory of electron waves but it holds true for all other waves in solid state. The periodic media in which elastic waves are studied in the aspect of the search for acoustic bandgaps are known as phononic crystals. Both bulk and surface waves are studied for over a decade, mainly in order to establish the full bandgap. This will allow to design more effective acoustic screens, waveguides for elastic waves,^1,2 acoustic filters and will allow for focusing of these waves.³ The phononic crystals of a lattice constant of an order of hundred nanometers can simultaneously be photonic crystals, which fosters the hope of getting interesting acousto - optical effects.

Because of the facility of their production, the most often studied phononic crystals are one- and two-dimensional ones, which means that the media are periodic in one or two dimensions,⁴ while in the other directions the medium is homogeneous. Literature gives information on investigation of bulk waves⁵ and surface waves⁶ in the media containing cylinders in a periodic lattices of different symmetries. Full energy gaps have been obtained thanks to a high acoustic contrast between the materials,⁷ rotation of the elasticity tensors of the cylinders,^5,6 taking into account piezoelectric properties^8–10 or by making use of local resonance effect.^1,12–14 The most often used simulation methods are that of finite elements, that of finite differences^2,15 and that of plane waves.^11,15 The results obtained by computation have been verified by experiments.^12,16–18

In this paper, propagation of surface waves in a two-dimensional phononic crystal is simulated by the finite elements method (FEM) and plane waves method (PWM). The medium, (Fig. 1), is an elastic half-space limited by a free surface x₃ = −h and is made of an anisotropic nickel film of the thickness h deposited on the (001) surface of silicon crystal. The nickel film is decorated with cylindrical holes of the depth equal to the nickel film thickness (h) and the radius r, arranged with the symmetry axis along the x₃ direction in a square lattice of the lattice constant a. The ratio of the volume taken by the holes to that taken by nickel is described by the filling factor f = πr²/a². The aim of the study was the check for the presence of acoustic bandgaps in such a medium at the border of the first Brillouin zone as a function of the filling factor.

The method of plane waves is based on the assumption that a solution of the wave equation in a periodic medium, the Bloch wave, can be found by a set of plane waves. For the surface waves in a periodic structure, this method has been described in detail by Tanaka and Tamura¹⁹ for cylinders of infinite depths. The system studied in our work is made of a silicon substrate coated with a thin nickel film of the height h with a periodic lattice of holes.

The solution of the equation of motion

where

$ρr⃗$ is the medium density at

$\vec{r}=(x_1,x_2,x_3)$

$r⃗=(x1,x2,x3)$⁠,

$c_{ijmn}\left(\vec{r}\right)$

$cijmnr⃗$ are the elasticity coefficients, and u_i are the components the displacement vector of the volume elements, is assumed to be a set of plane waves

The vector

$x⃗=(x1,x2)$⁠, and k₃ ≡ λ is a complex component of the wavevector in direction x₃, describing the decay of wave amplitude with depth.

In order to obtain the dispersion relation, the generalized eigenproblem, obtained from the equation of motion, should be solved (see Appendix  A):²⁰ 

with respect to λ and separately for the film and for the substrate. In this way we obtain a set of solutions in the form:

Where X_l = C_l for the film and X_l = S_l for the substrate, while

$x⃗=(x1,x2)$⁠. Then, selected solutions, ensuring that the wave is restricted to free surface, should be stitched with the help of border conditions requiring the equivalence of the displacement and stress at the interface. Together with the condition of disappearance of stress at the free surface of the layer, we get a set of equations (see Appendix  B), which can be expressed in the matrix form:

The set is satisfied if the matrix M determinant becomes zero. These solutions identify the pairs (⁠

$k⃗$⁠,ω) making the dispersion relation.

In order to check the dynamics of a continuous medium by the method of plane waves, a special application in C++ language was written. The calculations were performed using the libraries GSL and Linpack. In fact PWM is addressed to the systems of solid state- solid state and when applied to liquids it gives unreal results - transverse vibrations. However, it is possible to simulate the structures in which the cylinders are filled with vacuum by choosing the parameters to meet the condition ρ/c_ij → 0.^21,22

To carry out PWM calculations, the Fourier series should be limited to a finite number of terms. For our simulations we chose n = 169 terms to reproduce the medium, which was enough to get convergent results.

To perform simulation with the help of FEM, the acoustic module from the COMSOL Multiphysics program was used. The method is based on separation of the medium into small elements in which it is possible to approximate the solution of wave equation by a linear function, which permits a transformation of the differential equation into a set of algebraic equations.²³ To get a dispersion relation the following steps are taken:

• Bloch theorem is applied according to which it is possible to restrict the calculations to a single elementary cell,

• The depth of the elementary cell must be at least equal to the depth of the medium penetration by a surface wave,

• The bottom surface of the elementary cell is blocked as

  $\vec{u}=0$
  $u⃗=0$⁠,

• The surface waves are distinguished from other solutions by requirement that the geometric centre of vibrations is above 30% of the cell height, counted from the free surface.

As the depth of the surface wave penetration is proportional to its length, we have to omit the calculations near the centre of Brillouin zone as the cell must be very high near this centre, which extends the time of calculations and the results do not bring significant information. Moreover, the polarization of the surface waves is described by the parameter

Θ = 0 stands for the polarization in the sagittal plane (containing the wave vector and the vector normal to the free surface), while Θ = 1 stands for the polarization in the horizontal plane (Fig. 2).

The energy of surface wave can be permanently bound to the free surface only apart from the so-called sound cone,²⁴ defined by the curve of the slowest bulk wave on the dispersion relation. This curve corresponds to the frequency of wave propagating in the plane parallel to the surface: ω_|| = vk_||, where

$k‖=k12+k22$⁠. The bulk wave propagating out of the plane parallel to the surface has the frequency of

$\omega =v\sqrt{k^2_{||}+k^2_3}={\omega }_{||}{{\sqrt{k^2_{||}+k^2_3}}/{k_{||}}}$

$ω=vk‖2+k32=ω‖k‖2+k32/k‖$⁠, thus ω > ω_||, so ω > ω_|| and the dispersion curve is inside the sound cone. As the surface wave is a specific superposition of bulk waves, so to get the former from the latter the solution must have an imaginary component k₃. Then, the frequency of the surface wave is

${\omega }_s={\omega }_{||}{{\sqrt{k^2_{||}-k^2_3}}/{k_{||}}}$

$ωs=ω‖k‖2−k32/k‖$⁠, so ω_s < ω_||, and the dispersion relation for the surface waves is outside the sound cone.

For silicon without the nickel film and with the holes making a square lattice on the surface, no bandgap was obtained for surface waves.¹³ The same results was obtained for the structure of holes in the surface of crystalline nickel. Thus, we have decided to study the dispersion relations and acoustic bandgaps for the system of a square lattice of holes in a nickel layer on silicon crystal surface. The starting point of the simulations was the system of a uniform layer of nickel on silicon, and then the presence of holes was taken into account.

Fig. 3(a) presents the dispersion relation for surface waves for a homogeneous nickel film on silicon (001) surface with a sound cone, defined by the velocity of the slowest transverse bulk wave in silicon. On the surface of the homogeneous nickel film there is not only Rayleigh's wave but also Sezawa wave - a generalized Rayleigh's wave, both polarized in the sagittal plane, and a transverse horizontal Love's wave. The presence of these waves is characteristic for a thin film of acoustic impedance greater than that of the substrate.²⁵ As shown in the inset of Fig. 3(a), the velocity of Rayleigh's wave decreases with increasing thickness h of the nickel film (and increasing wavevector k), from the value characteristic of pure silicon to that characteristic of pure nickel. The velocities of Sezawa and Love waves also decrease and these waves occur in the limited range of kh, outside which they are coupled with bulk modes. Fig. 3(b) presents full dispersion relation in the reduced zone scheme. In the Γ → M, or [110] direction there are only two surface waves the Rayleigh's one and Love's one, along with a pseudo-surface wave (inside the sound cone). In the Γ → X and X → M directions all three surface modes propagate.

Because of the necessity of settling an appropriate height of the elementary cell in the FEM, the depth of surface waves penetration were determined for kh = 1 and kh = 0.34 (compare with Fig. 3(a)) in the [100] direction (Figs. 4 and 5). The Sezawa wave amplitude for kh = 0.34 (Fig. 5(c)) is nonzero even at a long distance from the surface because it is coupled with the bulk wave and its energy is emitted towards the inside of the medium. Apart from the above, at a distance of 1.5 of the wavelength, the amplitudes of surface waves almost disappear. Thus, the height of the elementary cell in the simulations was assumed as 1.5 of the wavelength at the smallest

$k⃗=π2ax̂1$ considered.

In order to determine the effect of a periodic lattice of cylindrical holes, the simulations of dispersion relations were made for different thicknesses of the nickel film h and different hole filling factors. For low filling factors as well as for a thick nickel film no acoustic bandgaps were found. Fig. 6 presents the dispersion relations for f = 0.6 and h/a = 0.3, h/a = 0.4 obtained by the method of plane waves and the method of finite elements. For a film with a periodic structure there are - among others - two surface modes polarized in the sagittal plane and a surface transverse horizontal mode. For a sufficiently high filling factor the full acoustic bandgap appears for surface waves (shaded area in Fig. 6).

The positions and widths of full bandgap versus the filling factor for h/a = 0.2, h/a = 0.3 and h/a = 0.4 are shown in Fig. 7. As follows from the figure, the position of the bandgap is shifted towards lower frequencies with increasing film thickness, maintaining the same lattice constant. The width of the bandgap in general increases with increasing filling factor, but irregularly because of the appearance of additional plane resonance bands for the large radius holes, as the regions between the holes act as resonance cavities for the waves of certain frequencies.²¹ For a particular lattice constant of 500[nm] the bandgap is in the range of hundreds of megahertz, e.g. it reaches 540[MHz] for thickness of a nickel layer of 150[nm] and a radius of holes of 240[nm].

For films of thicknesses from outside the range 0.2a ÷ 0.4a, bandgap also occurs but its width is much smaller than ones we showed. Moreover, it appear for the filling factor from a much narrower range.

Full bandgaps have been obtained for the surface waves propagating in the medium made of a thin nickel film with a periodic lattice of holes supported on silicon substrate. Such a structure is relatively simple to be obtained by lithographic methods. The width of the bandgaps had reached over 500[MHz] for the sample of the lattice constant 500[nm], which would permit experimental verification of the results by high-resolution Brillouin spectroscopy.

The results obtained by the two methods applied are consistent. On the other hand, in contrast to the finite element methods, the use of the plane wave method would not permit the investigation of more complex structures, e.g. multilayer systems in the aspect studied. This is a consequence of numerical technique limitation because of the necessity of work with large size matrices.¹⁵ 

We would like to thank Mr Bartlomiej Graczykowski for support in PWM and Comsol Multiphysics calculations.

The terms

$nG⃗,G⃗′im$ in the eigenproblem (3) for a square lattice are:

The set of equations is defined by the following expressions obtained on the basis of the assumed boundary conditions:

The upper indices S and B refer to the film and the substrate, while