Weyl points and topological surface states in a three-dimensional elastic lattice

Following the realization of Weyl semimetals in quantum electronic materials, classical wave analogues of Weyl materials have also been theorized and experimentally demonstrated in photonics and acoustics. Weyl points in elastic systems, however, have been a much more recent discovery. In this study, we report on the design of an elastic fully-continuum three-dimensional material that, while offering structural and load-bearing functionalities, is also capable of Weyl degeneracies and surface topologically-protected modes in a way completely analogous to the quantum mechanical counterpart. The topological characteristics of the lattice are obtained by \textit{ab initio} numerical calculations without employing any further simplifications. The results clearly characterize the topological structure of the Weyl points and are in full agreement with the expectations of surface topological modes. Finally, full field numerical simulations are used to confirm the existence of surface states and to illustrate their extreme robustness towards lattice disorder and defects.


I. INTRODUCTION
The search for novel topological states of matter has recently seen many exciting breakthroughs either in quantum or classical wave physics [1][2][3][4][5] . Phenomena such as the quantum Hall effect (QHE), the quantum spin Hall effect (QSHE), and the quantum valley Hall effect (QVHE) have been studied extensively and their nontrivial topological characteristics have been clearly connected to their ability to support backscattering-immune topological edge states 1,2 . Although the concept of topological material was discovered and developed in the field of condensed matter physics, the many conceptual and mathematical similarities between different fields of wave physics have led to the development of analogue concepts of topological materials in electromagnetic, acoustic, and elastic systems  . Recent studies have identified the existence of the so-called Weyl semimetals as three dimensional nontrivial topological materials capable of unidirectional backscattering-protected surface states 26 . Weyl semimetals are quantum materials whose behavior is characterized in terms of the Weyl Hamiltonian H(k) = ν x k x σ x + ν y k y σ y + ν z k z σ z , where ν i , k i , and σ i represent the components of the group velocity, momentum, and Pauli matrix, respectively. These materials can be thought of as a threedimensional extension of two-dimensional nontrivial topological materials characterized by Dirac degeneracies such as, for example, systems exhibiting QHE, QSHE, and QVHE. Similar to 2D Dirac materials, Weyl semimetals possess a degenerate nodal point formed by linear intersection of two bands in the three coordinate directions of the three dimensional reciprocal space. This degenerate point, called the Weyl point, carries a nonzero topological charge which confers the material nontrivial topological properties. The charge (i.e. the integral of the Berry curvature flux in a 2D manifold enclosing the Weyl point) can have either a positive or negative sign, hence determining if the degeneracy acts as a quantized source or a sink of berry curvature flux. The total topological charge of a Brillouin zone (BZ) in reciprocal space should always be zero 27 and therefore Weyl points always occur in pairs of opposite charge. Nontrivial topological surface states connecting the Weyl points of opposite charge exist on the boundary of the system. Apart from topological surface states, Weyl points have also been linked to other unusual effects such as chiral anomaly 28 and quantum anomalous Hall effect 29 .
Weyl points can be conceptually interpreted as the extension of Dirac points to three-dimensional lattices. However, some fundamental differences exist that have important implications on the robustness of these degeneracies. The existence of Dirac points can be guaranteed based on symmetry arguments (e.g. graphene-like lattices), provided that both parityinversion (P) and time-reversal (T ) symmetries are preserved. Breaking either P or T -symmetry would lift the degeneracy leaving behind a bandgap with topological significance. Unfortunately, the existence of Weyl points is much more ambiguous and it cannot be guaranteed a priori. However, some useful indications can come from the analysis of symmetries. In the case of Weyl points, P and T -symmetries impose contradicting requirements. In fact, considering a lattice having a Weyl point at the k high-symmetry point in reciprocal space, T -symmetry would require the topological charges at k and −k to be of the same sign while P-symmetry would require them to be opposite. Due to these contradicting requirements, it follows that Weyl points can only be obtained in systems with broken PT-symmetry. Thus, unlike for Dirac points in which symmetry breaking is a necessary condition to lift the degeneracy, for Weyl points symmetry breaking is a requirement for the existence of the isolated degeneracy. These more stringent requirements on symmetry conditions are at the foundation of the higher robustness of Weyl points against external perturbations. In fact, previous studies showed that they can only be annihilated by combining pair of points with opposite charges. Also, following the requirements dictated by both P and T -symmetries, the minimum number of Weyl points in a BZ is 4 (2 pairs) when P-symmetry is broken and 2 (1 pair) when T -symmetry is broken 3 .
Note also that, although in the current study only Weyl points with unit charge are considered, other studies [30][31][32][33][34][35][36][37][38] have reported systems possessing Weyl points having higher charge typically caused by the overlapping of multiple Weyl points of unit charge. Such points can be readily identified from the band structure because the modes do not intersect linearly at these points in the k x -k y plane. However, the intersection along k z would still be linear. Recent studies have also reported a new kind of Weyl point labeled type-II both in pho-tonic and phononic systems [39][40][41][42][43][44] . The degeneracies discussed in the present study are of type-I and represent the more direct analogue to the degenaracies originally studied by Weyl in quantum mechanical systems 45 .
Despite the successful implementation of Weyl points in photonic and acoustic systems, the realization of Weyl points in elastic systems has proven to be more challenging than other classical analogues. This increased complexity stemmed from wave coupling and mode conversion between different polarizations as well as from the difficulty in fabricating the necessary 3D structures. To-date, successful realizations were reported only in a few studies 35,61 . One study made use of beams and thin plates with a particular geometry 35 purposely selected to be manufactured by additive manufacturing. Another study utilized a truss-like construction made of beam elements 61 . Although elegant and simple to fabricate, these designs are discontinuous in nature therefore not well suited for applications where the load-bearing capabilities play a critical role (e.g. applications to structural materials). Also, in both the above mentioned studies, the topological characteristics were obtained based on a tight-binding (TB) formulation of the Hamiltonian. This approach works well only when the fundamental lattice is composed of weakly coupled resonant elements. In addition, the hopping parameters used in the TB model usually cannot be clearly connected to the actual geometric and design parameters. These limitations raise the more important question of how accurately TB models can describe continuous elastic lattices. In this regard, we propose an elastic analogue of a Weyl semimetal material based on a fully continuous design resulting in a solid load-bearing structure. Also, due to the strongly coupled nature of the elemnts in this design, we present a detailed study of the topological characteristics of such medium based on ab initio calculations, hence bypassing the limitations of TB approaches. This paper is organized as follows: Sec. II will introduce the proposed design and the corresponding dispersion properties in connection with Weyl points. Sec. III will present ab initio calculations of the topological invariant. Sec. IV will focus on the dispersion behavior of a supercell and the corresponding dispersion structure, in order to determine the occurrence of surface states. Finally, Sec. V will show full field simulations of the elastic Weyl material in order to illustrate the one-way, backscattering-immune nature of the surface states.

II. SYNTHESIS OF THE UNIT CELL AND DISPERSION PROPERTIES
From a general perspective, in order to develop the fundamental unit cell, we select a 2D triangular lattice (assumed in the xy-plane) and build the 3D geometry by periodically repeating this unit along the z-direction. The final result is a layered 3D structure having intact P and T -symmetries. Note that, as previously mentioned, starting from a 2D triangular lattice guarantees (due to symmetry arguments) the existence of Dirac degeneracy at the high symmetry points in the k x -k y momentum space. Periodically repeating this unit in the z-direction extends the Dirac degeneracy at virtually any k z value, thereby generating a line node degeneracy 3,55 along k z . In order to lift this three-dimensional degeneracy, P-symmetry can be broken by the proper introduction of additional structural element that do not respect inversion symmetry. For this purpose, we select slanted beam elements connecting the vertices of the hexagonal structure on two adjacent layers as shown in Fig. 3. The result is the introduction of a chiral coupling between the layers that breaks P-symmetry and all existing mirror symmetry. The most direct consequence is that the line degeneracy along k z is reduced to a point, hence giving rise to the Weyl point.
In the following section, we first describe the geometry and the dynamic properties of the parity-preserving fundamental unit. Then, we introduce the chiral coupling and investigate its effect compared with the P-symmetric unit.
A. 3D lattice with intact P-symmetry: line node degeneracy show the fundamental unit cell used to create the lattice. The cell is made of a vertical cylinder included between two homogeneous, hexagonal-shaped, thin plates. The cylinder has its longitudinal axis aligned with the z-direction. The cylinder is made of iron while the plates are made of aluminum. The selection of these materials was motivated by the intent of obtaining marked degeneracies at the high symmetry points in the momentum space. Clearly, several other combinations of structural materials leading to similar conditions could potentially be identified. The dimensions of the unit cell are: a = 30 mm, h 0 = 2 mm, r 1 = 5 mm, h 1 = 20 mm, where a is the lattice constant, h 0 is the thickness of each homogeneous flat plate, r 1 and h 1 are the radius and height of the cylinder, respectively. A 2D lattice in the x-y plane can be assembled by periodically repeating the unit cell. The final result is a lattice of cylinders in a hexagonal configuration sandwiched between two thin plates. This 2D lattice can be periodically repeated in the z-direction to form the final 3D lattice.
The analysis of the band structure of the 3D material can be conveniently performed by applying periodic boundary con-ditions to the fundamental unit cell. More specifically, periodic boundary conditions can be applied along the sideward faces of the plate in order to create the 2D lattice in the xyplane, while periodic boundary conditions applied on the top and bottom faces of the top and bottom plates will yield the complete 3D lattice. The resulting system is characterized by intact P and T -symmetries and the corresponding 3D Brillouin zone is shown in Fig. 1 (c). For a fixed value of the momentum component k z , the 3D BZ simplifies into a 2D BZ typical of a triangular lattice.
The dispersion curves along the boundary of the 2D BZ in the k x -k y plane at k z = 0, and in the k x -k z plane at k y = 0 have been calculated using the commercial finite element (FE) package COMSOL Multiphysics and are shown in Fig. 2. As visible in Fig. 2 (a), four modes intersect linearly at a frequency of approximately 40 kHz and give rise to two degenerate points at K. These degenerate points are indicated by the red and green boxes in the inset of Fig. 2 (a).
The resulting degenerate nodes are Dirac points which are protected by the lattice symmetry of the hexagonal (graphenelike) arrangement of vertical cylinders. As previously mentioned, this degeneracy exists at all k z , hence resulting in the line degeneracies marked in red and green colors in Fig. 2 (b). This behavior is a direct consequence of both P and Tsymmetries being preserved. Breaking P-symmetry would result in lifting the line degeneracy along the k z direction hence resulting in the formation of a Weyl points.

B. P-symmetry breaking and Weyl points
In order to break P-symmetry, a chiral coupling can be introduced by means of straight slanted cylinders connecting top and bottom layers in proximity of the vertices of the hexagonal layer. The resulting unit cell is shown in Fig. 3 while its 3D BZ remains the same as the one shown in Fig. 1 (c). The slanted cylinders are made of aluminum (as the plates) and have a radius of r 2 = 2 mm. The dispersion curves along the boundary of 2D BZ in the k x -k y plane at k z = 0 and in the k x -k z plane at k y = 0 are calculated again and shown in Fig. 4. The degenerate points resulting from the linear intersection of two modes in Fig. 2 (a) are also present in Fig. 4 (a) (marked by the same red and green boxes). However, the line degeneracies previously observed in Fig. 2 (b) are now lifted, due to P-breaking, and reduced to nodal degeneracies at the high-symmetry points K and H (see red and green circles Fig. 4 (b)). These degenaracies are the Weyl points that are formed by the linear intersection of the two (initially degenerate) modes along all three directions. Fig. 4 (a) shows the linear intersection in the k x -k y plane while Fig. 4 (b) shows the linear intersection in k x -k z plane.
As a result, the degeneracy marked in Fig. 4 (a) is lifted for nonzero values of k z (with k z = nπ/a z for integer n) hence giving rise to a topological bandgap as shown in Fig. 4 (c). The band structure shown in Fig. 4 (c) corresponds to k z = ± 0.1π a z , where a z = h 1 + 2h 0 = 24 mm is the total height of unit cell. The subscripts ± of the high symmetry points' labels indicate k z = ± 0.1π a z , as also depicted in Fig. 1 (c). The bandgap size varies with k z as shown in Fig. 4 (b). Thus, k z can be treated as a parameter that controls either the opening or closing of the topologically nontrivial bandgap, provided that the system remains periodic along the z-direction.

III. AB INITIO CALCULATIONS OF THE TOPOLOGICAL PROPERTIES
A. Topological charge In order to assess the topological significance of the degeneracies identified in the chiral 3D lattice and, consequently, the existence of Weyl points, this section presents a numerical investigation into the calculation of the topological charge. The charge can be calculated by integrating the Berry curvature flux on a 2D manifold S enclosing the Weyl point in the reciprocal space, that is where C is the topological monopole charge, i.e., the Chern number, dS is the vector surface element aligned with its local normal direction, and Ω Ω Ω is the Berry curvature (vector field) given by where k is the wavevector, and u n (k) is the displacement eigenstate as a function of k. Consider the Weyl point at K, at a frequency of approximately 40 kHz, as indicated by the red dashed box in Fig. 4 (a). This point corresponds to the intersection of the 16 th and the 17 th bands. To calculate the corresponding topological charge, we can integrate the Berry curvature flux threading the outer surface area (on the top, bottom, and side faces) of a fictitious prismatic manifold enclosing the K point (see Fig. 5). We anticipate, and show later, that the flux on the side faces is zero if a certain triangular prism is chosen as prototypical 2D manifold. Without loss of generality, we can select k z = ± 0.1π a z as the top and bottom faces of the prism. At these planes, we have K + and K − with the same (k x , k y ) coordinate as the K point. The z-component of the Berry curvature on the top plane is then calculated using Eq. 2 (see Fig. 6 (a) and (b)) for the 16 th and 17 th bands, respectively. Note that, at the two distinct points (i.e. they cannot be transformed into one another by translation symmetry operations), K + and K + points exhibit the same Berry curvature pattern. This is a direct result of the C 6 symmetry of the lattice. Therefore, if the 2D manifold is chosen to be the triangular prism with the right red dashed triangle in Fig. 6 (a) as its top face, we can assure that the flux integral on the side faces of the prism is zero due to symmetry in the reciprocal space. In fact, the C 6 symmetry guarantees that the outward flux on the side face x equals the inward flux on the side face y; the same cyclic symmetry along with the translation symmetry (from the left to the right, for example) forces the flux integral on the entire side surface z to vanish.
It follows that we only need to integrate the Berry curvature flux on the top and bottom faces. The bottom face which  is around the K − point, is actually the time reversed counterpart of the triangular face around the K + point (inversion of (k x , k y , k z )) as shown in the left red dashed triangle in Fig. 6 (a). T -symmetry requires that Ω Ω Ω(−k) = −Ω Ω Ω(k) which results in the z-component of the Berry curvature on the bottom face to have the same strength but opposite sign. Also, given that the dS element on the bottom face is oriented in the −k z direction, it still contributes a positive outward flux. As a result, the topological charge is the sum of the integrals of the Berry curvature calculated over the two red dashed triangles in Fig. 6 (a). While solving for the Bloch eigenmodes, we performed a parametric analysis across the k x and k y grid points by using an adaptive resolution (finer near K + point). The Berry curvature was calculated based on Eq. 2 using a finite difference formulation on the two triangular regions, and the results are shown in Fig. 6 Fig. 4 (a) is now lifted, due to the nonzero value of k z , hence leaving behind a topological bandgap. The degeneracy marked by the green box in Fig. 4 (a) is also lifted due to the nonzero value of k z , which opens a second (although narrower) bandgap. with opposite charge. Fig. 7 shows all the Weyl points within the half BZ (in the k z direction) and their corresponding topological monopole nature (either sink or source). On the k z = 0 plane, there are two Weyl points (K and K ) both with charge +1 (that is sources of Berry curvature), while on the k z = π/a z plane, each of the two Weyl points (H and H ) are of charge −1 (that is sinks of Berry curvature).

B. Bandgap-Chern number
An alternative way to obtain the topological charge of the Weyl point is to use a topological invariant called the bandgap-Chern number C g . This invariant is obtained by summing the Chern numbers of all the bands below the gap of interest. In the following, we consider k z as a parameter and focus on the bandgap that forms in the k x -k y plane. C g remains constant with varying k z except at K and H where the bandgap closes and reopens and C g experiences a sudden change in its value. This change is equal to the topological charge of the Weyl points on that plane because Weyl points act as either sources or sinks of Berry curvature flux 56 . This principle can be used to calculate the charge of the Weyl points in the band structure. Particularly, we focus on the degeneracies occurring at approximately 40 kHz. It is found that the Berry curvature integrals of the first 15 bands cancel each others, therefore C g reduces to the integral of Berry curvature in the 2D BZ for the 16 th mode only. Figure 8 shows the plot of C g as a function of k z for the Weyl points at 40kHz. It clearly shows that C g has a value of −1 for k z < 0 and +1 for k z > 0 which results in a net difference ∆C g = +2. This is consistent with the previous analysis that showed the k z = 0 plane to host two Weyl points of equal charge at K and K due to intact T -symmetry. The two Weyl points at K and K therefore have a topological charge of +1 each. In the same figure, there is also a change in the value of C g at k z = ± π a z . This indicates that the charge of the Weyl point at H (H ) is −1. Once again, this is expected and consistent with the fact that Weyl points always occur in pairs with opposite topological charge. The nonzero value of charge also confirms that the nodal degenerate points under consideration are in fact Weyl points with nontrivial topological significance. In light of the bulk-surface correspondence 2 , a nonzero C g indicates that topological surface modes can exist within the target bandgap while the magnitude of C g indicates how many surface modes should be expected. In our specific example, C g = ±1 for k z ∈ (0, π/a) or k z ∈ (−π/a, 0), respectively. Compared with trivial materials (such as vacuum), there is a k z -locked net difference of ±1 in C g . Therefore, one surface mode should be expected either in case of positive or negative k z . The corresponding surface states propagate unidirectionally either in the clockwise or counter-clockwise direction along the structure border (depending on the sign of C g ), as shown later on in the numerical simulations. This result further confirms that the selected geometry leads to the formation of Weyl points and it is capable of supporting topological surface states.
We note that the nonzero Chern number and the Berry curvature distribution for, as an example, k z = 0.1π/a z (Fig. 6) is reminiscent of 2D topological materials manifesting quantum Hall effect (QHE) 62 . In QHE materials, T -symmetry is broken due to the use of an external field. In the present material, the chiral structure breaks the z-mirror symmetry therefore a  8. The bandgap-Chern number C g associated with the gap around 40 kHz as a function of k z . As the bandgap closes and reopens at k z = 0, C g experiences a +2 discontinuous jump in its value, which indicates the sum of topological monopole charges on the k z = 0 plane. mode propagating in the z−direction in the infinite lattice is expected to exhibit k z -locked unidirectional surface states. In fact, at steady state, the role of temporal and the spatial (in this case z) variables can be interchanged. A similar idea was also presented in studies of 3D Floquet topological insulators 63 . This latter aspect will be further clarified in the next section.

IV. SUPERCELL DISPERSION AND SURFACE MODES
Weyl points with opposite topological charges are connected by nontrivial topological surface states that exist only on the boundary of the medium. In this section, we consider the physical response of a periodic supercell made out of the 3D unit cell described above. Consider a supercell having a finite dimension of 38 units in the y direction and infinite dimensions (obtained by means of periodic boundary conditions) in the x and z directions. Fig. 9 (a) shows the actual domain used for the numerical calculations.
We focus the following analysis on the topological bandgap that develops at the approximate frequency of 40kHz and that is marked by the red box in Fig. 4 (c). Fig. 9 (b) shows the dispersion curves for the supercell at k z = 0.1π a z . It is evident from the visual inspection of the dispersion that the supercell admits two additional modes, the surface modes, that exist within the topological bandgap. Bulk modes are in grey color while the surface modes are indicated in green and red. The green modes correspond to the lower edge of the supercell while red modes correspond to upper edge. This correspondence is also reflected in the plot of the eigenvectors on the supercell shown in Fig. 9 (c) and 9 (d). The dispersion of the surface modes corresponding to the upper edge has a positive slope (i.e. positive group velocity), which means that these modes are expected to travel in the positive x-direction. Similarly, the surface modes defined on the lower edge have negative slope and therefore are expected to travel in the negative x-direction.
Further, it can be observed that the supercell in Fig. 9 is  Fig. 4 (c). Topological surface modes exist in the bandgap and are indicated in red and green color. The bulk modes are indicated in grey color. The green colored modes occur at the lower edge of the supercell while the red colored modes occur at the upper edge as indicated in Fig. 9 (c) and Fig. 9 (d), respectively. The eigenvectors corresponding to the points indicated by (c) the green star, and (d) the red star in the dispersion curves and plotted on the supercell. obtained by terminating the triangular lattice along the zigzag edge. However, triangular lattices exhibit also another edge type named the armchair edge. We explored the possible propagation and existence of topological modes also on this edge type. For this purpose, a finite size supercell terminated along the x-direction was developed, as shown in Fig. 10 (a). The supercell extended indefinitely in the y and z directions (periodic boundary conditions were used to simulate the infinite size). The corresponding band structure at k z = 0.1π a z is shown in Fig. 10 (b). The bulk modes are marked in grey color, the surface modes corresponding to the right edge of the supercell are indicated in green color, and those corresponding to the left edge are indicated in red color. The eigenvectors corresponding to the surface states on both the right and left edges are shown in Fig. 10 (c) and 10 (d), respectively. Similar to the zigzag edge, the dispersion of the surface modes corresponding to the left edge has a positive slope (i.e. positive group velocity), which means that these modes are expected to travel in the positive y-direction. The surface modes defined on the right edge have negative slope and therefore are expected to travel in the negative y-direction. Therefore, if we assemble an extended 2D domain in the plane x-y we would expect to see the surface mode to propagate in a clockwise direction along the boundary of the lattice.
For k z = − 0.1π a z , the dispersion curves remain unaltered but the direction of propagation of the surface modes is re-versed. The surface mode corresponding to the upper edge acquires a negative group velocity (propagation in the negative x-direction) while the surface mode corresponding to the lower edge acquires a positive group velocity (propagation in the positive x-direction). Similarly, the surface mode corresponding to the left edge acquires a negative group velocity (propagation in the negative y-direction) while the surface mode corresponding to the right edge acquires a positive group velocity (propagation in the positive y-direction). An extended 2D domain in the x-y plane would then show a wave propagating in the anti-clockwise sense along the boundaries.
The combined behavior illustrated above is well consistent with those of dynamical systems characterized by Weyl points and corresponding topologically nontrivial modes. Results also suggest that, upon injecting a wave into the system, the direction of propagation can be controlled by properly choosing the sign of k z . This latter aspect will be further clarified by means of full field simulations in the next section.
Another interesting feature of the proposed geometry is that the group velocity of the surface modes can be controlled simply by changing the radius of the center vertical cylinder r 1 . Fig. 11 shows the dispersion curve of the supercell (a) for various values of r 1 . The group velocity of the surface modes decreases with increasing radius r 1 of the center cylinder. This reduction of the group velocity can be explained by observing that by increasing r 1 the homogenized mass density of the  Fig. 4 (c). Topological surface modes exist in the bandgap and are indicated in red and green color, while the bulk modes are indicated in grey color. The green colored modes occur at the right edge of supercell while red colored modes occur at the left edge of the supercell as indicated in Fig. 10 (c) and Fig. 10 (d) respectively. The eigenvectors corresponding to the points indicated by (c) the green star (d) the red star in the dispersion curves and plotted on the supercell. lattice increases, but without significantly increasing the rigidity. Hence, given that the eigenmode does not involve significant compressive or bending motion of the center cylinder, the surface modes tend to slow down. By tuning the radius r 1 , the surface states can gradually evolve from the single valley mechanism (similar to the 2D model discussed by Raghu 62 ) to the intervalley mechanisms (similar to the model discussed in Kane 64 ). In either case, the unidirectional gapless surface states connecting the lower and upper bulk bands are protected by the topological nature of the band structure and are guaranteed to be present. Note that, while a similar variation in the group velocity of the topological edge modes was earlier observed in 2D chiral materials 63 , this tuning ability of the surface states in 3D lattice was never observed and documented in earlier studies.

V. FULL FIELD NUMERICAL SIMULATIONS
Based on the 3D unit cell developed above, a complete domain in the x-y plane was also simulated in order to observe both the direction of propagation and the scattering behavior in presence of defects. The domain used is shown in Fig. 12 (a). The domain is rectangular and finite along both the xand y-directions. Both edges were set to free boundaries. The boundaries normal to the z-axis were assigned periodic boundary conditions so to render the lattice infinite in this direction. Perfectly matched layers (PML) were used at the left edge to absorb the incoming wave and allow a clear assessment of the direction of propagation of the surface state when performing a steady state analysis (i.e. it avoids the surface mode to propagate all around the boundary and get back to the source).
The domain was excited by a harmonic point force acting in the z-direction at the location indicated by the red star marker in Fig. 12 (a). A harmonic excitation at a frequency of 40 kHz was selected because right within the target topological bandgap where surface states exist. At k z = 0.1π a z , the source excites the surface state with positive group velocity which is observed propagating in the clockwise direction (see Fig. 12 (b)). At k z = − 0.1π a z , the situation is inverted and the source excites an counter-clockwise propagating surface state (Fig. 12  (c)). These results confirm previous observations concerning the one-way propagation property of the surface states.
These results also allow an additional observation on the back scattering immune properties of the surface states. As the wave propagates beyond the corner of the rectangular domain, the wave amplitude on both sides is exactly comparable. This suggests that there are no significant reflections taking place at the sharp corner and the wave is capable of propagating unidirectionally.
It should be noted that while the states propagate with the same intensity along the xand y-aligned edges, they are associated with different eigenstates. This can be clearly observed in Fig. 12 (b) and 12 (c), and it is due to the different structure of the edges (i.e. zigzag versus armchair).
The immunity to back scattering was also assessed by introducing a defect in the rectangular domain represented by a FIG. 11. Dispersion curves of the supercell shown in Fig. 9(a) for different values of the radius of the center cylinder. Values of r 1 equal to (a) 4mm (b) 5mm (c) 6mm (d) 7mm (e) 8mm (e) 9mm are considered. It is seen that the group velocity of the topological surface modes decreases with increasing r 1 and that the surface states evolve from a single-valley to an inter-valley dominated behavior. cutout on the top edge, as shown in Fig. 13 (a). The source was located on the left edge as indicated by the star marker. The excitation was kept identical to the previous simulation scenario. By exciting at k z = 0.1π a z , a clockwise propagating surface state is clearly obtained as shown in Fig. 13 (b). As the state encounters the defect on the top edge, it propagates around its perimeter and it continues along the original direction. Observing that the intensity of the wave amplitudes both before and after the defect are exactly comparable, and that there are no waves traveling back towards the source (reaching the section of the edge before the source), one can conclude that no appreciable back-scattering takes place due to the defect. The latter is a clear hallmark of topologically protected states.

VI. CONCLUSIONS
Weyl points and topological surface modes in elastic solid media are still in their early stages of formulation and design. The very few implementations proposed to-date are based on discrete designs that are not conducive to use in practical applications, particularly in those requiring load-bearing mate-rials. This study presented a fully continuous, load-bearing, elastic system capable of unidirectionally-propagating and topologically-protected surface states. The design of the 3D unit cell consisted in a layered prismatic lattice with hexagonal cross section in which the layers were spaced by solid cylindrical elements and by slanted circular beams connecting consecutive faces of the prismatic unit cell. The cylinders were used to set the necessary in-plane symmetry conditions and to provide high load carrying capacity. The slanted beams were used to achieve chirality and allowed achieving the necessary P-symmetry breaking conditions. Note that the layered structure resulted in a very feasible and practical design that has the potential to greatly facilitate the fabrication phase. The analysis of the lattice dynamics highlighted the existence of Weyl points following the breaking of the z-mirrorsymmetry and the P-symmetry of the lattice. To gain insight into the mechanism leading to the formation of these degeneracy points, we evaluated the topological invariants using ab initio calculations and without employing any further simplifications. Results were very well aligned with the expected quantized values. The relationship between the surface states and the topological invariants was also elaborated upon from different perspectives. k z -locked unidirectional propagating surface elastic states were predicted on the external surface of the 3D lattice and their existence was further confirmed by full field numerical simulations. Numerical results also confirmed the extreme robustness of these states against strong lattice defects. This study may serve as a basis to develop structures having surface-elastic-wave-guiding capabilities for applications in fields such as vibration control, energy harvesting, structural health monitoring, and on-chip telecommunication signal processing.