A finite element method for the quasiclassical theory of superconductivity

The Eilenberger-Larkin-Ovchinnikov-Eliashberg quasiclassical theory of superconductivity is a powerful method enabling studies of a wide range of equilibrium and non-equilibrium phenomena in conventional and unconventional superconductors. We introduce here a finite element method, based on a discontinuous Galerkin approach, to self-consistently solve the underlying transport equations for general device geometries, arbitrary mean free path and symmetry of the superconducting order parameter. We present results on i) the influence of scalar impurity scattering on phase crystals in $d$-wave superconducting grains at low temperatures and ii) the current flow and focusing in $d$-wave superconducting weak links, modeling recent experimental realizations of grooved high-temperature superconducting Dayem bridges. The high adaptability of this finite element method for quasiclassical theory paves the way for future investigations of superconducting devices and new physical phenomena in unconventional superconductors.


I. INTRODUCTION
The quasiclassical theory of superconductivity [1][2][3] allows the study of a wide range of phenomena in both conventional and unconventional superconductors on a mesoscopic scale. Solving the underlying transport equations is, however, a highly involved numerical problem, especially when a self-consistent determination of the self-energies is required. We propose here a finite element method (FEM), specifically a discontinuous Galerkin (DG) approach, that can be used to solve the underlying transport equations of the full Eilenberger quasiclassical theory in realistic superconducting systems in two or three dimensions. The DG method of solving transport equations in dimension D ≥ 2 was first proposed in the 1970s for the neutron transport equation in the context of nuclear reactors [4][5][6] but has since been applied to (in)compressible flow dynamics, chemical transport, and many other areas of physics [7].
Such a FEM offers several advantages over finite difference methods used in the past [8,9]. Firstly, the solution strategy introduced here is directly applicable to systems in one, two, or three dimensions, up to a slight adaption of the boundary conditions depending on the dimension. Secondly, it provides better adaptability of the discretization to the geometry under investigation as well as the degree of approximating functions which improves convergence [10].
For dirty conventional superconductors, where the mean free path is much shorter than the superconducting coherence length ξ 0 , the diffusive approximation due to Usadel is valid. An implementation of a FEM for the underlying Usadel diffusion equations has been reported [11], and used to analyze experimental results [12] as well as to investigate new physics in higher-dimensional structures [13]. The method we present here is valid for superconductors with arbitrary mean free path and orderparameter symmetry, and thus extends the possibility of finite-element analysis to a wide range of superconducting materials and phenomena. This paper is organized as follows. In Sec. II, we review the underlying equations of the quasiclassical theory of superconductivity and present the reformulation in terms of a discontinous Galerkin method. This is followed by results on two example problems in Sec. III. Firstly, we investigate the effect of scalar impurites on phase crystals in a closed system. Secondly, we model a current biased superconducting Dayem bridge as an example for current flow and focusing in a geometry with open boundaries. We conclude with an outlook on possible further applications of the method in Sec. IV.

A. Quasiclassical theory
The core of the quasiclassical theory of superconductivity is the Eilenberger equation, i v F · ∇ǧ(p F , R, ε) + ετ 31 −ȟ(p F , R, ε),ǧ(p F , R, ε) = 0, (1) for the quasiclassical Green's functionǧ(p F , R, ε) together with a normalization conditioň g(p F , R, ε) 2 = −π 2 . (2) The above form is for the time-independent steady state that we assume in this paper. In this case, the propagators only depend on momentum direction on the Fermi surface, p F , spatial coordinate, R, and energy, ε. In Eq. (1), v F is the Fermi velocity, the [A, B] denotes a commutator between matrices A and B, aˇmarks Keldysh-space matrices,τ 3 is the third Pauli matrix in Nambu (particle-hole) space, indicated by aˆ, andȟ is the self-energy matrix. The three elements of the Keldysh-space matrixǧ arě where the spectrum of the system is determined by the retarded (advanced) componentĝ R (ĝ A ), while information about the occupation of states is contained in the Keldysh componentĝ K . All three components are in turn matrices in Nambu space that we will detail below.
To shorten the notation, we often drop the explicit dependences on p F , R, and ε. The selfenergy matrixȟ has the same form asǧ in Keldysh space with the Nambuspace componentŝ This set of equations was first derived by Eilenberger [1], and separately Larkin and Ovchinnikov [2], and generalized to non-equilibrium by Eliashberg [3]. Instead of solving Eq. (1) directly, it is advantageous to use parametrizing functions forǧ that guarantee that the normalization condition, Eq. (2), is satisfied. One well-established choice is a parametrization in terms of coherence amplitudes γ,γ [14][15][16][17] and generalized distribution functions x,x [18,19]. The coherence amplitudes determine the retarded component where The expression for the advanced elementĝ A is analogous toĝ R . The Keldysh component also involves the distribution functions x andx, Generally all elements ofĝ R,A,K are matrices in a space of internal degrees of freedom such as spin, but in the following we assume a spin-singlet superconductor for simplicity. Then, γ can be written as γ = γ singlet iσ 2 -where γ singlet is a scalar function and σ 2 is the second Pauli matrix in spin space -while x is proportional to the unit matrix in spin space. One central symmetry of the theory is particle-hole conjugation, expressed as Starting from Eq. (1), a set of coupled equations for the parametrizing functions can be derived. In the steadystate, the equations for γ(p F , R, ε) and x(p F , R, ε) read Equations for the remaining amplitudesγ R/A and distributionx can be obtained via the tilde symmetry in Eq. (7). Note that both equations are transport equations along a transport direction v F , indicated by the directional derivative v F · ∇. Solutions of these differential equations have to be found along semi-classical trajectories determined by v F , a starting point, and an end point. Along such a trajectory, Eq. (8) is a Riccati equation that is well-studied in literature [20]. Observables are then determined by an average over all possible momentum orientations v F on the Fermi surface. For our twodimensional case we use a circular Fermi surface, where the orientation is fully determined by the angle to the x axis, ϕ F . The average of an observable or self-energy A is then After solving Eq. (1), all self-energies have to be recalculated until self-consistency is achieved. In the present paper, we usě whereȟ mf is the mean-field order parameter, given by (12) Here, η d is the basis function for the respective orderparameter symmetry. For the d x 2 −y 2 order parameter, which we will consider in this paper, we can choose η d = cos [2 (ϕ F − α)] where the angle α specifies the misalignment of the main crystal axis to the geometric axis. The self-energyȟ s describes scalar impurity scattering. Assuming an average dilute impurity concentration n i , the impurity self-energy is found from the t-matrix equation in the non-crossing approximation [21]: For scattering that is isotropic in momentum space with an s-wave scattering potential u 0 the elements ofť satisfy the equationst The two parameters of the model, the isotropic impurity potential u 0 and impurity concentration n i , can equivalently be written as a scattering energy Γ u = n imp /πN F and a scattering phase shift, δ 0 = arctan(πN F u 0 ). They can be combined into Γ ≡ Γ u sin 2 δ 0 , the so-called pair-breaking energy, which determines the normal-state mean free path = v F /(2Γ). We restrict ourselves here to the weak-scattering Born limit, or the opposite strong-scattering unitary limit, Solving Eqs. (8)-(9) allows to construct the full quasiclassical Green's functionǧ that solves Eq. (1), and can be used to calculate physical observables. For example, we find the momentum-resolved density of states which determines the full density of states, Here, N F is the normal-state density of states per spin at the Fermi level. The charge current density is obtained via We use j 0 = ev F N F k B T c as the unit for charge current. If both the system and all connected reservoirs are in equilibrium the Keldysh Green's function at temperature T simply readsĝ By replacing the energy ε with imaginary-axis Matsubara frequencies iε n = iπk B T (2n + 1) in Eq. (8) we can then solve for the Matsubara coherence function γ M . For better numerical performance alternatives such as the Ozaki poles of a continued-fraction expansion [22] can be used. This is the case for the results presented in this paper.

B. Trajectory method
Self-consistent numerical solutions of Eqs. (8) and (9) require a discretization of the self-energies on a grid, and similarly discretizing the trajectories in "steps". In (quasi-) one-dimensional systems it is always possible to choose the two to be commensurate. One possible choice is to model the energy landscape and the coherence functions as piecewise constant in between gridpoints. In this case, analytic solutions to Eq. (8)-(9) in a region of constant selfenergies can be found. This leads to a stepping method where the solution for the equation of motion, Eq. (8), is found along the trajectory by stepping along the trajectory step-by-step from one region to another.
In contrast, the steps cannot be made commensurate with the self-energy grid for all momentum orientations at once in two or three dimensions, as sketched Fig. 1. A stepping method is still possible and has been implemented for two-dimensional systems [9]. The discrepancy  between the self-energy grid and the trajectory steps require frequent interpolation from the self-energy grid to the "stepping grid". Here, we describe an alternative solution strategy that relies on a finite element method, more specifically a discontinuous Galerkin method. The solution to the differential equation is then defined on the same grid as the self-energies which avoids grid interpolation. The underlying assumption of finite element methods is that a given geometry can be decomposed into a collection of elements or cells, as indicated in Fig. 2. Note that we use a two-dimensional geometry and triangular cells for illustration, but the method introduced here can be directly applied to three-dimensional models as well. The original differential equation then has to be translated to its weak form. We will focus on the Riccati equation for γ R but will drop the superscript R in the following. The derivation is largely identical for the advanced function γ A . For ease of notation, we focus on the spin-degenerate case now which will result in a scalar equation for the coherence amplitude γ. The generalization to the full two-by-two spin structure follows along identical lines but will produce a coupled system of equations for four unknowns, with equally many independent test functions.
We start from Eq. (8) that can, in the scalar case, be rearranged to read Multiplying both sides by a, for now unspecified, test function φ(R) ≡ φ, and integrating over the entire domain Ω gives We now assume to have a triangulation T , meaning a collection of triangles T j that satisfy Ω = ∪ Tj ∈T T j . The integration is then split into a sum of integrals over each triangle, Performing a partial integration on the first term on the left hand side in Eq. (24), we arrive at where n j is the outer normal for ∂Ω j , the boundary of a given cell. For triangular cells each element boundary ∂Ω j consists of three edges. Each edge, in turn, is either part of the geometric boundary ∂Ω, or part of the collection of internal edges τ . The geometric boundary ∂Ω can be further decomposed into the inflow boundary ∂Ω − and the outflow boundary ∂Ω + , defined via The three sets are sketched in Fig. 2(b). In a discontinuous Galerkin method the unknown function γ as well as the test function φ can have different values on the two sides of a given edge. Each internal edge τ j ∈ τ is shared by two cells and thus integrated over twice in the sum over all the boundary integrals. It can be shown [23] that summing this double integration over each edge allows to rewrite the first term on the left-hand side in Eq. (25) as where the summations on the right-hand side are over individual edges rather than closed boundaries of each triangle. The different brackets in the first term on the right-hand side denote the jump [. . . ] and average {. . . } of a function over an edge shared by two cells 1 and 2. They are defined for vectors a and scalars φ as with the index denoting the value in the respective cell, and the outward-pointing normal vector n i in each of the two cells, see Fig. 3. We note that in the literature the jump operator is sometimes written including the normal vector, i.e., [a] = [a · n]. For edges on either the inflow or outflow boundary the function only has values on one side so no average or jump of the function values appears in the other two terms on the right-hand side of Eq. (28). Up until this point the derivation is identical for the retarded and advanced functions, γ R and γ A , but start to be different in the following. For the retarded function, γ R , a boundary condition γ B has to be provided on the inflow boundary ∂Ω − . We address how these starting values are obtained in Sect.II D. On the outflow boundary ∂Ω + and on internal edges τ , the values of γ R are unknown and determined in the solution step. This corresponds to solving Eq. (8) for the retarded function in the direction of v F with a boundary value at the start point of a given trajectory. For the advanced function γ A , ∂Ω − and ∂Ω + are swapped, corresponding to the stable integration direction opposite to v F in the stepping method. The remaining functionsγ R,A are similarly swapped with respect to their non-tilde counterpart. The naming of ∂Ω − (∂Ω + ) as inflow (outflow) boundary is thus only descriptive for γ R andγ A but follows the naming conventions established in literature.
Lastly, we note that for the imaginary-frequency Matsubara coherence function γ M , the propagation directions and flow boundaries correspond to the one for γ R (γ A ) for positive (negative) Matsubara poles.
As discussed in literature [24], better convergence and stability can be achieved for this weak form by replacing the average operator in the first term in Eq. (28) by the so-called upwind value, denoted by a subscript u and given by The upwind value enforces a propagation of the function values across cell edges in the transport direction v F , corresponding to the propagation along trajectories in the stepping method. A corresponding downwind value with swapped signs in the inequalities in Eq. (31) has to be used for the advanced function γ A . This is similar to the exchange of ∂Ω + and ∂Ω − for the boundary values. The final weak form for the retarded function γ R thus becomes Eq. (32) is of the form where f does not depend on the unknown function γ while L does. In the FEM language, L is the bilinear form and f the linear form of the weak formulation of Eq. (8).
An approximate weak solution γ R w to Eq. (32) can be found using a discontinuous Galerkin method. Such a method assumes that the approximate weak solution γ R w can, within each cell, be written as a sum of polynomials of finite order k. Depending on k, there are N j independent degrees of freedom, or nodes, for a given cell T j . In two dimensions, the common choice of k = 1 (linear functions) or k = 2 (quadratic functions) requires three and six nodes, respectively. For each of the N j nodes, there is an associated polynomial φ i . In this paper we use simple Lagrange polynomials that satisfy Using this cell-wise basis we can write where a j,i are expansion coefficients or weights. The function values on the nodes then fully determine a j,i and by Eq. (35) the function everywhere in the individual cell. The main difference to a continuous Galerkin method is that the values of a j,i on nodes shared between neighboring cells are independent, the global function can thus be discontinuous across cell edges. Such a discontinuous Galerkin method has been found to give better convergence for the transport equations and allows for higher mesh adaptability, see also the discussion in Ref. 6, 10, and 25. We expand the self-energies in the same basis with weights ∆ j,i and Σ j,i , while the test function φ used in Eq. (32) can be chosen to have the same form but with unit weights, φ = i φ i (x).
Inserting Eq. (35) and the basis functions into Eq. (32) gives an equation system for the unknown coefficients a j,i for each cell. Renumerating a j,i → a i and combining all the local equations, we get a collection of scalar equations of the form where i = 1, . . . , N is an index over the all N nodes in our domain. In two dimensions, a triangulation with N T triangles and N j nodes per cell will have a total number of N = N T N j nodes. Here, the tensor Q ijk is due to the quadratic term γ∆γ in Eq. (32), while the remaining linear parts of L(γ, φ) are combined into a matrix P ik . The nonlinearity prevents a solution of Eq. (36) by matrix inversion. Instead, the vector a is determined as an (approximate) zero of the residual vector r, with components where again i = 1, . . . , N as in Eq. (36). This approximate zero is found by iterative minimization. We find that a numerically efficient strategy is Newton iteration using the Jacobian matrix J ij (a, φ) of Eq. (32), which can be calculated from the residual via the Gateaux derivate Starting from a guess a 0 , one iterates the solution vector via Newton steps a n+1 = a n + (δa) n , where (δa) n is the solution to a linear equation system obtained by combining all J ij (a n , φ)(δa j ) n = −r i (a n , φ) for i = 1, . . . , N .

D. Boundary values
As outlined in Sect. II C, the weak form in Eq. (32) contains boundary values γ B , in the case of the retarded function on the inflow boundary ∂Ω − . Note that the inflow boundary will be different for different orientations of v F . These boundary values correspond to the starting value of the coherence function at the start point of a trajectory in a stepping method.
We treat the entire boundary ∂Ω, including the inflow boundary, as a collection of planar interfaces of finite length. Scattering at such atomically sharp interfaces cannot be described within quasiclassical theory itself, the theory has to be supplemented by external boundary conditions. Using scattering theory, such conditions have been derived for specularly scattering interfaces [18,19,[26][27][28][29] where the momentum parallel to the surface is conserved in the scattering process. The boundary conditions in the spin scalar case is at a given interface are expressed in terms of the normal-state scattering matrix where R and D are, respectively, the reflectivity and transmittivity of the interface satisfying R + D = 1. The boundary condition for γ R B is then of the form where γ R in,i=1,2 are the incoming function inside and outside the system. Within the FEM nomenclature such boundary conditions are referred to as Dirichlet boundary conditions. In a discontinuous Galerkin method, they are weakly enforced through the last term on the righthand side in Eq. (32). On a given segment the entire function is determined by the values on the nodes, hence it is sufficient to apply the boundary condition at the nodes only.
At a fully reflective segment, R = 1, specular scattering relates the outgoing function for an angle ϕ F to the function incoming from within the system for a different angle ϕ F . On a boundary segment with an outwardpointing normal n that spans an angle β n to the x-axis, we find This relation holds in systems with a two-dimensional Fermi surface where the orientation of v F can be parametrized by one scalar parameter only. In three dimensions an analogous relation for two scalar parameters has to be used. The weak form in Eq. (32), however, is valid in two and three dimensions if one replaces cells and edges with appropriately chosen polyhedra and surface planes.
Starting from a bulk guess as an incoming function, we use the incoming functions at an iteration n to obtain the outgoing function in the next, n + 1, iteration. At a fully transparent segment, D = 1, the outgoing function solely depends an incoming function from a reservoir that is not directly simulated. This reservoir can for example be a superconducting contact, or a normal-metal reservoir. We model both cases by associating a virtual "reservoir" point with each node on the given boundary segment, see Fig. 4.
We model a superconducting contact by using a bulk coherence function γ bulk as the incoming coherence function from the reservoir. On each virtual point, the order parameter and superfluid momentum p S are then iterated such that we enforce a certain current density across the respective boundary segment. In the case of a normal-metal reservoir, a voltage bias (temperature bias) can be modelled by fixed electrochemical potential (temperature) of the reservoir so no update is required. This induces a nonequilibrium occupation at the system edge that propagates in the system. For further details, we refer to our recent publications on nonequilibrium setups in (quasi) one-dimensional models [30,31]. For intermediate values of the transparency, D ∈ (0, 1), the boundary condition Eq. (40) gives a weighted combination of the two scenarios described above.

E. The non-equilibrium distribution x
A weak form for the non-equilibrium distribution x can be obtained in a largely identical procedure as outlined for the coherence amplitude γ in Sect. II C. After partial integration over the term containing the derivative v F · ∇x of splitting of the resulting boundary integrals has to be performed. The specification of boundary values and flow direction for x (x) then follows that of γ R (γ A ). The main difference to the coherence amplitude is that Eq. (9) is clearly linear in x. Hence the solution of the resulting matrix equation can be obtained through matrix inversion rather than the iterative solution of a nonlinear problem.

F. Selfconsistency & Numerical aspects
A self-consistent solution of the Eilenberger equation requires an iterative procedure of solving the underlying transport equation followed by an update of all selfenergies. In equilibrium, this leads to the following recipe for the retarded functions: 1. Provide a current guess of all self-energies inĥ R , see Eq. (4), and boundary values γ R B andγ R B .
2. Find weak solutions γ R w to Eq. (32), and similarlỹ γ R w , for all energies ε and momentum angles ϕ F .
3. Update all self-energies through Eqs. (12) and (16) or (17). 4. Update all boundary values γ R B andγ R B as described in Section II D Non-equilibrium situations also require weak solutions x w andx w and updates of the Keldysh self-energiesĥ K . The above procedure has to be repeated until a self-consistent solution is found, signalized for example by convergence of the self-energies up to a desired accuracy. Only then is charge-current conservation guaranteed everywhere in the system. For a given guess of all self-energies, the solution of Eq. (32) is an independent problem for each energy ε and momentum orientation ϕ F . This makes the problem highly parallelizable.
The solutions for all energies ε and angles ϕ F only need to be combined to obtain a new guess of the self-energies.
For the FEM solution step we use the open-source package Gridap [32], while the meshes are created in gmsh [33]. Detail on the specific meshes used for the results presented in this paper can be found in a supplementary material online.

III. RESULTS
To illustrate the strategies and the advantages of solving the Riccati equations with a FEM together with the self-consistency equations for the order parameter and scalar impurity self-energies, we present two examples. Both involve two-dimensional d-wave superconductors modeling for instance a single superconducting plane of a high-temperature superconductor. The first example is a closed system, a square island (see Fig. 5), where the superconducting crystal axes are rotated 45 • relative to the main axes of the square. For this orientation of the order parameter relative to the edges, below a phase transition temperature T * ≈ 0.17T c , time-reversal symmetry has been predicted to be broken in an unusual way [9,[34][35][36][37][38]. The characteristic of this state, referred to as a phase crystal [36], is a non-trivial, structured ordering of the superconducting phase χ(R) resulting in patterns of current flow consisting of loops near the edges. In the present example we show the correctness of the FEM solution by studying the detrimental influence of scalar impurity scattering within a self-consistent homogeneous scattering model. In the second example we study an open system, a current biased d-wave superconducting bridge with an inhomogeneous density of impurities, modeling experiments on grooved Dayem bridges [39,40]. In this case we also demonstrate how to self-consistently enforce a current boundary condition at the source and drain leads and current conservation across the bridge. This nicely illustrates the effect of current focusing and we predict the spatial variation of the superconducting phase over the bridge. In both example setups we assume that scattering at non-transparent boundaries is specular, as discussed in Sect. II D.

A. Effect of scalar impurities on phase crystals
Consider the square d-wave superconducting island in Fig. 5. When the crystal axes are misaligned by 45 • relative to the surface norms, the d-wave order parameter is suppressed to zero at the edges, see Fig. 5(a). This reflects the formation of zero-energy surface Andreev bound states [41][42][43]. At temperatures below T * ≈ 0.17T c a second-order phase transition was recently predicted, where time-reversal symmetry is broken and a non-uniform state appears with circulating currents near the edges [9]. The influence of mesoscopic roughness of the boundaries [34] and magnetic field [35] have been investigated earlier. Within a tight-binding model, also the influence of band structure was investigated [37]. Recently, an investigation of strong correlations within a Gutzwiller approximation preventing double occupancy on each tight-binding site was reported [38]. In the latter study it was shown that Anderson disorder suppresses T * , but surprisingly that strong correlations weakens the influnce of this disorder. In our example, we study the suppression of T * within a differ- ent model for the scalar impurities. In our quasiclassical treatment we self-consistently compute the impurity self-energy in a homogeneous scattering model for two types of disorder: the Born limiting of weak scattering, Eq. (16), and the unitary limit of diverging scattering potential strength, Eq. (17).
In both cases, larger impurity concentrations suppress the current flow up to a critical value Γ * at which the phase crystal is no longer energetically favorable and a uniform superconducting state is formed instead. Examples of current flow patterns for two temperatures, T = 0.03T c and T = 0.07T c , in the unitary limit for a larger value of the impurity concentration is shown in Fig. 6. As compared to the almost clean case in Fig. 5(c), for such impurity concentrations, such as Fig. 6 (normal state mean free path ≈ 22ξ 0 ), the pattern remains largely unchanged at the lower temperature while the absolute value of the current is reduced. At the higher temperature for the same impurity concentration, a pattern with a reduced number of loops is observed, see Fig. 6 As a measure of the total amount of current flow in the system we use the integrated quantity This average measure depends on temperature, impurity concentration, and type of impurities (Born or unitary), as summarized in Fig. 7. For the relatively clean case (blue circles), the suppression of j 2 with temperature leads to an estimate of the transition temperature T * . For T < T * the phase crystal is energetically favorable [9]. The corresponding temperature dependencies of j 2 for increasing impurity concentration Γ is shown in Fig. 7(a) for the Born limit and Fig. 7(b) for the unitary limit. Generally, for the same temperature and impurity concentration the currents are smaller for scattering in the Born limit. It is well-known that the surface Andreev bound states are less broadened in the unitary limit than in the Born limit [44]. As a result of the larger broadening of the Andreev bound states in the Born limit the spontaneous currents are energetically less favourable. For a temperature T = 0.1T c the currents are well developed in the clean case, with a typical wavelength for the current loops of ∼ 12ξ 0 (two loops with opposite circulations [9]). With increasing Γ, first the pattern changes to a single loop and then vanishes at this temperature at a mean free path of order ∼ 30ξ 0 in both the unitary and Born limits. These results indicate that the phase diagram in Γ − T -space can be rich and should be investigated in further detail, but this is beyond the scope of the present paper.
B. Focussing of current flow in a Dayem bridge As an example of an open system with current flow we consider a d-wave superconducting bridge, inspired by recent experimental realizations of Dayem bridges [39,40]. Fig. 8 shows the model for such a bridge. It consists of two superconducting reservoirs smoothly connected to a narrow channel of reduced width. The d-wave order parameter is in this case aligned with the bridge, such that the lobes are along x-and y-directions. Additionally, a groove can be etched into the channel. To model this groove we include a position-dependent impurity concentration of Gaussian shape in the center of the channel, In the results presented here we use x center = 45ξ 0 , at the center of the channel, and w Groove = √ 0.2ξ 0 . With these parameters the impurity concentration has a smooth Gaussian profile over a characteristic scale of The reduced mean free path in the groove results in a suppression of the order parameter by around twenty percent in the dirty region in the center of the channel already in equilibrium, as compared to the case without an impurity concentration profile, see in Fig. 9. Note also that the order parameter is, in both cases, partially suppressed at the tilted edges in the transition area between the wide contact and the narrow channel a result of the partially pair-breaking nature of the walls there.
We now enforce a constant boundary current density of j x = j b = 0.01j 0 at the left and right edges of the system. Upon a self-consistent calculation, the current that is flowing in the x direction from the left reservoir gets spatially focused to flow into the narrow channel.
By integration of Eq. (42) and application of Gauss' theorem in two-dimensions we can reinterpret the condi-tion of local current conservation as int C ∇ · j dΩ = C n · j ds = 0, (45) where C is a closed contour within our system and int C the enclosed area. An example for such a contour is shown in Fig. 9(a). We always include the left system edge as well as the upper and lower edges of the geometry. The last part of the contour is then a line parallel to the y-axis at different distances to the left edge. We can then integrate j x (x, y), the x component of the current, over the local height W (x) of the structure, Since no current can flow out of the reflective walls, the integrated current flow through the left and right contour edge, indicated as dashed lines in Fig. 9(a), has to be conserved throughout the system. In the results presented here, this is the case up to a relative error of δI x < 2 · 10 −2 I x (0). A self-consistent result of the spatially redistributed current flow is shown Fig. 10. The focusing of the current from the wide lead into the narrow channel is clearly visible. For the same imposed boundary current the flow pattern is basically identical for both grooved and nongrooved channels since the y-integrated current, Eq. (46), is conserved across the weak link in both cases. However, the suppression of the order parameter in the grooved Dayem bridge leads to an increased phase gradient in the channel. This can be more clearly seen by plotting the phase χ and x component of the superfluid momentum along a line cut in the x direction for constant y = W l /2, in the middle of the channel. Both quantities are shown in Fig. 11 for a non-grooved (blue dashed lines) and a grooved Dayem bridge (orange solid lines). In order to carry the same current in the x direction the phase gradient is enhanced in the groove region to compensate for the reduced order parameter. For increasing values of the boundary current, the suppression of the order parameter amplitude and the enhancement of the superfluid momentum become more and more pronounced, see Fig. 12.
Beyond the highest boundary current j b > 0.025j 0 superconductivity breaks down locally in the channel in the sense that the order parameter is locally suppressed to zero during the self-consistency iteration. After that we find no self-consistent stationary solution for the order parameter in the weak link.

IV. DISCUSSION AND OUTLOOK
In this paper, we have presented a reformulation of Eilenberger's quasiclassical theory of superconductivity  11. Variation of the superconducting order parameter along a line cut in the center of the channel, y ≈ W l /2, with Γ peak = 0 (dashed blue) and Γ peak = 5 · 10 −2 πkBTc (solid orange).
in terms of a discontinuous Galerkin method. We applied the method, firstly, to study the influence of scalar impurity scattering on phase crystals in d-wave superconductors, and, secondly, to investigate the current flow through a geometric constriction in the form of a Dayem bridge.
In the first case, we find that scalar impurity scattering suppresses the spontaneous flow patterns of phase crystals. As a result the average flow gets suppressed compared to a clean d-wave superconductor and the characteristic temperature T * , where the phase crystal is destroyed, is reduced. For the same impurity concentration the spontaneous flow is lower for impurities in the Born limit compared to the limit of unitary scattering. We attribute this to the increased broadening of the surface Andreev bound states in the Born limit which makes spontaneous flow energetically less favourable.
Our second example shows the strength of the FEM method in the application to open superconducting systems in the presence of current flow. Generally, the constriction in the Dayem bridge geometry leads to a increased phase gradient in the narrow channel. In the case of a grooved Dayem bridge, the reduced mean free path in the Groove requires a strong increase in the phase gradient within the groove region in order to carry the same current. As a result, the groove can lead to a quick suppression of superconductivity with increasing current although the flow itself is non-dissipative.
The method presented here gives a powerful and versatile numerical technique to study conventional and unconventional superconductors with arbitrary mean free path in two-and three-dimensional geometries. The adaptability of a FEM to arbitrary geometries allows the investigation of realistic device geometries and experimental setups. An extension to the case of full spinstructure of the quasiclassical propagators follows identical steps as discussed here and would allow to study for example the effects of spin-orbit interactions and other unconventional superconductors. Another extension, relevant for, e.g., the situation in the Dayem bridge when we can not find self-consistent stationary solutions at high imposed currents, would be to include time-dependence. ACKNOWLEDGMENTS We thank J. Wiman and O. Shevtsov for early discussions on the possibility of applying the finite element method within quasiclassical theory, and acknowledge the Swedish research council for financial support. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

APPENDIX: COMPARISON OF NUMERICAL PERFORMANCE
To benchmark the numerical performance of the DG method presented in this paper, we use a simple example in one dimension. Specifically, we assume a clean system with an order-parameter profile of the form where ∆ 0 = 1.5k B T c , x mid = 7.5ξ 0 , and let x ∈ [0, 15ξ 0 ]. This profile is similar to that close to a fully reflective, pair-breaking interface to a d-wave superconductor. For this scenario, we compare 1. the "stepping"method, see Sect. II B, 2. an implicit midpoint method, a well-known finitedifference scheme, see also Appendix C in Ref. [45], 3. the DG method presented here, for different order k of the approximating polynomials.
In all cases, we assume a uniform grid of N points (or nodes), which gives rise to a step size or cell size h = L/N = 15ξ 0 /N . As a starting value, we assume γ(x = 0) = 0 at the start of the trajectory.
In the absence of an analytic solution for the orderparameter profile in Eq. (48), we take the numerical solution of the stepping method for a very dense grid, using N max = 51200 points, as the reference solution. As a measure of the numerical error, we use the summed difference where we choose the x i such that they are included in the grid for all values of N . A plot of the error scaling as function of N , normalized to δγ(N min ) the error for the grid with the lowest number of points N min = 25, is shown in Fig. 13. For our benchmark problem, we find that the error δγ(N ) is of order O(1/N ) for both the stepping method and a corresponding DG method with piecewise constant polynomials (k = 0) as approximating functions. Once we use linear approximating functions (k = 1), we get an improved scaling behavior of order O(1/N 2 ), which we also find in the midpoint method. In both approaches, this improved error scaling is only found if the self-energy profile is also modeled as a linearly varying function in a given step length or cell. For a DG method with quadratic polynomials (k = 2), we find an error scaling of order O(1/N 3 ) (not shown in Figure 13). From the literature on DG methods we expect that for approximating polynomials up to order k, the error scaling is of order O((1/N ) k+1 ) for linear hyperbolic problems on uniform grids [10]. While some literature exists on DG methods for non-linear problems [46], it appears that a more thorough analysis of the error scaling for the nonlinear Riccati equation is an open question.
The computational time needed to solve this example in one dimension scales linearly with the number of nodes N and is thus of the same order of magnitude for all presented methods. In two dimensions, the midpoint method should scale as N 2 on a square geometry while the equivalent FEM is predicted to scale as N d 2 , where d is the so-called band width of entries around the diagonal [25]. Since d is determined by the upwind or downwind coupling of neighboring cells, d can be made small compare to N by numerical optimization.