Observing a scale anomaly and a universal quantum phase transition in graphene

One of the most interesting predictions resulting from quantum physics, is the violation of classical symmetries, collectively referred to as anomalies. A remarkable class of anomalies occurs when the continuous scale symmetry of a scale-free quantum system is broken into a discrete scale symmetry for a critical value of a control parameter. This is an example of a (zero temperature) quantum phase transition. Such an anomaly takes place for the quantum inverse square potential known to describe ‘Efimov physics’. Broken continuous scale symmetry into discrete scale symmetry also appears for a charged and massless Dirac fermion in an attractive 1/r Coulomb potential. The purpose of this article is to demonstrate the universality of this quantum phase transition and to present convincing experimental evidence of its existence for a charged and massless fermion in an attractive Coulomb potential as realized in graphene.

C ontinuous scale symmetry (CS)-a common property of physical systems-expresses the invariance of a physical quantity f(x) (e.g., the mass) when changing a control parameter x (e.g., the length). This property is expressed by a simple scaling relation, f(ax) = bf(x), satisfied ∀a > 0 and corresponding b(a), whose general solution is the power law f(x) = Cx α with α = ln b/ln a. Other physical systems possess the weaker discrete scale symmetry (DS) expressed by the same aforementioned scaling relation but now satisfied for fixed values (a, b) and whose solution becomes f(x) = x α G(ln x/ln a), where G(u + 1) = G(u) is a periodic function. Physical systems having a DS are also known as self-similar fractals 1 (Fig. 1a). It is possible to break CS into DS at the quantum level, a result which constitutes the basis of a special kind of scale anomaly 2,3 .
A well-studied example is provided by the problem of a particle of mass μ in an attractive inverse square potential 4,5 , which plays a role in various systems [6][7][8][9] and more importantly in Efimov physics 10,11 . Although well defined classically, the quantum mechanics of the scale-and conformal 12 -invariant Hamiltonian H = −Δ/2μ − ξ/r 2 (with ħ = 1) is well posed, but for large enough values of ξ, H is no longer self-adjoint 13,14 . The corresponding Schrödinger equation for a normalisable wave function ψ(r) of energy k 2 = −2μE is, where ζ ≡ 2μξ − l(l + d − 2) is a dimensionless parameter, d the space dimensionality and l the orbital angular momentum.
Equation (1) is invariant under the transformation r → λr and k → k/λ, ∀λ (CS), namely to every normalisable wave function of energy k 2 corresponds a continuous family of states with energies (λk) 2 , so that the bound spectrum is a continuum unbounded from below. Various ways exist to cure this problem, based on cutoff regularisation and renormalisation group [15][16][17][18][19][20][21] , and all lead for the low-energy spectrum to a quantum phase transition (QPT) monitored by ζ, between a single bound state for ζ < ζ c to an infinite and discrete energy spectrum for ζ > ζ c , independent of the regularisation procedure and given by which clearly displays DS. The critical value ζ c = (d − 2) 2 /4 depends on the space dimensionality only, and ϵ 0 is a regularization dependent energy scale. In the overcritical phase ζ > ζ c , the corresponding renormalization group solution provides a rare example of a limit cycle 15,16,22 . Building on the previous example, it can be anticipated that the problem of a massless Dirac fermion in an attractive Coulomb potential [23][24][25] , −Zα/r, is also scale invariant (CS) and that the spectrum of resonant quasibound states presents similar features and a corresponding QPT.
In this work, we demonstrate the existence of such a universal QPT for arbitrary space dimension d ≥ 2 and independently of the short distance regularisation. We obtain an explicit formula for the low-energy fractal spectrum in the overcritical regime. In contrast to the Schrödinger case equation (1), the massless Dirac Hamiltonian displays an additional parity symmetry which may be broken by the regularisation. In that case, the degeneracy of the overcritical fractal spectrum is removed and two intertwined geometric ladders of quasi-bound states appear in the s-wave channel. All these features are experimentally demonstrated using a charged vacancy in graphene. We observe the overcritical spectrum and we obtain an experimental value for the universal geometric ladder factor in full agreement with the theoretical prediction. We also explain the observation of two intertwined ladders of quasi-bound states as resulting from the breaking of parity symmetry. Finally, we relate our findings to Efimov physics as measured in cold atomic gases.

Results
The Dirac model. The Dirac equation of a massless fermion in the presence of a −Zα/r potential is obtained from the Hamiltonian (with ħ = c = 1), where (γ 0 , γ j ) are Dirac matrices. Here the dimensionless parameter monitoring the transition is β = Zα, where Z is the Coulomb charge and α the fine structure constant. The QPT occurs at the critical value β c = (d − 1)/2 (Supplementary Note 1) (A related anomalous behavior in the Dirac Coulomb problem has been identified long ago 26 but its physical relevance was marginal since it required non existent heavy-nuclei Coulomb charges Z ≃ 1/α ≃ 137 to be observed. Moreover, the problem of a massive Dirac particle is different due to the existence of a finite gap which breaks CS.). For resonant quasi-bound states, we look for scattering solutions of the form ψ in + e 2iη ψ sc , where η(E) is the energy-dependent scattering phase shift and ψ in,sc (r, E) are two component objects representing the radial part of the Dirac spinor which behave asymptotically as, for E j jr ) 1 and, using γ for E j jr ( 1 and for the lowest angular momentum channels. The two component objects V in,sc and U ± in;sc in Eqs. (4) and (5) are constants. It is easy to infer from (5) that β = β c plays a special role. Indeed for β > β c , there exists a family of normalisable solutions that admit complex eigenvalues E = −iϵ, hence the Hamiltonian (3) is not self-adjoint (H ≠ H † ). To properly define this quantum problem, a regularisation is thus needed for the too strong potential at overcritical values of β = Zα. This is achieved by introducing a cutoff length L and a boundary condition at r = L, which is equivalent to replacing the Coulomb potential at short distances by a well-behaved potential whose exact form is irrelevant in the low-energy regime EL ( 1.
The resulting mixed boundary condition can be written as where Ψ 1,2 represent the two components of the aforementioned radial part of the Dirac spinor. The resulting scattering phase shift η(E, L, h), which contains all the information about the regularisation, thus becomes a function of L and of the parameter h. The quasi-bound states energy spectrum is obtained from the scattering phase shift by means of the Krein-Schwinger relation 27,28 which relates the change of density of states δρ to the energy derivative of η, (This is also related to the Wigner time delay 29 and to the Friedel sum rule) Theoretical structure of quasi-bound spectrum. From now on, and to compare to experimental results further discussed, we consider the case d = 2, for which there is a single orbital angular momentum quantum number m 2 Z. The corresponding critical coupling becomes β c = |m + 1/2| ≥ 1/2, giving rise to the s-wave channels, m = 0, −1, for which β c = 1/2. Depending on the choice of boundary condition h, δρ(E) can be degenerate or non-degenerate over these two s-wave channels. This degeneracy originates from the symmetry of the (2 + 1) Dirac Hamiltonian (3) under parity, (x, y) → (−x, y), and its existence is equivalent to whether or not the boundary condition breaks parity (Supplementary Note 2). In what follows, we will consider the generic case in which there is no degeneracy. In the undercritical, β < β c , and low-energy regime EL ( 1, we observe (Figs. 1b, 2a) a single quasi-bound state originating from only one of the s-wave channels and which broadens as β increases. In the overcritical regime β > β c , this picture changes dramatically. (We emphasize that this picture remains valid for all values of β > β c and not only in the vicinity of β c ). The low-energy (EL ( 1) scattering phase shift displays two intertwined, infinite geometric ladders of quasi-bound states (Figs. 1c, 3) at energies E n still given by (2) but with ζ − ζ c now replaced by β 2 À β 2 c . (Moreover, note that the energy scale ϵ 0 for the Dirac case is different from the inverse square Schrödinger case defined in equation (1)). This sharp transition at β c belongs to the same universality class as presented for the inverse square Schrödinger problem, namely CS Experimental realization in graphene. A particularly interesting condensed matter system, where the previous considerations seem to be relevant is graphene in the presence of implanted Coulomb charges in conveniently created vacancies 30 . It is indeed known that low-energy excitations in graphene behave as a massless Dirac fermion field with a linear dispersion ϵ = ±v F |p| and a Fermi velocity v F ≃ 10 6 m/s 31 . These characteristics have been extensively exploited to make graphene a very useful platform to emulate specific features of quantum field theory, topology and especially QED 23 , since an effective fine structure constant α G ≡ e 2 /ħv F of order unity is obtained by replacing the velocity of light c by v F . It has been recently shown that single-atom vacancies in graphene can stably host local charge 30 . Density functional theory calculations have shown that, when a carbon atom is removed from the honeycomb lattice, the atoms around the vacancy site rearrange into a lower energy configuration 32 . The resulting lattice reconstruction causes a charge redistribution, which in the ground state has an effective local charge of ≈+1. Recent Kelvin probe force microscopy measurements of the local charge at the vacancy sites are in good agreement with the Density functional theory predictions. Vacancies are generated by sputtering graphene with He + ions 33,34 . Charge is modified and measured at the vacancy site by means of scanning tunnelling spectroscopy and Landau level spectroscopy as detailed in ref. 30 . Applying multiple pulses allows for a gradual increase in the vacancy charge, which in turn acts as an effective tunable Coulomb source. Moreover, the size of the source inside the vacancy is small (≈1 nm) as compared to the method of deposited metal clusters 35 . Using this method, we are able to observe the transition expected to occur at β = 1/2 and to measure and analyze three resonances for a broad range of β values.
To establish a relation between the measured differential conductance and the spectrum of quasi-bound states, we recall that the tunnel current I(V) is proportional to both the density of states ρ t (ϵ) of the STM tip and ρ(ϵ) of massless electronic excitations in graphene at the vacancy location. We also assume that the tunnel matrix element |t| 2 depends only weakly on energy and that both voltage and temperature are small compared to the Fermi energy and height of the tunnelling potential, so that the current I(V) = G t V is linear with V thus defining the tunnel conductance G t = 2π(e 2 /ħ)|t| 2 ρ t ρ(ϵ). Assuming that ρ t of the reference electrode (the tip) is energy independent, a variation δρ(ϵ) of the local density of states at the vacancy leads to a variation δI(V) of the current and thus to a variation δG t (V) of the tunnel conductance so that, at zero temperature, we obtain 36 where ρ 0 is the density of states in the absence of vacancy. By considering the vacancy as a local perturbation, each quasiparticle state is characterized by its scattering 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 obtained from equation (6) and combining together with equation (7) leads to the relation, between the differential tunnel conductance and the scattering phase shift. The measurements and the data analysis presented here were carried out as follows: positive charges are gradually injected into an initially prepared single atom vacancy and the differential conductance δG t (V) is measured at each step as a function of voltage. Since we are looking at the positions of resonant quasibound states, both quantities displayed in Figs. 2, 4 give the same set of resonant energies, independently of the energyindependent factor G t /πρ 0 . For low enough values of the charge, the differential conductance displayed in Fig. 2b, shows the existence of a single quasi-bound state resonance. The behavior close to the Dirac point, namely in the low-energy regime independent on the short distance regularization, is very similar to the theoretical prediction of Fig. 2a. When the build up charge exceeds a certain value, we note the appearance of three resonances, emerging out of the Dirac point. We interpret these resonances as the lowest overcritical (β > 1/2) resonances, which we denote E 1 , E ′ 1 , E 2, respectively. The corresponding theoretical and experimental behaviors displayed in Figs. 3, 4, show a very good qualitative agreement. To achieve a quantitative comparison solely based on the previous Dirac Hamiltonian equation (3), we fix L and the boundary condition h and deduce the theoretical β values corresponding to the respective positions of the lowest overcritical resonance E 1 (as demonstrated in Fig. 4). This allows to determine the lowest branch E 1 (β) for n = 1 represented in Fig. 5. Then, the experimental points E ′ 1 , E 2 are directly compared to their corresponding theoretical branch as seen in Fig. 5. We determine L and h, according to the ansatz h = a(m + 1), and obtain the best correspondence for L ≃ 0.2 nm, a ≃ −0.85. We compare the experimental E 2 /E 1 ratio with the universal prediction E nþ1 =E n ¼ e Àπ= p is fitted to the ratios E 2 /E 1 , yielding a statistical value of b = 3.145 with standard error of Δb = 0.06 consistent with the predicted value π. An error of ±1 mV is assumed for the position of the energy resonances.
A few comments are appropriate: (i) The points on the E 2 (β) curve follow very closely the theoretical prediction . This result is insensitive to the choice  of h, thus manifesting the universality of the ratio E n+1 /E n . (ii) In contrast, the correspondence between the E ′ 1 points and the theoretical branch is sensitive to the choice of h. This reflects the fact that while each geometric ladder is of the form equation (2) (with the appropriate ζ → β change), the energy scale ϵ 0 is different between the two thus leading to a shifted relative position of the two geometric ladders in Fig. 3. The ansatz taken for h is phenomenological (Supplementary Note 2), however, we find that in order to get reasonable correspondence to theory, the explicit dependence on m is needed. More importantly, it is necessary to use a degeneracy breaking boundary condition to describe the E ′ 1 (β) points. For instance, if the Coulomb potential is regularized by a constant potential for r ≤ L 37 , then both angular momentum channels (i.e., the E ′ 1 and E 1 points) become degenerate. The existence of the experimental E ′ 1 branch is therefore a distinct signal that parity symmetry in the corresponding Dirac description equation (3) is broken. In graphene, exchanging the triangular sublattices is equivalent to a parity transformation. Creating a vacancy breaks the symmetry between the two sub-lattices and is therefore at the origin of broken parity in the Dirac model. (iii) The value L ≃ 0.2 nm is fully consistent with the low-energy requirement E 1 L= hv F ' 0:03 ( 1 necessary to be in the regime relevant to observe the βdriven QPT.

Discussion
A further argument in support of the universality of this QPT is achieved by comparing the experimental results obtained in graphene with those deduced from a completely different physical problem. To that purpose, we dwell for a short while recalling the basics underlying Efimov physics 38 . Back to 1970, Efimov 10 studied the quantum problem of three identical nucleons of mass m interacting through a short range (r 0 ) potential. He pointed out that when the scattering length a of the two-body interaction becomes very large, a ) r 0 , there exists a scale-free regime for the low-energy spectrum, h 2 =ma 2 ( E ( h 2 =mr 2 0 , where the corresponding bound-states energies follow the geometric series ffiffiffiffiffiffiffiffi ffi ÀE n p ¼ Àε 0 e Àπn=s 0 À Á ; where s 0 ' 1:00624 is a dimensionless number andε 0 a problem-dependent energy scale. Efimov deduced these results from an effective Schrödinger equation in d = 3 with the radial (l = 0) attractive potential VðrÞ ¼ À s 2 0 þ 1=4 À Á =r 2 . Using Eqs. (1) and (2) and the critical value ζ c = (d − 2) 2 /4 = 1/4 for this Schrödinger problem, we deduce the ζ value for the Efimov effect to be s 2 0 þ 1=4 > ζ c corresponding to the overcritical regime of the QPT. The value of β matching to the Efimov geometric series factor e π=s0 is β E ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s 2 0 þ 1=4 p ¼ 1:1236, referred to as the fixed Efimov value. Despite being initially controversial, Efimov physics has turned into an active field especially in atomic and molecular physics where the universal spectrum has been studied experimentally [39][40][41][42][43][44][45][46] and theoretically 38 . The first two Efimov states E n (n = 1, 2) have been recently determined using an ultracold gas of caesium atoms 47 . Although the Efimov spectrum always lies at a fixed and overcritical value of the coupling, unlike the case of graphene where β can be tuned, the universal character of the overcritical regime allows nevertheless for a direct comparison of these two extremely remote physical systems. To that purpose, we include the Efimov value β E in the expression obtained for the massless Dirac fermion in a Coulomb potential and insert the corresponding data points obtained for cold atomic caesium in the graphene plot (Fig. 5) up to an appropriate scaling ofε 0 . The results are fully consistent thus showing in another way the universality presented.
There are other remote examples of systems displaying this universal QPT, e.g., flavoured QED3 48 , and the XY model (Kosterlitz-Thouless 8 and roughening transitions 22 ). Our results provide a useful and original probe of characteristic features of this universal QPT and motivate a more thorough study of this transition.

Methods
Our sample is stacked two layers of graphene on top of a thin BN flake (see Fig. 1f). The standard dry transfer procedure is followed to get this heterostructure. A large twisted angle between the two layers graphene is selected in order to weaken the coupling. The free-standing like feature for the top layer graphene is checked by the Landau levels spectroscopy. To achieve the diluted single vacancies, the sample is exposed to the helium ion beam for short time (100 eV for 5 s) followed by the high temperature annealing. The experiment is performed at 4.2 K with a home-built STM. The dI/dV (I is the current, V is the bias) is recorded by the standard lock-in technique, with a small AC modulation 2 mV at 473.1 Hz added on the DC bias. To tune the effective charge on the 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 findings of this study are available from the corresponding author upon request.