New solutions for graphene with scalar potentials by means of generalized intertwining

TL;DR

This paper develops a generalized intertwining method to analytically solve the two-dimensional Dirac equation with scalar potentials, relevant for modeling electrons in graphene under external electrostatic fields.

Contribution

It introduces a novel generalized intertwining approach to find exact solutions for the Dirac equation with scalar potentials in graphene.

Findings

Analytical solutions for wave functions in scalar potential fields.

Explicit examples demonstrating the method's application.

Potential for modeling electronic properties in graphene.

Abstract

The intertwining relations between superpartner Hamiltonians are the main ingredients of well known Supersymmetrical Quantum Mechanics (SUSY QM). In the present paper, the generalized form of intertwining is used for investigation of a massless (zero energy) two-dimensional Dirac equation with scalar external potential. This equation is related to the description of graphene and some other materials in the field of external electrostatic potential. The use of modified intertwining relations allows to find analytically solutions for the wave functions in the field of some external scalar potentials which depend on both space coordinates. A few examples of this construction are given explicitly.

Equations

$$\{Q^{+},Q^{-}{\}}_{+}=\hat{H};[\hat{H},Q^{\pm }]=0.$$

$$\hat{H}=\left(\begin{array}[]{cc}h^{(0)}&0\\ 0&h^{(1)}\\ \end{array}\right);\quad Q^{+}=\left(\begin{array}[]{cc}0&0\\ q^{-}&0\\ \end{array}\right);\quad Q^{-}=(Q^{+})^{\dagger}=\left(\begin{array}[]{cc}0&q^{+}\\ 0&0\\ \end{array}\right).$$

$$h^{(0)}$$

$$h^{(0)}{q}^{+}$$

$$h^{(0)}{q}^{+}={\tilde{q}}^{+}{h}^{(1)};h^{(1)}{\tilde{q}}^{-}=q^{-}{h}^{(0)},$$

$$D_1{\Psi }^{(1)}(\vec{x})=0;D_1\equiv (-i{\sigma }_1{\partial }_1-i{\sigma }_2{\partial }_2+V_1(\vec{x})),{\partial }_i\equiv \frac{\partial }{\partial x_i},i=1,2,$$

$$D_{1}N=M{D}_2;D_2=-i{\sigma }_1{\partial }_1-i{\sigma }_2{\partial }_2+V_2(\vec{x})$$

$$N(\vec{x})\equiv A_n{\partial }_n+A(\vec{x});M(\vec{x})\equiv B_k{\partial }_k+B(\vec{x}),$$

$${\Psi }^{(1)}(\vec{x})=N(\vec{x}){\Psi }^{(2)}(\vec{x}).$$

$$A_1=a_0^{(1)}+{\vec{a}}^{(1)}\cdot \vec{\sigma };$$

$$A_2=a_0^{(2)}-a_2^{(1)}\cdot {\sigma }_1+a_1^{(1)}\cdot {\sigma }_2+a_3^{(2)}\cdot {\sigma }_3;$$

$$B_1=a_0^{(1)}+a_1^{(1)}\cdot {\sigma }_1-a_2^{(1)}\cdot {\sigma }_2-a_3^{(1)}\cdot {\sigma }_3;$$

$$B_2=a_0^{(2)}+a_2^{(1)}\cdot {\sigma }_1+a_1^{(1)}\cdot {\sigma }_2-a_3^{(2)}\cdot {\sigma }_3,$$

$$A(\vec{x})\equiv a_0(\vec{x})+\vec{a}(\vec{x})\cdot \vec{\sigma };B(\vec{x})\equiv b_0(\vec{x})+\vec{b}(\vec{x})\cdot \vec{\sigma }$$

$$V_{\pm }(\vec{x})\equiv V_1(\vec{x})\pm V_2(\vec{x}).$$

$$2b_1(\vec{x})$$

$$A(\vec{x})$$

$$f_0(\vec{x})\equiv i{b}_0+a_1^{(1)}{V}_1(\vec{x});f_3(\vec{x})\equiv b_3-a_2^{(1)}{V}_1(\vec{x}).$$

$$y_2\equiv \alpha x_1+\beta x_2;y_1\equiv -\alpha x_2+\beta x_1;{\partial }_1=\alpha {\partial }_{y_2}+\beta {\partial }_{y_1};{\partial }_2=\beta {\partial }_{y_2}-\alpha {\partial }_{y_1},$$

$$\alpha \equiv a_3^{(1)}+i{a}_0^{(2)};\beta \equiv a_3^{(2)}-i{a}_0^{(1)},$$

$${\partial }_{y_2}{f}_0(\vec{x})-{\partial }_{y_1}{f}_3(\vec{x})=0.$$

$$f_0(\vec{x})={\partial }_{y_1}f(\vec{x});f_3(\vec{x})={\partial }_{y_2}f(\vec{x}).$$

$$({\alpha }^2+{\beta }^2)(\alpha {\partial }_1+\beta {\partial }_2)V_{-}(\vec{x})$$

$$({\alpha }^2+{\beta }^2)(\alpha {\partial }_2-\beta {\partial }_1)V_{+}(\vec{x})$$

$$({\alpha }^2+{\beta }^2)V_{+}(\vec{x})V_{-}(\vec{x})$$

$$a_0^{(1)}=a_0^{(2)}=0;{\alpha }^2+{\beta }^2=(a_3^{(1)}{)}^2+(a_3^{(2)}{)}^2=2.$$

$${\Psi }_A={\Psi }_A(z);{\Psi }_B={\Psi }_B(\bar{z}).$$

$$V_2(\vec{x})=0;V_{+}(\vec{x})=V_{-}(\vec{x})=V_1(\vec{x})$$

$$V_1(\vec{x})=c\cdot \exp (f(\vec{x})),$$

$$({\partial }_1^2+{\partial }_2^2)f(\vec{x})=-c^2\cdot \exp (2f(\vec{x})).$$

$$f(\vec{x})=-\ln{\biggl{(}(\vec{x}+\vec{\tau})^{2}+\gamma^{2}\biggr{)}}+\ln{(\frac{2\gamma}{c})},$$

$$V_1(\vec{x})=2\gamma \frac{1}{(\vec{x}+\vec{\tau }{)}^2+{\gamma }^2}.$$

$$\displaystyle(A_{1}\partial_{1}+A_{2}\partial_{2})=\left(\begin{array}[]{cc}(\alpha+i\beta)\partial+(\alpha-i\beta)\bar{\partial}&2\bar{\mu}\partial\\ 2\mu\bar{\partial}&-(\alpha+i\beta)\partial-(\alpha-i\beta)\bar{\partial}\\ \end{array}\right);$$

$$\displaystyle A(\vec{x})=\left(\begin{array}[]{cc}-(if_{0}+f_{3})&\frac{1}{2}V_{1}(\vec{x})(\beta+i\alpha)\\ \frac{1}{2}V_{1}(\vec{x})(\beta-i\alpha)&-(if_{0}-f_{3})\\ \end{array}\right);\,\mu\equiv a_{1}^{(1)}+ia_{2}^{(1)};\,\bar{\mu}\equiv a_{1}^{(1)}-ia_{2}^{(1)}.$$

Full text

**New solutions for graphene with scalar potentials by means of generalized intertwining

** M. V. Ioffe*1,aaaE-mail: ioffe000@gmail.com; corresponding author, D. N. Nishnianidze2,bbbE-mail: cutaisi@yahoo.com, E. V. Prokhvatilov1,*cccE-mail: evgeniprokhvat@yandex.ru

1 Saint Petersburg State University, 7/9 Universitetskaya nab., St.Petersburg, 199034 Russia.

2 Akaki Tsereteli State University, 4600 Kutaisi, Georgia.

The intertwining relations between superpartner Hamiltonians are the main ingredients of well known Supersymmetrical Quantum Mechanics (SUSY QM). In the present paper, the generalized form of intertwining is used for the investigation of a massless (zero energy) two-dimensional Dirac equation with scalar external potential. This equation is related to the description of graphene and some other materials in the field of external electrostatic potential. The use of modified intertwining relations allows to find analytically solutions for the wave functions in the field of some external scalar potentials which depend on both space coordinates. A few examples of this construction are given explicitly.

Keywords: intertwining relations, two-dimensional Dirac equation, graphene, supersymmetrical Quantum Mechanics.

PACS: 03.65.Fd; 73.22.Pr; 11.30.Pb; 03.65.-w;

1 Introduction.

During last decades, starting from the pioneer paper of Witten in 1981 [1], the supersymmetrical approach became a new interesting and effective study method in investigation of various classical and modern problems of Quantum Mechanics. This approach combined the characteristic features of two old well known methods: Darboux Transformations in Mathematical Physics [2], [3] and Factorization Method of Schrödinger [4], [5] in Quantum Mechanics (see details in [6] - [10]). Their combination was expressed in the algebraic form of SUSY QM where both commutation and anticommutation relations connect three superoperators - Superhamiltonian $\hat{H},$ and Supercharges $Q^{+},\,Q^{-}:$

[TABLE]

The superalgebra (1) admits a lot of different realizations, and each of them corresponds to a specific kind of quantum systems such as models with different space dimensionality, stationary/non-stationary models, models with scalar/matrix interaction, systems with Hamiltonians of second/first order in momenta, discrete/continuous spectra, one/many-particle systems etc. [11] - [20]. The simplest possible realization of the algebra (1) has the form:

[TABLE]

In particular, such realization was used in the framework of standard one-dimensional stationary scalar Quantum Mechanics, where the partner Hermitian Hamiltonians $h^{(0)},\,h^{(1)}$ are the second order differential operators and the components $q^{\pm}$ of supercharge are of first order in momenta. The latter ones play a dual role: they simultaneously factorize Hamiltonians and intertwine them,

[TABLE]

Just the intertwining relations (4) provide the relations between wave functions and the corresponding spectra of two spectral problems - with Hamiltonians $h^{(0)}$ and $h^{(1)}$ [6] - [10].

It is important that the intertwining relations (4) be considered separately when no factorization similar to (3) is fulfilled. For example, when both Hamiltonians and intertwining operators are differential operators of the same degree in derivatives [21], [22]. Exactly in this sense, the intertwining relations play a significantly more important role in SUSY QM.

In the present paper, we will go further by modifying essentially the SUSY intertwining themselves. The intertwining operators $q^{+},\,q^{-}$ in (4) will not be now the same in l.h.s and r.h.s., namely:

[TABLE]

with different $q^{\prime}s$ and $\widetilde{q}^{\prime}s$. From mathematical point of view, relations of similar form were studied in [23], [24] but in different context. Specifically, we will study solutions of relations (5) with a pair of first order $2\times 2$ two-dimensional matrix differential operators $D_{1},\,D_{2}$ of Dirac form (instead of Schrödinger Hamiltonians $h$ above) and also first order matrix intertwining operators $q^{\pm},\,\widetilde{q}^{\pm}.$

The motivation to use such a generalization is the following. While for the solution of the standard form of spectral problem $H\Psi_{n}=E_{n}\Psi_{n}$ with unknown eigenvalues $E_{n}$ and wave functions $\Psi_{n},$ the usual intertwining relations (4) are adequate, in the case of searching for solutions of $D\Psi=0,$ the use of generalized intertwining of the form (5) is enough. Do such homogeneous equations describe any known physical system? Yes, just such equations - two-dimensional Dirac equations with zero mass in the field of external (scalar or matrix) potential - have to be studied in the framework of theory of graphene and some other materials [25] - [28]. The two-component solutions of these Dirac equations provide the probability distribution of electron carriers. Here, the up and down elements of the Dirac two-component spinor are the wave functions corresponding to two different sublattices in graphene. Such form of Dirac equation was studied in many papers, mainly with external potentials depending only on one spatial variable [29] - [40].

Among others, the methods of (standard) SUSY Quantum Mechanics were also effectively used [41], [42] for the case of one-dimensional potential but recently in [43] for a class of external potentials, solutions with actual dependence on both coordinates were built analytically. In the present paper, the problem with external potentials $V(\vec{x})$ of a scalar (electrostatic) nature will be studied by means of generalized intertwining relations which should significantly increase the number of pairs of partner potentials. Among these potentials one can look for so simple that they are amenable to analytic solution of the corresponding Dirac equation. Then the partner Dirac equation can be also solved analytically due to generalized intertwining relations. The structure of the paper is the following. The next Section contains a detailed description of the generalized variant of SUSY intertwining relations for two-dimensional Dirac operators by means of matrix intertwining of first order in momenta. Several examples are collected in Section 3 where the solutions of both partner Dirac equations can be found analytically. In the conclusions, the results are summarized, and the role of zero modes of intertwining operator is discussed.

2 Method of generalized intertwining relations.

Let us consider the problem mentioned above of the analytical solution of two-dimensional massless (or, equivalently, zero energy) Dirac equation with a scalar potential:

[TABLE]

where the Fermi velocity was taken to unity, $\sigma_{1},\sigma_{2}$ - standard Pauli matrices, $x_{1},x_{2}$ - spatial coordinates, $V_{1}(\vec{x})$ - real scalar potential, and $\Psi^{(1)}(\vec{x})$ is a two-component ”spinor” with components $\Psi_{A}^{(1)}(\vec{x}),\Psi_{B}^{(1)}(\vec{x}).$ The main idea is to look for the partner Dirac operator $D_{2}$ which is intertwined with $D_{1}$ as

[TABLE]

with two different first order differential matrix operators:

[TABLE]

with summation over repeated indices $n,k=1,2.$ In (8), $A_{n},\,B_{k}$ are constant $2\times 2$ matrices but matrices $A(\vec{x}),\,B(\vec{x})$ may depend on coordinates. Let us suppose that the partner real potential $V_{2}(\vec{x})$ is simple enough so that the two-component solutions $\Psi^{(2)}(\vec{x})$ of equation $D_{2}\Psi^{(2)}(\vec{x})=0$ can be found. In such a case, due to generalized intertwining (7), solutions of Eq.(6) can be built as well:

[TABLE]

Thus, the problem can be reformulated as follows. One has to solve the generalized intertwining relations (7), i.e. to find both potentials $V_{1}(\vec{x}),\,V_{2}(\vec{x}),$ all constant matrices $A_{n},B_{k}$ and matrix functions $A(\vec{x}),B(\vec{x}).$ To find solutions, one has to decompose (7) step by step over Pauli matrices and over different partial derivatives.

Equating the coefficients with second derivatives in (7), we obtain:

[TABLE]

where all $a^{\prime}s$ are constant coefficients.

Analogously, equating the coefficients in two relations of intertwining (7) with first derivatives and using (10) and (11), one derives the relations between coefficients of expansions of (8),

[TABLE]

and combinations

[TABLE]

Namely, these relations are:

[TABLE]

The last step is to derive from (8) the relation between terms without derivatives. During this procedure, it is convenient to introduce two functions:

[TABLE]

In terms of linearly transformed coordinates,

[TABLE]

with constant complex parameters $\alpha,\,\beta$ defined by coefficients from (10), (11):

[TABLE]

one can check that these functions satisfy the condition:

[TABLE]

For this reason, $f_{0},\,f_{3}$ are expressed in terms of one function:

[TABLE]

After these preparations, simple but rather long calculations of the terms without derivatives in (8) lead to a compact system of nonlinear differential equations:

[TABLE]

Thus, from the intertwining relations (7), we have the relations (12) and the system of equations (18) - (20) together with (14). Arbitrary solution for potentials $V_{\pm}$ and function $f$ provides the pair of potentials $V_{1},\,V_{2}$ and intertwining operators $N(\vec{x}),\,M(\vec{x})$ which satisfy (7). Nonlinearity of Eq.(20) and the condition of reality of both potentials for (in general) complex constant parameters made the task nontrivial. To solve this system, one can approach by means of some suitable ansatzes - in particular, by choosing some special values of constant parameters above. For example, it is convenient to deal below with real normalized constants $\alpha,\,\beta,$ i.e. with

[TABLE]

3 Examples.

Example 1.

There is a special case when the massless Dirac equation (6) can be solved analytically without problems: the case with zero external potential $V(\vec{x})=0.$ The two-component solution of the corresponding Dirac equation is expressed by arbitrary functions of variables $z=x_{1}+ix_{2}$ or $\bar{z}=x_{1}-ix_{2}:$

[TABLE]

Due to the well known Liouville theorem [44], normalizability of such solution is possible only for the restricted regions of the plane with suitably chosen boundary conditions. Physically, this corresponds to the well known Klein paradox [25] - [28], [43].

Let us consider just the case:

[TABLE]

with parameters from (21). It follows from (18), (19) that:

[TABLE]

with real constant $c,$ and therefore, from (20) - that $f(\vec{x})$ must satisfy the nonlinear equation:

[TABLE]

This equation is familiar in the burning theory, and for real values of $c$ two different solutions are known [45].

First solution

For the first of them,

[TABLE]

with arbitrary constants $\tau_{1},\tau_{2},\gamma$, one obtains:

[TABLE]

Due to generalized intertwining relations (7), solutions of the Dirac equation (6) with this potential can be obtained according to (9). It is necessary to use the expressions derived above for coefficients $A_{1},\,A_{2}$ and for function $A(\vec{x}):$

[TABLE]

Taking into account the simplifying ansatz (21) and explicit analytical expressions (25), (26), after the straightforward calculations, one obtains from (22) both components of solution with potential $V_{1}(\vec{x})$:

[TABLE]

Here, $z,\bar{z}$ were defined above, the constants $\tau\equiv\tau_{1}+i\tau_{2},\,\bar{\tau}\equiv\tau_{1}-i\tau_{2}$ can be made equal zero by means of translation of $\vec{x}$, and the constants $\alpha,\,\beta,\mu,\,\bar{\mu}$ are still arbitrary. The initial solutions for zero potential $\Psi_{A}(z)$ and $\Psi_{B}(\bar{z})$ can be also chosen as arbitrary functions.

Let us illustrate this example by formulation of the possible boundary problem for this model on an upper half-plane $x_{2}\geq 0.$ The general boundary conditions [46] - [48] for this domain have the form of linear combination of the components $\Psi_{A}(x_{1},x_{2}=0)$ and $\Psi_{B}(x_{1},x_{2}=0):$

[TABLE]

where $a_{A},a_{B}$ are complex constants. It is necessary to find such functions $\Psi_{A}(z),\,\Psi_{B}(\bar{z})$ in (22), that after their substitution into the r.h.s. of (33), (34), the condition (35) will be fulfilled on the line $x_{2}=0.$ Direct substitution leads to:

[TABLE]

which can be solved by the following choice in (22):

[TABLE]

Second solution

The second solution of (24) is:

[TABLE]

where both components of the constant vector $\vec{\tau}\equiv(\tau_{1},\,\tau_{2})$ and $\gamma$ are arbitrary real constants, and the expression for the partner potential is given by (23):

[TABLE]

In this case, the solution actually depends on the projection of $\vec{x}$ along direction $\vec{\tau}.$ Calculation of the components of solutions of (6) gives:

[TABLE]

Example 2.

Let us choose the case when from the very beginning, the difference $V_{-}=V_{1}-V_{2}$ actually depends only on one direction in the plane:

[TABLE]

Then, (18) provides that also $f=f(y_{1})$ depends only on $y_{1},$ and from (20), one obtains that not only $V_{-},$ but also $V_{+}$ depend only on $y_{1}:$

[TABLE]

as well. Finally, the differential equation (19) can be integrated providing:

[TABLE]

It is convenient now to define new function:

[TABLE]

so that:

[TABLE]

By means of the change of function, the latter equation takes the form of the well known Riccati equation:

[TABLE]

which is solvable [49], [50] for some kinds of coefficient function $V_{2}(y_{1}).$ If one has solution $\chi(y_{1})$ for some specific $V_{2}(y_{1}),$ the partner potential can be built as:

[TABLE]

If one knows not only a pair of partner potentials $V_{1,2}(y_{1}),$ but also the solutions of Dirac equation for one of these potentials $V_{2}(y_{1}),$ the solutions of Dirac equation with the partner potential (48) can be built according to prescription (9).

First solution

The initial massless Dirac equation is easily solved not only for vanishing scalar potential as in Example 1 but also for a constant potential $V(\vec{x})=C.$ In this case it is reduced to the homogeneous Helmholtz equation for components of the wave function:

[TABLE]

Solutions of the homogeneous Helmholtz equation are well known [51], they depend crucially on the boundary value problem for some domain in the plane and on the chosen system of coordinates.

We use the fact that Eq.(47) is solvable for $V_{2}(y_{1})=\sqrt{2}c=const.$ The analytic expressions for solutions $\chi(y_{1})$ are different depending on the sign of the constant $(\nu^{2}-4c^{2}).$

For the positive $\lambda^{2}\equiv(\nu^{2}-4c^{2})>0:$

[TABLE]

for the opposite sign, $\kappa^{2}\equiv(4c^{2}-\nu^{2})>0:$

[TABLE]

and for $\nu^{2}=4c^{2}:$

[TABLE]

By means of translations along $y_{1},$ these solutions can be simplified correspondingly as:

[TABLE]

with the partner potentials defined explicitly by Eq.(48).

Second solution

Let us again consider the case when both potentials $V_{1,2}(\vec{x})$ and $f(\vec{x})$ depend only on $y_{1}$ (see Eqs.(41) – (48)) but with another choice for $V_{2}(y_{1})$ which admits solvability of (47). Two such solutions have relatively simple form:

[TABLE]

and

[TABLE]

Other solutions of Eq.(47) with the same form of $V_{2}$ exist but they lead to much more complicated expressions for $V_{1}$.

By means of generalized intertwining relations (7) and according to (9), solutions $\Psi^{(1)}_{A,B}(\vec{x})$ of the Dirac equation (6) with potential (51) can be obtained from the solutions $\Psi^{(2)}_{A,B}(\vec{x})$ of the Dirac equation with potential (50). The matrix differential operator which transforms $\Psi^{(2)}_{A,B}(\vec{x})$ to $\Psi^{(1)}_{A,B}(\vec{x})$ is:

[TABLE]

whose constant $\mu$ was defined in (32). Fortunately, the components $\Psi^{(2)}_{A,B}$ can be extracted from the papers [33], [52], [36], where solutions of Dirac equation for potential $V_{2}(x_{1})=\frac{(\nu-1)}{\sqrt{2}\cosh(x_{1})}$ of the form analogous to (50) (but depending on $x_{1}$) were obtained explicitly. The Dirac operator $D_{1}$ in (6) with potential $V_{1}(y_{1})$ depending only on $y_{1}$ and the Dirac operator in [52] with potential $V_{1}(x_{1})$ depending only on $x_{1}$ are connected by an unitary matrix. Therefore, the same unitary operator transforms solutions of [52] into solutions $\Psi^{(2)}_{A,B}(y_{1}),$ the latter are necessary to obtain $\Psi^{(1)}_{A,B}(y_{1})$ according to (9).

4 Conclusions.

It is known that the two-dimensional massless Dirac equation with potentials of different nature plays an important role in the description of graphene and some other materials [25] - [28]. Up to now, analytical solutions of such Dirac equation were obtained [29] - [40] for the restricted class of scalar (electrostatic) potentials by different approaches including methods of SUSY Quantum Mechanics [41], [42]. As a rule, these potentials depend only on one coordinate on the plane, and the corresponding one-dimensional part of the wave function is chosen to be normalizable. Meanwhile, the full wave function is non-normalizable due to plane wave multiplier in the second coordinate, this fact is in accordance with the well known Klein paradox [25] - [28]. In [43], the case of potentials with non-trivial dependence on both coordinates was studied, and solutions (also non-normalizable on the whole plane) were found for a class of such potentials.

In the present paper, the new class of two-dimensional potentials in the massless Dirac equation was considered. These potentials were obtained as solutions of modified intertwining relations with first order matrix intertwining operators. For several specific ansatzes, the solutions of Dirac equation were built analytically. In general, these solutions are also non-normalizable on the whole plane. In particular, for anzatses with dependence on variable along the single line of the plane, the wave functions include the non-normalizable plane wave multiplier along the orthogonal direction like in [41] - [52]. The wider class of model potentials amenable to solvability of the corresponding Dirac equation will be useful for further study of physically reproducible systems.

In conclusion, one more point must be discussed. Someone may be interested in the question whether there are other solutions of Dirac equation (6) besides those constructed above from generalized intertwining relations (7). In the standard SUSY QM the answer depends on existence of zero modes of the intertwining supercharge operator. If such zero modes exist, they provide the difference between discrete spectra of partner Hamiltonians and give some additional wave functions (see details in [21] and [6] - [10]). In the case of present paper, the answer is quite definite. Let us consider the equation Hermitian conjugate to (7):

[TABLE]

It is clear that generally speaking, it can produce additional solutions of the Dirac equation (6): an arbitrary zero mode of intertwining operator $M^{\dagger}$ may provide a new solution. This is just analogous to situation in standard SUSY QM. Thus, we must look for solutions $\Psi(\vec{x})$ which satisfy simultaneously:

[TABLE]

We skip for brevity the detailed analysis of the intertwining relations (7) together with relations of (56). The result is the following. For the existence of such solutions $\Psi(\vec{x}),$ among other relations, the functions $f_{0}(\vec{x}),\,f_{3}(\vec{x})$ given in (14) and (17) must satisfy the following equation:

[TABLE]

One can check that (57) is not satisfied in all examples of Section 3. Therefore, no zero modes of $M^{\dagger}$ give new solutions of Dirac equation, and the variety of solutions of (6) is exhausted by solutions built in Section 3.

5 Acknowledgments

The work of M.V.I. was supported by RFBR Grant No. 18-02-00264-a. We are grateful to Dr. D.Candido for information about paper [48].