Mechanical graphene

Topology [1–3] has become an important tool in advancing our understanding of electronic properties of solids. It plays an important role, for example, in determining the nature of surface states in a wide variety of systems including polyacetylene [4, 5], quantum Hall systems [6, 7], and topological insulators [8–13]. The success in applying topological ideas to electronic systems has recently inspired their generalization to certain classes of classical mechanical systems [14–32] by establishing a correspondence between quantum electronic Hamiltonians and the 'square root' of the mechanical dynamical matrix. This correspondence is exact for a class of Maxwell lattices for which there is a balance between the number of degrees of freedom and the number of constraints per unit cell. A prototype example is the one-dimensional mechanical model [14, 18, 33] whose excitation spectrum precisely matches the spectrum of the Su–Schrieffer–Heeger (SSH) model for polyacetylene. In the one-dimensional SSH model, electrons move on a bipartite lattice with different hopping matrix elements on alternating bonds that connect the sites on the A- and B-sublattice. The mechanical model consists of rigid rotors with pivot points fixed on the A-sublattice and connected by central-force springs residing on the the B-sublattice. There is a one-to-one correspondence between electrons on the A-sublattice in the SSH model and the rotors in mechanical model and between electrons on the B-sublattice and the springs of the mechanical model. We note that in contrast to other models of topological mechanics [23–31, 34] Maxwell lattices exhibit topologically protected zero frequency modes, along with an intrinsic particle–hole symmetry in the analog quantum Hamiltonian.

Here, we introduce and explore the properties of a model mechanical system whose bulk vibrational excitations are in correspondence with the electronic excitations of the particle–hole symmetric two-band tight-binding model of graphene. We consider a generalized graphene model describing nearest neighbor hopping on a honeycomb lattice in which the hopping matrix elements for the three bonds emanating from a given site are in general different. This model can represent strained graphene [35] for modest variations in the hopping amplitudes, and also arises in the analysis of Kitaev's honeycomb lattice model [36]. Our mechanical analog is constructed by closely following the paradigm used in the construction of the SSH analog. The honeycomb lattice of graphene, like the 1D SSH lattice, is bipartite with A- and B-sublattices, shown as red and blue disks in figure 1. Since there is only one degree of freedom per site in the graphene model, we assign a scalar variable $z({{\boldsymbol{R}}}_{A})$, which we identify as vertical height displacement, to each site on the A-sublattice, and we assign a kinetic energy ${\dot{z}}^{2}({{\boldsymbol{R}}}_{A})/2$ to that site. On each site in the B-sublattice, we need to assign the analog of a bond connecting sites on the A-sublattice, but this 'bond' is connected to three rather than the usual two sites. To each of these triangular bonds, which we call tri-bonds, we assign a kind of 'stretching' energy that depends quadratically on a certain linear combination of the heights at the three sites. It is important to keep straight that $z({{\boldsymbol{R}}}_{A})$ denotes the vertex heights of tri-bonds arranged on a rigid lattice, not the heights of the actual sites of a honeycomb lattice. Figure 1 shows our model mechanical graphene lattice, with red sites representing the locations of the height variables and gray triangles representing the tri-bonds.

Zoom In Zoom Out Reset image size

The tri-bonds are the analogs of springs in a standard ball-and-spring lattice, and each imposes a constraint that leads to a generalization of the Maxwell–Calladine theorem [37] relating the numbers of zero modes and states of self stress (SSS) to the number of degrees of freedom and number of constraints. Under periodic boundary conditions, our mechanical graphene lattice is a generalized Maxwell lattice [15, 37] in which the number of constraints equals the number of degrees of freedom, and it exhibits all of the properties of ball-and-spring Maxwell lattices: (1) Each zero mode in the bulk spectrum, which again matches the electronic spectrum of graphene, is accompanied by a SSS. (2) Each lattice is characterized by a topological polarization or by Weyl modes. (3) Periodic strips or finite lattices cut from a periodic lattice whose spectrum is fully gapped have a number of zero-frequency surface modes equal to the number of tri-bonds cut, i.e., at least one zero mode per surface wavenumber, residing on one or the other of the opposite free edges. (4) The number of zero modes at a given wavenumber on a free surface is determined by the topological polarization or by the positions of Weyl modes and by a local surface polarization. (5) Domain walls connecting lattices of different topological polarization harbor either zero modes or self stress for each wavenumber along the wall; those connecting different Weyl lattices have zero modes or SSS for some wavenumbers but not for others.

Though the bulk spectrum of our model and graphene are identical, their surfaces modes are not. The top and bottom edges of a horizontal strip of the mechanical lattice, which can be created by removing the row of horizontal of dashed tri-bonds in figure 1 from a periodic lattice, are different: the top surface exposes tri-bond vertices and the bottom a straight continuous line of tri-bond edges. These two surfaces correspond, respectively, to the bearded (with dangling bonds) and zigzag edges of graphene [38]. It is not possible to create a strip in the mechanical lattice like that in graphene in which both edges are zigzagged. Both edges of vertical strips exhibit a two-tri-bond zigzag pattern. The corresponding graphene edges correspond to armchair edges with an extra dangling bond at every second row.

This paper is divided into five sections of which this is the first. Section 2 presents details of our model and defines the equilibrium and compatibility matrices that establish the Maxwell–Calladine theorem. Section 3 treats the excitation spectrum and establishes a phase diagram with a region harboring Weyl points and regions that carry a topological polarization. Section 4 discusses zero modes at free edges or in domain walls. Section 5 presents physical models for the tri-bonds.

Figure 1 provides a visual image of our model. We take the Bravais lattice constant to be a and define primitive reciprocal lattice vectors

$$\begin{eqnarray}&&{{\bf{a}}}_{1}=a\hat{x},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\bf{a}}}_{2}=a\left(-\displaystyle \frac{1}{2}\hat{x}+\displaystyle \frac{\sqrt{3}}{2}\hat{y}\right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\bf{a}}}_{3}=a\left(-\displaystyle \frac{1}{2}\hat{x}-\displaystyle \frac{\sqrt{3}}{2}\hat{y}\right).\end{eqnarray} \tag{ 3 }$$

We also define the three vectors from a tri-bond centroid to its vertices:

$$\begin{eqnarray}&&{{\bf{d}}}_{1}=\displaystyle \frac{a}{\sqrt{3}}\hat{y},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\bf{d}}}_{2}=\displaystyle \frac{a}{\sqrt{3}}\left(-\displaystyle \frac{\sqrt{3}}{2}\hat{x}-\displaystyle \frac{1}{2}\hat{y}\right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\bf{d}}}_{3}=\displaystyle \frac{a}{\sqrt{3}}\left(\displaystyle \frac{\sqrt{3}}{2}\hat{x}-\displaystyle \frac{1}{2}\hat{y}\right).\end{eqnarray} \tag{ 6 }$$

The centroids of the tri-bonds are located on B sites of the honeycomb and their vertices lie at A sites. Each site on the A-sublattice is occupied by a unit mass that is constrained (say by a frictionless rod) to move in the vertical direction. We assume these vertical displacements $z({{\boldsymbol{R}}}_{A})$ are small enough that linear approximations can be used. Each tri-bond is pinned at a pivot point that is displaced from its centroid by

$$\begin{eqnarray}&&{{\boldsymbol{r}}}_{p}={x}_{1}{{\boldsymbol{d}}}_{1}+{x}_{2}{{\boldsymbol{d}}}_{2}+{x}_{3}{{\boldsymbol{d}}}_{3},\end{eqnarray} \tag{ 7 }$$

and the three tri-bonds meeting at a given site are connected in a way that allows each to freely rotate about an axis passing through its pivot point and the site in question.

Any location of the pivot point (within or outside the triangle spanned by the tri-bond vertices) can be specified by a unique triple $x=({x}_{1},{x}_{2},{x}_{3})$ with ${x}_{1}+{x}_{2}+{x}_{3}=1$. With this parametrization, the condition

$$\begin{eqnarray}&&e({{\boldsymbol{R}}}_{B})\equiv \displaystyle \sum _{i=1}^{3}{x}_{i}z({{\boldsymbol{R}}}_{B}+{{\boldsymbol{d}}}_{i})=0,\end{eqnarray} \tag{ 8 }$$

is satisfied for any rigid rotation of the tri-bond, where ${{\boldsymbol{R}}}_{B}$ is the B-sublattice location of the tri-bond centroid and ${{\boldsymbol{R}}}_{B}+{{\boldsymbol{d}}}_{i}$ are the A-sublattice positions of the vertices of the tri-bond. Violations of this condition necessarily cause a distortion of the tri-bond that costs some energy, which to lowest order in the displacements must be quadratic in $e({{\boldsymbol{R}}}_{B})$, giving us a total elastic energy

$$\begin{eqnarray}&&U=\displaystyle \frac{1}{2}k\displaystyle \sum _{B}{e}^{2}({{\boldsymbol{R}}}_{B}),\end{eqnarray} \tag{ 9 }$$

where the sum is over all sites ${{\boldsymbol{R}}}_{B}$ in the B-sublattice. This energy has exactly the same form as that of a lattice of harmonic springs where the sum is over bonds and e is the elongation of a bond. Following this analogy, we call $e({{\boldsymbol{R}}}_{B})$ the stretch of the tri-bond at ${{\boldsymbol{R}}}_{B}$. Of course, with the aid of equation (8) we can express U as a function of the displacements $z({{\boldsymbol{R}}}_{A})$ instead of stretches $e({{\boldsymbol{R}}}_{B})$.

In section 5, we will discuss two physical models for the tri-bond that produce the energy of equation (9) to harmonic order, one in which the energy arises from the stretching of springs attached to rigid tri-bond plates and the other in which tri-bonds are thin elastic sheets and the energy arises from bending them out of the plane. We will obtain equation (9) as an exact expression for the first model. For the second model we will argue that k is a function of the x_i's, which distorts the spectrum but does not affect the zero modes.

Given U, we can now construct expressions for the N_B tri-bond tensions $t({{\boldsymbol{R}}}_{B})$ conjugate to stretches of the tri-bonds and the N site forces $f({{\boldsymbol{R}}}_{A})$ conjugate to vertical displacements of the masses. As in the case of harmonic springs, $t({R}_{B})$ is simply the derivative of the energy with respect to $e({{\boldsymbol{R}}}_{B})$ and $f({{\boldsymbol{R}}}_{A})$ is the derivative of the energy with respect to the dynamical variable $z({{\boldsymbol{R}}}_{A})$:

$$\begin{eqnarray}&&t({{\boldsymbol{R}}}_{B})=\displaystyle \frac{\partial U}{\partial e({{\boldsymbol{R}}}_{B})}={ke}({{\boldsymbol{R}}}_{B}),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&f({{\boldsymbol{R}}}_{A})=\displaystyle \frac{\partial U}{\partial z({{\boldsymbol{R}}}_{A})}=\displaystyle \sum _{i}{x}_{i}\,t({{\boldsymbol{R}}}_{A}-{{\boldsymbol{d}}}_{i})\end{eqnarray} \tag{ 11 }$$

which may be written in terms of the displacements as

$$\begin{eqnarray}&&t({{\boldsymbol{R}}}_{B})=k\displaystyle \sum _{i}{x}_{i}\,z({{\boldsymbol{R}}}_{B}+{{\boldsymbol{d}}}_{i}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&f({{\boldsymbol{R}}}_{A})=k\displaystyle \sum _{i,j}{x}_{i}\,{x}_{j}\,z({{\boldsymbol{R}}}_{A}-{{\boldsymbol{d}}}_{i}+{{\boldsymbol{d}}}_{j}).\end{eqnarray} \tag{ 13 }$$

The tension t is one that induces 'stretching' of the tri-bond in a manner exactly analogous to the tension in a spring bond inducing stretching of the spring, and the force $f({{\boldsymbol{R}}}_{A})$ is a linear combination of the tension of bonds with endpoints sharing R_A, as is the case in systems of harmonic springs.

Introducing the N-dimensional vectors ${\boldsymbol{z}}$ and ${\boldsymbol{f}}$ of site-displacements and forces and the N_B dimensional vectors of ${\boldsymbol{e}}$ and ${\boldsymbol{t}}$ of stretches and tensions, we can write equations (8), (10), and (11) as

$$\begin{eqnarray}&&{\boldsymbol{e}}=C{\boldsymbol{z}};\qquad {\boldsymbol{t}}=k{\boldsymbol{e}};\qquad {\boldsymbol{f}}=Q{\boldsymbol{t}},\end{eqnarray} \tag{ 14 }$$

where C is the compatibility matrix with components

$$\begin{eqnarray}&&C({{\boldsymbol{R}}}_{B},{{\boldsymbol{R}}}_{A})=\displaystyle \sum _{i}{x}_{i}{\delta }_{{{\boldsymbol{R}}}_{B}+{{\boldsymbol{d}}}_{i},{{\boldsymbol{R}}}_{A}},\end{eqnarray} \tag{ 15 }$$

and Q is the equilibrium matrix with components

$$\begin{eqnarray}&&Q({{\boldsymbol{R}}}_{A},{{\boldsymbol{R}}}_{B})=\displaystyle \sum _{i}{x}_{i}{\delta }_{{{\boldsymbol{R}}}_{A}-{{\boldsymbol{d}}}_{i},{{\boldsymbol{R}}}_{B}}.\end{eqnarray} \tag{ 16 }$$

$Q={C}^{T}$ as required. In Fourier space,

$$\begin{eqnarray}&&e({\boldsymbol{q}})=C({\boldsymbol{q}})z({\boldsymbol{q}});\qquad f({\boldsymbol{q}})=Q({\boldsymbol{q}})t({\boldsymbol{q}}),\end{eqnarray} \tag{ 17 }$$

where $C({\boldsymbol{q}})$ and $Q({\boldsymbol{q}})$ in this case are 1 × 1 matrices with

$$\begin{eqnarray}&&C({\boldsymbol{q}})=\displaystyle \sum _{i}{x}_{i}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{q}}\cdot {{\boldsymbol{d}}}_{i}}={Q}^{* }({\boldsymbol{q}}).\end{eqnarray} \tag{ 18 }$$

As in the case of central-force springs, the null space of C consists of zero modes and that of Q of SSSs. The global and wave-number specific Maxwell–Calladine index theorems [15] follow immediately:

$$\begin{eqnarray}&&{N}_{0}-S=N-{N}_{B}\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}&&{n}_{0}({\boldsymbol{q}})-s({\boldsymbol{q}})=n-{n}_{B},\end{eqnarray} \tag{ 20 }$$

where N₀ and S are, respectively, the total number of zero modes and the total number of SSSs, ${n}_{0}({\boldsymbol{q}})$ and $s({\boldsymbol{q}})$ are the numbers of zero modes and SSS at wavenumber ${\boldsymbol{q}}$, and n and n_B are the number of sites (1 under PBC) and number of tri-bonds (1 under PBC) per unit cell.

The energy can be written in various ways in terms of these variables:

$$\begin{eqnarray}&&E=\displaystyle \frac{k}{2N}\displaystyle \sum _{{\boldsymbol{q}}}| e({\boldsymbol{q}}){| }^{2}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1}{2N}\displaystyle \sum _{{\boldsymbol{q}}}z(-{\boldsymbol{q}})D({\boldsymbol{q}})z({\boldsymbol{q}})\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1}{2{kN}}\displaystyle \sum _{{\boldsymbol{q}}}| t({\boldsymbol{q}}){| }^{2},\end{eqnarray} \tag{ 23 }$$

where $D({\boldsymbol{q}})=k| C({\boldsymbol{q}}){| }^{2}$ is the (one-dimensional) dynamical matrix. For the system with PBCs, the corresponding quantum Hamiltonian is block diagonal with 2 × 2 blocks of the form

$$\begin{eqnarray}H({\boldsymbol{q}})={\omega }_{0}\left(\begin{array}{cc}0 & {C}^{* }({\boldsymbol{q}})\\ C({\boldsymbol{q}}) & 0\end{array}\right),\end{eqnarray} \tag{ 24 }$$

where we define the normal mode frequency scale ${\omega }_{0}=\sqrt{k}$. The square of this Hamiltonian is diagonal with lower entry ${{kCC}}^{* }$ and upper entry ${{kC}}^{* }C=D$ (which are identical as C is one-dimensional). The eigenvalues of $H({\boldsymbol{q}})$, given by $\pm \omega ({\boldsymbol{q}})$, are simply the two square roots of D and thus specify the normal mode dispersion $\omega ({\boldsymbol{q}})$, with $C(0)=1$ implying ${\omega }_{0}=\omega ({\boldsymbol{q}}=0)$.

If we identify the bond hopping amplitudes ${t}_{i}={\omega }_{0}{x}_{i}$, then $H({\boldsymbol{q}})$ is formally identical to the Hamiltonian of the nearest neighbor tight binding model of graphene [35]. Equivalently, the constraint $\sum {x}_{i}=1$ corresponds to setting ${\omega }_{0}={\sum }_{i}{t}_{i}$.

The spectrum and topological properties of mechanical graphene depend on the location of the pivot point in the tri-bond. The features of the different models can thus be represented on a ternary phase diagram in which each point $({x}_{1},{x}_{2},{x}_{3})$ corresponds precisely to the placement of the pivot point. Unlike a typical ternary phase diagram, however, negative values of x_i (corresponding to pivot points outside the triangle) are allowed here.

In this section, we will derive the phase diagram shown in figure 2(a), in which there are three gapped regions with different topological polarizations ${{\boldsymbol{R}}}_{T}$ and four regions, ${{W}}_{\mathrm{0,1},\mathrm{2,3}}$, with Weyl points. The point $x=(1/3,1/3,1/3)$ corresponds to undistorted graphene. Other points correspond to strained graphene, though the regions outside ${{W}}_{0}$ correspond to degrees of strain that are probably too large to be physically realizable.

Zoom In Zoom Out Reset image size

Weyl points arise when $C({\boldsymbol{q}})=0$, which generically occurs for pairs of points $\pm {{\boldsymbol{q}}}^{* }$, where ${{\boldsymbol{q}}}^{* }$ is a solution of

$$\begin{eqnarray}&&C({{\boldsymbol{q}}}^{* })={x}_{1}{{\rm{e}}}^{{\rm{i}}{{\boldsymbol{q}}}^{* }\cdot {{\boldsymbol{d}}}_{1}}+{x}_{2}{{\rm{e}}}^{{\rm{i}}{{\boldsymbol{q}}}^{* }\cdot {{\boldsymbol{d}}}_{2}}+{x}_{3}{{\rm{e}}}^{{\rm{i}}{{\boldsymbol{q}}}^{* }\cdot {{\boldsymbol{d}}}_{3}}=0.\end{eqnarray} \tag{ 25 }$$

When $x=(1/3,1/3,1/3)$ (at the middle of ${{W}}_{0}$ in figure 2(a)) equation (25) is satisfied at the $\pm {K}$ points ${{\boldsymbol{q}}}^{* }=\pm {{\boldsymbol{q}}}_{{K}}$, where

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{{K}}=\displaystyle \frac{4\pi }{3a}\hat{x}\quad \Rightarrow \quad {{\rm{e}}}^{{\rm{i}}{{\boldsymbol{q}}}_{{K}}\cdot {{\boldsymbol{d}}}_{j}}={{\rm{e}}}^{2\pi {\rm{i}}(j-1)/3}.\end{eqnarray} \tag{ 26 }$$

This corresponds to the well known Weyl point at the Brillouin zone corner in unstrained graphene. In figure 3(a) we show the displacements of one of the zero frequency modes at ${{\boldsymbol{q}}}_{{K}}$. It is straightforward to see in figure 3(a) that this is indeed a 'floppy mode' in which, to linear order, each tri-bond undergoes a rigid rotation about its pivot point and hence is not stretched.

Zoom In Zoom Out Reset image size

Because the phase of $C({\boldsymbol{q}})$ advances by $2\pi$ as ${\boldsymbol{q}}$ wraps around ${{\boldsymbol{q}}}^{* }$, the Weyl points are locally protected and cannot be removed by a smooth deformation. Therefore, there must be a Weyl phase in a finite region around $(1/3,1/3,1/3)$. However, it is also clear from (25) that for $x=(1,0,0)$, $(0,1,0)$ or $(0,0,1)$ there are no solutions to $C({{\boldsymbol{q}}}^{* })=0$. Thus, there must be a gapped phase in the vicinity of the three corners of the dashed triangle in figure 2(a).

To determine the phase boundaries, we note that Weyl points can only only disappear if they meet, which must occur at a time-reversal-invariant point ${{\boldsymbol{q}}}^{* }=-{{\boldsymbol{q}}}^{* }+{\boldsymbol{G}}$, where ${\boldsymbol{G}}$ is a reciprocal lattice vector. ${{\boldsymbol{q}}}^{* }=0$ is ruled out for finite x because $C(0)={x}_{1}+{x}_{2}+{x}_{3}=1$, so this must occur at one of the three M points

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{{M}j}=\displaystyle \frac{2\pi }{\sqrt{3}a}{\hat{{\boldsymbol{d}}}}_{j},\end{eqnarray} \tag{ 27 }$$

which satisfy

$$\begin{eqnarray}{{\boldsymbol{q}}}_{{M}j}\cdot {{\boldsymbol{d}}}_{k}=\left\{\begin{array}{ll}2\pi /3 & j=k\\ -\pi /3 & j\ne k\end{array}\right..\end{eqnarray} \tag{ 28 }$$

Then equation (25) requires

$$\begin{eqnarray}&&{x}_{k}-{x}_{k-1}-{x}_{k+1}=0,\end{eqnarray} \tag{ 29 }$$

where the subscript is defined cyclically. Together with $\sum {x}_{i}=1$, this implies ${x}_{k}=1/2$ and ${x}_{k-1}+{x}_{k+1}=1/2$, which define the three red lines bounding the triangle inscribed in the dashed region of figure 2(a). The Weyl phase in the vicinity of $x=(1/3,1/3,1/3)$ thus corresponds to the region ${{W}}_{0}$. On the boundary of ${{W}}_{0}$ the Weyl points meet and annihilate.

The points outside the dashed triangle of figure 2(a) have the pivot outside the tri-bond and have one or two of the x_i negative. This corresponds in the graphene model to having bond(s) with a negative hopping amplitude. Systems with negative hopping amplitudes are closely related to systems with positive hopping amplitudes. The sign of the hopping amplitude on one of the three bonds (say along ${{\bf{d}}}_{1}$) can be changed by a nonuniform gauge transformation that changes the signs of all sites on every other horizontal (zig-zag) line of bonds on the honeycomb lattice. This transformation takes bond hopping amplitudes $({t}_{1},{t}_{2},{t}_{3})\to (-{t}_{1},{t}_{2},{t}_{3})$, which in our mechanical model takes $x=({x}_{1},{x}_{2},{x}_{3})\to x^{\prime} =(-{x}_{1},{x}_{2},{x}_{3})/({x}_{2}+{x}_{3}-{x}_{1})$. This gauge transformation does not change the normal mode frequencies except for an overall constant factor due to the relation between x_i and t_i. However, since the gauge transformation is at a nonzero wavevector ${{\boldsymbol{q}}}_{{M}1}$, it leads to a shift in the wavevector of the normal modes and hence a transformation of the dispersion relation:

$$\begin{eqnarray}&&{\omega }_{x^{\prime} }({\boldsymbol{q}})={({x}_{2}+{x}_{3}-{x}_{1})}^{-1}{\omega }_{x}({\boldsymbol{q}}+{{\boldsymbol{q}}}_{{M}1}).\end{eqnarray} \tag{ 30 }$$

This transformation maps $x=(1/3,1/3,1/3)$ to $x^{\prime} =(-1,1,1)$, so the normal mode spectra at these two pivot point locations are related by equation (30). In particular, for $x^{\prime}$ there are Weyl points at ${{\boldsymbol{q}}}_{{K}}+{{\boldsymbol{q}}}_{{M}1}$. Indeed, equation (25) is satisfied for ${{\boldsymbol{q}}}^{* }=\pm {{\boldsymbol{q}}}_{{L}}$, where

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{{L}}=\displaystyle \frac{2\pi }{3a}\hat{x}\quad \Rightarrow \quad {{\rm{e}}}^{{\rm{i}}{{\boldsymbol{q}}}_{{L}}\cdot {{\boldsymbol{d}}}_{(\mathrm{1,2,3})}}=({{\rm{e}}}^{0},{{\rm{e}}}^{\pi {\rm{i}}/3},{{\rm{e}}}^{-\pi {\rm{i}}/3}).\end{eqnarray} \tag{ 31 }$$

It can be checked that ${{\boldsymbol{q}}}_{{L}}=-{{\boldsymbol{q}}}_{{K}}+{{\boldsymbol{q}}}_{{M}1}$ up to a reciprocal lattice vector. More generally, the Weyl phase ${{W}}_{0}$ maps to the region ${{W}}_{1}$ in figure 2(a). Similar transformations identify the Weyl phases ${{W}}_{2}$ and ${{W}}_{3}$, whose boundaries are again given by the red lines defined by equation (29).

When x is in W₀ in figure 2(a), the Weyl points are at $\pm {{\boldsymbol{q}}}^{* }$, which reside in the dark gray regions of figure 2(b). When x is in ${{W}}_{1}$, the Weyl points are at $\pm {{\boldsymbol{q}}}^{* }$ residing in the light tan regions of figure 2(b), and there are symmetry related regions corresponding to ${{W}}_{2}$ and ${{W}}_{3}$. When a path through ${{W}}_{0}$ is traversed, beginning on one edge and ending another, a pair of Weyl points created at ${{M}}_{j}$ pass through the diametrically opposite dark gray regions of figure 2(b) and annihilate at ${{M}}_{j\pm 1}$. Similarly, for a path that begins on right boundary of ${{W}}_{1}$, passes through ${{W}}_{1}$, and terminates on the left boundary, the Weyl points are born at $\pm {{M}}_{2}$, pass through the opposite tan regions of figure 2(b), and annihilate at $\pm {{M}}_{3}$. For $x=(-\alpha ,(\alpha +1)/2,(\alpha +1)/2)$ with $\alpha \gg 1$, Weyl points occur at ${{\boldsymbol{q}}}^{* }\approx \pm (2\sqrt{2/\alpha }\,/a)\hat{x}$, converging to the Γ point as $\alpha \to \infty$.

Outside the ${{W}}_{n}$ regions there is a gap everywhere in the Brillouin zone. There are three disconnected phases that are topologically distinct. The topological polarization can be most easily evaluated at the three simple points (indexed by $j=1,2,3$) ${x}_{i}={\delta }_{{ij}}$. Then,

$$\begin{eqnarray}&&C({\boldsymbol{q}})={{\rm{e}}}^{{\rm{i}}{\boldsymbol{q}}\cdot {{\boldsymbol{d}}}_{j}}.\end{eqnarray} \tag{ 32 }$$

The topological polarization ${{\boldsymbol{R}}}_{T}$ is determined by the winding numbers n_i of this phase over the independent cycles of the Brillouin zone:

$$\begin{eqnarray}&&{n}_{i}=-\displaystyle \frac{1}{2\pi {\rm{i}}}{\oint }_{{C}_{i}}{C}^{-1}{\rm{d}}{C}=-\displaystyle \frac{1}{2\pi }{{\boldsymbol{b}}}_{i}\cdot {{\boldsymbol{d}}}_{j}\end{eqnarray} \tag{ 33 }$$

where C_i is the cycle along reciprocal lattice generator ${{\boldsymbol{b}}}_{i}$, which satisfies ${{\boldsymbol{b}}}_{i}\cdot {{\boldsymbol{a}}}_{j}=2\pi {\delta }_{{ij}}$ for Bravais lattice generators ${{\boldsymbol{a}}}_{j}$. (Note the minus sign in this equation, which appears because it is defined as an integral over $C({\boldsymbol{q}})$ rather than $Q({\boldsymbol{q}})$ as in [14].) Writing

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{T}={n}_{1}{{\boldsymbol{a}}}_{1}+{n}_{2}{{\boldsymbol{a}}}_{2}\end{eqnarray} \tag{ 34 }$$

then gives

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{T}=-{{\boldsymbol{d}}}_{j},\end{eqnarray} \tag{ 35 }$$

which follows from the 'completeness' relation ${\sum }_{i=1}^{2}{{\boldsymbol{a}}}_{i}{{\boldsymbol{b}}}_{i}=2\pi I$, where I is the unit matrix. The value of ${{\boldsymbol{R}}}_{T}$ depends on the real-space positions assigned to the A and B lattice sites, i.e., on our gauge choice. In our current symmetric choice, in which the B sites do not sit at a center of inversion, ${{\boldsymbol{R}}}_{T}$ is not a lattice vector, but the differences between of ${{\boldsymbol{R}}}_{T}$'s in different phases are. Arrows indicating ${{\bf{R}}}_{T}$ are indicated in figure 2(a). (In the gauge where the origin lies at an A site, we subtract a particular ${{\boldsymbol{d}}}_{k}$ from each ${{\boldsymbol{d}}}_{j}$ in equation (32), and the n_i are all 0 or ±1.)

Finally, we note that the line in the phase diagram corresponding to ${x}_{1}=0$ (where the pivot point is along the bottom edge of the tri-bond) corresponds to a one-dimensional limit, in which the system consists of decoupled horizontal lines that are similar to the SSH model. In this case there is a direct transition between topologically distinct gapped phases, which occurs at the blue point $x=(0,1/2,1/2)$. Here

$$\begin{eqnarray}&&C({\boldsymbol{q}})=\cos ({q}_{x}a/2){{\rm{e}}}^{-{{\rm{i}}{q}}_{y}a/(2\sqrt{3})},\end{eqnarray} \tag{ 36 }$$

so that there is a line of zero modes at ${q}_{x}=\pm \pi /a$ along the vertical line that connects ${{M}}_{2}$ and ${{M}}_{3}$, indicated by a blue line in figure 2(b). This situation is analogous to that in critical kagome lattices when there are parallel straight lines of bonds [39]. The zero modes are easy to visualize: as we have discussed, rotation of a tri-bond about any axis passing through its pivot point produces no stretch and costs no energy. Consider a horizontal line of edges containing the ${{\boldsymbol{d}}}_{2}$ and ${{\boldsymbol{d}}}_{3}$ vertices of the row of tri-bonds above it and the ${{\boldsymbol{d}}}_{1}$ vertex of the row of tri-bonds below it as shown in figure 3(b). Rotating neighboring tri-bonds in the upper row by $\delta \theta$ in opposite directions about the axis passing through the pivot point and ${{\boldsymbol{d}}}_{1}$ (the top vertex) while rotating neighboring tri-bonds in the lower row in opposite directions by $\delta \theta /\sqrt{3}$ about the axis along their bottom edges produces a zero mode. This operation only affects the given rows, and there is a zero mode for each line of bonds.

Associated with each zero mode, there must be a state of self stress. Alternating equal-amplitude stresses on tri-bonds along any row produces the desired zero-force state. In the case with ${x}_{1}=0$ and ${x}_{2}={x}_{3}=1/2$, the stress tends to bend the tri-bonds symmetrically either upward or downward about a line passing through the pivot point and vertex 1. Alternation of the sign of the stresses causes vertices 2 and 3 to experience equal and opposite forces from the two tri-bonds each shares along the x axis. When ${x}_{2}\ne {x}_{3}$, neighboring rows are still decoupled, and each is equivalent to the SSH model, whose critical point occurs when ${x}_{2}={x}_{3}$. Thus the lines defined by x_i = 0 for some i, which include the perimeter of the region where all x_i are positive, correspond to the 1D limit.

Strips with periodic boundary conditions in one direction or samples with free sides on all boundaries can be produced by removing lines of tri-bonds from the system under full periodic boundary conditions. Each line of cut tri-bonds creates two free edges. Since the number of sites and tri-bonds are equal under periodic boundary conditions, the index theorem reduces to

$$\begin{eqnarray}&&{N}_{0}-S={\rm{\Delta }}{N}_{B},\qquad \end{eqnarray} \tag{ 37 }$$

where ${\rm{\Delta }}{N}_{B}$ is the total number of tri-bonds cut to produce the free edges. A similar equation applies to each wavenumber q along the cut producing a strip:

$$\begin{eqnarray}&&{n}_{0}(q)-s(q)={\rm{\Delta }}{n}_{B},\end{eqnarray} \tag{ 38 }$$

where ${\rm{\Delta }}{n}_{B}$ is the number of bonds cut per unit cell of one of the exposed edges. How zero modes are distributed on the free edges depends on the topological polarization ${{\boldsymbol{R}}}_{T}$ and a local surface polarization ${{\boldsymbol{R}}}_{L}$ [14, 15] according to the same formula derived for central-force Maxwell lattices. The number of zero modes per edge unit cell (or equivalently per edge wavenumber q) for a given edge corresponding to a lattice 'plane' indexed by the reciprocal lattice vector ${\boldsymbol{G}}$ pointing along the edge's outer normal is

$$\begin{eqnarray}&&{n}_{0}=({{\boldsymbol{R}}}_{T}+{{\boldsymbol{R}}}_{L})\cdot {\boldsymbol{G}}/(2\pi ).\end{eqnarray} \tag{ 39 }$$

The local polarization ${{\boldsymbol{R}}}_{L}$ is simply the electric polarization at the given edge that arises from assigning a charge +1 to sites on the A-sublattice and a charge −1 to sites on the B sublattice. Of course only components of ${{\boldsymbol{R}}}_{L}$ parallel to ${\boldsymbol{G}}$ contribute to n₀ so we are free to add arbitrary components to ${{\boldsymbol{R}}}_{L}$ parallel to the edge.

It is instructive to look at a couple of examples. Consider the strip with edges parallel to the x-axis as shown in figure 3(b). To produce this strip, one tri-bond per surface unit cell had to be removed, so there is a total of one zero mode per wavenumber on the two exposed edges. The local polarization ${{\boldsymbol{R}}}_{L}^{L}$ on the lower edge with outer surface normal, ${\boldsymbol{G}}=-4\pi /\sqrt{3}a\hat{y}$, is equally well represented by $-{{\boldsymbol{d}}}_{1}/2$ , ${{\boldsymbol{d}}}_{2}$, or ${{\boldsymbol{d}}}_{3}$, giving a local contribution to the number of edge zero modes of 1/3. On the upper surface ${{\boldsymbol{R}}}_{L}^{U}=+{{\boldsymbol{d}}}_{1}$, for a contribution of 2/3 to edge-mode count. The topological count for the bottom surface is respectively 2/3, $-1/3$, and $-1/3$ for ${{\boldsymbol{R}}}_{T}$ equal to $-{{\boldsymbol{d}}}_{1}$, $-{{\boldsymbol{d}}}_{2}$, and $-{{\boldsymbol{d}}}_{3}$ for a total of one zero mode on the bottom surface and none on the top surface for ${{\boldsymbol{R}}}_{T}=-{{\boldsymbol{d}}}_{1}$ and no zero mode on the bottom and one the top surface for ${{\boldsymbol{R}}}_{T}=-{{\boldsymbol{d}}}_{2},-{{\boldsymbol{d}}}_{3}$. A similar analysis for a strip parallel to the y-axis yields for the number of zero modes on the left and right surfaces $({n}_{0}^{R},{n}_{0}^{L})=(1,1)$, $(0,2)$, and $(2,0)$ for ${{\boldsymbol{R}}}_{T}$ equal to $-{{\boldsymbol{d}}}_{1}$, $-{{\boldsymbol{d}}}_{2}$, and $-{{\boldsymbol{d}}}_{3}$, respectively. When there are Weyl points, zero modes shift from one side of a strip to the other at edge wavenumbers equal to the projections of the wavenumber of the Weyl points onto the edge.

Insight into equation (39) can be gained by considering what happens to C when the sites and tri-bonds are indexed at different positions (without changing the lattice itself). Let

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{A}^{\prime }={{\boldsymbol{R}}}_{A}+{\rm{\Delta }}{{\boldsymbol{R}}}_{A}\quad \mathrm{and}\quad {{\boldsymbol{R}}}_{B}^{\prime }={{\boldsymbol{R}}}_{B}+{\rm{\Delta }}{{\boldsymbol{R}}}_{B},\end{eqnarray} \tag{ 40 }$$

and define the 'gauge-transformed' compatibility matrix,

$$\begin{eqnarray}&&C^{\prime} ({{\boldsymbol{R}}}_{A}^{\prime },{{\boldsymbol{R}}}_{B}^{\prime })=C({{\boldsymbol{R}}}_{A}^{\prime }-{\rm{\Delta }}{{\boldsymbol{R}}}_{A},{{\boldsymbol{R}}}_{B}^{\prime }-{\rm{\Delta }}{{\boldsymbol{R}}}_{B}).\end{eqnarray} \tag{ 41 }$$

Then

$$\begin{eqnarray}&&C^{\prime} ({\boldsymbol{q}})={{\rm{e}}}^{-{\rm{i}}{\boldsymbol{q}}\cdot ({\rm{\Delta }}{{\boldsymbol{R}}}_{A}-{\rm{\Delta }}{{\boldsymbol{R}}}_{B})}C({\boldsymbol{q}})| C({\boldsymbol{q}})| \,{{\rm{e}}}^{-{\rm{i}}{\boldsymbol{q}}\cdot ({\rm{\Delta }}{{\boldsymbol{R}}}_{A}-{\rm{\Delta }}{{\boldsymbol{R}}}_{B}+{{\boldsymbol{R}}}_{T})}.\end{eqnarray} \tag{ 42 }$$

Thus if we chose ${\rm{\Delta }}{{\boldsymbol{R}}}_{A}-{\rm{\Delta }}{{\boldsymbol{R}}}_{B}={{\boldsymbol{R}}}_{L}$, we find that the total polarization of $C^{\prime} ({\boldsymbol{q}})$ is ${\boldsymbol{R}}^{\prime} ={{\boldsymbol{R}}}_{L}+{{\boldsymbol{R}}}_{T}$. Let ${\boldsymbol{q}}=({q}_{\perp },{q}_{\parallel })$, where q_⊥ and q_∥ are, respectively, the components of ${\boldsymbol{q}}$ perpendicular (positive toward the sample interior) and parallel to the lattice plane in question, and define $\eta =\exp ({{\rm{i}}{q}}_{\perp }{a}_{\perp })$, where a_⊥ is the depth of the surface unit cell, and $0\lt {q}_{\perp }\leqslant | {\boldsymbol{G}}|$, where ${\boldsymbol{G}}$ is the inner normal reciprocal lattice vector associated with the lattice plane. $C^{\prime} (\eta ,{q}_{\parallel })$ contains only positive powers of η, and thus no poles in η, and as a result, the integral

$$\begin{eqnarray}&&{n}_{0}({\boldsymbol{G}})=\displaystyle \frac{1}{2\pi {\rm{i}}}\oint {\rm{d}}\eta \displaystyle \frac{{\rm{d}}}{{\rm{d}}\eta }\mathrm{ln}C^{\prime} (\eta ,{q}_{\parallel })=\displaystyle \frac{{\boldsymbol{G}}\cdot ({{\boldsymbol{R}}}_{L}+{{\boldsymbol{R}}}_{T})}{2\pi }\end{eqnarray} \tag{ 43 }$$

counts the number of zero modes at the surface determined by ${\boldsymbol{G}}$. As particular examples, consider the bottom and left edges in figure 1. In the first case, ${{\boldsymbol{R}}}_{L}=-{{\boldsymbol{d}}}_{1}/2$ and

$$\begin{eqnarray}&&C^{\prime} (\eta ,{q}_{x})={x}_{1}\eta +{x}_{2}{{\rm{e}}}^{-{{\rm{i}}{q}}_{x}a/2}+{x}_{3}{{\rm{e}}}^{{{\rm{i}}{q}}_{x}a/2},\end{eqnarray} \tag{ 44 }$$

where $\eta =\exp ({{\rm{i}}{q}}_{y}\sqrt{3}a/2)$, has at most one zero and no poles, in agreement with our result that the top and bottom surfaces can have either one zero mode or none. In the second case, ${{\boldsymbol{R}}}_{L}={{\boldsymbol{d}}}_{2}$, and

$$\begin{eqnarray}&&C^{\prime} (\eta ,{q}_{x})={x}_{3}{\eta }^{2}+{x}_{1}\eta +{x}_{2}{{\rm{e}}}^{{{\rm{i}}{q}}_{y}\sqrt{3}a/2},\end{eqnarray} \tag{ 45 }$$

where $\eta =\exp ({{\rm{i}}{q}}_{x}a/2)$. Again, there are no poles, but the highest power of $\eta$ is 2, and according to equation (43), there can be 0, 1, or 2 zero modes in agreement with our previous results.

When there are Weyl modes, the number of zero modes on a free edge of a strip will change when ${q}_{| | }$ passes through the projection of a Weyl point onto that edge [32, 40]. The total number of zero modes does not change at this transition, so there must be a change of the opposite sign in the number of zero modes on the opposite edge. In other words, zero modes move from one side of the sample to the opposite at a projection of a Weyl point. A similar phenomenon occurs at domain walls in systems under periodic boundary conditions, in which the number of zero modes equals the number of SSSs and is given by

$$\begin{eqnarray}&&{\nu }_{T}={\boldsymbol{G}}\cdot ({{\boldsymbol{R}}}_{T}^{1}-{{\boldsymbol{R}}}_{T}^{2})/2\pi ,\end{eqnarray} \tag{ 46 }$$

where ${\nu }_{T}$ is equal to the number of zero modes per wavenumber if it is positive and minus the number of SSS if it is negative. Thus a change in the number of zero modes on a zero-mode domain wall, which occurs at q_∥ equal to the projected Weyl wavenumber, must be accompanied by an equal change in the number of Weyl states of self-stress on a self-stress domain wall.

Physical realization of the mechanical graphene model poses some technical challenges, but some straightforward approaches are possible. The simplest version of the tri-bond uses a spring that directly measures the stretch. The tri-bond consists of a rigid plate suspended on a vertical spring that is attached to a rigid ceiling, as shown in figure 4. We let ${\boldsymbol{p}}$ denote the point in three dimensions where the spring is attached to the plate. In the equilibrium configuration, ${\boldsymbol{p}}$ is in the z = 0 plane at a location specified by the two-dimensional pivot point vector r_p. The plates are connected at each vertex to a ball of mass m through a mechanism that allows free rotations about the axes passing through the center of the ball and ${\boldsymbol{p}}$.

Zoom In Zoom Out Reset image size

A possible design of the couplings of a single ball to three plates uses nested universal ball joints, as shown at left in figure 5. To allow for more general placements of the pivot point, the triangular plate may be replaced by a rigid frame with a lever arm that extends arbitrarily far from its centroid, with the spring attached to its end, as shown at right in figure 5. The pivot point can thus be placed at any location in the plane, allowing the realization of the complete phase diagram of figure 2.

Zoom In Zoom Out Reset image size

The derivation of the tri-bond energy of equation (9) for small displacements z_i is straightforward. We take the zero of the coordinate system to be at the centroid of the tri-bond, and the point at the end of the spring, about which the tri-bond can freely rotate, to be at position ${\boldsymbol{p}}={{\boldsymbol{r}}}_{p}+{z}_{p}\hat{{\boldsymbol{z}}}$, where z_p is the amount by which the spring is stretched from its equilibrium position z = 0. Let the $\delta {\rm{\Omega }}$ be a small rotation vector in a horizontal plane. Then under rotations and vertical stretch of the spring we have the vertex heights

$$\begin{eqnarray}{z}_{i} & = & {z}_{p}+\delta {\rm{\Omega }}\cdot [({{\boldsymbol{d}}}_{i}-{\boldsymbol{p}})\times \hat{{\boldsymbol{z}}}]\\ & = & {z}_{p}+\delta {\rm{\Omega }}\cdot [({{\boldsymbol{d}}}_{i}-{{\boldsymbol{r}}}_{p})\times \hat{{\boldsymbol{z}}}].\end{eqnarray} \tag{ 47 }$$

This equation immediately shows us that the rotational contribution to ${\sum }_{i}{x}_{i}{z}_{i}$ is zero because:

$$\begin{eqnarray}\displaystyle \sum _{i}{x}_{i}({{\boldsymbol{d}}}_{i}-{{\boldsymbol{r}}}_{p}) & = & \displaystyle \sum _{i}{x}_{i}\left[{{\boldsymbol{d}}}_{i}-\displaystyle \sum _{j}({x}_{j}{{\boldsymbol{d}}}_{j})\right]\\ & = & \displaystyle \sum _{i}{x}_{i}[1-{x}_{1}-{x}_{2}-{x}_{3}]{{\boldsymbol{d}}}_{i}\\ & = & 0.\end{eqnarray} \tag{ 48 }$$

Thus the tri-bond stretch e is

$$\begin{eqnarray}&&e\equiv \displaystyle \sum _{i}{x}_{i}{z}_{i}=\displaystyle \sum _{i}{x}_{i}{z}_{p}={z}_{p}.\end{eqnarray} \tag{ 49 }$$

The spring energy is ${{kz}}_{p}^{2}/2$, which immediately implies equation (9).

To mimic the spectrum of graphene, we need the kinetic energy matrix expressed in terms of the corner height variables $z({{\boldsymbol{R}}}_{A})$ to be a multiple of the identity: ${E}_{K}=(1/2)m{\sum }_{A}\dot{z}{({{\boldsymbol{R}}}_{A})}^{2}$. Thus we take the mass of a plate to be negligible compared to the mass of a corner ball. Note also that because the tri-bond plates are taken to be perfectly rigid, horizontal motions of the balls are not allowed to linear order, so there is no need for the frictionless vertical rods mentioned at the beginning of section 2.

In the limit of small deviations from equilibrium, this model directly mimics the generic discussion of the tri-bond above. Because the tri-bond stretch is directly encoded in a single spring and the height variables are literally encoded as heights of the ball joints, both ${\boldsymbol{e}}$ and ${\boldsymbol{z}}$ are immediately visible. For example, the self-stress state of the one-dimensional row of figure 3 is simply an alternating pattern of tension and compression in the springs along that row.

Other physical models could exhibit the mechanical graphene phase diagram, if not the exact same spectrum. All that is required is a potential energy of the form

$$\begin{eqnarray}&&U=\displaystyle \sum _{B}\displaystyle \frac{1}{2}k({x}_{1},{x}_{2},{x}_{3}){({x}_{1}{z}_{1}+{x}_{2}{z}_{2}+{x}_{3}{z}_{3})}^{2},\end{eqnarray} \tag{ 50 }$$

where the effective stiffness k may depend on the location of the pivot. Because the features of greatest interest are zero modes, the form of the kinetic energy does not affect their existence or locations in the Brillouin zone, though it will affect the details of the spectrum at finite frequency if the kinetic energy is not simply a multiple of ${\dot{z}}_{1}^{2}+{\dot{z}}_{2}^{2}+{\dot{z}}_{3}^{2}$.

A natural example realizes the tri-bond energy as the bending energy of a triangular elastic plate that is pinned at the pivot point and attached at each vertex to two other plates. While the precise form of the lowest potential energy of the sheet for arbitrary (small) choices of the corner heights is difficult to calculate due to the boundary conditions at the corners and on the free edges, it must vanish for all rigid rotations of the plate, for which ${x}_{1}{z}_{1}+{x}_{2}{z}_{2}+{x}_{3}{z}_{3}=0$, and thus cannot depend on any other linear combination of the corner heights.

Another possibility consists of rigid plates coupled pairwise by springs at their corners. In this model, there are two degrees of freedom per plate (the two rocking angles) and thus two vibrational bands. One band is fully gapped, consisting of modes in which the net displacement of the three plate corners at any given A site is zero. The other band exhibits the mechanical graphene phases, with each plate effectively acting as a tri-bond that couples the average displacements at each A site.

This paper has introduced a mechanical model that is a precise analog of the tight-binding model of graphene, and defines an appropriate two-dimensional generalization of the SSH analog introduced in [14]. This model system exhibits a rich phase diagram of Weyl phases, along with gapped phases with distinct topological polarizations. Our proposed structures are amenable to physical implementation, and it will be interesting to construct them and to probe their mechanical mode structures.

In addition, our construction introduces the tri-bond, which opens a new avenue for studies of Maxwell lattices. It will be interesting to use this approach, and generalizations of it, to construct new classes of two and three-dimensional mechanical systems.

The authors are grateful for the hospitality of the Lorentz Institute of for Theoretical Physics at the University of Leiden, where work on this project began. This work was supported in part by NSF grant DMR-1120901 (TCL), by two Simons Investigator grants (TCL and CLK).