Tunable Acoustic Valley-Hall Edge States in Reconfigurable Phononic Elastic Waveguides
This study investigates the occurrence of acoustic topological edge states in a 2D phononic elastic waveguide due to a phenomenon that is the acoustic analogue of the quantum valley Hall effect. We show that a topological transition takes place between two lattices having broken space inversion symmetry due to the application of a tunable strain field. This condition leads to the formation of gapless edge states at the domain walls, as further illustrated by the analysis of the bulk-edge correspondence and of the associated topological invariants. Although time reversal symmetry is still intact in these systems, the edge states are topologically protected when inter-valley mixing is either weak or negligible. Interestingly, topological edge states can also be triggered at the boundary of a single domain if boundary conditions are properly selected. We also show that the static modulation of the strain field allows tuning the response of the material between the different supported edge states.
Topological acoustics has rapidly emerged as a new and fascinating branch of physical acoustics. Following the groundbreaking research in solid state physics, which showed the existence of topological states of matter Hasan and Kane (2010); Ren et al. (2016), researchers have been able to formulate the acoustic analogue of selected mechanisms Yang et al. (2015); Süsstrunk and Huber (2015); Wang et al. (2015); Khanikaev et al. (2015); Mousavi et al. (2015); Nasha et al. (2015); Lu et al. (2016a); Chen and Wu (2016); He et al. (2016); Lu et al. (2016b, b). During the last decade, several studies have investigated topological materials based on broken space-inversion symmetry (SIS). In electronic materials, SIS was broken using different methodologies including graphene-like lattices with staggered sublattice potentials Xiao et al. (2007); Semenoff et al. (2008); Yao et al. (2009); Drummond et al. (2012); Kim et al. (2014), strained graphene Guinea et al. (2009), and multi-layered graphene under electric fields Jung et al. (2011); Vaezi et al. (2013); Ju et al. (2015); Zhang et al. (2013); Martin et al. (2008); Qiao et al. (2011); Li et al. (2016). The effect of the asymmetry is that of opening a gap at the original Dirac cone associated with the hexagonal lattice in which edge states are now supported. These edge states cannot be explained by the quantum Hall effect (QHE) mechanism. In fact, as TRS is intact the lattice still possesses a trivial topology within the context of QHE Hasan and Kane (2010); Raghu and Haldane (2008); Ochiai and Onoda (2009). However, due to the large separation in -space of the two valleys, valley-dependent topological invariants can be defined and used to classify the topological states of the different lattices. This approach, usually referred to as quantum valley Hall effect (QVHE), was recently investigated for application to fluidic acoustic waveguides Lu et al. (2016a, b), as well as elastic plates with local resonators Pal and Ruzzene (2017); Vila et al. (2017).
In the present study, we consider the dispersion and propagation behavior of a topological elastic phononic waveguide assembled from a truss-like unit cell (Fig. 1(a)). Compared to the previous studies Pal and Ruzzene (2017); Vila et al. (2017), we take a fully continuum modeling approach which provides a general methodology of analysis and allows mapping the topological behavior all the way back to the massless (or massive when SIS is broken) Dirac equation. Also, an in-depth study of the occurrence of edge states, either at domain walls or at the lattice boundaries, is presented. The analysis addresses both the topological significance of the different edge states as well as their lossless behavior. Further, despite the weak topological nature of the QVHE, we show that selective valley injection can be obtained by a dual-excitation approach able to target specific eigenstates. In addition, we show tunable and reconfigurable capabilities of the topological medium. Starting from an initial hexagonal lattice symmetry, the symmetry level can be tuned using an external static actuation which allows real-time tuning of the edge states (i.e. of the corresponding topological bandgap) as well as a complete reconfiguration of the lattice from a regular phononic material to a tunable topological medium, and viceversa. Note that, in this study, we are mostly concerned with the analysis of the physical behavior of the medium and of its tuning capability and we do not concentrate on the details of implementation of the tuning mechanism. Nevertheless, we note that the local perturbation could be practically implemented by using, as an example, mechanical or thermal loads. In the following, we assume that an isostatic pressure could be applied on the side walls to produce the deformation necessary for the topological phase transition.
The individual trusses forming the lattice are assumed made out of aluminum, and having width , and thickness where is the lattice constant. The fundamental lattice structure clearly exhibits SIS with inversion axes located at the nodal points. The lattice is then deformed by applying an internal pressure that lowers the symmetry to . The equivalent force produced by the is locally normal to the individual truss element as shown in Fig. 1(d). The application of either a positive or a negative pressure results in two different deformed states of the lattice that for simplicity we label as - and -states (Fig. 1(b),(c)), respectively. In these states, the SIS is broken and the states and are inverted images of each other, that is they turn into each other as switches sign.
We will use a combination of theoretical and numerical tools to show in the following that 1) the evolution from - to -state (and viceversa) is accompanied by a topological phase transition, and that 2) when a phononic system is obtained by assembling the two phases, an edge state is supported along the domain wall (DW, i.e. the interface between the - and -state). We note that these edge states can be created and modulated at the boundary of a single lattice by selecting proper boundary conditions and by tuning the internal pressure.
Ii Topological band structure analysis
ii.1 Phononic band structure
For a unit cell in prestrained conditions, the numerical analysis can be divided in two parts: 1) a nonlinear static analysis to calculate the state of displacement , stress , and strain induced by the applied pressure load, and 2) a linear wave propagation analysis around the pre-stressed equilibrium state.
In our analysis, the pre-stressed state was calculated using finite element analysis and used as input for the linear wave analysis. The wave propagation in the pre-stressed medium can be studied as a linear small oscillations problem around the new static equilibrium . The total displacement can then be expressed as . Under these conditions the linearized wave equation is given by
which is readily rewritten as:
where is the mass density, is the first Piola-Kirchhoff stress tensor, is the second Piola-Kirchhoff stress tensor with , where is a order elasticity tensor, and are the initial stress and strain, respectively, and the differential operator is taken with respect to the material frame. Solving the -dependent Bloch eigenvalue problem yields the phononic band structure and the eigenstates. The above described model was assembled and solved using the commercial finite element software COMSOL Multiphysics.
The frequency dispersion is normalized by the bulk shear wave speed in aluminum over the lattice constant . Note that the current system is effectively a flat waveguide that admits Lamb guided modes (symmetric and antisymmetric ), as well as shear horizontal modes. In this study, we concentrate only on modes (flexural modes) that are those more naturally excited in plate-like structures under external excitation. Other modes are filtered out based on particle motion polarization of the mode shapes. The value of the pressure perturbation was chosen in order to achieve large deflections without inducing plastic deformations. Note that, because the lattice was designed based on slender members, it can easily accumulate large deflections while maintaining small strain levels.
Fig. 2 (a) and (b) show the band structure of the flexural modes of the lattice with (intact SIS) and MPa (broken SIS), respectively. Note that the two lattices in states and have identical band structure (Fig. 2 (b)) since they are simply spatially inverted versions of each other and with intact TRS.
The analysis of the dispersion properties in Fig. 2 (a) reveals the existence of a degeneracy at the point from which locally linear dispersion curves emanate. This dispersion structure, known as the Dirac cone (DC) Neto et al. (2009), is protected by the lattice configuration. However, when SIS is broken (Fig. 2 (b)), the degeneracy is lifted and the initially degenerate modes can separate and give rise to a local bandgap. It is the adiabatic evolution of this specific dispersion structure under small perturbations that enables the generation of topological edge modes. We will show that the dispersion near the valleys in either the reference or the perturbed configuration can be mapped into a massless or a massive Dirac equation, respectively. Before proceeding, we should also note that the current configuration produces an incomplete bandgap (Fig. 2b). In principle, a complete bandgap could be obtained by optimizing the lattice geometric parameters Mousavi et al. (2015), however we will show that the valley-dependent design is fairly robust and does not strictly require a full bandgap.
As established above, the process of lifting the degeneracy and opening the bandgap is connected to the applied pressure perturbation within the cell. Fig. 3 (a) and (b) show the evolution of the gap bounds as well as of the width of the gap as a function of the applied pressure . An average upward shift of the frequency is observable in the gap bounds as increases. This behavior is directly related to the stress stiffening effect produced by the geometric nonlinear deformation of the cell. The analysis of the gap width indicates that the gap vanishes and reopens as crosses zero, therefore suggesting a topological phase transition in the -space.
ii.2 Berry curvature and valley Chern number
The evolution of the bandgap as a function of suggests that, a lattice obtained by connecting and domains should experience a topological phase transition associated with a vanishing gap occurring exactly at the DW. The topological nature of this transition can be characterized using a topological invariant, that is the Chern number (module ). The parameter is obtained by integrating the Berry curvature of the mode throughout the first Brillouin zone. For our system, is expected to be zero due to an odd distribution of the Berry curvature in -space, which should be expected given TRS is preserved Hasan and Kane (2010); Raghu and Haldane (2008); Ochiai and Onoda (2009). It follows that these lattices are classified as trivially gapped materials in the context of QHE systems. Nevertheless, for small SIS breaking, the Berry curvature is highly localized at the valleys, and the local integral of the Berry curvature converges quickly to a non-zero quantized value Xiao et al. (2007); Martin et al. (2008). This local integral is often referred to as the valley Chern number of the band and it is defined as , where the integral bounds extend to a local area around the valley. The right hand side of this equation is also referred to as topological charge. Previous studies in electronic systems Yao et al. (2009) have shown that this quantized value has important implications because it characterizes the bulk–edge correspondence. In fact, the difference between the valley Chern numbers of the upper (or lower) bands of two adjacent lattices indicates the number of gapless edge states expected at the DW.
Numerical integration of the Berry curvature shows that each valley of either the - or -lattices carries a topological charge of magnitude . The two lattices clearly exhibit Berry curvatures of opposite signs. It follows that, for the upper mode of the valley , the valley Chern numbers are for the -lattice and for the -lattice. The difference indicates the existence of a single gapless edge state at the DW between the - and -lattices. Fig. 4(a) shows the calculated dispersion surfaces of the flexural modes of the -lattice near the valley . Fig. 4(b) shows the Berry curvature corresponding to the upper mode. Equivalently, the lower band carries a Berry curvature of opposite sign.
Iii Domain wall edge state
There are two possible configurations of the DWs between the - and -lattices, that is above () or viceversa (see Fig. 5 (a) and (b)). These two configurations are not equivalent, and are dominated by a mirror symmetry with respect to either DW1 or DW2. As a direct consequence of the mirror symmetry, the edge states propagating along the DWs can be either symmetric or anti-symmetric (with respect to the DW interface).
This symmetry was exploited in order to reduce the size of the finite element model. In particular, we built an -state superlattice ribbon made of 37 triangular subcells stacked along the -direction (see the red part in Fig. 5 (a)) and having a length of about . Then, either symmetric or anti-symmetric boundary conditions along DW1 or DW2 were applied to simulate the presence of the -state ribbon. The model was solved to get the Bloch eigenvalues and eigenvectors. Results show that the edge state at DW1 is anti-symmetric, while at DW2 it is symmetric. Fig. 5 (c) shows the dispersion relation of such a ribbon structure.
The colorbar denotes the -position of the centroid of the superlattice weighted by the strain energy density. The blue color indicates edge states at DW1 while the red color indicates edge states at DW2. The two modes have equal and opposite group velocities but since they are supported by the two different edges of the lattice they cannot couple. In addition, if inter-valley mixing is neglected (due to large separation in -space between the and points), these edge states are immune to back-scattering. Note that this assumption is well verified in the absence of short range disorder. We highlight that this dynamic behavior is effectively equivalent to the quantum spin Hall effect in topological insulators Hasan and Kane (2010) if the spin index is replaced by the valley index.
In order to further characterize the egde modes, we performed full field numerical simulations. In particular, we simulated the steady-state response of two different DW shapes (i.e. straight and arbitrary) under a point excitation. Low-reflecting boundaries were used all around the model to suppress reflections. In both cases, results show that the edge states are well concentrated near the DWs and are guided along the wall itself (Fig. 5 (d),(e)).
Iv Semi-analytical approach
Although the band structure and the valley Chern number were accurately calculated via a numerical approach, it is still useful to introduce a semi-analytical model to further analyze the effects of SIS-breaking on the phononic waveguide. Starting from the (undeformed) reference lattice configuration which still preserves SIS, the perturbation approach Zhu and Semperlotti (2017) can be used to show that the linear dispersions at the valleys can be mapped to a massless Dirac equation. The method allows obtaining approximate solutions to the governing equations by using the degenerate modes as fundamental basis of an expansion process. In this method, the origin of the -axes is shifted to the degeneracy by letting (or ) and . Hence, the problem of determining the dispersion relations – is cast in the form of a eigenvalue problem associated to the following Hamiltonian:
where is the group velocity, and are Pauli matrices. The Hamiltonian maps to the massless Dirac equation associated with locally linear dispersion. By extracting the eigenstates from the previous simulations, we can numerically confirm that the group velocity obtained from the method matches well with the tangent slopes of the exact dispersion data obtained from numerical calculations (inset Fig. 2 (a)).
When SIS is broken by the application of a pressure perturbation , the degeneracy at (or ) is lifted, the modes become non-degenerate, and a gap opens up. It is well-known that breaking SIS introduces a -component into the Hamiltonian that can be expressed as Hasan and Kane (2010); Raghu and Haldane (2008):
It represents the Hamiltonian of a massive Dirac equation which we show applies also to our specific lattice structure:
Since the SIS-breaking perturbation is small, it can be expressed as a perturbation of the massless Hamiltonian,
where we assumed a general perturbation that potentially contains all the four Pauli matrix components () as well as the identity matrix . At the K (or ) point (), the Hamiltonian reduces to
Let be the two degenerate modes of the undeformed lattice and be the two non-degenerate modes of the deformed lattice, at . These modes are available from the numerical calculation of the band structure. By performing the inner product
it is seen that matrix is diagonal (or anti-diagonal, depending on the numbering of the modes), which implies that is also diagonal. This means that the perturbation term can only contain and components, at most. In particular, the splitting of the eigenfrequencies is connected to the component while the average upward shift of the bands (as a function of ) is connected to the component. Note that the term does not affect either the eigenvectors or the topological properties, hence it can be omitted in the Hamiltonian. The above arguments confirm that the Bloch eigenvalue problem of the SIS-broken lattice can be mapped to the massive Dirac equation (Eq. (4)).
The bulk dispersion is then given by with a gap equal to . From the numerical results in Fig. 3 (b) we find that , where is expressed in MPa.
The corresponding Berry curvature, associated with the evolution of the eigenvectors in -space, has a distribution sharply peaked at the two Dirac points,
where allows simple labeling of the two valleys and . Integrating Eq. (8) gives the valley Chern number . Under TRS, for either a lattice or , and at the same valley. Thus is always true. Hence, we can conclude that there is one gapless edge state living on the DW connecting the two lattice and .
V Edge states at external boundaries
The topological phase transition illustrated above and taking place between two lattice with broken SIS can be obtained also at the boundary of an individual lattice when proper boundary conditions are applied. The two specific boundary conditions considered in this study are traction-free and fixed. Each boundary condition can be applied on the top and bottom boundaries of the ribbon (either in - or -state), hence resulting in a total of eight possible configurations (Fig. 6 (a)-(d)). However, since the - and -lattices are inverted images of each other, a given boundary condition on the edge of an -ribbon corresponds to the inverted image of the same boundary condition on the opposite edge of a -ribbon, therefore yielding only four different edge dispersions (for example, see the pairs (a)/(c) and (b)/(d)) in Fig. 6. Note that, although they have the same edge dispersion, the valley indexes interchange.
Not all the four groups of boundary conditions yield edge modes. Only the two edges of the ribbon in Fig. 6 (a) (or the edges of ribbon in Fig. 6 (c)) have propagating edge states, as shown in Fig. 6 (e). The other two groups (Fig. 6 (b) and (d)) do not support edge states (Fig. 6 (f)). Note that in both cases (Fig. 6 (e) and (f)) there are localized modes near the free edge whose dispersion curves run along the bulk bands corresponding to the faster flexural branch. These modes have no topological significance and they are effectively surface-like waves with depth of the same order of the wavelength, which will be discussed in detail in a later section. These results can be explained by identifying the similarity between the dynamic behavior imposed by the boundary conditions and the DWs, respectively.
Recall that DW1 supports an anti-symmetric edge state. Anti-symmetry requires on the interface, where is any tangent vector to the interface, therefore having the direction of either or . To satisfy the antisymmetry condition and must be zero, while the component are non-zero. The modes, under consideration in this study, exhibit particle displacement mostly dominated by the component. When a fixed boundary condition is imposed along the edge, the particle displacement u is set to zero therefore resembling the behavior of the anti-symmetric interface. On the other hand, DW2 supports symmetric edge states. At the interface we have , without restriction on . Similarly to the discussion above, the traction free boundary condition will allow therefore behaving more closely to a symmetric interface.
As a result, each of the fixed boundaries in Fig. 6 (a) and (c) (note that this figure is equivalent to Fig. 5 in the main text; here reported for convenience) is dynamically equivalent to half of the DW1. This boundary configuration supports edge states similar to those observed along DW1, but without connecting the two continuous bulk bands (i.e. the red modes in Fig. 6 (e)). Similarly, each of the free boundaries in Fig. 6 (a) and (c) supports an edge state that does not cross the bandgap. On the contrary, the other two groups of boundary conditions (Fig. 6 (b) and (d)) break the similarity with the DWs, therefore the edge states are strongly suppressed.
Full field numerical simulations confirmed this behavior and showed that edge states can be excited on the fixed edge on the top of an -ribbon, but they cannot be excited under equivalent conditions in a -ribbon (see Fig. 7 (a) and (b)). These two conditions correspond to the fixed edges in Fig. 6 (c) and (b), respectively.
Vi lossless edge states
The results above deserve some further discussion to clarify why the edge states can be excited independently from the bulk states existing at the same frequency. Due to the lack of a full bandgap, when a point excitation is placed at the domain wall, both bulk and edge states can potentially be excited. However, two factors contribute to make the edge states predominant compared to the bulk states. First, the bulk states amplitude decays as from the source, while the edge state amplitude remains constant. Second, since the lattice is inhomogeneous, the mode shapes of either the bulk or the edge states have non-uniform distribution over the unit cell. If a point source is located at the point where the edge state has the largest displacement amplitude, most of the input energy will be injected into the edge rather than the bulk state. In Fig. 7 (b) if the edge state was suppressed, the long wavelength ripples of the bulk state would be visible. In summary, both states exists but their large separation in space allows selective excitation of the edge states. In a similar way, if a distributed source was used other than a point source, the valley injection could still be achieved by targeting a selected wavenumber that supports strong edge modes. This results could be achieved by using excitation methods such as comb transducers.
In addition, it can be shown that the edge state propagates without loss. At the frequency of the edge state, we calculate the equi-frequency contour (with the -axis is aligned with either the edge, the DW, or the boundary) of the faster bulk state which is indicated by the solid curve in Fig. 8. The curve shows an almost circular profile for the propagating bulk state while the dashed line indicates the imaginary component corresponding to an arbitrary large beyond the circle. This is equivalent to a conventional slowness diagram, but different by a constant factor . When the edge state propagates and induces bulk wave scattering, the component must be conserved, that is . Clearly since the bulk state is faster, there exists no on the solid curve satisfying this condition. The scattered bulk state then has a complex wavevector with purely imaginary component as illustrated graphically in Fig. 8. This means that the scattered bulk wave is evanescent along -direction and travels with the same speed as the edge state along -direction. As a result, the edge state is lossless and it does not necessarily require a bulk band gap as long as the bulk state is faster (in terms of phase velocity) than the edge state at the same frequency.
Vii Tuning the boundary edge states
In addition to the effect of the boundary conditions, adjusting the pressure of the outermost cells (those close to the boundary) can significantly tune the edge state dispersion. Yao et al. Yao et al. (2009) showed that on a graphene system with staggered sublattice potential, the edge states dispersion can be controlled by the on-site energy on the boundary cells. Similarly, in our system the phononic edge states dispersion can be manipulated by tuning the pressure of the outermost row of lattice. From the dispersion curve (Fig. 6 (e)), we observe that the edge state develops from the bulk band near the valley and as increases the edge state localizes more sharply on the edge. At , the energy is mostly confined in the outermost row of the lattice, therefore controlling the edge state requires tuning the outermost cell properties. For our phononic waveguide, a general trend can be established observing that a positive pressure on the outer cells tends to increase the edge state frequency at , while a negative pressure tends to decrease it. This suggests that the edge state dispersion can be tailored by controlling the local pressure in the edge cells.
Numerical results in Fig. 9 show that a combination of boundary conditions and local pressure allows tuning the edge states to either a partially gapped band, a gapless band, or a flat mode. As an example, by setting the outer cell pressure to near the fixed edge (cf. Fig. 6 (c)) the dispersion curve of the edge state bends up towards the top bulk band (Fig. 9 (a)). Similarly, by setting the outer cell pressure to near the free edge (cf. Fig. 6 (a)) the dispersion curve bends down towards the bottom bulk band (Fig. 9 (b)). Flat bands can also be obtained (Fig. 9 (c) and (d)) by applying and in the previous two cases, respectively. For the two cases supporting the gapless edge states crossing the gap, the pressure applied on the outer cells has different sign from the bulk lattice. This situation is somewhat equivalent to introducing a DW on the boundary and it results again in the occurrence of gapless edge states due to the contrast between bulk topological charges.
vii.1 Discussion of the non-topological edge state
These results deserve some additional discussion concerning the dispersion of the traction-free edge (Fig. 6 (e), (f) and Fig. 9 (b), (d)). Note that the dispersion curves highlight the existence of a localized mode near the free edge, which develops along the bulk band corresponding to the faster flexural branch. These parts have no topological significance and they are effectively surface-like waves with depth of the same order of the wavelength. To clarify this aspect further, consider a fin-like elastic waveguide with thickness and width . The fin can be thought as a continuum version of our phononic ribbon without the lattice structure. It is well-known that as the fundamental flexural mode develops into a surface wave (Fig. 10). This is similar to the and Lamb modes evolving into Rayleigh surface modes as . Therefore, this edge mode has no topological significance as confirmed also by the zero Berry curvature being associated with the bulk bands. Nevertheless, there is strong coupling between this mode and the topological edge states propagating along the same traction-free edge. To clarify the interaction between these two edge modes, we can plot again the results of Fig. 9 (b) with an emphasized nonlinear color scheme in Fig. 11. The gapless topological edge state and the non-topological edge state meet at the upper left of the cone, and there is a strong level repulsion indicating hybridization. When they enter the continuous bulk band of the other, they become strongly leaky due to the coupling with bulk modes.
Viii Selective valley-injection and one-way propagating edge states
Finally we show that, despite the TRS remains intact, the two valley-dependent edge states can be excited individually, therefore giving rise to one way propagation. Since the two valley edge modes are largely separated in momentum space, the coupling between them is vanishingly small. This means that valley-injection can be achieved by using a source having stronger coupling with one of the two valleys. As an example, referring to the Bloch modes, two point sources distanced by and having a phase difference can selectively excite a specific mode. Full field numerical results are shown in Fig. 12 (a, b) for the DW and in Fig. 12 (c, d) for the fixed boundary. Again, low reflecting boundaries were used all around the model to suppress reflections.
In conclusion, we presented the design of a tunable topological elastic phononic waveguide based on the acoustic analog of the QVHE. The lattice structure exploited SIS-breaking while preserving TRS, hence resulting in a weak topological acoustic material. SIS-breaking was induced by tunable elastic deformation of the original lattice structure while the topological phase transition was achieved by contrasting ribbons having different topological charges. This configuration could produce clear edge states at the domain wall interface. In a similar way, the proposed design was shown to be able to achieve topological edge states even in a single lattice configuration when in presence of specific boundary conditions, and the edge states could be efficiently tailored by tuning the pressure in the outermost edge cells. We also showed that, despite TRS being preserved, selective valley injection could be effectively achieved by a two-point excitation therefore giving rise to uni-directional edge modes. We also note that, in presence of smooth disorder, the inter-valley mixing is fairly weak which allows for back-scattering immune edge states.
Acknowledgements.The authors gratefully acknowledge the financial support of the Air Force Office of Scientific Research under the grant YIP FA9550-15-1-0133.
- M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
- Y. Ren, Z. Qiao, and Q. Niu, Rep. Prog. Phys. 79, 066501 (2016).
- Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong, and B. Zhang, Phys. Rev. Lett. 114, 114301 (2015).
- R. Süsstrunk and S. D. Huber, Science 349, 47 (2015).
- P. Wang, L. Lu, and K. Bertoldi, Phys. Rev. Lett. 115, 104302 (2015).
- A. B. Khanikaev, R. Fleury, S. H. Mousavi, and A. Alu, Nat. Comms. 6, 8260 (2015).
- S. H. Mousavi, A. B. Khanikaev, and Z. Wang, Nat. Commun. 6, 8682 (2015).
- L. M. Nasha, D. Klecknera, A. Reada, V. Vitellib, A. M. Turnerc, and W. T. M. Irvine, Proc. Natl. Acad. Sci. U.S.A. 112, 14495 (2015).
- J. Lu, C. Qiu, M. Ke, and Z. Liu, Phys. Rev. Lett. 116, 093901 (2016a).
- Z.-G. Chen and Y. Wu, Phys. Rev. Appl. 5, 054021 (2016).
- C. He, X. Ni, H. Ge, X.-C. Sun, Y.-B. Chen, M.-H. Lu, X.-P. Liu, and Y.-F. Chen, Nat. Phys. 12, 1124 (2016).
- J. Lu, C. Qiu, L. Ye, X. Fan, M. Ke, F. Zhang, and Z. Liu, Nat. Phys. 13, 369 (2016b).
- D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007).
- G. Semenoff, V. Semenoff, and F. Zhou, Phys. Rev. Lett. 101, 087204 (2008).
- W. Yao, S. A. Yang, and Q. Niu, Phys. Rev. Lett. 102, 096801 (2009).
- N. D. Drummond, V. ZÂ´olyomi, and V. I. Fal’ko, Phys. Rev. B 85, 075423 (2012).
- Y. Kim, K. Choi, and J. Ihm, Phys. Rev. B 89, 085429 (2014).
- F. Guinea, M. I. Katsnelson, and A. K. Geim, Nat. Phys. 6, 30 (2009).
- J. Jung, F. Zhang, Z. Qiao, and A. H. MacDonald, Phys. Rev. B 84, 075418 (2011).
- A. Vaezi, Y. Liang, D. H. Ngai, L. Yang, and E.-A. Kim, Phys. Rev. X 3, 021018 (2013).
- L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. V. Jr, C. Ojeda-Aristizabal, H. A. Bechtel, M. C. Martin, A. Zettl, J. Analytis, and F. Wang, Nature 520, 650 (2015).
- F. Zhang, A. H. MacDonald, and E. J. Mele, Proc. Natl. Acad. Sci. U.S.A. 110, 10546 (2013).
- I. Martin, Y. Blanter, and A. F. Morpurgo, Phys. Rev. Lett. 100, 036804 (2008).
- Z. Qiao, J. Jung, Q. Niu, and A. H. MacDonald, Nano Letters 11, 3453 (2011), pMID: 21766817.
- J. Li, K. Wang, K. J. McFaul, Z. Zern, Y. Ren, K. Watanabe, T. Taniguchi, Z. Qiao, and J. Zhu, Nat. Nanotechnol. 11, 1060â1065 (2016).
- S. Raghu and F. D. M. Haldane, Phys. Rev. A 78, 033834 (2008).
- T. Ochiai and M. Onoda, Phys. Rev. B 80, 155103 (2009).
- R. K. Pal and M. Ruzzene, New J. Phys. 19, 025001 (2017).
- J. Vila, R. K. Pal, and M. Ruzzene, (2017), arXiv:1705.08247v1 [cond-mat.mtrl-sci] .
- A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- H. Zhu and F. Semperlotti, (2017), arXiv:1706.01536v1 [cond-mat.mtrl-sci] .