# Bielectron vortices in two-dimensional Dirac semimetals

Charles A. Downing, M. E Portnoi

► To cite this version:

Charles A. Downing, M. E Portnoi. Bielectron vortices in two-dimensional Dirac semimetals. Nature Communications, 2017, 8, pp.897. 10.1038/s41467-017-00949-y. hal-03405974

HAL Id: hal-03405974

https://hal.science/hal-03405974

Submitted on 27 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

DOI: 10.1038/s41467-017-00949-y

OPEN

Bielectron vortices in two-dimensional Dirac semimetals

C.A. Downing<fnref n="1" />,² & M.E. Portnoi ID ²,³

Searching for new states of matter and unusual quasi-particles in emerging materials and especially low-dimensional systems is one of the major trends in contemporary condensed matter physics. Dirac materials, which host quasi-particles which are described by ultra-relativistic Dirac-like equations, are of a significant current interest from both a fundamental and applied physics perspective. Here we show that a pair of two-dimensional massless Dirac-Weyl fermions can form a bound state independently of the sign of the inter-particle interaction potential, as long as this potential decays at large distances faster than Kepler's inverse distance law. This leads to the emergence of a new type of energetically favorable quasiparticle: bielectron vortices, which are double-charged and reside at zero-energy. Their bosonic nature allows for condensation and may give rise to Majorana physics without invoking a superconductor. These novel quasi-particles arguably explain a range of poorly understood experiments in gated graphene structures at low doping.

<footnote n="1">¹Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504, Strasbourg F-67000, France. ²School of Physics, University of Exeter, Stocker Road, Exeter EX4 4QL, UK. ³International Institute of Physics, Universidade Federal do Rio Grande do Norte, Natal—RN 59078-970, Brazil. C.A. Downing and M.E. Portnoi contributed equally to this work. Correspondence and requests for materials should be addressed to C.A.D. (email: downing@ipcms.unistra.fr) or to M.E.P. (email: m.e.portnoi@exeter.ac.uk)</footnote>
potential barrier. It should be emphasized that the existence of these fully confined states does not require introducing any effective mass for the quasi-particles via either imposing sublattice asymmetry, which results in a finite band gap, or considering trigonal warping terms. The only requirement for the existence of zero-energy vortices in a strong enough (beyond a critical strength) radially symmetric potential is that it has a long-distance asymptotic decay faster than Coulombic. In practice, the latter condition always takes place in realistic quasi-2D Dirac semimetals due to either screening or, for the case of graphene, the presence of a metallic gate in close proximity to the 2D electron gas (which is necessary to control the carrier density).

Fully confined zero-energy vortices should be clearly distinguished from the widely discussed 'atomic collapse' peculiarities in the graphene density of states in a supercritical attractive Coulomb potential, since the potential decaying as $1/r$ cannot support square-integrable solutions. Notably, the experimentally observed maximum in graphene's density of states in the presence of supercritical impurities$^{40}$, which is attributed to the wavefunction collapse, may be also explained using the zero-energy vortices picture in conjunction with optimal screening. Indeed, the observed peaks in the density of states are too close to the Dirac points, and the spatial extent of the measured induced charge density around the impurities is of the order of tens of the graphene lattice constant, which is easier to explain in terms of the large-size vortices rather than the short-scale wavefunction collapse at the impurity center. Furthermore, there has been a recent glut of experiments on electrostatic confinement in graphene$^{41–46}$ which may, due to the long lifetimes found, be fingerprints of zero-energy bound states.

In this article, we generalize the principles behind the aforementioned single-particle picture of confinement to the two-body problem. We show that electrostatic binding of same charge particles into bielectron vortices is both possible and energetically favorable, the effects of which will be apparent in local density of states measurements.

## Results

**Model Hamiltonian.** A consideration of two particles with an interaction potential, in the framework of a four-by-four Dirac–Weyl Hamiltonian, shows that at zero energy the sign of the potential is irrelevant for confinement. This is because the interaction potential only appears as a logarithmic derivative or as squared. Thus, forming bielectron vortices is as much a possibility as binding electrons with holes to construct excitons. The binding of repelling particles is a consequence of the symmetric gapless band structure of graphene, such that the negative kinetic energy can fully compensate electrostatic repulsion. The considered bound pairs have to be static, since two particles may only bind if they have a zero total wavevector $K$; thus we deal with 'pinned' vortex pairs. This is because for a nonzero $K$ the angular momentum $m$ is no longer a good quantum number, and necessitates one to seek a solution as a linear combination of relative motion wavefunctions with all possible values of $m$. However, this expansion includes the non-square-integrable component corresponding to $m=0$ which acts to deconfine the whole quantum state.

It is important to consider either screened systems or gated structures, which modifies the interaction from a purely Coulombic potential$^{47}$ for which no square-integrable solutions exist. The presence of metallic gates inevitably leads to image charges resulting in a fast interaction decay$^{48,49}$ at large distances and it is reasonable to introduce a cutoff at short range to avoid a Coulombic singularity. Of course, in this setup the dielectric environment is still of great importance$^{50}$, as is the geometry of

Dirac materials have low-energy fermionic excitations described by a Dirac (or Dirac–Weyl in the case of vanishing mass) Hamiltonian. This intriguing property is found in a variety of condensed matter systems, from graphene to d-wave superconductors to the surface states of topological insulators$^1$. It follows that analogs to peculiar phenomena previously studied in high-energy physics have now entered the domain of mesoscopic physics. Typical examples include Klein tunnelling$^2$, Zitterbewegung$^3$ and atomic collapse$^4, 5$. While the most recent development is a hunt for three-dimensional Weyl fermions$^6$, arguably two-dimensional (2D) Dirac semimetals$^7$, studied extensively since the exfoliation of graphene, are even more interesting due to their topological non-trivialities and, in the case of graphene, the possibility of manipulating its properties with the help of electrostatic gates. The current focus on quasi-particles in these systems is on fermionic modes, namely: Dirac, Weyl and Majorana fermions. In this article, we consider another type of quasiparticle: charged bosons formed by the pairing of two Dirac–Weyl fermions at the apex of the Dirac cone. At this point the essential chirality of the 2D Weyl quasi-particles is suppressed and quite remarkably their binding becomes possible.

Electrons in 2D Dirac semimetals can be described by the rather exotic single-particle Hamiltonian $H_1 = v_F \boldsymbol{\sigma} \cdot \mathbf{p}$, where $v_F$ is the Fermi velocity and $\boldsymbol{\sigma} = (\sigma_x, \sigma_y)$ are the spin matrices of Pauli. Recently, special attention has been paid to the apexes of the resultant conical dispersion (Dirac points) in the band structures of 2D Dirac–Weyl systems, with studies ranging from defect-induced zero-energy states$^{8–10}$, which are speculated to be a source of magnetism in graphene, to Majorana zero modes in topological insulators in proximity to superconductors$^{11}$. In this work, we predict theoretically 2D Dirac semimetals as host materials for another type of hitherto overlooked quasiparticle associated with the Dirac point: stationary (zero center-of-mass motion) bielectron vortices. Intriguingly, these bosonic quasi-particles may be a source of a new type of condensate.

Pair formation in graphene was considered in connection to an excitonic insulator, discussed well before the isolation of graphene$^{12}$ and revisited thereafter$^{13–15}$. There has been significant interest in spatially-separated double-layer graphene$^{16–20}$, where several groups have studied Bose–Einstein condensation$^{21–23}$ and superfluidity$^{24, 25}$ in this gapped system. However, no gap has as yet been observed experimentally in monolayer graphene structures$^{26}$ and a question remains: can two charge carriers bind together in an ideal 2D Dirac–Weyl system? It has previously been claimed that excitons do not exist in gapless graphene$^{27}$, and so considerations of trigonal warping, which effectively introduces an angular-dependent single-particle mass, have been suggested as a route toward pair formation$^{28–30}$.

It is a commonly held belief that electrostatic confinement of 2D massless Dirac fermions is impossible as a result of the Klein paradox, where there is a perfect transmission for normally incident particles—this prompted proposals for localization via a variety of other means$^{31–35}$. An argument is usually made that conservation of pseudospin $\boldsymbol{\sigma} \cdot \mathbf{p} = \pm 1$ forbids bound states in purely electrostatic problems$^2$. However, at the Dirac point pseudospin is not well defined, a fact we exploit in this work. Indeed, unlike the case of finite energy, zero-energy bound states may form at the apex of the Dirac cone$^{36–38}$. Mathematically, this is because at finite energy the effective Schrödinger equation at long-range maps on to the problem of scattering states in a non-relativistic system$^{39}$; whereas at zero-energy solutions exist which decay algebraically, depending on the angular momentum quantum number $m$. When $m$ is nonzero the solutions are fully square-integrable, such that they are rotating ring-like states (vortices) avoiding the Klein tunneling due to their vorticity, which results in a nonzero momentum component along the
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Radial probability densities of bielectron vortices. A plot of the radial probability densities for the first two bielectron states, with the quantum numbers **a** (*m*, *n*) = (1, 0) and **b** (*m*, *n*) = (2, 0), as a function of position in two-dimensions. The spatial coordinates are measured in units of the length scale *d*. The color bar measures the dimensionless number associated with the probability density. The interaction strengths are given by Eq. (6)</figcaption></figure>

the device, which both contribute to the effective strength of the interaction. As we demonstrate below, the seemingly rigid conditions on the strength and spatial extent of the inter-particle potential, required to maintain the total energy at zero, are in fact easily satisfied for large-size vortices by linear screening provided by a small number of residual free carriers.

Previous theoretical works on excitonic effects in Dirac materials have approached the problem via either exact diagonalization⁵¹, the Bethe-Salpeter formalism¹⁵, ⁵² or in the language of a two-body matrix Hamiltonian⁵³⁻⁵⁵, which we will utilize here. The two-body Hamiltonian can be written as the Kronecker sum of the single-particle Hamiltonians $H = H_1 \oplus H_2$, or explicitly (as there are two sublattices and two particles) as the 4 × 4 matrix

$$H = v_F \begin{bmatrix} 0 & p_{x2} - ip_{y2} & p_{x1} - ip_{y1} & 0 \\ p_{x2} + ip_{y2} & 0 & 0 & p_{x1} - ip_{y1} \\ p_{x1} + ip_{y1} & 0 & 0 & p_{x2} - ip_{y2} \\ 0 & p_{x1} + ip_{y1} & p_{x2} + ip_{y2} & 0 \end{bmatrix}, \quad (1)$$

where the subscripts 1 and 2 refer to the two particles. The matrix Hamiltonian given by Eq. (1) is written for two electrons belonging to the same Dirac valley. It can be modified for the particles of different charge (electron and hole) and for two particles belonging to different valleys. Here and in what follows we also neglect spin, which is in principle important as it governs the parity of the relative motion function for the same-valley electrons. However, our immediate aim is to demonstrate the existence of bound states leaving classification of all possible pairs for a future work.

We expect our two-particle continuum theory to be a good approximation to the experimental reality, since in the single-particle picture the theory of zero-energy states³⁶⁻³⁸ has successfully predicted confinement effects seen in some recent experiments⁴¹, ⁴⁴, ⁴⁶. In addition, these toy model results from Dirac equations have been shown to be robust to sophisticated numerical experiments on finite sized flakes⁵⁶.

**Bielectronic solutions of the model.** The Hamiltonian (1) acts upon a two-particle wavefunction constructed via the Kronecker product $\Psi(\mathbf{r}_1, \mathbf{r}_2) = \psi_i(\mathbf{r}_1) \otimes \psi_j(\mathbf{r}_2)$, where $i, j = (A, B)$. In the absence of an interaction potential $U(\mathbf{r}_1 - \mathbf{r}_2)$, diagonalization of Eq. (1) yields four eigenenergies: $E = \pm v_F (p_{x1}^2 + p_{y1}^2)^{1/2} \pm v_F (p_{x2}^2 + p_{y2}^2)^{1/2}$. As is usual with two-body problems, we utilize the center-of-mass and relative motion coordinates:

$X = (x_1 + x_2)/2$, $Y = (y_1 + y_2)/2$, $x = x_1 - x_2$, $y = y_1 - y_2$. Upon assuming a translationally invariant system, such that the center-of-mass momentum $\hbar\mathbf{K}$ is a constant of motion, one can employ the ansatz $\Psi_i(\mathbf{R}, \mathbf{r}) = \exp(i\mathbf{K} \cdot \mathbf{R})\psi_i(\mathbf{r})$, where the index $i = (1, 2, 3, 4)$ numerates the four components of the wavefunction, which span the two sublattices and two particles. As shown in refs. ⁵³⁻⁵⁵ when $\mathbf{K} = 0$ one can rewrite the relative motion Cartesian coordinates $(x, y)$ in polar coordinates $(r, \theta)$, eventually reducing Eq. (1) to a system of three equations only for the transformed radial wavefunction components $\phi_i(r)$,

$$\begin{bmatrix} \frac{U(r)-E}{\hbar v_F} & \partial_r + \frac{m}{r} & 0 \\ 2(-\partial_r + \frac{m-1}{r}) & \frac{U(r)-E}{\hbar v_F} & 2(\partial_r + \frac{m+1}{r}) \\ 0 & -\partial_r + \frac{m}{r} & \frac{U(r)-E}{\hbar v_F} \end{bmatrix} \begin{bmatrix} \phi_1(r) \\ \phi_2(r) \\ \phi_3(r) \end{bmatrix} = 0, \quad (2)$$

with $m = 0, \pm 1, \pm 2, \dots$ and where one can take $\phi_4 = 0$.

Let us now consider a model interaction given by $U(r) = U_0/(1 + (r/d)^2)$, with an on-site energy $U_0$ and the long-range cutoff parameter *d*, which may be related to the separation between the 2D semimetal and the back-gate or to the screening length⁴⁸, ⁴⁹. This model potential provides a reasonable approximation to the more realistic potential decaying at large distances as $1/r^3$, for details see the Supplemental Information. Notably this functional form is well known in optics as the spatially inhomogeneous Maxwell's fish-eye lens⁵⁷, and remarkably is the simplest exactly solvable model, as the square well does not admit a nontrivial solution.

The system of equations (2) can be reduced to a second order differential equation for $\phi_2$ only, which admits an analytical solution for the chosen interaction. This solution is square-integrable only at the Dirac point ($E=0$). The same is true for any potential decaying faster than the Coulomb potential, so from now on we consider zero-energy states only. Now, when $r \sim 0$, one finds the usual short-range behavior $\phi_2 \sim r^{|m|}$. Meanwhile the asymptotic behavior as $r \to \infty$ is given by the decay $\phi_2 \sim r^{|m|-2\eta_m}$, where $\eta_m = (|m| + 1 + \sqrt{m^2 + 1})/2$. Thus, we are motivated to seek a solution of Eq. (2) with the ansatz

$$\phi_2(r) = \frac{(r/d)^{|m|}}{(1 + (r/d)^2)^{\eta_m}} f(r), \quad (3)$$

where $f(r)$ is a polynomial in $r$ that does not affect the short- and long-range behavior. Upon substituting Eq. (3) into Eq. (2), eliminating $\phi_{1,3}(r)$, and using the new variable $\xi = (r/d)^2$, we
arrive at the following equation for $f(\xi)$:

$$ \xi(1+\xi)^2 f''(\xi) + (1+\xi)[m+1+(m+2-2\eta_m)\xi]f'(\xi) + \left[ \left(\frac{1}{4}\frac{U_0d}{\hbar v_F}\right)^2 - \eta_m^2 \right] f(\xi) = 0, \quad (4) $$

which is a form of the Gauss hypergeometric equation⁵⁸. Its solution, regular at $\xi=0$, is given by

$$ f(\xi) = {}_2F_1\left(-n, -n + \frac{|U_0|d}{2\hbar v_F}; |m| + 1; \frac{\xi}{1+\xi}\right), \quad (5) $$

where we have terminated the power series in the Gauss hypergeometric function ${}_2F_1(a, b; c; x)$ to ensure decaying solutions at infinity. This termination leads to the following quantization condition for the formation of bound bielectron pairs

$$ \frac{|U_0|d}{\hbar v_F} = 4(n + \eta_m), \quad n = 0, 1, 2, \dots \quad (6) $$

The other wavefunction components $\phi_{1,3}(r)$ are readily obtainable from Eq. (2), and their long-range behavior $r \to \infty$ tells us that the $m=0$ state is non-square-integrable, since $(\phi_1, \phi_2, \phi_3) \to r^{-\sqrt{1+|m|^2}}(1, r^{-1}, 1)$. Thus, the pair states are rotating ring-like modes (vortices). Probability density plots are displayed in Fig. 1 for lowest node ($n=0$) states with $m=1, 2$. Most noticeable from the figure is the characteristic vortex-like shape of the bielectron states.

Notably, Eq. (6) displays two regimes of interest. In the subcritical case, the threshold value for the first confined state to appear is not met, $|U_0|d/\hbar v_F < \alpha_c$. The critical strength $\alpha_c \approx 6.83$ is found from Eq. (6) with $n=0, 1, 2, \dots, n=\infty$ ($|m|=1$). However, in the opposing (supercritical) domain $|U_0|d/\hbar v_F \ge \alpha_c$, and pairs may indeed form. Weak screening by a small number of mobile uncoupled carriers allows the system to adjust the inter-particle interaction potential so that it satisfies the strength condition given by Eq. (6) to support bound states, resulting in an energetically favorable drastic reduction in the chemical potential of the many-electron system, accompanied by a narrow spike in the density of states at zero energy. Indeed, pair formation due to doping is a well-known mechanism⁵⁹.

**Relation to experiments.** A more accurate treatment of the interaction potential, taking into account both an image charge necessarily present in gated structures and a regularization of the interaction as $r \to 0$, may be tackled numerically (see the Supplemental Information for details). The main difference is that the realistic potential falls at large distances as $1/r^3$, which is faster than the exactly solvable Lorentzian potential. As a result, the critical strength required for binding two electrons depends mostly on the Dirac semimetal fine structure constant (dimensionless interaction strength) $\alpha = e^2/(\kappa\hbar v_F)$ adjusted by a numerical factor of the order of unity, which depends on the ratio of the short-range potential cutoff and its long-range scale. Typical values are $\alpha \approx 2.19/\kappa, 4.38/\kappa$ for graphene¹ or surface states of three-dimensional topological insulators⁶⁰, respectively, where $\kappa$ is the relative permittivity of the material. According to ref.⁶¹, for gapless versions of silicene and germanene $\alpha \approx 4.06/\kappa, 4.13/\kappa$, respectively. Our numerical estimates (see the Supplemental Information for details) show that the parameter $\alpha_c$ required for forming the first ($|m|=1$) bielectron vortices is $\alpha_c \approx 2.5$. This condition is not satisfied for the case of suspended graphene ($\kappa=1$). However, the discrepancy is not very large and can be compensated by the moderate decrease of the Fermi velocity $v_F$ due to local stretching. Namely, in graphene the local expansion of the honeycomb lattice acts to decrease the

Fermi velocity and thus the effective potential strength may indeed enter the supercritical regime even for this system. In fact, strain-induced corrugations in real graphene samples have been shown to give rise to well-defined regions of electron-hole puddles⁶², ⁶³. Furthermore, the inclusion of static screening⁶⁴⁻⁶⁶ alone gives the interaction strength for suspended graphene almost sufficient for observing the vortices (see the Supplemental Information for details), and a small additional stretching will help their formation.

Most of the other gapless 2D Dirac systems¹ have Fermi velocities significantly smaller than that of graphene, so the critical strength condition can be easily satisfied. This suggests bielectron vortices should be present at moderate carrier densities in topological insulators, single-valley gapless mercury telluride quantum wells and silicene⁶⁷. Where the Fermi velocity cannot be locally adjusted by stretching, there should be a local pinning of the Fermi level in order to provide optimal screening which maintains the critical interaction strength until the vortices start overlapping (in analogy to a Mott transition). With further carrier density increase, screening effects will lead to the eventual disappearance of vortices when the long-range scale of the potential diminishes. Observing Fermi level pinning with moderate changes of carrier density in low-density ‘rigid’ Dirac–Weyl systems will be the most unambiguous proof of the existence of bielectron vortices.

One may speculate that fingerprints of bielectron vortices have already been observed in the range of experiments on graphene. Indeed, a reservoir of stationary, zero-energy bielectron vortices may offer a contribution to the experimentally seen Fermi velocity renormalization in gated graphene structures⁶⁸⁻⁷¹, which is observed instead of the widely theorized gap. According to this picture, the observed Fermi velocity renormalization could be an artifact of overestimating the number of charge carriers defining the position of the Fermi level; since a large number of them disappear into a many-body ground state of bosonic vortices. The best-known experiment⁶⁸ on Fermi velocity renormalization ($v_F \to v_F^*$) is based on measuring the cyclotron mass, given by $m_c = \hbar(\pi n)^{1/2}/v_F$, where $n$ is the carrier density. However, if a large amount of the carriers condense into a reservoir of zero-energy bosonic vortices, the corrected lower density $n \to n^*$ of remaining free fermions should be substituted into the cyclotron mass formula, then the smaller observed cyclotron mass may be explained without the need of renormalizing $v_F \to v_F^*$. The same is true for the quantum capacitance measurements⁷⁰, since the presence of the charged boson reservoir changes drastically the Fermi energy dependence on the total carrier density from the expected relation, which is used to estimate the renormalized $v_F^*$. Notably, both the original theory of Fermi velocity renormalization in free-standing graphene⁷² and its later refinement⁷³ are based on the long-range behavior of the unscreened Coulomb potential resulting in logarithmically divergent corrections at small $n$. Therefore, we believe that the applicability of these theoretical results should be taken with caution for screened and/or gated structures, in particular graphene on graphite⁷¹.

## Discussion

In conclusion, we have demonstrated that a gapless Dirac–Weyl 2D system with a short-range inter-particle interaction favors the existence of zero-energy charged bound pairs. The associated peak in the local density of states at the Dirac point, which is sensitive to the carrier density, should be taken into account for the interpretation of the scanning tunneling microscopy results. This peak could also serve as a source of carriers with energies corresponding to the strong nonlinear electromagnetic response⁷⁴ making low-doped graphene better suited for relevant applications⁷⁵.
Arguably, the reservoir of bosonic vortices can play a similar role to that of a superconductor in proximity to a Weyl semimetal by enforcing electron–hole symmetry. Indeed, adding an electron to the considered system is equivalent to adding a hole and another zero-energy vortex which makes this system a promising candidate in the on-going search of Majorana modes in solids. Notably, the particle–hole symmetry provided by the condensate is a necessary, rather than a sufficient, condition for creating Majorana modes. Searching for the most suitable Dirac–Weyl system, which will involve an appropriate Chern number analysis, is one of the avenues for future work.

The observed puddles of charged carriers in graphene in the case of long-range disorder⁷⁶ can be treated as many-body mesoscopic domains containing condensates of bosonic bipartite vortices, thus removing the controversy of having carrier puddles despite the absence of single-particle localization in smooth potentials due to the Klein phenomenon. Investigations of this new and unconventional many-body state, with special regard to possible occurrences of quantum critical phase transitions, will form part of a future work. The effect of puddles on the system is controlled by the tiny balance between the electrostatic energy from the positively and negatively charged droplets and the energy of the separating domain walls. The problem is therefore similar to the formation of Landau–Kittel domain structures in ferroelectric materials⁷⁷, which was solved recently for domains of an arbitrary shape⁷⁸. The phase diagram of the system with puddle decomposition could even be similar to that discovered in strained dioxyde vanadium, VO₂ nanoplatelets with metallic and insulator domain separation, controlled by long-range elastic forces⁷⁹.

## Methods

In this theoretical paper, all methods used are fully described in the Results section.

**Data availability.** The authors declare that all of the data supporting the findings of this theoretical study are available within the paper and its supplementary information files.

Received: 3 May 2017 Accepted: 8 August 2017

Published online: 12 October 2017

## References

1. Wehling, T. O., Black-Schaffer, A. M. & Balatsky, A. V. Dirac materials. *Adv. Phys.* **63**, 1–76 (2014).

2. Katsnelson, M. I., Novoselov, K. S. & Geim, A. K. Chiral tunnelling and the Klein paradox in graphene. *Nat. Phys.* **2**, 620–625 (2006).

3. Katsnelson, M. I. Zitterbewegung, chirality, and minimal conductivity in graphene. *Eur. Phys. J. B* **51**, 157–160 (2006).

4. Shytov, A. V., Katsnelson, M. I. & Levitov, L. S. Atomic collapse and quasi–Rydberg states in graphene. *Phys. Rev. Lett.* **99**, 246802 (2007).

5. Downing, C. A. & Portnoi, M. E. One-dimensional Coulomb problem in Dirac materials. *Phys. Rev. A* **90**, 052116–1–052116–5 (2014).

6. Hills, D. Y., Kusmartseva, A. & Kusmartsev, F. V. Current-voltage characteristics of Weyl semimetal semiconducting devices, Veselago lenses, and hyperbolic Dirac phase. *Phys. Rev. B* **95**, 214103 (2017).

7. Young, S. M. & Kane, C. L. Dirac semimetals in two dimensions. *Phys. Rev. Lett.* **115**, 126803 (2015).

8. Pereira, V. M., Guinea, F., Lopes dos Santos, J. M. B., Peres, N. M. R. & Castro Neto, A. H. Disorder induced localized states in graphene. *Phys. Rev. Lett.* **96**, 036801 (2006).

9. Yazyev, O. V. & Helm, L. Defect-induced magnetism in graphene. *Phys. Rev. B* **75**, 125408 (2007).

10. Espinosa-Ortega, T., Luk'yanchuk, I. A. & Rubo, Y. G. Magnetic properties of graphene quantum dots. *Phys. Rev. B* **87**, 205434 (2013).

11. Rokhinson, L. P., Liu, X. & Furdyna, J. K. The fractional ac Josephson effect in a semiconductor–superconductor nanowire as a signature of Majorana particles. *Nat. Phys.* **8**, 795–799 (2012).

12. Khveshchenko, D. V. Ghost excitonic insulator transition in layered graphite. *Phys. Rev. Lett.* **87**, 246802 (2011).

13. Drut, J. E. & Lahde, T. A. Is graphene in vacuum an insulator? *Phys. Rev. Lett.* **102**, 026802 (2009).

14. Stroucken, T., Gronqvist, J. H. & Koch, S. W. Optical response and ground state of graphene. *Phys. Rev. B* **84**, 205445 (2011).

15. Gamayun, O. V., Gorbar, E. V. & Gusynin, V. P. Supercritical Coulomb center and excitonic instability in graphene. *Phys. Rev. B* **80**, 165429 (2009).

16. Yu, Kharitonov, M. & Efetov, K. B. Electron screening and excitonic condensation in double-layer graphene systems. *Phys. Rev. B* **78**, 241401 (2008).

17. Min, H., Bistritzer, R., Su, J.-J. & MacDonald, A. H. Room-temperature superfluidity in graphene bilayers. *Phys. Rev. B* **78**, 121401 (2008).

18. Zhang, C.-H. & Joglekar, Y. N. Excitonic condensation of massless fermions in graphene bilayers. *Phys. Rev. B* **77**, 233405 (2008).

19. Sodemann, I., Pesin, D. A. & MacDonald, A. H. Interaction-enhanced coherence between two-dimensional Dirac layers. *Phys. Rev. B* **85**, 195136 (2012).

20. Abergel, D. S. L., Rodriguez-Vega, M., Rossi, E. & Das Sarma, S. Interlayer excitonic superfluidity in graphene. *Phys. Rev. B* **88**, 235402 (2013).

21. Lozovik, Yu. E. & Sokolik, A. A. Electron-hole pair condensation in a graphene bilayer. *JETP Lett.* **87**, 55–59 (2008).

22. Berman, O. L., Lozovik, Y. E. & Gumbs, G. Bose-Einstein condensation and superfluidity of magnetoexcitons in bilayer graphene. *Phys. Rev. B* **77**, 155433 (2008).

23. Berman, O. L., Kezerashvili, R. Y. & Ziegler, K. Coupling of two Dirac particles. *Phys. Rev. A* **87**, 042513 (2013).

24. Berman, O. L., Kezerashvili, R. Y. & Ziegler, K. Superfluidity of dipole excitons in the presence of band gaps in two-layer graphene. *Phys. Rev. B* **85**, 0035418 (2012).

25. Berman, O. L., Kezerashvili, R. Y. & Ziegler, K. Superfluidity and collective properties of excitonic polaritons in gapped graphene in a microcavity. *Phys. Rev. B* **86**, 235404 (2012).

26. Mayorov, A. S. et al. How close can one approach the Dirac point in graphene experimentally? *Nano Lett.* **12**, 4629–4634 (2012).

27. Ratnikov, P. V. & Silin, A. P. Size quantization in planar graphene-based heterostructures: pseudospin splitting, interface states, and excitons. *JETP Lett.* **114**, 512–528 (2012).

28. Mahmoodian, M. M. & Entin, M. V. Moving zero-gap Wannier-Mott excitons in graphene. *Europhys. Lett.* **102**, 37012 (2013).

29. Marnham, L. L. & Shytov, A. V. Metastable electron-electron states in double-layer graphene structures. *Phys. Rev. B* **92**, 085409 (2015).

30. Marnham, L. L. & Shytov, A. V. Bielectrons in the Dirac sea in graphene: the role of many-body effects, preprint at https://arxiv.org/abs/1512.02953 (2015).

31. Rozhkov, A. V., Giavaras, G., Bliokh, Y. P., Freilikher, V. & Nori, F. Electronic properties of mesoscopic graphene structures: charge confinement and control of spin and charge transport. *Phys. Rep.* **503**, 77–114 (2011).

32. Giavaras, G., Maksym, P. A. & Roy, M. Magnetic field induced confinement–deconfinement transition in graphene quantum dots. *J. Phys. Condens. Matter* **21**, 102201 (2009).

33. Downing, C. A. & Portnoi, M. E. Massless Dirac fermions in two dimensions: confinement in nonuniform magnetic fields. *Phys. Rev. B* **94**, 165407 (2016).

34. Raoux, A. et al. Velocity-modulation control of electron-wave propagation in graphene. *Phys. Rev. B* **81**, 073407 (2010).

35. Downing, C. A. & Portnoi, M. E. Localization of massless Dirac particles via spatial modulations of the Fermi velocity. *J. Phys. Condens. Matter.* **29**, 315301 (2017).

36. Bardarson, J. H., Titov, M. & Brouwer, P. W. Electrostatic confinement of electrons in an integrable graphene quantum dot. *Phys. Rev. Lett.* **102**, 226803 (2009).

37. Downing, C. A., Stone, D. A. & Portnoi, M. E. Zero-energy states in graphene quantum dots and rings. *Phys. Rev. B.* **84**, 155437 (2011).

38. Downing, C. A., Pearce, A. R., Churchill, R. J. & Portnoi, M. E. Optimal traps in graphene. *Phys. Rev. B.* **92**, 165401 (2015).

39. Tudorovskiy, T. Y. & Chaplik, A. V. Spatially inhomogeneous states of charge carriers in graphene. *JETP Lett.* **84**, 619–623 (2007).

40. Wang, Y. et al. Observing atomic collapse resonances in artificial nuclei on graphene. *Science* **340**, 734–737 (2013).

41. Margapoti, E. et al. Emergence of photoswitchable states in a graphene–azobenzene–Au platform. *Nano Lett.* **14**, 6823–6827 (2014).

42. Zhao, Y. et al. Creating and probing electron whispering-gallery modes in graphene. *Science* **348**, 672–675 (2015).

43. Mao, J. et al. Realization of a tunable artificial atom at a supercritically charged vacancy in graphene. *Nat. Phys.* **12**, 545–549 (2016).

44. Gutierrez, C., Brown, L., Kim, C.-J., Park, J. & Pasupathy, A. N. Klein tunnelling and electron trapping in nanometre-scale graphene quantum dots. *Nat. Phys.* **12**, 1069–1075 (2016).

45. Lee, J. et al. Imaging electrostatically confined Dirac fermions in graphene quantum dots. *Nat. Phys.* **12**, 1032–1036 (2016).
46. Bai, K.-K. et al. Generating nanoscale and atomically-sharp p-n junctions in graphene via monolayer-vacancy-island engineering of Cu surface, preprint at https://arxiv.org/abs/1705.10952 (2017).

47. Shytov, A. V., Katsnelson, M. I. & Levitov, L. S. Vacuum polarization and screening of supercritical impurities in graphene. *Phys. Rev. Lett.* **99**, 236801 (2007).

48. Fogler, M. M., Novikov, D. S. & Shklovskii, B. I. Screening of a hypercritical charge in graphene. *Phys. Rev. B* **76**, 233402 (2007).

49. Asgari, R., Katsnelson, M. I. & Polini, M. Quantum capacitance and Landau parameters of massless Dirac fermions in graphene. *Ann. Phys. (Berlin)* **526**, 359–365 (2014).

50. Perebeinos, V., Tersoff, J. & Avouris, P. Scaling of excitons in carbon nanotubes. *Phys. Rev. Lett.* **92**, 257402 (2004).

51. Paananen, T. & Egger, R. Finite-size version of the excitonic instability in graphene quantum dots. *Phys. Rev. B* **84**, 155456 (2011).

52. Wang, J., Fertig, H. A., Murthy, G. & Brey, L. Excitonic effects in two-dimensional massless Dirac fermions. *Phys. Rev. B* **83**, 035404-1-035404-15 (2011).

53. Sabio, J., Sols, F. & Guinea, F. Two-body problem in graphene. *Phys. Rev. B* **81**, 045428-1-045428-12 (2010).

54. Hartmann, R. R., Shelykh, I. A. & Portnoi, M. E. Excitons in narrow-gap carbon nanotubes. *Phys. Rev. B* **84**, 035437 (2011).

55. Lee, R. N., Milstein, A. I. & Terekhov, I. S. Quasilocalized states in a model of electron-electron interaction in graphene. *Phys. Rev. B* **86**, 035425 (2012).

56. Pieper, A., Heinisch, R. L., Wellein, G. & Fehske, H. Dot-bound and dispersive states in graphene quantum dot superlattices. *Phys. Rev. B* **89**, 165121 (2014).

57. Born, M. & Wolf, E. Principles of Optics (Cambridge University Press, 1999).

58. Gradshteyn, I. S. & Ryzhik, I. M. Table of Integrals, Series and Products (Academic, 1980).

59. Babichenko, V. S. & Kiselev, M. N. Doping-induced exciton transition. *JETP Lett.* **57**, 179–184 (1993).

60. Zhang, H. et al. Topological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surface. *Nat. Phys.* **5**, 438–442 (2009).

61. Yan, J.-A., Stein, R., Schaefer, D. M., Wang, X.-Q. & Chou, M. Y. Electron-phonon coupling in two-dimensional silicene and germanene. *Phys. Rev. B* **88**, 121403 (2013).

62. Gibertini, M., Tomadin, A., Guinea, F., Katsnelson, M. I. & Polini, M. Electron-hole puddles in the absence of charged impurities. *Phys. Rev. B* **85**, 201405 (2012).

63. Couto, N. J. G. et al. Random strain fluctuations as dominant disorder source for high-quality on-substrate graphene devices. *Phys. Rev. X* **4**, 041019 (2014).

64. Stern, F. & Howard, W. E. Properties of semiconductor surface inversion layers in the electric quantum limit. *Phys. Rev.* **163**, 816–835 (1967).

65. Portnoi, M. E. & Galbraith, I. Variable-phase method and Levinson's theorem in two dimensions: application to a screened Coulomb potential. *Solid State Commun.* **103**, 325–329 (1997).

66. Das Sarma, S., Adam, S., Hwang, E. H. & Rossi, E. Electronic transport in two-dimensional graphene. *Rev. Mod. Phys.* **83**, 407–470 (2011).

67. Drummond, N. D., Zolyomi, V. & Fal'ko, V. I. Electrically tunable band gap in silicene. *Phys. Rev. B* **85**, 075423 (2012).

68. Elias, D. C. et al. Dirac cones reshaped by interaction effects in suspended graphene. *Nat. Phys.* **7**, 701–704 (2011).

69. Chae, J. et al. Renormalization of the graphene dispersion velocity determined from scanning tunneling spectroscopy. *Phys. Rev. Lett.* **109**, 116802 (2012).

70. Yu, G. L. et al. Interaction phenomena in graphene seen through quantum capacitance. *Proc. Natl Acad. Sci. USA* **110**, 3282–3286 (2013).

71. Faugeras, C. et al. Landau level spectroscopy of electron-electron interactions in graphene. *Phys. Rev. Lett.* **114**, 126804 (2015).

72. Gonzalez, J., Guinea, F. & Vozmediano, M. A. H. Marginal-Fermi-liquid behavior from two-dimensional Coulomb interaction. *Phys. Rev. B* **59**, 2474–2477 (1999).

73. Hofmann, J., Barnes, E. & Das Sarma, S. Why does graphene behave as a weakly interacting system? *Phys. Rev. Lett.* **113**, 105502 (2014).

74. Mikhailov, S. A. & Ziegler, K. Nonlinear electromagnetic response of graphene: frequency multiplication and the self-consistent-field effects. *J. Phys. Condens. Matter* **20**, 384204 (2008).

75. Khurgin, J. B. Graphene: a rather ordinary nonlinear optical material. *Appl. Phys. Lett.* **104**, 161116 (2014).

76. Martin, J. et al. Observation of electron-hole puddles in graphene using a scanning single-electron transistor. *Nat. Phys.* **4**, 144–148 (2008).

77. Landau, L. D. & Lifshitz, E. M. Electrodynamics of Continuous Media (Elsevier, 1985).

78. Lukyanchuk, I. et al. High-symmetry polarization domains in low-symmetry ferroelectrics. *Nano Lett.* **14**, 6931–6935 (2014).

79. Tselev, A., Strelcov, E. & Lukyanchuk, I. A. et al. Interplay between ferroelastic and metal-insulator phase transitions in strained quasi-two-dimensional VO2 nanoplatelets. *Nano Lett.* **10**, 2003–2011 (2010).

## Acknowledgements

C.A.D. recognizes financial support from the EPSRC DTP (Award reference 1080089).

We also acknowledge support from the CNRS, the EU H2020 RISE project CoExAN (Grant No. H2020-644076), EU FP7 ITN NOTEDEV (Grant No. FP7-607521),

and the FP7 IRSES projects CANTOR (Grant No. FP7-612285), QOCaN (Grant No.

FP7-316432) and InterNoM (Grant No. FP7-612624). We thank E. Mariani, R.J.

Nicholas, L.A. Ponomarenko and B.I. Shklovskii for fruitful discussions.

## Additional information

Supplementary Information accompanies this paper at doi:10.1038/s41467-017-00949-y.

**Competing interests:** The authors declare no competing financial interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017