Confined states in graphene quantum blisters

In figures 1(a) and (b), we show schematically the atomic structure of a circular GQB with radius R. In the continuum limit, if the size of the GQB is much larger than the inter-atomic distance, one can describe charge carriers by a $4\times4$ tight-binding Hamiltonian written in the basis of orbital eigenfunctions of the four atoms making up the crystallographic unit cell of bilayer graphene [41]. The labelling of the different atoms is shown in figure 1(b). In the basis $\mathbf{\Psi}=(\Psi_{A 1}, \Psi_{B 1}, \Psi_{B 2}, \Psi_{A 2}){}^{T},$ the Hamiltonian in position representation in valley $K^{\tau}$ reads

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \newcommand{\re}{{\rm Re}} \displaystyle \label{starting_Hamiltonian} \hat{H}^\tau(\vec {r})=\left( \begin{array}{@{}cccc@{}} \tau \delta & v_{\rm F}\hat{\pi}_{+}^{\tau} & \gamma_{1}(\vec {r}) & 0 \nonumber \\ v_{\rm F}\hat{\pi}^{\tau}_{-} & \tau \delta & 0 & 0 \nonumber \\ \gamma_{1}(\vec {r}) & 0& - \tau \delta & v_{\rm F}\hat{\pi}^{\tau}_{-} \nonumber \\ 0 & 0& v_{\rm F}\hat{\pi}^{\tau}_{+} & - \tau \delta \nonumber \\ \end{array} \right). \nonumber \end{align} \tag{ 1 }$$

In equation (1), $v_{\rm F}\sim 10^6~{\rm m}~{\rm s}^{-1}$ is the graphene Fermi velocity [44] and $\tau=(+1, -1)$ for the K  −  and $K'-$valley, respectively. Furthermore, we have neglected skew hopping parameters that only affect the energy spectrum for very small energies [47] and only retain $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \gamma_{1}(\vec {r})$, the inter-layer coupling through van der Waals forces [48]. The quantity δ denotes the potential bias between the two layers and $\hat{\pi}_{\pm}^{\tau}$ are the momentum operators, which in polar coordinates become [29]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{momentum_operator} \hat{\pi}_{\pm}^{\tau}=\frac{\hbar}{\rm i} {\rm e}^{\pm {\rm i}\tau \phi}\left[ \frac{\partial}{\partial \rho} \pm \frac{{\rm i} \tau}{\rho}\frac{\partial}{\partial \phi} \right]. \nonumber \end{align} \tag{ 2 }$$

In equation (2), ρ is the radial distance to the center of the blister and ϕ is the azimuthal angle. Notice that the momentum operator depends on the valley in which the charge carriers reside [49, 50].

In equation (1) the function $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \gamma_{1}(\vec {r})$ describes the coupling between the two graphene layers. If we consider the GQB as in figure 1(a) where the layers depart from each other in the form of a kink, one can describe the inter-layer distance $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} c(\vec {r})$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \newcommand{\re}{{\rm Re}} \displaystyle c(\vec {r}) = \frac{(c_{\rm M}-c_0)}{2} \left[ \tanh\left(\frac{-\rho+R_{\rm QB}}{\xi}\right)+3 \right], \label{inter_layer_distance} \nonumber \end{align} \tag{ 3 }$$

where $c_{0} \sim 0.33~{\rm nm}$ is the equilibrium inter-layer distance, c_M is the inter-layer distance at the center of the blister, R_QB is the radius of the blister, ξ is the width of the interface between delaminated and AB-stacked bilayer graphene, and ρ is the radial component. Because the inter-layer coupling strength $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \gamma_{1}(\vec {r})$ arises from the overlap of two orbital eigenfunctions in the tight-binding formalism, its value decreases exponentially with increasing inter-layer distance. Following standard practice [5, 6, 51], we can write the inter-layer coupling function as

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \newcommand{\re}{{\rm Re}} \displaystyle \gamma_{1}(\vec {r}) = \gamma^{0}_{1} \exp\left(-\beta \frac{c(\vec {r})-c_{0}}{c_{0}}\right). \label{gamma1_coupling} \nonumber \end{align} \tag{ 4 }$$

In equation (4) we have introduced $\gamma_{1}^{0} = 0.38~{\rm eV}$ [48, 52, 53] the equilibrium inter-layer coupling. The quantity $\beta/c_{0}$ is the inverse inter-layer coupling decay length. For the calculations in this paper we choose $\beta \sim 13.3$, as was used before to match with the values for the skew hopping parameters in twisted bilayer graphene [51]. However, this value can be even larger when the blister is formed by insertion of nanoclusters in-between two graphene layers as these clusters screen the van der Waals interaction between the layers.

By plugging equation (3) into equation (4), one can calculate the radial dependence of the inter-layer coupling for a GQB. In figure 1(c) we compare the lateral profile of the GQB to the strength of the inter-layer coupling at the same positions. Because the magnitude of the inter-layer coupling decreases exponentially with inter-layer distance, already for a very small change in inter-layer distance, the coupling is diminished. Outside the GQB, however, the coupling still attains the value $\gamma_{1}^{0}$. As a result, to a very good approximation it is safe to assume that the two graphene layers are decoupled for the entire size of GQB, while just outside the blister the layers are coupled. Hence, we assume in this work that the inter-layer coupling profile has an abrupt transition at position R as shown in figure 1(c) by the black-dashed line. However, the results were shown to be robust against a smoothed blister [21].

In figures 1(d) and (e) we show the energy spectra, respectively, inside and outside the GQB under the application of a finite inter-layer bias δ. Because of the lack of inter-layer coupling the energy spectrum is linear and gapless inside the GQB. The application of a different potential to both layers, therefore, shifts the Dirac point in energy. As a result, for every energy there are electron or hole states available. Outside the GQB, figure 1(e) shows that the situation is substantially different. Because here the inter-layer coupling is strong, the inter-layer bias δ opens up a gap in the energy spectrum. In this region, only evanescent states are allowed and, therefore, the energy spectrum inside the GQB will be discrete and the corresponding modes are confined. The energy range where confinement appears is given by the range $\left[-\delta_{\rm G}, \delta_{\rm G}\right]$, where $\delta_{\rm G}$ is related to the inter-layer bias as [54]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta_{\rm G} = \delta \left(1+4\frac{\delta^{2}}{\gamma_{1}^{2}}\right)^{-1/2}. \label{eq_energy_range} \nonumber \end{align} \tag{ 5 }$$

In the following section we calculate the wavefunctions and the respective energy spectrum of the confined states in a GQB with an abrupt interface. For this, we first need to obtain the wavefunctions inside and outside the GQB and require continuity for each of the spinor components at the interface $\rho = R$ to find the energy eigenstates of the GQB in the presence of an inter-layer bias. In all calculations and results, the energy is scaled with the equilibrium inter-layer hopping parameter, $\gamma_1^{0}$, while $l = \hbar v_{\rm F} / \gamma_{1}^{0} \sim 1.65~{\rm nm}$ is the measure for the length scales.

In order to obtain the wavefunction outside the GQB, we solve the Schrödinger equation $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \hat{H}^{\tau}(\vec {r}) \Phi^{\tau}(\vec {r}) = E \Phi^{\tau}(\vec {r})$ for the Hamiltonian given in equation (1) with $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \gamma_{1}(\vec {r}) = \gamma_{1}^{0}$. The equation for the angle ϕ directly yields a relation between the phases of each spinor component. This means that the four-component wave function $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \Phi^{\tau}(\vec {r})$ in the τ valley can be written as [29]

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \newcommand{\re}{{\rm Re}} \displaystyle \label{wave_spinor_angular} \mathbf{\Phi^{\tau}}(\vec {r})=\left( \begin{array}{@{}cccc@{}} \phi_{A1}^{\tau}(\rho){\rm e}^{{\rm i}m\phi} \nonumber \\ {\rm i}\phi_{B1}^{\tau}(\rho){\rm e}^{{\rm i}(m-\tau)\phi} \nonumber \\ \phi_{B2}^{\tau}(\rho){\rm e}^{{\rm i}m\phi} \nonumber \\ {\rm i}\phi_{A2}^{\tau}(\rho){\rm e}^{{\rm i}(m+\tau)\phi} \end{array} \right). \nonumber \end{align} \tag{ 6 }$$

Solving the Schrödinger equation for the radial functions $\phi^{\tau}_{i}(\rho)$, we obtain the following set of coupled equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}-\frac{(\tau m-1)}{\rho}\right]\phi_{B1}^{\tau}=(E-\tau\delta)\phi_{A1}^{\tau}-\phi_{B2}^{\tau},\label{eq-06} \nonumber \end{align} \tag{ 7a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}+\frac{\tau m}{\rho}\right]\phi_{A1}^{\tau}=-(E-\tau\delta)\phi_{B1}^{\tau},\label{eq-07} \nonumber \end{align} \tag{ 7b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}+\frac{(\tau m+1)}{\rho}\right]\phi_{A2}^{\tau}=(E+\tau\delta)\phi_{B2}^{\tau}-\phi_{A1}^{\tau},\label{eq-08} \nonumber \end{align} \tag{ 7c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}-\frac{\tau m}{\rho}\right]\phi_{B2}^{\tau}=-(E+\tau\delta)\phi_{A2}^{\tau}.\label{eq-09} \nonumber \end{align} \tag{ 7d }$$

We remind the reader that in this set of equations, the energetic quantities are scaled with $\gamma_{1}^{0}$ and the radial component ρ by l, yielding dimensionless equations. The set of first-order differential equations can be written as a single fourth-order differential equation. As explained previously [55], this fourth-order differential equation has two sets of orthogonal solutions given by the solutions of the following second-order differential equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-10} \left[\frac{{\rm d}^{2}}{{\rm d}\rho^{2}}+\frac{1}{\rho}\frac{\rm d}{{\rm d}\rho}-\left( \frac{m^{2}}{\rho^{2}}+\alpha_{\pm}^{2} \right)\right]\phi_{A1}^{\tau}(\rho)=0. \nonumber \end{align} \tag{ 8 }$$

The two equations only differ by the value of

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_{\pm}^{2}=-(E^{2}+\delta^{2})\pm\sqrt{(E^{2}-\delta^{2})+4E^{2}\delta^{2}}. \label{alpha_equation} \nonumber \end{align} \tag{ 9 }$$

In the energy range where confinement is expected, the square root of equation (9) is imaginary. As a consequence, the solutions to equation (8) are Bessel functions with a complex argument [56]. Because we are outside of the GQB, the spinor components need to be finite in the limit $\rho\rightarrow \infty$, so we choose the modified Bessel function of the second kind $K_{m}(\alpha_{\pm} \rho)$ as solutions. Finally, notice that $\alpha_{+} = \alpha_{-}^{*}$, such that the two independent solutions of equation (8) can be written as a superposition of the real and imaginary part of $K_{m}(\alpha_{\pm} \rho)$, and we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-11} \phi_{A1}^{\tau}(\rho)=C_{1}^{\tau}\Re \left[ K_{m}(\alpha_{+} \rho) \right]+C_{2}^{\tau}\Im \left[ K_{m}(\alpha_{-} \rho) \right]. \nonumber \end{align} \tag{ 10a }$$

Using equations (7a)–(7d) we can obtain the other components explicitly as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-12} \phi_{B1}^{\tau}(\rho)=&~\frac{1}{(E-\tau \delta)}\left(C_{1}^{\tau}\Re \left[\alpha_{+} K_{m-\tau}(\alpha_{+} \rho) \right] \right. \nonumber \\ & \left. +C_{2}^{\tau}\Im \left[\alpha_{-} K_{m-\tau}(\alpha_{-} \rho) \right] \right), \nonumber \end{align} \tag{ 10b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-13} \phi_{B2}^{\tau}(\rho)=\frac{1}{(E-\tau \delta)}\left(C_{1}^{\tau}\Re \left[\eta^{+} K_{m}(\alpha_{+} \rho) \right]+C_{2}^{\tau}\Im \left[\eta^{-} K_{m}(\alpha_{-} \rho) \right] \right), \nonumber \end{align} \tag{ 10c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-14} \phi_{A2}^{\tau}(\rho)=&\frac{1}{(E^{2}- \delta^{2})}\left(C_{1}^{\tau}\Re \left[\eta^{+}\alpha_{+} K_{m+\tau}(\alpha_{+} \rho) \right] \right. \nonumber \\ & \left. + C_{2}^{\tau}\Im \left[\eta^{-}\alpha_{-} K_{m+\tau}(\alpha_{-} \rho) \right] \right). \nonumber \end{align} \tag{ 10d }$$

In these equations, we have introduced the compact notation $\newcommand{\e}{{\rm e}} \eta^{\pm}=\alpha_{\pm}^{2}+(E-\tau \delta){}^{2}$.

Inside the GQB the inter-layer coupling vanishes and, therefore, in equation (1) we have to put $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \gamma_{1}(\vec {r})=0$. Although the angular solution of the Schrödinger equation stays the same as in equation (6), the set of radial equations changes to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}-\frac{(\tau m-1)}{\rho}\right]\phi_{B1}^{\tau}=(E-\tau\delta)\phi_{A1}^{\tau}\label{eq-16}, \nonumber \end{align} \tag{ 11a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}+\frac{\tau m}{\rho}\right]\phi_{A1}^{\tau}=-(E-\tau\delta)\phi_{B1}^{\tau}\label{eq-17}, \nonumber \end{align} \tag{ 11b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}+\frac{(\tau m+1)}{\rho}\right]\phi_{A2}^{\tau}=(E+\tau\delta)\phi_{B2}^{\tau}\label{eq-18}, \nonumber \end{align} \tag{ 11c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\frac{\rm d}{{\rm d}\rho}-\frac{\tau m}{\rho}\right]\phi_{B2}^{\tau}=-(E+\tau\delta)\phi_{A2}^{\tau}\label{eq-19}. \nonumber \end{align} \tag{ 11d }$$

In this case, the set of equations is already decoupled for each layer. This allows to find a second-order differential equation for each layer as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-20} \left[\frac{{\rm d}^{2}}{{\rm d}\rho^{2}}+\frac{1}{\rho}\frac{\rm d}{{\rm d}\rho}-\left(\frac{m^{2}}{\rho^{2}}-\mu_{\pm}^{2} \right)\right]\phi_{B2/A1}^{\tau}(\rho)=0. \nonumber \end{align} \tag{ 12 }$$

In equation (12), $\mu_{\pm}=E\pm \tau\delta$ and the subscript of the function $\phi_{i}^{\tau}(\rho)$ refers to B1 for $\mu_{+}$ and to A1 for $\mu_{-}$. The solutions of equation (12) are Bessel functions. Dropping the ones singular at the origin, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-21} \phi_{A1}^{\tau}(\rho)=D_{1}^{\tau} J_{m}(\mu_{-} \rho), \nonumber \end{align} \tag{ 13a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-23} \phi_{B2}^{\tau}(\rho)=D_{2}^{\tau} J_{m}(\mu_{+} \rho). \nonumber \end{align} \tag{ 13b }$$

The other two components can then be found from equations (11b) and (11d) and yield

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-22} \phi_{B1}^{\tau}(\rho)=-\tau D_{1}^{\tau} J_{m-\tau}(\mu_{-} \rho), \nonumber \end{align} \tag{ 14a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-24} \phi_{A2}^{\tau}(\rho)=\tau D_{2}^{\tau} J_{m+\tau}(\mu_{+} \rho). \nonumber \end{align} \tag{ 14b }$$

We are now in a position to find the eigenstates and energylevels of a GQB. For this, we need to equate the spinor components inside the GQB with those outside at position $\rho = R$. Doing so, one obtains a set of four equations for four unknowns that can be written in a matrix formalism as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq-05} \mathbf{M^{\tau}}\left(\begin{array}{@{}cccc@{}} D_{1}^{\tau} \nonumber \\ D_{2}^{\tau} \nonumber \\ C_{1}^{\tau} \nonumber \\ C_{2}^{\tau} \end{array}\!\! \right)=\left(\begin{array}{@{}cccc@{}} -J_{m}(R\mu_{-}) & 0 & \Re[K_{m}(R\alpha_{+})] & \Im[K_{m}(R\alpha_{-})] \nonumber \\ \tau J_{m-\tau}(R\mu_{-}) & 0 & \Re[b_+K_{m-\tau}(R\alpha_{+})] & \Im[b_{-}K_{m-\tau}(R\alpha_{-})] \nonumber \\ 0 & -J_{m}(R\mu_{+})& \Re[c_+K_{m}(R\alpha_{+})]& \Im[c_{-}K_{m}(R\alpha_{-})] \nonumber \\ 0 & -\tau J_{m+\tau}(R\mu_{+})& \Re[d_+K_{m+\tau}(R\alpha_{+})]& \Im[d_{-}K_{m+\tau}(R\alpha_{-})] \nonumber \\ \end{array} \right)\left(\begin{array}{@{}cccc@{}} D_{1}^{\tau} \nonumber \\ D_{2}^{\tau} \nonumber \\ C_{1}^{\tau} \nonumber \\ C_{2}^{\tau} \end{array}\!\! \right)=0, \nonumber \end{align} \tag{ 15 }$$

where $b_{\pm}=\alpha_{\pm}/(E-\tau\delta)$, $c_{\pm}=\left[ (E-\tau \delta){}^2+\alpha_{\pm}^2 \right]/(E-\tau\delta)$, and $d_{\pm}=\alpha_{\pm}\left[ (E-\tau \delta){}^2+\alpha_{\pm}^2 \right]/(E^{2}-\delta^{2}).$ The energy levels E_m,n(R) of a GQB with radius R can be found through the roots of the determinant of the matrix $\mathbf{M^{\tau}}$. Here, n is the radial quantum number corresponding to $\vert n \vert$ modes in the radial direction that emerge with increasing the size of the blister. Subsequently, one can obtain the corresponding wavefunction by solving at the given energy and size R for the coefficients $C^{\tau}_{i}$ and $D^{\tau}_{i}$ and obtaining the eigenwavefunction $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \Phi^{\tau}_{m, n}(\vec {r})$. From this, the radial probability density (RPD) can be found as [57, 58]

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \newcommand{\re}{{\rm Re}} \displaystyle \mathcal{P}_{m,n}^{\tau}(\rho) = \rho \left\vert \Phi_{m,n}^{\tau}(\vec {r})\right\vert^{2}. \label{probability_density} \nonumber \end{align} \tag{ 16 }$$

Finally, the local density of states $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \mathcal{D}(\vec {r}, E)$ for a GQB with radius R can be derived from the eigenstates as

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\boldsymbol{#1}}} \newcommand{\re}{{\rm Re}} \displaystyle \mathcal{D}(\vec {r},E) = \sum_{m,n} \delta(E - E_{m,n}) \left\vert \Phi_{m,n}^{\tau}(\vec {r})\right\vert^{2}. \nonumber \end{align} \tag{ 17 }$$

In the numerical results displayed in the following section we will replace the Dirac function by a Gaussian profile with a finite spectral width Γ [47].

Now we turn to the discussion of the numerical results for the energy levels in three configurations of GQBs. Specifically, we consider the bias inside the blister to be the same as outside, opposite, and zero. The first configuration is the simplest example, i.e. the case where a homogeneous inter-layer bias potential δ is applied to the entire system. In figure 2 we show the energy levels as a function of the size R of the GQB for $m=0, \ldots, \pm 3$ with $\delta = 0.25~\gamma_{1}^{0}$. The results indicate that, indeed, for the energy range as defined in equation (5) the GQB has confined modes. Panel (a) of figure 2 shows that in the limit $R\rightarrow 0$, the GQB has two m  =  0 confined modes at energy $\pm \delta_{\rm G}$. As the radius of the GQB increases, the modes approach each other, anti-crossing one another around E  =  0. As the radius increases further, more m  =  0 modes are allowed inside the GQB. For a given radius, these modes are denoted by the radial quantum number n. The number of the confined modes crucially depends on the strength of the applied bias outside the blister $\delta_{>}$ and its size R. The energy spectrum of the different modes form anti-crossings with each other. As a consequence, the energy levels oscillate with the size of the GQB. These anti-crossings are a manifestation of coupling of states. This coupling is established through the coupled layers outside the blister, where the confined states on upper and lower layer can feel each other.

Zoom In Zoom Out Reset image size

In panels (b)–(d) of figure 2 we show the energy levels for non-zero angular quantum number m. These modes are only supported at larger radii R but their characteristic behavior is similar as the m  =  0 case shown in panel (a). Notice that the results are the same for positive as for negative m, i.e. $E_{m, n}=E_{-m, n}$. This is in contrast with previous studies where the symmetry between both signs of the angular quantum number is broken. [26, 28, 29]

To investigate the character of the different energy levels and the behavior at the anti-crossings, in figure 3 we show the RPD for m  =  0 at different configurations as indicated by the points in figure 2(a). We choose these points to be exactly at an anti-crossing, before, and after, and also near the continuum spectrum. Furthermore, we show the contribution of each sublattice to the probability density. Comparing for instance figures 3(a) and (b), we see that inside the GQB mainly holes reside on layer 1 while electrons reside on layer 2. In addition, we infer from figures 3(c) and (d) that at the anti-crossing electrons and holes are equally distributed on both layers. A transition in the residence of states takes place when passing an anti-crossing. The states with negative energy mainly reside on layer 2 instead of layer 1 before the anti-crossing point and vice versa for the positive one as shown in figures 3(e) and (f). This means that the modes $E_{0, \pm1}$ anti-cross each other in figure 2(a) and correspond mainly to states on one of the two layers before or after an anti-crossing. Peculiarly, however, in this configuration we find that states close to the continuum spectrum have also a large part of the probability located outside the GQB, see figures 3(g) and (h).

Zoom In Zoom Out Reset image size

We can find the radii R_j at which anti-crossing occur through the following partitioning relation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{condition_anticrossing_point} \int_{0}^{R_{j}}{\rm d}\rho\left[ \mathcal{P}_{m,n}^{\tau}(\rho) \right]_{\rm Layer~1}=\int_{0}^{R_{j}}{\rm d}\rho\left[ \mathcal{P}_{m,n}^{\tau}(\rho) \right]_{\rm Layer~2}, \nonumber \end{align} \tag{ 18a }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \mathcal{P}_{m,n}^{\tau}(\rho) \right]_{{\rm Layer}\ i}=\rho\left(\left\vert \phi^{m,n,\tau}_{Ai}(\rho) \right\vert^{2} +\left\vert \phi^{m,n,\tau}_{Bi}(\rho) \right\vert^{2}\right) \nonumber \end{align} \tag{ 18b }$$

where the radial part of the wave functions is normalized according to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Normalization} 2\pi\int_{0}^{\infty}{\rm d}\rho\sum_{i=1,2}\left[ \mathcal{P}_{m,n}^{\tau}(\rho) \right]_{{\rm Layer}\ i}=1. \nonumber \end{align} \tag{ 18c }$$

In other words, two energy levels form an anti-crossing when the probability of finding a state inside the GQB is the same for both layers. Note that a point in the middle between two subsequent anti-crossings associated with modes E_m,n also satisfies equation (18a). This point coincides with the anti-crossing in the second pair of energy branches, i.e. $E_{m, \left\vert n \right\vert +1}$, in the spectrum as can be inferred from figure 2(a).

Using equation (18a), we can find the radii R_j where the anti-crossings occur for any pair of energy branches E_m,n. For example, the first three anti-crossings of the first pair of energy branches $E_{0, \pm1}$ in figure 2(a) are located at R_j  =  (8.30, 20.74, 33.27)l. In between these three anti-crossings there are two points, where both layers also contribute the same to the RPD inside the GQB, located at R_j  =  (14.46, 26.98)l. For large GQB we notice that the strength of the anti-crossings becomes weaker. This is a result of leaking interaction between the two layers through the BLG outside the blister. In figure 4 we show the RPD outside the GQB at the first three anti-crossings labelled by the yellow dots (4, 9, 10) in figure 2(a). We see that the interaction between states on both layers becomes smaller with increasing the radius of GQB. Hence, in the limit $R\rightarrow \infty$ we expect the RPD to be zero outside the GQB and as result the anti-crossings will vanish and the states will be completely localized inside the blister. In this case, the GQB can be seen as a biased bilayer graphene nano-disk. We superimpose the first energy levels of a biased bilayer nano-disk with the respective angular momentum as gray dashed curves on figure 2. For a bilayer nano-disk we implement hard wall boundary conditions and the energy levels can be found by solving

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{m}(\mu_{+} R)=0, \ J_{m}(\mu_{-} R)=0. \nonumber \end{align} \tag{ 19 }$$

Zoom In Zoom Out Reset image size

The asymptotic behaviors of bessel function for small and large argument are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{asympt_Bessel0} J_{0}(x)=\left\{\begin{array}{@{}ll@{}}\sqrt{2/\pi x} \cos\left(x-\pi/4 \right)& x\rightarrow \infty \nonumber \\1-x^{2}/4 &x\rightarrow 0 \end{array}\right., \nonumber \end{align} \tag{ 20 }$$

using equation (20) one can show that $E\sim1/R$ for small and large size of bilayer nano-disk.

For $m\neq0$, we investigate the contribution of the two layers to the probability density for only $m=\pm1$ as shown in figure 5, and the findings also apply for $\left\vert m \right\vert>1$. We choose two points before the first anti-crossing of the modes $E_{1, \pm1}$ marked by red and green dots in figure 2(b). Because of the symmetry between m and  −m in this case ($E_{m, n}=E_{-m, n}$), these two points correspond to four modes as indicated in figure 5. We see that the states $m=\pm1$ with negative energy (green dot) mainly reside on the lower layer and vice versa for states with positive energy. This corresponds to the case with m  =  0 and it also holds here for the modes $E_{\pm1, \left\vert n \right\vert>1}$. Similar to the spectrum of m  =  0, the radii for which the energy forms anti-crossings can be also found using equation (18a). Notice that from the top panel of figure 5 the RPD acquires the layer symmetry

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \mathcal{P}_{m,n}^{\tau}(E,\rho) \right]_{\rm Layer~1} =\left[ \mathcal{P}_{-m,n}^{\tau}(-E,\rho) \right]_{\rm Layer~2}. \nonumber \end{align} \tag{ 21a }$$

Zoom In Zoom Out Reset image size

In the bottom panel of figure 5, we show the same results as in the top panel but in the vicinity of the $K'-$valley. Comparing top and bottom panels of figure 5, we find that the RPD also attains the following symmetry

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \mathcal{P}_{m,n}^{\tau}(E,\rho) \right]_{\rm Layer~1} =\left[ \mathcal{P}_{m,n}^{-\tau}(E,\rho) \right]_{\rm Layer~2}, \nonumber \end{align} \tag{ 21b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \mathcal{P}_{m,n}^{\tau}(E,\rho) \right]_{{\rm Layer}\ i} =\left[ \mathcal{P}_{-m,n}^{-\tau}(-E,\rho) \right]_{{\rm Layer}\ i}. \nonumber \end{align} \tag{ 21c }$$

Note that even though the RPD of each layer is different in each valley, the total RPD is the same in both valleys.

In figure 6 we show the local density of states for a GQB of fixed size $R=15~l$ as a function of the energy and distance from the origin for both layers. The results show that the layer selectivity of the modes is not only present for m  =  0, but also for the higher angular quantum numbers. Very pronounced is for example the $m=\pm 1$ mode that is strongly localized on layer 1 for negative energy and on layer 2 for positive energy. Such tunable layer localization was recently also observed on topological states in AB-BA domain walls in bilayer graphene [59]. Furthermore, the LDOS also shows that states with $\vert m \vert >1$ are not positioned at the center of the GQB, but more towards the edge or even outside the GQB in a classically forbidden region, specially, for those states close to the continuum spectrum. For example, the modes $\newcommand{\e}{{\rm e}} \epsilon_{1, \pm5}$ exist near the continuum spectrum of the AB-BLG and meanly localized outside the blister that is about R/3 far from the blister's edge. Note that the closed the modes to the continuum the far localized from the edge outside the blister.

Zoom In Zoom Out Reset image size

In the previous section we have considered the most straightforward case in which the inter-layer bias is the same in all parts of the sample. However, if the blister is formed by encapsulation of a metal colloid, the applied electric field also induces a dipole in the metallic nanoparticle. This will change the electrostatic potential on each layer. As a result, it can be strongly reduced inside the GQB with respect to outside it. To incorporate this difference, we investigate the case for which the inter-layer bias inside the GQB, $\delta_{<}$, is smaller than the bias $\delta_{>}$ outside.

In figure 7 we show the energy levels of a GQB with vanishing $\delta_{<}$ inside as a function of the radius of the blister for different values of the angular quantum number m. In this case, the energy levels do not show anti-crossings and approach each other as the size of the GQB increases in a monotonous way. In figure 7(a) we show the energy spectrum of the state with zero angular momentum. Here, contrary to the homogeneous bias case, each energy branch corresponds to states residing on a specific layer inside the GQB for any R. In figure 8 we show the RPD for different energy branches, labelled by yellow dots in figure 7(a). States on the first energy branch E_{0, −1} reside on the lower layer, see figure 8(a). While for the second branch E_{0, −2}, the states along it reside on the upper layer as shown in figure 8(b). Similarly, the third E_0,−3 and fourth E_{0, −4} branches, marked by points 3 and 4 in figure 7(a), the states reside mainly on the lower and upper layer respectively. This is illustrated in figures 8(c) and (d). Note that for the counterpart branches in the positive energy regime, the modes residence is opposite compared to the negative energy branches.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For non-zero values of the angular quantum number, however, the symmetry between positive and negative m is broken. This is most clear in the first pair of modes $E_{m>0, \pm1}$ in panels (b)–(d) in figure 7. The lowest of the two is only possible for negative m modes, while the highest is for positive m. These two modes, labelled by green and red dots in figures 7(b)–(d), are mainly localized at the interface of the blister as can be seen in figure 9. Moreover, inside the blister they mainly reside on the upper and lower layer for m  <  0 and m  >  0, respectively, as shown in figure 9. It turns out that they significantly reside on the disconnected sublattices A2 and B1 for the negative and positive angular momentum, respectively, as shown in figures 9(a)–(f). The appearance of these localized modes at the interface of the blister is one of its quintessential traits.

Zoom In Zoom Out Reset image size

Reminding ourselves that the lowest modes mainly reside on the upper layer while the upper modes reside on the lower layer, it follows that small GQBs only can support modes with a positive angular momentum on the lower layer while the negative angular momentum-modes reside on the upper layer. For large R, we notice that the broken symmetry between the negative and positive angular momentum is almost restored for modes whose radial quantum number $\left\vert n \right\vert >1$ such that $E_{m, n}(R) \approx E_{-m, n}(R)$ as shown in figures 7(b)–(d). The contribution of each layer to these modes is exactly the same as in the case of m  =  0 discussed in figure 7(a). For example, in the case of $m=\pm 1$, the first pair of modes labelled by the red circles in figure 7(b), mainly reside on the lower and upper layer for modes whose energy is negative and positive, respectively. The opposite occurs for the second pair of modes, labelled by the black circles, and such trend also holds for $\left\vert m \right\vert>1.$

In general, the energy levels still retain the following symmetry [55]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{m,n}(R) = -E_{-m,n}(R). \label{overall_symmetry} \nonumber \end{align} \tag{ 22 }$$

That this relation holds can be seen in figure 10 where the energy levels in a GQB with radius R  =  20l are plotted for a homogeneous inter-layer bias (panel (a)) compared with the case for a vanishing inter-layer bias (panel (b)) in the GQB. The results show that by changing the inter-layer bias inside the GQB, the m  =  0 modes are pushed away from each other while for the modes with finite angular momentum it even allows the total number of states in the GQB to be reduced.

Zoom In Zoom Out Reset image size

Note that considering inter-layer bias inside the blister that is smaller than outside and finite also allows confinement. The confined modes in this case still possess anti-crossings but at larger radius R. It also shows that for non-zero angular momentum the symmetry with respect to the sign of m remains strongly broken for small R, but that the modes with opposite m make a transition between different radial quantum number n.

Finally, in figure 11 we show the energy levels in the case when the bias inside the GQB is opposite to outside, i.e. $\delta_{<} = -\delta_{>}$. The m  =  0 result looks similar to the case of a homogeneous bias, but they are slightly different, i.e. $E_{0, n}(R, \delta_{<}) \neq E_{0, n}(R, -\delta_{<})$. For example, the first three anti-crossing points associated with the first pair of energy branches $E_{0, \pm1}$ are located at R_j  =  (10.93, 23.42, 35.96)l. The anti-crossings occur for slightly larger GQB in comparison with the homogeneous inter-layer bias case. Of particular importance is also the layer residency of states, where before the anti-crossing the states with negative energy mainly resided on the upper layer inside the GQB and vice versa for states with positive energy. This can be understood by considering the blister to be pure electron or hole doped [21]. Note that this is exactly the opposite to what happened in the case with homogeneous bias.

Zoom In Zoom Out Reset image size

For a finite angular momentum, the results in figure 11 show that the modes get pushed into each other, forming anti-crossings when the angular momentum quantum number is the same, while it crosses when m is opposite. As a result, for a given inter-layer bias, there is a radius R for which each of the non-zero m modes are degenerate at $E_{m, 1} = -E_{-m, 1} = 0$. We notice that the results of non-homogeneous bias also attains the symmetry $E_{m, n}(R) = -E_{-m, n}(R)$. In figure 12, we show the RPD of the modes $\newcommand{\e}{{\rm e}} \epsilon_{1, 1}$ and E_−1,1 at different radii indicated by the yellow dots in figure 11(b). It is evident from figures 12(a) and (b), that the degenerated modes $E_{\pm1, 1}=0$, labelled by point 1 in figure 11(b), are mainly localized at the interface of the blister. These modes are primarily confined on the upper and lower layer inside the GQB for positive and negative angular momentum, respectively, as can be seen from figures 12(a) and (b). While at points 2 and 3 the two layers contribute exactly the same to the confinement of the two states $E_{\pm1, 1}=\mp0.1 \gamma_{1}^{0}$ as shown in figures 12(c) and (d). Then, a transition occurs in the layer confinement, where at the points 4 and 5 the states $E_{\pm1, 1}=\mp0.06 \gamma_{1}^{0}$ principally reside on the upper and lower layer for positive and negative angular momentum, respectively. The behavior of the rest of modes $E_{\pm1, \left\vert n \right\vert>1}$ resembles that of the case m  =  0 in figure 11(a). This means that before and after anti-crossings, the states with negative and positive angular momentum reside on the upper and lower layer, respectively, while at the anti-crossings they are equally distributed. Analogously for $\left\vert m \right\vert>1$, we find that modes behave similarly to the case of $m=\pm1$.

Zoom In Zoom Out Reset image size

The inter-layer coupling $\gamma_{1}$ inside the GQB decreases very fast as the height of the blister increases However, if the blister is small and the layers in the blister remain loosely coupled, it is expected to show a band gap that is much smaller than outside the blister. Therefore, one also expects to find confined modes in this case. In this section, we investigate the energy levels of a non-zero inter-layer coupling in a GQB. In general, for a fixed gap outside the blister, the number of anti-crossings and their location mainly depend on the bias inside the blister. This allows to control the confinement to be mainly localized on a specific layer.

In figure 13 we show the energy levels of a GQB with a homogeneous inter-layer bias but different inter-layer coupling inside the blister as a function of the radius R. We see that for small radii, the energy levels are similar to the case of a completely decoupled blister, but that as the GQB grows, the energy levels do not cross the gap formed inside the GQB. As a consequence, for large R the oscillations of the lowest positive energy level are decreased and this level approaches the value of $\delta_{\rm G}$, dashed yellow lines, from equation (5) calculated with the inter-layer coupling inside the blister.

Zoom In Zoom Out Reset image size

Finally, the results in the current paper can be elegantly verified by measuring the conductance of electrons through the GQB. According to a recent experiment [34], a quantum dot with the size of few nanometers was realized and its electronic spectrum was probed by STM measurements. Thus, the same approach can be used to probe the electronic spectrum of the GQB. The ideal setup for this purpose is to apply a homogeneous bias to the GQB of strength that allows few modes to be confined within few nanometers size of the blister. The regime where the electronic properties of the GQB can be amenable to STM measurements restricted to a global bias of strength in the range $\backsim$$(30{{\rm {-}}}100)$ meV with the size R of the order $\backsim$$(15{{\rm {-}}}40)$ nm. Note that even the bias is not completely homogeneous or the inter layer coupling inside the blister is not strictly zero, confined states still exist.