Author: Ville Pyykkönen, Department of Applied Physics, Aalto University.
License: MIT License
If you use the code in your research, please cite V.A.J. Pyykkönen et al., Phys. Rev. Lett. 130, 216003 (2023).
Simulation of superconducting and normal transport through a two-terminal setup with a tight-binding scattering system. Interactions are considered in terms of self-consistent mean-field approximation of the BCS type. The code simulates single-particle observables at the non-equilibrium steady or stationary states. Specific quantities of interest are the particle current, particle number, and superconducting order parameter.
The simulated system consists of two reservoirs and a scattering system, which is connected to the reservoirs via leads. The scattering system and the leads are considered as tight-binding models. The leads are semi-infinite chains of lattice sites while the scattering system tight-binding geometry can be specified freely. The leads can be either normal or superconductors with specified order parameters. They have specified temperatures and chemical potentials. The scattering system is modeled as a tight-binding model with on-site interaction at the mean-field level. A gate potential is used to control the filling of the scattering system states. The mean fields are determined self-consistently with the non-equilibrium dynamics.
The system can be driven out of equilibrium by setting different chemical potentials, temperatures, or superconducting-order-parameter phases for the two leads. The observables such as current and particle number are obtained assuming that the system has reached a steady state. In other words, the transient effects and initial correlations are neglected. The observable expectation values and the self-consistent fields are determined by the non-equilibrium Green's function (NEGF) method.
Four case-specific methods are implemented to solve the expectation values of the observables.
Method 1: Closed System In case the leads are detached from the scattering system, one has a finite system. It can be directly diagonalized and the states filled with Fermi-Dirac statistics assuming the system is in equilibrium.
Method 2: Real-frequency NEGF with time-independent Hamiltonian In case the transport is assumed to be in a steady state and the system Hamiltonian is time-independent, one can solve the frequency-dependent Green's functions of the whole system using the NEGF method and obtain real-time expectation values for the observables by an inverse Fourier transform. This method works in and out of equilibrium.
Method 3: Real-frequency NEGF with time-periodic Hamiltonian In the case of two superconducting leads with respective chemical potential bias, the effective mean-field Hamiltonian is periodic in time even in the stationary state. To take this into account, the Fourier series components of the mean fields have to be considered.
Method 4: Padé summation at finite-temperature equilibrium If the two-terminal setup is at finite temperature and equilibrium, one can use the technique of imaginary-time Green's function, that is, Matsubara summation to obtain expectation values of the observables. Padé summation is a numerically faster converging algorithm to perform this task. The technique was introduced in T. Ozaki, Phys. Rev. B 75, 035123 (2007).
The methods used are described in detail in the following doctoral thesis. The package has also been used in the analyses of peer-reviewed publications Phys. Rev. B 103, 144519 (2021) and Phys. Rev. Lett. 130, 216003 (2023).
The simulations on the package are based on C++, with some features in C++17 standard. Eigen is used for linear algebra and hdf5 for storing the data. Also, example analysis scripts using Python, namely Numpy, Matplotlib, and h5py libraries, are provided.
Clone the source code
Install Eigen and hdf5 and add the appropriate links to the Makefile.
Eigen version 3.3.9 has been used in the code but probably newer versions also work. For Eigen, just download the software and link the location similarly to the example in the Makefile.
For hdf5 you can use the standard libhdf5-dev package if it is installed. Alternatively, you can use e.g. Anaconda package handling to install hdf5 and link the Anaconda libraries similarly to the given examples. Note in modifying the Makefile that the specific order of the -I, -rpath and -L are important.
Use the command 'make all' to compile the code. Optionally, you can add the handle -jN to parallelize the compilation, eg. 'make all -j4' to use 4 threads.
The code is used in two stages: preparation and execution.
Stage 1: the analysis is prepared by specifying the scattering system geometry and two-terminal setup parameters in parameters.cfg and geometry.cfg files. Examples of these configuration files are provided, where also various parameters and defining an analysis are illustrated with examples.
The analysis is prepared by the command ./prepare_analysis ./parameters.cfg ./geometry.cfg, which prepares the data point files. Each data point corresponds to a hdf5 file containing the parameters.
Stage 2: The simulations for data points are run by the process_datapoint executable by ./process_datapoints path_to_data.h5 parameters.cfg scf_params.cfg, where parameters.cfg contains the specific instructions and parameters for the simulation methods and scf_params.cfg specifies the self-consistent field algorithms used in the mean-field solution. The program will also generate a bash file to process the data point files and a sbatch file for cluster processing.
For practical use of the code, please take a look at the configuration files for instructions on how to set the various parameters. Namely, the geometry of the tight-binding model of the scattering system is specified in geometry.cfg and other parameters shown in the figure above are entered in parameters.cfg. How to apply the different numerical methods discussed above and how to modify their parameters is also shown in parameters.cfg.
The simulation preparation for analysis is described in the main file prepare_datapoints.cpp. The configuration files for the details are read by routines defined in config_parser.cpp and file_io.cpp files. The structure of the analysis preparation is shown in the figure below.
The simulation code structure behind the process_datapoint executable is demonstrated in the following figure.
The parts of the two-terminal setup correspond to classes Lead and ScatteringSystem specified in Lead.h and ScatteringSystem.h, respectively. The complete two-terminal setup corresponds to classes TwoTerminalSetup and SSJunction defined in TwoTerminalSetup.h and SSJunction.h. The former is used in most situations but for reasons concerning the solution algorithm, the SSJunction is used in non-equilibrium situations when both leads are superconducting. On top of their parameters, ScatteringSystemand Lead objects are responsible for their own single-particle Green's functions when detached. TwoTerminalSetup and SSJunctionobjects build their complete Green's functions using the Green's functions of the detached parts and the self-energies of the contacts between them and mean-field interaction. Green's functions are initially calculated in Fourier (frequency) representation and are subsequently inverse Fourier transformed to real-time to obtain observable expectation values. The inverse Fourier transform is handled using the da2glob algorithm Terje Espelid, ACM Trans. Math. Softw. 33, 21-es (2004). The algorithm is implemented in da2glob.h/cpp, Interval.h/cpp, Global.h/cpp and quadrature.h/cpp files.
The mean-field potentials are solved self-consistently using the methods defined in ScfSolver class, defined in ScfSolver.h. Specific methods, namely mixing, Good Broyden's method, Bad Broyden's method, and Anderson-Pulay mixing are defined in ScfMethod.h classes. These methods are discussed in the thesis mentioned above.
The equilibrium Green's functions at finite temperatures are often obtained via the method of Matsubara summation. It is based on the fact that instead of the continuous frequency, the discrete imaginary frequency contains the essential information near the origin of the complex frequency plane. However, the summation in its standard form converges slowly with respect to the number of elements. We have implemented the Pade summation algorithm which provides quadratic convergence instead of linear one. The frequencies are obtained by the functions in pade_frequencies.cpp.
In addition to the discussed core elements of the code, there are many utility functions implemented. config_parser.h/cpp implements a parser for the configuration files discussed above while file_io.h/cpp contain functions used to save and load data to hdf5 and csv format. The class defined in Observable.h handles the various Fourier components of observables in the case the system has a time-periodic solution due to bias between superconducting. Moreover, connectivity.h/cpp implement useful graph algorithms to find different disconnected parts from matrices to help speed up practical linear algebra operations. Finally, fd_dist.h/cpp impelement the Fermi-Dirac distribution function used in many analyses.