Zero-frequency and slow elastic modes in phononic monolayer granular membranes

Ultrasonics xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Ultrasonics journal homepage: www.elsevier.com/locate/ultras

Zero-frequency and slow elastic modes in phononic monolayer granular membranes Li-Yang Zheng ⇑, Hélène Pichard, Vincent Tournat, Georgios Theocharis, Vitalyi Gusev ⇑ LAUM, UMR-CNRS 6613, Université du Maine, Av. O. Messiaen, 72085 Le Mans, France

a r t i c l e

i n f o

Article history: Received 21 October 2015 Accepted 3 November 2015 Available online xxxx Keywords: Granular phononic crystal Zero-frequency modes

a b s t r a c t We theoretically study the dispersion properties of elastic waves in hexagonal and honeycomb monolayer granular membranes with either out-of-plane or in-plane particle motion. The particles interact predominantly via normal and transverse contact rigidities. When rotational degrees of freedom are taken into account, the bending and torsional rigidities of the intergrain contacts can control some of the phononic modes. The existence of zero-frequency modes, zero-group-velocity modes and their transformation into slow propagating phononic modes due to weak bending and torsional intergrain interactions are investigated. We also study the formation and manipulation of Dirac cones and multiple degenerated modes. This could motivate variety of potential applications in elastic waves control by manipulating the contact rigidities in granular phononic crystals. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction There has been a growing interest in investigating the propagation of elastic/acoustic waves in phononic crystals in the last decade [1,2]. Owing to Bragg scattering in spatially periodic media, phononic crystals possess exotic dispersion characteristics including, for example, frequency band gaps [3], negative refraction [4], subwavelength imaging [5], etc. Although studies of phononic crystals involving elastic behaviors have been reported in all dimensions, from 1D to 3D [6–8], most of the recent studies are focusing on different types of two-dimensional lattice that support bulk or edge modes [9–11]. For instance, in the nearly isostatic square/kagome lattices, by accounting for only the nearest neighbor central-force interactions, soft modes and zero-frequency bulk modes have been predicted [12,13]. More interestingly, when the kagome lattices are twisted, negative Poisson ratio and zerofrequency edge states could be achieved [14,15]. Other fascinating elastic properties, such as topological soft modes and topological edge modes, have also been reported in the kagome and honeycomb systems [16,17]. In non-consolidated granular crystals, the interactions between individual grains take place via local contacts, which are much smaller in size than the dimensions of the individual grains and ⇑ Corresponding authors. E-mail addresses: [email protected] (L.-Y. Zheng), vitali.goussev@ univ-lemans.fr (V. Gusev).

inherently much softer than the grains [18–22], e.g., Hertzian contacts. Even when granular crystals are consolidated via their curing, like opals, or by grain-connecting ligands, like in nanocrystal superlattices, the elastic links between the grains keep being significantly smaller and softer than the grains themselves. This induces propagation of elastic waves in granular structures at significantly slower velocities than in the individual grains [23–25] even if the rotational degrees of freedom of the individual beads are not strongly involved. In contrast to normal forces, which are central forces in most of the spring-mass systems [14–17], in granular crystals the shear forces due to transverse rigidity of the contacts are non-central and can initiate the rotation of the beads. Thus the rotational degrees of freedom of the individual beads, the particle dimensions and the interactions through non-central forces should be taken into consideration. It was theoretically predicted [6,18–20] and experimentally demonstrated [21,22] that, due to the rotational degrees of freedom of the particles, additional coupled rotational/transverse and pure rotational modes can propagate in granular crystals, while the pure transverse modes, predicted by the theoretical case of the frozen grain rotations, are modified into coupled transverse/rotational ones. Interestingly, accounting for the rotational degrees of freedom, in addition to translational ones, can also lead to the existence of zerofrequency (zero-energy) modes [26–28]. The additional rotational degrees of freedom provide extra flexibilities to dispersion engineering of phononic crystals and to the control of the propagation of elastic waves.

http://dx.doi.org/10.1016/j.ultras.2015.11.005 0041-624X/Ó 2015 Elsevier B.V. All rights reserved.

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

2

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

It should be mentioned that, a general theoretical approach for the analysis of acoustic waves in phononic crystals is known already for quite a long time [1–3], and the Cosserat continuum theory for the description of the long-wavelength acoustic wave propagating in the micropolar media with rotational degrees of freedom exists for more than a century [29–31]. In contrast, the analytical discretized models that study the rotational modes and their coupling to other modes, in particular to shear ones, have started to attract increasing interest only in recent years [6,18–20,26,32]. In this work, we exploit the Lagrangian method [33] to evaluate the intergrain interactions of granular phononic membranes. An important advantage of the method is that it is possible to find analytical solutions for modes at high symmetry points, like C, M and K. The analytical formulas are very useful by giving clear guidelines on how to control the phonons spectra and to design suitable phononic crystals. In particular, we theoretically study the dispersion relations of elastic waves in hexagonal and honeycomb monolayer granular membranes for both out-of plane and in-plane motion. We also demonstrate that rotational modes and their coupling to translational modes can provide more flexibilities and additional functionalities in the control of the elastic wave propagation. Specifically, the granular phononic structures are expected to be advantageous in the monitoring of bulk shear and surface Rayleigh acoustic waves [20,34]. Besides, the detailed analysis of zerofrequency modes is reported. In the honeycomb lattices, we demonstrate the existence of zero-group-velocity rotational modes with non-zero frequency, the propagation of which can be initiated by weak bending and torsional interactions between the beads. We also study how the number and parameters of these modes are changing in transition from hexagonal to honeycomb lattice. Finally, we predict the degenerated modes at C and K points that can be realized for particular values of bending and torsional rigidities. For example, when the bending/torsional rigidities have special values, the Dirac-like cones, triple degenerated points and even the double Dirac cones can be obtained. As reported in many previous publications, manipulation of Dirac cones could lead to many interesting effects [35–39]. For instance, by breaking the symmetry of the system, the opening of the gap in the Dirac K point in acoustic/elastic systems can give rise to the topological edge states propagating only along some particular directions [40–42]. Even for the Dirac cone at C point, one could expect to observe the pseudospin-resolved Berry curvatures of photonic bands and helical edge states characterized by Poynting vectors [43,44]. The analytical predictions of Dirac cones and double Dirac cones in this work could largely facilitate the study of the topological properties of elastic waves in more complex granular membranes with modified/broken symmetries. The study of these types of membranes is motivated by potential applications of the unidirectional waves propagating with reduced/avoided attenuation/scattering. In general, we believe that our theoretical analysis of the elastic waves in mechanically free membranes would be useful also in the studies of the interaction of the granular layers with the elastic substrates [45–48]. This paper is constructed as follows: in Section 2, the structures of the studied membranes and the interactions between beads are analyzed. The theoretical calculation and analysis of the modes in hexagonal monolayer membranes with out-of-plane motion is presented in Section 3. This analysis includes the phonon spectra, the zero-frequency modes, the degenerated modes and the Dirac cones at the high symmetry points. Then we are focusing on in-plane motion in hexagonal monolayer membrane in Section 4. In Sections 5 and 6, we turn our attention to the honeycomb monolayer membranes for both the case of out-of-plane and in-plane motions. In Section 7, we present the conclusions of this work.

2. Structures and intergrain interactions As shown in Fig. 1, the two-dimensional infinite monolayer membranes under consideration are composed of periodically ordered spherical particles with radius R, arranged in a hexagonal and a honeycomb lattice. The structures are characterized by the pffiffiffi lattice constant a ¼ 2R for the hexagonal lattice and a ¼ 2 3R for the honeycomb lattice. The corresponding first Brillouin zones are also depicted in Fig. 1. Considering different types of motion in these monolayer membranes, the following degrees of freedom are taken into account: (1) For the out-of-plane motion, the beads in the membrane exhibit out-of-plane displacement ðuÞ along z-axis and in-plane rotational angles u and / (u-rotation with the axis in the x-direction and /-rotation with the axis in the y-direction). (2) For the in-plane motion, the beads in the membrane possess out-of-plane rotation (u) along z-axis and in-plane displacements ux and uy along x-axis and y-axis, respectively. The dynamics and the coupling of these mechanical motions are controlled by the following forces and/or moments (see Fig. 2): (1) Shear forces, which are characterized by an effective shear rigidity ns (Fig. 2(a)). These forces are activated in the membrane due to a resistance of the contact to relative displacement of the beads in the direction transversal to the axis connecting their centers and due to in-phase rotation of the beads relative to the direction normal to the axis connecting their centers. (2) Torsional forces, which are characterized by an effective torsional (spin) rigidity nt (Fig. 2 (b)). The resistance of the contact to relative rotation of the beads along the axis connecting their centers can initiate these forces. (3) Bending forces, which are characterized by an effective bending rigidity nb (Fig. 2(c)). They originate from the resistance of the contact of beads to rolling. (4) Normal forces at the contact between two adjacent particles described by normal rigidity nn (Fig. 2(d)). This type of interaction can be excited when there is relative displacement between two adjacent beads along the axis connecting their centers. For the out-of-plane motion, the motions of beads in the membranes lead to the shear, torsional and bending interactions, while the normal forces are not initiated. For the in-plane motion, the normal, shear and bending interactions are activated, while the torsional interactions are not. In Fig. 2(a) the effective shear rigidity of the contact is represented by a single spring of an effective rigidity ns and the energy of interaction can be evaluated because the relative macroscopic displacements of the neighbor beads are known. The interaction energy is proportional to the product of the effective shear rigidity and the square of the relative displacement which is equal to the elongation of the effective spring. However, one should keep in mind that the shear interactions are in fact distributed at the complete surface of the contact. To appreciate the role of the finite dimensions of the contact, we split the effective shear spring into two independent springs of half rigidity separated from spatially not only themselves but also from the center of the contact. The characteristic distance of the separation is of the order of the contact radius d. This presentation of the shear interaction between the contact faces does not change the magnitude of shear interaction, however, it reveals the existence of the torsional rigidity of the contact. From Fig. 2(b) it is clear that, when the neighbor beads exhibit unequal rotations relative to the axis connecting their centers, one of the springs elongates while the other shrinks. Thus, it is the induced forces acting on the beads that are totally compensated being combined destructively, but not the moments. The moments due to the deformation of two beads will be added constructively, introducing resistance to torsional motion, which we describe as torsional interaction (Fig. 2(b)). It is worth noting that the elastic energy stored in shear interaction, presented in Fig. 2(a) is  ns ðu2  u1 Þ2 . The elastic energy stored in shear

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

3

Fig. 1. Structures of two granular membranes and their Brillouin zones. (a) Hexagonal membrane. (b) Honeycomb membrane. The directions of coordinate axes are also shown.

Fig. 2. Schematic presentation of shear and torsional interactions by resistance to relative motion of the beads of (a) a single central shear spring and (b) two separated noncentral shear springs. Schematic presentation of normal and bending interactions by resistance to relative motion of the beads of (c) two separated non-central normal springs and (d) a single central normal spring.

interaction caused by the rotations in Fig. 2(a) is  ns ½Rðw2 þ w1 Þ2 , where wi ði ¼ 1; 2Þ is the rotational angle shown in Fig. 2(a). Thus, it is the product of the bead radius and the angle that plays the role of the effective coordinate when we characterize the elongation of the effective springs caused by rotations. Then, the torsional elastic energy stored in the motion presented in Fig. 2(b), i.e.,  ns ½dðw2  w1 Þ2 , should be presented as  ns ðd=RÞ2 ½Rðw2  w1 Þ2 , demonstrating that the effective torsional rigidity of the contact is related to its shear rigidity by nt  ns ðd=RÞ2 . Thus for contact areas with small radius in comparison with the beads radius R,

the effective torsional rigidity is much weaker than the effective shear rigidity. From the physics point of view, this is a manifestation that torsional interactions (and also bending interactions described below) take place via the moments of forces and not through the forces. The moments can be small even for high amplitude forces if these are applied too close to the axis of rotation. In Fig. 2(c), in order to reveal the existence of bending interaction, we have splitted the normal rigidity of the contact (presented in Fig. 2(d)) into two separated normal springs of half rigidity. The elastic energy of the normal interactions is the same for both cases, i.e., those schematically illustrated in Fig. 2(c) and (d). However,

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

4

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

the spatially separated springs in Fig. 2(c) activate additionally the bending interaction as they counteract the rotation of the beads in opposite direction relative to the direction normal to the axis connecting their centers. Arguments similar to those presented above, where shear and torsional effective rigidities are compared, lead to the following relation between the effective bending and normal rigidities, nb  nn ðd=RÞ2 . The above derived estimates indicate that torsional and bending rigidities could be much smaller than shear and normal rigidities, respectively, i.e., nt  ns and nb  nn , in non-consolidated granular crystals, where commonly d  R and also if these small contacts are consolidated for example by curing. However, potentially it is possible to create interbead contacts with effective dimensions comparable to the dimensions of the beads by linking the beads with chemical ligands at micro/nano scale or by elastic rods at macro scale. Thus the granular crystals with nt P ns and nb P nn cannot be a priori excluded from the theoretical analysis. Because commonly normal rigidity is comparable with shear rigidity, i.e., nn P ns , then the granular crystals with nb P nt cannot be excluded from the theoretical analysis either.

2 Sð1;1Þ out ¼ sin ð2qx Þ þ sin ðqx þ qy Þ þ sin ðqx  qy Þ  X ; 2

Sð1;2Þ out

2

2

pffiffiffi

3 ¼ i sin 2ðqx þ qy Þ  sin 2ðqx  qy Þ ; 4

Sð1;3Þ out ¼  Sð2;2Þ out ¼

Sð2;3Þ out ¼ Sð3;3Þ out ¼

i 2 sinð4qx Þ þ sin 2ðqx þ qy Þ þ sin 2ðqx  qy Þ ; 4

i 3 3  3gb  gt h 2 2  sin ðqx þ qy Þ þ sin ðqx  qy Þ 2 4 1 2 þ gt sin ð2qx Þ  X2 ; P pffiffiffi

3 ð1  gb þ gt Þ cos2 ðqx  qy Þ  cos2 ðqx þ qy Þ ; 4 i 3 1  gb  3gt h 2 2  sin ðqx  qy Þ þ sin ðqx þ qy Þ 2 4 1 2  ð1  gb Þ sin ð2qx Þ  X2 ; P

ð6-1Þ ð6-2Þ ð6-3Þ

ð6-4Þ

ð6-5Þ

ð6-6Þ

3. Hexagonal monolayer membrane with out-of-plane motion

ð1;2Þ Sð2;1Þ out ¼ Sout ;

ð6-7Þ

3.1. Theory

ð1;3Þ Sð3;1Þ out ¼ Sout ;

ð6-8Þ

To analyze the out-of-plane motion in hexagonal membrane, we use the Lagrangian formalism [33]. Thus the Lagrangian of the central bead can be written as follows:

ð2;3Þ Sð3;2Þ out ¼ Sout ;

ð6-9Þ

L0 ¼ Ekin0  Epot0 ;

ð1Þ

where the kinetic energy Ekin0 and the potential energy Epot0 read as:

Ekin0 ¼

1 1 2 1 _2 _ þ I/ ; Mu_ 20 þ Iu 2 2 0 2 0

ð2Þ

Epot0 ¼

 1 X 2 ns Ds2i þ nb Dbi þ nt Dt2i ; 2 i

ð3Þ

3.2. Dispersion curves and complete band gaps

u0 is the displacement of the central bead, u0 and /0 are the inplane rotations along x- and y-direction, respectively. M is the mass of the beads and I is the momentum of inertia. For the case of homogeneous spheres, the momentum of inertia is I ¼ 2MR2 =5. Dsi ; Dbi and Dt i are the elongations in length of the effective springs, which are used to model the elastic response of the intergrain contact of the zero-th bead with its i-th neighbor and which activate the shear, bending and torsional forces, respectively. The explicit expressions of Dsi ; Dbi and Dt i are given in Appendix A. The Lagrange equation is [33]:

  @L d @L ¼ 0;  @q dt @ q_

ð4Þ

where q ¼ ðu; u; /Þ denotes the generalized coordinates and   _ u _ ; /_ the time-derivative of these coordinates. From the q_ ¼ u; Lagrange equation (Eq. (4)) one obtains the equations of motion of the central bead (see Appendix B). By looking for plane wave *

solutions of the form v ðt; x; yÞ ¼ v eixtikx xiky y , the dynamical equation can be obtained:

 * Sout qx ; qy v out ¼ 0; *

pffiffiffi where P ¼ MR2 =I, and qx ¼ kx a=4; qy ¼ 3ky a=4 are the normalized wave vectors of kx ; ky ; X ¼ x=2x0 is the normalized frequency pffiffiffiffiffiffiffiffiffiffiffi with x0 ¼ ns =M . gb ¼ nb =ns ðgt ¼ nt =ns Þ denotes the ratio of bending (torsional) to shear rigidity. The solution of this eigenvalue problem gives the X  k dispersion curves. It is worth noting that the description of the out-of-plane motion of the hexagonal membrane extends the theory developed in Ref. [19] by accounting for the role of the torsional (spin) rigidity of the contacts, which was neglected in Ref. [19].

ð5Þ

where v out ¼ ðu; U; WÞT is a column eigenvector with U ¼ Ru and W ¼ R/ and the superscript T represents the transposed vector.   The elements of the matrix Sout qx ; qy are given by:

In this section, the X  k dispersion curves for different values of the ratios gb and gt are derived. It will be shown that it is possible to create a complete band gap by a fine tuning of the contact rigidities. Without losing of generality, in Fig. 3 we show the dispersion curves, following the path KCMK of the first Brillouin zone, for the case of a fixed gt ¼ 0:1 and varying the gb . Each of the eigenmodes consists of three components, the displacement u, and the rotations U; W. The nature of the modes is labeled in different colors. Red curves correspond to the pure rotation modes, blue curves represent the coupled translational/rotational modes with a predominance of translation, while yellow curves depict the coupled translational/rotational modes with a predominance of rotation. All the modes at the high symmetry points C, M and K are decoupled to pure modes and they can be described analytically by the dynamical matrix. The frequencies of some of those are marked in the dispersion curves. Concerning the propagation along y-axis (CM segment), the eigenfrequency of the pure displacement u pffiffiffi mode at M point is fixed at Xu ¼ 2, while there is also a pure rotapffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tional U mode at XU ¼ 5ð3gb þ gt Þ=2. Since the latter depends on gb and gt , a forbidden band gap along this direction can exist and be tuned. In particular, as long as 3gb þ gt 6 8=5, this band gap exists pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi at 5ð3gb þ gt Þ=2 6 X 6 2 (see Fig. 3(a) and (b)). In the particular case of 3gb þ gt ¼ 8=5, the gap disappears (Fig. 3(c)), while for pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8=5 6 3gb þ gt 6 3, it reappears at 2 6 X 6 5ð3gb þ gt Þ=2 pffiffiffi pffiffiffiffiffiffi (Fig. 3(d) and (e)), and it stabilizes at 2 6 X 6 15=2 when

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

5

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

Fig. 3. Dispersion curves of the hexagonal membranes with out-of-plane motion when gb ¼ 0:1 while gt is varying. Red lines represent purely rotational modes. Blue (yellow) curves correspond to coupled translational/rotational pffiffiffiffiffiffi modes with the predominance of translation (rotation). A complete band gap appears above X ¼ 3=2 when 3gb þ gt P 9=5. Its width stabilizes at 3=2 6 X 6 15=2 when 3gb þ gt P 3. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3gb þ gt P 3 (Fig. 3(f)). Now, for waves propagating along the xdirection (CK segment), the eigenfrequency of the two degenerated pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi modes at K point is noted at X ¼ 15ð1 þ 3gb þ 3gt Þ=4 (a pure U mode degenerates with a pure W mode). Since this mode also depends on the values of gb and gt , a forbidden band gap along this direction can exist and be tuned. In particular, for gb þ gt P 7=15, as shown in Fig. 3(e) and (f), the band gap is located at pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3=2 6 X 6 15ð1 þ 3gb þ 3gt Þ=4. This gap has maximumwidth pffiffiffiffiffiffi 3=2 6 X 6 15=2 when the values of gb and gt are large ð3gb þ gt P 3Þ. Finally, for waves propagating in the MK direction, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a forbidden gap appears at 3=2 6 X 6 5ð3gb þ gt Þ=2 when 3gb þ gt P 9=5. Combining the information on modes at the C, M and K points, a complete band gap can be predicted, for example, in Fig. 3(e) and (f), above X ¼ 3=2 for 3gb þ gt P 9=5. Its width pffiffiffiffiffiffi reaches a maximum value 3=2 6 X 6 15=2 when 3gb þ gt P 3. 3.3. Zero-frequency modes For the extreme case where the impacts of bending and torsional forces are ignored (gb ¼ 0 and gt ¼ 0), for example, assuming ðd=RÞ2  1, the dynamical matrix reduces to:

23

 cos 2qx cos 2qy  12 cos 4qx  X2 pffiffi 6 0 S ¼6  23 i cos 2qx sin 2qy 4    2i sin 2qx cos 2qy þ sin 4qx 2

2

2 sin qy  X2 pffiffi 6 3 S0 ¼ 6 4 i 2 sin 2qy 0

i 3 2

pffiffi 3 2

sin 2qy   X2 3 cos 2q cos 2q x y 1  P 4 pffiffi  43 sin 2qx sin 2qy

sin 2qy 2

0

3

0

cos2 qy  XP

7 7: 5

0 1 ðcos2 2

Thus, the zero-frequency modes are u þ U modes, i.e., coupled translational/rotational modes with zero group velocity (nonpropagative modes). Although it has been predicted before in many systems, the physical explanation of the existence of the zerofrequency mode in granular crystal has seldom been discussed [18–20,26,32]. Fig. 4(b)–(d) demonstrate the movements of beads in three different positions along the CM direction. At C point, beads do not rotate while their displacements in z-axis are the same. It suggests that although beads could displace from their original places, the relative displacements of the neighbor beads are absent. Consequently, the displacements of the beads do not lead to the loading of the contacts. Thus the system has the same energy as the one with zero displacement and zero rotation. For the case of the point in the middle of CM, some beads rotate, others have only displacements, which implies that there could be loadings of the contacts between two adjacent layers. However, the rotations are compensated by the translations in a sense that the displacements along z-axis of the opposite sides of the interparticle contact are equal, even though the displacement of one side is the

pffiffi 3 i cos 2qx 2

The dispersion curve is presented in Fig. 4(a) [19]. In the case of waves propagating along the y-direction, a band of zero-frequency modes is noted along the CM direction ðqx ¼ 0; 0 < qy < p=2Þ. The   dynamical matrix S0 qx ; qy along this direction is reduced to the form:

2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 5 cos2 qy þ 2 =2 and two qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 u þ U coupled modes at X ¼ 0 and X ¼ 15ð1  7 sin qy =15Þ=2. Eq. (8) predicts a pure W mode at X ¼

qy þ 2Þ 

X2 P

ð8Þ

i 2

 1 4



sin 2qx cos 2qy þ sin 4qx pffiffi  43 sin 2qx sin 2qy

3

 

cos 2qx cos 2qy  1 þ 4 cos 2qx  2

7 7: 5

ð7Þ

X2 P

result of purely translational motion while the displacement of the other side is due to purely rotational motion. Thus, the contact remains unstrained. In the M point, the motion is dominated by rotations, where beads have the same translations with small amplitude and the contacts between the subsequent layers of beads in y-direction are unstrained because the subsequent layers of beads exhibit equal rotations but in opposite directions. For the other wave vectors along the y-axis the situation is similar, i.e., there is a combination of rotations and translations which keep the contacts between the particles unstrained and the energy of the

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

6

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

Fig. 4. (a) The zero-frequency mode along CM, propagating in y-axis, in hexagonal membrane with out-of-plane motion and corresponding movements of the beads at three different points. (b) C point: kx ¼ 0; ky ¼ 0; (c) The point of half CM: kx ¼ 0; ky ¼ p=4R; (d) M point: kx ¼ 0; ky ¼ p=2R.

system unmodified. Thus, to create this type of motion costs zero energy, i.e., zero-energy modes. In reality, in the considered granular membrane, it is the bending/torsional rigidities of the contacts that counteract the zero energy combinations of translations and rotations. Although the bending and torsional rigidities play negligible role in the description of modes whose propagation is mainly controlled by the shear rigidity of contacts, they play crucial role in activating the propagation of zero-energy modes where the interparticle contacts are unsheared due to the described above compensation of the translations and rotations (compare Fig. 4(a) to Fig. 3). It should be mentioned that the bending/torsional interactions, as shown in the following, can also initiate the propagation of non-propagating modes with non-zero frequency.

optic/acoustic effects [35–44]. At K point, namely by setting   qx ¼ p=6 and qy ¼ p=2, the dynamical matrix S qx ; qy gives: (1) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi an u mode at Xu ¼ 3=2, (2) a U mode at XU ¼ 15ð1 þ gb þ gt Þ=4, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and (3) a W mode at XW ¼ 15ð1 þ gb þ gt Þ=4. Thus, this Dirac cone originates from the degeneracy of two rotational modes at frepffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi quency X ¼ 15ð1 þ gb þ gt Þ=4. Basically, as shown in Fig. 5(a) and (b), the variation of the values of gb and gt cannot break the Dirac cone, but can only change its position in the dispersion diagram. As depicted in Fig. 5(c) and (d), by controlling gb and gt , the u mode could perfectly degenerate with the rotational modes when the relation gb þ gt ¼ 7=15 is satisfied, contributing to a triple degenerated mode at the K point. 4. Hexagonal monolayer membrane with in-plane motion

3.4. Degenerated modes and Dirac cone 4.1. Theory Due to symmetry [49,50], two bands are intersecting at K point forming a Dirac cone where the dispersion curves become linear around K point. Dirac cones have been reported in many systems and the wave physics around them has explained a lot of interesting

By analyzing the in-plane motion of the central bead in hexagonal membrane, the kinetic energy Ekin0 and the potential energy Epot0 are written in the forms:

Fig. 5. Dispersion curves in hexagonal membrane with out-of-plane motion for different values of bending and torsional rigidities. Dirac cones and triple degenerated modes at the K point noted by green arrows. The existence of Dirac cone cannot be broken by modifying the values of bending and torsional rigidities as shown in (a) and (b). Different combinations of bending and torsional rigidities can lead to triple degenerated modes as shown in (c) and (d). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

Ekin0 ¼

Epot0

1 1 1 2 _ ; Mu_ 2x0 þ M u_ 2y0 þ Iu 2 2 2 0

6   1X 2 ¼ nn Dn2i þ ns Ds2i þ nb Dbi ; 2 i¼1

ð9Þ

ð10Þ

where Dni is the elongation in length of the spring between the zeroth and the i-th bead that activates the normal forces. Substituting Eqs. (9) and (10) into Eqs. (1) and (4), one obtains the equations of motion of the central bead (see Appendix B). The dynamical equation for the in-plane motion could be found by using the solution form of plane waves:

 * Sin qx ; qy v in ¼ 0;

ð11Þ

 T   * with the eigenvector v in ¼ ux ; uy ; U . Sin qx ; qy is the dynamical matrix of the following form: ð1;1Þ

Sin

ð1;2Þ

Sin

ð1;3Þ Sin

ð2;2Þ

Sin

ð2;3Þ

Sin

¼

¼

g þ 3 4

 g 1  cos 2qx cos 2qy þ ð1  cos 4qx Þ  X2 ; 2

pffiffiffi 3ðg  1Þ sin 2qx sin 2qy ; 4

pffiffiffi 3 ¼ i cos 2qx sin 2qy ; 2 ¼

 1 3g þ 1  1  cos 2qx cos 2qy þ ð1  cos 4qx Þ  X2 ; 4 2

  1 ¼  i sin 2qx 2 cos 2qx þ cos 2qy ; 2  1  gb  X2 3 þ cos4qx þ 2cos 2qx cos2qy þ 3gb  ; 2 P

ð3;3Þ

¼

ð2;1Þ

¼ Sin ;

ð3;1Þ

¼ Sin ;

ð3;2Þ

¼ Sin ;

Sin Sin Sin Sin

ð1;2Þ

ð12-1Þ

ð12-2Þ

ð12-3Þ ð12-4Þ ð12-5Þ

ð12-6Þ ð12-7Þ

ð1;3Þ

ð12-8Þ

ð2;3Þ

ð12-9Þ

where g ¼ nn =ns denotes the ratio of normal to shear rigidity. Solving the eigenvalue problem of Eq. (11), one can obtain the dispersion curves for the in-plane motion.

7

4.2. Dispersion curves and complete band gaps To demonstrate the influence of rotation, Fig. 6 presents the dispersion curves when g ¼ 1 while the value of gb is increasing from 0 to 1.5. For the in-plane motion, the eigenmodes contain three components: two in-plane displacements ux and uy and one rotation U. The natures of modes are labeled in different colors accordingly. Blue curves correspond to the pure displacement modes, red curves to the coupled translational/rotational modes with a predominance of rotation while yellow curves to the coupled translational/rotational modes with a predominance of translation. Based   on the dynamical matrix Sin qx ; qy , all the modes at the high symmetry points C, M and K can be studied analytically. At the C point, the eigenmodes of displacement are fixed at Xu ¼ 0 and the U pffiffiffiffiffiffiffiffiffiffiffi mode appears at XW ¼ 15=2. For the case of M point, the eigenpffiffiffi frequencies of two u modes are degenerated at Xu ¼ 2, while the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U mode at XU ¼ 5ð1 þ 2gb Þ=2 depends on gb , which determines the width of the forbidden gap along the CM direction. For pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gb P 0, the forbidden gap 2 6 X 6 5ð1 þ 2gb Þ=2 grows with pffiffiffi pffiffiffiffiffiffiffiffiffiffiffi increasing bending rigidity. It stabilizes at 2 6 X 6 15=2 when gb P 1. For waves propagating along the CK direction, since the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi eigenfrequency of the U mode at XU ¼ 15ð3gb þ 1Þ=8 depends on the values of gb , changing the value of gb could open a forbidden band gap in this direction. For example, when gb P 1=15, there pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exists a forbidden band gap for 3=2 6 X 6 15ð1 þ 3gb Þ=8. Similarly, for waves propagating in the MK direction, the forbidden pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi band gap 3=2 6 X 6 15ð1 þ 3gb Þ=8 opens when gb P 1=15. Combining the information on the modes at C, M and K points, it can be predicted that a complete band gap appears above X P 3=2 for gb P 1=15. Its width grows up to a maximum value pffiffiffiffiffiffiffiffiffiffiffi 3=2 6 X 6 15=2 when gb P 1. Note that, as shown in Fig. 6, the zero-frequency modes do not exist in this case even when the bending rigidity comes to zero. Although, for compactness, we have analyzed above only the case of equal normal and shear rigidities, the dispersion curves can be obtained for their arbitrary ratio, when necessary. 4.3. Degenerated modes and Dirac cone The Dirac cone can be found at the K point, noted by green arrows in Fig. 7, but the nature of this Dirac cone is different from the case of out-of-plane motion. According to the dynamical matrix

Fig. 6. The dispersion curves of hexagonal membrane with in-plane motion when g ¼ 1 and gb is increasing from 0 to 1.5. Blue lines represent purely displacement modes. Red (yellow) curves correspond to coupled displacement-rotation modes with the pre-dominance of rotation (displacement). Forp0ffiffiffiffiffiffiffiffiffiffiffi 6 gb 6 1=15, there is no complete band gap. For 1=15 6 gb 6 1, the complete forbidden gap appears and its width grows with the increase of gb . It stabilizes at 3=2 6 X 6 15=2, when gb P 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

8

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

Fig. 7. Dispersion curves in hexagonal membrane with in-plane motion. Dirac cones and triple degenerated modes at K points noted by green arrows. In (a) and (b), changing bending/normal rigidities can modify the position of Dirac cone. Triple degenerated modes induced by different combinations of bending and normal rigidities in (c) and (d). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

  Sin qx ; qy at the K point, three modes exist: (1) an ux mode at pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Xux ¼ 9ðg þ 1Þ=8, (2) an uy mode at Xuy ¼ 9ðg þ 1Þ=8, and (3) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a U mode at XW ¼ 15ð3gb þ 1Þ=8. Different from the case of out-of-plane motion, the Dirac cone in this case is the result of the degeneracy of two translational modes. By setting the eigenfrequencies equal to each other, namely by assuming 3g  15gb ¼ 2, the displacement modes degenerate with the rotational mode, exhibiting a triple degenerated mode at K points, as depicted in Fig. 7(c) and (b).

5. Honeycomb monolayer membrane with out-of-plane motion 5.1. Dispersion curves The honeycomb monolayer membrane and its corresponding first Brillouin zone are shown in Fig. 1(b). The calculation of the dispersion curves is similar to the case of the hexagonal membrane with out-of-plane motion, but now the unit cell contains two sublattices A and B. Using the analytical Lagrangian method for each of them, one could finally get the dynamical equation by looking for plane wave solutions:

 * S0out q0x ; q0y v 0out ¼ 0; with the eigenvector

ð13Þ * 0 out

v

  ¼ ðuA ; UA ; WA ; uB ; UB ; WB ÞT . S0out q0x ; q0y is

the dynamical matrix (see the details in Appendix C). The dispersion curves when gt ¼ 0 while gb is changing are shown in Fig. 8(a)–(c), where a complete forbidden band gap can exist for a particular domain of values of gb . For 0 < gb < 0:8, the width of the complete pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi band gap 15ðgb þ gt Þ=2 6 X 6 6 decreases when gb is increased. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi For gb ¼ 0:8, the U modes and W modes ðX ¼ 15ðgb þ gt Þ=2Þ pffiffiffi degenerate with the u modes ðXu ¼ 6Þ, forming the triple degeneracy at C point, and the complete band gap disappears. For 0:8 < gb < 1, a complete band gap reappears, and its width grows when gb is increased. When gb P 1, the width of this gap stabilizes pffiffiffi pffiffiffiffiffiffiffiffiffiffiffi at 6 6 X 6 15=2. In Fig. 8(d)–(f), where gb ¼ 0 while gt is changing, the second band gap disappears when gt P 0:8. Nevertheless, a flat band, which supports the zero-frequency mode having zero group velocity as well, could be noticed when gb ¼ 0.

5.2. Zero-group-velocity and zero-frequency modes We note that two flat bands with zero-group-velocity appear in Fig. 8(a)–(c). As labeled in colors in the dispersion diagrams, the nature of these bands is purely rotational motion. In the absence of torsional rigidity, according to Eq. (13), the eigenfrequencies of the flat bands can be calculated. (1) The first flat band with the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi eigenfrequency at X ¼ 15gb =2. Changing the value of gb would only shift the position of the flat band, while the propagation of this zero-group-velocity mode cannot be initiated. For simplicity, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the motion of C point is analyzed by setting X ¼ 15gb =2 in   the dynamical matrix S0out q0x ; q0y . The eigenvectors are given by: uA ¼ uB ¼ 0; UA ¼ UB and WA ¼ WB . It implies that beads in two sublattices have only rotational motion with the same amplitudes but in opposite directions. As plotted in Fig. 2(c), beads rotating in opposite direction can activate only the bending interactions, without inducing normal interactions, thus the system exhibits non-zero oscillations and the eigenfrequency of the oscillation is determined only by bending rigidity. However we should notice that torsional rigidity, currently neglected, is expected to influence the considered oscillations, because the contacting beads are rotating in opposite directions. (2) The second flat band with pffiffiffiffiffiffiffiffiffiffiffi the eigenfrequency at X ¼ 15=2. We note that the position of this band keeps unchanged in the phonon spectra, when the bending rigidity is changing. By analyzing the eigenvectors for the fixed pffiffiffiffiffiffiffiffiffiffiffi eigenfrequency X ¼ 15=2 at C point, the following relations are expected: uA ¼ uB ¼ 0; UA ¼ UB and WA ¼ WB . This shows that beads have the same rotations in two sublattices. As plotted in Fig. 2(b), beads with the same rotations elongate the shear spring, thus the eigenfrequency relies on shear rigidity. Neither the bending interactions nor the currently neglected torsional interactions are involved because the contacting beads are rotating in the same direction. Analysis for other wavevectors can also be done by setting the flat bands frequencies in the dynamical matrix. One can acquire similar results where for the first flat band pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðX ¼ 15gb =2Þ, contacting beads, i.e., belonging to different sublattices A and B, have rotations with the same amplitude but in opposite direction, resulting in eigenfrequency controlled by bendpffiffiffiffiffiffiffiffiffiffiffi ing rigidity; while for the second flat band ðX ¼ 15=2Þ, all the beads exhibit the same rotations, leading to eigenfrequency depending only on shear rigidity. In both cases, propagation of wave is forbidden. The fact that the constant frequency modes

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

9

Fig. 8. Band structures for the out-of-plane motion of the honeycomb membrane with different torsional gt and bending rigidities gb . Red lines represent purely rotational modes. Blue (yellow) curves correspond to coupled displacement-rotation modes with the pre-dominance of displacement (rotation). In (a), (b) and (c), we set gt ¼ 0 and modify gb . For gb – 0, there is always a complete forbidden gap. When gb ¼ 0, two rotational modes degenerate with the u mode and the band gap disappears. In (d), (e) and (f), the zero-frequency mode exists when gb ¼ 0. The width of the lower complete band gap grows when increasing gt , while the band gap around X ¼ 3=2 disappears when gt > 0:8. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9. Dispersion curves of honeycomb lattices in the case of out-of-plane motion. (a) In the absence of bending and shear rigidity, gb ¼ gt ¼ 0, there are three degenerated zero-frequency modes and one constant frequency mode/oscillation. (b) Only gb – 0 destroys triple degeneracy by inducing two slow waves and transforming zero-frequency (zero-energy) mode into non-propagating oscillation of finite non-zero frequency. (c) Only gt – 0 transforms zero-group-velocity oscillation in slow mode and destroys triple degeneracy by inducing two slow modes and keeping a single zero frequency mode. (d) Both gb – 0 and gt – 0 transform all non-propagating modes represented in (a) into slow modes.

are composed of either collinear or anti-collinear rotations of the two sublattices A and B, with one mode controlled by bending rigidity only and another by shear only, can be proved in the general case of an arbitrary propagating wave vector and not only along the symmetry direction. However, the excitation of these non-propagating modes depends heavily on the torsional interaction. Even for weak torsional interaction, the zero-group-velocity modes could start to be transformed into modes propagating with slow velocity. When gb ¼ 0, zero-frequency modes exist for non-zero torsional rigidity. As shown in Fig. 9(a), for the case when both bending and torsional rigidities are ignored, the optical bands (three highest frequency bands) would transform into two propagative bands and one zero-group-velocity band, while the acoustical

bands (three lowest frequency bands) collapse into the flat bands that only support zero-frequency modes. As expected, the constant frequency mode/oscillation is a purely rotational mode, nonpropagative due to the lack of bending and torsional rigidities. For the case of gb – 0 only, shown in Fig. 9(b), the lower triple degeneracy bands are destroyed by transforming into two slow waves and another zero-group-velocity oscillation. As explained above, since torsional interaction is zero, the constant frequency bands are not transformed into propagating ones by controlling the bending rigidity only. When only gt – 0, Fig. 9(c), the upper zero-group-velocity oscillation is transformed into slow mode and two of the triple degeneracy bands will be transformed into slow modes and only a zero-frequency oscillation left. As depicted in Fig. 9(d), when both gb – 0 and gt – 0, all non-propagating

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

10

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

modes in Fig. 9(a) are transformed into slow modes. Intuitively, the beads in the honeycomb membrane interact with the three neighbor beads only. Compared with the ones in hexagonal membrane where each bead has six nearest neighbors, they undergo less constraints and possess more degrees of freedom. Ignoring one of the rotation-induced interactions (bending/torsion) would further release the bead from constraints. In the absence of bending and shear rigidity, gb ¼ gt ¼ 0, beads reach the maximum freedom, thus the number of non-propagating modes is maximum. Since there are more constraints to the beads in hexagonal membrane, one could expect less flat bands and less zero-frequency modes as shown in Fig. 5(a). The movements of beads in three different positions along CM when gb ¼ 0; gt ¼ 0:1 are depicted in Fig. 10(b)–(d). Similar to the case in Section 3.3, at C point, beads just translate in z-direction with zero displacement relative to their neighbors, leading to contacts kept unstrained. For the waves with nonzero wave vectors along the x-axis, the rotations and translations of the beads at each side of the interparticle contacts are combined in such a way that the resulting displacements of the opposite sides of the contact are equal and the contact is kept unstrained, as it has been discussed already in Section 3.3. In Fig. 10, we illustrate this general situation just by two particular examples. In the middle point of the CM (Fig. 10(c)), one could see non-equal translational motion of some neighboring beads in one direction. This could be compensated (in the sense of keeping the contact unstrained) by the rotation of one of the beads in the absence of the rotations of the other. In the M point (Fig. 10(d)), if the neighboring beads are rotating in the same direction, then they are exhibiting translational motion in the opposite directions and vice versa. In an arbitrary point along the CM direction all the beads are in general exhibiting both translations and rotations, but there are combinations in the amplitudes of these two motions that are keeping the intergrain contacts unstrained, thus creating the configurations supporting the zeroenergy modes. 5.3. Degenerated modes and Dirac cone Similar to the case of hexagonal membrane, we observe the Dirac cones and triple degenerated modes at K point. The difference is that in honeycomb membrane, there are two Dirac cones,

one for the acoustical type modes the other for optical type modes. Principally, two Dirac cones, which are based on rotational modes, lead to two possible ways for the triple degenerated modes: (1) the acoustical Dirac cone degenerates with the translational mode (u mode); (2) the optical one degenerates with another u mode. Fig. 11(a) and (b) show the results of these cases, respectively. Since it is impossible for the two Dirac cones to intersect together, double Dirac cone cannot be found at the K point. By carefully controling gb and gt , one could expect that there could be accidentally triple degenerated modes with linear dispersions at C point. Analytically, to get a triple degeneracy (Dirac-like cone) at the center of Brillouin zone, the only possible solution is that two rotational modes of the eigenfrequency pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X ¼ 15ðgb þ gt Þ=2 degenerate with the fixed u mode at X ¼ 6, giving the condition of gb þ gt ¼ 4=5 (see Fig. 11(c)). More interestingly, when gb and gt perfectly satisfy gb þ gt ¼ 1, a double Dirac cone (the degeneracy of four rotational modes), could be expected at the center of Brillouin zone in Fig. 11(d). Recently, it was reported that the dispersion in the vicinity of accidentally triple degeneracy at C point could lead to the zero-refractive-index in the photon and phonon systems [36]. A phononic crystal with the dispersion exhibiting double Dirac cone gives rise to the Talbot effect and the phenomenon of propagating waves becomes insensitive to defects [37]. Our theoretical results presented above indicate that the granular membranes, under conditions of some particular modifications, could support similar phenomena. 6. Honeycomb monolayer membrane with in-plane motion 6.1. Dispersion curves The structure of the membrane and the corresponding first Brillouin zone is shown in Fig. 1(b). The calculation of the dispersion curves is similar to the case of hexagonal membrane with in-plane motion, but for two coupled sublattices A and B, respectively. The dynamical equation has the following form: *

S0in ðq0 ; p0 Þv 0in ¼ 0;

ð14Þ

*

where v 0in ¼ ðuA ; UA ; WA ; uB ; UB ; WB ÞT is the eigenvector. S0in ðq0 ; p0 Þ is the dynamical matrix (see the details in Appendix C). The dispersion

Fig. 10. The zero-frequency mode along CM in honeycomb membrane with out-of-plane motion and corresponding movements of beads at three different points. (b) C point: kx ¼ 0; ky ¼ 0; (c) The point in the middle of CM: kx ¼ p=3a; ky ¼ 0; (d) M point: kx ¼ 2p=3a; ky ¼ 0.

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

11

Fig. 11. In (a) and (b), triple degenerated modes in honeycomb membrane with out-of-plane motion are induced by different combinations of bending and torsional rigidities. (c) Triple degenerated modes. (d) Double Dirac cone at C point.

Fig. 12. Band structures for in-plane motion of honeycomb membrane with different normal g and bending rigidities gb . (a), (b) and (c) are the cases of absence of bending rigidity. By changing g, the complete forbidden gap disappears when g > 1. In (d), (e) and (f), g is fixed at 1 while gb is increasing. The first complete forbidden band gap appears when gb > 2=15. For gb > 0:4, the second band gap is observed.

curves for gb ¼ 0 are shown in Fig. 12(a)–(c), where a complete forbidden band gap could exist for a particular range of values of g, the ratio of the normal to shear rigidity of the contact. For g < 1, the pffiffiffi width of the complete band gap X 6 3 shrinks when increasing g. For g P 1, the complete band gap disappears. We note that the zero-frequency mode exists only in CM direction. Changing the value of g cannot initiate the propagation of this mode. For the case when g ¼ 1 and gb – 0 illustrated in Fig. 12(d)–(f), the first band gap appears when gb P 2=15. For gb P 0:4, the second band gap pffiffiffi pffiffiffiffiffiffi is opened, and its width stabilizes at 6 6 X 6 15 when gb P 1. 6.2. Zero-frequency modes As shown in Fig. 12(a)–(c), the zero-frequency mode for inplane motion in honeycomb membrane exists only in the absence of bending rigidity. Different from the case of out-of-plane motion where zero-frequency modes could be expected in all directions, here a zero-frequency mode exists only in a single particular direction (CM). As we have described for the out-of-plane cases, motion of rotational type plays a major role for the formation of nonpropagating modes. For in-plane motion, only one rotational degree of freedom, i.e., rotation along z-axis, exists for the beads, thus only bending interaction can be activated. No bending rigidity

can still lead to zero-frequency mode, but only in CM because for other directions the motion in the system is dominated by translations so that rotations are too weak to compensate them and keep all the contacts unstrained. In addition, compared with the dispersions in hexagonal system for the in-plane motion (see Fig. 6), a zero-frequency mode could be found only in the honeycomb lattice. This is expected because the beads in the hexagonal membrane have more constraints than those in the honeycomb lattice. 6.3. Degenerated modes and Dirac cone As displayed in Fig. 13(a) and (b), one also could observe Dirac cones and triple degenerated modes at K point. The natures of the Dirac cones are translational modes, and the triple degeneracy is the consequence of the intersection of the Dirac cone degenerating with the third mode (U mode). By carefully choosing the combination of gb and g, accidentally triple degenerated modes with linear dispersions could be found at C point. Analytically, this triple degeneracy (Dirac-like cone) is formed by the degeneracy of the pffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi rotational mode ðXU ¼ 15 or XU ¼ 15gb Þ with two u modes pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðXu ¼ 3g þ 3Þ, giving the condition of g ¼ 4 or gb ¼ ðg þ 1Þ=5. When g ¼ 4 and gb ¼ 1; a double Dirac cone, i.e., the degeneracy of two Dirac cones, is observed at the center of Brillouin zone.

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

12

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

Fig. 13. Band diagrams for honeycomb membrane with in-plane motion, (a) and (b) Triple degenerated modes induced by different combinations of bending and torsional rigidities. (c) Triple degenerated mode. (d) Double Dirac cone at C points.

7. Conclusion

Appendix A. Length variations of the springs

To conclude, the analytical and numerical evaluations of both the out-of-plane and in-plane motions for the hexagonal/honeycomb monolayer membranes have been presented. Above a critical value of the ratio of normal/bending/torsional rigidities to shear rigidity of the interparticle contacts, a complete forbidden band gap for the elastic wave propagation exists. By introducing the rotational degrees of freedom, one could manipulate the dispersion of the phononic crystal in a flexible way. In the absence of bending/torsional rigidities, the non-propagating modes are observed in the honeycomb membranes or the hexagonal membranes for the case of out-of-plane motion. Rotational degrees of freedom are crucial for the existence of non-propagating modes. When torsional rigidity is ignored, the zero-group-velocity mode, whose nature is purely rotational, exists in the membrane by activating only the shearing or bending interactions. For the zero-frequency mode, at Brillouin center (zero wave vector), the relative translations of adjacent beads are zero. Thus, the contacts remain unloaded and the system costs zero energy. For non-zero wave vectors, unequal translations, in particular those of the same amplitude but opposite direction, could be compensated by the rotations of contacting beads in the same direction. This keeps the contacts between the particles unstrained and the energy of the system unmodified. We revealed the conditions in the considered membranes non-propagating modes at finite (non-zero) frequencies could also exist. These non-propagating modes have been identified as pure rotational oscillations. We demonstrated that even weak bending/torsional rigidities of the contacts can, under some conditions, transform the non-propagating zerofrequency modes and rotational oscillation into propagating modes of membrane. Thus the study of the non-propagating modes is useful for the design of the granular phononic crystals supporting the controllable propagation of the slowest waves. Based on the analytical derivations, we also discussed the formation and manipulation of the degenerated modes and Dirac cones. This could provide variety of possibilities for potential applications of granular membranes in the control of elastic waves.

A.1. Out-of-plane motion

Acknowledgments This work was supported by the fellowship of China Scholarship Council for L.-Y. Zheng. G.T. acknowledges financial support from FP7-CIG (Project 618322 ComGranSol). We also acknowledge fruitful discussions with Professor Ming-Hui Lu from College of Engineering and Applied Sciences, Nanjing University, Nanjing, Jiangsu, China.

In this work, Dsi ; Dbi and Dt i respectively correspond to the effective elongation of the effective springs accounting for the shear, bending and torsional contact interactions between the *

neighbor beads. We denote by e 0i the unit vector in the direction *

*

from the zero-th bead center to the i-th bead center. e x , e y and * ez *

represent the unit vector along x-, y- and z-axis, respectively. *

*

t 0i , the unit vector normal to e 0i and e z , has the form:

*

*

* e 0i .

t 0i ¼ e z  The effective elongations of the three types of springs are given by

* * * * a a Dsi ¼ ðui  u0 Þ  ðui þ u0 Þe x  t 0i  ð/i þ /0 Þe y  t 0i ; 2 2 * * * * a a Dri ¼ ðui  u0 Þe x  t 0i þ ð/i  /0 Þe y  t 0i ; 2 2 * * * * a a Dli ¼ ðui  u0 Þe x  e 0i þ ð/i  /0 Þe y  e 0i : 2 2

A.2. In-plane motion In the analysis of the in-plane motion, Dni is the elongation in length of the spring between the zeroth and the i-th bead that activates the normal forces. Dsi ; Dri still correspond to effective elongations of shear spring and bending spring, respectively. Considering the in-plane motion, they have the following forms *

*

*

*

Dni ¼ ðuxi  ux0 Þe x  e 0i þ ðuyi  uy0 Þe y  e 0i ; * * * * a Dsi ¼ ðuxi  ux0 Þe x  t 0i þ ðuyi  uy0 Þe y  t 0i  ðui þ u0 Þ; 2

a Dri ¼ ðui  u0 Þ: 2 Appendix B. Equations of motion for hexagonal membranes B.1. Out-of-plane motion According to Eqs. (1)–(4) in Section 3.1, one can obtain the equations of motion:

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

13

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

€0 ¼ Mu

6 X

ns Dsi ;

X

€1 ¼ Mu

i¼1 6 h * *  * *  * * i X €0 ¼ R Iu ns Dsi e x  t 0i þ nb Dbi e x  t 0i þ nt Dti e x  e 0i ;

€1 ¼ R Iu

6 h X

* *  * *  * * i ns Dsi e y  t 0i þ nb Dbi e y  t 0i þ nt Dt i e y  e 0i :

i¼1

Thus we can get the dynamical equation Eq. (5) in Section 3.1. B.2. In-plane motion

Xh

* *  * *  * * i ns Dsi e x  t 1i þ nb Dbi e x  t 1i þ nt Dti e x  e 1i ;

i¼0;2;4

i¼1

€0 ¼ R I/

ns Dsi ;

i¼0;2;4

€1 ¼ R I/

Xh

* *  * *  * * i ns Dsi e y  t 1i þ nb Dbi e y  t 1i þ nt Dt i e y  e 1i :

i¼0;2;4

Subtituting the expressions of elongation in Appendix A, the dynamical equation is obtained:

 * S0out q0x ; q0y v 0out ¼ 0; with the dynamical matrix in the form,

According to Eqs. (1), (4), (9) and (10), the equations of motion are given by: 6 h X

€ x0 ¼ Mu

* *  * * i nn Dni e x  e 0i þ ns Dsi e x  t 0i ;

i¼1 6 h X

€ y0 ¼ Mu

* *  * * i nn Dni e y  e 0i þ ns Dsi e y  t 0i ;

  Q S0out q0x ; q0y ¼ N

U



Q

;

where Q ; U; N are 3  3 matrices of the following forms,

0

B Q ¼B @

X2  3

0 X2

0

P

 3ð1þg2b þgt Þ

0

i¼1

0

0

0 X2 P

 3ð1þg2b þgt Þ

1

C C; A

  1 pffiffiffi iq0 0 0 0 0  3ie x sin q0y  e2iqx  eiqx cos q0y e2iqx þ 2eiqx cos q0y C B h i pffiffi pffiffiffi iq0 C B 0 33gb gt iq0x 3ð1gb þgt Þ iqx 0 0 2iq0x 0 x C; U¼B sin q  e cos q  g e ie sin q 3 ie t y y y 2 2 C B @  h iA pffiffi 0 0 0 0 0 3ð1gb þgt Þ iqx 1gb 3gt iqx 2iqx iqx 0 0 0 2iqx e  e cos qy ie sin qy  e cos qy þ ð1  gb Þe 2 2 0

 0  1 pffiffiffi iq0 0 0 0 e2iqx þ 2eiqx cos q0y  3ie x sin q0y e2iqx  eiqx cos q0y C B h i pffiffi pffiffiffi iq0 C B 0 33gb gt iq0x 3ð1gb þgt Þ iqx 0 0 2iq0x 0 x C: N¼B sin q  e cos q  g e ie sin q 3 ie  t y y y 2 2 C B @   h iA pffiffi 0 0 0 0 0 3ð1gb þgt Þ iqx 1gb 3gt iqx 2iqx iqx 0 0 0 2iqx   e e cos qy ie sin qy  e cos qy þ ð1  gb Þe 2 2 0

6 X € 0 ¼ R ½ns Dsi þ nb Dbi : Iu i¼1

Thus, we can get the dynamical equation Eq. (11) in Section 4.1.

Above, P ¼ MR2 =I. q0x ¼ kx a=2; q0y ¼

pffiffiffi 3ky a=2 are the normalized

wave vectors for kx ; ky . X ¼ x=x0 is the normalized frequency with pffiffiffiffiffiffiffiffiffiffiffi x0 ¼ ns =M. gb ¼ nb =ns and gt ¼ nt =ns denote the ratios of bending and torsional rigidities to shear rigidity, respectively.

Appendix C. Dynamical matrix for honeycomb membranes

C.2. In-plane motion

C.1. Out-of-plane motion

According to Eqs. (1), (4), (9) and (10), the equations of motion are given by: For sublattice A (beads marked by odd numbers):

According to Eqs. (1)–(4) in Section 3.1, one can write the equations of motion as follows: For sublattice A (beads marked by odd numbers):

X

€0 ¼ Mu

Xh

* *  * *  * * i ns Dsi e x  t 0i þ nb Dbi e x  t 0i þ nt Dt i e x  e 0i ;

i¼1;3;5

€0 ¼ R I/

Xh

* *  * * i nn Dni e x  e 0i þ ns Dsi e x  t 0i ;

i¼1;3;5

ns Dsi ;

i¼1;3;5

€0 ¼ R Iu

€ x0 ¼ Mu

Xh

€ y0 ¼ Mu

Xh

* *  * * i nn Dni e y  e 0i þ ns Dsi e y  t 0i ;

i¼1;3;5

€0 ¼ R Iu

X

½ns Dsi þ nb Dbi :

i¼1;3;5

* *  * *  * * i ns Dsi e y  t 0i þ nb Dbi e y  t 0i þ nt Dt i e y  e 0i :

i¼1;3;5

For sublattice B (beads marked by even numbers):

For sublattice B (beads marked by even numbers):

€ x1 ¼ Mu

Xh

* *  * * i nn Dni e x  e 1i þ ns Dsi e x  t 1i ;

i¼0;2;4

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005

14

L.-Y. Zheng et al. / Ultrasonics xxx (2015) xxx–xxx

€ y1 ¼ Mu

Xh

* *  * * i nn Dni e y  e 1i þ ns Dsi e y  t 1i ;

i¼0;2;4

€1 ¼ R Iu

X

½ns Dsi þ nb Dbi :

i¼0;2;4

Subtituting the expressions of elongation in Appendix A, the dynamical equation is obtained:

 * S0in q0x ; q0y v 0in ¼ 0;

with the dynamical matrix in the form,

S0in

  q0x ; q0y ¼

Q0 N0

! U0 ; Q0

where Q 0 ; U 0 ; N0 are 3  3 matrices of the following forms,

0 3ð1þgÞ

B Q0 ¼ B @

2

 X2

0 0

0 3ð1þgÞ 2

 X2

0

1

0

C C; A

0 2

3ð1 þ gb Þ  XP

[12] A. Souslov, A.J. Liu, T.C. Lubensky, Phys. Rev. Lett. 103 (2009) 205503. [13] X. Mao, N. Xu, T.C. Lubensky, Phys. Rev. Lett. 104 (2010) 085504. [14] K. Sun, A. Souslov, X. Mao, T.C. Lubensky, Proc. Natl. Acad. Sci. 109 (2012) 12369–12374. [15] C.L. Kane, T.C. Lubensky, Nature Phys. 10 (2014) 39–45. [16] J. Paulose, B.G. Chen, V. Vitelli, Nature Phys. 11 (2015) 153–156. [17] P. Wang, L. Lu, K. Bertoldi, Phys. Rev. Lett. 115 (2015) 104302. [18] A. Merkel, V. Tournat, V. Gusev, Phys. Rev. E 82 (2010) 031305. [19] V. Tournat, I. Perez-Arjona, A. Merkel, V. Sanchez-Morcillo, V. Gusev, New J. Phys. 13 (2011) 073042. [20] H. Pichard, A. Duclos, J.-P. Groby, V. Tournat, V. Gusev, Phys. Rev. B 86 (2012) 134307. [21] A. Merkel, V. Tournat, V. Gusev, Phys. Rev. Lett. 107 (2011) 225502. [22] J. Cabaret, P. Béquin, G. Theocharis, V. Andreev, V.E. Gusev, V. Tournat, Phys. Rev. Lett. 115 (2015) 054301. [23] A.V. Akimov, Y. Tanaka, A.B. Pevtsov, S.F. Kaplan, V.G. Golubev, S. Tamura, D.R. Yakovlev, M. Bayer, Phys. Rev. Lett. 101 (2008) 033902. [24] I. Lisiecki, V. Halté, C. Petit, M.-P. Pileni, J.-Y. Bigot, Adv. Mater. 20 (2008) 4176– 4179. [25] A. Ayouch, X. Dieudonné, G. Vaudel, H. Piombini, K. Vallé, V. Gusev, P. Belleville, P. Ruello, ACS Nano 6 (12) (2012) 10614–10621. [26] L.M. Schwartz, D.L. Johnson, S. Feng, Phys. Rev. Lett. 52 (1984) 831. [27] J.C. Maxwell, Phil. Mag. 27 (1865) 294. [28] C.R. Calladine, Int. J. Solids Struct. 14 (1978) 161–172. [29] C.S. Chang, L. Ma, Int. J. Solids Struct. 28 (1991) 67–86. [30] A.C. Eringen, Microcontinuum Field Theories. 5: Theory of Micropolar Elasticity, Springer, New York, 1999. [31] K.I. Kanatani, Int. J. Eng. Sci. 17 (1979) 419–432.

 pffiffi 1 0  pffiffiffi iq0 0 0 iq0  3þ2 g eiqx cos q0y þ ge2iqx  3ie x sin q0y  3ð2g1Þ ie x sin q0y C B     pffiffi C B 0 3ðg1Þ iqx 1þ3g iq0x 0 0 2iq0x 2iq0x iq0x 0 C; U0 ¼ B  ie sin q  e cos q þ e  e cos q  e y y y 2 2 C B @    A pffiffiffi iq0 0 0 0 0 0 2iqx iqx 0 iqx 0 2iqx x e  e cos qy ðgb  1Þ 2e cos qy þ e 3ie sin qy  pffiffi 1 0  pffiffiffi iq0 0 0 0 3ðg1Þ iqx ie sin q0y  3ie x sin q0y  3þ2 g eiqx cos q0y þ ge2iqx 2 C B    0  pffiffi C B 0 3ðg1Þ iqx 1þ3g iq0x 0 0 2iq0x 2iqx iq0x 0 C: N0 ¼ B ie sin q  e cos q þ e  e cos q e y y y 2 2 C B @  0   A pffiffiffi iq0 0 0 0 0 2iqx iqx 0 iqx 0 2iqx x sin qy  e e cos qy cos qy þ e 3ie ðgb  1Þ 2e

Above, P ¼ MR2 =I. q0x ¼ kx a=2; q0y ¼

pffiffiffi 3ky a=2 are the normalized

wave vectors for kx ; ky . X ¼ x=x0 is the normalized frequency with pffiffiffiffiffiffiffiffiffiffiffi x0 ¼ ns =M. gb ¼ nb =ns and g ¼ nn =ns denote the ratios of bending and normal rigidities to shear rigidity, respectively. References [1] R. Martinez-Salar, J. Sancho, J.V. Sanchez, V. Gomez, J. Llinares, F. Meseguer, Nature 378 (1995) 241. [2] M.H. Lu, L. Feng, Y.F. Chen, Mater. Today 12 (2009) 34. [3] A.S. Phani, J. Woodhouse, N.A. Fleck, J. Acoust. Soc. Am. 119 (4) (2006) 1995– 2005. [4] M.H. Lu, C. Zhang, L. Feng, J. Zhao, Y.F. Chen, Y.W. Mao, J. Zi, Y.Y. Zhu, S.N. Zhu, N.B. Ming, Nature Mater. 6 (2007) 744–748. [5] A. Sukhovich, B. Merheb, K. Muralidharan, J.O. Vasseur, Y. Pennec, P.A. Deymier, J.H. Page, Phys. Rev. Lett. 102 (2009) 154301. [6] H. Pichard, A. Duclos, J.-P. Groby, V. Tournat, V.E. Gusev, Phys. Rev. E 89 (2014) 013201. [7] G. Carta, M. Brun, A.B. Movchan, N.V. Movchan, I.S. Jones, Int. J. Solids Struct. 51 (2014) 2213–2225. [8] Z. Liu, X. Zhang, Y. Mao, Y.Y. Zhu, Z. Yang, C.T. Chan, P. Sheng, Science 289 (2000) 1734. [9] A. Spadoni, M. Ruzzene, S. Gonella, F. Scarpa, Wave Motion 46 (2009) 435–450. [10] D. Torrent, D. Mayou, J. Sánchez-Dehesa, Phys. Rev. B 87 (2013) 115143. [11] M. Brun, I.S. Jones, A.B. Movchan, Proc. R. Soc. A 468 (2012) 3027–3046.

[32] O. Mouraille, W.A. Mulder, S. Luding, J. Stat. Mech. (2006) P07023. [33] A.S.J. Suiker, A.V. Metrikine, R. de Borst, Int. J. Solids Struct. 38 (2001) 1563– 1583. [34] H. Pichard, A. Duclos, J.-P. Groby, V. Tournat, L.Y. Zheng, V. Gusev, 2015. Available from: <1510.05685>. [35] T. Kariyado, Y. Hatsugai, 2015. Available from: <1505.06679v1>. [36] X. Huang, Y. Lai, Z.H. Hang, H. Zheng, C.T. Chan, Nature Mater. 10 (2011) 582. [37] Z.G. Chen, X. Ni, Y. Wu, C. He, X.C. Sun, L.Y. Zheng, M.H. Lu, Y.F. Chen, Sci. Rep. 4 (2014) 4613. [38] J. Mei, Y. Wu, C.T. Chan, Z.Q. Zhang, Phys. Rev. B 86 (2012) 035141. [39] X. Zhang, Phys. Rev. Lett. 100 (2008) 113903. [40] X. Ni, C. He, X.C. Sun, X.P. Liu, M.H. Lu, L. Feng, Y.F. Chen, New J. Phys. 17 (2015) 053016. [41] Y.T. Wang, P.G. Luan, S. Zhang, 2014. Available from: <1411.2806>. [42] Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong, B. Zhang, Phys. Rev. Lett. 114 (2015) 114301. [43] R. Susstrunk, S.D. Huber, 2015. Available from: <1503.06808v1>. [44] L.-H. Wu, X. Hu, Phys. Rev. Lett. 114 (2015) 223901. [45] N. Boechler, J.K. Eliason, A. Kumar, A.A. Maznev, K.A. Nelson, N. Fang, Phys. Rev. Lett. 111 (2013) 036103. [46] A. Khanolkar, S. Wallen, M. Abi Ghanem, J. Jenks, N. Vogel, N. Boechler, Appl. Phys. Lett. 107 (2015) 071903. [47] A.A. Maznev, V.E. Gusev, Phys. Rev. B 00 (2015) 005400. [48] S.P. Wallen, A.A. Maznev, N. Boechler, 2015. Available from: <1506.08445>. [49] E. Kogan, V.U. Nazarov, Phys. Rev. B 85 (2012) 115418. [50] H. Tang, B.-S. Wang, Z.-B. Su, Symmetry and lattice dynamics, graphene simulation, in: Jian Gong (Ed.), InTech, 2011, ISBN: 978-953-307-556-3.

Please cite this article in press as: L.-Y. Zheng et al., Zero-frequency and slow elastic modes in phononic monolayer granular membranes, Ultrasonics (2015), http://dx.doi.org/10.1016/j.ultras.2015.11.005