We present SuperConga, an open-source framework for simulating equilibrium properties of unconventional and ballistic singlet superconductors, confined to two-dimensional (2D) mesoscopic grains in a perpendicular external magnetic field, at arbitrary low temperatures. It aims at being both fast and easy to use, enabling research without access to a computer cluster, and visualization in real-time with OpenGL. The core is written in C++ and CUDA, exploiting the embarrassingly parallel nature of the quasiclassical theory of superconductivity by utilizing the parallel computational power of modern graphics processing units. The framework self-consistently computes both the superconducting order-parameter and the induced vector potential and finds the current density, free energy, induced flux density, local density of states (LDOS), and the magnetic moment. A user-friendly Python frontend is provided, enabling simulation parameters to be defined via intuitive configuration files, or via the command-line interface, without requiring a deep understanding of implementation details. For example, complicated geometries can be created with relative ease. The framework ships with simple tools for analyzing and visualizing the results, including an interactive plotter for spectroscopy. An overview of the theory is presented, as well as examples showcasing the framework's capabilities and ease of use. The framework is free to download from https://gitlab.com/superconga/superconga, which also links to the extensive user manual, containing even more examples, tutorials, and guides. To demonstrate and benchmark SuperConga, we study the magnetostatics, thermodynamics, and spectroscopy of various phenomena. In particular, we study flux quantization in solenoids, vortex physics, surface Andreev bound-states, and a “phase crystal.” We compare our numeric results with analytics and present experimental observables, e.g., the magnetic moment and LDOS, measurable with, for example, scanning probes, STM, and magnetometry.

## I. INTRODUCTION

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),^{21} 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,^{36} Larkin and Ovchinnikov,^{37} and Eliashberg.^{38} 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, *ξ*_{0}. Quasiclassical theory is then a controlled expansion in small parameters like $\Delta /EF$, $Tc/TF$, or $1/(\xi 0kF)$, which are usually of order $10\u22122\u201310\u22123$ in metals. As a Green's function based theory,^{39} 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, $\u2113$, due to elastic impurity scattering is much smaller than the superconducting coherence length, *ξ*_{0}, 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^{53} and finite-element methods.^{54} 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.*^{55} Finally, we show that in a *d*-wave superconducting annulus, the free-energy parabolas are conventional^{33} 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\pi $ (*n* integer) around the annulus and integer flux $n\Phi 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.

## II. QUASICLASSICAL THEORY OF SUPERCONDUCTIVITY

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 $|\Delta (R)|$ and a phase $\chi (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=\u2212|e|$ is the charge of the electron, $\u210f$ 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 $\eta \Gamma (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.^{59} 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 $\Delta \alpha \beta (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\u2032)$ 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\Gamma $ 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 $\u27e8\cdots \u27e9pF$ 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\pi \u210fvFd$, where $kF,vF$ are the in-plane Fermi-wave number and Fermi velocity, respectively, and *d* is the inter plane distance.^{60} 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)=\u2207p\epsilon (p)|p=pF$, where $\epsilon (p)$ is the dispersion.

The two functions $g(pF,R;\epsilon n)$ and $f(pF,R;\epsilon 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 “ $\u2009\u0302$ ”-symbol. The propagator is determined from the quasiclassical counter part of the Gorkov–Dyson equation, the transport-like Eilenberger equation,

where $\tau \u03023$ is the third Pauli matrix in Nambu space. In addition to Eq. (10), the propagator obeys the normalization condition,

where $I\u0302$ is the identity matrix.

In our case of weak-coupling superconductivity in the clean limit, the self-energy $\Delta \u0302(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\u0303$) and un-tilded (*x*) quantities holds. The propagator may be evaluated at imaginary frequencies $\epsilon \u2192i\epsilon n=i\pi T(2n+1)$ giving the Matsubara propagator or at real frequencies giving either the Retarded $(\epsilon \u2192\epsilon +i\delta )$, or the Advanced propagator $(\epsilon \u2192\epsilon \u2212i\delta )$. 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\pi kBTc$, and a natural length scale, the coherence length $\xi 0=\u210fvF/2\pi kBTc$. We will use the natural units $\u210f=kB=1$. The normal state density of states at the Fermi level $NF$ is in units $[Energy\xd7unit\u2009\u2009cell\xd7spin]\u22121$. In addition to the coherence length, the length scale for screening of electromagnetic fields enters the theory through the penetration depth *λ*_{0}, defined as

We use the Ginzburg–Lanadu parameter *κ*, defined as $\kappa =\lambda 0/\xi 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 $\Phi 0=hc/2|e|$. The magnetic field,

is then given in units of $B0=\Phi 0/\pi \xi 02$ and the vector potential in units of $A0=\Phi 0/\pi \xi 0$. Finally, current densities will be given in units of $j0=2\pi 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 $g\u0302,\u2009\Delta \u0302$, 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 $\Delta (pF,R)$ and $A(R)$, recalculate Eqs. (3) and (7), and get back the same $\Delta (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,^{35} $\Omega S[T]\u2212\Omega N[T]$, and goes to zero as $\Delta \u21920$. To simplify notation, we denote this free energy difference $\Omega [T]$, and it has the form,

where $Bind(R)$ is the induced magnetic field, $|\Delta (R)|2=\u27e8|\Delta (pF,R)|2\u27e9pF$, and

The auxiliary propagator $g\u0302\lambda $ is obtained by solving the Eilenberger equation with scaled self-energy field $\lambda \Delta \u0302$ 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,^{36} 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).^{60}

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\Gamma $ and the Matsubara sum cutoff $\epsilon c$ have been replaced by the measurable transition temperature $Tc\Gamma $. 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, $\u27e8\eta i*(pF)\eta j(pF)\u27e9pF=\delta ij$, we can divide down the resulting prefactor of $\Delta \Gamma (R)$, obtaining the gap equation in a fix point form,

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

### A. Riccati formalism for inhomogeneous states

For inhomogeneous superconducting states we need to solve the first order differential equation for the propagator $g\u0302(pF,R;\epsilon )$ in Eq. (10), subject to the normalization condition in Eq. (11). The gradient term $ivF\xb7\u2207R\u2009g\u0302(pF,R;\epsilon )$ 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\u0302F=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 $\gamma (pF,R;\epsilon )$ and $\gamma \u0303(pF,R;\epsilon )$. 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,\u2009R$, 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 $\gamma (pF,R;\epsilon )$ is stable to integrate in the direction parallel to $vF$, while the equation for $\gamma \u0303(pF,R;\epsilon )$ is stable to be integrated in the opposite direction. After introducing the trajectory coordinate, we may let $vF\xb7\u2207\gamma \u2192vF\u2202s\gamma $.

## III. OVERVIEW OF THE IMPLEMENTATION

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.

### A. The main algorithm of SuperConga

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, $\gamma (Rmin)$, we solve Eq. (25) until we reach $Rmax$. Similarly, given an initial boundary value, $\gamma \u0303(Rmax)$, we solve Eq. (26) until we reach $Rmin$.

Clearly, the boundaries are special points in that we need a starting value for $\gamma (Rmin)$ and we obtain after stepping along the trajectory an end value $\gamma (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 $\Delta (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: $\Delta (R)=\Delta (x,y)\u2192\Delta (i,j)$ and $A(R)=A(x,y)\u2192A(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 $\phi \u2208[0,2\pi )$ and discretized according to $\phi k=k\Delta \phi $ (*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 $\phi k$ on the Fermi surface, we obtain a Fermi velocity $vF(\phi 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 $\phi 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 $\Delta (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 ${\phi k}$. Furthermore, due to the rotation of the grid, the trajectories hit the boundary at different points for each Fermi velocity $vF(\phi 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\u2208\u2202D$ and azimuthal angles.

Define domain $D$ and input parameters |

Initialize Δ and $Aind$ |

Initialize $\gamma \u2202D$ and $\gamma \u0303\u2202D$ ▹ The boundary |

while not converged do |

for all $pF$, ε _{n}do |

Rotate Δ and $A$ so that $y\u0302\u2225vF$ |

Compute γ and $\gamma \u0303$ along the y-axis ▹ See App. C |

Update $\gamma \u2202D$ and $\gamma \u0303\u2202D$ ▹ Write incoming $pF$ |

Update Δ, $j$, and Ω ▹ Eqs. (22), (8), (17) |

Compute $Aind$ ▹ See App. D |

Check convergence ▹ Eq. (27) |

$\Delta ,Aind\u2190$ Accelerator ▹ See App. E |

Update $\gamma \u2202D$ and $\gamma \u0303\u2202D$ ▹ Boundary condition |

Define domain $D$ and input parameters |

Initialize Δ and $Aind$ |

Initialize $\gamma \u2202D$ and $\gamma \u0303\u2202D$ ▹ The boundary |

while not converged do |

for all $pF$, ε _{n}do |

Rotate Δ and $A$ so that $y\u0302\u2225vF$ |

Compute γ and $\gamma \u0303$ along the y-axis ▹ See App. C |

Update $\gamma \u2202D$ and $\gamma \u0303\u2202D$ ▹ Write incoming $pF$ |

Update Δ, $j$, and Ω ▹ Eqs. (22), (8), (17) |

Compute $Aind$ ▹ See App. D |

Check convergence ▹ Eq. (27) |

$\Delta ,Aind\u2190$ Accelerator ▹ See App. E |

Update $\gamma \u2202D$ and $\gamma \u0303\u2202D$ ▹ Boundary condition |

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 $\Delta \phi $ 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^{77} 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\u2208{1,2,\u221e}$ is set by the user. The loop continues until the error $\epsilon 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, ${\Delta i,\Delta i\u22121,\Delta i\u22122,\u2026}$, and make a more educated guess for $\Delta 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.

### B. Framework overview

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,^{80} 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.

#### 1. Backend

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^{81}“small heavy sphere,” a variant of the Barzilai–Borwein (BB) method,^{82}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.^{80}

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.

#### 2. Frontend

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.

## IV. DEMONSTRATION OF SUPERCONGA

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.^{80} 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.^{83} 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 $\Phi ext=1.5\Phi 0$, the penetration depth $\lambda =5\xi 0$, and an *s*-wave order parameter with a phase winding of $\u22122\pi $ as a start guess.

“physics”: { |

“temperature”: 0.5, |

“external_flux_quanta”: 1.5, |

“penetration_depth”: 5.0, |

“crystal_axes_rotation”: 0.0, |

“gauge”: “symmetric”, |

“charge_sign”: −1, |

“order_parameter”: { |

“s”: { |

“critical_temperature”: 1.0, |

“initial_phase_shift”: 0.0, |

“initial_noise_stddev”: 0.0, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: -1.0 |

} |

] |

} |

} |

} |

“physics”: { |

“temperature”: 0.5, |

“external_flux_quanta”: 1.5, |

“penetration_depth”: 5.0, |

“crystal_axes_rotation”: 0.0, |

“gauge”: “symmetric”, |

“charge_sign”: −1, |

“order_parameter”: { |

“s”: { |

“critical_temperature”: 1.0, |

“initial_phase_shift”: 0.0, |

“initial_noise_stddev”: 0.0, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: -1.0 |

} |

] |

} |

} |

} |

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^{83}) 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\Phi 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.

python superconga.py simulate -C examples/swave_disc_vortex/simulation_config.json |

python superconga.py simulate -C examples/swave_disc_vortex/simulation_config.json |

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.

python superconga.py plot-simulation -L data/examples/swave_disc_vortex |

python superconga.py plot-simulation -L data/examples/swave_disc_vortex |

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,

python superconga.py postprocess -C examples/ |

swave_disc_vortex/postprocess_config.json |

python superconga.py postprocess -C examples/ |

swave_disc_vortex/postprocess_config.json |

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

{ |

“spectroscopy”: { |

“energy_max”: 0.5, |

“energy_min”: 0.0, |

“energy_broadening”: 0.008, |

“num_energies”: 128 |

}, |

“numerics”: { |

“convergence_criterion”: 1e-4, |

“norm”: “l2”, |

“num_energies_per_block”: 32, |

“num_fermi_momenta”: 720, |

“num_iterations_burnin”: −1, |

“num_iterations_max”: 1000, |

“num_iterations_min”: 0 |

}, |

“misc”: { |

“data_format”: “h5”, |

“load_path”: “data/examples/swave_disc_vortex”, |

“save_path”: “”, |

“verbose”: true |

} |

} |

{ |

“spectroscopy”: { |

“energy_max”: 0.5, |

“energy_min”: 0.0, |

“energy_broadening”: 0.008, |

“num_energies”: 128 |

}, |

“numerics”: { |

“convergence_criterion”: 1e-4, |

“norm”: “l2”, |

“num_energies_per_block”: 32, |

“num_fermi_momenta”: 720, |

“num_iterations_burnin”: −1, |

“num_iterations_max”: 1000, |

“num_iterations_min”: 0 |

}, |

“misc”: { |

“data_format”: “h5”, |

“load_path”: “data/examples/swave_disc_vortex”, |

“save_path”: “”, |

“verbose”: true |

} |

} |

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:

python superconga.py plot-postprocess -L data/ |

examples/swave_disc_vortex |

python superconga.py plot-postprocess -L data/ |

examples/swave_disc_vortex |

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 $\kappa \u2208[1,\u221e)$:

./examples/parameter_sweeps/ |

swave_disc_vortex_kappa_sweep.sh |

./examples/parameter_sweeps/ |

swave_disc_vortex_kappa_sweep.sh |

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

python examples/parameter_sweeps/ |

plot_vortex_kappa_sweep.py |

python examples/parameter_sweeps/ |

plot_vortex_kappa_sweep.py |

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:

“geometry”: [ |

{ |

“polygon”: { |

“add”: true, |

“vertices_x”: [0.0, 0.0, 100.0, 100.0], |

“vertices_y”: [5.0, −5.0, −5.0, 5.0] |

} |

}, |

{ |

“regular_polygon”: { |

“add”: true, |

“center_x”: 50.0, |

“center_y”: 0.0, |

“num_edges”: 8, |

“rotation”: 0, |

“side_length”: 15.0 |

} |

}, |

{ |

“disc”: { |

“add”: false, |

“center_x”: 50.0, |

“center_y”: 0.0, |

“radius”: 8.0 |

} |

} |

] |

“geometry”: [ |

{ |

“polygon”: { |

“add”: true, |

“vertices_x”: [0.0, 0.0, 100.0, 100.0], |

“vertices_y”: [5.0, −5.0, −5.0, 5.0] |

} |

}, |

{ |

“regular_polygon”: { |

“add”: true, |

“center_x”: 50.0, |

“center_y”: 0.0, |

“num_edges”: 8, |

“rotation”: 0, |

“side_length”: 15.0 |

} |

}, |

{ |

“disc”: { |

“add”: false, |

“center_x”: 50.0, |

“center_y”: 0.0, |

“radius”: 8.0 |

} |

} |

] |

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 *ξ*_{0}, 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\pi $), 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.^{84} 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.^{96} 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.

## V. EXAMPLES AND RESULTS

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,^{112} and in particular, to that of Pearl considering superconducting disks.^{113} 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.*^{55} on spectroscopy on superconducting nanostructures assembled by small squares of Mo_{79}Ge_{21}.

### A. Superconducting annulus in a magnetic field

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 $\Delta (pF,R)=\Delta (R)\eta (pF)\u2009exp\u2009[i\chi (R)]$, which, in general, is complex valued. The order-parameter matrix can then be decomposed as

with the transformation matrix

and $\Delta 0(R)$ a purely real amplitude (i.e., $\chi 0(R)\u22610$ while the basis function $\eta (pF)$ may still be complex valued). Applying the same transformation to the Eilenberger Eq. (10), suppressing the arguments of $g\u0302$ and $\Delta \u0302$, we see that if $g\u0302$ solves

with $\Delta \u0302=U\u0302\Delta \u03020U\u0302\u2020$, then $g\u03020=U\u0302\u2020g\u0302U$ 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\u226a\Delta /vF$. To do this, we write the propagator as a perturbation expansion $g\u0302=g\u03020+\delta g\u0302+O[ps2]$. Using the normalization condition on the propagator $g\u03022=\u2212\pi 2I\u0302$, we get

Since $g\u03020\u221di\epsilon n\tau \u03023\u2212\Delta \u03020(pF)$, we get via the Eilenberger equation the solution

defining $\Lambda n(pF)=|\Delta (pF)|2+\epsilon n2$. From the form of $\delta g\u0302$, we can directly see that the linear correction to the anomalous Green's function component, $f(pF;\epsilon n)$, is odd in frequency and will not give a correction in leading order in $ps$ to the order parameter $\Delta (pF)$. The diagonal part of $\delta g\u0302$ 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.^{114}

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,^{58}

and a magnetic contribution

where $Bext$ is the applied external magnetic field.

#### 1. The solenoid gauge

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=\pi (R>2\u2212R<2)$. Let us first consider the following text-book vector potential^{115} used to illustrate the Aharonov–Bohm effect,^{116}

where $\varphi \u0302$ is the unit vector along the azimuthal direction and $\Phi 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.^{117}

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\pi $, and we can write $\chi =n\phi $, where *n* can take any integer value. If we now look at Eq. (2) we get, reinstating $\u210f$ and *c*,

We see that $ps$ can be made to vanish for certain values of the external flux, $\Phi ext=\u2212n\Phi 0$, i.e., the superfluid momenta vanishes throughout the superconducting annulus for every external field that matches integer multiples of the flux quantum $\Phi 0$. At those matching fields, $\delta \Omega kin=0$. As a consequence of the gauge properties, we see that if we have a self-consistent solution for the pair $g\u0302n,\Delta \u0302n$ at a given magnetic flux $\Phi ext$ threading the hole, we can, via an appropriate gauge transformation $U\u0302$, construct a solution at any other external flux that differ from $\Phi 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=\u2212|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).

##### a. An s-wave superconducting annulus in the solenoid gauge

The simplest *s*-wave superconductor has an isotropic order parameter with $\Delta (pF)=\Delta $, where Δ is a complex scalar quantity and the basis function $\eta (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\u0302F(pF)$ and $p\u0302s$ are unit vectors. For this particular gauge, we let the penetration depth $\lambda \u2192\u221e$ and neglect screening effects. This leads to that the free-energy contribution $\delta \Omega magn$ is zero. The free energy comes in a family of parabolas as function of the external flux. Each parabola is centered around $\u2212n\Phi 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 $\Omega bulk(T)$ at zero external flux for the considered temperature.

We finally mention that at temperatures close to the transition temperature, $Tc,\u2009\delta \Omega kin$ can made to exceed the energy gain $\Omega bulk$ when the external flux is close to $(n+12)\Phi 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.

##### b. A d-wave superconducting annulus in the solenoid gauge

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\u2212y2$ 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 $\u223c(n+\Phi ext/\Phi 0)2$. However, as shown in Fig. 8(b), neither $\Omega bulk(T)$ nor $\delta \Omega kin(T,n,\Phi 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 $\Omega 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.^{126} 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.^{127} 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\pi $ around the annulus. Below $T\u223c0.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\xb7ps$, 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*\u22480.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 $\Phi ext\u22610$, 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.

#### 2. Solving London–Maxwell equation in a symmetric gauge

In contrast to the case with the solenoid gauge above, a more common situation experimentally is when the external magnetic field $Bext=Bextz\u0302$ is applied in a homogeneous fashion and penetrates also the superconductor. We incorporate this case through a symmetrical gauge $Aext(R)=12Bext(\u2212y\u2009x\u0302+x\u2009y\u0302)$ subject to the condition $\u2207\xb7A=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\pi n$ of the order parameter going around the annulus and include screening of the externally applied field. The temperature-dependent penetration depth $\lambda (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,

##### a. Analytic solution of the London–Maxwell equation

Equation (45) is a diffusion equation and can be solved analytically. We give the solution here for completeness. In a cylindrically symmetric case, the solution to Eq. (45) is given by the modified Bessel functions of first and second kind, $I\mu (x)$ and $K\mu (x)$, respectively. Writing $A=A(R)\phi \u0302$, the general solution, in the external magnetic field $Bext=\alpha \Phi 0z\u0302/A$, may be written as

where $A=\pi (R>2\u2212R<2)$ and $Ah=\pi R<2$ are the areas of the annulus and of the inner hole. The term $\u223c1/R$ is the particular solution needed when $n\u22600$. The unknown constants $(a,b,\beta )$ are determined so that $A$ is bounded for $R<R<$, and continuous at $R=R<$ and that the magnetic field $B=\u2207\xd7A$ is continuous at $R=R>$. We state them for completeness

with $M=I0>K2<\u2212I2<K0>$ and $I\mu >,<=I\mu (R>,<\lambda )$ and the same for $K\mu >,<$. The corresponding magnetic field, $B=B(R)z\u0302$, in and around the annulus is

The constant *β* gives the flux strength in the inner hole in a similar way as *α* does for the annulus.

##### b. An s-wave superconducting annulus

For the *s*-wave annulus, Eqs. (49) and (53) give a good account for the response to an external field. The free energy of the superconducting annulus can, thus, be evaluated using Eqs. (39) and (40), and the results compare well with results obtained by SuperConga. By inspection, one sees that the functional form of the free-energy correction due to an external field (*α*) and a possible phase winding (*n*) will be^{130}

We determine the factors $\delta \Omega 0,\omega \chi $, and *ω _{B}* numerically. They depend both on geometry and on temperature via the temperature dependence of the penetration depth.

The results for the free energy are shown in Fig. 10 for $T=0.2Tc$ and for a set of different penetration depths. The main purpose here is to verify SuperConga by comparison to the results obtained by linear response. The agreement is excellent for small phase windings and large penetration depths. However, the agreement worsens with decreasing *λ* and increasing *n*. This is a numerical issue that is remedied by a finer spatial resolution of the simulation. This highlights that one needs verify results from SuperConga by varying the parameters that are related to numerical accuracy. These are typically spatial resolution of the geometry, angular resolution of the Fermi-surface averages, and the cutoff in the frequency sums.

##### c. A d-wave superconducting annulus

For the *d*-wave superconductor, we solve Eqs. (7), (8), and (10) self-consistently. We set the zero-temperature penetration depth to $50\xi 0$ as most *d*-wave superconductors are extreme type-II superconductors. The results are shown in Fig. 11 for zero and one phase winding locked into the annulus. Well above $T\u22480.2Tc$ we find that the free energy has the parabolic functional dependence on applied field as for the *s*-wave superconductor. At lower temperatures, the low-energy surface Andreev states dominate, and their paramagnetic response modifies the free energy to host two minima at finite external flux, signaling that it is energetically favorable to form spontaneous currents as discussed above using the solenoid gauge.

##### d. Magnetic moment

The presence of a circulating current in an annulus can be detected as a magnetic moment,

Since the supercurrent, $j(R)$, we compute with SuperConga is confined to flow in a two-dimensional plane, the current density in a stack of layers is given as

i.e., the magnetic moment is given as

where $Nld$ is the thickness of the grain along the *z*-axis assuming *N _{l}* layers separated by the c-axis distance

*d*. The unit for the magnetic moment is $m0=j0\xi 0V$, where

with $B0=\Phi 0/(\pi \xi 02)$ and $V=NldA$. The parameters that determine *m*_{0} are, thus, the geometrical size (and shape) of the object and the magnetic penetration depth specific of the superconducting material. *M*_{0} can also be expressed in the perhaps more familiar form $m0=2\mu BnV$, where $\mu B=\u210f|e|/2me$ is the Bohr magneton and $n=vFpFNF$ the number density of charge carriers. Inserting numerical values for $\Phi 0$ and *μ*_{0} and setting the length scale to *nm*:s for *λ*_{0} and $V$ gives

Inserting the analytical form of the vector potential from Eq. (49) into Eq. (34), we can evaluate the magnetic field (*α*) and phase-winding (*n*) dependence of the magnetic moment of a stack of 2D *s*-wave annuli to be

where

How $d\alpha $ and *d _{n}* depend on the annulus radii and the penetration depth is shown inf Fig. 12. In Fig. 13, we show how the magnetic moment depends on the penetration depth and the forced phase winding as function of applied field. The computed values of the magnetic moment using SuperConga fall right on the analytical ones of Eq. (59).

While the magnetic moment of the *s*-wave annulus holds little surprise, the magnetic moment of the *d*-wave annulus shows trace of a low-temperature transition where screening currents are modified by the paramagnetic nature of the zero-energy Andreev states.^{57,107,131–133} The free energy dependence on temperature and magnetic field may be hard to directly measure whereas the magnetic moment of a *d*-wave annulus could be a possible route to detecting this low-temperature transition. In Fig. 14, we show the magnetic moment of a *d*-wave superconducting annulus as function of external flux at several temperatures. The dimensions are the same as for the s-wave case, i.e., $R>=25\xi 0$ and $R<=10\xi 0$, and the penetration depth is $\lambda 0=50\xi 0$. For low temperatures and small external fields, the signature of the Andreev states is that the magnetic moment of the annulus is reversed compared to that at high fields.

### B. A mesoscopic *s*-wave superconducting disk

As another example, we consider an *s*-wave superconductor in the form of a 2D-disk, with radius $R$, in an external magnetic field. We will allow for arbitrary strength of the applied magnetic field so that multiple vortices can enter the disk. The general solution can and will break the continuous rotational symmetry of the disk around its central axis. Investigating this problem requires that we solve the Eilenberger equation, Eq. (10), with a self-consistently determined order-parameter field, Eq. (3), and vector potential Eq. (7). The magnetic field where the first vortex enters is the first critical field of a superconductor and is referred to as $Bc1$. At this field, it becomes energetically favorable for the superconducting material to host one vortex located in its interior well away from its physical edges. This compared to generating large screening currents flowing at the sample edges. The first critical field will depend on several parameters. Parameters intrinsic to the superconducting material are the penetration depth, $\lambda (T)$, and the coherence length, $\xi (T)$. Both depend on temperature via the temperature dependence of the superconducting energy gap $\Delta (T)$. Extrinsic parameters are, e.g., the physical size and shape of the superconducting sample, which in our example here will simply be the radius $R$ of the disk. Also, temperature and magnetic-field strength are externally controlled. We will assume that our superconductor is a layered 2D-material with negligible inter-layer coupling. As the field is increased above $Bc1$, a cascade of mesoscopic vortex configurations are possible, which are separated by relatively small energy barriers. At high fields, a vortex lattice may be established.^{112} Reaching the second critical field, $Bc2$, bulk superconductivity ultimately vanishes. With our framework SuperConga, we can treat the three relevant length scales of this problem, i.e., $\xi 0,\lambda 0$, and $R$, on the same footing at general temperatures and at general magnetic fields up to $Bc2$ (and in principle $Bc3$). Throughout this section, we focus on a disk with a radius fixed to 25*ξ*_{0} at a temperature set to $0.2Tc$ if not otherwise stated. The analytic results that we give serve as an important tool to verify the results extracted from SuperConga.

Let us first focus on $Bc1$ which was nicely treated by Fetter in Ref. 129 for a thin-film superconducting disk. In the case that the penetration depth is the largest scale in the system, the back-coupling of the current can be omitted. The corresponding superfluid momentum is

Calculating the free-energy density correction due to the superfluid momentum above with $(n=\u22121)$ and without (*n* = 0), a single vortex gives

where $\Omega bulk$ is the bulk free-energy contribution Eq. (38) and

gives the free-energy contribution from the induced supercurrent due to the combined effect of the external field and to the phase-winding of $2\pi $ around the vortex core. The term $\epsilon core$, in principle, consists of two terms. The first is proportional to $ln\u2009(R/Rc)$, with the vortex-core radius $Rc$, due to the $1/R$ dependence of $ps$. The second, *ε*_{0}, is due to the gradient-energy terms that come from the suppression of the order-parameter amplitude needed to accommodate the singularity in the vortex-core center. This second term is not captured in the simple London–Maxwell theory as it requires a theory that includes also the coherence-length scale *ξ*_{0} in the calculations. In Fig. 15, we show calculations made with SuperConga for $\lambda 0=\u221e$. As seen, numerical results fall straight on the theoretical prediction for $(\Omega (\alpha )\u2212\Omega bulk)/A$ given by Eqs. (62) and (63). The main physical observable to extract is $Bc1$, defined as the field where having a vortex becomes energetically favorable compared to the vortex-free case. For this disk at this temperature, $Bc1$ occurs at $\alpha \u22483.8$.

When the penetration depth is comparable to or smaller than the size of the disk, we need to analyze Eqs. (49)–(53) in the limit that $R<\u21920$. The result is

with

The functional form of Eqs. (65) and (66) allows an explicit evaluation of the free-energy contributions $\delta \Omega kin$ and $\delta \Omega magn$ in the case of no vortictity (*n* = 0) and $\delta \omega kin$ is modified for finite *λ* in the following:

with

In the limit that $x\u21920$, corresponding to $\lambda \u2192\u221e,\u2009k(x)\u21921/16$, and $b(x)\u21920$, and we recover the prefactor $\delta \omega kin$ [Eq. (64)] as we should.

The vector potential, Eq. (65), is finite in the origin. This is not the case for the magnetic field, Eq. (66), which diverges as $ln\u2009(R)$ when $R\u21920$ in the origin.^{32,33} This divergence leads to that superconductivity must vanish locally in the vortex core (see, e.g., Fig. 4). The suppression of superconductivity happens over a couple of coherence lengths, and we determine the order-parameter profile using SuperConga. The functional form of the free energy, Eqs. (62) and (63), holds also for finite screening but with $\delta \omega kin(\lambda )$ and with a modified coefficient for linear-in-field term in Eq. (63) [see also Eq. (54)]. In Fig. 16, we show how the free-energy of a single vortex core, i.e., $\delta \omega kin(\lambda )\epsilon core$, depends on temperature. The penetration depth is considered for the values $\kappa =(2,20,\u221e)$, and the external field is zero. The temperature dependence of *ξ* and *λ* gives that close to $Tc$, in the regime of validity of Ginzburg–Landau theory, the physical size of the disk will be the smallest length-scale of the system. In this limit, the relative contribution of the vortex core to the free energy dominates, as seen in Fig. 16. Moving to lower temperatures Ginzburg–Landau theory is strictly speaking no longer applicable and we need to resort to the quasiclassical theory. Below $T=0.8Tc$, the vortex-core contribution for this geometry becomes less dominant and settles to be less than 5% of the total free energy. Instead, the free energy mainly depends on how well the supercurrents originating from the phase-winding around the vortex core are screened, i.e., focus is primarily on the ratio $\lambda 0/R$. We consider two cases of screening: a typical type-II superconductor with $\kappa =20.0$ for which the penetration depth compare to the size of the disk, and a marginal type II superconductor with $\kappa =2.0$ with $\lambda 0\u226aR$. In Fig. 16, we see that strong screening $(\kappa =2)$ reduces the energy cost of a single vortex compared to the case of moderate screening $(\kappa =20)$ so that we can extract $\delta \omega kin\epsilon core(\kappa =2)\u22480.5\delta \omega kin\epsilon core(\kappa =20)$.

In Fig. 17, we show a large set of computed free energies as a function of applied external magnetic field. The magnetic field $B=\alpha \Phi 0/A$ ranges from *α* = 0 to 200 so that we pass $Bc1$, where the first vortex enters, to $Bc2$ where the superconducting state is no longer stable and the disk becomes normal. Starting with the first critical field, $Bc1$, we lift results from the literature.^{113,129} For a superconducting disk of radius $R$ where $\lambda 0/R\u223c1$, the first critical field is given as

with *ϵ*_{0} due to the gradient terms. For a *pancake vortex* in the extreme 2D-case with $\lambda 0/R\u226a1$, the first critical field for a vortex in the 2D superconducting sheet is given as^{60,134}

The lower critical field, $Bc1$, that we extract correspond well to earlier theory if we associate the case *κ* = 20 to Eq. (71) and *κ* = 2 to Eq. (72). Extracting from the data in Fig. 17, we find

The numerical prefactor $4.4=ln\u2009(R/Rc)+\epsilon 0$ contains two unknown parameters $Rc$ and *ε*_{0} for *κ* = 20. On the other hand the prefactor in the case *κ* = 2, $1.62=ln\u2009\kappa +\epsilon 0$ contains only one unknown and we extract $\epsilon 0=1.62\u2212ln\u2009\kappa \u22480.93$. Taking the results from Fig. 16, we can determine the nominal vortex-core size $Rc$ as $ln\u2009\u2009Rc=ln\u2009R+2\epsilon 0(\kappa =2)\u22124.4\u22480.68$ or $Rc=2.7\xi 0$. This value is quite reasonable estimate for a vortex-core size if we compare with the order-parameter field displayed in Fig. 4 [or in the panels (b)–(e) in Fig. 17].

For large external fields, the upper critical field, $Bc2$, is given by the field scale when vortex cores start to overlap,^{33} i.e.,

We add a prefactor, *η*, as in our finite geometry, there is always a vortex-free surface layer around the perimeter of the sample that reduces the effective vortex-lattice area. The vortex-free ring at the perimeter of our disk has a width $\u223c5\xi 0$ and, thus, we can estimate $\eta =(20/25)2=0.64$. The upper critical field for our disk would be at $\alpha =12(R/\xi 0)2=312.5$ but we find that superconductivity vanishes at $\alpha \u2248200$. This reduction of $Bc2$ captured (by construction) by the smaller effective area of the vortex as we described above.

For intermediate fields, $Bc1<B<Bc2$, a vortex lattice emerges. In Fig. 17, we follow how the vortex lattice is formed for the two cases, $\kappa =2,20$. It is a very rich system with several local free-energy minima corresponding to different metastable configurations of vortices in the disk. For a stronger type-II superconductor, *κ* = 20, the lattice structure is less pronounced and only at high fields can we see a square lattice evolving in the center of the disk. On the route toward higher fields, one can find transitions between states comprising of distinct number of vortices. In the following example below, we come back to this and connect to recent experiments.

For the marginal type-II superconductor, the vortex lattice is established quite rapidly in increasing field above $Bc1$. The vortex lattice is triangular with the outer most layer of vortices creating an ever denser necklace of vortices following the perimeter.

### C. Spectroscopy on superconducting nanostructures

As we have seen above, the vortex physics in a finite geometry might be highly modified compared to a bulk system. In particular, the number of vortices and their separation do not only depend on the external field strength, but also on the geometry, the nature of boundaries, and the relation between system size and superconducting length scales, i.e., the coherence length and the penetration depth. This problem has attracted considerable interest both theoretically and experimentally during the last few years.^{3,19–21,55,88–96}

Figure 18 shows the arrangement of 1–6 vortices in a finite and irregular-shaped superconducting grain, inspired by the experiment of Timmermans *et al.*^{55} Finding the arrangement and number of vortices for a given magnetic field is a difficult task due to the complicated magnetic field dependence of the free energy. For this reason, Fig. 18(a) shows the local density of states at the minimum of the free energy. The magnetic field dependence of the free energy is shown in Fig. 18(b) for every number of a vortex. Each vortex follows a parabolic dependence on the magnetic field with a minimum marked as a triangle. The theoretically obtained vortex positions are in qualitative agreement with the experiment results in Ref. 55.

The transition between different vortices entering the grain can be visualized by measuring the zero-bias conductance close to the edge of the grain [Fig. 2(b) in Ref. 55]. Increasing the magnetic field starting from zero, screening currents prevent the magnetic field to penetrate the superconductor. The screening currents are responsible for a reduction of the superconducting gap and an increase in the density of states. The magnetic-field dependence of the density of states close to the edge is shown in Fig. 18(c) for each number of a vortex. The entrance of a vortex in the grain is accompanied by a reduction of the density of states at zero energy. These abrupt reduction of the local density of states are related to the reduction in the zero-bias conductance measurements in Ref. 55.

## VI. SUMMARY

SuperConga is an open-source framework for simulating mesoscopic superconducting grains within the quasiclassical theory of superconductivity. Its main strengths are ease-of-use, speed, and modularity. A Python frontend enables the user to quickly set up simulations of multi-component singlet superconducting grains in general 2D-geometries and in an applied magnetic field. The framework solves for the order parameters self-consistently and includes the back-coupling of the vector potential due to induced supercurrents. Real-time visualization is provided during simulations, and the Python frontend includes functionality for additional data analysis and visualization, as well as interactive spectroscopy of the local density of states.

The SuperConga framework is free to download under the GNU LGPL v3 license or higher, from its GitLab repository,^{83} https://gitlab.com/superconga/superconga. An extensive user manual has been published online,^{80} containing numerous pedagogical examples, tutorials, and guides. The framework has been developed and tested for Unix-based environments and generally runs on most modern laptops, desktops, as well as in cluster environments. The framework relies on the high-performance capabilities offered by GPU acceleration and CUDA^{78,79} and, therefore, requires a CUDA-capable device (i.e., NVIDIA GPU)^{135} to run. To the best of our knowledge, SuperConga is the only open-source code that uses quasiclassical theory to describe mesoscopic superconducting grains and ballistic devices. It therefore adds the missing “mesoscopic link” between more phenomenological methods, such as London–Maxwell^{136} and Ginzburg–Landau,^{117,137,138} and fully microscopic methods, such as density-functional theory^{139–143} and tight-binding approaches.^{144–147}

This paper outlines the functionality of version 1.0 of SuperConga. We foresee that the project website will evolve continuously and that a community of developers can be formed to enable further improvements that increase the functionality and scope of the framework to include, e.g., impurity scattering, spin degrees of freedom, multiband physics, more general boundary conditions, and non-equilibrium phenomena.

## ACKNOWLEDGMENTS

We thank Anton B. Vorontsov and Kevin Marc Seja for valuable discussions. We thank the Swedish Research Council (VR) and the European Research Council (ERC) for financial support. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at the computing centers C3SE, HPC2N, and NSC, funded by the Swedish Research Council through Grant Agreement No. 2022-06725.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Patric Holmvall and Niclas Wall Wennerdal contributed equally to this work.

**Patric Holmvall:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Niclas Wall Wennerdal:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Project administration (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Mikael Håkansson:** Conceptualization (equal); Methodology (equal); Software (lead); Visualization (equal). **Pascal Stadler:** Investigation (equal); Software (equal); Validation (equal); Visualization (equal). **Oleksii Shevtsov:** Conceptualization (equal); Formal analysis (supporting); Methodology (supporting); Software (supporting). **Tomas Löfwander:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Mikael Fogelström:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that supports the findings of this study are available within the article and its supplementary material.

### APPENDIX A: DEPENDENCIES

This appendix lists external dependencies and packages that SuperConga offloads functionality to. For a more detailed description, refer to the installation guide in the main README of the repository.^{83} SuperConga uses git as a version control system, with the central repository being hosted on GitLab.^{83} The SuperConga backend is written in C++, with most of the heavy lifting being done by CUDA^{78,79} kernels and Thrust^{148–150} algorithms running on the GPU. This functionality is provided by the CUDA toolkit, which also provides the NVCC compiler for compiling the CUDA code. Alternatively, Clang can be used for compilation. To run the compiled code, a CUDA-enabled device (an NVIDIA GPU) is required. Note that we use CMake^{151} as the build-system generator, and Ninja^{152} as the build system, but that it is possible to replace either of these with another choice. Linear algebra on the CPU uses Armadillo,^{153,154} and all tests are implemented in the doctest^{155} framework. The real-time plotting is done using OpenGL^{156,157} and ArrayFire-Forge.^{158} The package JsonCpp^{159} is used to handle JSON^{160} objects on the backend side, and HighFive^{161} for HDF5^{162} support. The SuperConga frontend requires a Python3 environment, in particular with the modules NumPy,^{163} Matplotlib,^{164} h5py,^{165} and pandas.^{166}

### APPENDIX B: THE SIMULATION-PARAMETER FILE

The simulation configuration is specified in a simulation_config.json file. It is a JSON^{160} file, which consists of a number of key-value pairs on the form “key”: value. The value of a key can be more key-value pairs, numbers, strings, etc. Below is the JSON file with the parameters used in the demonstration example in Sec. IV.

{ |

“physics”: { |

“temperature”: 0.5, |

“external_flux_quanta”: 1.5, |

“penetration_depth”: 5.0, |

“crystal_axes_rotation”: 0.0, |

“gauge”: “symmetric”, |

“charge_sign”: −1, |

“order_parameter”: { |

“s”: { |

“critical_temperature”: 1.0, |

“initial_phase_shift”: 0.0, |

“initial_noise_stddev”: 0.0, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: −1.0 |

} |

] |

} |

} |

}, |

“geometry”: [ |

{ |

“disc”: { |

“add”: true, |

“center_x”: 0.0, |

“center_y”: 0.0, |

“radius”: 15.0 |

} |

} |

] , |

“numerics”: { |

“convergence_criterion”: 1e-05, |

“energy_cutoff”: 16.0, |

“norm”: “l2”, |

“num_fermi_momenta”: 32, |

“num_iterations_burnin”: 0, |

“num_iterations_max”: 10000, |

“num_iterations_min”: 0, |

“points_per_coherence_length”: 20.0, |

“vector_potential_error”: 1e-06 |

}, |

“accelerator”: { |

“name”: “congacc” |

}, |

“misc”: { |

“data_format”: “h5”, |

“load_path”: “”, |

“save_frequency”: −1, |

“save_path”: “data/examples/swave_disc_vortex”, |

“verbose”: true, |

“visualize”: true |

} |

} |

{ |

“physics”: { |

“temperature”: 0.5, |

“external_flux_quanta”: 1.5, |

“penetration_depth”: 5.0, |

“crystal_axes_rotation”: 0.0, |

“gauge”: “symmetric”, |

“charge_sign”: −1, |

“order_parameter”: { |

“s”: { |

“critical_temperature”: 1.0, |

“initial_phase_shift”: 0.0, |

“initial_noise_stddev”: 0.0, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: −1.0 |

} |

] |

} |

} |

}, |

“geometry”: [ |

{ |

“disc”: { |

“add”: true, |

“center_x”: 0.0, |

“center_y”: 0.0, |

“radius”: 15.0 |

} |

} |

] , |

“numerics”: { |

“convergence_criterion”: 1e-05, |

“energy_cutoff”: 16.0, |

“norm”: “l2”, |

“num_fermi_momenta”: 32, |

“num_iterations_burnin”: 0, |

“num_iterations_max”: 10000, |

“num_iterations_min”: 0, |

“points_per_coherence_length”: 20.0, |

“vector_potential_error”: 1e-06 |

}, |

“accelerator”: { |

“name”: “congacc” |

}, |

“misc”: { |

“data_format”: “h5”, |

“load_path”: “”, |

“save_frequency”: −1, |

“save_path”: “data/examples/swave_disc_vortex”, |

“verbose”: true, |

“visualize”: true |

} |

} |

Under physics are all parameters directly related to the physics of the simulation, except the geometry which is specified separately. The temperature is given in units of the highest critical temperature of the order-parameter components, $Tc\u2261max{Tc,\Gamma}$. The external_flux_quanta is given in units of $\Phi 0$. The penetration_depth is given in units of *ξ*_{0}. If the penetration depth is zero or negative, it is considered to be infinite and the induced vector potential will be zero. The crystal_axes_rotation, i.e., the rotation of the crystal *ab*-axes relative to the simulation coordinate system, is given in units of a full turn, i.e., $2\pi $. The choice for the gauge is either symmetric, landau or solenoid. The sign of the carrier charge is controlled by charge_sign.

The field order_parameter is a list of components. The supported order-parameter components are $dx2\u2212y2$, *d _{xy}*, $gxy(x2\u2212y2)$, and

*s*. The critical_temperature is only meaningful with multiple components as it is the

*relative*$Tc,\Gamma $ that is relevant. Each component is initialized as

where $\delta \Gamma $ is controlled by initial_phase_shift, and *Z* is complex additive Gaussian noise with a standard deviation $\sigma \Gamma $ controlled by initial_noise_stddev in units of $2\pi kBTc$. The field vortices is list of vortices, with each vortex being specified by its center position, $C\Gamma =(center_x,center_y)$, and phase winding $w\Gamma =phase_winding$. Furthermore, $\epsilon (R)=atan2(y,x)$ is the polar angle of the coordinate $R$. We may add additional components by adding entries to the “order_parameter” field in the above JSON file. For example, we may add a subdominant *d _{xy}*-wave component initialized with a phase shift of $\pi /2$, some Gaussian noise, and no vortices.

“order_parameter”: { |

“s”: { |

“critical_temperature”: 1.0, |

“initial_phase_shift”: 0.0, |

“initial_noise_stddev”: 0.0, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: -1.0 |

} |

] |

}, |

“dxy”: { |

“critical_temperature”: 0.1, |

“initial_phase_shift”: 0.25, |

“initial_noise_stddev”: 0.01, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: -1.0 |

} |

] |

} |

} |

“order_parameter”: { |

“s”: { |

“critical_temperature”: 1.0, |

“initial_phase_shift”: 0.0, |

“initial_noise_stddev”: 0.0, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: -1.0 |

} |

] |

}, |

“dxy”: { |

“critical_temperature”: 0.1, |

“initial_phase_shift”: 0.25, |

“initial_noise_stddev”: 0.01, |

“vortices”: [ |

{ |

“center_x”: 0.0, |

“center_y”: 0.0, |

“winding_number”: -1.0 |

} |

] |

} |

} |

The geometry is built of three different kinds of components; disc, regular_polygon, and polygon. A disk is specified by its center which is set by center_x and center_y, and its radius. All parameters are in units of *ξ*_{0}. A regular_polygon is specified by its center coordinates center_x and center_y, the number of edges num_edges, its rotation, and the side_length. The side length and center is given in units of *ξ*_{0}, and the rotation in units of a full turn, i.e., $2\pi $. A polygon is specified by the *x*- and *y*-coordinates of its vertices, vertices_x and vertices_y respectively. The vertices should be given in counterclockwise order. All components have a parameter add, which determines if the area of the component should be added (true), or subtracted (false), from the geometry. For more general geometries it is possible to add and subtract several objects by adding to the geometry field. For example, we may put a square hole of side length $5\xi 0$ in the disk at $10\xi 0$ to the right of the center:

“geometry”: [ |

{ |

“disc”: { |

“add”: true, |

“center_x”: 0.0, |

“center_y”: 0.0, |

“radius”: 15.0 |

} |

}, |

{ |

“regular_polygon”: { |

“add”: false, |

“center_x”: 10.0, |

“center_y”: 0.0, |

“num_edges”: 4, |

“rotation”: 0, |

“side_length”: 5.0 |

} |

} |

], |

“geometry”: [ |

{ |

“disc”: { |

“add”: true, |

“center_x”: 0.0, |

“center_y”: 0.0, |

“radius”: 15.0 |

} |

}, |

{ |

“regular_polygon”: { |

“add”: false, |

“center_x”: 10.0, |

“center_y”: 0.0, |

“num_edges”: 4, |

“rotation”: 0, |

“side_length”: 5.0 |

} |

} |

], |

Note that the order of the components is important if any component is subtracted.

Under numerics are parameters controlling the numerics of the simulation. The convergence_criterion determines when the simulation is considered converged. When all computed quantities; each order-parameter component, the charge current density, the induced vector potential, and the free energy, has a residual [as given by Eq. (27)] smaller than convergence_criterion the simulation should stop. The energy_cutoff is the largest energy, in units of $2\pi kBTc$, to include in the energy sums. The *p*-norm to use when computing the residuals is controlled by norm, with $p\u2208{1,2,\u221e}$. The num_fermi_momenta is the number of discrete momenta to use when approximating the Fermi-surface averages $\u27e8\cdots \u27e9pF$. Note that the terms in the Fermi-surface averages are computed sequentially in contrast to the terms in the energy sum which are computed in parallel on the GPU.

The field num_iterations_burnin controls the number of iterations during which the order-parameter and vector potential are kept constant, while the coherence functions are converged and updated self-consistently, in particular on the boundary, i.e., $\gamma \u2202D$ and $\gamma \u0303\u2202D$. Once the burn-in period is over, all quantities are updated self-consistently again. This is particularly useful when resuming a previous simulation from file, since the order parameter and vector potentials are stored, but the coherence functions are not (due to occupying too much disk space). If the number of burn-in iterations is negative, it will run until convergence. The maximum and minimum of iterations to run is controlled by num_iterations_max and num_iterations_min, respectively. When num_iterations_max iterations has been run the simulation stops regardless of whether convergence has been achieved. Similarly, the simulation does *not* stop until num_iterations_min iterations has been run even if the simulation is considered converged. This is useful for avoiding endless simulations, or prematurely stopping a simulation that is moving slowly on a plateau, respectively. The spatial resolution of the simulation is set by points_per_coherence_length. The vector_potential_error is explained in Appendix D.

All parameters of the accelerator are explained in Appendix E.

Finally, under misc are miscellaneous parameters not affecting the simulation *per se*. data_format specifies which data format to use when saving files. The options are h5 and csv yielding compressed HDF5^{162} files or plain-text CSV files, respectively. Both data formats are standard and have excellent support in, e.g., MATLAB^{167} and Python (via h5py^{165} and pandas^{166}) By setting load_path SuperConga will read the files located there and use them as the initial values for the order parameter and the vector potential. How often to save the results are controlled by save_frequency. If save_frequency is negative the results will only be saved at the end of the simulation, and if it is zero the results will not be saved at all, otherwise it will save every *n*th iteration. The results are saved to save_path. If verbose is true the progress of the simulation will be written to the terminal, and if visualize is true the simulation will be visualized live.

### APPENDIX C: SOLVING THE RICCATI EQUATIONS

Using the trajectory coordinate, *s*, Eq. (25) in standard form is

with the dependence on $pF$ and *ε* dropped, and

for brevity. Instead of using the stepping method,^{70,168} i.e., using an analytical formula for $\gamma (s)$ but approximating $\Delta (s)$ and *z*(*s*) as piecewise constant, we solve Eq. (C1) along a trajectory using the implicit mid-point method, which we find to be both faster and more accurate. The implicit mid-point method has the following update rule:

with the discretized trajectory coordinate *s _{j}* =

*jh*, where

*j*is an integer and

*h*is the length between points. This yields a quadratic equation in $\gamma j+1$ to solve for each step along the trajectory;

with the solution,

Half-step values, $\Delta j+12$ and $zj+12$, are calculated by linear interpolation. Nearest-neighbor extrapolation is used at the boundary. If we have $\gamma (sj)$ at one point, we obtain $\gamma (sj+1)$ via Eq. (C9). The solution along the whole trajectory is obtained by starting with the initial boundary value $\gamma (s0=Rmin)$ and stepping along the trajectory to obtain $\gamma (sj)$ with *j* > 0. The stepping proceeds until a second boundary is reached (denoted $Rmax$, see Fig. 1). We compute $\gamma \u0303(s)$ by solving Eq. (26) in the same manner but stepping in the opposite direction from $Rmax$ to $Rmin$.

### APPENDIX D: COMPUTING THE VECTOR POTENTIAL

The induced vector potential is obtained by solving Eq. (7), in the Coulomb gauge, using the Green's function of the 2D Laplacian operator. The current is computed on a *N* × *N* grid, with every grid cell having the area *h*^{2}. We treat the current as being piecewise constant within each grid cell. With our choice of units, this yields

*G* is a normal matrix, so it can be eigendecomposed as $G=V\Lambda V\u22a4$, where *V* is a unitary matrix with (normalized) eigenvectors as its columns, and Λ is a diagonal matrix with the eigenvalues as diagonal elements. Equivalently,

where $Vm(k)$ is the eigenvector corresponding to the eigenvalue *λ _{k}*. Thus

and the cross correlation with the eigenvectors can be done separately. The eigendecomposition is done with Armadillo.^{153,154} In order to increase performance, only $N*$ eigenvalues are kept. It is computed as the smallest integer obeying

with the eigenvalues sorted $|\lambda 0|\u2265|\lambda 1|\u2265\cdots \u2265|\lambda N\u22121|$, and $\epsilon >0$ being a tolerance chosen by the user via the parameter vector_potential_error in the configuration file. If *ε* = 0, no eigendecomposition is done, and the current density is cross-correlated with the full 2D Green's function.

Finally, the mean of the charge-current density is, in general, not zero due to numerical errors or bad initialization. This has the undesired effect of introducing a constant term in the induced vector potential, which gives rise to an overall phase gradient through the grain and an associated current. For small penetration depths, this can prevent the simulation from converging. In order to remedy this, the mean of the charge-current density is subtracted before solving Eq. (7).

### APPENDIX E: CONVERGENCE ACCELERATORS

SuperConga solves Eqs. (7) and (20) self-consistently. With our choice of units, the solution to Eq. (7) is given by

in the Coloumb gauge, where $G(R,R\u2032)$ is the Green's function of the 2D Laplacian. In order to simplify the notation we concatenate the left- and right-hand sides of Eqs. (20) and (E1), yielding

Note that the vector potential is scaled by $\kappa 2$ in order to have similar magnitudes of all elements in $x$. SuperConga provides several different methods of solving Eq. (E2). Namely, basic Picard iterations, Polyak's^{81} “small heavy sphere,” a variant of the Barzilai–Borwein method,^{82} and a custom method. The user controls which method to use by setting the name in the accelerator in Listing 10. By only setting the name, all internal parameters of the accelerator will be set to default values. The user can change the internal parameters by specifying them in the configuration file, or via the CLI.

In the following discussion of different accelerators, it is convenient to define the difference between the right- and left-hand sides of Eq. (E2),

which enters all accelerator methods.

##### 1. Picard

The simplest way of solving Eq. (E2) is by using Picard iterations,

where $\alpha >0$ is the step size. Picard iterations can converge very slowly, or not at all. Hence, the need for a method accelerating the iterations, and preferably stabilizing them. Calling the Picard method an accelerator is a misnomer as it is the baseline. The other methods should improve on this method. How to use this accelerator is shown in Listing 13.

##### 2. Polyak

In Polyak's method, $x$ behaves similarly to the position of a particle moving through a viscous fluid, with a low Reynold's number, and a potential field [Eq. (E4)] taking the role of the negative potential gradient],

where $\beta \u2208(0,1)$ is the drag, and $\alpha >0$ is the step size. The velocity is initialized to zero, $v0=0$. How to use this accelerator is shown in Listing 14.

##### 3. Barzilai–Borwein

The Barzilai–Borwein (BB) method^{82} is similar to the Picard method but with an adaptive step size,

where $sBB1$ and $sBB2$ are two different variants of how to scale the step size. Note that we force the step size to be positive. Alternating between the variants gives better performance on the examples provided by SuperConga than using either one of them on its own. What step size to use during the first iteration is a free parameter and can be set by the user. How to use this accelerator is shown in Listing 15.

##### 4. CongAcc

CongAcc, short for (Super)Conga Accelerator, is our custom made accelerator. Introducing the *component* index, *c*, where a component is either the real or imaginary part of an order-parameter component, or the *x*- or *y*-component of the induced vector potential, in Eq. (E3), the update rule is as follows:

where *S _{C}* is the cosine similarity, $mi,c$ is an exponential moving average of the component difference $di,c$, and

*k*is a constant determining how much the step size should maximally increase (decrease), which can be set by the user. If $si,c<smin$, where $smin$ is the minimum tolerated similarity, the weight $w*$ is obtained by solving

ensuring that the step will be similar to $\delta i,c$. Equation (E17) is solved approximately using the bisection method. Note that we set $w0,c=0$ during the first iteration.

The reasoning is that if $mi,c$ and $di,c$ roughly point in the same direction, then we should have taken a larger step the previous iteration, and thus, the step size is increased. Analogously, if they roughly point in opposite directions, the step size should decrease.

How to use this accelerator is shown in Listing 16.

“accelerator”: { |

“name”: “congacc”, |

“step_size”: 0.5, |

“step_size_factor”: 1.234, |

“cos_similarity_min”: 0.7, |

“step_size_max”: 100.0, |

“step_size_min”: 0.001 |

} |

“accelerator”: { |

“name”: “congacc”, |

“step_size”: 0.5, |

“step_size_factor”: 1.234, |

“cos_similarity_min”: 0.7, |

“step_size_max”: 100.0, |

“step_size_min”: 0.001 |

} |

This accelerator should be regarded as experimental as no thorough analysis of it has been done.

##### 5. Comparing the accelerators

Table I shows the number of iterations needed to reach convergence in all the examples, (a)–(g), for each of the different accelerators available in SuperConga (using default parameters). The examples (a)–(c) are fairly easy as even Picard iterations converge quickly. It is examples (d)–(g) that really show the strength of the more advanced methods. In (d) Polyak, BB and CongAcc provide roughly a speedup of $3\u20137$ compared to Picard. Both Picard and Polyak fail to converge in example (f) and (g). This is due to the small penetration depth in those examples, making the default step sizes too large. They can be made to converge by manually changing their parameters, however. The adaptive methods, BB and CongAcc, have no problem with bad initial step sizes; they will quickly adjust the step size to something appropriate.

SuperConga examples and accelerator comparison . | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

. | Name . | $T/Tc$ . | $\Phi /\Phi 0$ . | κ
. | Symmetry . | Initialization . | Geometry . | Picard . | Polyak . | BB . | CongAcc . |

(a) | dwave_chiral | 0.5 | 0 | $\u221e$ | $dxy+idx2\u2212y2$ | Bulk | Disk: $R=12.5\xi 0$ | 80 | 34 | 25 | 37 |

(b) | dwave_octagon | 0.5 | 0 | $\u221e$ | $dx2\u2212y2$ | Bulk | Octagon: $S=10\xi 0$ | 71 | 33 | 25 | 32 |

(c) | dwave_plus_swave | 0.5 | 0 | $\u221e$ | $dxy+s$ | Bulk | Disk: $R=12.5\xi 0$ | 69 | 36 | 30 | 33 |

(d) | dwave_phase_crystal | 0.1^{a} | 0 | 100 | $dx2\u2212y2$ | Bulk + vortex^{b} | Irregular polygon^{c} | 3027 | 809 | 1054 | 439 |

(e) | swave_abrikosov_lattice | 0.2 | 20 | 10 | s | Bulk + giant vortex^{d} | Square: $S=25\xi 0$ | 1001 | 313 | 211 | 133 |

(f) | swave_disc_meissner | 0.5 | 0.5 | 5 | s | Bulk | Disk: $R=15\xi 0$ | F | F | 20 | 43 |

(g) | swave_disc_vortex | 0.5 | 1.5 | s | Bulk + vortex^{e} | Disk: $R=15\xi 0$ | F | F | 45 | 42 |

SuperConga examples and accelerator comparison . | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

. | Name . | $T/Tc$ . | $\Phi /\Phi 0$ . | κ
. | Symmetry . | Initialization . | Geometry . | Picard . | Polyak . | BB . | CongAcc . |

(a) | dwave_chiral | 0.5 | 0 | $\u221e$ | $dxy+idx2\u2212y2$ | Bulk | Disk: $R=12.5\xi 0$ | 80 | 34 | 25 | 37 |

(b) | dwave_octagon | 0.5 | 0 | $\u221e$ | $dx2\u2212y2$ | Bulk | Octagon: $S=10\xi 0$ | 71 | 33 | 25 | 32 |

(c) | dwave_plus_swave | 0.5 | 0 | $\u221e$ | $dxy+s$ | Bulk | Disk: $R=12.5\xi 0$ | 69 | 36 | 30 | 33 |

(d) | dwave_phase_crystal | 0.1^{a} | 0 | 100 | $dx2\u2212y2$ | Bulk + vortex^{b} | Irregular polygon^{c} | 3027 | 809 | 1054 | 439 |

(e) | swave_abrikosov_lattice | 0.2 | 20 | 10 | s | Bulk + giant vortex^{d} | Square: $S=25\xi 0$ | 1001 | 313 | 211 | 133 |

(f) | swave_disc_meissner | 0.5 | 0.5 | 5 | s | Bulk | Disk: $R=15\xi 0$ | F | F | 20 | 43 |

(g) | swave_disc_vortex | 0.5 | 1.5 | s | Bulk + vortex^{e} | Disk: $R=15\xi 0$ | F | F | 45 | 42 |

^{a}

The temperature is low, $T<T*$, so the solution is a phase crystal.^{57,58}

^{b}

The vortex has a phase winding *n* = −1 and is positioned outside of the grain, yielding a soft phase-gradient through the grain.

^{c}

The shape of the grain is a square, $S=25\xi 0$, with one corner removed, yielding four $[100]$-interfaces, and one pair-breaking $[110]$-interface.^{61}

^{d}

The giant vortex in the middle has a phase winding *n* = −13. It is unstable and decays into a lattice of singly quantized Abrikosov vortices.

^{e}

The vortex in the middle has a phase winding *n* = −1.

Listing 17 shows how to run one of the available examples with a specific accelerator, where <example> and <accelerator> should be replaced with the example and accelerator names, respectively.

### APPENDIX F: UNITS AND BASIS FUNCTIONS

We summarize in Table II our choice of units and the resulting natural scales for observables and other quantities. The available order parameter basis functions are listed in Table III.

Input to quasiclassical theory | ||

$vF$ | Fermi velocity | |

$NF$ | Normal-state density of states per spin at the Fermi level | |

$Tc$ | Superconducting transition temperature (critical temperature) | |

Geometry | ||

$A$ | Total superconducting computational area^{a} | |

$V$ | $=NldA$ | Volume of a stack of N layers with interlayer distance _{l}d |

Natural units | ||

$Tc$ | Temperature in units of critical temperature | |

$2\pi kBTc$ | Energy in units of critical temperature | |

ξ_{0} | $=\u210fvF2\pi kBTc$ | Length in units of superconducting coherence length |

$\Phi 0$ | $=hc2|e|$ | Magnetic flux in units of superconducting flux quantum |

TABLE II. (Continued.) | ||

Derived units and observables | ||

λ_{0} | $1\lambda 02=4\pi e2c2NFvF2$ | Penetration depth. Here we use Gaussian units. In SI units $4\pi /c2\u2192\mu 0$ with $\mu 0\u22484\pi \xd710\u22127J/A2m$ |

κ | $=\lambda 0\xi 0$ | Penetration depth in units of coherence length |

α | $=BextA\Phi 0$ | Symmetric gauge: external magnetic flux penetrating the superconducting computational area in units of the flux quantum solenoid gauge: external magnetic flux penetrating a circular hole centered at origo in units of the flux quantum |

A_{0} | $=\Phi 0\pi \xi 0$ | Vector potential |

B_{0} | $=\Phi 0\pi \xi 02$ | Magnetic flux density |

j_{0} | $=2\pi kBTc|e|vFNF$ | Charge-current density |

m_{0} | $=j0\xi 0V$ | Magnetic moment |

$\Omega A$ | $=(2\pi kBTc)2NFA$ | Free energy |

Input to quasiclassical theory | ||

$vF$ | Fermi velocity | |

$NF$ | Normal-state density of states per spin at the Fermi level | |

$Tc$ | Superconducting transition temperature (critical temperature) | |

Geometry | ||

$A$ | Total superconducting computational area^{a} | |

$V$ | $=NldA$ | Volume of a stack of N layers with interlayer distance _{l}d |

Natural units | ||

$Tc$ | Temperature in units of critical temperature | |

$2\pi kBTc$ | Energy in units of critical temperature | |

ξ_{0} | $=\u210fvF2\pi kBTc$ | Length in units of superconducting coherence length |

$\Phi 0$ | $=hc2|e|$ | Magnetic flux in units of superconducting flux quantum |

TABLE II. (Continued.) | ||

Derived units and observables | ||

λ_{0} | $1\lambda 02=4\pi e2c2NFvF2$ | Penetration depth. Here we use Gaussian units. In SI units $4\pi /c2\u2192\mu 0$ with $\mu 0\u22484\pi \xd710\u22127J/A2m$ |

κ | $=\lambda 0\xi 0$ | Penetration depth in units of coherence length |

α | $=BextA\Phi 0$ | Symmetric gauge: external magnetic flux penetrating the superconducting computational area in units of the flux quantum solenoid gauge: external magnetic flux penetrating a circular hole centered at origo in units of the flux quantum |

A_{0} | $=\Phi 0\pi \xi 0$ | Vector potential |

B_{0} | $=\Phi 0\pi \xi 02$ | Magnetic flux density |

j_{0} | $=2\pi kBTc|e|vFNF$ | Charge-current density |

m_{0} | $=j0\xi 0V$ | Magnetic moment |

$\Omega A$ | $=(2\pi kBTc)2NFA$ | Free energy |

^{a}

We use calligraphic capital letters for quantities related to geometry, such as area $A$ or disk radius $R$.

Order-parameter basis functions . | ||
---|---|---|

Symmetry class for $D4h$ . | Order parameter $\Delta (k)$ . | Basis function $\eta (\phi )$ for circular Fermi surface . |

$A1g$ (s) | 1 | 1 |

$A2g$ ($gxy(x2\u2212y2)$) | $kxky(kx2\u2212ky2)$ | $2\u2009sin\u2009(4\phi )$ |

$B1g$ ($dx2\u2212y2$) | $kx2\u2212ky2$ | $2\u2009cos\u2009(2\phi )$ |

$B2g$ (d) _{xy} | $kxky$ | $2\u2009sin\u2009(2\phi )$ |

Order-parameter basis functions . | ||
---|---|---|

Symmetry class for $D4h$ . | Order parameter $\Delta (k)$ . | Basis function $\eta (\phi )$ for circular Fermi surface . |

$A1g$ (s) | 1 | 1 |

$A2g$ ($gxy(x2\u2212y2)$) | $kxky(kx2\u2212ky2)$ | $2\u2009sin\u2009(4\phi )$ |

$B1g$ ($dx2\u2212y2$) | $kx2\u2212ky2$ | $2\u2009cos\u2009(2\phi )$ |

$B2g$ (d) _{xy} | $kxky$ | $2\u2009sin\u2009(2\phi )$ |

## References

_{2}Cu

_{3}O

_{7-δ}

_{3}/SrTiO

_{3}interface

_{2}Cu

_{3}O

_{7−δ}

_{2}Sr

_{2}CaCu

_{2}O

_{8+δ}

_{2}Sr

_{2}CaCu

_{2}O

_{8+δ}

_{2}Sr

_{2}CaCu

_{2}O

_{8+δ}

_{2}

_{2}RuO

_{4}

_{c}island

_{2}Cu

_{3}O

_{7−δ}island enhanced by a magnetic field

^{3}He

^{3}He-B

^{3}He-B

^{3}He-A