Bogoliubov-de Gennes equations in real space

Author: Yuki Nagai

This notebook is expensive. Running this notebook takes a few minutes.

We consider a superconductivity in a mean-field level.

Basic theory of superconductivity in real space

Let us consider a Hamiltonian for a fermion system given as \(H = \psi^{\dagger} \hat{\cal H} \psi/2\). The column vector \(\psi\) is composed of \(N\) fermionic annihilation \(c_i\) and creation operators \(c_i^{\dagger}\) \((i = 1,2,\cdots,N)\), \(\psi = (\{ c_i \}^T, \{ c_i^{\dagger} \}^T)\), where \(\{ c_i \} = (c_1,c_2,\cdots,c_N)^T\) and \(\{ c_i^{\dagger} \} = (c_1^{\dagger},c_2^{\dagger},\cdots,c_N^{\dagger})\). The row vector \(\psi^{\dagger}\) is also defined as \(\psi^{\dagger} = (\{c_i^{\dagger} \}^T, \{ c_i \}^T)\). The subscription \(i\) in \(c_i\) or \(c_i^{\dagger}\) indicates a quantum index depending on spatial site, spin, orbital, etc. The “Hamiltonian” matrix \(\hat{\cal H}\) is a \(2N \times 2N\) Hermitian matrix given as

\[\begin{split} \hat{\cal H} = \left( \begin{matrix} \hat{H} & \hat{\Delta} \\ \hat{\Delta}^{\dagger} & - \hat{H}^{\ast} \end{matrix} \right), \end{split}\]

where \(\hat{H}\) is a normal state Hamiltonian and \(\hat{\Delta}\) is a superconducting order parameter.

The BdG equations are regarded as the eigenvalue equation with respect to \(\hat{\cal H}\) expressed as

\[ \hat{\cal H} \vec{f}_{\gamma} = \epsilon_{\gamma} \vec{f}_{\gamma} \]

\[\begin{split} \vec{f}_{\gamma} = \left( \begin{matrix} \vec{u}_{\gamma} \\ \vec{v}_{\gamma} \end{matrix} \right) \end{split}\]

The column vectors \(\vec{u}_{\gamma}\) and \(\vec{v}_{\gamma}\) are \(N\)-component complex vetors. To solve the BdG equations is equivalent to diagonalization of \(\hat{\cal H}\) with a unitary matrix \(\hat{U}\),

\[ \hat{U}^{\dagger} \hat{\cal H} \hat{U} = \hat{D}, \: \: \: \hat{D} = {\rm diag} (\epsilon_1,\epsilon_2,\cdots,\epsilon_{2N}) \]

Mean-fields are calculated by one-particle Green’s functions. In the mean-field framework, the \(2N \times 2N\) matrix Green’s function with a complex frequency is defined as

\[ \hat{G}(z) = [z \hat{I} - \hat{\cal H}]^{-1}. \]

With the unitary matrix \(\hat{U}\), each component of \(\hat{G}(z)\) is epressed as

\[ G_{\alpha \beta}(z) = \sum_{\gamma=1}^{2N} U_{\alpha \gamma} U_{\beta \gamma}^{\ast} \frac{1}{z - \epsilon_{\gamma}} \]

If we set \(z = i \omega_n\) with the Matsubara frequency \(\omega_n = (2n+1)\pi T\), the above formula corresponds to Matsubara temperature Green’s function. The regarded and advanced Green’s functions are, respectively, defined as

\[ \hat{G}^{\rm R}(\omega) = \lim_{\tau \rightarrow 0+} \hat{G}(\omega + i \eta) \]

\[ \hat{G}^{\rm A}(\omega) = \lim_{\tau \rightarrow 0+} \hat{G}(\omega - i \eta) \]

In order to obtain physical quantities (e.g., density of states) from Green’s functions, we introduce the following useful \(2N\)-component unit-vectors \(\vec{e}(i)\) and \(\vec{h}(i)\) \((1 \le i \le N)\), which are, respectively, defined as

\[ [\vec{e}(i)]_{\gamma} = \delta_{i,\gamma}, \: [\vec{h}(i)]_{\gamma} = \delta_{i+N,\gamma} \]

For example, the local density of states with respect to the site \(i\) is given as

\[ N(\omega,i) = -\frac{1}{2\pi i} \vec{e}(i)^T \hat{d}(\omega) \vec{e}(i), \]

with the use of the spectral function \(\hat{d}(\omega)\) defined as

\[ [\hat{d}(\omega)]_{\alpha \beta} \equiv [\hat{G}^{\rm R}(\omega) - \hat{G}^{\rm A}(\omega)]_{\alpha \beta} = - 2\pi i \sum_{\gamma=1}^{2N} U_{\alpha \gamma} U_{\beta \gamma}^{\ast} \delta(\omega - \epsilon_{\gamma}) \]

Two types of mean-fields \(\langle c_i^{\dagger} c_j \rangle\) and \(\langle c_i c_j \rangle\) can be expressed as

\[ \langle c_i^{\dagger} c_j \rangle = - \frac{1}{2 \pi i} \int_{-\infty}^{\infty} d\omega f(\omega) \vec{e}(j)^T \hat{d}(\omega) \vec{e}(i), \]

\[ \langle c_i c_j \rangle = - \frac{1}{2 \pi i} \int_{-\infty}^{\infty} d\omega f(\omega) \vec{e}(j)^T \hat{d}(\omega) \vec{h}(i), \]

with \(f(x) \equiv 1/(e^{x/T} + 1)\).

Matsubara frequency representation

With the use of the analytic continuation, the mean-fields are rewritten as

\[ \langle c_i^{\dagger} c_j \rangle = T \sum_{n=-\infty}^{\infty} \vec{e}(j)^T \hat{G}(i \omega_n) \vec{e}(i) \]

\[ \langle c_i c_j \rangle = T \sum_{n=-\infty}^{\infty} \vec{e}(j)^T \hat{G}(i \omega_n) \vec{h}(i) \]

By solving the linear equations defined as

\[ (i \omega_n \hat{I} - \hat{\cal H}) \vec{x}(i,\omega_n) = \vec{h}(i), \]

the superconducting mean-field is expressed as

\[ \langle c_i c_j \rangle = T \sum_{n = - \infty}^{\infty} \hat{e}(j)^T \vec{x}(i,\omega_n). \]

For example, s-wave superconducting order parameter is defined as

\[ \Delta_i = U \langle c_i c_i \rangle, \]

The self-consistent cycle has the schematic form

\[ \hat{\Delta} \rightarrow \hat{G}(i \omega_n) \rightarrow \langle c_i c_i \rangle \rightarrow \hat{\Delta} \]

Reduced-shifted Conjugation Gradient method

The reduced-shifted conugate-gradient (RSCG) method is numerically efficient to calculate a matrix element of a Green’s function defined as a resolvent of a Hamiltonian operator, by solving linear equations with desired accuracy. The matrix elements with different frequencies are simultaneously obtained. The details are described in J. Phys. Soc. Jpn. 86, 014708 (2017). The RSCG method is provided in RSCG.jl package.

IR basis

The superconducting mean-fielads are expressed as

\[ \langle c_i c_j \rangle = \vec{e}(j)^T \hat{G}(\tau = \beta) \vec{h}(i) = \sum_{l=0}^{N_{\rm IR}-1} U_l(\beta) \vec{e}(j)^T \hat{G}_l \vec{h}(i) \]

Since the Matsubara Green’s function is expressed as

\[ G(i \omega_n) = \sum_{l=0}^{N_{\rm IR}-1} G_l U_l(i \omega_n), \]

one can obtain \(G_l\) by fitting the above equation.

Two-dimensional s-wave superconductor on the square lattice

The normal state Hamiltonian on the tight-binding model is expressed as

\[ \sum_{i,j} t_{ij} c_i^{\dagger} c_j - \mu \sum_i c_i^{\dagger} c_i \]

Here, \(t_{ij} = -t\) if \(j\) is a nearest neighbor for \(i\). The matrix form of the above Hamiltonain is expressed as

\[ \hat{H} = -\mu \hat{I} -t (\hat{T}_{+x} + \hat{T}_{-x} + \hat{T}_{+y} + \hat{T}_{-y}), \]

where \(\hat{T}_{+x}\) is a hopping matrix with respect to \(+x\) direction.

In Julia code, this Hamiltonian matrix is defined as

# Check if a given function called with given types is type stable
function typestable(@nospecialize(f), @nospecialize(t))
    v = code_typed(f, t)
    stable = true
    for vi in v
        for (name, ty) in zip(vi[1].slotnames, vi[1].slottypes)
            !(ty isa Type) && continue
            if ty === Any
                stable = false
                println("Type instability is detected! the variable is $(name) ::$ty")
            end
        end
    end
    return stable
end

using SparseArrays
using LinearAlgebra
function make_x_plushop(Nx,Ny,BC)
    N = Nx*Ny
    Txhop = spzeros(Int64,N,N)
    for ix=1:Nx
        for iy=1:Ny
            i = (iy-1)*Nx + ix
            jx = ix + 1
            jy = iy
            if BC == "PBC"
                jx += ifelse(jx > Nx,-Nx,0)
            elseif BC == "OBC"
            else
                error("BC = $BC is not supported")
            end
            if 1 <= jx <= Nx
                j = (jy-1)*Nx + jx
                Txhop[i,j] = 1
            end
        end
    end
    return Txhop
end

function make_x_minushop(Nx,Ny,BC)
    N = Nx*Ny
    Txhop = spzeros(Int64,N,N)
    for ix=1:Nx
        for iy=1:Ny
            i = (iy-1)*Nx + ix
            jx = ix - 1
            jy = iy
            if BC == "PBC"
                jx += ifelse(jx < 1,Nx,0)
            elseif BC == "OBC"
            else
                error("BC = $BC is not supported")
            end
            if 1 <= jx <= Nx
                j = (jy-1)*Nx + jx
                Txhop[i,j] = 1
            end
        end
    end
    return Txhop
end

function make_y_plushop(Nx,Ny,BC)
    N = Nx*Ny
    Tyhop = spzeros(Int64,N,N)
    for ix=1:Nx
        for iy=1:Ny
            i = (iy-1)*Nx + ix
            jx = ix 
            jy = iy + 1
            if BC == "PBC"
                jy += ifelse(jy > Ny,-Ny,0)
            elseif BC == "OBC"
            else
                error("BC = $BC is not supported")
            end
            if 1 <= jy <= Ny
                j = (jy-1)*Nx + jx
                Tyhop[i,j] = 1
            end
        end
    end
    return Tyhop
end

function make_y_minushop(Nx,Ny,BC)
    N = Nx*Ny
    Tyhop = spzeros(Int64,N,N)
    for ix=1:Nx
        for iy=1:Ny
            i = (iy-1)*Nx + ix
            jx = ix 
            jy = iy - 1
            if BC == "PBC"
                jy += ifelse(jy < 1,Ny,0)
            elseif BC == "OBC"
            else
                error("BC = $BC is not supported")
            end
            if 1 <= jy <= Ny
                j = (jy-1)*Nx + jx
                Tyhop[i,j] = 1
            end
        end
    end
    return Tyhop
end

function make_H_normal(Nx,Ny,μ,BC)
    N = Nx*Ny
    Tx_plushop = make_x_plushop(Nx,Ny,BC)
    Tx_minushop = make_x_minushop(Nx,Ny,BC)
    Ty_plushop = make_y_plushop(Nx,Ny,BC)
    Ty_minushop = make_y_minushop(Nx,Ny,BC)
    HN = sparse(I,N,N)*(-μ)
    t = 1.0
    
    HN += -t*(Tx_plushop + Tx_minushop + Ty_plushop + Ty_minushop)
    return HN
end

@assert typestable(make_x_plushop,(Int64,Int64,String))
@assert typestable(make_x_minushop,(Int64,Int64,String))
@assert typestable(make_y_plushop,(Int64,Int64,String))
@assert typestable(make_y_minushop,(Int64,Int64,String))

The onsite superconducting s-wave order parameter is defined as

function make_Δ(Δ)
    Nx,Ny = size(Δ)
    N = Nx*Ny
    Δmat = spzeros(ComplexF64,N,N)
    for ix=1:Nx
        for iy=1:Ny
            i = (iy-1)*Nx + ix
            Δmat[i,i] = Δ[ix,iy]
        end
    end
    return Δmat
end
@assert typestable(make_Δ,(Matrix{ComplexF64},))

The BdG Hamiltonian is defined as

function make_H_sc(Nx,Ny,μ,Δ,BC)
    HN = make_H_normal(Nx,Ny,μ,BC)
    matΔ = make_Δ(Δ)
    H = [
        HN matΔ
        matΔ' -conj.(HN)
    ]
    return H
end
@assert typestable(make_H_sc,(Int64,Int64,Float64,Matrix{ComplexF64},String))

We have to update the order parameters in the self-consistent cycle. The update method is implemented as

function update_H_sc!(H,Δ)
    matΔ = make_Δ(Δ)
    Nx,Ny = size(Δ)
    N = Nx*Ny
    H[1:N,1+N:2N] = matΔ
    H[1+N:2N,1:N] = matΔ'
end

@assert typestable(update_H_sc!,(AbstractMatrix{ComplexF64},Matrix{ComplexF64}))

Conventional method with Matsubara frequencies

If one wants to use the Matsubara frequencies, one has to introduce the cutoff of the frequencies. Then the mean-fields are expressed as

\[ \langle c_i c_j \rangle = T \sum_{n = - n_c}^{n_c} \hat{e}(j)^T \vec{x}(i,\omega_n). \]

Here, \(n_c\) is a cutoff value. With the use of the maximum frequency \(\omega_c\), we calculate the set of Matsubara frequencies defined as

function calc_ωn(T,ωc)
    M = Int((round(ωc/(T*π)))/2-1)
    println("num. of Matsubara freq: ",2M)
    ωn = zeros(ComplexF64,2M)
    for n=1:2M
        ωn[n] = π*T*(2.0*(n-M-1)+1)*im
    end
    return ωn
end
@assert typestable(calc_ωn,(Float64,Float64))

ωc = 100π
T = 0.01
ωn =  calc_ωn(T,ωc)

num. of Matsubara freq: 9998

9998-element Vector{ComplexF64}:
 -0.0 - 314.06501757937167im
 -0.0 - 314.00218572629984im
 -0.0 - 313.9393538732281im
 -0.0 - 313.87652202015624im
 -0.0 - 313.8136901670845im
 -0.0 - 313.75085831401265im
 -0.0 - 313.6880264609409im
 -0.0 - 313.62519460786905im
 -0.0 - 313.5623627547973im
 -0.0 - 313.4995309017255im
 -0.0 - 313.4366990486537im
 -0.0 - 313.3738671955819im
 -0.0 - 313.3110353425101im
      ⋮
  0.0 + 313.3738671955819im
  0.0 + 313.4366990486537im
  0.0 + 313.4995309017255im
  0.0 + 313.5623627547973im
  0.0 + 313.62519460786905im
  0.0 + 313.6880264609409im
  0.0 + 313.75085831401265im
  0.0 + 313.8136901670845im
  0.0 + 313.87652202015624im
  0.0 + 313.9393538732281im
  0.0 + 314.00218572629984im
  0.0 + 314.06501757937167im

With the use of RSCG.jl, the mean-fields are calculated by

using Pkg
Pkg.add(PackageSpec(name="RSCG", version = "0.1.2"))
using RSCG

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.9/Project.toml`
  No Changes to `~/.julia/environments/v1.9/Manifest.toml`

function calc_Δi!(i,N,H,Δold,T,U,ωn;mixratio = 0.5)
    j = i + N
    Gij = greensfunctions(i,j,ωn,H)
    Δi = U*T*sum(Gij)
    Δi = (1-mixratio)*Δold[i] + mixratio*Δi
    return Δi
end

@assert typestable(calc_Δi!,(Int64,Int64,AbstractMatrix{ComplexF64},Matrix{ComplexF64},Float64,Float64,Vector{ComplexF64}))

function calc_Δ!(Δnew,H,Δold,T,U,ωn;mixratio = 0.5)
    Nx,Ny = size(Δold)
    N = Nx*Ny
    map!(i -> calc_Δi!(i,N,H,Δold,T,U,ωn,mixratio = mixratio),Δnew,1:N) #If you use pmap! instead of map!, you can do the parallel computation.
    return
end

@assert typestable(calc_Δi!,(AbstractMatrix{ComplexF64},Matrix{ComplexF64},Float64,Float64,Vector{ComplexF64}))

Then, we can construct the BdG Hamiltonian matrix.

Nx = 8
Ny = 8
Δ = ones(ComplexF64,Nx,Ny)
Δold = copy(Δ)
Δnew = zero(Δ)
BC = "OBC" #open boundary condition
#BC = "PBC" #periodic boundary condition
U  =-2
μ = -0.2

H = make_H_sc(Nx,Ny,μ,Δ,BC)

128×128 SparseMatrixCSC{ComplexF64, Int64} with 704 stored entries:
⎡⠻⣦⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠢⡈⠱⣦⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠈⠢⡈⠛⣤⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠢⡈⠻⢆⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠢⡈⠻⣦⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠻⣦⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⢄⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠱⣦⡈⠦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⢄⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡌⠛⣤⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⢄⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠻⢆⡈⠢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⢄⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠻⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⣄⎥
⎢⠙⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠑⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠢⡈⠱⣦⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠑⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠛⣤⡘⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠑⢆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠲⡈⠻⢆⡈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠑⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠻⣦⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠻⣦⡈⠢⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠱⣦⡈⠢⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠛⣤⡈⠢⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠻⢆⡈⠢⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡈⠻⣦⎦

The self-consistent loop is defined as

itemax = 100
ix = Nx ÷ 2
iy = Ny ÷ 2
for ite = 1:itemax
    calc_Δ!(Δnew,H,Δold,T,U,ωn)
    update_H_sc!(H,Δnew)
    eps = sum(abs.(Δnew-Δold))/sum(abs.(Δold))
    println("$ite $eps ",Δnew[ix,iy])
    Δold .= Δnew
    if eps < 1e-3
        break
    end
end

1 0.19022941762340817 0.8018647273876987 + 0.0im

2 0.15123442982656063 0.674653022641663 + 0.0im

3 0.11957205212486025 0.589596797219806 - 5.293955920339377e-25im

4 0.09424713927918395 0.5308635372889943 + 3.970466940254533e-24im

5 0.07417361359379145 0.48926476258517093 + 1.9852334701272663e-24im

6 0.05834520068959513 0.4592081174486282 + 1.5220123270975708e-24im

7 0.0458979462622344 0.437146852174492 + 7.3619074517219455e-25im

8 0.03612064429598342 0.42075099653803383 + 4.177262093392789e-25im

9 0.028441883519312965 0.4084438556031036 - 8.261466367873364e-24im

10 0.022409027712510025 0.39913169403298343 - 4.1638204084388034e-24im

11 0.017666264771733595 0.3920398625330299 - 2.0901820103449318e-24im

12 0.013934920454239436 0.38661027271623244 - 1.0285473929214053e-24im

13 0.010997079557253684 0.38243510651602697 - 5.101377933979375e-25im

14 0.00868234585336176 0.37921283781072723 - 1.3138600807668442e-24im

15 0.0068573942739070855 0.37671838596315144 - 9.152074931303016e-24im

16 0.005417746377601394 0.374782341949813 - 4.857278873919537e-24im

17 0.004281510158206325 0.3732763752819832 - 2.4617266614618897e-24im

18 0.0033843513340726023 0.37210271002493217 - 1.2287953791995623e-24im

19 0.0026757083399357364 0.3711865098596343 - 2.0080738332326733e-25im

20 0.002115794136491162 0.3704702697211746 + 1.642941043553352e-25im

21 0.0016732885754822536 0.3699096269357205 + 6.560343992660704e-26im

22 0.001323477168596026 0.3694702916891107 + 4.336709142763973e-21im

23 0.0010469105457190362 0.3691256761125573 + 2.168222222483978e-21im

24 0.000828206739947293 0.3688551232542179 + 5.43040756280999e-21im

IR based method

IR based method for BdG equations is shown in arXiv:2111.13288. In this paper, we use the ultra-fast method (J. Phys. Soc. Jpn. 89, 074703 (2020)). In this notebook, we do not use the ultra-fast method for simplicity.

We use the SparseIR package. The sampling Matsubara frequencies are obtained as

using SparseIR
import SparseIR: valueim

wmax = 10.0
beta = 1/T
basis = FiniteTempBasis(Fermionic(), beta, wmax, 1e-7)
smpl = MatsubaraSampling(basis)
ωn_s = valueim.(smpl.sampling_points, beta)
println("num. of Matsubara freqs. ", length(ωn_s))
smpl_beta = TauSampling(basis; sampling_points=[beta])

num. of Matsubara freqs. 40

TauSampling64 with sampling points:
 100.0

The mean-fields are obtained by

function fit_ir(Gij,smpl_Matsubara,smpl_beta)
    gl = fit(smpl_Matsubara, Gij)
    G0 = evaluate(smpl_beta, gl)
    return -G0[1]
end

@assert typestable(fit_ir,(Vector{ComplexF64},typeof(smpl),typeof(smpl_beta)))

function calc_Δi_ir!(i,N,H,Δold,T,U,ωn,smpl_Matsubara,smpl_beta;mixratio = 0.5)
    j = i + N
    Gij = greensfunctions(i,j,ωn,H)
    G0 = fit_ir(Gij,smpl_Matsubara,smpl_beta)
    Δi = U*G0
    Δi = (1-mixratio)*Δold[i] + mixratio*Δi
    return Δi
end

@assert typestable(calc_Δi_ir!,(Int64,Int64,AbstractMatrix{ComplexF64},Matrix{ComplexF64},Float64,
            Float64,Vector{ComplexF64},typeof(smpl),typeof(smpl_beta)))


function calc_Δ_ir!(Δnew,H,Δold,T,U,ωn,smpl_Matsubara,smpl_beta;mixratio = 0.5)
    Nx,Ny = size(Δold)
    N = Nx*Ny
    map!(i -> calc_Δi_ir!(i,N,H,Δold,T,U,ωn,smpl_Matsubara,smpl_beta,mixratio = mixratio),Δnew,1:N) #If you use pmap! instead of map!, you can do the parallel computation.
    return
end

@assert typestable(calc_Δ_ir!,(Matrix{ComplexF64},AbstractMatrix{ComplexF64},Matrix{ComplexF64},Float64,
            Float64,Vector{ComplexF64},typeof(smpl),typeof(smpl_beta)))

Then, the self-consistent cycle is expressed as

Nx = 8
Ny = 8
Δ = ones(ComplexF64,Nx,Ny)
Δold = copy(Δ)
Δnew = zero(Δ)
BC = "OBC" #open boundary condition
#BC = "PBC" #periodic boundary condition
U  =-2
μ = -0.2

H = make_H_sc(Nx,Ny,μ,Δ,BC)

for ite = 1:itemax
    calc_Δ_ir!(Δnew,H,Δold,T,U,ωn_s,smpl,smpl_beta)
    update_H_sc!(H,Δnew)

    eps = sum(abs.(Δnew-Δold))/sum(abs.(Δold))
    println("$ite $eps ",Δnew[ix,iy])
    Δold .= Δnew
    if eps < 1e-4
        break
    end
end

1 0.1892161509461065 0.802878002123697 + 5.1557014979967645e-15im
2 0.150457979742895 0.6761326661696428 + 8.80983446052722e-15im

3 0.11898342486253458 0.5912904905536116 + 1.1488969651458835e-14im
4 0.0938032473327619 0.53265062730773 + 1.36657787201322e-14im

5 0.07383909072789409 0.4910866442832549 + 1.5569812951557572e-14im
6 0.05809241718246723 0.46103660085752246 + 1.7316895756868512e-14im

7 0.045706075191968086 0.438968929955149 + 1.911084338070148e-14im
8 0.03597424649498475 0.42256138584784747 + 2.0816526501213613e-14im

9 0.02832962425037164 0.41024116094347957 + 2.2474140822765472e-14im
10 0.022322530136773537 0.4009163619921634 + 2.419859729045251e-14im

11 0.017599332806191432 0.39381318906681373 + 2.596154136173753e-14im
12 0.013882932462044332 0.38837383352689137 + 2.765655138431156e-14im

13 0.010956570683617188 0.38419047851437443 + 2.943415894266075e-14im
14 0.008650706161258915 0.3809614754470959 + 3.1318193734686775e-14im

15 0.006832626218315066 0.37846154488904854 + 3.312705753896783e-14im
16 0.0053983225650163745 0.3765210964282115 + 3.498846581066755e-14im

17 0.0042662593924301375 0.3750116105891038 + 3.688766937567821e-14im
18 0.003372364638144916 0.3738351578380946 + 3.8794368679948055e-14im

19 0.0026662729485389196 0.37291675708548355 + 4.070300590892883e-14im
20 0.002108364299556271 0.37219877984003114 + 4.261968936617414e-14im

21 0.001667435271313639 0.3716367824732472 + 4.458369766355874e-14im
22 0.00131887096198135 0.37119637973635533 + 4.652715264787316e-14im

23 0.0010432770147928159 0.3708509406768796 + 4.8515735671392577e-14im
24 0.0008253431958382265 0.37057974003225014 + 5.054818747338997e-14im

25 0.0006529806704441708 0.37036665778602185 + 5.242669158398389e-14im
26 0.0005166496548387234 0.37019912241236475 + 5.4379927189697203e-14im

27 0.0004088023162402707 0.37006731470363485 + 5.630857374673937e-14im
28 0.0003234837873061463 0.36996355525718533 + 5.831334529613167e-14im

29 0.00025598222056606067 0.36988183592735135 + 6.026286444587446e-14im
30 0.00020257508040952168 0.3698174438375502 + 6.22325251854875e-14im

31 0.00016031110419138027 0.3697666830800834 + 6.426439743110256e-14im
32 0.0001268713459379015 0.3697266628159944 + 6.624610158104896e-14im

33 0.00010040920450751594 0.3696950766236306 + 6.82222249372604e-14im
34 7.946979714649064e-5 0.3696701575393463 + 7.018612390051859e-14im

Let us show the superconducting order parameter in real space. Since we consider the open boundary condition, the mean-field becomes large at four corner due to an interference effect.

using Plots
plot(1:Nx,1:Ny,abs.(Δnew),st=:contourf,
    xlabel = "ix",
    ylabel = "iy",
    zlabel = "|Delta|")

We plot the local density of states at the center.

M = 1000
σ = zeros(ComplexF64,M)
η = 0.05
σmin = -4.0 + im*η
σmax = 4.0+ im*η
for i=1:M
    σ[i] = (i-1)*(σmax-σmin)/(M-1) + σmin
end

i = (iy-1)*Nx + ix
j = i

Gij1 = greensfunctions(i,j,σ,H) 
plot(real.(σ),(-1/π)*imag.(Gij1),
    xlabel = "Energy [t]",
    ylabel = "Local DOS",)

Two-dimensional topological s-wave superconductor on the square lattice

We consider topological superconductivity. The s-wave superconductor becomes topologically non-trivial when there are Zeeman magnetic fields and Rashba spin-orbit coupling. The Hamiltonian in momentum space is given as

\[\begin{split} \hat{\cal H}(\vec{k}) = \left( \begin{matrix} h_0(\vec{k}) & i \Delta \sigma_y \\ - i \Delta^{\dagger} \sigma_y & - h_0^{\ast}(\vec{k}) \end{matrix} \right), \end{split}\]

with \( h_0(\vec{k}) = \epsilon(\vec{k}) - h \sigma_z + \alpha {\cal L}(\vec{k})\). The band dispersion is \(\epsilon(\vec{k}) = -2 t (\cos k_x + \cos k_y) - \mu\). The Zeeman magnetic field is \(h\). The Rashba spin-orbit term is described by \(\alpha {\cal L}(\vec{k}) = \alpha (\sigma_x \sin k_y - \sigma_y \sin k_x)\). Here, \(\sigma_{x,y,z}\) are the Pauli matrices.

By Fourier transforming this Hamiltonian, we have the Hamiltonian defined in real space.

We difine the Pauli matrices.

const σ0 = [1 0
0 1]
const σx = [0 1
1 0]
const σy = [0 -im
im 0]
const σz = [1 0
0 -1]

2×2 Matrix{Int64}:
 1   0
 0  -1

We can construct the Hamiltonian with the use of the kron:

function make_Htsc_normal(Nx,Ny,μ,BC,h,α)
    N = Nx*Ny
    Tx_plushop = make_x_plushop(Nx,Ny,BC)
    Tx_minushop = make_x_minushop(Nx,Ny,BC)
    Ty_plushop = make_y_plushop(Nx,Ny,BC)
    Ty_minushop = make_y_minushop(Nx,Ny,BC)
    HN = kron(sparse(I,N,N)*(-μ),σ0) 
    HN += kron(sparse(I,N,N)*(-h),σz) #Zeeman magnetic field
    t = 1.0
    
    HN += kron(-t*(Tx_plushop + Tx_minushop + Ty_plushop + Ty_minushop),σ0)
    
    Hax = kron((α/(2im))*(Tx_plushop - Tx_minushop ) ,σy)
    HN += Hax 
    Hay = kron((α/(2im))*(Ty_plushop - Ty_minushop ) ,σx)
    HN += Hay 
    
    return HN
end

@assert typestable(make_Htsc_normal,(Int64,Int64,Float64,String,Float64,Float64))

We also define the superconducting order parameter:

function make_Δtsc(Δ)
    Nx,Ny = size(Δ)
    N = Nx*Ny
    Δmat = spzeros(ComplexF64,N,N)
    for ix=1:Nx
        for iy=1:Ny
            i = (iy-1)*Nx + ix
            Δmat[i,i] = Δ[ix,iy]
        end
    end
    return kron(Δmat,im*σy)
end

@assert typestable(make_Δtsc,(Matrix{ComplexF64},))

Then, the superconducting BdG Hamiltonian is defined as

function make_Htsc_sc(Nx,Ny,μ,Δ,BC,h,α)
    HN = make_Htsc_normal(Nx,Ny,μ,BC,h,α)
    matΔ = make_Δtsc(Δ)
    H = [
        HN matΔ
        matΔ' -conj.(HN)
    ]
    return H
end

@assert typestable( make_Htsc_sc,(Int64,Int64,Float64,Matrix{ComplexF64},String,Float64,Float64))

The update method is implimanted as

function update_Htsc_sc!(H,Δ)
    matΔ = make_Δtsc(Δ)
    N,_ = size(matΔ)
    H[1:N,1+N:2N] = matΔ
    H[1+N:2N,1:N] = matΔ'
end

@assert typestable(update_Htsc_sc!,(AbstractMatrix{ComplexF64},Matrix{ComplexF64}))

The mean-fields are calculated as

function calc_Δitsc_ir!(i,N,H,Δold,T,U,ωn,smpl_Matsubara,smpl_beta;mixratio = 0.5)
    ispin = 1
    ii = (i-1)*2 + ispin
    jspin = 2
    jj = (i-1)*2 + jspin + N
    
    Gij = greensfunctions(ii,jj,ωn,H) 
    G0 = fit_ir(Gij,smpl_Matsubara,smpl_beta)            
    Δi = U*G0
    Δi = (1-mixratio)*Δold[i] + mixratio*Δi   
    return Δi
end

@assert typestable(calc_Δitsc_ir!,(Int64,Int64,AbstractMatrix{ComplexF64},Matrix{ComplexF64},Float64,
            Float64,Vector{ComplexF64},typeof(smpl),typeof(smpl_beta)))


function calc_Δtsc_ir!(Δnew,H,Δold,T,U,ωn,smpl_Matsubara,smpl_beta,;mixratio = 0.5)
    Nx,Ny = size(Δold)
    N = Nx*Ny*2
    map!(i -> calc_Δitsc_ir!(i,N,H,Δold,T,U,ωn,smpl_Matsubara,smpl_beta,mixratio = mixratio),Δnew,1:Nx*Ny) #If you use pmap! instead of map!, you can do the parallel computation.
end

@assert typestable(calc_Δtsc_ir!,(Matrix{ComplexF64},AbstractMatrix{ComplexF64},Matrix{ComplexF64},Float64,
            Float64,Vector{ComplexF64},typeof(smpl),typeof(smpl_beta)))

Then, we can solve the gap equations, self-consistently.

T = 0.01

beta = 1/T
wmax = 10.0

basis = FiniteTempBasis(Fermionic(), beta, wmax, 1e-5)
smpl = MatsubaraSampling(basis)
ωn = valueim.(smpl.sampling_points, beta)
println("num. of Matsubara freqs. ", length(ωn))
smpl_beta = TauSampling(basis; sampling_points=[beta])

U  =-5.6
itemax = 1000
μ = 3.5

Nx = 16
Ny = 16
Δ = 3*ones(ComplexF64,Nx,Ny)
Δold = copy(Δ)
Δnew = zero(Δ)
BC = "OBC"
h = 1
α = 1
Htsc =  make_Htsc_sc(Nx,Ny,μ,Δold,BC,h,α)

ix = Nx ÷ 2
iy = Ny ÷ 2
i = (iy-1)*Nx + ix
j = i


for ite = 1:itemax
    calc_Δtsc_ir!(Δnew,Htsc,Δold,T,U,ωn,smpl,smpl_beta,mixratio =1)
    update_Htsc_sc!(Htsc,Δnew)

    eps = sum(abs.(Δnew-Δold))/sum(abs.(Δold))
    println("$ite $eps ",Δnew[ix,iy]," ave: ",sum(abs.(Δnew))/length(Δnew))

    Δold .= Δnew
    if eps < 1e-2
        break
    end
end

num. of Matsubara freqs. 30
1 0.3728272031663908 

1.8802286062787428 - 2.0820223413844377e-14im ave: 1.881518390621739
2 0.209650742337318 

1.4861198287726276 - 4.1994906540563243e-14im ave: 1.487056674283118
3 0.13563820774942764 

1.2830170207309182 + 5.330995289932512e-11im ave: 1.285355106302483
4 0.0972523161393226 

1.1552868452410636 + 9.360232193763436e-11im ave: 1.160352015478455
5 0.07521015828299601 

1.0644744390740428 + 1.3505422426207966e-10im ave: 1.073083747489175
6 0.061407911146092134 

0.9944949237648685 + 3.187116676857662e-10im ave: 1.0071922114897494
7 0.05209959327112454 

0.9375697912624816 + 2.1840368945497596e-10im ave: 0.9547254445659928
8 0.045387786732905486 

0.8895959968259813 + 8.589252120964344e-10im ave: 0.9114040608700399
9 0.040235385603407746 

0.8482975309024207 + 6.149690686266491e-10im ave: 0.8747492182231559
10 0.03604857295301986 

0.8123534577627265 + 5.363451072088602e-10im ave: 0.8432360731227039
11 0.032483965845227475 

0.7809373993867752 + 1.3742306379620742e-9im ave: 0.8158690564865609
12 0.029346577429457656 

0.753468568363905 + 1.2249706269931202e-9im ave: 0.7919547202670784
13 0.026528633201899812 

0.7294847550366632 + 4.971410872983914e-10im ave: 0.7709774248320779
14 0.023971221111683563 

0.7085845368809598 + 3.994809620796865e-9im ave: 0.7525313875314211
15 0.02164101478045089 

0.6904049991209833 + 2.0900670310153102e-9im ave: 0.7362836093595521
16 0.019517048344906647 

0.6746151741872749 + 2.8943543563028365e-9im ave: 0.7219533098453342
17 0.017583787467650016 

0.6609149179468176 + 2.395281042460369e-9im ave: 0.7092999501530032
18 0.015827817555811634 

0.6490347990709716 + 1.8982173679248344e-9im ave: 0.6981156722522183
19 0.014236470218860176 

0.6387354469860635 + 1.6104222136288922e-9im ave: 0.6882200292928654
20 0.012797342842268314 

0.6298060229427608 + 1.5025395536408356e-9im ave: 0.6794560021896339
21 0.011498269632515516 

0.622062112513023 + 1.910165114608555e-9im ave: 0.671686771373637
22 0.010327447944612272 

0.6153432312811634 + 1.909253025728951e-9im ave: 0.6647929940233683
23 0.00927356128937015 

0.6095102715610458 + 1.2549056646787305e-9im ave: 0.6586704815714025

The superconducting order parameter is shown as

using Plots
plot(1:Nx,1:Ny,abs.(Δnew),st=:contourf,
    xlabel = "ix",
    ylabel = "iy",
    zlabel = "|Delta|")

We can see that the order parametere oscillates due to the Fridel oscillation.

The local density of states at the center is calculated as

using Plots
M = 1000
σ = zeros(ComplexF64,M)
η = 0.01
σmin = -4.0 + im*η
σmax = 4.0+ im*η
for i=1:M
    σ[i] = (i-1)*(σmax-σmin)/(M-1) + σmin
end


ix = Nx ÷ 2
iy = Ny ÷ 2
ispin  =1

i = (( iy-1)*Nx + ix-1)*2 + ispin
j = i
Gij1 = greensfunctions(i,j,σ,Htsc) 

ispin  =2
i = (( iy-1)*Nx + ix-1)*2 + ispin
j = i
Gij2 = greensfunctions(i,j,σ,Htsc) 
plot(real.(σ),(-1/π)*imag.(Gij1 .+ Gij2),
    xlabel = "Energy [t]",
    ylabel = "Local DOS",)

At the edge, there are mid-gap states since this is topologically non-trivial state. We can see the local density of states at the edge.

ix = 1
iy = Ny ÷ 2
ispin  =1

i = (( iy-1)*Nx + ix-1)*2 + ispin
j = i
Gij1 = greensfunctions(i,j,σ,Htsc) 

ispin  =2
i = (( iy-1)*Nx + ix-1)*2 + ispin
j = i
Gij2 = greensfunctions(i,j,σ,Htsc) 
plot(real.(σ),(-1/π)*imag.(Gij1 .+ Gij2),
    xlabel = "Energy [t]",
    ylabel = "Local DOS",)

We can see there are mid-gap states.