Keywords: cristal granulaire, ondes élastiques, modes de fréquence nulle granular crystal, wave dynamics, zero-frequency modes, topological edge waves, slow edge waves, Dirac point

ondes de bord topologiques, ondes lentes de bord, point de Dirac.

General introduction

The discovery of phononic crystals and metamaterials is one of the most exciting achievements in modern scientific research [1][2][3][4][5][6][7][START_REF] Khelif | Experimental study of guiding and filtering of acoustic waves in a two dimensional ultrasonic crystal[END_REF][START_REF] Laude | Surface acoustic wave trapping in a periodic array of mechanical resonators[END_REF][START_REF] Laude | Dispersion and polarization of surface waves trapped in high aspect ratio electrode arrays[END_REF][START_REF] Achaoui | Experimental observation of locally-resonant and Bragg band gaps for surface guided waves in a phononic crystal of pillars[END_REF][START_REF] Zigoneanu | Three-dimensional broadband omnidirectional acoustic ground cloak[END_REF][START_REF] Nicoletti | Acoustic diodes[END_REF]. The significance of those artificial materials is that they show exotic abilities in wave control that greatly compensate those missing by the conventional materials. Due to their extraordinary abilities in the control of acoustic/elastic wave propagation, many exotic phenomena that are impossible for the nature materials now can be achieved, such as negative reflection [6,7], energy band gap [START_REF] Khelif | Experimental study of guiding and filtering of acoustic waves in a two dimensional ultrasonic crystal[END_REF][START_REF] Laude | Surface acoustic wave trapping in a periodic array of mechanical resonators[END_REF], acoustic filter [START_REF] Laude | Dispersion and polarization of surface waves trapped in high aspect ratio electrode arrays[END_REF][START_REF] Achaoui | Experimental observation of locally-resonant and Bragg band gaps for surface guided waves in a phononic crystal of pillars[END_REF], acoustic cloaking [START_REF] Zigoneanu | Three-dimensional broadband omnidirectional acoustic ground cloak[END_REF], acoustic diode [START_REF] Nicoletti | Acoustic diodes[END_REF], etc.

In the last 10 years, the study of phononic crystals and metamaterials has been extended to discrete granular systems [START_REF] Job | Solitary wave trains in granular chains: Experiments, theory and simulations[END_REF][START_REF] Theocharis | Intrinsic energy localization through discrete gap breathers in one-dimensional diatomic granular crystals[END_REF][START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF][START_REF] Chong | Nonlinear coherent structures in granular crystals[END_REF][START_REF] Texier | Helical locomotion in granular media[END_REF][START_REF] Wallen | Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal[END_REF][START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF]. The concept of granular phononinc crystals or granular crystals has been widely introduced in those systems. In general, granular crystals are spatially periodic structures of particles, most often spherical homogeneous elastic beads, arranged in crystal lattices. In granular crystals, the beads are linked by interconnections, which are of much smaller dimensions and weights than the beads themselves. The interactions between beads take place predominantly via interparticles contact stiffness, e.g. normal and transverse rigidities of these elastic interconnections. These unique features lead to the propagation of elastic waves in granular structures at significantly slower velocities than in the individual grains [START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF].

In addition, due to the existence of non-central shear forces between the particles, rotations of individual grains can be initiated, resulting in the propagation of purely rotational and/or coupled rotational-translational waves in granular crystals [START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF][START_REF] Chong | Nonlinear coherent structures in granular crystals[END_REF][START_REF] Texier | Helical locomotion in granular media[END_REF][START_REF] Wallen | Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal[END_REF][START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF]. Those fascinating properties make granular crystals excellent media to study fundamental wave physics phenomena.

In parallel to the scientific activities in phononic crystals and metamaterials, in recent years the study of the topology in materials has also attracted a lot of attention [START_REF] Moore | The birth of topological insulators[END_REF][START_REF] Bliokh | Quantum spin Hall effect of light[END_REF][START_REF] Liang | Optical resonator analog of a two-dimensional topological insulator[END_REF][START_REF] Peng | Experimental demonstration of anomalous Floquet topological insulator for sound[END_REF][START_REF] Wang | Coriolis force induced topological order for classical mechanical vibrations[END_REF][START_REF] Huber | Topological mechanics[END_REF]. The topologically protected wave propagation is demonstrated in many different research fields including condensed matter physics [START_REF] Moore | The birth of topological insulators[END_REF], optics [START_REF] Bliokh | Quantum spin Hall effect of light[END_REF][START_REF] Liang | Optical resonator analog of a two-dimensional topological insulator[END_REF], acoustics [START_REF] Peng | Experimental demonstration of anomalous Floquet topological insulator for sound[END_REF] and elasticity [START_REF] Wang | Coriolis force induced topological order for classical mechanical vibrations[END_REF][START_REF] Huber | Topological mechanics[END_REF]. In the topological non-trivial materials/structures, the bulk of the materials is insulating for wave propagation, while the edge is metallic, supporting the transport of unidirectional edge waves.

The edge waves are topologically protected in a sense that they are immune to backscattering in the presence of imperfections and disorders. Interestingly, topological wave phenomena could be potentially found in granular crystals. As a peculiar type of phononic crystals for elastic wave, granular crystals not only exhibit those common properties such as allowed bands and forbidden gaps, but also show special advantages in other aspects such as the study of rotational waves, slow waves, etc. It is believed that the successful discovery of topological wave behaviour in granular crystals may pave the way to their practical applications in elastic wave devices and vibration isolation and confinement.

In this manuscript, we theoretically study the wave dynamics in two-dimensional granular crystals or granular membranes. The phononic properties in granular crystals are investigated [START_REF] Zheng | Zero-frequency and slow elastic modes in phononic monolayer granular membranes[END_REF], including the existence of Dirac points, zero-frequency modes, zero-group-velocity modes and their transformation into slow propagating phononic modes due to weak bending and torsional intergrain interactions. Furthermore, the wave dynamics on the edges/boundaries is also analyzed [START_REF] Zheng | Zero-frequency and extremely slow elastic edge waves in mechanical granular graphene[END_REF]. In the presence of edges/boundaries, zero-frequency and extremely slow elastic edge waves can be predicted in mechanical granular honeycomb crystals (granular graphene) when the torsional or/and bending rigidities are weak/vanishing. In addition, we found that by a proper design of the structure, the critical topological phase transition can take place in granular crystals. Topological properties of rotational edge waves in a granular graphene are theoretically demonstrated [START_REF] Zheng | Quasi-topological rotational waves in mechanical granular graphene[END_REF].

The organization of the manuscript is as follows: a brief review on the study of phononic crystals and metamaterials, granular crystals as well as the topological wave propagation is presented in the first Chapter. To begin with, we provide some fundamental concepts that are helpful for the understanding of wave dynamics in granular crystals. This includes the Brillouin zone, reciprocal space, lattice, band structure, contact mechanics and rotational degree of freedom. Some basic notions for the topological wave transport are also presented in the last two sections of Chapter 1.

In the second Chapter, we theoretically study the dispersion properties of elastic waves in hexagonal and honeycomb monolayer granular membranes with either out-of-plane or inplane particle motion. 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 degenerate modes. This could motivate variety of potential applications in elastic waves control by manipulating the contact rigidities in granular phononic crystals.

In the third Chapter, we report the theoretical description and analysis of edge elastic waves in a semi-infinite mechanical granular graphene structure. The dispersion curves of the edge waves are theoretically derived and analyzed for two configurations of boundary. Quasi-flat edge mode dispersion curves with near zero frequency are observed. These quasi-flat dispersion curves, supporting the propagation of waves with extremely slow group velocity, show a tendency to be perfect zero-frequency modes for zero torsional rigidity or vanish for zero bending rigidity, indicating that weak bending and torsional interbead interactions are critical in the transformation of zero-frequency modes into extremely slow propagating modes.

In the last Chapter, we theoretically demonstrate the topological properties of these mechanical rotational waves in a granular graphene, a two-dimensional monolayer honeycomb granular crystal with Dirac dispersion at the center of Brillouin zone. Around the Dirac point, effective spin, helicity and effective spin-orbit coupling are illustrated in the mechanical granular system. Finally, quasi-topological transport, where the rotational edge waves are nearly topological protected, is observed on the interface.

Chapter 1

Wave propagation in periodic structures and metamaterials

Introduction

In the last decades, the study of phononic crystals and metamaterials has become one of the most active research topics. Phononic crystals are usually constructed by scatterers arranged in periodic structures. The study of these periodic structures are actually motivated by the accumulated knowledge and advances in solid state physics. In the long-term study of solid state physics, researchers found that periodic structures can strongly influence/modulate the propagation of electrons. In 1987, these modulation effects were first extended to electromagnetic waves in periodic structures, which are called photonic crystals in analogy with the periodic structures of atoms (crystals) for propagation of electrons [START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF][START_REF] John | Strong localization of photons in certain disordered dielectric superlattices[END_REF]. The period of the spatial distribution of the scatterers in photonic crystals is comparable to the wavelength of the electromagnetic waves.

In parallel to the development of photonic crystals, the idea to control the propagation of acoustic waves by phononic crystals was also proposed [START_REF] Sigalas | Elastic and acoustic wave band structure[END_REF][START_REF] Kushwaha | Acoustic band structure of periodic elastic composites[END_REF]. Similar to the electromagnetic waves, acoustic waves are classical waves closely related to our daily life. Because of their important role in voice devices, communication, underwater detection, ultrasonic testing, seismic exploration, noise control, etc, effective manipulation/control of the propagation of acoustic waves has very strategic and practical significances. With the development of research in photonic/phononic crystals, the concept of metamaterials has also been proposed. Despite the fact that metamaterials are often designed as periodic structures, usually their properties do not depend on periodicity, but on local resonances of the individual scatterers. Metamaterials can be more efficient in influencing on wave propagation than the periodic structures.

Owing to their periodicity and/or local resonances, many exotic effects have been achieved with phononic crystals and metamaterials. For example, frequency band gaps [5][6][7][START_REF] Khelif | Experimental study of guiding and filtering of acoustic waves in a two dimensional ultrasonic crystal[END_REF][START_REF] Phani | Wave propagation in two-dimensional periodic lattices[END_REF], negative refraction [6,7,[START_REF] Lu | Negative birefraction of acoustic waves in a sonic crystal[END_REF], gradient-index lens [5], superlens [6], subwavelength imaging [START_REF] Sukhovich | Experimental and theoretical evidence for subwavelength imaging in phononic crystals[END_REF], cloaking [START_REF] Zhang | Broadband acoustic cloak for ultrasound waves[END_REF][START_REF] Stenger | Experiments on elastic cloaking in thin plates[END_REF], zero-reflective-index [START_REF] Zheng | Acoustic cloaking by a near-zero-index phononic crystal[END_REF][START_REF] Huang | Dirac cones induced by accidental degeneracy in photonic crystals and zero-refractive-index materials[END_REF], etc. The discovery of these artificial materials greatly complements the shortcomings of conventional materials.

In the latest few years, there exists another significant activity in evaluation of acoustic properties of granular phononic crystals, such as 1D chains [START_REF] Inserra | A method of controlling wave propagation in initially spatially periodic media[END_REF][START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF][START_REF] Sanchez-Morcillo | Second-harmonic generation for dispersive elastic waves in a discrete granular chain[END_REF], 2D membranes and packings [START_REF] Job | Solitary wave trains in granular chains: Experiments, theory and simulations[END_REF][START_REF] Theocharis | Intrinsic energy localization through discrete gap breathers in one-dimensional diatomic granular crystals[END_REF][START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF], and 3D assemblages of grains [START_REF] Chong | Nonlinear coherent structures in granular crystals[END_REF][START_REF] Texier | Helical locomotion in granular media[END_REF][START_REF] Wallen | Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal[END_REF][START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF]. Different from most of the phononic crystals, these granular materials have extraordinary ability to support the propagation of rotational waves in addition to common shear and longitudinal waves. These structures host extra acoustic modes, which do not exist in the absence of the rotational degrees of freedom. The rotational modes and their couplings to other modes, in particular to shear ones, could provide more flexibility and additional functionalities in control of the acoustic wave propagation. Specifically, the metamaterials based on granular phononic structures are expected to be advantageous in the monitoring of bulk shear and surface Rayleigh acoustic waves. This could be important for various applications, such as anti-seismic protection, attenuation and cloaking of ground noise, and signal processing with surface acoustic waves, for example. The phononic crystal for hypersound is constructed of air holes in an epoxy matrix with the lattice periodicity a = 1.36 µm and the one for heat is in a silicon matrix of air holes with a = 34 nm. [START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF] Recently, inspired by quantum systems, the study of topological phenomena in materials and metamaterials has also attracted increasing interest for scientific investigation. The developments include solid state physics [START_REF] Hasan | Colloquium: Topological insulators[END_REF][START_REF] Haldane | Model for a quantum Hall effect without landau levels: condensedmatter realization of the parity anomaly[END_REF][START_REF] Qi | Topological insulators and superconductors[END_REF], optics [START_REF] Haldane | Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry[END_REF][START_REF] Rechtsman | Photonic Floquet topological insulators[END_REF], acoustics [START_REF] Ni | Topologically protected one-way edge mode in networks of acoustic resonators with circulating air flow[END_REF][START_REF] Peano | Topological phases of sound and light[END_REF], mechanical systems [START_REF] Peano | Topological phase transitions and chiral inelastic transport induced by the squeezing of light[END_REF][START_REF] Chen | Topological mechanics of origami and kirigami[END_REF][START_REF] Paulose | Topological modes bound to dislocations in mechanical metamaterials[END_REF], etc. The topological properties of the band structure of waves can be exploited to achieve unique and exciting functionalities, such as the existence of edge waves that are immune to backscattering in the presence of imperfections and impurities, e.g. localized defects and corners. This concept can pave the way to structural configurations which may find practical application in acoustic waveguiding or vibrations isolation and confinement.

In this chapter, we firstly introduce some basic concepts of phononic crystals and metamaterials. In addition, fundamental theories and the developments in elasticity/acoustics of granular phononic crystals are demonstrated. Finally, basic topological phenomena in wave propagation are discussed.

Phononic crystals and metamaterials

In this section, we introduce the basic characteristics, classification and research progress in phononic crystals and metamaterials.

Phononic crystals

Phononic crystals are periodic composites or structures where the acoustic waves are exhibiting several dispersion/energy bands resulting from efficient multiple wave scattering. In particular, different bands can be separated by the energy/frequency gaps, where the propagation of the waves is forbidden. Fig. 1. 1(a) shows the schematics of phononic crystals of different dimensionality [START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF]. The concept of phononic crystal was proposed the first time quasi-simultanously by M. S. Kushwaha, et al. and Sigalas et al. in 1993 [36, 58]. They showed the existence of phonon band gap in a two-dimensional periodic elastic composite. In 1995, R. Martinez-Sala et al, carried out the first experimental demonstration of phononic crystal properties (left panel of Fig. 1.1(c)) on a minimalist sculpture in Madrid, proving the existence of phonon band gap [START_REF] Martinez-Salar | Sound attenuation by sculpture[END_REF].

Since then, the study of phononic crystal began to attract worldwide attention. Fig. 1.1(b) shows a typical band structure of a phononic crystal. Lines label the propagating branches of phonon modes, while the frequency ranges where the modes are absent correspond to phonon band gaps.

The orange area in Fig. 1.1(b) marks a complete band gap. In general, wave propagation in the band gap is suppressed, but in the pass band waves can be transmitted without attenuation.

It is the period of the structure which, through the condition of the Bragg scattering, controls the wavelength of the acoustic wave that is most influenced by the structure. Fig. 1.1(c) shows three phononic crystals corresponding to different frequency domains of acoustic waves.

In 2000, another mechanism leading to the existence of band gap in phononic crystals has been introduced by Z. Y. Liu et al [START_REF] Liu | Locally resonant sonic materials[END_REF], see Fig. 1.2. The phononic crystal has a simple microstructure unit consisting of a solid core material with relatively high density and a coating of elastically soft material. The highly contrast in the parameters of the two materials of the unit cell can lead to strong local resonances in the unit cell. The frequency of a local resonance can be much lower than the frequencies exhibiting Bragg scattering and, as a consequence, the band gaps could be realized for the acoustic wavelengths exceeding the lattice period by two orders of magnitude. Different from Bragg scattering mechanism where the modulation is due to non-resonant multiple scattering, in these phononic crystals the interaction of the wave with the near-resonant oscillation of the individual elastic elements of the unit cell is the origin for the existence of band gap. Thus in the phononic crystals studied in Ref. [START_REF] Liu | Locally resonant sonic materials[END_REF], the band gaps at lowest frequencies are due to their metamaterial properties (see Section 1.2.2) and not due to their phononic properties, i.e., spatial periodicity.

In addition to the band gap effect, phononic crystals also possess special wave guiding phenomena in the pass bands, for example, negative refraction [6,7,[START_REF] Lu | Negative birefraction of acoustic waves in a sonic crystal[END_REF] (Fig. 1.3(a)), negative bi-refraction [START_REF] Lu | Negative birefraction of acoustic waves in a sonic crystal[END_REF] (Fig. 1.3(b)), focusing [6,7,[START_REF] Lu | Negative birefraction of acoustic waves in a sonic crystal[END_REF][START_REF] Sukhovich | Experimental and theoretical evidence for subwavelength imaging in phononic crystals[END_REF] (Fig. 1.3(c)), surface wave trapping and guiding [START_REF] Khelif | Experimental study of guiding and filtering of acoustic waves in a two dimensional ultrasonic crystal[END_REF][START_REF] Laude | Surface acoustic wave trapping in a periodic array of mechanical resonators[END_REF][START_REF] Laude | Dispersion and polarization of surface waves trapped in high aspect ratio electrode arrays[END_REF][START_REF] Achaoui | Experimental observation of locally-resonant and Bragg band gaps for surface guided waves in a phononic crystal of pillars[END_REF], zero-refractive-index [START_REF] Zheng | Acoustic cloaking by a near-zero-index phononic crystal[END_REF][START_REF] Huang | Dirac cones induced by accidental degeneracy in photonic crystals and zero-refractive-index materials[END_REF][START_REF] Liu | Dirac cones at k=0 in acoustic crystals and zero refractive index acoustic materials[END_REF], etc.

Metamaterials

Metamaterials are artificial or engineered materials where the propagation of the waves is controlled by their interaction/hybridization with the spatially distributed resonant oscillators ("artificial atoms"). In the case of disordered distribution of the oscillators in space the dispersion/energy relations for the waves are mostly influenced near the resonance frequencies of the individual oscillators. In the case of the periodic distribution of the artificial atoms, the phononic effects related to Bragg scattering coexist with the interactions between the waves and artificial atoms. The modes of the phononic crystal can interact/hybridize with the oscillations of the artificial atoms. Thus, the single artificial atom not only can provide local resonance, but also can interact with the Bragg scattering caused by lattice periodicity, leading to much richer flexibilities in wave controlling. Thus, metamaterials, which consist of artificial atoms (b) Effective mass density as a function of frequency. [START_REF] Peano | Topological phase transitions and chiral inelastic transport induced by the squeezing of light[END_REF] in complex structures, have more features and advantages than phononic crystals composed of simple scatterers. A lot of novel effects can be achieved by metamaterials including acoustic hyperlens [START_REF] Li | Experimental demonstration of an acoustic magnifying hyperlens[END_REF], acoustic diodes [START_REF] Liang | Acoustic diode: rectification of acoustic energy flux in one-dimensional systems[END_REF][START_REF] Li | Tunable unidirectional sound propagation through a sonic-crystalbased acoustic diode[END_REF], cloaking [START_REF] Zhang | Broadband acoustic cloak for ultrasound waves[END_REF][START_REF] Stenger | Experiments on elastic cloaking in thin plates[END_REF][START_REF] Popa | Experimental acoustic ground cloak in air[END_REF][START_REF] Ergin | Three-dimensional invisibility cloak at optical wavelengths[END_REF][START_REF] Faure | Experiments on metasurface carpet cloaking for audible acoustics[END_REF], etc.

In the scientific research, the design of the artificial atom becomes the key issue that researchers need to solve for addressing a specific application of a metamaterial for wave control. In the early studies, a lot of attention was paid to the construction of the (meta)materials with effective negative material parameters, such as effective negative mass density [4] and effective negative bulk modulus [START_REF] Huang | Dirac cones induced by accidental degeneracy in photonic crystals and zero-refractive-index materials[END_REF]. As an example, a membrane-type artificial atom is shown in Fig. 1.4(a), where a flexible elastic membrane is fixed on a rigid rim [START_REF] Ma | Acoustic metasurface with hybrid resonances[END_REF][START_REF] Ma | Acoustic metamaterials: From local resonances to broad horizons[END_REF]. A rigid disk is attached to the center of the membrane. The resonance frequencies of this artificial atom can be controled by changing the dimensions of the membrane, the mass of the disk and the rigidity of the membrane (by prestress), for example. The plate including these artificial atoms (metaplate) can display effective mass negativity over certain frequency ranges, as shown in Fig. 1.4(b).

As the development of the study on metameterials, many novel wave phenomena can be realized by metameterials. For example, in electronics, a diode is a device that allows current to pass in one direction but not in the opposite one. This concept can be exploited in acoustic systems. The first idea for demonstration of acoustic asymetric propagation was to couple a one-dimensional (1D) phononic crystal made of alternating layers of glass and water to a nonlinear acoustic medium [START_REF] Liang | Acoustic diode: rectification of acoustic energy flux in one-dimensional systems[END_REF], as depicted in Fig. 1.5. A sound wave incident from the left, with frequency ω inside the band gap of the phononic crystal, is partially converted within the nonlinear medium to a secondary wave with frequency 2ω. Because ω lies in the band gap, the wave at frequency ω is reflected backward, but the secondary wave with frequency 2ω outside the band gap can freely pass through the structure. On the other hand, if the sound wave is incident from the right, the wave is completely reflected in the backward direction because ω lies inside the band gap. As a result, acoustic energy can pass through the system only from left to right, providing unidirectional sound transmission.

In recent years, the study of metamaterials has been extended to many different systems for different wavelengths. Many sophisticated designs have been implemented to achieve different novel phenomena, such as optical/acoustic black hole for absorption [START_REF] Genov | Mimicking celestial mechanics in metamaterials[END_REF][START_REF] Horstmann | Hawking radiation from an acoustic black hole on an ion ring[END_REF], zero refractive index materials for impedance matching [START_REF] Fleury | Extraordinary sound transmission through density-near-zero ultranarrow channels[END_REF], metamaterials for topological transport of waves [START_REF] Ni | Topologically protected one-way edge mode in networks of acoustic resonators with circulating air flow[END_REF], etc. At the same time, with the increasing interest of phononic crystals and metamaterials for surface/interface acoustic waves, it is believed that acoustic metamaterials with more exciting effects in different fields will continue to emerge.

Granular phononic crystals

Granular materials are conglomerations of discrete solid macroscopic particles, which can be made of different materials, size and shape. An example is a pile of sand, rice, coffee or corn flakes, etc. This kind of materials is commercially and academically important in applications as diverse as pharmaceutical industry, agriculture, energy production, geophysics, etc. The dynamics in granular materials is very complex, as they can behave as solids, liquids and gases [START_REF] Forterre | Flows of dense granular media[END_REF][START_REF] Luding | Towards dense, realistic granular media in 2D. Nonliearity[END_REF][START_REF] Herrmann | Modeling granular media on the computer[END_REF]. By assembling the particles in a periodic structure and defining the parameters (shape, size, elasticity...) of the particles, the granular materials which are called granular crystals can be created and studied. On one hand, due to the periodicity, granular crystals exhibit phonon band diagrams for elastic wave propagation as phononic crystals for acoustic waves. On the other hand, since the solid particles interact by forces and momenta applied to their geometrical contacts, granular materials also exhibit properties which are rather unique, such as rotational degrees of freedom [START_REF] Job | Solitary wave trains in granular chains: Experiments, theory and simulations[END_REF][START_REF] Theocharis | Intrinsic energy localization through discrete gap breathers in one-dimensional diatomic granular crystals[END_REF][START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF][START_REF] Chong | Nonlinear coherent structures in granular crystals[END_REF][START_REF] Texier | Helical locomotion in granular media[END_REF][START_REF] Wallen | Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal[END_REF][START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF] and strong nonlinearity [START_REF] Job | Nonlinear waves in dry and wet Hertzian granular chains[END_REF][START_REF] Daraio | Tunability of solitary wave properties in one-dimensional strongly nonlinear phononic crystals[END_REF].

Contacts between particles

In granular crystals, the interactions between elastic particles take place via their contacts. In order to study the wave dynamics in granular crystals, a description of the contact between two elastic beads is essential. The original work in contact mechanics by Heinrich Hertz dates back to 1881 [START_REF] Hertz | On the contact of rigid elastic solids[END_REF]. Hertz was attempting to understand the optical interference fringes (Newton's rings) between two lenses exhibiting an elastic deformation due to their contact interaction under pressure. Hertzian contact stress refers to the localized stresses that are developed as two In the Hertz-Mindlin theory, the contact area between two elastic spheres is small compared to the radius of the particles and planar which makes it possible to decouple the normal force-displacement relations from the tangential relations. In general, the contact area is elliptical. Figure 1.6 shows a schematics of contact between two particles. Consider the normal forces F N , the radius of the contact surface is written as,

δ = 3R * F N 4E * 1/3 , (1.1a)
with, 1

R * = 1 R 1 + 1 R 2 , 1 E * = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 , (1.1b)
where 1 and 2 represent the two particles. R is the radius of the particle, ν is the Poisson's ratio, and E is the Young's modulus of the particle. From Eqs. (1.1), the relation between the normal force F N and the deformation h of the particles is given by,

h = 3F N 4E * √ R * 2/3 . (1.2)
Equation (1.2) reveals that the deformation h is proportional to the 2/3 order of the normal force F N . The nonlinearity of the stress/strain relationship comes from the spherical geometry of the contacting particles and diminishes with the compressive force applied to the particles.

If a tangential force is applied to particles in additional to the normal force, deformation still happens in the vicinity of the contact surface. Considering that the normal force-displacement relation is decoupled from the tangential one, according to Coulomb's law, as long as the tangential force F T < µ f F N (µ f is the friction coefficient), there is no sliding on the contact surface.

While beyond this limit, slide can happen on the contact surface. The relation between the relative tangential displacement s and the tangential force F T is written as, is also nonlinear.

s = 3µ f F N 4δG * 1 -1 - F T µ f F N 2/3 , (1.3a) with, 1 G * = 2 -ν 2 1 4G 1 + 2 -ν 2 2 4G 2 , (1.3b) where G 1,2 = E 1,2 /[2(1+ν
In the case of a system in equilibrium, each contact can be considered as two springs, the first spring describing the normal interactions and the second characterising the transverse interactions. Considering the case when the two particles are identical, the normal rigidity is given by,

ξ n = 3RF N 4 1/3 E 2/3 (1 -ν 2 ) -2/3 , (1.4)
and the transverse rigidity (shear rigidity) is written, under the condition F T ≪ µ f F N , as,

ξ s = (6RF N ) 1/3 E 2/3 (1 -ν 2 1 ) 1/3 (2 -ν)(1 + ν)
.

(1.5)

Rotational degrees of freedom

When researchers study the elastic wave propagation in continuum media or in mechanical systems modelled by distribution of masses connected by the springs, usually rotational degrees of freedom are not taken into account. However, for discrete mechanical systems, the existence of auxetic materials [START_REF] Lakes | Deformation mechanisms in negative poisson's ratio materials : structural aspects[END_REF] and some chiral structures [START_REF] Spadoni | Numerical and experimental analysis of static compliance of chiral truss-core airfoils[END_REF].

In granular crystals, one important feature is that interactions between elastic particles can take place via shear forces. In contrast to There exists a number of theoretical and experimental works investigating the elastic waves in granular crystals accounting for rotational degrees of freedom of individual particles [16-19, 23, 30, 31, 81]. As one of the examples, the rotational modes have been theoretically studied in a hexagonal monolayer granular crystals [START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF], as shown in This work also proves that the Cosserat theory fails to predict correctly the dispersion of the elastic modes in granular crystals even in the long-wavelength limit.

Recently, the significance of rotations has been revealed experimentally for coupled rotationaltranslational waves in colloidal-based metamaterials [START_REF] Hiraiwa | Complex contact-based dynamics of microsphere monolayers revealed by resonant attenuation of surface acoustic waves[END_REF]. The sample is shown in Fig. 1.9(a), where a microscale granular crystal is composed of a two-dimensional monolayer of micro-sized spheres on a solid substrate. The experimental setup is presented in Fig. 1.9(b). The measurements of the resonant attenuation of laser-generated surface acoustic waves (SAWs) reveal three collective vibrational modes that involve displacements and rotations of the microspheres, as well as interparticle and particle-substrate interactions. To identify the modes, the authors tune the interparticle stiffness, which shifts the frequency of the rotational resonances while leaving the vertical normal resonance unaffected. From the measured contact resonance frequencies, they determine both interparticle and particle-substrate contact stiffnesses and find that the former is an order of magnitude larger than the latter. They found good agreement between experiment and theory.

In addition, there are other works studying elastic wave propagation in granular crystals and taking the rotational degrees of freedom into account [START_REF] Zheng | Zero-frequency and slow elastic modes in phononic monolayer granular membranes[END_REF][START_REF] Zheng | Zero-frequency and extremely slow elastic edge waves in mechanical granular graphene[END_REF][START_REF] Allein | Transversal-rotational and zero group velocity modes in tunable magneto-granular phononic crystals[END_REF][START_REF] Pichard | Surface waves in granular phononic crystals[END_REF]. Extra modes involving rotations are predicted, including rotational zero-frequency modes [START_REF] Zheng | Zero-frequency and slow elastic modes in phononic monolayer granular membranes[END_REF], rotational edge waves and rotational surface waves [START_REF] Zheng | Zero-frequency and extremely slow elastic edge waves in mechanical granular graphene[END_REF][START_REF] Pichard | Surface waves in granular phononic crystals[END_REF], etc.

Tunability and nonlinearity

In granular crystals, due to the Hertzian contact interactions between elastic particles, the dynamic vibration response of these structures is tunable from nearly linear to strongly nonlinear dynamical regimes by changing the ratio of static to dynamic interparticle displacements [START_REF] Allein | Transversal-rotational and zero group velocity modes in tunable magneto-granular phononic crystals[END_REF][START_REF] Allein | Tunable magneto-granular phononic crystals[END_REF].

Recently, the tunability of wave dynamics in granular crystals has been theoretically and experimentally studied in a one-dimensional magneto-granular chain, which is composed of spherical steel beads inside a properly designed magnetic field [START_REF] Allein | Tunable magneto-granular phononic crystals[END_REF]. This field is induced by an array of permanent magnets, located in a holder at a given distance from the chain, Fig The nonlinear effects in granular crystals have also been reported. In Ref. [START_REF] Cabaret | Nonlinear hysteretic torsional waves[END_REF], the propagation of nonlinear hysteretic torsional pulses in a vertical granular chain made of cm-scale, selfhanged magnetic beads is theoretically and experimentally studied. As predicted by contact mechanics, the torsional coupling between two beads is found to be nonlinear hysteretic. Considering a single contact between two spheres in a nonlinear resonance experiment as depicted in Fig. 1.11(d), the nonlinear hysteretic effect can be characterized at the contact. In Fig. 1.11(a), The tunablity of their dynamics from linear to strongly nonlinear make granular crystals a perfect medium for studying fundamental wave phenomena and engineering applications, in-Chapter 1: Wave propagation in periodic structures and metamaterials cluding solitary waves [START_REF] Nesterenko | Propagation of nonlinear compression pulses in granular media[END_REF][START_REF] Coste | Solitary waves in a chain of beads under Hertz contact[END_REF], discrete breathers [START_REF] Boechler | Discrete breathers in one-dimensional diatomic granular crystals[END_REF], tunable vibration filters [START_REF] Boechler | Tunable vibrational band gaps in one-dimensional diatomic granular crystals with three-particle unit cells[END_REF], acoustic diodes [START_REF] Devaux | Asymmetric acoustic propagation of wave packets via the self-demodulation effect[END_REF], acoustic rectifiers [START_REF] Boechler | Bifurcation-based acoustic switching and rectification[END_REF] and tunable functional switches [START_REF] Ganesh | From modal mixing to tunable functional switches in nonlinear phononic crystals[END_REF][START_REF] Li | Granular acoustic switches and logic elements[END_REF].

Lattice and band theory

Granular crystals are fabricated not from atoms or molecules but from macroscopic grains, which can be regarded as macroscopic analogues of atoms in a crystalline solid [START_REF] Porter | Granular crystals: Nonlinear dynamics meets materials engineering[END_REF]. Those similarities determine that the wave dynamics in granular crystals can be described by a similar method as in the atomic crystals in solid state physics [START_REF] Fjellvåg | Symmetry-operations, point groups, space groups and crystal structure[END_REF][START_REF] Galperin | Introduction to Modern Solid State Physics[END_REF]. In this section, we introduce some basic concepts in solid state physics which are also suitable for granular crystals.

Unit cell and lattice

An ideal crystal is a repetition of identical structural units in space. The periodicity is described by a lattice, which is the mathematical points at specific coordinates in space. The identical structural units (or base of the crystal) are placed at every lattice point. One can identify structural fragments that form the smallest repeating units. One can think of the smallest repeating unit as being inside a box. Such a box is called a unit cell. By stacking of the boxes in three dimensions a complete description of the structure of the solid material is obtained.

The unit cell contains one atom or more, depending on different crystal structures. In granular crystals, the atoms are the elastic spheres. The choice of a unit cell could be different, the most commonly used is the Wigner-Seitz primitive cell. Considering a lattice point is chosen, lines are drawn to all its nearest lattice points. At the midpoint of each line, another line is drawn normal to each of the first set of lines. The smallest area (or volume) enclosed by the lines is called the Wigner-Seitz primitive cell. Under this definition, a crystal is constructed by its primitive cells arranged in its lattice.

All possible lattices can be described by a set of three linearly independent vectors a 1 , a 2 , and a 3 , the unit vectors of the lattice. We can define a vector R n for the lattice,

R n = n 1 a 1 + n 2 a 2 + n 3 a 3 (1.6)
where n 1 , n 2 , and n 3 are integers. Eq. (1.6) implies the periodicity of the lattice. For a given lattice point, we could set it to be the original point with the coordinate (0, 0, 0), then any lattice point (n 1 , n 2 , n 3 ) can be reached by translation. R n is called the Bravais lattice.

In addition to the periodicity, crystals also exhibit symmetries, that is, crystals could maintain their structure after some symmetry operations such as, for example, identity, rotation, inversion, reflection. Those symmetries can be described by the point group of the lattice, and the crystal structures are classified according to their symmetry. For example, in two dimensions there are 5 distinct types of lattices as shown in Fig. 1.12. The lattice with the lowest symmetry is an oblique lattice, if |a 1 | = |a 2 | and the angle γ between a 1 and a 2 is not a rational fraction of π. crystals families [START_REF] Asbóth | A Short Course on Topological Insulators: Band Structure and Edge States in One and Two Dimensions[END_REF], see in Table 1.1. The centered lattice is special since it may also be considered as two sets of rectangular lattices packed together (shown with a dashed green and black rectangles in Fig. 1.12). More information about the lattices can be found in Ref. [START_REF] Fjellvåg | Symmetry-operations, point groups, space groups and crystal structure[END_REF] Figure 1.12: Two dimensional Bravais lattices of different symmetry. [START_REF] Fjellvåg | Symmetry-operations, point groups, space groups and crystal structure[END_REF] 1.4.2 Reciprocal space and band structure

Lattice and band theory

In crystals, any quantity, e.g. the electronic density in an atomic crystal, can be written as a periodic function,

f (r) = f (r + R n ), (1.7)
where R n is the Bravais lattice in the real space. Due to this periodicity, f (r) is possible to write in form of Fourier expansion,

f (r) = h f h (G h )e iG h r . (1.8)
Substituting into Eq. (1.7), we arrive at,

G h • R n = 2πN. (1.9)
where N is integer and the lattice constructed by G h is the reciprocal lattice. For the primitive vectors, (a 1 , a 2 , a 3 ), its reciprocal lattice vectors can be determined,

G h = h 1 • b 1 + h 2 • b 2 + h 3 • b 3 , (1.10) with b 1 = 2π a 2 × a 3 a 1 • (a 2 × a 3 ) , (1.11a) b 2 = 2π a 3 × a 1 a 2 • (a 3 × a 1 ) , (1.11b) b 3 = 2π a 1 × a 2 a 3 • (a 1 × a 2 ) , (1.11c)
where h 1 , h 2 and h 3 are integers. As R n is the Bravais lattice in the real space, G h is the Bravais lattice in the reciprocal space, which is also known as momentum space (or k-space).

The Wigner-Seitz primitive cell of the reciprocal lattice is called the first Brillioun zone. The importance of the Brillouin zone stems from the description of waves in a periodic medium, in which it is found that the solutions can be completely characterized by their behavior in a single Brillouin zone. Let us consider the vibrations in a one-dimensional diatomic granular chain. As shown in the red dash box in Fig. 1.13, the unit cell of the one dimensional (1D) diatomic granular chain contains two particles with different masses M (green) and m (black). The distance between two identical particles is 2a (lattice constant). Suppose the dynamics in the chain is linear and the particles displace only along the chain, which means that only normal forces exist in the chain. The normal rigidity of the contact between the particles is charaterized by an effective sping constant ξ n . The displacement of each bead is u i , thus the equations of motion for the 2n and 2n + 1 particles are given by,

M ü2n = ξ n (u 2n+1 + u 2n-1 -2u 2n ), (1.12a
)

mü 2n+1 = ξ n (u 2n+2 + u 2n -2u 2n+1 ). (1.12b)
Due to periodicity of the structure, we are looking for the solution of the form,

u 2n = Ae iωt-ik(2n)a , (1.13a
)

u 2n+1 = Be iωt-ik(2n+1)a , (1.13b)
where ω is the cyclic frequency, and k is the wave number. A and B are the amplitudes of the particle displacements. Substituting Eqs. (1.13) into Eqs. (1.12), the following relations can be obtained, 2ξ n -mω 2 -2ξ n cos ka

-2ξ n cos ka 2ξ n -M ω 2 A B = 0 (1.14)
For non-zero solutions, the determinant of Eq. (1.14) should be zero,

2ξ n -mω 2 -2ξ n cos ka -2ξ n cos ka 2ξ n -M ω 2 = 0. (1.15)
From Eq. (1.15), the dispersion relations (ω -k relations) can be derived, .16) where w = ω/ω c is the normalized frequency with ω c = (m + M )ξ n /mM . It predicts that there exist two branches of the elastic waves. The band structure (or dispersion curves) of the diatomic chain, revealing the ω -k relations, is shown in Fig. 1.14. Since the lattice constant is 2a, the first Brillouin zone (BZ) is [-π/2a, π/2a], which is exactly the range of k in the band structure. Actually, the dispersion curves are symmetric, even in k space. Thus, we only need to plot the part of k in the range [0, π/2a], then the band structure of the first BZ is known due to its symmetry. The corresponding range [0, π/2a] is called the irreducible Brillouin zone, which is the smallest area by reducing all the symmetries of the first BZ. In Fig. 1.14, two branches (blue and red curves) supporting the propagation of elastic modes along the chain are observed.

w 2 = 1 ± 1 m + M (m 2 + M 2 + 2mM cos 2ka) 1/2 . ( 1
Note that, there is a grey zone without any branches in it. This grey zone is a band gap for the elastic modes. It suggests that in this frequency range, waves can not propagate in the diatomic chain.

Topologically protected wave propagation

The band theory predicts that propagation of waves in a periodic structure leads to the allowed bands and forbidden gaps. Those predictions in solid state physics and in particular studying the electronic properties of materials give rise to a simple explanation of materials to be metals, semiconductors or insulators [START_REF] Galperin | Introduction to Modern Solid State Physics[END_REF]. In 1980, the discovery of the quantum Hall effect had shown that the simple division into band insulators and metals is not the end of the story, not even in band theory [START_REF] Asbóth | A Short Course on Topological Insulators: Band Structure and Edge States in One and Two Dimensions[END_REF][START_REF] Shen | Topological Insulators: Dirac Equation in Condensed Matters[END_REF][START_REF] Tkachov | Topological Insulators: The Physics of Spin Helicity in Quantum Transport[END_REF]. In the quantum Hall effect, a strong magnetic field confines the motion of electrons in the bulk, but the same field forces them into delocalized edge states on the surface.

A two-dimensional metal in strong magnetic field is thus an insulator in the bulk, but conducts along the surface, via a discrete number of edge states. The number of edge states was linked to the so-called Chern number, a topological invariant of the occupied bands. These materials are called Chern insulators [START_REF] Thouless | Quantized Hall conductance in a two-dimensional periodic potential[END_REF][START_REF] Berry | Quantal phase factors accompanying adiabatic changes[END_REF][START_REF] Hatsugai | Chern number and edge states in the integer quantum Hall effect[END_REF].

Over the last 20 years, theoretical progress over artificial systems has shown that the external magnetic field is not necessary for an insulator to have robust conducting edge states: instead, the nontrivial topology of the occupied bands is the crucial ingredient. The name topological insulator was coined for such systems, and their study became a blossoming branch of solid state physics. Since that time, the exploration of topological properties in materials and metamaterials has become an active area of research. The concept of topological insulators has also been extended to many other fields including optics [START_REF] Haldane | Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry[END_REF][START_REF] Rechtsman | Photonic Floquet topological insulators[END_REF], acoustics [START_REF] Ni | Topologically protected one-way edge mode in networks of acoustic resonators with circulating air flow[END_REF], optomechanics [START_REF] Peano | Topological phases of sound and light[END_REF][START_REF] Peano | Topological phase transitions and chiral inelastic transport induced by the squeezing of light[END_REF],

and elastic systems [START_REF] Chen | Topological mechanics of origami and kirigami[END_REF][START_REF] Paulose | Topological modes bound to dislocations in mechanical metamaterials[END_REF].

Systems supporting topologically protected wave propagation can be divided into two broad categories. The first one, analogous to the quantum Hall effect, relies on breaking time reversal symmetry (TRS) to produce chiral edge modes, that is edge waves in one dispersion branch only go one way and in the other branch go the other way [START_REF] Klitzing | New method for high-accuracy determination of the fine-structure constant based on quantized Hall resistance[END_REF][START_REF] Zhang | Experimental observation of the quantum Hall effect and Berry's phase in graphene[END_REF][START_REF] Novoselov | Two-dimensional gas of massless Dirac fermions in graphene[END_REF][START_REF] Chang | Experimental observation of the quantum anomalous Hall effect in a magnetic topological insulator[END_REF]. Generally this requires active components or the application of external fields. A second category, in analogy with the quantum spin Hall effect, preserves TRS in the system but exhibits a pair of helical edge modes on the boundaries or interfaces, that is the electrons on the edge are spin-polarized and their transport is unidirectional [START_REF] Hasan | Colloquium: Topological insulators[END_REF][START_REF] Qi | Topological insulators and superconductors[END_REF].

Time reversal symmetry breaking

The systems with non-zero Chern numbers belong to the first topological class with explicitly broken time reversal symmetry. Similar to the quantum Hall (QH) effect, in those systems, the bulk of the systems is insulating, and the transport of waves takes place only along the edge.

The flow of the edge waves is unidirectional and robust against defects. To achieve similarity to QH effect, the challenge is to efficiently break the TRS and induce non-zero Chern numbers to the bands. Note that, there are publications showing that even if TRS is preserved, QH effect (anomalous QH effect) still can be achieved [START_REF] Wei | Experimental demonstration of topologically protected efficient sound propagation in an acoustic waveguide network[END_REF][START_REF] Hafezi | Imaging topological edge states in silicon photonics[END_REF][START_REF] Hafezi1 | Robust optical delay lines with topological protection[END_REF]. For example, in Ref. [START_REF] Hafezi1 | Robust optical delay lines with topological protection[END_REF] coupled resonators supporting degenerate clockwise and counter-clockwise modes are needed in order to induce phase shifts for the degenerate modes and synthesize an effective magnetic field. We will discuss below only the case when TRS is broken.

The topologically protected wave propagation by breaking TRS has been demonstrated in different systems. The first experimental observation in photonic system was reported in 2009 by Wang et al [START_REF] Wang | Observation of unidirectional backscattering-immune topological electromagnetic states[END_REF]. In their experiment, a 2D photonic crystal is constructed by arranging gyromagnetic ferrite rods in a square lattice, Fig. 1.15(a). A uniform magnetic field is applied to break TRS. Two metallic plates and a metal wall are added to prevent radiation loss into air. In the presence of the magnetic field, a complete band gap appears between the second and third bands, which both have non-zero Chern numbers, Fig. 1.15(b). In addition, a gapless dispersion curve in red is observed. It supports the propagation of edge state inside the second band gap, which has only positive group velocity. The simulation of edge wave propagation on the edge with time is shown in Fig. 1.15(c), where the topologically protected edge mode can travel in only one direction and is very robust against reflections even when a large metallic scatterer is placed in the path of the propagation. Fig. 1.15(d) shows the transmission data of the edge waves for the forward and backward (red) cases, indicating that edge wave can not propagate in the backward direction. In general, the gyromagnetic response of materials is very weak in the optical regime. Thus, realization of topological propagation at optical frequencies remains a challenge [START_REF] Lu | Topological photonics[END_REF]. In 2015, topological edge states in an acoustic structure containing circulating fluids were reported by Z. Yang et al [START_REF] Yang | Topological acoustics[END_REF]. In their work, the acoustic structure is composed of a honeycomb phononic crystal, where the rigid solid cylinders are surrounded by the circulating fluids, Recently, topological edge waves in a different mechanical system with TRS breaking have been demonstrated [START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF][START_REF] Nasha | Topological mechanics of gyroscopic metamaterials[END_REF]. As shown in Fig. 1.17 (a) and (b), the designed structure is a hexagonal phononic crystal where the mass on the lattice point is connected to a gyroscope [122]. Due to the gyroscopic inertial effects, TRS of the system is broken, and the mechanical analogue of the electronic QH effect is realized. The band structure of the gyroscopic phononic crystal is shown in Fig. 1.17(c), where bands of non-zero Chern number are exhibited. The dispersion curves of edge wave in the bulk gaps are shown in Fig. 1.17(d) by the color lines.

The robustness of those edge states is confirmed by the simulation in Fig. 1.17(e), where the edge modes circumvent both the sharp corner and the line defect without any reflection. The experimental observation of the topological edge wave in the gyroscopic phononic crystal was later reported by L. M. Nash et al [START_REF] Nasha | Topological mechanics of gyroscopic metamaterials[END_REF].

Time reversal symmetry invariance

In 2005, Kane and Mele considered the spin orbit coupling in graphene [START_REF] Kane | Quantum spin Hall effect in graphene[END_REF][START_REF] Kane | Z 2 topological order and the quantum spin Hall effect[END_REF]. At sufficently low energy, they found that a single plane of graphene exhibits a quantum spin Hall (QSH) effect.

The QSH phase is a time reversal invariant electronic state with a bulk band gap that supports the transport of charge and spin in gapless edge states. They introduced a Z 2 topological index to characterize the system. In their work, the QSH phase exhibits a bulk energy gap with a pair of gapless spin polarized edge states in it. Since TRS invariant prevents the flip of spin states, those edge states are robust against small perturbations. However, because the spin-orbit coupling is extremely small, it is very difficult to observe the QSH effect in graphene.

Around the same year, B. A. Bernevig, T. L. Hughes and S. -C. Zhang proposed that the QSH effect can be realized in CdTe/HgTe semiconductor quantum wells (BHZ model) [START_REF] Bernevig | Quantum spin Hall effect and topological phase transition in HgTe quantum wells[END_REF]. The CdTe layer has a normal type of band structure, while the HgTe layer has inverted bulk bands.

By controlling the thickness of the well layer (HgTe), a critical phase transition happens tuning a conventional insulating phase to a topological one, leading to QSH effect in the quantum well. Later, König et al experimentally confirmed that a thin layer of HgTe is a topological insulator [START_REF] König | Quantum spin hall insulator state in HgTe quantum wells[END_REF]. Since then, a new topological class of materials exhibiting QSH effect are widely studied. Analogues of electronic QSH effect to achieve topological transport hence have been Unidirectional property and the robustness of edge waves against defects. [START_REF] Khanikaev | Photonic topological insulators[END_REF] extended to many other research fields. To achieve topological wave propagation with TRS, the key point is to construct the time reversal pair, for example, the spin-up and spin-down states in electronic systems.

The first photonic analogy of the QSH effect with TRS invariant was proposed in 2013 by A. B. Khanikaev et al [START_REF] Khanikaev | Photonic topological insulators[END_REF]. In their work, the time reversal pair is constructed by linear combinations of light with different polarizations: the TE+TM modes refers to as the spin '+' state, and the the TE-TM modes correspond to the spin '-' state, Fig. 1.18(a). Those modes can be achieved by metamaterials containing the split-ring resonators, Fig. 1.18(b). By arranging resonators in a honeycomb lattice, Dirac dispersion at K points is obtained, and lifted by inducing the bi-anisotropy to mimic the effective spin-orbit coupling, Fig. 1.18(c). On the interfaces of the metacrystals, a pair of spin polarized helical edge states is observed, Fig. 1.18(d).

This pair of edge states supports the one-way propagation of photons and robust against different kinds of disorder, Fig. 1.18(e). Later in 2014, this proposal of photonic topological insulator was experimentally realized by C. T. Chan's team [START_REF] Chen | Experimental realization of photonic topological insulator in a uniaxial metacrystal waveguide[END_REF].

In 2015, X. Hu's team proposed other way to achieve the spin analogy: the doubly degenerate Bloch modes due to the C 6 symmetry of the system [START_REF] Wu | Scheme for achieving a topological photonic crystal by using dielectric material[END_REF]. By deforming a honeycomb lattice of cylinders into a triangular lattice of cylinder hexagons, a double Dirac cone can be obtained and lifted at Gamma point by adjusting the lattice, leading to a band inversion and the topological phase transition. Later in 2016, H. Cheng et al demonstrated experimentally this design in acoustics [START_REF] He | Acoustic topological insulator and robust one-way sound transport[END_REF], see Fig. 1.19.

There are other designs of the spin analogy in photonic systems, in addition to using linear combinations of TE/TM modes [START_REF] Khanikaev | Photonic topological insulators[END_REF][START_REF] Chen | Experimental realization of photonic topological insulator in a uniaxial metacrystal waveguide[END_REF][START_REF] Ma | Guiding electromagnetic waves around sharp corners: Topologically protected photonic transport in metawaveguides[END_REF], for example, using the left circular polarization and right circular polarization of light to be the pseudo-spin states [START_REF] He | Photonic topological insulator with broken time-reversal symmetry[END_REF]. A pair of topological edge waves also can be observed in those systems.

In mechanical system, topological edge waves under TRS have also been reported. A dualscale phononic crystal slab as an elastic analogy of the QSH effect was demonstrated by S. H. Mousavi et al in 2015 [START_REF] Mousavi | Topologically protected elastic waves in phononic metamaterials[END_REF], Fig. 1.20. Two types of elastic modes are used to be the effective spins for phonons, and strong spin-orbit coupling is realized by breaking spatial mirror symmetry in z-direction. Thus topological gap is observed which supports the propagation of a pair of edge waves in it. In the same year, a mechanical topological insulator exhibiting QSH effect was experimentally demonstrated, Fig. 1.21. The mechanical system comprises of 1D pendula [START_REF] Ssstrunk | Observation of phononic helical edge states in a mechanical topological insulator[END_REF].

By carefully introducing couplings between pendula, the system exhibits Z 2 topological order and supports the topologically protected propagation of edge wave. During the recent studies of the topological effect in materials, many other topological phenomena have also been reported [START_REF] Mei | Pseudo-time-reversal symmetry and topological edge states in two-dimensional acoustic crystals[END_REF][START_REF] Deymier | Torsional topology and fermion-like behavior of elastic waves in phononic structures[END_REF][START_REF] Xiao | Synthetic gauge flux and Weyl points in acoustic systems[END_REF][START_REF] Lu | Valley vortex states in sonic crystals[END_REF][START_REF] Dong | Valley photonic crystals for control of spin and topology[END_REF][START_REF] Fu | Topological Crystalline Insulators[END_REF][START_REF] Sato | Topological superconductors: a review[END_REF], such as metamaterials with Weyl point [START_REF] Kane | Z 2 topological order and the quantum spin Hall effect[END_REF], analogues of valleytronic effect [START_REF] Lu | Valley vortex states in sonic crystals[END_REF][START_REF] Dong | Valley photonic crystals for control of spin and topology[END_REF], topological crystalline insulator [START_REF] Fu | Topological Crystalline Insulators[END_REF], topological superconductor [START_REF] Sato | Topological superconductors: a review[END_REF], etc. In addition, topological transition has been observed also in isostatic lattices [START_REF] Kane | Topological boundary modes in isostatic lattices[END_REF][START_REF] Sun | Surface phonons, elastic response, and conformal invariance in twisted kagome lattices[END_REF] and auxetic materials [START_REF] Kolken | Auxetic mechanical metamaterials[END_REF][START_REF] Zeb Rocklin | Transformable topological mechanical metamaterials[END_REF]. We believe that more and more phenomena associated with topology will continue to appear.

Topological band theory

As the classical band theory fails in prediction of the topological insulators, topological band theory has been exploited [START_REF] Bansil | Colloquium: Topological band theory[END_REF][START_REF] Kane | Topological band theory and Z 2 invariant[END_REF]. In this section, we will provide the foundations of topological band theory and explain how it can be used to characterize the topological phenomena.

Topological invariant

To describe the unique properties of wave propagation in topological non-trivial systems, topology has emerged as another essential degree of freedom. Mathematically, topology is introduced to classify different geometrical objects. Different topologies can be mathematically characterized by topological invariants, quantities that remain constant under arbitrary continuous deformations. For example, as shown in Fig. 1.22(a), the six objects have different geometries, but there are only three different topologies. The sphere can be smoothly deformed into the spoon without creating any hole. Topologically, they are equivalent. The torus and coffee cup are also topologically equivalent, so are the double torus and tea pot. For the above closed surfaces, the topological invariant is its genus, which corresponds to the number of holes within a closed surface. Objects with the same topological invariant are topologically equivalent; that is, they are in the same topological phase. Only when a hole is created or removed in the object does The magnitude of ∆C equals the number of gapless edge modes. [START_REF] Lu | Topological photonics[END_REF] the topological invariant change. This process is known as a topological phase transition.

The most fascinating and peculiar phenomena take place on the interface between two topological insulators with different topological invariants. The interface formed by these two topological insulators (Fig. 1.22(b), right) is topologically distinct from an ordinary interface, which is formed between topologically equivalent materials (Fig. 1.22(b), left). The difference lies in the dispersion spectra of the edge modes inside the bulk gap. On the left of Fig. 1.22c, the two bulk bands both have zero Chern numbers, so they can directly connect across the interface without closing the frequency gap. However, when the two materials have different Chern numbers, topology does not allow them to connect to each other directly. A topological phase transition must take place on the interface: this requires the closing of the frequency gap, neutralizing the Chern numbers, then reopening the gap. This phase transition (Fig. 1.22(c), right), ensures gapless edge states on the interface, edge states at all frequencies within the gap of the bulk. The gapless dispersion relations of the edge states are topologically protected, that is, their existence is guaranteed by the difference of the topologies of the bulk materials on the two sides. In general, the number of gapless edge modes equals the difference of the bulk topological invariants across the interface. This is known as the bulk-edge correspondence.

The topological protection of edge waves can be understood in reciprocal space. Fig. 1.22(d) shows the dispersion diagrams of both ordinary (left) and gapless (right) interfaces. On the left, the dispersion of the edge wave of the ordinary interface is disconnected from the bulk bands and can be continuously moved out of the frequency gap into the bulk bands by smoothly modifying the positions of the bulk bands. On the right, however, the gapless dispersion connects the bulk bands above and below the gap. It cannot be removed out of the gap by modifying the position of the lower and upper bulk bands. The only way to alter these connectivities is through a topological phase transition, that is, closing and reopening the bulk gaps.

The unidirectional property of the protected edge modes can also be explained from the slopes (group velocities) of the edge dispersions. An ordinary interface (Fig. 1.22(d), left) supports edge modes along both directions because it can be back-scattered (the same mode can be excited with opposite sign of wavevector). In contrast, a protected gapless edge mode on the topological interface (Fig. 1.22(d), right) is unidirectional as it has only positive (or only negative) group velocities (this mode can not be excited when the sign of the wavevector is opposite). It is obvious that there are no counter-propagating modes at the same frequencies for the one-way edge modes. This enables wave to flow around imperfections without backscattering since the wave can go only forward.

Berry phase and Chern number

Mathematically, the geometrical objects, undergoing smooth deformations into different shapes, are topologically equivalent unless holes are created in the deformation process. This notion can be extended to the topological classification of insulators [START_REF] Dauphin | Efficient algorithm to compute the Berry conductivity[END_REF][START_REF] Xiao | Berry phase effects on electronic properties[END_REF]. Suppose that the geometrical objects now are insulators, the smooth deformations correspond to the slow changes of the Hamiltonian, which is an adiabatic process. As long as the systems remain in the ground state, that is, no energy gap closing and reopening takes place, the insulators maintain topological equivalence in this process.

As we have discussed, the dynamical equation of a periodic structure usually can be expressed as an eigenvalue problem,

D(k)v(k) = Ωv(k), (1.17)
where D(k) is the dynamical matrix of the system. In electronic system, D(k) is the Hamiltonian, Ω is the energy and v(k) is the state of the system. The classical band theory cares only about the eigenvalues, Ω, which identify the band structure of the system. However, the information in the eigenvector v(k) is ignored. Eventually, it turns out that the topological property of the bands has close relations to the eigenvectors. To classify the topological invariant of bands, Berry phase plays a key role in topological band theory. The first work about Berry phase was published in 1984 by Michael Berry, when he was studying the adiabatic evolution of an eigenenergy state in a slowly changing quantum system. This notion can be extended to the band theory [START_REF] Boechler | Tunable vibrational band gaps in one-dimensional diatomic granular crystals with three-particle unit cells[END_REF]. Due to the periodicity in the crystals, under translation the Bloch states become,

v(k) → e iΘ(k) v(k). (1.18)
Θ(k) is a phase due to translation. According to the Bloch's theorem, Bloch states are invariant under the transformation in Eq. (1.18). We can define the Berry connection in this process,

A = i v(k)|∇ k |v(k) . (1.19) Under Eq. (1.18), it transforms as, A → A + ∇ k Θ(k). It is clear that A, depending on Θ(k),
is not gauge invariant. It is useful to define a gauge invariant vector, the Berry curvature, as

B = ∇ k × A.
Then, the Berry phase is defined,

Γ C = C dk • A = S dS • B, (1.20)
where C represents the close path of the BZ, and S is the surface of the BZ enclosed by the path C. Apparently, the Berry phase is gauge invariant.

To see the physical meaning of the Berry phase, let us consider a two-level system as an example [START_REF] Dauphin | Efficient algorithm to compute the Berry conductivity[END_REF][START_REF] Xiao | Berry phase effects on electronic properties[END_REF]. The Hamiltonian of the system can be given by

H(k) = e(k) • σ = e z e x -ie y e x + ie y -e z , (1.21)
where e(k) = (e x , e y , e z ) and σ represent the Pauli matrices,

σ x = 0 1 1 0 , σ y = 0 -i i 0 , σ z = 1 0 0 -1 . (1.22)
The corresponding eigenstates with the energy ±|e| are obtained, Let us consider the v -mode. According to Eq. (1.19), the Berry connection is given by,

v + =    cos θ 2 e -iφ sin θ 2    , v -=    sin θ 2 e -iφ -cos θ 2    , (1.23
A θ = i v -|∂ θ v -= 0, (1.25a) A φ = i v -|∂ φ v -= sin 2 θ 2 , (1.25b)
and the Berry curvature of v -mode in the spherical coordinates is, .26) where e = e/|e| is a unit vector and can be identified by a point on the sphere S 2 . According to Eq. (1.20), the Berry phase for this system reads,

B - θφ = (∂ θ A φ -∂ φ A θ ) e = 1 2 sin θ e. ( 1
Γ C = 1 2 S sin θ edθdφ = 1 2 S e |e| 3 dS.
(1.27)

Eq. (1.27) shows that the Berry curvature B = e/2|e| 3 can be regarded as a gauge field generated by a monopole at the origin point e = 0, which acts as source and drain of the Berry curvature flux. We find that the integral of the Berry curvature over a closed 2D space S,

1 2π S dθdφB - θφ = 1. (1.28)
It shows that the Berry curvature integrated over a closed manifold S is a multiple of 2π that is equal to the net number of monopoles inside. This number defines a topological invariant called the Chern number,

C = 1 2π S B • dS. (1.29)
For 2D crystals, the closed surface S is the 2D BZ. The topological band theory indicates that the eigenvectors need also to be computed and can indeed have a remarkable effect on the band structure.

Conclusion

As a new kind of artificial and functional materials/structures, phononic crystals and metamaterials exhibit fascinating abilities in wave control. Many novel phenomena can be achieved based on those functional materials/structures. In the recent years, the studies of granular crystals and metamaterials also show excellent abilities in the controlling of elastic waves [START_REF] Job | Solitary wave trains in granular chains: Experiments, theory and simulations[END_REF][START_REF] Theocharis | Intrinsic energy localization through discrete gap breathers in one-dimensional diatomic granular crystals[END_REF][START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF][START_REF] Chong | Nonlinear coherent structures in granular crystals[END_REF][START_REF] Texier | Helical locomotion in granular media[END_REF][START_REF] Wallen | Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal[END_REF][START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF].

In addition, due to interactions taking place via their interparticle contacts, they exhibit some advantages in the study of rotational waves [START_REF] Job | Solitary wave trains in granular chains: Experiments, theory and simulations[END_REF][START_REF] Theocharis | Intrinsic energy localization through discrete gap breathers in one-dimensional diatomic granular crystals[END_REF][START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF] and nonlinear waves [START_REF] Allein | Transversal-rotational and zero group velocity modes in tunable magneto-granular phononic crystals[END_REF][START_REF] Allein | Tunable magneto-granular phononic crystals[END_REF]. In this work, we continue to exploit the linear wave dynamics in the two-dimensional (2D) granular crystals.

Propagation of slow waves, zero-frequency modes and edge waves are studied in Chapters 2 and 3 [START_REF] Zheng | Zero-frequency and slow elastic modes in phononic monolayer granular membranes[END_REF][START_REF] Zheng | Zero-frequency and extremely slow elastic edge waves in mechanical granular graphene[END_REF]. In addition, in Chapter 4 topological wave propagation phenomena are demonstrated in granular crystals for rotational waves [START_REF] Zheng | Quasi-topological rotational waves in mechanical granular graphene[END_REF].

Chapter 2

Bulk wave propagation in granular crystals Contents

Introduction

There has been a growing interest in investigating the propagation of elastic/acoustic waves in phononic crystals in the last decade [1][2][3][4][5][6][7][START_REF] Khelif | Experimental study of guiding and filtering of acoustic waves in a two dimensional ultrasonic crystal[END_REF][START_REF] Laude | Surface acoustic wave trapping in a periodic array of mechanical resonators[END_REF][START_REF] Laude | Dispersion and polarization of surface waves trapped in high aspect ratio electrode arrays[END_REF][START_REF] Achaoui | Experimental observation of locally-resonant and Bragg band gaps for surface guided waves in a phononic crystal of pillars[END_REF][START_REF] Zigoneanu | Three-dimensional broadband omnidirectional acoustic ground cloak[END_REF][START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF][START_REF] Sigalas | Band structure of elastic waves in two dimensional systems[END_REF][START_REF] Martinez-Salar | Sound attenuation by sculpture[END_REF][START_REF] Liu | Locally resonant sonic materials[END_REF]. Owing to Bragg scattering in spatially periodic media, phononic crystals possess exotic dispersion characteristics including, for example, frequency band gaps [START_REF] Phani | Wave propagation in two-dimensional periodic lattices[END_REF], negative refraction [START_REF] Lu | Negative birefraction of acoustic waves in a sonic crystal[END_REF], subwavelength imaging [START_REF] Sukhovich | Experimental and theoretical evidence for subwavelength imaging in phononic crystals[END_REF] etc. Although studies of phononic crystals involving elastic behaviors have been reported in all dimensions, from 1D to 3D [START_REF] Liu | Locally resonant sonic materials[END_REF][START_REF] Pichard | Localized transverserotational modes in linear chains of equal masses[END_REF][START_REF] Carta | Dispersion properties of vortex-type monatomic lattices[END_REF], most of the recent studies are focusing on different types of twodimensional lattice that support bulk or edge modes [START_REF] Spadoni | Phononic properties of hexagonal chiral lattices[END_REF][START_REF] Torrent | Elastic analog of graphene: Dirac cones and edge states for flexural waves in thin plates[END_REF][START_REF] Brun | Vortex-type elastic structured media and dynamic shielding[END_REF]. 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 [START_REF] Souslov | Elasticity and response in nearly isostatic periodic lattices[END_REF][START_REF] Mao | Soft modes and elasticity of nearly isostatic lattices: randomness and dissipation[END_REF]. More interestingly, when the kagome lattices are twisted, negative Poisson ratio and zero-frequency edge states could be achieved [START_REF] Kane | Topological boundary modes in isostatic lattices[END_REF][START_REF] Sun | Surface phonons, elastic response, and conformal invariance in twisted kagome lattices[END_REF]. Other fascinating elastic properties, such as topological soft modes and topological edge modes, have also been reported in the kagome and honeycomb systems [START_REF] Paulose | Topological modes bound to dislocations in mechanical metamaterials[END_REF][START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF].

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 inherently much softer than the grains [START_REF] Nicoletti | Acoustic diodes[END_REF][START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF][START_REF] Chong | Nonlinear coherent structures in granular crystals[END_REF][START_REF] Texier | Helical locomotion in granular media[END_REF][START_REF] Wallen | Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal[END_REF][START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF], 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 [START_REF] Akimov | Hypersonic modulation of light in three-dimensional photonic and phononic band-gap materials[END_REF][START_REF] Lisiecki | Vibration dynamics of supra-crystals of cobalt nanocrystals studied with femtosecond laser pulses[END_REF][START_REF] Ayouch | Elasticity of an assembly of disordered nanoparticles interacting via either van der waals-bonded or covalent-bonded coating layers[END_REF] 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 [START_REF] Paulose | Topological modes bound to dislocations in mechanical metamaterials[END_REF][START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF][START_REF] Kane | Topological boundary modes in isostatic lattices[END_REF][START_REF] Sun | Surface phonons, elastic response, and conformal invariance in twisted kagome lattices[END_REF],

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 [START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF][START_REF] Chong | Nonlinear coherent structures in granular crystals[END_REF][START_REF] Texier | Helical locomotion in granular media[END_REF][START_REF] Wallen | Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal[END_REF][START_REF] Pichard | Localized transverserotational modes in linear chains of equal masses[END_REF] and experimentally demonstrated [START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF][START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF] 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 for 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 zero-frequency (zero-energy) modes [START_REF] Schwartz | Vibrational Modes in granular materials[END_REF][START_REF] Maxwell | On the calculation of the equilibrium and stiffness of frames[END_REF][START_REF] Calladine | Buckminster fuller's "tensegrity" structures and clerk Maxwell's rules for the construction of stiff frames[END_REF]. 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.

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][2][3][4][5], 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 [START_REF] Chang | A micromechanical-based micropolar theory for deformation of granular solids[END_REF][START_REF] Eringen | Microcontinuum Field Theories. 5: Theory of Micropolar Elasticity[END_REF][START_REF] Kanatani | A micropolar continuum theory for the flow of granular materials[END_REF].

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 [16-22, 81, 160, 166]. In this Chapter, we exploit the Lagrangian method [START_REF] Suiker | Comparison of wave propagation characteristics of the Cosserat continuum model and corresponding discrete lattice models[END_REF] 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 Γ, 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 [START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF][START_REF] Pichard | Surface waves in granular phononic crystals[END_REF]. Besides, the detailed analysis of zero-frequency modes is reported. In the honeycomb lattices, we demonstrate the existence of zero-groupvelocity 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 Γ 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 [START_REF] Kariyado | Manipulation of Dirac cones in mechanical graphene[END_REF][START_REF] Chen | Accidental degeneracy of double Dirac cones in a phononic crystal[END_REF][START_REF] Mei | First-principles study of Dirac and Dirac-like cones in phononic and photonic crystals[END_REF][START_REF] Zhang | Observing Zitterbewegung for photons near the Dirac point of a two-dimensional photonic crystal[END_REF]. 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 [START_REF] Wang | Coriolis force induced topological order for classical mechanical vibrations[END_REF][START_REF] Ni | Topologically protected one-way edge mode in networks of acoustic resonators with circulating air flow[END_REF][START_REF] Yang | Topological acoustics[END_REF]. Even for the Dirac cone at Γ point, one could expect to observe the pseudospinresolved Berry curvatures of photonic bands and helical edge states characterized by Poynting vectors [START_REF] Wu | Scheme for achieving a topological photonic crystal by using dielectric material[END_REF][START_REF] Ssstrunk | Observation of phononic helical edge states in a mechanical topological insulator[END_REF]. The analytical predictions of Dirac cones and double Dirac cones in this

Chapter 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 [START_REF] Boechler | Interaction of a contact resonance of microspheres with surface acoustic waves[END_REF][START_REF] Khanolkar | A selfassembled metamaterial for Lamb waves[END_REF][START_REF] Maznev | Waveguiding by a locally resonant metasurface[END_REF][START_REF] Wallen | Dynamics of a monolayer of microspheres on an elastic substrate[END_REF].

This Chapter is constructed as follows: in Sec. 2.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 Sec. 2.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 Sec. 2.4. In Sec. 2.5 and Sec. 2.6, we turn our attention to the honeycomb monolayer membranes for both the cases of out-of-plane and in-plane motions.

In Sec. 2.7, we present the conclusions of this chapter.

Interactions and dynamics in granular crystals

In this section, we discuss the individual movements of particle and the possible interactions between beads in two-dimensional infinite monolayer granular crystals. The equations of motion of hexagonal and honeycomb structures are derived. In two-dimensional infinite monolayer granular crystals, since rotation and translation of beads can take place, each individual particle possesses 6 degrees of freedom (DOF). However, when proper motions are imposed to the system, not all of them contribute to the elastic wave propagation. In this Chapter, we consider two different types of motion in granular structures as depicted in Fig. These forces are activated in the granular crystals 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.

Individual grain movements and intergrain interactions

(2) Torsional moments, which are characterized by an effective torsional (spin) rigidity . This type of interaction can be excited when there is relative displacement between two adjacent beads along the axis connecting their centers. (4) Bending moments, which are characterized by an effective bending rigidity (Fig. 2.2(c)) resulting from the distribution of the normal forces at the contact, when the beads are rotated relative to the axis parallel to the contact. They originate from the resistance of the contact of beads to rolling. 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.2(a) the effective shear rigidity of the contact is represented by a single spring of an effective rigidity ξ s 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 present the distributed shear interaction at the contact, caused by the relative rotations of the grains, as one supported by two separated springs of half rigidity separated 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 δ.

This presentation of the shear interaction between the contact faces does not change the magnitude of shear interaction, which will be the same as presented in Fig. 2.2(a), because the sum of these two forces is equal to zero. However, it reveals the existence of the torsional rigidity of the contact. From Fig. 2.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.2(b)). It is worth noting that the elastic energy stored in shear interaction,

presented in Fig. 2.2(a), is ∼ ξ s (u 2 -u 1 ) 2 .
The elastic energy stored in shear interaction caused by the rotations in Fig.

2.2(a) is ∼ ξ s [R(ψ 2 + ψ 1 )] 2
, where ψ i (i=1, 2) is the rotational angle shown in Fig. 2.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.

2(b), i.e., ∼ ξ s [δ(ψ 2 -ψ 1 )] 2 , should be presented as ∼ ξ s (δ/R) 2 [R(ψ 2 -ψ 1 )] 2
, demonstrating that the effective torsional rigidity of the contact is related to its shear rigidity by ξ t ∼ ξ s (δ/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.2(c), in order to reveal the existence of bending interaction, we have splitted the normal rigidity of the contact (presented in Fig. 2.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.2(c) and (d). However, the spatially separated springs in Fig. 2.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, ξ b ∼ ξ n (δ/R) 2 . The above derived estimates indicate that torsional and bending rigidities could be much smaller than shear and normal rigidities, respectively, i.e., ξ t ≪ ξ s and ξ b ≪ ξ n , in non-consolidated granular crystals, where commonly δ ≪ 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 ξ t ξ s and ξ b ξ n cannot be a priori excluded from the theoretical analysis. Because commonly normal rigidity is comparable with shear rigidity, i.e., ξ n ξ s , then the granular crystals with ξ b ξ t cannot be excluded from the theoretical analysis either.

Lagrangian formalism in granular crystals

The interactions in granular crystals can be theoretically analysed by the Lagrangian formalism [START_REF] Souslov | Elasticity and response in nearly isostatic periodic lattices[END_REF]. In general, there are connections only between beads with their nearest neighbors, thus we only consider the nearest couplings in granular crystals. The Lagrangian of beads can be written as,

L = E k -E p , (2.1)
where the kinetic energy E k and the potential energy E p read,

E k = 1 2 j M j χ2 j , (2.2)
E p = 1 2 l,i ξ l ∆l 2 i . (2.3)
where j = x, y, z is the coordinate index, l = s, n, b, t is the coupling index denoting the shear/normal/bending/torsional coupling, and i is the neighbor index of the particle labelling its i-th neighbor particle. M j is the mass/inertia of bead. χ j represents the movement component. • Out-of-plane motion

In this case, shown in Fig. 2.2, the shear, torsional and bending interactions take place in the granular crystals.

Firstly, we define several unit vectors as depicted in Fig. 2.3, where e i is the unit vector in the direction from the center of beads to the center of their i-th neighbor. e x , e y and e z represent the unit vectors along x-, y-and z-axis, respectively. d i is an unit vector in the direction normal to e i and e z . Thus, d i takes the form,

d i = e z × e i .
Under the definitions, the effective elongations ∆l i of those three types of springs are given by,

∆s i = (u i -u) - a 2 (ϕ i + ϕ)e x • d i - a 2 (ψ i + ψ)e y • d i , (2.5a
)

∆b i = a 2 (ϕ i -ϕ)e x • d i + a 2 (ψ i -ψ)e y • d i , (2.5b
)

∆t i = a 2 (ϕ i -ϕ)e x • e i + a 2 (ψ i -ψ)e y • e i , (2.5c)
where ∆s i , ∆b i and ∆t i correspond to the effective elongations of the effective springs accounting for the shear, bending and torsional contact interactions between beads with their i-th neighbor.

• In-plane motion

In this case, shown in Fig. 2.2, the normal, shear, and bending interactions are activated in the granular crystals. the effective elongations ∆l i of those three types of springs are given by,

∆n i = (u xi -u x )e x • e i + (u yi -u y )e y • e i , (2.6a
)

∆s i = (u xi -u x )e x • d i + (u yi -u y )e y • d i - a 2 (ϕ i + ϕ), (2.6b
)

∆b i = a 2 (ϕ i -ϕ), (2.6c)
where ∆n i , ∆s i and ∆b i correspond to the effective elongations of the effective springs accounting for the normal, shear, and bending contact interactions between beads.

As long as the structure of the granular crystal is given, the elongations of springs and the equations of motion can be obtained. In the next, we analyse separately the hexagonal and honeycomb structures with both out-of-plane and in-plane motions. A hexagonal granular membrane, a two-dimensional infinite monolayer granular crystal, is shown in Fig. 2.4, where periodically ordered spherical particles with radius R are arranged in a hexagonal lattice. The structure is characterized by the lattice constant a = 2R for the hexagonal lattice. The corresponding first Brillouin zones is also depicted in Fig. 2.4.

Equations of motion of hexagonal granular crystals

Considering the particle number "0" in Fig. 2.4, it has 6 nearest neighbors labeled by "1" to "6". Next, we discuss the cases of out-of-plane and in-plane motions, respectively.

• Out-of-plane motion According to Eqs. (2.1) -(2.4), the equations of motion are derived as,

I φ0 = R 6 i=1 (ξ s ∆s i e x • d i + ξ b ∆b i e x • d i + ξ t ∆t i e x • e i ), (2.7a
)

I φ0 = R 6 i=1 (ξ s ∆s i e y • d i + ξ b ∆b i e y • d i + ξ t ∆t i e y • e i ), (2.7b
)

M ü0 = 6 i=1
ξ s ∆s i .

(2.7c)

• In-plane motion
Similarly, the equations of motion are given by,

M üx0 = 6 i=1 (ξ n ∆n i e x • e i + ξ s ∆s i e x • d i ), (2.8a)
M üy0 = 6 i=1 (ξ n ∆n i e y • e i + ξ s ∆s i e y • d i ), (2.8b
)

I φ0 = R 6 i=1 (ξ s ∆s i + ξ b ∆b i ).
(2.8c)

In Eqs. (2.7) and (2.8), we need to calculate the inner produces of e x,y with e i and d i . Table 2.1 shows the results for the hexagonal structure. Thus, the elongations of effective springs and the explicit forms of equations of motion can be obtained. One should notice that, in the honeycomb structure, there are particles in two different positions marked as A and B. For example, as presented in Fig. 2.5, particles number "0", "2"

e x • e i 1/2 1 1/2 -1/2 -1 -1/2 e y • e i - √ 3/2 0 √ 3/2 √ 3/2 0 - √ 3/2 e x • d i √ 3/2 0 - √ 3/2 - √ 3/2 0 √ 3/2 e y • d i 1/2 1 1/2 -1/2 -1 -1/2

Equations of motion of honeycomb granular membranes

and "4" belong to particles in position A, while particles number "1", "3" and "5" are particles in position B. The unit cell of the honeycomb structure can be chosen as shown in the red box in Fig. 2.5. Considering the particle "0" and "1", each of them has three nearest neighbors.

To obtain the equations of motion, we analyse the cases of out-of-plane and in-plane motions, respectively.

• Out-of-plane motion Since there are two particles in one unit cell, for example, particles 0 and 1, we need to analyse all interactions including at least one of these particles. According to Eqs. (2.1) -(2.4), the equations of motion are derived as,

I φ0 = R i=1,3,5 (ξ s ∆s i e x • d i + ξ b ∆b i e x • d i + ξ t ∆t i e x • e i ), (2.9a)
I φ0 = R i=1,3,5 (ξ s ∆s i e y • d i + ξ b ∆b i e y • d i + ξ t ∆t i e y • e i ), (2.9b
)

M ü0 = i=1,3,5 ξ s ∆s i , (2.9c
)

I φ1 = R i=0,2,4 (ξ s ∆s i e x • d i + ξ b ∆b i e x • d i + ξ t ∆t i e x • e i ), (2.9d
)

I φ1 = R i=0,2,4 (ξ s ∆s i e y • d i + ξ b ∆b i e y • d i + ξ t ∆t i e y • e i ), (2.9e
)

M ü1 = i=0,2,4
ξ s ∆s i .

(2.9f)

• In-plane motion Similarly, the equations of motion of honeycomb structures for in-plane motion are given by,

M üx0 = i=1,3,5 (ξ n ∆n i e x • e i + ξ s ∆s i e x • d i ), (2.10a) M üy0 = i=1,3,5 (ξ n ∆n i e y • e i + ξ s ∆s i e y • d i ), (2.10b
)

I φ0 = R i=1,3,5 (ξ s ∆s i + ξ b ∆b i ), (2.10c
)

M üx1 = i=0,2,4 (ξ n ∆n i e x • e i + ξ s ∆s i e x • d i ), (2.10d
)

M üy1 = i=0,2,4 (ξ n ∆n i e y • e i + ξ s ∆s i e y • d i ), (2.10e
)

I φ1 = R i=0,2,4 (ξ s ∆s i + ξ b ∆b i ). (2.10f)
The inner produces of e x,y with e i and d i in Eqs. (2.9) and (2.10) can also be calculated.

Table 2.2 shows the results for the honeycomb structures. i-th 1 3 5

e x • e i 1 -1/2 -1/2 e y • e i 0 √ 3/2 - √ 3/2 e x • d i 0 - √ 3/2 √ 3/2 e y • d i 1 -1/2 -1/2 (b) Particle "1" i-th 0 2 4 e x • e i -1 1/2 1/2 e y • e i 0 - √ 3/2 √ 3/2 e x • d i 0 √ 3/2 - √ 3/2 e y • d i -1 1/2 1/2

Hexagonal monolayer membrane with out-of-plane motion

In this section, the out-of-plane motion is considered in the hexagonal granular membrane. The bulk wave properties are discussed.

Dispersion curves and complete band gaps

To analyze the out-of-plane motion in hexagonal membrane, the equations of motion are given in Eqs. (2.7). Suppose the solution of the equations of motion has a plane wave form,

v(t, x, y) = v out e iωt-ikxx-ikyy , (2.11)
where v out = [u; Φ; Ψ] is a vector of the motion component with Φ = Rϕ and Ψ = Rφ (R the radius of the particle), ω is the cyclic frequency, k x and k y are the wave vectors along x-and y-directions, respectively. According to Eqs. (2.7), the dynamical equation can be obtained,

S out v out = Ω 2 v out , (2.12)
where Ω = ω M/ξ s is the normalized frequency, M is the mass of beads, and S out is a 3 × 3 dynamical matrix, whose elements are of the following forms,

S 1,1 out = sin 2 2q x + sin 2 2q + + sin 2 2q -, (2.13a
)

S 1,2 out = √ 3i 4 (sin 2q + -sin 2q -), (2.13b
)

S 1,3 out = - i 4 (2 sin 4q x + sin 2q + + sin 2q -), (2.13c) S 2,2 out = 3P 2 - 3 -3η b -η t 4 P(sin 2 2q + + sin 2 2q -) + η t P sin 2 2q x , (2.13d)
S 2,3 out = √ 3P 4 (1 -η b + η t )(cos 2 2q --cos 2 2q + ), (2.13e)
S 3,3 out = 3P 2 - 1 -η b -3η t 4 P(sin 2 2q + + sin 2 2q -) -(1 -η b )P sin 2 2q x , (2.13f) S 2,1 out = -PS 1,2 out , (2.13g)
S 3,1 out = -PS 1,3 out , (2.13h)
S 3,2 out = S 2,3 out . (2.13i)
Above, q x = k x a/4 and q y = √ 3k y a/4 are the normalized wave vectors, and q ± = q x ± q y . P = M R 2 /I with I the moment of inertia of particle. n b,t = ξ b,t /ξ s is the bending/torsional rigidity normalized to the shear one. The solution of this eigenvalue problem gives the Ωk dispersion curves, which highly depends on the values of η b and η t .

In Fig. 2.6, the Ωk dispersion curves for different values of the ratios η b and η t 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. 2.6 we show the dispersion curves, following the path KΓMK of the first Brillouin zone (see Fig. 2.4), for the case of a fixed η b = 0.1 and varying the η t . Each of the eigenmodes consists of three components, the displacement u, and the rotations Φ, Ψ. 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 Γ, M and K are decoupled into 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 (ΓM segment), the eigenfrequency of the pure displacement u mode at M point is fixed at Ω u = √ 2, while there is also a pure rotational Φ mode at Ω Φ = 5(3η b + η t )/2. Since the latter depends on η b and η t , a forbidden band gap along this direction can exist and be tuned.

In particular, as long as 3η b + η t 8/5, this band gap exists at 5

(3η b + η t )/2 Ω √ 2 (see

Zero-frequency modes

For the extreme case where the impacts of bending and torsional forces are ignored (η b = η t = 0), for example, assuming (δ/R) 2 ≪ 1, the dynamical matrix S out reduces to, The dispersion curve is presented in Fig. 2.7(a). In the case of waves propagating along the y-direction, a band of zero-frequency modes is noted along the ΓM (q x = 0, 0 < q y < π/2) direction. The dynamical matrix S out along this direction is reduced further to the form,

       3 2 -cos 2q x cos 2q y - 1 2 cos4q x √ 3i 2 cos
S out =        2 sin 2 q y √ 3i 2 sin 2q y 0 - √ 3iP 2 sin 2q y 3P 2 cos 2 q y 0 0 0 P 2 (cos 2 q y + 2)        . (2.15)
Eq. (2.15) predicts a pure Ψ mode at Ω = 5(cos 2 q y + 2)/2 and two u+Φ coupled modes at Ω = 0 and Ω = 15 -7 sin 2 q y /2. Thus, the zero-frequency modes are u+Φ modes, i.e., coupled translational/ rotational modes with zero group velocity (non-propagative modes). Although it has been predicted before in many systems, the physical explanation of the existence of the zero-frequency mode in granular crystal has seldom been discussed. Γ 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 ΓM, 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 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 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. 2.7(a) to Fig. 2.6). It should be mentioned that the bending/torsional interactions, as shown in the following, can also initiate the propagation of modes with non-zero frequency which are non-propagating in their absence.

Degenerate modes and Dirac cone

Due to symmetry [START_REF] Kogan | Symmetry classification of energy bands in graphene[END_REF][START_REF] Tang | Symmetry and lattice dynamics, graphene simulation[END_REF], 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 acoustic effects [START_REF] Sigalas | Elastic and acoustic wave band structure[END_REF][START_REF] Kushwaha | Acoustic band structure of periodic elastic composites[END_REF][START_REF] Phani | Wave propagation in two-dimensional periodic lattices[END_REF][START_REF] Lu | Negative birefraction of acoustic waves in a sonic crystal[END_REF][START_REF] Sukhovich | Experimental and theoretical evidence for subwavelength imaging in phononic crystals[END_REF][START_REF] Zhang | Broadband acoustic cloak for ultrasound waves[END_REF][START_REF] Stenger | Experiments on elastic cloaking in thin plates[END_REF][START_REF] Zheng | Acoustic cloaking by a near-zero-index phononic crystal[END_REF][START_REF] Huang | Dirac cones induced by accidental degeneracy in photonic crystals and zero-refractive-index materials[END_REF][START_REF] Inserra | A method of controlling wave propagation in initially spatially periodic media[END_REF]. At K point, namely by setting q x = π/6 and q y = π/2, the dynamical matrix in Eqs.

Hexagonal monolayer membrane with in-plane motion

In this section, the in-plane motion is considered in the hexagonal granular membrane. The bulk wave dispersion, zero-frequency modes, degenerate modes are discussed.

Dispersion curves and complete band gaps

By analysing the in-plane motion in hexagonal membrane, the equations of motion are given in Eqs. (2.8). We still consider the plane wave solution for the equations of motion,

v(t, x, y) = v in e iωt-ikxx-ikyy , (2.16)
where v in = [u x ; u y ; Φ] is a vector of the moment component. Substituting into Eqs. (2.8), the dynamical equation for the in-plane motion can be found,

S in v in = Ω 2 v in , (2.17)
where v in = [u x ; u y ; Φ], and S in is a 3 × 3 dynamical matrix, whose elements are of the following forms,

S 1,1 in = η n + 3 4 (1 -cos 2q x cos 2q y ) + η n 2 (1 -cos 4q x ), (2.18a) S 1,2 in = √ 3(η n -1) 4 sin 2q x sin 2q y , (2.18b
)

S 1,3 in = - √ 3i
S 2,1 in = S 1,2 in , (2.18g)
S 3,1 in = -PS 1,3 in , (2.18h) S 3,2 in = -PS 2,3 in . (2.18i)
Above, η n = ξ n /ξ s is the normal rigidity normalized to the shear one. Solving the eigenvalue problem of Eq. (2.18), the Ωk dispersion curves for the in-plane motion can be obtained.

To demonstrate the influence of rotations, Fig. 2.9 presents the dispersion curves when ξ n = 1 while the value of ξ b is increasing from 0 to 1.5. For the in-plane motion, the eigenmodes contain Note that, as shown in Fig. 2.9, 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.

Degenerate modes and Dirac cone

The Dirac cone can be found at the K point, noted by green arrows in Fig. 2

Honeycomb monolayer membrane with out-of-plane motion

In this section, the out-of-plane motion is considered in the honeycomb granular membrane.

The bulk wave dispersion, zero-group-velocity, zero-frequency modes and degenerate points are discussed.

Dispersion curves

The honeycomb monolayer membrane and its corresponding first BZ are shown in Fig. 2.5. The calculation of the dispersion curves is similar to the case of the hexagonal membrane with out-ofplane motion, but now the unit cell contains two sublattices A and B. According to Eqs. (2.9), one could finally get the dynamical equation by looking for plane wave solutions:

S ′ out v ′ out = Ω 2 v ′ out , (2.19)
with the eigenvector

v ′ out = [u A ; Φ A ; Ψ A ; u B ; Φ B ; Ψ B ]. S ′
out is a 6 × 6 dynamical matrix, which is of the following form,

S ′ out = Q G N Q , (2.20)
where Q, G and N are 3 × 3 matrices of the forms,

Q =      3 0 0 0 3P(1 + η b + η t ) 2 0 0 0 3P(1 + η b + η t ) 2      , (2.21a)
G =       -ε -2ζ √ 3i̺ ε -ζ - √ 3iP̺ 3 -3η b -η t 2 Pζ -η t Pε - √ 3(1 -η b + η t ) 2 iP̺ -P(ε -ζ) - √ 3(1 -η b + η t ) 2 iP̺ 1 -η b -3η t 2 Pζ + (1 -η b )Pε       , (2.21b)
N =       -ε * -2ζ * √ 3i̺ * -ε * + ζ * - √ 3iP̺ * 3 -3η b -η t 2 Pζ * -η t Pε * √ 3(1 -η b + η t ) 2 iP̺ * P(ε * -ζ * ) √ 3(1 -η b + η t ) 2 iP̺ * 1 -η b -3η t 2 Pζ * + (1 -η b )Pε *       , (2.21c)
where ε = e -2iqx , ζ = e iqx cos q y and ̺ = e iqx sin q y . The dispersion curves when η t = 0 while η b is changing are shown in Figs. and (f). Nevertheless, when η b = 0, a flat band, which supports the zero-frequency mode having zero group velocity as well, could be noticed.

Zero-group-velocity and zero-frequency modes

We note that two flat bands with zero-group-velocity appear in Figs. Recently, it was reported that the dispersion in the vicinity of accidentally triple degeneracy at Γ point could lead to the zero-refractive-index in the photon and phonon systems [START_REF] Kushwaha | Acoustic band structure of periodic elastic composites[END_REF]. 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 [START_REF] Phani | Wave propagation in two-dimensional periodic lattices[END_REF]. Our theoretical results presented above indicate that the granular membranes, under conditions of some particular modifications, could support similar phenomena.

Honeycomb monolayer membrane with in-plane motion

In this section, the in-plane motion is taken into account in the honeycomb granular membrane.

The bulk wave properties, such as dispersion curves, zero-group-velocity and zero-frequency modes and degenerate points, are investigated.

Dispersion curves

The structure of the membrane and the corresponding first BZ is shown in Fig. 2.5. 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. According to Eqs. (2.10), the dynamical equation has the following form:

S ′ in v ′ in = Ω 2 v ′ in , (2.22)
with the eigenvector v ′ in = [u xA ; u yA ; Φ A ; u xB ; u yB ; Φ B ]. S ′ in is a 6 × 6 dynamical matrix, which takes the following form,

S ′ out = Q ′ G ′ N ′ Q ′ , (2.23)
where Q ′ , G ′ and N ′ are 3 × 3 matrices of the forms,

Q ′ =      3(1 + η n ) 2 0 0 0 3(1 + η n ) 2 0 0 0 3P(1 + η b )      , (2.24a)
G ′ =       - 3 + η n 2 ζ -η n ε - √ 3(η n -1) 2 i̺ - √ 3i̺ - √ 3(η n -1) 2 i̺ -( 1 + 3η n 2 ζ + ε) -(ε -ζ) √ 3iP̺ P(ε -ζ) (1 -η b )P(2ζ + ε)       , (2.24b
)

N ′ =       - 3 + η n 2 ζ * -η n ε * √ 3(η n -1) 2 i̺ * - √ 3i̺ * √ 3(η n -1) 2 i̺ * -( 1 + 3η n 2 ζ * + ε * ) (ε * -ζ * ) √ 3iP̺ * -P(ε * -ζ * ) (1 -η b )P(2ζ * + ε * )       . ( 2

Zero-frequency modes and degenerate modes

As shown in Figs. 2.15(a), (b) and (c), the zero-frequency mode for in-plane 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 (ΓM). As we have described for the out-of-plane cases, motion of rotational type plays a major role for the formation of non-propagating 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. The absence of bending rigidity can still lead to zero-frequency mode, but only in ΓM 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. 2.9), a zero-frequency mode could be found only in the honeycomb lattice. This is expected because the beads in the hexagonal membrane have more constrains than those in the honeycomb lattice.

Conclusion

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 does not modify its 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 under which the considered membranes non-propagating modes at finite (non-zero) frequencies could also exist. These nonpropagating 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 zero-frequency 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 the elastic waves. In the next Chapter, we will discuss the propagation of edge wave in granular crystals.

Introduction

In the previous Chapter, we have discussed the dispersion properties of elastic waves in hexagonal and honeycomb monolayer granular membranes with either out-of-plane or in-plane particle motion. When rotational degrees of freedom (DOF) are taken into account, the bending and torsional rigidities of the intergrain contacts can control some of the phononic modes. The zerofrequency modes, zero-group-velocity modes and their transformation into slow propagating phononic modes were predicted. We also studied the formation and manipulation of Dirac cones and multiple degenerate modes.

In the last two decades, graphene, a monolayer of carbon atoms arranged in a 2D honeycomb lattice, has become one of the most fascinating objects in material sciences and condensed-matter physics [START_REF] Geim | The rise of graphene[END_REF][START_REF] Castro Neto | The electronic properties of graphene[END_REF][START_REF] Sarma | Electronic transport in twodimensional graphene[END_REF]. Owing to its specific crystal symmetry, unusual energy band structure for charge carriers can be formed. An interesting consequence of a rich electronic spectrum is that the edge states may exist when boundaries are present. The existence of boundary, on one hand, breaks the crystal symmetry on the edges, but on the other hand, preserves some properties of the bulk, leading to the edge states having distinct propagation properties from the bulk modes and hence playing important roles in transport of electronic excitations [START_REF] Zak | Symmetry criterion for surface states in solids[END_REF][START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF][START_REF] Xiao | Surface impedance and bulk band geometric phases in one-dimensional systems[END_REF]. For example, the zero modes (zero-frequency edge states) were theoretically and experimentally studied in graphene ribbons [START_REF] Fujita | Peculiar localized state at zigzag graphite edge[END_REF][START_REF] Nakada | Edge state in graphene ribbons: Nanometer size effect and edge shape dependence[END_REF][START_REF] Wakabayashi | Electronic and magnetic properties of nanographite ribbons[END_REF][START_REF] Ryu | Topological origin of zero-energy edge states in particle-hole symmetric systems[END_REF][START_REF] Ryu | Zero-energy edge states and chiral symmetry breaking at edges of graphite sheets[END_REF][START_REF] Kohmoto | Zero modes and edge states of the honeycomb lattice[END_REF][START_REF] Esaki | Zero modes, energy gap, and edge states of anisotropic honeycomb lattice in a magnetic field[END_REF][START_REF] Checkelsky | Zero-energy state in graphene in a high magnetic field[END_REF]. More interestingly, the quantum Hall (QH) effect and quantum spin Hall (QSH) effect, where the topologically protected gapless edge states are robust channels for electronic transport, have been demonstrated [START_REF] Hatsugai | Chern number and edge states in the integer quantum Hall effect[END_REF][START_REF] Zhang | Experimental observation of the quantum Hall effect and Berry's phase in graphene[END_REF][START_REF] Novoselov | Two-dimensional gas of massless Dirac fermions in graphene[END_REF][START_REF] Kane | Z 2 topological order and the quantum spin Hall effect[END_REF].

Analogous to the above electron systems, the existence of edge state can emerge in completely different physical systems [START_REF] Wang | Coriolis force induced topological order for classical mechanical vibrations[END_REF][START_REF] Ni | Topologically protected one-way edge mode in networks of acoustic resonators with circulating air flow[END_REF][START_REF] Yang | Topological acoustics[END_REF][START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF][START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF][START_REF] Zhong | Acoustic analog of monolayer graphene and edge states[END_REF][START_REF] Rechtsman | Topological creation and destruction of edge states in photonic graphene[END_REF][START_REF] Bellec | Manipulation of edge states in microwave artificial graphene[END_REF][START_REF] Ochiai | Photonic analog of graphene model and its extension: Dirac cone, symmetry, and edge states[END_REF][START_REF] Ao | One-way edge mode in a magneto-optical honeycomb photonic crystal[END_REF][START_REF] Torrent | Acoustic analogue of graphene: observation of Dirac cones in acoustic surface waves[END_REF][START_REF] Kalesaki | Dirac cones, topological edge states, and nontrivial flat bands in two-dimensional semiconductors with a honeycomb nanogeometry[END_REF][START_REF] Savin | Vibrational Tamm states at the edges of graphene nanoribbons[END_REF]. In the photonic/phononic analog of graphene, novel phenomena, such as zero-energy mode on the edge [START_REF] Bellec | Manipulation of edge states in microwave artificial graphene[END_REF] and topological one-way edge states [START_REF] Ni | Topologically protected one-way edge mode in networks of acoustic resonators with circulating air flow[END_REF][START_REF] Hatsugai | Chern number and edge states in the integer quantum Hall effect[END_REF][START_REF] Ochiai | Photonic analog of graphene model and its extension: Dirac cone, symmetry, and edge states[END_REF][START_REF] Ao | One-way edge mode in a magneto-optical honeycomb photonic crystal[END_REF], have also been reported. However, in granular crystals, where most of the previous work addressed the structures of infinite size, only a few work have considered the influence of a boundary. For instance, in 1D semi-infinite monatomic granular chains, the existence of localized modes have been theoretically demonstrated [START_REF] Pichard | Localized transverserotational modes in linear chains of equal masses[END_REF]. In 3D cubic granular crystals, the Rayleigh-type surface acoustic waves (SAWs) and shear-horizontal-type SAWs have been theoretically analyzed, and the existence of ZGV SAWs of Rayleigh-type has been predicted [START_REF] Pichard | Surface waves in granular phononic crystals[END_REF]. Motivated by these recent results, the study of SAWs (or edge states) in 2D granular graphenes shows an interest as a necessary step towards their applications for elastic wave control, e.g. through the design of artificially modified granular graphenes (granular meta-graphenes) with non-trivial topological properties.

In this Chapter, we theoretically study the existence of edge states at the free boundaries of a semi-infinite mechanical granular graphene (GG), which is a 2D homogeneous monolayer of elastic beads packed into a honeycomb lattice of lattice constant a = 2 √ 3R as shown in Fig. 3.1.

An out-of-plane motion is imposed to the considered granular graphene with the zigzag (blue beads in Fig. 3.1) and armchair (orange beads in Fig. 3.1) boundaries. Consequently, each individual bead possesses one out-of-plane translational and two in-plane rotational DOFs, i.e. Fig. 2.1(a), and the mechanical interactions between the beads are provided by linear shear, bending and torsional rigidities, as depicted in Fig. 2.2(a)-(c), respectively. The dispersion curves of edge modes for varying bending/torsional rigidities are studied. Similar to the previous conclusions for the bulk mode [START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Kane | Topological band theory and Z 2 invariant[END_REF], due to the rotational DOFs, the dispersion branches of edge modes can be modified by tuning the bending and torsional rigidities. On the other hand, we demonstrate that nearly zero-frequency modes may appear on both the zigzag and armchair edges when bending/torsional rigidities are weak. This result is different from the cases in electronic/photonic systems where zero-frequency modes do not exist on the armchair edge when the sublattices A and B are the same [START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF][START_REF] Gusev | How acoustic waves are guided in buried subsurface channels in unconsolidated granular media[END_REF][START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: properties of rotational modes[END_REF][START_REF] Pavlov | A 2D granular medium with rotating particles[END_REF]. The physical origin of zero-frequency modes in the 2D granular monolayer honeycomb crystals has been interpreted in Ref. [START_REF] Zheng | Zero-frequency and slow elastic modes in phononic monolayer granular membranes[END_REF], where the zero-frequency mode results from the zero bending/torsional rigidities. In this Chapter, we

show that the existence of zero-frequency edge modes is also intimately related to the weakness of bending/torsional intergrain couplings.

The organization of this Chapter is as follows. In Sec. 3.2 the mechanical granular graphene is presented and the analytical model is described. The boundary conditions on the zigzag edge are introduced in Sec. 3.3, where the spectra of edge states and the zero-frequency modes are also evaluated, presented and discussed. For the armchair type of boundary, the boundary conditions, the existence of edge states and zero-frequency edge modes are systematically examined in Sec. 3.4. Finally, we present the conclusions of this work in Sec. 3.5.

Elastic structure and energy band projection

In this section, we show the structure of a granular graphene and discuss another band diagram, which is useful for the edge waves.

Granular graphene

The mechanical granular graphene presented in Fig. 3.1 exhibits two types of edges, zigzag and armchair, indicated in Fig. 3.1 by beads labeled in blue and orange, respectively. As discussed in Sec. 2.2, for the out-of-plane motion, each individual bead in the granular graphene exhibits the out-of-plane displacement (u) and the in-plane rotational angles ϕ and φ (Fig. 2.1(a)).

The dynamics and the couplings of these mechanical motions are controlled by the following forces and/or moments: (1) Shear forces (Fig. 2.2(a)). (2) Torsional moments (Fig. 2.

2(b)). (3)

Bending moments (Fig. 2.2(c)). Here, we would like to emphasise that the granular graphene is semi-infinite exhibiting a zigzag or an armchair edge configuration. These two boundaries are analyzed below separately.

To begin, we label the sublattice A (B) in a normalized coordinate (m, n) by A m,n (B m,n ), as shown in Fig. 3.1, where m and n are both integers. For the unit cell on site (m, n), the dynamical equation of bulk mode in Eq. ( 2. [START_REF] Pichard | Two-dimensional discrete granular phononic crystal for shear wave control[END_REF]) can be rewritten in the matrix form,

Q ′ G ′ N ′ Q ′ v A v B = Ω 2 v A v B , (3.1
v A = Q ′-1 (-G ′ )v B = S A v B , (3.2a)
(N ′ S A + Q ′ )v B = S B v B = 0, (3.2b)
where S A and S B are 3 × 3 matrices. Eqs. (3.2) suggest that, v A and v B are not independent.

Knowing one, the other can be simultaneously obtained from Eqs. (3.2). Let us consider that the displacement and rotations of particles A and B take the plane wave forms of solution,

v A = [u A ; Φ A ; Ψ A ] = K[α A ; β A ; γ A ]e i(ωt-mqx-nqy) , (3.3a) v B = [u B ; Φ B ; Ψ B ] = K[α B ; β B ; γ B ]e i(ωt-mqx-nqy) . (3.3b)
where K is a constant, α, β and γ are the coefficients of the motion components u, Φ and Ψ, respectively. The coefficients and their corresponding motion components satisfy the following relation,

u α = Φ β = Ψ γ = K. (3.4)
As long as the coefficients are given, the three motion components of an eigenmode can be settled. Since v A and v B are coupled by the interaction of the sublattices, we can substitute v B into Eq. (3.2), then the relations between those coefficients can be obtained,

α B = S (3,3) B S (2,2) B -S (3,2) B S (2,3) B , (3.5a
)

β B = S (3,1) B S (2,3) B -S (3,3) B S (2,1) B , (3.5b
)

γ B = S (3,2) B S (2,1) B -S (3,1) B S (2,2) B , (3.5c)
α A ; β A ; γ A = S A • α B ; β B ; γ B . (3.5d)
In Eqs. (3.5), the coefficients are given by S A and S B , which depend only on the normalized frequency Ω and wave vectors q x and q y . Consequently, an eigenmode of the system can be determined if Ω, q x and q y are known.

However, one should notice that, Ω is not independent from q x and q y . Their relation can be obtained from the dynamical equation shown in Eq. (2.19). For non-trivial solutions, Eq. (2.19) leads to the following condition,

det |S ′ out -Ω 2 I 6 | = 0. (3.6)
where I 6 is a 6 × 6 identity matrix. For the three variables Ω, q x and q y in Eq. (3.6), if the two of them are known, then the third one can be obtained. Note that, in this work, we analyze the granular system in which there is no loss or gain of energy, therefore Ω should be pure real.

Then the values of q x and q y determine the modes in the granular graphene to be propagative or evanescent.

For the propagating bulk modes, q x and q y are both real, leading to the pass bands for elastic wave in the structure. However, when modes are not in the pass band, i.e. in the band gap, the complex wave number reveals the fact that the amplitudes of modes are decaying when they penetrate into the structure. If the system has boundaries, there are possibilities for the existence of some special eigemodes that are localized near the boundaries [START_REF] Pichard | Localized transverserotational modes in linear chains of equal masses[END_REF]. Those modes support the edge waves which are propagating along the boundary but decaying into the bulk.

For example, suppose that an edge mode is propagating along the edge in y-direction, while decaying into the bulk. In this case, q y should be real, but q x must contain imaginary part so that edge wave is evanescent in the direction orthogonal to the surface. We will further discuss the edge solution in Sec. 3.3. those in the non-symmetry directions are taken into consideration, all the possible Ωk curves can be obtained and their combination in the k space would become the dispersion surfaces as shown in Fig. 3.2(c) by blue and gray surfaces. By displaying the eigenfrequency in the direction perpendicular to q x and q y plane, the three-dimensional (3D) band diagram is constructed, e.g.

Projection of bulk modes on a particular direction

Fig. 3.2(c).

In the 3D band diagram, all the propagating bulk modes which possess real q x and q y are included in the dispersion surfaces. However, actually the 3D band diagram is not very popular when the study is aiming at bulk wave dynamics of the system. Instead, the 2D band structure is more widely used since it is sufficient to predict the allowed bands and the forbidden bands for the bulk modes.

When studying the edge waves in the system, the 3D band diagram becomes a necessary step to obtain the projected band diagram, which is very efficient in determination of the regions in the Ωk space where the propagating edge waves could potentially exist. In the projected band structure, all the propagating modes in the 3D band diagram are projected into a particular direction along which the existence of the edge waves is studied. For example, suppose that the boundary we are interested in is the zigzag edge in a granular graphene as shown in Fig. 3.1.

In this case, the edge waves are propagating along the zigzag edge in y-direction, but their amplitudes are decaying when waves penetrate into the bulk. Physically, it implies that the wave vector q y of edge waves is real but q x is complex in such a way that edge waves can propagate along the edge, meanwhile being localized near the boundary. The projected band structure of the zigzag edge is obtained by projecting the modes in the 3D band diagram onto the q y direction as presented in Fig. 3.2(d). After the projection, the dispersion surfaces in the 3D band diagram are collapsed onto the blue and gray regions, in which the bulk modes could transport the energy from the edge into the bulk. The powerful advantages of the projected band structure can be exhibited through the scenario that all the possible edge modes must lie in the blank spaces of the projected band diagram. As depicted in Fig. 3.2(d), it reveals the dispersion relations between Ω and q y . For any point in the projected band structure, it corresponds to a real Ω and a real q y , indicating potential solutions for waves propagating along y-direction. Since the projected band (blue and gray areas) implies the existence of propagating bulk modes (q x is real), the existence of edge wave is only possible in the blank areas where q x is complex. When the boundary condition is also taken into account, the relation of Ω and q y can be derived, leading to the edge wave dispersion.

We can see that, although the band structure, e.g. and should be located between them. Otherwise edge wave will emit bulk modes and will not be evanescent [START_REF] Pichard | Localized transverserotational modes in linear chains of equal masses[END_REF][START_REF] Pichard | Surface waves in granular phononic crystals[END_REF].

Zigzag edge

In this Section, the waves dynamics on a free zigzag edge is discussed. The boundary condition, dispersion of the edge modes, and slow edge waves are analyzed.

Boundary condition

To get the edge wave dispersion, the boundary condition of the mechanically free zigzag edge of the granular graphene, as highlighted by blue beads in Fig. 3.1(a), has to be derived. Assuming the zigzag edge is introduced on site (0, n) by removing the beads A 0,n and those with m < 0, the boundary conditions are derived from the cancellation of the interactions between the removed grains A 0,n and the edge grains B 0,n . For example, on site (0, 0) the interactions between sublattice A 0,0 and B 0,0 should be zero, thus the boundary conditions are written as,

u A 0,0 -u B 0,0 + (Ψ A 0,0 + Ψ B 0,0 ) = 0, (3.7a) η t (Φ A 0,0 -Φ B 0,0 ) = 0, (3.7b) η b (Ψ A 0,0 -Ψ B 0,0 ) = 0. (3.7c)
where [u 0,0 ; Φ 0,0 ; Ψ 0,0 ] can be regarded as the motion components of the edge mode. If the relations between the edge wave components are given, the edge wave dispersion can be obtained by using Eqs. (3.7). For an edge solution with known Ω and q y , its evanescent properties into the bulk are determined by q x , which is governed by the bulk dispersion relation. As discussed in Sec. 3.2.1, when Ω and q y are given, possible values of q x can be calculated by Eq. (3.6).

Numerically, for fixed values Ω and q y , there are six solutions of q x . Considering that edge waves should be evanescent into the bulk, q x can not be real. This leads to q x having three solutions with negative imaginary part and other three with positive imaginary part. Under the zigzag edge configuration in Fig. 3.1, physically solutions with positive imaginary part can not exist in the system, otherwise the amplitudes of the corresponding eigenmodes are increasing into the bulk, which is definitely forbidden for edge waves. Therefore, for known Ω and q y , there are three corresponding evanescent eigenmodes in the bulk. They all contribute to the localization of the edge mode at the free boundary m = 0. Consequently, the solution of edge wave can be written in the following form which is the superpositions of the three evanescent bulk modes,

(u A ; Φ A ; Ψ A ) zig = 3 j=1 K j (α A j ; β A j ; γ A j )e i(ωt-mq xj -nqy) , (3.8a) (u B ; Φ B ; Ψ B ) zig = 3 j=1 K j (α B j ; β B j ; γ B j )e i(ωt-mq xj -nqy) , (3.8b)
where j is an index denoting different evanescent bulk modes, α j , β j and γ j are given by Eq. (3.5).

Substituting Eqs. (3.8) into the boundary conditions in Eqs. (3.7), we arrive at,

3 j=1 α A j -α B j + γ A j + γ B j K j = 0, (3.9a) 3 j=1 β A j -β B j K j = 0, (3.9b) 3 j=1 γ A j -γ B j K j = 0, (3.9c)
which can be rewritten in the following form, S zig (q x , q y )K zig = 0, (3.10) where

K zig = (K 1 , K 2 , K 3
) is the eigenvector of the boundary matrix, and S zig (q x , q y ) is the boundary condition matrix of zigzag edge, which is of the form,

S zig (q x , q y ) =     α A 1 -α B 1 + γ A 1 + γ B 1 α A 2 -α B 2 + γ A 2 + γ B 2 α A 3 -α B 3 + γ A 1 + γ B 3 β A 1 -β B 1 β A 2 -β B 2 β A 3 -β B 3 γ A 1 -γ B 1 γ A 2 -γ B 2 γ A 3 -γ B 3     . (3.11)
For the existence of the nontrivial solutions, the following equation must be satisfied, det |S zig (q x , q y )| = 0. (3.12)

Combining Eq. (3.6) with Eq. (3.12), q x can be eliminated, leading to the edge wave solution and Ω -q y dispersion curves of edge wave in the blank regions of the projected band diagram.

Edge modes

Following the theory above, the dispersion curves of the edge waves for different values of the ratios η b and η t are numerically derived in this part. Without loss of generality, in Fig. 3.

2
we present the projected band structures, where the bulk modes of the first BZ (in blue color)

have been projected to the q y direction, for the case of a fixed η t = 0.1 and a varying η b . As shown in Fig. 3.3, the bulk modes are concentrated mainly in two separated areas. The low frequency area, containing the three lowest bulk bands, is defined as lower-frequency phonon area (LPA) and the high frequency area, composed of three highest bulk bands, is noted as higher-frequency phonon area (HPA). There exists a large band gap between the LPA and HPA due to the fact that bending/torsional rigidities are much smaller than shear rigidities. For increasing bending/torsional rigidities, this band gap is squeezed. Dirac cones are observed at q y = ±π/3 in the LPA and HPA as shown in Fig. 3 the rotations should be predominantly induced by shear forces, Fig. 2.2(a), which have the maximal possible arm equal to R for creating the grain-rotating moment. However, this is not always the case because under some circumstances simultaneous displacements and rotations of the grains could keep the shear springs unstrained [START_REF] Bliokh | Quantum spin Hall effect of light[END_REF]. In these situations, the interactions between the beads through the bending and torsional moments, which are much weaker than the characteristic shear moments, could play an important role. As it is qualitatively illustrated in Fig. 2.2(b) and (c), the existence of both those moments is due to non-vanishing lateral size, δ = 0, of the contacts/links between the grains. The shear and normal forces distributed in the cross-section of the contact/link provide the resistance to rotation of the beads in opposite directions relative to the axis connecting their centers (Fig. 2.2(b)) and relative to the axis normal to the centers connecting axis (Fig. 2.2(c)), respectively. However, both shear and normal forces distributed across the contact have very small arms, ∼ δ ≪ R, for creating beadrotating moments. The torsional and bending moments induced by the distributions of shear and normal forces, respectively, appear to be much smaller than the nominal moment that could be created by the shear spring in Fig. 2.2(a). For the estimates, the interactions through torsional and bending moments could be characterized by the corresponding characteristic torsional and bending rigidities,

ξ t ∼ ξ s (δ/R) 2 and ξ b ∼ ξ n (δ/R) 2 (see Sec. 2.
2), respectively. Thus the waves, which propagation is due to torsional/bending moment-type interaction between the neighbor grains, will have the velocities ∼ δ ξ s,n /M ≪ R ξ s,n /M that are much slower than the velocities of classical acoustic waves, which propagation is due to shear/normal force-type interactions in granular crystals. Correspondingly, these waves caused by torsional/bending interactions between the grains, could be considered as extremely slow in comparison with the acoustic waves propagating in the material composing the grains. Note that, commonly normal rigidity is comparable with shear rigidity, i.e., ξ n ≥ ξ s , thus bending rigidity is comparable with torsional rigidity, ξ b ≥ ξ t .

As it follows from the discussion above, for the realization of extremely slow waves in mechanical graphene, the dimensions of the contacts between the spheres should be much smaller than the dimensions of the spheres. In non-cohesive granular crystals this condition could be achieved in the regime of weak pre-compression [START_REF] Allein | Transversal-rotational and zero group velocity modes in tunable magneto-granular phononic crystals[END_REF][START_REF] Allein | Tunable magneto-granular phononic crystals[END_REF][START_REF] Coste | Solitary waves in a chain of beads under Hertz contact[END_REF][START_REF] Nesterenko | Dynamics of heterogeneous materials[END_REF]. Otherwise the granular graphene could be potentially prepared by three-dimensional printing. If rod-type connections between the spheres are printed with dimensions much smaller than the sphere dimensions, then the effective torsional and bending rigidities of these artificially created contacts will be in general much weaker than their shear (normal) rigidities.

In the following, we discuss the influence of weak torsional or/and bending couplings on the modes propagating along the free zigzag edge being spatially localized in the vicinity of the edge.

Zero-frequency modes

In this part, we first discuss the case of weak torsional rigidity only, and later we turn to the influence of both weak bending and weak torsional rigidities under the condition that ξ b ≥ ξ t .

When torsional moments are weak, for example η t = 10 It suggests that the group velocity of edge states is extremely slow when torsional rigidity is weak.

The result presented in Fig. 3.5 confirms that elastic waves supported by the weak torsional rigidity could be of several orders of magnitude slower

Armchair edge

In this section, the boundary condition, dispersion relation and zero-frequency modes of on a free armchair edge are investigated.

Boundary condition

The mechanical granular graphene with the free armchair edge is illustrated in Fig. 3.1 highlighted in orange beads. Assuming the armchair edge is on site (m, 0) by removing all the beads with n < 0, the boundary conditions are derived from the absence of the interbead interactions between the cut layer (m, -1) and the edge layer (m, 0). For example, on site (0, 0), the interactions between the beads A 0,0 (B 0,0 ) and B -1,-1 (A 1,-1 ) should be zero, thus the boundary conditions are obtained, the Ω and q x are given, an edge solution is obtained. The evanescent nature of the edge mode into the bulk is determined by a complex q y , which can be evaluated from the bulk dispersion relation.

2(u B -1,-1 -u A 0,0 ) - √ 3(Φ B -1,-1 + Φ A 0,0 ) + (Ψ B -1,-1 + Ψ A 0,0 ) = 0, (3.13a) η t (Φ B -1,-1 -Φ A 0,0 ) + √ 3η t (Ψ B -1,-1 -Ψ A 0,0 ) = 0, (3.13b) - √ 3η b (Φ B -1,-1 -Φ A 0,0 ) + η b (Ψ B -1,-1 -Ψ A 0,0 ) = 0, (3.13c) 2(u A 1,-1 -u B 0,0 ) - √ 3(Φ A 1,-1 + Φ B 0,0 ) -(Ψ A 1,-1 + Ψ B 0,0 ) = 0, (3.13d) η t (Φ A 1,-1 -Φ B 0,0 ) - √ 3η t (Ψ A 1,-1 -Ψ B 0,0 ) = 0, (3.13e) √ 3η b (Φ A 1,-1 -Φ B 0,0 ) + η b (Ψ A 1,-1 -Ψ B 0,0 ) = 0. ( 3
When Ω and q x are fixed, one can use Eq. (3.6) to get the possible values of q y . For known Ω and q y , numerically 12 solutions of q y are expected. Considering that the edge wave should be in band gaps of the projected band structure, q y has six solutions with negative imaginary part and other six with positive imaginary part. Under the armchair edge configuration in Fig. 3.1, only those solutions with negative imaginary part have physical meaning in the system. Thus, for a solution of the armchair edge, there are six corresponding evanescent eigenmodes that contribute to the decay of amplitude when the edge mode penetrates into the bulk. Therefore, the solution of edge wave can be written in the form of the superpositions of the six evanescent bulk modes,

(u A , Φ A , Ψ A ) arm = 6 j=1 K j (α A j , β A j , γ A j )e i(ωt-mqx-nq yj ) , (3.14a)
(u B , Φ B , Ψ B ) arm = 6 j=1 K j (α B j , β B j , γ B j )e i(ωt-mqx-nq yj ) , (3.14b)
where j = 1, • • • , 6 denotes the six evanescent bulk modes. Substituting Eqs. 3.14 into Eqs. 3.13, the boundary conditions can be rewritten in the matrix form, S arm (q x , q y )K arm = 0, (

where K arm = [K 1 ; K 2 ; K 3 ; K 4 ; K 5 ; K 6 ] is the eigenvector of the boundary matrix and S arm is a 6 × 6 boundary condition matrix composed of the following elements,

S (1,j) arm = -2α A j - √ 3β A j + γ A j + (2α 2 j - √ 3β B j + γ B j )e i(qx+q yj ) , (3.16a)
S (2,j) arm = -β A j - √ 3γ A j + (β B j + √ 3γ B j )e i(qx+q yj ) , (3.16b)
S (3,j) arm = √ 3β A j -γ A j -( √ 3β B j -γ B j )e i(qx+q yj ) , (3.16c)
S (4,j) arm = (2α A j - √ 3β A j -γ A j )e i(-qx+q yj ) -2α 2 j - √ 3β B j -γ B j , (3.16d)
S (5,j) arm = (β A j - √ 3γ A j )e i(-qx+q yj ) -β B j + √ 3γ B j , (3.16e)
S (6,j) arm = ( √ 3β A j + γ A j )e i(-qx+q yj ) - √ 3β B j -γ B j , (3.16f)
For nontrivial solutions, the determinant of the boundary matrix must satisfy the following relation,

det|S arm (q x , q y )| = 0. (3.17)
The solutions for the edge mode of armchair edge are computed by combining Eq. (3.17) with Eq. (3.6).

Edge modes

In this part, the Ωk projected band structures for different values of the parameters η b and η t are derived. It is shown that it is possible to modify the dispersion diagram by a fine tuning of the contact rigidities. For example, the dispersion curves of the edge states for a fixed η b = 0.4 but changing η t are depicted in Fig. 3.7, where edge modes are labeled as red dash-dot curves and the blue areas represent the bulk modes. The particular number of existing edge bands depends on the value of η t . Explicitly, for modes in the LPA, Fig. 3.7(a) shows the dispersion curves for η b = 0.4 and η t = 0.1. Two edge bands appear between the second and third bulk bands. When the value of η t decreases to 0.05, one of these edge bands disappears, as shown in Fig. 3.7(b). Continuing to decrease η t to 0.03, one new edge band appears in each of the regions below the first bulk band as presented in Fig. 3.7(c). Reducing η t to 0.01, Fig. 3.7(d), one of the edge bands which existed between the first and second bulk bands vanishes. We can predict that, in the LPA, in the regions between the third and second bulk bands, between the second and the first bulk bands, and below the first bulk band there always exists at least one edge band each, while one additional edge band may appear in each of these regions depending on the bending and torsional rigidities. Similar phenomena of variation of the number of edge bands have also found in the HPA of the dispersion diagram. For example, when η t changes from 0.1

Zero-frequency modes

In this part, we first investigate the case of weak torsional coupling only, and later we discuss the influence of both weak bending and torsional couplings under the condition that ξ b ≥ ξ t .

For very small value of η t , the quasi-flat bands of edge modes at near zero frequency may appear in the projected dispersion diagram. Fig. 3.8 shows the dispersion curves of edge state for η b = 0.2 and η t = 10 -5 . There exist a near-zero-frequency band and a quasi-flat band of non-zero frequency just above it. We study the edge modes of q x ≈ π/4 (labeled by circles in Fig. 3.8(a)). The corresponding movements of beads are demonstrated in Fig. 3.8(b) and (c). It is clearly seen that the edge waves are propagating along x-axis being in-depth localized mostly on just a few layers near the edge. The insert in Fig. 3.8(a) shows the normalized group velocity (V g = v g • 2/(3ω 0 R)) of the first (near-zero-frequency) and the second (non-zerofrequency) quasi-flat modes as a function of η t . Note that the normalized group velocity in the first band is positive while of the second one is negative. Principally, both bands have very slow velocity when the torsional rigidity is weak. For the extreme case when torsional rigidity tends to zero, both bands have the tendency to be perfect flat bands at zero frequency supporting non-propagative edge modes. Similar to the results on zigzag edge, weak torsional rigidity can initiate the propagation of zero-frequency modes by transforming them into slow modes.

When both bending and torsional rigidity are weak, near-zero-frequency modes also can be found as shown in Fig. 3.9(a), where the dispersion curves of the LPA for η b = η t = 10 -5 are exhibited. Similar to the case of the zigzag edge, the small value of bending rigidity can lead to collapse of projected bulk bands into the region close to zero frequency, thus the associated edge states are confined to very narrow low frequency zone. In this case, the dispersion curves of edge modes are very close to the bulk bands, resulting in the edge states exhibiting weak localization in the vicinity of the edge. For example, we analyzed the edge waves with wave vector q x ≈ π/4 marked by the circles in Fig. 3.9(a). The corresponding movements of beads in those three bands are shown in Fig. 3.9(b)-(d). As explained in Sec.3.3.2, the closer the edge band is located to the bulk band, the weaker localization of the edge modes is. This prediction is valid in Fig. 3.9(b)-(d) where three edge states of q x ≈ π/4 are all weakly localized near the edge. Similar to the results on the zigzag edge, for zero bending rigidity, slow edge modes in the LPA disappear.

Conclusion

In this Chapter, we theoretically studied the existence of edge states with out-of-plane motion on the free zigzag and armchair boundaries of 2D semi-infinite mechanical granular graphenes.

Due to the additional rotational DOFs, the bending and torsional moments can be initiated between adjacent beads. By modifying the bending and torsional rigidities of contact, one could manipulate the dispersion of the edge states. Although the bending/torsional rigidities are weak in granular crystals, we demonstrate that they are crucial for the existence of slow propagating near zero-frequency modes on both zigzag and armchair edges. When torsional rigidity is weak, the extremely slow near-zero-frequency edge modes, whose group velocity depends on the tor-sional rigidity, exist on both zigzag and armchair edges. When bending rigidity also becomes weak in addition to torsional one, edge states in the LPA of the dispersion diagram are localized very close to zero frequency. Due to the fact that they are very close to the bulk modes in the dispersion diagram, edge modes in this case can be weakly localized near the edge. For weak torsional rigidity, the quasi-flat bands at near zero frequency transform into the perfect zero-frequency edge modes when torsional rigidity tends to zero, while for weak bending rigidity, those edge modes show the tendency of vanishing as bending rigidity is continuing to reduce.

The granular graphene membrane could be, in principle, prepared by three-dimensional printing.

In accordance with the theoretical estimations [START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF][START_REF] Haldane | Model for a quantum Hall effect without landau levels: condensedmatter realization of the parity anomaly[END_REF], the intersphere connections, with bending and torsional rigidities much weaker than shear one, could be prepared if the dimensions of the contacts are much smaller than the radius of the sphere. There are also possibilities to assemble a granular graphene supporting the propagation of slow elastic edge waves manually. The elastic spheres discussed in this manuscript could be in macro scale, which means that the particles can be displaced manually and the structure can be practically assembled nearly perfectly with almost zero disorder and missing contacts [START_REF] Paulose | Topological modes bound to dislocations in mechanical metamaterials[END_REF][START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF]. For example, in Ref. [START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF], when the external force (pre-compression) achieved through the magnetic interaction between the spheres is around 4N , the wave dynamics is linear and the magnitude of bending rigidity is around 5 orders of magnitude smaller than the shear rigidity, which is exactly the same order of magnitude that was assumed for some of the estimates in our work. Therefore, assembling magnetic beads nearly perfectly in the geometry of the granular graphene and achieving weak bending and torsional interactions is potentially possible. An adapted shaker attached to beads on the edge can be used as the source to initiate the different motions of the beads. Particle acceleration can be recorded by using small accelerometers glued to the side of several beads or by laser vibrometer measurements. The physical realizations of the described granular graphene are of interest in devices of signal processing using acoustic waves and in the design of acoustic delay lines. The investigation of edge modes in mechanical granular graphene could in turn provide considerable insights into their optical/acoustic graphene analogues. The study of the near-zero-frequency modes caused by the existence of weak rotational interactions is useful for the potential realistic design of granular crystals supporting controllable propagation of extremely slow edge waves.

In addition, it paves the way for the analysis of a variety of possible modifications of the above presented granular graphene with the goal of designing granular metamaterials supporting oneway propagating edge states. In the next Chapter, we will discuss the topological transport of edge wave in mechanical granular graphene.

Introduction

In the previous Chapter, we reported theoretically the description and analysis of edge elastic waves in a semi-infinite mechanical granular graphene structure. The dispersion curves of the edge waves were theoretically derived and numerically analyzed for various configurations of bead couplings, as well as the existence of edge states when the torsional or/and bending rigidities were weak/vanishing. Quasi-flat edge mode dispersion curves with near zero frequency were observed for both the zigzag and armchair edges. These results on edge waves in mechanical granular graphene structures with rotational degrees of freedom were the necessary preliminary step for the design of the granular meta-graphenes with topologically protected unidirectional edge states.

As discussed in Sec. 1.5, topological wave propagation has been observed in many different systems [115-143, 203, 204]. For mechanical systems, the topologically protected one-way edge modes can be realized by different designs. Examples include metamaterials with rotating gyroscopes [START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF][START_REF] Nasha | Topological mechanics of gyroscopic metamaterials[END_REF], mechanical structures subjected to Coriolis force caused by permanent rotation of the system [START_REF] Wang | Coriolis force induced topological order for classical mechanical vibrations[END_REF], elastic metaplate with effective spin-orbit coupling (SOC) [START_REF] Mousavi | Topologically protected elastic waves in phononic metamaterials[END_REF] and coupled pendula system exhibiting quantum spin Hall (QSH) effect [START_REF] Ssstrunk | Observation of phononic helical edge states in a mechanical topological insulator[END_REF]. In those proposals, the systems either require an external rotating bias to break time-reversal symmetry [START_REF] Peng | Experimental demonstration of anomalous Floquet topological insulator for sound[END_REF][START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF][START_REF] Nasha | Topological mechanics of gyroscopic metamaterials[END_REF] or need to be sophisticatedly design to induce Z 2 topological order [START_REF] Mousavi | Topologically protected elastic waves in phononic metamaterials[END_REF][START_REF] Ssstrunk | Observation of phononic helical edge states in a mechanical topological insulator[END_REF]. In addition, the role of rotational degrees of freedom of individual masses has seldom been evaluated [START_REF] Deymier | Torsional topology and fermion-like behavior of elastic waves in phononic structures[END_REF].

Due to the existence of non-central forces, which can initiate the rotation of individual mass particle, rotations become crucial in wave propagation in mechanical discrete systems [START_REF] Schwartz | Vibrational Modes in granular materials[END_REF][START_REF] Suiker | Comparison of wave propagation characteristics of the Cosserat continuum model and corresponding discrete lattice models[END_REF].

Especially in granular crystals, rotations can be easily excited by the transverse force between particles, leading to the propagation of rotational elastic waves. Therefore, rotation of individual particles should be taken into account when exploiting the topological properties of elastic wave in granular crystals.

To realize topologically protected transport of elastic waves in granular crystals, one can also design a granular system following the strategies either with time-reversal symmetry (TRS)

breaking or with TRS invariance, as discussed in Sec. 1.5. In mechanical systems, it has been reported that TRS can be broken by an external rotating bias, such as by rotating gyroscopes [START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF][START_REF] Nasha | Topological mechanics of gyroscopic metamaterials[END_REF] and by permanent rotation of the system [START_REF] Wang | Coriolis force induced topological order for classical mechanical vibrations[END_REF]. Nevertheless, the extension of these insights to the granular system remains a huge challenge, since the elastic particles themselves in granular crystals can exhibit rotations and out-of-plane displacements. So far, efficiently breaking TRS in granular system has not been reported. In this Chapter, we investigate the topological transport of rotational edge waves in the GG with TRS. As demonstrated in Sec. 1.5, one way to realize the topological wave propagation under TRS invariance is to construct a analogous system exhibiting the QSH effect [128, 130-132, 134, 136]. Following this idea, two degenerate spin states both having Dirac dispersion in the absence of SOC are required. In addition, a complete band gap supporting the propagation of topological edge waves is necessary in the presence of SOC. The existence of the Dirac degeneracy or band touching in the dispersion diagram is very crucial for the topological transition. Before and after the band touching, the corresponding phase of the two bands could be generally different in topology, resulting in nontrivial topological numbers for the bands [START_REF] Ochiai | broken symmetry and topology in photonic analog of graphene[END_REF]. Therefore, there are two important ingredients for the existence of topological transport with TRS: (1) A degenerate Dirac dispersion or a double Dirac cone when the system is unperturbed. (2) A topological band gap when topological transition is induced to the system. This Chapter is organized as follows: In Sec. 4.2, a quadruple degeneracy of the modes and wave dynamics near the center of Brillouin zone (BZ) are studied. In Sec. 4.3, a Dirac equation is derived to describe the wave dynamics when the degeneracy is achieved near the Γ point. The concepts of the effective spin and helicity are illustrated in the mechanical granular system. In Sec. 4.4, we show that the GG, in which wave propagation is described approximately by the Dirac equation, exhibits non-trivial topological properties and the quasi-topological transport of rotational edge can be achieved. In Sec. 4.5, the wave dynamics of edge wave on an interface is theoretically demonstrated. The physical explanation for the quasi-topological transport is suggested. Finally, we summarize the studies of this chapter in Sec. 4.6.

Dirac dispersion and wave dynamics

In this section, we focus on the GG that exhibits degenerate Dirac dispersion at the Γ point.

The wave dynamics around the Dirac point is studied.

Degenerate Dirac dispersion

The GG under consideration is shown in Fig. 3.1 (see Sec. 3.2.1). Considering the out-of-plane motion, the dispersion curves of the GG can be obtained by solving the eigenvalue problem of the dynamical equation in Eq. (2.12). By controlling the normalized bending and torsional rigidities, η b and η t , the Ω -k dispersion relation can be tuned. However, a degenerate Dirac cone is essential to exhibit non-trivial topological properties in the GG. As predicted in Sec. 2.5.3, when the condition η b +η t = 1 is satisfied, a double Dirac cone can be observed at the Γ point of the BZ.

Consequently, in this Chapter we discuss the GG under the condition η b + η t = 1. In addition, topological transition has to take place in the GG. Basically, this can be accomplished by adding a perturbation to the GG. In this Chapter, we only analyze the case when a perturbation η ′ is added to the normalized bending and torsional rigidities, which will lead to the topological transition in the GG (see Sec. 4.3.2). Thus, in the GG with the perturbation η ′ and η b + η t = 1, the dynamical equation in Eq. (2.12) is rewritten as,

Dv = Ω 2 v, (4.1)
with D of a 6 × 6 dynamical matrix taking the form,

D =             3 0 0 -α √ 3iγ -β 0 Ω 2 D + 3Pη ′ 0 - √ 3iPγ Ω 2 D 6 β -Pη ′ α - √ 3 Ω 2 D 6 iγ 0 0 Ω 2 D + 3Pη ′ Pβ - √ 3
Ω 2 D 6 iγ - Ω 2 D 6 β -Pη ′ α -α * √ 3iγ * β * 3 0 0 - √ 3iPγ * Ω 2 D 6 β * -Pη ′ α * √ 3 Ω 2 D 6 iγ * 0 Ω 2 D + 3Pη ′ 0 -Pβ * √ 3 Ω 2 D 6 iγ * - Ω 2 D 6 β * -Pη ′ α * 0 0 Ω 2 D + 3Pη ′             , (4.2)
where α = 2e iqx cos q y +e -2iqx , β = e iqx cos q y -e -2iqx , γ = e iqx sin q y , and * denotes the complex conjugate. Ω D = √ 3P is the Dirac frequency. In the previous Chapters, we assumed that the elastic spheres in the GG are homogeneously filled. In this case, the value of P is fixed to be 2.5.

However, if we consider the spheres are hollow or with a custom mass distribution along their radius, then P is changed. The extreme case is that all the mass of the sphere is at the sphere periphery. In this case, the value of P is 1.5. To illustrate the influence of P on the dispersion, from now on we analyze the GG with P = 1.55, and η b = η t = 0.5. Due to the existence of non-central transverse force, rotational-translational coupled modes exist in the GG. However, when we consider the wave dynamics in the region close to the degenerate point, the four branches marked as red and green are dominant by rotations. The contributions of u components in the eigenmodes of the four branches actually are very weak. To show the small influence of u components in the four rotational-dominant branches, the eigenmodes at the Γ point, under the basis v = [u A ; Φ A ; Ψ A ; u B ; Φ B ; Ψ B ], can be obtained by setting q x = q y = 0 in Eq. (4.1), which leads to the following forms of the eigenmodes,

Wave dynamics around the Dirac point

            0 1 0 0 0 0             ,             0 0 1 0 0 0             ,             0 0 0 0 1 0             ,             0 0 0 0 0 1             . ( 4.3)
Eq. ( 4.3) shows that at the Γ point, the rotational modes decouple completely from the translational modes, forming pure Φ and Ψ modes. The contribution of u components to the eigenmodes at the degenerate point is zero. Around the small region of the degenerate point, for example, at the point 0.07K in the ΓK direction as shown in Fig. 4.2(b), the eigenmodes of the red and green branches are obtained numerically as,

            0.0881i 0.5332 -0.4560i -0.0881i -0.5332 -0.4560i             ,             0.0258 -0.4857i 0.5133 0.0258 -0.4857i -0.5133             ,             0.0892 -0.4425i 0.5442 -0.0892 0.4425i -0.5442             ,             0.0273i -0.5122 0.4868i 0.0273i 0.5122 0.4868i            
.

This shows that away from the degenerate point the u, Φ and Ψ components couple to each other.

However, we can calculate the energy stored in u components for each mode (u 2 /(u 2 +Φ 2 +Ψ 2 )):

1.55%, 0.13%, 1.59%, 0.15%. It shows that the contribution of u component for the bead motion is no more than 2% compared to the Φ and Ψ modes. Therefore, we ignore the small influence of u components when discussing the wave dynamics around the Dirac point. Then the 6 × 6 dynamical equation in Eq.(4.1) can be approximated by a 4 × 4 one,

D(k)v 0 = Ω 2 v 0 , (4.4)
where v 0 = [Φ A ; Ψ A ; Φ B ; Ψ B ], and D(k) is a 4 × 4 dynamical matrix of the following form,

D(k) =        Ω 2 D + 3Pη ′ 0 Ω 2 D 6 β -Pη ′ α - √ 3 Ω 2 D 6 iγ 0 Ω 2 D + 3Pη ′ - √ 3
Ω 2 D 6 iγ - Ω 2 D 6 β -Pη ′ α Ω 2 D 6 β * -Pη ′ α * √ 3 Ω 2 D 6 iγ * Ω 2 D + 3Pη ′ 0 √ 3 Ω 2 D 6 iγ * - Ω 2 D 6 β * -Pη ′ α * 0 Ω 2 D + 3Pη ′        . (4.5)
In the vicinity of the degenerate point, namely for a small wave vector ∆k = (∆k x , ∆k y ) and for ∆Ω ≪ Ω D , Eq.(4.4) can be approximated using the Taylor expansion at k = 0. By ignoring the second and higher order terms and Ω 2 = (Ω D + ∆Ω) 2 ≈ Ω 2 D + 2Ω D ∆Ω, Eq.(4.4) can be rewritten as,

D(∆k)v 0 = Ωv 0 , (4.6)
with

D(∆k) =        Ω D + m 0 iV D ∆k x -m -iV D ∆k y 0 Ω D + m -iV D ∆k y -iV D ∆k x -m -iV D ∆k x -m iV D ∆k y Ω D + m 0 iV D ∆k y iV D ∆k x -m 0 Ω D + m        , (4.7)
where . This verifies the fact that the quadruple degeneracy at the Γ point exhibits a liner dispersion relation. In the next section, we will show that the wave dynamics around the degenerate point can be described by Dirac equation.

V D = Ω D R/

Effective spin and helicity

In this section, we connect the wave dynamics around the degenerate point to a Dirac equation and we establish the effective spin and helicity in the granular graphene.

Dirac Hamiltonian and effective spin

The rotational wave behavior in the vicinity of the Γ point can be approximately described by Eq. (4.6). We can transform the original basis v 0 = [Φ A ; Ψ A ; Φ B ; Ψ B ] to a circular polarized basis ψ = U v 0 by the unitary matrix,

U = 1 √ 2        1 i 0 0 0 0 i 1 0 0 -1 -i i 1 0 0        . (4.11)
The circular polarized basis can be presented as ψ = [p; d], where p and d take the following forms,

p = p ↑ ip ↓ = 1 √ 2 Φ A + iΨ A i(Φ B -iΨ B ) ; d = d ↑ id ↓ = 1 √ 2 -(Φ B + iΨ B ) i(Φ A -iΨ A ) , (4.12)
with "↑↓" denoting the left/right circular polarizations. By the proposed transformation, Eq. As can be seen, the granular system exhibits a pair of mechanical "spin" states which physically are the rotations of particle with left/right circular polarizations (the rotation vector is circulating left/right handed with time). For a small ∆k around the Dirac point, any eigenstate of the system can be expressed as the linear combination of p ↑↓ (d ↑↓ ) modes for the p (d) branches. This is similar to the case in the electronic system where any spin state of the system is the linear combination of spin up and spin down states.

Effective spin-orbit coupling and helicity

When the perturbation is imposed to the system, namely η ′ = 0 , the original granular system through the unitary transformation,

U ′ = 1 √ 2        1 0 1 0 0 1 0 1 1 0 -1 0 0 1 0 -1        . (4.15)
After the transformation U ′ HU ′-1 , Eq.(4.14) can be written as,

H + (∆k) 0 0 H -(∆k) S + S - = δΩ S + S - . (4.16)
S ± = p ± d is the new spin basis that does not depend on the perturbation. The explicit form of S ± is given by,

S + = S ↑ + S ↓ + = 1 √ 2 p ↑ + d ↑ i(p ↓ + d ↓ ) , S -= S ↑ - S ↓ - = 1 √ 2 p ↑ -d ↑ i(p ↓ -d ↓ ) , (4.17)
Eq.(4.16) predicts two decoupled massive Dirac equations, H ± (∆k)S ± = δΩS ± , with the Hamiltonian,

H ± (∆k) = V D ∆σ • k ± mσ z . (4.18)
This implies that the perturbed system supports a pair of the spin-polarized modes S ± . The second term in Eq. (4.18) describes the hybridization of p and d spinor subsystems, which can be regarded as the effective spin-orbit coupling (SOC) in the GG [START_REF] Nalitov | Spin-orbit coupling and the optical spin Hall effect in photonic graphene[END_REF][START_REF] Wang | Elastic spin-Hall effect in mechanical graphene[END_REF] When ∆k = 0, Eq. (4.16) reduces to be ±mσ z S ± = δΩS ± . Its solutions are the eigenstates of the Pauli matrix σ z . When ∆k = 0, the Hamiltonian H ± predicts the helicity-locking property of the bulk modes [START_REF] Boechler | Discrete breathers in one-dimensional diatomic granular crystals[END_REF]. To demonstrate the helicity, Eq. (4.16) can be rewritten as,

±λ(σ • e)S ± = δΩS ± , (4.19)
where λ = V

Topological transport of rotatinal edge wave

So far, we have discussed the wave dynamics around the Dirac point. Several important concepts, such as spin, SOC and helicity, have been extended to the granular system. In this section, we study the topological property of the GG.

Topological number

In general, the dispersion diagram contains most of the useful information about the system such as the allowed and forbidden bands for waves propagation. However, the information on topology of the system is hidden in the eigenstates [START_REF] Kane | Topological band theory and Z 2 invariant[END_REF][START_REF] Dauphin | Efficient algorithm to compute the Berry conductivity[END_REF][START_REF] Xiao | Berry phase effects on electronic properties[END_REF]. As discussed in Sec. 4.3.2, due to the appearance of perturbation, effective SOC is induced to the granular system. The existence of the effective SOC leads to the coupling of the original spin states p and d, resulting in the new spin states S ± to be the eigenmodes of the system. This changing of the eigenmodes indicates that the GG could have experienced a critical phase transition in topology of the system. To illustrate the topological transition in the GG, we need to classify the topological order of each branch based on their eigenstates around the Dirac point [START_REF] Mei | Pseudo-time-reversal symmetry and topological edge states in two-dimensional acoustic crystals[END_REF]. We first analyse the cases when η ′ = 0, then come to the case when η ′ = 0.

When η ′ = 0, namely m = 0, the wave dynamics around Dirac point is described by Eq. (4.13), which can be reduced to a pair of uncoupled equations. Here, we take p modes as an example to analyse. Similar conclusions are valid for the d modes. The wave dynamical equation of the p modes around the Dirac point is,

0 V D (∆k x -i∆k y ) V D (∆k x + i∆k y ) 0 p ↑ p ↓ = ∆Ω p ↑ p ↓ . (4.22)
It predicts the linear dispersion relation: ∆Ω ± = ±V D ∆k. From Eq. (4.22), the following relations can be obtained,

V D (∆k x -i∆k y )p ↓ = ∆Ωp ↑ , (4.23a)
V D (∆k x + i∆k y )p ↑ = ∆Ωp ↓ . (4.23b)
The relations above lead to two eigenwavefunctions near the Dirac point,

p + (∆k) = 1 √ 2   1 ∆k x + i∆k y ∆k   e ikr , p -(∆k) = 1 √ 2   ∆k x -i∆k y -∆k 1   e ikr . ( 4

.24)

The subscript ± labels the modes in the upper (+)/lower (-) Dirac cones. According to Eq. (1.20), we can calculate the Berry phase, a key parameter identifying/characterizing the topology of material, around the Dirac cone [START_REF] Mei | Pseudo-time-reversal symmetry and topological edge states in two-dimensional acoustic crystals[END_REF][START_REF] Chen | Accidental degeneracy of double Dirac cones in a phononic crystal[END_REF][START_REF] Zak | Berry's phase for energy bands in solids[END_REF],

Γ ± = i C p ± (∆k)|∇ ∆k |p ± (∆k) d∆k, (4.25)
where C represents the close path of the BZ. Let's take p + as an example. Substituting it into Eq. (4.25), we find,

Γ + = i C p + (∆k)|∇ ∆k |p + (∆k) d∆k = i 2 C 1, ∆k x -i∆k y ∆k      0 e x + ie y ∆k   +   1 ∆k x + i∆k y ∆k   ir    d∆k = i 2 C 1 |∆k| 2 (∆k x d∆k x + ∆k y d∆k y ) + i 2 C 2ird∆k. (4.26)
Above, e x,y are the unit vectors along x and y directions, respectively. For a close path C, the second term in the right side of Eq. (4.26) is obviously zero. Then Eq. (4.26) becomes,

Γ + = i 2|∆k| 2 C (∆k x d∆k x + ∆k y d∆k y ) = C d∆k • A, (4.27)
From Eq. (4.27), the Berry connection is defined as A = i(∆k x e x + ∆k y e y )/2∆k 2 . However, according to Eq. (4.25), the Berry connection can also be defined as

A ′ = i p + (∆k)|∇ ∆k |p + (∆k) =
Ar. This exactly reveals the fact that Berry connection, similar to the electromagnetic vector potential, is gauge dependent [START_REF] Xiao | Berry phase effects on electronic properties[END_REF]. Consequently, it is useful to define, in analogy to electrodynamics, a gauge field tensor derived from the Berry vector potential as discussed in Sec. 1.6.2, When η ′ = 0, the effective SOC is induced in the system, the eigenmodes of the GG can be obtained from Eq. (4.19). Take the S + modes as an example, the corresponding Hamiltonian is given by

B = ∇ k × A. ( 4
H + = V D σ • ∆k + mσ z = λ ± e • σ.
The corresponding eigenstates with energies

λ ± = ± m 2 + V 2 D ∆k 2 are obtained, v + =    cos θ 2 e -iφ sin θ 2    , v -=    sin θ 2 e -iφ -cos θ 2    , (4.30)
where θ and ψ are the azimuthal angle and the polar angle of vector e, respectively. They satisfy the following relations,

cos θ = m λ ± , sin θ = V D |∆k| λ ± , cos φ = ∆k x |∆k| , sin φ = ∆k y |∆k| . (4.31)
We can see that actually the S + modes system governed by the Hamiltonian H + is similar to the two-level system described in Sec. 1.6.2. Let us consider the v -mode, the Berry connection is given by,

A θ = i v -|∂ θ v -= 0, (4.32a)
A φ = i v -|∂ φ v -= sin 2 θ 2 , (4.32b)
and the Berry curvature is,

B -= ∂ θ A φ -∂ φ A θ = 1 2 sin θ. (4.33a)
Similarly, the Berry curvature of the v + mode can be obtained,

B + = - 1 2 sin θ. (4.33b)
The Berry curvature of the S -states can be calculated the same way. 1. Φ modes: Ω Φ1 = √ 3P and Ω Φ2 = 3P(1 + 2η ′ ).

Boundary condition

Ψ modes: Ω

Ψ1 = √ 3P and Ω Ψ2 = 3P(1 + 2η ′ ).
Suppose the two GGs have the same mass of beads, it is clear that when η ′ = 0.1, the band width of the forbidden gap is √ 3P < Ω < 3P(1 + 2η ′ ) for the GG in the left side. However, when η ′ = -0.1, the forbidden gap is 3P(1 -2η ′ ) < Ω < √ 3P for the GG in the right side.

This shows that the two GGs with the same mass of beads but opposite signs of perturbation

M R = f M = M/(1 + g) (g
is the tuning factor). On the interface, the unit cell contains two particles A L and B R (the superscripts L, R represent the GGs in the left or right sides), and particle A L is connected with particle B R through unperturbed bonds (black). In general, the total forces/moments on particle A L of the interface can be described by the equation,

M Äe = F ′ inner + F neighor , (4.34a)
where M denotes the mass/moment of inertia of particle A L , A e is the movement (displacement/rotation) of particle A L on the interface, F ′ inner represents the interactions within the unit cell of interface, i.e., with the particle B R , and F neighbor corresponds to the interactions with the adjacent unit cells. For the propagating edge waves, the boundary motion in Eq. (4.34)(a) should be ensured by the evanescent modes of the bulk. Meanwhile, the evanescent bulk modes in the GG of the left side material satisfy the following equation for the particle A L , M Äl = F inner + F neighor , (4.34b)

            0 0 0 1 0 -1 0 η ′ 0 0 -(η + η ′ ) 0 0 0 η ′ 0 0 -(η + η ′ ) 1 0 1 0 0 0 0 η 0 0 0 0 0 0 η 0 0 0             v L =             0 0 0 1 0 -1 0 0 0 0 -η 0 0 0 0 0 0 -η 1 0 1 0 0 0 0 (η -η ′ ) 0 0 η ′ 0 0 0 (η -η ′ ) 0 0 η ′             v R , (4.36)
where

v L,R = [v A e iq L,R x ; v B e -iq L,R
x ] with v A,B the motion components of the particles A or B, η = ξ b,t /ξ s and η ′ = ξ ′ /ξ s with ξ ′ denoting the perturbation of stiffness. Basically, in numerical calculations, the edge wave solution can be obtained by combining the dynamical equation of the bulk modes, Eq. (4.1), and the boundary condition, Eq. (4.36). However, the analytical formulas are too cumbersome. In fact, we are only interested in the wave dynamics of the interface near the Γ point. In this case the solution can be simplified significantly, see Sec. 4.5.1

One way propagating edge waves

Numerically, the projected band structure of the GGs with the interface in Fig. 4.7 is presented in Fig. 4.9(a), where the grey lines correspond to the bulk modes which are projected to the q y = √ 3k y R direction. Using the boundary condition, the dispersion curves of edge modes can be obtained, as labelled by color lines in Fig. 4.9(a). There are two branches of edge waves in the topological gap labelled as cyan and magenta lines. The cyan (magenta) lines correspond to upward S u (downward S d ) propagating rotational edge waves. The green lines mark the edge branches that support the propagation of the topological trivial waves (see Sec. In fact, the cyan and magenta branches located in the topological bad gap support the propagation of edge wave with different polarizations. To clear see this, we chose two edge It should be mentioned that a topological transport on the interface exists only when the edge wave dispersion is gapless. However, the dispersion of the edge waves in the GG is not gapless as there is a tiny gap at the Γ point, Fig. 4.9(b). In the bulk, the system has C 6 symmetry, while due to the mismatch of the masses in the interface unit cell, the C 6 symmetry is broken, resulting in a tiny gap. The existence of the tiny gap actually leads to the coupling of spin-polarized helical edge states, inducing the hybridized S u and S d rotational edge waves which are not perfectly topologically protected. This is confirmed by the spatio-temporal evolution simulations in Fig. 4.9(f), where the interface configuration is identical to Fig. 4.9(e) except that a hole is placed close to the interface. We monitor the particle at the position marked by the circle to see if there is back reflection. The rotational movements of the particle as a function of time are shown in Fig. 4.9(g) both for the cases with/without the hole. In the case of immune propagation and no backscattering, one would expect identical dynamics. However, a small difference between these two is observed which could originate from a small amount of back-reflection due to the existence of tiny gap and coupling of the edge waves.

Wave dynamics on the interface

In the previous section, we showed that a pair of unidirectional edge waves were observed on the interface constructed by two GGs with opposite signs of perturbation. Since the two GGs had different masses, the C 6 symmetry is broken on the interface, leading to a tiny gap in the dispersion relation of the edge wave. In this section, we theoretically investigate the origin and influence of this tiny gap.

Edge waves on the interface

For a small tuning on the mass in the right side, we can treat this tuning as perturbation on the system. Consequently, we can evaluate bulk motion equations at Ω D at the Γ point by using Taylor expansion. Under the spin basis in Eq. (4.17), the dynamical equation is given by,

H e        S ↑ + S ↓ + S ↑ - S ↓ -        = ∆Ω        S ↑ + S ↓ + S ↑ - S ↓ -        , (4.37)
with H e of the form,

H e = (1 + g)        m L,R ± m V D (∆k x -i∆k y ) 0 0 V D (∆k x + i∆k y ) m L,R ∓ m 0 0 0 0 m L,R ∓ m V D (∆k x -i∆k y ) 0 0 V D (∆k x + i∆k y ) m L,R ± m        . (4.38) Above, m L = m, m R = (1 + g)(Ω g -m)
with Ω g = gΩ D /2(1 + g). The sign + corresponds to the left side, while -to the right. From herein, in the left side, we have g = 0. The dispersion relations of bulk modes can be written as,

∆Ω = m L,R ± (1 + g) m 2 + V 2 D (∆k 2 x + ∆k 2 y ). ( 4

.39)

The edge wave should propagate along the interface (x = 0) and decay into the bulk (x = 0).

According to Eq. (4.39), the following relations are obtained from the dispersion relations,

∆k L x = iκ L = i m 2 V 2 D + ∆k 2 y - (∆Ω -m) 2 V 2 D , (4.40a)
∆k R x = -iκ R = -i m 2 V 2 D + ∆k 2 y - [∆Ω -(1 + g)(Ω g -m)] 2 (1 + g) 2 V 2 D ; (4.40b)
Eqs. (4.40) implies that for a given frequency, there are one S + and one S -modes. The corresponding eigenmodes of each side can be chosen as,

S L + =        α L 1 0 0        =         V D (∆k L x -i∆k y ) ∆Ω -2m 1 0 0         ; S L -=        0 0 1 β L        =         0 0 1 V D (∆k L x + i∆k y ) ∆Ω -2m         , (4.41a)
S R + =        α R 1 0 0        =         V D (∆k R x -i∆k y ) f ∆Ω + 2m 1 0 0         ; S R -=        0 0 1 β R        =         0 0 1 V D (∆k R x + i∆k y ) f ∆Ω + 2m         . ( 4

.41b)

On the interface (x = 0), the edge wave solution should be the combination of the bulk modes, which takes the form,

ϕ edge = Σ        S ↑ + S ↓ + S ↑ - S ↓ -        =        S ↑ + S ↓ + S ↑ - S ↓ -        L +        S ↑ + S ↓ + S ↑ - S ↓ -        R , (4.42) with 
      S ↑ + S ↓ + S ↑ - S ↓ -        L = (K L + S L + + K L -S L -)e -i∆kyy ;        S ↑ + S ↓ + S ↑ - S ↓ -        R = (K R + S R + + K R -S R -)e -i∆kyy , (4.43)
where K L,R ± are the coefficients, indicating the reltive contribution of different modes. Around the Γ point, the boundary condition in Eq. (4.36) can be approximated as (the u components are ignored),

       -η ′ 0 η 0 0 0 0 η + η ′ -η -η ′ 0 0 0 0 -η 0 -η ′               S ↑ + S ↓ + S ↑ - S ↓ -        L +        η ′ 0 -η 0 0 0 0 -η + η ′ η -η ′ 0 0 0 0 η 0 η ′               S ↑ + S ↓ + S ↑ - S ↓ -        R = 0. (4.44)
As shown in Fig. 4.8, we are looking for the edge wave solutions in the overlapping band gap

Ω 2 D < Ω 2 < Ω 2 D /f . Using the approximation Ω 2 ≈ Ω 2 D + 2Ω D ∆Ω, it is easy to obtain that 0 < ∆Ω < (1 + g)Ω g (
       -cα L cα R 1 -1 0 0 (1 + c)β L -(1 -c)β R -(1 + c)α L (1 -c)α R 0 0 -1 1 -cβ L cβ R               K L + K R + K L - K R -        = 0. (4.45)
Above, we set c = η ′ /η. Eq. (4.45) predicts the relations between components,

K R + = (1 + c)α L (1 -c)α R K L + , K R -= (1 + c)β L (1 -c)β R K L -. (4.46)
Thus, Eq. (4.45) can be reduced,

    2c 2 1 -c α L 1 - (1 + c)β L (1 -c)β R -1 + (1 + c)α L (1 -c)α R 2c 2 1 -c β L     K L + K L - = 0. (4.47)
Based on Eq. (4.47), it is possible to get analytical results for the existence of the tiny gap in the dispersion curves of edge waves.

Quasi-topological edge waves

The original of the tiny gap can be physically investigated. For simplicity, we only discuss the case when ∆k y = 0 below.

• The masses of both sides are equal

In this case, the scaling factor f = 1 (g = 0). For topological gapless edge wave, ∆Ω should be zero when ∆k y = 0. Then Eq. (4.45) leads to,

K L + = K R + , K L -= K R -. (4.48)
which predict two independent edge wave solutions:

ϕ ± = S L ± + S R ± . (4.49)
Eq. (4.49) shows that the edge states do not depend on perturbation η ′ , thus ϕ ± are the timereversal pair that support the edge modes with different spin polarization, which are perfectly topological-protected. However, from Eq. (4.40), we can obtain ∆k L x = ∆k R x = 0. Thus the modes in this case are not evanescent into the bulk. They are not the localized modes on the edge.

• The mass in the right side is modified In this case, the scaling factor f = 1 (g = 0). When ∆k y = 0, α L = β L and α R = β R . Suppose we still have ∆Ω = 0. From Eq. (4.40), it is clear that α L = β L = 0, while α R = β R = 0. Then Eq. (4.47) becomes, 0 1

-1 0 K L + K L - = 0. ( 4.50)
The determinant of Eq. (4.50) is nonzero, meaning that there is no edge wave solution. Consequently, ∆Ω must be non zero when g = 0. For non zero ∆Ω, from Eq. (4.40) it is apparent that α L = β L = 0 and α R = β R = 0. Thus, the zero determinant of Eq. (4.47) leads to the relation, where a = (1 -c)/(V D ∆k L x ) and b = (1 + c)/(V D ∆k R x ). Eq. (4.52) shows that there exists a gap for the edge wave at the Γ point. Thus, we can conclude that when f = 1, the mismatch of the masses on two sides of the interface leads to a gap for the edge wave.

1 -c α L - 1 + c α R = ±2c 2 . ( 4
The existence of a tiny gap has an influence on the nature of edge waves. Actually, in this case the uncoupled time-reversal pair of edge waves hybridize to each other, leading to nonprefect topological edge waves on the interface. From Eq. (4.47) it is clear that K L + and K L are not independent from each other when α L = β L = 0 and α R = β R = 0, which directly shows that the S + modes couple with S -modes. The coupled modes are the eigenstates of the system. The new edge waves can be theoretically predicted. According to Eq. (4.47), the following relation can be obtained,

K L + = 1 -c 2c 2 α L - 1 + c 2c 2 α R K L -= ±K L - (4.53)
The edge wave solution can take the form, and the modification on the mass f , which is not the case for non-trivial topological edge wave whose eigenstates are independent from perturbations [START_REF] Ma | Guiding electromagnetic waves around sharp corners: Topologically protected photonic transport in metawaveguides[END_REF]. According to the results above, we can conclude that when defects appear, the spin-polarized states can be scattered from one into another, turning the edge waves from being perfect immune to being non-perfect immune to the defects. Consequently, backscattering exists in the system as discussed in Fig. 4.9(g). However, as can be seen in Fig. 4.9(b), the width of the tiny gap is ∆ = 5e -5 when f = 0.92, suggesting that the hybridization of the spin-polarized helical edge states is very weak. Therefore, the backscattering is too small, and the propagation of rotational edge waves can be regarded as quasi-topologically protected. The rotational edge waves in the GG are still quite robust against defects. This can be confirmed by using the spatio-temporal evolution simulations. As an example, in Fig

ϕ edge = ±(1 -c)α R S L + + (1 -c)α R S L -± (1 + c)α L S R + + (1 + c)α L S R -e -i∆kyy . ( 4

Conclusion

The quasi-topological transport of rotational edge waves in the GG in a similar way to the QSHE have been demonstrated. The mechanical spin and helical states in the GG can be observed around the Γ point. By introducing perturbation on bending and torsional rigidities, effective SOC can be achieved, leading to the quasi-topological transport of rotational edge wave in the band gap. We show that the interface constructed by two GGs with different masses of the particles does not support the propagation of perfectly topological edge wave since it hybridizes the spin states. However, since the hybridization is very weak, the edge waves can be regarded as nearly topologically protected edge waves. The easy construction of granular crystals, combined with other features such as tunabilitiy upon external fields/mechanical loading and enhanced in a semi-infinite mechanical granular graphene structure. The band structures of the edge waves are theoretically derived and analyzed for zigzag and armchair edge configurations. Similar to the bulk bands, the dispersion of the edge states is also tunable by modifying bending and torsional rigidities of contacts. Since the bending/torsional rigidities are weak in granular crystals, quasiflat bands of edge wave with near zero frequency are observed. These quasi-flat dispersion curves, supporting the propagation of waves with extremely slow group velocity, show a tendency to be perfect zero-frequency modes for zero torsional rigidity or vanish for zero bending rigidity. In particular, when torsional rigidity is weak and bending rigidity is non-zero, the extremely slow edge modes can exist on both zigzag and armchair edges. When bending rigidity also becomes weak in addition to torsional one, all the branches of edge states in the low frequency area are shifted and concentrated to be near zero frequency and very close to the bulk modes, leading to the weak localization of edge waves on the boundaries.

The studies of bulk and edge waves made it possible to investigate the topologically protected edge wave propagation in granular crystals. We found that granular graphene exhibiting double Dirac cone at the Brillouin zone center is a potential structure to study the topological wave in granular system. Those topological wave effects, e.g. unidirectionality and insensibility to defects, can be explored in granular graphene. The notions such as effective spin, helicity and effective spin-orbit coupling are established. The existence quasi-topological rotational wave on the interface constructed by two topologically distinct granular graphenes with different masses of beads is predicted.

To conclude, the presented research results provide fundamental understandings of wave propagation in two-dimensional granular crystals and granular metamaterials and promote their potential applications in acoustic wave control and vibration isolation. We believe that the study of the zero-frequency modes and slow modes is useful for the potential realistic design of granular crystals supporting controllable propagation of extremely slow edge waves. In addition, it paves the way for the analysis of a variety of possible modifications of the granular graphene with the goal of designing granular metamaterials supporting one-way propagating edge states.

The easy construction of granular crystals, combined with other features such as tunabilitiy upon external fields/mechanical loading and enhanced nonlinear response make them a perfect testbed for further fundamental studies in the field of topological mechanics.

We hope that the study of the topological wave propagation in this work will motivate the enthusiasm for the scientific research on topological phenomena in granular system, and the experimental design for topological transport of edge wave based on granular crystals and metamaterials.

be clear, we take the out-of-plane motion as an example to discuss. For in-plane motion, the results can be obtained by following the same procedure. For out-of-plane motion, the green bonds represent the total couplings of three types of effective springs between the beads: (1) shear, (2) bending, and (3) torsional springs. Following analysis of interactions between beads in the main text in Sec. (m, n ∓ 1). Therefore, Eq. (A.2) can be regarded as the equation of motion of cell (m, n). For a GG of finite size, we can get its eigenmodes from its dynamical equation, which can be obtained motion components of the rectangular cells in (0, n) to be zero, i.e. X 0,n = 0 and to modifying the rigidities of the links connecting the beads with the wall, if required.

• Periodic boundary

As shown in In the discussion above, we consider that the boundaries were constructed along y-direction.

If we consider the boundaries along x-direction, similar boundary conditions can be formulated.

In general, the dynamical equation of the finite GG can be obtained by taking into account all the cell indices. As an example, let us consider the case in

Abstract

Granular crystals are spatially periodic structures of particles, most often spherical homogeneous elastic beads, arranged in crystal lattices. The beads are linked by interconnections, which are of much smaller dimensions and weights than the beads. The interactions between beads take place predominantly via normal and transverse rigidities of these elastic interconnections. These induce propagation of elastic waves in granular structures at significantly slower velocities than in the individual grains.

In addition, the existence of non-central shear forces between the particles can initiate the rotations of particles, leading to the rotational and rotational-translational coupled waves in the crystals. In the manuscript, the elastic waves dynamics in twodimensional monolayer granular crystals with either out-of-plane or in-plane particle motion is studied. The dispersion curves of the elastic waves, which can be tuned by controlling the rigidities of contacts, are theoretically derived and numerically analyzed for various configurations of bead couplings. The phononic properties of granular crystals are investigated, including the existence of Dirac points, zero-frequency modes, zero-group-velocity modes and their transformation into slow propagating phononic modes due to weak bending and torsional intergrain interactions. In addition, in the presence of edges/boundaries, zero-frequency and extremely slow elastic edge waves can be also predicted in mechanical granular honeycomb crystals (granular graphene) when the torsional or/and bending rigidities are weak/vanishing. In addition, by a proper design of the structure, the critical topological phase transition can take place in granular crystals. Topological properties of rotational edge waves in a granular graphene are theoretically demonstrated. Dirac dispersion of two degenerate uncoupled modes can be observed at the center of Brillouin zone. A perturbation, coupling these two modes, turns the topological order of granular graphene from trivial to nontrivial, leading to the possibility of topological transport of rotational edge waves in the mechanical granular structures. The developed theories could promote the potential applications of designed granular structures with novel elastic wave propagation properties.

Contents 1 . 1

 11 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Phononic crystals and metamaterials . . . . . . . . . . . . . . . . . . 1.2.1 Phononic crystals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Metamaterials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Granular phononic crystals . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Contacts between particles . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Rotational degrees of freedom . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Tunability and nonlinearity . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Lattice and band theory . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Unit cell and lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Reciprocal space and band structure . . . . . . . . . . . . . . . . . . . . 1.5 Topologically protected wave propagation . . . . . . . . . . . . . . . . 1.5.1 Time reversal symmetry breaking . . . . . . . . . . . . . . . . . . . . . 1.5.2 Time reversal symmetry invariance . . . . . . . . . . . . . . . . . . . . . 1.6 Topological band theory . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.1 Topological invariant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.2 Berry phase and Chern number . . . . . . . . . . . . . . . . . . . . . . . 1.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: (a) Schematics of phononic crystal of different dimensionalities. (b) An example of phonon band structure. The orange area marks the forbidden band gap. (c) Examples of phononic crystals efficient at different wavelength scales.The phononic crystal for hypersound is constructed of air holes in an epoxy matrix with the lattice periodicity a = 1.36 µm and the one for heat is in a silicon matrix of air holes with a = 34 nm.[START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF]

Figure 1 . 2 :

 12 Figure 1.2: Phononic crystal and band structure due to local resonance. [60] (a) Cross section of a coated lead sphere that forms the basic structure unit. (b) Schematics of the phononic crystals. (c) Calculated (solid line) and measured (circles) amplitude transmission coefficient. (d) Band structure.

Figure 1 . 3 :

 13 Figure 1.3: (a) Negative refraction. (b) Negative bi-refraction. (c) Focusing.[START_REF] Lu | Negative birefraction of acoustic waves in a sonic crystal[END_REF]

Figure 1 . 4 :

 14 Figure 1.4: (a) A schematic drawing of a membrane-type artificial atom to achieve negative effective mass density. (b) Effective mass density as a function of frequency. [54]

Figure 1 . 5 :

 15 Figure 1.5: Acoustic asymetric propagation based on non-linear effect [63].

Figure 1 . 6 :

 16 Figure 1.6: Schematics of contact between two elastic spheres.

Figure 1 . 7 :

 17 Figure 1.7: (a) Schematics of a hexagonal granular crystal. (b) Shear coupling on the contact between two elastic spheres. (c) Band structure of the granular crystal.[START_REF] Tournat | Elastic waves in phononic monolayer granular membranes[END_REF]

Figure 1 . 8 :

 18 Figure 1.8: (a) Experimental setup of and the three-dimensional granular crystal. (b) Dispersion curves. (b)Experimental results from a shear transducer (left) and a longitudinal one (right).[START_REF] Merkel | Experimental evidence of rotational elastic waves in granular phononic crystals[END_REF]

Fig. 1 .

 1 7(a). The predicted propagating modes involve an out-of-plane displacement and two rotations with axes in the crystal plane.Shear coupling at the contact is demonstrated inFig. 1.7(b). The pure rotational (red curves) and rotation-translation coupled modes (yellow and blue curves) are predicted in the dispersion relations, Fig. 1.7(c). The first experimental observation of the propagation of coupled rotational-translational elastic waves due to rotational degrees of freedom was achieved in a three-dimensional hexagonal closely packed granular crystal [23], as shown in Fig. 1.8(a). The dispersion relations evaluated including the rotational degrees of freedom are described in Fig. 1.8(b). The experimental results are depicted in Fig. 1.8(c), where coupled rotational-translational (T R/RT ) modes are observed.

Figure 1 . 9 :

 19 Figure 1.9: (a) Microscope image of the interface between the part of the surface covered by a monolayer of beads and uncovered part of the surface. The scale bar is 10 µm. (b) Schematic of the laser ultrasonic experimental setup.(c) Transmission spectra for SAWs propagating across the interface between the surface and the metasurface.[START_REF] Hiraiwa | Complex contact-based dynamics of microsphere monolayers revealed by resonant attenuation of surface acoustic waves[END_REF]

  . 1.10(a). By changing the strength of the magnetic field (changing the magnets B 1 to B 2 , Fig 1.10(b)), the contacts between particles are modified. Consequently, the dynamic response of the granular chain is tuned. The authors use the experimental setup in Fig. 1.10(c) to observe and evaluate elastic wave propagation in the chain. Experimental evidence of transverse-rotational modes of propagation is presented in Fig. 1.10(d), where band structure tunability is observed by changing the magnet B 1 (red curves) to B 2 (black curves).

Figure 1 .

 1 Figure 1.10: (a) Schematics of the magneto-granular chain. (b) Two different magnets in experiments. (c) Schematics of the laser experimental setup. (d) Linear spectrum of the velocity signal received for a chain of 15 beads under B1 (red) and B2 (black) permanent magnets. [92]

Figure 1 . 11 :

 111 Figure 1.11: Magnitude of the acceleration transfer function (detected acceleration A1 at bead 1 over the acceleration A0 at bead 0) of a "two-beads-one-contact" system for different excitation levels, experiment (a) and theory (b). (c) Relative resonance frequency shift as a function of the detected resonance amplitude. (d) Setupfor the single contact characterization in torsion.[START_REF] Cabaret | Nonlinear hysteretic torsional waves[END_REF]

Figure 1 . 13 :

 113 Figure 1.13: Schematics of a one-dimensional diatomic granular chain.

Figure 1 . 14 :

 114 Figure 1.14: Dispersion curves of vibrational modes in the diatomic chain when M/m = 2.

Figure 1 .

 1 Figure 1.15: (a) Experimental setup for measuring the one-way edge state. The inset is a picture of the photonic crystal. (b) The band structure of the gapless edge state (red) between the second and third bands of non-zero Chern numbers. (c) Simulation of the one-way mode and its robustness against a defect. (d) Experimental transmission of edge wave.[START_REF] Wang | Observation of unidirectional backscattering-immune topological electromagnetic states[END_REF]

Figure 1 .

 1 Figure 1.16: (a) A schematics of the phononic crystal with circulating fluid. (b) The band structure of the gapless edge state (red). (c) Simulation of the edge waves against different defects. (d) Simulation of the unidirectional property of the edge waves with different angular velocities of the circulating fluid. [121]

Fig. 1 .Fig. 1 .

 11 Fig. 1.16(a). Due to the existence of the circulating fluids, TRS is broken, leading to a nontrivial band gap, and the existence of topologically protected acoustic edge waves (red curves in Fig. 1.16(b)). To verify the robustness of the edge waves, different defects are introduced in the path of wave propagation in the simulations, including a cavity, a Z-shape and a 180-deg bend, Fig. 1.16(c) top to bottom. The unidirectional property of the edge waves is shown in Fig. 1.16(d), where a point sound source with the midgap frequency of the band gap is placed near the upper boundary. It shows that waves propagate to the left along the interface when the fluid is circulating with positive angular velocity. Changing the value of the angular velocity (Fig. 1.16(d) bottom) does not affect the unidirectional property of the edge wave. This is an example of acoustic analogue of the electronic edge states occurring in the QH effect.

Figure 1 .

 1 Figure 1.17: (a) A schematics of the hexagonal lattice. (b) A schematics of a gyroscope with the top tip pinned to a mass in the lattice. (c) Band structure of the gyroscope phononic crystal with the Chern numbers labeled on the bulk bands. (d) 1D band structure with bulk bands (black dots) and edge bands (color lines). (e) Snapshots of displacement field of the edge wave at different time.[START_REF] Wang | Topological phononic crystals with one-way elastic edge waves[END_REF]

Figure 1 .

 1 Figure 1.18: (a) Time reversal pair. (b) The bi-anisotropic photonic crystal. (c) Band structure with and without the bi-anisotropy of the cylinders in (b). (d) A pair of spin-polarized helical edge states in the gap. (e)Unidirectional property and the robustness of edge waves against defects.[START_REF] Khanikaev | Photonic topological insulators[END_REF]

Figure 1 . 19 :

 119 Figure 1.19: Experimental observation of acoustic analogue QSH effect. [131] (a) Sample of the acoustic metamaterial with interfaces. (b) Simulation of spin "-′′ state incident from port 1. (c) Simulation of spin "+ ′′ state incident from port 2. (d) The measured data on each port of (b). (e) The measured data on each port of (c).

Figure 1 .

 1 Figure 1.20: (a) Perspective view of the phononic crystal with broken z-mirror symmetry. (b) Band structure showing a topological gap. (c) Eigenmodes at the K point.[START_REF] Mousavi | Topologically protected elastic waves in phononic metamaterials[END_REF]

Figure 1 .

 1 Figure 1.21: (a) Illustration of two one-dimensional pendula. Each pendulum has three springs from the top to the bottom labeled by different colors. The top views of the unit cell are shown in the right of (a). The couplings between those springs from the top to the bottom are also depicted. (b) and (c) Steady states with a sequence of removed sites illustrating the stability of the edge states against boundary roughness. Steady states are displayed by colored disks representing pendula's polarization. The radius of the circle corresponds to their mean deflection; the black line indicates the position of the pendulum at a given fixed time. It shows that the motions are mostly confined on the edges. [135]

Figure 1 .

 1 Figure 1.22: (a) Geometric objects divided into into three different topologies. (b) Two interfaces formed by different (right) and same (left) topological materials. (c) Frequency bands of the structures. A topological phase transition takes place on the right, but not on the left. (d) Band structures of edge states for different interfaces. ∆C is the change in Chern number between the corresponding bulk bands on the right and left of the materials.The magnitude of ∆C equals the number of gapless edge modes.[START_REF] Lu | Topological photonics[END_REF]

  )where θ ∈ [0, 2π] is the azimuthal angle of e(k) and φ ∈ [0, π] is the polar angle. In these coordinates, e x = |e| sin θ cos φ, e y = |e| sin θ sin φ, e z = |e| cos θ.(1.24)

Figure 2 . 1 :

 21 Figure 2.1: Possible movements of individual particle. (a) Out-of-plane motion is accompanied by two inplane rotations ϕ, φ and one out-of-plane displacement u. (b) In-plane motion is accompanied by two in-plane translations ux, uy and one out-of-plane rotation ϕ.

2 . 1 . ( 1 )

 211 Out-of-plane motion. As shown in Fig.2.1(a), the beads in the granular structures exhibit out-of-plane translation (u) along z-axis and in-plane rotational angles ϕ and φ (ϕ-rotation with the axis in the x-direction and φ-rotation with the axis in the y-direction). (2) In-plane motion. As presented in Fig.2.1(b), the beads possess out-ofplane rotation (ϕ) along z-axis and in-plane translations u x and u y along x-axis and y-axis, respectively. The movements of individual particles lead to the following forces and/or moments (see Fig. 2.2): (1) Shear forces, which are characterized by an effective shear rigidity (Fig. 2.2(a)).

Figure 2 . 2 :

 22 Figure 2.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 non-central 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.

ξFigure 2 . 3 :

 23 Figure 2.3: Definitions of the unit vectors.

Figure 2 . 4 :

 24 Figure 2.4: Structures of a hexagonal granular membrane and its Brillouin zone.

A

  honeycomb granular membrane, a two-dimensional infinite monolayer honecomb granular crystal, is shown in Fig. 2.5. The lattice constant of the honeycomb structure is a = 2 √ 3R. The corresponding first Brillouin zones is also shown in Fig. 2.5.

Figure 2 . 5 :

 25 Figure 2.5: Structures of a honeycomb granular membrane and its Brillouin zone. The red box marks one unit cell of the structure.

Figs. 2

 2 Figs.2.6(a) and (b)). In the particular case of 3η b + η t = 8/5, the gap disappears (Fig.2.6(c)),

Figure 2 . 6 :

 26 Figure 2.6: Dispersion curves of the hexagonal membrane with out-of-plane motion when η b = 0.1 while ηt is varying. Red lines represent purely rotational modes. Blue (yellow) curves correspond to coupled translational/rotational modes with the predominance of translation (rotation). A complete band gap appears above Ω = 3/2 when 3η b + ηt 9/5. Its width stabilizes at 3/2 Ω √ 15/2 when 3η b + ηt 3.

Figure 2 . 7 :

 27 Figure 2.7: (a) The zero-frequency mode along ΓM, propagating in y-direction, in hexagonal membrane with out-of-plane motion and corresponding movements of the beads at three different points. (b) Γ point: kx = ky = 0; (c) The point of half ΓM: kx = 0, ky = π/4R; (d) M point: kx = 0, ky = π/2R.

  Figs. 2.7(b), (c) and (d) demonstrate the movements of beads in three different positions along the ΓM direction. At

  (2.13) give:(1) an u mode at Ω u = 3/2, (2) a Φ mode at Ω Φ = 15(1 + η b + η t )/4, and (3) a Ψ mode at Ω Ψ = 15(1 + η b + η t )/4. Thus, this Dirac cone originates from the degeneracy of two rotational modes at frequency Ω = 15(1 + η b + η t )/4. Basically, as shown in Figs. 2.8(a) and (b), the variation of the values of η b and η t cannot break the Dirac cone, but can only change its position in the dispersion diagram. As depicted in Figs. 2.8(c) and (d), by controlling η b and η t , the u mode could perfectly degenerate with the rotational modes when the relation η b + η t = 7/15 is satisfied, contributing to a triple degenerated mode at the K point.

Figure 2 . 8 :

 28 Figure 2.8: Dispersion curves in hexagonal membrane with out-of-plane motion for different values of the normalized bending and torsional rigidities. Dirac cones and triple degenerated modes at the K point are 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 degenerate modes as shown in (c) and (d).

2 Ω 15 /2 when η b 1 .

 2151 three components: two in-plane displacements u x and u y and one rotation Φ. The natures of modes are labelled 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 S in , all the modes at the high symmetry points Γ, M and K can be studied analytically. At the Γ point, the eigenmodes of displacement are fixed at Ω u = 0 and the Φ mode appears at Ω Φ = 15/2. For the case of M point, the eigenfrequencies of two u modes are degenerated at Ω u = √ 2, while the Φ mode at Ω Φ = 5(1 + 2η b )/2 depends on η b , which determines the width of the forbidden gap along the ΓM direction. For η b 0, the forbidden gap √ 2 Ω 5(1 + 2η b )/2 grows with increasing bending rigidity. It stabilizes at √ For waves propagating along the ΓK direction,

Figure 2 . 9 :

 29 Figure 2.9: The dispersion curves of hexagonal membrane with in-plane motion when ηn = 1 and η b 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). For 0 η b 1/15, there is no complete band gap. For 1/15 η b 1, the complete forbidden gap appears and its width grows with the increase of η b . It stabilizes at 3/2 Ω 15/2 when η b 1.

  .10, but the nature of this Dirac cone is different from the case of out-of-plane motion. According to the dynamical matrix S in at the K point, three modes exist: (1) an u x mode at Ω x = 9(η n + 1)/8, (2) an u y mode at Ω y = 9(η n + 1)/8, and (3) a Φ mode at Ω Φ = 15(3η b + 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 3η n -15η b = 2, the displacement modes degenerate with the rotational mode, exhibiting a triple degenerated mode at K points, as depicted in Fig. 2.10(c) and (b).

Figure 2 . 10 :

 210 Figure 2.10: Dispersion curves in hexagonal membrane with in-plane motion. Dirac cones and triple degenerated modes at K points are 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).

  2.11(a), (b) and (c), where a complete forbidden band gap can exist for a particular domain of values of η b . For 0 < η b < 0.8, the width of the complete band gap 15(η b + η t )/2 Ω √ 6 decreases when η b is increased. For η b = 0.8, i.e. Fig. 2.11(b), the Φ modes and Ψ modes (Ω = 15(η b + η t )/2) degenerate with the u modes (Ω = √ 6), forming the triple degeneracy at Γ point, and the complete band gap disappears. For 0.8 < η b < 1, a complete band gap reappears, and its width grows when η b is increased. When η b 1, the width

Figure 2 . 11 :

 211 Figure 2.11: Band structures for the out-of-plane motion of the honeycomb membrane with different torsional ηt and bending η b rigidities. 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 ηt = 0 and modify η b . For η b = 0.8, there is always a complete forbidden gap. When η b = 0.8, two rotational modes degenerate with the u mode and the band gap disappears. In (d), (e) and (f), the zero-frequency mode exists when η b = 0. The width of the lower complete band gap grows when increasing ηt while the band gap around Ω = 2 (blue rectangular in (d)) disappears when ηt 0.8.

  Fig.2.12(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, non-propagative due to the lack of bending and torsional rigidities. For the case of η b = 0 only, shown in Fig.2.12(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 η t = 0, Fig.2.12(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.2.12(d), when both η b = 0 and η t = 0, all non-propagating modes in Fig.2.12(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. Ignoring one of the rotation-induced interactions (bending/torsion) would further release the bead from constraints. In the absence of bending and shear rigidity, η b = η t = 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.2.8(a).

Figure 2 . 12 :

 212 Figure 2.12: Dispersion curves of honeycomb lattices in the case of out-of-plane motion. (a) In the absence of bending and shear rigidity, η b = ηt = 0, there are three degenerated zero-frequency modes and one constant frequency mode/oscillation. (b) Only η b = 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 ηt = 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 η b = 0 and ηt = 0 transform all non-propagating modes represented in (a) into slow modes.

Figure 2 . 13 :

 213 Figure 2.13: The zero-frequency mode along ΓM in honeycomb membrane with out-of-plane motion and corresponding movements of beads at three different points. (b) Γ point: kx = ky = 0; (c) The point in the middle of ΓM: kx = π/3a, ky = 0; (d) M point: kx = 2π/3a, ky = 0.

Figure 2 . 14 :

 214 Figure 2.14: 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 degenerate modes. (d) Double Dirac cone at Γ point.

  .24c) From Eq. (2.22), the dispersion curves for η b = 0 are shown in Figs. 2.15(a), (b) and (c), where a complete forbidden band gap could exist for a particular range of values of η n , the ratio of the normal to shear rigidity of the contact. For η n < 1, the width of the complete band gap Ω √ 3 shrinks when increasing η n . For η n 1, the complete band gap disappears. We note that the zero-frequency mode exists only in ΓM direction. Changing the value of η n cannot initiate the propagation of this mode. For the case when η n = 1 and η b = 0 illustrated in Figs. 2.15(d), (e) and (f), the first band gap appears when η b 2/15. For η b 0.4, the second band gap is opened, and its width stabilizes at √ 6 Ω √ 15 when η b 1.

Figure 2 . 15 :

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

Figure 2 . 16 :

 216 Figure 2.16: Band diagrams for honeycomb membrane with in-plane motion, (a) and (b) Triple degenerate modes induced in the K point by different combinations of bending and torsional rigidities. (c) Triple degenerate mode at the Γ point. (d) Double Dirac cone at Γ points.

Figure 3 . 1 :

 31 Figure 3.1: Structure of the granular graphene. Beads in deep blue label the zigzag edge and beads in orange for the armchair edge.

  ) where v A = [u A ; Φ A ; Ψ A ] and v B = [u B ; Φ B ; Ψ B ] are the motion components of particle A and B, respectively. Q ′ , G ′ and N ′ are the 3 × 3 matrices given by Eqs. (2.21). The dynamical equation in Eq. (3.1) leads to the relations,

  So far, we have only used a type of band diagram, e.g. Fig. 3.2(b). This band diagram shows the Ωk relations along the KΓMK directions, which are exactly the edges of the irreductible BZ highlighted by the gray triangle in Fig. 3.2(a). The red and green lines in Fig. 3.2(a) correspond to the symmetry directions of the BZ. Apparently, the band structure in Fig. 3.2(b) reveals only bulk modes propagating in the symmetry directions. If all the wave vectors in the BZ including

Figure 3 . 2 :

 32 Figure 3.2: (a) The first Brillouin zone (BZ) of the granular graphene. The color lines mark different symmetric directions of the BZ. (b) A two-dimensional dispersion curves of the granular graphene. The color lines correspond to bulk modes along the symmetric directions in (a). (c) A three-dimensional band structure of the granular structure. The color lines mark the symmetric directions in (a). (d) Projection of the three-dimensional band structure in (c) to the y-direction.

Fig. 3 .

 3 2(b), gives consistent definitions of the band gaps, the projected band diagram, e.g. Fig. 3.2(d), is much more insightful when boundaries are considered. In fact, it is necessary to use the projected diagram for the analysis of possible edge waves, because the edge waves solutions cannot overlap with propagative waves

Figure 3 . 3 :

 33 Figure 3.3: The projected band structures along qy direction of the graphene for ηt = 0.1 and η b = 0.1 (a), η b = 0.15 (b) and η b = 0.2 (c). The blue areas represent the projected bulk bands along y-direction, and the dash-dot red curves represent the edge states. The amplitudes of three components u, Φ and Ψ of the edge state at point qy = π/3 (labeled by circles in (a)-(c)) as a function of beads position along x-direction are shown for ηt = 0.1 and η b = 0.1 (d), η b = 0.15 (e) and η b = 0.2 (f). The blue curves correspond to sublattice A and the red dash curves to sublattice B.

  .3(a)-(c). In general, four branches of edge modes are obtained, labeled as red dash-dot lines in Fig.3.3(a)-(c), three of which are located in the LPA and another one in the HPA. Tuning the bending and torsional rigidities of contact is possible to control the propagation properties of the edge modes. For example, in the case of η b = 0.1 and η t = 0.1, the edge band below the first propagative bulk band is located very close to the bulk band, leading to edge modes possessing properties similar to those of the bulk modes. The amplitudes of the three components u, Φ and Ψ for wave vector q y = π/3 highlighted in Fig.3.3(a) are shown as a function of bead position along the bulk direction, as shown in Fig.3.3(d), where the blue curves correspond to beads A and the red dash curves to beads B. It implies that the edge mode exhibits weak in-depth localization as the oscillations of the beads spread into the bulk for a rather long distance. When increasing η b to 0.15, as shown in Fig.3.3(b), the edge mode of the wave vector q y = π/3 is easier to distinguish from the bulk modes than in the case of η b = 0.1 in the projected diagram, thus the edge state becomes more localized and its penetration depth into the bulk is shorter (Fig.3.3(e)). For η b = 0.2 in Fig.3.3(c), the edge mode with q y = π/3 has very distinct properties from the bulk mode, leading to the motion of beads localized only on a few layers near the free boundary as shown in Fig.3.3(f). These controllable properties of edge state suggest that the bending/torsional couplings are crucial in the propagation of elastic edge waves belonging to the LPA.It should be reminded that, as discussed in Sec. 2.2 the interactions between the grains through the bending and torsional moments usually are orders of magnitudes weaker than their interactions through the normal and shear forces (ξ b,t ≪ ξ n,s ), leading to the possible existence of extremely slow waves in granular crystals. To further uncover the slow wave nature, we continue modelling the granular crystals as spring-mass systems, Fig.2.2. When rotational DOFs are ignored, the spherical masses interact through forces transmitted by shear and normal springs with rigidties ξ s and ξ n respectively (see Fig.2.2(a) and (d)). The characteristic spatial scale in the granular crystal is of the order of grain radius, R, while characteristic frequencies of the vibrations can be estimated as being an order of ∼ ξ s,n /M . Thus the velocities of bulk waves are of the order, ∼ R ξ s,n /M . As it has been explained above, because of the weakness of contact rigidities ξ s and ξ n , those bulk waves are much slower than the shear and longitudinal waves in the material composing the beads. When the rotational DOFs are taken into account,

  -5 while η b = 0.1, the dispersion curves of edge state represented by red dash-dot curves in LPA are shown in Fig.3.4(a). It is clearly seen that the edge branch below the first propagating bulk area is a quasi-flat band with near zero frequency. We choose the point of q y = π/3 (labeled by a circle in Fig.3.4(a)) to analyze the near-zero-frequency mode. The movements of beads are shown in Fig. 3.4(b), where different colors of beads represent the different out-of-plane, i.e. along the z-axis, displacement of beads and the amplitude of displacement of beads is given in the color bar. The rotational directions of beads are presented by the arrows marked on beads and the amplitudes of rota-

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: (a) Dispersion curves of edge states for η b = 0.1 and ηt = 10 -5 in LPA. A quasi-flat band with near zero frequency can be observed. (b) Movement of beads of the near-zero-frequency mode at point qy = π/3 highlighted in (a) by a circle. The amplitude of the out-of-plane displacement of each bead is normalized. Thus the color marked on each bead represents the amplitude of displacement which can be read from the color bar. The orientation of black arrows on the beads represents the rotating directions of beads. The length of arrows represents the amplitude of rotations of beads, which is normalized to the out-of-plane displacement. The arrow in the right marks the largest amplitude of rotation 0.72.

Figure 3 . 6 :

 36 Figure 3.6: The projected band structure of granular graphene along the qy direction for η b = 2 × 10 -5 and ηt = 10 -5 . (a) One band of edge state exists in the HPA, while the bands of bulk modes in the LPA collapse into a narrow frequency zone. (b) The zoomed view of the LPA. Three bands of edge state can be observed in this area.

  .13f) Similar to the discussion in Sec. 3.3.1, the edge wave dispersion relation can be derived from Eqs. (3.13) and Eq. (3.6), if the motion components are given. For the armchair edge in Fig. 3.1, if

Figure 3 . 7 :

 37 Figure 3.7: Dispersion diagrams for the armchair boundary for a fixed η b = 0.4 but ηt = 0.1 (a), ηt = 0.05 (b), ηt = 0.03 (c) and ηt = 0.01 (d). Red dash-dot curves correspond to the edge states and the dark dots represent projected bulk modes.

Figure 3 . 8 :

 38 Figure 3.8: (a) The dispersion curves of edge modes with slow velocity for η b = 0.2 and ηt = 10 -5 . Two quasiflat bands appear when torsional rigidity is weak. The corresponding movements of beads for the wave vector qx ≈ π/4 noted by the circles are shown in (b) and (c). (b) The oscillation of beads on the edge for the second quasi-flat band. (c) The corresponding oscillation of edge state for the first quasi-flat band. The insert scheme in (a) illustrates the variations of group velocity of the quasi-flat bands as a function of ηt.

Figure 3 . 9 :

 39 Figure 3.9: The dispersion curves of edge modes in the LPA for η b = ηt = 10 -5 and the corresponding movements of the beads in the edge waves with vector qx ≈ π/4 highlighted by circles in (a). (a) Dispersion curves of edge states. Three branches, presented as red dash-dot curves, are found in the LPA. The oscillation for qx ≈ π/4 of the first (b), second (c) and third (d) edge bands are depicted. These edge states can penetrate into the bulk for a rather long distance.

  we fix η b = η t = 0.5 and we monitor the position shifting of the Dirac point by changing P, see the top panels of Fig.4.1. When the perturbation η ′ is non-zero, the Dirac point at the Γ point is lifted, see the bottom panels of of Fig.4.1. However, the existence of a complete band gad depends on the value of P. In this Chapter, we choose P = 1.55 to make sure that there is a complete band gap for non-zero η ′ and the band width is large. Following the discussion above,

Figure 4 . 1 :

 41 Figure 4.1: The influence of P on band structures when η b = ηt = 0.5. Top: perturbation η ′ = 0. Bottom: perturbation η ′ = 0.1.

Figure 4 . 2 :

 42 Figure 4.2: (a) Dirac dispersion and (b) the zoom of the black box in (a). The linear slopes of the dispersion relations near the degenerate point predicted by theory are shown by black lines.

4

 4 is the Dirac velocity, m = η ′ Ω D /2 is the effective mass coming from the perturbation on bending and torsional rigidities. When η ′ = 0, namely m = 0, Eq. (4.6) takes the following form,D 0 (∆k)v 0 (∆k) = Ωv 0 (∆k),(4.8)where D 0 (∆k) = H 0 + ∆H(∆k) with the Hamiltonian of the Dirac pointH 0 = diag(Ω D , Ω D , Ω D , Ω D ). ∆H(∆k) is given by, x -iV D ∆k y 0 0 -iV D ∆k y -iV D ∆k x -iV D ∆k x iV D ∆k y 0 0iV D ∆k y iV D ∆k x ) describes wave behaviors for a small ∆k. The determinant of Eq. (4.8) gives,det(D 0 -ΩI 4×4 ) = [(∆Ω) 2 -V 2 D (∆k 2 x + ∆k 2 y )] 2 . (4.10) For non-trivial solution, the determinant should be zero, leading to ∆Ω = ±V D |∆k|. Thus the eigenfrequency Ω = Ω D ± V D |∆k|, which exhibits the linear dispersion with the slopes ±V D . As presented in Fig. 4.2(b), the black lines mark the slopes ±V D . The degenerate black cone is in a good agreement with the dispersion near the degenerate point obtained from the equations of motion in Eq. (4.1)

13 )

 13 In fact, Eq. (4.13) can be decoupled into two massless Dirac equations:H(∆k)υ ′ = ∆Ωυ ′ , where υ ′ ∈ [p, d], H(∆k) = V D σ • ∆k is a 2 × 2Dirac Hamiltonian, and σ = (σ x , σ y , σ z ) are the Pauli matrices in Eq. (1.23). We can see that under the circular basis, the wave dynamics of the quadruple degeneracy can be described by two seperated massless Dirac equations, manifesting that the quadruple degerneracy corresponds to a double Dirac cone. Consequently, p = [p ↑ ; ip ↓ ] and d = [d ↑ ; id ↓ ] can be regarded as two uncoupled spinors, and the double Dirac cone is the degeneracy of the Dirac cones of p (red branches) and d (green branches) spinor subsystems at the Dirac frequency Ω D , as marked in Fig. 4.2. It is worth noting that the effective Hamiltonian in Eq. (4.13) can be also derived in the band theory of electronic topological insulators [124-126]. Since wave propagation in the GG around the Dirac point can be described by the Dirac equation, the eigenmodes at the Γ point, p ↑↓ and d ↑↓ can be regarded as the effective spin states in the GG. To see the spin nature of modes, the four states p ↑↓ and d ↑↓ are depicted in Fig. 4.3.

Figure 4 . 3 :

 43 Figure 4.3: Possible mechanical spin states in the GG. The left (right) circular polarization of rotation is the spin ↑ (↓) state. The red (blue) arrow represents rotation vector of the particle A (B).

  described by the massless Dirac equation is tuned into a different one characterized by the massive Dirac equation. According to Eq. (4.6), under the spin basis in Eq. (4.12), we can obtain the dynamical equation,H(∆k) mσ z mσ z H(∆k) ψ(∆k) = δΩψ(∆k),(4.14)where ∆Ω = Ω -Ω D -m, and m can be regarded as the "effective mass" due to perturbation η ′ . The existence of the mass term m tunes the dynamical equation of the granular system from a massless Dirac equation to a massive one. It is clear that the presence of the effective mass induces hybridization between the p and d spinors. The spin states p and d are no longer the eigenmodes of the system. The new eigenmodes (U ′ ψ) can be found by diagonalizing Eq. (4.13)

  . The signs "±" determine the nature of the eigenmodes to be S + or S -types. Due to the effective SOC, the Dirac point is lifted and the p and d modes couple to each other, leading to appearance of a gap and the spin-polarized modes S ± . To illustrate the appearance of a gap, first, the unperturbed (η ′ = 0) dispersion curves are shown in Fig. 4.4(a), where a double Dirac cone appears at the Dirac frequency Ω D = √ 3P in the center of the BZ. Modifying the bending and torsional rigidities, the Dirac cone is lifted and a complete band gap appears.

Figure 4 .

 4 4(b) shows the case when the bending and torsional rigidities are perturbed by η ′ = 0.1 (red bonds in the inset of Fig. 4.4(b)). The hybridization of modes are demonstrated in Fig. 4.4(c) and (d). The eigenmodes of the Γ point are shown in Fig. 4.4(c) when η ′ = 0 and in (d) when η ′ = 0.1. Each eigenmode is presented by its two rotational components Φ and Ψ separately. The amplitude of the rotational components can be evaluated from the color bar. It is clear that the eigenmodes in Fig. 4.4(d) are the combinations of the ones in Fig. 4.4(c), confirming the hybridization of modes in the GG.

Figure 4 . 4 :

 44 Figure 4.4: The band structures for η b = ηt = 0.5 and P = 1.55 are shown when perturbation η ′ = 0 in (a) and η ′ = +0.1 in (b). The four rotational branches around the Dirac point are marked from 1 to 4. The red branches in (a) represent the d modes, and the blue ones label p modes. The schematics of particles connected by black (red) bonds in the insets of (a) ((b)) indicates that the bending and torsional rigidities are unmodified (modified). The spin Chern number is shown in (b) for each band. The four eigenmodes of the Γ point when η ′ = 0 in (c) and when η ′ = 0.1 in (d).

2 D ∆k 2 + m 2

 222 and e is a unit vector defined by e = (V D ∆k x /λ, V D ∆k y /λ, m/λ). It is clear that σ • e commutes with H ± , namely [H ± , σ • e] = H ± (σ • e) -(σ • e)H ± = 0. Then the helicity h in the GG can be defined as,h = σ • e.(4.20)From Eq. (4.19), the helicity satisfies the following eigenvalue problem,(σ • e)S ± = ± δΩ λ S ± . (4.21)Note that, the Hamiltonian H ± predicts the dispersion relation δΩ = V 2 D ∆k 2 + m 2 . Thus the helicity h has the eigenvalues equal to ±1, and its eigenstates are the same as the ones of the Hamiltonian H ± . The eigenvalues ±1 determine the eigenstates to be right-handed (+1) or left-handed (-1) helical states. The vector ±(V D ∆k x /λ, V D ∆k y /λ, m/λ) is the helicity vector denoting the direction of the helicity[START_REF] Tkachov | Topological Insulators: The Physics of Spin Helicity in Quantum Transport[END_REF]. It should be mentioned that the definition in Eq. (4.20) is also valid when η ′ = 0, namely m = 0. For right-handed helical states, the helicity is aligned in the same direction as e, while it is the opposite for the left-handed helical states. The

Figure 4 . 5 :

 45 Figure 4.5: Projection of helical modes with perturbation (a) η ′ = 0 and (b) η ′ = 0.1 on the k space. The numbers 1,2,3,4 mark the modes corresponding to those in the branches labeled in the same number in Fig. 4.4. The arrow denotes the axis of helicity projection. The "currents" in the upper branches are outgoing (right-handed) while in the lower branches are incoming (left-handed). The color in the background represents the eigenfrequency.

. 28 )

 28 B is the Berry curvature, which is analogous to the magnetic field. Then according to Stokes's theorem the Berry phase can be written as a surface integralΓ + = S dS • B.(4.29) where S is the surface (BZ) enclosed by the close path C. It defines the Berry phase as the integral of the the Berry curvature over the BZ. Since the Berry connection A = i(∆k x e x + ∆k y e y )/2∆k 2 , it is clear that the Berry curvature B = ∇ k × A = 0 for the eigen function p + (∆k). Thus we can directly conclude that the Berry phase Γ + = 0. The same result is obtained if p -(∆k) is implied into Eq. (4.26). The results for p ± modes are also verified for the d ± modes. Thus we find that the Berry phase for the Dirac cone at the Γ point is zero. The Chern number, the integral of the Berry curvature over a closed manifold quantized in the units of 2π, is zero, indicating that the GG is topological trivial when η ′ = 0.

Figure 4 . 6 :

 46 Figure 4.6: Berry curvature near the Γ point. The numbers 1 to 4 label the four branches near the Dirac frequency from the top to the bottom as shown in Fig. 4.4(b).

Fig. 4 .

 4 6 shows the Berry curvature near the Γ point for the four branches from upper to lower subbands in Fig. 4.4(b).Berry phase, the integral of the Berry curvature over the BZ, is equal to ±π when the integral is performed on B ± . It is clear that the overall Berry phase is zero around the Γ point, which is guaranteed by time-reversal symmetry. However, the Berry phase for each of the branches at the Γ point is nonzero, then the topological number (more accurately, here it is spin Chern number characterizing the integral of the Berry curvature of the same spin over the BZ) is C s = +1, -1, +1, -1 from the upper subbands to the lower ones as shown in Fig.4.4(b). This shows that after the perturbation η ′ , the four branches around the Γ point exhibit non-trivial topological numbers. We can conclude that the crucial point is the Dirac point. When η ′ = 0, the Dirac point is lifted, and the GG undergoes a critical phase shift from the topological trivial to topological non-trivial. Consequently, the gap with band width 2m is topological non-trivial.

Figure 4 . 7 :

 47 Figure 4.7: Schematics of the interface. The unit cells in the bulk and at the interface have been marked by the circles.

η ′ = ±0. 1

 1 do not exhibit overlapping band gap as depicted in Fig. 4.8(a) and (b). In order to ensure that the two GGs exhibit overlapping band gaps, the mass of beads of the right GG has been modified by M R = f M with a scaling factor f = 0.92 (M is the mass of the left side). The dispersion curves in Fig. 4.8(b) and (c) confirm the existence of overlapping gap. Let us focus on the interface shown in Fig. 4.7 and derive the boundary condition. Suppose that the interface is constructed at position x = 0. The masses of beads on the two sides satisfy

Figure 4 . 8 :

 48 Figure 4.8: Dispersions for different parameters. (a) perturbation η ′ = -0.1 and the scaling factor f = 1. (b) η ′ = 0.1 and f = 1. (c) η ′ = -0.1 and f = 0.92. The gray regions mark the common band gap.

  4.5.2). The zooms of the boxes are shown in Fig. 4.9(b).

Figure 4 . 9 :

 49 Figure 4.9: (a) Dispersion of the edge wave along the y-direction. The grey lines correspond to the bulk modes and color lines to the edge modes. (b) The zooms of different edge branches. a tiny gap appears in the topological gap. However, The branches labelled by the blue and red lines still support upward and downward propagating edge modes, respectively. The green lines represent topological trivial edge modes. (c) Schematics revealing the opposite spin nature of modes in the blue and red branches in (b). The arrow in red (blue) in the unit cell represents the rotation direction of particle A (B). Time evolution simulations are implemented for the one-way upward propagating wave in (d) and downward propagating wave in (e). The green stars highlight the positions of source. (f) The same configuration as the one in (c) except that there is a hole close to the interface. The movement of a bead marked by the black circle is recorded to compare with the case in (e). (g) The rotational signals of the recorded bead. The green areas correspond to non-defect case and blue to the defect case.

. 51 )

 51 Finally, we can obtain the edge wave frequencies, ∆Ω = 2m(a -b) ± 2c 2 a -bf (4.52)

. 54 )

 54 It shows that the edge waves are the combinations of the S ± modes. Since α L,R rely on the values η ′ and f , Eq. (4.54) shows that the eigenstates of edge wave depend on perturbation η ′

Figure 4 . 10 :

 410 Figure 4.10: Robustness of the topological rotational wave against defect: a zigzag path containing two corners. (a) When downward propagating edge wave is launched from the position marked by the green star, it shows that rotational wave can turn the corners and keep propagating without noticeable backscattering. For comparison, a normal edge wave is excited from the same positions in (b). It can be seen that wave can not turn the corner and huge reflection can be observed.

  . 4.10(a) a zigzag path with two corners is constructed by combining two GGs with opposite signs of perturbation, η ′ = ±0.1. The same source (star) as in Fig. 4.9(e) is applied to launch the downward propagating rotational edge wave S d . It can be seen that the excited S d can turn the corners and keep propagating without significant backscattering. For comparison, in Fig. 4.10(b), a normal edge wave with q k = π/3 in the green branch of Fig. 4.9(b) is excited from the same position as Fig. 4.10(a),

Figure 4 . 11 :

 411 Figure 4.11: Schematics of the interface with defect.

2 . 2 . 4 ,

 224 the equations of motion for the four particles in the rectangular box can be obtained,-Ω 2 MA m,n = Saa m,n A m,n + Sab m,n B m+1,n + Sad m,n D m,n + Sd m,n D m,n-1 , (A.1a) -Ω 2 MB m,n = Sba m,n A m-1,n + Sbb m,n B m,n + Sc m,n C m,n-1 + Sbc m,n C m,n , (A.1b) -Ω 2 MC m,n = Scc m,n C m,n + Scd m,n D m,n + Sb m,n B m,n+1 + Scb m,n B m,n , (A.1c) -Ω 2 MD m,n = Sdc m,n C m,n + Sdd m,n D m+1,n + Sda m,n A m,n + Sa m,n A m,n+1, (A.1d) where A, B, C and D are the component vectors [u, Φ, Ψ] for beads 1, 2, 3, and 4, respectively. M = diag[1; 1/P; 1/P]. Saa, Sab, Sad, Sd, Sba, Sbb, Sc, Sbc, Scc, Scd, Sb, Scb, Sdc, Sdd,Sda, Sa are all 3 × 3 matrices. Their exact forms can be obtained by analyzing the interactions of each particles. As an example, let us consider particle 1 in the cell (m, n). The equation of motion of particle 1 can be obtained according to Eqs. (2.9)(a)-(c). It should be mentioned that, under the notations in Fig. A.1, particle 1 (the motion component vector is A m,n in Eq. (A.1a)) corresponds to the one labeled as 0 in Fig. 2.5. Its neighbour particle 4 (D m,n ) in cell (m, n) corresponds to the one marked as 3 in Fig. 2.5, particle 2 (B m+1,n ) in cell (m + 1, n) corresponds to the one noted as 1 in Fig. 2.5, and particle 4 (D m,n-1 ) in cell (m, n-1) corresponds to the one labeled as 5 in Fig. 2.5. In addition, the spring rigidity characterized by ξ in Eqs. (2.9)(a)-(c) is replaced by k in the numerical calculations. Therefore, from Eqs. (2.9)(a)-(c), we can easily obtain Saa, Sab, Sad, Sd of the following forms,Saa m,n = k n,m 1s + k n,m 1b -k n,m 1t ) 4 -k n,m 1s + k n,m 1b + 3k n,mFollowing the same procedure for particles 2, 3 and 4, all the matrices in Eqs. (A.1) can be obtained. Here, we will not list all of them but continue to discuss the eigenvalue problem in a finite size granular graphene. The four equations in Eqs. (A.1) describe the interactions involving the four particles in a rectangular cell. We can put them in a 12 × 12 matrix equation to characterize the total interactions of the rectangular cell in position (m, n), -Ω 2 X m,n = SA m,n X m,n + SB m,n X m-1,n + SC m,n X m,n-1 + SD m,n X m,n+1 + SE m,n X m+1,n , (A.2) where X = [B; C; D; A] is a 1×12 vector containing all the motion components of the rectangular cell. SA, SB, SC, SD and SE are 12 × 12 matrices of the following forms, diag[M; M; M; M]. The five matrices above have physical meanings: SA describes the interactions within the rectangular cell; SB and SE demonstrate the interactions of cell (m, n) with cell (m ∓ 1, n); SC and SD characterize the interactions of cell (m, n) with cell

  Fig. A.2(c), periodic boundaries are considered in the two ends of the structure (positions m = 1 and m = M ). In this case, the boundary conditions can be obtained by connecting the two ends together. Thus, the motion components of the cell in the two ends satisfy X 0,n = X M,n and X 1,n = X M +1,n . The springs on the edges satisfy k n,1 2 = k n,M +1 2 .

- 1 ,

 1 Fig. A.1, where the two boundaries at positions n = 1 and n = N (top and bottom edges) are free, while the two at positions m = 1 and m = M (left and right edges) are periodic. First of all, we consider the index n. According to Eq. (A.1), the following equation can be derived,-Ω 2 Y m = Qa m Y m + Qb m Y m-1 + Qc m Y m+1 ,(A.3)where Y m = [X m,1 ; X m,2 ; ...; X m,N ] is a 1 × 12N vector containing all the motion components of particles in position m. Qa m , Qb m and Qc m are 12N × 12N matrices of the forms, SA m,N -1 SD m,N -1 Qb m = diag[SB m,1 ; SB m,2 ; ...; SB m,N ], Qc m = diag[SE m,1 ; SE m,2 ; ...; SE m,N ]. Qa m describes the interactions of particles within the rectangular cells at position m, while Qb m (Qc m ) describes the interactions of cells at position m with those at position m -1 (m + 1). Next, we need to consider all the indices of m in Eq. (A.3). This leads to the final dynamical equation of the finite size GG, -Ω 2 V = QqV , (A.4) where V = [Y 1 ; Y 2 ; ...; Y M ] is a 1 × 12N M vector containing all the motion components of

Figure A. 3 :.

 3 Figure A.3: An eigenmode of a GG composed of 5 × 25 rectangular cells. Three motion components of beads are presented separately. The top and bottom edges are free boundaries, while the left and right edges are periodic boundaries.

  1,2 )] are the shear modulus of each particle. The relation in Eq. (1.3)(a)

Table 1 .

 1

	1: Two-dimensional lattice
	crystal family	Bravais lattice
	Monoclinic	Oblique
	Orthorhombic Rectangular, Centered
	Tetragonal	Square
	Hexagonal	Hexagonal

n is the common name for a distinct lattice type. Note that, in two dimensions, the 5 Bravais lattices can be grouped into 4

Table 2 .

 2

			1: Hexagonal granular membrane
	i-th	1	2	3	4	5	6

Table 2 .

 2

2: Honeycomb granular membrane (a) Particle "0"

  see the definition of Ω g in Eq. (4.38)). Substituting Eqs. (4.41) and (4.43) into Eq. (4.44), we arrive at,

1.2 Phononic crystals and metamaterials

Bibliography

Publications nonlinear response make them a perfect testbed for further fundamental studies in the field of topological mechanics. In addition, the study in this work could promote the potential applications of granular structures designed for the directional transport of elastic waves.

General conclusion

In this manuscript, we studied theoretically the linear wave dynamics and the topological property of elastic waves in two-dimensional (2D) granular crystals or granular membranes. Granular crystals as a particular type of phononic crystals exhibit excellent abilities in manipulation of elastic waves. Due to the rotational degrees of freedom, purely rotational and coupled rotationaltranslational modes can exist in granular systems. Regardless of these properties, we found that zero-frequency modes and slow modes can be observed in the two-dimensional granular systems.

Nevertheless, inspired by the study on topological wave behaviour in quantum/photonic/phononic systems, the possible existence of topological waves propagation in granular crystals was also investigated. Quasi-topological rotational waves on the interface are theoretically predicted and numerically confirmed in granular graphene, granular honeycomb crystals.

Considering the rotational degrees of freedom, we firstly studied the bulk wave behaviour in the hexagonal/honeycomb monolayer membranes when out-of-plane or in-plane motions are imposed. The tunability of the band structures for phononic modes as well as the complete band gaps are demonstrated. These properties are similar to those in other known phononic crystals.

However, distinguishing them from other phononic crystals, the existence of the non-propagating rotational modes and slow rotational modes is very common in granular crystals, especially when the bending/torsional rigidities are small. In the extreme case when bending/torsional rigidities is absent, the non-propagating modes are observed in both hexagonal and honeycomb membranes with out-of-plane motion.

To understand the physics of these non-propagating modes and slow modes, we analyzed the case when torsional rigidity is zero. It is found that due to the existence of rotational degrees of freedom, the non-propagating modes are predicted in the system when bending and/or torsional rigidities are zero. Their existences in the membrane are analyzed. For example, when wave vector is zero, rotation of the particles is not launched but the relative translations of adjacent beads are zero. When wave vectors is non-zero, the elongation of the effective springs between beads due to their unequal translations could be compensated by shrinking of the effective springs caused by their relative rotations. Consequently, both cases keeps the contacts between the particles unloaded and the energy of the system unmodified. For weak bending/torsional rigidities of the contacts, the balance of the compensation is broken, transforming the nonpropagating modes into propagating modes with extremely slow speed.

Based on the results obtained for bulk wave dynamics, we analyzed the edge wave dynamics (2) Fixed boundary; (3) Periodic boundary. We will discuss them below separately.

• Free boundary

As an example, Fig. A.2(a) shows a free boundary at positoin m = 1. In this case, the particles on the edge have zero interaction with those in the left side (m < 1). Mathematically, this is equivalent to setting the corresponding spring rigidities to be zero. For example, in In this part, the spatio-temporal simulation for elastic wave propagation in the structure is discussed. As discussed in Sec. A, Eq. (A.4) gives the eigenfrequencies and the eigenmodes of a finite size GG. However, in numerical analysis, if we also want to see the propagation of elastic wave in the structure, then the dynamical equation of the system should be time-dependent. which is a second order ordinary differential equation (ODE). Actually, in Eq. (B.1), V can be regarded as the "displacement" of the system, then V is the "acceleration". Therefore, we can introduce the "velocity" of the system U , which satisfies the following relations: 3). This equation can be rewritten as,

Let us still consider the case in

where U i is the velocity at time t i . From Eq. (B.4), the following relation can be obtained,

where ∆t = t i+1 -t i . Similarly, we can obtain another equation from Eqs. (B.3),

The RK4 method for the ODE problems in Eqs. (B.5) lead to [START_REF] Cellier | Continuous System Simulation[END_REF][START_REF] Carnahan | Applied numerical methods[END_REF],