Local sublattice symmetry breaking for graphene with a centrosymmetric deformation
Abstract
We calculate the local density of states (LDOS) for an infinite graphene sheet with a single centrosymmetric outofplane deformation, in order to investigate measurable strain signatures on graphene. We focus on the regime of small deformations and show that the straininduced pseudomagnetic field induces an imbalance of the LDOS between the two triangular graphene sublattices in the region of the deformation. Real space imaging reveals a characteristic sixfold symmetry pattern where the sublattice symmetry is broken within each fold, consistent with experimental and tightbinding observations. The open geometry we study allows us to make use of the usual continuum model of graphene and to obtain results independent of boundary conditions. We provide an analytic perturbative expression for the contrast between the LDOS of each sublattice, showing a scaling law as a function of the amplitude and width of the deformation. We confirm our results by a numerically exact iterative scattering matrix method.
pacs:
72.80.Vp,73.23.b,72.10.Fk,77.80.bnIntroduction.— Graphene under strain has been largely discussed in the literature and explored for different geometries, with particular features providing alternative routes to confine and control its charge carriers Pereira (); Vozmediano (); Peeters (). Significant development in the theoretical description of strained graphene elucidated how its electronic properties are modified on strained surfaces. At the microscopic level, a general deformation is described by modifications in the atomic positions which reflects in the Hamiltonian as local changes in the hopping parameter CastroNeto (); Papa (). In the continuum model these changes appear as an effective gauge field, and electrons with momentum around the Dirac valleys move in the deformed region as in the presence of a pseudomagnetic field Ando (). Strain also produces a deformation potential, i.e., a scalar field similar to a local chemical potential that can affect electron dynamics Ando (). Very recently, measurements in highquality graphene samples on particular substrates suggested a strong connection between random fluctuations in strain and transport properties Morpurgo ().
The use of strain effects to engineer graphene electronic properties has also been explored in several experiments in the last years Levy (); Klimov (); Mashoff (); Tomori (); JiongLu (); Georgiou (). As one of the most relevant findings, Levy et al. were able to show the presence of pseudo Landau levels generated by giant pseudomagnetic fields induced by homogeneous strain in graphene nanobubbles Levy (). This experimental confirmation that strain can have striking effects on the electronic properties of graphene has been followed by other experiments that explore the effect of different geometries as a path to control graphene electromechanically Tomori (); JiongLu (); Klimov (); Mashoff (). A generic deformation of a graphene sheet can cause inhomogeneous strain, which results in an effective nonuniform pseudomagnetic field and provides an experimental testbed to explore the interplay between highly tunable magnetic fields and Dirac fermions. For example, a scanning tunneling microscope (STM) tip has been used not only to probe samples, but also to continuously deform graphene nanomembranes, demonstrating electronic confinement due to nonuniform pseudomagnetic fields Klimov (). For a similar experimental setup, Mashoff et al. obtained atomically resolved STM images of stable and lifted regions of graphene Mashoff (). Whereas a hexagonal arrangement of the carbon atoms was found at unstrained regions, as expected for monolayer graphene Mashoff (), within the strained area a triangular pattern of bright spots was observed, signaling a symmetry breaking between A and B sublattices in some regions. At the time the authors speculated the effect to be caused by an instability in which the different sublattices acquire a zigzag configuration with respect to the substrate. However, a local sublattice rearrangement requires energies that are prohibitive within the regime of STM imaging making such scenario rather unlikely. Atomistic tightbinding modelsSloan (); Barraza (); Ramon () have predicted the development of such asymmetry but a continuum, symmetrybased description remains missing. Similar patterns were observed in earlier works but remained unexplained alsoXu (). These results indicate an incomplete understanding of the fundamental electronic properties of graphene samples where local manipulation produce effective inhomogeneous gauge and scalar fields.
In this work, we approach this problem by investigating the electronic properties of a graphene sheet in the presence of an axially symmetric outofplane deformation. The strain produced by such distortion is represented by a pseudomagnetic field with trigonal symmetry and embodies a good approximation to standard experimental configurations, while still allowing for analytical treatment. We use a scattering formalism based on the continuum description of graphene to address the question of confinement of electrons due to this deformation. In particular we calculate the LDOS and show that a noticeable imbalance in the distribution of charge density between the two graphene nonequivalent sublattices appears even for small deformations, providing a possible explanation for the experimentally observed sublattice asymmetry. We perform exact numerical calculations and show that these results are well described within an analytical perturbative approach for small deformations. We analyze the dependence of the maximum LDOS contrast between sublattices on energy and strain strength, providing a scaling dependence with the parameters of the deformation. Finally, we also show that the effective scalar field introduced by strain minimally modifies the predicted sublattice asymmetry.
Model.— The electronic properties of undeformed graphene are, for low energies and large system sizes, governed by two copies of a twodimensional (2D) Dirac Hamiltonian where is the velocity of graphene electrons, the electronic momentum around the K (K’) point, and are Pauli matrices reflecting the pseudospin degree of freedom associated with the sublattice structure of the honeycomb lattice Castro (). The strain is produced by a mechanical deformation modeled with a heightprofile that is centrosymmetric and is written generically as , where contains the radial profile, and the parameters and describe amplitude and effective radius of the deformation. In the following, to illustrate our results, we consider the case of a Gaussian bump with height profile . Note however that our results below are qualitatively valid for a generic profile with axial symmetry.
The effect of such deformation on the electronic properties in the continuum limit are described within the theory of elasticity. For an outofplane deformation the strain tensor of elasticity Landau () is derived from the height profile according to , which in polar coordinates () reads
(1) 
where sets the strength of the strain, while its spatial distribution is contained in the function . For the Gaussian profile, one has .
In the presence of the deformation, electrons experience the strain as a gauge field
(2) 
where we chose the zigzag direction to lie along the axis Vozmediano (). The coupling constant is CastroNeto (), being , and and the hopping parameter and the lattice constant of graphene.
For the radial symmetric deformation, the associated pseudomagnetic field shares the trigonal symmetry of the graphene lattice,
(3) 
where the spatial profile is given by the function . For the Gaussianshaped deformation, one has .
In addition to the gauge field, the electrons are also exposed to a scalar potential proportional to the trace of the strain tensor. In our model it is represented by with as a typical valueSloan (). The lowenergy electronic properties in the presence of the deformation are hence described by
(4) 
In this article we consider a bump that is smooth on the scale of the lattice constant, such that a coupling between the valleys can be neglected Sloan (); Ramon (); Blanter (); Bubaloo (). Moreover, we consider an infinite graphene sheet containing a single deformation, hence our results are independent of finite size effects and boundary conditionsBubaloo (); Blanter (); Barraza (); Ramon (); Sloan ().
Perturbation theory.— In this section we present analytic results for the change in the LDOS produced by the scattering of electrons off the deformation obtained with a perturbative approach in real space. We consider therefore small deformations that allow for an expansion in the parameter .
From now on we set , and work around the K valley. The effect of the K’ valley is discussed at the end of this section. We split the Hamiltonian in the kinetic part and the perturbation ,
(5) 
where we defined . We neglect the scalar potential in this part, we will include its effect in the next section.
Let us start with the states of the Dirac equation in the absence of the deformation. Here, we take circular waves as a set of basis states,
(6) 
where denotes the energy of the Dirac fermions (which we assume to be positive, for simplicity), and is a halfinteger index labeling the states according to their angular momentum. denotes the Bessel function of th order. We chose a normalization such that
(7) 
Our goal is to find the scattering state , that replaces when the bump is present. This is determined by the LippmannSchwinger equation
(8) 
which contains the Green’s function of graphene,
(9) 
where we defined
(10) 
are Hankel functions of first and second kind. A derivation of Eq. (9) is given in the Supplementary MaterialSupp ().
For small deformations, we solve the LippmannSchwinger equation in the Born approximation, and replace the scattering state on the righthand side of the equation by the unperturbed state . Note that our perturbative approach is valid for . The explicit form of the resulting scattering states is shown in the Supplementary MaterialSupp (). The trigonal symmetry of the pseudomagnetic field underlies the coupling between angular momentum states differing by .
The LDOS is obtained by calculating . The new states are properly normalized to leading order in , as the linear in correction is orthogonal to the unperturbed state. Since we are interested in the different sublattice occupations, we further introduce the sublatticeresolved LDOS
(11) 
where are projectors on the respective sublattice . For undeformed graphene, evaluating Eq. (11) with the free states produces the wellknown value of
(12) 
for the LDOS per sublattice.
We now want to discuss the effect of the deformation on the LDOS. Specifically we address the limit , which is the relevant case for experiments with a radius of a few lattice constants. In this case, one may simplify the results by using the asymptotic expressions of the Bessel and Hankel functions for small arguments Supp (). Upon retaining only the leading contribution for small energies, one finds for the corrections to leading order in
(13) 
with the function
(14) 
To leading order in , one can replace by in the denominator of Eq. (13). Notice that the relative LDOS correction has no dependence on energy. Thus, the deformation changes the local occupation in the different sublattices in opposite directions, and their spatial distribution shares the symmetry of the underlying pseudomagnetic field with a radial distribution governed by the function . Specifically for a Gaussian height profile, one finds
(15) 
The spatial distribution of the change in LDOS for sublattice A according to Eq. (13) is shown in Fig. 1 (c) for a Gaussian deformation of typical dimensions.
A quantity of experimental relevance is the LDOS contrast between sublattices, defined as
(16) 
which is plotted as a function of the radial distance in Fig. 2 (a). Note that, for a fixed width of the deformation, Eq (13) indicates that the contrast scales quadratically with the amplitude of a centrosymmetric deformation. This scaling is shown for a Gaussian deformation in Fig. 3 and compared with the exact numerical results presented in the next section.
To conclude this section let’s discuss the role of the two valleys K and K’ in these results. As mentioned above, the deformation is smooth enough that does not couple the valleys and their contributions add directly. To see that these are identical, note that the Dirac Hamiltonian takes the same form in both valleys when the spinors are written in the valley symmetric representationBeenakker (): around valley K, and around valley K’. Note that the components referring to and sublattices are interchanged between different valleys. On the other hand, the pseudomagnetic field enters the Dirac equation with opposite sign for each valley (in contrast to a real magnetic field that has the same sign in both valleys). These two effects ensure that their contributions to the sublattice occupancy contrast are identical.
Numerics.— In this section we discuss briefly our exact numerical approach for the continuum model, which is not restricted to the case of small amplitude deformations. The results obtained confirm our findings described in the previous section in the corresponding regimes. We use a slight modification of the method introduced in Ref. Martin, , which allows for the calculation of the scattering matrix for an arbitrarilyshaped scalar potential. The extension to include a vector potential is straightforward. To calculate the LDOS integrated over a certain (arbitrarily chosen) volume per sublattice, we include a fictitious additional scattering potential
(17) 
Such potential locally changes the electron’s energy in by a magnitude in the different sublattices . The LDOS per sublattice integrated over is then found from the scattering matrix via
(18) 
where specifies the sublattice. Fig. 2 (a) shows a comparison between analytical and numerical results for realistic parameters. Note that in the region of small amplitudes the contrast obtained with the expression from perturbation theory follows closely the exact solution given by the numerical approach. Numerical calculations in the presence of the scalar potential induced by the deformation were carried out for different values of the phenomenological parameter . Our results, as shown in Fig. 2 (b), confirm that its effect on the contrast is negligible (note that, to leading order in , affects the occupation of both sublattices in the same way, thus not affecting the contrast (16)). Furthermore, the energy independent value for the contrast predicted with the analytical approach (as long as the requirements and are met) is also verified in the numerical results as shown in Fig. 2 (b).
Conclusions.— We have shown that a centrosymmetric deformation, with its local breaking of the lattice translational symmetry, produces a local sublattice symmetry breaking on the electronic properties of a graphene sheet, and a consequent LDOS contrast between sublattices. Analytic expressions within the Born approximation predict the intensity of the LDOS contrast to be determined by the amplitude of the deformation and to be energy independent for the range of validity of the approximation. Exact numerics carried out with scattering matrix methods confirm the validity of these results, for experimentally realistic parameters. While our numerical approach allows us in principle to treat any size of deformations, we concentrated here on the study of small deformations and showed that there is a measurable LDOS contrast between sublattices even in the absence of Landau levels Castro (). The crossover to this regime will be published elsewhere. Our findings here provide an alternative interpretation for recent experimental observations on STM graphene nanomembranes Alex () and provide a quantitative way to guide the use of strain in the design of electronic properties of graphene flakes.
Acknowledgments.— We gratefully acknowledge support by CNPq, and DAAD (D.F.), DFG SPP 1459 and the A. v H. Foundation (M.S., S.V.K.), and NSF No. DMR1108285 (D.F., N.S.) as well as discussions with R. CarrilloBastos, A. Georgi, M. Morgenstern and P. NemesIncze.
Appendix A Greens functions
In this appendix, we show a derivation of the Green function for graphene given in Eq. 9 in the main text. The defining equation reads
(19) 
For the derivation, it is customary to introduce the auxiliary Green function via
(20) 
The advantage of the Green function is that it satisfies a scalar equation (i.e. the pseudospin degree of freedom has been eliminated).
(21) 
Moreover, is the solution to the Helmholtz equation in two dimensions. We would like to express the Green function in polar coordinates, thus we express the Laplacian as
(22) 
and we write for the function
(23) 
where the summation is running over integer . We then expand as
(24) 
where is determined by
(25) 
For , is thus satisfying a Bessel differential equation. We therefore write
(26) 
with , . In this way, we respect the conditions, that

is regular at ,

is behaving as an outgoing wave () for (this is the proper boundary condition for the retarded Green function ),

is continuous at .
The remaining coefficient is determined by the jump condition for the derivative of ,
(27) 
which yields . Thus, the result for the scalar Green function reads
(28) 
The Green function for graphene is then found in a straightforward albeit lengthy calculation from Eq. (20).
Appendix B Scattering states
The scattering states in the Born approximation are obtained by inserting Eq. 5, Eq. 6 and Eq. 9 into the following equation
(29) 
After insertion of the corresponding expressions, one finds
(30)  
The angular integration is nonzero for and . Calculating the matrix elements, one obtains the scattering states given by
(31) 
with
The leading order correction to the local density of states per sublattice up to first order on , when the deformation is present, is given by
(33)  
with for clean graphene. Using the asymptotic forms of the Bessel functions for small arguments
(34) 
(35) 
one finds the dominant contribution to from the second line, for m=1/2 (+c.c.):
(36) 
The LDOS correction for sublattice B is .
References
 (1) V. M. Pereira and A. H. Castro Neto, Phys. Rev. Lett. 103, 046801 (2009).
 (2) M. A. H. Vozmediano, M. I. Katsnelson, and F. Guinea, Phys. Rep. 496, 109 (2010).
 (3) M. Ramezani Masir, D. Moldovan, F. M. Peeters, Solid State Commun. 175, 76 (2013).
 (4) D. A. Papaconstantopoulos, M. J. Mehl, S. C. Erwin, and M. R. Pederson, in TightBinding Approach to Computational Materials Science, edited by P. Turchi, A. Gonis, and L. Colombo, Materials Research Society, Pittsburgh, 1998, p. 221.
 (5) V. M. Pereira, A. H. Castro Neto and N. M. R. Peres, Phys. Rev.B 80, 045401 (2009).
 (6) H. Suzura and T. Ando, Phys. Rev. B 65, 235412 (2002).
 (7) N. J. G. Couto, D. Costanzo, S. Engels, D. Ki, K. Watanabe, T. Taniguchi, C. Stampfer, F. Guinea, and A. F. Morpurgo, Phys. Rev. X 4, 041019 (2014).
 (8) N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. H. C. Neto, and M. F. Crommie, Science 329, 544 (2010).
 (9) T. Georgiou, L. Britnell, P. Blake, R. V. Gorbachev, A. Gholinia, A. K. Geim, C. Casiraghi, and K. S. Novoselov, Appl. Phys. Lett. 99, 093103 (2011).
 (10) N. N. Klimov, S. Jung, S. Zhu, T. Li, C. A. Wright, S. D. Solares, D. B. Newell, N. B. Zhitenev, and J. A. Stroscio, 336, 1557 (2012).
 (11) T. Mashoff, M. Pratzer, V. Geringer, T. J. Echtermeyer, M. C. Lemme, M. Liebmann, and M. Morgenstern, NanoLetters 10, 461 (2010).
 (12) H. Tomori, A. Kanda, H. Goto, Y. Ootuka, K. Tsukagoshi, S. Moriyama, E. Watanabe, and D. Tsuya, Appl. Phys. Express 4, 075102 (2011).
 (13) J. Lu, A. H. Castro Neto, and K. Ping Loh, Nature Communications 3, 823 (2012).
 (14) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
 (15) L. Landau and E. M. Lifshitz, Theory of Elasticity ( Volumen 7 of A Course of Theoretical Physics ) (Pergamon Press, Cambridge, 1970).
 (16) J. V. Sloan, A. A. P. Sanjuan, Z. Wang, C. Horvath, and S. BarrazaLópez, Phys. Rev. B 87, 155436 (2013).
 (17) A. A. Pacheco Sanjuan, Z. Wang, H. P. Imani, M. Vanevi’c, and S. BarrazaLópez. Phys. Rev. B 89, 121403 (2014).
 (18) R. CarrilloBastos, D. Faria, A. Latgé, F. Mireles, and N. Sandler, Phys. Rev. B 90, 041411(R) (2014).
 (19) K. Xu, P. Cao, and J. R. Heath, Nanolet. 9, 4446 (2009).
 (20) K.J. Kim, Ya. M. Blanter, and K.H. Ahn, Phys. Rev. B 84, 081401(R) (2011).
 (21) G. M. M. Wakker, R. P. Tiwari, and M. Blaauboer, Phys. Rev. B 84, 195427 (2011).
 (22) See Supplementary material for a detailed derivation of the Green’s function in App. A. In App. B, the explicit form of the resulting scattering states is shown.
 (23) C. W. J. Beenakker, Rev. Mod. Phys. 80, 1337 (2008).
 (24) M. Schneider and P. W. Brouwer, Phys. Rev. B 89, 205437 (2014).
 (25) A. Georgi, P. NemesIncze, and M. Moregenstern, Private Communication.