Guided modes and terahertz transitions for two-dimensional Dirac fermions in a smooth double-well potential

The double-well in quantum mechanics has been studied in relation to various physical phenomena, ranging from vibrations of polyatomic molecules Rosen and Morse (1932), through to applications in Bose–Einstein condensation Milburn et al. (1997). Solutions to the Schrödinger equation for smooth double-wells are equally broad, and have been analyzed using perturbative methods Muller-Kirsten (2012), instanton calculus Gildener and Patrascioiu (1977), the WKB approximation Landau and Lifshitz (2013) and other techniques. With the rise of two-dimensional (2D) Dirac materials Wehling et al. (2014), comes a fresh opportunity to revisit the double-well problem in a relativistic setting, and to conduct ultra-relativistic experiments without the need for powerful accelerators. Indeed, there has been significant progress in creating guiding potentials in Dirac materials Huard et al. (2007); Özyilmaz et al. (2007); Gorbachev et al. (2008); Liu et al. (2008); Williams et al. (2011); Rickhaus et al. (2015), most recently using carbon nanotubes as top gates Cheng et al. (2019). These top-gated structures allow the potential profile to be controlled by manipulating the top-gate voltage, allowing the creation of well-defined smooth double wells.

Several approaches have been considered to achieve the goal of confinement in Dirac materials. These range from utilizing magnetic fields De Martino et al. (2007); Masir et al. (2009); Roy et al. (2012); Downing and Portnoi (2016a, b), to implementing Fermi velocity engineering Peres (2009); Raoux et al. (2010); Concha and Tešanović (2010); Downing and Portnoi (2017), through to introducing a spatially-dependant mass term Downing and Portnoi (2019), and most commonly, using electrostatic potentials Pereira Jr et al. (2006); Cheianov et al. (2007); Silvestrov and Efetov (2007); Tudorovskiy and Chaplik (2007); Shytov et al. (2008); Beenakker et al. (2009); Zhang et al. (2009); Matulis and Peeters (2008); Bardarson et al. (2009); Zhao and Yelin (2010); He et al. (2014); Recher and Trauzettel (2010); Hartmann et al. (2010); Downing et al. (2011); Stone et al. (2012); Hartmann and Portnoi (2014); Downing and Portnoi (2014); Hasegawa (2014); Xu and Ang (2015, 2016); Downing et al. (2015); Lee et al. (2016); Hartmann and Portnoi (2017a); Bai et al. (2018). Confinement in massless Dirac materials is notoriously difficult due to the Klein tunneling effect Klein (1929); Katsnelson et al. (2006). However, total confinement can be achieved in such systems for zero-energy states within an electrostatically defined waveguide, whose potential vanishes at infinity. This is because the density of states vanishes outside of the waveguide Hartmann et al. (2010); Hartmann and Portnoi (2014); Downing and Portnoi (2014); Hartmann and Portnoi (2017a). For non-zero-energy states, the bound states within the potential can couple to the continuum states outside the channel, and are therefore poorly guided. However, for massive particles, bound states can occur in spite of the Klein phenomena.

The alternative geometry of transmission through potential barriers has also been a subject of extensive research, with the majority of studies utilizing “sharp but smooth potentials”, i.e., potentials which are step-like or have kinks but are assumed to be smooth on the scale of the lattice constant, so that the effects of inter-valley scattering are neglected. Supercritical transmission Adler (1971); Dong et al. (1998); Dombey and Calogeracos (1999); Dombey et al. (2000); Kennedy (2002); Villalba and Greiner (2003); Kennedy et al. (2004); Guo et al. (2009); Miserev (2016) and tunneling through barriers has been studied for a variety of one-dimensional (1D) model potentials in both massless and massive 2D Dirac systems  Roy and Khan (1993); Katsnelson et al. (2006); Cheianov and Fal’ko (2006); Jiang et al. (2006); Bai and Zhang (2007); Tudorovskiy and Chaplik (2007); Cheianov et al. (2007); Barbier et al. (2008); Nguyen et al. (2009); Barbier et al. (2009); Villalba and González-Árraga (2010); Pereira Jr et al. (2010); Xu et al. (2010); Miserev and Entin (2012); Moldovan et al. (2012); Reijnders et al. (2013); Fallahazad et al. (2015); Hsieh et al. (2018); Avishai and Band (2020); Brun et al. (2020) including square barrier structures such as the double barrier Pereira Jr et al. (2007, 2008); Wei-Yin et al. (2013), inverted double well Bahlouli et al. (2012); Alhaidari et al. (2012), asymmetric waveguides He et al. (2014); Xu and Ang (2016) and various other step-like structures Alhaidari et al. (2012); Hasegawa (2014); Xu and Ang (2015). A variety of approaches ranging from the transfer matrix method through to the WKB method have been used. However, there is a dearth of studies concerning potentials which span both positive and negative energies, i.e., contain both electron-like and hole-like guided modes. The exceptions are periodic potentials Brey and Fertig (2009) and sinusoidal multiple-quantum-well systems Xu et al. (2010). Multiple-quantum-well systems are shown to exhibit transmission gaps in the electrons and holes spectra at tilted angles of incidence Xu et al. (2010); Pham and Nguyen (2015), and the number of oscillations in the transmission window depends on the number of quantum wells. For a double well studied in this paper (see Fig. 1 for the Gedankenexperiment sketch), the width of the transmission gap simply corresponds to the energy difference between guided modes, whereas the oscillations are associated with the splitting of a guided mode into multiple modes due to tunneling between wells. In this work, we consider the regime where the potential results in an “inversion” of electron and hole states, i.e., the valence band states in the barrier are higher than the conduction band states in the two wells. In this regime, the repulsion of electron and hole states gives rise to interesting features in the eigenvalue spectrum, namely avoided crossings, which can be controlled via the applied potential. These level avoided crossings, also provide a clear physical picture behind the anisotropic energy gap opening near the Dirac point of graphene subjected to a periodic potential Brey and Fertig (2009); Wang et al. (2014).

The 2D Dirac equation, has been the subject of considerable interest due to the explosion of research in transition metal dichalcogenides (TMDs) Xiao et al. (2012), Weyl semimetals Young et al. (2012), topological insulators Hasan and Kane (2010); Qi and Zhang (2011); Zazunov et al. (2016), low-dimensional forms of carbon Neto et al. (2009); Hartmann et al. (2014) and their silicon analogs Liu et al. (2011), for which their low-energy excitations can be described by a Hamiltonian of the form

$$\hat{H}=\hbar v\left(\hat{k}_{x}\sigma_{x}+s_{\mathrm{K}}\hat{k}_{y}\sigma_{y}+k_{z}\sigma_{z}\right),$$

(1)

where $\hat{k}_{x}=-i\frac{\partial}{\partial x}$, $\hat{k}_{y}=-i\frac{\partial}{\partial y}$, $\sigma_{x,\,y,\,z}$ are the Pauli spin matrices, $v$ is the Fermi velocity which plays the role of the speed of light, $s_{\mathrm{K}}=\pm 1$ is the valley index number, and $k_{z}$ is proportional to the particle’s in-plane effective mass. For graphene, $v=v_{\mathrm{F}}\approx 10^{6}$ m/s, and $k_{z}=0$ Wallace (1947). For quasi-1D forms such as narrow-gap carbon nanotubes and certain types of graphene nanoribbons Dutta and Pati (2010); Chung et al. (2016), the operator $\hat{k}_{y}$ can be substituted by the number $k_{y}=E_{g}/(2\hbar v_{\mathrm{F}})$, where $E_{g}$ is the value of the bandgap, which can be manipulated via the application of external fields Portnoi et al. (2008); Hartmann (2011); Hartmann et al. (2011); Hartmann and Portnoi (2015); Hartmann et al. (2014, 2019). Examples of massive 2D Dirac systems include, but not limited to silicene, germanene, TMDs, and graphene on top of lattice-matched boron nitride Giovannetti et al. (2007).

The gapless spectrum of graphene, and the nearly gapless band structure of narrow-gap carbon nanotubes and ribbons caused a natural attraction to their optical properties in the terahertz (THz) spectral range, which have led to a menagerie of promising applications in the field of THz optoelectronics Hartmann et al. (2014). The main goal of this paper is to demonstrate that the gate-induced double-well geometry allows for tuneable THz transitions between the various guided modes. We approach this problem by calculating the dispersion of guided modes in several model potentials which are smooth at the atomic scale thus allowing us to disregard the problem of valley mixing caused by jumps and kinks in piecewise potentials. These models describe well the shape of an electrostatic potential created in plane by three differently charged nanowires placed above a metallic gate, as can be shown by the mirror charge method.

The detailed calculations of the energy spectrum for a smooth double well, with the middle barrier exceeding the side barriers are provided in Section II, whereas the model potential with a middle barrier below the side barriers is treated in Appendix A. The selection rules for dipole transitions between the guided modes in the potential described in Section II are analyzed in Section III followed by the summary of the results in Section IV.

In what follows, we shall consider a Dirac particle described by the Hamiltonian given by Eq. (1) subject to a double well one-dimensional potential $U\left(x\right)$. Hereafter, the valley index number, $s_{\mathrm{K}}$, is set to one. The other valley’s wave function can be readily obtained by performing a sign change on $k_{y}$. The Hamiltonian acts on the two-component Dirac wavefunction

$$\Psi=\left(\begin{array}[]{c}\psi_{A}\left(x\right)\\ \psi_{B}\left(x\right)\end{array}\right)e^{ik_{y}y}$$

(2)

to yield the coupled first-order differential equations

$$\left(V-E+\Delta_{z}\right)\psi_{A}-i\left(\frac{d}{d\tilde{x}}+\Delta_{y}\right)\psi_{B}=0$$

(3)

and

$$\left(V-E-\Delta_{z}\right)\psi_{B}-i\left(\frac{d}{d\tilde{x}}-\Delta_{y}\right)\psi_{A}=0,$$

(4)

where $\tilde{x}=x/L$ and $L$ is a constant with the dimension of length. $V=UL/\hbar v_{\mathrm{F}}$ and the charge carrier energy, $\varepsilon$, have been scaled such that $E=\varepsilon L/\hbar v_{\mathrm{F}}$. The charge carriers propagate along the y-direction with wave vector $k_{y}=\Delta_{y}/L$, which is measured relative to the Dirac point and $\Delta_{z}=k_{z}L$ represents their effective mass in dimensionless units. Finally, $\Psi_{A}\left(x\right)$ and $\Psi_{B}\left(x\right)$ are the wavefunctions associated with the $A$ and $B$ sub-lattices respectively. Substituting

$$\psi_{A}=\frac{1}{2}\left[\psi_{1}\exp\left(\frac{1}{2}i\phi\right)+\psi_{2}\exp\left(-\frac{1}{2}i\phi\right)\right]$$

and

$$\psi_{B}=\frac{1}{2}\left[\psi_{1}\exp\left(\frac{1}{2}i\phi\right)-\psi_{2}\exp\left(-\frac{1}{2}i\phi\right)\right],$$

where $\phi=\arctan\left(\Delta_{y}/\Delta_{z}\right)$, allows Eqs.(3,4) to be reduced to a single second-order differential equation in $\psi_{j}$:

$$-\frac{d^{2}\psi_{j}}{dZ^{2}}+\mathrm{V}_{s_{j}}\psi_{j}=M^{2}\psi_{j},$$

(5)

where

$$\mathrm{V}_{s_{j}}\left(Z\right)=W^{2}\left(Z\right)-s_{c}s_{j}\frac{dW\left(Z\right)}{dZ}$$

(6)

and the other spinor component is found by the relation

$$\psi_{j^{\prime}}=-\frac{1}{M}\left(V-E+s_{c}s_{j}\frac{d}{dZ}\right)\psi_{j},$$

(7)

where $\tilde{x}=s_{c}iZ$, $s_{c}=\pm 1,$ $W=V-E$, $s_{j}=\left(-1\right)^{j}$, and $j=1,\,2$ ($j^{\prime}=2,1$) corresponds to the spinor components $\psi_{1}$ ($\psi_{2}$) and $\psi_{2}$ ($\psi_{1}$) respectively. In this basis the particle’s transverse momentum $\Delta_{y}$ and in-plane effective mass $\Delta_{z}$, have been combined into a single effective mass, $M=\sqrt{\Delta_{y}^{2}+\Delta_{z}^{2}}$. Eq. (5) is of same form as the Schrodinger equation, i.e., a second order differential equation with no first order derivative. Therefore, if $\mathrm{V}_{s_{j}}$ is equal to a potential possessing known solutions to the Schrödinger equation, then for zero energy states, we can readily write down the bound state spectrum of the potential, $W$, which satisfies the non-linear Eq. (6) Dutt et al. (1988); Cooper et al. (1988); Ho and Roy (2014). It should also be noted that for zero energy modes, the left-hand side of Eq. (5) is of the form of the super-Hamiltonian, $W$ plays the role of the superpotential and the allowed $M$ plays the role of the energy eigenstates Dutt et al. (1988); Cooper et al. (1988). Thus there exists a plethora of potentials which admit solutions for zero-energy modes. For non-zero energy, the $W$ which satisfies the non-linear Eq. (6) for a given potential $\mathrm{V}_{s_{j}}$ will be eigenvalues of an energy-dependent potential Schulze-Halberg and Roy (2017). However, we are interested in potentials which are independent of energy and are suitable for the use of modelling double-wells in Dirac materials.

The solutions to many second order differential equations, such as the Schrödinger equation and the 2D Dirac equation, can be obtained in terms of hypergeometric functions Landau and Lifshitz (2013); Cook (1971); Moshinsky and Szczepaniak (1989); Kennedy (2002); Guo et al. (2009); Candemir and Bayrak (2013); Ishkhanyan (2018). Indeed any second order differential equation, possessing three regular singularities, can be re-expressed as Euler’s Hypergeometric differential equation. Terminating the resulting hypergeometric series, and/or utilizing their well-known connection formulae allows the eigenvalues of many potentials to be obtained. However, when considering the 2D Dirac equation for double-well potentials, the resulting second order differential equation’s may contain more than three regular singularities. For example, four regular singularities means the differential equation can be re-expressed as Heun’s Equation Heun (1888). For the family of potentials which result in a differential equation having two regular singularities and one irregular singularity of rank 1, the second order differential equation can be transformed into the confluent form of Heun’s equation. Indeed, reducing a system of coupled first-order differential equations to the confluent Heun equation has been exploited to solve various generalisations of the quantum Rabi model Koç et al. (2002); Xie et al. (2017). Despite the absence of a general connection formula, the confluent Heun functions can still serve as a powerful tool in studying confinement potentials, and has been extensively applied in the fields of General Relativity and Quantum Gravity Hortaçsu (2018).

In what follows, we consider the case where the 2D Dirac equation reduces to the confluent Heun equation for a family of potentials, some of which can be used to describe a double-well separated by a barrier. This quantum model is quasi-exactly solvable Turbiner (1988); Ushveridze (2017); Bender and Boettcher (1998); Koç et al. (2002); Downing (2013); Hartmann and Portnoi (2014); Hartmann (2014); Hartmann and Portnoi (2017a), i.e. only a subset of the eigenvalues can be found explicitly. We study bound states contained within double-well potentials and calculate their entire energy spectrum. Energy level splitting associated with quantum tunneling between wells, and avoided crossings associated with the inter-mixing of electron-hole states are discussed. Finally, terahertz (THz) applications utilizing the doublet states and avoided crossing points are considered.

Let us search for transformations which may be performed on the dependent and independent variables of Eq. (5), and the corresponding energy-independent potentials which allow Eq. (6) to be transformed into the confluent Heun equation. In some instances, the resulting confluent Heun series can be terminated Ronveaux and Arscott (1995), allowing a subset of the eigenvalues to be obtained exactly. In other instances, the entire energy spectrum can be obtained fully via the Wronskian method Zhong et al. (2013); Maciejewski et al. (2014); Hartmann (2014); Hartmann and Portnoi (2014); Xie et al. (2017); Hartmann and Portnoi (2017a). Similar approaches have been implemented in the non-relativistic regime, indeed, for the Schrödinger equation there exists 35 choices for the coordinate transformation, each leading to eleven independent potentials which are exactly solvable in terms of the general Heun functions Ishkhanyan et al. (2015); Ishkhanyan (2018). The Schrödinger equation also reduces to various forms of confluent Heun equations for potentials such as the Natanzon family Ishkhanyan (2016) and several others Ishkhanyan and Krainov (2016).

The B$\mathrm{\hat{o}}$cher symmetrical form of the Confluent Heun equation, also known as the generalized spheriodal equation, can be written using the notation from Ref. Ronveaux and Arscott (1995)

$$\frac{d}{d\xi}\left[\left(\xi^{2}-1\right)\frac{dy\left(\xi\right)}{d\xi}\right]+\left[-p^{2}\left(\xi^{2}-1\right)+2p\beta\xi-\lambda-\frac{m^{2}+s^{2}+2ms\xi}{\xi^{2}-1}\right]y\left(\xi\right)=0,$$

(8)

where the regular singular points are located at $\xi=1$ and $\xi=-1$. The first order derivative can be removed by transforming the independent variable to $v\left(\xi\right)=y\left(\xi\right)\sqrt{1-\xi^{2}}$, this allows Eq. (8) to be written in normal symmetrical form:

$$-\frac{d^{2}v\left(\xi\right)}{d\xi^{2}}+\left[p^{2}+\frac{\lambda-2p\beta\xi}{\xi^{2}-1}+\frac{m^{2}+s^{2}-1+2ms\xi}{\left(\xi^{2}-1\right)^{2}}\right]v\left(\xi\right)=0.$$

(9)

We shall now re-express Eq. (5), into the same form as Eq. (9). There are many transformations which may be performed on the dependent and independent variables of Eq. (5), which give rise to a second order differential equation possessing no first order derivative. Applying the transformation of the independent variable $\xi=\xi\left(Z\right)$ and $\psi_{j}=\exp\left[-\int\frac{1}{2}\left(\frac{1}{\Phi}\frac{d\Phi}{d\xi}\right)d\xi\right]\chi_{j}$ to the dependent variable yields

$$-\frac{d^{2}\chi_{j}}{d\xi^{2}}+\left[\frac{1}{4}\left(\frac{1}{\Phi}\frac{d\Phi}{d\xi}\right)^{2}+\frac{1}{2}\frac{d}{d\xi}\left(\frac{1}{\Phi}\frac{d\Phi}{d\xi}\right)-s_{c}s_{j}\frac{1}{\Phi}\frac{dV}{d\xi}+\frac{\left(V-E\right)^{2}-M^{2}}{\Phi^{2}}\right]\chi_{j}=0$$

(10)

where $\Phi=\frac{d\xi\left(Z\right)}{dZ}$. The potential

$$V=\frac{a_{2}\xi^{2}+a_{1}\xi+a_{0}}{1-\xi^{2}}\Phi_{\eta_{1},\,\eta_{2}}$$

(11)

represents a family of energy-independent potentials which for the quasi-one-dimensional Dirac problem admits wavefunctions in terms of the confluent Heun functions. Here $\Phi_{\eta_{1},\,\eta_{2}}=C_{\eta_{1},\,\eta_{2}}\left(1-\xi\right)^{\eta_{1}}\left(1+\xi\right)^{\eta_{2}}$, $\eta_{1}$ and $\eta_{2}$ can take the values of 1 and 0, and $C_{\eta_{1},\,\eta_{2}}$ is a constant. As mentioned in the introduction, this family of potentials belongs to the class of quantum models which are quasi-exactly solvable Turbiner (1988); Ushveridze (2017); Bender and Boettcher (1998); Koç et al. (2002); Downing (2013); Hartmann and Portnoi (2014); Hartmann (2014); Hartmann and Portnoi (2017a), where only some of the eigenvalues can be expressed explicitly.

When $\eta_{1}=\eta_{2}=1$, Eq. (11) can be expressed as the Rosen-Morse potential, $V=b_{1}\tanh^{2}\left(\tilde{x}\right)+b_{2}\tanh\left(\tilde{x}\right)$, whose eigenvalue spectrum, which is quasi-exact, was discussed in ref. Hartmann and Portnoi (2017a). The analytic bound state energy spectra of the Rosen-Morse Potential for the two-dimensional Dirac problem has also been obtained via the Nikiforov-Uvarov method Ikhdair (2010). When $\eta_{1}\neq\eta_{2}$, Eq. (11) can be expressed as a generalized Hulthen-like potential:

$$V=\frac{c_{1}\exp\left(-2\lambda Z\right)+c_{2}\exp\left(-\lambda Z\right)+c_{3}}{1+c_{4}\exp\left(-\lambda Z\right)},$$

(12)

where $\lambda$ and $c_{1,2,3,4}$ are consants. When $c_{1}=0$ and $c_{4}=-1$, Eq. (12) becomes a linear combination of the Hulthen potential Hulthen (1942) and its logarithmic derivative. The Hulthen potential has previously been investigated for the 2-D Dirac equation using an algebraic approach Roy and Roychoudhury (1990). The case where Eq. (12) reduces to a single exponential has also been studied in graphene waveguides  Stone et al. (2012). It should also be noted that when $\Phi=i$, the potential given by Eq. (11) can be reduced to the shifted 1D Coulomb potential, which has applications to charged impurities and excitons in carbon nanotubes Downing and Portnoi (2014).

Let us return to the more general form of the potential given by Eq. (11), one of the opportunities to describe a smooth double well is to choose the case of $\eta_{1}=\eta_{2}=0$. We require that our potential is non-singular and vanishes as $x\rightarrow\pm\infty$. Therefore we set, $a_{2}=0$ and the potential becomes

$$V=\frac{a_{1}Z+a_{0}}{1-Z^{2}}\equiv\frac{A_{1}\tilde{x}+A_{0}}{1+\tilde{x}^{2}},$$

(13)

where the potential parameters $A_{1}$ and $A_{0}$ are related to the parameters appearing in Eq. (9) via the relations:

$$p=s_{p}\sqrt{E^{2}-M^{2}},\qquad\beta=-is_{c}A_{1}E/p,$$

$$\lambda=-A_{1}\left(A_{1}+is_{j}\right)+2EA_{0},$$

$$\Lambda_{1,\,2}=s_{a}A_{0},\qquad\Lambda_{2,\,1}=is_{a}\left(A_{1}+is_{j}\right)s_{c},$$

where $\Lambda_{1}=m$, $\Lambda_{2}=s$ and $s_{a,\,p}=\pm 1$. Equation (13) is a linear combination of the Lorentzian and its logarithic derivative, which has known solutions for the Radial Dirac equation Downing et al. (2011).

Since the potential given by Eq. (13) is smooth and vanishes as $x\rightarrow\pm\infty$, it may be utilized for modelling top-gated structures in 2D Dirac materials. The Lorentzian has already been used to model the potential generated by a top-gate formed by a carbon nanotube Cheng et al. (2019). It should also be noted that when $p=0$, i.e. $\left|E\right|=M$, the Lorentzian admits exact energy eigenvalues (see Appendix B). Exponentially decaying potentials also play an important role in the modelling of heterostructure-devices based upon zero-energy modes in Dirac-materials, as these potentials are often quasi-exact, admitting some exact energy eigenvalues. As mentioned previously in the introduction, the majority of work on multiple-barrier structures focus on smooth but sharp potentials. However, realistic potential profiles vary slowly over the length scale of the Dirac material’s lattice constant, and discontinuous potentials have yet to be realized. Furthermore, many piecewise potentials do not result in smooth wavefunctions across the whole of configuration space due to the nontrivial nature of their boundary conditions Xu and Ang (2015, 2016); You et al. (2017). Smooth potentials do not suffer from this problem  Hartmann et al. (2010); Hartmann and Portnoi (2014, 2017a, 2017b), additionally they permit inter-valley scattering to be neglected.

The stationary points, $\tilde{x}_{s_{r}}$, of the potential Eq. (13) are located at

$$\tilde{x}_{s_{r}}=-r+s_{r}\sqrt{1+r^{2}},$$

(14)

where $r=A_{0}/A_{1}$ and $s_{r}=\pm 1$. Making the coordinate transformation $\tilde{x}=\left|x^{\prime}\right|+\tilde{x}_{s_{r}}$, allows the potential Eq. (13) to be utilized as a model for quasi-one-dimensional double-wells in realistic top-gated Dirac material heterostructures. When $A_{1}<0$ and $s_{r}=-1$ the potential is a double well (see Fig. 2), whose local minima are located a distance of $2\sqrt{1+r^{2}}$ apart. The two wells, of depth $A_{1}/(2\tilde{x}_{1})$ are separated by a barrier, of height $A_{1}/(2\tilde{x}_{-1})$, which acts as a single-well for hole-like particles, with the possible experimental set-up shown in Fig. (1). Although we are dealing with a single-particle Hamiltonian, it is convenient to call the states with the energy growing with increasing $|\Delta_{y}|$ when $|\Delta_{y}|\rightarrow\infty$ as electron-like states or electrons, whereas the hole-like states correspond to $E\rightarrow-\infty$ when $|\Delta_{y}|\rightarrow\infty$. In the limit that $r\gg 1$, the central barrier’s height becomes negligible in comparison to the depth of the two wells. Conversely, when $A_{1}>0$ and $s_{r}=-1$, the potential becomes a single electron well, lying between two barriers, which act as a double well for hole-like particles.

The non-symmetrical canonical form of the confluent Heun equation can be written as

$$\frac{\partial^{2}u\left(\tilde{\xi}\right)}{\partial\tilde{\xi}^{2}}+\left(4p+\frac{\gamma}{\tilde{\xi}}+\frac{\delta}{\tilde{\xi}-1}\right)\frac{\partial u\left(\tilde{\xi}\right)}{\partial\tilde{\xi}}+\frac{4p\alpha\tilde{\xi}-\sigma}{\tilde{\xi}\left(\tilde{\xi}-1\right)}u\left(\tilde{\xi}\right)=0$$

(15)

where

$$\tilde{\xi}=\frac{1-\xi}{2},\qquad\gamma=m+s+1,\qquad\delta=m-s+1,$$

$$\alpha=-\beta+m+1,\qquad\sigma=\lambda-2p\left(\beta-\gamma\right)-m\left(m+1\right)$$

and $u=u\left(p,\,\alpha,\,\gamma,\,\delta,\,\sigma;\,\frac{1-\xi}{2}\right)$ are the Heun Confluent functions Ronveaux and Arscott (1995), which are related to the solutions of $v\left(\xi\right)$ of Eq. (9) by the relation:

$$v\left(\xi\right)=\left(1-\xi\right)^{\frac{\gamma}{2}}\left(1+\xi\right)^{\frac{\delta}{2}}e^{-p\xi}u\left(\frac{1-\xi}{2}\right),$$

(16)

which from hereon we shall denote as $v\left(\xi\right)=\Psi_{s_{j},s_{c}}$. Currently, there are no universal expressions for arbitrary parameters connecting Heun functions about different singular points, and unlike for Gauss hypergeometric functions, there are no expressions relating the derivative of a confluent Heun function to another confluent Heun function, although particular instances have been obtained Fiziev (2009); Shahnazaryan et al. (2012); Hartmann and Portnoi (2017a). The dearth of connection formulae makes obtaining the analytic expressions for the complete eigenvalue spectrum of a wavefunction expressed in terms of confluent Heun function non-trival. However, for the potential, Eq. (13), the derivative can be expressed exactly in terms of Heun Confluent function with new parameters. Using the definition of the Heun Confluent function Ronveaux and Arscott (1995), $\Psi_{s_{j},s_{c}}$ can be expressed in terms $\Psi_{-s_{j},s_{c}}$ via the relation:

$$-\frac{1}{M}\left[\left(V-E\right)\Psi_{s_{j},s_{c}}+s_{j}i\frac{d\Psi_{s_{j},s_{c}}}{dx}\right]=\left(-\frac{s_{a}+A_{0}+is_{c}A_{1}}{2M}\right)^{-s_{j}s_{c}s_{a}}\Psi_{-s_{j},s_{c}},$$

(17)

which allows the solution to Eq. (6) for $\psi_{1}$ and its corresponding $\psi_{2}$ to be be written as

$$\psi_{1}=\sum_{s_{c},c_{p},s_{a}}C_{s_{c},c_{p},s_{a}}\left(-\frac{s_{a}+A_{0}+is_{c}A_{1}}{2M}\right)^{-\frac{s_{c}s_{a}}{2}}\Psi_{-1,s_{c}}$$

(18)

and

$$\psi_{2}=\sum_{s_{c},c_{p},s_{a}}C_{s_{c},c_{p},s_{a}}\left(-\frac{s_{a}+A_{0}+is_{c}A_{1}}{2M}\right)^{\frac{s_{c}s_{a}}{2}}\Psi_{1,s_{c}}$$

(19)

respectively, where $m=s_{a}A_{0}$, $s=is_{a}\left(A_{1}+is_{j}\right)s_{c}$, and $C_{s_{c},c_{p},s_{a}}$ are weighting coefficients. It should be noted that if $m$ and $s$ are exchanged, then the phase factor appearing in Eq. (17) and Eq. (19) must be multiplied by the factor $4^{s_{c}s_{a}}$. For bound states, we require that the term $e^{-p\xi}$ decays as $x\rightarrow\infty$, this imposes two conditions: First $s_{p}=s_{c}$, and second, the absolute value of the particle’s energy must be less than the reduced mass, i.e., $\left|E\right|<\left|M\right|$.

The functions $\psi_{1}$ and $\psi_{2}$ are neither even nor odd. For understanding optical selection rules and for a more clear mapping of our results to the conventional double well picture we shall move to the symmterized basis functions $\left|\psi_{I}\right\rangle=\left|\psi_{1}\right\rangle-\left|\psi_{2}\right\rangle$ and $\left|\psi_{II}\right\rangle=-i\left(\left|\psi_{1}\right\rangle+\left|\psi_{2}\right\rangle\right)$. It should be noted that in this basis, an exchange of the sign of $V$ and $E$ is formally equivalent to exchanging $\psi_{I}$ with $\psi_{II}$ for the original $V$ and $E$. Therefore, inverting the potential only results in a sign change of the eigenvalues. In the symmetrized basis one can construct a linear combination of functions Eq. (18) and Eq. (19) to form solutions which are electron- ($s_{a}=1$) and hole-like ($s_{a}=-1$):

$$\psi_{I}=\sum_{s_{a}}D_{s_{a}}\Im\left(\rho^{\star}\Psi_{1,1}+\rho^{-1}\Psi_{1,-1}\right)$$

(20)

$$\psi_{II}=\sum_{s_{a}}D_{s_{a}}\Re\left(\rho^{\star}\Psi_{1,1}+\rho^{-1}\Psi_{1,-1}\right)$$

(21)

where $\Re$ and $\Im$ are the standard notation for the real and imaginary parts of a complex number, $\rho=\left(-2M\right)^{-\frac{s_{a}}{2}}\left(s_{a}+A_{0}-iA_{1}\right)^{\frac{s_{a}}{2}}$ and $D_{s_{a}}$ are the weighting constants. For bound states we require that the spinor components, Eq. (20) and Eq. (21), vanish at infinity, i.e., $\psi_{I,II}(x^{\prime}\longrightarrow\pm\infty)=0$. While the coordinate transformation, $\tilde{x}=\left|x^{\prime}\right|+\tilde{x}_{s_{r}}$, imposes the continuity condition

$$\psi_{I}\left(\tilde{x}_{-}\right)=0,$$

(22)

for odd-states and

$$\left.\frac{\partial\psi_{I}}{\partial x}\right|_{\tilde{x}_{-}}=0,$$

(23)

for even-states. The coordinate transformation also requires that $\psi_{I,II}\left(x^{\prime}\right)=-\psi_{I,II}\left(-x^{\prime}\right)$ and $\psi_{I,II}\left(x^{\prime}\right)=\psi_{I,II}\left(-x^{\prime}\right)$ for odd and even states respectively. However, the radius of convergence, of the power series, is unity. Therefore, an iterative, analytic continuation method must be employed to evaluate the confluent Heun function beyond its radius of convergence Motygin (2018). First the power series, $u$ is evaluated at the point $\widetilde{\xi}_{1}$, which lies within the radius of convergence. A new series expansion is then performed about the point $\widetilde{\xi}_{1}$, which in turn is used to evaluate the point $\widetilde{\xi}_{2}$, and so on and so forth. This allows the energy eigenvalues to be found via a simple shooting method Press et al. (1992); Killingbeck (1987) which utilizes the wave-functions of Eq. (18) and Eq. (19).

In Fig. (3) we plot the energy spectrum of bound states (square-integrable along the x-direction) defined by the potential parameters $A_{1}=-6$ and $A_{0}=0$, as a function of effective mass, $M$, displaying only the four lowest lying electron-like doublets and all the hole-like states within the energy range. When $E\ll M+V_{min}$, (where $V_{min}$ is the minimum value of the potential) the spectrum resembles a single, hole-like particle, trapped within a positive potential barrier. The highest in energy hole-like branch is $s$-like (nodeless) for $\psi_{I}$, and the number of nodes increases by one as the branches progress towards more negative energies. For $E\gg-M+V_{max}$, (where $V_{max}$ is the maximum value of the potential), the spectrum is electron-like, and the energy splittings of the doublets (each shown by a red line with a blue line underneath in Fig. (3)) tend to zero with increasing $M$, as expected for non-relativistic particles when increasing mass suppresses tunneling. In this regime, it can be seen from Fig. (4) that the particles are localized in the region of the wells, and that the wavefunctions behave as a linear combination of the wavefunctions associated with the individual wells. This gives rise to two states where there would of been one if tunneling was forbidden. For $\psi_{I}$, the number of nodes increases by one as the branches progress to higher energies. However, since $\psi_{II}$ is related to the derivative of $\psi_{I}$ and the potential is even in $x$, i.e. $V(-x)=V(x)$, the asymmetric linear combination of the two single well functions forces a node in the barrier between them, while the symmetric combination does not. Which is reflected in the behavior of $\psi_{II}$ being precisely the same as a wavefuntion for a non-relativistic particle in a double well, whereas, the component $\psi_{I}$ has the opposite parity. Finally, in the energy zone $M-V_{min}<E<-M+V_{max}$, the solutions are a linear combination of electron-like and hole-like states. It can be seen from see Fig. (3) that a series of avoided crossings are opened in the energy spectrum due to the repulsion of barrier and well-doublet states. Also in this energy zone, the energy level splitting associated with quantum tunneling increases, until the bound states merge with the continuum states. This energy splitting between the two doublets can be controlled in realistic Dirac-material-heterostructures by either, adjusting the strength of the applied top-gate voltages, or, by the choice of geometry. This, coupled with the ability to change the position of the Fermi level via a back gate, gives rise to many possible device applications. Indeed, utilizing the dependence on the number of zero-energy modes on potential strength in single-well smooth electron waveguides, has been proposed as the basis switching devices in graphene Hartmann et al. (2010).

Within this section we provide the general formalism for calculating the dipole matrix element of THz transitions between guided modes of quasi-particles described by the Hamiltonian given by Eq. (1). For a Dirac particle subject to a 1D potential $U(x)$, the unperturbed Hamiltonian is given by $\hat{H}_{0}=\hat{H}+\mathrm{\bm{I}}U_{z}$, and the corresponding eigenfunctions of the unperturbed Hamiltonian are given by $\left(\psi_{A}\left(x\right),\psi_{B}\left(x\right)\right)^{T}e^{ik_{y}y}/\sqrt{N}$, where $N$ is a normalization factor given by the expression, $N=l\intop_{-\infty}^{\infty}\left(\left|\psi_{A}\right|^{2}+\left|\psi_{B}\right|^{2}\right)dx$, and $l$ is the sample length. In the presence of an electromagnetic field, the particle momentum operator, $\hat{\bm{p}}$, is modified such that $\hat{\bm{p}}\rightarrow\hat{\bm{p}}+e\bm{A}/c$, where $e$ is the elementary charge, and $\bm{A}$ is the magnetic vector potential, which is related to $\hat{\bm{e}}=\left(e_{x},e_{y}\right)$, the unit vector describing the polarization of the electromagnetic wave, via the relation $\bm{A}=A\hat{\bm{e}}$. For linearly polarized light, the polarization vector is expressed as $\left(\cos\left(\varphi_{0}\right),\sin\left(\varphi_{0}\right)\right)$, while for right- and left-handed polarized light it is $\left(1,-i\right)/\sqrt{2}$ and $\left(1,i\right)/\sqrt{2}$, respectively. The general form of the perturbation due to an electromagnetic wave impinging normally to a Dirac material is

$$\delta H=\frac{eAv_{\mathrm{F}}}{c}\left(\sigma_{x}e_{x}+s_{\mathrm{K}}\sigma_{y}e_{y}\right).$$

(24)

Using the wavefunctions given in Eq. (2), in the limit that $l\longrightarrow\infty$ the matrix element of transition becomes:

$$\left|\left\langle f\left|\delta H\right|i\right\rangle\right|=G_{1}\left|\intop_{-\infty}^{\infty}\left[\left(e_{x}-is_{\mathrm{K}}e_{y}\right)\psi_{A,f}^{\star}\psi_{B,i}+\left(e_{x}+is_{\mathrm{K}}e_{y}\right)\psi_{B,f}^{\star}\psi_{A,i}\right]dx\right|\delta_{k_{y,i},k_{y,f}},$$

(25)

where $G_{1}=eAv_{\mathrm{F}}/\left(c\sqrt{N_{i}N_{f}}\right)$, and the indices $i$ and $f$ correspond to the initial and final states respectively. For free two-dimensional massless Dirac fermions, Eq. (25) becomes valley independent Saroka et al. (2018). While in contrast, massive two-dimensional Dirac fermions have valley-dependent optical transition rules Hartmann et al. (2019). We shall now analyze the optical selection rules between guided modes contained within the electrostatically-controlled double well defined by Eq. (13). To do so it is more convenient to move from the original $\left|\psi_{A}\right\rangle$, $\left|\psi_{B}\right\rangle$ basis to the symmetrized one. The functions $\psi_{A}$ and $\psi_{B}$ can be expressed in terms $\psi_{I}$ and $\psi_{II}$ via the transformation $\left(\psi_{A},\psi_{B}\right)^{T}=U\left(\psi_{I},\psi_{II}\right)^{T}$, where

$$U=\frac{1}{2}\left(\begin{array}[]{cc}i\sin\left(\frac{1}{2}\phi\right)&i\cos\left(\frac{1}{2}\phi\right)\\ \cos\left(\frac{1}{2}\phi\right)&-\sin\left(\frac{1}{2}\phi\right)\end{array}\right).$$

(26)

To shift to the symmterized basis we make a unitary transformation to $\delta H$ with the unitary operator $U$. Under this change, the perturbation $\delta H$ transforms as $\delta\breve{H}=U^{\dagger}\delta HU$:

$$\delta\breve{H}=-\frac{ev_{\mathrm{F}}A}{4c}\left(e_{x}\sigma_{y}+\frac{s_{\mathrm{K}}e_{y}\Delta_{z}}{\sqrt{\Delta_{z}^{2}+\Delta_{y}^{2}}}\sigma_{x}+\frac{e_{y}\Delta_{y}}{\sqrt{\Delta_{z}^{2}+\Delta_{y}^{2}}}\sigma_{z}\right)$$

(27)

The optical transitions between different guided modes can be categorized into two distinct groups. First, when the parity of each of the symmetrized spinor components $\psi_{I}$ and $\psi_{II}$, is preserved during the transition, i.e., both $\psi_{I,i}$ and $\psi_{I,f}$ are even, or both odd. Second, when the parity of each component changes after the transition, i.e., $\psi_{I,i}$ is odd (even), while $\psi_{I,f}$ is even (odd). When a transition, which preserves the parity of each spinor component, occurs, its matrix element given by Eq. (25), can be expressed in terms of the symmetrized spinor components, Eq. (20) and Eq. (21) as

$$\left|\left\langle f\left|\delta H\right|i\right\rangle\right|=G_{2}\left|\intop_{-\infty}^{\infty}\,\psi_{f}^{\dagger}\left(\frac{e_{y}\Delta_{y}\sigma_{z}}{\sqrt{\Delta_{y}^{2}+\Delta_{z}^{2}}}\right)\psi_{i}\,dx\right|\delta_{k_{y,i},k_{y,f}},$$

(28)

where $\psi_{i}=\left(\psi_{I,i},\psi_{II,i}\right)$, $\psi_{f}=\left(\psi_{I,f},\psi_{II,f}\right)$, $G_{2}=G_{1}/4$, and the indices $i$ and $f$ correspond to the initial and final states respectively. When the transition occurs between states of differing parity, the matrix element of transition, Eq. (25), becomes

$$\left|\left\langle f\left|\delta H\right|i\right\rangle\right|=G_{2}\left|\intop_{-\infty}^{\infty}\,\psi_{f}^{\dagger}\left(e_{x}\sigma_{y}+\frac{s_{\mathrm{K}}e_{y}\Delta_{z}\sigma_{x}}{\sqrt{\Delta_{y}^{2}+\Delta_{z}^{2}}}\right)\psi_{i}\,dx\right|\delta_{k_{y,i},k_{y,f}}.$$

(29)

We will now analyze the optical selection rules in two distinct regimes. The first regime corresponds to the case where the dispersion lines are linear, here the electron-like branches have positive dispersion, whereas the hole-like branches have negative dispersion. The second regime corresponds to the case where there is a significant reconstruction of the eigenfunctions, i.e., a strong admixture of electron-like and hole-like eigenfunctions (see Eq. (20) and Eq. (21).

In the linear dispersion regime, one can clearly see from Fig. (4) that the electron-like states are highly localised in the well regions, whereas the hole-like states are highly localised in the barrier region. The overlap between electron-like and hole-like states becomes increasingly small with increasing $\Delta_{y}$, which leads to all transitions between branches of positive and negative dispersion being heavily suppressed. Similar to non-relativistic quantum wells, in our double well system, within the linear regime, transitions between electron-like (hole-like branches) branches are allowed for light-polarized normal to the direction of the waveguide, while for light-polarized along the direction of the waveguide these transitions rapidly become vanishingly small with increasing $\Delta_{y}$.

For a massless Dirac system, in the energy range far away from an avoided crossing, but in the region where the splitting between doublet states is maximal, transitions between the two guided modes comprising the doublet are only allowed for light linearly polarized normal to the direction of the waveguide. In stark contrast, in the vicinity of the avoided crossings, transitions strongly occur between the avoided crossing states for light linearly-polarized along the direction of the waveguide. It should also be noted that transitions are also allowed between an avoided crossing level and the state belonging to an “opposite” parity branch (opposite parity in the large effective mass limit) which lies in between, in this instance transitions are predominately polarized normally to the waveguide.

The energy-level splitting associated with quantum tunneling between wells and the avoided crossings associated with electron-hole repulsion can be utilized for terahertz (THz) applications. For the appropriate choice of parameters, the potential given by Eq. (13) can be used to create doublet states within the vicinity of graphene’s charge neutrality point with an energy-level splitting in the THz range. Similarly, the energy gap associated with avoided crossings can be engineered to fall within the THz regime. These avoided crossings would absorb in a narrow frequency range due to the presence of the van Hove singularity at the pseudo-gap edge. A detailed analysis of the optical selection rules shall be a topic of future study, and is beyond the scope of this paper. Our simple analysis serves to demonstrate the potentiality of double-well smooth electron waveguides as the basis of polarisation sensitive THz detectors.

We have found a class of 1D double-well potentials, for which the 2D Dirac equation can be solved in terms of confluent Heun functions, and calculated the corresponding energy spectra. The energy-level splitting associated with electron tunneling between wells as well the electron-hole avoided crossing gap, can be controlled via the parameters of the electrostatic potential. Dipole optical transitions between the two states comprising a doublet, as well as transitions across the avoided crossing gap are both allowed, but follow drastically different polarization selection rules. For the doublet transitions, light is strongly absorbed for linear polarizations oriented normal to the waveguide. For the avoided crossing transitions, absorption of radiation polarized along the direction of the waveguide dominates. The presence of the van Hove singularity at the bottom of the pseudo-gaps opened in the spectrum leads to a narrow absorption peak, which can be tuned via the applied top-gate voltages to occur in the THz range.

This work was supported by the EU H2020 RISE project TERASSE (H2020-823878). RRH acknowledges financial support from URCO (71 F U 3TAY18-3TAY19). The work of MEP was supported by the Ministry of Science and Higher Education of Russian Federation, Goszadanie no. 2019-1246.