Spectral analysis of Dirac materials in position-dependent magnetic and electric fields via Heun functions

TL;DR

This paper investigates the spectral properties of Dirac materials under position-dependent electromagnetic fields, utilizing Heun functions to solve the resulting differential equations and determine bound state solutions.

Contribution

It introduces a method to analyze Dirac materials with variable fields using Heun functions, providing new analytical solutions for the spectral problem.

Findings

Decoupled differential equations for Dirac Hamiltonian components.

Derived quantization relations for bound states.

Analytical solutions obtained using Heun functions.

Abstract

This work focuses on the study of the spectral problem for Dirac materials immersed in position-dependent magnetic and electric fields. To achieve this, the system of differential equations satisfied by the eigenfunction components of the Hamiltonian has been decoupled, and the solutions for some specific cases have been analyzed using Heun functions, which provide us with a quantization relation and allow us to determine the solutions for the bound states.

Equations

$${\mathcal{H}}_0=\nu \left[v_x{p}_x{\sigma }_x+v_y{p}_y{\sigma }_y+v_{\mathrm{t}}{p}_y{\sigma }_0\right],$$

$$\mathcal{H}=\nu \left[v_x{p}_x{\sigma }_x+v_y\left(p_y+{\mathcal{A}}_y(x)\right){\sigma }_y+v_{\mathrm{t}}\left(p_y+{\mathcal{A}}_y(x)\right){\sigma }_0\right]-\varphi (x){\sigma }_0,$$

$$\Psi (x,y)=e^{i{k}_{y}y}\bar{\Psi }(x),$$

$$\left[-i{v}_x\frac{\mathrm{d}}{\mathrm{d}x}{\sigma }_x+v_y\left(k_y+{\mathcal{A}}_y(x)\right){\sigma }_y+\left(v_{\mathrm{t}}{k}_y+v_{\mathrm{t}}{\mathcal{A}}_y(x)-\nu \varphi (x)-\nu E\right){\sigma }_0\right]\bar{\Psi }(x)=0.$$

$$\frac{\mathrm{d}}{\mathrm{d}x}\bar{\Psi }(x)=\left[i{\mathrm{F}}_1(x){\sigma }_x+{\mathrm{F}}_2(x){\sigma }_z\right]\bar{\Psi }(x),$$

$${\mathrm{F}}_1(x)=\nu \frac{E-\nu v_{\mathrm{t}}{k}_y+\varphi (x)-\nu v_{\mathrm{t}}\mathcal{A}(x)}{v_x},{\mathrm{F}}_2(x)=\frac{v_y}{v_x}\left(k_y+{\mathcal{A}}_y(x)\right).$$

$$\left(\mp \frac{\mathrm{d}}{\mathrm{d}x}+{\mathrm{F}}_2(x)\right){\psi }^{\pm }(x)={\mathrm{F}}_1(x){\psi }^{\mp }(x).$$

$${\mathcal{L}}^{\pm }=\left(\mp \frac{1}{{\mathrm{F}}_1(x)}\frac{\mathrm{d}}{\mathrm{d}x}+\frac{{\mathrm{F}}_2(x)}{{\mathrm{F}}_1(x)}\right),$$

$${\mathcal{L}}^{\pm }{\psi }^{\pm }(x)={\psi }^{\mp }(x).$$

$$\left[-\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}+\frac{{\mathrm{F}}_1^{{}^{\prime }}(x)}{{\mathrm{F}}_1(x)}\frac{\mathrm{d}}{\mathrm{d}x}+\left({\mathrm{F}}_2^2(x)-{\mathrm{F}}_1^2(x)\pm \frac{{\mathrm{F}}_1(x){\mathrm{F}}_2^{{}^{\prime }}(x)-{\mathrm{F}}_1^{{}^{\prime }}(x){\mathrm{F}}_2(x)}{{\mathrm{F}}_1(x)}\right)\right]{\psi }^{\pm }(x)=0.$$

$$\varphi (x)=-v_{\mathrm{d}}{\mathcal{A}}_y(x).$$

$${\beta }_{\nu }=\frac{v_{\mathrm{d}}+\nu v_{\mathrm{t}}}{v_y},\kappa =\frac{k{v}_y+{\beta }_{\nu }\left(E-\nu v_{\mathrm{t}}k\right)}{\sqrt{1-{\beta }_{\nu }^2}{v}_x},ϵ={\left(\frac{E+v_{\mathrm{d}}k}{\sqrt{1-{\beta }_{\nu }^2}{v}_x}\right)}^2,\mathrm{where}|{\beta }_{\nu }|<1.$$

$$\Omega =2\frac{v_y}{v_x}\sqrt{1-{\beta }_{\nu }^2}{\mathcal{B}}_0,\bar{\alpha }=-2,\bar{\beta }=\frac{2}{{\beta }_{\nu }}\sqrt{\frac{2ϵ}{\Omega }},\bar{\gamma }=\frac{2ϵ}{\Omega },{\bar{\delta }}^{\pm }=\mp \frac{2}{{\beta }_{\nu }}\sqrt{\frac{2(1-{\beta }_{\nu }^2)ϵ}{\Omega }},$$

$$z=\sqrt{\frac{\Omega }{2}}\left[x+2\sqrt{1-{\beta }_{\nu }^2}\frac{\nu v_{\mathrm{t}}k-E}{v_x{\beta }_{\nu }\Omega }\right],$$

$$\left\{-\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}+\frac{1}{z}\frac{\mathrm{d}}{\mathrm{d}z}+\left[{\left(z+\frac{\bar{\beta }}{2}\right)}^2-\bar{\gamma }+\frac{{\bar{\delta }}^{\pm }}{2z}\right]\right\}{\varphi }^{\pm }(z)=0,$$

$${\varphi }^{\pm }(z)=e^{-\frac{(z+\frac{\bar{\beta }}{2}{)}^2}{2}}{f}^{\pm }(z),$$

$$\left\{\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}+\left(\frac{1+\bar{\alpha }}{z}-\bar{\beta }-2z\right)\frac{\mathrm{d}}{\mathrm{d}z}+\left(\bar{\gamma }-\bar{\alpha }-2-\frac{{\bar{\delta }}^{\pm }+\bar{\beta }\left(1+\bar{\alpha }\right)}{2z}\right)\right\}{f}^{\pm }(z)=0.$$

$$\bar{\gamma }-\bar{\alpha }-2=2n.$$

$$E_n=-k{v}_{\mathrm{d}}+\mu v_x{\left(1-{\beta }_{\nu }^2\right)}^{\frac{1}{2}}\sqrt{n\Omega },n=0,1,....$$

$$f^{\pm }(z)=\sum _{k=0}^n{A}_k^{\pm }{z}^k=\sum _{k=0}^N{C}_k^{\pm }{u}_k(z),u_k(z)=\text{mbox}{H}_{{\tau }_0+k}\left(s_0(z+z_0)\right),$$

$$u_k^{\prime \prime }(z)-2s_0^2(z+z_0)u_k^{\prime }(z)+2s_0^2{\tau }_k{u}_k(z)=0,{\tau }_k={\tau }_0+k,$$

$$u_k^{\prime }(z)=2s_0{\tau }_k{u}_{k-1}(z),s_0(z+z_0)u_k(z)={\tau }_k{u}_{k-1}(z)+\frac{1}{2}{u}_{k+1}(z).$$

$$\sum _{k=0}^N{C}_k^{\pm }\left\{\left[\left(2s_0^2-2\right)z^2+\left(2s_0^2{z}_0-\bar{\beta }\right)z+1+\bar{\alpha }\right]u_k^{\prime }(z)+\left[\left(\bar{\gamma }-\bar{\alpha }-2-2s_0^2{\tau }_k\right)z-\frac{{\bar{\delta }}^{\pm }+\bar{\beta }(1+\bar{\alpha })}{2}\right]u_k(z)\right\}=0.$$

$$s_0=\pm 1,z_0=\frac{\bar{\beta }}{2}.$$

$$\sum _{k=0}^N{C}_k^{\pm }\left[R_k{u}_{k+1}(z)+Q_k^{\pm }{u}_k(z)+P_k{u}_{k-1}(z)\right]=0,$$

$$R_k=\left(n-{\tau }_k\right),Q_k^{\pm }=s_0\bar{\beta }\left({\tau }_k-n+\frac{1\pm {\left(1-{\beta }_{\nu }^2\right)}^{\frac{1}{2}}}{2}\right),P_k=2(n-{\tau }_k-1){\tau }_k.$$

$$P_0=0\text{mbox}and{R}_N=0,$$

$${\tau }_0=0\text{mbox}or{\tau }_0=n-1,$$

$${\tau }_0=n-N.$$

$${\tau }_0=n-1,R_k=\left(1-k\right),Q_k^{\pm }=s_0\frac{\bar{\beta }}{2}\left(2k-1\pm {\left(1-{\beta }_{\nu }^2\right)}^{\frac{1}{2}}\right),P_k=2k(1-n-k).$$

Full text

Keywords: Dirac materials, exact solutions, Heun funtion.

1 Introduction

In recent years, Dirac materials have captured the attention of numerous researchers, not only in condensed matter physics but also in other branches of physics. This interest stems from the fact that equations that describe their charge carriers in a low energy regime (near the Fermi level) can be used to analyze analogous physical systems or to explore the possibility of finding exact solutions applicable to various types of interactions [1, 2, 3, 4].

For example, a few years ago, supersymmetric quantum mechanics (SUSY QM) – a well-established formalism for analytically solving physical systems described by Schrödinger-like differential equations – was successfully applied to the study of interactions between monolayer and bilayer graphene electrons and perpendicular magnetic fields [5, 6, 7]. This approach enabled the derivation of exact solutions for configurations involving inhomogeneous magnetic fields, as the resulting equations exhibit a direct connection with the factorization method employed in SUSY QM. However, additional interactions, such as external electric fields, have since been incorporated into the physical models of graphene, making the analysis through exact solutions considerably more challenging [8, 9, 10]. Furthermore, the study of graphene has opened the door to the research of other two-dimensional Dirac materials, which have emerged as promising candidates for technological applications in fields such as valleytronics and spintronics [11, 12, 13, 14]. Consequently, the interaction of these materials with external electric and magnetic fields has regained attention from a mathematical physics perspective.

In this work, we revisit the problem of anisotropic Dirac materials with tilted cones, focusing on the interaction of their charge carriers with position-dependent external magnetic and electric fields. We aim to successfully implement a direct approach in which the energy spectrum and its corresponding eigenfunctions can be determined when more general configurations than those previously reported in the literature are used [15, 16, 17, 18, 19].

The organization of this paper is as follows. In Section 2, we present the effective Hamiltonian to be analyzed and an appropriate way to decouple the set of equations associated with the eigenvalue problem is discussed. Section 3 provides a list of solvable cases using the developed method, along with the polynomial solutions of the Heun equations. Finally, we conclude with remarks on the results obtained in Section 4.

2 The eigenvalue problem

At low energies, a wide variety of two-dimensional materials exhibit linear or quasilinear dispersion relations due to their lattice structure, which also induces anisotropies or Dirac cones tilting. For this reason, they are often referred to as Dirac materials, as their dynamics are described by a Dirac-like Hamiltonian given by

[TABLE]

where $\nu=\pm 1$ represents the valley index, $v_{x},v_{y}$ are the anisotropic Fermi velocities, $v_{\rm t}$ is the velocity term arising from the Dirac cone tilt, $\sigma_{x},\sigma_{y}$ are the Pauli matrices, and $\sigma_{0}$ is the $2\times 2$ identity matrix. When the system is affected by the presence of an external electric field $\mathbf{E}=\mathcal{E}(x)\hat{e}_{x}$ in the plane of the material (the $x$-$y$ plane) and a magnetic field $\mathbf{B}=\mathcal{B}(x)\hat{e}_{z}$ perpendicular to it, both with translational symmetry with respect to the $y$ coordinate, the Hamiltonian in Equation (1) must be modified to include the information of these fields through the minimal coupling rule as

[TABLE]

where $\phi(x)$ is the scalar potential that generates the electric field and $\mathcal{A}_{y}$ is the component of the vector potential in the Landau gauge ($\mathbf{A}=\mathcal{A}_{y}(x)\hat{e}_{y}$) that generates the magnetic field. Due to the translational invariance in the $y$-direction, the Hamiltonian given in Equation (2) commutes with the operator $p_{y}$; therefore, it is possible to express its eigenfunctions following the ansatz:

[TABLE]

with $\bar{\Psi}(x)=\left(\psi^{+}(x),i\psi^{-}(x)\right)^{\rm T}$. Thus, the eigenvalue problem $\mathcal{H}\Psi(x,y)=E\Psi(x,y)$ leads us to the following differential equation for $\bar{\Psi}(x)$:

[TABLE]

To determine the set of eigenfunctions and eigenvalues of the Hamiltonian, it is necessary to find the solutions corresponding to the previous equation. However, this leads us to a system of coupled differential equations for the components $\psi^{\pm}(x)$ of the spinor $\bar{\Psi}(x)$. Hence, in the next section, we will present an approach that allows us to find the solutions exactly.

2.1 Decoupling the differential system

In general, the solutions to Equation (4) strongly depend on the form of the fields. Therefore, we proceed as follows. First, by multiplying Equation (4) by $i\sigma_{x}$ on the left-hand side and algebraically manipulating, the system presented below is obtained:

[TABLE]

where

[TABLE]

The equation above leads us to a set of coupled differential equations for the components $\psi^{\pm}(x)$, given by:

[TABLE]

If we now define the operators $\mathcal{L}^{+}$ and $\mathcal{L}^{-}$ as:

[TABLE]

then the system of equations in (7) can be rewritten as:

[TABLE]

Therefore, by applying the operator $\mathcal{L}^{\mp}$ to both sides of Equation (9), the system of equations decouples, and after some algebra, two decoupled second-order homogeneous linear differential equations for the components $\psi^{\pm}(x)$ are obtained, which turn out to be:

[TABLE]

Note that Equation (10) represents a decoupled system of differential equations, which now allows us to directly find the components of the spinor in Equation (3) for arbitrary field profiles. However, in general, this system is not composed of eigenvalue equations that allow for the direct identification of the energy. Nevertheless, it is possible to obtain exact solutions for certain field profiles, as shown below.

3 Examples of exact solvability

In the following, we analyze some examples of external electric and magnetic field profiles for which it is possible to determine exact solutions of the system of equations in  (10), and thus obtain the energy spectrum and the corresponding eigenfunctions. To this end, we begin by assuming that the magnetic and electric field profiles are proportional to each other, that is, they differ by a multiplicative constant denoted by $-v_{\rm d}$, which is read as

[TABLE]

It can be observed that for $v_{\rm d}=0$, the electric field is zero. If $v_{\rm d}>0$, then both fields point in the same direction (either positive or negative), and if $v_{\rm d}<0$, then one field points in the positive direction while the other points in the negative direction. All three of these cases can be treated without loss of generality by assuming that $v_{\rm d}>0$ and subsequently taking the appropriate limit. For this reason, we will consider in the following discussion that both fields point in the positive direction.

Now, to simplify the expressions and calculations, we first define the following quantities:

[TABLE]

With these elements, we now proceed to analyze three examples of field profiles.

3.1 Constant fields

As a first example, we will consider the standard case on which both the magnetic field and the electric field are constant, namely, $\mathbf{B}=\mathcal{B}_{0}\hat{e}_{z}$ and $\mathbf{E}=\mathcal{E}_{0}\hat{e}_{x}$ with $\mathcal{B}_{0},\mathcal{E}_{0}>0$ (in this way, $v_{\rm d}=\mathcal{E}_{0}/\mathcal{B}_{0}$). Then, the $y$-component of the vector potential can be chosen as $\mathcal{A}_{y}(x)=\mathcal{B}_{0}x$ and $\phi(x)=-\mathcal{E}_{0}x$. Furthermore, we define the following quantities:

[TABLE]

and applying the change of variable given by

[TABLE]

Equation (10) transforms into

[TABLE]

where $\psi^{\pm}(x)=\phi^{\pm}(z)$. It can be observed that asymptotically, Equation (15) behaves like the one corresponding to a harmonic oscillator, which makes it natural to propose the following structure for $\phi^{\pm}$:

[TABLE]

being $f^{\pm}$ functions to be determined while the Gaussian form guarantees that $\phi^{\pm}$ are asymptotically correct solutions. By substituting (16) into (15), the following differential equation for $f^{\pm}(z)$ is obtained

[TABLE]

In the literature, Equation (15) is commonly referred to as the biconfluent Heun differential equation[20, 21], thus, $f^{\pm}(z)=\text{HeunB}(\bar{\alpha},\bar{\beta},\bar{\gamma},\bar{\delta}^{\pm};z)$. However, since the function $\psi^{\pm}$ must be square integrable, it is essential that $f^{\pm}(z)$ be a polynomial of finite order $n$, a condition that is satisfied if and only if

[TABLE]

We must emphasize that even though the differential equations in (10) satisfied by $\psi^{\pm}$ are not eigenvalue equations with energy as the eigenvalue. Actually, Equation (18) acts as a quantization condition that allows us to determine the energy spectrum of the system. Thus, from equations (12) and (13), it follows that

[TABLE]

being $\mu=\pm 1$, where the plus (minus) sign corresponds to the conduction (valence) band. Now, based on the information provided by the quantization relation in (18) and to determine the set of eigenfunctions corresponding to (19), we will show that the solutions $f^{\pm}$ to the Heun equation in (17) can be expressed as a linear combination of Hermite polynomials.

3.1.1 Polynomial solutions of the biconfluent Heun differential equation

By considering a finite-order polynomial $n$ as a solution to the biconfluent Heun differential equation in (17), the relation $\bar{\gamma}-\bar{\alpha}-2=2n$ is obtained. On the other hand, if we assume that such a polynomial solution can be written in terms of Hermite polynomials, then the following condition for $f^{\pm}$ must be satisfied:

[TABLE]

where $C^{\pm}_{k}$ are the coefficients to be determined and $\tau_{0}$, $s_{0}$, $z_{0}$ are constants that we can choose at will as long as the convergence of the function and equality with the sum of the powers are guaranteed.

In this regard, the Hermite polynomials in (20) satisfy the following differential equation

[TABLE]

in addition to the recurrence relations:

[TABLE]

Then, by substituting $f^{\pm}$ from (20) into Equation (17) and using (21) to replace the second-order derivative of $u_{k}$, it follows that

[TABLE]

In this manner, we can use the freedom we have over $s_{0}$ and $z_{0}$ to eliminate the terms proportional to $z^{2}u^{\prime}_{k}$, and $zu^{\prime}_{k}$ by choosing

[TABLE]

Thus, by making use of the relationship between $\bar{\gamma}$, $\bar{\alpha}$, and $n$, as well as equations in (24) and the recurrence relations in (22), Equation (23) leads us to the following equality:

[TABLE]

being

[TABLE]

By considering that the Hermite polynomials constitute a complete basis, the following conditions must be satisfied simultaneously:

[TABLE]

which allows to determine the values of $\tau_{0}$ and $N$, since the first one needs that

[TABLE]

while the second one holds if

[TABLE]

Equations (28) and (29) imply that $N=n$ or $N=1$. Although both choices lead to the same solution expressed in two different bases, the first one involves a more complicated calculation due to the following factors: 1) the upper limit of the sum depends on the degree of the polynomial considered for the solution; 2) regardless of the value of $n$, $P_{n-1}=0$, whereas with the other choice $P_{1}=0$ only occurs if $n=0$. Therefore, determining the coefficients $C^{\pm}_{k}$ is more laborious. For this reason, we will consider only the case $N=1$, then, it turns out that

[TABLE]

Thus, taking into account that $u_{k}(z)$ form a basis, Equation (25) can be expressed as:

[TABLE]

This system has a non-trivial solution if and only if $\mathbf{M}$ is non-invertible, i.e., its determinant vanishes:

[TABLE]

Then, by considering the expressions in (30) and manipulating algebraically, it can be shown that indeed $\det(\mathbf{M})=0$. Hence, $C^{\pm}_{0}$ and $C^{\pm}_{1}$ are not independent but satisfy the following relation:

[TABLE]

However, due to the quantization condition, it follows that $\bar{\beta}\propto\sqrt{n}$, which leads to the following two cases:

• 1

Case $n=0$. For this case, $\bar{\beta}=0$, thus $C^{\pm}_{0}=0$. In this way, the solution is given by

[TABLE]

where $C^{\pm}_{1}$ will be determined by the normalization of the final solution, and the argument was chosen simply as $z+z_{0}$, since for both possible values of $s_{0}$, the Hermite polynomial $\mbox{H}_{0}$ is equal to 1.

• 2

Case $n\geq 1$. For this case, $\bar{\beta}\neq 0$. Thus, according to Equation (33), it must hold that the solutions $f^{\pm}$ to the biconfluent Heun differential equation are given by:

[TABLE]

Notice that, regardless of the choice of $s_{0}=\pm 1$, both values will only generate solutions with a global phase difference. Therefore, without loss of generality, we can choose $s_{0}=1$.

Finally, since $z$ and $z_{0}$ depend implicitly on the energy, and therefore on $n$, it is convenient to define the variable

[TABLE]

In this manner, by taking into account equations (16), (34), and (35), the normalized eigenfunctions corresponding to the energies of Equation (19) are:

[TABLE]

where $\psi^{\pm}_{n}$ can be written as in [15, 17, 22]:

[TABLE]

with $\mathcal{N}_{n}=\frac{1}{\sqrt{2^{n}n!}}\left(\frac{\Omega}{8\pi}\right)^{\frac{1}{4}}$. Thus, knowing that the probability density and probability current for the Hamiltonian in Equation (1) are given by

[TABLE]

for the excited states we also have

[TABLE]

These expressions are also valid for the following cases. A brief discussion of these results is presented below.

3.1.2 Discussion

For the given configuration of electric and magnetic field profiles, the energy spectrum described in Equation (19) is discrete, infinite, and dispersive, exhibiting a linear behavior with respect to the wave number in the $y$-direction. This linearity presents a negative slope corresponding to $-v_{\rm d}$, i.e., the negative ratio of the field amplitudes (see Figure 1a). For that reason, if the electric field is null ($v_{\rm d}=0$), the energy levels become completely flat with respect to $k$ (non-dispersive), regardless of the anisotropic velocities or the value of $v_{\rm t}$. On the other hand, it is worth noting that as the parameter $\beta_{\nu}$ approaches $1$ (or $-1$)—that is, when the electric field amplitude approaches the critical value $\mathcal{E}_{0}=\mathcal{B}_{0}(v_{y}-\nu v_{\rm t})$—the energy levels tend to become degenerate at a single value of $-v_{\rm d}k$ (see Figure 1b).

Finally, while it is evident that the probability density depends on the valley and band indices (see Figure 1 c), a more interesting result arises from the $y$-component of the probability current. Due to the tilt of the Dirac cones, a non-zero current is generated in the ground state regardless of the values of the anisotropic velocities $v_{x}$ and $v_{y}$ (see Figure 1 d), in contrast to materials such as graphene[5].

3.2 Exponential decaying fields

The second example we consider involves both fields exhibiting an exponentially decaying profile, namely, $\mathbf{B}=\mathcal{B}_{0}e^{-\alpha x}\hat{e}_{z}$ and $\mathbf{E}=\mathcal{E}_{0}e^{-\alpha x}\hat{e}_{x}$, where $\mathcal{B}_{0},\mathcal{E}_{0},\alpha>0$. In this configuration $v_{\rm d}=\mathcal{E}_{0}/\mathcal{B}_{0}$, as before. The corresponding potentials can be chosen as $\mathcal{A}_{y}(x)=-\frac{\mathcal{B}_{0}}{\alpha}e^{-\alpha x}$ and $\phi(x)=\frac{\mathcal{E}_{0}}{\alpha}e^{-\alpha x}$. Moreover, if we define the following quantities:

[TABLE]

and considering the change of variable

[TABLE]

then the set of differential equations in (10) is transformed into

[TABLE]

where $\psi^{\pm}(x)=\phi^{\pm}(z)$. It is worth noting that, as $x\to-\infty$, the previous equation behaves like a modified Helmholtz equation in one dimension.Therefore, the solution is asymptotically of the form $\phi^{\pm}(z)\propto e^{\frac{\bar{\alpha}}{2}z}$. Similarly, as $x\to\infty$, the differential operator in Equation (44) behaves like $\frac{{\rm d}^{2}}{{\rm d}z^{2}}+\frac{1}{z}\frac{{\rm d}}{{\rm d}z}-\frac{\bar{\beta}^{2}}{4z^{2}}$, and hence the solution is asymptotically given by $\phi^{\pm}(z)\propto z^{\frac{\bar{\beta}}{2}}$. Based on these asymptotic behaviors, we propose that the solution takes the form

[TABLE]

By substituting (45) into (44), the following differential equation for $f^{\pm}(z)$ is obtained

[TABLE]

Equation (46) is commonly referred as the confluent Heun differential equation[23], so that $f^{\pm}(z)=\text{HeunC}(\bar{\alpha},\bar{\beta},\bar{\gamma},\bar{\delta},\bar{\eta}^{\pm};z)$. As in the previous case, for $\psi^{\pm}$ to be square-integrable, we must require that $f^{\pm}(z)$ be a polynomial of finite degree $n$. This condition is satisfied if

[TABLE]

Similarly to the first case, the previous equation acts as a quantization condition, allowing us to determine the energy values, which are found to be:

[TABLE]

where, once again, $\mu$ denotes the band parameter. By assuming that $f^{\pm}$ admits a polynomial solution, it can be shown that it is possible to express it in terms of known orthogonal polynomials. In this case, the associated Laguerre polynomials offer the best fit, as will be shown below.

3.2.1 Polynomial solutions of the confluent Heun differential equation

In order to demonstrate that the polynomial solution for $f^{\pm}$ can be written as a linear combination of associated Laguerre polynomials, we begin by defining the following parameters:

[TABLE]

In terms of the parameters $\bar{\xi}^{\pm}$ and $\bar{\zeta}^{\pm}$, the quantization condition is given by $\bar{\xi}^{\pm}+\bar{\zeta}^{\pm}=-\bar{\alpha}n$. On the other hand, if we assume that the polynomial solution can be expressed in terms of associated Laguerre polynomials, then the following condition must hold:

[TABLE]

where $C^{\pm}_{k}$ are coefficients to be determined, while $\tau_{0}$, $\rho_{0}$, $s_{0}$, and $z_{0}$ are constants that can be freely chosen, provided they ensure the convergence of the function and consistency with the polynomial expression.

Moreover, the differential equation for the Laguerre polynomials is well known, and therefore the following identity holds:

[TABLE]

and also the recurrence relations

[TABLE]

Thus, by substituting the ansatz in Equation (50) into the confluent Heun differential equation given in Equation (46), and replacing the second derivative of $u_{k}$ using the relation in Equation (51), we obtain—after some algebraic manipulations—that:

[TABLE]

By exploiting the freedom in the choice of parameters, we can set $s_{0}$ and $z_{0}$ in such a way that the terms proportional to $z^{3}u^{\prime}_{k}(z)$ and $u^{\prime}_{k}(z)$ vanish. Thus, we fix:

[TABLE]

With such a selection, a straightforward algebraic manipulation shows that Equation (53) becomes equivalent to:

[TABLE]

Now, using the equations in (52) to replace the terms $u_{k}^{\prime}$, $z\,u_{k}^{\prime}$, and subsequently making use of the freedom in choosing $\rho_{0}$ to eliminate the term proportional to $z\,u_{k}$, equation (55) leads to:

[TABLE]

with $\rho_{0}=n+\bar{\beta}-\tau_{0}-1$ and the terms $P_{k},\,Q^{\pm}_{k},R_{k}$ are given by

[TABLE]

However, since the polynomials $u_{k}$ form a basis, the conditions $P_{0}=0$ and $R_{N}=0$ must be satisfied simultaneously. The first condition allows us to determine the value of $\tau_{0}$, as it is fulfilled when

[TABLE]

while the second condition is satisfied if

[TABLE]

In this way, we observe that the value of $N$ must be either 1 or $-n$. Although both options lead to the same result, the calculation with $N=-n$ is more complicated, as it requires summing over non-positive integers, and the number of terms to be summed is not fixed but depends on the order of the polynomial in question. For this reason, in the present work we will consider only the case $N=1$, which leads to

[TABLE]

Once again, since $u_{k}$ form a basis, the following system is obtained:

[TABLE]

In order to avoid the trivial solution, the condition $\det(\mathbf{M})=R_{0}P_{1}-Q^{\pm}_{1}Q^{\pm}_{0}=0$ must be satisfied. It is straightforward to verify that, using the values defined in (42), this condition holds. Therefore, the coefficients $C_{0}^{\pm}$ and $C_{1}^{\pm}$ are not independent and satisfy the following relation:

[TABLE]

It is important to emphasize that the above condition does not include the case $n=0$, since for this value the condition $1\leq\tau_{k}$ is not satisfied for $k=0,1$. However, this is a trivial case that admits the representation shown in equation (50), provided that $C_{0}^{\pm}=0$. In this way, we once again obtain two cases:

• 1

Case $n=0$. For this case, we have $C^{\pm}_{0}=0$. Then, the solution is given by

[TABLE]

where $C^{\pm}_{1}$ will be determined by the normalization of the final solution.

• 2

Case $n\geq 1$. For this case, $C_{0}^{\pm}=\frac{\bar{\zeta}^{\pm}}{n}C_{1}^{\pm}$. Thus, $f^{\pm}$ is given by:

[TABLE]

It is important to emphasize that although $\bar{\beta}$ is constant with respect to $z$, its value depends on the energy and therefore on $n$, being defined as:

[TABLE]

Furthermore, since the argument of the associated Laguerre polynomials is given by $-\bar{\alpha}z$, it is convenient to define the variable $\theta$ as:

[TABLE]

Thus, by substituting the solutions from Equations (63) and (64), the eigenfunctions corresponding to the energy values given in Equation (48) are obtained as:

[TABLE]

with $\psi^{\pm}_{n}$ given by:

[TABLE]

where the constants $C_{n}$ and $\mathcal{N}_{n}$ are defined as:

[TABLE]

with the factor $I_{n}$ denotes the following integral

[TABLE]

We now present a brief discussion of these states and their corresponding energy levels.

3.2.2 Discussion

To begin, we must note that in this case it is necessary to impose a condition that ensures the energy levels in Equation (50) are real and that the corresponding wave functions in Equation (69) are square-integrable. This condition is given by

[TABLE]

where the dependence of $\kappa$ on $n$ arises from its relation to the energy (see Equation (12)). This condition causes the energy spectrum to become discrete but finite, with the maximum number of excited states determined by the fulfilment of the previous condition. The spectrum exhibits a linear behavior in $k$ with a slope of $-v_{d}$ for the ground state, whereas it is nonlinear for the excited states. Moreover, since the square-integrability condition defines a domain in $k$, it can be observed that, as the number of excited states increases, this domain becomes narrower. Likewise, since the number of excited states is finite, it is possible to define an upper bound for all energy levels whose behavior is linear with respect to $k$, as indicated by the gray dashed line in Figure 2a. On the other hand, as in the previous case, when $\beta_{\nu}$ approaches $1$ (or $-1$) the energy levels tend to become degenerate at a value of $-v_{d}k$ (see Figure 2b). Regarding the energy current, it is also observed that the $y$-component is nonzero due to the tilt of the Dirac cones (see Figure 2d).

Finally, we would like to mention that, by using the identity

[TABLE]

our result coincides with that derived in [16, 17], where a similarity transformation was employed to solve this problem.

3.3 Hyperbolic fields

For this example, we consider hyperbolic magnetic field profiles satisfying the relations $\mathbf{B}=\mathcal{B}_{0}\,\mathrm{sech}^{2}(\alpha x)\,\hat{e}_{z}$ and $\mathbf{E}=\mathcal{E}_{0}\,\mathrm{sech}^{2}(\alpha x)\,\hat{e}_{x}$, with $\mathcal{B}_{0},\mathcal{E}_{0},\alpha>0$. Then, the $y$-component of the vector potential and the scalar potential can be chosen as $\mathcal{A}_{y}(x)=\frac{\mathcal{B}_{0}}{\alpha}\tanh(\alpha x)$ and $\phi(x)=-\frac{\mathcal{E}_{0}}{\alpha}\tanh(\alpha x)$, respectively. Additionally, let us define the following quantities:

[TABLE]

[TABLE]

[TABLE]

Thus, by using the change of variable

[TABLE]

the system of differential equations in (10) becomes:

[TABLE]

where we take $\psi^{\pm}(x)=\phi^{\pm}(z)$ and the coefficients $A_{i}$ are given by

[TABLE]

[TABLE]

Note that $z\to 0$ as $x\to-\infty$. Furthermore, in this limit, the differential operator in (75) takes the form $\frac{{\rm d^{2}}}{{\rm d}z^{2}}+\frac{1}{z}\frac{{\rm d}}{{\rm d}z}-\frac{(\bar{\gamma}-1)^{2}}{4z^{2}}$, so asymptotically the solution must behave like $\phi^{\pm}(z)\propto z^{\frac{\gamma-1}{2}}$ in this region. On the other hand, when $x\to\infty$, we have $z\to 1$, and in this domain limit, the differential operator becomes $\frac{{\rm d^{2}}}{{\rm d}z^{2}}+\frac{1}{z-1}\frac{{\rm d}}{{\rm d}z}-\frac{(\bar{\delta}-1)^{2}}{4(z-1)^{2}}$. Therefore, in this region, the solution should asymptotically behave like $\phi^{\pm}(z)\propto(1-z)^{\frac{\delta-1}{2}}$. Based on all the above, we can propose that the solution to equation (75) takes the form

[TABLE]

Substituting this ansatz in (75), a new differential equation for $f^{\pm}$ is obtained:

[TABLE]

The above equation is known as the general form of the Heun differential equation [23]. Consequently, the function $f^{\pm}(z)=\text{HeunG}(\bar{\alpha},\bar{\beta},\bar{\gamma},\bar{\delta},\bar{q}^{\pm};z)$. As in previous cases, this solution $\psi^{\pm}$ will be square integrable if $f^{\pm}$ is a polynomial of finite degree $n$. For this to hold, the following condition must be satisfied:

[TABLE]

Note that this equation has solutions $\bar{\alpha}=-n$ or $\bar{\beta}=-n$. However, given the definitions of these parameters in (73), both conditions cannot occur simultaneously. Nevertheless, since these parameters can be freely interchanged without affecting any of the previous equations, we may assume, without loss of generality, that $\bar{\alpha}=-n$ and $\bar{\beta}\neq-n$ (or vice versa). Moreover, as in the previous cases, this condition acts as a quantization relation. In contrast to the quantization conditions for constant and exponential fields, this condition is not linear in $n$, but quadratic. As a result, there are four independent solutions for the energy, two of which are:

[TABLE]

with the band index $\mu=\pm 1$ and $n\in\mathbb{Z}_{\geq 0}$. For these solutions to be well-defined and for the spectrum to remain real, it is necessary to impose the following condition:

[TABLE]

Recalling that $\mathcal{F}$ is positive definite, for a given value of $\mathcal{F}$, this condition imposes an upper bound on the number of excited states. Now, returning to the other two possible energy spectra obtained by solving Equation (79), it can be observed that they result from implementing the transformation $n\alpha\rightarrow-n\alpha$ in Equation (80). However, this also modifies the condition required for the energy spectrum to be well-defined and real, demanding now that $\mathcal{F}<(v_{y}/v_{x})k$. To ensure that this condition is always satisfied, the fields would need to vanish or $k\rightarrow\infty$, which contradicts the goal of obtaining stationary states. For this reason, we will henceforth consider only the energy spectrum given in Equation (80), along with the condition specified in Equation (81). In this way, it is easy to verify that $\bar{\alpha}$ and $\bar{\beta}$ are given by

[TABLE]

On the other hand, similarly to the previous cases, the polynomial solution for $f^{\pm}$ can be written in terms of Jacobi polynomials, as we will show below.

3.3.1 Polynomial solutions of the general Heun differential equation

By considering that the polynomial solution to the general Heun differential equation can be expressed in a basis of Jacobi polynomials, the following condition must be satisfied:

[TABLE]

where $\tau_{0}$, $a_{0}$, $b_{0}$, $s_{0}$ and $z_{0}$ are constants that can be freely chosen, provided they ensure the convergence of the solution. Moreover, the functions $u_{k}$ satisfy the following differential equation:

[TABLE]

as well as the recurrence relations given by

[TABLE]

Thus, by inserting the ansatz from Equation (83) into the general Heun differential equation (78), and replacing the second derivative of $u_{k}$ by means of the Equation (84), one obtains, after some algebraic manipulation, the following expression:

[TABLE]

where the constants $A_{i}$ and $D_{i}$ are given by

[TABLE]

Using the freedom in the choice of $s_{0}$ and $z_{0}$ to eliminate the terms proportional to $u^{\prime}_{k}$ and $zu^{\prime}_{k}$ in Equation (86), and the freedom in $a_{0}$ and $b_{0}$ to ensure that the remaining terms share a common factor with $u_{k}$, we propose:

[TABLE]

It is worth mentioning that the previous conditions are not the only ones that lead to the required constraints. In fact, an equally valid set is given by $z_{0}=-\frac{1}{2}$, $s_{0}=\frac{1}{z_{0}}=-2$, $a_{0}=\bar{\gamma}-1$, $b_{0}=\bar{\delta}-1$. However, with these values, the final result will differ by a global phase, which is physically irrelevant. For this reason, and without loss of generality, we will only consider the case described by the conditions in Equation (88). Thus, Equation (86), after some algebraic manipulation, leads to

[TABLE]

Nevertheless, using the recurrence relations in (85) to replace the first derivative and the terms involving $zu_{k}(z)$, the previous equation transforms into

[TABLE]

where $R_{k}$, $P_{k}$, and $Q^{\pm}_{k}$ are given by

[TABLE]

Now, being that the functions $u_{k}(z)$ form a basis, it is necessary to demand that $P_{0}=0$ and $R_{N}=0$. Because the numerator of $P_{0}$ is a fourth-order polynomial in $\tau_{0}$, there will be several solutions; however, only one of them yields an integer value, corresponding to $\tau_{0}=-\bar{\alpha}-1=n-1$. Similarly, the condition $R_{N}=0$ yields four different solutions for $\tau_{0}$, but only two of them are integers, corresponding to

[TABLE]

Thus, $N$ can only take the values $1$ and $-n$. However, as in the previous cases, the value $-n$ leads to the same result but involves more complicated calculations. For this reason, we will only consider the value $N=1$. In this way, we obtain that

[TABLE]

As the functions $u_{k}$ form a basis, the following system is obtained:

[TABLE]

Once again, to avoid the trivial solution, it is necessary that $\det(\mathbf{M})=R_{0}P_{1}-Q^{\pm}_{1}Q^{\pm}_{0}=0$, which can be easily verified by a direct substitution of the parameters involved. In this way, $C^{\pm}_{1}$ and $C^{\pm}_{0}$ are not linearly independent and satisfy the following relation:

[TABLE]

Equations (84) and (95) lead to the following two cases:

• 1

Case $n=0$. In this case, the solution is given by a polynomial of degree zero, and therefore $C^{\pm}_{0}=0$. In this way, the functions $f^{\pm}$ are given by

[TABLE]

where $C_{1}^{\pm}$ is determined by the normalization procedure.

• 2

Case $n\geq 1$. For this case, $R_{0}\neq 0$, and thus we must have $C_{0}^{\pm}=\frac{\mathcal{F}\sqrt{\varepsilon}\left(1\pm\sqrt{1-\beta_{\nu}^{2}}\right)}{n\alpha\left(n\alpha-2\mathcal{F}\right)\beta_{\nu}}C_{1}^{\pm}$. Then, $f^{\pm}$ takes the form

[TABLE]

where $C_{1}^{\pm}$ is fixed by the final normalization of the function.

It is important to note that $\epsilon$, $\bar{\gamma}$ and $\bar{\delta}$ depend on the energy and therefore on $n$, i.e., $\epsilon=\epsilon_{n}$, $\bar{\gamma}=\bar{\gamma}_{n}$ and $\bar{\delta}=\bar{\delta}_{n}$. Moreover, if we define $\theta$ as the argument of the Laguerre polynomials, we have

[TABLE]

Therefore, the corresponding eigenfunctions for the energies in Equation (80) are given by

[TABLE]

where $\psi^{\pm}_{n}$ turn out to be:

[TABLE]

with the functions $w^{(i)}_{n}$ given by

[TABLE]

[TABLE]

being $M_{n}$ written as

[TABLE]

and $\mathcal{N}_{n}$ denotes the normalization constant, explicitly given by

[TABLE]

with $I_{n}$ being the integral defined by

[TABLE]

3.3.2 Discussion

For this case, it is necessary to impose an additional condition to guarantee the square-integrability of the states described in Equation (99), which is given by

[TABLE]

Thus, in addition to having a finite number of states due to the condition in Equation (81), we obtain a well-defined domain in $k$ for each energy level. Consequently, the resulting spectrum is discrete, finite, dispersive, and bounded. However, although the bounds remain linear in $k$, they do not exhibit the same slope on both sides. Furthermore, the domain tends to decrease as the number of excited states increases (see Figure 3a).

On the other hand, when the electric field approaches its critical value, the right-hand side of the inequality in Equation (106) is affected, thereby modifying the regions where the existence of excited states is allowed. Thus, when plotting the energy as a function of $\mathcal{E}_{0}$, the states tend to become degenerate again, approaching the value $-kv_{\mathrm{d}}$ (see Figures 3b). Finally, as in previous cases, both the probability densities and probability currents depend on the band and valley indices, with the latter showing a nonzero value in the ground state (see Figures 3c and 3d).

4 Conclusions

In this work, we analyzed the dynamics of Dirac materials subjected to position-dependent electric and magnetic fields by decoupling the corresponding system of differential equations. The resulting equations, although not adopting the standard eigenvalue form, exhibit deep connections with Heun equations. The quantization conditions derived from finite polynomial solutions of these equations allowed us to obtain discrete energy spectra for each field configuration considered.

Across the different scenarios studied—from constant and exponential profiles to hyperbolic field configurations—we observed a consistent set of behaviors and phenomena. In all cases, the energy spectrum features a term with linear dependence on the wave number $k$ in the $y$-direction. For the ground state, the slope of this linear term is universally given by $-v_{\mathrm{d}}$, corresponding to the ratio between the amplitudes of the electric and magnetic fields. For the excited states, the behavior varies depending on the specific field configuration. In the first two cases, the slope of the linear component remains constant across all energy levels. However, in hyperbolic configurations, the slope depends on the quantum number $n$, and only the ground state retains a slope equal to $-v_{\mathrm{d}}$. Therefore, we can conclude that the presence of the electric field generates dispersive energy levels, characterized by a linear dependence on $k$.

The presence of a critical electric field strength, defined by $\mathcal{E}_{0}=\mathcal{B}_{0}(v_{y}-\nu v_{\rm t})$, plays a crucial role. As the electric field approaches this critical value, all energy levels collapse to the unique value $E=-kv_{\mathrm{d}}$, indicating a degeneracy point consistently observed in each configuration (see Figures 1b, 2b, and 3b). In the exponential and hyperbolic cases, it is necessary to impose additional conditions to ensure the square-integrability of the wavefunctions. These constraints limit the number of allowed bound states, resulting in an energy spectrum that is finite, discrete, and bounded. Furthermore, the domain in $k$-space where these stationary states exist shrinks as the number of excited states increases.

Both the probability density and the current density exhibit explicit dependence on the valley and band indices. In particular, due to the tilt of the Dirac cones, the probability current remains nonzero even in the ground state—behavior that distinguishes these systems from ideal graphene. Importantly, the method developed in this work enables a direct solution to the eigenvalue problem without resorting to supersymmetric quantum mechanics, Lorentz or similarity (non-unitary) transformations, thereby offering a more transparent alternative to previously established approaches[19, 5, 17, 16, 22, 18].

Despite the methodological differences, our results for the first two cases are fully consistent with previous findings obtained. Furthermore, all cases successfully reproduce the known spectrum of graphene under magnetic fields in the appropriate limit, thus validating our theoretical framework and confirming its broader applicability.

Finally, we emphasize that although polynomial solutions to Heun equations are known, they are typically expressed as power series whose coefficients satisfy the so-called three-term recurrence relation [24, 23]. However, in the present work, we have not only found such polynomial solutions but also managed to express them in terms of well-known families of orthogonal polynomials.

Acknowledgments The authors acknowledge financial support from CONAHCYT Project FORDECYT-PRONACES/61533/2020. EDB also acknowledges the SIP-IPN research grant 20254000. DOC gratefully acknowledges the financial support from SECIHTI through the postdoctoral fellowship, CVU number 712322.

Data Availability Statement All data generated or analyzed during this study are included in this published article.