Magnetic Kronig-Penney model for Dirac electrons in single-layer graphene

The properties of Dirac electrons in a magnetic superlattice (SL) on graphene consisting of very high and thin (delta-function) barriers are investigated. We obtain the energy spectrum analytically and study the transmission through a finite number of barriers. The results are contrasted with those for electrons described by the Schrodinger equation. In addition, a collimation of an incident beam of electrons is obtained along the direction perpendicular to that of the SL. We also highlight the analogy with optical media in which the refractive index varies in space.


Introduction
During the last five years single-layer graphene (a monolayer of carbon atoms) has become a very active field of research in nanophysics [1,2]. It is expected that this material will serve as a base for new electronic and opto-electric devices. The reason is that graphene's electronic properties are drastically different from those, say, of conventional semiconductors. Charge carriers in a wide single-layer graphene behave like "relativistic", chiral, and massless particles with a "light speed" equal to the Fermi velocity and possess a gapless, linear spectrum close to the K and K ′ points. One major consequence is the perfect transmission through arbitrarily high and wide barriers, referred to as Klein tunneling.
One of the most challenging tasks is to learn how to control the electron behavior using electric fields in graphene. This task is made complicated precisely by the Klein tunneling according to which Dirac electrons in graphene can tunnel through arbitrarily wide and high electric barriers [3].
Alternatively, one can apply a magnetic field to control the electron motion. It was shown in numerous papers that an inhomogeneous magnetic field can confine standard electrons described by the Schrödinger equation [4,5,6]. The question then arises whether it can confine Dirac electrons in graphene. Up to now semiinfinite magnetic structures, that are homogeneous in one direction, were considered and made the task simpler by converting the problem into an one-dimensional (1D) one [7,8,9,10,11,12,13,14,15]. In particular, a magnetic confinement of Dirac electrons in graphene was reported in structures involving one [7] or several magnetic barriers [8,9] as well as in superlattices, without magnetic field for some very special values of the parameters involved [16]. In such structures standard electrons can remain close to the interface and move along so-called snake orbits [5] or in pure quantum mechanical unidirectional states [4].
Given the importance of graphene, it would be appropriate to study this magnetic confinement more systematically. We make such a study here by considering a magnetic Kronig-Penney (KP) model in graphene, i.e., a series of magnetic δ-function barriers that alternate in sign. This model can be realized experimentally in two different ways: 1) One can deposit ferromagnetic strips on top of a graphene layer but in a way that there is no electrical contact between graphene and these strips. When one magnetizes the strips along the x direction, cf. Fig. 1(a), by, e.g., applying an in-plane magnetic field, the charge carriers in the graphene layer feel an inhomogeneous magnetic field profile. This profile can be well approximated [22] by 2B 0 z 0 h/d(x 2 + z 2 0 ) on one edge of the strip and by −2B 0 z 0 h/(x 2 + z 2 0 ) on the other, where z 0 is the distance between the 2DEG and the strip, and d and h the width and height of the strip (see Fig. 1(b)). The resulting magnetic field profile will be modeled by two magnetic δ functions of height 2πB 0 h. Such ferromagnetic strips were deposited on top of a two-dimensional electron gas (2DEG) in a semiconductor heterostructure in Ref. [23].
2) It was recently shown that local strain in graphene induces an effective inhomogeneous magnetic field [24] (Fig. 1(c)). When one puts the graphene layer on a periodically structured substrate the graphene at the edges of the substrate becomes strained and the situation can be described by a magnetic δ-function profile such as that shown in Fig. 2.
In a quantum mechanical treatment of the above two systems the vector potential A(x) is the essential quantity and, within the Landau gauge, A(x) is nothing else than a periodic array of step functions. The Hamiltonian describing this system is periodic and consequently we expect the energy spectrum of the charge carriers in graphene to exhibit a band structure. The advantage of this magnetic Kronig-Penney (KP) model is mainly its analytical simplicity that provides some insight and allows a contrast with the same model for standard electrons [27]. To do that we adapt a method developed in optics, for a media with periodic in space refractive index. This optical method is clear and very well suited to the problem. Incidentally, there are many analogues of optical behavior in electronics, such as focusing [17,18,19], collimation or quasi-1D motion of electrons and photons [16,20,4,8], and interference [21] in a 2DEG.
The manuscript is organized as follows. In Sec. 2 we present the method and evaluate the spectrum and electron transmission through two antiparallel, δ-function magnetic barriers. In Sec. 3 we consider superlattices of such barriers and present numerical results for the energy spectrum. In Sec. 4 we consider a series of δ-function vector potentials and our concluding remarks are given in Sec. 5.

Characteristic matrix for Dirac electrons
An electron in a single-layer graphene, in the presence of a perpendicular magnetic field B(x), which depends on x, is adequately described by the Hamiltonian where p is the momentum operator, v F the Fermi velocity, and A(x) the vector potential.
To simplify the notation we introduce the dimensionless units: Here ℓ B is the magnetic length and t the tunneling strength. In these units Eq.
(1) takes the form Then the equation HΨ(x, y) = EΨ(x, y) admits solutions of the form with ψ I (x, y), ψ II (x, y) obeying the coupled equations Due to the translational invariance along the y direction we assume solutions of the form Ψ(x, y) = exp ik y y(U(x), V (x)) T , with the superscript T denoting the transpose of the row vector. For B(x) ∼ δ(x) the corresponding vector potential is a step function A(x) ∼ Θ(x). For A(x) = P constant, Eqs. (4) and (5) take the form Equations (3)-(7) correspond to those for an electromagnetic wave propagating through a medium in which the refractive index varies periodically. The two components of Ψ(x, y) correspond to those of the electric (or magnetic) field of the wave [28,29]. Equations (6) and Eq. (7) can be readily decoupled by substitution. The result is where Z = U, V . If E 2 → E ′ and (k y + P ) 2 → V ef f , Eq. (8) reduces to a Schrödinger equation for a standard electron where V ef f (k y , x) = (k y + P ) 2 can be considered as an effective potential. Taking θ 0 as the angle of incidence, we have k x = E cos θ 0 = [E 2 − k 2 y ] 1/2 and k y = E sin θ 0 are the wave vector components outside the medium and k ′ x = E cos θ = [E 2 − (k y + P ) 2 ] 1/2 is the electron wave vector inside the medium and θ = tan −1 (k y /k ′ x ) is the refraction angle. This renders Eq. (8) simpler with acceptable solutions for U and V U(x) = A cos (Ex cos θ) + B sin (Ex cos θ), For future purposes, we write U and V as a linear combination of U 1 , U 2 and V 1 , V 2 : We now multiply the equations of the first row by U 2 and U 1 , respectively, and those of the second by V 2 and V 1 . The resulting equations lead to where D = detD and Equation (12) shows that the determinant of the matrix (13) associated with any two arbitrary solutions of Eq. (8) is a constant, i.e, D is an invariant of the system of Eqs. (11). This also follows from the well-known property of the Wronskian of secondorder differential equations. For our purposes the most convenient choice of particular solutions is such that Then the solution with U(0) = U 0 , V (0) = V 0 , can be expressed as or, in matrix notation, as Since D is constant, the determinant of the square matrix N is a constant; its value, found by setting It is usually more convenient to express U 0 and V 0 as a function of U(x) and V (x). Solving for U 0 and V 0 we obtain This matrix M is unimodular, |M| = 1. Now we can find the characteristic matrix from Eqs. (9) and (10) as

Bound states
Regards to the average of vector potential we shall consider two different systems: one with zero average and the other with non zero average along the x-direction. First let us consider the magnetic field profile as shown in Fig. 2(a) for which the corresponding vector potential is where Θ(x) = 0(x < 0), 1(x > 0) is the theta function. This vector potential has a non-zero average, and the corresponding effective potential becomes (see Fig. 3(a)) as Here L is measured in the unit of magnetic length l B . There are two different cases which we have to consider. Case 1 for k y < −P/2: as shown in Fig. 3(a) by the full red curve, we have a 1D symmetric quantum well which, as is well-known, has at least one bound state (see also Fig. 3(c)). For E 2 < k 2 y the particle will be bound while for E 2 > k 2 y we have scattered states, or equivalently the electron tunnels through the magnetic barriers.
Case 2 for k y > −P/2: as shown in Fig. 3(a) by the dotted blue curve, the effective potential is like a barrier. We have a pure tunneling problem. With reference to Fig. 2, x 1 = 0 and x 2 = L, the solutions are as follows. For x < 0 the wave function is where κ = E cosh ξ and k y = E sinh ξ, while for x > L it is In the middle region, 0 < x < L, the wave function is given by With Matching the wave functions at x = 0 and x = L leads to a system of four equations relating the coefficients C, D, F , and Q. Setting the determinant of these coefficients equal to zero, we obtain the transcendental equation, the solution of it gives the energy spectrum cos θ cosh ξ cos k ′ L + sin θ sinh ξ sin k ′ L = 0.
For the special value of k y = −P , sin θ = 0, and we can rewrite Eq. (25) as cos (EL) = 0 or equally E n = n + 1 2 π L . The resulting bound states, as a function of k y , are shown by the red full curves in Fig. 4(a). The area of existence of bound states is delimited by the lines E = −k y and E = −(k y + P ). The number of bound states increases with |k y | which is also clear from the behavior of the min. and max. of the effective potential (see Fig. 3(c)). No bound states are found for k y > −P/2 as is also apparent from Fig.  3(c). For k y → −P/2 the potential is shallow and only one bound state exists. The average velocity v n (k y ) along the y direction is given by Fig. 4(a) it is clear that these bound states move along the y-direction, i.e. along the magnetic barriers. Their velocity v y > −v F is negative for k y < −P but as the electron is approaching k y → −P we have v y → 0. For k y > −P the velocity v y > v F is positive. This can be understood from the maximum and minimum of the effective potential which is shown in Fig. 3(c). The energy bound states can only exist between these two lines. Notice that the slope of min V ef f is negative for k y < −P while it turns positive for k y > −P which explains the k y dependence of the velocity. From Fig. 4(a) it is clear there are two different classes of bound states. The bound state which follows very closely the k y = −P curve and extends to the region −P < k y < 0 with energy close to zero has a wavefunction that is concentrated around the position of the two magnetic delta function and decays exponentially in the region 0 < x < L. The wave function of the other bound states are concentrated in a region between the two magnetic delta-functions (i.e. like in a standing wave fashion) and decays exponentially outside this region. Next, we consider a structure with zero-average vector potential as shown in Fig.  2(b), with corresponding effective potential shown in Fig. 3(b). The effective potential for k y < −P/2 and k y > P/2, consist of a potential well and a potential barrier and therefore has at least one bound state. Thus we expect bound states for all k y with energy between E = −(+)(k y + P ) and E = −(+)k y when k y < −P/2 (k y > P/2). The dispersion relation for those bound states are the solution of where M is the transfer matrix for the unit shown in Fig. 2(d). These bound states are shown by the red full curves in Fig. 4(b). Because of the spatial inversion symmetry of the vector potential the spectrum has the symmetry E(−k y ) = E(k y ). Notice that for −P < k y < P the lowest bound state has energy E ≈ 0. For −P/2 < k y < P/2 we have two potential barriers and therefore no bound states.

Reflection and transmission coefficients
Consider a plane wave incident upon a system of two δ-function magnetic barriers, identical in height but opposite in direction, placed at x = 0 and x = L, as shown schematically in Fig. 2(a). In this case the vector potential is constant for 0 ≤ x ≤ L, zero outside this region, and homogeneous in the y direction. Below we derive expressions for the amplitudes and intensities of the reflected and transmitted waves. Let A, R, and T denote the amplitudes of the incident, reflected, and transmitted waves, respectively. Further, let θ 0 be the angle of incidence and exit as shown in Fig.  2(b). The boundary conditions give The four quantities U 0 , V 0 , U, and V given by Eqs. (28) are connected by the basic relation Q 0 = MQ; hence, with J = m ′ 11 + m ′ 12 e iθ 0 and K = m ′ 21 + m ′ 22 e iθ 0 , we have where m ′ ij are the elements of the characteristic matrix of the medium, evaluated at x = L. From Eq. (29) we obtain the reflection and transmission amplitudes In terms of r and t the reflectivity and transmissivity are The characteristic matrix for a homogeneous vector potential is given by Eq. (19). Labeling with subscripts 1, 2, and 3 quantities which refer to the regions, respectively, I, II, and III of Fig. 2(a), and by L = x 2 −x 1 distance between the magnetic δ-functions, we have (β = EL cos θ i ) The reflection and transmission amplitudes r and t are obtained by substituting these expressions in those for J and K that appear in Eq. (30). The resulting formula can be expressed in terms of the amplitudes r 12 , t 12 and r 23 , t 23 associated with the reflection at and transmission through the first and second "interface", respectively. We have and similar expressions for r 23 and t 23 . In terms of these expressions r and t become r = r 12 + r 23 e 2iβ 1 + r 12 r 23 e 2iβ , t = t 12 t 23 e iβ 1 + r 12 r 23 e 2iβ .
The amplitude t of the transmission through the system is given by [8,9,12,24], where k y = E sin θ 0 , and k y + P = E sin θ i . This equation remains invariant under the changes E → −E , θ 0 → −θ 0 , θ i → −θ i . A contour plot of the transmission is shown in Fig. 4(a) and slices for constant k y and E are shown respectively in Fig. 4(c) and Fig. 4(d). By imposing the condition that the wave number k x be real for incident and transmitted waves, we find that the angles θ 0 and θ i are related by Equation (36) expresses the angular confinement of the transmission elaborated in Refs. [8], [12], [26], and [9]. Notice its formal similarity with Snell's law. Using Eq. (36) we obtain the range of incidence angles θ 0 for which transmission through the first magnetic barrier is possible For the special value of the energy E = P/2 and θ 0 in the range −π/2 ≤ θ 0 ≤ π/2, we have θ i = π/2 while for E = −P/2 the result is θ i = −π/2. Alternatively, we can put θ i = ±π/2 in Eq. (36) and obtain, for P > 0, the result where the +(−) sign corresponds to E > 0 (E < 0). A contour plot of the transmission as function of E and k y , obtained from Eq. (35), is shown in Fig. 4(a). In Fig. 4(a) we distinguish three different regions. In the region between E = −(k y + P/2) and E = −k y , the wave vector of the incident wave is imaginary and they are evanescent waves. In this region k ′ is real and it is possible to find localized states. The k and k ′ for the second region between E = k y + P/2 and E = −(k y − P/2) are real and the electron can tunnel through the magnetic δ-barriers. In the blue shadow region between E = k y + P/2 and E = k y − P/2, k is real but k ′ is imaginary and solutions inside the barrier are evanescent and there is very little tunneling which becomes very quickly zero. The transmission probability |T | = t · t * is equal to 1 for cos (2β) = 1. In this case the energy becomes E n = ± n 2 π 2 /L 2 + (k y + P ) 2 1/2 , n = 1, 2, ....
The condition cos (2β) = 1, or equivalently β = nπ = Ehcosθ 2 with n an integer, should be combined with that for the transmission to occur in the region delimited by the curves E = ±(k y + P ) and E = ±k y . For example, in Fig. 4(a) for k y = −0.05 and 0 < E < 0.2 we have 12 maxima. It is readily seen that with these parameters in Eq.
(39) we find 12 different energies as shown in Fig. 4(c). Fig. 4(b) shows a contour plot of the transmission for the structure shown in Fig. 2(d), which is symmetric around k y = 0. Notice that the number of resonances has increased substantially as compared to previous case which is due to the fact that we have twice as many magnetic barriers in our systems.

N units
We consider a system of N units, such as those shown in Fig. 2 with χ = 1 2 Tr M and u N the Chebyshev polynomials of the second kind: where Here ζ is the Bloch phase of the periodic system [33], which is related to the eigenfunctions of M. In the limit case of N → ∞, we have total reflection when ζ is outside the range (−1, 1).

Superlattice
Here we consider a finite number N of lattice unit shown in Fig. 2(c). We set β 2 = Eb cos θ 2 , β 1 = Ea cos θ 1 , p 2 = 1/ cos θ 2 , The characteristic matrix M 2 (L) for one period is readily obtained, in terms of these quantities, as in Sec. II, and from that the characteristic matrix M 2N (NL) of the multilayer system according to Eq. (41). Its elements are where s = p 2 p 1 , u N ≡ u N (χ), and The reflection and transmission coefficients of the multi-unit system are immediately obtained by substituting these expressions into Eq. (30). The numerical results are shown in Figs. 5, 6 for finite superlattice with N = 10 units. Two different type of structures are considered as shown in the insets to Figs. 6. The transmission doesn't have k y → −k y symmetry for the periodic system with magnetic delta up-down as is  apparent from Fig. 6(a). We contrast these results with the case in which we used an arrangement of magnetic delta-function as in previous structure plus another unit with opposite direction of magnetic delta function. As is clearly shown in Fig. 6(b), we  have k y → −k y symmetry for the transmission probability through this structure. The transmission resonances are more pronounced, i.e., the dips become deeper, when the number of barriers increases for both types of units. But the gaps occur when the wave is mostly reflected. The position of these gaps, which are especially pronounced as N increases, can also be found from the structure of the Bloch phase ζ, as shown in Figs. 5(c) and (d). In Fig. 5 the bound states are shown by the blue solid curves that are situated in the area −k y −P/2 < E < −k y +P/2 in case (a) and in −k y −P/2 < E < −k y for case (b) plus an area located symmetric with respect to k y . Notice in Fig. 5(a) that several bound states merge into a resonant states at E = −k y + P/2. This is different from Fig. 4(a) where each bound state becomes a resonant state at E = −k y .

Spectrum of a superlattice
Lets take N → ∞. We can find the energy-momentum relation from the previous standard calculation [32,33] by using where M is the characteristic matrix of one period, which results into With reference to the regions I and II shown in Fig. 2(a), we write k 1 = [E 2 − k 2 y ] 1/2 and k 2 = [E 2 − (k y + P ) 2 ] 1/2 and show the solution for E 2 > (k y + P ) 2 in Fig. 9. Differences of Eq. (49) from the corresponding result of Ref. [27], Eq. (7), for the case of Schrödinger electrons, is the term P 2 in the prefactor of the second term on the right-hand side and the linear E vs k spectrum instead of the quadratic one in Ref. [27]. If P is large the differences become more pronounced. Our numerical results for the energy spectrum are shown in Figs. 7, 8, 9, 10, and 11. The results for standard and Dirac electrons show similarities but also important differences. The first band shows a qualitative difference near k y ≈ 0, see Fig. 9. As Figs. 7, 9(a) and 9(b) show, the band behavior in the k = k x direction for fixed k y is constant and almost symmetric about k y = 0; the motion becomes nearly 1D for relatively large k y . From the contour plots of Figs. 9(b) and (d), as well as from Fig. 7(c), we infer a collimation along the k y direction, i.e., v y ∝ ∂E/∂k y ≈ v F and v x ≈ 0, which is similar to that found for a SL of electric potential barriers [16] for some specific values of the barrier heights. Also there are no gaps for k y ≈ 0 in Fig. 10(a) but there are for the case of Dirac electrons as seen in Fig. 10(b). This difference can be traced back to the presence of P 2 in the dispersion relation Eq. (49) when compare to the same equation for the standard electron. The even-number energy bands in Fig. 11(b) are wider than those in Fig. 11(a) and, as a function of the period, the energy decreases faster for Dirac electrons. This behavior of the bands for Dirac electrons is very similar to that for the frequency ω vs k y or L in media with a periodically varying refractive index [29]. This is clearly a consequence of the linear E − k relation. Notice the differences between the lowest bands shown in panels (a) and (b) in Fig. 12 and in particular the difference between the corresponding drift velocities as functions of k y .

A series of δ-function vector potentials
In the limit that the distance between the opposite directed magnetic barriers decrease to zero the vector potential approaches a δ-function [24]. We consider a series of magnetic δ-function vector potentials A(x) = ∞ n=−∞ A 0 δ(x − nL) as shown in Fig. 11(a). First we consider a single such potential, that is zero everywhere except at x = 0. We start with Eqs. (6) and (7) which becomes now The solutions are readily obtained in the form φ a = A cos (εx cos θ) + B sin (εx cos θ) , x < 0, C cos (εx cos θ) + D sin (εx cos θ) , x > 0, Integrating Eqs. (50), (51) around 0 gives and We now consider the entire series of δ-function vector potentials shown in Fig. 13(a) and use Eqs. (54) and the periodic boundary condition Ψ I (0) = e ikxL Ψ II (L). The resulting dispersion relation for the superlattice is where η = 1 + ℓ B A 0 . From Eq. (57) we can find the energy spectrum as (58) We can define (1/L) cos −1 2η(1 + η 2 ) −1 cos(k x L) = s, and obtain E = ±[2nπ + k 2 y + s 2 ] 1/2 .
The energy bands around the Dirac point are plotted in Fig. 12(b). Notice that: 1) there is an opening of a gap at the Dirac point, 2) the motion is strongly 1D, i.e. along the k y -direction, and 3) higher subbands have a smaller dispersion.

Concluding remarks
We developed a magnetic Kronig-Penney model for Dirac electrons in graphene. The model is essentially a series of very high and very narrow magnetic δ-function barriers alternating in signs. The treatment of the transmission through such a series of barriers followed closely the one developed in optics for media in which the refractive index varies in space [28,29]. We contrasted a few of the results with those for standard electrons described by the Schrödinger equation [27].
In several cases the energy spectrum or the dispersion relation were obtained analytically, cf. Eqs. (25), (39), (47), (48), and (57), largely due to the simplicity of the model and the adapted method from optics. For only two magnetic δ-function barriers, opposite in sign, we saw several bound states, whose number increases with |k y |, and a reduction of the wavevector range for which tunneling is possible, cf. Fig.  4(a). This is in line with that reported earlier for single [7] and multiple [9] barriers. The reduction becomes stronger as we increase the number of barriers, cf. Fig. 4(b). We also made contact with Snell's law in optics, cf. Eq. (36): the term P/E represents the deviation from this law.
An important feature of the superlattice results is a collimation of an incident electron beam normal to the superlattice direction at least for large wave vectors. As easily seen from Figs. 7 and 8, for |k y | ≥ 2 we have v x ∝ ∂E/∂k x ≈ 0 for the first three minibands in the middle panels and nearly five minibands in the right panels. This occurs for both standard electrons and Dirac electrons but notice an important difference for |k y | ≈ 0 shown clearly in Fig. 9. This collimation is similar to that reported in Ref. [16] for superlattices involving only electric barriers but with somewhat unrealistic large barrier heights.
It is also worth emphasizing the differences and similarities in the first two minibands and the corresponding drift velocities as functions of k y for different periods L and constant k x as shown in Fig. 12. Notice in particular the resemblance between the drift velocities in the lowest miniband for standard electrons and the second miniband for Dirac electrons.
Given that ferromagnetic strips were successfully deposited on top of a 2DEG in a semiconductor heterostructure [23], we hope they will be deposited on graphene too and that the results of this paper will be tested in a near future.