Bound states in the continuum: Localization of Dirac-like fermions

The electronic behavior of graphene has been the focus of a great amount of work since its isolation in 2004 [1]. Graphene has an outstanding mobility at room temperature, being an excellent metallic material at the nanoscale. Since its charge carriers behave as chiral massless Dirac fermions, it is not possible to confine them by purely electrostatic means, so other strategies, such as employing bilayer graphene [2], have to be considered for nanoelectronic applications, for which gapped materials are needed to fabricate transistors and logic gates.

The difficulty to confine charge carriers in graphene has a fundamental origin: Klein tunneling was predicted to occur in monolayer graphene [3], and it was subsequently detected in transport experiments [4,5]. Devising alternative ways to confine massless fermions is therefore a relevant topic in graphene research. One possible way to attain localized states in graphene is by means of bound states in the continuum (BICs). The remarkable advances in graphene experiments [6] in which their electronic properties are tuned by applied external potentials, suggest that BICs could be observed in these structures.

Bound states in the continuum were first predicted by von Neumann and Wigner in 1929, who constructed a potential that yielded a truly localized, square-integrable state completely embedded in a continuum [7]. Since then, a great number of theoretical works have explored the feasibility of BICs in diverse setups. In fact, they have been predicted to occur in atomic and molecular systems [8–10], as well as in mesoscopic structures [11,12]. Systems with quantum dots can be exploited to produce BICs in low-dimensional structures [13–17]. The formation of BICs is a result of the interference of resonant states via the continuum. Several mechanisms can give rise to these bound states: they can appear due to symmetry effects [17,18], or as a result of a destructive interference process of resonant states for certain values of the physical parameters of the system [19–21], or finally based in the Fabry-Perot interferometer, of application to photonic systems [22].

As the physics of BICs relies on interference, it is not limited to quantum systems. In fact, in ref. [23] ballistic transport through a quantum dot was studied, demonstrating the possibility of a classical analogue of BICs. Furthermore, exploiting the analogy between electronics and photonics, several proposals have predicted BICs in photonic materials [24–26]. And time-dependent fields can be employed to achieve such localized states, as has been recently proposed [27,28]. On the other hand, there are only two experimental works reporting the measurement of BICs; remarkably, it has been achieved in a photonic system [29] and in a semiconductor quantum well system [30]. The observation by Capasso et al. [31] of an electronic state above the barrier of a semiconductor heterostructure was considered by other authors to be a BIC, but in fact it turned out to lie in a minigap of the superlattice [29]. Therefore, the search for electronic systems which could reveal the existence of robust BICs with unambiguous features is an important field of research. As discussed above, one possibility for detection of BIC is to employ graphene-based devices.

In this work we show that BICs can be detected in graphene trilayer systems. Our findings indicate that this nanostructure is an excellent candidate to observe bound states in the continuum by combining the measurement of the local density of states (LDOS) and the conductance through the system as a function of the gate voltage. This gate potential produces the coupling of the BIC to the continuum, so the LDOS at the BIC energy will be reduced for increasing values of the gate potential; however, we find that this coupling yields a maximum of the transmission through the system that persists for large gate voltages. Thus, the combination of LDOS and transport measurements provides a way to identify BICs in electronic systems.

We employ a π-orbital tight-binding model, which gives an excellent description of the electronic properties of graphene systems around the Fermi level. The in-plane nearest-neighbor interaction given by a single hopping parameter γ, which we take as $-3\ \text{eV}$. One interlayer hopping parameter between atoms directly placed on top of each other couples the flakes, being $\gamma^{\prime}=0.1 \gamma$.

We consider a trilayer armchair flake with direct AAA stacking. The leads are also armchair nanoribbons connected to the central flake, as shown in fig. 1. This structure can be alternatively viewed as an infinite armchair nanoribbon with two flakes placed symmetrically above and below it. Note that, although Bernal (AB) stacking is more stable for graphite, the direct or AA stacking has been experimentally found in few-layer graphene [32]. We choose the direct stacking because the AAA-stacked system can be easily mapped into a one-dimensional chain with two hoppings, yielding analytic expressions of the transmission that perfectly fit the conductance obtained numerically in the first transmission channel, as we shall later show in this paper. The widths considered correspond to metallic armchair nanoribbons, which play the role of contacts in this system. We give the length of the finite flake in translational unit cells, which corresponds to the number of 4-atom units along it. For the width we use the standard notation, giving the number of dimer chains across the ribbon.

Zoom In Zoom Out Reset image size

The total Hamiltonian is described by the sum of each graphene flake (up, central and down) Hamiltonian $H_{f}^{\alpha}$, with ($\alpha=u, c$ and d), the interaction between the flakes β ($\beta=u$ and d), with the central flake $H_{\beta - \textit{cf}}$, the interaction between central flake with the contacts $H_{\textit{cf}-\textit{leads}}$ and the Hamiltonian of the contacts $H_{\textit{leads}}$:

$$\begin{equation} H=\sum_{\alpha=u,c,d}H_{f}^{\alpha}+\sum_{\beta=u,d} H_{\beta- \textit{cf}}+ H_{\textit{cf}- \textit{leads}}+H_{\textit{leads}}, \end{equation} \tag{ 1 }$$

with

$$\begin{eqnarray} H_{f}^{\alpha}&=\sum_{i =1}^{N}\varepsilon_{i}^{\alpha}c_{i}^{\alpha {\dag}} c_{i}^{\alpha}+\gamma\sum_{\substack{i,j=1\\ &&\langle i,j \rangle}}^{N}(c_{i}^{\alpha {\dag}}c_{j}^{\alpha}+\text{H.c.}), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} H_{\beta - \textit{cf}}&=\gamma^{\prime}\sum_{i=1}^N (c_{i}^{\beta {\dag}} c_{j}^{c}+ \text{H.c.}), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} H_{\textit{cf}-\textit{leads}}&=\gamma\left\{(c_{1}^{c {\dag}}c_{0}^{c}+c_{N}^{c {\dag}} c_{N+1}^{c})+\text{H.c.}\right\}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} H_{\textit{leads}}&=\sum_{\substack{i>N, \\ &&i<1}}\varepsilon_{i}^{c}c_{i}^{c{\dag}} c_{i}^{c}+\gamma\sum_{\substack{i,j>N,\\ &&i,j<1\\ &&\langle i,j \rangle}}(c_{i}^{c {\dag}} c_{j}^{c}+\text{H.c.}), \end{eqnarray} \tag{ 5 }$$

where $\varepsilon_{i}^{\alpha}$ is the site energy for atom i in layer α; $c_{i}^{\alpha}\ (c_{i}^{\alpha {\dag}})$ is the annihilation (creation) operator of one electron in atom i of layer α; γ is the nearest-neighbor hopping between atoms inside a layer and $\gamma^{\prime}$ is the hopping between atoms directly on top of each other in neighboring layers. Note that as the leads are smoothly connected to the central flake constituting a nanoribbon, we leave the superindex c for lead variables and operators.

This problem is often solved numerically by standard Green function techniques, but it also admits an analytical approach. To this end, a graphene monolayer is mapped onto a one-dimensional (1D) chain [33], so the trilayer is modeled as three coupled 1D chains.

The stationary state for one monolayer graphene can be written as $|\psi\rangle= \sum_{j,m}^{N}\left(\varphi_{j,m}^{A}|j,m\rangle^{A}+\varphi_{j,m}^{B}|j,m\rangle^{B} \right)$, where $\varphi_{j,m}^{A}$ and $\varphi_{j,m}^{B}$ represent the probability amplitudes to find one electron in the dimer $(j, m)$ in atoms A and B respectively, fig. 2. From the eigenvalue equation for one monolayer graphene, $H_{\textit{layer}}^{\alpha}|\psi\rangle=E|\psi\rangle$, two linear difference equations are obtained [33]:

$$\begin{eqnarray} \varepsilon\Psi_{j,m}^{A}&=&\gamma(\Psi_{j-1,m-1}^{B}+\Psi_{j-1,m+1}^{B}+\Psi_{j,m}^{B}), \\ \varepsilon\Psi_{j,m}^{B}&=&\gamma(\Psi_{j+1,m-1}^{A}+\Psi_{j+1,m+1}^{A}+\Psi_{j,m}^{A}). \end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

Considering plane wave solutions with q transverse wave number, $\varphi_{j,m}^{A,B}=e^{iqm}\phi_{j}^{A,B}$, the equations of motion are

$$\begin{eqnarray} &&E\phi_{j}^{A}=\gamma(2\cos{q}\;\phi_{j-1}^{B}+\phi_{j}^{B}), \\ &&E\phi_{j}^{B}=\gamma(2\cos{q}\;\phi_{j+1}^{A}+\phi_{j}^{A}). \end{eqnarray} \tag{ 7 }$$

This is equivalent to a one-dimensional chain with different hoppings between A-B and B-A atoms, namely, γ and $\eta=2 \gamma \cos q$. Incidentally, note that close to the Dirac point K, q takes the value $2\pi /3$, so the equations are written as

$$\begin{eqnarray} &&E\phi_{j}^{A}=\gamma(-\phi_{j-1}^{B}+\phi_{j}^{B}), \\ &&E\phi_{j}^{B}=\gamma(-\phi_{j+1}^{A}+\phi_{j}^{A}). \end{eqnarray} \tag{ 8 }$$

Thus, from the tight-binding equations of motion near the K point, we arrive in the continuum limit at the Dirac equation, $E\tilde{\Phi}(x) = v_F\sigma_2\partial_x \tilde{\Phi}(x)$, where $v_F=\gamma a/\hbar$ is the Fermi velocity and a is the lattice constant of graphene.

Now we consider the trilayer system as composed of three 1D chains with interchain hopping $\gamma^{\prime}$, as depicted in fig. 3. The equations of motion for trilayer graphene can be written in a basis that decouples the transversal modes. The Hamiltonian representing this system is given by

$$\begin{equation} H=\left(\begin{array}{c@{\quad}c@{\quad}c} \varepsilon_{1\text{D}}+\Delta & \gamma^{\prime} & 0 \\ \gamma^{\prime} & \varepsilon_{1\text{D}} & \gamma^{\prime} \\ 0 & \gamma^{\prime} & \varepsilon_{1\text{D}}-\Delta \end{array} \right), \end{equation} \tag{ 9 }$$

where $\varepsilon_{1\text{D}}$ is the energy of the one-dimensional chain, and Δ is the gate potential. The above Hamiltonian is written in the basis $\Psi=\left(\Psi^{u},\Psi^{c},\Psi^{d}\right)$, where $\Psi^{u}, \Psi^{c}$ and $\Psi^{d}$ represent the atomic wave functions of the chain up, central and down, respectively.

Zoom In Zoom Out Reset image size

Solving the eigenvalue equation $H|\Psi\rangle=\varepsilon|\Psi\rangle$, we obtain

$$\begin{eqnarray} \overline{\varepsilon}_{1}&=&\varepsilon_{1\text{D}}+\sqrt{\Delta^{2}+ 2\gamma^{\prime^{2}}}, \nonumber\\ \overline{\varepsilon}_{2}&=&\varepsilon_{1\text{D}}, \\ \overline{\varepsilon}_{3}&=&\varepsilon_{1\text{D}}-\sqrt{\Delta^{2}+ 2\gamma^{\prime^{2}}}. \nonumber \end{eqnarray} \tag{ 10 }$$

The energies $\overline{\varepsilon}_{1}, \overline{\varepsilon}_{2}\ \text{y}\ \overline{\varepsilon}_{3}$ are the renormalized energies for the atoms of the uncoupled chains up, central and down, respectively. The corresponding eigenvectors form a new basis which diagonalizes the Hamiltonian $\bar{H}=P^{-1}HP$, and therefore decouples the three chains, where the matrix P is

$$\begin{equation} P=\left(\begin{array}{c@{\quad}c@{\quad}c} \displaystyle\frac{\gamma^{\prime}}{\sqrt{\Delta^{2}+2\gamma^{\prime^{2}}}-\Delta}&1 &\displaystyle-\frac{\gamma^{\prime}}{\sqrt{\Delta^{2}+2\gamma^{\prime^{2}}}+\Delta} \\ 1 & -\displaystyle\frac{\Delta}{\gamma^{\prime}} & 1 \\ \displaystyle\frac{\gamma^{\prime}}{\sqrt{\Delta^{2}+2\gamma^{\prime^{2}}}+\Delta}&-1 &\displaystyle\frac{\gamma^{\prime}}{-\sqrt{\Delta^{2}+2\gamma^{\prime^{2}}}+\Delta} \end{array}\right), \end{equation} \tag{ 11 }$$

so that

$$\begin{equation} \bar{H}=\left(\begin{array}{c@{\quad}c@{\quad}c} \overline{\varepsilon}_{1} & 0 & 0 \\ 0 & \overline{\varepsilon}_{2} & 0 \\ 0 & 0 & \overline{\varepsilon}_{3} \end{array}\right). \end{equation} \tag{ 12 }$$

With the above consideration, it is possible to calculate the scattering amplitude through the trilayer by applying the standard boundary conditions for the resulting spinors in the three regions of the system, namely, left lead, trilayer flake, and right lead. Such a procedure leads to the following transmission probability:

$$\begin{equation} T=|t|^{2}=\frac{2\mu^{2}F^{2}\left(1-\cos{2k_{1\text{D}}}\right)}{\mu^{4}H^{2}+ 4\mu^{2}G^{2}-4\mu G\cos{k_{1\text{D}}}+2\mu^{2}H\left(2\mu G\cos{k_{1\text{D}}}-\cos{2k_{1\text{D}}}\right)+1}, \end{equation} \tag{ 13 }$$

here $\mu=\gamma^{\prime}/(\Delta^{2}+2\gamma^{\prime^{2}})$, where Δ is the gate potential, $k_{1\text{D}}$ is the wave vector of the 1D chain and $H=F^2-G^2$ with F, G given by

$$\begin{eqnarray*} F&=& \frac{\gamma^{\prime}}{U_{N}(q_{1})}-\frac{\Delta\nu}{U_{N}(q_{2})}+ \frac{\gamma^{\prime}}{U_{N}(q_{3})},\\ G&=&\gamma^{\prime}\frac{U_{N-1}(q_{1})}{U_{N}(q_{1})}-\Delta\nu\frac{U_{N-1}(q_{2})} {U_{N}(q_{2})}+\gamma^{\prime}\frac{U_{N-1}(q_{3})}{U_{N}(q_{3})}, \end{eqnarray*}$$

with $\nu=-\Delta/\gamma^{\prime}$, $U_n(x)$ are the Chebyshev polynomials of the second kind, and $x=q_1,q_2, q_3$ are the allowed wave vectors in the renormalized chains composing the trilayer, given by

$$\begin{eqnarray} q_{1}&=&\arccos{\left[\frac{\varepsilon_{1\text{D}}-\sqrt{\Delta^{2}+2 \gamma^{\prime^{2}}}}{2}\right]}, \nonumber\\ q_{2}&=&\arccos\left[\frac{\varepsilon_{1\text{D}}}{2}\right],\\ q_{3}&=&\arccos{\left[\frac{\varepsilon_{1\text{D}}+\sqrt{\Delta^{2}+2 \gamma^{\prime^{2}}}}{2}\right]}. \nonumber \end{eqnarray} \tag{ 14 }$$

Notice that these three wave vectors correspond to the three states of the trilayer, namely, the coupled bonding and antibonding solutions and the nonbonding state, which is independent of the interlayer coupling. This latter state gives rise to the BIC.

From the theoretical perspective, a BIC can be identified by a Dirac delta-like peak in the LDOS embedded in the continuum, being a square-integrable, localized state coexisting with nonnormalizable, propagating states. However, this does not constitute an unequivocal description from the experimental viewpoint, given that resonances or Fabry-Pérot states can also give rise to very sharp peaks in the density of states. We have studied the evolution of the LDOS and the conductance of these BIC states in graphene trilayers as a function of the gate voltage Δ, showing that they can be identified unambiguously considering both LDOS and transport measurements.

In order to check our results, we also obtain numerically the conductance and the LDOS for these systems employing a decimation technique to obtain the Green's function of the structure, as done previously for similar systems [34–36]. In the energy region with only one propagating state in the leads, the results agree perfectly. Within the coupled chains model, the LDOS is obtained as the sum of the squared moduli of the amplitudes over the trilayer flake.

In fig. 4 we present the transmission and LDOS for the N = 12 trilayer calculated with the coupled chain model for the symmetric case $(\Delta =0)$ and for a nonzero, albeit small, value of the gate potential, $\Delta=0.01 \gamma$. Energies are given in units of γ throughout the paper. The LDOS shows very sharp peaks superposed to broader peaks over a nonzero density for all energies. The sharp features correspond to BICs. The broader peaks correspond to antiresonances related to bonding and antibonding states of the trilayer, strongly coupled to the continuum. The shape of these superimposed peaks can be described by the sum of two Lorentzians with quite different line shapes: one corresponds to the strongly coupled state, and the other to the BIC, yielding a Dirac delta in the limit of zero gate voltage [17].

Zoom In Zoom Out Reset image size

The transmission probability, shown in the top panels of fig. 4, presents antiresonances at the energies of the strongly coupled states, but there is not any signature of BICs in the transport properties for $\Delta =0\gamma$. However, for $\Delta=0.01\gamma$ a sharp transmission peak appears for each BIC, reaching the maximum value T = 1. This change indicates that the BICs are beginning to couple to the continuum, being in fact quasi-bound states in the continuum (quasi-BICs). Due to this coupling, they contribute to the conductance of the system. The LDOS still shows the maxima at the energies of the quasi-BICs.

In fact, the gate potential plays the role of the asymmetry parameter, which controls the coupling of quasi-BICs to the continuum in this system. Remarkably, increasing the value of the gate potential does not destroy the characteristic maximum in the transmission. Figure 5 shows the LDOS and transmission probabilities for larger values of the gate voltage, $\Delta=0.05 \gamma$ and $\Delta=0.1\gamma$. The coupling of the quasi-BICs to the continuum is so large that they no longer show any peak in the LDOS; however, their presence can be inferred from the maxima in the transmission that persist for these voltages. Most notably, the maximum value of the transmission still reaches 1, while the other maxima appearing between resonances have a smaller value, dependent on the length of the flake, similarly to the results found for bilayer flakes [35,37]. Besides, the width of the maxima in the transmission corresponding to these quasi-BIC states increases with the value of the gate potential (the coupling parameter), thus yielding a robust signal that can be studied as a proof of the existence of BICs in trilayer graphene systems.

Zoom In Zoom Out Reset image size

When $\Delta\rightarrow 0$, the transmission around the energy of the quasi-BICs can be written as a superposition of Fano and Breit-Wigner line-shapes,

$$\begin{equation} T(\varepsilon)=\frac{(\varepsilon-\varepsilon_0)^2}{(\varepsilon-\varepsilon_0)^2+ \Gamma_{+}^2}+\frac{\Gamma_{-}^2} {(\varepsilon-\varepsilon_0)^2+\Gamma_{-}^2}, \end{equation} \tag{ 15 }$$

where $\varepsilon_0$ is the position of the quasi-BIC, $\Gamma_{+}=2\Gamma_0$ and $\Gamma_{-}=\Delta^2/2\Gamma_0$. The parameter $\Gamma_0$ is related to the interlayer coupling $\gamma^{\prime}$, so that T = 1 for zero interlayer coupling and it yields a peak of width $\propto \Delta^2$, as it can be appreciated in figs. 4 and 5. This dependence of the width of the transmission peak as a function of the gate voltage is another signature of the quasi-BIC states.

To show that the electronic behavior described above is robust under different conditions, we have numerically calculated the transmission probability for trilayer ribbons in ABA stacking and, for trilayer ribbons in AAA stacking with a small amount of disorder at their edges.

In the case of ABA trilayer system, fig. 6(a) and (b) show that the appearance of BICs in the conductance curves is independent of the stacking of the system, thus, for highly symmetric systems it is necessary to apply a gate potential to couple these bound states with the continuum spectra of the ribbon. As in the previous case, the conductance peaks reach the maximum of the transmission at the BICs energy levels.

Zoom In Zoom Out Reset image size

In the case of a defective AAA stacking, we have introduced a small amount of disorder to the trilayer flakes by removing few atoms at the ribbons edges. In this case, due to the loss of spatial symmetry of the system, quasi-BICs appear in the conductance curves even for $\Delta = 0$. As a gate voltage is applied to the system, quasi-BICs become broader and as a consequence, the associated conductance peaks reach the maximum value of the transmission.

The above behavior confirms that these interference phenomena are robust for different system configurations, even in the absence of symmetry. In the light of these results, we believe that the existence and visualization of BICs in the suggested systems is possible.

Although pure BICs do not give any contribution to the conductance when they are uncoupled, coupling them to the continuum by means of a gate voltage makes them contribute to the current with a maximum transmission, with a characteristic dependence on the external potential. The persistence of the maximum, even for large values of the applied voltage, enables their characterization by the analysis of the LDOS and electronic conductance. The analytic calculations presented in this letter are valid in the low-energy range (close to the Fermi level) for metallic nanoribbons and for pristine graphene. As the key point of our proposal is to detect these states by a combination of LDOS and transport measurements, metallic ribbons are the natural option. These states could also be detected in a symmetric stacking of metallic zigzag nanoribbons, although it cannot be treated analytically. On the other hand, in the literature it has been reported that BICs are robust under correlation effects [38,39], therefore, we expect that in this kind of symmetric configuration, independently of the electron correlation effects, BICs should exist and may appear in the conductance curves if the symmetry of the potentials is slightly broken by external manipulation.

This work has been partially supported by the Spanish DGES under grant FIS2012-33521, Chilean FONDECYT grants 1140571 (PO), 1140388 (LR), and DGIP/USM internal grant 11.14.68 (MP, PO), CONICYT ACT 1204 (LR, MP, PO). LC gratefully acknowledges the hospitality of the Universidad Técnica Federico Santa María (Chile).