Anomalous Induced Density of Supercritical Coulomb Impurities in Graphene Under Strong Magnetic Fields

Recently, the supercritical Coulomb impurity problem has been revived [1, 2, 3, 4, 5, 6, 7, 8, 9] in two-dimensional graphene [10]. (The problem of Coulomb impurity in three-dimensional systems was intensively investigated many years ago. For a comprehensive review, see Ref. [11].) It has been experimentally demonstrated that single-atom vacancies in two-dimensional graphene can stably host local charge. Using various experimental techniques, the supercritical regime can be achieved [12, 13]. In the absence of a magnetic field, the induced density [2, 3, 4, 5, 7] has several interesting properties. However, how electron-electron interactions would affect these results is not well-known. The purpose of this paper is to investigate the induced density of the Coulomb impurity in the presence of magnetic fields. The advantage of applying a magnetic field is that, in some cases, it reduces the effect of electron-electron interactions due to an excitation gap between the Landau levels. We find that impurity states exhibit several unusual properties and give rise to an anomalous induced density, in addition to anticrossings that lead to new impurity cyclotron resonances. Before we present our main findings, we provide a brief introduction of the Coulomb impurity problem both in the absence and in the presence of a magnetic field.

The continuum model Hamiltonian of the two-dimensional Coulomb impurity model in the absence of a magnetic field takes the following form

$\displaystyle\mathcal{H}_{0}=v_{F}\vec{\sigma}\cdot\vec{p}-\frac{Ze^{2}}{\kappa r}+\Delta\sigma_{z},$

(1)

where $v_{F}\approx 10^{6}\textrm{ m/s}$ is the Fermi velocity, $\vec{\sigma}=(\sigma_{x},\sigma_{y})$ represents the Pauli spin matrices, $Ze$ is the impurity charge, $\kappa$ is the effective dielectric constant, and $\sigma_{z}$ is the Pauli matrix in the $z$ direction. Additionally, $\Delta$ represents a finite mass gap. In the presence of a strong Coulomb potential, to avoid pathological oscillations of wavefunctions towards the impurity origin, the impurity charge is introduced with a size of $R$. This breaks continuous scale symmetry into discrete scale symmetry (see Appendix A for an explanation of this effect). The coupling strength is defined as the ratio between two energy scales,

$\displaystyle g=E_{C}/E_{D}=\frac{Ze^{2}}{\kappa\hbar v_{F}},$

(2)

where $E_{C}=Ze^{2}/\kappa R$ and $E_{D}=\hbar v_{F}/R.$ In the absence of a magnetic field and zero mass gap $\Delta=0$, subcritical and supercritical regimes separate at the critical coupling strength $g_{c}=1/2$ [2, 4, 1, 6].

Nishida [7] showed that, in the absence of a magnetic field, discrete scale invariance in the induced density in the supercritical regime $|g|>g_{c}=1/2$ has the following form

$\displaystyle\rho(\vec{r})=\sum_{|J|<g}\frac{F_{J}(r/r_{J}^{*})}{r^{2}}+N_{0}\delta(\vec{r}),$

(3)

where $J$ is the total angular momentum and $r_{J}^{*}$ is a $J$-dependent regularization parameter. In the subcritical region, $|g|<1/2$, only the scale-independent first $\delta$-function term is present [2, 4, 3], with the analytical form of $N_{0}$ given in Ref. [5]. The noteworthy feature is that the universal function $F_{J}(r/r_{J}^{*})$ displays log-periodic and discrete scale invariance [7], characterized by

$\displaystyle F_{J}(r/r_{J}^{*})=F_{J}\left(e^{n\pi/\sqrt{g^{2}-J^{2}}}r/r_{J}^{*}\right),$

(4)

where $n$ is an integer. The induced density exhibits a power-law tail, $\rho(r)\sim 1/r^{2}$, as $r\rightarrow\infty$. The role of screening [5] in this induced density in the presence of electron-electron interactions has not been well investigated.

In the Coulomb impurity problem in magnetic fields, regardless of the value of the magnetic field, the critical dimensionless coupling remains constant, $g=g_{c}$ [14]. The problem was investigated both below [15, 16, 14] and above [14, 17, 18] critical coupling $g_{c}$. The continuum model Hamiltonian of the Coulomb impurity problem in a magnetic field reads

$\displaystyle\mathcal{H}=v_{F}\vec{\sigma}\cdot\left(\vec{p}+\frac{e}{c}\vec{A}\right)-\frac{Ze^{2}}{\kappa r}+\Delta\sigma_{z}.$

(5)

The graphene sheet lies in the $xy$ plane, and we use a symmetric gauge with vector potential $\vec{A}=\frac{B}{2}(y,-x)$. The first term gives rise to graphene Landau levels with magnetic length $\ell\equiv\sqrt{\hbar c/eB}$. The Landau level energy in the absence of impurity and the gap $\Delta$ is $\mathcal{E}_{n}/E_{M}=\text{sgn}(n)\sqrt{2|n|}$, where the magnetic energy scale is $E_{M}=\hbar v_{F}/\ell$. Discrete scale invariance is broken because of this magnetic length scale. The probability densities of the Landau level states in the absence of the impurity potential form rings with width $\ell$. However, in the supercritical region, there is another length scale, $R\ll\ell$, as the Coulomb potential must be regularized for supercritical impurity potentials [14]. There are now two relevant dimensionless parameters,

$\displaystyle g=E^{\prime}_{C}/E_{M}=\frac{Ze^{2}}{\kappa\hbar v_{F}}\text{ and }R/\ell,$

(6)

where $E^{\prime}_{C}=Ze^{2}/\kappa\ell$ is the characteristic energy scale of the Coulomb impurity. The coupling strength $g$ is the ratio between the Coulomb energy and the Landau level energy spacing. Notice that despite a magnetic field being considered here, this coupling strength is identical to the one in Eq. (2). When $g$ is large, many Landau levels are coupled by the Coulomb potential. Note that $g$ is independent of $\ell$. The other parameter $R/\ell$ characterizes the regularization parameter of the Coulomb impurity.

In the supercritical region $g>g_{c}=1/2$, the following properties are found. (i) No complex energy solutions (resonances) are possible in the Coulomb impurity problem in magnetic fields: the effective potential does not allow resonant states since the vector potential diverges while the Coulomb potential goes to zero in the limit $r\rightarrow\infty$ [19]. (ii) Regardless of the size of the mass gap $\Delta$, the critical dimensionless coupling strength remains a constant $g=g_{c}$, unlike the case of zero magnetic field [6]. (iii) There can be different types of impurity bound states in a magnetic field: quasi-log-periodic [14] states for $g\gg 1$ and tightly bound states for $g\lesssim 1$, as depicted in Fig. 1.

So far, we have briefly reviewed the basic properties of electronic wave functions of the Coulomb impurity problem. Now, we present our main results for the induced density in the supercritical regime, particularly focusing on strong magnetic fields. We consider only values of the dimensionless coupling strength $g$ where Landau impurity bands do not overlap. In such a case, the Landau level mixing in a filled impurity Landau band due to many-body effects is weak [20]. We investigate how the properties of induced density in the presence of a magnetic field compare to those in its absence. We find that the dimensionless induced density of such a filled Landau impurity band $N$ has the following structure in the supercritical region:

	$\displaystyle\rho_{N}(\vec{r})$	$\displaystyle=\ell^{2}\sum_{J}|\Psi_{N,J}(\vec{r})|^{2}$		(7)
		$\displaystyle=\ell^{2}\sum_{|J|\leq g}|\Psi_{N,J}(\vec{r})|^{2}+\ell^{2}\sum_{|J|>g}|\Psi_{N,J}(\vec{r})|^{2}$

with $\vec{r}$ being the vector position from the impurity charge. The explicit form of $\Psi_{N,J}(\vec{r})$ will be presented below in Eq. (11). Here, the sum over $J$ is for all the states in the $N$th filled Landau impurity band. We find the following similarities and differences in comparison to the zero-field mathematical structure of discrete scale invariance:

1. 1.

   We find, as in the presence of discrete scale invariance, that states with angular momentum $J\leq g$ strongly contribute to the anomalous induced density near $r=0$. For $N=1$ the peak in the induced density is at $r=0$ and is most pronounced. However, for $N=-1$ the peak in the induced density is away from the impurity center. The induced density displays small oscillations for $r\gtrsim\ell$, but without log-periodic oscillations for $r\gg\ell$.

2. 2.

   There is no sharp change in the peak value of $\rho_{N}(r)$ near the critical strength $g_{c}=1/2$. The transition is smooth, but the peak value increases rapidly as $g$ exceeds $g_{c}$.

3. 3.

   The second term of Eq. (7) leads to a unique effect present in a magnetic field: the induced density approaches a constant value $1/(2\pi\ell^{2})$ for $r>d_{s}$, where $d_{s}$ is the screening length. (In the absence of an impurity, the density of a filled Landau level is independent of $r$ and equal to this constant value.)

In addition, Landau impurity band states $\Psi_{N,J}(\vec{r})$ display anticrossings that lead to anomalous impurity cyclotro resonances.

Our paper is organized as follows. In Section II, we explain how our numerical method is implemented using a Hamiltonian matrix. Its eigenvalues and eigenstates that are relevant to the induced density are also explained. The properties of Landau impurity bands are explained in Section III. The induced density is computed in Section IV, and its properties are elucidated. Section V explains some unusual features of impurity cyclotron resonances due to the anticrossing of Landau impurity states. Finally, discussion and a summary are given in Section VI.

We find the eigenstates and eigenvalues of the problem numerically by converting it into the diagonalization of the Hamiltonian matrix.

We introduce the following wavefunctions to construct the basis states of the Hamiltonian:

$$\psi_{n,m}(\vec{r})=c_{n}\begin{pmatrix}-\textrm{sgn}({n})i\phi_{|n|-1,m}({\vec{r}})\\ \phi_{|n|,m}({\vec{r}})\end{pmatrix},$$

(8)

where $n=\ldots,-2,-1,0,1,2,\ldots$ and $m=0,1,2,\ldots$ are, respectively, the inter-Landau-level index and intra-Landau-level index. These two-component states $\psi_{n,m}(\vec{r})$ are graphene Landau level states in the absence of an impurity, and their energy is given by $\mathcal{E}_{n}$ [see below Eq.(5)]. The wave function of each component $\phi_{p,q}(\vec{r})$ is defined in Appendix B. These wave functions are defined only for $p\geq 0$ and $q\geq 0$, and are widely used in ordinary two-dimensional gases in a magnetic field [21]. In Eq. (8), when $|n|-1<0$ (equivalently, $n=0$), $\phi_{|n|-1,m}(\vec{r})=0$ by definition. In this case, only the second component is non zero, and the wave function is chiral. The normalization condition of $\psi_{n,m}$ requires $c_{0}=1$ and $c_{n}=1/\sqrt{2}$ for $n\neq 0$. Note that $\text{sgn}(0)=0$.

In constructing the basis states of the impurity problem in a symmetric gauge, it is useful to utilize the $z$ component of the total angular momentum:

$$J=|n|-m-1/2,$$

(9)

as it is a good quantum number. Using $|n|=J+m+1/2$, we find that the possible values of $J$ are half integers: $\pm 1/2,\pm 3/2,\pm 5/2,\ldots$. (The $z$ component of the total angular momentum operator is $\hat{J}=-i\partial/\partial\theta+\sigma_{z}/2$, where $\theta$ is the polar angle.) Table I lists possible values of $n$ for a given value of $J$.

By relabeling $\psi_{n,m}(\vec{r})$ using index $J$ instead of $m$, we now introduce the basis states of the Hamiltonian matrix:

$\displaystyle\psi^{\prime}_{n,J}(\vec{r})=\psi_{n,m}(\vec{r})=\psi_{n,|n|-J-1/2}(\vec{r}).$

(10)

The eigenstate wave function of the Hamiltonian with eigenenergy $E_{N,J}$ can be expressed as a linear combination of graphene Landau level states as follows:

$\displaystyle\Psi_{N,J}(\vec{r})=\sum_{n}C_{n}^{N,J}\psi^{\prime}_{n,J}(\vec{r}),$

(11)

where $N$ is defined as the Landau impurity band index. Formally, it is defined by taking the limit $g\rightarrow 0$, where an impurity state reduces to a basis state: $\Psi_{N,J}=\psi^{\prime}_{n,J}$. In other words, only one term exists in Eq. (11), and $N=n$. The expansion coefficients $\left\{C^{N,J}_{n}\right\}$ are column eigenvectors.

For a given value of $J$, using these basis states $\psi^{\prime}_{n,J}(\vec{r})$, we form the total Hamiltonian matrix in the Hilbert subspace labeled by $J$. The relevant Hamiltonian consists of a Dirac term, the Coulomb potential, and a mass term:

$$\mathcal{H}_{n,n^{\prime}}=\mathcal{T}_{n,n^{\prime}}+\mathcal{V}_{n,n^{\prime}}+\mathcal{D}_{n,n^{\prime}}.$$

(12)

Since the graphene Landau level states are the basis, the matrix of the Dirac Hamiltonian in a magnetic field is a diagonal matrix with Landau level energies:

$$\mathcal{T}_{n,n^{\prime}}=\text{sgn}(n)\sqrt{2|n|}\delta_{n,n^{\prime}}.$$

(13)

Here, the energy is measured in units of the magnetic energy $E_{M}$. Note that employing the orthogonality in Eq. (30), the matrix elements of the mass term can be computed as

$$\mathcal{D}_{n,n^{\prime}}=\frac{\Delta}{E_{M}}\langle{\psi_{n,m}}|\sigma_{z}|{\psi_{n^{\prime},m^{\prime}}}\rangle=-\frac{\Delta}{E_{M}}\delta_{n,-n^{\prime}}.$$

(14)

Matrix elements of the Coulomb potential are written as

$\displaystyle\mathcal{V}_{n,n^{\prime}}$

$\displaystyle=\langle{\psi_{n,m}}|\frac{-Ze^{2}}{\kappa rE_{M}}|{\psi_{n^{\prime},m^{\prime}}}\rangle=-\langle{\psi_{n,m}}|g\frac{\ell}{r}|{\psi_{n^{\prime},m^{\prime}}}\rangle,$

(15)

which eventually can be simplified to the following form

	$\displaystyle\mathcal{V}_{n,n^{\prime}}=-2\pi gc_{n}c_{n^{\prime}}\Bigg{[}$	$\displaystyle\text{sgn}(nn^{\prime})2^{\alpha_{1}-1/2}A_{|n|-1,m}A_{|n^{\prime}|-1,m^{\prime}}I_{\beta_{1}-\alpha_{1}/2,\beta^{\prime}_{1}-\alpha_{1}/2}\left(\alpha_{1}-\frac{1}{2},\alpha_{1},\alpha_{1}\right)$		(16)
		$\displaystyle+2^{\alpha_{2}-1/2}A_{|n|,m}A_{|n^{\prime}|,m^{\prime}}I_{\beta_{2}-\alpha_{2}/2,\beta^{\prime}_{2}-\alpha_{2}/2}\left(\alpha_{2}-\frac{1}{2},\alpha_{2},\alpha_{2}\right)\Bigg{]},$

where $\alpha_{1}=\left|J-1/2\right|$, $\beta_{1}=\frac{2|n|-J-3/2}{2}$, $\beta^{\prime}_{1}=\frac{2|n^{\prime}|-J-3/2}{2}$, $\alpha_{2}=\left|J+1/2\right|$, $\beta_{2}=\frac{2|n|-J-1/2}{2}$, and $\beta^{\prime}_{2}=\frac{2|n^{\prime}|-J-1/2}{2}$. Here, $A_{n,m}$ is the normalization factor defined in Appendix B [see Eq. (31)].

To derive the above analytical form, we use the following identity of Laguerre polynomials [22]:

	$\displaystyle I_{n,m}(\mu,\alpha,\beta)$	$\displaystyle=\int_{0}^{\infty}t^{\mu}\exp(-t)L^{\alpha}_{m}(t)L^{\beta}_{n}(t)dt$		(17)
		$\displaystyle=\Gamma(\mu+1)\frac{(\beta-\mu)_{n}(\alpha+1)_{m}}{m!n!}\;{}_{3}F_{2}\left(-m,\mu+1,\mu-\beta+1;\mu-\beta+1-n,\alpha+1;1\right),$

where $\Gamma(x)$ is the gamma function, the Pochhammer symbol is defined as $(a)_{n}\equiv\Gamma(a+n)/\Gamma(a)$, and ${}_{3}F_{2}\left(a_{1},a_{2},a_{3};b_{1},b_{2};1\right)$ is the generalized hypergeometric function.

In the supercritical region we must regularize the Coulomb impurity potential by introducing the radius of the impurity charge $R$. For each value of $J$, inter-Landau level numbers within $|n|\leq(N_{c}-1)/2$ are included (we will call $N_{c}$ the Landau level cutoff). The regularization parameter is related to the matrix dimension $N_{c}$ as follows:

$$R\sim\ell\sqrt{2/N_{c}}.$$

(18)

This comes from the fact that the Landau level state with the highest index $N_{c}$ has this minimum length scale, namely, the distance between adjacent nodes in the wavefunction.

Some examples of the expansion coefficients $\left\{C^{N,J}_{n}\right\}$ are given in Fig. 2. Various plots of the probability densities of these eigenstates are shown in Appendix C.

Plots of eigenvalues of Hamiltonian $\mathcal{H}$ as a function of coupling strength $g$ are presented in Fig. 3. Figure 3(a) corresponds to a smaller value of $R/\ell$ compared to Fig. 3(b). The following points can be observed from the plot. Firstly, an impurity splits Landau level degeneracy. This splitting, measured in units of the magnetic energy $E_{M}$, increases as coupling strength $g$ increases. The Landau levels $n=0$ and $n=1$ are mostly affected. The small magnetic field $B$ limit is approached with $N_{c}\rightarrow\infty$, i.e., $R/\ell\rightarrow 0$ [see Eq. (18)]. (Our numerical approach is not suited for investigating this limit because it requires a prohibitively large value of $N_{c}$, meaning a prohibitively large Hamiltonian matrix.) Secondly, there are values of $g$ where Landau impurity bands do not overlap, and we will focus on these ranges of $g$, specifically to the left of the black arrows. Finally, in this range of $g$, the $n=0$ and $n=1$ Landau levels are strongly affected by the change in $R/\ell$. In addition, the Landau level splitting is smaller for a smaller value of $R/\ell$.

Suppose that states of a Landau impurity band are filled, and they do not overlap in energy with other Landau impurity band states. (There are values of $g$ for which Landau impurity bands do not overlap, positioned to the left of the black arrows in Fig. 3.) In such cases, mixing of Landau impurity band states with other band states due to many-body effects is weak, as demonstrated in Refs. [20, 23].

Impurity cyclotron resonance [24] may be used to detect the discrete energy levels in the energy spectrum. The optical matrix elements between the graphene Landau level states in the absence of impurity ($g=0$) are evaluated by using the formula $\vec{j}=v_{F}\vec{\sigma}$ [25]:

$\displaystyle\langle\psi_{n,m}|\sigma_{x}|\psi_{n^{\prime},m^{\prime}}\rangle=M_{nn^{\prime}}\delta_{mm^{\prime}}.$

(21)

Here, we consider the current along the $x$ axis. (We recall that graphene is in the $xy$ plane.) The explicit matrix elements for optical selection rules are given in Tables 2, which implies that $M_{nn^{\prime}}$ is non zero only for $\Delta n=n^{\prime}-n=\pm 1$. Combined with the implied rule of the Kronecker delta $\delta_{mm^{\prime}}$, we can infer that the allowed transitions must satisfy either $\Delta J=1$ or $\Delta J=-1$.

However, in the presence of a Coulomb potential, the optical matrix elements must be evaluated using Landau impurity band states. We find

			$\displaystyle\mathbf{T}_{(N,J)\rightarrow(N^{\prime},J^{\prime})}=\langle\Psi_{N,J}|\sigma_{x}|\Psi_{N^{\prime},J^{\prime}}\rangle$
		$\displaystyle=$	$\displaystyle\sum_{n,n^{\prime}}(C_{n}^{N,J})^{*}\langle\psi^{\prime}_{n,J}|\sigma_{x}|\psi^{\prime}_{n^{\prime},J^{\prime}}\rangle C_{n^{\prime}}^{N^{\prime},J^{\prime}}$
		$\displaystyle=$	$\displaystyle\sum_{n,n^{\prime}}(C_{n}^{N,J})^{*}M_{n,n^{\prime}}\delta_{|n|-J,|n^{\prime}|-J^{\prime}}C_{n^{\prime}}^{N^{\prime},J^{\prime}}$
		$\displaystyle=$	$\displaystyle\sum_{n,n^{\prime}}(C_{n}^{N,J})^{*}\left[M^{(+1)}_{n,n^{\prime}}\delta_{J+1,J^{\prime}}+M^{(-1)}_{n,n^{\prime}}\delta_{J-1,J^{\prime}}\right]C_{n^{\prime}}^{N^{\prime},J^{\prime}},$

where $M^{(\pm 1)}_{n,n^{\prime}}=\pm ic_{n}c_{n^{\prime}}\delta_{|n|\mp 1,|n^{\prime}|}\text{sgn}(n)$ correspond to angular momentum $J$ increasing or decreasing by $1$.

The form $\mathbf{T}_{(N,J)\rightarrow(N^{\prime},J^{\prime})}$ above suggests the possibility of anomalous transitions within the same impurity band, i.e., $\Delta N=0$. In the limit $g=0$ the condition $\Delta N=0$ changes into $\Delta n=0$, which is not optically allowed, as mentioned above. However, $\Delta N=0$ is possible at numerous finite values of $g$ because two impurity Landau levels may cross each other [17], as shown in Fig. 11. We compared our numerical results with the eigenspectrum obtained using the shooting method to solve the Dirac equation as described in Ref. [17], and obtained similar results using $N_{c}=201$.

Let us analyze an example of an anomalous optical transition to gain a better understanding. Its matrix element is given by

$$\mathbf{T}_{(1,-1/2)\rightarrow(1,1/2)}=\langle\Psi_{1,-1/2}|\sigma_{x}|\Psi_{1,1/2}\rangle.$$

(23)

This transition is depicted in Fig. 12(a), and the dependence of $|\mathbf{T}_{(1,-1/2)\rightarrow(1,1/2)}|$ on $g$ is plotted in Fig. 12(b). We observe the following properties. (i) For small values of coupling strength $g$, this optical matrix element is small. It can be explained by noting that with a small coupling strength, this optical matrix element is approximated by transition between graphene Landau level states, $|\psi^{\prime}_{nJ}\rangle=|\psi^{\prime}_{1,1/2}\rangle\rightarrow|\psi^{\prime}_{1,-1/2}\rangle$, which is forbidden because $\Delta n=0$. (ii) For a strong coupling strength, such as with $g>1.5$, the transition $|\psi^{\prime}_{-1,1/2}\rangle\rightarrow|\psi^{\prime}_{0,-1/2}\rangle$, which is allowed because $\Delta n=-1$, contributes significantly to $\mathbf{T}_{(1,-1/2)\rightarrow(1,1/2)}$. (iii) There is a crossover in $\mathbf{T}_{(1,-1/2)\rightarrow(1,1/2)}$ as a function of $g$, that occurs around $g=0.8$, which is a consequence of strong Landau level mixing.

We investigated the induced density of the supercritical Coulomb impurity in the regime where filled Landau impurity bands do not overlap and the effect of electron-electron interactions is significantly reduced. The strong coupling between graphene Landau level states by the impurity potential is non trivial and can lead to several anomalous effects. The induced density of a filled Landau impurity band can exhibit a sharp peak near the impurity center, much narrower than the magnetic length. However, due to strong coupling between graphene Landau levels, this peak can be located away from the center of the impurity, depending on the properties of different Landau impurity bands. We found, like in the presence of discrete scale invariance, that states with angular momentum $J\leq g$ strongly contribute to the induced density near $r=0$. In addition, the impurity charge is screened despite the Landau impurity band being completely filled. We also showed that additional impurity cyclotron resonances exist that involve the anticrossing of Landau impurity band states.

While it is desirable to conduct a Hartree-Fock calculation [23], we do not anticipate qualitative changes in the induced density of a filled Landau impurity band, although some quantitative adjustments in the results are expected [26]. A scanning tunneling microscope [27] may be useful in investigating the anomalous induced density. Impurity cyclotron measurements [24] may also prove useful.