General solution of 2D and 3D superconducting quasiclassical systems: coalescing vortices and nanoisland geometries

An extension of quasiclassical Keldysh-Usadel theory to higher spatial dimensions than one is crucial in order to describe physical phenomena like charge/spin Hall effects and topological excitations like vortices and skyrmions, none of which are captured in one-dimensional models. We here present a numerical finite element method which solves the non-linearized 2D and 3D quasiclassical Usadel equation relevant for the diffusive regime. We show the application of this on three model systems with non-trivial geometries: (i) a bottlenecked Josephson junction with external flux, (ii) a nanodisk ferromagnet deposited on top of a superconductor and (iii) superconducting islands in contact with a ferromagnet. In case (i), we demonstrate that one may control externally not only the geometrical array in which superconducting vortices arrange themselves, but also to cause coalescence and tune the number of vortices. In case (iii), we show that the supercurrent path can be tailored by incorporating magnetic elements in planar Josephson junctions which also lead to a strong modulation of the density of states. The finite element method presented herein paves the way for gaining insight in physical phenomena which have remained largely unexplored due to the complexity of solving the full quasiclassical equations in higher dimensions.

Scientific RepoRts | 6:22765 | DOI: 10.1038/srep22765 limit) and must be supplemented by appropriate boundary conditions. Focusing on the diffusive limit, as it is often the experimentally relevant one, a variety of options are available depending on the physical situation at hand. In the simplest case of perfectly transparent interfaces, the Green function is taken as continuous across the interface. This is clearly an idealized scenario and the more realistic Kupriyanov-Lukichev 13 boundary conditions describe an interface in the tunneling limit where there exists a substantial interface resistance. Boundary conditions for an arbitrary interface transparency were developed in ref. 14. When the interface has magnetic properties, either because of an intrinsically thin magnetic layer inserted between e.g. two metals or superconductors or if one of the regions separated by the interface is magnetic on its own, one must use spin-dependent boundary conditions. Pioneered in refs 15,16, these were brought to a more tractable form by Cottet and co-workers in the diffusive limit 17 . However, up until recently there existed a knowledge gap in terms of how to describe strongly polarized magnetic interfaces in quasiclassical theory. Eschrig et al. solved this problem in ref. 9.
It is clear that the development of a numerical routine that is able to solve the quasiclassical Keldysh equations in higher dimensions than 1D will be of great value in terms of studying a vast number of physical phenomena, including various types of Hall effects, spin swapping, and topological excitations such as magnetic skyrmions and vortices. None of these phenomena can be captured in an effective 1D model. Furthermore, the ability to handle complex higher dimensional geometries numerically allows for the modeling of systems which are more closely related to experiments. For instance, superconducting nanoisland systems and vortices in mesoscopic structures have received much attention experimentally [18][19][20][21][22] . These systems require not only solution in 2D or 3D, but also the description of non-trivial geometries within the numerical framework. Such solutions have been investigated using the Ginzburg-Landau formalism in the context of flux patterns and vortex states in superconductors [23][24][25] . The ability to aid experiments with numerical routines that are both geometry and dimension independent would be highly beneficial to their study. Nevertheless, explicit solutions of the full quasiclassical equations in two dimensions have rarely been reported 26,27 . In the linearized regime, corresponding to a weak proximity effect, several works have considered the 2D solution of the Usadel equation [28][29][30][31][32] . Motivated by this, we report as the main result of this paper the description of a finite element method that we have developed which is capable describing mesoscopic systems in 2D and 3D using quasiclassical theory without any linearization. As far as the authors are aware, this is the first work to solve the Usadel equations in 3D. After going through the details of this method, we show its application to three model systems. One of our main findings is that in a 2D Josephson junction exposed to a magnetic flux, it is possible to control not only the geometrical array in which superconducting vortices arrange themselves, but it is also possible to cause coalescence and thus tune the number of vortices. In addition, we show that the supercurrent flow through planar junction geometries can be tailored by the magnetization pattern and strength and also spatially modulates the proximity-induced density of states, which can be probed by STM-measurements. We organize our presentation as follows. First, we introduce the system of coupled NLDEs that define the central equations in quasiclassical theory. The finite element method solving these equations in 2D and 3D is described in detail in the next section. We proceed to show the application of this method to three different hybrid structures where a superconducting material is coupled to a normal metal with external flux, and to a ferromagnet respectively. Finally, we provide a discussion of our results and concluding remarks.

Theory
In this section, we write down the quasiclassical equation of motion for ǧ in the diffusive limit and its belonging boundary conditions. The task at hand is then to solve this numerically in 2D and 3D, and we demonstrate how this can be accomplished using a finite element method in the next section. g is an 8 × 8 matrix satisfying =ǧ 1 2 with the following structure: holds in equilibrium, so that in this scenario one only needs to determine ĝ R in order to completely specify ĝ. The structure of the retarded Green function looks as follows: R where = ε g g( ) and = ε f f ( ) denote the 2 × 2 normal and anomalous Green function matrices in spin space, respectively. The … ∼ operation means complex conjugation and reversal of the energy argument ε → (− ε ).
The Usadel equation reads: 3 where D is the diffusion coefficient, ρ = ρ ρˆdiag( , ) 3 Scientific RepoRts | 6:22765 | DOI: 10.1038/srep22765 where  h describes the magnitude and direction of the magnetic exchange field while σ = σ σ σ  ̲ ̲ ̲ ̲ ( , , ) x y z is the vector of Pauli matrices. In the presence of gauge fields, such as a U(1) magnetic vector potential A describing an external magnetic field one has to replace the gradient operator with its covariant equivalent: 3 where q is the charge of the fermion field. A similar substitution is also made if one wishes to include an SU(2) gauge field  →  that describes antisymmetric spin-orbit coupling of Rashba or Dresselhaus type. In this work, we will use the standard Kupriyanov-Lukichev 13 boundary conditions as a realistic description of the interface regions. While originally derived for the tunneling regime, these boundary conditions have been shown to give good results also for moderately to highly transparent interfaces 33 , which are considered herein. For an interface separating a material 1 on the left side from a material 2 on the right side, they read: j j j j 1 2 Here, ζ j = R B /R j describes the ratio between the interface resistance and the bulk resistance of region j while L j is the length of region j. Here,  n is the unit vector normal to the interface pointing from region 1 to 2. At interfaces to air, no current is allowed to flow and the boundary condition is where  n again represents the unit vector normal to the air interface. Equations (3), (6), and (7) define a system of coupled differential equations with belonging boundary conditions and the task is to find the solution ǧ . For concreteness, we restrict our attention to an equilibrium scenario where only the retarded Green function matrix ĝ R must be found. Even with this restriction, the equations are capable of describing a variety of different mesoscopic systems. The equation system for ĝ R is identical to the one for ǧ , as can be verified by direct insertion of Eq. (1) in the place of ĝ, by replacing all …̌ matrices with their . equivalents. Before proceeding to a description of the finite element method we have used to solve this equation set in 2D and 3D, it is useful to introduce a Ricatti parametrization 34 of ≡ĝ g R . This parametrization, first applied in ref. 35 in the context of the Usadel equation, simplifies the numerical implementation of the equations by exploiting the symmetries and normalization of ǧ . One introduces two matrices in spin-space, γ and γ , which define ĝ as follows: This parametrization satisfies both the proper symmetry relations between the elements of ĝ as well as the normalization condition =ĝ 1 2 . Equations (3), (6), and (7) comprise a set of second-order coupled partial nonlinear differential equations which, when solved, determine the Green function ǧ of the system. Various physical quantities of interest may then be computed, such as the charge current density J Q and the density of states (DOS), given as: Another physical quantity that may be computed is the pair correlation function, Ψ , indicating the degree to which superconducting correlations exist in the system. It is given as: refers to the element in column i and row j of the Keldysh Green function matrix. A general analytical solution of equations (3), (6), and (7) is impossible. Some progress can be made by linearizing the equations, as is often done when considering a superconducting proximity effect. However, this approximation limits the validity of the obtained results and may cause the loss of novel physical phenomena that are only captured when the full equations are used. To do so, one must use a numerical approach. So far, only a handful of works have managed to solve the 2D Usadel equation numerically. This has been done in the full proximity effect regime for a superconductor/normal metal/superconductor junction in refs 26,27. To the best of our knowledge, no work has ever reported a solution of the Usadel equations in 3D.

Implementation of the finite element method.
We here present a way to solve the quasiclassical equations in 2D and 3D using a finite element method. Its detailed description follows below. After its presentation, we show its application to 2D and 3D model systems by solving the equations without any approximations.
Inserting 8 into equation 3 results in the following: Scientific RepoRts | 6:22765 | DOI: 10.1038/srep22765 where δ models the effect of inelastic quasiparticle scattering (the so-called Dynes parameter 36 ). We set δ /Δ = 10 −3 in this paper where Δ is the bulk superconducting gap. In Eq. (12), it is also possible to include self-energies corresponding to spin-flip and spin-orbit scattering on impurities which act pair-breaking on superconducting correlations. This typically amounts to a reduction of the magnitude of the superconducting proximity effect and we omit these terms in the present work. We also note that the effect of Rashba and Dresselhaus spin-orbit interactions were derived in Ricatti-parametrized form very recently 37,38 . As γ and γ  are 2 × 2 matrices, thus containing 4 elements each, it is clear that the solution of equation 12 involves solving a system of 8 coupled NLDEs. For brevity, we introduce the notation where γ ij and γ  ij are elements of γ and γ  respectively. Equation 12 may then be written as where α is an element of equation 13 and F (α) is a function that performs the matrix multiplications of equation 12 and extracts the appropriate element. Similarly, the boundary conditions become in the Riccati parametrization: where the negative sign should be used for a boundary where region j is to the right of region i, and the positive sign for a boundary where region j is to the left of region i. A similar expression is found for ⋅ ∇γ n i by applying the … ∼ operation to equation 15. These are Neumann boundary conditions of the type where B (α) works in a similar manner as F (α) . By multiplying equation 14 by a test function η  r ( ) and integrating over the domain Ω in which the equations are defined, one gets what is called the weak formulation of the NLDEs (not to be confused with the weak proximity effect approximation): where the divergence theorem has been used and ∂ Ω is the boundary of Ω. The unit vector ν  is an outward pointing surface normal, and is either parallel or antiparallel with the normal vector  n as defined in the Kupriyanov-Lukichev boundary conditions. It may thus be expressed as ν = ν ⋅     n n ( ) . It is assumed that the domain Ω can be discretized into a mesh of N el elements, i.e., N el subdomains Ω n , so that equation 17 becomes So far, no approximations have been made, and provided it is continuous in Ω n , the exact solution of equation 18 exists in the infinite space of polynomials P(Ω n ). To progress further, we will use the Galerkin method, a common finite element formulation technique treated in most books on the subject, e.g 39 . The method consists of restricting the space in which solutions are sought, from P(Ω n ) to a finite dimensional space of polynomials P N (Ω n ) consisting of all polynomials of degree N or lower. Normally, N is equal to 1 or 2.
On each element there are defined N n nodes, containing the degrees of freedom of the system -in this case the solution of the Usadel equation at the location of the node -and it is possible to define N n polynomials, φ  r ( ) j , that interpolate between them. These interpolation functions span the space of P N (Ω n ) and are used as a basis for the approximate solution of equation 18: We now consider the boundary term. Having meshed the domain Ω, it is obvious that some of the element domains Ω n intersect with the boundary ∂ Ω. In fact, the boundary is the union of all these intersections. It follows that the nodes associated with these intersections also lie on the boundary, and so there are defined interpolation functions also here. With the dimensionality of ∂ Ω being one less than Ω, the surface interpolation functions φ j S , which are zero everywhere but on the boundary, are found by evaluating the element interpolation functions at the surface, i.e., φ = φ  r ( ) j S j S where  r S is a surface coordinate. With the approximation given in 19, equation 18 is in general not satisfied, so that for every element the right hand side becomes equal to a residual, α R j , however due to the nonlinearities introduced by F (α) and B (α) this needs to be done iteratively by Newton-Raphson iterations: with J ij the Jacobian matrix in the 8 dimensional parameter space, given as Finally, 22 needs to be assembled into a global system of equations by summing over all elements, taking element connectivity into account. This involves restructuring and expanding the element matrices into a global system matrix: where  is an 8M × 8M matrix, and M is the number of nodes in the system. The integrals over the element domains are performed by changing coordinates to a reference element, and integrating numerically by means of a Gauss quadrature. This puts restrictions on how distorted a mesh can be, as the Jacobian for the coordinate transformation has to exist. In general, a structured mesh where the deviation from the geometry of the reference element is small will often give higher accuracy and reduce the computation time as the sparsity of the assembled matrices is increased.

Results
Application: 2D and 3D superconductor/ferromagnet junctions. The main advantage of the finite element method over the finite difference method, a method commonly used to solve partial differential equations numerically, is that it is formulated entirely without specifying element type, interpolation functions, spatial dimension or the geometry. This gives it the flexibility to solve PDEs on geometries which would be challenging to solve with the finite difference method. Here, we have used second order Lagrange polynomials as interpolation functions with quadrilateral (QUAD9) and hexagonal (HEX27) elements in 2D and 3D respectively. We illustrate this in the following. For the numerical implementation, we use the finite element library libMesh 40 and its integration with the PETSc library of numerical equation solvers 41,42 . In the following, the superconducting regions will be treated as reservoirs such that the bulk expression for the Green function =ǧ g BCS will be used. The superconductors thus effectively enter the problem as boundary conditions. 2D Josephson junction with external magnetic flux. It is well known that for a Josephson junction where an external flux is applied to the intermediate region, the supercurrent exhibits a Fraunhofer interference pattern. In refs 26,27 a 2D superconductor/normal/superconductor Josephson junction was studied in the presence of an external magnetic flux. The authors revealed that the Fraunhofer interference pattern would qualitatively change its dependence on the external flux depending on the width of the junction W relative to its length L. When W ≫ L a conventional Fraunhofer pattern was found, when W ≪ L the supercurrent was monotonically decaying. Moreover, it was shown that the Fraunhofer interference pattern was accompanied by a regular array of proximity-induced vortices in the transversal direction of the normal metal region. The vortices are not present in the narrow width limit. Experimental verification of the appearance of proximity-induced vortices was recently reported in ref. 43 which considers Josephson junctions generated by a network of superconducting nanocrystals.
Here, we explore how the vortices disappear from the system as the width is reduced and the system transitions to a vortex-less state. We also determine how a change in the phase difference between the superconducting leads affect the vortices, and will show that this does not always correspond to a shift of the vortex array along the transverse direction. To illustrate the ease with which the finite element method handles non-trivial geometries, we consider a Josephson junction with a bottleneck in the normal metal region, as shown in Fig. 1. We assume that the currents in the system are small, so that the magnetic field remains unaffected. As will be shown, it turns out to be possible to tune the geometry of the array along which the superconducting vortices align, swapping from a vertical necklace to a horizontal row of vortices and vice versa. Moreover, we demonstrate that changing superconducting phase difference, tunable e.g. via a current-bias, causes vortices to merge. This offers an interesting route to exerting external control over topological excitations in superconducting hybrid structures. It is seen that with no bottleneck, and with W ≫ L, a linear array of vortices along the y-axis is found. This is shown in Fig. 2a and is in agreement with refs 26,27. The number of vortices is simply equal to the number of flux quanta in the system. Furthermore, the fact that the vortices align themselves in an array implies that they repel, a feature they share with the Abrikosov vortices found in type II superconductors. Decreasing the width, pushes the vortices closer together which is energetically less favorable. With no phase difference between the superconducting leads, the phase correlation function is symmetric about both the x-and y-axis. This means that when the system becomes too narrow to sustain four vortices, two vortices must simultaneously translate vertically out of the system in a way which maintains this symmetry, as seen in Fig. 2b. The two remaining vortices are seen to be forced closer together until they eventually meet at the origin, from which it may be inferred that for the given flux and geometry, the presence of two vortices is energetically favorable regardless of their separation. In particular, it is observed that within numerical precision, the vortices are found to completely overlap in Fig. 2c, resulting in a single vortex. The winding of the phase of the pair correlation function along a contour around this vortex is found to be 4π , implying a topological charge of 2, see inset of Fig. 2c. As the bottleneck is introduced, and the width further decreased, the vortices split symmetrically along the x-axis, as shown in Fig. 2d,e. This behavior may also be explained by the symmetry of the system, which restrains the positions of the two vortices to be symmetric about the origin, on either the x-axis or the y-axis. As the vortices evidently feel a stronger repulsion from the edges than from each other, they are pushed together. However once they meet at the origin, they are free to separate along the x-axis and thus reduce the energy in the system. By continuing to decrease the bottleneck width, W b , a point where even two vortices may not be sustained is eventually reached. The boundary conditions that constitute the superconducting leads enforce a constant pair correlation, and so it becomes increasingly difficult to maintain the curvature necessary for the vortices to exist as one approaches the superconductors. In other words, the vortices may not in a continuous fashion exit the system along the x-axis. Instead, the vortices are seen to return to the y-axis, and be expelled vertically.
While the vortices separate along the x-axis for decreasing bottleneck width, the length of the narrowing area is large enough to contain them, and so the system behaves as if the width is uniformly decreased. It has been verified that by reducing the horizontal extent of the bottleneck, it is possible to create a situation where the vortices are pushed to the wide regions of the junction, at which point they become virtually independent of the bottleneck width W b .
We also show the results for varying phase difference between the superconductors, for a geometry with W b = 0.6W, shown in Fig. 3. We find that not only are the positions of the vortices changed by varying the phase difference φ , but also the number of vortices is altered. Figure 3a-c show the absolute value of the pair correlation function. With no phase difference, two vortices are located symmetrically along the x-axis. As φ increases from 0, the two vortices coalesce at the origin. Further increase translates one of the vortices in the negative y-direction, until only a single vortex remains. We have confirmed that the spatial rearrangement of the vortex pattern and the merging of vortices also takes place even without the bottleneck geometry, i.e. for a rectangular N region. The dependence of the vortex positions on the phase difference may be explained by the magnetic field, which in the small current approximation, permeates the normal metal unhindered. The current generated by the phase difference is altered by the field which in turn influences which locations that are energetically favorable for the vortices. The current density for each of the cases considered are shown in Fig. 3d-f. Close to a vortex, where pair correlation is low, currents are induced by the magnetic field and circulate counter-clockwise. Due to the omnipresent magnetic field, screening currents circulating clockwise are generated which dominate when pair correlation is high. As the pair correlation function is weakened upon approaching the vortex core, one observes an abrupt change in the current density pattern at a certain distance from the vortex.
The phase of the pair correlation function is given as θ = Ψ Ψ I R arctan( / ). By integrating ∇ θ along a contour going around a point where the pair correlation function vanishes, a value of 2π is found. This can be seen directly from the phase plots in Fig. 3g-i, as any curve around a zero of the pair correlation function has to traverse two discontinuous jumps of value π . In other words, these points have a topological charge of one, showing that they are indeed vortices. With the approximation of weak currents we do not however find flux quantization, as this requires a self consistent calculation of the magnetic field.
3D ferromagnetic nanoisland. We also demonstrate how the finite element method is capable of dealing with fully three-dimensional structures with non-rectangular geometry by considering a superconductor/ferromagnet bilayer as depicted in Fig. 4. The ferromagnet is cylindrical with a radius of R = 2ξ and a height of L z = 0.4R, and is placed atop an assumed infinite superconductor. Such a geometry is inspired by ref. 44 which experimentally explores the appearance of magnetic field induced superconductivity in a lattice of ferromagnetic islands placed on top of a superconductor. While the experimental setup is far too sophisticated for their results to be recreated by the example considered herein, it does demonstrate the relevance of the model.
We use Kupriyanov-Lukichev boundary conditions with a resistance ratio of ζ = 3. We compute the density of states (DOS) for this structure with an exchange field h equal to 0.3Δ, 0.5Δ and 0.7Δ in the vertical direction, as shown in Fig. 5a. The results are identical with the one-dimensional solution to the S/F bilayer, displaying an enhanced DOS at the Fermi level (ε = 0) and a spin-split minigap structure [45][46][47][48][49][50][51] . The spatial distribution of the DOS is nearly constant for this particular parameter set choice in F as illustrated by Fig. 5b, thus proving the correctness of the method.
3D ferromagnet with superconducting islands. To illustrate the 3D capabilities of the method developed on a system which cannot be described by an effective 1D model, we consider two variations of a system where two superconducting islands are placed on a ferromagnet with dimensions L x × L y × L z = 10ξ × 7ξ × ξ , as shown in Fig. 6. To avoid self-consistency iterations, the islands are approximated by bulk BCS superconductors, and are included as Kupriyanov-Lukichev boundary conditions with a resistance ratio of ζ = 1.5. The dimensions of the islands are 2.5ξ × 2.5ξ . The motivation for these analyses is to study the current flow between the superconducting islands and the spatial modulation of the density of states due the proximity effect in the presence of a supercurrent. To this end, the islands are given a phase difference of φ = π 2 . The configurations considered are:   Configuration A may be realized experimentally by growing the ferromagnetic film on a substrate and placing superconducting electrodes on top of it. Configuration B may be created by placing the lower superconducting electrode on the substrate and subsequently grow a normal metal film on top of it (a similar type of geometry was considered in the context of Fraunhofer patterns in ref. 52). The upper superconducting electrode is then placed atop the film. The spatial distribution of the magnetization used in Configuration B can be generated by placing a strong ferromagnet on top of the normal film. This will magnetize the film across the thickness, creating a cross section approximately equal to the ferromagnet within which the magnetization is constant. The ability of the ferromagnet to induce magnetization within the normal film abates quickly as the distance from it increases, thus generating the distribution shown in Fig. 6c.  The results from both configurations are given in Fig. 7. For Configuration A it is seen that the current is largely confined to the region between the islands, passing from one to the other, shown in Fig. 7a,g. Due to the magnetization being uniform, the supercurrent travels along a path that minimizes the distance, and thus has no component pointing along the transversal direction (the y-axis of Fig. 7). However, as the current enters and exits the ferromagnet vertically, it is seen to arc into the thickness of the film, as seen in Fig. 7c,e. Configuration B has a magnetization of = ∆  h r ( ) 15 within a horizontal radius  r of one superconducting coherence length ξ from the center of the ferromagnet. This has a significant effect on the supercurrent. The supercurrent flows vertically into the system from the lower superconductor, as seen in Fig. 7d. However, rather than flowing directly to the upper superconductor, as would have been the case for a homogeneous magnetization, the current is seen to avoid the area of highest magnetization by following a semicircular path, shown in Fig. 7b,d,h. The exchange field has a detrimental effect on the superconducting correlations as it breaks up the Cooper pairs. For this reason it is natural that the path selected by the supercurrent eventually transitions from the shortest route, to a path where the central area is avoided as h increases. In this sense, the exchange field is seen to influence the supercurrent in a way which is analogous to the way a resistance influences a normal current. It is interesting that this transition has occurred already for h = 15Δ, which is to be considered a somewhat weak magnetization, and may provide means for customizing the supercurrent path.
In Fig. 8 we show the density of states (DOS) along two different lines on the surface of the ferromagnet in Configuration B, thus simulating the measurement of ∝ DOS by scanning tunneling microscopy. The DOS is seen to feature a strong spatial modulation. Along the x-axis, as seen in Fig. 8a, the probed line passes directly underneath the upper superconductor. Here the characteristic peaks associated with the superconducting DOS are observed at ε = ± Δ. Similar peaks are also created on the surface above the lower superconductor. Furthermore, a slight suppression of the DOS is found in the same regions, at the level of ε = ±  h r ( ), which is typical for superconductor/ferromagnet hybrid structures 49 . The second line is placed opposite the lower superconductor, in the y-direction as shown in Fig. 8b. Also here, the characteristic peaks and split gap is found. The proximity effect is seen to decay as one moves away from the position of the superconductor, so that the DOS approaches that of a normal metal, which is reasonable.

Conclusion
We have demonstrated how the full, spin dependent, Usadel equation may be solved by the finite element method. The method excels in solving differential equations for non-trivial geometries and may find use in solving a wide range of problems which have not been manageable with other methods. A natural development of the finite element method presented herein would be to incorporate the kinetic equations coming from the Keldysh part of the quasiclassical equations in non-equilibrium situations. The methodology may also be generalized to handle time dependent problems such as domain wall motion. Work is currently ongoing on these subjects which may find interesting applications in the field of superconducting spintronics 53,54 .