Symmetry breaking and (pseudo)spin polarization in Veselago lenses for massless Dirac fermions
K. J. A. Reijnders, M. I. Katsnelson
TL;DR
This paper investigates how polarization in Veselago lenses for massless Dirac fermions causes symmetry breaking and spin polarization effects, with analytical and semiclassical methods revealing displacement and current differences.
Contribution
It introduces a semiclassical approach to analyze symmetry breaking and spin polarization effects in Veselago lensing for Dirac fermions, providing analytical formulas and identifying scattering regimes.
Findings
Polarization causes spatial symmetry breaking in Veselago lensing.
Analytical formula for vertical displacement of the main focus.
Identification of two scattering regimes in current injection.
Abstract
We study Veselago lensing of massless Dirac fermions by n-p junctions for electron sources with a certain polarization. This polarization corresponds to pseudospin for graphene and to real spin for topological insulators. Both for a point source and for injection into a sample through a narrow lead, we find that polarization leads to spatial symmetry breaking. For the Green's function, this results in a vertical displacement, or even complete vanishing of the main focus, depending on the exact polarization. For injection through a lead, it leads to a difference between the amounts of current emitted with positive and negative transversal momenta. We study both systems in detail using the semiclassical approximation. By comparing the results to the exact solutions, we establish that semiclassical methods provide a very effective way to study these systems. For the Green's function, we…
Click any figure to enlarge with its caption.
Figure 1
Figure 10
Figure 11
Figure 2
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 11Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Symmetry breaking and (pseudo)spin polarization in Veselago lenses for massless Dirac fermions
K. J. A. Reijnders, M. I. Katsnelson
Radboud University, Institute for Molecules and Materials, Heyendaalseweg 135, 6525 AJ Nijmegen, The Netherlands
(March 7, 2024)
Abstract
We study Veselago lensing of massless Dirac fermions by n-p junctions for electron sources with a certain polarization. This polarization corresponds to pseudospin for graphene and to real spin for topological insulators. Both for a point source and for injection into a sample through a narrow lead, we find that polarization leads to spatial symmetry breaking. For the Green’s function, this results in a vertical displacement, or even complete vanishing of the main focus, depending on the exact polarization. For injection through a lead, it leads to a difference between the amounts of current emitted with positive and negative transversal momenta. We study both systems in detail using the semiclassical approximation. By comparing the results to the exact solutions, we establish that semiclassical methods provide a very effective way to study these systems. For the Green’s function, we derive an easy-to-use analytical formula for the vertical displacement of the main focus. For current injection through a lead, we use semiclassical methods to identify two different scattering regimes.
I Introduction
Focussing is an effect well-known in optics, where light rays are refracted by a lens to create spots of high intensity. A particular kind of lens was proposed by Veselago, Veselago (1968) who investigated lenses made of materials with a negative refractive index. Such lenses have already been realized in metamaterials, Smith et al. (2000); Houck et al. (2003); Grbic and Eleftheriades (2004) chiral metamaterials Tretyakov et al. (2003); Pendry (2004); Zhang et al. (2009); Xiong et al. (2010) and photonic crystals, Parimi et al. (2004); Cubukcu et al. (2003) and can be used to produce an image with subwavelength resolution. Pendry (2000); Grbic and Eleftheriades (2004) Soon after the discovery of graphene, Cheianov et al. Cheianov et al. (2007) realized that n-p junctions in this material would be ideally suited to create an electronic analog of a Veselago lens.
Graphene is a two-dimensional gapless semiconductor, whose low-energy charge carriers are governed by the Dirac equation. Wallace (1947); McClure (1956); Slonczewski and Weiss (1958); Semenoff (1984); Castro Neto et al. (2009); Katsnelson (2013) This results in peculiar behavior of its electrons, most notably Klein tunneling Klein (1929); Katsnelson et al. (2006); Cheianov and Fal’ko (2006); Shytov et al. (2008); Tudorovskiy et al. (2012); Reijnders et al. (2013): an electron normally incident on a potential barrier is transmitted with unit probability. The transmission probability decreases as the angle of incidence increases, meaning that an electron beam is collimated. Cheianov and Fal’ko (2006); Shytov et al. (2008); Tudorovskiy et al. (2012); Reijnders et al. (2013) A few years after its discovery, Klein tunneling was shown experimentally, Young and Kim (2009); Stander et al. (2009) and recent experiments show that it is still in center of attention today. Chen et al. (2016); Lee et al. (2016a); Gutiérrez et al. (2016) Notably, the angular dependence of the transmission coefficient has recently been measured. Chen et al. (2016)
A graphene n-p junction exhibits (Veselago) lensing, because for electrons the group velocity is parallel to the phase velocity, whereas for holes the group velocity is opposite to the phase velocity. Cheianov et al. (2007) Klein tunneling is crucial in this process, because it makes the n-p interface highly transparent to electrons. Recently, two experimental groups have demonstrated Veselago lensing in graphene samples. In the first experiment, Lee et al. (2015) the authors measured ballistic transport accross a graphene device, and found an increase in the background-subtracted current in the bipolar regime. In the second experiment, Chen et al. (2016) transverse magnetic focussing was employed to show Veselago lensing, allowing the authors to simultaneously measure the angle-dependent transmission coefficient.
Theoretical papers on the subject have considered the Green’s function for a quasi-one-dimensional n-p junction, Cheianov et al. (2007) or have looked at circular n-p junctions. Cserti et al. (2007); Péterfalvi et al. (2010); Wu and Fogler (2014) In the latter case, semiclassical considerations were also presented, Péterfalvi et al. (2010); Wu and Fogler (2014) though the semiclassical approximation to the wavefunction near the main focus was not computed. Another study Choi et al. (2014) considered a Veselago lens in a graphene nanoribbon, and showed that the geometrical phase that is acquired when scattering off a zigzag edge influences the interference pattern. Finally, a numerical study Milovanović et al. (2015) was conducted where the authors considered n-p junctions in graphene samples of realistic size, with current entering from a narrow lead. The authors compared their findings to a semiclassical billiard model, Beenakker and van Houten (1989); Milovanović et al. (2013) and generally found good agreement.
In this paper, we perform a theoretical study of Veselago lenses formed by quasi-one-dimensional n-p junctions. We mainly consider the Green’s function, although at the end of the paper we also briefly consider the situation where current flows into a sample through a lead that is attached on one of its sides. Our emphasis is not so much on the classical focussing, but rather on the matrix character of the Dirac Hamiltonian and on how it influences the interference pattern. In particular, we consider the case where the point source or incoming wave has a certain sublattice or pseudospin polarization, meaning that the current is not equally distributed among the two graphene sublattices. The fact that we are dealing with spinors makes this problem different from optical problems, where one is usually concerned with the Helmholtz equation.
We study these interference effects using the semiclassical approximation, which is valid when we have a small parameter in our problem. Earlier studies Tudorovskiy et al. (2012); Reijnders et al. (2013) have shown that this requires either large length scales or high energies. The first step of our analysis is to carefully review the classical problem, for which we need a few elements of the general theory of caustics and wave fronts, a theory known as catastrophe theory. Berry and Upstill (1980); Poston and Stewart (1978); Arnold et al. (1982); Arnold (1975) We then apply the stationary phase approximation Fedoryuk (1977); Guillemin and Sternberg (1977); Maslov and Fedoryuk (1981) to our solution. However, the simplest form of this approximation fails near the main focus, where our primary interest lies. In order to quantitatively study interference effects near the main focus, we therefore employ the Pearcey approximation. Pearcey (1946); Connor and Farrelly (1981a); S. Yu. Dobrokhotov et al. (2014) We also briefly consider the uniform approximation, Ursell (1972); Connor and Farrelly (1981a) which, in a different form, was successfully applied to graphene for a rather large semiclassical parameter. Reijnders et al. (2013) Because we compare the various approximations with the exact solution, our study can also be considered as a benchmark for the application of various semiclassical methods to graphene.
One of our interests is to see if pseudospin polarization could lead to symmetry breaking between the and -valleys in graphene. If this is the case, then it may provide another way of creating valley polarization in graphene. Rycerz et al. (2007) Since charge carriers in both valleys obey the same classical Hamiltonian, Katsnelson (2013) it is clear that the valley polarization we are looking for can only result from quantum interference. Therefore, it is unlikely that a polarization of 100% could be realized in our system. Such a polarization could for instance be detected using the valley Hall effect Xiao et al. (2007); Gorbachev et al. (2014) or second harmonic generation. Wehling et al. (2015)
We believe that there may be ways to realize such a sublattice polarization in graphene experimentally. Firstly, one could inject electrons on a single site using a scanning tunneling microscope (STM) with an atomically sharp tip. Secondly, one could consider a device where electrons tunnel into a graphene layer through hexagonal boron nitride (h–BN). Because the strengths of the carbon–nitrogen and carbon–boron interactions differ, Sachs et al. (2011); Bokdam et al. (2014); van Wijk et al. (2014) this could lead to an asymmetry between graphene’s sublattices. In this context, we note that it has recently been shown experimentally that a device with a few layers of h–BN between two layers of (bilayer) graphene can be used to manipulate the valley and pseudospin state of Dirac electrons. Wallbank et al. (2016) We believe that a graphene sample with current flowing in through a graphene lead at one of its sides would be easier to realize. Here, one could create an initial (i.e. in the lead) sublattice polarization by using a substrate that acts differently on both sublattices, giving rise to a mass term in the Dirac equation. Katsnelson (2013)
Regarding the experimental realization of high-energy states in graphene, which would be needed to study the deep semiclassical limit, we note that it is possible to create hole-doped states with energies around 0.5 – 0.6 eV. This can be achieved by molecular doping Bae et al. (2010); Wehling et al. (2008); Nair et al. (2013) with HNO3 or NO2, but can also be reached on a SiO2/Si substrate after proton irradiation. Lee et al. (2016b). Electron doping can for instance be achieved using aniline Liu et al. (2011), with which one can reach energies of about 0.25 eV. Doping graphene with alkali metals, such as lithium, very high electron doping above 1 eV can be achieved, Khademi et al. (2016) which can also induce superconductivity. Profeta et al. (2012); Ludbrook et al. (2015)
In our theoretical considerations, we consider a sharp n-p junction. Although semiclassical tunneling has been extensively studied for smooth n-p junctions, Shytov et al. (2008); Tudorovskiy et al. (2012); Reijnders et al. (2013) studying the Green’s function for such a junction is far from straightforward. Although considering a sharp barrier is not very realistic from an experimental point of view, Chen et al. (2016) we do not believe that this will significantly influence the main results. An indication for this is given by the aforementioned numerical study, Milovanović et al. (2015) where the authors found that in going from a sharp barrier to a smooth barrier the main features of the results were preserved. One notable effect should be the broadening of the main focus. Phong and Kong (2016); Milovanović et al. (2015)
Finally, let us briefly discuss the relation between Veselago lensing in (chiral) metamaterials and in graphene. Whereas in metamaterials negative refraction typically occurs in a narrow frequency band around a resonance, in graphene it occurs for the full range of energies for which the Dirac equation is applicable, which means for energies until about 1 eV. Katsnelson (2013) Since, within the Dirac approximation, the classical Hamiltonians of the two valleys in graphene are equal, the classical trajectories in both valleys coincide. Therefore, if (pseudo)spin polarization leads to valley polarization, this has to happen because of quantum interference. The situation is quite different in chiral metamaterials, where the refractive index is different for left-handed and right-handed circularly polarized light. Tretyakov et al. (2003); Pendry (2004); Zhang et al. (2009); Xiong et al. (2010) Hence, the rays, which are the analogs of the classical trajectories, are different for both types of handedness. Furthermore, since one refractive index is typically negative, whilst the other one is positive, the classical rays are focussed for only one handedness and a well-defined polarization can be created.
Although graphene will be our main example in this paper, we stress that the behavior of its charge carriers is not unique. Another class of materials whose electrons follow the massless Dirac equation is formed by the two-dimensional surfaces of three-dimensional topological insulators. Hasan and Kane (2010); Qi and Zhang (2011); Moore (2010); Hasan and Moore (2011); Bansil et al. (2016) We are then dealing with real spin instead of pseudospin and using a spin-polarized STM one can inject a single spin. Therefore, we will generally refer to charge carriers governed by the Dirac equation as massless Dirac fermions and clearly indicate it when we specialize to the case of graphene.
The paper is organized in the following way: in section II, we discuss the basic equations that describe the Green’s function of an electronic Veselago lens for massless Dirac fermions. Subsequently, we discuss classical focussing and caustics in section III, and quantum interference and symmetry breaking in section IV. In section V, we discuss the semiclassical evaluation of the Green’s function, and compare various approximations with the exact solution. A semiclassical derivation of the vertical displacement of the maximum that results from (pseudo)spin polarization is presented in section VI. The resulting formula is tested for the case of graphene. In section VII, we briefly consider the case where current enters a graphene sample through a narrow graphene lead. We successively discuss the wavefunction, symmetry breaking and the semiclassical evaluation of the wavefunction. Finally, we present our conclusions in section VIII.
II Veselago lenses
In this section, we introduce the equations that describe an electronic Veselago lens for massless Dirac fermions and review the key results from the literature. We only consider the case of the Green’s function here, postponing the case where an electronic current enters the sample from one of its sides to section VII. We split our considerations into three parts. In the first subsection, we define the Green’s function and introduce the proper dimensionless parameters. Subsequently, we briefly review the classical focussing that was discussed in Ref. Cheianov et al., 2007. Finally, we write down the wavefunction for a Veselago lens formed by an n-p junction with a (pseudo)spin polarized source.
II.1 The Green’s function and dimensionless parameters
The Hamiltonian for two-dimensional massless Dirac fermions is equal to Katsnelson (2013)
[TABLE]
where is the two-dimensional unit matrix, the two-dimensional vector consists of the Pauli matrices and is the momentum operator. All position vectors are two dimensional, i.e. . The function represents the potential to which the charge carriers are subject and the quantity is the Fermi velocity. For the specific case of graphene in the nearest neighbor approximation, it is defined by , where eV is the hopping parameter and nm is the distance between two carbon atoms. Katsnelson (2013)
The Green’s function for the Hamiltonian (1) is defined by
[TABLE]
where indicates the source from which the particles are emitted with energy . For an arbitrary electron source , the equation of motion reads
[TABLE]
and the solution is given in terms of the Green’s function as
[TABLE]
In most of this paper, we will assume that we are dealing with a point source that has a certain polarization, which is pseudospin (sublattice) for the case of graphene Katsnelson (2013) and true spin for the case of the two-dimensional surfaces of three-dimensional topological insulators, Qi and Zhang (2011) i.e.
[TABLE]
For convenience, we will assume that the constants are dimensionless and that they form a vector that is normalized, i.e. . In practice, depending on the normalization of the source, these constants will however have a dimensionality, which can easily be incorporated into the description. It should also be noted that in the case of graphene one cannot use the continuum approximation when atomically sharp features are present. Hence, the notion of a point source implies that the diameter of the source is much larger than the interatomic distance , yet much smaller than the electronic wavelength : . Inserting the source (5) into Eq. (4), we obtain the wavefunction for our problem as
[TABLE]
Since we want to perform a semiclassical analysis later on, for which we need to know the true, dimensionless, semiclassical parameter, we restate the problem and its solution in dimensionless parameters. The intrinsic length scale of the problem is given by the distance from the point source to the junction, which will be introduced more precisely in the next subsection. Furthermore, let be the typical energy scale of the problem, which we take to be in this paper. Alternatively, one could set , with the typical value of , without any essential difference. These definitions allow us to define the dimensionless small parameter , and the dimensionless quantities , , and . Furthermore, we define . Taking into account that , we find that Eq. (2) becomes
[TABLE]
Defining and , we find that Eqs. (3)-(6) remain valid when we replace all quantities by their dimensionless counterparts. Therefore, the source and the wavefunction are given by
[TABLE]
In the following sections, we will work almost exclusively with these redefined (dimensionless) quantities and omit the tildes. Unless explicitly stated, we will always be referring to the dimensionless quantities defined here, rather than their original counterparts.
Briefly returning to the case of graphene, we remark that the Hamiltonian (1) is only valid near one of the two conical points in the Brillouin zone, namely at the -point. Katsnelson (2013) Near the other conical point, the so-called -point, the Hamiltonian reads
[TABLE]
Therefore, the Green’s function near the -point is related to the Green’s function near the -point by
[TABLE]
II.2 Classical focussing
Before considering the focussing of the electrons, let us first consider the classical motion of massless Dirac fermions. The matrix Hamiltonian
[TABLE]
describes both electrons () and holes () within the same equation. One can extract the classical Hamiltonian functions that correspond to this matrix Hamiltonian by replacing the momentum operators by -numbers and computing the eigenvalues. Berlyand and S. Yu. Dobrokhotov (1987); Belov et al. (2006); Tudorovskiy et al. (2012) We obtain two classical Hamiltonian functions,
[TABLE]
corresponding to electrons () and holes (). This readily shows that the group velocity of electrons is parallel to their momentum, , whereas for holes the group velocity is opposite to the momentum, .
Let us now, following Ref. Cheianov et al., 2007, consider electrons emitted by a point source at position incident on a one-dimensional n-p junction. We assume that the potential consists of a single step at :
[TABLE]
where is the Heaviside step function and . As before, in practice this means that the length scale of the potential increase satisfies . Then it is clear that the characteristic length scale of the system, which we used to define dimensionless parameters in the previous section, is equal to .
Now consider an electron incident on this potential from the left under an angle , with momentum , . At the interface, part of this electron is reflected, whilst another part is transmitted, with momentum . Since we consider scattering to right-moving holes, and , we have for . As the potential does not depend on , the transversal momentum is conserved and we find the relation
[TABLE]
which is nothing but Snell’s law for an electronic system. Cheianov et al. (2007) However, a very important characteristic of this system is that the refractive index is negative, which means that and have opposite signs. Therefore, the junction has the ability to focus electrons emitted by a source on the left-hand side, as can be seen in Fig. 1. We will discuss this focussing in more detail in section III.
Finally, we note that the maximal angle under which electrons can be classically transmitted is . For , this means that electrons that are incident on the barrier under an angle larger than
[TABLE]
will not be transmitted. This is related to the concept of a boundary angle in optics. For , all electrons that are incident on the boundary can be classically transmitted, and we can set .
II.3 The wavefunction for a polarized source
Now that we have reviewed classical focussing by an n-p junction, let us consider the wavefunction induced by the source (8). In appendix A, we solve Eq. (7) and obtain the Green’s function (91). Combining this result with Eq. (8), we find that the wavefunction induced by a (pseudo)spin polarized source equals
[TABLE]
where
[TABLE]
is the classical action. The limits of integration in Eq. (16) are determined by , where was defined in the previous subsection.
Of course, one can also use different source terms than (8), as was done in Ref. Cheianov et al., 2007.
III Caustics
In the limit where the dimensionless parameter , which we introduced in section II.1, is small, the main contribution to the integral in the wavefunction (16) is given by the stationary points of the action, Fedoryuk (1977); Guillemin and Sternberg (1977); Maslov and Fedoryuk (1981) i.e. the points where vanishes. This means that the main contribution is given by the points that are on the classical trajectories of the system. Arnold (1989) We find that they are given by
[TABLE]
Naturally, these are equivalent to the trajectories that were obtained before in Ref. Cheianov et al., 2007, as reviewed in section II.2.
There are also singular points, at which the second derivative vanishes. These points form a curve that separates the region where each point lies on three trajectories (and hence interference takes place) from the region where each point lies on a single trajectory, as can be seen in Fig. 1. Focussing takes places on such curves, which are known as caustics. Berry and Upstill (1980); Poston and Stewart (1978); Arnold et al. (1982); Arnold (1975) Some calculus yields that these points are defined by
[TABLE]
with the corresponding -value given by Eq. (18). Alternatively, Eqs. (19) and (18) can be cast into the form Cheianov et al. (2007)
[TABLE]
We can also look at the caustic from the point of view of the trajectories. If we parametrize them as , then the caustic is the set of points where the Jacobian vanishes. Indeed, some algebra shows that the Jacobian is proportional to the second derivative of the action:
[TABLE]
Hence the second derivative vanishes if and only if the Jacobian does.
Let us now consider the caustic in somewhat more detail. We first note that the transformation that sends to reflects a trajectory in the -axis, which implies that the set of trajectories is symmetric with respect to the line . Therefore, the caustic should be symmetric with respect to the -axis as well. Indeed, we see that is invariant under reflection of . Alternatively, we can also see directly from Eq. (20) that the caustic is symmetric.
Second, let us consider the shape of the caustic. For general , it consists of two so-called fold lines Berry and Upstill (1980); Poston and Stewart (1978); Arnold et al. (1982) meeting into a cusp point, see Fig. 1(a). From the symmetry considerations presented above, we conclude that this cusp has to lie on the -axis and therefore corresponds to . Equation (19) then implies that it is located at .
According to catastrophe theory, Arnold et al. (1982); Poston and Stewart (1978) a smooth change of variables can bring the action near the caustic into a certain normal form, which is a polynomial with its degree depending on the type of caustic. For points on the fold lines, the third derivative of does not vanish, and this normal form is a third order polynomial without a quadratic term. Arnold et al. (1982); Poston and Stewart (1978) In the Arnold classification, Arnold et al. (1982); Arnold (1975) this type of caustic is denoted by . At the cusp point, denoted by in the Arnold classification, the third derivative vanishes as well, but the fourth derivative is nonzero. It turns out that we can therefore express the action near this point as a fourth order polynomial without cubic term. In Fig. 1, we see that we can have two types of cusp catastrophes, depending on the value of the potential . For , we see that the cusp is the rightmost point of the caustic, whereas for , it is the leftmost point. The difference between the two types is the sign of the fourth derivative, which carries over to a plus or minus one in front of the quartic term of the normal form. For , this is a plus one, for , this is a minus one.
The theory of Lagrangian singularities Arnold et al. (1982); Poston and Stewart (1978) shows that in Hamiltonian systems in two dimensions the only generic singularities that can occur are folds and cusps. According to this theory, any other singularity will turn into one of these cases when an arbitarily small change is made to the system. However, the system that we are considering has an additional parameter that can be tuned, namely the potential strength . As we have seen, we can change the sign of the fourth derivative of the action from positive to negative by changing the potential. In doing so, we will inevitably pass through the point where the fourth derivative vanishes, and hence through a higher order singularity. By symmetry, this higher order singularity is again located on the -axis, and therefore corresponds to vanishing . At , the action (17) is an even function of , which means that its Taylor expansion in only contains terms of even order. In a generic setting, we can only expect the coefficients in front of the quadratic and the quartic terms in the expansion to vanish at this higher order singularity, since we have just two parameters, and . This would imply a singularity corresponding to a sixth order polynomial, i.e. a two-dimensional section of the so-called butterfly catastrophe . Arnold et al. (1982); Poston and Stewart (1978) However, looking at the action (17), we see that when , and not only the second and the fourth derivative vanish, but that in fact all derivatives of with respect to vanish. In this very special case, the n-p junction acts as an ideal lens and focusses all trajectories in a single point, as shown in Ref. Cheianov et al., 2007 and depicted in Fig. 1.
We wish to emphasize that this behavior is not generic, and is a special feature of the system under consideration. In fact, arbitarily small changes to the spatial setup, such as a non-straight barrier interface, or arbitarily small changes to the dispersion will ruin the perfect focus. In real graphene samples, we expect at least two corrections to the Hamiltonian (11) to contribute to the breaking down of the perfect focus. The first of these is the presence of next-nearest neighbor hopping, Katsnelson (2013) which will slightly change the classical trajectories of the system. Furthermore, it destroys Klein tunneling, although its influence on the transmission through an n-p junction was shown to be small. Kretinin et al. (2013) The second important correction to the Hamiltonian is trigonal warping, Katsnelson (2013) the influence of which will become stronger as the energy increases. As with next-nearest neighbor hopping, trigonal warping will change the classical trajectories of the system. Furthermore, it also destroys Klein tunneling for almost all orientations. Ando et al. (1998); Logemann et al. (2015)
IV Quantum interference and symmetry breaking
In the previous section, we saw that the classical trajectories and the caustic are symmetric with respect to the -axis. Let us now consider the symmetry of the Green’s function. First, we note that the classical action (17) satisfies . Then, making the change of variables in the integral, it is easy to show that
[TABLE]
Now let us consider the wavefunction (8) induced by a (pseudo)spin polarized point source. For its norm, , we obtain the equality
[TABLE]
This equals only whenever . Therefore, the wavefunction will in general not be symmetric, even though the classical trajectories are.
In Fig. 2, we have plotted the density , given by Eq. (16), near the cusp point for three different polarizations. For the polarizations and , the intensity is symmetric about the -axis, in accordance with what we just showed. For , the symmetry is broken and we see that the maximum of the wavefunction is displaced. This shift is due to quantum interference and is an effect of the (pseudo)spin polarization of the source. In Fig. 3, we show sections of the wavefunction along a line parallel to the -axis and through for various polarizations. We see that as the ratio decreases, the position of the maximum shifts more and more towards negative , while the intensity at the maximum decreases. When , the situation is once again symmetric, but the main focus has disappeared completely. Therefore, we conclude that we can markedly change the position of and the intensity at the central focus by changing the polarization.
Briefly returning to the case of graphene, we see from Eqs. (10), (22) and (8) that
[TABLE]
which means that the densities for the two valleys in graphene are related to each other by a reflection in the -axis. In particular, changes sign, which means that the maxima for the two valleys are on opposite sides of the -axis. In the following two sections, we investigate how large this asymmetry is and whether this may provide another way of realizing a valley filter in graphene.
V Semiclassical evaluation of the wavefunction
To gain a better understanding of this asymmetry and the factors that influence it, we investigate the wavefunction (16) with the semiclassical approximation. This will also give us more insight in the intensity at the central focus and in the way the size of the focus scales.
Central to the semiclassical approximation is the dimensionless small parameter that we introduced in section II.1. In section III, we already saw that in the limit , the main contribution to the integral (16) is given by the stationary points of the action, which give rise to the classical trajectories (hence the name semiclassical approximation, as we are in a situation that is ‘almost’ classical). In this limit, we can expand the wavefunction (16) as an asymptotic series in powers of .
In the simplest case, we consider points that are not on the caustic, which means that does not vanish at any of the stationary points . Such stationary points are called nondegenerate. Looking at Fig. 1, we see that we can distinguish two regions. In the first region, each point lies on a single trajectory, and hence the action only has one stationary point. In the second region, each point lies on three trajectories, and the action has three stationary points. In appendix B.1, we discuss how the leading order contribution of a nondegenerate stationary point to the integral (16) can be obtained by the conventional stationary phase approximation, Fedoryuk (1977); Guillemin and Sternberg (1977); Maslov and Fedoryuk (1981) with the result given by Eq. (96). In the first region, this directly gives us the leading order term of the wavefunction. In the second region, we need to compute the contribution of each of the three stationary points, and then add these contributions to find the correct approximation to the wavefunction (16). We will henceforth refer to these results as the Wentzel-Kramers-Brillouin (WKB) approximation.
In section III, we discussed that at the caustic the second derivative vanishes. Therefore, the result (96) diverges and we need to obtain the main contribution to the integral (16) in a different way. In appendix B, we show that the simplest approximation for the wavefunction near a caustic can be obtained by making a Taylor expansion of the action in up to the first nonvanishing term.
For the fold caustic, discussed in appendix B.2, this means that we have to expand up to third order, from which one obtains an expression in terms of the Airy function, Airy (1838) see also e.g. Ref. Connor and Farrelly, 1981a. The final result, presented in Eq. (105), is valid in an neighborhood of the fold. We remark that expression (105) was derived under the assumption that the limits of integration are infinite, whereas in Eq. (16) they are finite. This is, however, not a problem, since the main contribution to the integral comes from a narrow vicinity of the stationary points. Fedoryuk (1977); Guillemin and Sternberg (1977); Maslov and Fedoryuk (1981) Since all of the latter lie between the finite limits of integration in the integral, we can extend the limits of integration to infinity without changing the leading order term. Furthermore, we note that when we expand the action to even higher orders near the fold caustic, we will only get corrections beyond the leading order, i.e. terms in higher powers of . Finally, a more accurate result can be obtained by using the uniform Airy approximation, Chester et al. (1957) see also e.g. Ref. Connor and Farrelly, 1981a, but we will not consider this approximation in this paper.
The leading order approximation to Eq. (16) for a point near the fold caustic then consists of two terms. The first term is the one with the Airy function that we just discussed. The second term is a WKB term that comes from the third stationary point. In terms of the trajectories plotted in Fig. 1, this term originates from the trajectory that is not tangent to the caustic near the point , but rather “crosses” the caustic. We henceforth refer to the sum of these two terms as the Airy approximation.
Since our main interest in this paper is the asymmetry that is induced near the main focus, we now concentrate on the wavefunction near the cusp. Since the third derivative vanishes at the cusp, one has to expand the action up to fourth order. In appendix B.3.1, we review how this leads to an expression for the wavefunction near the cusp caustic that involves the Pearcey function, Pearcey (1946); Connor and Farrelly (1981b); Connor and Curtis (1982) see also e.g. Ref. Connor and Farrelly, 1981a, which is defined in Eq. (108). The result, presented in Eq. (116), contains the coefficients and , defined in Eqs. (107) and (99), which can be obtained by taking derivatives of the action . As we already saw in section III, the cusp corresponds to , which considerably simplifies the calculations. After some calculus, we find that the nonzero coefficients , which are the -th derivatives of the action at the cusp point, are given by
[TABLE]
As we already discussed in section III, we see that is positive for , and negative for . Furthermore, we obtain the coefficients as
[TABLE]
Comparing Eqs. (16) and (93), we see that the amplitude does not depend on , and that
[TABLE]
Finally, combining the above results with the general result (116), we find that the leading order approximation of the wavefunction in an neighborhood of the cusp is given by
[TABLE]
The first thing that should be noted about the result (28) is that, regardless of the polarization, it is symmetric with respect to the -axis, because of the fact that the Pearcey function is even in its second argument. Therefore, this approximation is insufficient if we want to understand the asymmetry. Second, we note that this approximation is not valid when the the potential is equal or close to , i.e. when we are close to the ideal focus, since in that case the coefficient vanishes or becomes very small, and the result (28) diverges.
Before we take a closer look at the asymmetry, we first want to see how well the approximation (28) works for the symmetric polarization . To this end, we compare it with the exact wavefunction, Eq. (16), which is evaluated by numerical integration. We also compare it with the result of the uniform approximation, Ursell (1972); Connor and Farrelly (1981a) which is discussed in appendix B.4. In this approximation, we do not perform a Taylor expansion of the action, but instead bring the action to its normal form near the cusp by an exact change of variables. The final result, shown in Eq. (128), is given as a sum of the Pearcey function and its derivatives. In order to make the comparison complete, we also include the WKB approximation and the Airy approximation that we discussed before. These are not expected to work well in the vicinity of the cusp.
In Fig. 4, we compare these five approximations for three different values of the small parameter . We see that the uniform approximation perfectly coincides with the exact wavefunction (16). Furthermore, we observe that the Pearcey approximation is not very accurate for large values of , but becomes much better when we decrease . Note that although it typically overestimates the magnitude of the wavefunction, it correctly predicts the position of the maximum for all three values of . This implies that we may be able to find the position of the asymmetry by including higher order corrections, even for rather large . As predicted, the WKB approximation works well far away from the cusp, but diverges as we come close to it, as does the Airy approximation. Near the fold caustic, the Airy approximation performs well for a large range of distances and for all three values of , whereas the WKB approximation only works far away from the fold.
As we just saw, the leading order Pearcey approximation (28) is not enough to reproduce the asymmetry that we found in section IV. Therefore, let us look at higher order corrections to the Pearcey approximation, which are discussed in appendix B.3.2. The corrections can come from two different sources, namely from higher order terms in the Taylor expansion of the action and from higher order terms in the expansion of the amplitude, i.e. the part of the integrand in Eq. (16) that precedes the exponent with the action. As we discussed in section III, the cusp point lies on the line , which implies that the action (17) is symmetric with respect to . Therefore, all terms in its Taylor expansion that are odd with respect to vanish at the cusp, and in particular the fifth order term vanishes. This means that the second term of in Eq. (119) is irrelevant, as . Hence the only correction of is given by Eq. (120). Using the results (25), (26) and (27), we obtain the first-order correction to the leading order term (28) as
[TABLE]
where represents the derivative of the Pearcey function with respect to its second argument, as defined in Eq. (121). Since is odd in its second argument, the sum of the leading order term (28) and the first-order correction (29) does not necessarily have its maximum at .
In Fig. 5, we compare the Pearcey approximation including the first-order correction with the exact solution, the uniform approximation and the Airy and WKB approximations on the line that goes through the cusp point and is parallel to the -axis. Comparing Fig. 5 with the middle panels of Fig. 4, we see that for polarization the result does not qualitatively differ from the leading order approximation, although the numerical values are slightly different. For polarization , we see that with the first-order correction (29) we correctly reproduce the position of the maximum, even though it is no longer at . This holds for both the large and the small value of .
When the polarization equals , we see from Eq. (28) that the term proportional to vanishes. This makes sense, since we already saw in section IV that the central resonance vanishes in this case. Hence, the leading order term for this case is given by Eq. (29), which correctly reproduces the position of the two maxima that lie symmetrically on both sides of . However, since vanishes at , this approximation predicts that the wavefunction also vanishes on the -axis, which is incorrect, as can be seen from Fig. 5(c). Therefore, we have also included the second correction, which is of , in the Pearcey approximation plotted in Fig. 5(c). Looking at Eq. (119), and remembering that both and vanish in our case, we easily see that this correction is given by Eq. (123). Taking the results (25), (26) and (27) into account, we obtain the second order correction as
[TABLE]
where is the derivative of the Pearcey function with respect to its first argument, as defined in Eq. (122). With this second-order correction, we see from Fig. 5(c) that we have a reasonable approximation for the value of the wavefunction at . We remark that this correction does not substantially influence our prediction for the position of the maximum for this polarization.
Figure 5 clearly shows that we can greatly decrease the intensity at the central focus by changing the polarization from to , as discussed in the previous section. Let us now use the Pearcey approximation that we have developed to derive an equation for the ratio between the intensities for these two polarizations. In order to arrive at a simple expression, we will use the value of the various Pearcey approximations at the cusp point . Although this is not the position of the main focus, the wavefunction at this point gives us a good indication of its value at the maximum. In Fig. 4, we see that for the largest , the exact value of the maximal intensity is approximately equal to the value of the Pearcey approximation at the cusp, due the fact that the latter gives too large values. For the smallest , we see that the maximal intensity is about a factor of two larger than the Pearcey approximation at the cusp. For the polarization , we use the leading order Pearcey approximation, given by Eq. (28). At the origin, the Pearcey function takes a particularly simple form, as , where is the gamma function. Korn and Korn (1961) This identity can easily be proven directly using the definition (108) of the Pearcey function and the definition of the gamma function. Since we also have , see also Eq. (109), we obtain
[TABLE]
For the polarization we use the second order correction (30), since both the leading order term (28) and the first order correction (29) vanish on the -axis. One can show that and . Therefore, we find that
[TABLE]
For the ratio between the two intensities at the cusp point, we then obtain
[TABLE]
which shows that the relative decrease of the intensity at the main focus is proportional to the small semiclassical parameter .
We finish this section by showing the densities that the various approximation schemes give near the cusp point for two different values of . Looking at the comparisons in Fig. 4(a), we see that at it is difficult to construct a global approximation for by combining the different local approximations, since there is no region where the Pearcey approximation smoothly joins the stationary phase approximation. Therefore, we conclude that for a global approximation only the uniform approximation is adequate. In Fig. 6, we show the results of this approximation . Comparing Figs. 2 and 6, we see that the agreement is excellent, as we already inferred from the comparisons along the various sections.
A few words about the implementation of the uniform approximation are in place here. In the region where each point lies on three trajectories, the equation has three real roots . The values of these roots are restricted, since for we have and for we have . We can obtain these roots numerically, and subsequently determine the action , its second derivative and the amplitude at these points, from which we can obtain the parameters for the uniform approximation, as explained in detail in appendix B.4. When only lies on a single trajectory, the equation still has three roots, but this time only one of them is real and two of them are complex. However, the absolute value of these complex roots is not necessarily restricted. We have found that when the complex roots become too large in absolute value, the performance of the uniform approximation becomes rather poor. When we impose on the same demands that hold for the case when is real, i.e. when and when , we obtain good agreement. However, this means that we cannot use the uniform approximation far away from the caustic, which gives rise to the large white area in Fig. 6. Note in particular the strange situation that occurs for , and , where we have three real roots, two of which have an absolute value larger than .
Coming back to the densities that the various approximations predict near the cusp point, we see that for , it is possible to construct a global approximation for by combining various local approximations. From the comparisons in Fig. 4, we see that for polarization , we can use the Pearcey approximation in an area around the cusp that is described by an ellipse with semi-major axis (along the -direction) and semi-minor axis (along the -direction). Furthermore, we can use the Airy approximation along a distance from the fold caustic in the outward direction (where there is only a single trajectory) and along a distance in the inward direction (where there are three trajectories). Outside of both these regions, we can use the WKB approximation. We remark that we do not need to patch the different approximations together by determining certain constants, since all of the different approximations are simplifications of the same wavefunction (16) that are appropriate for a certain region.
In Fig. 7, we show the combination of the various approximations, as well as the region for each of the approximations. We see that the final result nicely coincides with the exact wavefunction (16), which was evaluated numerically. Since the result of the uniform approximation perfectly coincides with the exact wavefunction (16), it is not shown separately.
VI Derivation of the displacement from semiclassical considerations
In the previous section, we saw that, for a certain set of parameters, we could reproduce the vertical position of the maximum, even though it was displaced from the -axis, by using the first order correction (29) to the leading order term (28). In this section, we perform a more systematic study of the vertical displacement of the maximum, which is caused by the (pseudo)spin polarization, and obtain a simple formula for the shift.
Let us therefore try to derive a formula for the -coordinate of the maximum when equals . To this end, we consider the sum of the leading order Pearcey approximation (28) and its first correction (29), i.e. , at . In order to find the maximum, we need to find the points where the first derivative vanishes. Unfortunately, this equation cannot easily be solved, as it involves the Pearcey function and one of its partial derivatives. Therefore, let us approximate the Pearcey function by its Taylor expansion. First of all, we note that it is even in its second argument, see Eq. (109), which means that when we perform a Taylor expansion in the second argument, all terms of odd order vanish. In the previous section, we already determined the zeroth order coefficient of the Taylor expansion, which is equal to . Furthermore, from the definitions (108) and (122), it is easy to see that . Using the result for from the previous section, we can then obtain the second order coefficient . Combining our results, we find that
[TABLE]
where the last equality can be obtained by direct computation or by using the differential equation that is satisfied by the Pearcey function, see Ref. Connor and Farrelly, 1981b. From these equalities, we see that the fourth order coefficient is much smaller than both and . Therefore, we obtain a rather accurate approximation by replacing the Pearcey function by its second order Taylor expansion, that is,
[TABLE]
where the second relation is a consequence of the symmetries of the Pearcey function, see Eq. (109). We can use the same approximation for the derivative, which naturally gives
[TABLE]
Using the approximations (35) and (36) for the Pearcey function and its derivative, we find, after some algebra, that
[TABLE]
where , see also the equivalent definition in appendix B.3.1, and the amplitude and its derivative , see Eq. (27), are to be evaluated at . We remark that this equation is valid for both and , since .
At this point we recall that the vertical displacement of the maximum is not a classical effect, since the classical focus lies at . Therefore, the -coordinate of the maximum cannot be of order unity. Instead, we expect to be of order , since the effect is caused by (quantum) interference. Indeed, we see that Eq. (37) is a cubic equation in , which can be solved analytically to find the maximum . However, we can also do one additional approximation, using the main assumption of the semiclassical approximation, which is that the semiclassical parameter is small. Compared to the first two terms, the last three terms in Eq. (37) have an additional factor of . Therefore, in a crude approximation, we can neglect them. The resulting linear equation in can easily be solved and we obtain
[TABLE]
Interestingly, this result does not depend on the coefficients and of the Taylor expansion. Returning from dimensionless units to regular units, we find that
[TABLE]
Note that the factors of that are present in both and have cancelled, yielding a result that does not depend on the length scale of the system.
In the remainder of this section, we compare several results for the position of the maximum. The first of these is the exact value, obtained by numerically determining the maximum of the wavefunction (16) at . The second is the value obtained by numerically determining the maximum of the Pearcey approximation, composed of the leading order term (28) and the first correction (29). The third and fourth results are the result (39) and the solution of the cubic equation (37), respectively. To gain understanding of the numbers involved, we now specialize to the case of graphene, for which , see section II.1.
Figure 8(a) shows the dependence of on the length scale of the system. We see that is indeed largely independent of length, as predicted by Eq. (39), with the exact solution showing only a slight variation. In Fig. 8(b), we consider the dependence of the ratio on length. For this purpose, we define the width of the peak as the full width at half maximum (FWHM) of . Looking at the leading-order Pearcey approximation (28), we expect the peak width to scale as , since both and contain a factor , and the coefficient does not depend on . In Fig. 8(b), we indeed see a clear power law scaling of , and from a fit we find a value close to for the power, as predicted. So although the position of the maximum roughly stays the same with increasing length, the displacement from the -axis will be harder to see because the width of the maximum increases. This effect can be clearly seen when comparing the middle panels of Fig. 5.
In Fig. 8(c) and (d), we consider the dependence of and on the electron energy. It is clear that our result (39) performs very well. In fact, all results are quite close to each other and show a clear power law behavior , with power . Looking at Eq. (28), we expect the width to scale as . Indeed, we see from Fig. 8(d) that shows a clear power law behavior, and from a fit we find that the power is close to . This implies that although we can increase by lowering the energy, this will also increase the width of the peak, yielding only a small increase in the ratio .
The dependence of on the potential is shown in Figs. 8(e)–(h). From the exact value, we see that is largely independent of the potential, being only slightly larger at the point . Interestingly, our result (39) outperforms the other two approximations for , with the difference becoming smaller as increases.
In Fig. 8(i), the dependence of on the polarization is shown. As already noted in section IV, the displacement from the -axis becomes larger as the ratio decreases. It is clearly seen that our approximation (39) and the solution of the cubic equation (37) give good results when is larger than approximately , but fail for smaller ratios. The reason for this is that for large values of , we can no longer approximate the Pearcey function by its second order Taylor expansion around . This is indicated by the fact that the Pearcey approximation, consisting of the leading order term (28) and the first correction (29), gives good results for all polarizations. Adding the second correction (30) does not substantially change the result. Note that when we invert the polarization, that is, when we consider in the range minus one to one, the position of the maximum changes sign with respect to Fig. 8. This can be seen directly from Eq. (23) and is particularly clear from our result (39).
Finally, we show the dependence of and on the coordinate in Fig. 8(j) and (k). We see that is roughly constant, and that the width varies somewhat. Interestingly, the exact solution and the Pearcey approximation follow slightly different trends and intersect around . The fact that the position of the maximum does not show a large variation with means that we can safely use our results to obtain an estimate for the asymmetry at the main focus, which is generally not located at the cusp point.
Looking at all the different dependencies in Fig. 8, we conclude that our approximation (39) gives quite accurate predictions for the position of the maximum, even though it was derived using several approximations. It only fails when the ratio comes close to minus one, which is the point where the central resonance disappears completely.
In section IV, we showed that, for the case of graphene, the displacement of the maximum in the -valley is opposite to the displacement in the -valley. Since the effect is rather large, on the order of a few nanometers for energies around 100 meV, we believe that it would be possible to observe it experimentally by measuring the spatial profile of the wavefunction with the help of an STM. Another possibility would be to try to place a tiny contact near the predicted maximum and to measure the valley composition of the current using the valley Hall effect. Xiao et al. (2007); Gorbachev et al. (2014) Since a typical laser beam is larger than a few nanometers in size, we believe that it would not be possible to measure the effect using second harmonic generation. Wehling et al. (2015)
An important remark is that, for typical energy and length scales, does not exceed 0.5, as can be seen in Fig. 8(b) and (d). So although the peak displacement is rather large, the peaks are also rather broad, making it much harder to identify them. This ratio improves as the energy and length of the device become smaller, although it should be noted that both the peak displacement and its width increase as the energy decreases. Because of the rather small value of , we do not think that the effect is large enough to create an effective valley filter in graphene.
VII Current entering from a lead
In this section, we no longer consider the Green’s function, but discuss the related problem where current flows into the sample through a lead on one of its sides. We first construct the wavefunction for the general case and subsequently specialize to the case of a graphene sample with graphene leads. In the second subsection, we consider the symmetries of the wavefunction. In the final subsection, we consider the semiclassical evaluation of the wavefunction.
VII.1 Derivation of the wavefunction
For definiteness, we henceforth assume that current enters the sample from the left, through a lead of width which is located between the points and . As before, we consider a sample with a potential that consists of a single step, see Eq. (13). This gives rise to a setup that is qualitatively similar to the Green’s function, but instead of considering a single point source, we now consider a lead with a finite width, which, following Huygens’ principle, Arnold (1989) may be regarded as a collection of point sources.
Defining the characteristic length scale of the system by , we can define the dimensionless quantities , , , and in the same way as in section II.1. Omitting tildes, this leads to the Hamiltonian (11). A new dimensionless parameter in the problem is the lead width, which is naturally defined as . Because of the translational symmetry of the lead, the wavefunction of each mode in the lead can naturally be decomposed as the product of a phase factor and a transversal wavefunction . We henceforth assume that the each of the modes is normalized in such a way that it carries unit current, which means that in dimensionless units
[TABLE]
In order for this equality to hold in units with dimensions as well, we set . As before, we omit the tildes from here on and deal exclusively with these newly defined quantities, unless otherwise indicated.
In mathematical terms, we can now formulate the problem at hand as an initial value problem, namely
[TABLE]
where the dimensionless equals minus one. In appendix C, we solve this problem for an arbitrary initial wavefunction . The general solution (140) is a linear combination of the independent solutions and , defined in Eqs. (71) and (72), which correspond to waves coming in from minus infinity and infinity, respectively. When we consider the case where no current flows into the sample from the right, the coefficient in front of should be zero, and it is sufficient to consider only the term proportional to . Note that in a realistic sample, the former coefficient is not necessarily zero, as the finite length and width of the sample will introduce scattering between various modes. Nevertheless, we expect the approximation to hold in reasonably sized samples, as the induced scattering will be small. In the appendix, we show that in that case the wavefunction is approximately given by
[TABLE]
where is a row vector and is the Fourier transform of , defined in Eq. (75).
In the remainder of this section, we consider the specific example of a graphene sample, with current entering through a graphene lead. In particular, we consider a graphene lead with zigzag edges, which do not mix the two valleys and . Within the continuum approximation, which is valid for sufficiently broad leads, we can then obtain the wavefunction in the lead by setting the boundary conditions Brey and Fertig (2006)
[TABLE]
which are valid for both valleys. Within the graphene lead, we allow for the presence of a constant mass, which can for instance arise in the context of chemical functionalization, Katsnelson (2013) or for graphene on a substrate, such as h–BN. Sachs et al. (2011); Bokdam et al. (2014); van Wijk et al. (2014) In the Hamiltonian (11) for the -valley, it manifests itself as an additional term . Solving the eigenvalue equation for a constant potential and a constant mass, and imposing the boundary conditions (43), we find that the wavefunction within the lead equals Brey and Fertig (2006)
[TABLE]
where is the so-called boxcar function:
[TABLE]
The momenta and are defined by the relations and
[TABLE]
Furthermore, the factor is defined by
[TABLE]
When , it is easy to show that , and that its value alternates between successive bands. Finally, the normalization factor ensures that the mode carries unit current, i.e. that Eq. (40) is satisfied. We remark that in the above computation we have disregarded the surface states, Brey and Fertig (2006) and hence do not consider very low energies. The computations for the -valley are entirely analogous.
The last issue that we need to consider is the relation between and . Let us first look at the case where and consider a single incoming mode. Since the lead and the left side of the sample have the same potential and the same mass, we can expect that there will be very little backreflection into the lead. Though small, this reflection will in reality however not be zero, because of the finite width of the lead. Furthermore, note that when we completely neglect the backreflection, the coefficient in front of no longer vanishes, as a computation using Eq. (140) shows. Nevertheless, this coefficient is still small, and we consider the approximation of by feasible.
When the mass inside the lead does not vanish, i.e. , the situation is rather different. In this case, the dispersion in the lead, differs from the dispersion in the sample, . We therefore expect significant backreflection into the lead, which increases as the mass increases. This means that we can no longer approximate by , but that we should include multiple left-moving modes with appropriate reflection coefficients.
VII.2 Symmetries of the wavefunction
Let us first consider the symmetry of the wavefunction in the absence of a mass term, i.e. for . In this case and it is easy to show that
[TABLE]
As we discussed in the previous subsection, we can approximate by in this case. When we consider its Fourier transform, we see that it has the same symmetry, i.e
[TABLE]
Using this identity, we can show that the wavefunction (42) also possesses this symmetry:
[TABLE]
The main ingredient of the calculation is the change of variables in the integral. Under this transformation, becomes , as can be seen from Eq. (71). We conclude from Eq. (50) that is symmetric in the -axis when .
Before considering the case of a nonzero mass, let us first consider a second symmetry of the lead wavefunction. From Eq. (44), we see that is real, irrespective of the mass as long as . Because of the properties of the Fourier transform (75), this means that
[TABLE]
This symmetry implies that , which, roughly speaking, means that the amount of current that has positive is the same as the amount of current that has negative . For , the latter equality is also implied by Eq. (50). However, the more general statement (51) is also true for nonzero masses.
When , the symmetry (48) is clearly broken, since one sees from Eqs. (44) and (47) that a mass term creates a difference in the amplitudes on the two sublattices. However, something more fundamental is going on when , since at the point where the lead and the sample join, the dispersion relation changes. As discussed in the previous subsection, we therefore expect significant backreflection into the lead and we cannot simply approximate by the incoming mode . Instead, we need to take a linear combination of the incoming mode and several reflected modes, with appropriate coefficients. These reflection coefficients are in general complex, meaning that is no longer a real function. Hence, the symmetry (51) is broken, and we can expect the amount of current that is emitted with positive transversal momentum to be different from the amount of current that is emitted with negative transversal momentum. Therefore, we expect the effect of sublattice polarization for this case to be quite different from the effect for the case of the Green’s function, which we discussed elaborately in the previous sections. Since the determination of the reflection coefficients is in general not an easy task, we do not pursue this problem further in this paper. However, from our previous considerations it is clear that a sublattice polarization, originating from a mass term within the lead, should lead to an asymmetry.
VII.3 Semiclassical evaluation
In this final subsection, we consider the semiclassical evaluation of the wavefunction (42) for a sample where the current enters from a lead with zero mass, i.e. is given by , Eq. (44), with . We make the dependence on the lead mode explicit by including the mode number in the notation, i.e. we write and .
Let us first consider the “deep” semiclassical limit, where both and the dimensionless parameter . When we want to apply the stationary phase approximation to the solution (42), we should be aware of the dependence of on . Explicitly writing down the Fourier transform, we obtain
[TABLE]
We insert this Fourier transform into the wavefunction (42) and specialize to the case , whence is given by Eq. (71). We then see that we are dealing with a sum of two two-dimensional integrals, that should be considered separately. The actions for these two integrals are given by
[TABLE]
The stationary points correspond to those points where the partial derivatives with respect to and vanish. The condition yields the condition
[TABLE]
which is very similar to the condition (18) that we had for the Green’s function. Furthermore, yields
[TABLE]
Together, these two conditions determine the classical trajectories of the system. We see that, as in the case of the Green’s function, the trajectories are straight lines and are focussed by the n-p junction. However, this time they are not emitted from a single point, but from a line, parametrized by the variable . This can be seen as an illustration of Huygens’ principle: each point of the lead acts as a point source. However, in this situation the transversal momenta of the trajectories are strongly constrained and can take only two values, namely , with the transversal momentum of the mode in the lead. We also note that, because of this constraint, the set of trajectories covers only a limited region of space.
Caustics in the system arise when the Hessian matrix , the matrix of second derivatives of the action, is degenerate, i.e. , see appendix B.1. Computing the second derivatives of the action (53), we see that and that . Hence , and we always have one positive and one negative eigenvalue. We therefore conclude that, as long as we are in the deep semiclassical limit, there are no caustics in the system.
In this limit, we can therefore construct an approximation for the wavefunction, given by Eqs. (42) and (52), by employing the WKB approximation, as explained in appendix B.1. The calculation is rather involved, and in particular requires a careful analysis of the transversal momenta , defined by Eq. (46). Numbering the modes from the lowest value of up and starting at one, one can show that . The final result is then given by
[TABLE]
where is the transmission coefficient (73). We note that is a function of , and (which equals ) through the condition (54) for a stationary point.
From the equations for the stationary point, we see that the two trajectories that emerge from the point meet each other in the point . Around this point there is a rhombus shaped region where interference occurs. Furthermore, we note that the point is different for each mode , unless , in which case the point is the same for all modes, as discussed in section III. We also remark that the WKB approximation (56) has the symmetry (50). This can be easily seen when one uses the identity . Note that this symmetry implies that within the interference region, both components of the wavefunction are given by a cosine for modes with , whereas they are given by a sine for modes with .
Secondly, let us consider the case where is small, and is rather large, which means that is small, i.e. we are dealing with a relatively narrow lead. In that case, the physical situation is slightly different from the one sketched in the previous paragraphs. In order to understandy why, one needs to consider the width of the Fourier transform of the modes in the lead. To determine it, we first note that each of the components of the wavefunction is the product of a trigonometric function and a boxcar function. Therefore, the Fourier transform of each of these components is the convolution of the Fourier transforms of the two functions that make up the product. Since the Fourier transform of a trigonometric function is the sum of two delta functions, it is the Fourier transform of the boxcar function that determines the width of the Fourier transform of the lead wavefunction. This Fourier transform is easy to compute, and we obtain
[TABLE]
To obtain an estimate of the full width at half maximum (FWHM) of this Fourier transform, we expand the sine up to third order in its argument. Solving for the point where the absolute value of the function is half of its maximal value and subtracting the two solutions, we find that the FWHM is approximately given by
[TABLE]
where all units are dimensionless. Going back to units with dimensions, we see that the width is determined by the dimensionless parameter .
When this dimensionless semiclassical parameter in the lead is rather large, the Fourier transform is broad. In that case, it is less appropriate to consider the solution (42) as a double integral to which one should apply the WKB approximation. Instead, one should rather think of it as a single integral and consider as a function that does not depend on . From a physical perspective, one might say that, from each point in the lead, trajectories come out at all angles, instead of at just two. This means that we enter a regime for which the classical picture is qualitatively more similar to the one discussed in section III, albeit with a lead as the source of electrons instead of a single point. In particular, we expect the formation of caustics when .
To illustrate this, let us consider a particular situation. In Fig. 9, we show the norm of the Fourier transform of the wavefunction for various modes in a graphene lead of width nm. In Fig. 9(a), the energy of the electrons is meV, which means that . One indeed sees that the Fourier transform of the only mode in the lead is very broad and that it does not make much sense to speak about characteristic momenta. When we raise the energy of the electrons to meV, the dimensionless parameter and there are six modes in the lead. Now the Fourier transform of the various modes is much narrower, as can be seen in Fig. 9(b), and it does make sense to speak about characteristic momenta.
In Fig. 10, we show the norm of the wavefunction (42) for various electron energies and potential heights for fixed length scales nm and nm. Comparing Figs. 10(a) and (b), we see that for meV, we are indeed in a situation that is qualitatively similar to the Green’s function, since trajectories come out of the lead at all angles. When , they are focussed in a single point, and when , a caustic occurs and we see the characteristic interference pattern. When meV, this pattern is already much less pronounced, see Fig. 10(c). For an electron energy meV, the various modes in the lead carry designated momenta, as is particularly clear from Fig. 10(e), where the wavefunction that results from the second mode in the lead is shown. For , Fig. 10(d), all trajectories are once again focussed in a single point, whereas for this is not the case. However, in Fig. 10(f), where we show the intensity , i.e. the sum of the intensities that result from the separate modes in the lead, we do not see a clear interference pattern characteristic of a caustic. Instead, we have a rather sharp focussing spot, indicating that we are in the regime where we can approximate the wavefunction (42) by its WKB approximation (56).
Comparing Fig. 10(f) with the numerical results from Ref. Milovanović et al., 2015, we see that there is qualitative agreement between the two approaches. Unfortunately, their numerical simulations use rather high energies, typically 0.4 eV, in combination with rather wide leads, typically 50 nm, so their pictures only show one of the two regimes that we have identified, i.e. the WKB regime. However, going to lower energies or narrower leads, we believe that it should be straightforward to observe the other regime, in which caustics occur, as well. In their paper, they also observe a lowered transmission for leads with zigzag edges as compared to leads with armchair edges, which they attribute to a poorer lensing ability of an armchair p-n interface as compared to a zigzag p-n interface. As expected, our study, which stays within the framework of the continuum approximation, does not offer any alternative explanations for this effect and we believe that further research would be necessary to elucidate its nature.
In Fig. 11, we compare the total density for the wavefunction (42) and its WKB approximation (56) for meV, nm and nm. We see that the WKB-approximation captures the essential behavior of the wavefunction, but that the discrepancy is rather large, especially away from the maximum. We ascribe this discrepancy, which is notably smaller for than for , to the fact that the width of the Fourier transform of the lead wavefunction is still rather large. For meV and nm, our estimate (58) gives in dimensionless units, which means that for each mode, apart from , there are still a lot of other values of that contribute to the scattering.
VIII Conclusion
In this paper, we have studied two realizations of electronic Veselago lenses for massless Dirac fermions. We have found that in both cases the presence of (pseudo)spin polarization leads to symmetry breaking. By comparing the exact solutions with various semiclassical approximations, we have established that the semiclassical approximation is an effective tool to study focussing in graphene.
For the case of the Green’s function, we have demonstrated that, depending on the (pseudo)spin polarization, the main focus can either be vertically displaced, or can vanish completely. When the polarization equals , the main focus is displaced from the -axis, on which it lies in the absence of (pseudo)spin polarization. Specializing to the case of graphene, the size of this effect is typically on the order of several nanometers and the effect is opposite for electrons in the -valley and in the -valley. However, since the ratio of the displacement and the peak width is typically smaller than 0.25, we believe that this effect is not strong enough to create an effective valley filter in graphene. Nonetheless, we think that the effect could be measured experimentally. An initial sublattice polarization of could be realized by injecting electrons onto a single site with an STM, and a smaller amount of symmetry breaking could perhaps be attained by considering tunneling through hexagonal boron nitride. Sachs et al. (2011); Bokdam et al. (2014); van Wijk et al. (2014) Subsequently, one could for instance try to measure the displacement using an STM, or with the valley Hall effect. Xiao et al. (2007); Gorbachev et al. (2014) When the polarization equals , the main focus vanishes completely. Although we believe that such a polarization would be hard to realize in graphene, it is likely to be attainable in topological insulators, where we are dealing with real spin instead of pseudospin.
To gain more insight into the effect of different polarizations, we have studied the Green’s function using various semiclassical approximations. We have demonstrated that the vertical position of the main focus can be well predicted using the Pearcey approximation, provided that we include its first correction. This approximation also gives good results for the horizontal position of the maximum, even for large values of the semiclassical parameter. Furthermore, we have shown that the uniform approximation shows very good agreement with the exact solution, making it the preferred approximation when one is not only interested in the position of the maximum, but also in its value. Using the Pearcey approximation with various corrections, we have derived Eq. (33), which shows that the ratio between the peak intensities for polarizations and is proportional to the dimensionless semiclassical parameter of the system. Finally, we have derived Eq. (39), which reveals how the displacement of the main focus depends on the different system parameters. We have demonstrated that it shows excellent agreement with the exact solution.
For the case of current entering a graphene sample through a narrow graphene lead, we have used the semiclassical approximation to identify two different regimes. When the dimensionless semiclassical parameter in the lead is rather large, while the semiclassical parameter of the system is small, which happens for instance for low energies or narrow leads, we expect caustics to be formed in the system. On the other hand, when both semiclassical parameters are small, we have shown that a rather sharp focussing spot will occur. We believe that the transition between these two regimes should be visible both in experiment and in numerical simulations. The effects of symmetry breaking in this system are less clear, since a mass term in the lead not only breaks the reflection symmetry in the -axis, but also the symmetry of the Fourier transform of the total wavefunction in the lead, due to the presence of reflected waves. Because of this, we can generally expect the amount of current with positive transversal momentum to differ from the amount of current with negative transversal momentum.
Acknowledgements
We are grateful to Timur Tudorovskiy for stimulating discussions that initiated this work, and to Sergey Dobrokhotov for valuable discussions on the semiclassical evaluation of the wavefunction. Furthermore, we are grateful to Kostya Novoselov, Alex Khajetoorians, Misha Titov and Erik van Loon for helpful discussions about the project.
The authors acknowledge support from the ERC Advanced Grant 338957 FEMTO/NANO and from the NWO via the Spinoza Prize.
Appendix A Green’s function
In this appendix, we derive the Green’s function defined by Eq. (7) with the potential (13). In the first subsection, we define the scattering states for the problem without a source term, which we will need in the construction later on. In appendix A.2, we present a rigorous derivation of the Green’s function, based on the method of variation of constants.
A.1 Scattering states
Let us consider the eigenstates of the matrix Hamiltonian (11). Since the potential step (13) is independent of , the transversal momentum is constant. Therefore, we can define . In the electron region (), we define the right-moving and left-moving states by
[TABLE]
where is defined by , . In the hole region (), we define the right-moving and left-moving states by
[TABLE]
where is defined by , . After some calculus, we find that for both right-moving states, the probability current equals 1, whereas for the left-moving states it equals -1. Furthermore, we find that vanishes.
With these definitions, we can define two independent solutions for scattering by a potential step. One solution incoming from the left,
[TABLE]
and one solution incoming from the right,
[TABLE]
Matching each of these solutions at the barrier interface gives
[TABLE]
Note that the conservation of probability current automatically ensures that .
A.2 Derivation of the Green’s function
In order to obtain the solution to Eq. (7), we need the Fourier transform and its inverse, which we define as
[TABLE]
Performing the Fourier transform with respect to , we obtain
[TABLE]
We now construct the solution to
[TABLE]
for arbitrary , coming back to Eq. (76) only at the very end. Following the method of variation of constants, Korn and Korn (1961) we can seek the solution as
[TABLE]
where and are the solutions to the homogeneous equation with , see Eqs. (71) and (72). Since we are looking for the Green’s function, we demand that there are no waves incoming from . Hence, we pose the boundary conditions
[TABLE]
Inserting the trial solution (78) into (77), we find
[TABLE]
We now confine ourselves to the situation where the source term vanishes at . Then, because of the linear independence of the solutions, for positive , and we only need to consider the region , i.e. we only have a source in the electron region. Multiplying Eq. (80) by and using our previous results for the probability current, we find that
[TABLE]
where the last equality defines for . Multiplying Eq. (80) by , we find that
[TABLE]
Using the result (81), we obtain
[TABLE]
where we have defined for in the last equality.
Having obtained the derivatives and , we can find the coefficients themselves by integrating. Taking the boundary conditions (79) into account, we have
[TABLE]
Inserting this into Eq. (78), using the results (81) and (83), and using that for , we obtain
[TABLE]
where
[TABLE]
is the Green’s function. Note that the zero upper boundary for is not a fundamental limitation. Rather, it is a result of the fact that we have confined our attention to the situation where vanishes for . If desired, one can expand the description and determine the derivatives and for positive by multiplying Eq. (80) by and . This gives a natural way of expanding the definitions of and to the region .
Now let us come back to the original problem. Comparing Eqs. (76) and (77), we see that
[TABLE]
We are mainly interested in the Green’s function in the hole region, i.e. the region . Using the definition of and of for , we obtain
[TABLE]
Finally, applying the inverse Fourier transform to the Green’s function (90), we find the solution to Eq. (7) as
[TABLE]
where the classical action is given by
[TABLE]
Note that we have set the integration limits in Eq. (91) to , which is defined in the last paragraph of section II.2, since the action becomes imaginary for larger values of . Classically, corresponds to the maximal angle for which an electron emitted by the source can propagate to the hole region. For , we have and larger momenta do not give rise to propagating waves, whence we can ignore their contribution far from the source. For , is defined by Eq. (15). Modes with momentum larger than will (classically) be reflected by the barrier, since they cannot propagate in the hole region, and can therefore be ignored sufficiently far away from the barrier.
Appendix B Evaluation of oscillatory integrals
In this appendix, we consider the prototype integral
[TABLE]
in the limit . The vector is two-dimensional and we assume that vanishes for sufficiently large . In the first subsection, we allow to be a vector; in the second and third subsections, we only consider scalar . The scalar function is henceforth referred to as the action. When becomes small, we can apply the so-called stationary phase approximation Fedoryuk (1977); Guillemin and Sternberg (1977); Maslov and Fedoryuk (1981) to the integral (93). In this appendix, we briefly discuss this method for regular points and near caustics. In particular, we will discuss in detail how we can obtain good results near a cusp caustic. Although we will also give some derivations, the emphasis will be on the results.
B.1 WKB approximation
In this subsection, we allow the variable of integration to be an -dimensional vector. In the limit , the main contribution to the integral (93) is given by the critical points, Fedoryuk (1977); Guillemin and Sternberg (1977); Maslov and Fedoryuk (1981) where the gradient of the phase function with respect to vanishes, i.e.
[TABLE]
Let us start by considering the simplest case, where the action has a nondegenerate critical point :
[TABLE]
which means that the Hessian matrix is invertible at the critical point. In this case, the implicit function theorem states that there exists a neighborhood of and a function such that Eq. (95) holds for all points in this neighborhood. One can then show that, for a nondegenerate critical point , one has Fedoryuk (1977); Guillemin and Sternberg (1977); Maslov and Fedoryuk (1981)
[TABLE]
In the physical literature, these kind of approximations, notably for the one-dimensional case, are usually referred to as the Wentzel-Kramers-Brillouin (WKB) approximation. Although we do not discuss it here, we remark that the need to make a consistent choice for the sign of naturally leads to the notion of the Maslov index, see Refs. Maslov, 1973; Maslov and Fedoryuk, 1981. Furthermore, because of the aforementioned implicit function theorem, the result (96) can be extended to a neighborhood of the point .
When there are multiple critical points for a given value of , one has to compute the right-hand side of Eq. (96) for each of them. The integral (93) then equals the sum of these results. From a physical perspective, this means that we have interference between multiple trajectories.
B.2 Fold caustic: Airy approximation
From here on, we consider only scalar . When the critical point is degenerate, that is,
[TABLE]
the approximation (96) diverges and is no longer valid. In this appendix, and in appendix B.3.1, we show how the leading order term of the asymptotic expansion of for near can be obtained near a fold and a cusp caustic. We remark that there is an extensive body of literature on approximating the integral near a caustic, and that most approximations were gradually developed. A good overview of the various approximations in the context of semiclassical collission theory, a subject in which they have been used extensively, is given in Ref. Connor and Farrelly, 1981a. The derivations that we present in this appendix and in appendix B.3.1 closely follow appendix 2 of Ref. S. Yu. Dobrokhotov et al., 2014, only extending some of their arguments. In turn, the derivation presented there makes extensive use of the ideas of catastrophe theory and the stationary phase approximation, as presented in Refs. Fedoryuk, 1977; Arnold et al., 1982. We have nevertheless chosen to include these derivations in this appendix, in order to make the paper self-contained and to make appendix B.3.2 more accessible to the reader.
When the degenerate stationary point lies on a fold caustic, we know Arnold et al. (1982) that the third derivative of the action does not vanish. Let us therefore consider the Taylor expansion of the action up to third order in around , i.e.
[TABLE]
where and . We note that and are related through Eq. (97). Subsequently, we expand the coefficients up to first order in , that is,
[TABLE]
Note that the constant parts of and vanish due to Eq. (97).
We now show how to express the leading order term of the asymptotic expansion of the integral (93) near the fold caustic in terms of the Airy-function Airy (1838) , which has the integral representation
[TABLE]
In the integral , we first make the substitution , and subsequently , where . We also make a Taylor expansion of in around . We then see that the leading order term in the asymptotic expansion is , whereas the terms of in the action give a contribution of . Similarly, the first order term in the expansion of gives a contribution of . Using Eq. (98) and making the above substitutions, we therefore arrive at
[TABLE]
We now need to determine in which neighborhood of the fold caustic we can use this formula. To this end, we note that
[TABLE]
The integral is heavily oscillating for small and has different asymptotic expansions for different values of . In particular, we cannot use Eq. (102) when we are far away from the point on the fold caustic, since the coefficient will be too large to justify the equality (101). More specifically, we could say that the argument of the Airy function should not be large. If it were large, we could expand the Airy function for large arguments and we would be in the regime of the WKB approximation. Therefore, a safe estimate seems to be to demand that the argument of the Airy function is , with . Setting for instance , we find from Eq. (104) that we can use Eq. (102) in an neighborhood of the fold caustic. Of course, this is an estimate and it may be possible to use the approximation in a larger neighborhood. This is however dependent on the details of the problem, for instance on the values of the coefficients and .
Since for , we have and , we see that the errors that are introduced by neglecting the second order terms in the Taylor expansions of the coefficients are smaller than those introduced in Eq. (101). Using Eqs. (103) and (104) and keeping only the zeroth order term of the Taylor expansion of around , we can then simplify Eq. (102) to
[TABLE]
B.3 Cusp caustic: Pearcey approximation
B.3.1 Leading order approximation
For a point on the cusp caustic, the third derivative of the action vanishes as well, but the fourth derivative does not. We therefore expand the action up to fourth order in , similar to Eq. (98):
[TABLE]
where . As in the previous section, we expand the coefficients up to first order in . The expansions of , and are equal to those in Eq. (99). For the other coefficients, we have
[TABLE]
For near , we can then express the leading order term of the asymptotic expansion of the integral (93) in terms of the Pearcey function Pearcey (1946) , which is defined by the integral
[TABLE]
where the superscript plus or minus corresponds to the sign in front of the term. This function has two important symmetries, namely Connor and Farrelly (1981b); Connor and Curtis (1982)
[TABLE]
These can be easily verified using the definition (108). Furthermore, its two partial derivatives satisfy
[TABLE]
Although the Pearcey function is not implemented in most computer algebra systems, it can be efficiently computed using the methods in Refs. Connor and Farrelly, 1981b; Connor and Curtis, 1982.
As in the previous section, we first make the substitution , and subsequently , where , in the integral . Making a Taylor expansion of in around , we see that the leading order term in the asymptotic expansion is . The terms of in the action give a contribution of , as does the first order term in the Taylor expansion of . Therefore,
[TABLE]
After performing the aforementioned substitution, we obtain the following expression for :
[TABLE]
The sign in is taken as the sign of , and hence as the sign of .
Making use of the expansions of the , given in Eqs. (107) and (99), we find that
[TABLE]
Following the reasoning in the previous section, we then demand that both arguments of the Pearcey function are , with . Setting for example , we see from Eqs. (114) and (115) that a safe estimate for the neighborhood in which we can use Eq. (112) is .
Assuming that , we find that , and . Therefore, for this neighborhood, the largest error that we introduce in making a first order Taylor expansion of the coefficients is , which is much smaller than the errors of that are introduced in Eq. (111).
Using the above results, and replacing by its zeroth order Taylor approximation around , we can simplify Eq. (112) to
[TABLE]
We can use this approximation in an -neighborhood of the point .
B.3.2 Higher order corrections
Let us now look at higher order corrections to Eq. (116). In the previous section, we identified two sources of corrections of , namely the terms of in the action, and the first order term in the Taylor expansion of in . We also saw that corrections that come from the Taylor approximations in are of . Therefore, let us look at higher order terms in the Taylor expansion of the action, i.e. S. Yu. Dobrokhotov
[TABLE]
where was defined in Eq. (106). Let us also consider higher order terms in the Taylor expansion of the amplitude in , that is, S. Yu. Dobrokhotov
[TABLE]
When we make the substitutions and , where , in Eq. (93), we see that is and that is . Therefore, these terms are small compared to , which is of order one, and we can make a Taylor expansion of the exponent. Gathering all terms of the same order, we obtain
[TABLE]
We know from the discussion in the previous section that the first term becomes Eq. (116). Let us now look at the first term of . After we have performed the substitutions and have discarded all terms of and higher, we find that it equals
[TABLE]
One can prove the last equality by using the same arguments as in the previous subsection. By we mean the derivative of the Pearcey function with respect to its second argument, given by
[TABLE]
In a similar way, one can express the first term of in Eq. (119) in terms of the derivative of the Pearcey function with respect to its first argument, given by
[TABLE]
After some algebra, we obtain
[TABLE]
We will not go into the other terms in Eq. (119) at this point, as they will prove to be irrelevant for the problem that we discuss in the main text.
B.4 Uniform approximation near the cusp
In the previous sections, we performed a Taylor expansion of the action until the first nonvanishing term, and constructed an approximation for the integral (93) based on this expansion. Using the theorems of catastrophe theory, a uniform approximation of can be constructed. For points near the fold caustic this was first done in Ref. Chester et al., 1957, and for the cusp caustic in Refs. Ursell, 1972; Connor and Farrelly, 1981a. In this appendix, we summarize the construction of the uniform approximation near a cusp caustic as presented in Refs. Ursell, 1972; Connor and Farrelly, 1981a, using slightly different conventions.
Consider a cusp point , at which the first three derivatives of the action with respect to vanish. Then, for points in the vicinity of this cusp point , a transformation exists, with inverse transformation , such that the action in the new variable has the form Poston and Stewart (1978); Arnold et al. (1982); Ursell (1972)
[TABLE]
where the sign in front of equals the sign of and both and . Note that this is an exact transformation and that we are no longer using a truncated Taylor expansion here.
Changing our integration variable in the integral (93) from to , we obtain
[TABLE]
where we have introcuded the new amplitude function
[TABLE]
Subsequently, we expand up to second order in , i.e.
[TABLE]
One can then show that the following equality holds Ursell (1972)
[TABLE]
where the derivatives of the Pearcey function were defined in Eqs. (121) and (122) and the sign in the definition of the Pearcey function corresponds to the sign in front of in Eq. (124). When higher order terms in the expansion (127) are taken into account, the constants , and in Eq. (128) are replaced by series in integer powers of , as proven in Ref. Ursell, 1972.
In order to compute the solution (128), we need to obtain the parameters , and , and , and for a given point . In order for the mapping to be one-to-one, the stationary points on the left-hand side of Eq. (124) should correspond to those on the right-hand side. Concerning the left-hand side of the equation, let us assume that we know the action at the three stationary points of for a given . On the right-hand side, the stationary points are defined by the equation
[TABLE]
When the discriminant
[TABLE]
is positive, this cubic equation has three distinct real roots. When is negative, the equation has one real root and two complex conjugate roots and when the discriminant vanishes, all roots are real, but there is a multiple root. We call these three roots of Eq. (129) . Note that in a practical implementation, it is important that both sets of stationary points are ordered in the same way, e.g. from small to large when all numbers are real. When this is not the case, one could for instance take the first stationary point to be the real one, followed by the two complex ones ordered by their imaginary part.
Requiring that the stationary points on both sides of Eq. (124) coincide, we find that the following set of equalities has to hold.
[TABLE]
Let us, for simplicity, disregard the case of a degenerate critical point for a moment. Then we can subtract the second equation in (131) from the first and the third from the first to eliminate , which gives
[TABLE]
Given initial guesses for the parameters and , one can find the three stationary points from Eq. (129). These can be inserted into Eq. (132) to find new values for and and this process can be iterated until self-consistency is reached. As initial guesses for the parameters and one can use the result of the Taylor expansion of the action,
[TABLE]
cf. Eq. (116). For a degenerate critical point, we can use the fact that the discriminant vanishes to obtain a relation between and . Inserting this into Eq. (129), one obtains expressions Connor and Farrelly (1981a) for in terms of . We can then obtain from Eq. (132). For more details, we refer to Ref. Connor and Farrelly, 1981a. Alternatively, one can obtain , and by using an algebraic method, that was described in Ref. Connor and Farrelly, 1981a. However, this method requires tracing certain solutions across the caustic, and we therefore find it less convenient.
Finally, one needs to determine the parameters , and in the expansion (127). By combining Eqs. (126) and (127) and evaluating the result at the three critical points , we arrive at a system of three linear equations
[TABLE]
from which , and can easily be found. However, this requires the computation of the derivative at the stationary points. By taking the second derivative with respect to on both sides of Eq. (124), we find that
[TABLE]
We are interested in the value of this derivative at the critical points, where the second term on the left-hand side vanishes by definition. When the stationary point is nondegenerate, the first term on the left-hand side is nonzero and we obtain
[TABLE]
When we are dealing with a fold point, where vanishes as well, we can take the derivative of expression (135) with respect to once more to obtain an equation for . For the cusp point, we can obtain by differentiating Eq. (135) twice. Combining Eqs. (134) and (136), we see that we can obtain , and when we know the values of and at the critical points.
Appendix C Initial value problem
In this appendix, we construct the solution of the initial value problem (41), considered in section VII. First, we take the Fourier transform of both the equation and the initial condition with respect to , which gives
[TABLE]
In appendix A.1, we constructed two linearly independent solutions of the eigenvalue problem, and , see Eqs. (71) and (72). The solution of Eq. (137) is a linear combination of these two solutions, i.e.
[TABLE]
that satisfies the initial condition. Here, \big{(}\overline{\Psi}_{>}(x)\;\;\;\overline{\Psi}_{<}(x)\big{)} denotes the matrix with columns and . Inserting the initial condition into Eq. (138) and multiplying by the inverse of the matrix on the right-hand side, we find that
[TABLE]
Combining the previous results and taking the inverse Fourier transform with respect to , we find the solution of the initial value problem (41) as
[TABLE]
where the integration limits are , as states with higher transversal momentum do not propagate in the sample, see also the discussion at the end of appendix A.
When the electronic current only flows into the sample from the left-hand side, we can confine our attention to the term proportional to (see the main text). This means that we are only interested in the first component of the vector \big{(}\overline{\Psi}_{>}(x_{s})\;\;\;\overline{\Psi}_{<}(x_{s})\big{)}^{-1}\;\overline{\Psi}_{0}(p_{y}). After some calculus, we find that this first component equals
[TABLE]
Therefore, we can approximate the full solution (140) by
[TABLE]
When one considers a situation with , and is interested only in the hole region, the integration limits should be reduced to , as discussed at the end of appendix A.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Veselago (1968) V. G. Veselago, Sov. Phys. Usp. 10 , 509 (1968).
- 2Smith et al. (2000) D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, Phys. Rev. Lett. 84 , 4184 (2000).
- 3Houck et al. (2003) A. A. Houck, J. B. Brock, and I. L. Chuang, Phys. Rev. Lett. 90 , 137401 (2003).
- 4Grbic and Eleftheriades (2004) A. Grbic and G. V. Eleftheriades, Phys. Rev. Lett. 92 , 117403 (2004).
- 5Tretyakov et al. (2003) S. Tretyakov, I. Nefedov, A. Sihvola, S. Maslovski, and C. Simovski, J. Electromagn. Waves. Appl. 17 , 695 (2003).
- 6Pendry (2004) J. B. Pendry, Science 306 , 1353 (2004).
- 7Zhang et al. (2009) S. Zhang, Y.-S. Park, J. Li, X. Lu, W. Zhang, and X. Zhang, Phys. Rev. Lett. 102 , 023901 (2009).
- 8Xiong et al. (2010) X. Xiong, W.-H. Sun, Y.-J. Bao, M. Wang, R.-W. Peng, C. Sun, X. Lu, J. Shao, Z.-F. Li, and N.-B. Ming, Phys. Rev. B 81 , 075119 (2010).