Tilted double Dirac cone and anisotropic quantum-spin-Hall topological insulator in mechanical granular graphene

Dirac degeneracies are essential ingredients to control topological charge exchanges between bands and trigger the unique edge transport properties of topological materials. In addition, when Dirac cones are tilted, exotic phenomena can emerge such as anomalous Hall effect or unconventional Klein tunneling. However, the unique topological transport properties arising from the opening of tilted Dirac cone degeneracies have been left completely uncharted. Here, we demonstrate a new form of Dirac degeneracy that occurs in mechanical granular graphene (MGG): a tilted double Dirac cone, composed of two counter-tilted type-I Dirac cones. Different from the reported C6 systems, we show that the tilted double Dirac cone is present in a C2 granular graphene. Remarkably, a pair of anisotropic helical edge waves appears when the degeneracy is lifted. This leads to an anisotropic quantum spin-Hall topological insulator that possesses unique wave propagation properties, including anisotropic edge dispersion and direction-dependent edge-bulk mode conversion.


Introduction
In topological materials, the existence of edge/surface states is induced by the topology of bulk bands [1][2][3][4][5]. A non-trivial topology is typically acquired through the lifting of a point degeneracy between bands, such as the one stemming from a linear band crossing (Dirac cone) in momentum space. The crossing point is termed Dirac point (or Weyl point in three-dimensional spin-non-degenerate systems) [6][7][8]. The presence of Dirac cones, as commonly seen in graphene or graphene-like structures, is protected by time-reversal symmetry (TRS) and inversion symmetry. Breaking the TRS leads to quantum Hall effect [9] while breaking the inversion symmetry results in quantum valley Hall phase [10]. Moreover, keeping the TRS and inversion symmetry, Dirac cones can be lifted by the intrinsic spin-orbit coupling (SOC). In this case, the topological materials, called Z 2 topological insulators, exhibit the so-called quantum spin Hall (QSH) phase [11][12][13] and robust helical edge states. Similar phenomena have also been investigated in classical photonic/phononic/mechanical systems. There, effective spin degeneracy has been emulated using double Dirac cones (overlap of two Dirac cones), allowing for the observation of classical analogs to QSH effects [14][15][16][17][18][19][20][21].
Dirac cones can be tilted due to, for example, mechanical deformation [22], external magnetic fields [23] or anisotropic hoppings [24]. In this case, the zero-energy Fermi surface can take different shapes, depending on the strength of the tilt [25][26][27]. In type-I tilted Dirac cone, the strength of the tilt is relatively weak, thus the Fermi surface remains a point (the Dirac point). Examples include the hydrogenated graphene [28], organic semiconductor [29] and type-I Weyl semimetals [30]. When the strength of the tilt is strong, a type-II Dirac cone emerges with the Fermi surface turning into two lines crossing at the Dirac point. This dispersion has been reported in type-II Dirac and Weyl semimetals [31][32][33][34], screw symmetry photonic crystals [35], graphyne-like lattices [36] and metasurfaces [37]. Recently, a type-III Dirac cone is demonstrated in different systems [38][39][40], where the Fermi surface is a single line at the critical transition from a point-like (type-I) to a lines-crossing (type-II) Fermi surface. The successful implementation of different kinds of Dirac cones could open the door for advanced control of waves [41][42][43]. In particular, manipulation of tilted Dirac cones could bring drastically different transport phenomena. For example, in type-I tilted Weyl metals with broken TRS, anomalous Hall effect can be observed [44,45]. Moreover, in some Dirac materials, the anisotropic property due to tilted energy dispersions leads to the unconventional Klein tunneling [46,47].
However, the physical properties of tilted double Dirac cones, i.e. a fourfold degeneracy obtained at the band touching points between two counter-tilted (opposite tilting angle) Dirac cones, have been left completely unexplored. This is due to the difficulty in obtaining such degeneracies. The requirement of extra hopping/deformation lowers the symmetry of structures, typically preventing these complex forms of symmetry-protected Dirac cones to occur in low symmetry lattices [48]. Yet, tilted double Dirac cones may bring novel wave physics and an unique form of Dirac/Weyl semimetals. Exploring this new area rises fundamental questions that must be answered. For example, how to obtain a tilted double Dirac cone in a simple lattice? Once the tilted Dirac degeneracy is lifted by SOC, would a novel form of QSH effect appear?
And if yes what is the influence of the tilt on the properties of the helical edge states?
In this work, we answer the above questions, demonstrating a novel form of anisotropic QSH effect. We first demonstrate how to form a tilted double Dirac cone in a modified mechanical granular graphene (MGG) under out-of-plane motion excitation [figure 1(a)]. One of the key ingredients is that the MGG, in which each bead possesses one translational and two rotational degrees of freedom (DOFs) [49][50][51][52], provides sufficient number of modes for a fourfold degeneracy. This is the key for a tilted double Dirac cone in a C 2 symmetry lattice, which is completely distinguished from the band-folding technique on C 6 structures [19,53]. Another key ingredient is the existence of rotations which play a significant role as providing two pairs of rotational degenerate modes at the Brillouin zone (BZ) center. In the presence of two in-plane rotations, a double Dirac cone can be achieved by tuning the rigidities. Interestingly, when modifications are introduced which lower the symmetry of the MGG from C 6 to C 2 , the existence of the double Dirac dispersion is not affected due to the accidental nature of the degeneracy. However, the modifications lead to a novel type of Dirac dispersion: a tilted double Dirac cone. We show that lifting the Dirac degeneracy can introduce a Z 2 band gap, which supports a pair of topological helical edge modes. We further demonstrate the transport properties of the anisotropic spin-locked helical edge waves. Our results as well as the recent experimental studies confirm the relevance of MGG as a reliable platform to study novel Dirac and topological wave physics reported here [54,55].

Model
The MGG, a two-dimensional (2D) monolayer granular crystal, is shown in figure 1(a). Considering the out-of-plane motion, each bead exhibits one out-of-plane translation u and two in-plane rotations ϕ, φ (with respect to the x-and y-axes, respectively) i.e., figure 1(b). The movements of individual particles lead to the following forces and/or interactions: (1) shear forces with the rigidity ξ s (figure 1(c)); (2) torsional moments with the rigidity ξ t (figure 1(d)); (3) bending moments with the rigidity ξ b (figure 1(e)). For simplicity, we set ξ t = ξ b , and denote them as η b = η t = ε e , where η b,t = ξ b,t /ξ s are the ratios of bending/torsional rigidities to shear rigidity. Each unit cell (highlighted by black dash box) can be viewed as a molecule with internal-couplings ε i (red bonds) within the molecule, and with external-couplings ε e (blue bonds) to its neighbors. By analyzing the motions between beads related to the unit cell and applying the Bloch wave solution, a dynamical equation of the MGG is given by (the full derivation can be found in references [54,55]), where with Φ = Rϕ, Ψ = Rφ (R is the radius of the beads). * denotes the complex conjugate. Ω = ω M/ξ s is the normalized frequency, ω is the cyclic frequency. D 1 and D 2 are 3 × 3 matrices of the following forms, Above, l = (1 − 4ε e )/2, and α = P e −2iq x , β = P e iq x cos q y , γ = √ 3P e iq x sin q y . q x = k x R, q y = √ 3k y R are the normalized wave vectors for k x , k y , respectively. P = MR 2 /I with I the moment of inertia. Normally, P ∈ [1.5 2.5] corresponding to a hollow sphere where all the mass is at the sphere periphery (P = 1.5) and a homogeneously filled sphere (P = 2.5), respectively. To this end, we set P = 1.55.
Equation (1) leads to the Ω − k dispersion as shown in figure 2(a) when ε i = ε e = 0.5. It can be seen that a double Dirac cone appears at the Dirac frequency Ω D = √ 3P at the Γ point (the BZ center). The four degenerate modes at the Γ point are purely rotational modes. When ε i = ε e , the graphene structure holds a C 6v symmetry. An accidental degeneracy of the wave functions corresponding to E 1 and E 2 irreducible representations of the C 6v group occurs when 2ε e + ε i = 3/2, leading to a double Dirac cone of purely rotational modes at the Γ point. Interestingly, even when the symmetry is reduced to a C 2 by setting ε i = 0.6, ε e = 0.45 [figure 2(b)], the double Dirac cone is preserved at the Γ point since 2ε e + ε i = 3/2 is still valid. Although the irreducible representations of the C 2 group are one-dimensional (symmetry or asymmetry), the tilted double Dirac cone is a degeneracy of twofold symmetric and asymmetric modes, see appendix A. In the MGG, the number of pure rotational modes is twice the number of the irreducible representations of C 2 , the extra modes of the same nature (rotations) would double the number of wave functions of the same irreducible representations, leading to twofold symmetric and asymmetric modes in the C 2 MGG at the Γ point. By keeping 2ε e + ε i = 3/2, accidental degeneracy of the twofold modes occurs again, forming the double Dirac cone. However, the influence of symmetry reduction to C 2 is to tilt each individual Dirac cone in opposite angle with respect to the Dirac point, leading to the tilted double Dirac cone. The full gap by tuning ε e and ε i [ε e = ε i , figure 2(c)] does not break C 2 symmetry, but destroys the accidental degeneracy, giving rise to a Z 2 topological phase transition of bulk bands accompanying with the band inversion as discussed appendix A.

Tilted double Dirac cone
To unveil the bulk properties of the tilted double Dirac cone and derive the effective Hamiltonian in the vicinity of the Dirac point (δk = (δk x , δk y )), we follow an asymptotic model developed in reference [56]  (see appendix B). The following dynamical equation is obtained, where are the Pauli matrices and σ 0 is the identity matrix.
denotes the 'effective mass' due to SOC that breaks the degeneracy [54]. m s = 4/3 * ε V D δk x is a symmetry-breaking mass with the modified parameter ε = ε i − ε e . For nonzero δk, it suggests that m s relies on the propagating direction. The four components spinor ψ = [p; d] consists two circular polarized where  The tilt parameter defined asṽ ≡ determines the tilted type of Dirac cone [23]. According to In this work, we considerṽ < 1, namely the tilted Dirac cones belong to type-I category. However, by tuning ε , other types of tilted double Dirac cone can also be achieved, see appendix C. Due to the tilting effect, wave behavior around the Dirac point exhibits anisotropic properties. To confirm that the anisotropic properties stem from the tilt, we further simplify the Hamiltonian in equation (5) Above, V ± = (1 ± 4/3 * ε cos θ)V D is the anisotropic Dirac velocity, and θ ∈ (−π, π] is the angle of wave vector to the x-axis. The tilted Hamiltonian in equation (6) shows that the slope of the linear dispersion (V ± ) around the Dirac point is anisotropic, depending on the propagation angle θ. When m = 0, the elegance of H θ ± to capture the wave physics around the Dirac point is shown in figure 3(f), where the cross section of the counter-tilted double Dirac cone [obtained from equation (1)] at Ω = 2.16 is presented in black curves. The red and green dashed ellipses correspond to the prediction of H θ + and H θ − , respectively. The perfect match of the two colored ellipses with the black curves verifies our theory, that is, the Dirac cone of H θ + (red ellipse) is tilted to the positive x-direction, while H θ − (green ellipse) to the negative, resulting in the anisotropy around the Dirac point.
It should be mentioned that the Hamiltonian in equation (6) can be rewritten as H θ ± = V D δk ± · σ ± mσ z where δk ± = δk + A ± with A ± = ±4/3 * ε δk cos θ. This indicates that the symmetry reduction to C 2 equally introduces a pseudomagnetic field, defined as B ± = ∇ × A ± as discussed in references [57][58][59]. However, this reduction does not break TRS because the total pseudomagnetic field of the counter-tilted double Dirac cones, B = B + + B − , is zero. See the discussion of TRS in appendix D. Under TRS, the nonzero m due to SOC induces a Z 2 topological phase, leading to a nontrivial topological gap around the Dirac point as discussed in appendix A. Thus, equation (6) is an anisotropic QSH Hamiltonian, formed by two copies of tilted Dirac fermion Hamiltonians H θ + and H θ − . Consequently, equation (6) predicts the existence of a pair of time-reversal edge modes (helical edge waves) with exotic transport properties.

Anisotropic quantum spin-Hall effect
The anisotropic QSH Hamiltonian H θ ± supports the unidirectional propagation of topological rotational waves along different interface configurations. As shown in the inserts of figures 4(a) and (b), a vertical zigzag interface (VZI) with θ = −π/2 and an inclined zigzag interface (IZI) with θ = −π/6 are constructed. We set cyan bonds e = 0.5 for both sides, while the red bonds ε i = 0.6 (ε = 0.1, gray beads), and the black bonds ε i = 0.4 (ε = −0.1, green beads). On the interface, the bonds is equal to 0.5. We also scale down the mass of beads in green by a factor f = 0.92 to ensure that a full gap appears around Ω D . Figures 4(a) and (b) show the edge dispersion. The bulk modes are labeled in gray, and the edge modes are marked by red and blue. More details about the calculation of edge mode spectra can be found in appendix E. It can be seen that the two zigzag interfaces support a pairs of topological edge waves, since the edge waves (red and blue lines) are spin-locked, they support unidirectional transport along the interfaces. As examples, the spatiotemporal evolution of edge waves (propagation of rotation-dominated waves) on the VZI are implemented in figures 4(c) and (d). An interface unit cell is used to be a harmonic source (green star) to excite the eigenmodes in the red or blue branches. The color scale represents the total magnitude of rotation (|Φ| + |Ψ|) in each bead. It shows that the edge wave transport is either upward [figure 4(c)] or downward [ figure 4(d)], manifesting the one-way propagation property. The unidirectionality of edge transport along the IZI is demonstrated in figure 4(e), where the helical edge waves propagate solely to the θ = −π/6 direction.
As can be seen from figures 4(a) and (b), there is a difference in the edge dispersions along the VZI and the IZI direction. This is translated into distinct wave propagation along these directions. It is reasonable to relate this direction-dependent edge transport with the resulting anisotropic property of the bulk waves around the Dirac point when the condition 2ε e + ε i = 3/2 is kept. Once the tilted double Dirac cone is lifted, the edge waves, that appear in the bulk gap, also exhibit such a direction-dependent transport. Thus, we observe an interesting bulk-edge correspondence that 'transfers' the anisotropy character of the underlying bulk waves to the resulting edge waves.
Remarkably, the direction-dependent property of edge wave transport leads to the topological edge-bulk mode conversion in the MGG. To demonstrate, an interface path is constructed by combining the VZI and IZI [inserts of figures 4(a) and (b)]. Thus, a corner appears to connect the VZI and IZI as shown in figure 5. When edge waves are excited from the source (stars), the turning of waves on the corner depends highly on the frequency location of different mode types. As shown in figure 5(a), there are three cases of mode types corresponding to the three colored areas in figures 4(a) and (b). The spatiotemporal evolution of edge waves for the three cases are shown in figures 5(b)-(d), respectively. In case I where topological edge modes do not exist in the IZI [ figure 5(b)], the edge waves of the VZI, which are topologically protected, can only transfer to the bulk modes when meet the corner. Therefore, the topological edge waves spread into the bulk and reflect back when bulk waves come to the free edge on the bottom, leading to zero transport of waves along the IZI. In case II [ figure 5(c)], since all the three types of modes can exist in the MGG, the edge waves of the VZI can transfer into the bulk and the IZI at the corner, resulting in the spreading of waves into both the bulk and the IZI. In case III [figure 5(d)], due to the fact that the MGG does not support the propagation of bulk modes, i.e., frequencies in the band gap, the edge waves in the VZI transfer only to the IZI at the corner. Figure A1. Dispersion curves and the corresponding rotational modes around the tilted Dirac cone for 2ε e + ε i < 3/2 (left panel) and for 2ε e + ε i > 3/2 (right panel).

Conclusions
In this work, we report a tilted double Dirac cone in a C 2 MGG under out-of-plane motion excitation. We found that the bulk dynamics around the Dirac cone is described by an effective anisotropic QSH Hamiltonian, which in contrast with the conventional case, is formed by two copies of the counter-tilted Dirac fermion Hamiltonians. The tilting effect results in a pair of helical edge waves with rich anisotropic transport properties. Our results can stimulate the study of novel Dirac/Weyl materials, and far-reaching applications for beam-splitters, switches and filters.
(Ω D + δΩ) 2 ≈ Ω 2 D + 2Ω D δΩ, we arrive at the following equation, by ignoring the second and higher order corrections, By left multiplying the projector X † to equation (A4), finally we can get the following equation, where the effective Hamiltonian H is given by, This circularly polarized spinor basis defines the pseudo-spin in the granular graphene. Under the new basis ψ, equation (A5) reads, where H 0 = V D (k x σ x + k y σ y ) is the Dirac Hamiltonian and the spin basis, The subscripts L, R represent the left or right circular polarization defined by the unitary matrix U (i.e. equation (A7)) acting on the linear basis μ. The tilting of the Dirac cones can be investigated by applying the following unitary matrix, Equation (A8) becomes, with s = U ψ and H ± = ±V A δk x σ 0 + V D (δk x σ x + δk y σ y ) the general form of type-I tilted Dirac Hamiltonian. When m = 0, H t describes two counter-tilted Dirac cones H ± overlapping at the Dirac point. By applying the unitary matrix, where θ satisfies δk x = |δk| cos θ, δk y = |δk| sin θ, H t can be rewritten as follow, or in a compact form,

Appendix C. Type-II and type-III double Dirac cones
According to the Hamiltonian H θ ± , the tilt parameterṽ = 4|ε |/3 determines the tilted type of the double Dirac cone. Figure A2 shows the type-II and type-III double Dirac cones by tuning ε .

Appendix D. Pseudo-TRS and Berry curvatures
In a spin-1/2 system, the QSH topological phase is protected by TRS where T 2 = −1 guarantees the existence of the Kramers' pair. And a T-invariant Bloch Hamiltonian must satisfy, Due to its bosonic nature where T 2 = 1, the analogies of QSH phase in spinless systems are not protected by the real TRS but a pseudo-TRS is constructed usually by other symmetries, e.g., C 6 in many reported photonic/phononic systems. Next, we show that the effective Hamiltonian of the MGG around the Γ point satisfies the same property as equation (A14) under the C 2 symmetry. Let us consider the rotation axis along z-axis passing through the midpoint of the center distance of A and B beads. The C 2 operation can be described by a rotation matrix (we ignore z component as its contribution to the Dirac physics around the Γ point is negligible), In addition, sublattices A and B change positions after a 180 • rotation, Thus, a C 2 operator is expressed as C π = S ⊗ R π . It is easy to prove that the Hamiltonian in equation (A6) satisfies, Compared equations (A14) and (A17), it leads to an important conclusion that the C π acting on the effective Hamiltonian of the MGG has the same property as T on the Bloch Hamiltonian of a spin-1/2 system. Under the C 2 rotation, the Hamiltonian is invariant. In this way, pseudo-TRS T p can be constructed based on the C π as long as it holds T 2 P = −1 and commutes with C π . Therefore, one pseudo-TRS is constructed as T P = −iσ 0 ⊗ σ y K (K is complex conjugation). It is possible to check that T 2 P = −1 and [C π , T P ] = 0. One should also notice that the C 2 rotation holds (δk x , δk y ) → (−δk x , −δk y ), H(δk x , δk y ) → H(−δk x , −δk y ) with the rotating axis must passing through the Γ point of the BZ center. Thus, the discussion here is valid around the Dirac point at the high symmetry Γ point.
The Berry curvature around the Dirac point can be calculated using the eigenmodes of H θ ± . Taking H θ + as an example, its eigenvalues are λ ± = ± V 2 + (δk 2 x + δk 2 y ) + m 2 , and the corresponding eigenmodes can The definition of the Berry connection is given by, A = i v ± |∇ k |v ± , and the Berry curvature by B = ∇ k × A. Therefore, based on equation (A18), we can numerically calculated the Berry curvatures around the Dirac frequency. In figure A3, we show the Berry curvatures of the four branches around the Γ point. The Chern number is defined as the integral of the Berry curvature of the same spin over the Brillouin zone. Since our granular system exhibits pseudo-spin and pseudo-time-reversal symmetry, then the pseudo-spins up and down have independent Chern integers N + and N − , which are −1, +1 for the two-fold degeneracy of lower cone, and +1, −1 for the two-fold degeneracy of upper cone based on the Berry curvatures in figure A3. Due to the real TRS, the overall Chern number (N = N + + N − ) is always zero for each two-fold degeneracy. However, the spin Chern number, defined as N s = N + − N − , is nonzero due to pseudo-time-reversal symmetry. This leads to a Z 2 invariant in the granular graphene.

Appendix E. Edge wave spectra
As shown in figure A4, we consider a strip consisting of 60 cells with the zigzag interface in the middle. The direction perpendicular to the strip is treated as periodic. Since each particle has three DOFs, we can write down in total 3 × 2 × 60 equations of motion including the interface configuration. Combining the boundary condition of free edges on the two ends of the strip, the spectra including bulk, interface and free edge modes can be derived, see the bottom panel of figure A4. In main text, figures 4(a) and (b) are the projections of bulk modes to the VZI and IZI based on the bulk Hamiltonian. To be more specific, all the bulk modes at each point of the BZ can be obtained from the bulk Hamiltonian, and they are projected to the directions along VZI or IZI when we focus on the waves along the interfaces. The edge wave dispersions for the two interfaces are obtained by further solving the boundary condition of the interfaces. By comparing the results of the VZI in figure 4(a) with the one in figure A4, it is clear that they are consistent to each other. There, we highlight the interface branches as blue and red, while we mark the bulk region in grey.