This is the submitted version of the article: Morresi T., Binosi D., Simonucci S., Piergallini R., Roche S., Pugno N.M., Simone T.. Exploring event horizons and Hawking radiation through deformed graphene membranes. 2D Materials, (2020). 7. 041006: - . 10.1088/2053-1583/aba448. Available at: https://dx.doi.org/10.1088/2053-1583/aba448 Forging graphene pseudospheres to mimic curved spacetimes Tommaso Morresi,1,2 Daniele Binosi,1 Stefano Simonucci,3 Riccardo Piergallini,3 Stephan Roche,4∗ Nicola M. Pugno,2,5,6 Simone Taioli1,7∗ 1 European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT*-FBK), Trento, Italy 2 Laboratory of Bio-Inspired & Graphene Nanomechanics - Department of Civil, Environmental and Mechanical Engineering, University of Trento, Italy 3 School of Science and Technology, University of Camerino, Camerino, Italy 4 Institut Catal de Nanociència i Nanotecnologia - ICN2, Barcelona, Spain 5 Ket-Lab, Edoardo Amalfi Foundation, Rome, Italy 6 School of Engineering and Materials Science, Materials Research Institute - Queen Mary University of London, London, UK 7 Trento Institute for Fundamental Physics and Applications (TIFPA-INFN), Trento, Italy Quantum mechanics and general relativity are among density of state (LDOS) for these extended curved struc- the most successful theories of modern physics. Each and tures. Comparison between the numerically evaluated every phenomena that they have predicted, from entan- LDOS and the theoretically predicted one shows, within glement to black holes, has been experimentally tested uncertainties, its thermal nature, establishing the pres- and verified. On the other hand, some tension exists ence of a black hole type horizon in the device. when the two theories are unified in a single framework, which, in turn, makes it extremely challenging to obtain Beltrami’s pseudosphere represents the hyperbolic firm theoretical predictions. Perhaps, the only exception counterpart of the regular sphere: it is a surface of rev- is the discovery by Hawking that, from a quantum me- olution characterized by a constant negative Gaussian chanical point of view, black holes are not completely curvature κ = −1/Rp2 , with Rp the pseudosphere radius. black1 : they emit ‘Hawking radiation’ consisting of pho- Gauss Bonnet’s theorem shows immediately that a car- tons, neutrinos and, to a lesser extent, all sorts of mas- bon allotrope arranged in a pseudosphere shape requires sive particles. However, direct detection of this radiation, the presence of Stone-Wales (SW) defects with heptago- which is thermal in nature, is unlikely: the latter is pro- nal defects exceeding the pentagonal ones by exactly six portional to the inverse of the black hole mass, which, units. In addition, Hilbert’s theorem states that no an- for the smallest observed black hole, implies T = 60 nK alytic complete surfaces of constant negative Gaussian which is 9 orders of magnitude smaller than the cur- curvature can be embedded in R3 , so that a graphene rent cosmic microwave background temperature. On the pseudosphere will be neither compact nor finite. Insights other hand, so-called black hole analogues, first proposed for finding a (local) minimum energy tiling of the Bel- by Unruh2 , are rapidly turning from promising to consol- trami’s pseudosphere has been offered in20 ; but this study idated avenues in the study of various thermodynamics was inconclusive in: i) delivering a general approach to aspects. This is particularly true for sonic analogues built the tessellation of hyperbolic surfaces; ii) scaling-up the from ultracold gases3–12 , for which not only spontaneous graphene pseudosphere size; iii) measuring the surface’s Hawking radiation has been experimentally observed13 , electronic structure. And properties ii) and iii) will be of but also its correlation spectrum shown to be thermal paramount importance in ascertain the capacity of the and with a temperature given by the system’s surface structure to act as an analogue gravity model (see be- gravity14 , thus vindicating Hawking’s predictions. The low). We start the structure generation from a planar state-of-the-art for solid-state black hole analogues is at graphene sheet, in which we force the presence of six a less advanced state15 : while all current experimental heptagonal faces in the center (see Fig. 1a panel i). The approaches face major challenges mainly related to ma- initial configuration of the pseudosphere (Fig. 1a panel terial synthesis and device fabrication, in the last cou- ii) is then obtained by simply projecting the graphene ple of years key conceptual advances have been achieved net on the Beltrami’s surface along the z-axis (see Meth- that show hopes for some of the fundamental questions ods), and shows carbon-to-carbon bond lengths in the to be addressed in these systems too, especially in con- flat (curved) regions which are shorter (longer) than the junction to the implementation of the Sachdev-Ye-Kitaev typical bond distances in graphene (aCC = 0.142 nm). model16,17 . Here we construct a solid-state black-hole A sequence of bond-switching trial moves and structural analogue consisting of a graphene crystal with a three- optimization steps with a modified Keating potential to connected tessellation engineered to shape it in the form favour the formation of hexagonal cells, is then applied of a constant negative curvature surface, known as a Bel- (see Methods), and accepted/rejected according to an trami’s pseudosphere18,19 . We develop a novel computa- energy-minimization criterion (Fig. 1a, panels iii through tional method to build realistic negative curvature car- v). After O(104 ) moves the algorithm efficiency drasti- bon structures comprising several million atoms; and de- cally drops, which limits the radius size of the minimized vise a tight-binding (TB) approach to calculate the local structures (panel vi) to few nm and the number of car- bon atoms to O(103 ). Scaling up of the numbers of atoms 2 a) a) b) b) I)i) ii) II) iii) III) Rp~1.58 Rp ≈ 1.58 nmnm Rp ≈Rp2.74 ~2.74 nm nm N=1,202 N = 1,202 N =N=3,600 3,600 iv)IV) v)V) vi) VI) ... Rp~4.75 nm …. Rp ≈N=10,782 4.75 nm Rp73.96 Rp ≈ ~73.96 nm nm N = 10,782 N =N=2,615,976 2,615,976 FIG. 1. a) Optimization of a carbon pseudosphere containing N = 1626 atoms with radius Rp = 18.26 Å. Panels from i) to vi) represent different optimization stages, consisting of trial bond switch–twist moves. Carbon atoms that do not belong to the hexagonal faces are marked in red color. Starting from an initial configuration (panel i) almost entirely tiled with hexagonal polygons, atoms rearrange to fill uniformly the surface reaching a local minimum after a few thousands steps (panel iv). After reaching this configuration, trial moves are rejected at an increasingly high rate and the Metropolis algorithm becomes inefficient. b) A dualization sequence (described in Methods), leading to a thousandfold increase in the number of carbon atoms and close to a hundredfold increase in the pseudosphere radius. is next implemented through a custom dualization algo- three different structures differing by the number and lo- rithm (Fig. 1b), by which the √ pseudosphere radius and cation of the SW defects. In the energy range E ∈ [−6, 6] number of atoms scale like ∼ 3 and ∼ 3 respectively, eV, the LDOS shows a graphene-like shape for all the while conserving both the bond distance as well as the pseudospheres independently of the radius and defect dis- number of defects (see Methods). Each dualization step tribution. With respect to the pristine graphene region is next followed by a bond switching optimization run i), region ii), shows Van-Hove singularities associated to to counteract the former tendency of splitting apart the the π band peaks which are broadened and shifted; this SW defects of the original structure and thus artificially is due to the slightly elongated carbon bonds character- increasing its total energy. Repeated application of this izing this region, which represents the would-be Hilbert procedure allows to reach a thousandfold increase in the horizon (where the pseudosphere ends as a consequence number of carbon atoms (N = 2, 615, 976) and a pseudo- of Hilbert theorem) of the pseudosphere. After that, the sphere radius Rp = 73.96 nm. LDOS stays the same at a qualitatively level indepen- Next, the electronic properties of these structures can dently from the depth inside the pseudosphere’s ‘throath’ be characterized in terms of the LDOS near the Dirac at which it is evaluated. A similar behaviour (Fig. 3 points21 where electrons behave as relativistic massless a) persists in the biggest structures studied, which also pseudo-particles. Given the O(106 ) carbon atoms of the confirms that the bulge happening in the vicinity of the structures realized, the LDOS will be evaluated through Fermi energy seen in the blown up region in Fig. 3 a), a multi-orbital Tight Binding (TB) approach implement- and which could be also already seen in the central region ing the Kernel Polynomial Method (KPM) to avoid the of Fig. 2), is a genuine curvature effect. Projecting the diagonalization of the Hamiltonian. Due to curvature, LDOS over single atomic sites both inside and outside the the pz orbitals forming the π band are not anymore or- pseudosphere (sites s1−5 in Fig. 3 b)) and separating out thogonal to the in-plane direction; the same happens the nearest neighbours contributions (that would corre- for the sp2 -hybridized orbitals, which do not lay in the spond to the A and B sublattices in pristine graphene), graphene plane. Thus, a multi-orbital TB approach, in we find that the LDOS spectrum around the Fermi en- which all four valence orbitals (2s, 2px , 2py , 2pz ) are in- ergy is significantly asymmetric for the two nonequivalent cluded in the simulation as opposed to the pz orbital sublattices while for E  EF it is practically indistin- alone, has been employed (see Methods, where details guishable. concerning the parametrization of the Hamiltonian are The largest simulated pseudosphere has a ratio reported and well-established results on graphene and aCC /Rp ∼ 2×10−3 , which quantifies how well the Hilbert carbon nanotube structures reproduced). The LDOS horizon and the Rindler-type event horizon emerging projected onto longitudinal circles in regions located at a when treating the pseudosphere as a 2+1 dimensional different z-depth along pseudospheres obtained at various space-time coincide19 (0 representing coincidence). In stages of the dualization procedure (and thus character- the ideal case18,19 , in the vicinity of the horizon Hawking ized by varying number of atoms N and radius Rp ), is radiation takes place with the massless electron’s action plotted in Fig. 2. In each case we evaluate the LDOS for before (respectively after) tunnelling has taken place be- 3 a) a) Rp ~ 73.97 nm Rp ~R14.26 p ≈ 14.26 nm nm N = 96,930 (a) Rp ≈ 73.96 nm N=2,615,976 N = 2,615,976 v N=96,930 iv v ii iii iii ii iii iv i ii iv i i iii ii iv v i v b) Rp ~R24.66 p ≈ 24.66 nm nm b) Rp ~ 73.97 nm v N = 290,736 N=290,736 (b) Rp ≈ 73.96 nm N=2,615,976 Rp ~ 73.97 N =nm 2,615,976 iv N=2,615,976 s5 iii s5 ii iii s4 s3 s5 ii s3 s4 s1 s2 i i s3 s4 s1 s2 iv s2 v s1 c) Rp R ≈ 42.70 ~ p42.7 nm nm v N = 876,046 N=876,046 (c) FIG. 3. a) Same as before but for the three biggest pseudo- iv spheres (N = 2, 615, 976 and Rp = 73.96 nm, right). We also iii ii iii zoom on the Fermi energy to expose the LDOS bulge in this ii region. b) Pseudosphere’s LDOS projected over the shown i i sites si . Projection over the two would-be inequivalent sub- iv lattices A and B of graphene is shown; the insets zoom near v the Fermi energy to show the LDOS asymmetry. FIG. 2. Panels a) through c): LDOS projected onto five re- event horizon is present, then it would be revealed by gions at different depth along the revolution axis z (i through the the thermal nature of the power spectrum, since18,19 , v) for various pseudospheres characterized by three different neglecting boundary terms, values of the number of atoms N , radius Rp . For each case of N and Rp we report the LDOS for three pseudosphere re- 4 1 Rp2 −2u/Rp E ~vF aCC u/r alizations B1−3 differing by the configurations of SW defects, ρ= e ; Θ= e , π (~vF )2 a2CC E e κΘ − 1 κ 2πr2 of which a representative is shown on the right of each panel. (2) The Fermi energy is set equal to zero in all cases. where κ is the Boltzman constant and Θ the temperature, which reaches its maximum at the horizon where r = Rp , ing u = Rp ln Rp /aCC , and Θ = ~vF /2πκRp ∼ 16 K for our largest pseudosphere. Notice in addition that only √ Z S = ivF d3 x g ψ̄γ µ Dµ ψ, (1) electrons with an energy E = ~vF /Rp ∼ 9 meV have a long enough wavelength to experience the effects of the curvatures and thus 2. Given the exponential nature with: vF ∼ c/300 the Fermi velocity; γ the Dirac of the power spectrum (2), the effect of the presence of matrices; Dµ the SO(2,1) (respectively, flat) covariant the horizon should manifest in a marked asymmetry of derivative; and, finally, g the determinant of the pseu- the LDOS around the Fermi energy, as measured by the dosphere metric (respectively 1) with (conformal) met- contrast ric ds2B = φ2 (u)ds2R , where φ(u) = `/r eu/r (with ` a constant that in the physical case is to be indentified LDOSE<0 − LDOSE>0 with aCC ) and ds2R is the Rindler-type metric ds2R = ν= . (3) LDOSE<0 + LDOSE>0 φ−2 (u)(dt2 − du2 ) − r2 dv 2 . (u and v are the curvilin- ear coordinates spanning the pseudosphere). Hawking Results for the LDOS and its contrast projected on sites radiation is then characterized by a LDOS obtained by located near the Hilbert/Rindler horizon are shown in evaluating the power spectrum ρ of the 2-point function Fig. 4. The behavior of the curves is remarkably sim- SB = h0M |ψB ψ̄B |0M i with |0iM the flat vacuum. If an ilar; recall, in fact, that the continuum approximation 4 a)a) Theoretical Numerical sublattices can be neither avoided nor disentangled. In particular notice that the contrast while being off in mag- nitude is qualitatively the same. In summary, the Beltrami’s pseudosphere tiled by car- bon atoms arranged in a defected graphene net may rep- resent a viable analogue model of a quantum field theory in curved space-time in general, and a black-hole hori- zon in particular. The success of our numerical effort in verifying the analytical predictions obtained within bb)) a continuum representation, suggests that this analogue system, if experimentally realized (for example through direct optical forging24 ), can allow measuring a Hawking temperature of the order of tens of K by measuring the surface LDOS through either low temperature Scanning Tunnlling Microscope or optical near-field spectroscopy. ACKNOWLEDGMENTS FIG. 4. a) LDOS lineshape and contrast projected over the pseudosphere sites s4,5 shown in Fig. 3 b), according to analytical predictions22,23 and numerical results. N.M.P. is supported by the European Commission un- der the Graphene Flagship Core 2 grant No. 785219 (WP14, “Composites”) and FET Proactive (“Neurofi- formula (2) is valid when neglecting both local elastic bres”) grant No. 732344 as well as by the Italian Ministry effects induced by the curvature as well as the presence of Education, University and Research (MIUR) under the of SW defects; and that in our computational model, “Departments of Excellence” grant L.232/2016. The au- the interplay between curvature and defects (which the thors acknowledge Bruno Kessler Foundation (FBK) for Gauss-Bonnet theorem makes inseparable) is inherently providing unlimited access to the KORE computing fa- present and, thus, its effect on the LDOS of the graphene cility. 1 11 S. W. Hawking, Black hole explosions, Nature 248, 30 P. E. Larre, A. Recati, I. Carusotto, N. Pavloff, Quantum (1974). fluctuations around black hole horizons in Bose-Einstein 2 W. G. Unruh, Experimental black hole evaporation, Phys. condensates, Phys. Rev. A85, 013621 (2012). 12 Rev. Lett. 46, 1351 (1981). J. Steinhauer, Measuring the entanglement of analogue 3 T. A. Jacobson, G. E. Volovik, Event horizons and ergore- Hawking radiation by the density-density correlation func- gions in He-3, Phys. Rev. D58, 064021 (1998). tion, Phys. Rev. D92, 024043 (2015). 4 13 L. J. Garay, J. R. Anglin, J. I. Cirac, P. Zoller, Black holes J. Steinhauer, Observation of quantum Hawking radiation in Bose-Einstein condensates, Phys. Rev. Lett. 85, 4643 and its entanglement in an analogue black hole, Nature (2000). Phys. 12, 959 (2016). 5 14 C. Barcelo, S. Liberati, M. Visser, Towards the observation J. R. M. de Nova, K. Golubkov, V. I. Kolobov, J. Stein- of Hawking radiation in Bose-Einstein condensates, Int. J. hauer, Observation of thermal Hawking radiation and its Mod. Phys. A18, 3735 (2003). temperature in an analogue black hole, Nature 569, 688 6 S. Giovanazzi, Experimental black hole evaporation, Phys. (2019). 15 Rev. Lett. 94, 061302 (2005). M. Franz, M. Rozali, Mimicking black hole event horizons 7 R. Balbinot, A. Fabbri, S. Fagnocchi, A. Recati, I. Caru- in atomic and solid-state systems, Nature Rev. Mater. 3, sotto, Non-local density correlations as signal of Hawking 491 (2018). 16 radiation in BEC acoustic black holes, Phys. Rev. A78, S. Sachdev, J. Ye, Gapless spin fluid ground state in a 021603 (2008). random, quantum Heisenberg magnet, Phys. Rev. Lett. 70, 8 I. Carusotto, S. Fagnocchi, A. Recati, R. Balbinot, A. Fab- 3339 (1993). 17 bri, Numerical observation of Hawking radiation from A. Kitaev, A simple model of quantum holography, KITP acoustic black holes in atomic BECs, New J. Phys. 10, Strings Seminar and Entanglement Program . 18 103001 (2008). A. Iorio, G. Lambiase, The Hawking-Unruh phenomenon 9 J. Macher, R. Parentani, Black hole radiation in Bose- on graphene, Phys. Lett. B716, 334 (2012). 19 Einstein condensates, Phys. Rev. A80, 043601 (2009). A. Iorio, G. Lambiase, Quantum field theory in curved 10 P. N. Recati, A., I. Carusotto, Bogoliubov Theory of acous- graphene spacetimes, Lobachevsky geometry, Weyl sym- tic Hawking radiation in Bose-Einstein Condensates, Phys. metry, Hawking effect, and all that, Phys. Rev. D90, Rev. A80, 043603 (2009). 025006 (2014). 5 20 S. Taioli, R. Gabbrielli, S. Simonucci, N. M. Pugno, A. Io- 97, 170201 (2006). 28 rio, Lobachevsky crystallography made real through car- H. Gould, J. Tobochnik, W. Christian, An introduction to bon pseudospheres, Journal of Physics: Condensed Matter computer simulation methods (1988). 29 28, 13LT01 (2016). S. Alfthan, Computational studies of silicon interfaces and 21 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. amorphous silica, G5 artikkeliväitöskirja (2006). 30 Novoselov, A. K. Geim, The electronic properties of J. C. Slater, G. F. Koster, Simplified lcao method for the graphene, Rev. Mod. Phys. 81, 109 (2009). periodic potential problem, Phys. Rev. 94, 1498 (1954). 22 31 A. Iorio, G. Lambiase, The hawking-unruh phenomenon S. Yuan, M. Rösner, A. Schulz, T. O. Wehling, M. I. Kat- on graphene, Physics Letters B 716, 334 (2012). snelson, Electronic structures and optical properties of par- 23 A. Iorio, G. Lambiase, Quantum field theory in curved tially and fully fluorinated graphene, Phys. Rev. Lett. 114, graphene spacetimes, lobachevsky geometry, weyl symme- 047403 (2015). 32 try, hawking effect, and all that, Phys. Rev. D 90, 025006 T. Stauber, J. I. Beltran, J. Schliemann, Tight-binding (2014). approach to penta-graphene, Scientific Reports 6, 22672 24 A. Johansson, et al., Optical Forging of Graphene into (2016). 33 Three-Dimensional Shapes, Nano Lett. 17, 6469 (2017). P. Giannozzi, et al., Quantum espresso: a modular and 25 A. Kumar, M. Wilson, M. F. Thorpe, Amorphous open-source software project for quantum simulations of graphene: a realization of zachariasen’s glass, Journal of materials, Journal of Physics: Condensed Matter 21, Physics: Condensed Matter 24, 485003 (2012). 395502 (2009). 26 34 F. Wooten, K. Winer, D. Weaire, Computer generation of A. Weiße, G. Wellein, A. Alvermann, H. Fehske, The kernel structural models of amorphous si and ge, Phys. Rev. Lett. polynomial method, Rev. Mod. Phys. 78, 275 (2006). 35 54, 1392 (1985). J. Garcia, The kernel polynomial method for quantum 27 E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumb- transport in disordered systems, Ph.D. thesis (2015). sch, Structural relaxation made simple, Phys. Rev. Lett. METHODS Tiling the pseudosphere The tiling of the Beltrami’s pseudosphere by carbon atoms, which represents an interesting geometrical problem on its own right, has been achieved through the following steps: 1. Set the length of the pseudosphere by fixing the maximum value of the coordinate along the axis of revolution (z). 2. Determine the number of carbon atoms N that are needed if one were to tile the surface of the Beltrami’s pseudosphere with the same density of planar graphene (0.379 atoms/Å2 ). Periodic boundary conditions are also introduced to saturate the outer carbon atom bonds belonging to r = Rp and z = 0 (the Hilbert horizon); finally, a rectangular supercell periodically repeated along the x and y directions, respectively. 3. Construct a planar graph (N, F, E) consisting of N vertices, F faces and E edges. The N vertices represent compressed carbon atoms with shortened carbon-to-carbon bond lengths, aCC <1.42 Å; each vertex is linked to three nearest neighbours by edges (representing bonds) and is shared by three faces. 4. Map the initial graph onto the Beltrami’s pseudosphere surface via a one-to-one p transformation, by which the revolution axis coordinate z of the vertices is unambiguously determined ∀ x + y 2 < Rp by fixing 2 "s s # x2 + y 2 x2 + y 2 z = z(x, y) = R 1− − atanh 1 − ; Rp2 Rp2 5. Find the atomic arrangements with N ∼ 1000-2000 that minimize a surface potential energy of the Keating type25 : 3 α X 2 2 2 3 β X  a2CC 2 X E= 2 rij − aCC + 2 rij · rik + + ch (|Fi | − 6)2 (4) 16 aCC i,j 8 aCC 2 i,j,k Fi where α = 25.88 eV Å −2 is the bond stretching force constant, aCC = 1.42 Å is the carbon-to-carbon bond length in planar graphene, rij is the distance between atoms i and j, and β ∼ α5 is the bond–bending force 6 a) A 6 A A 5 a) C 6 B C 6 B C 7 7 B 6 D D 5 D A 6 A A 5 a) A A A C 6 B C 6 B B C 7 7 B (i) 6 DD (ii) D D (iii) 5 D D b) (i) (ii) (iii) b) b) FIG. 5. a) An example of a bond-switch trial move in a graphene lattice. Starting from an initial configuration showing 6-fold rings (dashed circles), the bonds between the carbon atoms A–B, and C–D are cut, and new bonds A–C and B–D are formed. Twisting then transforms the four 6-fold rings into two 5-fold rings and two 7-fold rings (adapted from Ref.29 ). b) Dualization sequence based on the three-connectivity of the graph. The parent parent geometry is shown on the left, while the resulting (daughter) geometry is on the right. constant. The last term favours the formation of hexagonal faces; as Fi labels the polygons of the net |Fi | is the number of vertices of the polygons (ch = 0.35 empirically). To reach the energy minimum repeat the following steps, typically O(104 ) times: • perform random switchings/twists of atomic bonds, based on the Wooten, Winer and Weaire (WWW) method26 (Fig. 5 a); • let the geometry relax through molecular dynamics simulations based on the Fast Inertial Relaxation Engine (FIRE) approach27 ; • accept the move only if it lowers the total energy of the system according to a Metropolis algorithm28 . 6. Execute on the minimized surfaces a dualization sequence, to increase the number N of atoms and correspond- ingly the radius of the pseudosphere ((Fig. 5 b). By using the three-connectivity of the graph√one creates a hexagon around each vertex of the initial optimized structure; rescale the bond lengths with a 3 factor and repeat from 5. Tight Binding parameter estimate Low energy electronic properties of geometries containing millions of atoms, have been evaluated using a TB approach, which is well known to describe correctly the dispersion of graphene around the six Dirac K-points in the first Brillouin zone (IBZ)21 . Due to the pseudosphere curvature, a multi-orbital TB approach has been developed, in which all four valence orbitals (2s, 2px , 2py , 2pz ) are included in the simulation through the Hamiltonian: iξ a†i,ξ ai,ξ + † X X H= tij ξ,γ ai,ξ aj,γ , (5) ξi ξ,γhiji 7 a) a/acc=0.96 a/acc=1 9 9 6 6 3 3 Energy [eV] Energy [eV] 0 0 -3 -3 -6 -6 -9 -9 -12 DFT -12 DFT -15 TB -15 TB -18 -18 -21 -21 Γ M K Γ Γ M K Γ a/acc=1.04 a/acc=1.08 9 9 6 6 3 3 Energy [eV] Energy [eV] 0 0 -3 -3 -6 -6 -9 -9 -12 DFT -12 DFT -15 TB -15 TB -18 -18 -21 -21 Γ M K Γ Γ M K Γ FIG. 6. TB fit of graphene bands obtained from DFT electronic structure calculations. From left to right and from top to bottom we report the bands for a biaxial compression of the graphene cell equal to 4 %, unstrained graphene, and a biaxial tensile strain of the cell equal to 4 %and 8%. where ξ, γ are orbital label indices while i, j are site indices; tij ξ,γ indicates the hopping parameters; and the symbol hiji means that the nearest neighbours approximation is adopted. The parameters tij ξ,γ describing the hopping between orbitals in different sites were computed within the Slater-Koster formulation30 , which provides a scheme to relate the orbital symmetry, distances and directions of neighbour atoms. Owing to the non-planarity of our geometry we cannot make use of the multi-orbital parametrization typically used for graphene31,32 where the onsite energy of the pz -symmetry orbitals are treated differently from the x, y orbital cartesian components along the in-plane directions (that is px = py 6= pz ). Therefore, we derive the TB parameters by fitting ab-initio Density Functional Theory (DFT) simulations of the graphene bands by further imposing that the onsite energies for the p orbitals are the same (px = py = pz ). DFT simulations of equilibrium and strained configurations of graphene were carried out by using the Quantum Espresso code suite33 ; in particular we use a norm–conserving PBE pseudopotential (C.pbe-mt gipaw.UPF) and an energy cut-off for the wavefuntion expansion on plane-waves set equal to 100 Ry. The k-point mesh is a 40401 grid for the calculation of both the ground state density and the band structures. Convergence of the integrals over the Brillouin zone was improved by smearing the occupancy with a 0.136 eV width Gaussian function. The TB parameters that we obtained using Eq. (5) for unstrained (a ≡ aCC ) and strained graphene a−aCC (a 6= aCC ) are (see Fig. 6): s = −2.8 eV and p = 0 eV as onsite energies; Vssσ (a) = −5.6 · a 6= aCC · e− 0.55 eV, a−aCC a−aCC a−aCC Vspσ (a) = 5.2 · a/aCC · e− 0.75 eV, Vppσ (a) = 4.6 · a/aCC · e− 0.55 eV and Vppπ (a) = −2.44 · a/aCC · e− 0.41 eV as hopping parameters between different orbitals. Kernel Polynomial Method For the evaluation of the Local Density of States (LDOS) we resorted to The Kernel Polynomial Method (KPM), which is a numerical approach useful to access spectral quantities of extended systems for which a direct diagonalization of the full Hamiltonian matrix is computationally unfeasible. It consists in the expansion of the sought quantity in terms of a set of orthogonal polynomials, and then in improving the convergence of the expansion with a kernel to avoid spurious Gibbs oscillations34 . In particular, we used Chebyshev polynomials for the expansion, and Jackson kernel to increase convergence, resolution, and accuracy34 . Under these conditions a generic function can be expanded according to c −1 " NX # 1 f (x) = √ µ0 g0 + µn gn Tn (x) , (6) π 1 − x2 n=1 8 a) b) b) bc)) Nc=2400 Nc=3600 Nc=4800 a) LDOS [arb. units] Nc=6000 -0.03 -0.015 0 0.015 0.03 Energy [eV] a) zoom around EF=0 -20 -15 -10 -5 0 5 10 15 Energy [eV] b) rad=0.61 nm rad=3.32 nm rad=6.71 nm rad=13.49 nm TDOS [arb. units] rad=20.27 nm zoom around EF=0 -3 -2 -1 0 1 2 3 Energy [eV] FIG. 7. a) (top) LDOS for different values of the cut-off parameter Nc obtained by projection on a site of a pristine graphene rectangular cell with sides equal to 199.97 nm × 199.81 nm, containing 1, 525, 188 carbon atoms; The inset represents a zoom near the Fermi energy. (bottom) total DOS of (n,n) nanotubes for n=9 (r = 0.61 nm), n=49 (r = 3.32 nm), n=99 (r = 6.71 nm), n=199 (r = 13.49 nm) and n=299 (r = 20.27 nm). b) The LDOS angular and radial dependence obtained by projecting over sites belonging to a SW defect (colored atoms in the right panel). c) The LDOS radial dependence. R1 where Tn (x) are Chebyshev polynomials of the first kind, µn = −1 dx f (x)Tn (x) are the coefficients of the expansion and the gn are the Jackson kernel coefficients defined as   1 πn πn π gn = (Nc − n + 1) cos + sin cot , (7) Nc + 1 Nc + 1 Nc + 1 Nc + 1 with Nc the truncation number related to the maximum momentum. The best achievable resolution through this kernel is r J π ∆ = 1 − cos . (8) Nc + 1 We refer to34 for the details about the calculation of µn ; here it suffices to emphasize that it is based on the Stochastic Evaluation of Traces, which requires a certain number R of random initial states. As expected, the bigger is R, the more accurate becomes the evaluation of the coefficients; we found that R = 100 was enough for all calculations carried out. Tests of the LDOS calculation The convergence with respect to the Nc parameter can be tested in the calculation of the LDOS for the benchmarks cases of planar graphene and armchair carbon nanotubes. The LDOS for four different values of Nc ranging from 2400 to 6000, is shown in Fig. 7 a) top. While at a wide energy scale the curves are indistinguishable, zooming near the Fermi energy (set to zero as usual), shows higher truncation values for Nc captures more faithfully the expected linear dispersion relation; on the other hand, there is a threshold to the number of terms in the summation after which spurious oscillations sets is, spoiling convergence. This can be understood by noticing that the energy separation between levels in periodic graphene is infinitesimal and the DOS is a continuous function. Then, since the pseudosphere in our simulations is a large but finite system and the energy separation of the levels increases with respect to infinite periodic structures, a too big value of Nc may result in a KPM energy resolution marginally above the finite energy separation between levels of our finite system, thus leading to poor convergence35 . For non-planar systems, we have computed the total DOS of (n, n) nanotubes for n = 9, 99, 199 and 299 (radius r = 0.61, 3.32, 6.71 9 and 13.49 nm). The total DOS is reported in the bottom panel of Fig. 7 a), where we observe that the DOS lineshape of these armchair nanotubes is reproduced surprisingly well already for the moderate value of Nc = 2000 and that, as expected, the confinement effects become less important upon increasing the radius size. Since a method for estimating the value of Nc that trades-off between accuracy and computational efficiency exists only for pristine structures that do not have any defect35 , selecting the best Nc is a trial and error process. For the pseudosphere case we found Nc = 8000 to be the optimal value (with Nc = 16000 used when resolving the Fermi energy region in Fig. 4). LDOS of a single SW defect in planar graphene Let us finally investigate the effect of the presence of a single Stone-Wales (SW) defect, like the ones present in the pseudosphere, on the LDOS of a graphene net (N = 823, 860). In particular, we study both the LDOS projected over different symmetry sites of the SW defect (with Nc = 5000), thus obtaining information on the angular dependence (Fig. 7 b), as well as the LDOS projected over sites increasingly far from the SW defect, thus obtaining insights on the radial dependence (Fig. 7 c). On top of a marked angular dependence, we see that the shape of the LDOS is dramatically modified near the defect site, while far from it the planar graphene shape is recovered; the presence of the SW defect still affects the LDOS projected at distances of ≈ 80 Å with small oscillations in the spectrum. Near the Fermi energy one observes in addition a marked asymmetry persisting again up to a distance of ≈ 80 Å .