Elastic waves in phononic monolayer granular membranes

The vibrational properties of out-of-plane elastic waves in hexagonal monolayer granular membranes were studied theoretically. The predicted propagation modes involve an out-of-plane displacement and two rotations with axes in the membrane plane. Shear and bending rigidities at the contact between beads were considered. Both the cases of freely suspended membranes and membranes coupled to a rigid substrate were analyzed. Dispersion relations and the existence of band gaps are presented and discussed for various contact properties. For freely suspended membranes with sufficient contact bending rigidity, it is shown that complete band gaps exist. The results obtained may be of interest for testing with acoustic waves the elasticity of recently developed granular membranes composed of nanoparticles (of interest because of their phoxonic properties) and more generally for the control of designing devices for membrane wave propagation.


Introduction
The study of the properties of phononic crystals and of the control of elastic/acoustic wave propagation by metamaterials has become a very active area of research, owing to their numerous interesting applications in, for example, high-frequency surface acoustic wave filtering [1,2]. There is also much interest among researchers in the study of acoustic metamaterials, related to new or enhanced wave processes including, for example, acoustic cloaking, negative refraction or subwavelength imaging [3][4][5][6]. Although devices involving plates or membranes have been widely used, most of the recent studies on metamaterials concern three-dimensional (3D) materials supporting bulk or surface acoustic waves [2], whereas very little has been reported on plate or membrane waves [7][8][9]. Recently, even at the smallest scale suitable for gigahertz acoustic applications, nanogranular media were produced, allowing simultaneously both photonic and phononic (phoxonic) functionalities [10]. For instance, technologies have been developed for building crystals consisting of nanoparticles [11][12][13][14][15], ordered monolayers of particles [16][17][18] and freely suspended membranes composed of a single layer of periodically ordered particles [19][20][21].
To describe the high-frequency (short-wavelength) vibrational behavior of such periodic structures composed of interacting discrete particles, the continuous approximation of the elasticity theory becomes inadequate and the discrete nature of the medium should be explicitly taken into account [22,23]. For most of the discrete structures, including atomic crystals, the particles can be considered as point masses interacting through central forces [22]. However, to model the elastic behavior of granular crystals composed of finite-size particles interacting through forces located at the contact points between particles [24][25][26][27], it becomes necessary to take into account the particle dimensions, the finite rotational inertia (rotational degrees of freedom) and the interactions through non-central forces. In particular, in addition to shear and longitudinal acoustic modes, rotational (micro-rotational) modes are predicted in such structures [28,29].
For more than a century, there exist theories for the description of long-wavelength acoustic wave propagation in the so-called micropolar media, the Cosserat continuum theory (micropolar elasticity) [30][31][32]. However, when wavelength is shortened and becomes comparable to the particle size, these theories are unable to describe the dispersive phononic properties. Nevertheless, the predictions of the Cosserat continuum theory (micropolar elasticity) [30][31][32], which incorporates in continuous elasticity the possible rotations of micro-mechanical elements 3 composing micro-inhomogeneous materials such as microcrystalline materials, can be used for a comparison to the low frequency limit of the theories for elastic wave propagation in ordered granular structures. This comparison provides an opportunity to identify micro-mechanical interactions contributing to particular elastic constants appearing in Cosserat-type theories of micro-inhomogeneous but macroscopically homogeneous and isotropic materials in addition to two constants of classical linear elasticity of micro-homogeneous isotropic materials [33][34][35].
It should be mentioned that, although a general theoretical approach for analysis of acoustic waves in discrete periodic systems has been presented earlier [36], the analysis of 2D-3D structures composed of finite size particles has been limited mostly to predictions of the properties of bulk longitudinal and shear acoustic-type elastic modes [37][38][39][40] and guided shear acoustic modes in macroscopically inhomogeneous granular crystals [41][42][43]. Only in recent years the rotational, coupled rotational/transversal and transversal/rotational modes have been included in analyses [44][45][46]. The analysis of 2D structures, which includes rotational degrees of freedom in particle motion, has been reported only for in-plane motion [33][34][35]. Out-ofplane motion has been analyzed only for surface acoustic waves on discrete structures of point masses [23] and in the framework of the Cosserat continuum models [47,48].
Here, we present the results of a theoretical investigation of the elastic modes in freely suspended granular membranes and in membranes interacting with a substrate. The out-of-plane motion incorporating rotations of particles is analyzed for the first time. The developed theory provides the basis for the application of elastic waves in the non-destructive evaluation and testing of granular membranes and in the testing of the adhesion of a particle monolayer to a substrate.

Model
We consider an infinite monolayer membrane composed of periodically ordered spherical particles, arranged in a hexagonal lattice ( figure 1(a)). The structure is supposed to consist of homogeneous spheres with radius R and is characterized by a lattice constant a = 2R. The corresponding first Brillouin zone is shown in figure 1(b), and the modulus of the reciprocal vector k is 4π/( √ 3a).

Shear rigidity of the contacts
In order to analyze the out-of-plane modes of the membrane, we model the shear forces at the contacts as springs, characterized by a constant shear rigidity ξ at the contact. The elongation of the springs introduces forces and momenta that induce the motion of the particles: out-ofplane displacement (u) and in-plane rotation angles ϕ and ψ (ϕ-rotation with the axis in the x-direction and ψ-rotation with the axis in the y-direction; figure 1(b)). We present the equations of motion for the central particle, which are given by the equation for the displacement u 0 of the center of the central particle, and the equations for the rotations ϕ 0 and ψ 0 : where m is the mass of the particle and, for the particular case of homogeneous spheres, the momentum of inertia is I = 2 5 m R 2 . The term δu i denotes the elongation of the spring at the contact between the central and i-particle, i.e. the relative displacement between 0 and i-particle at the contact point. The form of relative displacements δu i can be given in terms of particle displacements in the z-direction and rotations. For the six particles surrounding the central one in a cell of the lattice, relative displacements are explicitly given by the expressions where the expressions for the displacements of particle i in the contact point with particle j (u i/j ) are given in the appendix, and ϕ i and ψ i are the rotations of the i particle with the axes collinear to the xand y-directions, respectively. The motion equations can be written for the new variables = Rϕ and = Rψ (the arcs on the particle surface) and solved in the form of plane waves where (x 0 , y 0 ) are the coordinates of the central particle. We note that only the relative coordinates ( x, y) between the central particle and the neighbor particles are involved in the out-of-plane modes of the membrane. After the substitution of equations (4)-(9) into the system of equations of motion (1)-(3) and plane-wave substitution with the help of equation (10), the equation of motion can be written in the form Normalizing the distances to the lattice constant a, the elements of the matrix S are then given by where we have defined 2 = ω 2 /(4ω 2 0 ) with ω 2 0 = ξ/m and α = 1 4 k x a, β= The -k dispersion curves shown in figure 2(a) are obtained from det(S) = 0. The dispersion curves are depicted, following the KOMK path in the first Brillouin zone [2], in figure 1(b). In the new normalization (a = 1) these points correspond, in the reciprocal space, to O ≡ (0, 0), M ≡ (0, π/2) and K ≡ (π/6, π/2). The segments OM and OK correspond to waves propagating along symmetry directions. The segment OM corresponds to the y-direction. Due to the symmetry of the structure, any rotation of the propagation direction of nπ/3 (with n being an integer) relative to the z-axis provides no modification of the propagation properties. Accordingly, the segment OK is a propagation direction equivalent to the x-direction, i.e. it exhibits the same propagation properties. In figure 2(a), different propagation directions are separated by dashed lines and the nature of the modes is labeled. The plotted eigenvalues have been colored relatively to the eigenvectors that have been classified. Red curves correspond to pure rotation modes, blue curves to coupled displacement-rotation modes (with a predominance of displacement) and yellow curves to coupled displacement-rotation modes (with a predominance of rotation). The eigenmodes of the infinite membrane motion are consisting, in general, of three different components (the out-of-plane displacement u and two in-plane rotation angles and ). Nevertheless, for particular propagation directions one of the components is decoupled from the others. In the case of a propagative wave in the k x (respectively k y ) direction, a mode called -mode (respectively -mode) and shown in red in figure 2 appears uncoupled from the mixed modes (u + ) (respectively (u + )). In the case of waves propagating along the y-direction (0 < β < π/2), we note a zero-frequency (u + )-mode, whereas two other modes corresponding to u + and -modes propagate with frequencies 2 = 15 4 (1 − 7 15 sin 2 β) and 2 = 15 4 (1 − 1 3 sin 2 β), respectively. For the x-direction, a compact presentation of the dispersion relation can be obtained for the uncoupled -mode: 2 = 15 4 (1 − sin 2 α), where 0 < α π/3. For the particular plane vibrations with k x = 0, k y = 0, the three components are decoupled from each other: the two modes corresponding to rotations have the same eigenfrequency given by 2 = 15/4, whereas the displacement mode has zero frequency.
The frequencies satisfying det(S) = 0 for any propagation direction in the plane (α, β) are shown in figure 2(b). For the three obtained eigenvalues of S, the sixth-order symmetry around the z-axis is observed on these iso-surface plots.
The zero frequency u + -mode (non-propagative mode of zero group velocity) exists because there are such combinations of displacements and rotations of the particles which do not lead to the loading of the contact (to the elongation of the shear spring). In other words, in the model presented above, there are configurations of the displaced and rotated particles which have the same energy as the background configuration with zero displacement and zero rotation. The analysis demonstrates that the zero mode is allowed, particularly because in the model developed above (figure 1(c)), the rotations of two neighboring particles with the same angle but opposite directions do not cause any loading of the contact. In reality, it is the bending rigidity of the contacts that counteracts this type of motion ( figure 3(c)). Whereas the bending rigidity plays a negligible role in comparison with shear rigidity in unconsolidated granular crystals [44][45][46], it plays an important role when the contacts are consolidated (sintered).

Bending rigidity of the contacts
We model, in the following, the bending rigidity of sintered contacts of radius r Rθ R by including in the contact description (additionally to the shear spring) two springs with normal 8 rigidities ξ B (see figure 3(c)). They are oriented in the membrane plane and located at the lowest and highest points of the contact along the z-direction.
In this extremely simplified presentation of the physical nature of bending rigidity, even rotations of the neighboring particles in opposite directions but with the same angle will cause elongation of one of the springs and contraction of the other. The corresponding forces have non-zero momenta relative to the rotation axes and will try to make the particles return to their equilibrium state. For example, an additional (bending-related) momentum of forces acting on particle 0 because of the rotations of particles 0 and 2 relative to the y-axis can be described by In the proposed simple model, bending-related momenta are effective not only for contacts located on the rotation axis (e.g. for the contact between particles 0 and 2 in the case of their rotation along the x-axis). In the latter case, spin (torsional) rigidity of the contact [49] could play a role. As a consequence of the contact bending rigidity, which is a function of both ξ B and the dimension Rθ of the contact along the z-axis (where θ is the contact angular dimension; see figure 3(a)), there is a moment along the t i j -axis tangent to the contact between neighbor beads i and j in the x-y-plane. This moment exists even if the beads rotate relative to the direction t i j to exactly opposite angles χ i j = −χ ji ( figure 3(c)).
Here, the angle χ i j of rotation of the bead i relative to the axis t i j is equal to where α i j = α ji is the angle between the t i j -axis and the y-axis (α 01 = α 04 = π/3, α 03 = α 06 = 2π/3, α 02 = α 05 = 0). Due to bending rigidity, the sum of the bending momenta provides the following contributions to the momenta rotating the central beads ( j = 0) along the xand y-axes: The matrix S is modified in the following additional terms: where the non-dimensional bending parameter p B has been defined as p B = θ 2 2 ξ B ξ . By considering bending rigidity of the contact, the rotation around the x-axis influences the rotation around the y-axis and vice versa. The dispersion relations of the membrane are modified by the presence of bending rigidity, as can be noted in the modes corresponding to the propagation along the y-axis (OM segment), where a second coupled mode propagates even for small values of the bending parameter ( figure 4).
The frequencies of all the modes in the points M, K and O can be found analytically. Their analysis demonstrates that a forbidden frequency gap for wave propagation along the y-direction ( Combining the information about the mode frequencies in the critical points M, K and O, it can be predicted that a complete band gap opens above = 3/2 for p B > 3/5. Its width grows up to a maximum value 3/2 √ 15/2 when p B = 1 and retains the same value on increasing p B above 1. These theoretical predictions agree with the numerical results presented in figure 4. For small values of the parameter p B , the frequency of the lowest mode is proportional to the square root of p B , 2 2 p B sin 4 β/[1 − (7/15) sin 2 β], demonstrating that propagation of this mode is due to the bending rigidity of the contacts. Measurements of the velocity of this mode provide direct access to the evaluation of bending rigidity. Analytical expressions for some symmetry directions of propagation can be obtained. In particular, for waves propagating along the y-direction the rotation mode frequency is given by and for waves propagating along the x-direction, the frequency of the rotation mode is Dispersion curves for different bending parameter values are shown in figure 4 and verify the expressed analytical results. In figure 4 (and in the following), the color code (blue, green, red) corresponds to the order of appearance of the eigenvalues along the diagonal of the dynamical matrix S after a Schur decomposition. The Schur decomposition provides ordering of the eigenmodes (and consequently ordering of the associated eigenvalues) in accordance with the dominance of the component u, as a first criterion. If, for several modes, components u are equal, the ordering is obtained in accordance with the dominance of , and finally, a comparison of the component is used if necessary. This ordering is used to attribute the successive colors blue, green and red to the three ordered eigenvalues in figure 4.

Normal contact with a rigid surface
We have analyzed the out-of-plane modes of a free suspended membrane considering both the rotational degrees of freedom of the particles and the finite size of the contact between beads. We include now the interaction with an infinitely rigid substrate supporting the membrane. The interaction with the substrate gives rise to the existence of additional forces and momenta at  the contact between an individual particle and the substrate, which are modeled as springs with their corresponding stiffness constant.
An additional returning force appears due to the motion along the z-axis, related to the normal rigidity of the contact of a sphere with the substrate ξ N (the symbol refers to the interaction with the rigid substrate). This introduces modifications in only one term of S, given explicitly by where we have defined the adimensional parameter p N as p N = ξ N /ξ . The contact of the membrane with a rigid substrate has an influence on the dispersion relations. The effect of normal rigidity in the absence of bending rigidity is shown in figure 5: when the membrane is on a rigid substrate it is composed rather of interacting oscillators than of oscillating masses. This manifests, in particular, in the fact that the wave with k = 0, which corresponds to pure displacements without rotations, has now a non-zero frequency.
The growth of the frequency of the displacement-dominated lowest mode in the point O is described by 2 = p N and is independent of the bending rigidity of the inter-bead contacts.
As follows from our analytical evaluations and figure 5 above, the displacement-dominated mode in the case p N = 0 becomes the lowest (in frequency) mode for p B 7/15 in the point K and for p B 8/15 in the point M. The analysis demonstrates that at larger values of p B and non-zero values of p N the frequency of the lowest mode in these points depends only on the bead-substrate normal interactions: 2 = (3/2) 2 + p N and 2 = 2 + p N at the points K and M, respectively. Thus the continuous increase of p N for a fixed value of p B > 9/15, where the complete band gap exists in the free-standing ( p N = 0) membrane, will lead to closing of the complete gap. For 3/5 < p B < 1 the complete gap closes because the frequency of the first (lowest) mode at the point K becomes higher than the frequency of the second mode at the point M. This happens when p N = 3 4 [5 p B − 3]. For p B 1 the complete band gap associated with bending rigidity (see figure 4) closes when p N becomes larger than 3/2 because the frequency of the first (lowest) mode in the point K becomes higher than the coincident frequencies of the second and the third modes at the point O. In figure 6, we illustrate for p B = 1.5 the closing of this complete gap, which is caused by the increased normal rigidity of the contacts between the beads constituting the membrane and the substrate. While this complete gap closes, the opening of a low-frequency complete band gap is also observed. This low-frequency complete band gap exists only if both bending rigidity of contacts between beads and normal rigidity of bead-substrate contacts are simultaneously nonzero. The bending rigidity is necessary to open a local low-frequency band gap at the point M, while the normal rigidity is necessary to open a local low-frequency band gap at the point O. For small values of the rigidities p B and p N , these local gaps grow proportionally to corresponding rigidities (as 0 15 p B /4 and 0 p N , respectively). Thus, for p r 1 and p N 1, the complete low-frequency band gap is controlled by max(15 p B /4, p N ) and experimental determination of the complete gap could provide a measure either of the bending rigidity of the contacts between the beads or of the normal rigidity of the bead-substrate contacts.

Conclusions
The above-presented results on the analytical and numerical evaluation of the out-of-plane motion of free-standing membranes composed of a single layer of ordered spherical particles demonstrate that, above a critical value of the ratio of bending to shear rigidity of the interparticle contacts, a complete forbidden band gap for the elastic wave propagation exists. This complete band gap may disappear in the membranes supported by a substrate. Closing of the high-frequency ('optical') complete gap with increasing normal rigidity of the contacts between the particles and the rigid substrate has been predicted. Also, a low-frequency ('acoustical') complete band gap has been predicted to exist when the bending rigidity of contacts between beads and the normal rigidity of bead-substrate contacts are both important. More complex phenomena are expected if the shear rigidity of the contacts between the particles and the substrate is taken into consideration, because the shear rigidity of the particle-substrate contact couples the out-of-plane motion in the supported membranes to the in-plane motions. It would also be interesting to investigate in future the role of the spin (torsional) rigidity of the inter-particle contacts. The theoretical results presented above will guide the development of opto-acoustical methods for the non-destructive non-contact evaluation of the elastic properties 14 of free-standing and supported monolayer particle membranes [13-18, 20, 21], which are based on the all-optical monitoring of the elastic waves through both their generation and detection by lasers [50][51][52].