Two-body problem for two-dimensional electrons in Bernervig-Hughes-Zhang model

TL;DR

This paper investigates how two interacting electrons can form bound states in a two-dimensional topological insulator model, revealing that band inversion can enhance electron pairing and affect the bound state spectrum.

Contribution

It provides a detailed analysis of two-electron bound states in the Bernevig-Hughes-Zhang model, highlighting the influence of topological phases and band inversion on pairing mechanisms.

Findings

Bound states can form within the band gap due to electron interactions.

Band inversion in topological phases can promote electron pairing.

The bound state spectrum varies significantly with the coupling parameter.

Abstract

We study the two-body problem for two-dimensional electron systems in a symmetrized Bernevig-Hughes-Zhang model which is widely used to describe topological and conventional insulators. The main result is that two interacting electrons can form bound states with the energy in the gap of the band spectrum. The pairing mechanism can be interpreted as the formation of a negative reduced effective mass of two electrons. The problem is complicated because the relative motion of the electrons is coupled to the center-of-mass motion. We consider the case of zero total momentum. Detail calculations are carried out for the repulsive interaction potential of steplike form. The states are classified according to their spin structure and two-particle basis functions that form a given bound state. We analyze the spectra and electronic structure of the bound states in the case of both topological and trivial phases and especially focus on effects originating from the band inversion and the coupling of the electron and hole bands. In the trivial phase and the topological phase with the large coupling parameter $a$, the bound state spectra are qualitatively similar. However, when $a$ is less a certain value, the situation changes dramatically. In the topological phase, new states arise with a higher binding energy at lower interaction potential, which evidences that the band inversion can favor pairing the electrons.

Equations

$${\hat{H}}_0(\hat{\mathbf{k}})=\left(\begin{array}{cc}\hat{h}(\hat{\mathbf{k}})& 0\\ 0& {\hat{h}}^{\ast }(-\hat{\mathbf{k}})\end{array}\right)$$

$$\hat{h}(\hat{\mathbf{k}})=\left(\begin{array}{cc}M-B{\hat{k}}^2& A({\hat{k}}_x+i{\hat{k}}_y)\\ A({\hat{k}}_x-i{\hat{k}}_y)& -M+B{\hat{k}}^2\end{array}\right),$$

$$\begin{gathered} (|E\uparrow ,E\uparrow ⟩,|E\uparrow ,H\uparrow ⟩,|E\uparrow ,E\downarrow ⟩,|E\uparrow ,H\downarrow ⟩, \\ |E\downarrow ,H\uparrow ⟩,\dots ,{|H\downarrow ,E\downarrow ⟩,|H\downarrow ,H\downarrow ⟩)}^T. \end{gathered}$$

$$\hat{H}(1,2)={\hat{H}}_0({\hat{\mathbf{k}}}_1)\oplus {\hat{H}}_0({\hat{\mathbf{k}}}_2)+V({\mathbf{r}}_1-{\mathbf{r}}_2)\cdot {\hat{\mathbf{I}}}_{16\times 16},$$

$$\hat{H}(1,2)\Psi ({\mathbf{r}}_1,{\mathbf{r}}_2)=E\Psi ({\mathbf{r}}_1,{\mathbf{r}}_2).$$

$$\begin{array}{rl}& {\Psi }_1(1,2)=\left(\begin{array}{c}{\psi }_1\cdot |E\uparrow E\uparrow ⟩\\ {\psi }_2\cdot |E\uparrow H\uparrow ⟩\\ {\psi }_5\cdot |H\uparrow E\uparrow ⟩\\ {\psi }_6\cdot |H\uparrow H\uparrow ⟩\end{array}\right),{\Psi }_2(1,2)=\left(\begin{array}{c}{\psi }_3\cdot |E\uparrow E\downarrow ⟩\\ {\psi }_4\cdot |E\uparrow H\downarrow ⟩\\ {\psi }_7\cdot |H\uparrow E\downarrow ⟩\\ {\psi }_8\cdot |H\uparrow H\downarrow ⟩\end{array}\right),\\ & {\Psi }_3(1,2)=\left(\begin{array}{c}{\psi }_9\cdot |E\downarrow E\uparrow ⟩\\ {\psi }_{10}\cdot |E\downarrow H\uparrow ⟩\\ {\psi }_{13}\cdot |H\downarrow E\uparrow ⟩\\ {\psi }_{14}\cdot |H\downarrow H\uparrow ⟩\end{array}\right),{\Psi }_4(1,2)=\left(\begin{array}{c}{\psi }_{11}\cdot |E\downarrow E\downarrow ⟩\\ {\psi }_{12}\cdot |E\downarrow H\downarrow ⟩\\ {\psi }_{15}\cdot |H\downarrow E\downarrow ⟩\\ {\psi }_{16}\cdot |H\downarrow H\downarrow ⟩\end{array}\right).\end{array}$$

$$\left\{\hat{H}_{j}-\bigl{[}\varepsilon-2v(\mathbf{r})\bigr{]}\mathbf{I}_{4\times 4}\right\}\Psi_{j}(\mathbf{r},\mathbf{R})=0.$$

$$\epsilon =\frac{E}{|M|},{\mathbf{r}}^{\prime }=\mathbf{r}\sqrt{\frac{|M|}{|B|}},a=\frac{A}{\sqrt{|MB|}},v({\mathbf{r}}^{\prime })=\frac{V(\mathbf{r})}{2|M|}.$$

$${\Psi }_j(\mathbf{R},\mathbf{r})={\Psi }_{j,\mathbf{K}}(\mathbf{r})e^{i\mathbf{KR}},$$

$$\left\{\begin{array}[]{rl}\left(-\dfrac{\varepsilon}{2}+\lambda+\hat{k}^{2}+\dfrac{K^{2}}{4}+v(r)\right)\psi_{1}(\mathbf{r})-\dfrac{a}{2}\left(\hat{k}_{+}-\dfrac{K_{+}}{2}\right)\psi_{2}(\mathbf{r})+\dfrac{a}{2}\left(\hat{k}_{+}+\dfrac{K_{+}}{2}\right)\psi_{5}(\mathbf{r})&=0\\ -\dfrac{a}{2}\left(\hat{k}_{-}-\dfrac{K_{-}}{2}\right)\psi_{1}(\mathbf{r})+\left(-\dfrac{\varepsilon}{2}+\mathbf{\hat{k}K}+v(r)\right)\psi_{2}(\mathbf{r})+\dfrac{a}{2}\left(\hat{k}_{+}+\dfrac{K_{+}}{2}\right)\psi_{6}(\mathbf{r})&=0\\ \dfrac{a}{2}\left(\hat{k}_{-}+\dfrac{K_{-}}{2}\right)\psi_{1}(\mathbf{r})+\left(-\dfrac{\varepsilon}{2}-\mathbf{\hat{k}K}+v(r)\right)\psi_{5}(\mathbf{r})-\dfrac{a}{2}\left(\hat{k}_{+}-\dfrac{K_{+}}{2}\right)\psi_{6}(\mathbf{r})&=0\\ \dfrac{a}{2}\left(\hat{k}_{-}+\dfrac{K_{-}}{2}\right)\psi_{2}(\mathbf{r})-\dfrac{a}{2}\left(\hat{k}_{-}-\dfrac{K_{-}}{2}\right)\psi_{5}(\mathbf{r})+\left(-\dfrac{\varepsilon}{2}-\lambda-\hat{k}^{2}-\dfrac{K^{2}}{4}+v(r)\right)\psi_{6}(\mathbf{r})&=0\,.\end{array}\right.$$

$$\left\{\begin{array}[]{rl}\left(-\dfrac{\varepsilon}{2}+\lambda+\hat{k}^{2}+\dfrac{K^{2}}{4}+v(r)\right)\psi_{3}(\mathbf{r})+\dfrac{a}{2}\left(\hat{k}_{-}-\dfrac{K_{-}}{2}\right)\psi_{4}(\mathbf{r})+\dfrac{a}{2}\left(\hat{k}_{+}+\dfrac{K_{+}}{2}\right)\psi_{7}(\mathbf{r})&=0\\ \dfrac{a}{2}\left(\hat{k}_{+}-\dfrac{K_{+}}{2}\right)\psi_{3}(\mathbf{r})+\left(-\dfrac{\varepsilon}{2}+\mathbf{\hat{k}K}+v(r)\right)\psi_{4}(\mathbf{r})+\dfrac{a}{2}\left(\hat{k}_{+}+\dfrac{K_{+}}{2}\right)\psi_{8}(\mathbf{r})&=0\\ \dfrac{a}{2}\left(\hat{k}_{-}+\dfrac{K_{-}}{2}\right)\psi_{3}(\mathbf{r})+\left(-\dfrac{\varepsilon}{2}-\mathbf{\hat{k}K}+v(r)\right)\psi_{7}(\mathbf{r})+\dfrac{a}{2}\left(\hat{k}_{-}-\dfrac{K_{-}}{2}\right)\psi_{8}(\mathbf{r})&=0\\ \dfrac{a}{2}\left(\hat{k}_{-}+\dfrac{K_{-}}{2}\right)\psi_{4}(\mathbf{r})+\dfrac{a}{2}\left(\hat{k}_{+}-\dfrac{K_{+}}{2}\right)\psi_{7}(\mathbf{r})+\left(-\dfrac{\varepsilon}{2}-\lambda-\hat{k}^{2}-\dfrac{K^{2}}{4}+v(r)\right)\psi_{8}(\mathbf{r})&=0\,.\end{array}\right.$$

$$\left\{\begin{array}[]{rl}\left[2v(r)\!-\!\varepsilon\!+\!2\lambda\!+\!2\hat{k}^{2}\right]\psi_{3}(\mathbf{r})\!+\!a\hat{k}_{-}\psi_{4}(\mathbf{r})\!+\!a\hat{k}_{+}\psi_{7}(\mathbf{r})&\!=0\\ a\hat{k}_{+}\psi_{3}(\mathbf{r})+\left[2v(r)-\varepsilon\right]\psi_{4}(\mathbf{r})+a\hat{k}_{+}\psi_{8}(\mathbf{r})&\!=0\\ a\hat{k}_{-}\psi_{3}(\mathbf{r})+\left[2v(r)-\varepsilon\right]\psi_{7}(\mathbf{r})+a\hat{k}_{-}\psi_{8}(\mathbf{r})&\!=0\\ a\hat{k}_{-}\psi_{4}(\mathbf{r})\!+\!a\hat{k}_{+}\psi_{7}(\mathbf{r})\!+\!\left[2v(r)\!-\!\varepsilon\!-\!2\lambda\!-\!2\hat{k}^{2}\right]\psi_{8}(\mathbf{r})&\!=0\,.\end{array}\right.$$

$${\Psi }_{2,K=0}(\mathbf{r})=\sum _m{\Psi }_{2m}{e}^{im\phi }=\sum _m\left(\begin{array}{c}{\psi }_{3m}(r)\\ {\psi }_{4m}(r)e^{i\phi }\\ {\psi }_{7m}(r)e^{-i\phi }\\ {\psi }_{8m}(r)\end{array}\right)e^{im\phi }.$$

$$\begin{gathered} {\hat{k}}_m^4{\Phi }_m+[2\lambda +a^2-g_2(\epsilon ,r)]{\hat{k}}_m^2{\Phi }_m-2g_1(\epsilon ,r)\frac{d}{dr}\left({\hat{k}}_m^2{\Phi }_m\right) \\ -(2\lambda +a^2)g_1(\epsilon ,r)\frac{d}{dr}{\Phi }_m+[1-\tilde{\epsilon }(r{)}^2-\lambda g_2(\epsilon ,r)]{\Phi }_m=0, \end{gathered}$$

$${\hat{k}}_m^2=-\frac{d^2}{d{r}^2}-\frac{1}{r}\frac{d}{dr}+\frac{m^2}{r^2},\tilde{\epsilon }(r)=\frac{\epsilon }{2}-v(r),$$

$$g_1(\epsilon ,r)=\frac{v^{\prime }}{\tilde{\epsilon }},g_2(\epsilon ,r)=\frac{2{v^{\prime }}^2}{{\tilde{\epsilon }}^2}+\frac{v^{\prime \prime }}{\tilde{\epsilon }}+\frac{v^{\prime }}{r\tilde{\epsilon }},$$

$${\psi }_{3m,8m}(r)$$

$${\psi }_{4m,7m}(r)$$

$$v(r)=\left\{\begin{array}[]{ll}v_{0},&r<r_{0},\\ 0,&r>r_{0},\end{array}\right.$$

$$\begin{array}[]{rl}\psi_{3m}\bigm{|}_{-}^{+}&=0,\\ \psi_{8m}\bigm{|}_{-}^{+}&=0,\\ 2\frac{d\psi_{3m}}{dr}+ia(\psi_{4m}+\psi_{7m})\biggm{|}_{-}^{+}&=0,\\ \frac{d\psi_{3m}}{dr}+\frac{d\psi_{8m}}{dr}\biggm{|}_{-}^{+}&=0.\end{array}$$

$$\left\{\begin{array}[]{rl}\left[2v(r)\!-\!\varepsilon\!-\!2\!+\!2\hat{k}^{2}\right]\psi_{1}(\mathbf{r})\!-\!a\hat{k}_{+}\psi_{2}(\mathbf{r})\!+\!a\hat{k}_{+}\psi_{5}(\mathbf{r})&\!=0\\ -a\hat{k}_{-}\psi_{1}(\mathbf{r})+\left[2v(r)-\varepsilon\right]\psi_{2}(\mathbf{r})+a\hat{k}_{+}\psi_{6}(\mathbf{r})&\!=0\\ a\hat{k}_{-}\psi_{1}(\mathbf{r})+\left[2v(r)-\varepsilon\right]\psi_{5}(\mathbf{r})-a\hat{k}_{+}\psi_{6}(\mathbf{r})&\!=0\\ a\hat{k}_{-}\psi_{2}(\mathbf{r})\!-\!a\hat{k}_{-}\psi_{5}(\mathbf{r})\!+\!\left[2v(r)\!-\!\varepsilon\!+\!2\!-\!2\hat{k}^{2}\right]\psi_{6}(\mathbf{r})&\!=0\,.\end{array}\right.$$

$${\Psi }_{1,K=0}(\mathbf{r})=\sum _m{\Psi }_{1m}{e}^{im\phi }=\sum _m\left(\begin{array}{c}{\psi }_{1m}(r)e^{i\phi }\\ {\psi }_{2m}(r)\\ {\psi }_{5m}(r)\\ {\psi }_{6m}(r)e^{-i\phi }\end{array}\right)e^{im\phi }.$$

$$\begin{array}[]{rl}\psi_{1m}\bigm{|}_{-}^{+}&=0,\\ \psi_{6m}\bigm{|}_{-}^{+}&=0,\\ \frac{d\psi_{1m}}{dr}-ia\psi_{2m}\biggm{|}_{-}^{+}&=0,\\ \frac{d\psi_{6m}}{dr}-ia\psi_{2m}\biggm{|}_{-}^{+}&=0.\end{array}$$

$$\epsilon =$$

$$\begin{gathered} A^2{\hat{k}}^2{\Phi }_m-A^2\frac{2V^{\prime }(r)}{E-2V(r)}\frac{d{\Phi }_m}{dr}+[EV(r)-V^2(r)]{\Phi }_m \\ +\left(M^2-\frac{E^2}{4}\right){\Phi }_m=0. \end{gathered}$$

$$-\frac{A^2}{2|M|}{\hat{k}}^2{\Phi }_m-\frac{A^2}{2M^2}{V}^{\prime }(r){\Phi }_m+V(r){\Phi }_m=\Delta E{\Phi }_m.$$

$${\hat{k}}^2{\Phi }_m+\frac{1}{r-r_0}\frac{d{\Phi }_m}{dr}+\frac{M^2}{A^2}{\Phi }_m=0.$$

$$\Phi (r)\simeq \mathrm{const}|r-r_0|K_1\left(\left|\frac{M}{A}(r-r_0)\right|\right),$$

$$\begin{array}[]{ll}\psi_{3m}(r)=A^{\pm}_{m}\mathcal{F}_{m}(Q_{\pm}r),&\psi_{4m}(r)=B^{\pm}_{m}\mathcal{F}_{m+1}(Q_{\pm}r),\\ \psi_{7m}(r)=C^{\pm}_{m}\mathcal{F}_{m-1}(Q_{\pm}jr),&\psi_{8m}(r)=D^{\pm}_{m}\mathcal{F}_{m}(Q_{\pm}r),\end{array}$$

Full text

I Introduction

The electron-electron interactions in topological insulators (TIs) currently are one of the most challenging problems in which one can expect the emergence of new and nontrivial properties of electronic systems. The role of the electron-electron interactions in TIs is still poorly understood, but it is already clear that a lot of highly unusual effects arises due to the interplay between Coulomb interactions and topological aspects of matter [1]. Recent experiments clearly demonstrate that in many cases the inclusion of the electron-electron interaction is crucially important for understanding the electron transport in TIs (see, e.g., Refs. [2, 3]). The theoretical researches are focused mainly on the influence of the electron-electron interactions on the topological phase transitions and on the possibility of a topological phase to be formed due to the interaction (for a review of recent studies in this field on the two-dimensional (2D) TIs see, e.g., Ref. [4]).

There is another aspect of the many-body problem associated with the formation of stable or metastable complexes of two or more electrons, such as Cooper pairs, excitons, exciton complexes, etc. Investigations in this direction are carried out very intensively for conventional materials for many years and currently continue to attract great interest [5]. However, such states are still insufficiently studied for topologically nontrivial electronic systems, although it is clear that new properties of the bound states can appear thanks to the unusual band structure of TIs.

Recent studies were focused on electron-hole bound states. It was found that chiral excitons arise on a surface of TIs with a magnetically induced gap in the surface state spectrum [6]. They differ from conventional excitons by a chiral structure and a modified spectrum. Chiral excitons give rise to resonant manifestations in Faraday and Kerr effects [7]. Similarly, edge excitons were found to be formed in 2D TIs in the presence of an in-plain magnetic field [8]. In the bulk of the crystal, the exciton states are affected by the geometrical properties of the Bloch bands. The Berry curvature and quantum geometric tensor essentially modify the exciton states and their spectrum [9, 10]. The studies of excitons in topologically trivial narrow-gap materials revealed a substantial dependence of the exciton properties on the electron dispersion in the bands. Such investigations were carried out in recent years for the quasirelativistic dispersion of electrons and holes in a gapped graphene and carbon nanotubes [11, 12, 13, 14].

In the present paper, we address to the problem of two interacting electrons. In this case, the formation of bond states does not seem obvious because of the Coulomb repulsion. Nevertheless, bound states can be formed because of peculiarities of the band structure. This conclusion can be drawn from recent studies of two-electron states in graphene. Sabio, Sols, and Guinea [15] investigated the problem of two interacting Dirac fermions and revealed a singular behavior of the two-particle wave function at a definite distance between the particles, which evidenced a partial localization of the interacting electrons. The existence of quasilocalized two-particle states in this case was demonstrated in Ref. [16]. Further studies showed that the formation of stationary two-particle states with localized wave function becomes possible when the quasiparticle dispersion substantially deviates from the linear one. First, Mahmoodian and Entin found that the trigonal warping of the spectrum results in the formation of the excitonlike states in some regions of the momentum space [17]. Then, Marnham and Shytov introduced quadratic momentum terms into the single-particle kinetic energy and came to the conclusion that the bound, Cooper-pairlike, states could appear in double-layered structures, but they were metastable [18].

The situation of TIs is obviously more interesting, because in this case there is a gap in the single-particle spectrum and the single-particle states have a more complex orbital structure. The electronic states are a superposition of the states of the electron and hole bands, which are characterized by effective masses of opposite sign. One can therefore expect a nontrivial dynamics of the particles under the action of the Coulomb forces since the relative motion of the particles is determined by a reduced effective mass, the sign of which is not obvious in advance, i.e., without knowledge of the orbital composition of the two-particle wave function which in its turn is determined by the solution of corresponding Schrödinger equation.

This conjecture is supported by the results of recent studies of the bound states localized at impurities with a short-range potential [19, 20, 21]. It turns out that the potential of any sign produces bound states of two kinds in the energy gap of the 2D TIs, in contrast to the topologically trivial case where only one bound state exists. For example, in the case of an impurity with negative potential, one state is formed as a result of the attraction of the electronlike quasiparticle. The captured particle is localized in the center. Other state, on the contrary, arises as a result of the repulsion of the holelike quasiparticle. In this state, the particle is localized around the impurity similarly to edge states. In other words, an impurity produces a bound state in both cases: when the impurity attracts a particle or repels it.

In this connection a natural question arises whether two electrons form a bound state in 2D TIs when the Coulomb force acts between them? To answer this question in the present work we study two-particle states within the model proposed by Bernevig, Hughes, and Zhang [22] (BHZ). The model is widely used for 2D TIs, but it describes also a trivial phase under appropriately chosen parameters, so that we can compare the results obtained in both cases to reveal effects that arise only in a topological phase. Two-electron states have not yet been studied in the BHZ model. We solve this problem and show that bound states indeed arise. It is found that the bound states can be formed in both the topological and trivial phases at any sign of the pair interaction potential. However, in the topological phase, new bound states appear in addition to those in the trivial phase. They arise at lower interaction potential and have a higher bounding energy. We study general properties of the bound state spectra and classify the states according to their spin structure as singletlike and tripletlike ones. With respect to orbital degrees of freedom, the bound states are well classified only in the case of small interaction potential where the states are separated into two groups. It the one group, the states are mainly formed by basis states in which both electrons are in the same (electron or hole) band. It the other group, the pairing electrons are in the different bands.

The paper is organized as follows. Section II presents general equations. Here, we also classify the two-particle states and simplify the problem by addressing to the case of zero center-of-mass momentum and to a model potential. In Sec. III, we present qualitative arguments explaining the bound state formation. Section IV is devoted to singletlike bound states. The tripletlike states are considered in Sec. V. In Sec. VI, a specific case is studied to show that new bound states arise due to the band inversion in the topological phase. In Sec. VII the topologically trivial case is considered and the band-inversion effect is discussed. Finally, in Sec. VIII, we summarize the results.

II General equations

We start with a statement of the two-body problem in the BHZ model. The BHZ model presents single-particle electronic states in the frame of the $\mathbf{kp}$ theory with using four-band basis $(|E\uparrow\rangle,|H\uparrow\rangle,|E\downarrow\rangle,|H\downarrow\rangle)^{T}$, where $|E\uparrow\rangle$ and $|E\downarrow\rangle$ are a superposition of the electron- and light-hole states with the moment projection $m_{J}=\pm 1/2$; $|H\uparrow\rangle$ and $|H\downarrow\rangle$ are the heavy-hole states with $m_{J}=\pm 3/2$. The single-particle Hamiltonian that determines the spinor of the envelope functions reads

[TABLE]

[TABLE]

where $\mathbf{\hat{k}}$ is the quasimomentum operator, $A$, $B$, and $M$ are the parameters of the BHZ model. Here, for simplicity, we do not take into account the spin-orbit interaction, which can actually be present due to structural inversion asymmetry and bulk inversion asymmetry. The terms describing the asymmetry of the electron and hole bands also are dropped for simplicity. These assumptions do not have a decisive impact on the results but greatly simplify the calculations. The BHZ model describes both topological and trivial phases of a 2D electron system in a crystal. Trivial phase is realized at the ordinary arrangement of the electron and hole bands, when $M/B<0$. In the topological phase the band structure is inverted, $M/B>0$.

Two-particle wave functions are represented by a spinor of 16th order, $\Psi(\mathbf{r}_{1},\mathbf{r}_{2})=\left(\psi_{1},\psi_{2},\psi_{3},...,\psi_{16}\right)^{T}$, which defines the envelope functions in the basis:

[TABLE]

The Hamiltonian of two interacting electrons has the form

[TABLE]

where $V(\mathbf{r})$ is the pair interaction potential which is supposed to be a given function.

The wave function $\Psi(\mathbf{r}_{1},\mathbf{r}_{2})$ is determined by the Schrödinger equation

[TABLE]

Due to the block-diagonal structure of the single-particle Hamiltonian (1), the Schrödinger equation (5) splits into four uncoupled equations for the following wave functions:

[TABLE]

Here, for clarity, we have written both the envelope functions $\psi_{1},\psi_{2},\dots$ and the corresponding basis functions. The states described by $\Psi_{1}(1,2)$ and $\Psi_{4}(1,2)$ are composed of the spin-up and spin-down orbitals. Therefore, they can be conventionally classified as tripletlike states. Similarly, the states $\Psi_{2}(1,2)$ and $\Psi_{3}(1,2)$ can be called singletlike ones. These terms are not strict here. In Sec. IV it will be shown that the wave functions $\Psi_{2}(1,2)$ and $\Psi_{3}(1,2)$ describe the same bound state.

The wave functions $\Psi_{j}(1,2)$ are determined by equations of the following form:

[TABLE]

Here and in what follows we use dimensionless notations:

[TABLE]

For convenience, the prime in the variable $\mathbf{r^{\prime}}$ will be omitted. To separate the topological and trivial phases we introduce a parameter $\lambda=M/|M|=\pm 1$ and assume that $B<0$. In this case, $\lambda=1$ corresponds to the trivial phase and $\lambda=-1$ corresponds to the topological phase.

The operators $\hat{H}_{j}$ are $4\times 4$ matrices, the elements of which are expressed via the operators $\hat{\mathbf{k}}_{1}$ and $\hat{\mathbf{k}}_{2}$. Before we present the equations of motion in an explicit form, it is meaningful to modify the wave functions taking into account that the system is translationally invariant.

In order to study the bound states it would be natural to try to separate the relative motion of the particles from their movement as a whole. Therefore we switch to the center-of-mass frame, defining the new coordinates: $\mathbf{R}=(\mathbf{r}_{1}+\mathbf{r}_{2})/2$ and $\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$. However, within the BHZ model the relative motion and the motion of the center of mass are not separated because the nondiagonal terms in Eq. (2), which determine the coupling of the electron and hole bands, depend on the momenta of each particle. Nevertheless, since the system is translationally invariant, the wave function can be represented in the form:

[TABLE]

where $\mathbf{K}$ is the total momentum of the pair.

Of most interest are the functions $\Psi_{j,\mathbf{K}}(\mathbf{r})$ that determine the relative motion. The fact that $\Psi_{j,\mathbf{K}}(\mathbf{r})$ is explicitly dependent on $\mathbf{K}$ means that the spectrum of the bound states and their structure depend on the total momentum.

The equations defining the functions $\Psi_{j,\mathbf{K}}(\mathbf{r})$ are different for all four states given by Eq. 6. In the case of the tripletlike state with the moments up, the function $\Psi_{1,\mathbf{K}}(\mathbf{r})=\left(\psi_{1},\psi_{2},\psi_{5},\psi_{6}\right)^{T}$ is determined by the following equation system:

[TABLE]

The wave function $\Psi_{2,\mathbf{K}}(\mathbf{r})=\left(\psi_{3},\psi_{4},\psi_{7},\psi_{8}\right)^{T}$ of the singletlike state is defined by the following equations:

[TABLE]

The wave functions $\Psi_{4,\mathbf{K}}(\mathbf{r})$ and $\Psi_{3,\mathbf{K}}(\mathbf{r})$ are described by similar equations which are not presented here for the sake of brevity.

In Eqs. (10) and (11), $\hat{\mathbf{k}}=(\hat{\mathbf{k}}_{1}-\hat{\mathbf{k}}_{2})/2$ is the operator of the relative momentum, $\hat{k}_{\pm}=\hat{k}_{x}\pm i\hat{k}_{y}$ and $K_{\pm}=K_{x}\pm iK_{y}$.

II.1 The case of zero total momentum

In order to understand the essential features of the problem, of most interest is the case of zero total center-of-mass momentum. One can expect that in this case the pairing effect is particularly important since the kinetic energy of the pair is minimal. This is also the simplest case for calculations since Eqs (10) and (11) are considerably simplified. From these equations, one can suppose that the finite value of the total momentum does not strongly change the solutions as long as $K$ is small enough. This is why we will focus on the case where $\mathbf{K}=0$.

To be specific, consider first the state described by the wave function $\Psi_{2,\mathbf{K}}(\mathbf{r})$. In the case of $\mathbf{K}=0$, Eq. (11) is simplified as follows:

[TABLE]

For further analysis, it is convenient to go to polar coordinates ($r,\varphi$) and expand the wave functions in the Fourier series:

[TABLE]

In this way, the system of Eqs (12) is reduced to independent systems of four equations defining the Fourier components $\psi_{3m},\psi_{4m},\psi_{7m}$, and $\psi_{8m}$ for each $m$. We do not write them explicitly, so as not to clutter the paper.

Generally speaking, it is possible to transform the problem further by reducing it to an equation for a single function. This turns out to be useful for the further analysis. If we introduce the function $\Phi_{m}(r)=\psi_{3m}(r)+\psi_{8m}(r)$, one can exclude all functions $\psi_{3m},\psi_{4m},\psi_{7m},\psi_{8m}$ and obtain a single equation for $\Phi_{m}(r)$,

[TABLE]

where the following designations are used

[TABLE]

[TABLE]

$v^{\prime}$ and $v^{\prime\prime}$ are the potential derivatives.

The functions $\psi_{3m}$, $\psi_{4m}$, $\psi_{7m}$, and $\psi_{8m}$ are expressed through $\Phi_{m}(r)$ as follows:

[TABLE]

Equation (14) can be quite simply analyzed. It is seen that the equation has a singular point in which $\varepsilon-2v(r_{c})=0$. If the interaction is repulsive, $v(r)>0$, the singularity exists only for $\varepsilon>0$. In this case, the solution can be analyzed by the expansion of the function $\Phi_{m}(r)$ near the singular point: $\Phi_{m}(r)=|r-r_{c}|^{\lambda}\sum_{l}a_{l}(r-r_{c})^{l}$. In this way we come to the conclusion that $\Phi_{m}(r)$ and the components of the spinor $\Psi_{2m}(r)$ do not diverge at the point $r=r_{c}$. This fact allows one to further simplify the problem by using a model potential.

II.2 Model steplike potential

The physical understanding of the structure of two-particle states and their spectrum can be obtained by considering a model potential $v(r)$, which has the basic properties of the real potential of the pair interaction. As a model potential we choose a steplike function,

[TABLE]

which is widely used and usually gives a good effective description of a more general class of short-range potentials.

When using the steplike potential, an important point is to obtain matching conditions for the wave functions at the radius $r=r_{0}$. They should be obtained by integrating the full equations defining $\Psi_{j,\mathbf{K}}$ over the transition region, $|r-r_{0}|<\delta$, assuming that $v(r)$ is a finite value. Finally, the limit $\delta\to 0$ should be taken.

In this way, we arrive at the following matching equations in the case $\mathbf{K}=0$. For the singletlike states, one obtains

[TABLE]

It is interesting to note that the function $\Phi_{m}(r)$ is continuous at $r=r_{0}$.

The same approach can be used for the tripletlike states. To be specific, we consider the state $\Psi_{1,\mathbf{K}}(\mathbf{r})$ at $\mathbf{K}=0$. The components $\psi_{1}(\mathbf{r})$, $\psi_{2}(\mathbf{r})$, $\psi_{5}(\mathbf{r})$, and $\psi_{6}(\mathbf{r})$ of the envelope function spinor are defined by the following equations:

[TABLE]

It seen that the components $\psi_{5}(\mathbf{r})$ and $\psi_{2}(\mathbf{r})$ are connected by a simple relation: $\psi_{2}(\mathbf{r})=-\psi_{5}(\mathbf{r})$.

In the polar coordinates, $\Psi_{1,K=0}$ is presented in the form of the Fourier series:

[TABLE]

Equations defining the components $\psi_{1m}(r)$, $\psi_{2m}(r)$, $\psi_{5m}(r)$ and $\psi_{6m}(r)$ are easily obtained from Eq. (21). The matching equations have the form:

[TABLE]

Thus, Eqs. (13), (12), and (20) fully define the singletlike state $\Psi_{2,K=0}(\mathbf{r})$. Correspondingly, Eqs. (21), (22), and (23) define the tripletlike state $\Psi_{1,K=0}(\mathbf{r})$. Equations defining $\Psi_{3,K=0}(\mathbf{r})$ and $\Psi_{4,K=0}(\mathbf{r})$ can be obtained in a similar way. These equation are straightforwardly solved in Secs. IV–VII, but before presenting the results of the calculations in detail, it is reasonable to stay on a qualitative picture of the bound state formation based on simplified models.

III A qualitative picture

In this section, we provide physical arguments, which qualitatively explain the mechanism of the bound state formation. They allow one also to better understand the main types of the bound states, which are obtained by solving the equations presented in the previous section. These arguments are derived from simplified models with using different additional assumptions.

III.1 Step potential

First, consider the case of the interaction potential of the step form. Let us divide the space of the relative coordinate $\mathbf{r}$ into two regions: the interaction region, $r<r_{0}$, where $v(r)=v_{0}$, and the outer region, $r>r_{0}$, where the interaction is absent.

In the outer region, the particles move freely, so that the spectrum of two particles with the zero total momentum contains two bands and zero-energy level:

[TABLE]

The bands correspond to the particle configuration in which both particles have the energy in the conduction band or in the valence band. The zero-energy level is infinitely degenerate. This is a zero-energy mode, which corresponds to the case where the particles are in different bands with opposed momenta.

In the interaction region, the situation is very similar. Equation (5) shows that the only effect of the inter-particle interaction on the two-particle spectrum is the shift of the energy by $2v_{0}$:

[TABLE]

Now imagine the energy diagram of the two-particle system in the space of the relative coordinate, see Fig. 1. Here, the areas in which there are propagating solutions of the Schrödinger equation, are colored in light green (darkened). Uncolored areas are classically inaccessible for the particles. Inside them, the solutions decay. It is obvious that in the energy interval $-2<\varepsilon<2(-1+v_{0})$, the propagating solutions exist in the interaction region. In the outer region, the wave function decays if $\varepsilon<2$. The propagating solutions can interfere within the interaction region to form a two-particle bound state.

Another possibility to realize propagating solutions in the interaction region appears near the energy level of the zero-energy mode $\varepsilon\approx 2v_{0}$ in the interaction region, $r<r_{0}$. At this energy level, the wave functions decay outside the interaction region when $v_{0}<1$. In this way, localized two-particle states can also arise.

One can say that an effective quantum dot is formed in the interaction region where the bound states can be formed in the energy interval $-2<\varepsilon<2(-1+v_{0})$ and near the energy $\varepsilon\approx 2v_{0}$.

Of course, these arguments are very qualitative, but in many respects they are true. The qualitative results agree with the rigorous solutions given in the following sections.

III.2 A simplified model

Another qualitative approach that illustrates the possibility of the bound state formation is based on the simplified model, which was used in early works on electron systems with inverted band structure [23, 24]. In terms of the BHZ model, the simplified model neglects the term $Bk^{2}$ in the diagonal elements of the $\hat{h(k)}$ matrix in the Hamiltonian $\hat{H}_{0}$, see Eqs (1) and (2). This simplification is usually justified by the fact that the calculations are carried out within the $\mathbf{kp}$ theory where $k$ is supposed to be small. The system, reduced in such a way, in many cases leads to qualitatively correct results, but something is lost in it. The idea to qualitatively analyze localized states within two-band model by reducing it to an effectively single-band model goes back to the old work of L.V. Keldysh [25], where this idea is applied to single-particle states localized at an impurity.

This simplification reduces the problem to the two-particle Dirac equation, which is well known in the quantum field theory. In recent years, the two-particle Dirac equation with a simplified electron-interaction interaction was adapted to the narrow-gap and gapless electronic systems in graphene and carbon nanotubes [15, 12, 16, 13, 14].

In this way, we arrive at the following results. To be specific we consider only the singletlike states and turn to Eq. (14) for the function $\Phi_{m}(r)$. Turning to the reduced model, we have to put $B=0$ in the BHZ Hamiltonian. Since $B$ is used in defining the dimensionless variables, introduced in Eq. (8), one needs to go back to the dimensional quantities, which will be used only in Eqs. (28)–(31) below. Finally, in the reduced model we arrive at the following equation instead of Eq. (14):

[TABLE]

Note that a similar equation was used for the two-particle systems in Refs. [15, 12, 16, 13, 14].

One can consider this equation as a single-particle problem, where the first term is the kinetic energy and the third term plays a role of an effective potential. Its sign depends on the energy and the magnitude of $V(r)$. So, the effective potential is negative (i.e. attractive) when $E<0$. It can be negative also at $E>0$, if the real potential $V$ is high enough. Since the effective potential is attractive in a wide range of $E$, one can expect that bound states can be formed, at least in the case where the second term in Eq. (28) is not large [for example, when $V(r)$ is a slow varying function].

It is interesting to consider the same equation from other point of view. We redefine the potential so that it becomes positive (for example, by multiplying the equation by a number). In this case, the first term, which plays the role of kinetic energy, becomes negative and therefore the effective reduced mass of the two particles is negative.

This is well illustrated by considering Eq. (28) in the limiting case where $V(r)\ll|M|$, $E=-2|M|+2\Delta E$ and $\Delta E\ll|M|$. Equation (28) takes the form

[TABLE]

It is seen that the effective reduced mass is negative and can be defined as $m^{*}=-\hbar^{2}|M|/A^{2}$. Thus, a bound state is formed by the positive potential. This is obvious, at least, if one neglects the second term, which really can be dropped if the characteristic length $l$ of the potential change is large, $l^{2}\gg A^{2}/(2|M|)$.

Another possibility for a bound state to appear arises because of the singularity of the second term in Eq. (28) in the point $r=r_{0}$ where $E=2V(r_{0})$. Again, consider a simplified case where $E/2,V(r)\ll|M|$ and $E/2$ is smaller than the maximum value of the interaction potential. If $V(r)$ is a monotonic function, there is one singular point. We focus on the solution of Eq. (28) near the point $r=r_{0}$ by expanding the potential: $V(r)=E/2+V^{\prime}(r_{0})(r-r_{0})+\dots$. In this case, Eq (28) takes the form

[TABLE]

This equation is easily solved in the vicinity of the point $r=r_{0}$. For $m=0$, one obtains the following solution

[TABLE]

where $K_{1}(z)$ is the modified Bessel function of the second kind.

This fact argues that there can be a solution localized near the point $r=r_{0}$, but it fails to determine the eigenenergy $E$ since the complete solution satisfying boundary conditions is not found. Nevertheless, the qualitative behavior of the wave function agrees with the results of Refs. [15, 16] where such a singularity was studied in the case of massless Dirac fermions and it was found that a quasibound state appears with the wave function effectively localized near $r=r_{0}$. The estimate (31) qualitatively agrees also with the total solution of the problem, which will be presented in the next sections.

IV Singletlike bound states

In this section, the singletlike states are studied by the direct solution of Eq. (11) in the case where the potential has the step form and $\mathbf{K}=0$ . First, we consider of the states $\Psi_{2,\mathbf{K=0}}(\mathbf{r})$. We find their spectrum and the spatial distribution of all components of the envelope function spinor. Then the results are generalized to the states $\Psi_{3,\mathbf{K=0}}(\mathbf{r})$, and finally we obtain the two-particle wave function which is antisymmetric with respect to the permutation of the particles. To be specific we consider below the case of the topological insulator ($\lambda=-1$). The topologically trivial case will be presented in Sec. VII.

IV.1 Spectrum

The states described by the wave function $\Psi_{2,\mathbf{K}=0}(\mathbf{r})$ (the index $\mathbf{K}=0$ is dropped hereinafter) are determined by Eqs. (12). We solve these equations in the regions $r<r_{0}$ and $r>r_{0}$, and match the found functions at $r=r_{0}$ with using Eqs. (20).

In the case of the step potential, Eqs. (12) are easily solved in terms of the Bessel functions. The fundamental set of solutions for the components of the spinor $\Psi_{2}$ has the form:

[TABLE]

where the wave numbers $Q_{\pm}$ are the roots of the dispersion equation, which has a unified form in both regions:

[TABLE]

where $\widetilde{\varepsilon}$ takes different values for the interaction region and the outer region,

[TABLE]

For convenience we have denoted here $\varepsilon_{0}\equiv\varepsilon/2$, which is the energy of an electron pair per particle.

The explicit expression for $Q_{\pm}$ reads

[TABLE]

In what follows it is important that $Q_{\pm}$ can be real, imaginary or complex, depending on the parameter $a$ and the energy $\widetilde{\varepsilon}$. The map of possible values of $Q_{\pm}$ on the plane $(a^{2},\widetilde{\varepsilon})$ is shown in Fig. 2.

In Eq. (32), $\mathcal{F}_{m}(Q_{\pm}r)$ is a fundamental solution of the Bessel equation. $\mathcal{F}_{m}(Q_{\pm}r)$ can be written as any pair of the Bessel functions: $J_{m}(Q_{\pm}r)$ and $Y_{m}(Q_{\pm}r)$; $I_{m}(Q_{\pm}r)$ and $K_{m}(Q_{\pm}r)$; $H_{m}^{(1),(2)}(Q_{\pm}r)$; etc.. The choice of the pair of Bessel functions in a specific case is determined by the values of $Q_{\pm}$ at given $a$ and $\widetilde{\varepsilon}$ (in accordance with the map in Fig. 2), and by the behavior of the Bessel function at $r\to 0$ and $r\to\infty$.

To be specific consider the case where $a^{2}>4$.

In the energy interval $-1<\varepsilon_{0}<-1+v_{0}$, the solution of Eqs (12) can be presented in the following form:

(i) at $r<r_{0}$,

[TABLE]

where

[TABLE]

and

[TABLE]

(ii) at $r>r_{0}$,

[TABLE]

where

[TABLE]

and

[TABLE]

Now the functions (36) and (39) should be matched at the boundary $r=r_{0}$. Using Eqs. (20), we get a homogeneous equation system for the coefficients $A_{+},A_{-},B_{+},B_{-}$. The equations are very cumbersome, so we do not give them and subsequent equations in an explicit form. The determinant $\mathfrak{D}$ of this equation system is a function of the energy $\varepsilon_{0}$ and the parameters $a,v_{0},r_{0},m$. The eigenenergies are determined by the equation

[TABLE]

It turns out that this equation has several solutions:

[TABLE]

where $n$ is a root number at given parameters $a$, $v_{0}$, $r_{0}$, and the angular number $m$. One can say that $n$ is the radial quantum number. The upper symbol indicates that this is a singletlike state.

In the energy interval $-1+v_{0}<\varepsilon_{0}<1$, the solution of Eqs (12) differs from that considered above since in the interaction region both roots $Q_{\pm}$, see Eq. (35), are imaginary. Therefore the solution in the region $r<r_{0}$ should be composed of the Bessel functions $I_{m}(|Q_{\pm}|)$. This is the only difference from Eqs (36) and (39). Moreover, it is clear that Eqs (36) and (39) are formally correct in the interval $-1+v_{0}<\varepsilon_{0}<1$, if one considers $k_{+}$ as a complex number.

Equation (42) for the eigenvalues of the energy is very cumbersome and complicated. In this paper, we solve it numerically. This approach enables us to find solutions for a finite value of fundamentally important parameters $v_{0}$ and $r_{0}$. As a result we demonstrate the presence of bound states and the main features of their spectrum. The main result of these studies is that two-particle bound states exist in a wide range of the parameters $v_{0}$, $r_{0}$ and $a$. Energy levels of the bound states lie in the gap of the band spectrum.

The analysis has shown that the spectrum of the two-particle states is more complicate than one could expect from the qualitative arguments of Sec. III. In this section we have restrict ourselves to the bound states with zero angular number and the region of the parameter $|a|>2$. In this case, the calculations turn out to be more simple and the results seem to be quite general. Qualitatively new features of the bound states are expected when the parameter $|a|<2$. This case will be studied in Sec. VI.

There are two groups of the bound states in accordance with the qualitative arguments of Sec. III. They are classified by the energy that a bound state has at low interaction potential. The energy levels of the first group appear at the bottom of the gap of the two-particle band spectrum and then rise with increasing $v_{0}$. The bound states of the second group have the energy near the center of the band gap at low $v_{0}$. The behavior of the energy levels of the first and second groups with the increase of the interaction potential is illustrated in Fig. 3(a) in the case where the interaction radius $r_{0}$ is not large, so that only one state of the first group exists, when $m=0$. In contract, the second group contains two states with $m=0$. Of course, the states with $m\neq 0$ also exist in both groups, but the dependence of their energy on the parameters $v_{0}$ and $r_{0}$ is more complicated than one might expect at first glance.

It is obvious that this classification is justified only at low interaction potential. When $v$ is comparable with the band gap, this classification is very conventional and little constructive. Nevertheless, we will stick to it to trace the evolution of the bound states with increasing the interaction potential.

The dependence of the bound state energy on the interaction radius $r_{0}$ is illustrated in Fig. 3(b). The states of the first group behave as follows. With increasing $r_{0}$, new roots of the determinant $\mathfrak{D}(\varepsilon;a,v_{0},r_{0},m)$ successively appear at the bottom of the gap. This remembers the usual picture of quantization in a quantum dot. In the case we are studying, such a quantum dot is effectively formed by the interaction potential as illustrated in Fig. 1. Of course, the quantization conditions are very different from those in ordinary quantum dots in one-band model with a quadratic dispersion.

The states of the second group show a completely different behavior. With the increase of $r_{0}$, no additional roots appear with zero angular number. This feature could be understood as a result of the fact that the quantum state is localized along the perimeter of the effective quantum dot similarly to an edge state, rather than inside it. In this case, only the angular motion is quantized.

Thus, in order to elucidate the mechanism of the bound state formation it is interesting to analyze the spatial distribution of the electron density and the density of all components of the envelope function spinor.

IV.2 Electronic structure of the bound states

The envelope functions $\psi_{3}(r)$, $\psi_{4}(r)$, $\psi_{7}(r)$ and $\psi_{8}(r)$ in the state $\Psi_{2}(\mathbf{r})$ can be calculated straightforwardly with using Eqs (36), (39) and coefficients $A_{+}$, $A_{-}$, $B_{+}$, and $B_{-}$.

We begin with the states of the first group.

IV.2.1 First group of bound states

The radial distribution of the density of all spinor components [$\psi_{3}(r)$, $\psi_{4}(r)$, $\psi_{7}(r)$, and $\psi_{8}(r)$] is shown in Figs. 4(a)–(c) for the bound state corresponding to the line 1.1 in Fig. 3(a). Here it should be noted that in the case of zero angular number, $m=0$, the spinor components $\psi_{4}(r)$ and $\psi_{7}(r)$ coincide though, in general, $\psi_{4}(r)\neq\psi_{7}(r)$. The discontinuity of some envelope functions or their derivatives at $r=r_{0}$ originates from the singularity of the potential. A separate investigation of the solutions in the vicinity of the point $r=r_{0}$ in the case where the potential is approximated by a linear function with large gradient, shows that the wave function is also continuous but sharply changes.

It is seen that the spinor components $\psi_{3}(r)$ and $\psi_{8}(r)$ have the largest amplitude. They describe the contribution of the two-particle basis states $|E\uparrow E\downarrow\rangle$ and $|H\uparrow H\downarrow\rangle$ into the total wave function, respectively. Thus, the states of this group are formed mainly by those orbital components, in which both particles are in the electron band or in the hole band. The contribution of the mixed components $|E\uparrow H\downarrow\rangle$ and $|H\uparrow E\downarrow\rangle$ is small. In addition, the amplitude of the mixed components strongly decreases with decreasing $v_{0}$. Another conclusion is that the electron density is distributed in the volume of the effective quantum dot, though there is also a small density located at the edge.

The radial distribution of the total density $|\Psi_{2}(r)|^{2}=|\psi_{3}(r)|^{2}+|\psi_{4}(r)|^{2}+|\psi_{7}(r)|^{2}+|\psi_{8}(r)|^{2}$ is shown in Fig. 4(d).

Now we turn to the second group of the bound states.

IV.2.2 Second group of bound states

The radial distribution of the densities of the spinor components is shown in Fig. 5(a,b,c) for the bound state shown by the line 2.1 in Fig. 3.

In this group of states, the amplitude of the components $\psi_{4}(r)$ and $\psi_{7}(r)$ noticeably increases as compared with the first group states. These components represent the contribution of the mixed states of the electron and hole bands ($|E\uparrow H\downarrow\rangle$ and $|H\uparrow E\downarrow\rangle$) to the total wave function. However, the main distinction from the first group states is that the particle density is concentrated at the edge of the effective quantum dot. Hence, this state can be considered as a kind of edge states.

IV.3 Two-electron wave function

The true wave function of the two electrons is to be antisymmetric with respect to the permutation of particles. The wave function described by the spinor $\Psi_{2,\mathbf{K}}(\mathbf{r})$ does not satisfy this requirement. Therefore the two-electron wave function should be presented in the form:

[TABLE]

where the arguments $1,2$ denote the coordinates of two electrons and $\Psi_{2}(1,2)$ is the wave function of the state described by the envelope function spinor $\Psi_{2}(\mathbf{r})$. The permutation of the particles includes both the replacement $\mathbf{r}\to-\mathbf{r}$ and the interchange of the particle coordinates in the two-particle basis functions, Eq. (3). Taking into account this fact, we arrive at the following wave function:

[TABLE]

As it is seen, the wave function (45) can not be factorized into orbital and spin functions. Therefore, this state can not be called a singlet state in the usual sense. Nevertheless we continue to use this nonstrict term.

We finish this section by considering another singletlike state $\Psi_{3}(\mathbf{r})$. Straightforward calculations show that the components of this spinor [$\psi_{9}(\mathbf{r})$, $\psi_{10}(\mathbf{r})$, $\psi_{13}(\mathbf{r})$, and $\psi_{14}(\mathbf{r})$] are determined by the same equations as the components of $\Psi_{2}(\mathbf{r})$. One can show that the components of $\Psi_{3}(\mathbf{r})$ are connected with those of $\Psi_{2}(\mathbf{r})$ by the following replacement: $\psi_{9}\to\psi_{3}$, $\psi_{10}\to-\psi_{7}$, $\psi_{13}\to-\psi_{4}$ and $\psi_{14}\to\psi_{8}$. Taking into account this replacement together with the replacements in the basis functions, it is easy to see that the spinor wave function $\Psi_{3}(\mathbf{r})$ differs from the wave function $\Psi_{2}(\mathbf{r})$ simply by the permutation of the particles. Thus, the quantum state described by the spinor $\Psi_{3}(\mathbf{r})$ coincides with the already-studied state, $\Psi_{m}^{(s)}(1,2)$.

V Tripletlike bound states

The tripletlike states are studied similarly to the singletlike ones with using the same simplifications. Therefore, we do not go into the details and only present main results.

There are two tripletlike states: $\Psi_{1,\mathbf{K}}(\mathbf{r})$ and $\Psi_{4,\mathbf{K}}(\mathbf{r})$ which differ only in the direction of the spins. Since the system studied here has $S_{z}$ symmetry, the other properties of these states are the same. Below we consider only the state $\Psi_{1,\mathbf{K}}(\mathbf{r})$ at $\mathbf{K}=0$ and focus on the topologically nontrivial case as in the previous section.

V.0.1 The spectrum

The energy spectrum of the tripletlike bound states is generally similar to the spectrum of the singletlike states, but there are some differences in details. The tripletlike bound states can also be divided into two groups which differ in the energy at the low interaction potential. The states of the first group have the energy at the bottom of the band gap, while the energy of the second-group states lies near the center of the gap. The bound-state energy is determined by three parameters of the model ($a$, $v_{0}$ and $r_{0}$) and two quantum numbers: the radial quantum number $n$ and the angular quantum number $m$. The dependence of the energy on the potential amplitude is illustrated in Fig. 6(a) for $m=0$. The interaction radius is chosen so small that there is only one energy level of the states of the first group. In contrast, the second group contains two states even if $v_{0}$ is small. The evolution of the spectrum with increasing $r_{0}$ is shown in Fig. 6(b). It is seen that new states with $m=0$ arise only in the first group. They are characterized by the radial quantum number.

V.0.2 Electronic structure of the bound states

Electronic structure of the tripletlike bound states in many respects is also similar to that of the singletlike state, but there are many significant differences in the spatial distribution of the densities of the spinor components related to the electron and hole bands.

The radial distribution of the spinor-component densities is shown in Figs 7–9. To compare the results with those for the singletlike states, the parameters $a$, $v_{0}$ and $r_{0}$ are chosen the same as in Figs 4 and 5.

In the case of the first group states, Fig. 7, the main feature is that the components representing the configuration where the particles are in the different bands, such as $|E\uparrow H\uparrow\rangle$, strongly increase in the tripletlike states in comparison with the corresponding singletlike state. Another peculiarity is that the density of the components corresponding the configuration in which both particles are in the same band, such as $|E\uparrow E\uparrow\rangle$, turns to zero in the center.

The second group of the bound states contains two states with $m=0$, both states being present at $v_{0}=2$. First, consider the states with lower energy (see the line 2.1 in Fig. 6). The spatial distribution of the spinor-component densities in this state is shown in Fig. 8. Of largest value are the components in which the particles are in the different bands. Their density is distributed mainly in the bulk of the effective quantum dot, as it is seen in Fig. 8(b). This density distribution strongly differs from that in the case of the singletlike states, where the density in the low-energy branch of the second group is located near the edge of the effective quantum dot, Fig. 5.

The higher-energy states of the second group (see the line 2.2 in Fig. 6) are in contrast located at the edge of the effective quantum dot, Fig. 9(a-d). One should note that in the case of the singletlike states, the higher-energy states of the second group have strongly different distribution of the spinor-component densities. When $v_{0}$ is not small, the predominant components are those in which the particles are in the same bands.

The tripletlike wave functions are antisymmetrized in the same way as described above, so that the antisymmetric wave function reads

[TABLE]

VI Bound states in topological phase with nearly flat bands

We turn to the question of how the coupling of the electron and hole bands affect the bound states in the topological phase where the bands are inverted. The band coupling is characterized by the parameter $a$. When $|a|>2$, the spectrum and the electronic structure of the bound states are little changed qualitatively with varying $a$. However at $|a|<2$ the situation changes in two aspects.

First, the single-particle wave functions with the energy in the gap are changed radically since the wave vector becomes complex. So that the wave functions not only decay with the distance but also oscillate. The two-particle wave functions behave similarly, since their wave vectors $Q_{\pm}$, see Eq. (35), are complex. Therefore additional oscillating components appear in the fundamental solutions, such as Eq. (32), that form the bound-state wave function. In this case one can expect the appearance of new solutions.

The second aspect is that the single-particle spectrum also changes essentially with $a$. The spectrum shape changes from nearly parabolic one at $|a|\gg 1$ to that of a mexican-hat form at $|a|<\sqrt{2}$. Correspondingly, the effective mass of electrons also changes very strongly and even changes its sign. The effective mass near the band boundaries is known to play an important role. It is usually supposed that the two-electron bound state is formed due to a negative single-particle energy dispersion near the top band boundary [26, 27].

In this section, we consider the two-particle bound states in the case where $a=\sqrt{2}$. This case is very interesting for two reasons. First, at $a=\sqrt{2}$ the real and imaginary parts of the wave vectors of the states with the energy in the gap are of the same magnitude. Therefore evanescent states in the gap are described by the wave functions which have an oscillating component. Because of this, one can expect that nontrivial interference effects appear in the presence of a spatially inhomogeneous potential.

Other reason in that the effective mass goes to $\pm\infty$, respectively, at the bottom of the conduction band and the top of the valence band and does not change the sign with changing the energy in the bands. The single-particle energy dispersion is described by the equation $\varepsilon_{\pm}=\pm\sqrt{1+k^{4}}$.

The calculations are carried out in the same way as in Sec. IV. The two-particle energy in the interaction region and outside of it has the form

[TABLE]

The characteristic wave vectors $Q_{\pm}$ defined by Eq. (35) take the following values

[TABLE]

for $\widetilde{\varepsilon}^{\,2}>1$ and

[TABLE]

for $\widetilde{\varepsilon}^{\,2}<1$.

The fundamental solutions $\mathcal{F}_{m}(Qr)$ of the equation system describing the spinor components in the case of both singletlike and tripletlike states are as follows:

(i) For $\widetilde{\varepsilon}^{\,2}>1$, the functions $\mathcal{F}_{m}(Qr)$ are the Bessel functions of the first and second kinds: $J_{m}(qr)$, $Y_{m}(qr)$, $I_{m}(qr)$, and $K_{m}(qr)$, where $q=\sqrt[4]{\widetilde{\varepsilon}^{\,2}-1}$.

(ii) For $\widetilde{\varepsilon}^{2}<1$, the fundamental solutions are $\mathcal{V}_{m}^{\pm}=\mathrm{ber}_{m}(\gamma r)\pm i\mathrm{bei}_{m}(\gamma r)$ and $\mathcal{W}_{m}^{\pm}=\mathrm{ker}_{m}(\gamma r)\pm i\mathrm{kei}_{m}(\gamma r)$, where $\mathrm{ber}_{m}(z)$, $\mathrm{bei}_{m}(z)$, $\mathrm{ker}_{m}(z)$, and $\mathrm{kei}_{m}(z)$ are the Kelvin functions, and $\gamma=\sqrt[4]{1-\widetilde{\varepsilon}^{\,2}}$.

The calculations lead to the following results for the singletlike states. The spectrum of the bound states significantly changes as compared to the case of $a>2$. The main difference is that new states appear. Fig. 10 presents the bound-state spectrum for the same parameters ($r_{0}$ and $m$) as in Fig. 3. New branches are seen to appear in both groups of states.

In the first group, there are two states with $m=0$ at small $v_{0}$ (the lines 1.1 and 1.2) in contrast to the case of $a>2$ where there is only one state. One of the states has a very small binding energy, while the other state has a large energy and arises at a much smaller interaction potential.

Let us discuss the nature of these states. According to the mechanism of the bound state formation due to the negative single-particle dispersion near the valence-band top [26], one could expect that the energy of the bound state will be close to the amplitude of the interaction potential $v_{0}$, since the effective mass tends to infinity. This really happens with the large-energy state (line 1.2). Obviously, the mechanism of formation of the second state is not directly related to the effective mass in the valence band, since its energy is very close to the bottom of the gap, which indicates an anomalously small effective mass. Therefore we conclude, that the low-energy state (line 1.1) arises because of an interference effect of the evanescent states which results in the appearance of new roots of Eq. (42).

The states 1.1 and 1.2 differ greatly in the spatial distribution of the density of the spinor components and the magnitude of different components. This is illustrated in Figs 11 and 12. It can be seen that, in state 1.1, the electron density is mainly localized around the interaction region, which indicates the similarity of this state with the edge state. In state 1.2, on the contrary, the electron density is localized mainly in the effective quantum dot and goes slightly beyond its limits.

Another new state (line 1.3 in Fig. 10) can not be strictly attributed to any group since it arises at large interaction potential, when $v_{0}>1.725$. We conventionally classify it to the first group since its energy is lower than $2v_{0}$.

The spectrum of the bound states of the second group is also changed compared with the case of $a>2$. Two states (branches 2.1 and 2.2 in Fig. 10) slightly increase their energy at small $v_{0}$. With increasing $v_{0}$, the energy of the state 2.1 becomes a non-monotonic function of $v_{0}$ as it is shown in Fig. 10. In addition, a new branch (line 2.3) appears when $v_{0}>1$ with the energy close to $2v_{0}$.

VII Bound states in the topologically trivial phase: band inversion effect

To complete the picture of the two-particle bound states in the BHZ model we present here the results of the study in the case of a topologically trivial phase. In this case, the parameter $\lambda$ in Eqs (11), (10) should be set equal to $\lambda=+1$ and the calculations are carried out similarly to those described in Secs II, IV, V. When $\lambda=1$, the fundamental solutions $\mathcal{F}_{m}(Qr)$ that define the components of the spinors Eq. (6) are expressed via the Bessel functions of real arguments for any value of the band coupling parameter $a$. Therefore the results do not dramatically depend on $a$. The main result is that two-electron bound states can also exist in topologically trivial phase. Their spectra are not very different from those in the topological phase at $a>2$, but the composition of the spinor components is very different in some cases. Below we demonstrate this for the singletlike states.

The spectrum of the singletlike states is presented in Fig. 13 for the same parameters ($a=2.1$, $r_{0}=2.0$ and $m=0$) as in the case of the topological phase in Fig. 3. It is seen that the spectrum of the bound states and the dependence of their energy on the potential amplitude in both cases are qualitatively similar. The energy of the first group state (line 1.1) is somewhat larger than that in the topological phase (line 1.1 in Fig. 3). The energies of the states of the second group (lines 2.1 and 2.2) also differ not strongly.

The spatial distribution of the spinor component densities in the state of the first group is shown in Fig. 14. The main difference from the topologically nontrivial case is that the spinor component $\psi_{8}$ greatly predominates over the others.

The states of the second group do not so strongly differ in the ratio of the spinor components from the topological case. In Fig. 15 we present only the spatial distribution of the total electron density for both branches (2.1 and 2.2) of the spectrum.

VII.1 Band inversion effect

Here we discuss why the energies of the bound states of the first group in the trivial and topological phases are different and what role the band inversion plays in their formation.

Let us turn again to the mechanism which is commonly used to explain the formation of two-electron bound states in a periodic potential. Two-electron bound states were first discovered in connection with inverse hydrogen absorption spectra observed experimentally [26]. From the very beginning it was supposed that the bound states are formed due to a negative single-particle energy dispersion near the top band boundary [26]. This idea was developed within several simple models of periodic potential for one-dimensional [27] and 2D systems [28, 29]. This mechanism was used to explain the pairing of repulsive ultracold rubidium atoms in an optical lattice [30] and the observations of two-electron peaks in the photo-double-ionization spectra of aromatic hydrocarbons [31]. According to this mechanism the binding energy is determined by the effective mass near the top band boundary of the single-particle spectrum.

It is evident that the results of our calculations do not agree with this idea at least in the case of a topologically nontrivial phase where electron and hole bands are inverted. Indeed, according to the BHZ model the effective mass close to the top of the valence top in the topological phase is larger than in the trivial one. When $|a|>\sqrt{2}$, the ratio of effective masses in topological and trivial phases is equal to $m_{top}/m_{triv}=(a^{2}+2)/(a^{2}-2)$, which equals $\approx$2.66 at $a$=2.1. Therefore, it can be expected that the energies of bound states of the first group in the trivial and topological phases will differ significantly. On the contrary, the calculations show that these energies are unexpectedly close (cf. lines 1.1 in Figs. 3 and 13). We argue below that this contradiction indicates that the reduced effective mass, which determines the energy of the bound state, differs significantly from the reduced effective mass, which is determined by the single-particle dispersion, especially in the topological phase where the electron and hole bands are inverted.

First, consider the topologically trivial case. Figure 14 shows that the bound state is mainly formed by the single-particle states of the hole band, $|H\uparrow,H\downarrow\rangle$. This is consistent with the structure of the band states near the valence band top. It is well known that in the trivial case these states are formed also by the hole band states. Hence our calculations confirm the known point of view that the two-electron bound states are formed by the single-particle states near the top of the valence band.

Now let us turn to the topologically nontrivial phase. From Fig. 4 it is seen that the bound state is mainly formed by the basis states of both the electron and hole bands, $|E\uparrow,E\downarrow\rangle$ and $|H\uparrow,H\downarrow\rangle$. Furthermore, both bands contribute to the total density almost equally. Direct calculations for a variety of the interaction potential amplitudes (not shown here) confirm this result even in if the binding energy is very small. This is easy to understand, since the energy of the bound state lies in the gap where the electron and hole bands overlap, and therefore the electron-band states largely contribute to the total wave function. In contrast, in the trivial case the bands do not overlap.

It should be noted that the spinor structure of the bound state is not consistent with that of the band states at the top of the valence band. In the case of inverted bands, the valence band states are well known to be formed by the electron band states. This is easy to see directly from the Hamiltonian (2) at $a>\sqrt{2}$ and $k\to 0$.

Since the spinor structures of the bound states and valence band states are very different, there are no arguments to think that the reduced effective masses are the same in both cases. Because the bound states are formed by a mixture of the states of the electron and hole bands and their products (specifically, $|E\uparrow,E\downarrow\rangle$, $|E\uparrow,H\downarrow\rangle$, $|H\uparrow,E\downarrow\rangle$, and $|H\uparrow,H\downarrow\rangle$), the reduced effective mass depends on the weights of the components, which should be found by a direct solution of the two-particle Schrödinger equation with a given interaction potential. It is evident that these weights are dependent on the profile and amplitude of the interaction potential.

VIII Conclusions

In this paper we have investigated the two-body problem for 2D electron system described by the symmetric BHZ model in the case of both topologically nontrivial and trivial phases. The main conclusion is that the interaction between the electrons leads to the formation of two-particle bound states at any sign of the pair interaction potential. The pairing of electrons under the action of a repulsive potential becomes possible due to the formation of a negative reduced effective mass of two electrons. The two-electron bound states have the charge $2e$ like the Cooper pair, but their energy lies in the gap of the band spectrum. In this respect they are akin to excitons.

In the case where the spin $S_{z}$ in conserved in the single-particle states, the two-particle states are classified according to the moments of their constituent electrons as a singletlike state (with opposed moments of the electrons) and two tripletlike states (with parallel electron moments).

The bound state spectrum has been studied in a simplified case of zero total momentum of the pair. In general, the spectrum is dependent on the total momentum since the relative motion of the electrons is coupled to the motion of the center of mass. The bound state energy lies in the gap of the two-electron band spectrum. Since the interaction potential is a function only of the distance between electrons, the bound states are specified by an angular quantum number $m$. General properties of the bound state spectrum are found with using a steplike model potential. The states are well classified into two groups in the case of small interaction potential amplitude when they noticeably differ in the energy. The states of the first group have the energy close to the bottom of the gap of the band spectrum, and the energy of the second group states lies near the middle of the gap.

In the trivial phase, the states of the first group are mainly formed by the basis states in which both electrons are mainly in the hole-band states, such as $|H\uparrow,H\downarrow\rangle$. In this case, it is obvious that the reduced effective mass is negative. At a given angular quantum number the bound states are specified by the radial quantum number, so that there is a series of these states. The states of the second group are mainly formed by the basis states composed of different bands such as $|E\uparrow,H\downarrow\rangle$. In these states, one of the electrons is the electron band and the other is in the hole band. It is clear that in this case the reduced effective mass also can be negative. Interestingly, in this group of states there are two states at a given $m$ and no other states appear with increasing the interaction radius $r_{0}$ in the range we have studied, though the bound state energies depend on $r_{0}$. Only one of the state can disappear with increasing $r_{0}$.

In the topologically nontrivial phase, the situation is more complicated. It depends on the coupling of the electron and hole bands. When the coupling parameter is large, $|a|>2$, the bound state spectrum at small $v_{0}$ is qualitatively similar to that in the trivial case, though the energies are noticeably changed. However, the electronic structure of the first group states is changed dramatically. In contrast to the trivial case, these states are formed by the basis states of both bands, to be exact, by the states $|E\uparrow,E\downarrow\rangle$ and $|H\downarrow,H\uparrow\rangle$. That is the bound states are a superposition of the states in which both electrons are in the hole and electron band states. Due to this fact the binding energy turns out to substantially increase as compared with the trivial case.

When the coupling of the bands is not strong, $|a|<2$, the situation changes radically as a consequence of the fact that the evanescent states forming the bound state in the gap contain an oscillating component. We demonstrate this by a detailed study of the case where $a=\sqrt{2}$, which allows one to find the solution exactly. In this case the single-particle spectrum is nearly flat at the band boundaries. The bound state spectra are strongly changed in both groups. The main effect is that new bound states arise in the spectrum in addition to the states of the same type as in the case of $|a|>2$. Of particular interest is the new state appearing in the first group. The new state has a much higher binding energy and arises at much lower interaction potential than other states. This fact shows that the band inversion can favor pairing the electrons when the band coupling is not strong.

The mechanism of the band-inversion impact on the formation of two-particle bound states is caused by two factors: a strong change in the composition of the basis states, which mainly form a given bound state because of the band inversion, and the appearance of oscillating evanescent states.

In the trivial phase, the states of the first group are formed mainly by the basis states of the hole band, such as $|H\uparrow,H\downarrow\rangle$. In contract, in the inverted-band case with strong coupling, $|a|>2$, the bound states of the first group are mainly formed by the basis states of both the electron and hole bands, such as $|H\uparrow,H\downarrow\rangle$ and $|E\uparrow,E\downarrow\rangle$, even if the interaction potential is small.

In the case of nearly flat bands, $|a|=\sqrt{2}$, new addition states arise in the bound state spectrum. In the new state, the weight of the mixed basis states, such as $|E\uparrow,H\downarrow\rangle$, is noticeably increased. Since the bound states are formed by the mixture of the two-particle basis states, which strongly differs from that forming the conduction and valence bands, the reduced effective mass, which appears in the bound state formation, can be essentially different from the reduced effective mass determined by the band spectrum. The weights of the basis states in the bound state are determined by the solution of the Schrödinger equations, like Eqs. (10) and (11), with a given interaction potential. Unfortunately, we failed to find any general relationships for the reduced effective mass in the two-particle bound state.

To complete the picture it is interesting to study the case when the single-particle spectrum in the bands is of a mexican-hat shape. However, this situation requires a separate careful study because of strong singularity of density of states, which appears in this case at the boundaries of the single-particle band spectrum. It is well known that such singularity remarkably facilitates the pairing of electrons [32, 33, 34, 35, 36].

Another interesting point refers to the life time of the two-electron bound states. Since the bound states have the energy in the band gap, they can decay into the states of noninteracting electrons under the action of external disturbances. However, the probability of this decay is greatly reduced if the lower band is filled by electrons. In this connection an important question arises about the life time of the bound states in the presence of many electrons. This question also requires a separate study.

Acknowledgements.

This work was supported by Russian Science Foundation under the grant No. 16-12-10335 in the part related to the general properties of two-particle states (Sec. II) and numerical calculations, and by Russian Foundation for Basic Research under the grant No. 17-02-00309.