Edge states of mechanical diamond and its topological origin

A mechanical diamond, a classical mechanics of a spring-mass model arrayed on a diamond lattice, is discussed topologically. Its frequency dispersion possesses an intrinsic nodal structure in the three-dimensional Brillouin zone (BZ) protected by the chiral symmetry. Topological changes of the line nodes are demonstrated associated with modification of the tension. The line nodes projected into two-dimensional BZ form loops which are characterized by the quantized Zak phases by 0 and $\pi$. With boundaries, edge states are discussed in relation to the Zak phases and winding numbers. It establishes a bulk-edge correspondence of the mechanical diamond.


Introduction
Topological semimetal is a system in which its band gap is finite almost everywhere in the Brillouin zone, except on some sets of isolated points. In short, it is a system with singular gapless points [1]. Generally, existence of a singular gapless point leads to nontrivial topology. That is, a singularity often serves as a source of a "twist" of Bloch wave functions captured by the Berry curvature, which gives rise to nontrivial topology. When a system is characterized by a point-like singularity associated with linear dispersion, i.e., a Dirac cone, it is Dirac/Weyl semimetal, which is a representative topological semimetal [2,3,4,5,6]. Very recently, the other kinds of the topological semimetal, nodal line semimetal where the gapless points form a line, typically a closed loop, begins to attract attention as a new stage to play with topology [7]. Just as in the case of fully gapped topological insulators, topological semimetals are characterized by topologically protected edge modes as a consequence of the bulk-edge correspondence [8].
On the other hand, there are intensive efforts to export the idea of topological insulators and topological semimetals to classical world [9,10,11,12,13]. Especially, classical mechanical systems are interesting playgrounds because of its simplicity and flexibility to tune parameters [11,12,13,14,15,16,17,18,19,20,21]. In fact, we have already seen topological edge modes in various mechanical systems [11,12,13,15,16,17,18]. One specific example to show the easy-totune feature is mechanical graphene, which is a honeycomb spring-mass model [22,16,18,23,24]. In mechanical graphene, Dirac cones are known to exist, and their number and positions in the Brillouin zone can be controlled by simply changing the tension of springs in equilibrium [23]. It is worth noting that the effect of the equilibrium tension can be interpreted as a spin-orbit coupling if clockwise and counterclockwise motion of the mass point is mapped to spin [25,26].
In this paper, we make an analysis on a spring-mass model with the diamond structure, namely, mechanical diamond, which is a natural extension of mechanical graphene to threedimension [23]. Mechanical diamond is a typical classical counterpart of a nodal line semimetal. It is revealed that the equilibrium tension induces unique evolution of the structure of the gapless line nodes. Topological properties of mechanical diamond are also investigated by relating edge modes, the quantized Zak phase, and the winding number. In accordance with the unique line node structure, the edge modes also shows unique distribution on the surface Brillouin zone, including the situation with the edge mode multiplicity 2 or 3. We confirm that these features are well captured by the quantized Zak phase and the winding number, which establishes the bulk-edge correspondence in mechanical diamond.
The paper is organized as follows. In Sec. 2, basic notions of mechanical diamond are explained. Then, the frequency dispersion is shown in Sec. 3. Section 4 is devoted to the topological arguments, and the paper is summarized in Sec. 5.

Formulation
Let us introduce our model, mechanical diamond, which consists of mass points aligned in the diamond structure and springs connecting the nearest neighbor pairs of the mass points [See Fig. 1(a)]. This three-dimensional model is a natural extension of two-dimensional mechanical graphene [23], a spring-mass model with honeycomb structure. As in the case of 2D mechanical graphene, parameters characterizing our model are mass of mass points m (m is fixed to unit for simplicity), spring constant κ, natural length of springs l 0 and distance between the neighboring mass points R 0 . Note that R 0 and l 0 does not necessarily match with each other, namely, R 0 can be larger than l 0 if we apply a proper boundary condition to exert uniform outward tension to the system. In order to investigate the dynamics of the system, we introduce a quantity x Ra = t (x Ra , y Ra , z Ra ) describing displacement of the each mass point from the equilibrium position. Here, R designates a lattice point and a is a sublattice index.  We assume that the elastic energy U s of a specific spring can be expressed as U s = 1 2 κ(l − l 0 ) 2 with l = |R 0 + x − x |, where x , x , and R 0 are two displacement vectors for two mass points, and a vector connecting the two equilibrium positions, respectively. By expanding U s up to the second order in δx = x − x , we obtain where γ µν R 0 = (1 − η)δ µν + ηR µ 0R ν 0 ,R 0 = R 0 /|R 0 | and η ≡ l 0 /R 0 . Here, we implicitly take the summation over µ and ν, which run through the directions of displacement x, y, and z, i.e., we apply the Einstein summation convention. The second term that is linear in δx has no  The third term is characterized by the parameter η. For η = 1, at which the springs have the natural length in equilibrium, the angle between R 0 and δx has significant impact on the elastic energy because of the factorR µ 0R ν 0 . On the other hand for η = 0, the limiting case of the stretched springs in equilibrium, U s has no dependence on the direction of δx, in contrast to the case of η = 1. Therefore, the parameter η enables us to tune the property of the spring-mass model.
For mechanical diamond, the equation of motion can be derived by evaluating the elastic energy substituting R (i) 0 in Fig. 1(b) into Eq. (1). If we further assume that the system is periodic in space and time, namely, by introducing φ aµ (k) as u µ ka = e iωt φ aµ (k) and where Γ µν . 1(a) and 1(b).] The eigenfrequency and the eigenmodes are calculated by diagonalizinĝ Γ(k). Note thatΓ (k) =Γ(k) − κ(4 − 8 3 η)1 anticommutes withΥ = diag(1, 1, 1, −1, −1, −1), i.e.,Γ (k) has the chiral symmetry. This chiral symmetry can be used to make topological characterization ofΓ(k), since the shift proportional to1 does not modify the eigenmodes. Figure 2 shows the frequency band dispersion of the system for several values of η. Here, κ is scaled by a factor 1/(1 − 2 3 η) to remove the η dependence of the total band width. Since we have six degrees of freedom, two from the sublattices and three from the directions x, y, and z, per a unit cell, we find six bands. As in the case of the single orbital tight-binding model with the diamond structure, the gap between the third and fourth bands is closed on the X-W line regardless of the value of η. (We employ the standard notation for the high symmetry points in the fcc Brillouin zone.) Around these gapless points, the gap grows linearly in the direction perpendicular to the X-W line. Actually, in the 3D Brillouin zone, the gapless points form line nodes, which are protected by the chiral symmetry. This point will be discussed in detail soon later. Interestingly, a new gapless point, which again forms a line node, is identified on the W-L line for η < 3/4, namely, we can generate a line node that is absent in the single orbital tight-binding model on the diamond lattice simply by applying tension to control η.

Frequency dispersion and nodal line
In order to obtain a global mapping of the geometry of the line nodes in the 3D Brillouin zone, we employ the Berry phase. This is because the direct assessment of degeneracy in the eigenfrequency spectrum requires some care, and it is safer to detect a "twist" in the eigenmodes associated with the singular degeneracy, which is captured by the Berry phase. In the numerical calculation of the Berry phase, we follow the idea in Ref. [27]. In the following, we describe the procedure to map the line nodes in order.
First, the Brillouin zone is decomposed into small cubes whose corners are specified by k l = (k l 1 , k l 2 , k l 3 ) where l = (l 1 , l 2 , l 3 ), k lµ = 2πl µ /N B , and l µ = 0, ..., N B − 1. We also introduce a triplet ψ(k) = (|n 1 (k) , |n 2 (k) , |n 3 (k) ), where |n i (k) is the eigenvector of the ith band at k. Here, "triplet" because we are focusing on the gap between the third and fourth bands. Then, we define a U (1) link variable by Here, we have introduced a shorthand notationê µ = (2π/N B )(δ 1µ , δ 2µ , δ 3µ ), and N ±êµ (k) ≡ | det ψ † (k)ψ(k ±ê µ N D )|. Now, we can assign the Berry phase to the each square surface of a small cube by taking the edge of the square as an integration path to define the Berry phase. For instance, for the square spanned by k l , k l +ê µ , k l +ê ν , and k l +ê µ +ê ν (µ = ν), the assigned Berry phase is computed as Arg Uê ν (k l +ê µ + a N Dê ν ) In our case, the Berry phaseυ µν (k l ) is quantized into 0 or π (modulo 2π), owing to the chiral symmetry. Furthermore,υ µν (k l ) = π implies that there exist odd number of Dirac points on the associated square [23,28,29]. In the 3D view point, a Dirac cone on a square patch means that a line node threads through that patch. Now, if we set N B large enough so that the possibility of multiple line nodes threading one square patch is excluded, we can follow each line node by following Berry phase π. Namely, as far as the square patches are sufficiently small, we can draw a reasonably smooth line node by connecting the central points of the nearby patches with Berry phase π. Note that we can always find a nearby patch with Berry phase π when we have a certain patch with Berry phase π, since a line node cannot be terminated abruptly.  Before moving onto the next topic, we explain why the gapless points in our model form a line, 1D object. Although our model is a six band model, here we focus on a two band effective model that is valid in the vicinity of the gap closing point. The most generic form of the effective Hamiltonian can be written as where R(k) is a three-dimensional vector consists of real quantities R 1 , R 2 and R 3 , and σ is a vector with the Pauli matrices as the components. In our case, the chiral symmetry induces a constraint on R(k). For simplicity, we choose σ 3 as a chiral operator, or assume {H, σ 3 } = 0, which restricts the hamiltonian to the form H(k) = R 1 (k)σ 1 + R 2 (k)σ 2 . The condition to have a gap closing point is |R(k)| = 0, which in fact gives two conditions R 1 (k) = 0 and R 2 (k) = 0. That is, we have two conditions and three parameters, k 1 , k 2 , and k 3 . Therefore, a manifold satisfying the zero gap condition has a dimension 3 − 2 = 1, which supports the existence of the line nodes [30]. Note that the choice of σ 3 as a chiral operator does not mean loss of generality since it is merely a matter of the basis choice. (a)

Zak phase and topological edge modes
In order to relate nontrivial topology induced by the bulk band singularity and edge modes, we calculate the Zak phase, whose definition is Here, A is a non-Abelian Berry connection A = ψ † dψ, which is a 3×3 matrix-valued one-form associated with a multiplet ψ = (|n 1 , |n 2 , |n 3 ). The integration path L is taken to be a periodic path of fixed (k 1 , k 2 ) and 0 ≤ k 3 < 2π, a line along k 3 -direction. Just as in the case of the Berry phase in the previous section, the Zak phase is quantized into 0 or π due to the chiral symmetry [31]. Figures 5(a)-5(e) show the Zak phase as a function of (k 1 , k 2 ), which is obtained using the same method as in the previous section by modifying the integration path.
Blue regions indicate υ(k 1 , k 2 ) = 0, while red regions indicate υ(k 1 , k 2 ) = π. The obtained results are consistent with the three-fold symmetry and the reflection symmetry of the diamond lattice. The lines dividing the blue and red regions are the line nodes (See Fig. 4) projected on to the surface under consideration, that is, the jumps in υ(k 1 , k 2 ) are associated with the bulk singularity. The complicated structures for intermediate η can be understood by relating each of the loop on the surface to each of the line node in the bulk. This point will be discussed later. If we consider a surface parallel to the plane spanned by a 1 and a 2 , on which k 1 and k 2 are good quantum numbers, we expect at least one topologically protected edge mode for (k 1 , k 2 ) with v(k 1 , k 2 ) = π. This statement can be confirmed by investigating edge mode explicitly. For that purpose, we consider a system periodic in the a 1 and a 2 direction, but finite in the a 3 direction, so as to have surface parallel to the a 1 -a 2 plane as in Fig. 6(a). Here, we apply a fixed boundary condition, where the last springs at the surface are connected to the wall. This choice of the fixed boundary condition is important to keep "the chiral symmetry" of the system [23]. Then, the frequency dispersion as a function of (k 1 , k 2 ) is obtained by diagonalizingΓ edge (k), whose explicit form iŝ whereΓ 1 (k) = −κ(γ 4 + e −ik·a 1γ 1 + e −ik·a 2γ 2 ),Ẑ = κ(γ 1 +γ 2 +γ 3 +γ 4 ) andΓ 2 = −κγ 3 . Importantly, due to the fixed boundary condition, the diagonal part is uniform, which means that the eigenvectors are unaffected by the diagonal terms and we still have a chiral symmetry in the same sense as the previous analysis. Figures 6(b)-6(e) show the frequency dispersion on a line in the surface Brillouin zone specified in Fig. 6(b), obtained using a system with 300 layers in z-direction. We find flat in-gap modes that are localized to the surface. As η reduces from 1 to 0, the location of the flat mode changes associated with the changes in the bulk spectrum. Notably, there is only one pair of edge modes aroundΓ point for η = 1, while there are three pairs of edge modes aroundB point for η = 0. This means that, if we remove the degeneracy originating from the existence of two surfaces (top and bottom), i.e., focusing on one of the surfaces, the flat edge mode for η = 1 is nondegenerate while the one for η = 0 is three-fold degenerate.
Let us relate the edge mode and the Zak phase. Figure 7(a) shows the Zak phase and the multiplicity of the pair of the edge modes for η = 0.5. As in the case of Fig. 5, the blue and red regions correspond υ(k 1 , k 2 ) = 0 and π, respectively. From this picture, it is found that even number of pairs exist for the region with υ(k 1 , k 2 ) = 0, while odd number of pairs for υ(k 1 , k 2 ) = π. It is relatively easy to understand the region with multiplicity 0 and 1. For those regions, an established argument on the relation between the Zak phase and the edge mode applies as it is. Then, the question is how to understand the region with multiplicity 2 and 3. Intuitively, it is understood by pulling the projected loops in the 2D surface Brillouin zone back into the loops in the 3D bulk Brillouin zone. As we can see from Fig. 7(b), even when two projected loops are overlapping and crossing with each other, there is an adiabatic way to resolve the overlap in 3D. Assuming that each loop carries one edge mode as in the standard argument, and any adiabatic change keeps the topological properties intact, the multiplicity larger than one can be attributed to the overlapping loops.
So far, we have been focusing on the Zak phase so as to emphasize the role of the band singularity as a source of the twist in the wave function, or eigenvector in our case. However, in order to capture the number of the edge modes beyond the parity of the number of the edge modes, we have to introduce the other topological invariant, winding number. The definition of the winding number is made possible by the chiral symmetry. Therefore, the constant diagonal term should be subtracted before applying the following argument to our model. For a system with the chiral symmetry, an appropriate choice of the basis set leads to the hamiltonian of the form Then, the winding number is evaluated as [32] where the path L is taken to be the same as the path used to define the Zak phase that is (k 1 , k 2 ) fixed and 0 ≤ k 3 < 2π. With this choice of the path, the winding number becomes a function of (k 1 , k 2 ). Then, it is expected that N w (k 1 , k 2 ) captures the number of edge modes at (k 1 , k 2 ). This expectation is confirmed by our numerical calculation of N w . Note that we have a relation υ = πN w (mod 2π), indicating that the Zak phase has an ability to capture the parity of the number of edge modes. Generically speaking, existence of more than one zero modes near the same boundary implies gap opening (deviation from zero energy) due to the interaction among them. However in the present model of nearest neighbor spring mass model, the edge states of each boundary have fixed chirality, χ, (Υ |ψ χ edge = χ |ψ χ edge , χ = ±) [33]. Then due to the selection rule, ψ ± edge | (Γ edge − κ(4 − 8 3 η)1) |ψ ± edge = 0, all of the multiple zero mode edge states remain at the zero energy. It justifies the winding number itself specifies the number of edge states.

Summary
To summarize, we have investigated the topological properties of mechanical diamond, a threedimensional spring-mass model with the diamond structure. We have shown interesting evolution of the gapless line nodes as a function of η, which is a parameter representing tension of the springs in equilibrium. The structure of line nodes is especially complicated for intermediate η.
In accordance with the complicated structure of the line nodes, the multiplicity of the edge modes show the complicated distribution in the surface Brillouin zone. We have also established the bulk-edge correspondence in mechanical diamond by relating the edge modes and two kinds of bulk topological number, the quantized Zak phase and the winding number, where the chiral symmetry plays an essential role.