Manipulation of Dirac Cones in Mechanical Graphene

Recently, quantum Hall state analogs in classical mechanics attract much attention from topological points of view. Topology is not only for mathematicians but also quite useful in a quantum world. Further it even governs the Newton’s law of motion. One of the advantages of classical systems over solid state materials is its clear controllability. Here we investigate mechanical graphene, which is a spring-mass model with the honeycomb structure as a typical mechanical model with nontrivial topological phenomena. The vibration spectrum of mechanical graphene is characterized by Dirac cones serving as sources of topological nontriviality. We find that the spectrum has dramatic dependence on the spring tension at equilibrium as a natural control parameter, i.e., creation and annihilation of the Dirac particles are realized as the tension increases. Just by rotating the system, the manipulated Dirac particles lead to topological transition, i.e., a jump of the “Chern number” occurs associated with flipping of propagating direction of chiral edge modes. This is a bulk-edge correspondence governed by the Newton’s law. A simple observation that in-gap edge modes exist only at the fixed boundary, but not at the free one, is attributed to the symmetry protection of topological phases.

Scientific RepoRts | 5:18107 | DOI: 10.1038/srep18107 = R l 0 0 . When > R l 0 0 , each spring gives finite force even in equilibrium, which means that we should make the system to avoid shrinkage as a whole. For instance, we need to support the system by appropriate choice of the boundary condition. Dynamical variables for this model are deviations of mass points from the equilibrium positions, written as where R denotes a lattice point, and a is a sublattice index.

Results
Formulation. In mechanical graphene, the elastic energy cost for finite x Ra has an special importance. Here, we see it for a pair of mass points, introducing parameters R 0 , ′ x , and ″ x as in Fig. 1b. Throughout this paper, we adopt an assumption that all springs have good linearity, i.e., the energy of a spring is always written as . Note that if η = 1, the system is not stretched, while if η < 1, the system is stretched. Summation over indices appearing as a pair is assumed µ ν ( , = , ) x y . This formula indicates that U s strongly depends on the angle between R 0 and δx for η = 1, while U s is independent of the direction of δx for η = 0. For η = 1, no force is given by a spring if δx is normal to R 0 , since such motion only leads to the rotation of the spring without elastic energy cost. On the other hand, for η = 0, there is finite force even for δ = x 0, and then the motion with δ ⊥ x R 0 can increase the energy since there is a finite component of force normal to R 0 . In short, η controls the relation between the restoring forces for δx R 0 and δ ⊥ x R 0 , or in other words, the relation between the longitudinal and the transverse wave modes.
The physical meaning of the last statement is illustrated in Fig. 1d. There, we consider the case that a mass point is supported by two linearly arranged springs. For η = 1, the illustrated sidewise shift does not lead to harmonic oscillation since the potential energy cost for this shift is not quadratic but higher order in the shift size. On the other hand, for η<1, the sidewise shift leads to harmonic oscillation, since the restoring force given as a partial component of the force by the springs is linear in the shift size. So, the motion in perpendicular to the spring direction, which corresponds to a transverse mode, is significantly affected by η. In contrast, the motion in parallel with the spring direction, which corresponds to a longitudinal mode, is not affected by η as far as good linearity of the springs is assumed, since with good linearity, the elastic energy has a constant second derivative with respect to the spring deformation.
Having this in mind and assuming the periodicity of the system, we introduce new dynamical variables µ u ka and in order to obtain normal modes. Then the equation to be solved reduces to  Fig. 1c.] Practically, the ω-spectrum is obtained by diagonalizing Γ k . It is worth noting that except the constant diagonal elements, Γ k has no term connecting the same kind of sublattices. Since the constant diagonal elements only add a constant contribution to ω 2 , we say that Γ k has a "chiral" symmetry, since Γ k anticommute with ϒ = ( , , − , − ) diag 1 1 1 1 if the constant diagonal part is subtracted. In fermionic systems, the chiral symmetry plays important roles in various situations. For instance, it stabilizes Dirac cones in 2D cases 29 . Now, the problem is quite similar to the phonon problem in graphene [30][31][32][33][34] . However, it is difficult to control η in a wide range in real graphene. Furthermore, the effective model for graphene phonon does not respect the "chiral symmetry", which plays important roles in topological characterization. This is because if the vibration in graphene is modeled by a spring-mass model, it requires springs connecting the next nearest neighbor sites and more.
Manipulation of Dirac Cones. Now, we follow the η dependence of the ω-spectrum [ Fig. 1e-h]. Throughout this paper, the spring constant is scaled as κ η = /( − / ) 1 1 2 in order to make the total band width constant. For η = 1 (not shown), the second and third bands of ω 2 -dispersion is same as the one in the nearest neighbor tight-binding model for graphene (NN-graphene) in which Dirac cones are there at K-and K'-points, but the first and the fourth bands are exactly flat and stick to the top and bottom of the other bands. (See Ref. 18.) Note that Γ k at η = 1 is essentially same as the Hamiltonian of p-orbital honeycomb lattice model 35 . The ω-dispersion for η = / 5 6 is shown in Fig. 1e. It is similar to the one for η = 1 and the Dirac cone is still there at the K-and K'-points, but the first and the fourth bands are no longer flat. At η = / 2 3 [ Fig. 1f], the gap at the M-point between the second and the third bands closes. In the vicinity of the M-point, the dispersion is linear in Γ -M direction but quadratic in M-K direction. Such a linear vs quadratic dispersion is characteristic for merging of a Dirac cone pair 24 . In fact, soon after η becomes less than 2/3, two new Dirac cones other than those at the K-and K'-points appear on M-K lines. [See Fig. 1g.] Due to the three fold rotational symmetry, there are six new Dirac cones in the whole Brillouin zone. The new Dirac cones move from the M-point to the K-point as η approaches to zero. For small η where the new Dirac cones get close to the K-and K'-points, the dispersion relation looks like the one for bilayer graphene with trigonal warping 36 . Finally, at η = 0 [ Fig. 1h], the new Dirac cones are absorbed to the original Dirac cones at the K-and K'-points, and the dispersion becomes like the one of NN-graphene with extra double degeneracy. In η = 0 limit, there is no distinction between the longitudinal and transverse waves, and then, oscillations in x-and y-direction decouple, which leads to the two-fold degeneracy.
The transition at η = / 2 3 is characterized by parity of the "Herring number", which is defined for electronic systems as 37 where d is spatial dimension, k r s are the time reversal invariant momenta, and ( ) ± k N r is the number of occupied states with parity eigenvalue ± 1 at k r . This number is also used to obtain Z 2 topological number of the inversion symmetric topological insulator, or to check the existence of Dirac cones 38 . For 2D cases, parity of N H determines the parity of the number of Dirac cones in half of the Brillouin zone as explained below. In Ref. 39, in order to detect Dirac cones in generic 2D systems, we have used the Berry phase defined as where k and ⊥ k are two momenta in independent directions, and ⊥ u nk k is a Bloch wave function 39,40 . When the system respects both of the time reversal and spatial inversion symmetry, θ ( ) k is quantized into 0 or π for the spinless case 39 . As far as the quantization is kept, θ ( ) k has to show a jump when it changes as a function of k , but the jump should be associated with a singularity of the band structure that is nothing more than a Dirac cone. Having this in mind, we can say that if θ θ π π ( ) ≠ ( ) 0 m od 2 θ θ π π ( ( ) = ( ) ) 0 m od 2 , odd (even) number of Dirac cones exist in the region π < < k 0 , i.e., half of the Brillouin zone. On the other hand, the inversion symmetry gives us a relation 38,41 . Therefore, the Herring number has an ability to detect the number of Dirac cones. This idea has been applied to the 2D organic Dirac fermion systems 41,42 . It is also possible to apply the idea to "mechanical" graphene. We Edge Modes and Symmetry Protection. Next, we see edge states using ribbon geometry with zigzag edges. In this case, the equation to be solved is similar to equation (2), but with larger matrices, since the Fourier transformation is applicable only in one direction. Here, we adapt the fixed boundary condition, that is, the mass points at the boundary are connected to the springs fixed to the wall. [See Fig. 2f,g.] The obtained frequency spectra for several values of η as functions of k , momentum parallel to the edge, are shown in Fig. 2a-d. For η = / 5 6, we observe in-gap edge modes near = k 0 as in the case of η = 1 18 . The edge modes forms a flat band, since the "chiral symmetry" is preserved at the fixed boundary. At η = 0, the system exactly inherits the property of NN-graphene, i.e., the edge modes are found near π = k instead of = k 0. In short, the position of the edge modes in the edge Brillouin zone is switched between near = k 0 and π = k as η changes. Namely the new Dirac cones emerged from the M-points bring new edge modes near π = k and take away the edge modes near = k 0. The difference in momenta along the edge is reflected in the oscillation pattern in the real space. In order to convince this expectation we perform a calculation using a rectangular system with long zigzag edges and short armchair edges. In Fig. 3, the typical eigenmodes whose frequency is close to 3 , a frequency of the bulk Dirac points, is shown for η = / 1 3 and η = 1. For η = 1, the mass points at the boundary basically oscillate in phase with neighboring mass points on the edge with large wave-length background, reflecting the fact that the edge modes are existing near = k 0 in the momentum space. On the other hand, for η = / 1 3, the neighboring mass points on the edge tend to oscillate anti-phase, reflecting the edge mode position in the momentum space. For η = 1, the mass points at the boundary oscillate approximately normal to the boundary, while for η = / 1 3, oscillation in parallel to the boundary is also allowed. In Fig. 3, only the typical eigenmodes are shown, but we have confirmed that the above statements are basically valid for the mode near ω = 3. Therefore, the oscillation pattern of the eigenmode localized at the boundary can be used to detect the phase transition associated with the Dirac cone merging experimentally.
In order to see the importance of the boundary condition, the ω-spectrum for η = 1 obtained with the free boundary condition, where the springs at the boundary are simply removed, is shown in Fig. 2e. (Note that η ≠ 1 does not go well with the free boundary condition since the finite tension results in the deformation at the boundary.) We find no edge mode in the gap. This difference between the two boundary conditions is deeply related to the symmetry protection of topological phases, which is one of the most important concept in modern theory for topological matters. In the real graphene case, it is possible to assign a topological origin of the edge modes 43 . There, the chiral symmetry plays an important role in defining the bulk topological number, the quantized Berry phase, and preserving in-gap edge states. The fixed boundary respects the "chiral symmetry" of Γ k . On the other hand, the free boundary does not respects it since the absence of the springs at the boundary leads to nonuniform diagonal elements in Γ k . Therefore, even the bulk topological property is same, the topological edge states are not necessarily preserved at the free boundary in contrast to the fixed boundary, which reflects the manifestation of the symmetry protection.

Time Reversal Symmetry Breaking and Chiral Edge Modes.
Lastly, we break the time reversal symmetry by uniform rotation as a whole with angular frequency Ω , which brings two new forces, the centrifugal and the Coriolis force. For simplicity, only the Coriolis force, which is first order in Ω , is taken account of since the centrifugal force, which is second order in Ω , is negligible for small Ω . The equation of motion becomes 18 .
It is no longer possible to reduce the problem to diagonalization of Γ k . Instead, the ω-spectrum is obtained by ] 0 k 2 analytically, and solving a quartic equation of ω 2 . After obtaining ω, φ k n is derived as a nontrivial solution for equation (8). For η = 1, the finite Ω induces a gap at the Dirac point 18 . By solving equation (8), it is confirmed that such a gap opening remains at work also for η ≠ 1. At some critical point η η = c , the gap is closed at the M-point. Then, for η less than η c , the system is again fully gapped except η = 0 where the gap is again absent. η c depends on Ω , but it is close to 2/3 in the small Ω limit.
For the topological analysis of this gap closing and opening, we define the Chern number C as The summation is taken over = , n 1 2, since we are now focusing on the gap between the second and third bands. C can be evaluated with a well established numerical method 44 . With Ω = .
0 02, a representative value for small enough Ω , = C 1 for η = 1, which is consistent with the previous work 18 . As η is reduced, the Chern number Before closing, we make a comment on the effects of dissipation. Assuming that the frictional force is proportional to the velocity, the effects of friction are easily incorporated by slightly modifying Â . Then, although the friction prevents the long distance propagation of the chiral edge modes, we can still observe the unidirectional nature of the chiral edge modes as far as the friction is enough small. Very recently, an interesting aspect of Dirac cones in photonic system with radiational leakage is reported 45 . Then, addressing the relation between radiational loss in photonic systems and frictional dissipation in mechanical systems is an interesting future subject.

Discussion
To summarize, we first propose a simple and feasible way, uniform stretch of the system, to control the vibration spectrum of mechanical graphene. An interesting point of the proposed procedure is it does not involve symmetry breaking, but we observe generation, migration, and creation of Dirac cones. It is also shown that the transitions associated with the motion of Dirac cones are detected by observing real space pattern of the edge modes. For the case without the time reversal symmetry, an intimate relation between the bulk topological number and chiral edge modes in a mechanical system is established, which indicates universality of the bulk-edge correspondence, which is one of the most important concepts in the topological approach. In addition, the boundary condition dependence of the edge modes is explained in terms of the symmetry protection of the topological state, which reveals new aspect of the symmetry in classical mechanical systems.