Mapping borophene onto graphene: Quasi-exact solutions for guiding potentials in tilted Dirac cones

We show that if the solutions to the (2+1)-dimensional massless Dirac equation for a given 1D potential are known, then they can be used to obtain the eigenvalues and eigenfunctions for the same potential, orientated at an arbitrary angle, in a tilted anisotropic 2D Dirac material. This simple set of transformations enables all the exact and quasi-exact solutions associated with 1D quantum wells in graphene to be applied to the confinement problem in tilted Dirac materials such as borophene. We also show that smooth electron waveguides in tilted Dirac materials can be used to manipulate the degree of valley polarization of quasiparticles travelling along a particular direction of the channel. We examine the particular case of the hyperbolic secant potential to model realistic top-gated structures for valleytronic applications.


Introduction
It can be shown using supersymmetric methods that whenever the Schrödinger equation can be solved exactly for a onedimensional (1D) potential, there exists a corresponding potential for which the two-dimensional (2D) Dirac equation admits exact eigenvalues and eigenfunctions 1 . A broad class of quasi-1D potentials can also be solved quasi-exactly by transforming the 2D Dirac equation into the Heun equation 2 or one of its confluent forms [3][4][5] , or via the application of Darboux transformations [6][7][8][9][10][11] . These exact and quasi-exact bound-state solutions have direct applications to electronic waveguides in 2D Dirac materials 2-5, 12, 13 , such as graphene, where the low-energy spectrum of the charge carriers can be described by a Dirac Hamiltonian 14 , and the guiding potential can be generated via a top gate [15][16][17][18][19][20] . Recent advances in device fabrication, utilizing carbon nanotubes as top gates, has enabled the detection of individual guided modes 21 , opening the door to several new classes of devices such as THz emitters 5,13 , transistors 12 , and ultrafast electronic switching devices 22 . These advances in electron waveguide fabrication technology make the need for analytic solutions all the more important, since they are highly useful in: determining device geometry, finding the threshold voltage required to observe a zero-energy mode 12 , calculating the size of the THz pseudo-gap in bipolar waveguides 13 , as well as ascertaining the optical selection rules 5,13 in graphene heterostructures.
In extension to the well-known case of graphene, Dirac cones can in general possess valley-dependent tilt 23 . There are only a handful of 2D electronic systems that have been predicted to host these tilted cones 24-33 , one of which is 8-Pmmn borophene 23, 34 , which has attracted considerable attention. In general, boron-based nanomaterials are a growing field of interest [35][36][37] ; indeed, exploring the tilt of 2D Dirac cones in the context of 8-Pmmn borophene has recently led to a plethora of theoretical works spanning many fields of research 38-84 from optics to transport, and many more. The spectacular rise of borophene, and the growing interest in tilted Dirac materials, has led to the revisiting of several well-known problems in graphene, e.g., Klein tunneling 47 and transport across quasi-1D heterostructures 43,52,73,85 , within the context of tilted Dirac materials. As mentioned previously, methods such as supersymmetry and reducing the Dirac equation to the Heun equation, utilize solutions of known problems, to generate solutions to new ones. This begs the question: does a simple mapping exist which would allow us to harness the large body of exact and quasi-exact solutions for 1D waveguides in graphene and then apply them to materials with tilted Dirac cones?
In what follows, we show that the differential equations governing guided modes in an anisotropic tilted Dirac waveguide (orientated at an arbitrary angle) can, via a simple transformation, be mapped onto the graphene problem, i.e., transformed into the massless 2D Dirac equation 14 , for the same potential, but of modified strength, effective momentum, and modified energy scale. After outlining the transformation, we study the particular case of the hyperbolic secant potential, which in graphene is known to admit quasi-exact solutions to the eigenvalue problem; but nevertheless, the whole spectrum can be obtained via a semi-analytic approach 2, 12 . We use this known graphene waveguide spectrum to generate the corresponding tilted waveguide Figure 1. A schematic diagram of an electrostatic potential, U(x ), created by an applied top-gate voltage in a tilted 2D Dirac material. The waveguide is orientated at an angle of θ , relative to the x-axis of the crystal. The potential is invariant along the y -axis, and varies in strength along the x -axis. The x − y axes are denoted by the solid black arrows, whereas the crystallographic axes x − y are shown by the light gray arrows. spectrum, which is verified using a transfer matrix method. Finally, we discuss valleytronic applications.

Transformation
The Hamiltonian describing the guided modes contained within a smooth electron waveguide in a tilted Dirac material can be written aŝ wherek x = −i∂ x ,k y = −i∂ y , σ x,y are the Pauli spin matrices, σ 0 is the identity matrix, v x and v y are the anisotropic velocities, v t is the tilt velocity and s = ±1 is the valley index number; here s = 1 and s = −1 are analogous to the K and K valley respectively. In what follows, we set s = 1, but it should be noted that the other valley's eigenvalues can be obtained by replacing v y and v t with −v y and −v t . In general, the crystallographic orientation is not known, nor is it currently possible to deposit the top gate at a selected angle relative to the crystallographic axis. Therefore, we shall solve for the case of a waveguide at an arbitrary angle relative to the crystallographic axis (x − y). The electrostatic potential, U (x, y), is 1D, directed along the y -axis, and varies along the x −axis (see Fig. 1 for geometry), i.e., U = U (x ). We rotate the x − y axes counterclockwise through an angle θ . The new axes x − y are defined by the original coordinates via the transformation: x = x cos θ + y sin θ , Hence, the wave vector operators in the non-rotated coordinate systemk = k x ,k y are expressed in the rotated coordinate frame,k = k x ,k y , via the relations: The Hamiltonian, Eq. (1), acts on the two-component Dirac wavefunction Ψ = (ψ A (x ) , ψ B (x )) e ik y y to yield the coupled first-order differential equationsĤΨ = εΨ, where ψ A and ψ B are the wavefunctions associated with the A and B sublattices of the tilted Dirac material. These coupled first-order differential equations can be recast into the same equations used to describe guided modes propagating along a smooth electron waveguide in graphene:

2/12
where the effective potential V , energy E, and momentum ∆ are obtained from the original tilted case via the relations: where V (x ) = U (x ) L/hv x , E = εL/hv x and ∆ = k y L, and L is a constant, associated with the effective width of the potential.
We define the tilt and anisotropy parameters as t = v t /v x and T = v y /v x , respectively, and l = 1 − (1 − T 2 ) sin 2 θ . The eigenfunctions of the guided modes in the effective graphene sheet, Φ, can be mapped onto the tilted Dirac spinor components, ψ A and ψ B , via the expression: where ϕ = arctan (T tan θ ) and µ = (l − t sin θ ) 1 2 (l + t sin θ ) − 1 2 . It then follows that if the eigenfunctions and eigenvalues are known for the potential V in graphene, one can immediately write down the eigenfunctions and eigenvalues of a 1D confining potential of the same form, orientated at an arbitrary angle, in a tilted Dirac material. Conversely, if a quasi-1D potential readily admits exact or quasi-exact solutions for the tilted case, and no solutions are known for the graphene problem, then our mapping method can be used to obtain the eigenfunctions and eigenvalues for the case of graphene. This mapping also reveals the angular dependence of the number of guided modes contained within the waveguide. Namely, it can be seen from Eq. (5) that rotating the orientation of the waveguide is equivalent to varying the effective depth of the potential (see Fig. 2b). Indeed, the effective potential's depth, V 0 , is equal to the actual potential's depth, V 0 at θ = 0 and rises to a maximum value of It should be noted that to perform the same transformations for the other chirality (i.e., the s = −1, or graphene K valley analog), one must exchange t with −t, and T with −T . Therefore, the eigenvalue spectrum of one valley can be obtained by a reflection of the other valley's eigenvalue spectrum about the k y -axis. It can be seen from Eq. (5) that in the absence of a tilt term, i.e., t = 0, the eigenvalue spectrum of a given valley is symmetric with respect to k y . Thus, both chiralities have the same band structure. Similarly for t = 0, if the waveguide is orientated such that cos θ = 0, the eigenvalue spectrum is chirality-independent. In all other cases, the energy spectrum for a given valley lacks E(k y ) = E(−k y ) symmetry. This gives rise to the possibility of utilizing smooth electron waveguides in tilted Dirac materials as the basis of valleytronic devices. This will be discussed in the penultimate section.

Quasi-exact solution to the tilted Dirac equation for the hyperbolic secant potential
In this section, we shall apply our simple transformations, given in Eq. (5), to generate the energy spectrum of a smooth electron waveguide in a tilted Dirac material for a potential which has been studied in depth in graphene: This potential, shown in Fig. 2a, belongs to the class of quantum models which are quasi-exactly solvable 2, 86-92 , where only some of the eigenfunctions and eigenvalues are found explicitly. The depth of the well is given by V 0 , and the potential width is characterized by the parameter L. Here V 0 and L are taken to be positive parameters. In graphene, the wavefunctions can be solved in terms of Heun polynomials, which reduce to hypergeometric functions for the case of zero energy 2, 12 . For zero-energy modes ( E = 0), the permissible values of ∆ are given by the simple relation ∆ = V 0 − n − 1 2 , where V 0 is the depth of the effective potential and n is a non-negative integer. For non-zero energies, exact energy eigenvalues can be obtained when the Heun polynomials are terminated 2 . To illustrate the power of our mapping method, we apply the transformations given in Eq. (5) to the well-known zero-energy solutions to the 2D Dirac equation for the 1D hyperbolic secant potential. The corresponding tilted Dirac equation solutions become: The relative strength of the effective potential's depth, V 0 , compared to the actual potential depth V 0 .
and their corresponding n = 0 eigenfunctions for two different waveguide orientations are shown in Fig. 3, for the case of borophene, i.e., v Fig. 4 we plot the borophene eigenvalue spectrum for the potential defined by V 0 = 3.2 for two orientations, θ = 0 and θ = π/2, as well as the graphene waveguide spectra (quasi-analytically determined 2 ) used in the mapping. In the same figure, we also plot the numerical solutions to the tilted Dirac problem obtained via a transfer matrix method (see supplementary material). We show in blue crosses the exact solutions given in Eq. (8) together with the complete set of mapped quasi-exact solutions given in ref. 2 . It can be seen from Fig. 4 that the waveguide orientated at θ = π/2 contains more bound states than the waveguide orientated at θ = 0. This is a result of the effective graphene potential being deeper, and thus supporting more guided modes. It should be noted that for potentials which vanish at infinity, i.e., V (±∞) = E = 0, only the zero-energy modes are truly confined, since the density of states vanishes outside of the well. Guided modes occurring at non-zero energies can always couple to continuum states outside of the well, thus having a finite lifetime.

Valleytronic applications
It has been suggested that the valley quantum number can be used as a basis for carrying information in graphene-based devices 93 in an analogous manner to spin in semiconductor spintronics 94,95 . Unlike in the case of graphene (in the conical regime), smooth electrostatic potentials in tilted Dirac materials can be utilized as a means to achieve valley polarization. The majority of studies have focused on tunneling across electrostatically-induced potential barriers, and valley filtering and beam splitting have been demonstrated 43, 47, 52, 73, 85 . It has also been shown that the allowed transmission angles through a potential can be controlled using magnetic barriers 85 . We propose a change in geometry: rather than studying chirality-dependent transmission across barriers, we shift the focus to studying guided modes along quasi-1D confining potentials. The conductance along such a channel can be measured by placing one terminal at each end. According to the Landauer formula, when the Fermi level is set to energy E (by modulating the back-gate voltage 21 ), the conductance along the waveguide is simply 2(n K + n K )e 2 /h, where n K and n K are the number of modes belonging to the s = 1 and s = −1 chirality, respectively (or the K and K valley in the effective graphene sheet), at that particular energy. In a 2D Dirac material subject to a quasi-1D potential, the introduction of the tilt parameter breaks the E(∆) = E(−∆) symmetry for a given valley. Indeed, for a given valley, the additional tilt term increases the particle velocity along one direction of the barrier and decreases it along the opposite direction; and vice versa for the other valley. For the case of type-I Dirac materials, i.e., t < T , it can be seen from Fig. 5 that for a given sign of ∆, the eigenvalues of the critical solutions (sometimes referred to as the zero-momentum solutions: i.e., bound states with energy |E| = |∆|) belonging to the two valleys are different. Thus, providing that t = 0 and cos θ = 0, there exists a range of energies for which there will be more bound states propagating along a particular direction belonging to one valley than the other, i.e., valley polarization. The degree of valley polarization can be controlled by varying the strength of the electrostatic potential and by changing the position of the Fermi level, which in practical devices is achieved by modulating the back-gate voltage 21 . Full valley polarization can be achieved for energies less than the lowest lying supercritical state (defined as a bound state with energy E = −∆) belonging to the valley where the tilt term enhances the particle's velocity, indicated by the shaded region in Fig. 5.   For type-III tilted Dirac materials, i.e., t = T , full valley polarization occurs for all energies and all orientations of the waveguide. For such materials the infinite number of positive-energy critical solutions of graphene map onto the zero-energy modes of a type-III Dirac material. Consequently, the infinitely many zero-energy modes will give rise to a sharp peak in the conductance along the channel when the Fermi energy is in the proximity of E = 0. This is in stark contrast to the case of graphene, where the creation of an infinite number of zero-energy states requires an infinitely deep and wide potential 4, 96-99 . Since the potential possesses infinitely many bound states, which infinitely accumulate at E = 0, a type-III Dirac material could be used as a THz emitter; namely, the Fermi level could be set below E = 0, and optical photons can be absorbed from low-lying energy levels to E = 0. Then the photo-excited carriers can relax back down to the Fermi level via the emission of THz photons through the closely spaced energy levels in the proximity of zero energy.
Lastly, for type-II tilted Dirac materials, i.e., t > T , full valley polarization occurs for all energies; however, bound states occur only for orientations within the range −1/ √ t 2 − T 2 < tan (θ + nπ) < 1/ √ t 2 − T 2 . In the limit at which θ becomes imaginary, i.e., the boundary at which the equi-energy surfaces become unbounded, the effective potential required for mapping diverges.

Conclusion
We have shown that if the eigenvalues and eigenfunctions of a quasi-1D potential in graphene are known, then they can be used to obtain the corresponding results for the same potential (with modified strength, orientated at arbitrary angle), for an anisotropic, tilted 2D Dirac material. Therefore, all the rich physics associated with guiding potentials in graphene, e.g., THz pseudo-gaps in bipolar waveguides, can be revisited in the context of tilted Dirac materials, but with the distinct advantage of knowing the eigenvalues and eigenfunctions. We have also shown that in stark contrast to smooth electron waveguides in graphene, valley degeneracy can be broken in tilted 2D Dirac materials for a broad range of waveguide orientations, anisotropy and tilt parameters. The degree of valley polarization along the waveguide can be controlled by varying the potential strength of the top gate, and also by changing the back-gate voltage. Tilted 2D Dirac materials, such as borophene, are therefore promising building blocks for tunable valleytronic devices.

Methods
The Supplementary Information contains a full description of the transfer matrix method used to calculate the band structure of guiding potentials in 2D Dirac materials with tilted Dirac cones.

Data Availability
This study did not generate any new data.

Numerical Model of a Waveguide in a 2D Dirac Material with Tilted Dirac Cones
In order to confirm the results of the main text, we discretize smooth waveguides and numerically calculate their band structure using the transfer matrix method. We begin by determining the wavefunction in a single 1D square well for a Dirac material with tilted cones, as an extension to the graphene case 1 . An arbitrary number of square wells are placed one after another to approximate the smooth waveguides discussed in the main text. For the system of sequential wells we set up a transfer matrix which links the wavefunction components in the leftmost and rightmost well. Bound states yielding the band structure are found numerically by placing boundary conditions on the transfer matrix. The boundary conditions ensure the wavefunction has plane-wave solutions inside the waveguide and decaying solutions at the edges of the system.

Solution to the Dirac Equation with Tilted Dirac Cones in a Square Well
The Dirac equation for electrons in a 2D semimetal with tilted Dirac cones subjected to the potential U(x, y) has the form v x stσ 0py + sT σ ypy + σ xpx Ψ(x, y) = ε −U(x, y) Ψ(x, y), (S1) with energy ε, momentump i where i = x, y, Fermi velocity v x , tilt parameter t, anisotropy parameter T and valley index s = ±1. The identity matrix is σ 0 , the Pauli matrices are σ x and σ y , and the spinor wavefunction is Ψ(x, y) = (φ A (x, y), φ B (x, y)) which is written in the basis of the Bloch sums of the Dirac material. We wish to find the wavefunction of the j th well where j = 1, 2, ...N. Each well has finite depth U j with arbitrary width along the x -axis (x j−1 ≤ x ≤ x j ) and is invariant along the perpendicular y -axis (see Fig. 1 of the main text). Next we perform a coordinate transformation on the momentum operatorŝ p x =p x cos θ −p y sin θ andp y =p x sin θ +p y cos θ before inserting their definitionp i = −ih∂ i , where i = x , y . Due to the invariance along the y -axis we can insert a plane-wave solution along this axis Ψ j (x , y ) = Ψ j (x )e ik y y . It is at this point that we introduce an arbitrary length scale L allowing us to define the dimensionless parameters ξ = x /L, V j = U j L/hv x , E = εL/hv x and ∆ = k y L. Let the wavefunction in the j th well be of the form Ψ j (ξ ) = (φ A j , φ B j ) e ib j ξ , where b j is in general a complex number. All preceding steps yield the eigenvalue problem where 0 is the null vector. We define the above eigenvalue problem as where for non-trivial solutions we must ensure |H j − σ 0 (E −V j )|= 0. This step yields two solutions (±) for the complex parameter where l = 1 − (1 − T 2 ) sin 2 θ . The first term of Eq. (S3) is always real whilst the second term can be either real or imaginary depending on the parameters. We shall now introduce the subscript ± corresponding to the two roots of b j,± . Inserting these solutions back into Eq. (S2) allows us to define φ B j, The wavefunction in the region ξ j−1 ≤ ξ ≤ ξ j , where ξ j = x j /L can be written as a linear combination of the forward and backward propagating plane-wave solutions: where C j is a normalization factor. This can be written as Ψ j (ξ ) = Ω j (ξ )(α j , β j ) where Ω j (ξ ) = e ib j,+ ξ e ib j,− ξ Λ j,+ e ib j,+ ξ Λ j,− e ib j,− ξ , where α j = C j φ A j,+ and β j = C j φ A j,− .

Transfer Matrix
We now consider a smooth waveguide approximated by N square wells as sketched in Fig. S1a. The analytics will be limited to a simple waveguide that terminates at the same maximum potential to the left and right, i.e., V 1 = V N = V max , but the analysis can be readily extended to the general case. We must first consider our conditions for a bound state. From the perspective of single wells, our bound state must have a plane-wave solution in the deepest well, V min , which means that in this region, Im(b j,± ) = 0. Additionally, we must ensure that the wavefunction decays in the leftmost and rightmost wells, V max , meaning Im(b 1,± ) = Im(b N,± ) = 0 so that it disappears as ξ → ±∞. These two criteria are satisfied by the condition 2/3 and f 2 = − |T | l 2 − t 2 sin 2 θ l 2 . (S9) For the reader's convenience these regions have been sketched for two waveguide orientations (see Fig. S1b and c). Now that we have established the requirements for bound states in the waveguide, we need to match the wavefunctions of each square well at their boundaries. If we begin on the right side of the waveguide, we must match the wavefunctions Ψ j (ξ ) in regions j = N − 1 and j = N at the interface ξ N−1 . This will allow us to write the wavefunction components of region N as a function of the components in region N − 1: (S10) Iterating this process will allow us to write the components of the wavefunctions in the right region ( j = N) as a function of the components of the left region ( j = 1) as (α 1 , β 1 ) = T(α N , β N ) , where the transfer matrix is defined as To ensure the wavefunction decays at ξ → ±∞ we set α 1 = 0 and β N = 0, hence This places a constraint on the transfer matrix T αα = 0. Therefore, bound states can be found as energies E and wavevectors ∆ that satisfy the criteria defined in Eq. (S7) and that ensure T αα = 0, which can be found numerically.