Granular Graphene: direct observation of zigzag and armchair edge waves

TL;DR

This paper introduces a mechanical granular graphene model using steel beads, demonstrating Dirac points and edge wave phenomena, providing a new platform for studying advanced wave behaviors.

Contribution

It presents the first experimental realization of granular graphene with observed Dirac points and edge wave dynamics, highlighting its potential for exploring topological and nonlinear wave phenomena.

Findings

Observation of Dirac points in granular graphene

Experimental evidence of edge wave propagation and turning effects

Granular graphene as a platform for novel wave phenomena

Abstract

We propose a mechanical granular graphene obtained by replacing the carbon atoms with macroscopic spherical stainless steel beads in contact. The experimental measured dispersion relation is presented, in conjunction with evidence of the Dirac points. In addition, wave propagation along the zigzag and a robust turning effect of edge waves from the zigzag to the armchair boundary is experimentally revealed, even in the absence of a full band gap for bulk modes. Our work shows that mechanical granular graphene can serve as an excellent experimental platform to study novel Dirac, topological and nonlinear wave phenomena.

Figures

Equations

$$n_{\beta }$$

$$s_{\beta }$$

$$b_{\beta }$$

$$M{\ddot{u}}_{\alpha ,\mathrm{m},\mathrm{n}}$$

$$M{\ddot{v}}_{\alpha ,\mathrm{m},\mathrm{n}}$$

$$I{\ddot{\phi }}_{\alpha ,\mathrm{m},\mathrm{n}}$$

$${\ddot{\mathbf{U}}}_{\mathrm{m},\mathrm{n}}^A=S_0{\mathbf{U}}_{\mathrm{m},\mathrm{n}}^A+S_1{\mathbf{U}}_{\mathrm{m},\mathrm{n}}^B+S_2{\mathbf{U}}_{\mathrm{m}-1,\mathrm{n}+1}^B+S_3{\mathbf{U}}_{\mathrm{m}-1,\mathrm{n}-1}^B,$$

$${\ddot{\mathbf{U}}}_{\mathrm{m},\mathrm{n}}^B=D_0{\mathbf{U}}_{\mathrm{m},\mathrm{n}}^B+D_1{\mathbf{U}}_{\mathrm{m},\mathrm{n}}^A+D_2{\mathbf{U}}_{\mathrm{m}+1,\mathrm{n}+1}^A+D_3{\mathbf{U}}_{\mathrm{m}+1,\mathrm{n}-1}^A,$$

$$M_0{\mathbf{U}}_{1,\mathrm{n}}^B+D_1{\mathbf{U}}_{1,\mathrm{n}}^A=0,$$

$$M_1{\mathbf{U}}_{21,\mathrm{n}}^A+S_1{\mathbf{U}}_{21,\mathrm{n}}^B=0,$$

$$M_2{\mathbf{U}}_{\mathrm{m},1}^A+S_3{\mathbf{U}}_{\mathrm{m}-1,0}^B=0,$$

$$M_3{\mathbf{U}}_{\mathrm{m},1}^B+D_3{\mathbf{U}}_{\mathrm{m}+1,0}^A=0,$$

$$M_4{\mathbf{U}}_{\mathrm{m},41}^A+S_2{\mathbf{U}}_{\mathrm{m}-1,42}^B=0,$$

$$M_5{\mathbf{U}}_{\mathrm{m},41}^B+D_2{\mathbf{U}}_{\mathrm{m}+1,40}^A=0,$$

$${\xi }_n$$

$${\xi }_s$$

$${\xi }_b$$

$$\delta$$

$${\ddot{\mathbf{U}}}_{\mathrm{m},\mathrm{n}}^A$$

$${\ddot{\mathbf{U}}}_{\mathrm{m},\mathrm{n}}^B$$

Full text

††preprint: APS/123-QED

I Introduction

Graphene, a single-layer of carbon atoms in honeycomb lattice, has recently emerged as an appealing system for conducting fundamental studies in condensed matter physics and in particular Dirac physics phenomena. Neto ; Sarma ; Geim ; Neto ; Waka . The great advances in experiments in graphene were awarded with the 2010 Nobel prize. Despite the enormous progress, there are still great difficulties in designing/modifying graphene at will at the nanoscale. This leads many researchers to propose and study other artificial microscopic and even macroscopic graphene analogues for further fundamental studies. These settings include the use of molecules review , ultracold atoms ultracold , photons RechtsmanNature ; DuboisNC ; ZhangPRL ; AblowitzPRA or phonons TorrentPRL ; YuNM in honeycomb lattice.

The study of edge wave in finite crystals has been a long studied topic in condensed matter physics Shockley ; Tamm . In genuine graphene, zero-energy electronic states are predicted in nanoribbons with zigzag/beared boundaries Fujita ; Nakada and confirmed by means of scanning-tunneling microscopy Kobayashi ; Niimi . However, armchair boundaries do not support electronic edge states unless defects appear on the edges Kobayashi or the system is anisotropic Kohmoto . The interest on edge states has been significantly renewed by recent advances in the study of topological physics. It has been shown that robust edge states/modes can appear in topological insulators Hasan ; Qi , crystallineFu and higher orderBenalcazar topological insulators. In these systems, the existence of edge states is directly connected with the topological properties of bulk bands. This is also the case for the electronic edge states of graphene, since they are related to the winding number of bulk eigenmodes Ryu ; Mong ; Delplace .

In addition to genuine graphene, an extensive body of works has been published over the last few years, studying edge waves in different artificial graphene structures, particularly in photonics. A number of edge modes has been observed experimentally, such as conventional flat bands as well as unconventional edge branches in zigzag and beared boundaries Kuhl ; Plotnik ; Milicevic . Regarding armchair boundaries, the previous reports of armchair edge waves are either in anistropic microwave artificial graphene BellecNJP , or in a photonic graphene-like structure of coupled micropillars MilićevićPRL , where the existence of edge states is due to the coupling of $p_{x,y}$ photonic orbitals. Considering the study of vibrational edge waves, the existence of edge modes has been predicted in genuine SavinPRB and granular graphene ZhengEML . Moreover, it has also been shown that in mechanical graphene analogues both flat and dispersive unconventional edge modes can be found at zigzag edges under fixed boundary conditions Coriolis ; Kariyado . However, to the best of our knowledge, there is no experimental observation of edge waves in mechanical graphene till now.

In this work, we propose another type of artificial graphene, that is, granular graphene, which can be thought of as a mechanical analogue of graphene whose carbon atoms are replaced by macroscopic elastic beads and chemical bonds are substituted by contact interactions via various stiffnesses. Compared to genuine graphene or other mechanical graphene structures, granular graphene possesses extra physical features that make it very appealing from a fundamental point of view. One of those is the existence of multiple degrees of freedom (translations and rotations) MerkelPRL ; HiraiwaPRL ; PichardPRB ; AlleinAPL . This, in combination with the honeycomb lattice geometry, leads to Dirac cones in the dispersion relation TournatNJP ; ZhengUltrasonics and topological helical edge waves ZhengPRB . From an experimental standpoint, previous reports on two-dimensional (2D) granular crystals usually focus on closely-packed hexagonal or square lattices with mechanical constraints located on the borders GillesPRL ; CostePRE2008 ; PrlChiara . Thus far however, no direct observation has been made on Dirac cones or edge wave propagation on granular graphene. The obstacles blocking further experimental investigations on granular systems include difficulties in constructing different structures and stability issues, particularly for looser packings like the honeycomb structure.

In this article, we present how such structural difficulties have been overcome using external periodic magnetic fields. The proposed magneto-granular graphene (MGG) is structurally stable and acts as a nearly free-standing granular structure. We then experimentally show the dispersion curves, evidence of the Dirac point, and for the first time, direct observations of the edge wave propagation in the MGG. We theoretically show that in the MGG, edge waves can exist not only on free zigzag, but also on free armchair boundaries. More importantly, for a range of frequencies we show the existence of novel bulk modes where translations are extended in the bulk of the structure while rotations are localized at the edges. In the same frequency range, when translations are also constrained at the edges, i.e., in partial band gaps of the bulk modes, edge modes can also appear on the free zigzag or both the free zigzag and armchair edges. This leads to an interesting turning effect of edge waves from zigzag to armchair boundary in the frequency range where edge modes appear on both free zigzag and armchair edges. Aside from the topological wave mechanism, where edge transport occurs in the full gap for bulk wave and is protected by the bulk topology ZhuPRB , the turning effect demonstrated here originates from the coexistence of wave solutions on the zigzag and armchair edges over a certain frequency range. The role of the topology in MGG, like in the recent works of higher order topological insulators, is a remaining intriguing question.

II Experimental set-up and modeling

The MGG is depicted in Fig. 1B, where $820$ stainless steel beads (diameter $d=7.95$ mm, density $\rho=7678$ kg/m3, Young’s modulus $E=190$ GPa and Poisson’s ratio $\nu=0.3$) are precisely placed in a honeycomb lattice, in contact with one another. This layout stems from a properly designed external magnetic field that is induced by permanent cylindrical NdFeB magnets (remanent magnetization $1.37$ T, diameter $6$ mm, and length $13$ mm) placed in a honeycomb configuration within the wood matrix, Fig. 1A. The external periodic magnetic field magnetizes the elastic beads, resulting in equivalent pre-compression forces between beads and thus a mechanically stable structure. Between the elastic beads and the substrate, a thin layer of rubber (thickness $2$ mm) has been set to minimize the mechanical coupling of the granular graphene with the substrate, and to damp the transmission of elastic waves into the wood matrix. The experimental set-up is shown in Figs. 1C. In-plane motion is excited by the driver connected to a piezoelectric transducer (Panametrics V3052). Each bead in the structure exhibits one out-of-plane rotation $\varphi$ around the $z-$axis and two in-plane translations $u$ and $v$ along the $x-$ and $y-$axes, respectively. The $u$ and $v$ components of each bead can be monitored separately by two laser vibrometers, which are sensitive to changes in the optical path length along the beam direction.

Regarding the mechanical contact interactions between adjacent beads, we consider normal, shear and bending interactions, as characterized by the contact rigidities $\xi_{n}$, $\xi_{s}$, and $\xi_{b}$ respectively, Fig. 1E. Once pre-compression has been determined (around $\sim 1.55$ N by means of measurement), $\xi_{n}$, $\xi_{s}$ and $\xi_{b}$ can be obtained from Hertzian contact mechanics Johnson ; Mindlin , see Methods. For the three types of interactions between adjacent beads, i.e., Fig. 1E, the elongations corresponding to the effective normal $n_{\beta}$, shear $s_{\beta}$, and bending $b_{\beta}$ contact springs can be expressed as,

[TABLE]

where $\alpha=A,B$ is the sublattice index. Considering the honeycomb structure, each sublattice bead is in contact with three other beads as denoted by a neighboring index $\beta=1,2,3$. We define $\bm{e}_{\beta}$ as unit vectors in the directions from the center of $\alpha$ bead to the center of its $\beta$-th neighbors. $\bm{e}_{x}$, $\bm{e}_{y}$ and $\bm{e}_{z}$ represent the unit vectors along $x-$, $y-$ and $z-$axes, respectively. $\bm{l}_{\beta}$ are unit vectors normal to $\bm{e}_{\beta}$ and $\bm{e}_{z}$ with the form $\bm{l}_{\beta}=\bm{e}_{z}\times\bm{e}_{\beta}$. As displayed in Fig. 1D, we can label the sublattice $\alpha$ in a normalized coordinate $(\mathrm{m},\mathrm{n})$ (with the bead center positions serving as the coordinate) by $\alpha_{\mathrm{m},\mathrm{n}}$, where $\mathrm{m}$, $\mathrm{n}$ are both integers representing the normalized center positions of beads in the $x-$ and $y-$axes, respectively. On site $(\mathrm{m},\mathrm{n})$, the equations of motion can be expressed as follows,

[TABLE]

Above, $M$ is the mass of a bead and $I$ is its moment of inertia. The dots on the top represent derivation over time. It can be seen from Eqs. (2) that bending interactions can not lead to the displacement of beads, i.e., Eqs. (2a), (2b), while normal interactions do not contribute to the rotation of beads, i.e., Eq. (2c). Based on the equations of motion in Eqs. (2), wave dynamics in the MGG can be described by,

[TABLE]

where $\bm{U}_{\mathrm{m},\mathrm{n}}^{\alpha}=[u_{\alpha};v_{\alpha};\Phi_{\alpha}]_{\mathrm{m},\mathrm{n}}$ with $\Phi=\varphi d/2$ are the motion of particle $\alpha$ in the normalized coordinates. $S_{i}$ and $D_{i}$ ($i=0,1,2,3$) are $3\times 3$ matrices, see supplemental materials (SM). By applying the Bloch periodic boundary conditions in both $x-$ and $y-$axes, i.e., $\bm{U}_{\mathrm{m},\mathrm{n}}^{\alpha}=\bm{U}^{\alpha}e^{i\omega t-iq_{x}\mathrm{m}-iq_{y}\mathrm{n}}$ with the normalized wave vectors $q_{x}=k_{x}d/2$, $q_{y}=\sqrt{3}k_{y}d/2$, Eqs. (3) can be mapped into an eigenvalue problem which leads to the dispersion curves of an infinite MGG as shown in Figs. 2A and  2B.

Considering that the MGG in experiments is of a finite size $21\times 41$, there are free zigzag edges at positions $(\mathrm{m},\mathrm{n})=(1,\mathrm{n})$, $(\mathrm{m},\mathrm{n})=(21,\mathrm{n})$ and free armchair edges at $(\mathrm{m},\mathrm{n})=(\mathrm{m},1)$, $(\mathrm{m},\mathrm{n})=(\mathrm{m},41)$. At the mechanically free boundaries, which can be obtained by removing parts of the neighbors of edge beads, the beads are interacting with a smaller number of neighboring beads than in the volume. Therefore, the boundary conditions are derived from the cancellation of interactions between the removed beads and the edge beads, which leads to the following boundary conditions:

[TABLE]

for the zigzag edges, and,

[TABLE]

for the armchair edges with $M_{j}$ ($j=0,1,2,3,4,5$) $3\times 3$ matrices, see SM. To account for dissipation, a phenomenological on-site damping term BoechlerNat has also been introduced into the right-hand side of Eqs. (3), $-1/\tau\bm{\dot{U}}_{\mathrm{m},\mathrm{n}}^{\alpha}$, with $\tau$ characterizing the decay time of waves. This coefficient has been chosen to fit the experimental results. More details about the dissipation implementation can be found in the Methods section.

III Dispersion curves and Dirac point

To measure the MGG dispersion, in-plane motion has been activated using a frequency sweep excitation from $500$ Hz to $35$ kHz by the bead-driver located at position $(1,22)$. The $u$, $v$ components of particle $B$ in each unit cell are collected by the laser vibrometers. By scanning all particles $B$, the spatial frequency signals of translation can be obtained, which in turn yields the dispersion curves by applying a double Fourier transform. Figures 2A and 2B present the dispersion curves of an infinite granular graphene without dissipation. The color scale level reflects the weights of $u$ (red curves) and $v$ (green curves) components in each mode. The corresponding numerical dispersion curves, mimicking the experimental process, are displayed in Figs. 2C-2D, while the experimental ones are shown in Figs. 2E-2F for the $u$ and $v$ components respectively. Figures 2C-2F indicate that up to $\sim 20$ kHz, the experimental dispersion curves are in good agreement with both the theoretical and numerical curves since the branches are translation-dominated. As expected, the branches with frequencies above $\sim 20$ kHz are absent due to the fact that these modes are rotation-dominated, which are not easily detected by the laser vibrometers.

Interestingly, Figs. 2E and 2F reveal the band crossing at the K point around frequency $10$ kHz. The observation of this crossing provides evidence of the Dirac cone in MGG, originating from the honeycomb lattice symmetry. As shown in Methods, there are theoretically two Dirac cones with the Dirac frequencies $\omega_{\pm}$ in granular graphene considering the in-plane motion. The band crossing around $10$ kHz corresponds to the Dirac point $\omega_{-}$ at the K point of the Brillouin zone (BZ). Note that, another Dirac cone is also predicted around $\omega_{+}\sim 24$ kHz, Figs. 2A-2B. However, this Dirac point is not visible in Figs. 2E and 2F due to the fact that the translation signals around $24$ kHz are weak and thus hidden in the color scale.

In order to observe the $\omega_{+}$ Dirac point in the MGG, another set of experiment has been performed around the target frequency. Experimentally we choose the source to be a frequency sweep excitation from $18$ kHz to $26$ kHz. However, there are still two main difficulties to be overcome: (1) Collection of the weak translational signals. Since the modes are dominated by rotation over this frequency region, the signal of translational components is consequently weak. Thus, this rotation-dominated modes are not easily detected by the laser vibrometers since they are only sensitive to changes due to displacements of beads. When dissipation is also taken into account, the weak translational motion becomes weaker due to the attenuation during propagation. (2) The resolution of the Dirac point. Since the number of the eigenmodes is related to the size of the MGG, larger size of the structure results in a number of eigenmodes, which in turn lead to better resolution of dispersion around the Dirac point. However, as explained in (1), large size of the sample can lead to disadvantages for measurement due to the fact that the translational signals of the particles far away from the source can be too weak to be measured by the vibrometers. As a compromise, in this experiment for the $\omega_{+}$ Dirac point, we chose the size of sample to be $11\times 41$. This keeps the resolution point along the $q_{y}$ unchanged, but decreases the length along $q_{x}$ to reduce the influence of attenuation on the translational signal. By reconstructing a new sample of size $11\times 41$ with stainless steel beads (diameter $8$ mm, density $7650$ kg/m3, Young’s modulus $210$ GPa and Poisson’s ratio $0.3$), the pre-compression between beads in this sample is measured to be around $F_{0}=1$ N. Therefore, we calculate the Dirac frequency to be at $\omega_{+}\sim 22.585$ kHz. By scanning all the particles $B$ and doing the 2D real-reciprocal space Fourier transformation, iso-frequency contour of a given frequency can be obtained. In Fig. 2G, we show the iso-frequency contours at the Dirac frequency $\omega_{+}$ obtained experimentally and numerically considering the same size of sample and dissipation. We observe that the mode at $\omega_{+}$, displayed by the high value in the iso-frequency contours, is only present around the K points which reveals an evidence of Dirac point in the MGG.

IV Edge waves

Another interesting feature that appears in the dispersion of the finite-sized MGG is the existence of branches in the $\Gamma$K and MK directions around $20$ kHz (green ellipses in Fig. 2C-2F. These branches correspond to edge waves and in this section, we will study these first theoretically and then experimentally.

By considering free boundaries in Eqs. (4) and (5), the edge wave dispersion for the zigzag and armchair edges is calculated, see Fig. 3A. In the calculations of the edge dispersion, we assume that the free zigzag (armchair) edges are located at $\mathrm{m}=1$ and $\mathrm{m}=21$ ($\mathrm{n}=1$ and $\mathrm{n}=41$), while along the $y-$ ($x-$) axis are infinite. Therefore, based on Eqs. (3) and the boundary conditions in Eqs. (4) (Eqs. (5)) along with the Bloch periodic conditions in the $y-$ ($x-$) axis, the edge wave dispersions in Fig. 3A are obtained. The gray regions correspond to bulk, while the red (blue) curves to edge wave solutions. In total, $2$ edge branches for the zigzag and $3$ for the armchair are present. This increased number of edge states, especially the existence of edge states at the armchair edge, is not encountered in the genuine graphene. Similar edge modes have been only reported in photonic lattices with orbital bands MilićevićPRL , or mechanical honeycomb lattices with purely in-plane translations Coriolis ; Kariyado , namely 2D mass-spring honeycomb systems with two translational degrees of freedom per sublattice. However, in granular graphene, there are three degrees of freedom per sublattice (two translations $u$, $v$, one rotation $\phi$). As we will see, this significantly enriches the edge physics of this system.

Considering the edge wave in the rotation-dominated region (above $\sim 22$ kHz), as shown in Fig. 3A, it can be seen that no edge modes are found in the zigzag boundary. This is different from the conventional graphene where a flat band of zero edge modes exists in the zigzag edges Kohmoto ; Ryu . In mechanical graphene, similar flat band edge states can be found under fixed boundaries conditions. As it was commented in Ref. [Kariyado ], free boundary conditions, like the ones used in this work, break the chiral symmetry on the free zigzag edges, leading to the absence of a flat band of edge modes. However, a novel armchair edge branch appears, whose topological origin can be further discussed in future investigation. In Fig. 3B, we present a close-up of the edge wave dispersion around $20$ kHz for the zigzag and armchair cases, respectively. Interestingly, there is an overlapping region from $\sim 19.84$ kHz to $\sim 20.07$ kHz, where edge modes can be found on both the zigzag and armchair edges, while due to the absence of a full bulk gap, bulk modes also exist.

To shed more light on the edge physics of the MGG, we carry on eigenmode analysis of the structure. Since the MGG is of a finite size $21\times 41$, a dynamical equation describing wave behavior of the MGG can be derived from Eqs. (3)$-$(5) by taking into account all the coordinate indices $\mathrm{m}$ and $\mathrm{n}$. Consequently, the eigenvalue problem from the dynamical equation of the MGG can be thoroughly solved. In Fig. 4, we choose to show several eigenmodes around the edge wave frequencies. Starting from the rotation-dominated region, we show two eigenmodes with the eigen-frequencies $24.80$ kHz in Fig. 4A, and $24.82$ kHz in Fig. 4D. The color scales of the three components suggest that in this frequency region, modes are dominated by rotation. This is consistent with the measured dispersion curves in Figs. 2E and  2F, where modes are not detected above $\sim 22$ kHz since rotational signals can not be recorded by the laser vibrometers. It can be also seen that the eigenmode in Fig. 4A shows an extended property since all parts of the structure are involved in the motion. However, the eigenmode in Fig. 4D exhibits a property of localization as the motion mostly is confined only on the free armchair boundaries. Note that, due to the boundaries, the eigenmodes of a finite-size graphene can be viewed as the contributions of those bulk and edge modes of the infinite MGG. Vibrational modes strongly localized at the edges of the structure can be called as edge modes of the MGG. Therefore, the extended mode at $24.80$ kHz could be viewed as a mode dominated by the contribution from the bulk modes of infinite MGG of the grey region of Fig. 3, while the eigenmode in Fig. 4D could be dominated by the mode marked by the yellow dot in Fig. 3A in the armchair edge mode branch. Regarding the overlapping region, two eigenmodes close to each other at $19.89$ kHz and $19.88$ kHz are presented in Fig. 4B and 4E, respectively. It shows that the eigenmode in Fig. 4B also exhibits an extended property with translations involving most part of the structure, while the eigenmode in Fig. 4E has a localized mode property similar to the one in Fig. 4D but with motion constrained on both the zigzag and armchair boundaries. Thus, this eigenmode at $19.88$ kHz has a dominant contribution from the edge modes of infinite granular graphene in the zigzag and armchair edge branches labelled by green dots in Fig. 3B. The structure of this eigenmode confirms the prediction of the existence of edge branches in both zigzag and armchair ribbon (blue, red curves in Fig. 3B). Finally, in the region that lies below the overlapping region ($\sim 19.61$ kHz to $\sim 19.84$ kHz), Fig. 4C and 4F, the behavior of eigemodes is quite similar to those in Fig. 4B and 4E but with the absence of motion in the armchair boundaries. The eigenmode at $19.62$ kHz, Fig. 4F, is confined only on the zigzag edges, indicating the dominant contribution from the edge mode in the zigzag branch of the infinite graphene marked by the black dot in Fig. 3B.

Another intriguing property observed in Fig. 4 is that the extended modes in the region from $\sim 19.61$ kHz to $\sim 20.58$ kHz manifest a very interesting behavior. As shown in Fig. 4B, the translational components $u,v$ of the eigenmodes are spread in the whole finite structure but the rotation is localized only in the boundaries. In Fig. 4C, similar properties as the one in Fig. 4B are observed but now the rotational component is only confined in the zigzag boundaries as the rotation of the edge mode in Fig. 4F. This highlights a novel behavior of the dynamics in finite-size granular graphene, where one can find modes with extended translations in the structure while localized rotation on the boundaries. To the best of our knowledge such a behavior has not been reported before in other graphene structures, and the rich wave physics originates here from the extra rotational degree of freedom. As a result, the rotation of the beads can have a very distinct behavior compared to the translation of the beads in the MGG. This can lead to interesting potential applications like rotation isolation devices in a more general mechanical system and advanced wave control.

The existence of edge waves in the MGG can be confirmed directly in experiment. To observe the edge wave propagation, the experimental set-up is the same as the one shown in Fig. 1C, while a harmonic signal of duration $10$ ms with an initial linear ramping has been used as the source. All particles $B$ are still scanned to record the $u$, $v$ components. Figure 5 displays the measurements of total displacement amplitude ($\sqrt{u^{2}+v^{2}}$) for two separate times with a signal at $20$ kHz. As shown in Fig. 5A, when $t=1.7$ ms, the displacements are primarily localized on the zigzag edge while decaying into the bulk. Despite the excitation of bulk waves at this frequency, decay of the bulk wave is expected due to both dissipation and 2D geometrical spreading which provides a better observation of the elastic edge wave at $20$ kHz. The numerical simulation of the experimental process is shown in Fig. 5C, where the translational components of just the particles $B$ are shown. A good agreement between experiment and simulation is achieved.

V Novel turning effect of edge waves

We now turn our attention to the frequency range at which edge modes coexist in both zigzag and armchair boundaries. For example, since $20$ kHz is located in this frequency range, one should expect that when the zigzag edge wave of $20$ kHz reaches the corner, this wave can be mode-converted into an armchair edge wave. To observe this phenomenon, the spatial pattern at $t=3.7$ ms is depicted in Fig. 5B and 5D. Indeed, wave motions are seen to be localized on both the zigzag and armchair boundaries. This can be further confirmed by the close-up of the experimental spatial pattern of motion at the lower MGG corner as presented in Fig. 5E. To demonstrate this turning effect more clearly, we have chosen two particles marked by black and blue hexagons in Fig. 5E and we plot their time evolution in Fig. 5F. In addition, Fig. 5G provides the spatial distribution of edge waves, which are obtained from experiments and simulations by focusing on rows $\mathrm{n}=22$, $\mathrm{n}=32$ and column $\mathrm{m}=10$, as labeled in Fig. 5B by arrows. In rows $\mathrm{n}=22$ and $\mathrm{n}=32$ the motion distribution of the mode shows a similar profile (the amplitude is normalized to the first bead on the left zigzag edge), confirming the edge mode property as the motion propagates along the edge while decays very fast into the bulk (a distance of around $x=9d$). For $\mathrm{m}=10$, i.e. bottom panel of Fig. 5G, the translational signal also reveals a localized profile close to the armchair boundary confirming that the movement of beads in the armchair edge is due to the turning effect but not from the bulk modes. Note that, as indicated in Figs. 5B and 5D, due to dissipation and slow propagation velocity, the edge wave at $20$ kHz on the armchair edge is damped before propagating a long distance, e.g. $15d$. Further investigation of zigzag and armchair edge wave dynamics both with and without losses can be found in the SM. Note also that the spatial pattern in Fig. 5B shows a small asymmetry of wave propagation in the upward and downward directions in the experiments. This is most likely due to uncertainties in the pre-compression forces and asymmetric excitation of motion due to small misalignment between the driving bead and the set-up.

We now investigate the edge wave propagation along a corner of an angle of 120$\degree$, which connects a zigzag to another zigzag boundary. Considering that the experimental MGG is still of a finite size $21\times 41$, there is a free zigzag boundary at $(\mathrm{m},\mathrm{n})=(1,\mathrm{n})$. In order to build the second zigzag boundary, the particles above the line connecting the position $(\mathrm{m},\mathrm{n})=(1,11)$ and $(\mathrm{m},\mathrm{n})=(12,1)$ are removed, see Fig. 6.

The experimental set-up for zigzag to zigzag edge wave measurement is the same as that shown in Fig. 1C, where a harmonic signal of duration $10$ ms with an initial linear ramping and a frequency of $20$ kHz has been used as the source at position $(\mathrm{m},\mathrm{n})=(1,26)$. All particles $B$ are still scanned to record the $u$, $v$ components. Since edge modes are found on the zigzag edges at around $20$ kHz, one expects that when the zigzag edge wave reaches the 120$\degree$ corner, this wave can turn to the other zigzag boundary. To observe this phenomenon, Fig. 6A displays the measurements of total displacement amplitude ($\sqrt{u^{2}+v^{2}}$) at a given time after reaching the steady state. Numerically, the simulation of the experimental process is shown in Fig. 6B, where the translational components of only the particles $B$ are shown.

Indeed, wave motions are seen to be localized on the new zigzag edge and a good agreement between experiment and simulation is achieved. We can observe that bulk waves are also present because there is no gap for bulk mode at this frequency. In addition, Figure 6C provides greater detail for the spatial distribution of zigzag edge waves at a given point of time by focusing on rows $\mathrm{n}=18$, and cut line from $(\mathrm{m},\mathrm{n})=(5,8)$ to $(\mathrm{m},\mathrm{n})=(9,20)$, as labeled in Fig. 6(B) by arrows. As shown in Fig. 6C, the two profiles have a similar form with the translational signals becoming very weak (the amplitude is normalized to the first bead of each zigzag edge) after a distance of around $x=9d$. This confirms that bead movement in the new zigzag boundary is due to the turning effect from the zigzag to zigzag edges but not from the bulk modes.

VI conclusions

In this work, we propose a new artificial graphene called magneto granular graphene. This structure is composed of stainless steel beads in contact and placed in a properly designed magnetic field. The latter magnetizes the beads resulting in equivalent pre-compression forces between beads, and thus a mechanically stable structure, free of mechanical borders. The MGG proposed in this work can be used as a perfect experimental benchmark for fundamental study of Dirac and edge wave physics in mechanical systems.

Considering the wave behavior in the MGG, first we obtain the dispersion relation and the Dirac points. Then, we turn our attention to the edge physics of the structure. We show that the MGG supports unconventional edge waves that can exist also in armchair free boundaries, in contrast with genuine or other artificial graphene. In addition we show that for a range of frequencies, the structure supports edge vibrations both on the zigzag and armchair boundaries. Interesting enough, in this region the bulk modes are extended in their translation motions but localized at the edges regarding their rotation motion. Such a unique behavior has not been reported before to the best of our knowledge. Novel applications such as rotational isolators could be then designed using MGG or other flexible mechanical metamaterials with rotational elements Qian ; Babaee ; Deng .

Moreover, we also demonstrated that the coexistence of edge wave solutions in both zigzag and armchair boundaries lead to a turning effect from zigzag to armchair/zigzag free boundaries. This does not require a full bulk gap, which normally is necessary in the scenario of pseudospin topologically-protected wave propagation, like the case of helical edge waves. The role of the topology in MGG, like in the recent works of higher order topological insulators, is a remaining intriguing question and might lead to potential study of novel topological phase in mechanical systems. Finally, taking advantage of the intrinsic nonlinearities of the granular crystals, the MGG proposed herein offers a perfect platform to explore a wide array of novel nonlinear bulk, edge waves in mechanical graphene, similarly to solitons Nesterenko ; Chong , nonlinear waves CostePRE ; JobPRL ; LeonardExpMech ; CabaretPRL and breathers BoechlerPRL ; TheocharisPRE in simpler crystal structures.

Acknowledgements.

This work has been funded by RFI Le Mans Acoustique in the framework of the APAMAS and Sine City LMac projects and by the project CS.MICRO funded under the program Etoiles Montantes of the Region Pays de la Loire and partly funded by the Acoustic Hub project.

Appendix A Contact description in granular graphene

Considering the in-plane motion in the magneto-granular graphene (MGG), there are normal, shear and bending interactions characterized by contact rigidities $\xi_{n}$, $\xi_{s}$ and $\xi_{b}$, respectively, as represented in Fig. 1E. For stiffness of macroscopic elastic spheres in the MGG, the contact mechanism between the beads can be modeled by the Hertzian contact Johnson ; Mindlin . This leads to the expressions of the rigidities:

[TABLE]

As long as the physical parameters of the beads are known and the precompression $F_{0}$ is measured, $\xi_{n}$, $\xi_{s}$ and $\xi_{b}$ can be obtained. In our experiment, the stainless steel beads have the parameters: Young’s modulus $E=190$ GPa, Poisson’s ratio $\nu=0.3$, diameter $d=7.95$ mm, and density $\rho=7678$ kg/m3. The procompression can be determined around $F_{0}\sim 1.55$ N by means of measurement. This leads to the rigidities: $\xi_{n}\simeq 6.19\times 10^{6}$ N/m, $\xi_{s}\simeq 5.09\times 10^{6}$ N/m, and $\xi_{b}\simeq 3.04\times 10^{2}$ N/m.

Appendix B Calculation of Dirac points

Let us consider the modes at the corner (K point) of the BZ, by applying the periodic boundary condition, i.e., $\bm{U}_{\mathrm{m},\mathrm{n}}^{\alpha}=\bm{U}^{\alpha}e^{i\omega t-iq_{x}\mathrm{m}-iq_{y}\mathrm{n}}$, Eqs. (3) lead to two degenerate modes at the K point:

$\omega_{D_{\pm}}=\sqrt{\dfrac{g\pm\sqrt{g^{2}-h}}{4M}}$.

Above, $g=3[\xi_{n}+\xi_{s}+2\mathrm{P}(\xi_{b}+\xi_{s})]$, and $h=72\mathrm{P}(\xi_{n}\xi_{b}+\xi_{n}\xi_{s}+\xi_{b}\xi_{s})$ with $\mathrm{P}=MR^{2}/I$. These degenerate modes originate from the symmetry of honeycomb lattice, which are two Dirac points with frequencies $\omega_{D_{\pm}}$ at the K point.

Appendix C Dissipation

In order to compare the experimental results with the numerical simulations, the attenuation of the wave during the propagation has to be considered. In our model, the attenuation is implemented by a phenomenological linear viscous on-site dissipation BoechlerNat considering a time of decay, $\tau$, for elastic waves which can take different values as a function of displacement polarization. This leads to extra terms in the right hand side of Eqs (3),

[TABLE]

Therefore, wave dynamics of the MGG considering dissipation can be described by combining the boundary conditions in Eqs. (4) and Eqs. (5) with Eqs. (7). We can notice that Eqs. (7) are second order ordinary differential equations of time. As a consequence, we numerically obtain the time evolution of elastic wave propagation in the structure by solving Eqs. (7) using Runge Kutta fourth order method. More details about the time evolution of elastic wave propagation in two-dimensional granular crystals can be found in Ref.ZhengPRB . By fitting the experimental results with the numerical ones, we can estimate that $\tau$ is about 1 ms for both polarizations of displacement in our experiment.