{'fraction_non_alphanumeric': 0.035629008848728876, 'fraction_numerical': 0.010346926354230066, 'mean_word_length': 4.318600597609562, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quasiparticle scattering from vortices in d-wave superconductors I: Superflow contribution

In the vortex state of a d-wave superconductor, massless Dirac quasiparticles are scattered from magnetic vortices via a combination of two basic mechanisms: effective potential scattering due to the superflow swirling about the vortices and Aharonov-Bohm scattering due to the Berry phase acquired by a quasiparticle upon circling a vortex. In this paper, we study the superflow contribution by calculating the differential cross section for a quasiparticle scattering from the effective non-central potential of a single vortex. We solve the massless Dirac equation in polar coordinates and obtain the cross section via a partial wave analysis. We also present a more transparent Born-limit calculation and in this approximation we provide an analytic expression for the differential cross section. The Berry phase contribution to the quasiparticle scattering is considered in a separate paper.


I. INTRODUCTION
The low energy excitations of a d-wave superconductor are the Bogoliubov quasiparticles that reside in the vicinity of the four gap nodes at which the order parameter vanishes. The physics of these quasiparticles can be probed experimentally via a variety of low temperature electrical and thermal transport measurements. Since quasiparticles are part electron and part hole, their energy is well defined but their charge is not. Thus, it is thermal current that follows quasiparticle current. Since T T c , we are well within the superconducting state. Since T ∆ 0 (the gap maximum), transport is dominated by quasiparticles excited in the vicinity of the gap nodes. (Quasiparticle dispersion is therefore given by the anisotropic Dirac spectrum, E = (v 2 f k 2 1 + v 2 2 k 2 2 ) 1/2 ,where v f is the Fermi velocity, v 2 is the slope of the gap, and k 1 and k 2 are defined locally about each node. We shall choose our axes such that gap nodes are located at ±p Fx and ±p Fŷ in momentum space.) Furthermore, the temperatures of interest are low enough that sources of inelastic scattering are frozen out. This is what we mean by low temperature quasiparticle transport. In the mixed state, the remaining energy scales are the impurity scattering rate, Γ 0 , the vortex scattering rate, Γ v , and the temperature. We will focus on the weak-field high-temperature regime where Γ v T and T is the dominant energy scale. In this regime, the quasiparticles responsible for transport are thermally generated rather than impurity-induced 1-3 or magnetic field-induced [4][5][6] .(Note that thermal transport in the opposite, low T , regime has been discussed frequently in the literature [7][8][9][10][11] .) In this high T regime, the physical situation is relatively simple. Thermally excited quasiparticles carry the heat current. To understand the thermal transport, we need only understand how they scatter from magnetic vortices. Furthermore, since H H c2 , the vortices are dilute, separated by distances large compared to the quasiparticle de Broglie wavelength. We can therefore learn a lot by considering the scattering of quasiparticles from a vortex and assuming that all such scattering events are independent. We take the following path. In Sec. II, we develop our picture of a nodal quasiparticle scattering from a vortex. Then in Sec. III, we consider the Bogoliubov-de Gennes equation, make a singular gauge transformation, and shift the origin of momentum space to the location of one of the gap nodes. The resulting problem is one of an (anisotropic) Dirac fermion scattering from an effective non-central potential (due to the superflow) in the presence of antiperiodic boundary conditions (a consequence of our gauge choice). After obtaining the quasiparticle current functional in Sec. IV, we define our vortex model in Sec. V and make several important approximations. As far as scattering is concerned, the vortex has two parts: a circulating superflow and a Berry phase factor of (-1) acquired upon circling a vortex. In this paper, we consider scattering due to superflow current in a single vortex, without the Berry phase effect. (We consider the Berry phase contribution in a separate paper 12 .) Furthermore, we neglect the anisotropy of the Dirac dispersion and take v f = v 2 . Although these approximations will prevent our results from being quantitatively accurate, we argue that our qualitative results reflect the essential physics of the problem. This notion is supported by Ref. 13, wherein this model, with these approximations, was used to derive a thermal Hall conductivity in good agreement with the qualitative features of the experimental results for Y Ba 2 Cu 3 O 6.99 by Zhang et al 22 . In Secs. VI and VII, we expand the quasiparticle wave function in angular momentum eigenstates and obtain the resulting radial equations. We solve these equations outside the vortex in Sec. VIII and inside the vortex in Sec. IX. In Sec. X, we match solutions and sum over all four gap nodes in order to obtain the single vortex scattering cross section. Results of our numerical arXiv:1006.2818v3 [cond-mat.supr-con] 13 Aug 2011 Figure 1: Schematic depiction of vortices in the mixed state of an extreme type II superconductor. Vortices are pinned to randomly distributed pinning centers and are separated by an average inter-vortex distance, 2R. Since ξ R λ, vortex cores are point-like while the magnetic field profile of each vortex is quite extended. Therefore, the magnetic field is approximately uniform across the sample even though the vortices are well separated.
calculations are presented in Sec. XI. In the same section (Sec. XI) we write a closed-form expression for the differential cross section using the more transparent Born approximation details of which are given in Appendix A. Conclusions are discussed in Sec XII. Note that the results of this calculation were used in Ref. 13 (without the detailed derivation presented here) to calculate thermal transport coefficients.

II. SINGLE VORTEX SCATTERING OF BOGOLIUBOV QUASIPARTICLES
Consider the mixed state of a d-wave superconductor. Our picture is that of Fig. 1. In the presence of a magnetic field (H > H c1 ), magnetic vortices penetrate the sample (a 2D CuO 2 layer). Vortices are distributed randomly, pinned to local defects. The cuprates are extreme type II superconductors in which the coherence length, ξ, is much smaller than the penetration depth, λ. As a result, while the vortex cores may be well separated, the magnetic field profiles overlap significantly such that there is little variation in the magnetic field across the sample. We therefore adopt the extreme type II limit of ξ → 0 and λ → ∞ and take the magnetic field to be constant, H = Hẑ. In this limit, there are only two remaining length scales. The first, 1/k, is set by the temperature such that Note that in the isotropic limit, this reduces to the magnitude of the quasiparticle momentum (measured from a node). The second length scale, R, is half of the average distance between vortices. With one flux quantum per vortex, HπR 2 = Φ 0 = hc/2e, so we define R ≡ c/eH. In terms of R, we can define the (2D) density of vortices to be n v = H/Φ 0 = 1/πR 2 . Note that the ratio of these two length scales yields which is the inverse of the argument of the Simon-Lee scaling functions 14 .
We wish to consider a quasiparticle, in a state with a particular energy and current, scattering from a magnetic vortex into another state, of the same energy, with a different current. At low temperatures, quasiparticles are excited in the vicinity of the four gap nodes in momentum space. The nodal structure of the Brillouin zone is depicted in Fig.2(b). Following Simon and Lee 14 , we have defined coordinate axes, rotated by 45 • from the usual crystal axes, such that the gap nodes are located on-axis at ±p Fx and ±p Fŷ . In the neighborhood of a node, the quasiparticle excitation energy is given by the anisotropic Dirac spectrum, E k = v 2 f k 2 1 + v 2 2 k 2 2 , where the k 1 and k 2 axes are defined, respectively, to be perpendicular and parallel to the local Fermi surface. The surfaces of equal energy are therefore ellipses centered about each of the nodes. The quasiparticle current is given by the group velocity v G = ∇ k E k , directed outward from the node centers, perpendicular to the ellipses of constant energy.
In coordinate space, the quasiparticle is represented as an incident plane wave. The single (hc/2e) vortex is centered at the origin and surrounded by a circulating supercurrent. Due to the interaction between the quasiparticle and the vortex, the quasiparticle is scattered into an outgoing radial wave. The situation is depicted schematically in Fig. 2(a). By solving the Bogoliubov-de Gennes equation for the quasiparticle wave functions and considering the nature of the quasiparticle current before and after such a scattering event, we will compute the quasiparticle cross section for scattering from a vortex in a d-wave superconductor. Our analysis shall be similar in spirit to that conducted for the s-wave case by Cleary 15,16 , who extended the work of Caroli, de Gennes, and Matricon 17 (on vortex core bound states) to the problem of quasiparticle scattering from a vortex in an s-wave superconductor.

III. BOGOLIUBOV-DE GENNES EQUATION
Consider the Bogoliubov-de Gennes (BdG) equation for a d-wave superconductor in the presence of a constant perpendicular magnetic field, A = 1 2 Hrφ, and with an order parameter that winds once about the origin, ∆(r) = ∆ 0 e iφ : Here p = −i ∇, {a, b} = (ab + ba)/2, E is the quasiparticle energy, and the form of the gap operator is that required to provide d-wave symmetry 10,14 . Entangled within the complex differential form of the gap operator,∆, is the fact that upon circling an hc/2e vortex, the quasiparticle acquires a Berry phase factor of (-1). The Hamiltonian would be simplified if we could effectively strip the gap function, ∆(r), of its phase. This is accomplished via the application of the singular gauge transformation In this gauge, known as the Anderson gauge, where is the gauge invariant superfluid momentum (superflow), v f = p F /m, v 2 = ∆ 0 /p F , and R ≡ c/eH. Although this form of the Bogoliubov-de Gennes equation is much simplified, the boundary conditions have become more complicated. While the original wave function was defined with periodic boundary conditions, Ψ(r, φ) = Ψ(r, φ + 2π), the transformed wave function is not single-valued and has antiperiodic boundary conditions, Φ(r, φ) = −Φ(r, φ + 2π). Hence, we have introduced a branch cut such that with each trip around the origin, the wave function changes sign. In effect, the Berry phase contribution has been extracted from the Hamiltonian and encoded in the antiperiodic boundary conditions imposed on the wave function. Note that when considering a system of many vortices, the resulting sea of branch cuts can be problematic. In response, Franz and Tesanovic introduced an alternate gauge transformation better suited to the vortex lattice 18,19 . However, for the single vortex case that we consider, the Anderson gauge described above is most convenient. All observables, such as differential cross section or transport coefficients, are independent of gauge choice. We would like to emphasize that the Berry phase effect is not an artifact of our gauge choice. It is physical, and an intrinsic part of the original Bogoliubov de Gennes equations. In the untransformed Hamiltonian, which has periodic boundary conditions, it resides in the complex differential form of the gap operator. The gauge transformation drastically simplifies the Hamiltonian at the cost of introducing antiperiodic boundary conditions. Thus, effectively, the Berry phase contribution is extracted from the Hamiltonian and moved to the boundary conditions.
Recall that the low energy excitations of a d-wave superconductor are concentrated about the four nodal points in momentum space where the gap vanishes. We can therefore further simplify our Hamiltonian by shifting the origin of momentum space to the location of one of the nodes. Shifting to node 1 we find that where α = v 2 /v f and we have used the fact that P s = P s (r)φ to commute p with P s . Here H is written as the sum of a linear (Dirac) Hamiltonian, H D , and a quadratic (curvature) Hamiltonian, H C , to emphasize that the second term is smaller than the first by a factor of E/E F . We will focus on the dominant term, H D , in the remainder of this paper. In what follows we shall seek solutions to the Bogoliubov-de Gennes equation developed above. However, we must first consider the nature of quasiparticle current in a d-wave superconductor.

IV. QUASIPARTICLE CURRENT
In order to calculate the cross section for a quasiparticle scattering from a vortex, we must be able to write down the currents corresponding to the incident and scattered wave functions. Since the incident and scattered currents will be considered in the far field where the quasiparticles are free, we wish to determine the quasiparticle current as a functional of Φ for P s = 0. Setting P s = 0 in Eq. (7) we find that the BdG Hamiltonian becomes Following Ref. 20, we can write down a mean-field Hamiltonian in terms of the real-space electron creation/annihilations operators, Ψ † α and Ψ α .
Since quasiparticles carry well-defined spin (and heat) the quasiparticle current is equal to the spin current divided by the quasiparticle spin. The spin density operator is expressed aŝ where s ≡ 1/2 and the spin current operator is, in turn, determined via Making use of the fermionic anticommutation relations, we find that where η α ≡ ±1. Noting that each of these terms can be manipulated into the form of a divergence and making use of Eq. (16) we can write the spin current operator aŝ where ∇ = ∂ xx + ∂ yŷ and ∇ ⊥ ≡ ∂ yx + ∂ xŷ . The real-space creation/annihilation operators, Ψ † α and Ψ α , can be expressed as a weighted sum of Bogoliubov operators, γ † nα and γ nα , via where u n and v n are particle and hole wave functions for state n. Plugging these forms into the spin current operator and evaluating for state and spin σ yields an expression for the spin current as a functional of u and v .
Dropping the eigenstate label, dividing by the spin sη σ , and writing the result in terms of the particle-hole 2-vector we obtain the quasiparticle current As in the preceding section, it is convenient to shift the origin of momentum space to a nodal point. Shifting to node 1 yields where α = v 2 /v f . Once again, the second (curvature) term is smaller than the first (Dirac) term by a factor of E/E F and shall be neglected in what follows.

V. MODEL AND APPROXIMATIONS FOR SINGLE VORTEX
Given the BdG Hamiltonian and the quasiparticle current functional, we are ready to consider a model of quasiparticle scattering from a single vortex. In reality, each quasiparticle encounters a sea of vortices separated by distances on the order of 2R = 4 c/eH. The superflow circulating around each vortex is peaked near the vortex centers but falls off slowly and overlaps in the regions between vortices. The total flux through the sample is equal to one (hc/2e) flux quantum per vortex but is distributed smoothly since the penetration depth is much larger than the inter-vortex distance. In order to model this situation via single vortex scattering, we approximate the effect of neighboring vortices by cutting off the superflow distribution about our single vortex at a distance R from its center. By construction, the flux through this circle is exactly one flux quantum. Thus, the superfluid momentum takes the form of Eq. (8) inside the vortex (r < R) and is zero outside the vortex (r > R) Note that at the vortex edge (r = R), the phase gradient and vector potential terms cancel and the superflow goes continuously to zero. We shall consider a quasiparticle plane wave incident upon such a vortex and calculate the scattering cross section. Before going forward, however, we make the following additional approximations.
First, recall from Sec. III that, with our gauge choice, the wave function must obey antiperiodic boundary conditions such that it changes sign each time the quasiparticle winds around the vortex. This is the Berry phase contribution to the problem and it has the interference effect of an Aharonov-Bohm π-flux at the origin. Thus, even in the absence of a superflow, it is possible for quasiparticles to scatter, due solely to the antiperiodic boundary conditions. The effect of this Berry phase is considered in a companion paper.
Second, recall that the quasiparticle excitation spectrum for a d-wave superconductor is that of an anisotropic Dirac cone with two characteristic velocities, v f and v 2 , which can be quite different. This anisotropy enters the BdG equation and current functional in the form of the parameter α = v 2 /v f . The fact that α is not equal to one has the effect of complicating the form of the BdG equation in polar coordinates. Thus, in order to make this partial differential equation separable, we should scale out the anisotropy. This can be accomplished by scaling the y-coordinate by a factor of α. In the scaled coordinates, the elliptical gap nodes become circular and we obtain the isotropic Dirac equation. Unfortunately, this scaling also has a negative consequence. Whereas the vortex was originally circular in real space, the scaling makes it elliptical. Thus, in the scaled coordinates, the superfluid momentum (and the vortex boundary) become much more complicated. While such complications can be dealt with (via a significantly more involved computation), they were not considered in this investigation. Possible consequences of anisotropy have been studied in Ref. ? , where it was noted that there are no drastic experimental signatures due to the anisotropy of the Dirac cone. Thus, for simplicity, rather than scaling out the anisotropy, we consider the more straightforward case where the gap nodes are isotropic and v f = v 2 . For this isotropic case, both the gap nodes and the vortex are circular and we can separate the BdG equation in polar coordinates. However, since v f exceeds v 2 by a factor of 10 to 20 in the cuprates 21 , this is clearly an approximation. As a result, we expect only qualitative agreement with reality. That this is the case was demonstrated by Durst, Vishwanath, and Lee 13 who used this approximation to provide a clear, though qualitative, explanation of the thermal Hall conductivity measurements of Ong and co-workers 22 Finally, recall that the BdG Hamiltonian can be expressed as the sum of a linearized (Dirac) part, H D , and a quadratic (curvature) part, H C . We see that H C is small compared to H D as long as we are sufficiently far from the vortex center (r > 1/p F ). Hence we must cutoff our model at the scale of the vortex core (ξ ∼ 10/p F ). As our final assumption, we select a reasonable core size and model the vortex core as a region with vanishing superflow. We now have a well-defined scattering problem, which is solved considering the linearized Hamiltonian. Since we shall only consider the isotropic case from this point forward, the Hamiltonian becomes the quasiparticle energy is and the quasiparticle current functional takes the form Note that the effective potential induced by the superflow is smooth on the scale of 1/p F . The Fourier components for large-momentum inter-node scattering are therefore small. Hence, we expect that the contribution to the scattering cross section of inter-node scattering will be subdominant to that of intra-node scattering. Thus, for simplicity, we include only intra-node scattering in what follows.

VI. ANGULAR MOMENTUM EIGENSTATES
Consider the linearized Hamiltonian, H D , in the isotropic limit, Eq. (28). Note that for P s = 0, this is just the Dirac Hamiltonian for massless spin-1/2 fermions in two dimensions with the γ-matrix representation For the physical Dirac equation we know that the total angular momentum operator takes the form, J = L + 1 2 Σ where Σ = i 2 3jk γ j γ k . Therefore, in the above γ-matrix representation Evaluating the commutator with H 0 D we find (as expected) Therefore, there exists a complete set of simultaneous eigenstates of H 0 D and J. The eigenstates of L = −i ∂ ∂φ take the form and the eigenstates of Σ = τ 2 take the form Therefore, the simultaneous eigenstates of J and H 0 D have the form where the radial functions, f n (r) and g n (r), are determined from the solution of a pair of coupled radial equations.
Since we have neglected the Berry phase contribution and adopted periodic boundary conditions for this case, the requirement of single-valued wave functions demands that n is an integer.
For P s = 0, H D looks like the massless Dirac equation in the presence of an effective scalar potential, where (τ ·r) ≡ τ 3 cos φ + τ 1 sin φ and we note that (τ ·r) 2 = 1. It is also useful to write the wave function as where we have defined Note that (τ ·r) transforms χ + n into χ − n and vice versa, and also that cosines and sines shift the χ's up and down in angular momentum, Making use of these relations, we see that both HΦ and EΦ can be written as sums over χ ± n weighted by coefficients that are functions only of the radial coordinate. Equating coefficients of χ ± n on both sides of the BdG equation and dividing out an overall factor of −iv f k yields a set of fully coupled differential equations for f n and g n , where we have defined a dimensionless radial coordinate, ρ ≡ kr, and a dimensionless superfluid momentum, P s (ρ) ≡ (1/ρ − ρ/(kR) 2 )/2. In this manner, the solution of the BdG equation is reduced to the that of a system of coupled ordinary differential equations.

VIII. OUTSIDE THE VORTEX
According to our model of the single vortex, the superfluid momentum vanishes at a distance, R, from the origin. Thus, for ρ > kR, we sit outside the vortex. Here P s = 0 and the quasiparticles are free. Since the superfluid momentum (Eq. 8) includes the effects of both the order parameter phase gradient and the magnetic field, quasiparticles outside the vortex are subject to neither. Thus, for ρ > kR, we consider quasiparticles subject only to the free Dirac Hamiltonian, H 0 D , and with current defined via the linearized current functional, j D [Φ]. In this regime, it is possible to obtain explicit solutions to the BdG equation, define incident and scattered wave functions, and, by constructing these incident and scattered waves from the free basis functions, write down an expression for the scattering cross section in terms of the coefficients of the basis functions. Then we need only match solutions with those inside the vortex to obtain the cross section.

A. Free Solutions
For P s = 0 , the radial equations, (44) and (45), take the dramatically simpler form Note in particular that, while f n is coupled to g n , the equations for functions of different n are independent. Eliminating g n from the equations above yields which is the defining equation for the Bessel functions. Therefore, where A n and B n are complex constants and J n (ρ) and Y n (ρ) are Bessel functions of the first and second kind. Substituting back for g n then yields where we have used the Bessel function identity (∂/∂ρ − ν/ρ)Z ν = −Z ν+1 . Hence, we can write down the free wave function where the coefficients shall remain undetermined until we match with solutions inside the vortex.

B. Incident and Scattered Waves
We wish to obtain the cross section for a plane wave, with quasiparticle current in the incident direction, scattering off a vortex as a radial wave, with quasiparticle current in the scattered direction. If the incident momentum is k = (k, θ) and the final momentum is k = (k, φ), then the incident direction is the direction of the group velocity at momentum k and the scattered direction is the direction of the group velocity at momentum k . For general, anisotropic nodes, the group velocity need not be parallel to the momentum. However, for the isotropic case that we consider and the group velocity and momentum are parallel. Therefore, if Φ i denotes the incident wave function and Φ s denotes the scattered wave function, then we require Recall that outside the vortex, quasiparticles are subject to neither an order parameter phase gradient nor a magnetic field. Thus, the incident wave function is a plane wave. This is consistent with the well-known results of Franz and Tesanovic (Franz and Tesanovic, PRL 84, 554 (2000)) who showed that the low-energy quasiparticle states of a d-wave superconductor in the vortex state are Bloch waves of massless Dirac fermions rather than Landau Levels.
(For a discussion of the analyses that led to this important result, the reader is referred to Refs. 25, 26, 27, and 19.). Inspection of the form of the current functional, Note also that this form solves the BdG equation, as it must in the absence of the vortex. The appropriate scattered radial wave is then Here f (φ − θ) is the scattering amplitude and the e iφ/2 prefactor has been added to make the wave function singlevalued.

C. Constructing the Cross Section
We shall now construct, from our free solution basis functions, an asymptotic wave function containing the correct incident and scattered waves. The factor in front of the scattered wave will then be our scattering amplitude.
We begin with the free wave function obtained in Eq. (50). Note that in the asymptotic limit, the integer-index Bessel functions can be expressed as J n (ρ) = η n 2 πρ cos(ρ − π/4 − |n|π/2) Y n (ρ) = η n 2 πρ sin(ρ − π/4 − |n|π/2) where η n = 1 for n ≥ 0 and η n = (−1) n for n < 0. Therefore, if we shift n → n − 1 in the second term of Eq. (50), plug in the asymptotic forms, and decompose the sines and cosines into exponentials, we obtain We can replace our two complex constants, A n and B n , with two new complex constants, a n and b n , by defining A n + iB n ≡ i n e −i(n+ 1 2 )θ 1/2 + e i π 4 (−1) n a n .
Making use of these definitions, noting that η n i n = i |n| , and reorganizing terms, we find Noting the Bessel function expansion of a plane wave Since we have yet to restrict a n and b n , this asymptotic wave function is still totally general. However, we have succeeded in rearranging it into a suggestive form. The three terms above are easily understood. The first is our incident plane wave, the second is the outgoing radial wave, and the third is an incoming radial wave. By construction, we require an incident plane wave and an outgoing radial wave. To realize this scenario, our asymptotic boundary conditions require that there be no additional incoming wave. Thus, we require g(ϕ) = 0 and must therefore set a n = 0 for all n. In terms of our original constants, this restriction requires that Furthermore, we obtain a simple form for the asymptotic wave function where we recall that k = (k, θ), r = (r, φ), and ϕ = φ − θ.
Applying the current functional, j D , to both the incident and scattered parts of this wave function yields the incident and scattered current density Then the differential cross section is Integrating over ϕ yields the total cross section By weighting with a factor of (1 − cos ϕ), we obtain the transport cross section and by weighting with a factor of sin ϕ we find the the skew cross section All that remains is to calculate the coefficients b n . These, of course, are determined by the details of the quasiparticle scattering. Hence, we must now look inside the vortex.

IX. INSIDE THE VORTEX
Inside the vortex, the situation is more complicated. For ρ < kR, the superfluid momentum is nonzero and takes the form, P s = (1/ρ−ρ/(kR) 2 )/2. Therefore, even the linearized Hamiltonian, H D , mixes angular momentum eigenstates. As discussed by Simon and Lee 14 , Ye 11 , and Vishwanath 9 , since the linearized Hamiltonian is time-reversal invariant, the resulting skew cross section (after summing over all four nodes) is zero. We solve for the radial functions, f n (ρ) and g n (ρ). The radial equations look as follows: Since the radial equations are fully coupled by the superfluid momentum terms, the equations for all n should be solved simultaneously via an infinite-dimensional matrix equation. By numerical necessity, we shall cut off the coupling at some large n such that we consider a total of N angular momentum eigenstates where − N 2 ≤ n ≤ N 2 − 1. This is physically reasonable as we do not expect very large angular momenta to have a significant effect on the low energy physics. The result is a 2N × 2N matrix differential equation where z(ρ) is a 2N -component vector containing the f n 's and the g n 's, M(ρ) is a 2N × 2N matrix. Since only neighboring angular momenta are coupled, M(ρ) is a rather sparse matrix. Furthermore, due to the form of Eqs. (75) and (76) and the simple ρ-dependence of the superfluid momentum, this matrix can be written as where B, A 0 , and A 1 are constant 2N × 2N matrices. In terms of this matrix notation, our procedure will be to solve the homogenous equation.

A. Homogeneous Solutions via Method of Frobenius
Consider the homogeneous equation Since M(ρ) diverges as 1/ρ at the origin, the equation has a regular singular point at ρ = 0. Thus, its solution requires local analysis about the origin. In particular, we shall employ a matrix generalization of what is known as the Method of Frobenius 28 .
Consider a solution of the (Frobenius) form where α is a complex number and a m is a vector coefficient for each integer m. Plugging this into our differential equation yields Shifting m → m − 1 for the A 0 term, shifting m → m − 2 for the A 1 term, and noting that the total coefficient of each power of ρ must equal zero, we obtain the following equations: Note that Eq. (82) is just the eigenvalue equation for B. It is solved if α is an eigenvalue of B and a 0 is the corresponding eigenvector. This eigensystem has 2N solutions which we can label α k and a k 0 for k = 1 . . . 2N . The other two equations, (83) and (84), define matrix recursion relations. If B−(α k +m)1 1 is non-singular for m = 1, 2, 3, . . ., then we can solve these equations to obtain a k m given a k 0 and α k . This condition will always be satisfied for the largest eigenvalue of B but will fail for any α k for which α k +m is also an eigenvalue of B. In general, this circumstance would require additional solutions not of the Frobenius form. Fortunately, for P s = 0, the eigenvalues of B are separated by non-integer numbers. Therefore, we can obtain 2N solutions of the Frobenius form and use the equations above to calculate the vector coefficients recursively. Due to the constraints of our numerics, we cutoff each series at m = m max . As such series represent expansions about the origin, the more terms we include the larger the ρ at which the solutions will be valid. Since the differential equation contains no singular points besides the origin, the solution can be pushed as far from the origin as necessary by including additional terms. Since our solutions need only be valid out to the edge of the vortex, ρ = kR, we shall cutoff each series at the m max for which the solutions are valid out to kR.

B. Boundary Condition at Origin
Although all 2N solutions, z k , satisfy the homogeneous equation (79), we must reject solutions that show unphysical behavior at the origin. While it is permissible for the wave function or the current density to diverge at the origin, we must require that observable quantities like the total probability at the origin and the total current passing through the origin be well behaved. Since solutions of the Frobenius form yield contributions to the wave function of the form Φ(ρ → 0) ∼ ρ α for small ρ, we seek conditions on the allowed values of α.
Consider a circle of radius → 0 about the origin. The total quasiparticle probability within the circle is If the origin is to be just like any other point in space, then this probability must vanish as → 0. Hence, we require α > −1. While this restricts the allowed solutions, it turns out that the quasiparticle current provides a more restrictive condition.
Now consider a semicircle of radius → 0 about the origin, oriented about theθ direction. The total current passing though such a semicircle in theθ direction is For the origin to be physical, this current must vanish as → 0. Hence, we require α ≥ −1/2.
For P s = 0, the eigenvalues of B are always such that N of them are strictly larger than −1/2 and N of them are strictly smaller than −1/2. Thus, in order to satisfy these boundary conditions at the origin, we shall keep only the N solutions, z k , for which α k > −1/2. The full solution to the homogeneous equation is then where the c k are N complex coefficients to be determined via solution matching at the edge of the vortex.

X. FOUR-NODE CROSS SECTION
For ρ > kR, we found 2N radial solutions (Bessel functions) and reduced the 2N coefficients (A n ,B n ) to N coefficients with the condition that there be no incident radial wave in the asymptotic limit. (This required B n = i(A n − A 0 n ).) For ρ < kR, we found 2N radial solutions (Frobenius series) and reduced the 2N coefficients to N coefficients with the condition that the solutions be well-behaved at the origin. (This required c k = 0 for α k < − 1 2 .) At ρ = kR, all 2N radial solutions must be continuous. (Note that since the Dirac equation is 1st order, only continuity is required, while for the 2nd order Schrodinger equation, differentiability would also be required.) Thus, we have 2N equations which can be solved for the remaining 2N coefficients by matching solutions at ρ = kR. We solve the resulting 2N × 2N matrix equation via LU decomposition to obtain c k for 1 ≤ k ≤ N and A n for − N 2 ≤ n ≤ N 2 − 1. As discussed in prior sections, solutions are matched to obtain to obtain the coefficients to get scattering cross sections. Recall that at the outset of this calculation, we shifted the origin of momentum space to the center of node 1. Thus, in the discussions that followed, we have been considering quasiparticles scattered from one state in the vicinity of node 1 to another state in the vicinity of node 1. The resulting cross section is therefore only the cross section for these node-1 quasiparticles. However, given a quasiparticle current in any particular direction, quasiparticles from all four nodes will contribute equally. Thus to obtain the physical cross section, we must average over the cross sections for quasiparticles at each of the four nodes.
Our results for node 1 can be easily generalized to node j = {1, 2, 3, 4} by transforming coordinates to those appropriate to node j. In accordance with the d-wave structure of the gap, we can define a local coordinate system at each of the four nodes with ak 1 axis pointing along the direction of increasing k and ak 2 axis pointing along the direction of increasing ∆ k . Note that while nodes 1 and 3 define right-handed coordinate systems, nodes 2 and 4 define left-handed coordinate systems. We can therefore transform from node 1 to node j simply by rotating our incident and scattered angles (θ and φ) and then changing the sign of these angles to account for the handedness of the local coordinate system.
Thus, to obtain results for quasiparticles about node j, we need only input each θ j and take the output as a function of (−1) j+1 ϕ. Then the physical cross sections are Note that we do not account for quasiparticles that are scattered from one node to another. However, since the effective potential (induced by the superflow) is smooth on the scale of 1/p F , the Fourier components for large-momentum inter-node scattering will be small. Hence, it is reasonable to neglect such contributions.

XI. RESULTS
The procedure outlined in the preceding sections was implemented numerically for a range of intervortex distances, kR = γT / √ H, from 0.5 to 15. The larger the kR, the more angular momentum eigenstates that contribute to the vortex scattering cross section, and the more runtime-intensive the computation. For kR = 15, we set N = 46 and considered total angular momenta, n + 1/2, ranging from -22.5 to 22.5. Calculating solutions to the BdG equation both inside and outside of the vortex, the program matches solutions at ρ = kR to obtain the coefficients, b n , first for node 1 and then for nodes 2, 3, and 4. For each of the four nodes, Eqs. (65) and (71) are then used to calculate the differential cross section for quasiparticles in the vicinity of that node. The coefficients and differential cross sections for kR = 8 are shown in Fig. 3 (Left). Polar plots of the four-node differential cross sections for integer kR from 1 to 15 are shown in Fig. 3 (Right). In all cases, dσ/dϕ has a two peak structure and vanishes for ϕ = π. For small kR, the peaks are centered about ±π/2. As kR increases, the magnitude of dσ/dϕ increases while the peaks sweep closer to the forward direction.
The single vortex cross section discussed above was calculated to all orders in the linearized potential . However, to gain some physical insight regarding the calculated form of our results, it is instructive to consider the more transparent Born-limit calculation (valid to first order in the potential) whereby dσ/dϕ is proportional to the square of the Fourier transform of the scattering potential. For the linearized Hamiltonian, our effective potential is V D = −v f (1/r−r/R 2 )θ(R−r) sin φ where the θ-function imposes a cutoff at r = R. Given this input, the Born approximation yields an analytical expression for the differential cross section. Details of this calculation is given in Appendix A. We Figure 3: (Left) Matched coefficients and differential cross sections for kR = 8. The coefficients, bn, denote the contributions from eigenstates with total angular momentum n + 1/2. These coefficients and the corresponding differential cross sections are plotted for quasiparticles about nodes 1, 2, 3, and 4. The final row contains the four-node average dσ/dϕ plotted in both Cartesian and polar form. The schematic above indicates the directions from which quasiparticles about each node approach the non-central effective potential. (Right) Fixed-scale polar plots of the four-node differential cross section for integer kR from 1 to 15. Note that the cross section magnitude grows with increasing kR.
where J 1 (x) is the first order Bessel function of the first kind. This form comes about because the two-dimensional Fourier transform of 1/r yields 1/q where q 2 = |k − k | 2 = 2k(1 − cos ϕ). Hence the square is proportional to 1/(1 − cos ϕ) and goes like 1/ϕ 2 for small angles. The rest of the expression results from our cutoff at the vortex edge (r = R) and serves to suppress dσ/dϕ for ϕ < 1/kR. Integrating over angles therefore yields a total cross section that is linear in kR for large kR. Recall that the transport cross section includes an additional weighting factor of (1 − cos ϕ). This precisely cancels the leading functional dependence of dσ/dϕ and results in a transport cross section that approaches a constant for large kR. Our numerical results for σ (see Fig. 4) are therefore quite logical. Thus, the Born approximation is sufficient to capture the qualitative features of our transport cross sections.

XII. CONCLUSIONS
In this paper, we considered the scattering of thermally-excited quasiparticles from a single magnetic vortex in a d-wave superconductor. The scattering effect of the vortex can be divided into two contributions, one due to the superflow circulating about the vortex and another due to the Berry phase acquired by the quasiparticle upon circling the vortex. Both effects are present in the Bogoliubov-de Gennes equation for a single vortex. By applying a singular gauge transformation, we have explicitly separated them, the superflow resulting in an effective non-central potential and the Berry phase generating antiperiodic boundary conditions for the quasiparticle wave function. In a separate paper 12 , we have studied the Berry phase contribution in the absence of the superflow effect, essentially the Aharonov-Bohm scattering of a massless Dirac fermion. Here, by adopting periodic (rather than antiperiodic) boundary conditions, we have considered the superflow contribution in the absence of the Berry phase effect. The resulting problem is that of a massless Dirac fermion scattering from a non-central potential (non-central because the gap node about which we linearize is shifted away from the origin of momentum space). We solved this problem by calculating the eigenstates inside and outside the vortex, building an incoming plane wave and an outgoing radial wave, and computing the resulting differential cross section. We plotted the size and shape evolution of the differential cross section as a function of the ratio of the magnetic length R to the de Broglie wavelength 1/k. Intregrating, we computed the transport cross section and showed that it saturates to a constant value for large kR. This result was previously used in Ref. 13 to compute thermal transport coefficients. Here, we have provided the details of its derivation. We have also provided an alternate derivation, and a closed-form solution, valid within the Born approximation. These results, taken together with the results of Ref. 12, shed significant light on the nature of the scattering of massless Dirac quasiparticles from magnetic vortices. What remains is a complete analysis, taking into account both the superflow and Berry phase contributions and the interference between them.

XIII. ACKNOWLEDGMENTS
We would all like to thank Sasha Abanov, Patrick Lee, Zlatko Tesanovic, and Ashvin Vishwanath for very helpful discussions. This work was supported by the NSF under grant No. DMR-0605919. M.K. was also supported by the NSF under grant No. DMR-0906866. S.G. was also supported by the DOE under grant no. DE-FG02-09ER16052.
The desired solution has the form In the position basis one could write it conveniently as We will express the operator (E − H 0 ) −1 in the position basis as where < x | p >= 1 2π e i p· x (A8) Therefore, simplifying Eq. (A44) further we can write the full cross section as dσ dϕ = π 32k 1 sin 2 ( ϕ 2 ) (1 − J 1 (2kR | sin ϕ/2 |) kR | sin ϕ/2 | ) 2 (cos(ϕ + θ) + cos(θ)) 2 (A47) The above quasiparticle scattering cross section is the contribution from one node only. We will now calculate the contributions from the other three nodes see (Eq. (90)). θ and ϕ are the Node 1 angles and we calculate the rest of the node contribution with respect to this node.