SuperConga: An open-source framework for mesoscopic superconductivity

Over the last few decades, there has been a tremendous improvement in the fabrication and controllable downscaling of samples and devices. The physical size of a sample can now be comparable to characteristic quantum length scales of the composing materials. For a superconducting material, which is the topic of this article, the relevant length scales are the superconducting coherence length, the length scale over which the superconducting order parameter typically varies, and the magnetic penetration depth, which determines the length scale of screening. In direct connection with the ability to make smaller samples and devices, new measurement and sensing tools have been developed. Current scanning probes can spatially resolve minute variations in magnetic fields using scanning SQUIDS^1–5 or make detailed spatial maps of low-lying quasiparticle states and coherences using scanning tunneling spectroscopy (STS).^6–15 By fabricating nano-scale cantilevers, the intrinsic state of a superconducting island can be monitored as changes in the cantilever-oscillation frequency.^16–18 Other techniques of detecting and following the superconducting state in nanoscale devices include Hall magnetometry,^19,20 scanning electron microscopy (SEM),²¹ and scanning Hall probe microscopy.^22–24 Scanning single electron transistors (SET) are local probes that can map out charging, or parity, effects that become prominent for ultra-small superconductors where the energy-level spacing and the superconducting pairing are of the same order of magnitude.^25–29 Experimental data gathered using the above-mentioned methods give rich and complex information about the strongly correlated material system under study. To fully understand these data, one needs to compare and contrast them to theoretical predictions and modeling.

Superconductivity in metals is explained by the Bardeen–Cooper–Schrieffer (BCS) theory.^30–33 Its extension and generalization to unconventional superconductivity and superfluids^34,35 makes the theory a corner stone of condensed matter physics. One method to model superconducting devices and to explore fundamental problems of superconductivity with the BCS theory as a starting point is the quasiclassical theory of superconductivity. It is the extension of Landau's Fermi-liquid theory to include superconducting phenomena and was pioneered by Eilenberger,³⁶ Larkin and Ovchinnikov,³⁷ and Eliashberg.³⁸ The quasiclassical approximation relies on a separation between relevant energy, temperature, and length scales in the normal-metal and superconducting states. The normal-metal state is characterized by the Fermi energy, $EF$⁠, Fermi temperature, $TF$⁠, and inverse Fermi wave number, $1/kF$⁠, while the superconducting state is characterized by the superconducting order parameter, Δ, transition temperature, $Tc$⁠, and coherence length, ξ₀. Quasiclassical theory is then a controlled expansion in small parameters like $Δ/EF$⁠, $Tc/TF$⁠, or $1/(ξ0kF)$⁠, which are usually of order $10−2–10−3$ in metals. As a Green's function based theory,³⁹ it is very powerful and include in its most general form material effects such as impurity scattering, Fermi-liquid effects, electron–phonon interaction, and non-equilibrium situations for instance imposed by external fields or potentials.^40,41 At the same time, it is versatile enough that it can be adapted to realistic device sizes and geometries.^42–45 In its simplest form, quasiclassical theory is equivalent to the Andreev approximation of the Bogoliubov–deGennes equation for the ballistic case.^32,46,47

When the mean free path, $ℓ$⁠, due to elastic impurity scattering is much smaller than the superconducting coherence length, ξ₀, it is possible to derive diffusion equations from the more general Eilenberger–Larkin–Ovchinnikov equations. These diffusion equations, first derived by Usadel,^48,49 are more easy to handle and is widely used to describe conventional diffusive s-wave samples and devices.^50–52 Solution methods in the diffusive case include Nazarov's circuit theory⁵³ and finite-element methods.⁵⁴ Here, we focus instead on the more general equations in the ballistic regime.

A complication when solving the equations of quasiclassical theory is that one often needs to resort to numerical methods. Developing the necessary codes is technically demanding and time-consuming. Today, there is no open-source code general enough to describe ballistic devices or superconducting grains of mesoscopic size with the quasiclassical theory. This lack of an open-source code forces all researchers and students in the community to re-implement the theory for each individual problem.

To remedy this, in this paper, we present an open source framework for studies of two-dimensional superconducting grains. The application program interface (API) is written in C++ and CUDA and takes advantage of the speedup made possible by running the code on Graphics Processing Units (GPUs). At the same time, the Python-based frontend is sufficiently easy to use that any user interested only in the physics and results never has to dwell on the technical details of the implementation. The first version of the API that we present here can be used to study conventional and unconventional singlet superconducting grains of general geometry in two dimensions with an applied external magnetic field. The framework is sufficiently modular that in the future it can be extended to include more aspects of quasiclassical theory.

This paper is organized in the following: In Sec. II, we give an introduction to quasiclassical theory with the aim to provide a self-contained account of the theory, which the framework is based on. In this section, we also show some simple model calculations, serving as a background to the examples worked out in later sections. In Sec. III, we give an overview of the main algorithm of SuperConga and give sufficient background to understand the parameters that have to be set while running the simulations. In Sec. IV, we show by a simple example how to run the framework using the Python frontend. Section V contains a few more involved examples that benchmark and highlight the capabilities of SuperConga. We study the rather complicated physics and energetics of vortices in a superconducting grain. In particular, we find $Bc1$ where it becomes energetically favorable to have one vortex in a finite size disk. By comparison to formulas from literature, we are also able to extract the vortex-core energy. At high fields, we demonstrate that many different vortex configurations can be stabilized with almost the same free energy. For lower fields, such configurations might be more separated in energy but can be studied as metastable states. In another example, we find the vortex configurations and the field-dependent local density of states (LDOS) in a grain studied experimentally by Timmermans et al.⁵⁵ Finally, we show that in a d-wave superconducting annulus, the free-energy parabolas are conventional³³ at high temperatures, while at low temperatures zero-energy Andreev bound states at the edges modifies the energetics through their paramagnetic response to the external field. Interestingly, spontaneous circular surface currents^56–58 can coexist at high external magnetic fields with superconducting phase windings $n2π$ (n integer) around the annulus and integer flux $nΦ0$ threading the inner circular hole. Section VI contains a summary of the paper.  Appendixes A–F contain additional information on the SuperConga external dependencies and parameter configuration file, details about self-consistency convergence acceleration, the technique to solve for the vector potential and the equations of motion, as well as tables summarizing our choice of units and order-parameter basis functions.

The superconducting state is characterized by an order parameter,

that is zero for temperatures above the superconducting transition (⁠$T>Tc$⁠) and non-zero below it. All superconductors break U(1)-gauge symmetry, and hence, the minimal order parameter is described by an amplitude $|Δ(R)|$ and a phase $χ(R)$⁠, both, in general, dependent on the spatial coordinates $R$⁠. The phase is directly coupled to the superfluid momentum, $ps(R)$⁠, and the electromagnetic gauge field, $A(R)$⁠, via the gauge invariant expression,

where $e=−|e|$ is the charge of the electron, $ℏ$ is Planck's constant, and c is the speed of light. In addition to breaking the U(1) symmetry, some superconducting compounds also break symmetries of the crystal lattice. This possibility is encoded in the complex-valued basis function $ηΓ(pF)$⁠, which depends on the Fermi momentum, $pF$⁠. The index Γ denotes the irreducible representation of the crystallographic point group that the basis function belongs to.⁵⁹ In the general case, we also need to account for the spin-degree of freedom, and the order-parameter field should, therefore, be written as a spin-matrix $Δαβ(pF,R)$⁠. However, we will focus on spin-singlet superconductivity, and the form in Eq. (1) is adequate for this.

The singlet order parameter satisfies the following gap equation:

where $V(pF,pF′)$ is the effective superconducting pairing interaction. We use the approximation that the interaction can be separated into symmetry channels of the crystal point group,

where $VΓ$ is the pairing strength in the symmetry channel Γ. In this case, the dependencies on position and momentum direction of the order parameter separate as in Eq. (1). The angle bracket $⟨⋯⟩pF$ denotes a Fermi-surface average, which in 2D is a line integral around the Fermi surface according to^60,61

where the total normal-state density of states per spin at the Fermi level is

We consider the extreme layered superconductor with conduction in stacked two-dimensional planes. In this case, $NF=kF/2πℏvFd$⁠, where $kF,vF$ are the in-plane Fermi-wave number and Fermi velocity, respectively, and d is the inter plane distance.⁶⁰ The vector potential is determined through Ampère's law,

where the charge current density is defined as

where the factor 2 accounts for spin degeneracy. The Fermi velocity at $pF$ on the Fermi surface is $vF(pF)=∇pε(p)|p=pF$⁠, where $ε(p)$ is the dispersion.

The two functions $g(pF,R;εn)$ and $f(pF,R;εn)$ appearing in Eqs. (3) and (8) are components of the quasiclassical Green's function, or propagator,

It is a 2 × 2 matrix in electron–hole (Nambu) space, as indicated by the “ $̂$ ”-symbol. The propagator is determined from the quasiclassical counter part of the Gorkov–Dyson equation, the transport-like Eilenberger equation,

where $τ̂3$ is the third Pauli matrix in Nambu space. In addition to Eq. (10), the propagator obeys the normalization condition,

where $Î$ is the identity matrix.

In our case of weak-coupling superconductivity in the clean limit, the self-energy $Δ̂(pF,R)$ entering in Eq. (10) simply consists of the order parameter,

There is some redundancy in the parametrization of the propagator, and the following symmetry^35,62

between tilded (⁠$x̃$⁠) and un-tilded (x) quantities holds. The propagator may be evaluated at imaginary frequencies $ε→iεn=iπT(2n+1)$ giving the Matsubara propagator or at real frequencies giving either the Retarded $(ε→ε+iδ)$⁠, or the Advanced propagator $(ε→ε−iδ)$⁠. Here, δ is a small and positive energy broadening. To compute the order parameter, the current, and the back-coupling through the gauge field, it is enough to work with the Matsubara propagator. For the local density of states, defined as

we need to also compute the retarded Green's function with the order parameter and gauge field as input.

The superconducting state has a characteristic energy scale given by the transition temperature $Tc$⁠, namely, $2πkBTc$⁠, and a natural length scale, the coherence length $ξ0=ℏvF/2πkBTc$⁠. We will use the natural units $ℏ=kB=1$⁠. The normal state density of states at the Fermi level $NF$ is in units $[Energy×unit cell×spin]−1$⁠. In addition to the coherence length, the length scale for screening of electromagnetic fields enters the theory through the penetration depth λ₀, defined as

We use the Ginzburg–Lanadu parameter κ, defined as $κ=λ0/ξ0$⁠, to give the value of the penetration depth in relation to the coherence length. To define the magnetic units we start from the superconducting magnetic flux quantum $Φ0=hc/2|e|$⁠. The magnetic field,

is then given in units of $B0=Φ0/πξ02$ and the vector potential in units of $A0=Φ0/πξ0$⁠. Finally, current densities will be given in units of $j0=2πkBTc|e|vFNF$⁠. A summary of these units and natural scales are given in Table II in  Appendix F.

The quasiclassical theory is a conserving theory in the sense that the equations for $ĝ, Δ̂$⁠, and $A$ above may be derived as stationarity conditions from the Luttinger–Ward free-energy functional.^35,63 For the theory to be conserving, the stationarity equations are to be solved to self-consistency. At convergence, we solve Eq. (10) with $Δ(pF,R)$ and $A(R)$⁠, recalculate Eqs. (3) and (7), and get back the same $Δ(pF,R)$ and $A(R)$ within some given tolerance. Current conservation is then fulfilled to the same tolerance. At self-consistency, we have found an extremum of the quasiclassical version of the Luttinger–Ward functional.^35,56,64 This functional always give a measure of the free energy with respect to the normal state,³⁵ $ΩS[T]−ΩN[T]$⁠, and goes to zero as $Δ→0$⁠. To simplify notation, we denote this free energy difference $Ω[T]$⁠, and it has the form,

where $Bind(R)$ is the induced magnetic field, $|Δ(R)|2=⟨|Δ(pF,R)|2⟩pF$⁠, and

The auxiliary propagator $ĝλ$ is obtained by solving the Eilenberger equation with scaled self-energy field $λΔ̂$ so the integral over the dummy variable λ in Eq. (18) may be evaluated. In the framework SuperConga, we also apply the original Eilenberger-form of the free-energy functional,³⁶ in which case the term with the λ-integration is simplified to

The advantage of the Eilenberger-form is that it can be evaluated without the λ-integration at little computational cost. However, it must be used with some caution as it has been shown to not be generally valid and final results for the free energy needs to be confirmed by Eq. (17).⁶⁰ 

The above form of the free energy is consistent with the following regularized gap equation:

for superconductivity with the symmetry given by the representation Γ. The dependencies on the pairing interaction strength $VΓ$ and the Matsubara sum cutoff $εc$ have been replaced by the measurable transition temperature $TcΓ$⁠. After this regularization, the Matsubara sum is not cutoff dependent. The connection between Eqs. (3) and (20) can be seen as replacing the coupling strength by

By moving the second term on the right-hand side of Eq. (20) over to the left-hand side, and using the fact that the basis functions are orthonormal, $⟨ηi*(pF)ηj(pF)⟩pF=δij$⁠, we can divide down the resulting prefactor of $ΔΓ(R)$⁠, obtaining the gap equation in a fix point form,

The currently available symmetries and the corresponding basis functions $ηΓ(pF)$ are listed in Table III in  Appendix F.

For inhomogeneous superconducting states we need to solve the first order differential equation for the propagator $ĝ(pF,R;ε)$ in Eq. (10), subject to the normalization condition in Eq. (11). The gradient term $ivF·∇R ĝ(pF,R;ε)$ couples the momentum direction and spatial coordinates. Quasiparticle trajectories are then defined by this directional derivative and we may introduce a trajectory coordinate s as

where $v̂F=vF/vF$ is a unit vector. Starting with an initial value, for instance at a boundary point $Rmin$⁠, the propagator should be found by integrating along the trajectory until another boundary point $Rmax$ is met. At that point, a boundary condition must be supplemented to the theory. In general, boundaries present strong perturbations and their treatment lay outside the strict validity of the quasiclassical approximation. Nevertheless, one can model scattering off surfaces or at semitransparent interfaces via effective boundary conditions derived and carried over from a more microscopic theory.^{35,47,62,65–70} In certain cases, boundary and interface effects may arise, which lay outside the standard quasiclassical approximation.^71–73 In SuperConga a simple specular boundary condition is used at all boundaries, see illustration in Fig. 1.

In practice, the most convenient and stable formulation is based on a parametrization^62,74–76 of the Green's function in terms of coherence functions $γ(pF,R;ε)$ and $γ̃(pF,R;ε)$⁠. They are related by the symmetry in Eq. (13) and quantify electron–hole coherence along the trajectory in the superconducting state. The quasiclassical Green's function can be written in terms of them as

where we for brevity dropped the explicit dependencies on $pF, R$⁠, and ε. The normalization of the Green's function in Eq. (11) is now automatically satisfied. The coherence functions satisfy a set of Riccati equations,

The equation for $γ(pF,R;ε)$ is stable to integrate in the direction parallel to $vF$⁠, while the equation for $γ̃(pF,R;ε)$ is stable to be integrated in the opposite direction. After introducing the trajectory coordinate, we may let $vF·∇γ→vF∂sγ$⁠.

In this section, we give a brief overview of the main algorithm of SuperConga and describe a few key features of the framework. The aim of this section is to explain why certain parameters related to the implementation have to be set and how the parameter choices may influence the simulations.

We solve the Riccati equations, Eqs. (25) and (26), numerically (see  Appendix C) along trajectories through the grain. An example of a trajectory is shown in Fig. 1. Given an initial boundary value, $γ(Rmin)$⁠, we solve Eq. (25) until we reach $Rmax$⁠. Similarly, given an initial boundary value, $γ̃(Rmax)$⁠, we solve Eq. (26) until we reach $Rmin$⁠.

Clearly, the boundaries are special points in that we need a starting value for $γ(Rmin)$ and we obtain after stepping along the trajectory an end value $γ(sN)$ at another boundary, where N is the number of steps taken along the trajectory s. All boundary values are related to each other through the boundary condition. The complication is then that the end values become, after the boundary condition has been applied, a starting value for a different trajectory. This can be solved by having an initial guess for all boundary starting values and self-consistently computing the boundary values by repeatedly stepping along all trajectories and applying the boundary condition. Since the order parameter $Δ(R)$ is also solved self-consistently by iteration, the update of the boundary values can be done in parallel. At the end of the computation, when the gap equation is satisfied, the boundary condition is also self-consistently satisfied. A caveat in SuperConga is that while the order parameter and vector potential are saved to file, the coherence functions are not, due to being too disk-space intensive. When restarting the calculation using a nominally self-consistent order parameter, the self-consistency of the boundary conditions must be reached again. This is solved via a “burn-in” period, where the order parameter and vector potential are kept constant for a few iterations when restarting a calculation, until the coherence functions and boundary conditions are also self-consistent.

In SuperConga, the order parameter and the vector potential have been discretized in two-dimensional (2D) real space: $Δ(R)=Δ(x,y)→Δ(i,j)$ and $A(R)=A(x,y)→A(i,j)$⁠, where i and j are integers. The grid is a simple square grid with spacing h, i.e., x = ih and y = jh, but we still take into account the exact location and shape of boundaries, such that there is no aliasing or similar artifacts. Figure 2 visualizes the implementation and discretization. The Fermi surface is parameterized by an azimuthal angle $φ∈[0,2π)$ and discretized according to $φk=kΔφ$ (k integer). The Fermi surface averages are computed through numeric integration, via one of the following methods: Boole's, Simpson's 3/8, Simpson's 1/3, or the trapezoidal rule, depending on if the number of momenta is divisible by 4, 3, 2, otherwise, respectively. For a particular point $φk$ on the Fermi surface, we obtain a Fermi velocity $vF(φk)$⁠, which determines a direction for a set of trajectories in real space. Each trajectory coordinate has been discretized as outlined above. Considering all parallel trajectories for a certain point on the Fermi surface, the discrete grid in real space then forms a simple square grid with the same lattice spacing h, but rotated relative to the reference grid where the order parameter and vector potential are defined. Bilinear interpolation is then performed on the order parameter and the vector potential from the reference frame to the rotated frame. After solving the Riccati equations, the Riccati amplitudes are known in the rotated frame. The $φk$ terms in the Fermi-surface averages for the order parameter, and the charge-current density are computed in the rotated frame by constructing the relevant matrix element of the Green's function, summing over energies, and then rotating back to the reference frame. In this way, the order parameter $Δ(i,j)$⁠, observables such as the current $j(i,j)$⁠, and also the vector potential $A(i,j)$ can be updated through Eqs. (3)–(7).

The boundary condition is problematic as the azimuthal angle of the specularly reflected momentum at a given surface might not exist in the set ${φk}$⁠. Furthermore, due to the rotation of the grid, the trajectories hit the boundary at different points for each Fermi velocity $vF(φk)$⁠. In order to determine the initial conditions for the coherence functions, we once again perform bilinear interpolation. Here, the interpolation is between spatial coordinates on the boundary $R∈∂D$ and azimuthal angles.

An overview of the algorithm of SuperConga is shown in Algorithm 1. The user defines, e.g., the system domain, physical parameters, the discretization h in real space, the discretization $Δφ$ in momentum space, and an energy cutoff for the frequency summations [see, e.g., Eq. (3)]. The frequency summations are not over the usual Matsubara frequencies, but instead a more efficient summation due to Ozaki⁷⁷ has been implemented. In addition, a convergence criterion for the self-consistency must be set. The global error is defined as

where i here denotes the iteration number in the self-consistency loop, and $p∈{1,2,∞}$ is set by the user. The loop continues until the error $εG$ is smaller than the convergence criterion, or the given maximum number of iterations have been performed. The gap equation can be solved by direct fixed-point iteration, but in SuperConga a few methods for convergence acceleration have been implemented. In simplified terms, the accelerators may save a number of computed order parameters, ${Δi,Δi−1,Δi−2,…}$⁠, and make a more educated guess for $Δi+1$⁠. When observables, the vector potential and free energy are also computed, the convergence criterion is the same as for the order parameter. All the details of the stepping along trajectories, interpolation, as well as the treatment of the boundaries are automatically taken care of by the API. The user must be aware of, however, that the discretization in real and momentum space can be a source of error, and it is necessary to check convergence by running the same problem with finer grids. Also, depending on the problem being solved, the path to self-consistency may not be direct. For the simplest examples in this paper, this is not a problem. However, for more advanced problems, it is advisable to experiment with different starting guesses for the order parameter and also carefully choose a suitable convergence acceleration strategy.

The main algorithm and functionality of SuperConga is implemented as a fairly modular framework, in the form of a backend written in C++ and CUDA.^78,79 A user can in principle choose which modules to use and in what manner, to set up highly customized simulations, or even extend the existing functionality by exploiting the modularity. However, this requires writing C++ code that imports or adds to the desired features. While this is certainly possible and a valid approach, a main goal of SuperConga is to also enable users to perform advanced simulations without having to write any code, by instead only specifying simulation parameters, but without losing crucial functionality. To make this possible, SuperConga comes with a user friendly frontend, and a set of backend interfaces. The frontend manages and validates the user input, and then sends it to the appropriate backend interface. The backend interface takes care of importing and calling the necessary backend modules with the settings provided by the user. See Fig. 3 for a schematic overview. SuperConga also comes with a set of ready-to-run examples (see Sec. IV for a demonstration), a detailed user manual full of tutorials and guides,⁸⁰ a set of helpful tools, a Singularity container which simplifies setup and improves portability, and an extensive testing framework to test the validity of the code. The latter ranges from smaller unit tests of basic functionality to larger end-to-end physics test of, e.g., current conservation and well-known analytic results. In the following, we give further insight into the backend and frontend, and list the most important modules.

The backend modules are implemented as a set of classes and functions via header files. The modularity is typically built up by an inheritance structure, from virtual base classes to derived classes. The virtual base classes describe the general signature or pattern for a module but does not necessarily contain any functionality. This is instead taken care of by the derived classes, which inherits the signature from the corresponding virtual class, and implements the actual functionality (or offloads it to external libraries, see App. A for a description of dependencies). For example, BoundaryCondition.h contains a virtual base class defining the general properties of a boundary condition. The derived class in BoundaryConditionSpecular.h inherits from BoundaryCondition.h and implements a specular boundary condition. Similarly, there are virtual base classes for Fermi surfaces, geometric shapes, and even Riccati solvers, with specific implementations of, e.g., a circular Fermi surface, disk- and polygon-geometries, and a Riccati solver for finite grains. The idea is that the modularity should allow for straightforward extension of anything from boundary conditions to new observables. Below is a list of the most important modules:

• Boundary conditions: Implementation of the boundary conditions to handle interface scattering of quasiparticles, the latter which follow the Eilenberger–Riccati trajectories. Currently only implements specular boundary conditions for superconductor-vacuum interfaces.

• Fermi surface: A class that takes care of the momentum-averaging over the Fermi surface. Can either treat circular Fermi-surfaces, or parametrize a more general Fermi-surface shape based on a tight-binding hopping model.

• Geometry: Classes and functionality for creating mesoscopic grains from basic shape primitives (i.e., disks and polygons).

• Observables: Implementation of various observables and quantities, specifically the superconducting gap, charge current density, magnetic moment, LDOS, free energy, vector potential, and the magnetic induction. Also, takes care of computing the residual for the self-consistency iteration.

• Accelerator: Implementation of self-consistency accelerators. The user can currently choose from basic Picard iterations, Polyak's⁸¹ “small heavy sphere,” a variant of the Barzilai–Borwein (BB) method,⁸² and CongAcc which is an adaptive method developed for SuperConga. See  Appendix E for details and comparisons.

• Order parameter: Class for setting up the order parameter, possibly with multiple components like d + is and chiral d-wave. Note that the implementation is currently limited to spin-singlet with a single band and spin degenerate Fermi-surfaces.

• Riccati solver: Implements a solver for the Riccati equations, currently based on the mid-point method for confined (finite) geometries. To give an idea, this class could more generally be extended with, e.g., a higher order Runge–Kutta solver, or other types of systems like bulk or semi-infinite superconductors.

• Self-consistency solver: A collection of classes that take care of adding the contribution of the Riccati trajectories for each degree of freedom (energy, Fermi-surface angle, coordinate) to self-consistently solve the order parameter from the superconducting gap equation, together with the vector potential from the corresponding Maxwell equation.

See the API documentation for more modules and information.⁸⁰ 

SuperConga currently comes with two backend interfaces, that expose the functionality of the above modules to the user (via the frontend). The first backend interface is simulation.cu, which acts as the “main” program, essentially implementing Algorithm 1 to self-consistently solve the gap equation and Maxwell's equation. More specifically, simulation.cu starts by generating a superconducting grain and Fermi surface based on the user-provided geometry and model. The superconducting order parameter and its components are initialized in this grain, according to the specified basis functions, transition temperatures, and start guess. The program then continues to set up the self-consistency accelerator, Riccati solver, boundary conditions, observables, and vector potential. An OpenGL instance is initialized if live visualization is enabled. Finally, the self-consistency loop is run until reaching either the convergence criterion or maximum number of iterations, upon which data are saved to disk. The second backend interface is postprocess.cu, which computes the LDOS based on results from a previous simulation (e.g., computed with simulation.cu).

Note that both backend interfaces take as command-line input the location of a configuration file, which should contain all parameters necessary to completely specify the simulation. The simulation will not run if the configuration file is not found, or if the contents are invalid. The purpose of the frontend is to make it easier to set up all simulation parameters correctly, with validation and proper help messages to guide the user if anything goes wrong.

The SuperConga frontend consists of a main run-file, superconga.py, which draws on an assortment of functionality implemented in a modular Python library, the latter located in the folder frontend/. The goal of the frontend is to facilitate setting up and analyzing simulations. For example, the frontend consists of Python modules for parsing command-line arguments, validating them and providing help for the user, reading and writing configuration files used by the backend, calling the correct backend interfaces, reading and visualizing simulation data, converting data between different formats, and modifying data by, e.g., adding noise. The functionality is divided into a set of subcommands, used via python superconga.py <subcommand>. A list of available frontend subcommands is obtained by calling python superconga.py –help, and each subcommand also has its own help message. A short summary of each subcommand will now be given, and Sec. IV demonstrates how to use them.

The first subcommand is setup, which is used to setup and generate the SuperConga build system via CMake.

The second subcommand is compile, which is used after the setup to compile the SuperConga framework. These two steps have to be performed before any simulations can be run.

The third subcommand is simulate, which calls the binary generated from simulation.cu to run self-consistent simulations. It takes as input the location of a configuration file, defining all simulation parameters. Optionally, each parameter can be set or overridden via the command-line interface (CLI). The parameters will then be validated for errors, attempting to provide the user with helpful messages if failing. If valid, the binary is called with the user-defined settings.

The fourth subcommand is plot-simulation, which takes the location of data from a self-consistent simulation, and plots computed quantities (see Fig. 4 for example). The user can control what to plot, in what order, and which settings to use, for example, colormaps and fonts. If the user has specified that data should be saved when running simulate, then plot-simulation is automatically called, using the default settings, and a PDF is saved together with the data.

The fifth subcommand is postprocess, which calls the binary generated from postprocess.cu to calculate the LDOS from self-consistently converged data. The runner takes as command-line input the location of the data, and parameters specifying which resolution to compute the LDOS with. The latter parameters can either be set directly via the CLI, or specified in an external configuration file.

The sixth subcommand is plot-postprocess, which is a spectroscopy tool. It takes as input the path to a directory containing LDOS data computed with postprocess. If found, the data are read and visualized in an interactive plot (see Fig. 5 for example).

The seventh subcommand is plot-convergence, which plots the residuals of computed quantities vs iteration number. This tool can be run on convergence data from either simulate or postprocess subcommands.

The eighth subcommand is convert, which is used to convert data back and forth between the file formats h5 and csv.

Having given an overview of the SuperConga framework and its frontend, we will now demonstrate how to use it to set up simulations and do data analysis.

The SuperConga framework ships with a few ready-to-run examples that illustrate how to set up and run simulations. These are also explained in the online documentation, together with detailed tutorials and research guides.⁸⁰ One of the examples will now be shown, namely, the simulation of a single Abrikosov vortex in a conventional superconducting disk. The demonstration can be split into the following generic steps:

• Pre-processing: Setting up and defining the simulation,

• Processing: Running the self-consistency simulations, optionally with live visualization,

• Post-processing: Data analysis, visualization and spectroscopy.

The simulation configuration for the vortex example can be found in examples/swave_disc_vortex/in the SuperConga main folder.⁸³  Appendix B lists the full contents of this file, with a detailed explanation of the format. In short, each configuration file contains all necessary parameters needed to completely specify a simulation. The parameters are grouped into different sections: Listing 1 shows the group of “physics” parameters for the Abrikosov example, specifying the temperature $T=0.5Tc$⁠, the external flux $Φext=1.5Φ0$⁠, the penetration depth $λ=5ξ0$⁠, and an s-wave order parameter with a phase winding of $−2π$ as a start guess.

Other parameter groups deal with, for example, the geometry, numerical accuracy, convergence criteria, convergence acceleration, and so on. The parameters can also be set or changed directly via command line, effectively overriding the settings in the configuration file, as shown below. Assuming that the framework has been properly installed and compiled (see the online documentation and installation guide⁸³) the configuration file can be used to start the simulation by running the following from the root directory of SuperConga,

which should run a few tens of iterations until the convergence criterion in simulation_config.json is fulfilled, printing simulation status for each iteration to terminal. Here, the flag -C is used to specify the relative path to the configuration file. Adding the arguments “-T 0.1 -B 10.0” will change the temperature to $0.1Tc$⁠, and the external flux to $10Φ0$⁠. To start with an order parameter from a file, corresponding to, e.g., a previous simulation or an arbitrary computed guess, the -L argument can be used. Using the flag --help gives further information about which parameters are available and their usage.

We note that live visualization can be turned on or off with the flags --visualize and --no-visualize, respectively. The main purpose of the live-visualization is to get an overview of the simulation progress. Thus, it is geared toward speed rather than producing production-ready plots. For the latter, we instead recommend to visualize the fully converged results, via the data files generated by the program. The user can either use their favorite plotting tool, or use SuperConga as described in the following. Having followed the example in Listing 2, the Abrikosov vortex results can be visualized by running the following from the root directory.

This plots all the computed quantities, as shown in Fig. 4. Various properties of the plot can be controlled via command-line arguments (as described by the help message), like fonts, colormaps, which quantities are plotted, and if saving the plot directly to file rather than drawing in a window. The latter makes it easy to automatically generate plots from a large number of simulations.

The local density of states (LDOS) can be calculated from converged data as a post-processing step, by calling,

where the contents of postprocess_config.json is shown in Listing 5.

This file specifies, e.g., which energies and resolution to use in the LDOS calculation, and where the order-parameter data are located. The empty save path indicates that the LDOS data will be saved to the same directory as the order parameter data.

The LDOS can be particularly tricky to visualize due to the high dimensionality. To aid with this, the framework provides a spectroscopy tool in the form of an interactive LDOS visualizer. Assuming that the vortex LDOS has been calculated with the above post-processing example in Listing 4, vortex spectroscopy can be done by running the following from the root directory:

Figure 5 shows a screenshot of the interactive LDOS plotter for this example, illustrating that the user can click in the window to choose which energy to plot the LDOS at (as a 2D-heatmap vs coordinates), and which point in space to plot the LDOS vs energy at (as a 1D curve).

To extract the temperature dependence, or the dependence on any other parameter for that matter, the most straightforward approach is to write a script which calls the program multiple times, but with unique values of the parameter. A few such examples are included in SuperConga, and the following one illustrates how to simulate an Abrikosov vortex for different values of Ginzburg–Landau coefficient κ, throughout the whole type-II range $κ∈[1,∞)$⁠:

The results, seen in Fig. 6, can subsequently be visualized via

So far, the demonstration has focused on a superconducting disk. SuperConga allows the user to specify more general composite geometries by successively adding and removing disks and polygons. Both regular and free-form polygons are implemented, enabling quite intricate systems that, e.g., has holes, are multiply connected, or even disconnected but coupled via induction. For example, Fig. 1 was created in this way by setting the “geometry” parameter section of the configuration file. Listing 9 shows a SQUID-like geometry, created by adding a central region shaped like an octagon, two arms in the form of a rectangle, and then removing a disk from the central region:

In Listing 9, the entry “geometry” consists of a list which specifies in which order (from top to bottom) the shapes should be added or removed, the latter indicated by “add”: true or “add”: false, respectively. All coordinates and lengths are given in units of ξ₀, and the user can choose whichever origin they prefer, as it does not influence the physics. The disk is uniquely defined via its center coordinates and radius, a regular polygon by center coordinates, number of edges, rotation (in units of $2π$⁠), and side-length. For free-range polygons, a list of x- and y-coordinates for the vertices are specified (in counterclockwise order). Currently, the free-range polygons are limited to convex hulls. Note that an arbitrary number of each shape in principle can be listed, which was previously used to, e.g., study mesoscopic roughness by generating boundaries with random removal of polygons.⁸⁴ Note, however, that for extremely large number of added/removed shapes, the performance might eventually be reduced due to having to deal with a large set of complicated trajectories and boundary conditions.

It is straightforward to extend the above example of a single Abrikosov vortex, to vortex lattices, by increasing the external field. As an example, Fig. 7 plots the zero-energy LDOS in nanoscale grains of various geometries exposed to high external fields. Vortex cores show up as dot-like dark spots, illustrating vortex lattices with various mesoscopic properties. In contrast to bulk lattices, the number of vortices and their separation in a finite system does not only depend on the external field strength, but also on, e.g., the boundary conditions, the system size relative to the penetration depth, and the shape of the system. While bulk lattices are usually triangular, a mesoscopic vortex lattice can mimic the shape and symmetry of the superconducting grain. In addition, rather than entering one by one as the external field is increased, several vortices might suddenly enter at a time, sometimes as concentric shells (hence referred to as the vortex-shell effect). These effects have been observed, e.g., in grains shaped like disks,^3,19,21 squares,^{21,55,88–91} triangles,^21,92–95 and pentagons.⁹⁶ Such scenarios are great examples of when a sufficiently sophisticated simulation tool like SuperConga is crucial, since one needs to capture the combined effects of a confined geometry and a finite penetration depth with back-coupling of the induction. In particular, one needs to use a powerful convergence accelerator that can properly traverse the phase space. For example, we note that with the proper choice of start guess and accelerator settings, it is in principle possible to study giant vortices,^97–102 as well as vortex-antivortex pairs.^{90,91,103–106} Careful treatment and use of the adaptive convergence accelerators show, however, that these structures are generally unstable and decay to more traditional vortex structures with a lower free energy.

In this section, we present several examples of studies that are quite straightforward to perform with the framework SuperConga. In a few special cases, it is possible to find analytic results, which enable us to check the overall correctness of the framework. First, we study a 2D superconducting annulus. In the case of an s-wave superconductor, we may compare with analytics, while in the d-wave superconducting case, we demonstrate the corrections due to suppression of the order parameter at the edges and the formation of spontaneous currents at low temperature.^{57,58,84,107–111} Second, we study a superconducting disk also subjected to an external magnetic field. Here, we return to the well-studied problem of vortex lattice formation, which has a long history going back to the original work of Abrikosov,¹¹² and in particular, to that of Pearl considering superconducting disks.¹¹³ For the cases of zero or one vortex at the center of the disk, we compare with analytic results. Then, we continue with higher fields, where we have to resort to a numerical treatment. Finally, we simulate irregular superconducting islands in an external magnetic field. The study is inspired by experimental work done by Timmermans et al.⁵⁵ on spectroscopy on superconducting nanostructures assembled by small squares of Mo₇₉Ge₂₁.

We begin by studying homogeneous superflow in a superconductor. Superflow implies that there is a finite superfluid momentum caused by phase gradients and/or a vector potential, as defined in Eq. (2). The fact that the phase gradient enters in the superfluid momentum becomes clear when doing a gauge transformation. In quasiclassical theory this is done as follows. We start with the order-parameter $Δ(pF,R)=Δ(R)η(pF) exp [iχ(R)]$⁠, which, in general, is complex valued. The order-parameter matrix can then be decomposed as

with the transformation matrix

and $Δ0(R)$ a purely real amplitude (i.e., $χ0(R)≡0$ while the basis function $η(pF)$ may still be complex valued). Applying the same transformation to the Eilenberger Eq. (10), suppressing the arguments of $ĝ$ and $Δ̂$⁠, we see that if $ĝ$ solves

with $Δ̂=ÛΔ̂0Û†$⁠, then $ĝ0=Û†ĝU$ is a solution to

The superfluid momentum $ps$ in Eq. (2) is now naturally formed, and we get the Eilenberger equation in the form

with the purely real order-parameter matrix.

Let us look at the linear response to a small and spatially homogeneous superfluid momentum, $ps≪Δ/vF$⁠. To do this, we write the propagator as a perturbation expansion $ĝ=ĝ0+δĝ+O[ps2]$⁠. Using the normalization condition on the propagator $ĝ2=−π2Î$⁠, we get

Since $ĝ0∝iεnτ̂3−Δ̂0(pF)$⁠, we get via the Eilenberger equation the solution

defining $Λn(pF)=|Δ(pF)|2+εn2$⁠. From the form of $δĝ$⁠, we can directly see that the linear correction to the anomalous Green's function component, $f(pF;εn)$⁠, is odd in frequency and will not give a correction in leading order in $ps$ to the order parameter $Δ(pF)$⁠. The diagonal part of $δĝ$ allows us to calculate the current to linear order in $ps$ as

with the angle-dependent Yoshida function,

Equation (34) is usually written in the well-known form,

defining

the angle-dependent superfluid-density, or superfluid stiffness tensor.¹¹⁴ 

As a final step, we write the change in the free energy in linear response as

where

is the bulk superconducting free energy (see, e.g., Ref. 35), and $A$ is the area of the superconducting sample considered. The change in free energy in Eq. (37) due to superflow consists of a kinetic contribution,⁵⁸ 

and a magnetic contribution

where $Bext$ is the applied external magnetic field.

With the above introduction to superflow in a superconductor, we are ready to study a superconducting annulus with outer (inner) radius $R>$ (⁠$R<$⁠), as illustrated in Fig. 8(c). The area of the annulus is $A=π(R>2−R<2)$⁠. Let us first consider the following text-book vector potential¹¹⁵ used to illustrate the Aharonov–Bohm effect,¹¹⁶ 

where $ϕ̂$ is the unit vector along the azimuthal direction and $Φext$ is an externally applied magnetic flux threading only the hole. We call this particular gauge the solenoid gauge. It can be seen as a mathematical idealization to guarantee zero external magnetic flux through the superconducting annulus, while the vector potential remains finite everywhere. One possible physical realization of this gauge would be to let the superconducting annulus encircle a strong magnet.¹¹⁷ 

Using the gauge transformation Eq. (28), we can impose a phase winding in the order-parameter field going around the annulus. One of the constraints on an order parameter is that it is single valued. This condition leads to that χ in Eq. (28) must be chosen so that the total winding going around the void is a multiple n of $2π$⁠, and we can write $χ=nφ$⁠, where n can take any integer value. If we now look at Eq. (2) we get, reinstating $ℏ$ and c,

We see that $ps$ can be made to vanish for certain values of the external flux, $Φext=−nΦ0$⁠, i.e., the superfluid momenta vanishes throughout the superconducting annulus for every external field that matches integer multiples of the flux quantum $Φ0$⁠. At those matching fields, $δΩkin=0$⁠. As a consequence of the gauge properties, we see that if we have a self-consistent solution for the pair $ĝn,Δ̂n$ at a given magnetic flux $Φext$ threading the hole, we can, via an appropriate gauge transformation $Û$⁠, construct a solution at any other external flux that differ from $Φext$ by an integer number of flux quanta. This setup nicely demonstrates the fundamental concept of flux quantization.^118–122 Note that the negative sign of the phase winding n reflects our choice of unit, $e=−|e|$⁠, where the particle probability current flows with the gradient of the phase, while the associated charge current is in the opposite direction, see Eq. (35).

The simplest s-wave superconductor has an isotropic order parameter with $Δ(pF)=Δ$⁠, where Δ is a complex scalar quantity and the basis function $η(pF)=1$⁠. For such an order parameter, scattering off sample surfaces leaves the superconducting state unaffected. As a consequence, linear-response theory will be sufficient to describe the annulus up to quite large superfluid momenta. We can, thus, calculate the kinetic contribution to the free energy Eq. (39), using the homogeneous solution of Eq. (20), to get

after integrating over the superconducting area and setting

where $v̂F(pF)$ and $p̂s$ are unit vectors. For this particular gauge, we let the penetration depth $λ→∞$ and neglect screening effects. This leads to that the free-energy contribution $δΩmagn$ is zero. The free energy comes in a family of parabolas as function of the external flux. Each parabola is centered around $−nΦ0$ where n is the integer flux quantization, or phase winding, enforced on the order parameter. This is clearly shown in Fig. 8(a), where we also show results obtained by direct calculation using SuperConga, for phase windings n = 0 and n = −1. The bottom of each parabola corresponds to the free energy $Ωbulk(T)$ at zero external flux for the considered temperature.

We finally mention that at temperatures close to the transition temperature, $Tc, δΩkin$ can made to exceed the energy gain $Ωbulk$ when the external flux is close to $(n+12)Φ0$⁠. This is a manifestation of the Little–Parks effect,^33,123–125 where one can induce a field-modulation of the critical temperature of a superconducting film enclosing a hole.

Next, we change the superconductor to be a d-wave superconductor. The order-parameter field for the d-wave depends on position on the Fermi surface, $pF$⁠, via the basis function, either the $dx2−y2$ or d_xy as listed in Table III. This momentum dependence leads to that scattering off a surface may be pair breaking and one must solve for the spatial dependence of the order parameter using Eq. (20).

As for the s-wave superconductor, the free energy in an external magnetic field come in a family of parabolas $∼(n+Φext/Φ0)2$⁠. However, as shown in Fig. 8(b), neither $Ωbulk(T)$ nor $δΩkin(T,n,Φext)$⁠, are captured correctly by a calculation with a homogeneous order parameter magnitude and superflow treated in linear response. The order-parameter magnitude of a d_xy-wave superconducting annulus is shown in Fig. 8(d). As seen there are regions near the surfaces where the magnitude is severely reduced, which leads to correction to $Ωbulk(T)$⁠.

The pair-breaking at interfaces that can be found in unconventional superconductors is due to low-energy Andreev states that are localized in a region of a few coherence lengths width at the surface.¹²⁶ These states modify low-temperature properties of the superconductor. A disk-shaped unconventional superconductor was studied by Suzuki and Asano finding a paramagnetic instability at low temperatures as a response to a small magnetic field.¹²⁷ We find similar paramagnetic behavior for the annular superconductor at low temperatures and relate the effects to spontaneous generation of phase gradients,^57,107 in a state known as a phase crystal.^58,61,128

In Fig. 9, the free-energy is shown for a set of temperatures for the two cases where the initial guess for the order-parameter field is homogeneous, and where we lock a phase winding of $2π$ around the annulus. Below $T∼0.2Tc$⁠, the existence of zero-energy Andreev states start to drive the superconductor to generate a finite $ps$ by spontaneous phase gradients.^57,107 This finite superfluid momentum gives rise to a Doppler shift, $vF·ps$⁠, that lifts the zero-energy states away from the Fermi level and, thus, lowers the free energy of the system. The free-energy minimum below $T*≈0.18Tc$ is shifted to occur at finite magnetic flux. At the lowest temperature, we show here, $T=0.1Tc$⁠, we can find several metastable states at small fluxes. The orange rhombus in Fig. 9 has a purely real order parameter and the zero-energy states are not shifted away from the Fermi level. This state is only stable at $Φext≡0$⁠, and any small seed of a phase gradient in the order parameter, or a minute magnetic perturbation, will generate a configuration with spontaneous currents. The two low-energy states we find are combinations of two possible states as showed in Fig. 9.

In contrast to the case with the solenoid gauge above, a more common situation experimentally is when the external magnetic field $Bext=Bextẑ$ is applied in a homogeneous fashion and penetrates also the superconductor. We incorporate this case through a symmetrical gauge $Aext(R)=12Bext(−y x̂+x ŷ)$ subject to the condition $∇·A=0$⁠.

Let us study the London–Maxwell equation for the vector potential $A$ in the symmetrical gauge,^{32,33,113,129,130}

for the same annulus geometry as in Sec. V A 1. We allow for a phase winding $2πn$ of the order parameter going around the annulus and include screening of the externally applied field. The temperature-dependent penetration depth $λ(T)$ entering Eq. (45) is defined as

As long as the inner radius is much larger than the coherence length and the external magnetic field is moderate in strength, the superfluid momentum will be small, and we can disregard any reduction of the order-parameter amplitude for an s-wave annulus. In this case, linear-response theory remains valid and gives a platform to verify more involved calculations produced by SuperConga. Below, we introduce the dimensionless applied magnetic flux,