Time evolution of electron waves in graphene superlattices

Graphene is a carbon-based material where the atoms are arranged in a honeycomb lattice.^1–8 This genuinely two-dimensional material is characterized by unusual and remarkable electronic properties, including a “relativistic”-type spectrum, such that the propagation of low-energy electrons in graphene is described by the massless Dirac equation.³ These properties pushed graphene into the frontline of condensed matter physics research.^1–7 

Interestingly, it may be feasible to gain additional control over the transport properties of electrons in graphene by artificially introducing a new length scale into the system in the form of a periodic electrostatic potential.^9–16 These structures, known as graphene superlattices, may be realized using different techniques, such as with periodically patterned gates, with the deposition of adatoms on graphene’s surface, or using a crystalline substrate.^17–24 

Rooted in the superlattice concept, it was recently proposed that it may be possible to extend to electronics some phenomena and devices originally discussed in the context of electromagnetism,^25,26 such as the “perfect lens”^27,28 or an electron “wormhole”.²⁹ In these works, the propagation of the electrons in the superlattices was studied using an effective medium approach, wherein the granular details of the superlattices are homogenized.²⁷ Within this framework the structure is regarded as a continuum and the dynamics of the wave function envelope is described by an effective Hamiltonian.

The main objective of the present work is to demonstrate that the effective medium theory proposed in Ref. 27 can be used to determine the time evolution of electronic states in graphene superlattices. To this end, we develop a finite-difference time-domain (FDTD) algorithm to characterize the propagation of electron waves in superlattices using both microscopic and macroscopic formalisms. We present a detailed comparison of the physical response predicted by the two approaches. It is important to highlight that the application of numerical methods to the Schrödinger and Dirac equations is well known and has been widely reported in the open literature (e.g. Refs. 30–39). In particular, Refs. 30 and 31 study the propagation of electron waves in graphene heterojunctions using the FDTD method. However, the key novelty of our work is the verification that the effective medium formalism developed in Ref. 27 can be used to predict the physical response in the frame of time evolution problems as well as in the frame of stationary state problems in complex propagation scenarios. We are unaware of similar studies in related physical platforms. To this end, a FDTD algorithm based on a leapfrog update scheme⁴⁰ is developed and applied to the characterization of graphene superlattices using both the macroscopic and microscopic models. It is underlined that the theory developed in Refs. 30 and 31 cannot be directly applied to the propagation of electron waves in the context of the effective medium model considered here, and this is the reason why we develop our own numerical scheme to solve the modified Dirac equation. It is important to make clear that what we designate here by “microscopic model” corresponds to the description of the electron wave propagation in graphene using the two-dimensional Dirac equation, which is itself an effective medium theory. Our effective Hamiltonian corresponds thus to a second level of homogenization. The Dirac equation is a valid starting point for a microscopic theory when the period of the superlattice is much larger than the atomic constant of graphene. Furthermore, the Dirac equation has been the basis of several successful theoretical predictions supported by experiments, such as the negative refraction of Dirac fermions and the connection between the optical conductivity of graphene and the fine structure constant.^41,42 One of the nontrivial aspects of an effective medium approach is the correct formulation of boundary conditions at interfaces. Here, the boundary condition introduced in Ref. 29 is implemented in the FDTD code, and its validity is numerically confirmed in complex propagation scenarios.

The article is organized as follows. In Sec. II the effective medium model is succinctly reviewed. Section III describes the FDTD algorithm used to characterize the time evolution of electron waves. In Secs. IV and V, the theory is applied to the propagation of stationary states and to the time evolution of initial electronic states, respectively. Finally, in Sec. VI conclusions are drawn.

The propagation of the charge carriers in graphene may be described using the massless Dirac equation³:

being $H ˆ =−iħ v F σ x ∂ ∂ x + σ y ∂ ∂ y +V r$ the “microscopic” Hamiltonian operator near the K point, where V is the microscopic potential, $ψ= Ψ 1 , Ψ 2$ is the two-component pseudospinor, v_F ≈ 10⁶m/s is the Fermi velocity and σ_x, σ_y are the Pauli matrices. It is relevant to mention that for superlattices made from adatoms or for graphene-boron nitride superlattices the Hamiltonian gains an additional σ_z component.⁴³ Such a term can be easily incorporated into our FDTD discretization, but for simplicity in what follows we focus on superlattices described by a microscopic electrostatic potential with a one-dimensional spatial variation of the form $V x = V av + V osc sin 2 π x / a$ (see Fig. 1(a)). Here, V_av is the average potential and V_osc is the peak amplitude of the oscillating part of the potential.

By solving the Dirac equation (1) it is possible to completely characterize the wave function ψ in both spatial and time domains. However, since the microscopic potential V has a complex spatial dependence, this approach may be computationally demanding and provides limited physical insights.

A solution to reduce the complexity of the problem is to use effective medium methods. It was recently shown that the electronic states with the pseudo-momentum near the Dirac K point can be accurately modeled using an effective medium framework.^27,29 In this approach, the microscopic potential is homogenized and the effective Hamiltonian treats the superlattice as a continuum characterized by some effective parameters.^27,29 For an unbounded superlattice the effective Hamiltonian is of the form:

where $σ χ = σ x x ˆ +χ σ y y ˆ$ and χ is an effective medium parameter designated by anisotropy ratio. The value of χ depends on the amplitude of the fluctuating part of the microscopic potential V_osc and is numerically determined using the “first principles” homogenization approach described in Ref. 27. The explicit dependence of χ on the peak amplitude of the oscillating part of the potential V_osc is shown in Fig. 1(b), and reveals that the proper tuning of V_osc may enable a regime characterized by an extreme anisotropy, where χ = 0. It is well known that graphene superlattices may have strongly anisotropic Dirac cones and particle velocities, and may allow for the propagation of electron waves with virtually no diffraction.^18,19,29 The stationary states energy dispersion obtained with the macroscopic framework is given by²⁷:

where E is the energy of the electrons and $k= k x , k y$ is the wave vector associated with the electronic state, measured with respect to the K point. In pristine graphene the anisotropy ratio is unity (χ = 1) and the energy dispersion is determined by the usual Dirac cone, as depicted in Fig. 1(ci)). This linear energy dispersion corresponds to an isotropic propagation velocity. Quite differently, the energy dispersion for a superlattice with extreme anisotropy (χ = 0) corresponds to a Dirac cone stretched along the y-direction, as shown in Fig. 1(cii), so that the electrons are allowed to propagate only along the x-direction. It is important to note that for sufficiently strong modulations of the electric potential new Dirac points can emerge in the energy spectrum.^17,21–24 This effect is fully described by the microscopic theory [Eq. (1)] but not by the effective Hamiltonian which predicts a unique Dirac point.

The minimum value of V_osc that leads to an extreme anisotropy is V_osca/ħ v_F ≈ 7.55. Interestingly, this amplitude for V_osc is precisely coincident with an analytical solution of Ref. 12 for new zero-energy states (due to the zero-crossing of the energy spectrum for large k_y) and for a strong enhancement of the conductance in superlattices with a sinusoidal profile. Since in an extreme anisotropy regime there is a strong enhancement of the electric response,⁴⁴ it follows that the continuum approximation captures the essence of the findings of Ref. 12. Furthermore, it will be shown that despite the emergence of new Dirac points the effective Hamiltonian describes extremely well the electron wave propagation in graphene superlattices with a strong potential modulation when the electron state has a characteristic width a few times larger than the lattice period.

The knowledge of the time dynamics of electrons is essential to predict the response of graphene-based electronic devices. As a starting point, it is convenient to consider a generalized Schrödinger equation with a fictitious source term:

The term $j= j 1 j 2 T$ may be regarded as an external source that injects carriers into the system. As explained later with detail, this fictitious source is useful to characterize extended (non-localized) stationary energy states. In case of heterojunctions with a spatially varying anisotropy ratio χ, the macroscopic Hamiltonian (2) should be written as:

We replaced $χ σ y i ∂ ∂ y → 1 2 χ σ y i ∂ ∂ y + 1 2 ∂ ∂ y χ σ y i$ to ensure that the Hamiltonian remains Hermitian when the anisotropy ratio χ depends explicitly on the y coordinate. In addition, a spatially dependent χ requires that the macroscopic potential V_av is transformed as $V a v → V a v + ∂ u ∂ x ħ v F$ (Eq. (5b)), where $u=u χ$ is defined as in Eq. (A3) of  Appendix A. This transformation is required for the correct modeling of the wave propagation at an interface between distinct superlattices, as discussed in  Appendix A.

For conciseness, next we present a unified description of the FDTD method that applies to both the microscopic and macroscopic (effective Hamiltonian) approaches. It is implicit that in the microscopic approach χ = 1 and V_ef = V. The system (4) with the Hamiltonian (5) can be spelled out as:

To discretize this system and obtain the time update equations in an explicit form, the spatial domain is discretized into a rectangular grid, such that the node distances along the x- and y-directions are taken equal to Δ_x and Δ_y, as shown in Fig. 2. The pseudospinor is sampled at time instants separated by the time step Δ_t.

As usual, the partial derivatives are replaced by finite differences such that for a generic physical entity F:

where ∂_l ≡ ∂/∂l for l = x, y, t, and $F x , y , t =F p Δ x , q Δ y , n Δ t ≡F p , q , n$⁠. In our algorithm the components of the pseudospinor Ψ₁ and Ψ₂ are sampled at different points of space-time such that:

Therefore, the discretization of Ψ₁ and Ψ₂ ensures that the partial derivatives of Ψ₁ (Ψ₂) in the spatial domain are defined over the same spatial subgrid as Ψ₂ (Ψ₁). The time derivatives of Ψ₁ and Ψ₂ are also defined in staggered time grids. Applying these principles to the system (6) leads to the following time update equations:

where $V p , q = V e f p Δ x , q Δ y$⁠, $χ p , q =χ p Δ x , q Δ y$⁠, etc. In the derivation of the above equations, we use interpolation formulas such as $Ψ 1 p , q , n + 1 2 ≈ 1 2 Ψ 1 , p , q n + Ψ 1 , p , q n − 1$ to evaluate the wave function at points not lying in the pertinent subgrid. By sequentially using the explicit update equations (9) in a leapfrog scheme,⁴⁰ it is possible to determine $Ψ 1 , p , q n$ and $Ψ 2 , p + 1 2 , q + 1 2 n + 1 2$ at a generic instant of time n from the knowledge of the initial state of the system (⁠$Ψ 1 , p , q n$ and $Ψ 2 , p + 1 2 , q + 1 2 n + 1 2$ calculated at n = 0). This update scheme is completely analogous to the FDTD algorithm for electromagnetic waves.⁴⁰ In Appendix B, we formally demonstrate that our algorithm is stable provided the time step satisfies:

Moreover, in  Appendix C, it is shown that unbounded regions (with the graphene sheet infinitely extended) can be numerically emulated using a “perfectly matched layer” (PML). For complex heterostructures the FDTD algorithm may require substantial computational resources. Given the processing power provided by current graphics processing units (GPUs), the algorithm was implemented based on parallel computing methods. A validation of the FDTD algorithm for electron waves propagating in simple graphene heterostructures is reported in the Supplementary Materials.⁴⁵ In the next sections, the algorithm is applied to graphene superlattices.

In the following, we investigate the propagation of extended stationary states (with a definite energy) in graphene superlattices using both the microscopic and effective medium formalisms.

First, we consider the propagation problem in an unbounded graphene superlattice. Within the continuum approximation described in Sect. II, the wave packet (group) velocity is given by [see Eq. (3)]:

Thus, unlike in pristine graphene (χ = 1) in general the group velocity is not parallel to the quasi-momentum k. In particular, in the extreme anisotropy limit (χ = 0) the group velocity satisfies $v=± v F x ˆ$⁠, and hence in this case all the electron states flow along the x-direction and the superlattice supports diffractionless propagation.^11–16 To illustrate this effect, we consider the propagation of a Gaussian electron wave with initial beamwidth radius R_G = 14.14a and normalized energy E₀a/ħ v_F = 0.2 in an unbounded graphene superlattice. The superlattice is characterized by χ = 0, $V av − E 0 a/ħ v F =0.1$⁠, being a = 10nm the lattice period. In a first stage, the superlattice is treated as a continuum.

An (extended) stationary state can be characterized with the FDTD method using a fictitious “electron source”, i.e. with a suitable j in Eq. (4). The role of the source is to imitate the continuous flow of the incoming stationary wave packet. The explicit formulas that are used to generate an incoming Gaussian wave packet with energy E = E₀ are given in Appendix D. The time dependence of the fictitious source is of the form e^−iω₀t with ω₀ = E₀/ħ, and the source is turned on at time t = 0 with an initial state ψ_t=0 = 0. In order to imitate an unbounded structure the computational domain is truncated with a PML. After a sufficiently long time, the wave function will reach a steady state such that the time variation of ψ is also of the form e^−iω₀t in all points of space and the integral $∫ ψ 2 dxdy$ becomes time independent. In all the calculations of the article, we used a time step $Δ t =0.35 Δ v F 2 ≈0.62as$ for a spacing between nodes Δ = Δ_x = Δ_y = 0.25nm, consistent with the stability criterion defined in Eq. (10).

In the present example, the stationary regime in the FDTD method is reached after 16 × 10³ time steps so that the propagation time is t ≈ 0.99ps. The spatial distribution of the probability density function is depicted in Fig. 3(a) and reveals that the beamwidth of the electron wave is unaffected by the propagation in the superlattice. Note that in our plots the x- (stratification) direction is the vertical axis. This confirms that within the continuum approximation a superlattice with χ = 0 is insensitive to diffraction effects.

To validate these results we applied the FDTD algorithm to the corresponding microscopic structure taking into account that the microscopic potential has a sinusoidal-type spatial variation (see Fig. 1(a)) with V_osca/ħ v_F ≈ 7.55. This value of V_osc corresponds to a vanishing χ in the continuum model (see Fig. 1(b)).²⁷ The probability density function obtained with the microscopic approach is depicted in Fig. 3(b), and agrees extremely well with the effective medium theory results. This is corroborated by Figs. 3(c) and 3(d), which show the longitudinal and transverse profiles of the probability density function determined using the continuum and the microscopic models. Moreover, in Figs. 4(a)-4(b) we also represent the amplitude and phase of the two components of the pseudospinor $ψ= Ψ 1 , Ψ 2$ along the y = 0 line. As is well known, each component of the pseudospinor is associated with a different sublattice of graphene.^3,46 Thus, the FDTD results show that the continuum model accurately characterizes the effective response of both sublattices. It should be noted that the region x > 52a of the computational domain corresponds to the PML region, whose boundary is marked by the dashed vertical line in Fig. 3(c) and in Fig. 4. As seen, the PML region effectively “absorbs” the electron wave without reflections, mimicking an open boundary. We numerically verified (not shown) that the results are virtually unchanged for a graphene strip with a finite transverse width W along the y-direction, being W a few times larger than R_G. Thus, in the extreme anisotropy limit the edges play no role on the electron wave propagation.

It is interesting to note that due to the granularity of the superlattice the microscopic model results have considerable fluctuations on the scale of a, particularly the phase of the pseudospinor oscillates wildly in each period as seen in Figs. 4(a)-4(b). To filter out these oscillations, the microscopic results are spatially averaged, so that for a physical entity F we calculate $F a v x , y = 1 a ∫ − a / 2 a / 2 F x + x ′ , y d x ′$⁠. The transverse profiles of the amplitude and phase of the pseudospinor $ψ= Ψ 1 , Ψ 2$ calculated using the effective medium theory and the microscopic model with spatial averaging are shown in Fig. 4(c)-4(d). The results demonstrate that with the spatial averaging the microscopic and the effective medium model results are virtually coincident, which is consistent with the theory of Ref. 27. In the rest of the article, the spatial averaging is applied to all the microscopic model results.

To determine the limits of validity of the continuum approximation, we calculated the stationary regime results for Gaussian wave packets with increasingly small radial width R_G (Fig. 5). The simulations show that for R_G = 1.9a the continuum model still captures the physical response of the superlattice, but for smaller beamwidths (e.g. R_G = 0.4a) it gives results diverging from the microscopic theory. Indeed, for R_G = 0.4a the wave is diffracted by the structure leading to a significant broadening of the wave packet, as illustrated in Fig. 5(f). This finding is consistent with the general limits of validity of effective medium methods which are expected to break down when the wave becomes more localized than the lattice period.²⁷ Indeed, when R_G < a, electronic states with large values of the wave vector (e.g. states associated with new Dirac points) can influence the wave propagation. These states are not described by the continuum model.

Next, we study the wave propagation in complex graphene heterostructures. Specifically, in a first step we investigate the electron wave scattering by a superlattice nanostrip encapsulated in pristine graphene. The superlattice nanostrip has the same structural parameters as in Fig. 3 and thickness W = 20a. Figure 6 shows the transmission coefficient T for the pseudo-spinor amplitude calculated with the FDTD code for E = E₀, with E₀a/ħ v_F = 0.2, as a function of the pseudo-momentum k_y. These results are compared with “exact” analytical results for plane wave incidence obtained using the theory of Ref. 29. Note that the analytical results depend if the superlattice is regarded as a granular structure or as a continuum.²⁹ 

As seen in Fig. 6, there is an excellent agreement between the FDTD results and the analytical theory for both the microscopic (Fig. 6(b)) and the continuum formulations (Fig. 6(a)). Note that as k_ya increases from 0 to 0.2 the incidence angle varies from 0^o to 90^o (grazing incidence). For k_ya > 0.2 the incident wave is evanescent and the x-component of the incident wave vector becomes imaginary so that the states are not normalizable.

A crucial aspect is that within the macroscopic framework the wave function is not continuous at the interface, but it rather satisfies a boundary condition derived in Ref. 29 and further discussed in  Appendix A. In the FDTD algorithm the effect of the nontrivial boundary condition is described by the function $u χ$ in Eq. (5b). When $u χ$ is taken equal to zero the boundary condition reduces to the continuity of the pseudospinor.

To highlight the importance of using the correct boundary condition in the continuum approximation, we consider the scenario wherein a Gaussian electron wave with R_G = 5a, normalized energy E₀a/ħ v_F = 0.2 impinges on a superlattice nanostrip with an incidence angle θ_i = 53^o. The input interface is at x = 10a. The Gaussian wave is created by a fictitious source placed at the x = 0 plane. The graphene nanostrip is a superlattice with thickness W = 20a, average potential $V av − E 0 a/ħ v F =0.1$ and anisotropy ratio χ = 0.25 (in the microscopic model V_osca/ħ v_F ≈ 6.13). Figure 7 shows the calculated probability density function for the microscopic theory, the effective medium theory with the correct boundary condition, and the effective medium theory assuming the continuity of the pseudospinor. The last case corresponds to setting u = 0 in the update equations.

As seen in Fig. 7(a)-7(c) the continuum model only concurs with the exact microscopic theory when the correct parameter u is used. This is made clear in Fig. 7(d), which shows the longitudinal profiles of the probability density function. The red dotted curve, which corresponds to the incorrect boundary condition (i.e. to the continuity of the macroscopic pseudospinor), underestimates the value of the wave function inside the superlattice. To have a clearer idea of the discrepancy introduced by the incorrect boundary condition we show in Fig. 7(e) the transmission coefficient for a plane wave that impinges on the same graphene nanostrip.²⁹ Consistent with the results of Fig. 7(d), the incorrect boundary condition strongly underestimates the transmission across the nanostrip, especially for k_y > E₀/ħ v_F.

In the following, the problem of time evolution of a given initial electronic state is considered. This problem involves solving the FDTD equations subject to a given initial condition $ψ t = 0$⁠. Of course, in this context one does not need to consider a fictitious source and hence we set j = 0 in Eq. (4). In this work, the initial state is assumed to be of the form (the normalization of the wave function is arbitrary):

This initial state corresponds to a Gaussian wave packet initially centered at $x 0 , y 0$ and with a radial width R_G. The parameters E₀ and $k 0 = k x 0 , k y 0$ (with $E 0 − V =ħ v F k x 0 2 + χ 2 k y 0 2$⁠) play the role of “energy” and quasi-momentum of the wave packet, even though, strictly speaking the considered state is not an eigenstate of the energy operator. The vector k₀ determines the direction of propagation. Similar to the previous sections, the computational domain is surrounded by a PML to avoid reflections from the sidewalls. Thus, after a sufficiently long time the initial state will be absorbed by the PML and the probability of finding the electron inside the computational domain approaches zero.

In the first example, the initial electronic state propagates in an unbounded superlattice characterized by an anisotropy ratio χ₀ and average potential V₀ = 0, as shown in Fig. 8(a). The parameters that characterize the initial state are R_G = 2.82a, $E 0 a ħ v F =1.9$ and $k 0 = 1 . 9 , 0 /a$⁠, with a = 10nm; the initial wave packet is centered at $x 0 , y 0 = − 15 a , 0$⁠.

First, we consider the electron wave propagation in pristine graphene (χ₀ = 1). Using the FDTD algorithm the wave function was sampled at three different time instants, t = 0, t = t₁ = 2000Δ_t and t = t₂ = 4000Δ_t, with the time step Δ_t defined as in Sect. IV. Figure 8(b) represents the transverse profiles of probability density function at the x=const. lines wherein the amplitude of $ψ x , t y = 0 2$ is maximal for a fixed t. As expected, the time evolution of the initial electronic state causes the Gaussian wave packet to broaden and increase its characteristic size. Indeed, in pristine graphene there is no preferred direction of motion because the group velocity is independent of the direction of propagation, and hence the wave is diffracted. Note that in this example the microscopic and continuum results are coincident because there is no microscopic potential.

We did similar studies for the case of superlattices characterized by anisotropy ratios χ₀ = 0.7 and χ₀ = 0 (V_osca/ħ v_F ≈ 3.58 and V_osca/ħ v_F ≈ 7.55, respectively, in the microscopic model). The corresponding results are depicted in Figs. 8(c) and 8(d), respectively, and the associated time animations can be found in Ref. 45. As seen, consistent with the findings of Sect. IV, as the value of χ₀ approaches zero the broadening of the wave packet becomes insignificant. In particular, for χ₀ = 0 (Fig. 8(d)) the electronic state is unaffected by diffraction and the shape of the wave front does not change with time. Importantly, the results obtained with the exact microscopic model are largely coincident with the results obtained with our effective medium theory, confirming that the effective Hamiltonian can be used - when the initial state is less localized than the period a  − to predict the time evolution of the electronic states in a superlattice.²⁷ 

Similar to the previous section, the continuum model breaks down when the characteristic size of the initial electronic state is comparable to the period of the superlattice. This property is made clear by the results of Fig. 9, which shows that the effective medium theory still concurs extremely well with the microscopic theory for an initial state with R_G = 4a (Figs. 9(a) and 9(b)). However, when the width of the initial state is decreased to R_G = a or R_G = 0.25a (Figs. 9(c)-9(f)) the spatial spectrum of the initial state is formed by harmonics with large wave vectors that are not described by the continuum approximation, and hence the two formulations predict quite different time evolutions.

In the final example, we consider an electron wave propagating in a cascade of several superlattices embedded in pristine graphene (χ₀ = 1 and V₀ = 0) (Fig. 10(a)). With reference to Fig. 10(a), the parameters of the superlattices are: anisotropy ratio χ_1,2 = ∓0.2 (V_osc,1a/ħ v_F ≈ 8.90 and V_osc,2a/ħ v_F ≈ 6.41 in the microscopic model), average potential V_1,2a/ħ v_F = ∓0.5 and thicknesses W₁ = 7.5a and W₂ = 15a. The initial state propagates along the x-direction and has a radial width R_G = 2.82a and E₀a/ħ v_F = 1.9. Similar to the previous example, the time evolution problem is solved using the FDTD algorithm and the longitudinal profile of probability density distribution is sampled at specific time instants: t = 0, t₁ = 2000Δ_t, t₂ = 3800Δ_t, t₃ = 5600Δ_t and t₄ = 7600Δ_t. As depicted in Fig. 10(b), the sampling times are such that the maximum of $ψ 2$ lies at the mid-plane of each superlattice nanostrip in Fig. 10(a). Remarkably, despite the complexity of the microscopic potential, the results obtained with the continuum model are nearly overlapped with the microscopic theory results.

To highlight the agreement between both models, we also calculated the transverse profiles of the wave function at the same time instants. As seen in Figs. 10(c) and 10(d), there is an almost exact agreement between the two formalisms, and only in the case t = t₁, the profiles are not exactly overlapped. Figures 10(e)) and 10(f) also reveal that the agreement is similarly good for the individual components of the pseudo-spinor. Interestingly, one may detect in these profiles the presence of some trembling motion in the wave function. Such a feature is perfectly captured by the continuum approximation. The trembling motion trailing the main beam is plausibly a manifestation of the Zitterbewegung effect,^47–51 caused by the interference between positive and negative energy states. In our problem, the Zitterbewegung is due to the fact that states with energies E − V < 0 may be excited in the graphene superlattice regions and hence coexist with states with E − V > 0. Indeed, the presence of interfaces may imply a redistribution of the relative weight of the waves with E − V > 0 (forward waves) and the waves E − V < 0 (backward waves, with opposite phase direction). Note that due to the Klein tunneling^3–5 the reflections are rather weak in our problem, and hence the pulse broadening should not be attributed to them.

In this article it was numerically demonstrated that an effective Hamiltonian proposed in an earlier work (Ref. 27) can be used to characterize the electron transport and the time evolution of a given initial state in complex graphene heterostructures. To this end, a FDTD algorithm was developed to study the electron wave propagation in the frame of the effective medium theory. We derived the conditions of stability of the algorithm and proposed a way to mimic open boundaries based on an absorbing boundary condition. Importantly, the numerical results reveal that the effective Hamiltonian accurately describes both the extended stationary states (continuous spectrum) and the time dynamics of an initial wave packet, provided the initial state is not more localized than the characteristic spatial period of the superlattice. Surprisingly, this property holds even for strong modulations of the electric potential that lead to the emergence of new Dirac points (e.g., for V_osca/ħ v_F ≈ 7.55¹²). Moreover, the numerical results confirm that within the effective Hamiltonian framework the correct boundary condition does not correspond to the continuity of the pseudospinor. We envision that the developed computational tools can be useful to characterize complex graphene electronic devices and may be generalized to determine the optical response (e.g. the conductivity) of graphene superlattices.

This work was funded by Fundação para Ciência e a Tecnologia under projects PTDC/EEI-TEL/2764/2012 and UID/EEA/50008/2013. D. E. Fernandes acknowledges support by Fundação para a Ciência e a Tecnologia, Programa Operacional Potencial Humano/POPH, and the cofinancing of Fundo Social Europeu under the fellowship SFRH/BD/70893/2010.

As discussed in the main text, a graphene superlattice may be regarded as a continuum characterized by the Hamiltonian (2). This Hamiltonian can be used to describe the electron wave propagation in a homogeneous unbounded region. Crucially, it was demonstrated in Ref. 29 that for graphene-based heterostructures, e.g. at an interface between a graphene superlattice (regarded as a continuum) and pristine graphene, the boundary conditions at the interface are not trivial and do not imply the continuity of the pseudospinor. Indeed, it is possible to ensure the continuity of the probability density current j_x = v_Fψ^* ⋅ σ_x ⋅ ψ at an interface (x = x₀) between two distinct graphene based nanomaterials without having a continuous pseudospinor.²⁹ In particular in Ref. 29, a generalized boundary condition of the type $U x 0 − ⋅ψ x 0 − = U x 0 + ⋅ψ x 0 +$ was put forward to characterize a step-like discontinuity of the anisotropy ratio parameter, where U is a unitary operator of the form $U= e i u σ x = cos u i sin u i sin u cos u$⁠, with $u χ$ some function that depends on the microscopic potential. Within this theory, it is clear that the term $−iħ v F σ x ∂ ∂ x$ in Eq. (2) needs to be replaced by $−iħ v F σ x U − 1 ∂ ∂ x U$⁠. Hence, we should write the Schrödinger equation as:

Importantly, this modified equation holds at all points of space, including the interfaces. Using the formulas $U − 1 ∂ ∂ x U= ∂ ∂ x + U − 1 ⋅ ∂ U ∂ x = ∂ ∂ x + ∂ u d x U − 1 ⋅ d U d u$ and $U − 1 ⋅ d U d u =i σ x$⁠, we may simplify Eq. (A1) as follows:

The above equation is consistent with the effective Hamiltonian in Eq. (5) of the main text. The function $u χ$ may be heuristically determined with a numerical fitting as explained in Ref. 29. In case the microscopic electrostatic potential is characterized by a sinusoidal spatial variation (Fig. 1), we found that the function u is given by:

Here, it is shown that the update equations (9) are numerically stable provided the time step is sufficiently small. In the following, it is assumed that the structure is spatially homogeneous (χ and V_ef are independent of the spatial coordinates) and that there is no electron source (j = 0). Our aim is to characterize the stationary states of the system. To this end, we look for plane-wave type solutions of Eq. (9) with:

where ξ_p = e^iθ_p and ξ_q = e^iθ_q are the spatial phase-shifts between adjacent nodes. In addition, we assume a time variation such that $Ψ 1 , p , q n + 1 =λ Ψ 1 , p , q n$ and $Ψ 2 , p , q n + 1 / 2 =λ Ψ 2 , p , q n − 1 / 2$ with $λ=λ ξ p , ξ q$⁠. The algorithm is stable provided $λ ≤1$ for arbitrary values of the spatial phase-shifts ξ_p, ξ_q with $ξ p = ξ q =1$⁠. Substituting Eq. (B1) into Eq. (9) it is readily found that the plane wave states must satisfy:

The above system of equations may be written in a compact matrix form as follows:

To have nontrivial solutions, λ is required to satisfy the characteristic equation:

Using ξ_p = e^iθ_p and ξ_q = e^iθ_q it can be shown that:

The solutions of characteristic equation are:

with the parameters $A= V 2 ħ Δ t$ and $B 2 =− v F Δ t 2 D − D + >0$ real-valued. It is straightforward to check that if $B 2 2 − A 2 −1<0$ then

Thus, we can conclude that the FDTD algorithm is stable when $B 2 2 − A 2 −1<0$⁠. This condition is equivalent to

The above inequality should hold for arbitrary values of the spatial phase-shifts ξ_p = e^iθ_p and ξ_q = e^iθ_q. To satisfy this stability criterion it is sufficient to ensure that $χ 2 Δ x 2 Δ y 2 Δ x 2 + Δ y 2 χ 2 v F Δ t 2 <1$⁠, which is the same as Eq. (10).

Here, it is explained how an unbounded graphene region can be mimicked using the FDTD method. This can be done with a “perfectly matched layer” (PML) that ideally should “absorb” any incoming electron wave. We propose to do this by tailoring the electrostatic potential V such that in the PML region V is complex-valued and the wave propagation is damped. Analogous ideas are used in electromagnetics to imitate the wave propagation in unbounded regions.⁵² It should be noted that a Hamiltonian with a complex-valued V is not Hermitian and hence it is compatible with damping in the wave propagation. It can be easily verified that to have a damped wave propagation the electric potential V must have a negative imaginary part: $Im V <0$⁠. For example, for pristine graphene the relevant Hamiltonian is $H ˆ =−iħ v F σ⋅∇+V$ and hence from the Schrödinger equation one obtains a modified probability current conservation law $∇⋅ j e + ∂ W ∂ t + p loss =0$⁠, being $j e = v F ψ * ⋅ σ x ⋅ ψ x ˆ + ψ * ⋅ σ y ⋅ ψ y ˆ$ the probability current, $W= ψ 2$ the probability density and $p loss =−Im V 2 ħ ψ 2$⁠. The term p_loss is the instantaneous rate of “probability absorption” per unit of volume, and hence should be positive to have a damped propagation. Thus, it is required that $Im V <0$⁠.

In our numerical implementation the PML completely surrounds the material region of interest with dimensions L_x × L_y (Fig. 11). In a time evolution problem of a certain initial localized electronic state, the wave will eventually reach and be absorbed by the PML. Hence, after a sufficiently long time it will eventually abandon the region of interest (confined by the PML boundary) so that $∫ ψ 2 dxdy→0$⁠. The thicknesses of the vertical and horizontal PML regions are d_PML,x and d_PML,y, respectively (Fig. 11).

In our numerical implementation it is assumed that χ in the PML regions is the same as in the material adjacent to the PML. Similarly, the real part of the potential V is the same as for the neighboring material (V₀). We adopted the following potential distribution in the regions 1 and 2 of the PML:

where 0 < d_y < d_y,PML is the distance to the PML boundary (e.g. for −d_PML,y < y < 0 one has $d y = y$ and for L_y < y < L_y + d_PML,y one has d_y = y − L_y). Note that $Im V <0$ to ensure the wave absorption in the PML. Similarly, in regions 2 and 3 of the PML we use:

where 0 < d_x < d_x,PML is the distance to the relevant PML boundary. Note that in regions 2 and 3 the imaginary part of the potential is continuously increased (in absolute value) until it reaches a constant value. This ensures a good absorption of the electron waves. In our simulations we used a PML with thicknesses d_PML,x = 15a and d_PML,y = 15a.

As discussed in the main text, it is possible to study the extended electronic states (the continuous spectrum) with a given energy E = E₀ (e.g. the scattering of a plane wave by a graphene heterostructure) using the FDTD method. More specifically, an incoming incident wave can be emulated by introducing a fictitious electron source $j= j 1 j 2$ in the massless Dirac equation (4). For plane wave incidence, the following current distribution is adopted,

where ω₀ = E₀/ħ, $r 0 = x 0 , y 0$ defines the location of the source, and j₀ is the vector

In the above $k 0 = k x 0 , k y 0$ is the wave vector associated with the plane wave, which satisfies $E 0 − V =ħ v F k x 0 2 + χ 2 k y 0 2$⁠. If E₀ > V one should choose $k x 0 > 0$⁠, otherwise $k x 0 < 0$⁠.

Next it is shown that in a steady-state (t → ∞) this source emits a plane wave of the form $ψ= 1 ħ v F E 0 − V k x 0 + i χ k y 0 e i k 0 ⋅ r − r 0$ into the region x > x₀. Indeed, in a steady-state (t → ∞) the solution of (4) under the excitation (D1) satisfies: