Electronic properties of curved few-layers graphene: a geometrical approach

We show the presence of non-relativistic L\'evy-Leblond fermions in flat three- and four-layers graphene with AB stacking, extending the results obtained in [Curvatronics2017] for bilayer graphene. When the layer is curved we obtain a set of equations for Galilean fermions that are a variation of those of L\'evy-Leblond with a well defined combination of pseudospin, and that admit L\'evy-Leblond spinors as solutions in an approriate limit. The local energy of such Galilean fermions is sensitive to the intrinsic curvature of the surface. We discuss the relationship between two-dimensional pseudospin, labelling layer degrees of freedom, and the different energy bands. For L\'evy-Leblond fermions an interpretation is given in terms of massless fermions in an effective 4D spacetime, and in this case the pseudospin is related to four dimensional chirality. A non-zero energy band gap between conduction and valence electronic bands is obtained for surfaces with positive curvature.


INTRODUCTION
Graphene is a two-dimensional (2D) material which enables the realization of different heterostructures and electronic devices with novel quantum properties [1][2][3]. The electronic and transport properties of graphene in a multilayer stacking can be modified and tuned by changing the number of layers of graphene in the system. Remarkably, the low energy electronic band structure of multilayered graphene evolves from Dirac cones in monolayer graphene with massless Weyl excitations, to parabolic bands with massive excitations in the case of bilayer graphene, to more complex band structures and peaked density of states when a few-layer graphene system is realized with different stacking orders [4][5][6][7]. Few-layer graphene can be fabricated from graphite by mechanical exfoliation [8,9] and by chemical techniques [10][11][12] controlling the stacking order. Experimental characterizations of electronic and transport properties of trilayer graphene have been reported in Ref. [13][14][15]. An electric field perpendicular to the graphene layers has been shown experimentally to open a band gap in the single-particle electronic energy spectrum in bilayer [16] and trilayer [17,18] graphene, being this effect of crucial importance in realizing semiconducting behavior for electronics. Theoretical predictions of a similar band gap opening in the energy spectrum of few-layer graphene have been reported [19]. The electronic band structure of graphene and many other 2D systems can be also tailored by strain of their lattice in different directions, inducing several possible quantum effects. Strain plays a key role in iron-based superconductors [20], in cuprate superconductors [21], and in superconducting diborides [22]. In cuprate superconductors the displacement of the Cu ions along the transversal direction induces strain modulations and bending in the copper-oxide planes in the form of lattice stripes leading to amplification effects of the superconducting parameters [23,24]. Curved 2D systems are a source of different kind of compressive and tensile strains, with a wide range of strain percentages of the lattice depending on the topology and on the curvature radius of the surface of the curved systems. Therefore we expect that curvature will be another possible way to engineer the electronic band structure of 2D systems, as we will discuss in this work. Recently the low energy electronic properties of bilayer graphene have been studied in the context of a geometrical approach, using the Lvy-Leblond equation [25] for massive particles [26] in a curved space. Although the study of graphene through relativistic approaches is a quite new field of research, a significant amount of literature already exists, in the study of deformed monolayer graphene employing the relativistic 2D Dirac equation in a curved metric [27][28][29][30][31][32][33]. This connection between graphene and relativity has been explored also in other contexts, for instance for the evolution of the free electron current density in graphene with defects [34,35]. On the other hand, the scientific challenge to connect Dirac theory of the spinning electron with gravitation dates back to the seminal works of Hermann Weyl [36]. Interestingly, Weyl fermions with zero mass and definitive chirality have not been identified in high-energy physics, but they have been clearly observed in novel topological materials in condensed matter systems [37]. The role of two-dimensionality in generating Weyl states in systems with parabolic electronic bands has been demonstrated in Ref. [38], while the influence of zero helicity states at the interfaces of layered heterostructures, leading to local magnetization and mass anisotropy, has been discussed in Ref. [39]. The above described example demonstrates that solid state systems with geometrical constraints and specific topological properties can work as model systems for physical phenomena with very different energy scales, as high-energy physics and gravitation in a curved D-dimensional space-time. The mapping proposed in Ref. [26] allowed to study the evolution of the band structure, both conduction and valence bands and their band-gap, as a function of geometrical curvature. Positive (spherical-like) curvatures have a local effect of opening a band gap in the spectrum, while negative (hyperbolic-like) curvatures tend to close the band gap. The band-gap energy is predicted to be tunable and proportional to the curvature radius of the deformation. This paves the way to "Curvatronics" with bilayer graphene, an interesting possibility for applications to control in a static way the electronic properties of layered systems with the curvature. In this work we extend the analysis of Ref. [26], analysing massive and massless fermionic excitations in the more complicated case of few-layer flat graphene systems, while for the case of the curved bilayer we refine the treatment, proposing a set of equations that take into account the exact combination of pseudospin states that defines physical states, i.e. eigenstates of the energy. Similar considerations apply to multi-layer systems and will be studied separately.
The work is organised as follows. In sec. we discuss Galilean fermions in few layers graphene: we begin in recalling the results obtained in [26] for bilayer graphene, and then in and we show that Galilean fermions are present in the case of three and four layers as well. We obtain the exact spinor solutions for all energies, and study them in the proximity of the Galilei points, where the electronic dispersions acquire their extreme values: the results is given by 2D spinors that satisfy the Galilei invariant Lévy-Leblond equation, and we show where these fermions are localised. Then in sec. we discuss various aspects of the Lévy-Leblond equation, both in the flat and curved case. We recall its Galilean invariance and discuss its two independent solutions, which correspond to states with different pseudo-spin. It has been discussed in [26] how solutions of the Lévy-Leblond equation in 2D can be lifted to solutions of the massless Dirac equation in 4D. We show here how states with a definite pseudospin lift to states with definite chirality, with a perfect match of degrees of freedom: the states corresponding to positive and negative isospin are mapped into the two solutions of the massless 4D Dirac equation, one with left and one with right chirality. Then we discuss in detail the case of the curved bilayer. We present a set of covariant equations that describe the generalization of the flat massive fermions, show how the equations predict that energy eigenstates are given by a well defined mixture of pseudospin states and infer that the scalar curvature of the surface alters the local energy density.

Bilayer graphene
In this section we review the results obtained in [26] for bilayer graphene with AB Bernal stacking, presenting them in a manner that is suitable to generalisation to a higher number of layers.
Bilayer graphene with AB stacking presents one set of metallic bands E (±) 1 , and another of isolating bands E (±) 2 . We define a quantity with the dimensions of mass where γ is the interlayer hopping parameter, γ ∼ 0.4eV, and v F ∼ 10 6 ms −1 is the Fermi velocity in graphene. The energy bands can be described in terms of m 0 and of the effective mass m = ± m0 2 for bilayer graphene, where the ± factor denotes either positive or negative bands: It is useful to define a dimensionless factor α by such that, in the bilayer case, α = 2. The k · p model that describes these bands is defined in terms of a spinor λ with 4 components. The Hamiltonian is given by where κ = τ k x + ik y is the wave number of the excitation, with τ = ±1 denoting the two inequivalent Fermi points 1 . These are the points in the Brillouin zone where the conduction and valence bands touch, and for graphene they are at the level of the Fermi energy. The k · p model can be obtained in three steps: 1) assuming that both the action of the full Hamiltonian, and the overlap between wavefunctions, is restricted to nearest neighbour contributions (the tightbinding model restricted to the nearest-neighbor hoppings gives a satisfactory description of the electronic structure of bilayer graphene at low energies, see e.g. [40]), 2) expanding the momentum close to the Fermi points, and 3) assuming that the most relevant modes for the expansion of the full wavefunction Ψ of a single electron are given by p z orbitals, where z is the direction transverse to the bilayer, as these are the orbitals that contribute the most to conduction. In formulas where A, B correspond to inequivalent types of carbon atoms in the same plane, and the index 1, 2 corresponds to different layers. Given the form (4) of the Hamiltonian, we can recognise that B1 and A2 are the two atoms that are directly aligned along the transversal direction. Solving the eigenvalue equation Hλ = Eλ we find that there are no non-trivial eigenspinors with λ 1 = 0 when |κ| = 0. The component λ 1 can then be set to 1 as the eigenvalue equation is homogeneous. The full solution follows from linear algebra where σ = ±1 is a factor labelling different bands that appear in the consistency condition satisfied by the energy In [26] we expanded the solution around the Fermi points in the dimensionless parameter ǫ = vF |κ| γ . For the metallic E (±) 1 bands one obtains For simplicity we set τ = +1 and focus on one of the Fermi points. To begin, it is useful noticing that for all values of the constant A the following spinors satisfy the Lévy-Leblond equations [25] i where D = iσ j k j is the 2-dimensional Dirac operator in phase space and the σ j are the Pauli matrices in the standard basis. As we will discuss in sec. , these are coupled, first order equations that are invariant under the Galilei group of non-relativistic mechanics and describe non-relativistic spin 1 2 fermions. In particular, the spinors with λ 3 = 0 = λ 4 correspond to a solution with definite pseudospin, and those with λ 1 = 0 = λ 2 to an independent solution with opposite pseudospin. We notice from the explicit form of the solution (8-10) that, if in the solution we exchange at the same time the energy E + 1 with E − 1 , keeping κ unchanged, then the components λ 1 and λ 3 of the solutions are unchanged, while λ 2 and λ 4 change sign. Thus, by taking a linear combination of such eigenstates of energies E + 1 and E − 1 it is possible to obtain 'checkered' states where either λ 2 = 0 = λ 4 , i.e. only sites of type A are excited, or viceversa states with λ 1 = 0 = λ 3 , i.e. only sites of type B are excited. Conversely, taking a linear combination of the latter it is possible to build eigenstates of the energy. We also notice that, as explicit electronic excitations are defined by their Fourier coefficients for each value of the momentum κ, these are free parameters in our analysis, which concerns only the energy spectrum. This is why we are able to set λ 1 = 1 at |κ| = 0. In sec. we will see how these results are generalized in the case of curved sheets. We will obtain equations for Galilean fermions that are a variation of the Lévy-Leblond equation, which involve a well defined combination of pseudospin. We will see that, at least in the axially symmetric case, it is still true that linear combinations of states with opposite energies give checkered states.
In analogy with the fact that for monolayer graphene the Fermi points are called Dirac points because of the relativistic emergent symmetry, for the bilayer we will use the wording Galilei points. At the Galilei point, i.e. for κ = 0, the solutions have λ 2 = 0, λ 3 = 0. As shown in Figure 1a, this corresponds physically to activating only the p z orbitals of the A 1 and B 2 atoms: these are not directly aligned along the z direction, thus giving a small overlap of the electron orbitals.
Conversely, when we perform a similar analysis for the insulating bands E find and This time, at the Galilei point λ 1 = 0, λ 4 = 0 and only the orbitals of the B 1 , A 1 atoms are activated, as shown in are out of phase, thus minimising it.

Three-layer graphene
Here we extend our approach to the three-layer graphene with the ABA stacking of the layers, which is represented in Figure 2a.
The effective low energy Hamiltonian for ABA stacked graphene around the K point can be found in [41]: Here we see that the atoms that are directly aligned along the transverse direction are B 1 and A 2 , as before, and A 2 and B 3 . Analysing the band structure one finds two touching bands linear in momentum and a set of metallic and insulating bands given by We begin by analising the metallic and insulating bands: we perform an analysis near the Fermi points in order to describe the detailed structure of the spinors. This will allow us to build spinors that satisfy the Lévy-Leblond equations. We expand the energy bands in the dimensionless parameter ǫ = vF |κ| γ , finding for the two metallic bands with effective mass m = ± γ √ 2v 2 F , or equivalently, recalling (3), α = √ 2, a smaller value than for AB bilayer graphene. For the two insulating bands This is the same gap, for given mass, found in the bilayer case. The energy is again in non-relativistic form and therefore we expect that non-relativistic fermions will appear. The full formula (22) for the energy levels can be rewritten in terms of m as It is a different formula from (2) but it has an identical low energy limit. The spinor solution of the equation Hλ = Eλ is given by where σ = ±1 is a factor labelling different bands that appears in the consistency condition satisfied by the energy Now we analyse the spinor solutions close to the Fermi point for the metallic bands: we substitute the low energy limit (23) into the spinor solution, and retain the lowest order in the parameter ǫ. The result is given by (8)-(10) plus If we build χ 1 and χ 2 as in (11), (12) then again these satisfy the Lévy-Leblond equations. There is also another set of Lévy-Leblond equations with spinors where B is another arbitrary constant. As for the bilayer, at the Fermi point κ = 0 the p z orbitals of the atoms directly aligned B1, A 2 , B 3 are turned off. So we see that the non-relativistic fermions for the metallic bands of three-layer graphene consist of states where interlayer hopping happens between all three layers, generalizing the physics occurring in bilayer graphene: interlayer hopping is a fundamental ingredient in order to generate a mass for the excitations. The same procedure, applied to the insulating bands yields the expansion (15)-(17), plus the low energy limit This time at the Fermi point the contributions of the atoms A 1 , B 2 , A 3 are turned off. The Lévy-Leblond spinors in this case are given by We conclude this section analysing the bands that are linear in momentum. The solution close to the bottom of such Dirac bands is given by Then the spinors satisfy the Weyl equations for massless spinors These solutions are localised on the layers 1 and 3, and they vanish in the intermediate layer 2. On the layers, they describe massless modes moving at the speed of sound v F . One can think of these solutions in this way: the electronic wavefunctions undergo totally destructive inteference in the middle layer, and the remaining electrons are confined to the top and bottom layers. Since interlayer hopping is forbidden in these localized states, then the excitations become massless again as for monolayer graphene.

Four-layer graphene
The case of four-layer graphene with the ABAB stacking of the layers is finally considered. Its schematic representation is reported in figure 2b. The ABAB Hamiltonian H (ABAB) k around the K point is given by This is the same as for the three-layer case, with the addition of atoms B 3 , A 4 directly aligned. We plain to repeat the analysis done in the previous section: first finding the shape of the bands, then looking for special types of excitations. In this case we will find only Lévy-Leblond spinors, and no massless excitations. Analysing the eigenvalues of H (ABAB) k we find 8 bands, expressed in terms of two masses: In our notation α 1 = 4 And for m 2 The components of the energy eigenspinors are given by where σ = ±1 is a factor labelling different bands that appears in the consistency condition satisfied by the energy At the bottom of the metallic bands E Similarly, at the bottom of the insulating bands E Therefore we see repeating the same pattern that we found for the bilayer and three-layer cases: at the Galilei point for the metallic bands the component related to the A 4 atom is turned off, together with those of all the atoms that are directly aligned, and at the Galilei point for the insulating bands conversely the B 4 component is zero.
The new Lévy-Leblond spinors for the metallic bands are while for the insulating bands where C is another arbitrary constant. We see here that with four layers, no configurations with totally destructive interference in the middle layers appear. All states present interlayer hopping and are massive.

THE CURVED CASE
The Lévy-Leblond equation The Lévy-Leblond equations (13), (14) are coupled, first order equations. Lévy-Leblond in his seminal paper [25] showed that these are obtained from the theory of representations of the Galilei group applied to particles of spin 1 2 . In general, the Lévy-Leblond equations can be solved in the following way: solving for χ 1 in (14) and plugging it into (13) we obtain the second order, uncoupled equation which in fact is the Pauli equation in the special case of the external magnetic field set to zero, see [26] for the case of a non-zero field. This equation admits the independent solutions which correspond to different values of 2D (pseudo)spin. Then using (14) χ 1 is obtained by differentiation of χ 2 as By linearity one can then construct the full solution: Let us compare this general formula with, for example, the Lévy-Leblond spinor obtained for the metallic bands E (±) 1 in the case of bilayer graphene: corresponds to states of definite value of the pseudospin. Similar formulae hold for Lévy-Leblond spinors related to the insulating bands, as well as for the three-and four-layer case. States with definite pseudospin in the case of bilayer graphene close to the Galilei points are schematically shown in Figure 3.
We now recall the relationship with the massless Dirac equation in 4 dimensions. It has been shown in [26] that given a solution (χ 1 , χ 2 ) of the Lévy-Leblond equations it is possible to construct a massless Dirac spinor in an effective 4-dimensional Minkowski spacetime where u, v are conjugate null coordinates. The massless 4D spinor is given bŷ In terms of degrees of freedom, the massless Dirac equation in 4D has two degrees of freedom, corresponding to two different allowed chiralities. In turn, the Lévy-Leblond equations admit two independent solutions, as just seen. This fact is well known in the literature [42][43][44][45], and is related to the concept of the Eisenhart-Duval lift of non-relativistic mechanics [46][47][48]. The relationship can be clearly seen in terms of the Gamma matrices adapted to the 4D lift: In terms of these the chirality matrix is calculated as showing that the spinorΨ built from the combination a = 1, b = 0 has a definite chirality, which by convention we can call right, and the one built using a = 0, b = 1 has opposite chirality, say left. That chirality in 4D can be related to pseudospin in 2D can be understood in the following way. First, we notice that in 4D we can introduce a timelike variable t and a spacelike variable z according to 2dudv = −dt 2 + dz 2 . Then, from eq.(67) chirality is decomposed into a boost in the z direction, times a rotation around the z axis. Spinors of the form (66) with either χ 1 = 0 or χ 2 = 0 are eigenspinors of boosts along the z direction: this is related to the fact that the metric (65) is written in null coordinates u, v and that the spinors are adapted to the coordinates. So for these spinors 4D chirality and 2D isospin are directly related.

The role of curvature
We now consider curved few layers of graphene. Our chief example will be the bilayer case, however similar considerations apply to the case of three or four layers. When the radius of curvature is small enough we consider a Hamiltonian that satisfies the following two requirements. First, the Hamiltonian should be self-adjoint, so that the energy eigenvalues are real. Second, it should be a covariant generalization of (4), which reduces to (4) when the metric is flat. In order to do so, we replace the momentum k with −i∇, the 2D spinorial covariant derivative that includes the spin connection of the curved metric g ij dq i dq j , with q i = {x, y}. In sec. we have seen how the Lévy-Leblond spinors appeared reorganizing the spinor components as per (11). Motivated by this insight we make the following change of basis The Hamiltonian changes accordingly into or which is (formally) self-adjoint since D † = −D. The electronic wavefunction is and similarly to what we have done in the flat case we regroup the components into The eigenvalue equation for H ′ κ becomes We can decouple the equations by obtaining χ 2 from the second, and inserting it into the first. An important case is that of metallic bands, for which the Eχ 1 term in the top equation is negligible with respect to the other two: in this limit Recalling that γ = 2|m|v 2 F , these equations are similar to the stationary Lévy-Leblond equations, with the difference that the term −2mv 2 F is replaced by the matrix 0 γ γ 0 . To understand the relationship with what we found in the flat case, we isolate χ 2 in (81), plug it into (80) and obtain This is the eigenvalue equation for the energy in the curved case. In flat space and for metallic bands we found so that (82) became and in fact our solution was Notice that in this case χ 1 is an eigenspinor of the matrix σ 1 with eigenvalues ∓1, so that we are justified in exchanging the matrix 0 γ γ 0 with −2mv 2 F . However, in general the operator on the right hand side of (82) does not commute with σ 1 , and as a consequence eqs.(80-81) are a variation of the Lévy-Leblond equations.
Our theory predicts that there will appear a gap in the spectrum if there are no solutions of (82) for the eigenvalue E = 0. In this particular case the analysis reduces to finding zero modes of D , and it is a known result that if the surface has positive curvature then the spectrum of D does not include the eigenvalue zero.
Eq.(82) is in general non-trivial to study, and to gain a better understanding we focus our attention on the axisymmetric case. That is, we consider a metric of the type For example for C(r) 2 = sin(r) 2 the geometry is that of a sphere, with constant positive curvature, and for C(r) 2 = r 2 + a 2 , with a a constant, the geometry describes two flat ends joined by a tube, and has negative, non-constant, curvature. In the generic axisymmetric case the Dirac operator is given by where so that We will also assume that the operator acts on a spinor of the form χ 1 = e ilϕ a(r) b(r) , l = ± 1 2 , ± 3 2 , . . . . From our calculation then eq.(82) reduces to the pair of equations where we defined the second order, self-adjoint, positive operators A = OO † , B = O † O. These equations can be decoupled, giving or in other words on the space of χ 1 spinors the operator γ 2 E 2 is represented by To summarize, we obtained a fourth order equation for the square of the energy eigenvalues. Our equation is consistent since it is a known fact that for positive operators A and B the spectrum of AB is the same of that of BA, and it is positive [49]. Being fourth order, our equations will be challenging to solve analytically in general: the solution will require a careful but feasible numerical approach, which we postpone to future works. There are however some properties that are easy to discuss. First, since they are quadratic in E, then the spectrum will be symmetric under E → −E. Second, as the operator in (93) commutes with the matrix σ 3 , then we can see that if a, b are eigenfunctions related to the eigenvalue E, then a, −b are eigenfunctions related to the eigenvalue −E. Thus, in analogy to the flat case, we are able to construct checkered solutions from linear combinations of solutions with opposite energy, and viceversa. However, in general it will be the case that a = b, since A = B. Therefore the solutions cannot be related to Lévy-Leblond fermions. There is however a regime where we reasonably expect that A ∼ B. One can calculate Suppose the function 2 C ′ C 2 has a maximum µ. We can substitute B in (89) for A + O † , O and, in those cases where we can ignore the extra term and reduce the eigenvalue equations to These admit two types of solutions: solutions with b = a give i.e. a is an eigenfunction of A and the energy is negative, or b = −a which gives i.e. a is an eigenfunction of A and the energy is positive. In this approximation of low angular momentum to energy ratio the problem simplifies considerably, as one needs to solve 2nd order differential equations. Summarising our results, we expect that deforming the few-layer graphene sheets provides an efficient way to increase or reduce an energy band gap between conduction and valence bands, since few layers graphene has several degrees of freedom with respect to which deformations can be introduced: both along the sheets and perpendicular to it. As an illustrative example, possible positive or negative curvatures of the three-layer graphene are depicted in Figure 4.

CONCLUSIONS
In this work we have extended the geometrical approach and the results of our previous work on curvatronics with bilayer graphene of Ref. [26] to the case of three and four layers of AB stacked graphene, and we have refined the role of isospin states in the case of the curved bilayer. We found that for flat layers Lévy-Leblond fermions are still present at the bottom of the parabolic bands, which we call Galilei points of the Brillouin zone, and we gave a geometrical interpretation of the solutions in terms of the p z orbitals of graphene associated to them, on the basis of their non-trivial spatial overlap. We have discussed in detail the effective electronic masses present for three and four layers, comparing them with the bilayer graphene case, as well as the explicit form of the solutions. We analysed the Lévy-Leblond equation and showed how energy eigenstates are given by the superpositions of checkered states. Using the relationship between the 2D Lévy-Leblond equation and the 4D massless Dirac equation, we were able to relate 2D pseudospin with chirality in 4D. Lastly, we discussed the fact that curving few layers graphene provides a tool to tune the energy band gap, which has implications for tunable electronics based on curvature, or curvatronics with few-layer graphene systems. For instance, a graphene-based device with alternating regions of positive or negative curvature, inducing insulating or conductive states, could be exploited to produce non-volatile memories. It would be interesting for future research to study in detail the consequences of having curved few layers, as well as generalize our results to a higher number of layers, and to study the fermionic excitations away from the Galilei points. Another avenue of possible future research is studying the physics of pseudospin in settings with non-trivial curvature or topology, which might lead to novel effects.
Finally, very recently the class of Weyl metamaterials has been introduced in Ref. [50]. In these systems, chiral Weyl fermions are moving in an artificial 3D curved-space geometry, with applications to tunable novel electronic devices through curvature engineering and with connections to cosmological theories. Therefore, the new field of quantum physics in curved spaces is likely to become of great and practical interest in the near future.