Zero-energy vortices in Dirac materials

TL;DR

This paper reviews methods for confining massless Dirac particles in scalar potentials, presenting exact solutions and implications for experiments in two-dimensional Dirac materials, including a two-body Dirac-Kepler problem analysis.

Contribution

It introduces a route to achieve discrete zero-energy states in Dirac particles and provides an analytic solution to the two-body Dirac-Kepler problem, advancing understanding of confinement and atomic collapse.

Findings

Discrete zero-energy states can be achieved in scalar potentials.

Analytic solution to the two-body Dirac-Kepler problem.

Insights into atomic collapse in mesoscopic Dirac systems.

Abstract

In this brief review, we survey the problem of electrostatic confinement of massless Dirac particles, via a number of exactly solvable one- and two-body models. By considering bound states at zero energy, we present a route to obtain truly discrete states of massless Dirac particles in scalar potentials, circumventing the celebrated Klein tunnelling phenomenon. We also show how the coupling of two ultrarelativistic particles can arise, and discuss its implications for cutting-edge experiments with two-dimensional Dirac materials. Finally, we report an analytic solution of the two-body Dirac-Kepler problem, which may be envisaged to bring about a deeper understanding of critical charge and atomic collapse in mesoscopic Dirac systems.

Figures

Equations

$$H=v_{\mathrm{F}}\mathbf{\sigma }\cdot \hat{\mathbf{p}}+V(r),$$

$$\Psi(r,\theta)=\frac{e^{\mathrm{i}m\theta}}{\sqrt{2\pi}}\left(\begin{array}[]{c}\chi_{A}(r)\\ \mathrm{i}e^{\mathrm{i}\theta}\chi_{B}(r)\end{array}\right),$$

$$\left({\partial }_r+\frac{m+1}{r}\right){\chi }_B$$

$$\left(-{\partial }_r+\frac{m}{r}\right){\chi }_A$$

$${\chi }_A^{\prime \prime }+\left(\frac{1}{r}+\frac{U^{\prime }}{\epsilon -U}\right){\chi }_A^{\prime }+\left({\left[\epsilon -U\right]}^2-\frac{m^2}{r^2}-\frac{m}{r}\frac{U^{\prime }}{\epsilon -U}\right){\chi }_A=0,$$

$${\chi }_A^{\prime \prime }+\left(\frac{1}{r}\right){\chi }_A^{\prime }+\left({\epsilon }^2-\frac{m^2}{r^2}\right){\chi }_A=0,$$

$${\chi }_A=c_1{J}_m(\epsilon r)+c_2{Y}_m(\epsilon r),$$

$$U(r)=-\frac{\alpha }{r},$$

$${\chi }_{A,B}=r^{{\gamma }_m-1/2}{e}^{-\mathrm{i}\epsilon r}{f}_{A,B}(\epsilon r),$$

$${\gamma }_m=\sqrt{{{\alpha }_m^{\mathrm{c}}}^2-{\alpha }^2},$$

$${\alpha }_m^{\mathrm{c}}=|m+1/2|,$$

$$U(r)=-U_0\Theta (a-r),$$

$$\left(\begin{array}[]{c}\chi_{A}^{\mathrm{in}}\\ \chi_{B}^{\mathrm{in}}\end{array}\right)=\frac{C_{nm}}{a}\left(\begin{array}[]{c}J_{m}(U_{0}r)\\ J_{m+1}(U_{0}r)\end{array}\right),$$

$$\left(\begin{array}[]{c}\chi_{A}^{\mathrm{out}}\\ \chi_{B}^{\mathrm{out}}\end{array}\right)=\frac{C_{nm}}{a}\left(\begin{array}[]{c}0\\ J_{m+1}(U_{0}a)\left(\frac{a}{r}\right)^{m+1}\end{array}\right),m=0,1,...$$

$$U_{0}a={\alpha }_{m,n},$$

$$\left(\begin{array}[]{c}\chi_{A}^{\mathrm{out}}\\ \chi_{B}^{\mathrm{out}}\end{array}\right)=\frac{C_{nm}}{a}\left(\begin{array}[]{c}J_{m}(U_{0}a)\left(\frac{r}{a}\right)^{m}\\ 0\end{array}\right),m=-1,-2,...$$

$$U(r)=\frac{-U_0}{1+r^2/d^2},$$

$$p_m=\frac{1}{2}\left(1+|m|+|1+m|\right).$$

$${\chi }_A=\frac{{\xi }^{\frac{|m|}{2}}}{(1+\xi {)}^{p_m}}w(\xi ),$$

$$\xi (1+\xi {)}^2{w}^{\prime \prime }(\xi )$$

$$+(1+\xi )\left[1+|m|+(2+|m|-2p_m)\xi \right]w^{\prime }(\xi )$$

$$+\left[(\frac{U_{0}d}{2}{)}^2-p_m^2\right]w(\xi )=0,$$

$$w(\xi )={}_2{F}_1\left(p_m+\frac{1}{2}{U}_{0}d,p_m-\frac{1}{2}{U}_{0}d;1+|m|;\frac{\xi }{\xi +1}\right).$$

$$\lim_{r\to\infty}\left(\begin{array}[]{c}\chi_{A}\\ \chi_{B}\end{array}\right)\propto\frac{1}{r}\left(\begin{array}[]{c}\frac{1}{r^{|}1+m|}\\ \frac{1}{r^{|}m|}\end{array}\right),$$

$$U_{0}d=2\left(n+p_m\right),n=0,1,\mathrm{2...}$$

$$U(r)=\{\begin{array}{ll}-\frac{\alpha }{R},& r\le R,\\ -\frac{\alpha }{r},& r>R,\end{array}$$

$$\left(\begin{array}[]{c}\chi_{A}^{\mathrm{in}}\\ \chi_{B}^{\mathrm{in}}\end{array}\right)=\frac{C}{R}\left(\begin{array}[]{c}J_{m}(\alpha r/R)\\ J_{m+1}(\alpha r/R)\end{array}\right),$$

$$\left(\begin{array}[]{c}\chi_{A}^{\mathrm{out}}\\ \chi_{B}^{\mathrm{out}}\end{array}\right)=\frac{C}{R}J_{m}(\alpha)\left(\begin{array}[]{c}1\\ \frac{m+\gamma_{m}+1/2}{\alpha}\end{array}\right)\left(\frac{R}{r}\right)^{\gamma_{m}+1/2},$$

$$\alpha J_{m+1}(\alpha )=\left(m+{\gamma }_m+1/2\right)J_m(\alpha ),$$

$$H=v_{\mathrm{F}}\mathbf{\sigma }\cdot {\hat{\mathbf{p}}}_1\oplus v_{\mathrm{F}}\mathbf{\sigma }\cdot {\hat{\mathbf{p}}}_2+V,$$

$$E=\pm v_{\mathrm{F}}\sqrt{p_{x_1}^2+p_{y_1}^2}\pm v_{\mathrm{F}}\sqrt{p_{x_2}^2+p_{y_2}^2}.$$

$$H_{K}=\scalebox{0.7}{\mbox{$\displaystyle\begin{bmatrix}0&-\partial_{z}+\frac{K}{2}e^{-\mathrm{i}\Theta_{K}}&\partial_{z}+\frac{K}{2}e^{-\mathrm{i}\Theta_{K}}&0\\ -\partial_{\bar{z}}+\frac{K}{2}e^{\mathrm{i}\Theta_{K}}&0&0&\partial_{z}+\frac{K}{2}e^{-\mathrm{i}\Theta_{K}}\\ \partial_{\bar{z}}+\frac{K}{2}e^{\mathrm{i}\Theta_{K}}&0&0&-\partial_{\bar{z}}+\frac{K}{2}e^{-\mathrm{i}\Theta_{K}}\\ 0&\partial_{\bar{z}}+\frac{K}{2}e^{\mathrm{i}\Theta_{K}}&-\partial_{\bar{z}}+\frac{K}{2}e^{\mathrm{i}\Theta_{K}}&0\end{bmatrix}$}},$$

Full text

1 Introduction

After the turn of the 21st century, the discovery of Dirac materials led to the Dirac equation becoming a cornerstone of modern mesoscopic physics [1], [2], [3], [4]. Starting with the experimental realization of graphene [5], and followed swiftly after by the unearthing of topological insulators [6], the study of massless Dirac fermions became increasingly important from both fundamental and applied perspectives.

One of the most famous phenomena associated with Dirac particles is Klein tunneling [7], an effect in relativistic quantum mechanics whereby a particle impinging on a barrier tunnels into the Dirac sea of antiparticles, and in doing so heavily suppresses the chances of reflection [8]. In Dirac materials, there is a direct analogy when electrons inter-band tunnel into hole states [9]. The most famous result in this area is the perfect transmission of Dirac electron normally incident on a high potential barrier [10] [11], which has been observed experimentally in graphene [12].

The appearance of Klein tunneling in Dirac materials raises the question of how to confine massless charge carriers. This problem is of crucial importance for future electronic devices which require the complete control of these elusive particles - most clearly for transistors, which require well-defined on-off digital logic. Various ideas have been proposed to induce bound states of massless Dirac fermions, including the application of magnetic fields, strain engineering, chemical doping and the opening up of a band gap, as reviewed in Ref. [13]. However, these methods in some sense blunt the key attribute which made Dirac materials so attractive initially, namely their high electron mobility. Therefore, a hope remains that electrostatic confinement of ultrarelativistic particles is possible with negligible cost to their superlative transport properties.

In this brief review, we recount the efforts made to trap Dirac particles electrostatically. In particular, we show how considerations of zero-energy single particle states in scalar potentials allows one to indeed uncover bound states, albeit with some restrictions [14], [15], [16], [17]. We go on to describe the ultra-relativistic two body problem [18], detail how the coupling of two massless Dirac fermions is indeed possible, and outline its potential consequences [19]. We also comment on atomic collapse in Dirac materials [20], and present an analytical solution for the two-body Dirac-Kepler problem, which has implications for ongoing research into the full many-body problem [21].

The rest of this work is organized as follows. In Sec. 2, we provide the continuum theory behind the absence of bound states in Dirac materials. We highlight how electrostatic confinement is indeed possible at zero-energy in Sec. 3, and discuss the latest experiential work in the area in Sec. 4. We introduce the formalism of the two-body problem in Sec. 5, where we provide the exact solution of the Dirac-Kepler problem and review the possible nontrivial pairing at zero-energy. In Sec. 6 we outline competing proposals to create bound states in Dirac materials, and Sec. 7 is devoted to a brief survey of zero-energy states in wider physics research. Finally, we draw some conclusions in Sec. 8, and sketch out some perspectives for this direction of research topic.

2 On bound states with the massless Dirac equation

Here we examine the possibility to obtain bound states from the addition of scalar potentials to the massless Dirac equation. Our modus operandi is asymptotic analysis of the Dirac equation for arbitrary potential wells, where we search for square-integrable solutions associated with bound states [Sec. 2.1]. We treat the Coulomb problem separately due to its long-range nature [Sec. 2.2].

The single particle 2D Dirac-Weyl Hamiltonian reads [5], [6]

[TABLE]

where $v_{\mathrm{F}}$ is the Fermi velocity of the material, the momenta $\bm{\hat{p}}=(\hat{p}_{x},\hat{p}_{y})$, $\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ are Pauli’s spin matrices and $V(r)$ is a scalar potential. We shall act on the Hamiltonian (1) with a two-component spinor wavefunction in the form

[TABLE]

where the quantum number $m=0,\pm 1,\pm 2,...$ is directly related to the total angular momentum quantum number $j_{z}=m+1/2$. Explicitly, the spinor (2) satisfies $J_{z}\Psi=j_{z}\Psi$, with $J_{z}=-\mathrm{i}\partial_{\theta}+\sigma_{z}/2$. The wavefunction (2) separates the variables in the Schrödinger equation $H\Psi=E\Psi$, leading to a pair of coupled equations for the radial components

[TABLE]

where the scaled eigenvalue $\varepsilon=E/\hbar v_{\mathrm{F}}$ and scaled potential function $U(r)=V(r)/\hbar v_{\mathrm{F}}$. Decoupling the two first-order equations (3) leads to the following second-order differential equation for the upper component of the radial wavefunction

[TABLE]

here ′ denotes taking a derivative with respect to $r$. The corresponding equation for $\chi_{B}$ is obtained with the replacements $A\to B$ and $m\to-(m+1)$ in Eq. (4). In the following two subsections, we consider the asymptotics of Eq. (4), as well as the special case of the Coulomb problem.

2.1 Long-range asymptotics

As pointed out by Tudorovskiy and Chaplik [22], the large-$r$ asymptotics of Eq. (4) for some (faster than Coulomb) decaying potential well is exactly Bessel’s differential equation

[TABLE]

since the terms with $U(r\to\infty)$ and $U(r\to\infty)^{\prime}$ can be safely neglected. However, unlike massive particles as governed by the 2D Schrödinger equation, the scaled energy appears here as squared ($\varepsilon^{2}$). Therefore the asymptotic solutions of Eq. (5) are independent of the sign of the energy, and thus behave like the scattering (and not the bound state) solutions of the equivalent problem for non-relativistic particles. Explicitly, the solutions of Eq. (5) are the superposition

[TABLE]

in terms of the Bessel functions of the first [$J_{m}(z)$] and second [$Y_{m}(z)$] kinds of order $m$, and where $c_{1,2}$ are some constants. The solution (6) suggests a slow asymptotic decay like $\chi_{A,B}\propto 1/\sqrt{r}$, which is characteristic of scattering states. Thus the question of how to create bound states, which are defined by square-integrable wavefunctions with the typical exponential decay $\chi_{A}\propto e^{-r/r_{0}}$, is raised.

2.2 One-particle Coulomb problem

The above asymptotic analysis does not hold for long-range potentials, although it can be shown it does not change the fact that discrete states are seemingly impossible to create.

We consider the 2D Dirac-Kepler problem [23] [20], [24], [25], [26], [27], [28], [29] with the Coulomb potential

[TABLE]

where $\alpha$ is the dimensionless strength parameter (and indeed the effective fine structure constant for 2D Dirac materials). One may find from Eqs. (3) the solutions

[TABLE]

where the somewhat complicated functions $f_{A,B}(\varepsilon r)$ are derived explicitly in Ref [24]. In Eq. (8), the ‘atomic collapse’ parameter

[TABLE]

defines two regimes of the problem: sub-critical ($\alpha<\alpha^{\mathrm{c}}_{m}$) with real $\gamma_{m}$; and super-critical ($\alpha>\alpha^{\mathrm{c}}_{m}$) with imaginary $\gamma_{m}$. The critical strength of the Coulomb potential is given by the $m$-dependent quantity

[TABLE]

with important threshold value $\alpha^{\mathrm{c}}_{0}=1/2$. In the super-critical regime, there are pathological oscillations in the wavefunction at short range ($r\to 0$), mimicking the celebrated atomic collapse phenomena of an electron falling into the nucleus [30], [31]. Furthermore, so-called quasi-bound states may also appear [20], however there are no truly bound states with square-integrable wavefunctions, as can be seen from the $\chi_{A,B}\propto 1/\sqrt{r}$ asymptotic decay in Eq. (8) with $r\to\infty$.

The hallmarks of atomic collapse are believed to have been seen in two landmark experiments. In Ref. [32], the supercritical Coulomb regime was reached by pushing clusters of calcium dimers together on the surface of gated graphene devices, with the aid of a scanning tunneling microscope (STM). STM spectroscopy was then used to probe the atomic collapse-like states as a function of the Fermi level, which is tunable via electrostatic gating. In Ref. [33], the authors exploit a single-atom vacancy to create a stable supercritical charge. Using an STM tip to apply voltage pulses to the vacancy, the positive charge can be increased so that nearby conduction electrons experience the supercritical phase. Resultant measurements through STM microscopy and Landau level spectroscopy suggest the signatures of atomic-collapse states in the density of states.

It should also be mentioned that results of the atomic cluster experiment of Ref. [32] leave some option questions. For example, the peak in the density of states for critically charged clusters appears precisely at the Dirac point, rather than a few electron volts below it, and with the further increase of charge it even moves in the ‘wrong’ direction above the Dirac point. In addition, the density distribution in the ‘collapsed’ states is crater-like and covers tens of lattice constants, which are all signatures of a fully-confined zero-energy state as discussed in the next section.

The atomic collapse observed in graphene is not directly analogous to the standard textbook case of atomic collapse [34], since the particles are massless and hence there is no bandgap. However, the characteristic feature of bound states diving into the continuum can be seen with the massive Dirac particles, as found in certain carbon nanotubes for example [35].

3 Zero energy states of the massless Dirac equation

The no-go theorem of the Sec. 2.1 in fact requires an important caveat, which we detail here. This stipulation opens up a route to create bound states electrostatically, a feat against popular and widely held beliefs.

The impossibility of supporting bound states of massless Dirac fermions was based on an asymptotic analysis of Eq. (5). However, so as to avoid a reduction to the scattering states of Eq. (6), one may consider the special case of zero-energy states ($\varepsilon=0$), which were previously implicitly ignored. Then the asymptotic solution of Eq. (5) is instead $\chi_{A}\propto 1/r^{|m|}$, such that bound states may indeed form for certain values of angular momentum as indexed by $m$ [36], [37]. The resultant algebraically decaying wavefunctions are a characteristic of this problem, in stark contrast to the usual exponential decay of non-relativistic problems. Furthermore, we note that the sign of the potential is not important for the quantization condition, since the potential enters Eq. (5) as either squared or as a logarithmic derivative for $\varepsilon=0$.

In what follows, we investigate two simple models which elucidate the features of such zero energy bound states [Secs. 3.1 and 3.2], before analyzing the zero-energy Coulomb problem [Sec. 3.3].

3.1 Finite circular well

Let us start by considering the simplest possible model, of bound states in the finite circular well [38], [39], [40]

[TABLE]

where $U_{0}$ and $a$ are the strength (in units of inverse length) and spatial range of the well respectively, and $\Theta(z)$ is the Heaviside step function.

Solving the coupled equations (3) inside the potential well ($r<a$) yields

[TABLE]

where $C_{nm}$ is a normalization constant to be determined. Outside the potential well ($r>a$), one finds

[TABLE]

The requirement of continuity of both of the radial wavefunction components leads to the condition $J_{m}(U_{0}a)=0$, such that bound states appear at

[TABLE]

where $\alpha_{m,n}$ is the $n$-th zero of the Bessel function of the first kind of order $m$. Therefore, there are a countably infinite number of zero-energy bound states, which occur at certain critical values of the potential strength. Crucially, the states with $m=0$ are extended states (appearing at $U_{0}a\simeq 2.41,5.52,8.65...$), since they are marginally non-integrable due to their slow algebraic decay [from Eq. (13), $\chi_{B}\propto 1/r$]. The first square-integrable bound state is for $m=1$ (where $\chi_{B}\propto 1/r^{2}$), at the threshold potential strength $U_{0}a\simeq 3.83$. Further bound states successively appear with increasing the potential strength $U_{0}a$. All of the bound-state supporting values of Eq. (14) are distinct, since the Bessel function of the first kind of integer order $J_{m+n}(z)$ does not have common zeros with $J_{m}(z)$ as follows from Bourget’s hypothesis [41].

Most notably, Eq. (13) displays a complete suppression of probability density on the upper wavefunction component ($\chi_{A}$), in a direct manifestation of the chirality of the system. When the model describes electrons in a honeycomb lattice, like in graphene, this means the electron density is only on the B sublattice of the honeycomb structure.

The analysis leading to Eq. (13) may be extended for negative values of $m$, yielding the solution

[TABLE]

Now the outer radial wavefunction demonstrates the absence of probability density on the lower wavefunction component ($\chi_{B}$), due to its zero-eigenvalue condition $J_{m+1}(U_{0}a)=0$, or $U_{0}a=\alpha_{m+1,n}$. The solution (15) shows that states with $m=-1$ are extended and those with $m\leq-2$ are bound states. All of these features highlight the intrinsic $A\to B$ and $m\to-(m+1)$ symmetry of the system, as follows from the governing Eq. (3).

It is instructive to now consider a more realistic, smooth potential well, which has the added benefit of an analytic expression for the supported zero-energy states.

3.2 Lorentzian well

We shall now consider fully-confined zero-modes trapped in the smooth confining potential [36] [88]

[TABLE]

where $U_{0}$ and $d$ parametrize the strength (in units of inverse length) and spatial extent of the potential respectively.

We wish to solve the coupled Eq. (3) by initially finding just the upper radial wavefunction component $\chi_{A}$. A short-range analysis of Eq. (4) suggests that at $r\to 0$ our solution should be of the form $\chi_{A}\sim r^{|m|}$; whilst a similar long-range analysis suggests a decay like $\chi_{A}\sim r^{|m|-2p_{m}}$ as $r\to\infty$. Here we have introduced the parameter

[TABLE]

We switch to a new variable $\xi=(r/d)^{2}$, and try an appropriate ansatz solution of Eq. (4) in the form

[TABLE]

where $w(\xi)$ is an unknown function, having no bearing on the small or large asymptotics of the solution. We find from Eq. (4) the following equation for $w(\xi)$

[TABLE]

The solution of Eq. (3.2) can be given in terms of the Gauss hypergeometric function $\,{}_{2}F_{1}\left(a,b;c;z\right)$ as

[TABLE]

Thus we have found, from Eq. (18) with Eq. (20), the explicit form of $\chi_{A}$, while $\chi_{B}$ is readily obtained from Eq. (3b). It is implicit that the power series in Eq. (20) must be terminated, so as to satisfy the required conditions on its limiting behavior determined at the outset [given above Eq. (17)]. Given this termination, the radial asymptotics are

[TABLE]

which highlights the algebraically decaying, but marginally non-square-integrable, extended states with $m=\{0,-1\}$ and the otherwise bound states with $m\neq\{0,-1\}$, as was the case of the finite circular well of the previous Subsection 3.1.

We obtain the conditions for bound states by assigning the second argument in Eq. (20) to be a non-positive integer, which results in

[TABLE]

where $p_{m}$ is defined in Eq. (17). Therefore a staircase of bound states is formed with increasing potential strength $U_{0}d$, starting from the threshold of $U_{0}d=4$, since $U_{0}d=2$ only supports extended states with $m=0$ or $m=-1$. Clearly, all of the modes of Eq. (22) appear at even integer values of $U_{0}d=2N$, where $N=1,2,3,...$ and possess an accidental degeneracy of $2N$. For example, when $U_{0}d=2$ there is a two-fold degeneracy with the states $(n,m)=(0,0)$ and $(0,-1)$, while for $U_{0}d=4$ there is a four-fold degeneracy, encompassing $(n,m)=(0,1)$, $(1,0)$, $(0,-2)$ and $(1,-2)$.

Our analyses of the toy models in this section allows us to gather together some general conclusions. In order to have bound states, we require: zero-energy states ($\varepsilon=0$); a scalar potential $U(r)$ decaying faster than $1/r$, a critical strength $U_{0}$ and spatial extent $d$ of the well, and a large enough angular momentum, or vorticity, $m\neq\{0,-1\}$. Therefore, we term these exotic bound states ‘zero-energy vortices’.

3.3 Coulomb problem

It is also worth considering zero-energy states with the Coulomb potential. At finite energy ($\varepsilon\neq 0$) the exact solution is in terms of non-square-integrable eigenfunctions, leading to the absence of bound states as was discussed in Sec 2.2. However, the special case of zero-energy states ($\varepsilon=0$) which interests us were not considered before. Indeed, zero-energy states were implicitly excluded by the variable choice $\varepsilon r$ in Eq. (8).

We study the regularized Coulomb potential

[TABLE]

with the (dimensionless) strength $\alpha$ and length scale $R$. Inside the cutoff region ($r<R$) we find

[TABLE]

where $C$ is the normalization constant. Outside the cutoff region ($r<R$), we find

[TABLE]

where we have used the atomic collapse parameter defined in Eq. (9). In order to describe a sufficiently fast decaying wavefunction, one requires $m(m+1)>\alpha^{2}$. The wavefunction continuity condition on $\chi_{B}$ suggests that bound states should appear at certain values of $\alpha$, which satisfy

[TABLE]

which is notably $R$ independent. However, this transcendental equation (26) has no solutions, as one can readily see graphically. Therefore we conclude that even at zero-energy, long-range potentials are not able to support bound states.

We should mention that the pure Coulomb potential is of questionable relevance for realistic low-dimensional materials since it does not account for dielectric screening, which leads to its modification as discussed by Rytova and Keldysh for thin films [42] [43] and by Cudazzo and co-workers for 2D materials [44]. Indeed, the problem of zero energy states in screened potentials is a particularly fruitful avenue for further research [45].

3.4 Adding a magnetic flux tube

The addition of a magnetic flux tube to the Hamiltonian (1) via the vector potential $\bm{A}=(0,\hbar f/er,0)$, where $f=\Phi/\Phi_{0}$ is the number of flux quanta $\Phi_{0}=h/e$, provides a useful extra degree of freedom [46]. After introducing a generalized momentum $\bm{\hat{p}}\to\bm{\hat{p}}+e\bm{A}$ to incorporate such a vector potential, one obtains same coupled equations as without a flux tube [c.f. Eqs. 3] after the transformation $m\to\tilde{m}=m+f$. Therefore, bound states are no longer determined by just the angular momentum quantum number, with the number of flux quanta inside the flux tube giving rise to a degree of tunability which may be exploited in future experiments.

3.5 Supersymmetric solutions

The experimental breakthroughs in realizing Dirac materials have led to enormous theoretical activity regarding exact solutions of Dirac-like equations. One neat method to obtain bound states for a variety of scalar potentials is through the techniques of supersymmetric quantum mechanics [47], [48], [49] as applied to the Dirac equation [50], [51]. Further analytic solutions have been generated by exploiting the relation between solitonic solutions of the Kortweg-de Vries equation and the Dirac equation [52].

4 Experimental signatures of confinement

The principle signature of bound or quasi-bound states in Dirac materials is through the local density of states data obtained by STM measurements [53], [54]. As mentioned previously, the results of the atomic cluster experiments of Ref. [32] can also be explained by full-confined zero-energy states. In particular, the movement of the apparent maximum in the density of states above the Dirac point for supercritical clusters can be explained by screening of the clusters by free electrons to precisely the critical potential strength value, producing the combined picture of a zero-energy states peak plus an additional screening electron density.

A landmark experiment on confinement in graphene focused on whispering-gallery modes [55], so-called because of how sound waves in the circular chamber walls of St. Paul’s Cathedral in London are confined by the principle of reflection. Resonators were made by circular pn junctions created by a scanning tunneling probe and Klein scattering of the electron waves at the edges of the formed cavity led to resonances in the local density of states.

In Ref. [56], the authors used STM to map the electronic structure of massless Dirac fermions confined by circular graphene pn junctions, allowing for the observation of the energy levels and nodal patterns inside the quantum dot. STM was also used in Ref. [57] to probe the transmission of Dirac electrons through a sharp circular potential well defined by substrate engineering. Later experimental work on confinement in circular pn junctions created by the tip-induced charge, has seen signatures of both atomic collapse states and whispering-gallery modes [58]. Notably, the theoretical model for the tip potential was based on an infinite harmonic well which is known to not have bound states, and the truly bound states at zero-energy were overlooked. Therefore a consistent treatment of whispering gallery modes with a potential reaching a plateau at long range, rather than an infinity, is ripe for investigation.

Monolayer graphene was also studied in Ref. [59], where the authors studied the interplay of a homogeneous magnetic field and electrostatic confinement. An STM tip was used to induce a confining potential in the Landau gaps of bulk graphene, thus circumventing Klein tunneling. A magnetic field was also exploited in Ref. [60], where a dramatic increase in the energy of states in circular graphene p-n junction resonators were observed when a small critical magnetic field was reached, in a manifestation of the topological behavior of Dirac particles.

In Ref. [61], the authors realized nanoscale p-n junctions with atomically sharp boundaries in graphene monolayers by creating a monolayer vacancy island of copper surface, and they saw quasibound states with a finite lifetime in the formed graphene quantum dots. The realization of truly bound states in nanoscale graphene quantum dots was recently reported in [62], where the quantum dots were electronically isolated by boundaries generated by strong coupling between the graphene and the substrate. The experimental quest for the trapping and controllable manipulation of electronic modes in massless Dirac materials remains an active and ongoing pursuit.

5 Two-body problem

We now consider the simplest few-body problem, that of two interacting Dirac particles, in order to gain some insight into the many-body problem. Our aim is to understand the nature of the coupling between massless particles, considering the difficulties for confinement alluded to in the previous sections.

The ultra-relativistic two-body problem in 2D has been studied thoroughly in recent years, mainly from matrix Hamiltonian approaches [18], [63], [64] but also from the Bethe-Salpeter equation [65], [66]. The analogous problem in one spatial dimension has also been extensively studied, for example in Refs. [67], [68], and [69].

A variety of exotic phenomena have been predicted at the few-body level, including: the superconducting coupling of electrons in graphene due to phonons [70]; metastable electron-electron states [71], [72] and exciton-hole states [73], [74] due to trigonal warping; the formation of excitons in double layer graphene [75]; two-particle scattering [76]; and two-body bound states in gapped graphene [77]. Furthermore, it is known that the self-trapping of Wannier-Mott excitons in polar crystals occurs at zero-energy, and the excitons experience a short-range self-trapping potential (decaying like $1/r^{4}$) due to electron-hole charge compensation [78].

5.1 Formalism

The two-particle 2D Dirac-Weyl Hamiltonian reads [18]

[TABLE]

where the momenta of particle $\{1,2\}$ is $\bm{\hat{p}}_{1,2}=(\hat{p}_{x_{1,2}},\hat{p}_{y_{1,2}})$, the interaction is $V=V(\sqrt{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}})$ and $\oplus$ denotes the direct sum. Neglecting for the moment the interaction, there are four eigenvalues from the Schrödinger equation $H\Psi=E\Psi$, which are

[TABLE]

Now with the interaction $V$ taken into account, let us move into center-of-mass [$X=(x_{1}+x_{2})/2$, $Y=(y_{1}+y_{2})/2$] and relative motion [$x=x_{1}-x_{2}$, $y=y_{1}-y_{2}$] coordinates, before transforming into polar coordinates [$(X,Y)\to(R,\Theta)$, $(x,y)\to(r,\theta)$]. We try the solution $\Psi(\mathbf{R},\mathbf{r})=e^{\mathrm{i}\mathbf{K}\cdot\mathbf{R}}\psi(r,\theta)$, where $\hbar\mathbf{K}$ is the center-of-mass momentum. This leads to the equation $H_{K}\psi=\mathrm{i}[\varepsilon-U(r)]\psi$, with

[TABLE]

where $K=\sqrt{K_{x}^{2}+K_{y}^{2}}$, $\tan(\Theta_{K})=K_{y}/K_{x}$, $\partial_{z}=e^{-\mathrm{i}\theta}\left(\partial_{r}-\tfrac{\mathrm{i}}{r}\partial_{\theta}\right)$ and $\partial_{\bar{z}}$ is its conjugate. Upon separating the variables with the following radial spinor wavefunction

[TABLE]

and changing the basis functions to $\chi_{1,2}=R_{1}\pm R_{4}$ and $\mathrm{i}\chi_{3,4}=R_{2}\pm R_{3}$ for simplicity, we obtain a system of four equations to be solved

[TABLE]

Most notably Eq. (31c) has the trivial solution $\chi_{3}=0$, while the remaining trio of equations can be manipuated to provide a single equation for $\chi_{4}$ only:

[TABLE]

while is analogous to the one-particle equation (4). In what follows, we seek bound state solutions of Eq. (32), but first we investigate the Coulomb interaction in the next subsection.

5.2 Coulomb interaction

Most famously for Cooper pairing [79], two-body problems in condensed matter physics are known to give unique insights into the full many-body problem. The ultra-relativistic two particle problem is highly important for a range of phenomena, from critical coupling to pair condensation to spontaneous gap generation [21]. Thus far, only the exact solution for the $m=0$ state has been reported [18], here we seek to derive the general solution for all values of $m$.

We consider the bare Coulomb potential (7) and attempt to first find the fourth wavefunction component $\chi_{4}$ from Eq. (32). An asymptotic analysis of Eq. (32) suggests the behavior $\chi_{4}\propto e^{-\mathrm{i}\varepsilon r/2}$ at large-range, and $\chi_{4}\propto r^{\gamma_{m}-1/2}$ at short-range. Here the two-body atomic collapse parameter [c.f. Eq. (9)] reads

[TABLE]

As in the single particle case of Sec. 2.2, there are two regimes of the interest: sub-critical ($\alpha<\alpha^{\mathrm{c}}_{m}$) and super-critical ($\alpha>\alpha^{\mathrm{c}}_{m}$). The critical strength of the Coulomb potential [c.f. Eq. (10)] is

[TABLE]

with the first few values $\alpha^{\mathrm{c}}_{0,1,2}\simeq\{1,1.41,2.24\}$. These supercritical values are comparable to those reported in the literature by other means. For example, by simulating the phase diagram of graphene as a function of the Coulomb coupling between quasiparticles via Monte Carlo calculations, the critical couplings $\alpha^{\mathrm{c}}_{0}=1.11\pm 0.06$ were found [80] [81], [82]. Renormalization group calculations carried out perturbatively in the interaction strength yield $\alpha^{\mathrm{c}}_{0}\simeq 0.833$ [83]. Analytics in the framework of the Schwinger-Dyson equation gives $\alpha^{\mathrm{c}}_{0}\simeq 1.13$ [84], and $\alpha^{\mathrm{c}}_{0}\simeq 1.62$ from the Bethe-Salpeter equation [65].

Taking into account the aforementioned asymptotics, we undertake a peeling-off procedure with the ansatz

[TABLE]

so that Eq. (32) provides the following equation for $g(r)$:

[TABLE]

Making the natural change of variable $\xi=-\varepsilon r/\alpha$ leads to the standard form of the confluent Heun equation [86] [87]

[TABLE]

where the parameters $\{\rho,\beta,\lambda,\mu,\nu\}$ are given by

[TABLE]

With the more familiar parameter parameterization $(\mu,\nu)\to(\delta,\eta)$ [85], we obtain

[TABLE]

The equation (37) has as its the solution the confluent Heun function $H_{\mathrm{C}}\left(\rho,\beta,\lambda,\delta,\eta,\xi\right)$ [86], or explicitly

[TABLE]

which, with Eq. (35), completes the formal solution of $\chi_{4}$. Therefore, with Eqs. (31) we obtain the full solution of the ultrarelativistic Dirac-Kepler two-body problem, which reads

[TABLE]

where $C$ is a normalization constant. Further analysis of the exact solution of Eq. (41), using relations for the asymptotics of the confluent Heun function, promises to explain the nature of two-body quasibound states, charge criticality and condensation in Dirac materials [21]. Additionally, calculations on the two-body problem beyond the pure Coulomb interaction, to take into account dielectric and even dynamical screening, is an important open problem.

The coupling of two massless Dirac fermions into a truly bound state remains an outstanding and nontrivial task. As can be expected from the form of the principle two-particle Eq. (32), which is highly reminiscent of the key one-particle Eq. (4), one is forced to consider zero-energy states in a bid to bind ultra-relativistic particles, which we now turn to.

5.3 Coupling at zero-energy

Now we consider zero-energy states in a faster than Coulomb potential at zero-energy ($\varepsilon=0$). We choose to work with the Lorentzian interaction of Eq. (16) to gain an analytic understanding of the system [19], and we note that the sign of the interaction is irrelevant for the quantization condition derived from Eq. (32), since at zero-energy the interaction only enters as square or as a logarithmic derivative.

When $r\sim 0$, one finds the typical short-range behavior $\chi_{4}\sim r^{|m|}$. Meanwhile, the asymptotic behavior as $r\to\infty$ is given by the decay $\chi_{4}\sim r^{|m|-2l_{m}}$, where

[TABLE]

is a key dimensionless quantity. Thus, we are motivated to seek a solution of Eq. (32) with the following ansatz and variable change

[TABLE]

Here $f(\xi)$ is a polynomial in $\xi$ that does not affect the short- and long-range behavior of $\chi_{4}$ and it is determinable from Eq. (32) via:

[TABLE]

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

[TABLE]

To ensure decaying solutions at infinity, one needs to terminate the power series in Eq. (45). Then the radial wavefunction asymptotics of the full wavefunction are

[TABLE]

so that states with $m=0$ are marginally non-square-integrable, due to the slowly decaying component $\chi_{2}\propto 1/r$. Therefore truly bound two-particle states have a non-vanishing $m$, and are known as zero-energy vortices.

The quantization condition for the formation of bound pairs follows from the second argument of Eq. (45) as

[TABLE]

where $l_{m}$ is defined in Eq. (42). The marginally non-square-integrable states (with $m=0$) appear at values of $U_{0}d$ divisible by four [$U_{0}d=4(n+1)$], while the bound states (with $m\neq 0$) occur at irrational values, starting from the threshold $U_{0}d\simeq 6.83$. We plot in Fig. 1 the lowest square-integrable two-particle bound state, characterized by the quantum numbers $(m,n)=(1,0)$, showing its vortex-like structure and axial symmetry.

This beautiful result of pairing of massless Dirac particles has some interesting implications [19]. Namely, we see that a pair of ultrarelativistic fermions can form a bound state at zero-energy, if the interaction decays faster than the Coulomb interaction. Since the binding is independent of the sign of the inter-particle interaction, it suggests electron-electron pairing into bielectrons is as possible as electron-hole coupling into excitons. Furthermore, the obeyed bosonic statistics allows for a novel type of condensation, whose characteristics require further study.

6 Confinement by other means

Here we briefly mention other proposals to bind massless Dirac particles, without recourse to zero-energy states, see also the review of Ref. [13].

The most popular method is to apply an external magnetic field perpendicular to the 2D material plane [89], [90]. It was shown theoretically that various magnetic quantum dots are able to sustain either quasibound [91] or fully bound states [92], [93], [94] depending on the spatial configuration. Fictitious or pseudo-magnetic fields, where certain strain configurations give rise to effective magnetic fields and even a type of quantum Hall effect, have also been extensively studied [95], [96], [98].

It was also suggested theoretically that a spatially inhomogeneous Fermi velocity could lead to bound states [101], [102], [103], [104], a situation that might occur due to strain or perhaps superlattice effects. Other schemes have focused on opening up a band gap via substrate engineering [99], or functionalizing the material to change its structure, for example the synthesis of fluorinated graphene [100].

7 Zero modes in physics

We note in passing that zero-energy states are of great importance across several subfields of physics.

The inherent mathematical interest in zero-modes in non-relativistic quantum mechanics has been studied by Makowski and co-workers in Refs. [105], [106], [107]. In scattering theory, bound states at zero-energy are intrinsically linked to scattering phase shifts via the Levinson theorem, which has been investigated with the variable phase method in a series of papers [108], [109], [110], [111].

It has been shown theoretically that zero energy states in superconducting vortices can be fixed by particle-hole symmetry and accomplished with the aid of local lattice distortions [112]. Such lattice deformations have been observed in a spin-polarized neutron scattering experiment on superconducting niobium [113].

Zero-modes of Dirac equations are highly exotic [114], [115], since they can be associated with fractional charge [116], [117], topologically protected edge and soliton-like states [118], [119], and celebrated Majorana modes [120]. Most recently, analogues of some of these zero-energy state phenomena with Dirac-like equations have been both predicted and observed in a diverse array of systems, from photonics [121] to plasmonics [122] to cold atoms [123].

8 Conclusions

In this short review, we have described the problem of creating bound states in Dirac materials. The principle solution to the problem discussed here was the electrostatic confinement of zero-energy states. We have discussed some simple toy models of confinement at the single-particle level, which unveils the unique characteristics of bound states of massless Dirac particles. We have noted that the latest experimental literature displays signatures of zero-energy states in their reported data.

We have also described the ultra-relativistic two-body problem, and provided the solutions for the Dirac-Kepler problem and for the prototypical model of two-particle pairing. In particular, we have described how electron-electron pairing into bielectrons is remarkably both possible and energetically favorable. As a by-product of research into this topic, the number of analytic solutions of Dirac-like equations has been unexpectedly increased [124], showcasing how fundamental theory and applied physics are entwined in Dirac material research.

{acknowledgement}

This work was supported by the Government of the Russian Federation through the ITMO Fellowship and Professorship Program. C.A.D. acknowledges support from the Juan de la Cierva program (MINECO, Spain). M.E.P. acknowledges support from the EU H2020 RISE project CoExAN (Grant No. H2020-644076). We grateful to K. S. Gupta for fruitful discussions and C.A.D. thanks the AVE 03093 from Puerta de Atocha for hospitality during the writing of this paper.