Observing a scale anomaly and a universal quantum phase transition in graphene
O. Ovdat, Jinhai Mao, Yuhang Jiang, E. Y. Andrei, E. Akkermans
OObserving a scale anomaly and a universal quantum phasetransition in graphene
O. Ovdat , Jinhai Mao , Yuhang Jiang , E. Y. Andrei , and E. Akkermans Department of Physics, Technion – Israel Institute of Technology, Haifa 3200003, Israel Department of Physics and Astronomy, Rutgers University,Piscataway, New Jersey 08854
Abstract
One of the most interesting predictions resulting fromquantum physics, is the violation of classical symmetries,collectively referred to as anomalies. A remarkable classof anomalies occurs when the continuous scale symmetryof a scale free quantum system is broken into a discretescale symmetry for a critical value of a control parame-ter. This is an example of a (zero temperature) quantumphase transition. Such an anomaly takes place for thequantum inverse square potential known to describe ’Efi-mov physics’. Broken continuous scale symmetry intodiscrete scale symmetry also appears for a charged andmassless Dirac fermion in an attractive 1 /r Coulombpotential. The purpose of this article is to demonstratethe universality of this quantum phase transition andto present convincing experimental evidence of its exis-tence for a charged and massless fermion in an attractiveCoulomb potential as realised in graphene.Continuous scale symmetry (CS) – a common prop-erty of physical systems – expresses the invariance of aphysical quantity f ( x ) (e.g., the mass) when changing acontrol parameter x (e.g., the length). This property isexpressed by a simple scaling relation, f ( ax ) = b f ( x ),satisfied ∀ a > b ( a ), whose gen-eral solution is the power law f ( x ) = C x α with α =ln b/ ln a . Other physical systems possess the weakerdiscrete scale symmetry (DS) expressed by the sameaforementioned scaling relation but now satisfied forfixed values ( a, b ) and whose solution becomes f ( x ) = x α G (ln x/ ln a ), where G ( u + 1) = G ( u ) is a periodicfunction. Physical systems having a DS are also knownas self-similar fractals [1] (Fig. 1a). It is possible to breakCS into DS at the quantum level, a result which consti-tutes the basis of a special kind of scale anomaly [2, 3].A well studied example is provided by the problem ofa particle of mass µ in an attractive inverse square po-tential [4, 5] which plays a role in various systems [6–9]and more importantly in Efimov physics [10, 11]. Al-though well defined classically, the quantum mechanicsof the scale – and conformal [12] – invariant Hamilto-nian H = − ∆ / µ − ξ/r (with (cid:126) = 1) is well posed, ∗ These authors contributed equally to this work.Correspondence should be addressed to E. Akkermans([email protected]) but for large enough values of ξ , H is no longer self-adjoint [13, 14]. The corresponding Schr¨odinger equa-tion for a normalisable wave function ψ ( r ) of energy k = − µE is, ψ (cid:48)(cid:48) ( r ) + d − r ψ (cid:48) ( r ) + ζr ψ ( r ) = k ψ ( r ) (1)where ζ ≡ µξ − l ( l + d −
2) is a dimensionless parameter, d the space dimensionality and l the orbital angular mo-mentum. Equation (1) is invariant under the transfor-mation r → λ r and k → k/λ , ∀ λ (CS), namely to everynormalisable wave function of energy k corresponds acontinuous family of states with energies ( λk ) , so thatthe bound spectrum is a continuum unbounded from be-low. Various ways exist to cure this problem, based oncutoff regularisation and renormalisation group [15–21],and all lead for the low energy spectrum to a quantumphase transition (QPT) monitored by ζ , between a singlebound state for ζ < ζ c to an infinite and discrete energyspectrum for ζ > ζ c , independent of the regularisationprocedure and given by k n ( ζ ) = (cid:15) e − πn √ ζ − ζ c , n ∈ Z (2)which clearly displays DS. The critical value ζ c = ( d − / (cid:15) is a regularization dependent energy scale. In the over-critical phase ζ > ζ c , the corresponding renormaliza-tion group solution provides a rare example of a limitcycle [15, 16, 22]. Building on the previous example,it can be anticipated that the problem of a masslessDirac fermion in an attractive Coulomb potential [23–25]), − Zα/r , is also scale invariant (CS) and that thespectrum of resonant quasi-bound states presents similarfeatures and a corresponding QPT.In this work, we demonstrate the existence of sucha universal QPT for arbitrary space dimension d ≥ s -wave channel. All these features are ex-perimentally demonstrated using a charged vacancy in1 a r X i v : . [ c ond - m a t . m e s - h a ll ] S e p 𝜆𝐸 = 𝜆 −1 𝜌 𝐸 (c) (b) 𝛽 < 𝛽 𝑐 𝑉 𝑟 = −𝑍𝛼 𝑟 𝒓𝑽 𝒓 V bias V gate G/GBN SiO Si 𝐸 𝑛 = −𝜀 𝑒 − 𝜋𝑛𝛽 −𝛽 𝑐2 𝛽 > 𝛽 𝑐 𝒓𝑽 𝒓 (a)(f)(d) (e) 𝑑𝐼/𝑑𝑉(a.u.)𝑑𝐼/𝑑𝑉(a.u.) Figure 1:
Schematic visualization of the purpose and main results of this paper . (a) Sierpinski gasketas typical featuring of such iterative fractal structures. This QPT is realized experimentally by creating single-atom vacancies in graphene. The function ρ ( E ) is the density of states and obeys a scaling relation characterisingthe existence of discrete scale symmetry. (b)-(c) Illustration of the universal Quantum Phase Transition (QPT)obtained by varying the dimensionless parameter β ≡ Zα (see text for precise definitions) in the low energyspectrum of a massless fermion in a Coulomb potential V = − Zα/r created by a charge Z . (b) For low values β < β c , there is a single quasi-bound state close to zero energy. (c) For overcritical values β > β c , the low energyspectrum is a ladder E n characterized by a discrete scale symmetry { E n } = { λE n } for λ = exp( π/ (cid:112) β − β ).(d)-(e) Experimental dI/dV maps of charged vacancy for fixed β < β c (d) and β > β c (e). The images illustratethe characteristic probability density of the resonances in (b) and (c). (f) Scanning tunnelling microscopy (STM)setup. Local charge Z is accumulated at the single vacancy in graphene by applying voltage pulses to the STMtip.graphene. We observe the overcritical spectrum and weobtain an experimental value for the universal geomet-ric ladder factor in full agreement with the theoreticalprediction. We also explain the observation of two inter-twined ladders of quasi bound states as resulting fromthe breaking of parity symmetry. Finally, we relate ourfindings to Efimov physics as measured in cold atomicgases. Results
The Dirac model.
The Dirac equation of a masslessfermion in the presence of a − Zα/r potential is obtainedfrom the Hamiltonian (with (cid:126) = c = 1), H = − iγ γ j ∂ j − βr (3)where ( γ , γ j ) are Dirac matrices. Here the dimension-less parameter monitoring the transition is β = Zα ,where Z is the Coulomb charge and α the fine struc-ture constant. The QPT occurs at the critical value β c = ( d − / (see Supplementary Note 1). For reso- A related anomalous behaviour in the Dirac Coulomb prob-lem has been identified long ago [26] but its physical relevance nant quasi-bound states, we look for scattering solutionsof the form ψ in + e iη ψ sc , where η ( E ) is the energy-dependent scattering phase shift and ψ in , sc ( r, E ) are twocomponent objects representing the radial part of theDirac spinor which behave asymptotically as, ψ in , sc ( r, E ) = r − d (cid:0) V in , sc (2 i | E | r ) ∓ iβ e ∓ iEr (cid:1) (4)for | E | r (cid:29) γ ≡ (cid:112) β − β , ψ in , sc ( r, E ) = r − d (cid:16) U − in , sc (2 iEr ) − iγ + U +in , sc (2 iEr ) iγ (cid:17) , (5)for | E | r (cid:28) V in , sc and U ± in , sc in equations (4),(5) are constants. It is easy to inferfrom (5) that β = β c plays a special role. Indeed for β > β c , there exists a family of normalisable solutionswhich admit complex eigenvalues E = − i(cid:15) , hence theHamiltonian (3) is not self-adjoint ( H (cid:54) = H † ). To prop-erly define this quantum problem, a regularisation is was marginal since it required non existent heavy-nuclei Coulombcharges Z (cid:39) /α (cid:39)
137 to be observed. Moreover, the problem ofa massive Dirac particle is different due to the existence of a finitegap which breaks CS. β = Zα . This is achieved by introducing acutoff length L and a boundary condition at r = L ,which is equivalent to replacing the Coulomb potentialat short distances by a well behaved potential whose ex-act form is irrelevant in the low energy regime EL (cid:28) h = Ψ ( r, E ) / Ψ ( r, E ) | r → L + , where Ψ , representthe two components of the aforementioned radial partof the Dirac spinor. The resulting scattering phase shift η ( E, L, h ), which contains all the information about theregularisation, thus becomes a function of L and of theparameter h . The quasi-bound states energy spectrumis obtained from the scattering phase shift by means ofthe Krein-Schwinger relation [27, 28] which relates thechange of density of states δρ to the energy derivative of η , δρ ( E ) = 1 π dη ( E ) dE . (6) Theoretical structure of quasi-bound spec-trum . From now on, and to compare to experimentalresults further discussed, we consider the case d = 2for which there is a single orbital angular momentumquantum number m ∈ Z . The corresponding criticalcoupling becomes β c = | m + 1 / | ≥ /
2, giving rise tothe s -wave channels, m = 0 , − β c = 1 /
2. De-pending on the choice of boundary condition h , δρ ( E )can be degenerate or non-degenerate over these two s -wave channels. This degeneracy originates from thesymmetry of the (2 + 1) Dirac Hamiltonian (3) underparity, ( x, y ) → ( − x, y ), and its existence is equivalentto whether or not the boundary condition breaks par-ity (see Supplementary Note 2). In what follows, we willconsider the generic case in which there is no degeneracy.In the undercritical, β < β c , and low energy regime EL (cid:28)
1, we observe (see Figs. 1b and 2a) a singlequasi-bound state originating from only one of the s -wave channels and which broadens as β increases. Inthe overcritical regime β > β c , this picture changesdramatically. The low energy ( EL (cid:28)
1) scatteringphase shift displays two intertwined, infinite geometricladders of quasi-bound states (Figs. 1c and 3) at ener-gies E n still given by (2) but with ζ − ζ c now replacedby β − β . This sharp transition at β c belongs tothe same universality class as presented for the inversesquare Schr¨odinger problem, namely CS of the quasi-bound states spectrum is broken for β > β c into a DSphase characterized by a fractal distribution of quasi-bound states. The QPT thus reflects the lack of self-adjointness of the Hamiltonian (3) and the necessary reg-ularisation procedure leads to a scale anomaly in whichCS is broken into DS. Experimental realization in graphene.
A par-ticularly interesting condensed matter system where theprevious considerations seem to be relevant is graphenein the presence of implanted Coulomb charges in con- This is also related to the Wigner time delay [29] and to theFriedel sum rule We emphasize that this picture remains valid for all values of β > β c and not only in the vicinity of β c4 Moreover, note that the energy scale (cid:15) for the Dirac case isdifferent from the inverse square Schr¨odinger case defined in (1). Π ¶ Η ¶ E - - Β ‡ Β ‡ Β ‡ H a L - - Β =
Β =
Β = H b L d I d V @ a . u . D E H eV L E - E D H eV L Figure 2:
Experimental and theoretical picturein the undercritical regime . (a) Theoretical be-haviour of π dη/dE for d = 2 showing quasi-bound statesof a massless Dirac fermion in the undercritical regime β < /
2. In the scale-free low energy EL (cid:28) m = − m = 0 (purple) branch shows no peak independentlyof the choice of boundary condition (see SupplementaryNote 2). While increasing β , the resonance shifts tolower energy and becomes broader. (b) Excitation spec-trum measured in graphene using STM as a function ofthe applied voltage V . The determination of the param-eter β is explained in the text. - - - - - Π ¶Η¶E E H eV L - - - - - - - - - - -
100 1001021041061081010 Π ¶Η¶E E H eV L Figure 3:
Theoretical behaviour of the low en-ergy and scale free part of the quasi-bound statesspectrum in the overcritical regime for d = 2 and β = 1 . > β c (= 1 / The lower plot displays the de-tailed structure of the infinite geometric ladders. Notethat the m = − m = 0 (purple) laddersare intertwined. These results are independent of theboundary condition.3eniently created vacancies [30]. It is indeed knownthat low energy excitations in graphene behave as amassless Dirac fermion field with a linear dispersion (cid:15) = ± v F | p | and a Fermi velocity v F (cid:39) m / s [31].These characteristics have been extensively exploited tomake graphene a very useful platform to emulate spe-cific features of quantum field theory, topology and espe-cially QED [23], since an effective fine structure constant α G ≡ e / (cid:126) v F of order unity is obtained by replacing thevelocity of light c by v F .It has been recently shown that single-atom vacanciesin graphene can stably host local charge [30]. Densityfunctional theory (DFT) calculations have shown thatwhen a carbon atom is removed from the honeycomb lat-tice, the atoms around the vacancy site rearrange into alower energy configuration [32]. The resulting lattice re-construction causes a charge redistribution which in theground state has an effective local charge of ≈ +1. Re-cent Kelvin probe force microscopy measurements of thelocal charge at the vacancy sites are in good agreementwith the DFT predictions. Vacancies are generated bysputtering graphene with He + ions [33, 34]. Charge ismodified and measured at the vacancy site by means ofscanning tunnelling spectroscopy and Landau level spec-troscopy as detailed in [30]. Applying multiple pulses al-lows for a gradual increase in the vacancy charge, whichin turn acts as an effective tunable Coulomb source.Moreover, the size of the source inside the vacancy issmall ( ≈ β = 1 / β values.To establish a relation between the measured differ-ential conductance and the spectrum of quasi-boundstates, we recall that the tunnel current I ( V ) is pro-portional to both the density of states ρ t ( (cid:15) ) of theSTM tip and ρ ( (cid:15) ) of massless electronic excitations ingraphene at the vacancy location. We also assumethat the tunnel matrix element | t | depends only weaklyon energy and that both voltage and temperature aresmall compared to the Fermi energy and height of thetunnelling potential, so that the current I ( V ) = G t V is linear with V thus defining the tunnel conductance G t = 2 π (cid:0) e / (cid:126) (cid:1) | t | ρ t ρ ( (cid:15) ). Assuming that ρ t of the ref-erence electrode (the tip) is energy independent, a vari-ation δρ ( (cid:15) ) of the local density of states at the vacancyleads to a variation δI ( V ) of the current and thus to avariation δG t ( V ) of the tunnel conductance so that, atzero temperature, we obtain [36] δG t ( V ) G t = δρ ( (cid:15) ) ρ , (7)where ρ is the density of states in the absence of va-cancy. By considering the vacancy as a local pertur-bation, each quasi-particle state is characterised by itsscattering phase shift taken to be the phase shift η ( E )of the quasi-bound Dirac states previously calculated.Then, the change of density of states δρ ( E ) is obtainedfrom (6) and combining together with (7) leads to the - - -
50 0246 246 E E E '1 l og H ¶ Η ¶ E L E H meV L - - -
50 013 13 E E '1 E d I (cid:144) d V E - E D H meV L Figure 4:
Experimental and theoretical picturein the overcritical regime.
Upper plot: Theoreticalbehaviour of the low energy and scale free part of theovercritical ( β = 1 .
33) quasi-bound states spectrum ob-tained from (6). The blue (purple) line corresponds to m = − m = 0). Lower plot: Experimental values ofthe (STM) tunnelling conductance measured at the po-sition of charged vacancies in graphene. The labelling E , E , E (cid:48) of the peaks is explained in the text.relation, dδIdV = G t πρ dη ( E ) dE (8)between the differential tunnel conductance and thescattering phase shift.The measurements and data analysis presented herewere carried out as follows: positive charges are grad-ually injected into an initially prepared single atom va-cancy and the differential conductance δG t ( V ) is mea-sured at each step as a function of voltage. Since we arelooking at the positions of resonant quasi-bound states,both quantities displayed in Figs. 2 and 4 give the sameset of resonant energies, independently of the energy-independent factor G t /πρ . For low enough values of thecharge, the differential conductance displayed in Fig. 2b,shows the existence of a single quasi-bound state reso-nance. The behaviour close to the Dirac point, namelyin the low energy regime independent on the short dis-tance regularization, is very similar to the theoreticalprediction of Fig. 2a. When the build up charge exceedsa certain value, we note the appearance of three reso-nances, emerging out of the Dirac point. We interpretthese resonances as the lowest overcritical ( β > /
2) res-onances which we denote E , E (cid:48) , E respectively. Thecorresponding theoretical and experimental behavioursdisplayed in Figs. 3, 4, show a very good qualitativeagreement. To achieve a quantitative comparison solelybased on the previous Dirac Hamiltonian (3), we fix L and the boundary condition h and deduce the theoret-ical β values corresponding to the respective positionsof the lowest overcritical resonance E (as demonstratedin Fig. 4). This allows to determine the lowest branch E ( β ) for n = 1 represented in Fig. 5. Then, the ex-perimental points E (cid:48) , E are directly compared to their4 ‰ - - - E - E D H e V L Β E H Β L E H Β L E '1 H Β L ‰‰ ‰‰ Figure 5:
Behaviour of the energies E n ( β ) of thequasi-bound state spectrum. The curves are ob-tained from (2) for E ( β ) , E (cid:48) ( β ) , E ( β ) as adapted tothe massless Dirac case. The black and cyan dots cor-respond to the values measured in graphene. The twopink x’s are the values of Efimov energies measured inCaesium atoms [37, 38] which corresponds to the (over-critical) fixed Efimov value β E = 1 . L and h , according to the ansatz h = a ( m +1),and obtain the best correspondence for L (cid:39) . a (cid:39) − .
85. We compare the experimental E /E ratiowith the universal prediction E n +1 /E n = e − π/ √ β − / as seen in Fig. 6. A trend-line of the form e − b/ √ β − / is fitted to the ratios E /E yielding a statistical valueof b = 3 .
145 with standard error of ∆ b = 0 .
06 consis-tent with the predicted value π . An error of ± mV isassumed for the position of the energy resonancesA few comments are appropriate: (i) The points onthe E ( β ) curve follow very closely the theoretical pre-diction E n +1 /E n = e − π/ √ β − / . This result is insensi-tive to the choice of h , thus manifesting the universalityof the ratio E n +1 /E n . (ii) In contrast, the correspon-dence between the E (cid:48) points and the theoretical branchis sensitive to the choice of h . This reflects the factthat while each geometric ladder is of the form (2) (withthe appropriate ζ → β change), the energy scale (cid:15) isdifferent between the two thus leading to a shifted rel-ative position of the two geometric ladders in Fig. 3.The ansatz taken for h is phenomenological (see Sup-plementary Note 2), however, we find that in order toget reasonable correspondence to theory, the explicit de-pendence on m is needed. More importantly, it is nec-essary to use a degeneracy breaking boundary condi-tion to describe the E (cid:48) ( β ) points. For instance, if theCoulomb potential is regularised by a constant potentialfor r ≤ L [41], then both angular momentum channels(i.e., the E (cid:48) and E points) become degenerate. The ex-istence of the experimental E (cid:48) branch is therefore a dis-tinct signal that parity symmetry in the correspondingDirac description (3) is broken. In graphene, exchang-ing the triangular sublattices is equivalent to a paritytransformation. Creating a vacancy breaks the symme-try between the two sub-lattices and is therefore at the origin of broken parity in the Dirac model. (iii) Thevalue L (cid:39) . E L/ (cid:126) v F (cid:39) . (cid:28) β -driven QPT. Β E (cid:144) E Figure 6:
Comparison between the experimen-tally obtained E /E ratio and the universal fac-tor e − π/ √ β − / . Blue points: the ratio E /E ob-tained from the position of the points in Fig. 5. Cyanpoint: Universal Efimov energy ratio as measured inCaesium atoms [37, 38]. Blue line (dashed): the corre-sponding optimized curve, fitted according to the model e − b/ √ β − / and corresponding to b = 3 .
145 with stan-dard error of ∆ b = 0 .
06 consistent with the predictedvalue π . The shaded pink region is the ± b confi-dence interval of the curve. Cyan line: universal lowenergy factor e − π/ √ β − / . Purple line: theoretical ra-tio E /E obtained from the exact solution of the Diracequation. As β → . | E n | becomes smaller thereforethe green and purple curves coincide for low β . The errorbar on the resonance energies is ± mV . Discussion
A further argument in support of the universality of thisQPT is achieved by comparing the experimental resultsobtained in graphene with those deduced from a com-pletely different physical problem. To that purpose, wedwell for a short while recalling the basics underlyingEfimov physics [42]. Back to 1970, Efimov [10] stud-ied the quantum problem of three identical nucleons ofmass m interacting through a short range ( r ) poten-tial. He pointed out that when the scattering length a of the two-body interaction becomes very large, a (cid:29) r ,there exists a scale free regime for the low energy spec-trum, (cid:126) /ma (cid:28) E (cid:28) (cid:126) /mr , where the correspond-ing bound states energies follow the geometric series( √− E n = − ˜ (cid:15) e − πn/s ) where s (cid:39) . (cid:15) a problem-dependent energyscale. Efimov deduced these results from an effectiveSchr¨odinger equation in d = 3 with the radial ( l = 0) at-tractive potential V ( r ) = − (cid:0) s + 1 / (cid:1) /r . Using equa-tions (1), (2) and the critical value ζ c = ( d − / / ζ value forthe Efimov effect to be s + 1 / > ζ c correspondingto the overcritical regime of the QPT. The value of β e π/s is β E = (cid:112) s + 1 / . E n ( n = 1 ,
2) have been recently determined using an ul-tracold gas of caesium atoms [38]. Although the Efimovspectrum always lies at a fixed and overcritical value ofthe coupling, unlike the case of graphene where β can betuned, the universal character of the overcritical regimeallows nevertheless for a direct comparison of these twoextremely remote physical systems. To that purpose, weinclude the Efimov value β E in the expression obtainedfor the massless Dirac fermion in a Coulomb potentialand insert the corresponding data points obtained forcold atomic caesium in the graphene plot (Fig. 5) upto an appropriate scaling of ˜ (cid:15) . The results are fullyconsistent thus showing in another way the universalitypresented.There are other remote examples of systems display-ing this universal QPT, e.g., flavoured QED3 [48], andthe XY model (Kosterlitz-Thouless [8] and rougheningtransitions [22]). Our results provide a useful and origi-nal probe of characteristic features of this universal QPTand motivate a more thorough study of this transition. Methods
Our sample is stacked two layers of graphene on top ofa thin BN flake (see Fig. 1f). The standard dry transferprocedure is followed to get this heterostructure. A largetwisted angle between the two layers graphene is selectedin order to weaken the coupling. The free-standing likefeature for the top layer graphene is checked by the Lan-dau levels spectroscopy. To achieve the diluted single va-cancies, the sample is exposed to the helium ion beam forshort time (100 eV for 5 s) followed by the high temper-ature annealing. The experiment is performed at 4.2 Kwith a home-built STM. The dI/dV ( I is the current, V is the bias) is recorded by the standard lock-in tech-nique, with a small AC modulation 2 mV at 473.1 Hzadded on the DC bias. To tune the effective charge onthe vacancy, we apply the voltage pulse (-2 V, 100 ms)with the STM tip directly locating on top of the vacancy. Data availability
The data that support the find-ings of this study are available from the correspondingauthor upon request.
Acknowledgments
This work was supported by the Israel Science Founda-tion Grant No. 924/09. Funding for the experimentalwork provided by DOE-FG02-99ER45742 (STM/STS),NSF DMR 1207108 (fabrication and characterization).
Author contributions
O. Ovdat and E. Akkermans proposed observing theaforementioned quantum phase transition in graphenewith a charged vacancy. They have contributed to in-terpreting and solving the theoretical model as well asanalysing the experimental data and making contactwith the theory. J. Mao, Y. Jiang and E. Y. Andrei con-ceived of and designed the experiment, as well as per-formed the measurements and analysed the data. Allauthors contributed to discussions and preparation ofthe manuscript.
Competing Financial Interests
The authors de-clare no competing financial interests.
References [1] E. Akkermans, in
Fractal Geometry and Dynami-cal Systems in Pure and Applied Mathematics II:Fractals in Applied Mathematics , Vol. 601, editedby D. Carf, M. L. Lapidus, E. P. J. Pearse, andM. van Frankenhuijsen (American MathematicalSociety (AMS), 2013) pp. 1–21.[2] S. L. Adler, Phys. Rev. , 2426 (1969).[3] J. S. Bell and R. Jackiw, Il Nuovo Cimento A (1965-1970) , 47 (1969).[4] K. M. Case, Phys. Rev. , 797 (1950).[5] L. D. Landau, Quantum mechanics : non-relativistic theory (Butterworth-Heinemann, Ox-ford Boston, 1991).[6] J.-M. L´evy-Leblond, Phys. Rev. , 1 (1967).[7] H. E. Camblong, L. N. Epele, H. Fanchiotti, andC. A. Garc´ıa Canal, Phys. Rev. Lett. , 220402(2001).[8] D. B. Kaplan, J.-W. Lee, D. T. Son, and M. A.Stephanov, Phys. Rev. D , 125005 (2009).[9] C. Nisoli and A. R. Bishop, Phys. Rev. Lett. ,070401 (2014).[10] V. Efimov, Physics Letters B , 563 (1970).[11] V. Efimov, Sov. J. Nucl. Phys , 589 (1971).[12] R. W. Jackiw, Diverse topics in theoretical andmathematical physics (World Scientific, 1995).[13] K. Meetz, Il Nuovo Cimento (1955-1965) , 690(1964).[14] D. M. Gitman, I. Tyutin, and B. L. Voronov, Self-adjoint Extensions in Quantum Mechanics: Gen-eral Theory and Applications to Schr¨odinger andDirac Equations with Singular Potentials , Vol. 62(Springer, 2012).[15] S. Albeverio, R. Høegh-Krohn, and T. T. Wu,Phys. Lett. A , 105 (1981).616] S. R. Beane, P. F. Bedaque, L. Childress, A. Kry-jevski, J. McGuire, and U. van Kolck, Phys. Rev.A , 042103 (2001).[17] E. J. Mueller and T.-L. Ho, arXiv:cond-mat/0403283 (2004).[18] E. Braaten and D. Phillips, Phys. Rev. A , 052111(2004).[19] H.-W. Hammer and B. G. Swingle, Annals ofPhysics , 306 (2006).[20] S. Moroz and R. Schmidt, Annals of Physics ,491 (2010).[21] A. De Martino, D. Kl¨opfer, D. Matrasulov, andR. Egger, Phys. Rev. Lett. , 186603 (2014).[22] E. B. Kolomeisky and J. P. Straley, Phys. Rev. B , 12664 (1992).[23] A. V. Shytov, M. I. Katsnelson, and L. S. Levitov,Phys. Rev. Lett. , 246802 (2007).[24] V. M. Pereira, J. Nilsson, and A. H. Castro Neto,Phys. Rev. Lett. , 166802 (2007).[25] Y. Nishida, Phys. Rev. B , 165414 (2014).[26] I. Pomeranchuk and J. Smorodinsky, J. Phys.(Moscow) , 100 (1945).[27] E. Akkermans, Journal of Mathematical Physics , 1781 (1997).[28] E. Akkermans, G. Dunne, and E. Levy, in Opticsof Aperiodic Structures: Fundamentals and DeviceApplications , edited by L. D. Negro (Pan StanfordPublishing, 2013).[29] F. T. Smith, Phys. Rev. , 349 (1960).[30] J. Mao, Y. Jiang, D. Moldovan, G. Li, K. Watan-abe, T. Taniguchi, M. R. Masir, F. M. Peeters, andE. Y. Andrei, Nat Phys , 545 (2016).[31] E. Y. Andrei, G. Li, and X. Du, Reports onProgress in Physics , 056501 (2012).[32] Y. Liu, M. Weinert, and L. Li, Nanotechnology ,035702 (2015).[33] O. Lehtinen, J. Kotakoski, A. V. Krasheninnikov,A. Tolvanen, K. Nordlund, and J. Keinonen, Phys.Rev. B , 153401 (2010).[34] J.-H. Chen, L. Li, W. G. Cullen, E. D. Williams,and M. S. Fuhrer, Nat Phys , 535 (2011).[35] Y. Wang, D. Wong, A. V. Shytov, V. W. Brar,S. Choi, Q. Wu, H.-Z. Tsai, W. Regan, A. Zettl,R. K. Kawakami, S. G. Louie, L. S. Levitov, andM. F. Crommie, Science , 734 (2013).[36] E. Akkermans and G. Montambaux, in MesoscopicPhysics of Electrons and Photons (Cambridge Uni-versity Press, 2007) Chap. 7. [37] T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl,C. Chin, B. Engeser, A. D. Lange, K. Pilch,A. Jaakkola, H.-C. Ngerl, and R. Grimm, Nature , 315 (2006).[38] B. Huang, L. A. Sidorenkov, R. Grimm, and J. M.Hutson, Phys. Rev. Lett. , 190401 (2014).[39] S.-K. Tung, K. Jim´enez-Garc´ıa, J. Johansen, C. V.Parker, and C. Chin, Phys. Rev. Lett. , 240402(2014).[40] R. Pires, J. Ulmanis, S. H¨afner, M. Repp, A. Arias,E. D. Kuhnle, and M. Weidem¨uller, Phys. Rev.Lett. , 250404 (2014).[41] V. M. Pereira, V. N. Kotov, and A. H. Castro Neto,Phys. Rev. B , 085101 (2008).[42] E. Braaten and H.-W. Hammer, Physics Reports , 259 (2006).[43] S. E. Pollack, D. Dries, and R. G. Hulet, Science , 1683 (2009).[44] N. Gross, Z. Shotan, S. Kokkelmans, andL. Khaykovich, Phys. Rev. Lett. , 163202(2009).[45] T. Lompe, T. B. Ottenstein, F. Serwane, A. N.Wenz, G. Z¨urn, and S. Jochim, Science , 940(2010).[46] S. Nakajima, M. Horikoshi, T. Mukaiyama,P. Naidon, and M. Ueda, Phys. Rev. Lett. ,143201 (2011).[47] M. Kunitski, S. Zeller, J. Voigtsberger, A. Kalinin,L. P. H. Schmidt, M. Sch¨offler, A. Czasch,W. Sch¨ollkopf, R. E. Grisenti, T. Jahnke, D. Blume,and R. D¨orner, Science , 551 (2015).[48] T. Appelquist, D. Nash, and L. C. R. Wijeward-hana, Phys. Rev. Lett. , 2575 (1988).7e present in supplementary figures 1 and 2 all overcritical and undercritical measurements, and their corre-sponding theoretical plots. In supplementary figure 3 we present the STM topography image of an isolated vacancyin graphene.Supplementary Figure 1: Experimental and theoretical picture in the undercritical regime. a. Theo-retical behaviour of the undercritical ( β < /
2) quasi-bound states spectrum of massless Dirac fermions. b. STMmeasurement of the quasi-particle spectrum at the position of the charged vacancy in the undercritical regime. Seemain text for more details.
Supplementary Note 1: Dirac Coulomb problem in d +1 dimensions – critical coupling,phase shift and quasi bound states In what follows, we study a system described by a massless Dirac particle in the presence of an electric potentialthat has an inverse radial tail in d + 1 dimensions. We show that beyond a critical coupling value, an anomalousbreaking of conformal symmetry occurs and the system is described by a discrete scale invariant spectrum. Weobtain for a general over critical coupling and short range behaviour of the potential, an expression for an infiniteseries of geometrically spaced quasi bound states.The massless Dirac equation with an attractive potential V ( r ) = − βr , β ≡ Zα in d + 1 dimensions is iγ µ ( ∂ µ + ieA µ ) ψ ( x ν ) = 0 (1)where µ = 0 . . . d , A µ is the electromagnetic potential (EM) eA = − β/rA i = 0 i = 1 , . . . , d. (2)1upplementary Figure 2: Experimental and theoretical picture in the overcritical regime. a. Theoreticalbehaviour of the overcritical ( β > /
2) quasi-bound states spectrum of massless Dirac fermions. b. STM measure-ment of the quasiparticle spectrum at the position of the charged vacancy in the overcritical regime. See main textfor more details.and γ µ are d + 1 matrices satisfying the anti-commutation relation { γ µ , γ ν } = 2 η µν (3)with η µν being the d +1 Minkowski metric with a ’mostly minus’ sign convention. The Hamiltonian of the systemis expressed as H = γ γ j p j − β/r (4)where j = 1 . . . d and corresponds to the scale invariant eigenvalue equation Hψ = Eψ equivalent to (1).Utilizing rotational symmetry, the angular part of supplementary equation (1) can be solved and the radialdependence of ψ ( x ν ) is given in terms of two functions Ψ ( r ) , Ψ ( r ) [1] determined by the following set of equationsΨ ( r ) + ( d − K )2 r Ψ ( r ) = (cid:18) E + βr (cid:19) Ψ − Ψ ( r ) − ( d − − K )2 r Ψ ( r ) = (cid:18) E + βr (cid:19) Ψ (5)where K ≡ ( ± (cid:0) l + d − (cid:1) d > m + 1 / d = 2 , (6) l = 0 , , . . . and m ∈ Z are orbital angular momentum quantum numbers. In terms of these radial functions, thescalar product of two eigenfunctions ψ, ˜ ψ is given by Z dV ψ † ˜ ψ = Z dr r d − (cid:16) Ψ ∗ ( r ) ˜Ψ ( r ) + Ψ ∗ ( r ) ˜Ψ ( r ) (cid:17) . We introduce a short distance radial cut-off L and assume that there exist an electric Coulomb potential V ( r ) = − βr for r > L and some interaction at r < L that can be modelled by a BC at r = L . The equivalent mixedboundary condition of (5) can be written as follows [2] h = Ψ ( r )Ψ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r → L + (7)where h is determined by the short range physics and in general can depend on E, L and K . Specification of thecut-off L , coupling β , boundary condition h (and the already determined angular dependence) determine a specificset of solutions to equation (1).Two independent solutions to equations (5) are given by S ± γ with2upplementary Figure 3: Characteristic topography signature of an isolated vacancy in graphene.
Thetriangular interference pattern arises due to the local crystal distortion and corresponding electronic state recon-struction. This is the feature that was used to identify single atom vacancies in this work. S γ ( ρ ) ≡ r − d e − ρ/ ρ iγ (cid:20) F ( i ( γ + β ); 1 + 2 iγ ; ρ ) (cid:18) i (cid:19) + γ + βK F (1 + i ( γ + β ); 1 + 2 iγ ; ρ ) (cid:18) i (cid:19)(cid:21) (8)where F ( a, b, z ) is Kummer’s function [3], γ ≡ p β − K and ρ ≡ iEr . The | E | r (cid:28) S ± γ is givenby S γ = r − d ρ iγ (cid:18)(cid:18) i γ + βK i + γ + βK (cid:19) + O ( | E | r ) (cid:19) . (9)Solutions corresponding to outgoing and ingoing radial waves for r → ∞ are given by the combinations ψ in = (cid:0) C γ in S γ + C − γ in S − γ (cid:1) = r − d (cid:18)(cid:18) − i (cid:19) e − iEr − iβ log(2 | E | r ) + O (cid:18) | E | r (cid:19)(cid:19) ψ sc = (cid:0) C γ sc S γ + C − γ sc S − γ (cid:1) = r − d (cid:18)(cid:18) i (cid:19) e + iEr + iβ log(2 | E | r ) + O (cid:18) | E | r (cid:19)(cid:19) (10)where C γ sc ≡ Ke − πβ (cid:0) e πβ − e − πγ (cid:1) Γ( i ( γ + β ))( e πγ − e − πγ ) Γ(1 + 2 iγ ) C γ in ≡ e πγ e πβ (cid:0) e − πβ − e − πγ (cid:1) Γ(1 + i ( γ − β ))( e πγ − e − πγ ) Γ(1 + 2 iγ ) (11)For the case E <
0. The log dependence in (10) is characteristic of the long range Coulomb tail and is irrelevant tothe physics of the scattering problem [4]. Thus, a general solution to (5) can be written as (cid:18) Ψ Ψ (cid:19) ∝ ψ in ( ρ ) + e iη ψ sc ( ρ ) (12)where the scattering phase shift η is determined by the boundary condition (7) at r = Le iη ( EL,h ) = − ψ in, ( ρ ) − hψ in, ( ρ ) ψ sc, ( ρ ) − hψ sc, ( ρ ) (cid:12)(cid:12)(cid:12)(cid:12) ρ =2 iEL (13)3he energy derivative of η is given by the logarithmic derivative of the RHS of (13) dηdE = 12 i ddE ln − ψ in, ( ρ ) − hψ in, ( ρ ) ψ sc, ( ρ ) − hψ sc, ( ρ ) (cid:12)(cid:12)(cid:12)(cid:12) ρ =2 iEL ! (14)The explicit expression is complicated and is therefore omitted. Resonant quasi-bound states appear for specificvalues of the energy E at which dη/dE exhibits a sharp maxima. We would like to obtain an analytic expressionfor the position of the resonances in the regime | E | L (cid:28) E → ε ≡ E R − i W [5] and look for solutions of (5) and (7) with no e − iEr plane wavesolution for r → ∞ . The lifetime of the resonance is given by W − and is required to be positive. A solution with noingoing wave is obtained by the requirement that e − iη ( EL,h ) = 0 or alternatively the vanishing of the denominatorin (13) such that h = ψ sc, ( ρ ) ψ sc, ( ρ ) (cid:12)(cid:12)(cid:12)(cid:12) ρ =2 iεL . (15)Focusing our attention at the regime | ε | L (cid:28)
1, one finds very different results depending on the value of β . For β < β c ≡ | K | , γ is pure imaginary. Thus, from (9) and (10) we get that both ψ sc, , ∝ ρ − √ K − β and thereforeequation (15) is independent of ε to leading order in | ε | L . This inconsistency means that for fixed L and β there areno quasi bound states arbitrarily close to zero energy. In contrary, for β > β c and to leading order in | ε | L , equation(15) is given by h = C γs ρ iγ (cid:16) i + β + γK (cid:17) + ( γ → − γ ) C γs ρ iγ (cid:16) i β + γK (cid:17) + ( γ → − γ ) (16)where h ≡ h | E → . Solving supplementary equation (16) for ρ iγ gives ρ iγ = z where z ≡ C − γs ((1 − ih )( β − γ ) − ( h − i ) K ) C γs (( h − i ) K − (1 − ih )( β + γ )) .Inserting ρ ≡ iεL and solving for ε yields ε n = − (cid:15) e i Θ e − πnγ n ∈ Z (17)where ε ≡ L | z iγ | > ≡ arg (cid:18) iz iγ (cid:19) . The regime of validity of this result is for n > n min , where n min issuch that | ε n min | L ∼
1. Using Γ − function identities, the phase Θ of − ε n can be simplified toΘ = arg (cid:18) iz iγ (cid:19) = π − γ log( | z | )= − π γ log (cid:18) sinh ( π ( β + γ ))sinh ( π ( β − γ )) (cid:19) . (18)Note that Θ is independent on the boundary condition h . Recall that γ ≡ p β − K , where K is given in (6) for d ≥
2. For β > | K | ≥ , arg ( − ε n ) is a very slowly varying function bounded in the region 0 < arg ( − ε n ) < . π ,as can be seen in supplementary figure 4. This is consistent with negative resonances E n = Re ( ε n ) < W n ≡ − ε n ) >
0. Therefore, in the regime | ε | L (cid:28) ε n describes a geometrically spaced set ofinfinitely many negative quasi bound states for all β > | K | , h, d ≥
2. The width of the resonance is given by W n ≡ − ε n ) < . π | ε n | which describes resonances getting sharper indefinitely as n → ∞ . In the vicinity ofthe transition point β − | K | &
0, arg ( − ε n ) ∼ = − π − e | K | π + O ( ξ − | K | ) / which for d = 2 , m = 0( ⇔ | K | = ) agreeswith the result of reference [6] in the main text. Supplamentary Note 2: Dirac Coulomb system as a model for a charged vacancy ingraphene
As explained in the main text, we model excitations around the vacancy as massless Dirac particles and neglect anyinteractions between them. Therefore, for r > L , L being the characteristic size of the vacancy, they are describedby the Dirac equation (5) with d = 2. The full spinor in this case is given in polar coordinates by ψ ( r, φ ) = e imφ (cid:18) Ψ ( r ) i Ψ ( r ) e iφ (cid:19) (19)in a representation where H = ~σ · ~p − β/r . The components of the spinor correspond to amplitudes of the wavefunction on each of the two graphene sublattices. In the near vacancy region, r < L we assume some unknowninteraction that can be taken into account by the boundary condition (7) at r = L . Conventionally used choices are:4upplementary Figure 4: The phase of − ε n as a function of the coupling β and | K | for β > | K | . The boldblack line represents ξ = | K | . It can be seen that the phase of − ε n takes values in the range [0 , . π ] for d ≥ V < = − β/L [7] corresponding to h = J m +1 ( β + EL ) /J m ( β + EL ) , where J n ( x ) is Bessel’s function; (ii) zero wavefunction on one of the graphene lattice sites [8] corresponding to h = 0; (iii) infinite mass term on boundary term [9] corresponding to h = 1. Generically, h can depend on E, L and m . From (6), the critical coupling is β c ≡ | m + 1 / | ≤ /
2, giving rise to two angular momentum s-wave channels, m = 0 , − β c = 1 / P y , in which ( x, y ) → ( − x, y ) [10]. The action of the parity operators on Ψ( r, φ ) is defined as ψ ( r, φ ) ≡ P y ψ ( r, π − φ )= σ y ψ ( r, π − φ )= − ie i ( − m − φ (cid:18) Ψ ( r ) − i Ψ ( r ) e iφ (cid:19) . (20)This transformation can also be accounted for (up to an unimportant overall phase) byΨ ( r ) → Ψ ( r ) , Ψ ( r ) → − Ψ ( r ) , m → − m − K ≡ m + 1 /
2, thus m → − m − ⇔ K → − K . Indeed, the Dirac equation (5), is invariant under (21).From (21) it is apparent that solutions corresponding to angular momentum m and − m − m and − m −
1. The answer depends on the short distance region. Forexample, if we describe the r < L regime by condition (i) then the quasi bound states will necessarily be degenerateover the m and − m − r > L , r < L respectsparity symmetry. As a result, h in (i) transforms like Ψ ( r ) / Ψ ( r ) under (21), i.e., h → − h − . Thus, by applying(21) on both side of (14) it is straightforward to obtain that dη m /dE = dη − m − /dE . However, if the potential inthe r < L regime break parity, that is, the corresponding boundary condition does not transform as h → − h − under (21) (for example h = 0 , dη m /dE = dη − m − /dE and the degeneracy will be broken. Specifically, for m = 0 , − β > β c = 1 / h and therefore is sensitiveto the detail of the short range physics. 5upplementary Figure 5: The Landau levels spectrum for twisted bilayer graphene under . (a) dI/dVcurve on graphene at B = 10T with V b = − I = 20pA. (b) Fit of the Landau level sequence in (a) usedto extract the Fermi velocity.All the figures in the main text describing the behaviour of dη/dE as a function of E and the values of the quasibound states E n ( β ) are extracted from the exact relation appearing in equation (14). Supplementary methods
Sample fabrication . In this work, we use G/G/BN on SiO to perform the experiment. The hBN thin flakes wereexfoliated onto the SiO surface followed by a dry transfer process using a sacrificial PMMA thin film to stack thefirst graphene layer on the hBN flake. Before stacking the top layer graphene, the PMMA was removed with acetoneand IPA, followed by furnace annealing in forming gas (10% H and 90% Ar) at 230 ◦ C for 3 hours. The secondlayer graphene was stacked by using the same procedure as the first layer. Au/Ti electrodes were deposited by thestandard SEM lithography for the STM contact. After the liftoff process, the sample was annealed again in furnacewith forming gas to remove the PMMA residues. Subsequently, the sample is loaded in the UHV chamber for furtherannealing at 230 ◦ C overnight. To generate the single vacancies, the sample is exposed to a 100 − + ionbeam followed by high temperature in situ annealing. The other stacked samples, G/BN/SiO and G/G/SiO , wereprepared by the same procedure. Characterization by STM topography and Landau level (LL) spectroscopy . The STM experiment isperformed at 4 .
2K using a cut PtIr tip. The dI/dV spectroscopy is performed using the standard lockin method withbias modulation typically 2mV at 473 . substrate by using an intermediate graphene layerand an hBN buffer underneath the layers. When the twist angle between the two stacked graphene layers exceeds10 ◦ the two layers are electronically decoupled at the experimentally relevant energies. Therefore a large twist anglewas chosen to ensure a linear dispersion near the Dirac point. LL spectroscopy provides a direct way to prove thelayer decoupling. For single layer graphene, the energy level sequence is given by: E N = sign( N ) v F p e ~ | N | B ,where N is the LL index, v F is the Fermi velocity, ~ is the reduced Plank constant. For fixed magnetic field, theroot N dependence for the sequence is the fingerprint of single layer graphene. supplementary figure 5(a) shows theLL spectrum of the G/G/BN sample at 10T. By fitting the LLs sequence (supplementary figure 5(b)), we obtainthe value of the Fermi velocity v F = (1 . ± . × m/s .6 upplementary References [1] S.-H. Dong, Wave Equations in Higher Dimensions (Springer, 2011).[2] C. N. Yang, Comm. Math. Phys. , 205 (1987).[3] M. Abramowitz and I. A. Stegun,
Handbook of mathematical functions: with formulas, graphs, and mathematicaltables , 55 (Courier Corporation, 1964).[4] J. R. Walter Greiner, B. Mller,
Quantum Electrodynamics of Strong Fields , edited by W. Greiner (Springer-Verlag Berlin Heidelberg, 1985).[5] H. Friedrich,
Scattering Theory (Springer-Verlag Berlin Heidelberg, 2013).[6] A. V. Shytov, M. I. Katsnelson, and L. S. Levitov, Phys. Rev. Lett. , 246802 (2007).[7] V. M. Pereira, V. N. Kotov, and A. H. Castro Neto, Phys. Rev. B , 085101 (2008).[8] A. V. Shytov, M. I. Katsnelson, and L. S. Levitov, Phys. Rev. Lett. , 236801 (2007).[9] V. M. Pereira, J. Nilsson, and A. H. Castro Neto, Phys. Rev. Lett. , 166802 (2007).[10] R. WINKLER and U. ZLICKE, The ANZIAM Journal57