Algebraic Theory of Crystal Vibrations: Localization Properties of Wave Functions in Two-Dimensional Lattices

Abstract

The localization properties of the wave functions of vibrations in two-dimensional (2D) crystals are studied numerically for square and hexagonal lattices within the framework of an algebraic model. The wave functions of 2D lattices have remarkable localization properties, especially at the van Hove singularities (vHs). Finite-size sheets with a hexagonal lattice (graphene-like materials), in addition, exhibit at zero energy a localization of the wave functions at zigzag edges, so-called edge states. The striped structure of the wave functions at a vHs is particularly noteworthy. We have investigated its stability and that of the edge states with respect to perturbations in the lattice structure, and the effect of the boundary shape on the localization properties. We find that the stripes disappear instantaneously at the vHs in a square lattice when turning on the perturbation, whereas they broaden but persist at the vHss in a hexagonal lattice. For one of them, they eventually merge into edge states with increasing coupling, which, in contrast to the zero-energy edge states, are localized at armchair edges. The results are corroborated based on participation ratios, obtained under various conditions.

1. Introduction

The isolation of graphene [1], a monoatomic layer of carbon atoms arranged in a honeycomb (hexagonal) lattice, has stimulated recently considerable interest in two-dimensional (2D) structures. Most of the interest has come from the extraordinary electronic band structure [2,3] of graphene. It stems solely from the symmetry properties of the honeycomb lattice, which is formed by two interpenetrating triangular lattices [4]. Another important feature of graphene is its soft structure. Graphene shares many properties of soft membranes [5] due to the fact that it has out-of-plane vibrational modes [6,7,8] that are not found in 3D crystals. These modes, also called flexural modes [2], play a crucial role in many applications. In addition, if one of the dimensions is shrunk to zero, one can have graphene ribbons, or carbon nanotubes [9] obtained by rolling graphene along a given direction and reconnecting the carbon bonds. These structures are one-dimensional (1D). Finally, the carbon atoms can be arranged spherically to form fullerene molecules [10]. These are zero-dimensional (0D) structures with discrete energy states. They can be obtained from graphene by wrapping the sheet into a sphere [11].

Recently [12], we have introduced an algebraic approach to lattice systems particularly well suited to study 2D lattices, which is a generalization of that to molecules (0D) [13] applied in [14] and to linear lattices [15]. It was used in [12] to calculate energy dispersion relations (EDR) and the density of states (DOS) of flexural one-dimensional (1d) out-of-plane vibrations, both analytically and numerically. The numerical simulations were done on linear lattices with n sites and on square and hexagonal lattices with $n\times m$ sites. In this article, we study unusual localization properties of the wave functions of the hexagonal lattice (graphene), which are not present in other types of 2D lattices. The wave functions are obtained by means of the algebraic model.

One peculiarity of graphene sheets is the existence of zero-energy electronic states that are nonvanishing only at the sheet boundary and accordingly are called edge states. Such states were predicted to be localized at clean zigzag edges [16,17,18,19], or more generally, at parts of the boundary where the particle-hole sublattice symmetry is not broken [20]. Previous studies focused on these edge states [21]. They occur in the energy region where the dispersion relation of graphene is linear [2,3]. The band structure, however, becomes more complex with increasing energy and eventually the dispersion relation becomes quadratic [22]. The transition takes place at the van Hove singularities in the DOS, which generally occur in two-dimensional crystals with a periodic structure [23]. There, the wave functions exhibit a further peculiar feature both in square and in hexagonal lattices, namely striped flexural modes, which, when frozen, produce the so-called ripples. Various edge geometries and edge states have been studied experimentally in graphene [24,25,26,27]. Furthermore, van Hove singularities have been observed in measurements [28,29,30,31,32]. In the vicinity of a van Hove singularity, the DOS diverges logarithmically. Accordingly, arbitrarily weak distortions in the graphene structure might produce large effects in the electronic properties. Similarly, edge states might be affected by distortions of the edge structure, which, actually, were detected in the experiments. We present results on the robustness of the localization exhibited by the wave functions at zero energy [21] and at the van Hove singularities with respect to perturbations in the lattice structure and to the choice of the boundary shape. In addition, the connection with soft membranes will become apparent when analyzing the wave functions.

The algebraic method can be used both for bosons and fermions (electrons). The structure of the equations is identical and depends only on the symmetry of the unit cell of the lattice, $D_{6h}$ for the honeycomb lattice. The problem of out-of-plane vibrations of the lattice is therefore formally identical to that of the electronic band structure of the lattice, except for the replacement of boson operators with fermion operators. The only difference is that, while for bosons, we can put any number at each site, for fermions, we can put only one at each site (or two if spin is taken into account). This analogy will be also briefly mentioned.

First, we briefly outline the main features of the algebraic model in Section 2. The associated algebraic Hamiltonian, actually, has a structure similar to that of the corresponding tight-binding Hamiltonian [2]. In Section 3 and Section 4, we review the salient features of the wave functions of finite-size square and hexagonal lattices, respectively. Finally, we investigate in Section 5 the unusual localization properties of the wave functions in terms of participation ratios. Here, we consider hexagonal lattices of rectangular and of Africa shape. Classical billiards of corresponding shapes exhibit integrable and chaotic dynamics, respectively. In References [33,34], the spectral properties of artificial graphene billiards of these shapes were studied experimentally and shown to depend on the properties of the classical dynamics. Similarly, we will demonstrate that the choice of a hexagonal lattice with the shape of the billiard with chaotic dynamics leads to stronger localization of the edge states and the states in the vicinity of the van Hove singularities than that with integrable dynamics.

2. Algebraic Model

The algebraic model of crystal vibrations, denoted in short here AMC, has been introduced in detail in Reference [12]. Therefore, we will review only briefly its main ideas. The AMC consists in inserting an algebra g at each site i and coupling the algebra at site i with that at site j. The model is a generalization of the algebraic approach to molecules (0D) [13]. For one-dimensional (1d) vibrations, the algebra is $g\equiv u\left(2\right)$, for two-dimensional (2d) vibrations, it is $g\equiv u\left(3\right)$, and for three-dimensional (3d) vibrations, it is $g\equiv u\left(4\right)$. One-dimensional vibrations have been studied both in molecules (0D) [14], and in polymers (1D) [15,35,36]. Applications to 2D lattices were introduced in [12], where both harmonic and anharmonic vibrations were studied.

The $u\left(2\right)$ algebra at each site can be realized with two boson operators $b,s$ [37]. The elements of $u\left(2\right)$ are

$$\begin{array}{c}u\left(2\right)\doteqdot {\widehat{J}}_{+}=b^{†}s,{\widehat{J}}_{-}=s^{†}b,{\widehat{J}}_z=\frac{1}{2}\left(b^{†}b-s^{†}s\right),\\ \widehat{N}=\left(b^{†}b+s^{†}s\right).\end{array}$$

(1)

For harmonic vibrations, we need to consider only the Heisenberg–Weyl algebra

$$hw\left(2\right)\doteqdot b^{†},b,b^{†}b,1,$$

(2)

obtained from Equation (1) by the process of contraction [38].

The algebraic Hamiltonian of the lattice takes then the general form

$$\begin{array}{c}H=\epsilon {\sum }_{i=1}^n{b}_i^{†}{b}_i+\frac{U}{2}{\sum }_{i=1}^n{b}_i^{†}{b}_i\left(b_i^{†}{b}_i-1\right)+V{\sum }_{i\ne j}^n{b}_i^{†}{b}_i{b}_j^{†}{b}_j\hfill \\ +{\sum }_{i\ne j}^n{\lambda }_{ij}\left(b_i^{†}{b}_j+b_j^{†}{b}_i\right),\hfill \end{array}$$

(3)

where the sum goes over the n lattice sites. The last term is also called “hopping” or “exchange” Majorana operator ${\widehat{M}}_{ij}$ because it was introduced by Majorana in 1937 [39] in the context of nuclear physics. Within the algebraic theory of molecules [13], it was introduced in 1991 [14]. This Hamiltonian is identical to that of the Bose–Hubbard model [40], except that, for applications in crystal vibrations, $\epsilon$ is large and positive, and U, V are small, $U,V\approx 0$. In this article, we take $U=V=0$ and consider the simplified form

$$\begin{array}{c}H=\epsilon {\sum }_{i=1}^n{b}_i^{†}{b}_i-{\lambda }^{\left(I\right)}{\sum }_{\langle i,j\rangle }\left(b_i^{†}{b}_j+b_j^{†}{b}_i\right)-{\lambda }^{\left(II\right)}{\sum }_{\langle \langle i,j\rangle \rangle }\left(b_i^{†}{b}_j+b_j^{†}{b}_i\right)+...,\hfill \end{array}$$

(4)

where we have expanded the Majorana interaction into one between nearest neighbors with interaction coefficient ${\lambda }^{\left(I\right)}$, and one between next-to-nearest neighbors with coefficient ${\lambda }^{\left(II\right)}$, etc. The coefficients ${\lambda }^{\left(I\right)},{\lambda }^{\left(II\right)},...$ are symmetry adapted to reflect the symmetry of the unit cell. Our sign convention for these coefficients is that of Reference [12], and they are assumed to be positive.

In a previous publication [12], we have studied the energy eigenvalues of the Hamiltonian Equation (4) and, from those, we have obtained the EDR and the DOS for linear (1D), square and hexagonal (2D) lattices, with symmetry $C_{\infty v}$, $D_{4h}$, and $D_{6h}$, respectively. In this article, we study the corresponding eigenfunctions. In particular, we consider here the fundamental vibrations of the lattice, obtained by diagonalizing the Hamiltonian Equation (4) in the n-dimensional basis

$$\left|1_i\right.〉\equiv \left|0,0,...,1_i,...,0,0\right.〉,$$

(5)

where $\left|1_i\right.〉$ denotes the state with one quantum of vibration (phonon) at location i. The wave functions in this basis, called the local basis, are written as

$$\left|{\psi }_{\nu }\right.〉=\sum _i{c}_i^{\left(\nu \right)}\left|1_i\right.〉,i=1,...,n\nu =1,...,n.$$

(6)

3. Square Lattice

The symmetry adaptation of this lattice, Figure 1, was discussed in [12]. The basis vectors ${\mathit{a}}_i$ and ${\mathit{b}}_i$ generating the lattice and the reciprocal lattice, respectively, are given as

$${\mathit{a}}_1=a(1,0),{\mathit{a}}_2=a(0,1),{\mathit{b}}_1=\frac{2\pi }{a}(1,0){\mathit{b}}_2=\frac{2\pi }{a}(0,1).$$

(7)

Within a unit cell, with symmetry $D_{4h}$, we have two types of interactions, nearest-neighbor, ${\lambda }^{\left(I\right)}$, and next-to-nearest neighbor, ${\lambda }^{\left(II\right)}$. The coordinates of the four nearest neighbors of the atom at lattice site i in Figure 1 are given as

$$\frac{{\mathit{\delta }}_1}{a}=-{\mathit{a}}_1,\frac{{\mathit{\delta }}_3}{a}={\mathit{a}}_2,\frac{{\mathit{\delta }}_5}{a}={\mathit{a}}_1,\frac{{\mathit{\delta }}_7}{a}=-{\mathit{a}}_2,$$

(8)

those of the next-to-nearest neighbors by

$$\frac{{\mathit{\delta }}_2}{a}=-{\mathit{a}}_1+{\mathit{a}}_2,\frac{{\mathit{\delta }}_4}{a}={\mathit{a}}_1+{\mathit{a}}_2,\frac{{\mathit{\delta }}_6}{a}={\mathit{a}}_1-{\mathit{a}}_2,\frac{{\mathit{\delta }}_8}{a}=-{\mathit{a}}_1-{\mathit{a}}_2.$$

(9)

The EDR including both interactions is given by [9]

$$\begin{array}{cc}\hfill E\left(\mathit{k}\right)& =\epsilon -{\lambda }^{\left(I\right)}{\sum }_{\mathrm{odd}j}{e}^{i\mathit{k}·{\mathit{\delta }}_j}-{\lambda }^{\left(II\right)}{\sum }_{\mathrm{even}j}{e}^{i\mathit{k}·{\mathit{\delta }}_j}\hfill \\ \hfill & =\epsilon -{\lambda }^{\left(I\right)}2\left[cos2\pi k_{1}a+cos2\pi k_{2}a\right]\hfill \\ \hfill & -{\lambda }^{\left(II\right)}2\left[cos2\pi (k_1+k_2)a+cos2\pi (k_1-k_2)a\right].\hfill \end{array}$$

(10)

Here, $k_1$ and $k_2$ are the components of $\mathit{k}$ in the reciprocal lattice basis, $\mathit{k}=k_1{\mathit{b}}_{\mathbf{1}}+k_2{\mathit{b}}_{\mathbf{2}}$, and the sign convention for ${\lambda }^{\left(I\right)},{\lambda }^{\left(II\right)}$ is that of [12,41]. Accordingly, in terms of the components $k_x$ and $k_y$, the EDR is given by [12,41,42]

$$\begin{array}{ccc}\hfill E\left(k_x,k_y\right)=\epsilon & -& {\lambda }^{\left(I\right)}2\left(cos{k}_{x}a+cos{k}_{y}a\right)\hfill \\ \hfill & -& {\lambda }^{\left(II\right)}4cos{k}_{x}acos{k}_{y}a,\hfill \\ \hfill 0& \le & k_{x}a,k_{y}a\le \pi .\hfill \end{array}$$

(11)

In Figure 2, we show the corresponding energy surface (left) and the plot of isofrequency contours (right) in the quasimomentum plane $(k_x,k_y)$.

The DOS is obtained from Equation (11) by integration,

$$g\left(E\right)=\frac{a^2}{{\pi }^2}{\int }_0^{\pi /2}\mathrm{d}{k}_x{\int }_0^{\pi /2}\mathrm{d}{k}_y\delta (E-E\left(\mathit{k}\right)),$$

(12)

where $\delta \left(\right)$ is Dirac’s $\delta$-function. For ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\le 1/2$ it is given by [12]

$$\begin{array}{ccc}\hfill g\left(\tilde{E}\right)& =& \frac{4K\left(k_1^2\right)}{{\pi }^2\sqrt{{\left(1+2{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\right)}^2-8\left({\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\right)\tilde{E}}},\hfill \\ \hfill k_1^2& =& \frac{4(1-\tilde{E})\tilde{E}}{{\left(1+2{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\right)}^2-8\left({\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\right)\tilde{E}},\hfill \\ \hfill 0& <& \tilde{E}<1,\hfill \end{array}$$

(13)

where $\tilde{E}=E/8{\lambda }^{\left(I\right)}+1/2$ and $K\left(k_1^2\right)$ is a complete elliptic integral of the first kind. For ${\lambda }^{\left(II\right)}=0$, i.e., only nearest-neighbor interactions, the DOS takes the simple form

$$g\left(\tilde{E}\right)=\frac{4K\left(k_1^2\right)}{{\pi }^2},k_1^2=4(1-\tilde{E})\tilde{E}.$$

(14)

Both in Equations (13) and (14), the DOS has a logarithmic singularity, called a van Hove singularity [23] (vHs). The upper panels in Figure 3 show the DOS for six different values of ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}$ with ${\lambda }^{\left(I\right)}=1$ fixed for a finite-size square shaped square lattice with $n_1\times n_2=42\times 42$ sites and Dirichlet boundary conditions. Equation (11) provides the EDR of an infinitely extended square lattice. For a finite-size sheet, the values of $(k_1,k_2)$ and $(k_x,k_y)$ are discretized, depending on the shape of the sheet. Numerically, the determination of the eigenvalues and wave functions of a lattice with N atoms involves the diagonalization of an N-dimensional matrix. The eigenvector component $c_{i,i^{\prime }}^{(\kappa ,{\kappa }^{\prime })}$, with $(i,i^{\prime })$ and $(\kappa ,{\kappa }^{\prime })$ denoting the discretized coordinates $(x,y)$ and quasimomenta $(k_x,k_y)$, respectively, yields the amplitude of the wave function associated with the eigenvalue at the corresponding lattice site [43].

We also determined the quasimomenta corresponding to the eigenstates. For this, we computed the momentum distributions, which were obtained from the Fourier transform of the wave functions ${\psi }_n(x,y)$,

$${\tilde{\psi }}_n(k_x,k_y)={\int }_{\Omega }\mathrm{d}x\mathrm{d}y{\psi }_n(x,y)e^{-i(k_{x}x+k_{y}y)}.$$

(15)

It is peaked at the quasimomenta $(k_x,k_y)$ associated with the state ${\psi }_n(x,y)$.

We display examples for these wave functions and the associated momentum distributions in Figure 4. We consider to this end the case of a square shaped lattice with $n_x=n_y=42$ sites, which was cut out parallel to the basis vectors ${\mathit{a}}_{1,2}$ from an infinitely extended square lattice. The total number of states is here $n_x\times n_y=1764$. In Figure 4, we show the squared wave functions, i.e., the intensities ${\left|c_{i,i^{\prime }}^{(\kappa ,{\kappa }^{\prime })}\right|}^2$, obtained numerically by diagonalizing the Hamiltonian Equation (4). At the lower band edge, shown in the first row of Figure 4, the wave functions represent vibrations of the lattice as a membrane. Around the vHs at $\nu =882$, shown in the second row of Figure 4, we have the occurrence of some unusual features consisting of stripes, particularly evident for the wave function $\nu =880$. This spatial localization will be discussed further in Section 5. At the upper the band edge the wave functions represent again membrane-like vibrations of the lattice. These wave functions, with the index counted from the upper band edge down, differ from those counted from the lower band edge up only in the signs of their components. The figures also include momentum distributions. They are peaked around the minimum and the maximum of the EDR for eigenvalues close to the band edges, and along the equipotential contour connecting the saddle points of the EDR around the van Hove singularities.

The results in Figure 4 are for a square lattice with square-shape and Dirichlet boundary conditions. Similar studies have been done for square lattices with a rectangular shape of the boundary. Furthermore, we considered square and rectangular shaped sheets, which were cut out of a square lattice not along the rows of the square lattice but along diagonal lines, i.e., along next-to-nearest neighbor atoms. For a square-shaped boundary, the wave function patterns are similar to those shown in Figure 4 at the band edges. This is in accordance with the findings of Reference [33], that, at the band edges, the spectral properties of a lattice only depend on the shape of the sheet. At the van Hove singularities, we again observe the occurrence of stripes in the wave function patterns; however, here they are parallel to the sides of the sheet, not to the diagonal as in Figure 4. In conclusion, the striped structure of the wave functions observed at the van Hove singularities persists independently of the shape of the sheet.

4. Hexagonal Lattice

The spectral properties of the hexagonal lattice are of particular interest in applications to graphene [2] and artificial graphene [22,33,44]. The symmetry adaptation of this lattice, Figure 5, was discussed in [12]. The unit cell can also be viewed as composed of two triangular subcells A and B. The basis vectors ${\mathit{a}}_i$ and ${\mathit{b}}_i$ generating each triangular lattice and the reciprocal lattice, respectively, are given as

$$\begin{array}{c}\hfill {\mathit{a}}_1=\frac{a}{2}(3,\sqrt{3}),{\mathit{a}}_2=\frac{a}{2}(3,-\sqrt{3}),\end{array}$$

(16a)

$$\begin{array}{c}\hfill {\mathit{b}}_1=\frac{2\pi }{3a}(1,\sqrt{3}),{\mathit{b}}_2=\frac{2\pi }{3a}(1,-\sqrt{3}),\end{array}$$

(16b)

with a denoting the distance between neighboring atoms. Within the unit cell, with symmetry $D_{6h}$, we have three types of interactions, nearest neighbor, ${\lambda }^{\left(I\right)}$, next-to-nearest neighbor, ${\lambda }^{\left(II\right)}$, and third neighbor, ${\lambda }^{\left(III\right)}$. The interaction ${\lambda }^{\left(II\right)}$ couples atoms within the same subcell, while ${\lambda }^{\left(I\right)}$ and ${\lambda }^{\left(III\right)}$ couple atoms in subcell A with those in subcell B. The coordinates of the three nearest neighbors of the atom at lattice site i in Figure 5 are given as

$${\mathit{\delta }}_1=-\frac{1}{3}({\mathit{a}}_1+{\mathit{a}}_2),{\mathit{\delta }}_5={\mathit{a}}_1+{\mathit{\delta }}_1,{\mathit{\delta }}_9={\mathit{a}}_2+{\mathit{\delta }}_1.$$

(17)

Since the hexagonal lattice is composed of two interpenetrating triangular lattices, one has to diagonalize a $2\times 2$ matrix

$$\left|\begin{array}{cc}{H}_{AA}-E& H_{AB}\\ H_{BA}& H_{BB}-E\end{array}\right|=0,$$

(18)

in order to obtain the EDR for an infinite-size hexagonal lattice [4]. Taking into account only nearest-neighbor interactions, the off-diagonal elements are given as

$$\begin{array}{ccc}\hfill H_{AB}& =& {\lambda }^{\left(I\right)}{\sum }_{j=1,5,9}{e}^{i\mathit{k}·{\mathit{\delta }}_j}\hfill \\ & =& e^{i\mathit{k}·{\mathit{\delta }}_1}\left(1+e^{i\mathit{k}·{\mathit{a}}_1}+e^{i\mathit{k}·{\mathit{a}}_2}\right)\hfill \\ \hfill H_{BA}& =& H_{AB}^{\star }.\hfill \end{array}$$

(19)

Choosing $H_{AA}=H_{BB}=0$, the diagonalization of Equation (18) yields for the EDR in terms of the components $k_1$ and $k_2$ of $\mathit{k}$ in the reciprocal lattice basis

$$\begin{array}{c}E\left(k_1,k_2\right)=\mp {\lambda }^{\left(I\right)}\sqrt{3+2\left[cos2\pi a{k}_1+cos2\pi a{k}_2+cos2\pi a(k_1-k_2)\right]},\hfill \end{array}$$

(20)

and with $k_1=\left(3k_x+\sqrt{3}{k}_y\right)/4\pi$, $k_2=\left(3k_x-\sqrt{3}{k}_y\right)/4\pi$ in the $(k_x,k_y)$ plane

$$\begin{array}{ccc}\hfill E\left(k_x,k_y\right)& =& \mp {\lambda }^{\left(I\right)}\sqrt{3+u(k_x,k_y)},\hfill \\ \hfill u(k_x,k_y)& =& 2cos\sqrt{3}a{k}_y+4cos\frac{\sqrt{3}}{2}a{k}_{y}cos\frac{3}{2}a{k}_x.\hfill \end{array}$$

(21)

The EDR including the next-to-nearest neighbor interactions is given by [2,12]

$$E\left(k_x,k_y\right)=\mp {\lambda }^{\left(I\right)}\sqrt{3+u(k_x,k_y)}-{\lambda }^{\left(II\right)}u(k_x,k_y).$$

(22)

Figure 6, left panel, shows the two sheets of the energy surface, Equation (21), of an infinite size hexagonal lattice, in the quasimomentum plane $(k_x,k_y)$, where $k_{x}a$ and $k_{y}a$ are varied within the first Brillouin zone. The energy surfaces touch each other conically at the corners of the first Brillouin zone. The right panel shows the corresponding plot of isofrequency contours in the quasimomentum plane $(k_x,k_y)$.

The DOS associated with Equations (20) and (21) can be written in explicit analytic form [45] in terms of an elliptic integral

$$\begin{array}{ccc}\hfill g\left(\omega \right)& =& \frac{9}{{\pi }^2}\sqrt{\frac{1-2\omega }{3}}K\left(\sqrt{\frac{\omega {\left(2-3\omega \right)}^3}{\left(1-2\omega \right)}}\right),\hfill \\ \hfill 0& <& \omega <\frac{1}{3},\hfill \\ \hfill g\left(\omega \right)& =& \frac{9}{{\pi }^2}\frac{\left(1-2\omega \right)}{\sqrt{3\omega {\left(2-3\omega \right)}^3}}K\left(\sqrt{\frac{\left(1-2\omega \right)}{\omega {\left(2-3\omega \right)}^3}}\right),\hfill \\ \hfill \frac{1}{3}& \le & \omega \le \frac{1}{2},\hfill \end{array}$$

(23a)

where $\omega =\left(\tilde{E}+1\right)/2$, $\tilde{E}=E/3{\lambda }^{\left(I\right)}$, $-1<\tilde{E}<1$. This function has two van Hove singularities (vHs${}^{\pm }$) at $\tilde{E}=\pm 1/3$ and a Dirac point (DP) at $\tilde{E}=0$. Equations (21) and (22) provide the EDR of an infinite-size hexagonal lattice. We are interested in wave function properties of finite sheets. The corresponding quasimomentum components $(k_x,k_y)$ take discrete values $(\kappa ,{\kappa }^{\prime })$, which depend on the shape of the sheet. Furthermore, for finite-size lattices, the DOS has an additional peak at $E=0$. The upper panels in Figure 7 show the DOS of a hexagonal lattice with Africa shape [22,34] and Dirichlet boundary conditions for six different values of ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}$ with ${\lambda }^{\left(I\right)}=1$ fixed and ${\lambda }^{\left(III\right)}=-{\lambda }^{\left(II\right)}/2$. The peak at $E\simeq 0$ is attributed to the edge states that are localized along the zigzag edges in the lattice structure. Note that the positions of the van Hove singularities vHs${}^{\pm }$ and also of the DP are slightly shifted with increasing ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}$ as expected from tight-binding model calculations [46].

For the general case of a sheet with N sites, the energies and wave functions can be obtained by diagonalizing an N-dimensional matrix. We show in Figure 8 the wave functions of a hexagonal lattice with rectangular shape with $N=1656$ sites and Dirichlet boundary conditions, including both space and momentum distributions. The zigzag edges are parallel to the longer sides of the rectangle and the lattice consists of in total 24 rows of 34 and 35 lattice sites, respectively. At the lower band edge, see the first row of Figure 8, the wave functions represent vibrations of the lattice as a membrane, similarly as in the square lattice. Around the vHs${}^{-}$, shown in the second row of Figure 8, the wave functions exhibit a striking structure consisting of stripes along interior zigzag lines. This behavior, already seen in the square lattice, is even more pronounced here. Around the DP, see the third row of Figure 8, we observe the above mentioned edge states. Their wave functions are concentrated at the zigzag edges of the lattice, and are vanishingly small everywhere else. The situation repeats at the vHs${}^{+}$, with clear stripes, and at the upper band edge, with membrane-like vibrations. The only difference between the wave functions at the lower and upper band edges is the sign of wave-function components at the sites of the lattice, $c_{i,i^{\prime }}^{(\kappa ,{\kappa }^{\prime })}$. The fourth to sixth rows in Figure 8 show the corresponding momentum distributions. Near the lower band edge (fourth row), they are peaked around the $\Gamma$ point of the EDR, on equienergy lines of the M points close to the vHs (see fifth row) and around the K points in the vicinity of the DP (sixth row).

The results in Figure 8 are for a rectangular shape of the lattice. In view of the remarkable properties of the wave functions at the van Hove singularities, we have also investigated the wave functions for more complicated shapes, in particular the Africa shape [22,34,47,48,49]. We show in Figure 9 again the wave functions at the lower band edge (first row), the vHs${}^{-}$ (second row), and the DP (third row) for Dirichlet boundary conditions. For such shapes, we also observe the extraordinary features described in the paragraphs above. In addition, at both van Hove singularities, the stripes appear along all three equivalent zigzag directions. Two of these, inclined at angles $\pm 2\pi /3$ to the horizontal axis, were not seen in the rectangular shape due to boundary effects. In addition, the momentum distributions are extended to the whole hexagonal equienergy contour corresponding to the vHs, in contrast to the rectangular shape, fifth row in Figure 8.

5. Localization Properties

We illustrated in the previous sections that the wave functions of the square and hexagonal lattice are localized into stripes at the van Hove singularities and for the latter in addition at the zigzag edges around the DP. In this section, we quantify the localization in terms of participation ratios [50,51] in the local basis, defined for 2D lattices through the wave function components $c_{i,i^{\prime }}^{(\kappa ,{\kappa }^{\prime })}$ at sites $(i,i^{\prime })$ with $(\kappa ,{\kappa }^{\prime })$ denoting the discretized quasimomenta $(k_x,k_y)$ as

$${{\mathrm{PR}}_{\mathrm{loc}}}^{(\kappa ,{\kappa }^{\prime })}=\frac{1}{{\sum }_{i{i}^{\prime }}{\left|c_{i,i^{\prime }}^{(\kappa ,{\kappa }^{\prime })}\right|}^4}.$$

(24)

Similar results are obtained for the Shannon entropy [52]. Note that, in order to reproduce the DOS of real graphene, nonzero second- and third-neighbor interactions had to be introduced [9,53]. Accordingly, we are particularly interested here in how localization is influenced by (i) the addition of second- and third-neighbor interactions, and (ii) noise affecting the hopping coefficients ${\lambda }^{\left(I\right)}$, ${\lambda }^{\left(II\right)}$ and ${\lambda }^{\left(III\right)}$. For this, we keep ${\lambda }^{\left(I\right)}=1$ fixed and vary ${\lambda }^{\left(II\right)}$ and ${\lambda }^{\left(III\right)}$. Since, in general, in a square lattice ${\lambda }^{\left(III\right)}\ll {\lambda }^{\left(I\right)}$, we can neglect it in the study of their localization properties. In hexagonal lattices, on the other hand, ${\lambda }^{\left(III\right)}\simeq -{\lambda }^{\left(II\right)}/2$ [33], so, in the corresponding studies, we simply chose ${\lambda }^{\left(III\right)}=-{\lambda }^{\left(II\right)}/2$ and ${\lambda }^{\left(I\right)}$ and ${\lambda }^{\left(II\right)}$ with opposite sign [54].

Figure 3 shows the DOS (upper panels) and the corresponding participation ratio of the fundamental vibration of a square lattice with Dirichlet boundary conditions and square shape for six different values of ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}$. Here, the values of the hopping coefficients ${\lambda }^{\left(I\right)},{\lambda }^{\left(II\right)}$ are the same for each lattice site. We see from this figure that around the vHs-energy $E\approx 0$, states of increased localization persist even for large values of the relative strength ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\le 1.0$ of the second-neighbor interactions, identified by small values of ${\mathrm{PR}}_{\mathrm{loc}}$ around the singularity. However, a close inspection reveals that the minimum ${\mathrm{PR}}_{\mathrm{loc}}$ jumps abruptly from ${\mathrm{PR}}_{\mathrm{loc}}\approx 76$ at ${\lambda }^{\left(II\right)}=0$ to a value saturated at ${\mathrm{PR}}_{\mathrm{loc}}\approx 410$, for an infinitesimally small value ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}>0$, as checked numerically down to ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}={10}^{-10}$. The corresponding striped eigenstate for ${\lambda }^{\left(II\right)}=0$ is nonvanishing only along one of the diagomals and the two neighboring ones. For ${\lambda }^{\left(II\right)}>0$, it spreads into states resembling the eigenstates of a lattice with pure second-neighbor interaction, i.e., ${\lambda }^{\left(II\right)}\ne 0,{\lambda }^{\left(I\right)}=0$ (not shown here), in which the vibration involves only one half of the lattice sites. This is reflected in the fact that the saturation value of ${\mathrm{PR}}_{\mathrm{loc}}$ at the vHs for ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}>0$ equals one half of the value corresponding to the ground state, which is ${\mathrm{PR}}_{\mathrm{loc}}^{\mathrm{g}.\mathrm{s}.}\approx 822$, indicating that the wave function components are nonzero on one half of the sites. Accordingly, stripes are not present in the square lattice for ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}>0$ independently of the shape of the boundary.

Figure 7 shows the participation ratios of the fundamental vibrations of the hexagonal lattice with an Africa shape. The corresponding results for a rectangular shape are not shown because they are qualitatively similar to those for the Africa shape. As stated above, we included third-nearest neighbor interactions in this case where we chose ${\lambda }^{\left(III\right)}=-{\lambda }^{\left(II\right)}/2$. Similarly to the situation in the square lattice discussed above, this figure indicates remarkable localization properties. In particular, the minimum value of ${\mathrm{PR}}_{\mathrm{loc}}$ at the DP ($E_{\mathrm{DP}}=3/2{\lambda }^{\left(II\right)}$) is practically independent of $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|$, underlining the robustness of the edge states. However, the minimum values of ${\mathrm{PR}}_{\mathrm{loc}}$ near the lower and the upper van Hove singularities ($E_{\mathrm{vHs}}^{\pm }=\pm {\lambda }^{\left(I\right)}+{\lambda }^{\left(II\right)}$) show a jump at small values of $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\lesssim 0.001$, from ${\mathrm{PR}}_{\mathrm{loc}}\approx 124$ to ${\mathrm{PR}}_{\mathrm{loc}}\approx 238$ for the Africa shape and from ${\mathrm{PR}}_{\mathrm{loc}}\approx 140$ to ${\mathrm{PR}}_{\mathrm{loc}}\approx 500$ for the rectangular one. These values for ${\mathrm{PR}}_{\mathrm{loc}}$ already indicate that the localization as compared to the ground state, where ${\mathrm{PR}}_{\mathrm{loc}}$ equals half the number of sites, is stronger in the former case than in the latter. In contrast to the square lattice, the peculiar striped states persist for the van Hove singularities even though they are broader above $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\simeq 0.001$ as indicated by the respective participation ratio. Above $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\simeq 0.2$, a third peak appears in the DOS near the upper band edge. There, at the vHs${}^{+}$, the stripes disappear with increasing $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|$ and the wave functions become more and more localized at the armchair edges, that is, e.g., in the hexagonal lattice with rectangular shaped at the shorter sides. This is illustrated in the upper rows of Figure 10 and Figure 11.

Examples are shown for wave functions at minimal values of the participation ratio close to the vHs${}^{-}$, the vHs${}^{+}$ and the peak in the DOS appearing above $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\simeq 0.2$ for $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|=0.25$. In both cases, the wave functions are localized in broadened stripes along interior zigzag edges at the vHs${}^{-}$, strongly localized at armchair edges at the vHs${}^{+}$ and within the lattice at the third peak in the DOS. Note that this feature at the vHs${}^{+}$ is opposite to that near the DP, where the wave functions are localized at the zigzag edges. This persistence of the striped structure at vHs${}^{-}$ and its disappearance at the vHs${}^{+}$ and the simultaneous increase of the localization at the armchair edges is reflected in the corresponding participation ratios.

Substantially different behaviors between the square and the hexagonal lattices are observed under the influence of noise. We introduce the noise by generating individual hopping coefficients ${\lambda }_{ij}^{\left(n\right)}$, $n=I,II,III$ at each site $i,j$ in the general Hamiltonian of Equation (3) randomly with a uniform distribution, yielding ${\lambda }_{ij}^{\left(n\right)}={\lambda }^{\left(n\right)}(1+{\delta }_{ij}^{\left(n\right)})$, where we again set ${\lambda }^{\left(III\right)}=-{\lambda }^{\left(II\right)}/2$. Here, ${\lambda }^{\left(n\right)}$ denotes the average value across the whole lattice and $|{\delta }_{ij}^{\left(n\right)}|$ denotes the noise strength with respect to ${\lambda }^{\left(n\right)}$ at site $(i,j)$.

In Figure 12, we show the evolution of the participation ratios of the vHs-states compared to that of the ground states (g.s.) of the square lattice as a function of ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}$ without noise (left) and with noise of relative strength $|{\delta }_{ij}^{\left(n\right)}|\le {10}^{-2}$ (right). Note that the state with ${\mathrm{PR}}_{\mathrm{loc}}\approx 76$ in the noiseless case at ${\lambda }^{\left(II\right)}=0$ is due to the formation of stripes in the wave function structure (see Section 5). Evidently, noise leads to an enhancement of localization at the vHs for ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\lesssim 0.03$. For ${\lambda }^{\left(II\right)}=0$, the wave function is nonvanishingly only along one of the diagonals, as reflected in the corresponding participation ratio, ${\mathrm{PR}}_{\mathrm{loc}}\approx 51$. With increasing ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}$, the eigenstate at the vHs delocalizes not abruptly, but smoothly until it saturates at ${\mathrm{PR}}_{\mathrm{loc}}\approx 411$ for ${\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}\gtrsim 0.03$.

Similarly to the square, the noise in the hexagonal lattice acts against, and surpasses, the effect of delocalization caused by the second- and third-neighbor interaction, as can be seen in Figure 13 and Figure 14. Compared to the noiseless situation with an abrupt jump for $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|>0$ (left panels), in the noisy case, increasing $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|$ induces instead a gradual rise in ${\mathrm{PR}}_{\mathrm{loc}}$ (right panels). In contrast to the square lattice, the localization of the striped states at the vHs${}^{-}$ is enhanced in the noiseless case even for $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\gtrsim 0.03$, where saturation of the participation ratio takes place. The minimal participation ratio at the vHs${}^{+}$ gradually decreases with increasing $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|$ for the rectangular shape both with and without noise until it saturates at around $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\simeq 0.24$ similar to that at the third peak appearing in the DOS above $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\simeq 0.2$. The localization of the latter is considerably enhanced when turning on the noise. Similar results are obtained for the Africa shape; however, there the minimal participation ratio at the vHs${}^{+}$ starts to decrease only at $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|\simeq 0.2$, simultaneously with the occurrence of the third peak in the DOS at the upper band edge.

In all, the localization is stronger in the hexagonal lattice with the shape of a classically chaotic Africa billiard than in that with the shape of a rectangle with classically integrable dynamics. From these observations, we may conclude that both the introduction of noise and the choice of the shape of a billiard with classically chaotic dynamics enhance the localization at the peaks in the DOS. The comparison of the wave functions in the upper and lower rows of Figure 10 and Figure 11 clearly corroborates this behavior.

Finally, in order to illustrate the main object of study of this article, i.e., the striped structure in the wave functions, we show in Figure 15, the smoothed probability density $|c^{(\kappa ,{\kappa }^{\prime })}{(x,y)|}^2$ corresponding to such vibrations at the van Hove singularities of hexagonal lattices with rectangular (left) and Africa-shaped (right) boundaries. This figure shows clearly the ripples of flexural modes in graphene-like materials.

We also note that physically relevant values of $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|$ are $1/27\simeq 0.04$ for a Lennard–Jones potential, $\approx 1/22$ for artificial graphene [12] and $\approx 0.22$ for graphene [53], thus we expect that the effects discussed above can be realized experimentally. Although we do not report it, these localization properties persist also for other shapes of the boundary, like sectors from a circle [55].

In our previous papers [12,22], a connection was made between van Hove singularities and excited state quantum phase transitions (ESQPT) whose definition is given in [56,57]. In relation to this connection, we note that the localization properties of the $U\left(2\right)\otimes U\left(2\right)\otimes ...\otimes U\left(2\right)$ model described here appear to be the same as those of $U\left(n\right)$ models, described in [58]. The explicit relation between van Hove singularities (and their associated striped wave functions) and ESQPT will be discussed in a forthcoming publication.

6. Discussion

In this paper, the wave functions of square and hexagonal lattices (2D) have been studied. The wave functions in 2D show remarkable properties at the van Hove singularities and at the Dirac point. Particularly interesting are the properties at the van Hove singularities consisting of stripes.

In case of pure nearest-neighbor interactions, the striped states occur at van Hove singularities in both the square and the honeycomb lattices. However, while in the square lattice, the striped states disappear rapidly, when turning on second-nearest neighbor coupling ${\lambda }^{\left(II\right)}$, and/or noise, in the hexagonal lattice, the stripes—oriented parallel to the zigzag edges—persist to realistic values of $|{\lambda }^{\left(II\right)}/{\lambda }^{\left(I\right)}|$ and become even more localized under the influence of weak noise. The orientations of the stripes in finite hexagonal systems depend on the shape of the boundary, as we have demonstrated by the rectangular and Africa shapes with regular and chaotic spectral properties, respectively. Localization may be enhanced by either introducing noise or chaoticity or, even better, both. The results of Section 5 explicitly show that graphene shares many properties of soft membranes [5] (Figure 8).

7. Conclusions

Although here we have concentrated on systems with Dirichlet boundary conditions and 1d vibrations, the AMC [12] is quite general and can be used also for systems with cyclic boundary conditions. These are of interest in graphene ribbons [9] and carbon nanotubes (1D systems), and finally in fullerene [10] (0D). The AMC can also be extended to 2d vibrations both in 1D and 2D systems. This is done by inserting at each site the algebra $g\equiv u\left(3\right)$ instead of $g\equiv u\left(2\right)$.

In addition, we mention that, since the Hamiltonian (3) is formally identical for bosons and fermions — the only difference being that for boson systems, we can put any number at each site, while, for fermion systems, we can put only one (or two if spin is included)—the same properties of the wave functions discussed here will appear in the equivalent fermion case, i.e., in the electronic band structure of graphene. Finally, because of its generality, the AMC can be used for other two-dimensional materials beyond graphene by adapting the boson operators of Section 2 to the appropriate symmetry of the lattice.