On the density of states of circular graphene quantum dots

TL;DR

This paper presents a straightforward method to compute the local density of states in axially symmetric graphene nanostructures, demonstrating excellent agreement with recent experimental measurements of circular graphene quantum dots.

Contribution

It introduces a simple, effective approach for calculating the local density of states in axially symmetric graphene structures, applicable to quantum dots and rings.

Findings

Method accurately reproduces experimental data

Approach is applicable to various axially symmetric potentials

Results show strong agreement with STM measurements

Abstract

We suggest a simple approach to calculate the local density of states that effectively applies to any structure created by an axially symmetric potential on a continuous graphene sheet such as circular graphene quantum dots or rings. Calculations performed for the graphene quantum dot studied in a recent scanning tunneling microscopy measurement [{\sl Gutierrez et al. Nat. Phys. \textbf{12}, 1069--1075 (2016)}] show an excellent experimental-theoretical agreement.

Figures

Equations

$$\mathcal{H}=\vec{\sigma }\cdot \vec{p}+U(r),$$

$${\Psi }^{(E)}(r,\varphi )=e^{ij\varphi }\left(\begin{array}{c}{e}^{-i\varphi /2}{\chi }_A^{(E,j)}(r)\\ e^{+i\varphi /2}{\chi }_B^{(E,j)}(r)\end{array}\right),$$

$$i\frac{\partial {\chi }^{(E,j)}(r)}{\partial r}={\mathcal{H}}_r{\chi }^{(E,j)}(r),$$

$${\mathcal{H}}_r=\left(\begin{array}{cc}i\frac{j-\frac{1}{2}}{r}& U(r)-E\\ U(r)-E& -i\frac{j+\frac{1}{2}}{r}\end{array}\right).$$

$$\rho (E,r)=\sum _{j=-\infty }^{+\infty }{\rho }^{(j)}(E,r),$$

$${\rho }^{(j)}(E,r)\propto \frac{1}{\Delta E}\Vert {\chi }^{(E,j)}(r){\Vert }^2,$$

$${\chi }^{(E,j)}(r)=W_f(r)C_f,$$

$$W_f(r)=\left(\begin{array}{cc}{J}_{j-\frac{1}{2}}(q_{f}r)& Y_{j-\frac{1}{2}}(q_{f}r)\\ i{\tau }_f{J}_{j+\frac{1}{2}}(q_{f}r)& i{\tau }_f{Y}_{j+\frac{1}{2}}(q_{f}r)\end{array}\right),$$

$$\Delta E=\frac{\pi }{L}.$$

$$\frac{4L\Vert C_f{\Vert }^2}{|E-U_f|}=1.$$

$${\chi }^{(E,j)}(r)=\mathcal{N}\left(\begin{array}{c}{J}_{j-\frac{1}{2}}(q_{i}r)\\ i{\tau }_i{J}_{j+\frac{1}{2}}(q_{i}r)\end{array}\right),$$

$$i\frac{\partial {\tilde{\chi }}^{(E,j)}}{\partial r}={\tilde{\mathcal{H}}}_r{\tilde{\chi }}^{(E,j)}(r)$$

$${\tilde{\mathcal{H}}}_r={\mathcal{H}}_r-\frac{i\beta }{r}.$$

$$\rho (E)={\int }_0^{R_{\max }}4\pi dr\sum _{j=\frac{1}{2}}^{+\infty }{\rho }^{(j)}(E,r),$$

Full text

Quantum dots are among the most intensively studied nano-structures. From the application point of view, it is desirable to create quantum dots by feasible and controllable confinement potentials. For conventional semiconductors, such a confinement potential can be easily realized experimentally, e.g., using an appropriate system of gates. The gate-induced electrostatic potentials can be tuned externally to confine electrons to localized states with some desired properties Chakraborty (1999). As for the mono-layer graphene, due to the Klein tunneling, it has been challenging to experimentally realize the potentials that can induce strictly localized electronic states Katsnelson et al. (2006); Neto et al. (2009). Fortunately, though electrostatic potentials fail to create truly bound electronic states, they can trap the charge carriers in quasi-bound states (QBSs) with a trapping time long enough to satisfy application requirements Matulis and Peeters (2008). Thus, various confinement potential models have been probed to seek for appropriate QBS-structures Matulis and Peeters (2008); Chen et al. (2007); Hewageegana and Apalkov (2008); Recher et al. (2009); Bardarson et al. (2009); Downing et al. (2011); Schneider and Brouwer (2014); Wu and Fogler (2014); Giavaras et al. (2009); Ramezani Masir et al. (2011); Giavaras and Nori (2012); Kim et al. (2012); Gutiérrez-Rubio and Stauber (2015). Notably, most of the potentials probed Matulis and Peeters (2008); Chen et al. (2007); Hewageegana and Apalkov (2008); Recher et al. (2009); Bardarson et al. (2009); Downing et al. (2011); Schneider and Brouwer (2014); Wu and Fogler (2014); Giavaras et al. (2009); Ramezani Masir et al. (2011); Giavaras and Nori (2012); Kim et al. (2012); Gutiérrez-Rubio and Stauber (2015) are axially symmetric, implying that the examined graphene quantum dots are circular in shape (circular graphene quantum dots - CGQDs).

Each QBS is characterized by its energy and trapping time, expressing respectively as the energy and the width of a resonance emerging in the local density of states (LDOS). One can therefore identify QBSs by analyzing the structure of the LDOS Matulis and Peeters (2008). Alternatively, one can also directly find the energy spectrum of QBSs by solving the Dirac equation with an outgoing wave boundary condition. Generally, the QBS spectrum is then complex: while the real parts give the energy positions of QBSs, the imaginary parts give the inverse of their trapping times Chen et al. (2007); Hewageegana and Apalkov (2008); Nguyen et al. (2009).

The interest in CGQDs has particularly raised in the recent tunneling spectroscopy measurements Zhao et al. (2015); Freitag et al. (2016); Gutiérrez et al. (2016); Lee et al. (2016). It was suggested that the tip of a scanning tunneling microscopy (STM) can be finely adjusted to create a quantum dot on the continuous graphene sheet Zhao et al. (2015); Freitag et al. (2016); Lee et al. (2016). It was claimed that all the graphene quantum dots realized in these experiments are practically circular Zhao et al. (2015); Freitag et al. (2016); Lee et al. (2016). Impressively, STM is also the tool to detect the LDOS of graphene quantum dots with high precision. In fact, it has been used to explore the electron whispering-gallery mode resonators Zhao et al. (2015) and to directly image the wave functions of QBSs Zhao et al. (2015); Gutiérrez et al. (2016); Lee et al. (2016).

To theoretically describe the aforementioned experimental data, one has to calculate the LDOS for the CGQD of interest. In Refs. Gutiérrez et al. (2016) and  Lee et al. (2016) the LDOS was calculated using the scattering and the the finite different methods, respectively. It was claimed in Ref. Gutiérrez et al. (2016) that the experimental data agree well with the calculated LDOS, except for the QBS of lowest angular momentum. This QBS has made a puzzle by experimentally appearing at the energy considerably higher than theoretically predicted.

In the present work, stimulated by the beautiful STM measurements, we suggest an approach to efficiently calculate the LDOS of any realistic CGQD created by an axially symmetric electrostatic potential, avoiding the indirect calculation of scattering coefficients or the computationally expensive finite difference method. As illustrations, we calculate the LDOS in two typical cases of step and smooth confinement potentials. In the former case, the LDOS was calculated for the CGQD measured/calculated in Ref. Gutiérrez et al. (2016). Our results describe very well the whole experimental QBS spectrum reported in Ref. Gutiérrez et al. (2016). Notably, although not as high as the observed value, our results suggest that the QBS of lowest angular momentum is actually expected to be at the energy higher than that calculated in Ref. Gutiérrez et al. (2016); some part of the puzzle is therefore resolved. The factor that makes our approach different from that used in Ref. Gutiérrez et al. (2016) is clarified. For the studied CGQDs, we also show that the resonance widths (RWs) extracted from our calculated LDOS practically coincide with those obtained from the corresponding complex QBS spectrum of the Dirac equation and qualitatively describe the experimental data. In the case of smooth potentials, we calculate the LDOS for CGQDs created by the Lorentzian shape potential, which is believed to describe the potential induced by a charged STM-tip Downing et al. (2011).

For a general CGQD, the Hamiltonian that describes the low energy properties of trapped electrons has the Dirac-Weyl form:

[TABLE]

where $\vec{\sigma}=(\sigma_{x},\sigma_{y})$ are Pauli matrices, $\vec{p}=-i(\partial_{x},\partial_{y})$ is the $2$-dimensional momentum operator, and $U(r)$ is an axially symmetric potential. For the simplicity, we restrict to the case of experiments Zhao et al. (2015); Freitag et al. (2016); Gutiérrez et al. (2016); Lee et al. (2016) where the valley scattering can be neglected and use units such that $\hbar=1$ and the Fermi velocity $v_{F}=1$.

To calculate the LDOS for the studied CGQDs, one has to solve the eigenvalue equation of Hamiltonian (1) with a proper normalization. Suppose $E$ and $\Psi^{(E)}(r,\phi)$ are the associated eigenvalue and eigenfunction of this Hamiltonian. Since the potential $U(r)$ is axially symmetric, the eigenfunction $\Psi^{(E)}(r,\phi)$ can be found in the form

[TABLE]

where the total angular momentum $j$ takes half-integer values and $\chi_{A/B}^{(E,j)}(r)$ are the radial wave functions on the graphene $A/B$-sublattices. The radial wave function $\chi^{{(E,j)}}(r)=(\chi^{{(E,j)}}_{A}(r),\chi^{{(E,j)}}_{B}(r))^{T}$ follows the equation

[TABLE]

where the formal radial Hamiltonian $\mathcal{H}_{r}$ is defined by

[TABLE]

Certainly, because of the circular symmetry of the structure, the LDOS also depends only on the radial coordinate $r$ and can be found as

[TABLE]

with

[TABLE]

where $\Delta E$ is the level spacing at the energy $E$ and $\chi^{{(E,j)}}(r)$ has to be subjected to a proper normalization condition. However, for the considered quantum dots, states are only quasi-bound; strictly speaking, the energy spectrum is continuous and the wave function cannot be normalized. To introduce the level spacing $\Delta E$ and the normalization condition for $\chi^{{(E,j)}}(r)$, we follow the approach suggested in Ref. Matulis and Peeters (2008). In this approach, the quantum dot is imagined to be embedded in a fictitious large graphene disk of radius ${L}$. This effectively replaces the continuous energy spectrum by dense discrete levels. Note that these discrete levels are independent of the local potential applied to the graphene disk to create a quantum dot. The applied potential however changes the wave functions, and thus the electronic density locally. The LDOS describes this perturbation of the electronic density around the quantum dot (relative to the uniform density away from the potential).

As the disk is so large that for much of its area, the potential $U(r)$ is practically flat. Therefore one can assume there exits some distance $r_{f}\ll{L}$ such that for $r\geq r_{f}$, the potential could be considered constant, $U(r)\equiv U_{f}$. Consequently, for $r\geq r_{f}$, the wave function can be expressed in terms of two integral constants $C_{f}=(C_{f}^{(1)},C_{f}^{(2)})^{T}$:

[TABLE]

where

[TABLE]

with $q_{f}=|E-U_{f}|$, $\tau_{f}=\operatorname{sign}(E-U_{f})$, and $J_{j\pm\frac{1}{2}}$ and $Y_{j\pm\frac{1}{2}}$ are Bessel functions of the first and the second kind, respectively. The two columns of this $W_{f}(r)$-matrix are just the two independent basic solutions to the radial Hamiltonian $\mathcal{H}_{r}$ in the region considered (see Ref. Nguyen et al. (2016) for the details). Then, using the fact that the wave function vanishes at $r={L}$, one finds the level spacing to be Matulis and Peeters (2008)

[TABLE]

Next, the normalization condition for the wave function can be found by requiring that the integration of the electronic probability density over the whole disk must be $1$. Note that for much of the large disk area outside the quantum dot, the wave function is of the form (7). Although the electronic density in this area is small, it spans the whole (fictitious) macroscopic disk. Therefore, the electronic density outside the quantum dot gives the main contribution to the integration taken over the disk, while the contribution from the relatively small area inside the quantum dot can be ignored Matulis and Peeters (2008). This ultimately results in the following normalization condition Matulis and Peeters (2008):

[TABLE]

The remaining problem is to find an appropriate initial condition so that the differential eq. (3) can be solved. The case when the potential can be considered to be flat near the origin, namely, $U(r)=U_{i}$ for $r\leq r_{i}$, has been studied using the $T$-matrix method Nguyen et al. (2016). In this case, the eigenfunction of eq. (3) near the origin has the simple form

[TABLE]

with $q_{i}=|E-U_{i}|$, $\tau_{i}=\operatorname{sign}(E-U_{i})$ and $\mathcal{N}$ being the normalization coefficient. One then can just take the solution (11) at $r_{i}$ as the initial values and solve eq. (3) for $\chi^{{(E,j)}}(r)$. The normalization coefficient is found by imposing the condition (10), where $C_{f}$ is related to $\chi^{{(E,j)}}(r_{f})$ by eq. (7). With the wave function normalized, the LDOS can be now calculated using eqs. (5) and (6).

Though the assumption that the radial potential $U(r)$ is flat near the origin ($r<r_{i}$) is really observed in different CGQD-models Matulis and Peeters (2008); Hewageegana and Apalkov (2008); Recher et al. (2009); Bardarson et al. (2009); Gutiérrez et al. (2016), with regard to the confinement potentials as those induced by the STM-tip in experiments reported in Ref. Zhao et al. (2015); Freitag et al. (2016); Lee et al. (2016), there is a need to relax this assumption. Note that, on the other hand, any electrostatic confinement potential should tend to be constant at large distances.

When the potential $U(r)$ is not flat near the origin, the wave function (11) is no longer an exact solution of eq. (3). However, if $U(r)$ is continuous near the origin, the asymptotic form of (11) still correctly describes the asymptotic behaviour of the solution, namely $\chi^{{(E,j)}}(r)\sim(\propto r^{|j-\frac{1}{2}|},\propto r^{|j+\frac{1}{2}|})$ [Here, to avoid irrelevant factors we use the symbols $\propto$]. However, except for $j=\pm\frac{1}{2}$, this asymptotic solution vanishes at the origin ($r=0$), thus cannot be used as the initial value to solve the differential eq. (3) for $\chi^{{(E,j)}}(r)$. The simple trick to get around this problem is to introduce the regularized wave function $\tilde{\chi}^{{(E,j)}}(r)=r^{-\beta}\chi^{(E,j)}(r)$, where $\beta=\min\{|j-\frac{1}{2}|,|j+\frac{1}{2}|\}$. This regularized wave function $\tilde{\chi}^{{(E,j)}}(r)$ follows the evolution

[TABLE]

with

[TABLE]

The initial condition for the regularized wave function $\tilde{\chi}^{{(E,j)}}(r)$ is now regular, namely, $\tilde{\chi}^{{(E,j)}}(0)=(\tilde{\mathcal{N}},0)^{T}$ if $|j-\frac{1}{2}|<|j+\frac{1}{2}|$ and $\tilde{\chi}^{{(E,j)}}(0)=(0,\tilde{\mathcal{N}})^{T}$ otherwise. Here $\tilde{\mathcal{N}}$ is the normalization coefficient, which is found using eqs. (7) and (10) as described above. With the wave function $\chi^{(E,j)}(r)$ determined from $\tilde{\chi}^{{(E,j)}}(r)$, the LDOS can again be calculated using eqs. (5) and (6).

As a typical illustration for the suggested approach, we calculate the LDOS for the CGQD studied in Ref. Gutiérrez et al. (2016). This CGQD is believed to exhibit a sharp boundary so that the radial confinement potential can be modeled as a step one, $U(r)=V_{0}\Theta({R_{0}}-r)$, where $V_{0}=\operatorname{constant}$, $\Theta(x)$ is the Heaviside step function, and ${R_{0}}$ is the dot radius. In Ref. Gutiérrez et al. (2016), $R_{0}=5.93$ nm, $V_{0}=0.43$ eV and the potential is applied to a piece of graphene with dimension that can be estimated to be $\approx 40$ nm. There it was also shown that this piece can be well treated as an infinite graphene sheet. For this step potential $U(r)$ the LDOS can be calculated from eqs. (5) and (6) with the wave function $\chi^{(j,E)}(r)$ determined directly from eq. (3). Note that without an external magnetic field, the energy spectra are identical for $j>0$ and $j<0$ Chen et al. (2007); Hewageegana and Apalkov (2008), we therefore consider only the case $j>0$.

In Fig. 1$(a)$ and $(b)$, we show respectively the experimental data provided by the authors of Ref. Gutiérrez et al. (2016) and the LDOS calculated using our approach for the same CGQD. Both figures were plotted in the same format, giving a clear view of the full spectrum through a cross-sectional slice of the CGQD. To access to the resonances in the LDOS spectra, we calculate the corresponding experimental and theoretical total density of states (TDOS) via

[TABLE]

where $R_{\mathrm{max}}$ is the maximal radius probed in the experiment (here $R_{\mathrm{max}}=8\ \mathrm{nm}$), which encircles the major maxima of the LDOS. For the theoretical LDOS, the summation over angular momenta is truncated at $j_{\mathrm{max}}=\frac{31}{2}$ as higher momenta do not contribute significantly to the LDOS at this energy scale. The obtained theoretical and experimental LDOSs are presented in Fig. 1(c) by the solid and dashed lines, respectively.

Comparing Fig 1(a) and Fig 1(b), obviously, on the whole, there is a very good agreement on the resonance energies between the experimental data and our theoretical calculation. Particularly, our calculation gives a better agreement for the state of lowest angular momentum $j=\frac{1}{2}$, compared to the calculation reported in Ref. Gutiérrez et al. (2016). A closer analysis shows that the difference between the two calculated results mainly stems from the normalization eq. (10). This normalization was absent from the scattering calculations in Ref. Gutiérrez et al. (2016), but naturally appears in the direct calculation described above. As a particular verification of this intuitive assessment, we compare the partial LDOSs for just the state of $j=\frac{1}{2}$ those calculated: from the scattering approach in Ref. Gutiérrez et al. (2016) [Fig. 1$(d)$]; from eq. (6) discarding the normalization (10) [Fig. 1$(e)$]; and from eq. (6) with the normalization (10) [Fig. 1$(f)$]. Obviously, while the two figures Fig. 1$(d)$ and Fig. 1$(e)$ are practically identical, the normalization pushes the state in Fig. 1$(f)$ closer to the experimental position. It should be noted that, despite this improvement, a small discrepancy between the theoretical prediction and the experimental position of this state still persists. We speculate that, having the largest level width Chen et al. (2007); Matulis and Peeters (2008), the $j=\frac{1}{2}$ state is more susceptible to various fluctuations such as the imperfection of the dot boundary (as suggested in Ref. Gutiérrez et al. (2016)) or the thermal noise, which further affect its position.

Another level that is also particularly interesting is that of $j=\frac{3}{2}$. It shows very small level width in comparison to the other levels. With very long wave length and energy near the zero point, it is in fact resembles the so-called zero-energy bound state Downing et al. (2011).

Although the experimental and theoretical TDOS-lines in Fig. 1$(c)$ show the resonances at almost the same energies, we also notice that the theoretical widths are noticeably smaller than the experimental ones. To assess this experimental-theoretical discrepancy, we determine the widths of both the experimental ($\delta_{\mathrm{exp.}}E$) and theoretical ($\delta_{\mathrm{th.}}E$) resonances in Fig. 1$(c)$ by fitting the TDOS around each resonance to a Lorentzian peak Matulis and Peeters (2008); Gutiérrez et al. (2016). The resonance of $j=\frac{9}{2}$ is excluded due to the low quality of the experimental data. For the rest, the widths obtained from the two TDOSs are compared in Fig. 2. It seems that for all the resonances examined, the experimental widths (vertical axis) are in the same amount of $\approx 0.028\ \mathrm{eV}$ larger than the theoretical ones (horizontal axis). This systematic smearing of resonances (that makes peaks wider and lower) may be caused by, as was already noted in Ref. Gutiérrez et al. (2016), the fact that the electrons in graphene have a non-zero probability of transition into the surface of the copper substrate, which leads to a decrease of the trapping times in QBSs. The thermal noise might be an additional reason for this resonance smearing.

For the same QBSs of $j=\frac{1}{2}$, $\frac{3}{2}$, $\frac{5}{2}$, and $\frac{7}{2}$ of the studied CGQD, we also calculate $(-\operatorname{Im}E)$ of the QBS complex energies $E$ using our $T$-matrix approach suggested in Ref. Nguyen et al. (2016). For a given QBS, the quantity $(-\operatorname{Im}E)$ should provide a direct measure of the resonance width. In the inset in Fig. 2, we show the obtained $(-\operatorname{Im}E)$ in comparison with the resonance widths extracted from the calculated TDOS. We find that for all the resonances examined, the resonance widths extracted from TDOS (horizontal axis) and the corresponding imaginary parts of the QBS complex energies (vertical axis) are in agreement with the relative accuracy of at least $92\%$. This gives an additional confidence to the current discussion.

As another illustration, we calculate the LDOS for the CGQDs of the type that is created by a charged STM-tip in the experiments reported in Refs. Zhao et al. (2015); Lee et al. (2016). For such the CGQDs we follow Ref. Downing et al. (2011) and model the smooth confinement potentials as: $U(r)=V_{0}/[1+(r/\bar{R})^{2}]$ (Lorentzian potentials), where $V_{0}$ and $\bar{R}$ measure the strength and the width of the potential, respectively. Although the LDOS can still be calculated from eqs. (5) and (6), since the potential is not flat in the vicinity of the origin, the wave function $\chi^{(E,j)}(r)$ in these equations should be determined from the regularized one $\tilde{\chi}^{{(E,j)}}(r)$ of eq. (12). We use the ODE45 Shampine and Reichelt (1997) solver to solve this differential equation (with the maximal step size $10^{-2}\bar{R}$). Note that the differential equation for $\tilde{\chi}$ is solved successively from $r_{i}=10^{-8}\bar{R}$ towards infinity; so there is no approximation in this step when one stops the solver at some $r=r_{f}$. In the present LDOS calculation, we stop the solver at $r_{f}=10\bar{R}$ to calculate $C_{f}$ from eq. (7). This is equivalent to (approximately) regarding $U(r)=0$ for $r\geq r_{f}=10\bar{R}$. The TDOS is calculated by integrating eq. (14) in the same area, $R_{\mathrm{max}}=10\bar{R}$ (with the angular momentum series also terminated at $j_{\mathrm{max}}=\frac{31}{2}$). We show in Fig. 3 the LDOS $(a)$ and the corresponding TDOS $(b)$ calculated for the CGQD with the Lorentzian confinement potential of $\bar{R}=30\ \mathrm{nm}$ and $V_{0}=0.2$ eV. Certainly, from this TDOS, we can extract the resonance widths in the same way as presented above. Owing to the lack of detailed experimental data available for comparison, we would like simply to note the rather dense and narrow resonances emerged in Fig. 3(b). The whole spectrum is also very sensitive to both parameters $V_{0}$ and $\bar{R}$.

Thus we have presented an approach to calculate the LDOS of CGQDs. This approach equally applies to practically any structure created by axially symmetric electrostatic potentials on a continuous graphene sheet. It can be easily extended to include a mass term in the Hamiltonian (1) Nguyen et al. (2016). Under an external magnetic field, the current formulation does not however apply directly and further studies are needed; see Recher et al. (2009); Giavaras and Nori (2011, 2012); Kim et al. (2012) for alternative approaches.

Acknowledgement: We would like to thank Abhay Narayan Pasupathy and Christopher Gutierrez for sharing their excellent experimental data and for helpful discussions. We also thank C. Huy Pham for his comments on the draft. This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.02-2015.48.