A Polynomial Approach to the Spectrum of Dirac-Weyl Polygonal Billiards
aa r X i v : . [ phy s i c s . c o m p - ph ] J u l A Polynomial Approach to the Spectrum ofDirac-Weyl Polygonal Billiards
M. F. C. Martins Quintela ∗ and J. M. B. Lopes dos Santos3rd August 2020 Centro de Física das Universidades do Minho e do Porto, CF-UM-UP e Departmento de Física eAstronomia, Universidade do Porto, Rua do Campo Alegre 4169-007, Porto, Portugal.The Schrödinger equation in a square or rectangle with hard walls is solved in everyintroductory quantum mechanics course. Solutions for other polygonal enclosures onlyexist in a very restricted class of polygons, and are all based on a result obtained byLamé in 1852. Any enclosure can, of course, be addressed by finite element methodsfor partial differential equations. In this paper, we present a variational method toapproximate the low-energy spectrum and wave-functions for arbitrary convex polygonalenclosures, developed initially for the study of vibrational modes of plates. In view of therecent interest in the spectrum of quantum dots of two dimensional materials, describedby effective models with massless electrons, we extend the method to the Dirac-Weylequation for a spin-1/2 fermion confined in a quantum billiard of polygonal shape, withdifferent types of boundary conditions. We illustrate the method’s convergence in caseswhere the spectrum in known exactly and apply it to cases where no exact solution exists.
One of the first problems solved in introductory quantum mechanics courses is that of particleenclosed in a square or rectangular box with hard walls—Dirichlet boundary conditions, ψ ( r ) = 0 at the boundary. Separation of variables easily turns the problem into two one-dimensional problemsand the solutions are a finite sum of plane waves. The circle and the ellipse are also amenable to anexact treatment, again with separation of variables. But the vast majority of students are never ledto consider the case of general polygonal enclosures. ∗ [email protected] The essential idea of the polynomial method [18, 19, 20] is to generate orthonormal polynomials ofincreasing order that satisfy, from the start, the required boundary conditions (BC). In the
Oxy x and y , P ( x, y ) = a + a x + a y will be zero in a straight line, a + a x + a y = 0 . It follows that a product of n such polynomials trivially satisfies the condition P n ( x, y ) = 0 at the boundary of a polygonal convex region. After normalization, such a function isa zero-th order approximation to the ground state.As an example, for a right triangle with vertexes (0 , , (1 , , (0 , we would choose ψ ( x, y ) = N xy (1 − x − y ) (1)with N = (cid:20) ˆ dx ˆ − x dy [ xy (1 − x − y )] (cid:21) − / = 12 √ If we multiply ψ ( x, y ) by any polynomial, ψ ( x, y ) = P ( x, y ) ψ ( x, y ) , we end up with a functionthat still satisfies the BC. We therefore proceed to generate a set of linearly independent functions ψ n ( x, y ) = P n ( x, y ) ψ ( x, y ) (2)where {P n ( x, y ) : n = 0 , , . . . } is a family of ascending-order monomials/polynomials. These canbe generated by taking into account the symmetries of the enclosure in question, as a simple orderingof monomials, or as a generalization of Legendre polynomials discussed by Larcher [21]. In [20] theauthors resort to the sorting (cid:8) , x, y, xy, x , y , x y, y x, x y , x , ... (cid:9) . (3)We truncate the series at a finite n and orthonormalize the functions, typically through Gram-Schmidt orthogonalization. The Hamiltonian matrix in this truncated basis is calculated and di-agonalized. For regular polygons, one can take into account the symmetry group of the enclosureby replacing the monomials of Eq.(3) by increasing order polynomials that transform according toirreducible representations of the symmetry group [22]. The Hamiltonian matrix becomes blockdiagonal, and the resulting eigenstates end up classified according to their symmetry properties.The problem of the Schrödinger equation for a particle confined to a square well, provides asuitable example to assess convergence. The exact spectrum of the Schrödinger equation for a particle confined to a square well is E n ,n = ~ m π L (cid:0) n + n (cid:1) , (4)where L is the side length of the enclosure, and n positive integers.3entering the enclosure at the origin, the starting wavefunction is ψ ( x, y ) = N "(cid:18) L (cid:19) − x L (cid:19) − y . (5)The higher order basis functions are generated as described previously, by multiplication by x m y n monomials. Both the orthogonalization and the calculation of the matrix elements of the kineticenergy operator reduce to the calculation of integrals of x m y n monomials in the domain of the en-closure; these are easily calculated symbolically and stored. Keeping the computation symbolic upuntil the diagonalization of the Hamiltonian avoids stability problems of the Gram-Schmidt proced-ure. We computed the Hamiltonian matrix for different basis sizes, and diagonalized it numericallyusing standard packages. The results are shown in Fig.1, overlayed with the exact eigenvalues.The method is seen to converge towards the exact eigenvalues as we increase the basis size. Asthe index of the eigenvalue increases, we require a larger basis, as expected. The numerical resultsindicate, roughly, that to the get the lowest n eigenvalues to within 1% accuracy requires of theorder of n basis functions. Eigenvalue Index E n e r g y ( h _2 / m . π / L ) Polyn. Method (6 polyns.)Polyn. Method (15 polyns.)Polyn. Method (28 polyns.)Polyn. Method (45 polyns.)Polyn. Method (66 polyns.)Polyn. Method (91 polyns.)Exact Solution
Figure 1: Convergence of the eigenvalues of the Hamiltonian in a square enclosure (polynomialmethod versus exact solution). 4 .3 The hexagonal enclosure
To apply this procedure to the Schrödinger equation for a particle confined to an hexagonal well ofside length L , we define the initial normalized polynomial as ψ ( x, y ) = N (cid:18) L − y (cid:19) " L − (cid:18) x − y √ (cid:19) × " L − (cid:18) x + y √ (cid:19) . (6)We again calculated the Hamiltonian matrix elements for different basis sizes, and diagonalized itnumerically using standard packages. The results are shown in Fig.2, overlayed with the eigenvaluesobtained via direct numerical integration of the partial differential equation (PDE), by finite elementsmethod. The convergence is similar to the one found in the case of the square. Eigenvalue Index E n e r g y ( h _2 / m . π / L ) Polyn. Method. (6 polyns.)Polyn. Mehod (15 polyns)Polyn. Method (28 polyns.)Polyn. Method (45 polyns.)Polyn. Method (66 polyns.)Polyn. Method (91 polyns.)Numerical Integration
Figure 2: Convergence of the eigenvalues of the Hamiltonian in an hexagonal enclosure (polynomialmethod versus direct numerical integration of PDE).As an illustration of this method, we applied it to different polygons with increasing number ofsides and plot the lowest eigenvalues against the inverse number of sides (squared) in Fig.3, wherethe result for the circle, ǫ ∞ i , i = 1 , , , is the analytical solution.5 (Number of Sides) -2 E n e r g y ( h _2 / m . / L ) First EigenvalueSecond Eigenvalue (2x degenerate)Third Eigenvalue
Figure 3: Comparison of the first three distinct energy levels as a function of the inverse of thenumber of sides of the polygon squared.The dashed lines are fits of the form ǫ i ( n ) = ǫ ∞ i + a n − + a n − , where n is the number of sides.A suitable fit could not be found for the third eigenvalue, and the point for n = 3 was excluded.The remaining data points are well fitted with a = 0 . We now turn to the extension of these results to the Dirac-Weyl equation.
Berry and Mondragon [8] were the first to address the issue of the boundary conditions of a two-component spinor solution of the Dirac-Weyl equation in a confining planar enclosure. The eigenvalueequation is − i ~ v F σ · ∇ Ψ = E Ψ . (7)where σ = ( σ x , σ y ) . We assume a confining boundary at y = 0 , excluding the region y < . Theprobability current is given by j = v F Ψ † σ Ψ = v F h σ i , (8)which means that a necessary condition for confined states is h σ y i boundary = 0 . (9)6his condition implies that the spinor at the boundary must be an eigenstate of σ u = σ · u , where u = cos θ e z + sin θ e x . i.e ., (cid:20) cos θ sin θ sin θ − cos θ (cid:21) (cid:20) ψ A ψ B (cid:21) = η (cid:20) ψ A ψ B (cid:21) , η = ± , (10)or simply ψ B ψ A = η − cos θ sin θ . (11)Denoting t := tan ( θ/ , simple trigonometric identities lead us to ψ B ψ A = t, (12)with a spinor, at the boundary, given by (cid:20) ψ A ψ B (cid:21) = 1 √ t (cid:20) t (cid:21) (13)for η = +1 . The case for η = − is obtained by the substitution t → − /t , or by changing θ → π − θ .As such, this is not a different solution, and we keep only η = +1 .Some special cases are:• t = 0 or t → ∞ (Dirichlet boundaries) (cid:26) t = 0 , ψ B = 0; t → ∞ , ψ A = 0 ; (14)so that, at the boundary, h σ z i = ± .• t = ± ψ A = ± ψ B . (15)In this case, the spin points along the x -direction.The latter case corresponds to the infinite-mass confinement by Berry and Mondragon [8]. Theyconsidered adding a spatially dependent mass (gap) term, − i ~ v F σ · ∇ Ψ + m ( r ) σ z = E Ψ . (16)and showed that continuity of the spinor components at the boundary, with m ( r ) = 0 inside theenclosure, and m ( r ) → ∞ outside, leads to BC of the type of Eq.(15). Similarly, one can show thata mass term of the form [22] m ( r ) ( σ z − cos θ ) (17)7eads to the more general BC of Eq.(12) with t = tan ( θ/ . The cases t = 0 or t → ±∞ correspondto Dirichlet boundaries in one of the spinor components, as in the case of zigzag boundaries ingraphene.If the boundary is at an angle φ with the x -axis, one must rotate the spinor in Eq. 13, resultingfinally in Ψ = e − iφσ z / √ t (cid:20) t (cid:21) = 1 √ t (cid:20) e − iφ/ e iφ/ t (cid:21) . (18)We will now show how to extend the polynomial method defined in Section 2 to this new type ofBC. As an illustration we begin with 1D cases, where one easily computes the exact spectrum. Consider the simple case where t = 1 and the confined region is the interval y ∈ [ − L/ , L/ . TheBC are ψ B ( − L/ /ψ A ( − L/
2) = 1 at one end, and ψ B ( L/ /ψ A ( L/
2) = − , at the other ( φ = π in Eq.(18)). Measuring energies in units of ~ v F /L , ǫ = E/ ( ~ v F /L ) , − ∂ y ψ B = ǫψ A ∂ y ψ A = ǫψ B . (19)Squaring the Hamiltonian, one sees that the solutions have the form ψ B = f e iqy + f e − iqy ψ A = f e iqy + f e − iqy , (20)with ǫ = sq , s = ± . With no loss of generality we assume q > . The Dirac-Weyl equation implies f = isf .f = − isf The BC then imply (1 − ist ) f + (1 + ist ) f e − iqL = 0(1 + ist ) f + (1 − ist ) f e iqL = 0 . (21)A non-zero solution to this homogeneous system of 2 equations, for t = 1 , requires (1 − is ) e iqL = (1 + is ) e − iqL (22)or q = π L (2 n + 1) (23)and the energies are ǫ s ( n ) = s ~ v F π L (2 n + 1) n ≥ , s = ± (24)8he Dirac-Weyl equation in free space has electron-hole symmetry; it is easily seen that for anysolution Ψ = [ ψ A , ψ B ] T of energy E , the state Ψ = [ ψ A , − ψ B ] T is a solution of energy − E . Butthe BC we are choosing, ψ B /ψ A = ± , breaks this symmetry, because Ψ and Ψ cannot both obeythese BC. This can be traced to the mass term we introduced outside the confining region to derivethe BC: it explicitly breaks this particle-hole symmetry. Nevertheless, the spectrum of eigenvaluesremains symmetrical (Eq.24), because it is still true for this complete 1D Hamiltonian, inside andoutside the enclosure, that σ x is a chiral symmetry operator, σ x ˆ Hσ x = − ˆ H . -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 Eigenvalue index -20-18-16-14-12-10-8-6-4-202468101214161820 E n e r g y ( h _ v F π / L ) Polynomial Method (16 polynomials)Exact Solution
Basis Size E n e r g y ( h _ v F π / L ) First Eigenval.Second Eigenval.Third Eigenval.
Figure 4: Exact solution vs polynomial approximation for a 1-D strip with t = 1 . Convergence ofthe first three eigenvalues in the inset.To implement the polynomial method, we must first realize that the parities of ψ A and ψ B areopposite. We therefore choose for the spinor components as the lowest order polynomials that satisfythe BC. For the conduction band ( s = 1) we have, Ψ (1)0 = N (cid:20) − yL/ (cid:21) (25)where N is a normalization constant. To obtain the valence band solution, s = − , it is enoughto swap ψ A ↔ ψ B . The basis is now generated by functions of the form Ψ ( s ) n = P n ( y )Ψ ( s )0 , or-thonormalized by the Gram-Schmidt process. One must take care, however, to orthogonalize statesbetween the valence and conduction bands, as the states generated this way for s = 1 and s = − are not orthogonal to start with.Having generated the basis, one easily computes and diagonalizes the Hamiltonian matrix. Weobtained the approximate spectrum visible in Fig.4; the inset shows the convergence analysis of thefirst three eigenvalue of the conduction band. An almost exact match is also present in the first9igenfunction of the conduction band, as shown in Fig.5.Having checked the validity of this method against a simple solution of the Dirac-Weyl equation,we will now proceed to consider planar enclosures. -0,5 -0,4 -0,3 -0,2 -0,1 0 0,1 0,2 0,3 0,4 0,5 Adimensional Position ( y/L ) -0,8-0,6-0,4-0,200,20,40,60,81 Sp i no r P r obab ilit y A m p lit ud e s ( Ψ ) Ψ A (Exact Solution) Ψ B (Exact Solution) Ψ A (Approximated) Ψ B (Approximated) Figure 5: Comparison between the exact solution (full lines) and the lowest-energy polynomial ei-genfunction (dots) for polynomials. To generalize this procedure to two-dimensional enclosures, we must be aware of a difference betweenDirichlet boundaries and boundaries with finite t . If one has a Dirichlet BC, one can square theDirac-Weyl Hamiltonian and easily prove that each spinor component obeys a Schrödinger equationinside the enclosure where the mass term is zero. As such, one can use the method we presentedbefore for this equation to calculate the square of the eigenvalues and the spinor component thatis zero at the boundary. Since the Dirac-Weyl equation has particle-hole symmetry—if one changesthe sign of one of the spinor components in a state with energy E , one obtains a state of energy − E ,which still preserves the Dirichlet BC—, one effectively obtains all the eigenvalues. Applying theDirac operator to the component which is zero at the Boundary, one generates the other componentof the spinor and completes the solution [14]. The same procedure cannot be used for finite t , when the BC only fixes the ratio between the spinorcomponents, ψ B /ψ A = te iφ . The circle is the only 2D enclosure where an exact solution is known10or the infinite-mass BC we are using [22]. The solution is given by [23] Ψ ǫ,m ( r, θ ) = (cid:20) e imθ J m ( q ǫ r ) ise i ( m +1) θ J m +1 ( q ǫ r ) (cid:21) (26)where s = ± is the band index, J m is the Bessel function of the first kind of order m and ǫ = s ~ v F q ǫ is the energy of the state. Imposing the boundary condition, the allowed energy levels are those inwhich J m +1 ( q ǫ R ) J m ( q ǫ R ) = s. (27)To apply the polynomial method to this enclosure, we must build the initial functions carefully.The rotation symmetry is continuous: ˆ U φ Ψ = e − i ( φ/ σ z R φ Ψ = (cid:20) e − iφ/ ψ A ( r, θ − φ ) e iφ/ ψ B ( r, θ − φ ) (cid:21) . Provided (cid:20) ψ A ( r, θ − φ ) ψ B ( r, θ − φ ) (cid:21) = (cid:20) e − imθ ψ A ( r, θ ) e − i ( m +1) θ ψ B ( r, θ ) (cid:21) we get ˆ U φ Ψ = e − ikφ Ψ , (28)with k = m + 1 / . The only functions which are invariant under this symmetry are linear combin-ations of | z | n = r n , for n ≥ . For each m , the initial states that respect the symmetries and theBC, ψ B /ψ A = ie iθ are (with z = x + iy and ¯ z the complex conjugate of z ) Ψ ( m ≥ ( z, ¯ z ) = z m (cid:20) − iz (cid:21) , Ψ ( m< ( z, ¯ z ) = (¯ z ) − ( m +1) (cid:20) ¯ zi (cid:21) . (29)With no loss of generality, we discuss the states with m = 0 . The procedure is analogous for theremaining representations. A polynomial state can have the form | P ( r ) , Q ( r ) i = (cid:20) − iP ( r ) zQ ( r ) (cid:21) (30)if P and Q are even polynomials in r such that P (1) = Q (1) = 1 , thus respecting both the BCand the symmetries in question. Simple monomials P n ( r ) = r n satisfy both conditions and anon-orthogonal basis is {| P , P i , | P , Q i , | P , P i , | P , P i , | P , P i , . . . } (31)11fter applying the G-S process, the Hamiltonian is straightforwardly computed and diagonalized fordifferent basis sizes.The comparison with the exact results is presented in Table 1. As the method quickly convergestowards the exact results, we will now apply it to the square enclosure with the same BC.Basis Size ǫ ǫ ǫ ǫ . − . . − . . . − . . − . . − . . − . Exact . − . . − . Table 1: Convergence of polynomial method for the circular enclosure.
The case of square with edges at x, y = ± L/ and t = 1 has no known analytic solution. As theseboundary conditions correspond to an uniform gap outside the enclosure, the Hamiltonian will beinvariant under π/ rotations. ˆ U = e − i π σ z R π/ The one-dimensional representations of this symmetry are given by ˆ U Ψ = e iα Ψ , where α = π n π ∈ [ − π, π ] . (32)The other possible symmetries of the Hamiltonian would be reflections. These, however, change thesign of the mass gap outside the enclosure and are also not compatible with the imposed boundaryconditions.The following are low order polynomial functions, belonging to each of the representations of thesymmetry group. Φ π/ = (cid:20) ¯ z (cid:21) , Φ − π/ = (cid:20) zz (cid:21) , Φ π/ = (cid:20) ¯ z ¯ z (cid:21) , Φ − π/ = (cid:20) z (cid:21) . (33)Each time we multiply one of these states by z m ¯ z n , α changes as α → α + ( n − m ) π . (34)12e stay in the same irreducible representation, provided α is unchanged modulo π , i.e. , if n − m = 0 mod 4 . (35)Invariant monomials of z and ¯ z are | z | , z , (¯ z ) and their products. There are two issues that needto be addressed at this point: (a) the states of Eq.(33) do not satisfy the BC, and are therefore not suitable as starting statesupon which to build a basis. Taking the α = − π/ as an example, we require states of theform Φ = P (cid:16) | z | , z , (¯ z ) (cid:17) zQ (cid:16) | z | , z , (¯ z ) (cid:17) (36)where the invariant polynomials, satisfy P − Q z = 1 for y = − / . If this is verified at thisedge, the symmetry will ensure the BC are also verified at the other edges. The choice P = − i − iz + i z , Q = 1 + z − | z | −
15 ¯ z (37)proved convenient in the sense that the behaviour of the ket Φ ( x, y ) near the origin resemblesthat of the lowest eigenvalue of the circular enclosure, (b) the basis must be built by multiplying each component of the spinor by increasing order invariantpolynomials which are required to be unity at the edges of the enclosure so that the BC arepreserved. Such polynomials can be generated order by order, by imposing P n (cid:16) | z | , z , (¯ z ) (cid:17) =1 for y = − / .Once these steps are performed, the calculation proceeds as for the circle, by orthogonalizing thebasis and calculating of the Hamiltonian matrix, separately for each of the 4 symmetry representa-tions. The necessary work is halved because the transformation σ x K ( K is the complex conjugationoperator) changes the sign of the Hamiltonian and preserves the BC; the states of the representationswith α = 1 / , π/ can be obtained from those of α = − / , − π/ by this transformation; theeigenvalue spectrum remains symmetrical.The spectrum is displayed in Fig.6, with the low-energy regime in the inset. The number ofpolynomials in the legend refers to the basis size in each irreducible representation.13
30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30
Eigenvalue Index -50-40-30-20-1001020304050 E n e r g y ( h_ v F / L ) Figure 6: Comparison of the spectrum of the uniform-BCs square billiard with increasing number ofpolynomials.This method again converges, as it did in case of the circular enclosure.As stated before, the method also provides the eigenstates as polynomials. The valence band statewith energy closest to zero is plotted in Fig.7. At the center of the square, the behaviour of the twospinor components is not unlike that of the corresponding state of the circle, but is quite differentat the edges. The two components have the same absolute value at the edges, as required by theBC, and both vanish at the corners, as that is the only way the BC of the edges that meet at thecorner can be satisfied.To finalize the analysis of the method, we tested it in a problem with BCs analogous to thosefound in zigzag-terminated hexagonal graphene flakes.Figure 7: Absolute values of the spinor wave-functions for the eigenstate of energy ǫ = − . ~ v F /L .Left panel, | ψ A | ; middle, | ψ B | ; right, both components, showing that | ψ A | = | ψ B | at theedges as required by the BC. Near the center of the square the behaviour resembles thatof the corresponding state for the circle, but here the spinor wavefunction is zero at thevertexes, as this is the only way to satisfy the BC for the two edges that meet at the vertex(colour online). 14 .5 Dirichlet Hexagonal Enclosure We applied this method to the study of an hexagonal enclosure with the BCs equivalent to those ofan hexagonal graphene flake whose edges are all zigzag-terminated. This specific set of BCs can beexpressed by alternating t between t = 0 and t → ∞ , depending on the termination sublattice.As we are working with Dirichlet boundaries in alternating sides for each spinor component, wecan construct the initial polynomial as in equation 5. Considering a regular hexagon centered at theorigin of side-length L , the two spinor components for the starting state will be defined as ψ ,A ( x, y ) = " √ L y L + (cid:18) x − y √ (cid:19)(cid:21) × (cid:20) L − (cid:18) x + y √ (cid:19)(cid:21) (38) ψ ,B ( x, y ) = ψ A ( x, − y ) . With this, we diagonalize the Hamiltonian h Ψ i | H † H | Ψ j i = h ψ i,A | (cid:0) − ~ v F ∇ (cid:1) | ψ j,A i ++ h ψ i,B | (cid:0) − ~ v F ∇ (cid:1) | ψ j,B i . (39) Eigenvalue Index E n e r g y ( h _ v F / L ) Figure 8: Comparison of the spectrum of the zigzag-like hexagonal enclosure with increasing numberof polynomials.The resulting spectrum is visible in Fig.8 for different basis sizes. The numerical solution was15btained via
Wolfram Mathematica ® by imposing Dirichlet BCs on three alternating sides of theenclosure in question, leaving the remaining free.Finally, we will compare these results with the ones obtained by Gaddah [14] for the problemof the triangular billiards with boundary condition ψ A = 0 (in our notation, t → ∞ ). The exactspectrum described by the author follows the relation ( L is the side-length of the equilateral trianglein question) E n ,n = ± πL ~ v F q n + n n + n , (40)where n ≥ n > . As described in the article, the states with n > n are (at least) × degenerate.The BCs we imposed on a hexagon, ψ A = 0 on three alternating sides (Eq.(38)), also enforce ψ A = 0 on the edges of an equilateral triangle, that contains the hexagon region in question. The conditionfor ψ B is the same, as it is a simple inversion of this same triangle. As a consequence, after re-scalingthe exact spectrum for the triangle by p A hex /A triang , we obtain a match between the two as shownin Fig.10. We also compare the density plots for the first three approximate eigenfunctions to thosepresent in Gaddah’s work. This comparison (using 10 polynomials) is presented in Fig.9 for | ψ A | ,with the corresponding regions highlighted.Figure 9: | ψ A,j | for the hexagon (top) and for the triangle (bottom) j ∈ { , , } (color online).16 Eigenvalue Index E n e r g y ( h _ v F / L ) Polynomial Method (hexagon, 66 polynomials)Exact Solution (triangle, rescaled)
Figure 10: Comparison of the spectrum of the zigzag-like hexagonal billiard against the exact solutionfor the triangle.
With the polynomial method we have been able to construct approximate solutions to Schrödingerand Dirac-Weyl equations for planar convex polygonal enclosures. The method has been appliedbefore to the determination of frequencies of plate vibrations in the Mechanical Engineering com-munity, but as far as we know its generalization to the Dirac billiards has not been reported before.In the this case we considered different types of boundary conditions, including situations whereonly the ratio of spinors components is specified.The method is based on constructing a truncated basis of states, all of which satisfy the prescribedBC from the start. By comparing with cases where the exact solution is known (most often byseparation of variables) we were able to ascertain the convergence of the method. Typically, a basissize of order n is required to obtain n eigenvalues with close to 1% accuracy. The method providesnot only the lowest energies (or lowest absolute energies for the Dirac -Weyl case) but also eigenstatesas finite order polynomials in the coordinates. Acknowledgments
The authors acknowledge financing of Fundação da Ciência e Tecnologia, of COMPETE 2020 pro-gram in FEDER component (European Union), through projects POCI-01-0145-FEDER-028887and UID/FIS/04650/2019. The authors also acknowledge financial support from Fundação para aCiência e Tecnologia, Portugal, through national funds, co-financed by COMPETE-FEDER (grant17-ERA-NET2/0002/2016 – UltraGraf) under the Partnership Agreement PT2020.
References [1] G. Lamé. Leçons sur la Théorie Mathématique de l’Elasticité.
Bachelier, 1852 .[2] M.A. Pinsky. The Eigenvalues of An Equilateral Triangle.
Siam Journal On MathematicalAnalysis , 11(5):819–827, 1980.[3] BJ McCartin. Eigenstructure of the equilateral triangle, Part I: The Dirichlet problem.
SiamReview , 45(2):267–287, JUN 2003.[4] V Amar, M Pauri, and A Scotti. Schrödinger-equation For Convex Plane Polygons - A TilingMethod For The Derivation Of Eigenvalues And Eigenfunctions.
Journal of MathematicalPhysics , 32(9):2442–2449, Sep 1991.[5] Wai-Kee Li and S. M. Blinder. Particle in an equilateral triangle: Exact solution of a nonsep-arable problem.
Journal of Chemical Education , 64(2):130, February 1987.[6] W. A. Gaddah. A Lie group approach to the Schrödinger equation for a particle in an equilateraltriangular infinite well.
European Journal of Physics , 34(5):1175–1186, July 2013.[7] V Amar, M Pauri, and A Scotti. Schrödinger-equation For Convex Plane Polygons .2. A No-go Theorem For Plane-waves Representation Of Solutions.
Journal Of Mathematical Physics ,34(8):3343–3350, Aug 1993.[8] Berry Michael Victor and Mondragon R. J. Neutrino billiards: time-reversal symmetry-breakingwithout magnetic fields.
Proceedings of the Royal Society of London. A. Mathematical andPhysical Sciences , 412(1842):53–74, July 1987.[9] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva,S. V. Dubonos, and A. A. Firsov. Two-dimensional gas of massless Dirac fermions in graphene.
Nature , 438(7065):197–200, November 2005.[10] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim. The electronicproperties of graphene.
Reviews of Modern Physics , 81(1):109–162, January 2009.[11] L. A. Ponomarenko, F. Schedin, M. I. Katsnelson, R. Yang, E. W. Hill, K. S. Novoselov, andA. K. Geim. Chaotic Dirac Billiard in Graphene Quantum Dots.
Science , 320(5874):356–358,April 2008.[12] Florian Libisch, Christoph Stampfer, and Joachim Burgdörfer. Graphene quantum dots: Bey-ond a Dirac billiard.
Physical Review B , 79(11):115423, March 2009.1813] M. Zarenia, A. Chaves, G. A. Farias, and F. M. Peeters. Energy levels of triangular andhexagonal graphene quantum dots: a comparative study between the tight-binding and theDirac approach.
Physical Review B , 84(24):245403, December 2011.[14] W. A. Gaddah. Exact solutions to the Dirac equation for equilateral triangular billiard systems.
Journal of Physics A: Mathematical and Theoretical , 51(38):385304, September 2018.[15] Liang Huang, Ying-Cheng Lai, and Celso Grebogi. Characteristics of level-spacing statistics inchaotic graphene billiards.
CHAOS , 21(1), MAR 2011.[16] Liang Huang and Ying-Cheng Lai. Perspectives on relativistic quantum chaos.
Communicationsin Theoretical Physics , 72(4), APR 1 2020.[17] Y Shimizu and A Shudo. Polygonal Billiards - Correspondence Between Classical TrajectoriesAnd Quantum Eigenstates.
Chaos Solitons & Fractals , 5(7):1337–1362, Jul 1995.[18] R.B. Bhat. Flexural vibration of polygonal plates using characteristic orthogonal polynomialsin two variables.
Journal of Sound and Vibration , 114(1):65–71, January 1987.[19] K. M. Liew, K. Y. Lam, and S. T. Chow. Free vibration analysis of rectangular plates usingorthogonal plate function.
Computers & Structures , 34(1):79–85, January 1990.[20] K. M. Liew and K. Y. Lam. A Set of Orthogonal Plate Functions for Flexural Vibration ofRegular Polygonal Plates.
Journal of Vibration and Acoustics , 113(2):182–186, April 1991.[21] H. Larcher. Notes on orthogonal polynomials in two variables.
Proceedings of the AmericanMathematical Society , 10(3):417–423, 1959.[22] M. Quintela. From the 1D Schrödinger infinite well to Dirac-Weyl graphene flakes. Master’sthesis, Faculdade de Ciências, Universidade do Porto, 2019.[23] B. Wunsch, T. Stauber, and F. Guinea. Electron-electron interactions and charging effects ingraphene quantum dots.
Phys. Rev. B , 77:035316, Jan 2008.19, 77:035316, Jan 2008.19