Infinite bound states and hydrogen atom-like energy spectrum induced by a flat band

In the past decades, the physics induced by flat bands have attracted great interest [1–15]. In contrast with ordinary dispersion bands, a lot of novel phenomena, for example, the ferromagnetism transition [16], localization [17], super-Klein tunneling [18–22], quantum Hall-like states [23, 24], zitterbewegung [25], preformed pairs [26], strange metal [27], high-T_c superconductivity/superfluidity [28–31], etc, appear. Flat bands have some prominent noteworthy features. First of all, the states of the flat band have localized properties. For example, the states only occupy a finite number of unit cells in the lattice model [32]. Secondly, the flat band is very sensitive to small perturbations, e.g. the magnetism transitions driven by an interaction [33, 34]. Thirdly, the infinitely large density of states (DOS) can give rise to a linear dependence of interaction strength of superconductor/superfluid order parameters [35–39].

It is well known that the behaviors of DOS near the threshold of the continuous energy spectrum play crucial roles in the formations of bound states [40]. For example, in the three-dimensional case, the DOS near the threshold (E_th = 0) follows the law ${\rho }_{0}(E)\propto \sqrt{E}\to 0$ as E → 0, then only when the potential strength is sufficiently large [for a square well potential with width a > 0, the depth |V| > π² ℏ²/(8ma²)] [41], a bound state can exist. Meanwhile, for a two-dimensional counterpart, due to a finite DOS near the band edge, there exist bound states no matter how weak the attractive potential is. For the one-dimensional case, owing to the divergence of DOS near the band edge [${\rho }_{0}(E)\propto 1/\sqrt{E}\to \infty$ as E goes to zero], the bound states also exist for arbitrarily weak attractive potential. In addition, in the presence of a weak potential well, the shallow bound-state energy is proportional to the square of potential strength. Due to the singular density of DOS of the flat band, it is expected that the existence conditions of bound states will be affected significantly.

The bound-state problem in the two-dimensional flat band model with a short-ranged potential well and a long-ranged Coulomb potential has been investigated by Gorbar et al [42]. It was found that for arbitrary weak potential, there exist bound states which split from the flat band. In the presence of Coulomb potential, the flat band cannot survive. Furthermore, Van Pottelberge found that the flat band gradually evolves into a continuous band (uncountable set) with the increase in Coulomb potential strength [43].

In this work, we investigate the bound states in a one-dimensional spin-1 Dirac-type Hamiltonian. It is found that for an arbitrarily weak delta potential, similar to that in the two-dimensional case [42], there exists a bound state which is generated from the flat band. Differently from the ordinary one-dimensional case, the bound-state energy is linearly dependent on the potential strength as the strength goes to zero. It is believed that the existence of a hydrogen atom-like energy spectrum E_n ∝ 1/n² [41] usually requires a long-ranged Coulomb potential 1/r. Surprisingly for the flat band system, even a short-ranged square well potential can result in an infinite number of bound states (countable set). The existence of infinite bound states originates from an arbitrarily strong effective potential induced by the flat band. In addition, when the bound-state energy is much smaller than the typical energy scale of the system, it displays a hydrogen atom-like energy spectrum, i.e. E_n ∝ 1/n².

The work is organized as follows. In section 2, the three energy bands, free-particle Green function and DOS are investigated. Next, we solve the bound-state problems for several typical short-ranged potentials in section 3. At the end, a summary is given in section 4.

In this work, we consider a three-component spin-1 Dirac-type Hamiltonian [22] in one dimension, i.e.

$$\begin{align}\hfill H& ={H}_{0}+{V}_{\text{p}}(x)\hfill \\ \hfill {H}_{0}& =-i{v}_{\text{F}}\hslash {S}_{x}{\partial }_{x}+m{S}_{z},\hfill \end{align} \tag{ 1 }$$

where V_p(x) is the potential energy, H₀ is the free-particle Hamiltonian, v_F > 0 is the Fermi velocity, and m > 0 is the energy gap parameter. S_x and S_z are spin operators for spin-1 particles [44, 45], i.e.

$$\begin{equation}{S}_{x}=\left[\begin{matrix}\hfill 0\hfill & \hfill \frac{1}{\sqrt{2}}\hfill & \hfill 0\hfill \\ \hfill \frac{1}{\sqrt{2}}\hfill & \hfill 0\hfill & \hfill \frac{1}{\sqrt{2}}\hfill \\ \hfill 0\hfill & \hfill \frac{1}{\sqrt{2}}\hfill & \hfill 0\hfill \end{matrix}\right];\qquad {S}_{z}=\left[\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \end{matrix}\right],\end{equation} \tag{ 2 }$$

in the usual basis |i⟩ with i = 1, 2, 3. In the whole manuscript, we use the units of v_F = ℏ = 1. The above Hamiltonian may be realized in photonic systems [46, 47]. When V_p(x) = 0, for a given momentum k, the free-particle Hamiltonian H₀ has three eigenstates and the eigenenergies, i.e.

$$\begin{align}\hfill \langle x\vert -,k\rangle & ={\psi }_{-,k}(x)=\frac{1}{2\sqrt{{k}^{2}+{m}^{2}}}\left(\begin{matrix}\hfill \sqrt{{k}^{2}+{m}^{2}}-m\hfill \\ \hfill -\sqrt{2}k\hfill \\ \hfill \sqrt{{k}^{2}+{m}^{2}}+m\hfill \end{matrix}\right){\text{e}}^{\text{i}kx},\hfill \\ \hfill {E}_{-,k}& =-\sqrt{{k}^{2}+{m}^{2}};\hfill \\ \hfill \langle x\vert 0,k\rangle & ={\psi }_{0,k}(x)=\frac{1}{\sqrt{2({k}^{2}+{m}^{2})}}\left(\begin{matrix}\hfill -k\hfill \\ \hfill \sqrt{2}m\hfill \\ \hfill k\hfill \end{matrix}\right){\text{e}}^{\text{i}kx},\hfill \\ \hfill {E}_{0,k}& =0;\hfill \\ \hfill \langle x\vert +,k\rangle & ={\psi }_{+,k}(x)=\frac{1}{2\sqrt{{k}^{2}+{m}^{2}}}\left(\begin{matrix}\hfill \sqrt{{k}^{2}+{m}^{2}}+m\hfill \\ \hfill \sqrt{2}k\hfill \\ \hfill \sqrt{{k}^{2}+{m}^{2}}-m\hfill \end{matrix}\right){\text{e}}^{\text{i}kx},\hfill \\ \hfill {E}_{+,k}& =\sqrt{{k}^{2}+{m}^{2}}.\hfill \end{align} \tag{ 3 }$$

|−, 0, +; k⟩ denote the states for the lower, middle (flat band) and upper bands, respectively, and E_−,k , E_0,k , E_+,k represent the three corresponding energy bands. It is found that a flat band with zero energy (E_0,k = 0) appears between the upper and lower bands (see figure 1). The possible bound states (if any) can only exist in the gaps among the three bands, i.e. 0 < E < m and −m < E < 0 (the regions A and B in figure 1).

Zoom In Zoom Out Reset image size

In order to explore the properties of the flat band, we calculate the free-particle Green function G₀ (corresponding to H₀) and DOS ρ₀(x, E) [40]. In momentum space, the free-particle Green function can be expressed as

$$\begin{align}\hfill & {G}_{0}(k,{k}^{\prime },z) \equiv \langle k\vert \frac{1}{z-{H}_{0}}\vert {k}^{\prime }\rangle \hfill \\ \hfill & \qquad =\left[\begin{matrix}\hfill {G}_{0,11}(k,{k}^{\prime },z)\hfill & \hfill {G}_{0,12}(k,{k}^{\prime },z)\hfill & \hfill {G}_{0,13}(k,{k}^{\prime },z)\hfill \\ \hfill {G}_{0,21}(k,{k}^{\prime },z)\hfill & \hfill {G}_{0,22}(k,{k}^{\prime },z)\hfill & \hfill {G}_{0,23}(k,{k}^{\prime },z)\hfill \\ \hfill {G}_{0,31}(k,{k}^{\prime },z)\hfill & \hfill {G}_{0,32}(k,{k}^{\prime },z)\hfill & \hfill {G}_{0,33}(k,{k}^{\prime },z)\hfill \end{matrix}\right],\hfill \\ \hfill & \qquad =\frac{\delta (k-{k}^{\prime })}{2({m}^{2}z+{k}^{2}z-{z}^{3})}\hfill \\ \hfill & \quad \qquad \times \left[\begin{matrix}\hfill {k}^{2}-2z(m+z)\hfill & \hfill -\sqrt{2}k(m+z)\hfill & \hfill -{k}^{2}\hfill \\ \hfill -\sqrt{2}k(m+z)\hfill & \hfill 2({m}^{2}-{z}^{2})\hfill & \hfill \sqrt{2}k(m-z)\hfill \\ \hfill -{k}^{2}\hfill & \hfill \sqrt{2}k(m-z)\hfill & \hfill {k}^{2}+2z(m-z)\hfill \end{matrix}\right].\hfill \end{align} \tag{ 4 }$$

In coordinate space, the matrix elements of the free-particle Green function are [48]

$$\begin{equation}{G}_{0,ij}(x,{x}^{\prime },z)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}k\,\mathrm{d}{k}^{\prime }\,{\text{e}}^{\text{i}(kx-{k}^{\prime }{x}^{\prime })}{G}_{0,ij}(k,{k}^{\prime },z),\end{equation} \tag{ 5 }$$

or explicitly written down

$$\begin{align}\hfill {G}_{0,11}(x,{x}^{\prime },z)& =\frac{-(m+z)}{2\sqrt{{m}^{2}-{z}^{2}}}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert }\hfill \\ \hfill & \quad -\frac{\sqrt{{m}^{2}-{z}^{2}}}{4z}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert }+\frac{\delta (x-{x}^{\prime })}{2z},\hfill \\ \hfill {G}_{0,12}(x,{x}^{\prime },z)& =\frac{-i(m+z)\mathrm{sign}(x-{x}^{\prime })}{2\sqrt{2}z}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert },\hfill \\ \hfill {G}_{0,13}(x,{x}^{\prime },z)& =\frac{\sqrt{{m}^{2}-{z}^{2}}}{4z}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert }-\frac{\delta (x-{x}^{\prime })}{2z},\hfill \\ \hfill {G}_{0,21}(x,{x}^{\prime },z)& ={G}_{0,12}(x,{x}^{\prime },z),\hfill \\ \hfill {G}_{0,22}(x,{x}^{\prime },z)& =\frac{\sqrt{{m}^{2}-{z}^{2}}}{2z}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert },\hfill \\ \hfill {G}_{0,23}(x,{x}^{\prime },z)& =\frac{i(m-z)\mathrm{sign}(x-{x}^{\prime })}{2\sqrt{2}z}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert },\hfill \\ \hfill {G}_{0,31}(x,{x}^{\prime },z)& ={G}_{0,13}(x,{x}^{\prime },z),\hfill \\ \hfill {G}_{0,32}(x,{x}^{\prime },z)& ={G}_{0,23}(x,{x}^{\prime },z),\hfill \\ \hfill {G}_{0,33}(x,{x}^{\prime },z)& =\frac{(m-z)}{2\sqrt{{m}^{2}-{z}^{2}}}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert }\hfill \\ \hfill & \quad -\frac{\sqrt{{m}^{2}-{z}^{2}}}{4z}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert }+\frac{\delta (x-{x}^{\prime })}{2z},\hfill \end{align} \tag{ 6 }$$

where the sign function sign(x) = 1 for x > 0, sign(x) = −1 for x < 0 and δ(x) is a Dirac delta function. The Dirac delta term δ(x − x')/(2z) in the Green function arises from the flat band, which reflects the localized properties of states in the flat band [32].

The DOS at energy E (per unit length) is given by the imaginary part of the Green function [40], i.e.

$$\begin{align}\hfill {\rho }_{0}(x,E)& =-\underset{{x}^{\prime }\to x}{\mathrm{lim}}\,\frac{\text{Im}\left\{\mathrm{Tr}[{G}_{0}(x,{x}^{\prime },z=E+i{0}_{+})]\right\}}{\pi }\hfill \\ \hfill & =-\frac{\text{Im}\left\{{\sum }_{i=1,2,3}{G}_{0,ii}(x,x,z=E+i{0}_{+})\right\}}{\pi }\hfill \\ \hfill & =-\frac{\text{Im}\left[\frac{\delta (0)}{E+i{0}_{+}}+\frac{-(E+i{0}_{+})}{\sqrt{{m}^{2}-{(E+i{0}_{+})}^{2}}}\right]}{\pi }\hfill \\ \hfill & =\frac{\theta (E-m)\vert E\vert }{\pi \sqrt{{E}^{2}-{m}^{2}}}+\delta (0)\delta (E)+\frac{\theta (-E-m)\vert E\vert }{\pi \sqrt{{E}^{2}-{m}^{2}}},\hfill \end{align} \tag{ 7 }$$

where θ(x) = 1 if x > 0, otherwise θ(x) = 0 and δ(0) = lim_x'→x   δ(x − x') → ∞.

It shows that the flat band of zero energy gives rise to an infinitely large DOS [δ(0)] [42], which reflects the fact that there are infinite states of the flat band in the continuous model. Furthermore, the δ(E) singularity in DOS ρ₀(x, E) causes a 1/z singularity in the Green function (see equation (6) and appendix A). In addition, near the thresholds of the continuous spectrum (E_th = ±m), the DOS shows a similar divergence to that in the ordinary one-dimensional case, e.g. ${\rho }_{0}(x,E)\propto 1/\sqrt{\vert E-{E}_{\text{th}}\vert }\to \infty$ as E → E_th (see figure 1).

In the following, we will see that such an infinitely large DOS and 1/z singularity of the Green function have important effects on bound states (see appendix B). We may naively think that, due to the non-dispersion of the flat band, an arbitrarily weak potential may pull some zero-energy states out of the flat band, and then these states form bound states. It is shown that, indeed, for some kinds of square well potentials, there exist infinite bound states for arbitrarily weak potential (see the next section).

In the following manuscript, we assume the potential energy V_p has the following diagonal form in the usual basis |i = 1, 2, 3⟩, namely,

$$\begin{align}\hfill {V}_{\text{p}}(x)& ={V}_{11}(x)\otimes \vert 1\rangle \langle 1\vert +{V}_{22}(x)\otimes \vert 2\rangle \langle 2\vert \hfill \\ \hfill & \quad +{V}_{33}(x)\otimes \vert 3\rangle \langle 3\vert =\left[\begin{matrix}\hfill {V}_{11}(x)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {V}_{22}(x)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {V}_{33}(x)\hfill \end{matrix}\right].\hfill \end{align} \tag{ 8 }$$

3.1. Delta potential

First of all, we assume the potential energy has the form of the delta potential and satisfies V₁₁(x) = V₃₃(x) ≡ 0, and V₂₂(x) = gδ(x) with potential strength g. The Schrödinger equation (Hψ = Eψ) can be written in terms of a three-component wave function, i.e.

$$\begin{align}\hfill \begin{aligned}\hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E-m]{\psi }_{1}(x),\hfill \\ \hfill -i{\partial }_{x}[{\psi }_{1}(x)+{\psi }_{3}(x)]/\sqrt{2}& =[E-g\delta (x)]{\psi }_{2}(x),\hfill \\ \hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E+m]{\psi }_{3}(x).\hfill \end{aligned}\end{align} \tag{ 9 }$$

Furthermore, using equation (9) to eliminate wave functions for the first and third components, we get an effective Schrödinger equation (second-order differential equation) for ψ₂(x):

$$\begin{equation}{\partial }_{x}^{2}{\psi }_{2}(x)+\frac{[E-g\delta (x)]({E}^{2}-{m}^{2})}{E}{\psi }_{2}(x)=0.\end{equation} \tag{ 10 }$$

For bound states and x ≠ 0, the wave function ψ₂(x) can be written as

$$\begin{align}\hfill & {\psi }_{2}(x)=C{e}^{-\lambda x}\quad x > 0,\hfill \\ \hfill & {\psi }_{2}(x)=D{e}^{\lambda x}\quad x< 0,\hfill \end{align} \tag{ 11 }$$

where $\lambda =\sqrt{{m}^{2}-{E}^{2}} > 0$, C and D are two coefficients for the decayed part of the wave function. The continuity of ψ₂(x) at origin x = 0 requires C = D. Integrating both sides of equation (10) near the origin x = 0, it is shown that the derivative of ψ₂(x) satisfies

$$\begin{align}\hfill {\psi }_{2}^{\prime }({0}^{+})-{\psi }_{2}^{\prime }({0}^{-})& =\frac{g({E}^{2}-{m}^{2})}{E}{\psi }_{2}(0),\hfill \\ \hfill & \to -2\lambda C=\frac{g({E}^{2}-{m}^{2})}{E}C.\hfill \end{align} \tag{ 12 }$$

Based on equation (12), the bound-state energy E_B is obtained:

$$\begin{equation}{E}_{\text{B}}\equiv E=\frac{mg}{\sqrt{4+{g}^{2}}}.\end{equation} \tag{ 13 }$$

First of all, there always exists one bound state whether the potential is repulsive (g > 0) or attractive (g < 0) (see figure 2). However, for the ordinary one-dimensional case, only when the potential is attractive (g < 0), there is one bound state. Secondly, when the potential strength g goes to zero, the bound-state energy E_B is linearly proportional to the strength, i.e. E_B ∝ g as g → 0 [see panel (b) of figure 2]. Such behavior is very different from the ordinary one-dimensional bound-state energy, which is proportional to g². The linear dependence on the strength is also consistent with the fact that the superfluid/superconductor order parameter is linearly proportional to the two-body interaction strength in a flat band superfluid system [29, 35, 36]. Most of all these peculiar behaviors are due to the infinitely large DOS of the flat band (see appendix B).

Zoom In Zoom Out Reset image size

Thirdly, due to the restrictions of the thresholds of the upper and lower continuums, when the potential strength is very large (g → ∞), the bound-state energy E_B approaches the thresholds of continuums ±m asymptotically. This is also very different from the usual one-dimensional bound-state energy, which goes to infinity (∝g² → ∞).

3.2. Square well potential of type I

In this subsection, we assume the potential satisfies V₁₁(x) = V₂₂(x) = V₃₃(x) and

$$\begin{align}\hfill & {V}_{11}(x)=0\quad x > a/2,\hfill \\ \hfill & {V}_{11}(x)=V\quad -a/2< x< a/2,\hfill \\ \hfill & {V}_{11}(x)=0\quad x< -a/2,\hfill \end{align} \tag{ 14 }$$

where a > 0 and V are the width and depth (or height) of square well potential, respectively. Now the Schrödinger equation is

$$\begin{align}\hfill \begin{aligned}\hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E-m-{V}_{11}(x)]{\psi }_{1}(x),\hfill \\ \hfill -i{\partial }_{x}[{\psi }_{1}(x)+{\psi }_{3}(x)]/\sqrt{2}& =[E-{V}_{11}(x)]{\psi }_{2}(x),\hfill \\ \hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E+m-{V}_{11}(x)]\left.\!\right]{\psi }_{3}(x).\hfill \end{aligned}\end{align} \tag{ 15 }$$

Similarly to the above subsection, we get an effective second-order differential equation

$$\begin{align}\hfill & {\partial }_{x}^{2}{\psi }_{2}(x)+{\partial }_{x}\left[\mathrm{ln}\left(\frac{E-{V}_{11}(x)}{{(E-{V}_{11}(x))}^{2}-{m}^{2}}\right)\right]{\partial }_{x}{\psi }_{2}(x)\hfill \\ \hfill & \qquad +[{(E-{V}_{11}(x))}^{2}-{m}^{2}]{\psi }_{2}(x)=0.\hfill \end{align} \tag{ 16 }$$

In the above derivations, we assume that E − V(x) ≠ 0. For bound states and |x| > a/2 (V₁₁(x) ≡ 0), the wave function ψ₂(x) is exponentially decayed, i.e.

$$\begin{align}\hfill & {\psi }_{2}(x)=C{e}^{-\lambda x}\quad x > a/2,\hfill \\ \hfill & {\psi }_{2}(x)=D{e}^{\lambda x}\quad x< -a/2,\hfill \end{align} \tag{ 17 }$$

where $\lambda =\sqrt{{m}^{2}-{E}^{2}} > 0$. For |x| < a/2, the wave function ψ₂(x) can be

$$\begin{equation}{\psi }_{2}(x)=A\,{\text{e}}^{\text{i}kx}+B\,{\text{e}}^{-\text{i}kx},\end{equation} \tag{ 18 }$$

where $k=\sqrt{{(E-V)}^{2}-{m}^{2}} > 0$, A and B are two coefficients for the oscillation parts of the wave function.

In the following, we assume that the wave functions for all three components have finite values (they are not infinitely large). Integrating both sides of equation (15) near x = ±a/2, we find ψ₂(x) and ψ₁(x) + ψ₃(x) are continuous functions at x = ±a/2 [21, 45]. Furthermore, taking the relations of

$$\begin{align}\hfill & \frac{-i{\partial }_{x}{\psi }_{2}(x)}{\sqrt{2}[E-m-{V}_{11}(x)]}={\psi }_{1}(x),\hfill \\ \hfill & \frac{-i{\partial }_{x}{\psi }_{2}(x)}{\sqrt{2}[E+m-{V}_{11}(x)]}={\psi }_{3}(x),\hfill \end{align} \tag{ 19 }$$

into account, we get equations for coefficients A, B, C and D, i.e.

$$\begin{align}\hfill \begin{aligned}\hfill A\,{\text{e}}^{\text{i}ka/2}+B\,{\text{e}}^{-\text{i}ka/2}& =C{e}^{-\lambda a/2},\hfill \\ \hfill A\,{\text{e}}^{-\text{i}ka/2}+B\,{\text{e}}^{\text{i}ka/2}& =D{e}^{-\lambda a/2},\hfill \\ \hfill \frac{2(E-V)(ikA\,{\text{e}}^{\text{i}ka/2}-ikB\,{\text{e}}^{-\text{i}ka/2})}{(E-m-V)(E+m-V)}& =\frac{-2E\lambda C{e}^{-\lambda a/2}}{{E}^{2}-{m}^{2}},\hfill \\ \hfill \frac{2(E-V)(ikA\,{\text{e}}^{-\text{i}ka/2}-ikB\,{\text{e}}^{\mathrm{i}ka/2})}{(E-m-V)(E+m-V)}& =\frac{2E\lambda D{e}^{-\lambda a/2}}{{E}^{2}-{m}^{2}}.\hfill \end{aligned}\end{align} \tag{ 20 }$$

Letting the corresponding determinant vanish, we get an equation for the bound-state energy:

$$\begin{align}\hfill & [2{E}^{4}-4{E}^{3}V+2E{m}^{2}V-{m}^{2}{V}^{2}+2{E}^{2}(-{m}^{2}+{V}^{2})]\mathrm{sin}(ka)\hfill \\ \hfill & \quad +2E(V-E)\sqrt{({m}^{2}-{E}^{2})[{(E-V)}^{2}-{m}^{2}]}\mathrm{cos}(ka)=0.\hfill \end{align} \tag{ 21 }$$

Panel (b) of figure 2 reports the bound-state energy E_B as a function of potential strength V. In all the figures, we set the width of the potential well a = 1/m.

When |E − m| ≪ m, and |V| ≪ m, the resulting bound-state energy can be approximated by

$$\begin{equation}{E}_{\text{B}}=E\simeq \mp m\pm \frac{m{a}^{2}{V}^{2}}{2}.\end{equation} \tag{ 22 }$$

It is shown that near the thresholds of continuums ±m, and in the case of weak potential, the bound-state energy (relative to the thresholds) is proportional to the square of strength (E_B ± m ∝ V²), which is consistent with the ordinary one-dimensional bound-state energies (see figure 2).

This is because for the square well potential of type I, the three-diagonal matrix elements of potential are equal (in the usual basis |i = 1, 2, 3⟩) [see equation (14)]. One can also view the potential energy in the energy band basis, i.e. |−, k⟩, |0, k⟩ and |+, k⟩; consequently the diagonal form is also not changed approximately. Therefore, for the three energy bands (flat band, upper and lower bands), they have the same potentials V₁₁(x).

For the flat band, the states are spatially localized. Once the potential is turned on, the energy of the states within the range of the potential well only shift by a value of V, i.e. E = 0 → E = V [see the black dashed line in panel (b) of figure 2]. These resulting new states basically inherit the localized properties of the original flat band states. Therefore, we call them trivial bound states (localized flat band states), which we are not interested in.

In addition, near the thresholds of the upper and lower continuums, the DOS has a similar divergence behavior to that in the ordinary one-dimensional case, i.e. $1/\sqrt{\vert E-{E}_{\text{th}}\vert }\to \infty$ as E → E_th [40] [see figure 1 and equation (7)]. So, in such cases, the bound-state problem can be reduced to an ordinary one-dimensional bound-state problem near the thresholds (especially for shallow bound states). Consequently, the resulting bound-state energies are very similar to the usual one-dimensional results, e.g. E_B ± m ∝ V².

For a given potential strength V, the number of (non-trivial) bound states is always finite for a square well potential of type I (with three equal diagonal matrix elements).

3.3. Square well potential of type II

In this subsection, we assume the potential energy satisfies V₁₁(x) = V₃₃(x) ≡ 0 and

$$\begin{align}\hfill & {V}_{22}(x)=0\quad x > a/2,\hfill \\ \hfill & {V}_{22}(x)=V\quad -a/2< x< a/2,\hfill \\ \hfill & {V}_{22}(x)=0\quad x< -a/2.\hfill \end{align} \tag{ 23 }$$

Now the Schrödinger equation can be rewritten as

$$\begin{align}\hfill \begin{aligned}\hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E-m]{\psi }_{1}(x),\hfill \\ \hfill -i{\partial }_{x}[{\psi }_{1}(x)+{\psi }_{3}(x)]/\sqrt{2}& =[E-{V}_{22}(x)]{\psi }_{2}(x),\hfill \\ \hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E+m]{\psi }_{3}(x).\hfill \end{aligned}\end{align} \tag{ 24 }$$

Similarly to the above subsections, we get an effective Schrodinger equation (second-order differential equation):

$$\begin{equation}{\partial }_{x}^{2}{\psi }_{2}(x)+\frac{[E-{V}_{22}(x)]({E}^{2}-{m}^{2})}{E}{\psi }_{2}(x)=0.\end{equation} \tag{ 25 }$$

In the above derivation, we assume that E ≠ 0.

For bound states and |x| > a/2, the wave function ψ₂(x) can be written as

$$\begin{align}\hfill & {\psi }_{2}(x)=C{e}^{-\lambda x}\quad x > a/2,\hfill \\ \hfill & {\psi }_{2}(x)=D{e}^{\lambda x}\quad x< -a/2,\hfill \end{align} \tag{ 26 }$$

where $\lambda =\sqrt{{m}^{2}-{E}^{2}} > 0$. For |x| < a/2, the wave function ψ₂(x) is

$$\begin{equation}{\psi }_{2}(x)=A\,{\text{e}}^{\text{i}kx}+B\,{\text{e}}^{-\text{i}kx},\end{equation} \tag{ 27 }$$

where $k=\sqrt{\frac{[E-V]({E}^{2}-{m}^{2})}{E}} > 0$.

Similarly, integrating the both sides of equation (24) near x = ±a/2, we find ψ₂(x) and ψ₁(x) + ψ₃(x) are continuous functions at x = ±a/2. Furthermore, taking the relations

$$\begin{align}\hfill & \frac{-i{\partial }_{x}{\psi }_{2}(x)}{\sqrt{2}(E-m)}={\psi }_{1}(x),\hfill \\ \hfill & \frac{-i{\partial }_{x}{\psi }_{2}(x)}{\sqrt{2}(E+m)}={\psi }_{3}(x),\hfill \end{align} \tag{ 28 }$$

into account, it can be shown that the derivative ∂_x ψ₂(x) at x = ±a/2 is also continuous.

Using these continuity conditions of ψ₂(x) at x = ±a/2, we get linear equations for A, B, C and D

$$\begin{align}\hfill \begin{aligned}\hfill A\,{\text{e}}^{\mathrm{i}ka/2}+B\,{\text{e}}^{-\text{i}ka/2}& =C{e}^{-\lambda a/2},\hfill \\ \hfill A\,{\text{e}}^{-\text{i}ka/2}+B\,{\text{e}}^{\mathrm{i}ka/2}& =D{e}^{-\lambda a/2},\hfill \\ \hfill ikA\,{\text{e}}^{\mathrm{i}ka/2}-ikB\,{\text{e}}^{-\text{i}ka/2}& =-\lambda C{e}^{-\lambda a/2},\hfill \\ \hfill ikA\,{\text{e}}^{-\text{i}ka/2}-ikB\,{\text{e}}^{\mathrm{i}ka/2}& =\lambda D{e}^{-\lambda a/2}.\hfill \end{aligned}\end{align} \tag{ 29 }$$

In order to have non-trivial solutions to the above linear equations, the corresponding determinant should vanish, i.e.

$$\begin{align}\hfill Det\left\vert \begin{matrix}\hfill {\text{e}}^{\mathrm{i}ka/2}\hfill & \hfill {\text{e}}^{-\text{i}ka/2}\hfill & \hfill -{e}^{-\lambda a/2}\hfill & \hfill 0\hfill \\ \hfill {\text{e}}^{-\text{i}ka/2}\hfill & \hfill {\text{e}}^{\mathrm{i}ka/2}\hfill & \hfill 0\hfill & \hfill -{e}^{-\lambda a/2}\hfill \\ \hfill ik\,{\text{e}}^{\mathrm{i}ka/2}\hfill & \hfill -ik\,{\text{e}}^{-\text{i}ka/2}\hfill & \hfill \lambda {e}^{-\lambda a/2}\hfill & \hfill 0\hfill \\ \hfill ik\,{\text{e}}^{-\text{i}ka/2}\hfill & \hfill -ik\,{\text{e}}^{\mathrm{i}ka/2}\hfill & \hfill 0\hfill & \hfill -\lambda {e}^{-\lambda a/2}\hfill \end{matrix}\right\vert =0.\end{align} \tag{ 30 }$$

Then, we get an equation of bound-state energy

$$\begin{equation}2E\sqrt{V/E-1}\,\mathrm{cos}(ka)=(V-2E)\mathrm{sin}(ka),\end{equation} \tag{ 31 }$$

where $k=\sqrt{\frac{[E-V]({E}^{2}-{m}^{2})}{E}} > 0$. Taking a = 1/m, the resulting bound-state energies are reported in figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When |E| ≪ |V|, the above equation (31) is simplified to

$$\begin{equation}\mathrm{sin}\left(\sqrt{\frac{({m}^{2}-{E}^{2})V}{E}}a\right)=0.\end{equation} \tag{ 32 }$$

So the energy

$$\begin{equation}{E}_{n}\equiv {E}_{\text{B}}=E=\frac{-{n}^{2}{\pi }^{2}+\sqrt{{n}^{4}{\pi }^{4}+4{m}^{2}{a}^{4}{V}^{2}}}{2{a}^{2}V},\end{equation} \tag{ 33 }$$

where the quantum number n = 1, 2, 3, ... (see figure 3). It is shown that there exist infinite bound states, which are generated from the flat band. Due to the infinitely large DOS of the flat band, an arbitrarily small potential V indeed pulls some states out the flat band and then they form bound states. For a given n, when the potential is very strong (V → ±∞), the bound-state energy approaches the thresholds ±m, which is very similar to the delta potential case. For weak potential, the energy is proportional to the potential strength, e.g. E_n ∝ V.

We should emphasize that the existence of an infinite number of bound states is quite a universal phenomenon in the flat band system with a quite arbitrary potential well of type II (not limited to square well potential; for more explanations see appendix B). As long as the energy is small enough (very near the zero energy of the flat band), the resulting effective potential $(\tilde{V}\propto {V}_{22}/E\to \infty )$ will be arbitrarily strong, and then the system will have infinite bound states.

Furthermore, when |E| ≪ m, the above equation (32) becomes

$$\begin{equation}\mathrm{sin}\left(\sqrt{\frac{{m}^{2}V}{E}}a\right)=0,\end{equation} \tag{ 34 }$$

so we get the energy

$$\begin{equation}{E}_{n}={E}_{\text{B}}=E\simeq \frac{{m}^{2}{a}^{2}V}{{n}^{2}{\pi }^{2}},\end{equation} \tag{ 35 }$$

where the quantum number n = 1, 2, 3, ... (see figure 4). It is interesting that the bound-state energy displays a hydrogen atom-like energy spectrum. Outside the square potential well (|x| > a/2), the bound-state wave functions are all exponentially decayed [see equation (26)]. Within the potential well, i.e. |x| < a/2, the wave function oscillates wildly for large quantum n [see equation (27)]. When |E| → 0, the DOS can be approximately obtained:

$$\begin{equation}\text{DOS}(E)=\frac{{\Delta}n}{{\Delta}E}\simeq \vert \frac{\mathrm{d}n}{\mathrm{d}{E}_{n}}\vert =\frac{ma\vert V{\vert }^{1/2}}{2\pi \vert E{\vert }^{3/2}}.\end{equation} \tag{ 36 }$$

As |E| → 0, DOS(E) ∝ 1/|E|^3/2 → ∞.

The origin of existence of the infinite number of bound states in square well potential is due to the 1/z singularity in the Green function and the δ(E) singularity of the DOS induced by the flat band (see appendices A and B). For an ordinary bound-state problem in a quantum mechanics textbook, the energy E usually appears in the numerator of the square of the wave vector, e.g. k² ∼ E. Meanwhile, due to the 1/z of the Green function in the flat band model, the energy E can appear in the denominator, e.g. k² ∝ 1/E [see equations (34) and (25)]. This is why an infinite number of bound states and even the hydrogen atom-like energy spectrum for the square well potential of type II can appear, while for the square well potential of type I, the 1/z singularity does not appear in the bound-state equation [see equations (16) and (18)]. Consequently, the number of (non-trivial) bound states is always finite for a given potential strength V. In order to see the importance of the 1/z singularity, we give more explanations in the appendices A and B. For other types of square well potentials and general potential strengths, there also exist infinite bound states and a 1/n² energy spectrum (see appendix C).

In conclusion, we analytically calculate the free-particle Green function and DOS for the one-dimensional spin-1 Dirac model. It is found that the infinitely large DOS of the flat band results in a 1/z singularity of the Green function. Furthermore, we solve the bound-state problems in the flat band system for several typical potential wells. It is found that the flat band and its ensuing 1/z singularity of the Green function play crucial roles in the formations of bound states. Due to the 1/z singularity, the bound-state energy shows a linear dependence on the potential strength. In addition, there exists an infinite number of (non-trivial) bound states for some kinds of square well potential. Furthermore, the bound-state energy displays a hydrogen atom-like energy spectrum. Our findings provide alternative ways to get the hydrogen atom-like energy spectrum in condensed matter physics.

Finally, in the presence of long-ranged Coulomb potential, there also exist infinite bound states [55–57]. However, the bound-state energy is usually proportional to 1/n, which is different from 1/n² of the short-ranged square well potential.

The above findings provide some useful insights and aid our understanding of the many-body physics of the flat band. For example, the infinite bound states induced by small potential imply that an arbitrarily small interaction would dominate the physics, which also reflects the fact that the flat band is not unstable under interactions. Since the bound states can appear for both repulsive (positive strength) and attractive (negative strength) potentials, then one can expect that even a repulsive interaction may result in superfluid/superconductor pairing states in the flat band system [58].

There are also other interesting questions; for example, how to introduce phase shift for the states in the flat band, how the potential affects the phase shifts, what are the relationships between the phase shift and the number of bound states (Levinson theorem) [49–54], etc, need further investigation.

YC Zhang thanks Professor Shizhong Zhang for the useful discussions. This work was supported by the NSFC under Grant Nos. 11874127 and 11504095. YC Zhang is also thankful for the startup grant from Guangzhou University.

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

There exists a general relation between (the trace of) the Green function and DOS [40], i.e.

$$\begin{equation}\mathrm{Tr}[{G}_{0}(x,x,z)]={\int }_{-\infty }^{\infty }\mathrm{d}E\,\frac{{\rho }_{0}(x,E)}{z-E}.\end{equation} \tag{ A1 }$$

In addition, if the DOS has a δ(E) singularity, e.g. ρ₀(x, E) ∝ δ(E), the Green function behaves as

$$\begin{equation}\mathrm{Tr}[{G}_{0}(x,x,z)]={\int }_{-\infty }^{\infty }\mathrm{d}E\,\frac{{\rho }_{0}(x,E)}{z-E}\propto 1/z.\end{equation} \tag{ A2 }$$

Furthermore, for a general matrix element of the Green function, if a flat band appears, there also exists the 1/z singularity [see equation (6) in the main text]. It is known that the Green function can be represented as (spectral decomposition)

$$\begin{align}\hfill {G}_{0,ij}(x,{x}^{\prime },z)& =\langle x;i\vert \frac{1}{z-{H}_{0}}\vert {x}^{\prime };j\rangle \hfill \\ \hfill & =\sum\limits _{n,\lambda }\frac{{\psi }_{n,\lambda }(x;i){\psi }_{n,\lambda }^{\ast }({x}^{\prime };j)}{z-{E}_{n,\lambda }},\hfill \end{align} \tag{ A3 }$$

where n is the energy band index, λ denotes other quantum numbers for a given energy band n. Here the eigenfunction ψ_n,λ (x) satisfies H₀ ψ_n,λ (x) = E_n,λ ψ_n,λ (x). If a flat band appears, for example, E_n,λ ≡ 0 for some specific index n, G_0,ij (x, x', z) includes a term $\sim {\sum }_{\lambda }{\psi }_{n,\lambda }(x;i){\psi }_{n,\lambda }^{\ast }({x}^{\prime };j)/z$, then it would display a 1/z singularity in general. It should be remarked that the above conclusions not only apply to a free-particle Hamiltonian H₀, but also to a total Hamiltonian H = H₀ + V_p [for example, see equation (1)].

In order to see the importance of 1/z singularity of the Green function in bound-state equations, we reformulate the bound-state problem with the Green function method [40].

First of all, we introduce the full Green function G(z):

$$\begin{equation}G(z)=\frac{1}{z-H},\end{equation} \tag{ B1 }$$

which corresponds to the total Hamiltonian H = H₀ + V_p. In addition, G(z) satisfies the Lippmann–Shwinger equation, i.e.

$$\begin{equation}G(z)={G}_{0}(z)+{G}_{0}(z){V}_{\text{p}}G(z),\end{equation} \tag{ B2 }$$

where

$$\begin{equation}{G}_{0}(z)=\frac{1}{z-{H}_{0}},\end{equation} \tag{ B3 }$$

is the free-particle Green function in operator form. Taking a potential energy in the form of

$$\begin{equation}V(x)={V}_{22}(x)\otimes \vert 2\rangle \langle 2\vert =\left[\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {V}_{22}(x)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\end{equation} \tag{ B4 }$$

as an example (both the delta potential and square well potential of type II in section 3 belong to such a class), the Lippmann–Shwinger equation (in coordinate representation) takes the following form:

$$\begin{align}\hfill {G}_{ij}(x,{x}^{\prime },z)& ={G}_{0,ij}(x,{x}^{\prime },z)\hfill \\ \hfill & \quad +\int \mathrm{d}{x}_{1}\,{G}_{0,i2}(x,{x}_{1},z){V}_{22}({x}_{1}){G}_{2j}({x}_{1},{x}^{\prime },z).\hfill \end{align} \tag{ B5 }$$

Specifically, the matrix element

$$\begin{align}\hfill {G}_{22}(x,{x}^{\prime },z)& ={G}_{0,22}(x,{x}^{\prime },z)\hfill \\ \hfill & \quad +\int \mathrm{d}{x}_{1}\,{G}_{0,22}(x,{x}_{1},z){V}_{22}({x}_{1}){G}_{22}({x}_{1},{x}^{\prime },z).\hfill \end{align} \tag{ B6 }$$

The bound-state solutions are determined by the corresponding homogenous integral equation, namely

$$\begin{equation}{G}_{22}(x,{x}^{\prime },z)=\int \mathrm{d}{x}_{1}\,{G}_{0,22}(x,{x}_{1},z){V}_{22}({x}_{1}){G}_{22}({x}_{1},{x}^{\prime },z).\end{equation} \tag{ B7 }$$

For simplification, taking a specific value x' = 0, and setting z = E and G(x) ≡ G₂₂(x, 0, E), the resulting integral equation is

$$\begin{equation}G(x)=\int \mathrm{d}{x}_{1}\,{G}_{0,22}(x,{x}_{1},E){V}_{22}({x}_{1})G({x}_{1}).\end{equation} \tag{ B8 }$$

Inserting ${G}_{0,22}(x,{x}^{\prime },z)=\frac{\sqrt{{m}^{2}-{z}^{2}}}{2z}{e}^{-\sqrt{{m}^{2}-{z}^{2}}\vert x-{x}^{\prime }\vert }$ [see equation (6) in the main text] in equation (B8), we get

$$\begin{equation}G(x)=\frac{\sqrt{{m}^{2}-{E}^{2}}}{2E}\int \mathrm{d}{x}_{1}\,{e}^{-\sqrt{{m}^{2}-{E}^{2}}\vert x-{x}_{1}\vert }{V}_{22}({x}_{1})G({x}_{1}).\end{equation} \tag{ B9 }$$

It is shown that the 1/z singularity of the Green function matrix element G_0,22 introduces a factor 1/E, which results in the linear dependence on the potential strength in bound-state energy. For example, when the potential is a delta function, i.e. V₂₂(x) = gδ(x), then the integral equation becomes

$$\begin{equation}G(x)=\frac{g\sqrt{{m}^{2}-{E}^{2}}}{2E}{e}^{-\sqrt{{m}^{2}-{E}^{2}}\vert x\vert }G(0).\end{equation} \tag{ B10 }$$

Then

$$\begin{equation}G(0)=\frac{g\sqrt{{m}^{2}-{E}^{2}}}{2E}G(0),\end{equation} \tag{ B11 }$$

and

$$\begin{equation}1=\frac{g\sqrt{{m}^{2}-{E}^{2}}}{2E}.\end{equation} \tag{ B12 }$$

So we can get the bound-state energy

$$\begin{equation}{E}_{\text{B}}=E=\frac{mg}{\sqrt{4+{g}^{2}}},\end{equation} \tag{ B13 }$$

which is consistent with equation (13) in the main text.

In the case of square well potential, the above equation (B9) is essentially equivalent to the effective Schrödinger equation [equation (25) in the main text]. This is because the differential equation equation (25)

$$\begin{equation}{\partial }_{x}^{2}{\psi }_{2}(x)+\frac{[E-{V}_{22}(x)]({E}^{2}-{m}^{2})}{E}{\psi }_{2}(x)=0,\end{equation} \tag{ B14 }$$

can be also transformed into an corresponding integral equation for the bound-state problem. First of all, let us write it in the form of an ordinary second-order Schrödinger equation, e.g.

$$\begin{equation}-{\partial }_{x}^{2}{\psi }_{2}(x)+\tilde{V}{\psi }_{2}(x)=\tilde{E}{\psi }_{2}(x),\end{equation} \tag{ B15 }$$

where the effective total energy $\tilde{E}$, effective potential energy $\tilde{V}$

$$\begin{align}\hfill \tilde{E}& ={E}^{2}-{m}^{2},\hfill \\ \hfill \tilde{V}(x)& =\tilde{E}-\frac{[E-{V}_{22}(x)]({E}^{2}-{m}^{2})}{E}\hfill \\ \hfill & =\frac{{V}_{22}(x)({E}^{2}-{m}^{2})}{E},\hfill \end{align} \tag{ B16 }$$

and $\tilde{E}< 0$ for bound states. In the determining effective total energy, we assume the potential energy V₂₂(x) goes to zero as x → ±∞. It shows that when V/E > 0, the effective potential $\tilde{V}=V({E}^{2}-{m}^{2})/E< 0$ is attractive, and then bound states would exist in one dimension. Furthermore, if the energy is very small, e.g. 0 < E/V ≪ 1 and |E|/m ≪ 1, the effective potential strength can be arbitrarily large $(\vert \tilde{V}\vert \simeq \vert V\vert m/\vert E\vert \to \infty )$, so there would exist infinite bound states. We see that the factor 1/E plays crucial roles in the existence of an infinite number of bound states.

It should be remarked that the existence of infinite bound states is quite a universal phenomenon in the flat band system with a potential well of type II. This is because as long as the energy E is sufficiently small, the factor 1/E would result in an arbitrarily strong effective potential $\tilde{V}(x)$. Consequently, a strong effective potential would support infinite bound states.

The integral equation corresponding to the above equation (B15) is

$$\begin{equation}{\psi }_{2}(x)=\int \mathrm{d}{x}_{1}\,{\tilde{G}}_{0}(x,{x}_{1},\tilde{E})\tilde{V}({x}_{1}){\psi }_{2}({x}_{1}),\end{equation} \tag{ B17 }$$

where the effective free-particle Green function (corresponding to the effective Hamiltonian ${\tilde{H}}_{0}=-{\partial }_{x}^{2}$ in one dimension) [40] is

$$\begin{align}\hfill {\tilde{G}}_{0}(x,{x}_{1},\tilde{E})& =\frac{-{e}^{-\sqrt{-\tilde{E}}\vert x-{x}_{1}\vert }}{2\sqrt{-\tilde{E}}}\hfill \\ \hfill & =\frac{-{e}^{-\sqrt{{m}^{2}-{E}^{2}}\vert x-{x}_{1}\vert }}{2\sqrt{{m}^{2}-{E}^{2}}}.\hfill \end{align} \tag{ B18 }$$

Inserting it in equation (B17), then

$$\begin{equation}{\psi }_{2}(x)=\frac{\sqrt{{m}^{2}-{E}^{2}}}{2E}\int \mathrm{d}{x}_{1}\,{e}^{-\sqrt{{m}^{2}-{E}^{2}}\vert x-{x}_{1}\vert }{V}_{22}({x}_{1}){\psi }_{2}({x}_{1}).\end{equation} \tag{ B19 }$$

It is found that G(x) and ψ₂(x) satisfy the same integral equation for bound states [see equations (B9) and (B19)]. During the above derivations, we see that 1/z singularity in the Green function is very important in the bound-state equation. The 1/z singularity has close relations to the factor E in the denominator of the effective potential $\tilde{V}$ [see equations (B16), (B9) and (B19)]. Most of the peculiar behaviors of bound-state energies, for example linear dependence on the potential strength, an infinite number of bound states and a hydrogen atom-like spectrum, can be attributed to the 1/z singularity of the Green function.

In this appendix, we discuss the bound states for other types of square well potentials. First of all, we assume all the matrix elements V₁₁(x), V₂₂(x) and V₂₂(x) are square well potentials with same width a, i.e.

$$\begin{align}\hfill & {V}_{11}(x)=0\quad x > a/2,\hfill \\ \hfill & {V}_{11}(x)={V}_{11}\quad -a/2< x< a/2,\hfill \\ \hfill & {V}_{11}(x)=0\quad x< -a/2,\hfill \\ \hfill & {V}_{22}(x)=0\quad x > a/2,\hfill \\ \hfill & {V}_{22}(x)={V}_{22}\quad -a/2< x< a/2,\hfill \\ \hfill & {V}_{22}(x)=0\quad x< -a/2,\hfill \\ \hfill & {V}_{33}(x)=0\quad x > a/2,\hfill \\ \hfill & {V}_{33}(x)={V}_{33}\quad -a/2< x< a/2,\hfill \\ \hfill & {V}_{33}(x)=0\quad x< -a/2,\hfill \end{align} \tag{ C1 }$$

where V₁₁, V₂₂ and V₃₃ are the potential strengths of V₁₁(x), V₂₂(x) and V₂₂(x), respectively.

Now the Schrödinger equation can be rewritten as

$$\begin{align}\hfill \begin{aligned}\hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E-m-{V}_{11}(x)]{\psi }_{1}(x),\hfill \\ \hfill -i{\partial }_{x}[{\psi }_{1}(x)+{\psi }_{3}(x)]/\sqrt{2}& =[E-{V}_{22}(x)]{\psi }_{2}(x),\hfill \\ \hfill -i{\partial }_{x}{\psi }_{2}(x)/\sqrt{2}& =[E+m-{V}_{33}(x)]{\psi }_{3}(x).\hfill \end{aligned}\end{align} \tag{ C2 }$$

Similarly, eliminating ψ₁(x) and ψ₃(x), we get

$$\begin{align}\hfill & {\partial }_{x}^{2}{\psi }_{2}(x)+{\partial }_{x}{\psi }_{2}(x){\partial }_{x}\hfill \\ \hfill & \quad \times \left[\mathrm{ln}\left(\frac{2E-[{V}_{11}(x)+{V}_{33}(x)]}{[E-m-{V}_{11}(x)][E+m-{V}_{33}(x)]}\right)\right]\hfill \\ \hfill & \quad +\frac{2(E-{V}_{22}(x))(E+m-{V}_{33})[E-m-{V}_{11}(x)]}{2E-[{V}_{11}(x)+{V}_{33}(x)]}{\psi }_{2}(x)=0.\hfill \end{align} \tag{ C3 }$$

When |x| > a/2, the bound-state wave function ψ₂(x) can be written as

$$\begin{align}\hfill & {\psi }_{2}(x)=C{e}^{-\lambda x}\quad x > a/2,\hfill \\ \hfill & {\psi }_{2}(x)=D{e}^{\lambda x}\quad x< -a/2,\hfill \end{align} \tag{ C4 }$$

where $\lambda =\sqrt{{m}^{2}-{E}^{2}} > 0$. For |x| < a/2, the wave function ψ₂(x) is

$$\begin{equation}{\psi }_{2}(x)=A\,{\text{e}}^{\text{i}kx}+B\,{\text{e}}^{-\text{i}kx},\end{equation} \tag{ C5 }$$

where

$$\begin{equation}k=\sqrt{\frac{2(E-{V}_{22})(E-m-{V}_{11})[E-m-{V}_{33}]}{2E-({V}_{11}+{V}_{33})}} > 0.\end{equation} \tag{ C6 }$$

Integrating both sides of equation (C2) near x = ±a/2, we find ψ₂(x) and ψ₁(x) + ψ₃(x) are continuous functions near x = ±a/2. Furthermore, taking the relations

$$\begin{align}\hfill & \frac{-i{\partial }_{x}{\psi }_{2}(x)}{\sqrt{2}[E-m-{V}_{11}(x)]}={\psi }_{1}(x),\hfill \\ \hfill & \frac{-i{\partial }_{x}{\psi }_{2}(x)}{\sqrt{2}[E+m-{V}_{33}(x)]}={\psi }_{3}(x),\hfill \end{align} \tag{ C7 }$$

into account, we get

$$\begin{align}\hfill & A\,{\text{e}}^{\mathrm{i}ka/2}+B\,{\text{e}}^{-\text{i}ka/2}=C{e}^{-\lambda a/2},\hfill \\ \hfill & A\,{\text{e}}^{-\text{i}ka/2}+B\,{\text{e}}^{\mathrm{i}ka/2}=D{e}^{-\lambda a/2},\hfill \\ \hfill & \frac{[2E-({V}_{11}+{V}_{33})](ikA\,{\text{e}}^{\mathrm{i}ka/2}-ikB\,{\text{e}}^{-\text{i}ka/2})}{(E-m-{V}_{11})(E+m-{V}_{33})}\hfill \\ \hfill & \qquad =\frac{-2E\lambda C{e}^{-\lambda a/2}}{{E}^{2}-{m}^{2}},\hfill \\ \hfill & \frac{[2E-({V}_{11}+{V}_{33})](ikA\,{\text{e}}^{-\text{i}ka/2}-ikB\,{\text{e}}^{\mathrm{i}ka/2})}{(E-m-{V}_{11})(E+m-{V}_{33})}\hfill \\ \hfill & \qquad =\frac{2E\lambda D{e}^{-\lambda a/2}}{{E}^{2}-{m}^{2}}.\hfill \end{align} \tag{ C8 }$$

Letting the corresponding determinant vanish, we can get bound-state energies. When V₁₁ = V₂₂ = V₃₃ = V, the bound-state energy is reduced to equation (21) of potential of type I. Meanwhile, when V₁₁ = V₃₃ ≡ 0, V₂₂ = V, the bound-state energy is given by equation (31) of type II. We see for a given potential strength V, there exist finite bound states in the case of type I, while the system has infinite bound states in the case of type II.

From the discussions in appendix B, we know that a crucial requirement for the existence of infinite bound states is that the effective potential [or wave vector k in equation (C6)] can become infinitely large. So if the denominator of k² [see equation (C6)] cannot be canceled by the numerator, the wave vector k and effective potential can be arbitrarily large for some E. Then, there would exist infinite bound states. However, when the denominator of k² can be canceled by the its numerator, the wave vector k is always finite. In such a case, the system only has finite bound states, e.g. the case of V₁₁ = V₂₂ = V₃₃ = V (the case of type I).

In the following, we investigate the case of V₂₂ = V₃₃ ≡ 0, V₁₁ = V, where the denominator of k² cannot be canceled by its numerator. Then it is expected that there would exist infinite bound states. In the whole manuscript, we refer to such type of potential as type III. The bound-state energy equation is

$$\begin{align}\hfill & 2(-2E+V)\sqrt{\frac{m-E}{m+E}}k\,\text{cos}(ka)\hfill \\ \hfill & \qquad +(4{E}^{2}-4Em-3EV+mV)\mathrm{sin}(ka)=0,\hfill \end{align} \tag{ C9 }$$

where $k=\sqrt{\frac{2E(E+m)[E-m-V]}{2E-V}} > 0$. The results are reported in figures 5 and 6. When V < 0 and E → E_th = m, the resulting bound-state energy

$$\begin{equation}{E}_{\text{B}}\simeq m-\frac{m{a}^{2}{V}^{2}}{2},\end{equation} \tag{ C10 }$$

which is also consistent with result of one-dimensional bound states [see figures 2 and 5 and equation (22)].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, when the energy E approaches V/2, e.g. 2E − V → 0, the wave vector k → ∞. Consequently, there exist infinite bound states near the line E = V/2 (see black line of figure 5). When |E − V/2|/m ≪ 1, the above equation (C9) is simplified into

$$\begin{align}\hfill & \mathrm{sin}\left(\sqrt{\frac{2E(E+m)[E-m-V]}{2E-V}}a\right)\hfill \\ \hfill & \simeq \mathrm{sin}\left(\sqrt{\frac{-V{(m+V/2)}^{2}}{2E-V}}a\right)=0.\hfill \end{align} \tag{ C11 }$$

Then we get the energy

$$\begin{equation}{E}_{\text{B}}=E\simeq \frac{V}{2}-\frac{V{(m+V/2)}^{2}{a}^{2}}{2{n}^{2}{\pi }^{2}},\end{equation} \tag{ C12 }$$

where n = 1, 2, 3, ... There also exist infinite bound states for arbitrarily small potential strength V (see figure 5). The bound-state energy (relative to V/2) also displays the hydrogen atom-like energy spectrum (see figure 6).

The above discussions indicate that, for other types of square well potentials and general potential strengths, there also exist infinite bound states and a 1/n² energy spectrum.