Confinement in bilayer graphene via intra- and inter-layer interactions

In bilayer graphene, two flakes of graphene are close each other such that their electrons can mutually interact. The relative orientation of the two layers can vary. In case of Bernal (or AB-) stacking, the two layers are relatively shifted such that some bonds are parallel in the two lattices. Denoting the atoms in the two triangular sublattices of the jth layer as A_j and B_j , the atoms A₂ of the upper layer sit just above the B₁ atoms of the lower layer, whereas the atoms A₁ are below the centers of the hexagons of the upper lattice. The effective Hamiltonian (in the four-band approximation) for the low-energy particles can be written as [1]

$$\begin{equation}{H}_{blg}=\left(\begin{matrix}{{\epsilon}}_{{A}_{1}}& {\pi }^{{\dagger}}& 0& {v}_{3}\pi \\ \pi & {{\epsilon}}_{{B}_{1}}& \gamma & 0\\ 0& \gamma & {{\epsilon}}_{{A}_{2}}& {\pi }^{{\dagger}}\\ {v}_{3}{\pi }^{{\dagger}}& 0& \pi & {{\epsilon}}_{{B}_{2}}\end{matrix}\right),\quad \pi =-i{\partial }_{x}+{\partial }_{y}.\end{equation} \tag{ 1 }$$

The basis of the bispinors is (A₁, B₁, A₂, B₂). The quantities ${{\epsilon}}_{{A}_{j}({B}_{j})}$, j ∈ {1, 2}, correspond to the on-site energies that can originate from an external electric potential, spin–orbit interaction or an interaction with the substrate. The parameter γ represents interaction of the electrons on the sites A₂ and B₁. The term proportional to v₃ is related to the interlayer trigonal warping [1], which is frequently set to zero in the literature. The four-band Hamiltonian (1) was used e.g. in the analysis of strains and their effect on electronic [2–4] or topological properties of bilayer graphene [5]. It was used in the analysis of confinement of Dirac fermions in quantum dots formed by doping [6] or by local variation of the interlayer coupling related to local delamination of the bilayer graphene (graphene blisters) studied recently in [7, 8]. It serves well for description of other bilayer Dirac materials, see e.g. bilayer silicene [9–11]. Qualitatively the same operator with v₃ = 0 appears in description of spin–orbit interaction in graphene, see e.g. [12, 13].

Electrons on the binding sites A₂ and B₁ form dimers. For E ≪ γ, the dynamics on the non-dimer sites A₁ and B₂ gets dominant, and the effective Hamiltonian can be derived from (1), see [1, 14]. The latter case is known as the two-band approximation, in which the energy operator reads as $H=\left(\right.\enspace \enspace _{{({\pi }^{{\dagger}})}^{2}}^{0}\enspace \enspace {}_{0}{}^{{\pi }^{2}}\left.\right)$. This framework proved to be useful in the analysis of various situations where electric or magnetic fields are inhomogeneous, or the bilayer is subjected to inhomogeneous deformations [14–17]. Let us also mention that the exactly solvable two-band Hamiltonians with inhomogeneous magnetic field were constructed [18]. Coherent states for the system described by this operator were found in [19]. In most cases, the studied systems possessed translational or rotational symmetry so that effectively one-dimensional systems were analyzed.

In the current article, we are interested in systems where, besides the on-site interactions and (pseudo-)magnetic field, the inter-layer coupling γ can also be inhomogeneous. We are particularly interested in situations where fluctuations of the involved interactions can confine the Dirac fermions. In this quest, we prefer working within the framework of the four-band approximation. Therefore, we will consider the settings described by the following energy operator

$$\begin{equation}H=\left(\begin{matrix} {{\epsilon}}_{1}(x) & -i{\partial }_{x}+A(x) & 0 & i{v}_{3}{\partial }_{x} \\ -i{\partial }_{x}+{A}^{\ast }(x) & {{\epsilon}}_{2}(x) & \gamma (x) & 0 \\ 0 & \gamma (x) & {{\epsilon}}_{2}(x) & -i{\partial }_{x}+{A}^{\ast }(x) \\ i{v}_{3}{\partial }_{x} & 0 & -i{\partial }_{x}+A(x) & {{\epsilon}}_{1}(x) \end{matrix}\right).\end{equation} \tag{ 2 }$$

When compared to (1), this operator can be matched with the Hamiltonian of bilayer graphene for inhomogeneous ${{\epsilon}}_{{A}_{1}}={{\epsilon}}_{{B}_{2}}={{\epsilon}}_{1}(x)$, ${{\epsilon}}_{{B}_{1}}={{\epsilon}}_{{A}_{2}}={{\epsilon}}_{2}(x)$ and with the longitudinal momentum k_y = 0. Therefore, (2) can describe dynamics of the fermions that bounce on the potential in normal direction. The operator (2) also contains an additional potential term A_x σ₀ ⊗ σ₁ − A_y σ₃ ⊗ σ₂, A = A_x + iA_y . There exist numerical methods for solving the stationary equation associated with (2), such as the Fourier grid Hamiltonian [20] and its extension for matrix Hamiltonians [21]. They provide information in a finite region (local), and the intrinsic numerical error is carried on when computing other physically relevant quantities, such as expectation values. In this regard, analytic solutions prove to be more reliable as they provide global information, which becomes valid on the whole considered domain. Additionally, they are error-free by definition, which makes them outstanding candidates when it comes to implementing perturbative theories. Despite their advantages, it might be a challenging and sometimes impossible task to find them explicitly. Still, there are mathematical techniques that help to find the analytical solutions.

It was discussed recently in [22] that the Hamiltonian (2) belongs to the class of reducible operators where the solution of the associated stationary equation can be found via two, lower-dimensional, dynamical equations with auxiliary interactions. This substantially simplifies the determination of exact solutions. Indeed, let us make an ansatz for the eigenstates of (2) in the following manner,

$$\begin{equation}\mathbf{\Psi }={({\xi }_{1},{\xi }_{2},{\xi }_{2},{\xi }_{1})}^{\text{T}},\qquad \mathbf{\Xi }={({\chi }_{1},{\chi }_{2},-{\chi }_{2},-{\chi }_{1})}^{\text{T}}.\end{equation} \tag{ 3 }$$

Then, the bispinors Ψ and Ξ satisfy ³

$$\begin{equation}H\enspace \mathbf{\Psi }=E\enspace \mathbf{\Psi },\qquad H\enspace \mathbf{\Xi }=\bar{E}\enspace \mathbf{\Xi },\end{equation} \tag{ 4 }$$

provided that the spinors $\mathbf{\xi }={({\xi }_{1},{\xi }_{2})}^{\text{T}}$ and $\mathbf{\chi }={({\chi }_{1},{\chi }_{2})}^{\text{T}}$ are solutions of the following two equations,

$$\begin{equation}{h}_{1}\enspace \xi =\left[-i\left(\begin{matrix} {v}_{3} & 1 \\ 1 & 0 \end{matrix}\right){\partial }_{x}+\left(\begin{matrix} {{\epsilon}}_{1}(x) & A(x) \\ {A}^{\ast }(x) & {{\epsilon}}_{2}(x)+\gamma (x) \end{matrix}\right)\right]\xi =E\enspace \xi ,\end{equation} \tag{ 5 }$$

$$\begin{equation}{h}_{2}\enspace \chi =\left[-i\left(\begin{matrix} -{v}_{3} & 1 \\ 1 & 0 \end{matrix}\right){\partial }_{x}+\left(\begin{matrix} {{\epsilon}}_{1}(x) & A(x) \\ {A}^{\ast }(x) & {{\epsilon}}_{2}(x)-\gamma (x) \end{matrix}\right)\right]\chi =\bar{E}\enspace \chi ,\end{equation} \tag{ 6 }$$

with E and $\bar{E}$, in general, different from each other.

These equations resemble the 1 + 1 dimensional Dirac equation up to the matrix coefficient in the kinetic term that contains the nonvanishing constant v₃. Let us investigate how the equations (5) and (6) can be solved and what are the physical scenarios that can be studied in this way.

For the sake of generality, let us consider the matrix Hamiltonian

$$\begin{equation}h=-i\left(\begin{matrix} {v}_{3} & 1 \\ 1 & 0 \end{matrix}\right){\partial }_{x}+\left(\begin{matrix} {V}_{1}(x) & A(x) \\ {A}^{\ast }(x) & {V}_{2}(x) \end{matrix}\right)\enspace ,\quad {V}_{1}(x),{V}_{2}(x):\mathbb{R}\to \mathbb{R}\enspace ,\quad A(x):\mathbb{R}\to \mathbb{C}.\end{equation} \tag{ 7 }$$

The stationary equation

$$\begin{equation}h\left(\begin{matrix} {\psi }_{1} \\ {\psi }_{2} \end{matrix}\right)=E\left(\begin{matrix} {\psi }_{1} \\ {\psi }_{2} \end{matrix}\right),\end{equation} \tag{ 8 }$$

gives rise to the set of coupled equations for the spinor components ψ₁ and ψ₂,

$$\begin{equation}-i{v}_{3}{\psi }_{1}^{\prime }-i{\psi }_{2}^{\prime }+A(x){\psi }_{2}+({V}_{1}(x)-E){\psi }_{1}=0,\end{equation} \tag{ 9 }$$

$$\begin{equation}-i{\psi }_{1}^{\prime }+{A}^{\ast }(x){\psi }_{1}+({V}_{2}(x)-E){\psi }_{2}=0,\end{equation} \tag{ 10 }$$

with f' ≡ ∂f/∂_x . The equations (9) and (10) can be decoupled by fixing

$$\begin{equation}{\psi }_{2}=\frac{-i{\psi }_{1}^{\prime }+{A}^{\ast }(x){\psi }_{1}}{E-{V}_{2}(x)},\end{equation} \tag{ 11 }$$

so that (9) turns into a second-order differential equation for ψ₁. In order to bring it into the Sturm–Liouville form, we make an additional energy-dependent transformation

$$\begin{equation}{\psi }_{1}(x)={\text{e}}^{-\text{i}\phi (x)}{\tilde{\psi }}_{1}(x),\end{equation} \tag{ 12 }$$

where

$$\begin{equation}\phi (x)=\frac{i}{2}\enspace \mathrm{log}(E-{V}_{2}(x))+\frac{1}{2}{\int }^{x}\mathrm{d}{x}^{\prime }\left({v}_{3}(E-{V}_{2}(x))+2\enspace \mathrm{Re}\enspace A({x}^{\prime })\right).\end{equation} \tag{ 13 }$$

Notice that the first term of ϕ(x) is not a pure phase so that the normalization of ψ₁ and ${\tilde{\psi }}_{1}$ can differ. This should be kept in mind when imposing boundary conditions on the functions ${\tilde{\psi }}_{1}$ in general. The equation for ${\tilde{\psi }}_{1}$ reads as

$$\begin{align} -{\tilde{\psi }}_{1}^{{\prime\prime}}& +\left(-\mathrm{Im}\enspace {A}^{\prime }(x)+{(\mathrm{Im}\enspace A(x))}^{2}-{v}_{3}(E-{V}_{2}(x))\;\mathrm{Re}\enspace A(x)\right.\left.-(E-{V}_{1}(x))(E-{V}_{2}(x))\right. \\ & \left.-\frac{{v}_{3}^{2}}{4}{(E-{V}_{2}(x))}^{2}-\frac{\mathrm{Im}\enspace A\enspace {V}_{2}^{\prime }(x)}{E-{V}_{2}(x)}+\frac{3{({V}_{2}^{\prime }(x))}^{2}}{4{(E-{V}_{2}(x))}^{2}}+\frac{{V}_{2}^{{\prime\prime}}(x)}{2(E-{V}_{2}(x))}\right){\tilde{\psi }}_{1}=0. \end{align} \tag{ 14 }$$

Therefore, ${\tilde{\psi }}_{1}$ can be found as the zero-mode of the Schrödinger equation with a potential given in terms of A(x), V₁(x), and V₂(x). Although it is a rather complicated task to find its solutions in general, one can identify (14) with a Schrödinger equation that possesses a solvable potential V₀, $-{\tilde{\psi }}_{1}^{{\prime\prime}}+{V}_{0}{\tilde{\psi }}_{1}=0$. However, by doing so, one has to solve a nonlinear differential equation for either Im A or V₂. To overcome this issue, one can fix either Re A(x) or V₁(x) so that the potential in (14) coincides with V₀. This way, one of the two quantities can be fixed as

$$\begin{equation}{V}_{1}={V}_{1}({V}_{2},A,{V}_{0},E)\quad \text{or}\quad \mathrm{Re}\;A=\mathrm{Re}\enspace A({V}_{2},{V}_{1},\mathrm{Im}\enspace A,{V}_{0},E).\end{equation} \tag{ 15 }$$

Then we can acquire a solution ${\tilde{\psi }}_{1}$ for a single energy level E. Indeed, changing E would alter the interaction V₁(x) so that we would deal with a different physical setting. We will discuss this situation in the section 4.

For constant V₂, equation (14) simplifies considerably as the third line in (14) cancels out, leading to

$$\begin{equation}-{\tilde{\psi }}_{1}^{{\prime\prime}}+{V}_{E}(x){\tilde{\psi }}_{1}={\lambda }_{E}{\tilde{\psi }}_{1},\end{equation} \tag{ 16 }$$

where

$$\begin{align} {V}_{E}(x)& =-\mathrm{Im}\enspace {A}^{\prime }(x)+{(\mathrm{Im}\enspace A(x))}^{2}-{v}_{3}(E-{V}_{2})\;\mathrm{Re}\enspace A(x) \\ & \quad -(E-{V}_{2})(E-{V}_{1}(x)), \end{align} \tag{ 17 }$$

$$\begin{equation}{\lambda }_{E}=\frac{{v}_{3}^{2}}{4}{(E-{V}_{2})}^{2}.\end{equation} \tag{ 18 }$$

The equation (16) corresponds to a stationary Schrödinger equation with an energy-dependent potential term V_E (x). The energy-dependent part in (17) cancels out effectively when V₁ is constant and either v₃ vanishes or Re A is a constant. In the latter case, the term v₃(E − V₂)A₀ + (E − V₁)(E − V₂) adds a shift to the eigenvalue λ_E .

Notice that when E = V₂, the equations (9) and (10) can be solved directly with

$$\begin{equation}{\psi }_{1}={\text{e}}^{-\text{i}\int {A}^{\ast }(x)\mathrm{d}x},\qquad {\psi }_{2}=-\frac{i\int {\text{e}}^{\text{i}\int (A(x)-{A}^{\ast }(x))\mathrm{d}x}({V}_{1}-{V}_{2}-{v}_{3}{A}^{\ast }(x))\mathrm{d}x}{{\text{e}}^{\text{i}\int A(x)\mathrm{d}x}}.\end{equation} \tag{ 19 }$$

We can recover the 2 × 2 bilayer Hamiltonians h₁ and h₂ in (5) and (6), respectively, through the following identification,

$$\begin{equation}{h}_{1}\equiv h{\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)+\gamma (x)\end{subarray}}\enspace ,\qquad {h}_{2}\equiv h{\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)-\gamma (x)\\ {v}_{3}\to -{v}_{3}\end{subarray}}.\end{equation} \tag{ 20 }$$

Therefore, if we have

$$\begin{align} h\left(\begin{matrix} {\psi }_{1} \\ {\psi }_{2} \end{matrix}\right)=E\left(\begin{matrix} {\psi }_{1} \\ {\psi }_{2} \end{matrix}\right),\end{align} \tag{ 21 }$$

then there also holds

$$\begin{equation}{h}_{1}\xi =\varepsilon \enspace \xi ,\quad \xi =\left(\begin{matrix} {\xi }_{1} \\ {\xi }_{2} \end{matrix}\right)\equiv \left(\begin{matrix} {\psi }_{1} \\ {\psi }_{2} \end{matrix}\right){\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)+\gamma (x)\end{subarray}}\enspace ,\quad \varepsilon \equiv E{\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)+\gamma (x)\end{subarray}},\end{equation} \tag{ 22 }$$

$$\begin{equation}{h}_{2}\chi =\tilde{\varepsilon }\enspace \chi ,\quad \chi =\left(\begin{matrix} {\chi }_{1} \\ {\chi }_{2} \end{matrix}\right)\equiv \left(\begin{matrix} {\psi }_{1} \\ {\psi }_{2} \end{matrix}\right){\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)-\gamma (x)\\ {v}_{3}\to -{v}_{3}\end{subarray}}\enspace ,\quad \tilde{\varepsilon }\equiv E{\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)-\gamma (x)\\ {v}_{3}\to -{v}_{3}\end{subarray}}.\end{equation} \tag{ 23 }$$

The corresponding bispinor solutions of (4) are

$$\begin{align} \mathbf{\Psi }& =\frac{1}{\sqrt{2}}\left(\xi ,{\sigma }_{1}\xi \right)=\frac{1}{\sqrt{2}}{\left({\psi }_{1},\frac{-i{\psi }_{1}^{\prime }+{A}^{\ast }(x){\psi }_{1}}{E-{V}_{2}(x)},\frac{-i{\psi }_{1}^{\prime }+{A}^{\ast }(x){\psi }_{1}}{E-{V}_{2}(x)},{\psi }_{1}\right)}^{\text{T}}{\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)+\gamma (x)\end{subarray}}, \\ \mathbf{\Xi }& =\frac{1}{\sqrt{2}}\left(\chi ,-{\sigma }_{1}\chi \right)=\frac{1}{\sqrt{2}}{\left({\psi }_{1},\frac{-i{\psi }_{1}^{\prime }+{A}^{\ast }(x){\psi }_{1}}{E-{V}_{2}(x)},-\frac{-i{\psi }_{1}^{\prime }+{A}^{\ast }(x){\psi }_{1}}{E-{V}_{2}(x)},-{\psi }_{1}\right)}^{\text{T}}{\vert }_{\begin{subarray}{c}{V}_{1}(x)\equiv {{\epsilon}}_{1}(x)\\ {V}_{2}(x)\equiv {{\epsilon}}_{2}(x)-\gamma (x)\\ {v}_{3}\to -{v}_{3}\end{subarray}}.\quad \end{align} \tag{ 24 }$$

The bispinors are normalized provided that ψ and χ are normalized. Let us notice that if we can find a fundamental systems of solutions of (5) and (6) for $E=\varepsilon =\tilde{\varepsilon }$, then we can find the fundamental system for (4) for E = ɛ.

When ₂ and γ are constants, both h₁ and h₂ correspond to h with constant V₁ and V₂ (recall that the explicit form of V₂ differs in h₁ and h₂). In this case, the potential V_E and λ_E in (16) acquire this form

$$\begin{equation}{V}_{E}(x)=-\mathrm{Im}\enspace {A}^{\prime }(x)+{(\mathrm{Im}\enspace A(x))}^{2}-{v}_{3}(E-{V}_{2})\;\mathrm{Re}\enspace A(x),\end{equation} \tag{ 25 }$$

$$\begin{equation}{\lambda }_{E}=\frac{{v}_{3}^{2}}{4}{(E-{V}_{2})}^{2}+(E-{V}_{2})(E-{V}_{1})\enspace .\end{equation} \tag{ 26 }$$

For V₂ constant, the transformation (12) does not alter integrability, such that when ${\tilde{\psi }}_{1}$ is square-integrable, so is ${\psi }_{1}={\text{e}}^{-\text{i}\phi (x)}{\tilde{\psi }}_{1}$. Therefore, we get solutions of (5) and (6) as in (22) and (23), respectively. The two solutions (3) for the stationary equation (4) then read as in (24).

Now, if we allow ₁(x), ₂(x) and γ(x) to be inhomogeneous, we can still keep V₂ constant in either h₁ or h₂, provided either ₂(x) + γ(x) or ₂(x) − γ(x) is constant, respectively. For convenience, let us fix that there holds ⁴

$$\begin{equation}{V}_{2}={{\epsilon}}_{2}(x)+\gamma (x),\quad {V}_{2}\in \mathbb{R}.\end{equation} \tag{ 27 }$$

In this case, we get constant V₂ = ₂(x) + γ(x) in h₁ whereas h₂ corresponds to h with inhomogeneous V₂(x) = ₂(x) − γ(x), see (20). Notice that V₁ = ₁(x) can be inhomogeneous now. The equation (5) reduces into (16) whereas (6) leads to (14). It is feasible to solve analytically only one of the two equations by fixing ₁(x) appropriately. In either case, the analytical solutions of (4) are either Ψ or Ξ in (24) with V₁ → ₁(x). It is worth noticing that the solutions of (16) depend just on ₂(x) + γ(x), they do not 'feel' the explicit form of the inter-layer interaction γ(x) and of the on-site coupling ₂(x). These interactions can be changed without altering ${\tilde{\psi }}_{1}$ provided the change complies with (27). As we can find only part of the solutions of (4) analytically, we consider the models with inhomogeneous ₂(x) and γ(x) that satisfy (27) as quasi-exactly solvable.

Solution of (16) with the energy-dependent potential (17) is nontrivial. Identification of a complete set of solutions to an energy-dependent Schrödinger equation is, if possible, a challenging task in most cases [23, 24]. We will focus on the cases where confinement of Dirac fermions is caused either by A(x), ₁(x), or γ(x). The three situations will be discussed in the following three sections separately. Our strategy will be to identify the potential (17) with the potential of a known solvable system. This way, we will be able to identify a set of square-integrable solutions and the corresponding set of eigenvalues λ_E .

The potential term A_x σ₀ ⊗ σ₁ − A_y σ₃ ⊗ σ₂ in (1) resembles a magnetic vector potential. However, the corresponding magnetic field would have a different sign on the two layers, which might be physically unfeasible. Nevertheless, it is known that deformations of graphene layers are manifested in the form of the effective (pseudo-)magnetic vector potential [2, 25, 26]. Therefore, we can see A as a combination of the magnetic and pseudo-magnetic field that can acquire different values on the two layers, see e.g. [27].

In this section, we fix ₁, ₂ as well as γ in (2) to be constant,

$$\begin{equation}{{\epsilon}}_{1},\enspace {{\epsilon}}_{2},\enspace \gamma \in \mathbb{R}.\end{equation} \tag{ 28 }$$

It allows us to convert the task of solving both (5) and (6) into solution of the Schrödinger equation with energy-dependent potential (16). We will identify the latter equation with a stationary equation of known exactly solvable model. This way, we shall find the explicit solutions ${\tilde{\psi }}_{1}$ and energies E of (16). As a consequence, the corresponding solutions of (5) and (6) are found through the relationships (22) and (23), respectively. We will discuss two settings, the harmonic oscillator and Rosen–Morse system.

3.1. Harmonic oscillator case

Let us fix the vector potential A(x) as a complex-valued function linear in x in both its real and imaginary parts, so that

$$\begin{equation}A(x)={A}_{0}+m({w}_{r}+i\enspace {w}_{i})x\enspace ,\quad m,{w}_{r},{w}_{i}\in \mathbb{R}.\end{equation} \tag{ 29 }$$

We can interpret the vector potential term as the consequence of external homogeneous magnetic field and a mechanical strain that gives rise to homogeneous pseudo-magnetic field. Homogeneous pseudo-magnetic field in bilayer graphene was discussed in [26], whereas Dirac fermions in bilayer graphene in presence of homogeneous magnetic field were discussed in [4, 14].

The effective potential V_E (x) in (16) takes the form

$$\begin{equation}{V}_{E}(x)={m}^{2}{w}_{i}^{2}{x}^{2}-m{w}_{r}{v}_{3}(E-{V}_{2})x-m{w}_{i}-{v}_{3}{A}_{0}(E-{V}_{2}).\end{equation} \tag{ 30 }$$

Clearly, the latter implies that we must solve the eigenvalue equation of the well-known stationary oscillator. Its general solution can be found once we cast the eigenvalue equation into the confluent hypergeometric equation. It reads explicitly as (for details, see [28])

$$\begin{align} {\tilde{\psi }}_{1}& ={e}^{-\frac{m{w}_{i}}{2}{\left(x-\frac{{{\Omega}}_{E}}{{w}_{i}}\right)}^{2}}\left({\ell }_{1}\enspace {}_{1}F_{1}\left[\frac{1}{4}-\frac{{\lambda }_{E}}{4m{w}_{i}},\frac{1}{2};m{w}_{i}{\left(x-\frac{{{\Omega}}_{E}}{{w}_{i}}\right)}^{2}\right]\right. \\ & \quad +\left.{\ell }_{2}\left(x-\frac{{{\Omega}}_{E}}{{w}_{i}}\right)\enspace {}_{1}F_{1}\left[\frac{3}{4}-\frac{{\lambda }_{E}}{4m{w}_{i}},\frac{3}{2};m{w}_{i}{\left(x-\frac{{{\Omega}}_{E}}{{w}_{i}}\right)}^{2}\right]\right). \end{align} \tag{ 31 }$$

Here, ℓ₁ and ℓ₂ are constant coefficients, ${}_{1}F_{1}[a,b;z]$ stands for the Kummer or confluent hypergeometric function [29], and

$$\begin{align} {{\Omega}}_{E}& {:=}\frac{{w}_{r}{v}_{3}}{2m{w}_{i}}(E-{V}_{2})\enspace ,\quad {\lambda }_{E}=\frac{{v}_{3}^{2}}{4}{\left(1+\frac{{w}_{r}^{2}}{{w}_{i}^{2}}\right)}^{2}{(E-{V}_{2})}^{2}+(E-{V}_{2})(E-{V}_{1}+{v}_{3}{A}_{0})+m{w}_{i}. \end{align} \tag{ 32 }$$

It is worth noticing that the energy-dependent terms in the potential (30) cause just shifted decentering Ω_E of the harmonic oscillator and shift the energies λ_E .

Now, from the asymptotic behavior of the hypergeometric functions, it is straightforward to determine the physical values E for which the function ${\tilde{\psi }}_{1}$ becomes square integrable. This is achieved by imposing a polynomial behavior on the confluent hypergeometric function so that the Gaussian term (31) vanishes faster than the polynomial at x → ±∞. This leads to an exponentially vanishing function for large |x|. Such a polynomial behavior is achieved if a = −n in ${}_{1}F_{1}[a,b;z]$, with n being a non-negative integer or zero. To simplify the discussion we consider two cases. First, the conditions ℓ₂ = 0 and $\left(\frac{1}{4}-\frac{{\lambda }_{E}}{4m{w}_{i}}\right)=-n$ lead to λ_E = 2mw_i (2n + 1/2), from which we obtain ${\tilde{\psi }}_{1}$ in terms of even Hermite polynomials H_2n (z). Second, the conditions ℓ₁ = 0 and $\left(\frac{3}{4}-\frac{{\lambda }_{E}}{4m{w}_{i}}\right)=-n$ allow us to obtain λ_E = 2mw_i ((2n + 1) + 1/2), and solutions for ${\tilde{\psi }}_{1}$ in terms of odd Hermite polynomials H_2n+1(z).

We can thus unify both the even and odd solutions as

$$\begin{align} {\tilde{\psi }}_{1;n}\equiv {\tilde{\psi }}_{1}& ={e}^{-\frac{m{w}_{i}}{2}{\left(x-\frac{{{\Omega}}_{E}}{{w}_{i}}\right)}^{2}}\enspace {H}_{n}\left[\sqrt{m{w}_{i}}\left(x-\frac{{{\Omega}}_{E}}{{w}_{i}}\right)\right]\enspace ,\quad {\lambda }_{E}=2m{w}_{i}\left(n+\frac{1}{2}\right). \end{align} \tag{ 33 }$$

The physical energies E_n are determined after comparing λ_E in (32) with (33). We get

$$\begin{equation}{E}_{n}^{(\pm )}=\frac{2({V}_{1}+{V}_{2})-2{v}_{3}{A}_{0}+{V}_{2}{v}_{3}^{2}\left(1+\frac{{w}_{r}^{2}}{{w}_{i}^{2}}\right)\pm 2\sqrt{{({V}_{1}-{V}_{2}-{v}_{3}{A}_{0})}^{2}+2nm{w}_{i}\left(4+{v}_{3}^{2}\left(1+\frac{{w}_{r}^{2}}{{w}_{i}^{2}}\right)\right)}}{4+{v}_{3}^{2}\left(1+\frac{{w}_{r}^{2}}{{w}_{i}^{2}}\right)}.\end{equation} \tag{ 34 }$$

Here, it is worth remarking that the case n = 0 should be addressed with caution, as the exact value of ${E}_{0}^{(\pm )}$ depends on the sign of V₁ − V₂ − v₃ A₀. That is,

$$\begin{equation}\begin{aligned} & {E}_{0}^{(+)}=\begin{cases}{V}_{2}+4\frac{{V}_{1}-{V}_{2}-{v}_{3}{A}_{0}}{4+{v}_{3}^{2}\left(1+\frac{{w}_{r}^{2}}{{w}_{i}^{2}}\right)},\quad & {V}_{1}-{V}_{2}-{v}_{3}{A}_{0} > 0 \\ {V}_{2},\quad & {V}_{1}-{V}_{2}-{v}_{3}{A}_{0}< 0 \end{cases}, \\ & {E}_{0}^{(-)}=\begin{cases}{V}_{2},\quad & {V}_{1}-{V}_{2}-{v}_{3}{A}_{0} > 0 \\ {V}_{2}+4\frac{{V}_{1}-{V}_{2}-{v}_{3}{A}_{0}}{4+{v}_{3}^{2}\left(1+\frac{{w}_{r}^{2}}{{w}_{i}^{2}}\right)},\quad & {V}_{1}-{V}_{2}-{v}_{3}{A}_{0}< 0 \end{cases}, \end{aligned}\end{equation} \tag{ 35 }$$

from which we see that E = V₂ for either ${E}_{0}^{(+)}$ or ${E}_{0}^{(-)}$. We can find the spinor corresponding to this energy level from (19). It reads as

$$\begin{align} {\left.{\psi }_{0}^{(\pm )}\right\vert }_{E={V}_{2}}& ={\mathcal{N}}_{0}^{(\pm )}{\text{e}}^{-\text{i}\left({A}_{0}x+\frac{m{w}_{r}}{2}{x}^{2}\right)}{e}^{-\frac{m{w}_{i}}{2}{x}^{2}} \\ & \quad \times \left(\begin{matrix} 1 \\ \sqrt{\frac{\pi }{m{w}_{i}}}\frac{{V}_{1}-{V}_{2}-{v}_{3}{A}_{0}}{2}{e}^{m{w}_{i}{x}^{2}}Erf(\sqrt{m{w}_{i}}\enspace x)+\frac{{v}_{3}}{2}\left(\frac{{w}_{r}}{{w}_{i}}-i\right) \end{matrix}\right). \end{align} \tag{ 36 }$$

It is not square-integrable as its the second component diverges for x → ±∞. The latter means that E = V₂ is not a physical energy, and we thus identify the point spectrum of h as

$$\begin{equation}Sp(h)=\begin{cases}{\left\{{E}_{n}^{(+)}\right\}}_{n=0}^{\infty }\cup {\left\{{E}_{n+1}^{(-)}\right\}}_{n=0}^{\infty }\enspace ,\quad & {V}_{1}-{V}_{2}-{v}_{3}{A}_{0} > 0 \\ {\left\{{E}_{n+1}^{(+)}\right\}}_{n=0}^{\infty }\cup {\left\{{E}_{n}^{(-)}\right\}}_{n=0}^{\infty }\enspace ,\quad & {V}_{1}-{V}_{2}-{v}_{3}{A}_{0}< 0 \end{cases}.\end{equation} \tag{ 37 }$$

Now, we identify the corresponding spinor to each element in Sp(h). From the general solution (31), and after some calculations, we get

$$\begin{align} {\psi }_{n}^{(\pm )}& =\left(\begin{matrix} {\psi }_{1;n}^{(\pm )} \\ {\psi }_{2;n}^{(\pm )} \end{matrix}\right)={\mathcal{N}}_{n}^{(\pm )}{\text{e}}^{-\text{i}\left(\frac{{v}_{3}}{2}({E}_{n}^{(\pm )}-{V}_{2})+{A}_{0}\right)x}{\text{e}}^{-\text{i}\frac{m{w}_{r}}{2}{x}^{2}}{e}^{-\frac{{({z}_{n}^{(\pm )}(x))}^{2}}{2}} \\ & \quad \times \left(\begin{matrix} {H}_{n}({z}_{n}^{(\pm )}(x)) \\ -\frac{{v}_{3}}{2}\left(1+i\frac{{w}_{r}}{{w}_{i}}\right){H}_{n}({z}_{n}(x))-2in\frac{\sqrt{m{w}_{i}}}{{E}_{n}^{(\pm )}-{V}_{2}}{H}_{n-1}({z}_{n}^{(\pm )}(x)) \end{matrix}\right), \end{align} \tag{ 38 }$$

for n = 0, 1, ..., and ${E}_{0}^{(\pm )}\ne {V}_{2}$, where we have introduced the reparametrized coordinate and decentering shift

$$\begin{equation}{z}_{n}^{(\pm )}(x)=\sqrt{m{w}_{i}}\left(x-\frac{{{\Omega}}_{n}^{(\pm )}}{{w}_{i}}\right)\enspace ,\quad {{\Omega}}_{n}^{(\pm )}{:=}\frac{{w}_{r}{v}_{3}}{2m{w}_{i}}({E}_{n}^{(\pm )}-{V}_{2}),\end{equation} \tag{ 39 }$$

respectively.

Although the Hermite polynomials in (38) depend explicitly on the energy, it is still possible to compute the normalization factor for each spinor. This is done by exploiting the well-known properties of the Hermite polynomials, leading to, up to a global complex-phase,

$$\begin{equation}{\mathcal{N}}_{n}^{(\pm )}=\sqrt{\frac{1}{{2}^{n}n!}\sqrt{\frac{m{w}_{i}}{\pi }}}{\left(1+\frac{{v}_{3}^{2}}{4}{\left(1+\frac{{w}_{r}^{2}}{{w}_{i}^{2}}\right)}^{2}+\frac{2nm{w}_{i}}{{({E}_{n}^{(\pm )}-{V}_{2})}^{2}}\right)}^{-\frac{1}{2}}\enspace ,\quad \enspace n=0,1\dots ,\end{equation} \tag{ 40 }$$

which holds true only for the elements in (37).

Heretofore, we have determined the eigenvalue problem related to h, and now the corresponding information for the reduced Hamiltonians h₁ and h₂ may be extracted directly from the relationships given in (22) and (23). For clarity, we use the notation

$$\begin{equation}{\varepsilon }_{n}^{(\pm )}={\left.{E}_{n}^{(\pm )}\right\vert }_{\begin{subarray}{c}{V}_{1}\equiv {{\epsilon}}_{1}\\ {V}_{2}\equiv {{\epsilon}}_{2}+\gamma \end{subarray}}\enspace ,\qquad {\tilde{\varepsilon }}_{n}^{(\pm )}={\left.{E}_{n}^{(\pm )}\right\vert }_{\begin{subarray}{c}{V}_{1}\equiv {{\epsilon}}_{1}\\ {V}_{2}\equiv {{\epsilon}}_{2}-\gamma \\ {v}_{3}\to -{v}_{3}\end{subarray}},\end{equation} \tag{ 41 }$$

to denote the physical energies of h₁ and h₂, respectively. The corresponding bispinors Ψ_n and Ξ_n of (2) are similarly extracted via (24).

Particularly, we depict the energy levels structure of the bilayer Hamiltonian H in figure 1, for m = 1, A₀ = 0, ₁ = −1, w_r = 2, w_i = 1.8, ₂ = 1.5, v₃ = 0.3, and γ = 0.7. In this configuration, V₁ − V₂ − A₀ v₃ < 0, and so the energies ${\varepsilon }_{0}^{(+)}$ and ${\tilde{\varepsilon }}_{0}^{(+)}$ are both removed from the spectrum. Similarly, the corresponding bispinors are discarded. The behavior of the probability distributions related to the bispinors ${\mathbf{\Psi }}_{n+1}^{(+)}$, ${\mathbf{\Xi }}_{n+1}^{(+)}$, ${\mathbf{\Psi }}_{n}^{(-)}$, and ${\mathbf{\Xi }}_{n}^{(-)}$ are depicted in figures 1(a) and (b) for n = 0, 1.

Zoom In Zoom Out Reset image size

3.2. Morse–Rosen potential

Now, let us associate A(x) with a smooth step-like profile, defined in terms of a purely imaginary function of the form,

$$\begin{equation}A(x)=i{U}_{0}\left(\frac{\ell (\ell -1)}{2(\kappa -1)}+(\kappa -1)\mathrm{tanh}\enspace ({U}_{0}x)\right)\enspace ,\quad \ell ,\kappa \geqslant 1.\end{equation} \tag{ 42 }$$

In this case, the effective potential in (16) becomes

$$\begin{align} {V}_{E}(x)& =-{U}_{0}^{2}\kappa (\kappa -1){[sech({U}_{0}x)]}^{2}+{U}_{0}^{2}\ell (\ell -1)\mathrm{tanh}({U}_{0}x) \\ & \quad +{U}_{0}^{2}\left[\frac{{\ell }^{2}{(\ell -1)}^{2}}{4{(1-\kappa )}^{2}}+{(1-\kappa )}^{2}\right], \end{align} \tag{ 43 }$$

which corresponds to the Morse–Rosen interaction [30] (also known as the hyperbolic Rosen–Morse potential), one of the well-known exactly solvable models in quantum mechanics [28, 31–33]. Interestingly, for ℓ = 1, the potential (43) reduces to the Pöschl–Teller interaction, a particular case to be discussed in detail in the upcoming sections. Let us notice that Dirac fermions in (single-layer) graphene in presence of (42) were discussed in [34].

For the rest of this section, we focus on the case ℓ ⩾ 1 and κ ≠ 1. We will consider bound states that comply with the boundary condition ${\left.{\tilde{\psi }}_{1}\right\vert }_{x\to \pm \infty }\to 0$. Thus, taking the differential equation for ${\tilde{\psi }}_{1}$ into the hypergeometric form [28], and after some calculations, we get the eigenstates

$$\begin{equation}{\tilde{\psi }}_{1;n}={\mathcal{C}}_{n}{\left(1-\mathrm{tanh}({U}_{0}x)\right)}^{\frac{{\alpha }_{n}}{2}}{\left(1+\mathrm{tanh}({U}_{0}x)\right)}^{\frac{{\beta }_{n}}{2}}{P}_{n}^{({\alpha }_{n},{\beta }_{n})}(\mathrm{tanh}({U}_{0}x)),\end{equation} \tag{ 44 }$$

$$\begin{align} {\lambda }_{E}& =-{U}_{0}^{2}({\beta }_{n}^{2}-{\beta }_{0}^{2}) \\ & =-{U}_{0}^{2}\left(\frac{{\ell }^{2}{(\ell -1)}^{2}}{4{(n+1-\kappa )}^{2}}+{(n+1-\kappa )}^{2}-\frac{{\ell }^{2}{(\ell -1)}^{2}}{4{(1-\kappa )}^{2}}-{(1-\kappa )}^{2}\right), \end{align} \tag{ 45 }$$

where n = 0, 1, ..., n_max, and ${P}_{n}^{(\alpha ,\beta )}(z)$ stands for the Jacobi polynomials, with

$$\begin{equation}{\alpha }_{n}{:=}(\kappa -n-1)-\frac{\ell (\ell -1)}{2(n+1-\kappa )}\enspace ,\quad {\beta }_{n}{:=}(\kappa -n-1)+\frac{\ell (\ell -1)}{2(n+1-\kappa )},\end{equation} \tag{ 46 }$$

and the normalization constant [31]

$$\begin{equation}{\mathcal{C}}_{n}={2}^{n+1-\kappa }{\left(\frac{{U}_{0}n!{\alpha }_{n}{\beta }_{n}{\Gamma}(2\kappa -n-1)}{(\kappa -n-1){\Gamma}({\alpha }_{n}+n+1){\Gamma}({\beta }_{n}+n+1)}\right)}^{\frac{1}{2}}.\end{equation} \tag{ 47 }$$

From (44), it follows that ${\tilde{\psi }}_{1;n}$ is square-integrable only when both α_n , β_n > 0. This leads us to a condition for the existence of at least one bound state and an upper bound n_max given by

$$\begin{equation}{(\kappa -1)}^{2} > \frac{\ell (\ell -1)}{2}\enspace ,\quad {n}_{\mathrm{max}}=\lfloor \kappa -1-\sqrt{\frac{\ell (\ell -1)}{2}}\rfloor ,\end{equation} \tag{ 48 }$$

respectively.

The physical energies E_n are then determined by comparing λ_E in (45) with (26), from which we obtain a polynomial equation of second-order for E_n . We thus get the energies

$$\begin{equation}{E}_{n}^{(\pm )}=\frac{2({V}_{1}+{V}_{2})+{v}_{3}^{2}{V}_{2}\pm 2\sqrt{{({V}_{1}-{V}_{2})}^{2}+{U}_{0}^{2}(4+{v}_{3}^{2})({\beta }_{0}^{2}-{\beta }_{n}^{2})}}{4+{v}_{3}^{2}}.\end{equation} \tag{ 49 }$$

The second spinor component ψ_2;n is determined from (11), and the corresponding spinors take the form

$$\begin{align} {\psi }_{n}^{(\pm )}& =\left(\begin{matrix} {\psi }_{1;n}^{(\pm )} \\ {\psi }_{2;n}^{(\pm )} \end{matrix}\right)={\text{e}}^{-\text{i}\frac{{v}_{3}}{2}({E}_{n}^{(\pm )}-{V}_{2})x}{(1-z)}^{\frac{{\alpha }_{n}}{2}}{(1+z)}^{\frac{{\beta }_{n}}{2}} \\ & \quad \times \left(\begin{matrix} {P}_{n}^{({\alpha }_{n},{\beta }_{n})}(z) \\ -\left(\frac{{v}_{3}}{2}+\frac{i{U}_{0}}{({E}_{n}^{(\pm )}-{V}_{2})}\left(\frac{\ell (\ell -1)}{2(\kappa -1)}+\frac{\ell (\ell -1)}{2(n+1-\kappa )}+nz\right)\right){P}_{n}^{({\alpha }_{n},{\beta }_{n})}(z) \\ \qquad -i\frac{{U}_{0}(2\kappa -1-n)(1-{z}^{2})}{2({E}_{n}^{(\pm )}-{V}_{2})}{P}_{n-1}^{({\alpha }_{n}+1,{\beta }_{n}+1)}(z), \end{matrix}\right) \end{align} \tag{ 50 }$$

with z ≡ z(x) := tanh(U₀ x).

The exact value of ${E}_{0}^{(\pm )}$ depends on the sign of V₁ − V₂. In analogy to the oscillator-like interaction of the previous section, we have ${E}_{0}^{(+)}={V}_{2}$ for V₁ < V₂, and ${E}_{0}^{(-)}={V}_{2}$ for V₁ > V₂. The corresponding spinors (19) are not square integrable so that these values do not belong to the point spectrum of h.

We found the following set Sp(h) of discrete energies of h,

$$\begin{equation}Sp(h)=\begin{cases}{\left\{{E}_{n}^{(+)}\right\}}_{n=1}^{{n}_{\mathrm{max}}}\cup {\left\{{E}_{n}^{\left(\right.-}\right\}}_{n=0}^{{n}_{\mathrm{max}}}\enspace ,\quad {V}_{1}< {V}_{2} \\ {\left\{{E}_{n}^{(+)}\right\}}_{n=0}^{{n}_{\mathrm{max}}}\cup {\left\{{E}_{n}^{\left(\right.-}\right\}}_{n=1}^{{n}_{\mathrm{max}}}\enspace ,\quad {V}_{1} > {V}_{2} \end{cases}.\end{equation} \tag{ 51 }$$

Since α_n and β_n do not depend on the energy upper-index (±), the upper bound n_max is the same for both energies ${E}_{n}^{(\pm )}$. For n_max = 0, the set of discrete energies is just $\left\{{E}_{0}^{(-)}\right\}$ for V₁ < V₂, and $\left\{{E}_{0}^{(+)}\right\}$ for V₁ > V₂. Moreover, the point spectrum may be empty if the inequality in (48) is not fulfilled.

Now, the discrete energy levels associated with h₁ and h₂ are obtained via (22) and (23),

$$\begin{equation*}Sp({h}_{1})=Sp(h){\vert }_{\begin{subarray}{c}{V}_{1}\equiv {{\epsilon}}_{1}\\ {V}_{2}\equiv {{\epsilon}}_{2}+\gamma \end{subarray}},\qquad Sp({h}_{2})=Sp(h){\vert }_{\begin{subarray}{c}{V}_{1}\equiv {{\epsilon}}_{1}\\ {V}_{2}\equiv {{\epsilon}}_{2}-\gamma \\ {v}_{3}\to -v3\end{subarray}}.\end{equation*}$$

The bispinor solutions of (4) can be obtained via (24). In figure 2, we illustrate calculated energy levels of h₁, h₂ and H, and probability density related to the bispinors Ψ_n and Ξ_n . In particular, we have considered U₀ = 1, κ = 2.5, ℓ = 1.1, ₁ = −1, ₂ = 1.5, v₃ = 0.3, and γ = 0.7. In such a case, we obtain n_max = 1, so that we generate two physical solutions for each reduced Hamiltonian. On the other hand, for both h₁ and h₂, we get V₁ < V₂, which means that ${E}_{0}^{(+)}$ is discarded from the point spectrum. Each reduced Hamiltonian contributes with three physical energies, and so the bilayer Hamiltonian H contains six energy levels. See figure 2(c).

Zoom In Zoom Out Reset image size

In this section, we focus on the case where A(x) in (2) vanishes. The on-site interactions ₁(x), ₂(x) as well as the inter-layer coupling γ(x) can be inhomogeneous. Let us suppose that ₂(x) and γ(x) are related by (27), i.e. V₂ = ₂(x) + γ(x) is constant. It brings the equations (5) and (6) into

$$\begin{equation}{h}_{1}\xi =\left(-i\left(\begin{matrix} {v}_{3} & 1 \\ 1 & 0 \end{matrix}\right){\partial }_{x}+\left(\begin{matrix} {{\epsilon}}_{1}(x) & 0 \\ 0 & {V}_{2} \end{matrix}\right)\right)\xi =E\enspace \xi ,\end{equation} \tag{ 52 }$$

$$\begin{equation}{h}_{2}\chi =\left(-i\left(\begin{matrix} -{v}_{3} & 1 \\ 1 & 0 \end{matrix}\right){\partial }_{x}+\left(\begin{matrix} {{\epsilon}}_{1}(x) & 0 \\ 0 & {{\epsilon}}_{2}(x)-\gamma (x) \end{matrix}\right)\enspace \right)\chi =E\chi .\end{equation} \tag{ 53 }$$

As discussed in the section 2, the solutions of (52) can be found via the Schrödinger equation with energy-dependent potential (16). It acquires the following simple form

$$\begin{equation}-{\tilde{\psi }}_{1}^{{\prime\prime}}-(E-{V}_{2})(E-{{\epsilon}}_{1}(x)){\tilde{\psi }}_{1}=\frac{{v}_{3}^{2}}{4}{(E-{V}_{2})}^{2}{\tilde{\psi }}_{1}{\vert }_{{V}_{2}={{\epsilon}}_{2}+\gamma },\end{equation} \tag{ 54 }$$

where the spinor components of $\xi ={({\xi }_{1},{\xi }_{2})}^{\text{T}}$ are determined from ${\tilde{\psi }}_{1}$ through

$$\begin{equation}{\xi }_{1}={\text{e}}^{-\text{i}\frac{{v}_{3}}{2}(E-{V}_{2})x}{\tilde{\psi }}_{1}\enspace ,\qquad {\xi }_{2}=-\frac{i\enspace {\psi }_{1}^{\prime }}{({V}_{2}-E)},\end{equation} \tag{ 55 }$$

see (22) and (23). The equation (53) reduces into (14). With current fixing of the quantities and denoting $\mathcal{B}(x)={{\epsilon}}_{2}(x)-\gamma (x)$, the latter equation reads as

$$\begin{align} -{\tilde{\chi }}_{1}^{{\prime\prime}}& +\left(-\frac{{v}_{3}^{2}}{4}{(E-\mathcal{B}(x))}^{2}-(E-{{\epsilon}}_{1}(x))(E-\mathcal{B}(x))\right. \\ & \quad \left.+\frac{3{(\mathcal{B}{(x)}^{\prime }(x))}^{2}}{4{(E-\mathcal{B}(x))}^{2}}+\frac{\mathcal{B}{(x)}^{{\prime\prime}}(x)}{2(E-\mathcal{B}(x))}\right){\tilde{\chi }}_{1}=0, \end{align} \tag{ 56 }$$

where ${\chi }_{1}=\sqrt{E-\mathcal{B}(x)}{\text{e}}^{-\text{i}\frac{{v}_{3}}{2}\left(x-{\int }^{x}\mathcal{B}(s)\mathrm{d}s\right)}{\tilde{\chi }}_{1}$, see (12) and (23).

We find it physically reasonable to consider the systems where ₁(x) is bounded. We shall match either (54) or (56) with the stationary equation of a solvable quantum systems that meets these requirements. Let us consider the stationary equation of the Pöschl–Teller system

$$\begin{equation}-{\tilde{\psi }}_{1}^{{\prime\prime}}-\kappa (\kappa -1){U}_{0}^{2}\enspace {\text{sech}}^{2}\enspace {U}_{0}x\enspace {\tilde{\psi }}_{1}-\lambda {\tilde{\psi }}_{1}=0,\quad \kappa > 1.\end{equation} \tag{ 57 }$$

It is worth noticing that Dirac electrons in bilayer graphene were studied in presence of Pöschl–Teller electrostatic potential in [35], see also [36]. The equation (57) is a special case (ℓ = 1) of the Rosen–Morse equation (43) that was discussed in the previous section. Therefore, we can use (44) and write down the square integrable solutions ${\tilde{\psi }}_{1;n}$ of (57),

$$\begin{equation}{\tilde{\psi }}_{1;n}={\mathcal{C}}_{n}^{(\pm )}sech{({U}_{0}x)}^{{\alpha }_{n}}{P}_{n}^{({\alpha }_{n},{\alpha }_{n})}(\mathrm{tanh}({U}_{0}x)),\end{equation} \tag{ 58 }$$

where n = 0, 1, ..., n_max = ⌊κ − 1⌋, and α_n = (κ − n − 1). The corresponding eigenvalues λ_n are

$$\begin{equation}{\lambda }_{n}=-{U}_{0}^{2}{(1-\kappa +n)}^{2}=-{U}_{0}^{2}\enspace {\alpha }_{n}^{2}.\end{equation} \tag{ 59 }$$

Now, we shall identify either (54) or (56) with (57). Let us start with (54),

$$\begin{equation}(E-{V}_{2})(E-{{\epsilon}}_{1}(x))+\frac{{v}_{3}^{2}}{4}{(E-{V}_{2})}^{2}=\kappa (\kappa -1)\enspace {U}_{0}^{2}\enspace {\text{sech}}^{2}\enspace {U}_{0}x-{U}_{0}^{2}\enspace {(1-\kappa +n)}^{2}.\end{equation} \tag{ 60 }$$

There are different ways how to fix ₁(x), and each of them leads to different values of E. Let us discuss some of them.

4.1. Case I

To begin with, let us consider the inhomogeneous on-site interaction

$$\begin{equation}{{\epsilon}}_{1}(x)=\mathcal{A}\enspace {U}_{0}^{2}{sech}^{2}\enspace {U}_{0}x,\quad \mathcal{A} > 0.\end{equation} \tag{ 61 }$$

Then (60) is satisfied provided that we fix κ ≡ κ and E such that they solve the following two equations,

$$\begin{equation}{\kappa }_{{\epsilon}}({\kappa }_{{\epsilon}}-1)=\mathcal{A}\enspace ({V}_{2}-E)\enspace ,\quad {U}_{0}^{2}\enspace {(1-{\kappa }_{{\epsilon}}+n)}^{2}=({V}_{2}-E)\left(\frac{{v}_{3}^{2}}{4}({V}_{2}-E)-E\right).\end{equation} \tag{ 62 }$$

That is, the Pöschl–Teller amplitude κ depends explicitly on the energy through

$$\begin{equation}{\kappa }_{{\epsilon}}=\frac{1}{2}+\sqrt{\mathcal{A}({V}_{2}-E)+\frac{1}{4}}.\end{equation} \tag{ 63 }$$

On the other hand, the second equation in (62) yields to a fourth-order polynomial equation for E of the form

$$\begin{equation}({V}_{2}-E)\left(\left(1+\frac{{v}_{3}^{2}}{4}\right)({V}_{2}-E)-{V}_{2}\right)=-{U}_{0}^{2}{\left(n+\frac{1}{2}-\sqrt{\mathcal{A}\enspace ({V}_{2}-E)+\frac{1}{4}}\right)}^{2},\end{equation} \tag{ 64 }$$

the solutions of which become unfeasible to obtain in the general setup. Despite such complexity, we can proceed further and obtain some additional information.

The straightforward calculations show that the square-integrable condition (59) still holds in this case, with κ → κ given in (63). Thus, with the current choice of parameters, the solutions of (57)

$$\begin{equation}{\tilde{\psi }}_{1;n}={\mathcal{C}}_{n}sech{({U}_{0}x)}^{{\bar{\alpha }}_{n}}{P}_{n}^{({\bar{\alpha }}_{n},{\bar{\alpha }}_{n})}(\mathrm{tanh}({U}_{0}x)),\end{equation} \tag{ 65 }$$

$$\begin{equation}{\bar{\alpha }}_{n}=-n-1+{\kappa }_{{\epsilon}}=-n-\frac{1}{2}+\sqrt{\mathcal{A}({V}_{2}-E)+\frac{1}{4}},\end{equation} \tag{ 66 }$$

are square-integrable provided that ${\bar{\alpha }}_{n}$ is positive. In order to keep ${\bar{\alpha }}_{n}$ real, the term inside the square-root of ${\bar{\alpha }}_{n}$ must be positive. Since $\mathcal{A} > 0$, we get $E< {V}_{2}+\frac{1}{4\mathcal{A}}$. Additionally, the requirement ${\bar{\alpha }}_{n} > 0$ is satisfied provided that $E< {V}_{2}-\frac{n(n+1)}{\mathcal{A}}$. Now, the right-hand side of (64) is negative, as it a multiple of $-{\bar{\alpha }}_{n}^{2}$. Therefore, one obtains real solutions of (64) for (V₂ − E) only if the term on the left is negative as well, which is quadratic and convex on (V₂ − E). We obtain $0< ({V}_{2}-E)< \frac{4{V}_{2}}{4+{v}_{3}^{2}}$, for V₂ > 0, and $\frac{4{V}_{2}}{4+{v}_{3}^{2}}< ({V}_{2}-E)< 0$, for V₂ < 0. Therefore, any real solution of (64) has to lie inside the one of the following intervals

$$\begin{equation}\begin{aligned} & E\in \left(\frac{{v}_{3}^{2}}{4+{v}_{3}^{2}}{V}_{2},{V}_{2}\right)\cap \left(-\infty ,{V}_{2}-\frac{n(n+1)}{\mathcal{A}}\right)\enspace , & & for\enspace {V}_{2} > 0, \\ & E\in \left({V}_{2},\frac{{v}_{3}^{2}}{4+{v}_{3}^{2}}{V}_{2}\right)\cap \left(-\infty ,{V}_{2}-\frac{n(n+1)}{\mathcal{A}}\right)\enspace ,\quad & & for\enspace {V}_{2}< 0, \end{aligned}\enspace \end{equation} \tag{ 67 }$$

where in the latter is clear that $\frac{{v}_{3}^{2}}{4+{v}_{3}^{2}}< 1$, for ${v}_{3}\in \mathbb{R}$. The requirement that the intersections are non-empty sets an upper bound for possible values of n. Indeed, when V₂ > 0, the intersection is nonempty for n ⩽ n_max where

$$\begin{equation}{n}_{\mathrm{max}}=\lfloor \sqrt{\frac{4\mathcal{A}{V}_{2}}{4+{v}_{3}^{2}}+\frac{1}{4}}-\frac{1}{2}\rfloor \enspace ,\quad for\enspace {V}_{2} > 0.\end{equation} \tag{ 68 }$$

It provides us with an upper bound for the maximum number of physical solutions which is n_max + 1. It is worth noticing that the trigonal warping term acts against the confinement here; the larger is |v₃|, the smaller is n_max. When V₂ < 0, it is clear that only n = 0 leads to an non-empty intersection of the energy intervals. This yields to E = V₂. Nevertheless, the expression (19) suggests that the corresponding spinor is not square integrable. Thus, such a solution is discarded, and no physical solutions are produced for V₂ < 0.

Interestingly, even if E has to be found by numerical means, we have been able to extract general information about the spectrum and number of physically allowed solutions. Furthermore, the spinor may be computed explicitly from (55) and (65) as

$$\begin{align} {\psi }_{n}& =\left(\begin{matrix} {\psi }_{1;n} \\ {\psi }_{2;n} \end{matrix}\right)={\mathcal{N}}_{n}\enspace {\text{e}}^{\text{i}\frac{{v}_{3}}{2}({V}_{2}-{E}_{n})x}{\left(sech({U}_{0}x)\right)}^{{\bar{\alpha }}_{n}}\times \\ & \quad \left(\begin{matrix} {P}_{n}^{({\bar{\alpha }}_{n},{\bar{\alpha }}_{n})}(z(x)) \\ -\left(\frac{{v}_{3}}{2}+i\frac{{U}_{0}{\bar{\alpha }}_{n}}{{V}_{2}-{E}_{n}}\enspace \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}({U}_{0}x)\right){P}_{n}^{({\bar{\alpha }}_{n},{\bar{\alpha }}_{n})}(z(x))+i\frac{{U}_{0}(n+1+2{\bar{\alpha }}_{n})}{2({V}_{2}-{E}_{n})}\enspace {sech}^{2}({U}_{0}x){P}_{n-1}^{({\bar{\alpha }}_{n}+1,{\bar{\alpha }}_{n}+1)}(z(x)) \end{matrix}\right), \end{align} \tag{ 69 }$$

where z(x) := tanh(U₀ x), and n = 0, 1, ..., n_max. The corresponding set of bispinors follow from (3), which in this case are given through

$$\begin{equation}{\mathbf{\Psi }}_{n}=\frac{1}{\sqrt{2}}{({\psi }_{1;n},{\psi }_{2;n},{\psi }_{2;n},{\psi }_{1;n})}^{\text{T}}\enspace ,\quad n=0,1,\dots ,{n}_{\mathrm{max}}\enspace .\end{equation} \tag{ 70 }$$

To illustrate our results, we consider V₂ = 2.3, v₃ = 0.1, U₀ = 1, and $\mathcal{A}=2.7$, so that real energies lie inside the interval E ∈ (0.005 735, 2.3). Moreover, from (68), the maximum number of physical energies is n_max + 1 = 3. The numerical values for E_n and ${\bar{\alpha }}_{n}$ are shown in table 1. For each n, we obtain two energies E_n , which are all real for n = 0, 1, 2, and complex for n ⩾ 3. The physical energies are identified as E₀ = 1.425 92, E₁ = 0.368 25, and E₂ = 0.006 3725 as they render ${\bar{\alpha }}_{n}$ positive. The associated probability distributions are depicted in figure 3(a).

Table 1. Numerical solutions of the characteristic equation (64). We have fixed the parameters as $\mathcal{A}=2.7$, v₃ = 0.1, U₀ = 1. For h₁ we used V₂ = ₂ + γ = 2.3, whereas for h₂ we have V₂ = ₂ − γ = 1.7.

	h1		h2
n	En	${\bar{\alpha }}_{n}$	En	${\bar{\alpha }}_{n}$
0	2.3	0	1.7	0
	1.425 92	1.115 55	1.185 07	0.780 747
1	2.112 15	−0.629 828	1.475 78	−0.575 127
	0.368 25	0.837 884	0.217 073	0.562 499
2	1.615 15	−1.051 17	0.779 102	−0.845 786
	0.006 3725	0.038 266	0.079 2946	−0.349 208
3	0.462 606 + 0.648 775i	−1.186 08 − 0.378 511i	0.089 0298 + 1.035 22i	−1.265 95 − 0.625 57i
	0.462 606 − 0.648 775i	−1.186 08 + 0.378 511i	0.089 0298 − 1.035 22i	−1.265 95 + 0.625 57i

Zoom In Zoom Out Reset image size

4.1.1. Particular setup for h₂

As we have remarked, the solutions of (54) are insensitive with respect to the explicit form of γ and ₂, and thus we can impose some further restrictions in order to obtain information about the reduced Hamiltonian h₂. Particularly, if

$$\begin{equation}{{\epsilon}}_{2},\enspace \gamma \in \mathbb{R},\end{equation} \tag{ 71 }$$

the equation (56) brings us back to the system (62), where now ${V}_{2}\equiv \mathcal{B}={{\epsilon}}_{2}-\gamma$. Therefore, we can solve the equation in exactly the same manner as we did in table 1 for different values of V₂ now. This is illustrated in the fourth and fifth columns of table 1, where the energy levels are determined for V₂ = ₂ − γ = 1.6, v₃ = 0.1, U₀ = 1, and $\mathcal{A}=2.7$. The corresponding probability distributions associated with the bispinors Ξ_n are depicted in figure 3(b).

4.2. Case II

Fixing of ₁(x) in (61) allowed us to have the on-site interaction independent on n. The price we paid was that the equation (64) had to be solved numerically. Let us consider the other choice of parameters such that (60) is satisfied. We fix

$$\begin{equation}{{\epsilon}}_{1}(x)=\kappa \enspace \mathcal{A}\enspace {U}_{0}^{2}{sech}^{2}\enspace {U}_{0}x.\end{equation} \tag{ 72 }$$

Additionally, there must hold

$$\begin{equation}-\mathcal{A}(E-{V}_{2})=(\kappa -1){U}_{0}^{2},\quad (E-{V}_{2})\left(\frac{{v}_{3}^{2}}{4}(E-{V}_{2})+E\right)=-{U}_{0}{(1-\kappa +n)}^{2}.\end{equation} \tag{ 73 }$$

We can see that inclusion of κ into ₁(x) lowered the order of κ in the second equation in (73) when compared to (62). We can solve (73) for E and either κ, U₀ or V₂. As we require ₁(x) to be independent of n, we solve the equation for V₂, and E,

$$\begin{align} {V}_{2}(\mathcal{A},\kappa ,n)& =\left(1+\frac{{v}_{3}^{2}}{4}\right)\frac{(\kappa -1)}{\mathcal{A}}+\frac{\mathcal{A}\enspace {(1+n-\kappa )}^{2}{U}_{0}^{2}}{\kappa -1}, \\ E(\mathcal{A},\kappa ,n)& =\frac{\mathcal{A}\enspace {U}_{0}^{2}\enspace {(1+n-\kappa )}^{2}}{\kappa -1}+\frac{(\kappa -1){v}_{3}^{2}}{4\mathcal{A}}. \end{align} \tag{ 74 }$$

${V}_{2}(\mathcal{A},\kappa ,n)$ is a parabola in n with minimum ${V}_{2}({n}_{0})=\frac{(\kappa -1)}{\mathcal{A}}\left(1+\frac{{v}_{3}^{2}}{4}\right)$ at n₀ = κ − 1. For each n from the allowed interval n ∈ {0, ..., ⌊κ − 1⌋}, the on-site interaction ₁(x) remains the same. Nevertheless, the value of V₂ gets changed correspondingly.

For each of this specific configurations, we are able to find a localized solution ${\tilde{\psi }}_{1;n}$. The bispinor solution of (4) is then

$$\begin{align} \mathbf{\Psi }& =\left({\psi }_{1},\frac{-i{\psi }_{1}^{\prime }}{E-{V}_{2}(\mathcal{A},\kappa ,n)},\frac{-i{\psi }_{1}^{\prime }}{E-{V}_{2}(\mathcal{A},\kappa ,n)},{\psi }_{1}\right),\quad {\psi }_{1}={\text{e}}^{-\text{i}\frac{{v}_{3}}{2}(E-{V}_{2})x}{\tilde{\psi }}_{1;n}{\vert }_{\begin{subarray}{c}E=E(\mathcal{A},\kappa ,n)\\ {V}_{2}={V}_{2}(\mathcal{A},\kappa ,n)\end{subarray}}. \end{align} \tag{ 75 }$$

The solution is invariant with respect to the changes of ₂(x) and γ(x) that preserve V₂, including the case where both ₂(x) and γ(x) are constant. When this is the case, the equation (56) reduces into the equation that coincides with (54), yet for V₂ = ₂ − γ. If ₂ and γ are such that

$$\begin{equation}{{\epsilon}}_{2}+\gamma ={V}_{2}(\mathcal{A},\kappa ,n)\quad \text{and}\quad {{\epsilon}}_{2}-\gamma ={V}_{2}(\mathcal{A},\kappa ,\tilde{n}),\quad n,\tilde{n}\in \left\{0,\dots ,\lfloor \kappa -1\rfloor \right\},\end{equation} \tag{ 76 }$$

then we can get bound state solutions for each equations (52) and (53) with energies ${\epsilon}=E(\mathcal{A},\kappa ,n)$ and $\tilde{\varepsilon }=E(\mathcal{A},\kappa ,\tilde{n})$. The bispinor Ξ corresponding to the latter energy is

$$\begin{align} \mathbf{\Xi }& =\left({\psi }_{1},\frac{-i{\psi }_{1}^{\prime }}{E-{V}_{2}(\mathcal{A},\kappa ,\tilde{n})},\frac{i{\psi }_{1}^{\prime }}{E-{V}_{2}(\mathcal{A},\kappa ,\tilde{n})},-{\psi }_{1}\right),\quad {\psi }_{1}={\text{e}}^{\text{i}\frac{{v}_{3}}{2}(E-{V}_{2})x}{\tilde{\psi }}_{1;\tilde{n}}{\vert }_{\begin{subarray}{c}E=E(\mathcal{A},\kappa ,\tilde{n})\\ {V}_{2}={V}_{2}(\mathcal{A},\kappa ,\tilde{n})\end{subarray}}. \end{align} \tag{ 77 }$$

We show density of states of the corresponding bispinors Ψ and Ξ in figure 4.

Zoom In Zoom Out Reset image size

4.3. Case III

Let us relax the condition (27), i.e. both ₂ + γ and ₂ − γ can be inhomogeneous. In this case, the equations (5) and (6) lead to (14), yet with different form of V₂(x) in each case. As we mentioned in section 2, we can find configuration of ₁(x), ₂(x) and γ(x) such that (14) is partially solvable. Let us consider the case with V₂(x) = γ(x) + ₂(x). We fix V₁(x) ≡ ₁(x) as

$$\begin{equation}{{\epsilon}}_{1}(x)=E+\frac{{v}_{3}^{2}}{4}(E-{V}_{2}(x))-\frac{3{({V}_{2}^{\prime }(x))}^{2}}{4{(E-{V}_{2}(x))}^{3}}-\frac{{V}_{2}^{{\prime\prime}}(x)}{2{(E-{V}_{2}(x))}^{2}}+\frac{{V}_{0}(x)}{E-{V}_{2}(x)}.\end{equation} \tag{ 78 }$$

Then the equation (14) reduces into

$$\begin{equation}-{\tilde{\psi }}_{1}^{{\prime\prime}}+{V}_{0}(x){\tilde{\psi }}_{1}=0.\end{equation} \tag{ 79 }$$

Let us identify (79) with the stationary equation of the Pöschl–Teller system again, ${V}_{0}(x)=-\kappa (\kappa -1){U}_{0}^{2}\enspace {sech}^{2}\enspace {U}_{0}x+{U}_{0}^{2}{(1+n-\kappa )}^{2}$. When n is a positive integer n ∈ {0, ⌊κ − 1⌋}, the equation (79) has a localized solution ${\tilde{\psi }}_{1;n}$, see (58). Notice that this solution is independent on the explicit choice of V₂(x) and E as they do not appear in (79). Nevertheless, both V₂(x) = γ(x) + ₂(x) and E affect the form of the bispinor solution (24) via (12),

$$\begin{align} & \mathbf{\Psi }(x)=\left({\psi }_{1}(x),\frac{-i{\psi }_{1}^{\prime }(x)}{E-{{\epsilon}}_{2}(x)-\gamma (x)},\frac{-i{\psi }_{1}^{\prime }(x)}{E-{{\epsilon}}_{2}(x)-\gamma (x)},{\psi }_{1}\right), \\ & {\psi }_{1}(x)=\sqrt{E-{{\epsilon}}_{2}(x)-\gamma (x)}{\text{e}}^{-\text{i}\frac{{v}_{3}}{2}(E-{{\epsilon}}_{2}(x)-\gamma (x))x}{\tilde{\psi }}_{1;\tilde{n}}(x). \end{align} \tag{ 80 }$$

In the current setting, E plays rather the role of an interaction parameter. We can tune the interaction ₁(x) by changing E such that it confines a bound state with energy equal to E. In order to keep Ψ(x) square-integrable, we require that $\sqrt{E-{V}_{2}(x)}$ is a bounded function.

We fix V₂(x) = ₂(x) + γ, $\gamma \in \mathbb{R}$. When ₂(x) is periodic, ₁(x) shares its periodicity up to the last term in (78) that represents a periodicity defect. For instance, if we fix

$$\begin{equation}{{\epsilon}}_{2}(x)={V}_{2}(x)-\gamma =c\enspace \mathrm{cos}\enspace x,\quad 0< c< E,\end{equation} \tag{ 81 }$$

then ₁(x) reads explicitly

$$\begin{equation}{{\epsilon}}_{1}(x)=E+\frac{(E-c\enspace \mathrm{cos}\enspace x){v}_{3}^{2}}{4}+\frac{-2c\enspace \mathrm{cos}\enspace x(c\mathrm{cos}\enspace x-E)-3{c}^{2}\enspace {\mathrm{sin}}^{2}\enspace x}{4{(E-c\enspace \mathrm{cos}\enspace x)}^{3}},\end{equation} \tag{ 82 }$$

$$\begin{equation}\qquad +\frac{{(1+n-\kappa )}^{2}{U}_{0}^{2}}{E-c\enspace \mathrm{cos}\enspace x}+\frac{(1-\kappa )\kappa {U}_{0}^{2}\enspace {sech}^{2}\enspace {U}_{0}x}{E-c\enspace \mathrm{cos}\enspace x}.\end{equation} \tag{ 83 }$$

We illustrate the interactions for different values of parameters in figure 5 together with density of probability of the bound state.

Zoom In Zoom Out Reset image size

Up to now, the inter-layer coupling γ had rather implicit influence on considered solutions as it was 'hidden' in V₂. Let us see whether we can get an analytical solution of confined Dirac fermions by inhomogeneous γ. It is worth noticing in this context that confinement by inhomogeneous γ with rotational symmetry was analyzed numerically in [7, 8]. We focus on the situation where ₁ and ₂ are constant (in order to eliminate localization by on-site potentials) and γ is inhomogeneous. We fix

$$\begin{equation}{v}_{3}=0,\quad {V}_{1}={{\epsilon}}_{1}=\text{const.},\quad {{\epsilon}}_{2}=\text{const.},\quad {V}_{2}(x)={{\epsilon}}_{2}+\gamma (x).\end{equation} \tag{ 84 }$$

Then the equations (5) and (6) for the spinor components ξ_j and χ_j , for j = 1, 2, can be decoupled through the relationships

$$\begin{equation}-\frac{i{\xi }_{2}^{\prime }}{E-{{\epsilon}}_{1}}={\xi }_{1},\quad -\frac{i{\chi }_{2}^{\prime }}{E-{{\epsilon}}_{1}}={\chi }_{1},\end{equation} \tag{ 85 }$$

leading to the effective energy-dependent Schrödinger equation

$$\begin{align} -{\xi }_{2}^{{\prime\prime}}-({{\epsilon}}_{1}-E)\gamma (x){\xi }_{2}& =({{\epsilon}}_{1}-E)({{\epsilon}}_{2}-E){\xi }_{2}\enspace ,\quad -{\chi }_{2}^{{\prime\prime}}+({{\epsilon}}_{1}-E)\gamma (x){\chi }_{2} \\ & =({{\epsilon}}_{1}-E)({{\epsilon}}_{2}-E){\chi }_{2}. \end{align} \tag{ 86 }$$

Notice that the difference between the latter equations relies on the sign of the energy-dependent potential, which both coincide qualitatively with (54), therefore, we can follow the same steps as in the previous section. We shall identify both equations in (86) with (57). We set ⁵

$$\begin{equation}\gamma (x)=\mathcal{A}\enspace {U}_{0}^{2}{sech}^{2}\enspace {U}_{0}x+{\gamma }_{0}\enspace ,\quad {\gamma }_{0},\enspace \mathcal{A} > 0\enspace ,\quad {U}_{0}\in \mathbb{R}.\end{equation} \tag{ 87 }$$

It has similar form to ₁(x) in (61), however, it acquires nonvanishing constant value γ₀ asymptotically now. In analogy to the previous section, we identify the set of solutions as

$$\begin{equation}\begin{aligned} & {\psi }_{2;n}={\mathcal{C}}_{2}{(sech({U}_{0}x))}^{{\eta }_{n}}{P}_{n}^{({\eta }_{n},{\eta }_{n})}(\mathrm{tanh}({U}_{0}x)), \\ & {\chi }_{2;n}={\mathcal{D}}_{2}{(sech({U}_{0}x))}^{{\tilde{\eta }}_{n}}{P}_{n}^{({\tilde{\eta }}_{n},{\tilde{\eta }}_{n})}(\mathrm{tanh}({U}_{0}x)), \end{aligned}\end{equation} \tag{ 88 }$$

where ${\mathcal{C}}_{2}$ and ${\mathcal{D}}_{2}$ are the respective normalization factors, and

$$\begin{equation}{\eta }_{n}={\nu }_{E}-n-1\enspace ,\qquad {\tilde{\eta }}_{n}={\tilde{\nu }}_{E}-n-1,\end{equation} \tag{ 89 }$$

together with

$$\begin{equation}{\nu }_{E}=\frac{1}{2}+\sqrt{\mathcal{A}({{\epsilon}}_{1}-E)+\frac{1}{4}}\enspace ,\qquad {\tilde{\nu }}_{E}=\frac{1}{2}+\sqrt{-\mathcal{A}({{\epsilon}}_{1}-E)+\frac{1}{4}}.\end{equation} \tag{ 90 }$$

On the other hand, from the relationship

$$\begin{equation}({{\epsilon}}_{1}-E)({{\epsilon}}_{2}-E+\delta \enspace {\gamma }_{0})=-{U}_{0}^{2}{\left(n+\frac{1}{2}-\sqrt{\delta \mathcal{A}({{\epsilon}}_{1}-E)+\frac{1}{4}}\right)}^{2}\enspace ,\quad \delta =+1,-1,\end{equation} \tag{ 91 }$$

we extract the energies of h₁ and h₂ after choosing δ = +1 and δ = −1, respectively.

It is worth to recall that an immediate solution for the energy equation (91) can be found for n = 0 and E = ₁ in both cases δ = ±1. Nevertheless, such a solution is discarded as it is not square-integrable.

Likewise in (67), we can obtain the energy intervals in which E takes real values for h₁ and h₂. First, we should guarantee that η_n and ${\tilde{\eta }}_{n}$ in (89) are real and positive in order to get square-integrable solutions. Next, the left term in (91) should be negative as the equation would have no real solutions otherwise. Combining both results, we get the intervals where the real roots of (91) have to lie,

$$\begin{equation}{h}_{1}:\enspace E\in \left({{\epsilon}}_{2}+{\gamma }_{0},{{\epsilon}}_{1}-\frac{n(n+1)}{\mathcal{A}}\right),\quad {{\epsilon}}_{1} > {{\epsilon}}_{2}+{\gamma }_{0},\end{equation} \tag{ 92 }$$

$$\begin{equation}{h}_{2}:\enspace E\in \left({{\epsilon}}_{2}+{\gamma }_{0},{{\epsilon}}_{1}-\frac{n(n+1)}{\mathcal{A}}\right),\quad {{\epsilon}}_{1}< {{\epsilon}}_{2}-{\gamma }_{0}.\end{equation} \tag{ 93 }$$

One cannot get square-integrable eigenstates of either h₁ or h₂ corresponding to real energies for other values of ₁. Now, the intervals (93) are nonempty for some values of n only. This way, we get an upper bounds n_max and ${\tilde{n}}_{\mathrm{max}}\enspace$ for the number of bound states of h₁ and h₂, respectively, that we can obtain this way. They are

$$\begin{align} {n}_{\mathrm{max}}& =\lfloor \sqrt{\mathcal{A}({{\epsilon}}_{1}-{{\epsilon}}_{2}-{\gamma }_{0})+\frac{1}{4}}-\frac{1}{2}\rfloor \enspace ,\quad {\tilde{n}}_{\mathrm{max}}=\lfloor \sqrt{\mathcal{A}({{\epsilon}}_{2}-{{\epsilon}}_{1}-{\gamma }_{0})+\frac{1}{4}}-\frac{1}{2}\rfloor . \end{align} \tag{ 94 }$$

The corresponding bispinors are given by

$$\begin{align} {\mathbf{\Psi }}_{n}& =\frac{1}{\sqrt{2}}\left(-i\frac{{\psi }_{2;n}^{\prime }}{E-{{\epsilon}}_{1}},{\psi }_{2;n},{\psi }_{2;n},-i\frac{{\psi }_{2;n}^{\prime }}{E-{{\epsilon}}_{1}}\right)\enspace ,\enspace \enspace {\mathbf{\Xi }}_{n}=\frac{1}{\sqrt{2}}\left(-i\frac{{\xi }_{2;n}^{\prime }}{E-{{\epsilon}}_{1}},{\xi }_{2;n},-{\xi }_{2;n},i\frac{{\xi }_{2;n}^{\prime }}{E-{{\epsilon}}_{1}}\right), \end{align} \tag{ 95 }$$

for ₁ > ₂ + γ₀ and ₁ < ₂ − γ₀, respectively, with ψ_2;n and χ_2;n given in (88).

To illustrate the results presented in this section, let us fix the parameters as ${U}_{0}=\mathcal{A}=1$, ₁ = 1.5, ₂ = −1.5, and γ₀ = 0.3. Since ₁ > ₂ + γ₀, we would expect bound states only for h₁, besides the non-physical solution E = ₁ = 1.5. Moreover, from (94), one may see that only two bound states can be generated. Such an information may be verified in table 2, where we obtain two physical energies for h₁ as E₀ = −0.685 308 and E₁ = −1.183 15. Although there are more real energies, they do not satisfy the finite-norm condition η_n > 0. The behavior for the corresponding probability densities ${\mathcal{P}}_{n}$ is depicted in figure 6.

Zoom In Zoom Out Reset image size

In the article, we focused on the systems described by Dirac Hamiltonians of the form (2) that appear in the analysis of bilayer Dirac materials. We were interested in analytical treatment of confined states that can appear due to local fluctuations (61), (72), (87) or periodicity defects (83) of the involved interactions.

We have made use of the fact that equation (4) is reducible in terms of equations (5) and (6) with 2 × 2 Hamiltonians. In section 2, we showed that both equations could be solved through the solutions of the Schrödinger equation (14), whose potential is a nonlinear function of the interactions and their derivatives. We focused on the specific case where (14) can be significantly simplified into a Schrödinger equation with energy-dependent potential (16).

We considered confinement by a combination of external magnetic field and mechanical deformations in section 3 where the energy-dependent Schrödinger equation was identified with the stationary equation of the harmonic oscillator or the Rosen–Morse system. In section 4, we focused on confinement by inhomogeneities of the on-site and inter-layer interactions. We showed that Dirac fermions can be confined by a local fluctuation or periodicity defect of the on-site interaction ₁. We demonstrated this fact on the systems with Pöschl–Teller-type interactions (61), (72), or periodic interactions with a localized defect (83). Finally, in section 5, we considered a situation where only the interlayer interaction was inhomogeneous. Here we fixed the trigonal warping term vanishing. It allowed us analytical treatment of decoupled equations with energy-dependent potential (86).

In all the scenarios, we faced the need to solve a Schrödinger equation with energy-dependent potentials. In section 3, the latter equation occurred due to the presence of the trigonal warping term, v₃ ≠ 0. If such a term were absent, decoupling of (9) would produce a Schrödinger equation with energy-independent potential. In sections 4 and 5, ₁(x), ₂(x) or γ(x) were considered inhomogeneous functions to account for the presence of an electrostatic component in the potential term of the reduced equations (5) and (6). It is known [37] that decoupling of components in the stationary equation for 2 × 2 Dirac Hamiltonian with electric potential leads to Schrödinger equation with energy-dependent potential. The problems related to the solution of such an equation are avoided when bound states with zero energy are of interest, see e.g. [38–41]. The zero modes in presence of inhomogeneous electric potential and an effective mass were discussed recently in [42]. In this context, the Hamiltonian (2) can be understood as two, coupled 2 × 2 Dirac Hamiltonians with electrostatic potential accompanied by an effective mass term. For these systems, we found localized states with energies distinct from zero.

The solvability of the presented models relies on the specific form of the interactions. As much as it is rather impossible to prepare such fields in the experiments, we can see a broader application of the analytically solvable models. Besides representing a technically feasible setting that demonstrates a specific physical phenomenon (confinement of Dirac fermions in our case), they can be used to analyze physically more relevant interactions via perturbative approach where they can serve as the initial (unperturbed) system. They can also provide reference models for application of numerical methods.

In the light of the computational complexities mentioned, complete analysis of the spectral properties of the presented models would necessarily involve the use of numerical methods. Our aim in this article was different. We demonstrated confinement of Dirac fermions in bilayer graphene by different types of interactions in terms of analytical solutions. Our analytic approach is not restricted to the problems addressed in this paper and can be extended to a broader class of interactions. Indeed, in the presented models, we used exactly solvable models (e.g. harmonic oscillator, Rosen–Morse model) to find localized solutions of (16). However, to demonstrate the confinement, it would be sufficient to consider broader class of systems, where one bound state is known analytically at least. It would be also possible to analyze existence of the confined states via qualitative (e.g. asymptotic) properties of the interactions, following the same steps. Nevertheless, it goes beyond the scope of the present article and it should be addressed in the future.

MC-C thanks Department of Physics of the Nuclear Physics Institute of CAS for hospitality. MC-C acknowledges the support of CONACYT, Project FORDECYT-PRONACES/61533/2020. MC-C also acknowledges the Conacyt fellowship 301117. VJ was supported by GAČR Grant No. 19-07117S. KZ acknowledges the support from the project 'Physicists on the move II' (KINEÓ II), Czech Republic, Grant No. CZ.02.2.69/0.0/0.0/18_053/0017163.

All data that support the findings of this study are included within the article (and any supplementary files).