Phase space analysis of quantum transport in electronic nanodevices

Electronic transport in nanodevices is commonly studied theoretically and numerically within the Landauer-Büttiker formalism: a device is characterized by its scattering properties to and from reservoirs connected by perfect semi-infinite leads, and transport quantities are derived from the scattering matrix. In some respects, however, the device becomes a ‘black box’ as one only analyses what goes in and out. Here we use the Husimi function as a complementary tool for quantitatively understanding transport in graphene nanodevices. It is a phase space representation of the scattering wavefunctions that allows to link the scattering matrix to a more semiclassical and intuitive description and gain additional insight in to the transport process. In this article we use the Husimi function to analyze some of the fascinating electronic transport properties of graphene, Klein tunneling and intervalley scattering, in two exemplary graphene nanodevices. By this we demonstrate the usefulness of the Husimi function in electronic nanodevices and present novel results e.g. on Klein tunneling outside the Dirac regime and intervalley scattering at a pn-junction and a tilted graphene edge.


Introduction
Typically quantum transport simulations of electronic nanodevices are based on the Landauer-Büttiker formalism. There, a nanodevice is regarded as a scattering region that is connected to electron reservoirs by semiinfinite leads (see figures 1(b), (c) for examples). The central quantity of this formalism is the scattering matrix S, from which one can obtain transmission probabilities, electric and thermal conductivities, as well as other useful quantities [1]. Even though this approach gives a wealth of information on transport through the device, in some respects the device appears to be a 'black box'.
This becomes an apparent weakness when one wants to connect the quantitative results the formalism produces to the physical intuition obtained by the semiclassical picture, or when trying to understand the role played by the different components of a complex (not easy to compartmentalize) device. In that case one wants to analyse the scattering wavefunctions inside the device, and how they populate position and momentum space. For example, if one wants to study and understand Klein tunneling (briefly reviewed in section 1.1) in a graphene nanodevice, then information about the momentum orientation ('angle of incidence') of the wavefunction inside the device (and specifically before a pn-junction) is important. In order to complement the scattering matrix information, and to get an intuitive connection with the semiclassical picture, here we will use the Husimi function Q which transforms a wavefunction into a phase space (quasi-)distribution. We note that most quantum transport simulators of course also offer the possibility of computing the quantities like the quantum current inside the nanodevice, which also allow looking into the device. It has, however, already been shown that Q is the generalization of the quantum current [2] and has numerous advantages over it 3 .
Husimi functions have been introduced to quantum mechanics a long time ago [3] and since then have been used in various areas of physics, like quantum optics [4,5] and ocean acoustics [6]. Husimi functions also play a prominent role in the field of quantum chaos, which tries to unravel the properties of complex quantum systems. For example, they have been used to understand the structure of the eigenfunctions in paradigmatic chaotic systems like quantum maps and billiards [7][8][9][10][11], transport in quantum ratchets [12], the properties of optical microdisc lasers [13][14][15] and even electronic transport in disordered systems [16]. Recently Mason et al [2,17,18] have introduced a processed Husimi map in tight-binding models of nanodevices allowing to recover and visualize classical paths in coordinate space. But in general solid state physics does not yet take much advantage of this useful tool.
In this article we will use the Husimi function to analyse graphene nanodevices. Graphene is a fascinating material for studying quantum transport, due to the abundance of new physics it brought into light quickly after its discovery [19] e.g. weak (anti-)localization effects connected with the existence of two inequivalent valleys [20] or the Klein tunnel effect and its potential impact on technological applications [21]. Motivated by the intention to understand the impact of Klein tunneling (section 1.1) and intervalley scattering on transport in arbitrary graphene nanodevices, in this article we will use the Husimi function as a distribution in position and momentum space to analyse two exemplary devices. We find that the the combination of Husimi functions and semiclassical considerations is a powerful tool to understand transport phenomena in graphene nanodevices. The observations we report also include the mode dependence of intervalley scattering at pn-junctions and the behaviour of Klein tunnelling at trigonal warping (section 2) as well as the quantification of intervalley scattering at a tilted graphene edge (section 3).

Klein tunneling in graphene
Klein tunneling is a paradox of relativistic quantum mechanics first discovered by Klein in 1929 [22, 23], but it has become a reality in graphene [21]. This is due to the fact that low-energy excitations in graphene are well described by the Dirac equation for massless relativistic particles with linear energy dispersion, | |  k , and the wave function on the honey comb lattice of graphene, separated in two sublattices, is represented by a spinor, corresponding to the Dirac spinor. In this formalism chirality-conservation leads to a fascinating transport anomaly (for more details and a pedagogical review see [24]): conduction electrons are able to penetrate (almost) arbitrarily high potential barriers! If a conduction electron impinges normally on a potential step higher than the Fermi energy it will not get reflected but fully transmitted into a symmetric state in the valence band inside the high potential region instead.
More quantitatively, consider a plane wave (an electronic eigenstate of the Dirac equation [24]) with wavevector k incident on a pn-junction (i.e. a potential step), like the one of figure 1(b). The transmission probability depends on the angle of the wavevector (the angle of incidence) and on the width of the junction w [21,25], where for a very steep junction (  w 0) one expects Step in in out Here λ F is the Fermi wavelength of the electrons, f out the angle of the wavevector of the transmitted wave (the angle of refraction) and k 1,2 are the wavenumbers of the incident and the transmitted wave, respectively. In both expressions, at normal incidence, i.e. for f in =0, the transmission is perfect, T=1, independent of the height and shape of the potential step! The transmission probability in the Klein tunneling process crucially depends on the angle of incidence of the electron on the junction. This angle is clearly defined for a plane wave: the wavevector angle. But what is the 'angle of incidence' of a more complex electron wavefunction in a realistic scenario, like the scattering wavefunction in a nanodevice of arbitrary spatial complexity? Basic examples of such wavefunctions are shown in figures 1(d), (e). In general a plane wave will not always be a good representation of these wavefunctions. Yet, if one wants to understand the impact of Klein tunneling on transport in such nanodevices, one needs information about an angle of incidence. Because such information is inherently present in the (classical) phase space, below we will use the Husimi function to obtain this information. But first let us introduce the nanodevices we will be using as examples.

Graphene devices
In our work we will study transport in tight-binding models of the two graphene-based devices shown in figure 1. Notice however that our methodology can be applied to tight-binding based quantum transport in any device (and in 3 dimensions just a well), as it will become clear below. For the graphene devices, Device A is the conceptually simplest device in which one can study Klein tunneling in a realistic scenario (i.e. a finite nanodevice): a graphene nanoribbon (GNR) of constant width with a pn-junction in its middle. In device A the boundary conditions are chosen such that it forms a 'zigzag' nanoribbon. These have been studied by Bray and Fertig in detail within the Dirac approximation [26] and many of their properties are known analytically (for small Fermi energies). Analytical descriptions in this case are possible because of the many symmetries that are present. Device B however breaks the conservation of k y as well as the reflection symmetry along the x axis. Note also that for ω=π/6 the 'scattering edge' in device B (highlighted in green in figure 1(c)) exactly is an armchair boundary. In both devices we create pn-junctions via a linear increment of the potential energy from the n region with bias voltage -V 2 0 to the p region with bias +V 2 0 over a range w (see figure 1). The kinetic energy E of the incoming electrons is connected to the Fermi energy E F by = + E E V 2 F 0 . Zigzag GNRs have a dispersion relation shown in figure 1(a) and discussed in detail in [26,27]. For a given Fermi energy M bands of the dispersion intersect the energy level at positive slope, thus having positive group velocity. This results in M incoming modes (M is always odd and scales linearly with the width of the GNR). We order the modes by decreasing k x , as shown in figure 1(a). Importantly, for small energies the two (inequivalent) Dirac valleys ¢ K K , are well separated in momentum space, which leads to the incoming modes being valleypolarized. This means that modes 1 to⌊ ⌋ M 2 (where⌊·⌋denotes the integer part) come from valley ¢ K , while modes⌊ ⌋ + M 2 1to M come from valley K. We also stress that K has one additional incoming mode, see figure 1(a).
All of our quantum transport simulations are tight binding calculations performed with the software Kwant [28]. The devices are finite scattering regions that are coupled to semi-infinite leads (which are also GNR). The modes (eigenfunctions) of the leads enter the device and are subsequently scattered, defining the scattering wavefunctions y m for each mode. As we will consider transport always from the left to the right lead, we will only use a part of the scattering matrix which we define as the N×M transmission matrix  , where M and N are the total number of modes in the left and right lead, respectively. The element  nm is the transmission amplitude from the m-th (incoming) mode of the left lead to the n-th (outgoing) mode of the right lead. The total transmission probability T m of each individual incoming mode is given by  (with d =r r r 0 and D spatial dimensions) which is a Gaussian envelope in space with origin r 0 multiplying a plane wave with wavevector k 0 . 1 is the normalization factor in the case of continuous space, so that | á ñ = W W 1 and that Δ x=Δ y=σ. The key property of these wavepackets is that they minimize the uncertainty relation between position and momentum. Here σ is the spatial uncertainty and thus is a parameter that controls the trade-off between the uncertainty in position (σ) or momentum space (1/(2σ)).
The Husimi function Q is defined as the magnitude of a projection of a wavefunction onto | ñ where for continuous space systems we have , , ; 6 0 0 0 0 * where the integration extends over the full spatial domain of the device (in our case in two dimensions). For a tight-binding system the projection is turned into a sum due to the discrete nature of the lattice j being the wavefunction at lattice site j with position r j 4 . The normalization factors here depend on the lattice and thus in the following we will omit normalizations.
We stress that Q is a rigorous method for transforming a wavefunction into a phase space (quasi-probability) distribution, entirely consistent with the framework of quantum mechanics. Q is the Weierstrass transform of the Wigner function [29] (which is also entirely consistent with quantum mechanics) and have been used as a versatile tool for understanding complex quantum and other wave dynamics in several scenarios [2,[6][7][8][9][10][11][12][13][14][15][16][17][18][29][30][31][32]. Q also respects the uncertainty principle, which is why it is a quasi-probability distribution: it drops Kolmogorov's third axiom because all points in a phase space box of size~ D are indistinguishable from each other. While the marginals of Q are not the true position and momentum distributions (i.e. the amplitudes of the position and momentum wavefunctions), they are just smoothened versions of them with smoothening factor σ or ( ) s 1 2 so the information loss in sufficiently semiclassical systems is not only small, but also entirely controllable.

Husimi function in a lead
Before moving on to the numeric applications of this paper, it is helpful to obtain some intuition about Q in analytically treatable examples. Let us first simply consider a plane wave ( ) exp . It is straightforward to calculate its Husimi function with We thus find with ¢ =k k k 0 . As one would expect the Husimi function Q does not depend on r 0 , since a plane wave is spatially homogenous, and Q is a Gaussian in momentum space since the integral above is the Fourier transform of a Gaussian. In addition, Q only depends on the difference between the wavevectors of the Gaussian wave packet and the plane wave.
Let us now examine the case of a lead (or waveguide). Here the wavefunction X is a plane wave in the longitudinal direction, while the transverse part is the quantum well wavefunction, i.e.
Notice that the following analytic result applies to both square lattice leads (whose low-energy excitations follow a dispersion k 2 and we know explicitly for a lead of width W), but also zigzag graphene 4 For each r 0 we use only lattice sites that are within | | s - r r 3 j 0 , to reduce computation time. Notice the complex conjugation y* in equations (6) and (7), which sometimes is omitted in the literature. While for closed systems with time-reversal symmetry it can be omitted, it is crucial for open systems, which we explore here, and for systems with broken time reversal symmetry like those in magnetic fields.
nanoribbons (GNRs) because the transverse wavefunctions there are as well sine modes (see end of section 1.2). The only difference is that for GNRs the expression for k m is given by equation (13).
As we know that in the Husimi function the longitudinal component will result to = ¢ s -Q e x k 2 x 2 2 , we now tend to the integral of the sine function with the Gaussian function in finite limits. Although analytically solvable, its expression is not so trivial In our later analysis we will use , with an appropriate k F , as discussed in further detail in section 1.5. To obtain an intuition around equation (10) we are plotting it in figure 2. Notice that ). Based on this Q y we can define a marginal distribution over f as Here there are two properties to point out. Q is symmetric over k m =0 (and equivalently f = 0). This is already true from equation (10), but it can also be understood simply from the fact that ( ) k y sin m decomposes to , which is a superposition of two plane waves with opposite directions. Furthermore, as y 0 moves further away from the center of the lead W 2, the accuracy of the Husimi function drops significantly ( figure 2(d)). This has implications for calculating Husimi functions exactly at the edges of a tight-binding device, which we will discuss in section 2.2.  (11)) for different m. (d) Cuts through Q y for different y 0 (also shown in panel (b)).

Calculating the Husimi function
For each scattering wavefunction ψ m we compute the Husimi function Q. To reduce the dimensionality of (the parameters of) Q we only evaluate it at well chosen transverse cuts (x 0 =const.), e.g. just in front or behind the pn-junction. We will thus obtain a distribution of incoming and outgoing wavevectors that 'pass through' these cuts as a function of the transverse coordinate y 0 .
We further reduce the dimensionality of Q by exploiting energy conservation. In figure 3(a) we show the two-dimensional dispersion of graphene with l = 1 the band index, t=2.8 eV the hopping constant and » a 0.142 nm the carbon-carbon distance. Measuring the Husimi function we see that Q localizes strongly on the two-dimensional energy contour corresponding to the Fermi energy (even though in reality the dispersion relation of a finite GNR is in fact onedimensional). This allows us to reduce the dimensionality of Q by using the 2D dispersion. In the following we measure Q for wavevectors that populate the Fermi energy contour at equally spaced angles f.
A difference with the honeycomb lattice versus the square lattice is that there are six valleys in the twodimensional dispersion, as seen in figure 3. Thus we compute Q in all six of them, and the angle is measured with respect to the respective Dirac point K ξ of each valley, i.e. f = ¢ ¢ k k arctan y x with ¢ =x k k K . We denote this by ( ) . This reduces Q from depending on both k k , x y to be only a function of f. The parameter σ we will choose such that the wavevector uncertainty satisfies where k is the (average) magnitude of the wavevector with respect to Dirac point (see appendix). This yields typical values of s » 8 nm for small energies, while for higher energies σ can be smaller than 4 nm. In the following and for device A the notation = x n will denote a cut in the n region of the device, 3σ before the pn-junction, while x=p will denote a cut 3σ after the junction. For device B the slice location is given explicitly (in the rest of the text we measure space in nm and energy in eV).
For the zigzag nanoribbon we can compare the numerically computed Q with the analytical expression(10) (because we keep x fixed, we in principle calculate Q y with our numerical scheme). The theory of Brey and Fertig shows that in the Dirac regime the transverse wavefunctions of a GNR are sine modes, ( ) k y sin m [26] and thus their Husimi function is given by. (10). To obtain the transverse wavevectors k m for zigzag GNRs one needs to solve [26] (14) (averaged only over the three appropriate valleys instead of all six) with red color.

Klein tunneling
We first want to test the usage of Q in a well studied situation where much can be inferred analytically: Klein tunneling in device A at small energies [26] (see [24] for a review on Klein tunneling in graphene). Figure 3 The top panels show Q for valley K 2 , the bottom for ¢ K 2 . We show Q both before (incoming & reflected) and after (outgoing) the pn-junction for two modes. What we have seen is that for device A before the junction, Q in valley ( ) ¢ K 2 is the mirror reflection of Q in valley ( ) ¢ K 3 while in valley ( ) ¢ K 1 we find an almost exact superposition of the Qs in ( ) ¢ K 1 and ( ) ¢ K 2 . Figure 3 shows that for all modes the incoming Q nicely localizes at a single angle. We also show in red the marginal distributions ; , . 14 W 0 Because in this setup the incoming Q is highly localized, we do not need the entire distribution and can simply choose the maximum location of ( ) f Q , Φ, to represent the 'incoming angle' for each mode ; p , for 0, 2 15 (we use Q of valley ( ) ¢ K 2 exclusively for this, and we also know which of the two valleys is the incoming one). As discussed in the end of section 1.5, we can obtain numerically the wavevector of the transverse component of the wavefunction k m , based on the theory of Bray and Fertig [26] (notice that this is invalid outside the Dirac regime). Kwant also provides k x , the wavevector of the longitudinal component. We then compare Φ with ( ) n = k k arctan m x in figure 5(a). We see that only for the 'highest' (meaning high energy band) modes of each cone Φ does not have a perfect agreement with ν. We now want to use Φ to compare the results of the tightbinding calculations with theoretical result for the Klein tunneling at a pn-junction, equations (2) and (1). In figure 5(c) we plot the theoretical curves and the values of T m versus Φ for each mode, for two different pnjunction widths w, and find very good agreement. This does not only hold for the case of a symmetric pnjunction, i.e. E F =0, but also for higher and lower Fermi energies, as shown in figure 5(d) for w=10 nm (for other parameter values we also find excellent agreement).
It is clear that through Q we can find the parameter f in . Now we want to show that we can even obtain the transmission probabilities from the Husimi function using the theoretical transmission formulas. Using the marginal distribution of equation (14) we can compute the transmission of a mode as the average where T si either equations (1) or (2). Notice that in principle equation (16) could be resolved analytically, since we know the expressions for both T as well as Q (see equation (10)) for the simple device A. Unfortunately, we were not able to indeed resolve the integral analytically, but numeric integration is always possible. In figure 5(b) we compare á ñ T with T m and again we find a near perfect match (also for many more parameters than the ones shown). Equation (16) will also give a good estimate of the transmission value in cases where the distribution is not strongly localized at a single angle, allowing us to use the integrated transmission in more complicated cases like those in section 3.

Intervalley scattering
We now turn to study intervalley scattering, which describes the scattering of a wavefunction from one valley to another (inequivalent) one, e.g. from K to ¢ K . We discussed in section 1.2 that for zigzag GNRs and low energies every incoming mode is valley-polarized [26]. Intervalley scattering has found considerable interest in the literature, and was first discussed in the context of weak localization [20,[35][36][37]. Later work focused on valley filters and valley 'spintronics', see [38][39][40] and references therein. The discussions in the literature so far have been qualitative and mostly theoretical.
The Husimi function is an excellent tool to study intervalley scattering, because it directly provides information in momentum space at different positions in the device. In fact, Mason et al have used a processed Husimi projection technique in [2] to study intervalley scattering in graphene billiards. Here we will use a simpler approach directly using the Husimi function. As one can already see from figure 3 p p Î -2, 2 ) has most weight in one valley (the 'incoming valley') V i , while the other (the 'complementary') valley V c contains either just noise or only the reflected wave (compare the scales of the colorbars). In panels (c, e) it is evident there exist modes that undergo intervalley scattering, as for panel (e) the outgoing valley ¢ K has significantly more weight than what it had in the incoming case of panel (d). We want to define two intuitive measures for intervalley scattering. We first define the following weights (the sums are over all equivalent valleys) n, 17 α, is used for the normalization to the incoming mode. β and γ measure the weights of the transmitted wave that are localized in the same valley as the incoming mode and its complement, respectively. With these quantities we define  2 come from ¢ K while the higher modes come from the K valley which has an additional incoming band (see figure 1(a) or figure 5(g)). The perplexing result of figure 5(e) can be qualitatively explained based on this extra mode and the unitarity of the scattering matrix S [1] (i.e. current conservation). To aid the following argument, in figure 5(g) we show a sketch of where is each incoming mode transmitted. The lines connecting incoming and outgoing modes have widths directly proportional to the transmission amplitude | |  im 2 . After transmission, each mode 'tries' to scatter into a the same valley at negative energy to conserve the valley pseudospin (green dots in figure 5(g)). Likewise should the reflected part scatter into modes in the same valley at the same energy level but with negative group velocity. Modes 1 to⌊ ⌋ M 2 have no problem achieving this, as within their valley the outgoing channels are more than the incoming ones and thus available channels always exist. This is not the case however for modes⌊ ⌋ + M 2 1to M, since the number of outgoing channels within the same valley is one less, both for transmission and reflection. As the mode number increases the outgoing channels are filled and the higher modes have to move some of their weight to other channels (as a specific outgoing channel cannot be filled with more than total transmission of 1, see [1]). The only remaining channels that can accommodate these modes exist in the ¢ K valley (right valley of figure 1(a)) which leads to intervalley scattering.

Trigonal warping and Klein tunneling
Klein tunneling applies to graphene because for small energies the Dirac equation is a valid approximation. In Klein tunneling the important angle is the wavevector angle (with respect to the Dirac points), see equations (1), (2). The group velocity angle θ coincides with f for small energies, however as the energy increases and trigonal warping effects begin to be significant, this is not the case anymore and q f ¹ [19]. As there is no theoretical result on the tunneling behavior of graphene for energies beyond the Dirac regime, one is left to wonder: for higher energies is the Klein tunneling picture still relevant? And if yes, are the tunneling properties still dictated by f? This is an interesting question since the physical propagation direction is governed by θ. We can answer this using the Husimi function. We significantly increase the energies in device A, setting = V 5 0 and keeping E F =0, yielding incoming energy of = » E t 2.5 0.9 which shows strong trigonal warping. Once again we compute incoming angles using Φ as in equation (15)  , due to the warping of the energy contour, see below. This regime is not covered by [26] and so the transverse wavefunctions are not necessarily sines (and we found Q to perform much better than attempting to anyway fit sines to the transverse wavefunctions).
At higher energies the two valleys provide very different incoming modes, as seen in figure 6(c). For valley ¢ K there is a 'flat' front, greatly limiting the possible group velocities. The contrary is happening in valley K where the contour with positive group velocity spans more angles. In addition, in the ¢ K case the incoming wavevector angle is limited to (k k , x y are measured with respect to the center of the BZ here). Klein tunneling assumes equivalence between the two valleys as it depends on the wavevector angle. To see whether some remnant of Klein tunneling exists at higher energies, we have to look for some tunneling property that not only decays exponentially with increasing angle of incidence, but also stays 'as similar' as possible between the two valleys. In figures 6(d), (e) we compare the transmission probability of each mode T m versus the wavevector angle Φ and group velocity angle θ.
The result surprised us, since we find a Klein tunneling-like behaviour in T m versus f. We were rather expecting T m versus θ to show similar behaviour at the two valleys, because θ corresponds to the physical propagation direction. We do not suggest that Klein tunneling straightforwardly applies to higher energies. In figure 6(f) the characteristic perfect transmission at normal incidence (f = 0 in ) is lost, nevertheless, it is clear that the tunneling probability as a function of the wavevector angle is quite similar to what would be expected for Klein tunneling.

Applications of the Husimi function: device B
In this section we study transport through the asymmetric device B (see figure 1(c)) in which the incoming modes are scattered both from the boundary ('scattering edge', highlighted in green) and the pn-junction. There are two main questions we want to address. First, to what extend can we use the existing expressions describing Klein tunneling to understand the transmission properties of such a device? These expressions are derived for plane waves, which have infinite spatial extend and are characterized by a single angle f in . Due to the boundary induced scattering the wavefunction in device B cannot be well approximated by a single plane wave. Can we use the Husimi technique to connect the transmission through the device to Klein tunneling? And also, how much can we push this technique, with respect to the physical size of the configurations we can examine?
Second, we want to understand how the type of the scattering edge affects intervalley scattering. There is strong theoretical evidence that the armchair termination is in some way unique, while a random termination behaves like zigzag [27,41,42]. In addition, in the theoretical treatment of graphene nanoribbons in [26], the authors showed that the armchair termination mixes valleys while the zigzag keeps them separated. These  generally enhanced at armchair boundaries. Here we quantify this effect by using the Husimi function, similarly as in 2.2 and we will show that intervalley scattering is indeed enhanced drastically at armchair edges.
Let us stress that in device B the lead modes and the angles ν are not of much use, since the waves are deflected by the titled boundary of device B and also because k y is not conserved until the pn-junction. On the other hand, ( ) f Q is just as valid here as it was in section 2.1. It also becomes clear from figure 7(a) that many of the scattering waves inside L 2 cannot be approximated using a single angle, which means that one needs the entire distribution.

Tunneling
We find that we can apply the Klein tunneling formulas 'locally' even in small devices and when the incoming waves are not single plane waves. We show this numerically using the integrated transmission formula, equation (16) with Q measured at location = + x L L 2 1 2 (which is 3σ before the pn-junction). However, now we can't compare á ñ T with T m directly, because T m also accounts for the back-scattering from the boundary inside L 1 . To compensate for that, we compute the transmission of equation (3) once without any pn-junction at all. We call this quantity T 0 . We now have to compare · á ñ T T 0 with T m , which we do in figures 7(c)-(f) for various orientations of the boundary.
We see that the integrated transmission matches the transmission obtained through the pn-junction (using the scattering matrix) very well. This good agreement means that the Klein tunneling formula still locally describes the tunneling properties at the pn-junction, even when the nanodevice is small (we found good results for W 2 as small as 20 nm) and the incoming wave is not a simple plane wave. In addition, this also means that the Husimi function accurately decomposed the incoming scattering wave into a representative distribution of angles of incidence. From this we can see that Q allows us to separate the contributions of Klein tunneling (or any other transmission function ( ) ¼ T ) in T m , which could be useful in other scenarios as well.

Intervalley scattering
We now want to explore the intervalley scattering induced by the scattering edge and not the pn-junction. Therefore we first obtain the scattering wavefunctions y m in device B without a pn-junction (i.e. ). We measure Q using a slice at x=L 1 (exactly where the scattering boundary ends) and we compute I 2 from Q there. The results are shown in figures 7(g)-(i). We note that the results we display below do not have significant differences if one uses slices at x>L 1 .
An important benefit of using I 2 (over e.g. measures used in [17]) is that it does not depend on, or demands measuring Q for r 0 exactly at the boundaries of the nanodevice. This is crucial as the accuracy of the Husimi function dramatically drops at the boundaries, since most lattice sites around a circle of 3σ from r 0 do not even exist. We have observed in our simulations that this leads to numeric artifacts and should be avoided, and we have also established this to be true in the analytic treatment of Q in section 1.4.
There are two interesting observations to be made. First, the intervalley scattering from a lattice termination is fundamentally different from that seen in section 2.2 which results from a pn-junction. In the present case both valleys always undergo intervalley scattering. The second observation is what we expected from existing theory and now quantified using a well-defined measure: armchair lattice terminations induce much more intervalley scattering than any other termination orientation. This can be seen firstly in figures 7(g), (h) where I 2 has clearly higher values, but most prominently in panel i where we plot the average intervalley scattering per mode, i.e. 2 I 2 has a very sharp peak at ω=π/6, where the boundary termination is exactly armchair.

Conclusion
In this paper we have used the Husimi function both numerically and analytically to analyse quantum transport through tight-binding nanodevices and have demonstrated that it is a very useful tool. We have for example shown that even in situations where the angle of incidence on a tunnel barrier is not easily discernible we can use the Husimi distribution to evaluate Klein tunneling at this barrier. For higher Fermi energies the Husimi function allowed us to analyse the tunneling behavior in the regime of triangular warped Dirac cones. We have also shown how Q can be straightforwardly used to accurately define and measure intervalley scattering. Through this we have shown that pn-junctions not only introduce intervalley scattering, but that is unexpectedly strongly valley-asymmetric. We also confirmed quantitatively that the strongest geometric intervalley scatterer in graphene is indeed the armchair termination.
A main goal of this work was to show that the Husimi function is a helpful tool to augment the Landauer-Büttiker transport formalism and is useful to have in the toolbox of condensed matter physicists. The Husimi function complements, and not competes with, the scattering matrix approach. Notice that to compute Q numerically one needs the scattering wavefunctions. What the Husimi function is able to do is to offer an additional level of depth that allows one to look directly into the device and even specific parts of the device (in both coordinate and momentum space at the same time, not possible just from the scattering wavefunctions ψ m ). Whether this level of detail is necessary or useful depends of course on the exact problem one wants to study, and thus cannot be discussed generally. What is important is that if such a level of detail is sought after, the Husimi function can provide it. We also point out that even though the current work was applied to graphene, the methodology involving the Husimi function is in no way limited to it.

Acknowledgments
The authors would like to thank Theo Geisel for helpful contributions.

Appendix. Angle uncertainty
For a given value of the parameter σ, the wavepacket has a known uncertainty in both position and momentum ≔ ( ) s D = D x k What we are interested about is the uncertainty in the propagation angle. For small energies the propagation angle is the same for the wavevector and the group velocity defined as . For any nonlinear function, uncertainty propagation is given by therefore we see that If we want to have a constant s f for measurements at different energies, then we will use σ such that (assuming also D D = x k 1 2)