Density of States in Graphene with Vacancies: Midgap Power Law and Frozen Multifractality

The density of states $\varrho (E)$ of graphene is investigated numerically and within the self-consistent $T$-matrix approximation in the presence of vacancies within the tight binding model. The focus is on compensated disorder, where the concentration of vacancies $n_A$ and $n_B$ in both sublattices is the same. Formally, this model belongs to the chiral symmetry class BDI. The nonlinear sigma model predicts for BDI a Gade-type singularity $\varrho (E)\sim |E{|}^{-1}\exp [-|\log (E){|}^{-1/x}]$. Our numerical data are comparable to this result in a preasymptotic regime that gives way, however, at even lower energies to $\varrho (E)\sim E^{-1}|\log (E){|}^{-\tilde{x}}$, $1\le \tilde{x}<2$. We take this finding as evidence that, similar to the case of dirty $d$-wave superconductors, generic bipartite random hopping models may also exhibit unconventional (strong-coupling) fixed points for certain kinds of randomly placed scatterers if these are strong enough. Our research suggests that graphene with (effective) vacancy disorder is a physical representative of such systems.

Figure 1

Density of states of graphene with $\overline{n}=0.1\%–8\%$ vacancies in either sublattice. Comparison of SCTMA and tight-binding simulation.

Figure 2

Data collapse of the inverse time series on a master curve for different vacancy concentrations: $\overline{n}=0.05\%,0.5\%,3\%$. It is seen that at preasymptotic times $1\ll \mathfrak{D}\tau \ll \mathfrak{D}{\tau }_{\overline{n}}^{*}$, Gade-Wegner scaling, Eq. (1), works well within the simulation time window. However, at very long times $\tau \gg {\tau }_{\overline{n}}^{*}$, the increase of $\overline{n}/\varrho (\tau )$ becomes sublinear so that Gade-Wegner scaling is violated. This time regime is numerically accessible at vacancy concentration $\overline{n}\gtrsim 8\%$. Motivated by Ref. [23], we fit $\sqrt{\ln (\tau /{\tau }_0+{\mathfrak{a}}_0)}/A_0$, $\tilde{x}=3/2$, but alternatives [dashed lines; $\ln (\tau )$, $\tilde{x}=2$ and $\ln \ln (\tau )$, $\tilde{x}=1$] are also shown for comparison. The fluctuations in the raw data reflect the stochastic nature of the methodology. (Fitting parameters for fits shown here and in all subsequent figures are listed in the Supplemental Material [24].)

Figure 3

Evolution of $\overline{n}/\varrho (\tau )$ at larger vacancy concentrations ($\overline{n}=3\%–10\%$) where very long (effective) times $\tau \gg {\tau }_{\overline{n}}$ fall within our simulation window. Fits correspond to $\tilde{x}=3/2$ (solid lines) and $\tilde{x}=1$ (dashed lines), defined in Eq. (6).

Figure 4

Top panels: master curves for different $q$ values as obtained after rescaling of $x,y$ axes with energy-dependent scale factors $\xi (E)$ and $Q_q=\xi (E{)}^{{\tau }_q}$ (not shown). Parameters: $\overline{n}=4\%$, $L=64–2048$, $\epsilon ={10}^{-3},5\times {10}^{-5},{10}^{-6},{10}^{-7}$. IPR-distribution functions are given in the Supplemental Material [24]. Bottom panel: multifractal spectrum as estimated from fitting to $Q_q$; it displays frozen multifractality.

Figure 5

Localization length as obtained from the DoS (triangles, data Fig. 1) and from the generalized multifractal analysis (symbols) at vacancy concentration 4%. The conversion of the DoS has been done following Eq. (3) assuming $\tilde{y}=1$. Four fits display Fourier transformation of Eq. (5) (red dashed line, Gade-Wegner-form) and Eq. (6) (blue dot-dashed line: $\tilde{x}=2$; green solid line: $\tilde{x}=3/2$; yellow line: $\tilde{x}=1$). As seen here, Eq. (6) with $\tilde{x}=3/2$ fits all regimes best (using three fitting parameters). Inset: Conversion of DoS into $\xi (E)$ assuming $\tilde{y}=0$. Comparison illustrates that $\tilde{y}$, indeed, enters the data interpretation in an important way, since for $\tilde{y}=0$ only $\tilde{x}=1$ would provide an acceptable fit.