Bethe-Slapter-type equation solver

Bethe-Slapter-type equation solver and the application to Cooper-pair linear response approach.

function BSeq_solver(param, G2::GreenFunc.MeshArray, kernel, kernel_ins, qgrids::Vector{CompositeGrid.Composite},
    Euv; Ntherm=120, rtol=1e-10, α=0.7, source::Union{Nothing,GreenFunc.MeshArray}=nothing,
    source_ins::GreenFunc.MeshArray=GreenFunc.MeshArray([1], G2.mesh[2]; dtype=Float64, data=ones(1, G2.mesh[2].size))
)

Bethe-Slapter equation solver by self-consistent iterations.

\[ R(\omega_n, k) = \eta(\omega_n, k) - \frac{1}{\beta} \sum_m \frac{d^dp}{(2\pi)^d} \Gamma(\omega_n,k;\omega_m,p) G^{(2)}(\omega_m,p)R(\omega_m,p) ,\]

where $\eta(\omega_n, k)$ is the sourced term, $\Gamma(k,\omega_n;p,\omega_m)$ is the particle-particle four-point vertex with zero incoming momentum and frequency, and $G^{(2)}(p,\omega_m)$ is the product of two single-particle Green's function.

Arguments

• param: parameters of ElectronGas.
• G2: product of two single-particle Green's function (::GreenFunc.MeshArray).
• kernel: dynamical kernel of the Legendre decomposed effective interaction.
• kernel_ins: instant part of the the Legendre decomposed effective interaction.
• qgrids: momentum grid of kernel (::Vector{CompositeGrid.Composite}).
• Euv: the UV energy scale of the spectral density.
• Ntherm: thermalization step. By defalut, Ntherm=120.
• rtol: tolerance absolute error. By defalut, rtol=1e-10.
• α: mixing parameter in the self-consistent iteration. By default, α=0.7.
• source: dynamical part of sourced term in imaginary-time space. By default, source=nothing.
• source_ins: instant part of sourced term. By default, source_ins=1.

Return

• Inverse of low-energy linear response 1/R₀ ($R_0=R(\omega_0, k_F)$)
• Bethe-Slapter amplitude F in imaginary-frequency space

\[ F(\omega_m, p) = G^{(2)}(\omega_m,p)R(\omega_m,p)\]

• dynamical part of R in imaginary-time space
• instant part of R in imaginary-time space

function G02wrapped(Euv, rtol, sgrid, param)

Returns the product of two bare single-particle Green's function.

\[ G_0^{(2)}(\omega_n, k) = 1/(\omega_n^2+\omega^2)\]

where $\omega= k^2/(2m)-\mu$.

function G2wrapped(Σ::GreenFunc.MeshArray, Σ_ins::GreenFunc.MeshArray, param)

Returns the product of two single-particle Green's function from given dynamical and instant parts of self-energy (Σ and Σ_ins).

\[ G_0^{(2)}(\omega_n, k) = 1/\left([\omega_n - \mathrm{Im} \Sigma(\omega_n, k)]^2+ [\omega+ \mathrm{Re} \Sigma(\omega_n, k) - \Sigma_{\mathrm{shift}}]^2\right)\]

where $\omega= k^2/(2m)-\mu$ and $\Sigma_{\mathrm{shift}}=\mathrm{Re} \Sigma(\omega_0, k_F)$.

function calcF!(F::GreenFunc.MeshArray, R::GreenFunc.MeshArray, R_ins::GreenFunc.MeshArray, G2::GreenFunc.MeshArray)

Calculation of the Bethe-Slapter amplitude F from the product of single-particle Green's function G2 and the dynamical and instant parts of R, R_ins. Compute in frequency space to avoid \tau integration.

\[ F(\omega_n, k) = G^{(2)}(\omega_n, k) [R(\omega_n, k)+R_{\mathrm{ins}}(k)]\]

function calcR!(F::GreenFunc.MeshArray, R::GreenFunc.MeshArray, R_ins::GreenFunc.MeshArray,
    source::Union{Nothing,GreenFunc.MeshArray}, kernel, kernel_ins, qgrids::Vector{CompositeGrid.Composite})

Calculate $kR(\tau, k)$ in three dimensions by given pF(\tau, p) and kernel. Compute in imaginary time space to aviod frequency convolution.

\[ kR(\tau, k) = k\eta(\tau, k) - \int \frac{dp}{4\pi^2} H(\tau,k,p) pF(\tau,p),\]

where $H(\tau,k,p)\equiv kp W(\tau,k,p)$ is the helper function of interaction (see Legendre Decomposition of Interaction) and the kernel argument. The dynamical source $k\eta(\tau, k)$ will be added if it is given as GreenFunc.MeshArray.

function calcR_2d!(F::GreenFunc.MeshArray, R::GreenFunc.MeshArray, R_ins::GreenFunc.MeshArray,  
    source::Union{Nothing,GreenFunc.MeshArray}, kernel, kernel_ins, qgrids::Vector{CompositeGrid.Composite})

Calculate $kR(\tau, k)$ in two dimensions by given pF(p) and kernel Compute in imaginary time space to aviod frequency convolution.

\[ kR(\tau, k) = k\eta(\tau, k) - \int \frac{pdp}{4\pi^2} W(\tau,k,p) pF(\tau,p),\]

where $W$ is the interaction and the kernel argument. The dynamical source $k\eta(\tau, k)$ will be added if it is given as GreenFunc.MeshArray.

function initFR(Euv, rtol, sgrid, param)

Initalize the Bethe-Slapter amplitude F in imaginary-frequency space and R ≡ F(GG)⁻¹ in imaginary-time space.

Arguments:

• Euv: the UV energy scale of the spectral density. parameter for DLR grids.
• rtol: tolerance absolute error. parameter for DLR grids.
• sgrid: momentum grid of F and R
• param: parameters of ElectronGas.

Return

• $F(\omega_n, k)=0$ as a GreenFunc.MeshArray
• the dynamical part of R in the imaginary-time space as a GreenFunc.MeshArray. In the imaginary-frequency space,

\[ R(\omega_n, k) = \frac{1}{\Omega_c^2+e^2}\left(1- \frac{2\omega_n^2}{\omega_n^2+\Omega_c^2} \right)\]

where $e=k^2/(2m)-\mu$ and $\Omega_c =0.01$.

• the instant part of R as a GreenFunc.MeshArray, $R_{\mathrm{ins}}(k)=0$.

function linearResponse(param, channel::Int; Euv=100 * param.EF, rtol=1e-10,
    maxK=10param.kF, minK=1e-7param.kF, Nk=8, order=8, sigmatype=:none, int_type=:rpa, α=0.7)

Implmentation of Cooper-pair linear response approach.

Arguments:

• param: parameters of ElectronGas.
• channel::Int: orbital angular channel (0: s-wave, 1: p-wave, ...)
• Euv: the UV energy scale of the spectral density. By default, Euv=100EF.
• rtol: tolerance absolute error. By defalut, rtol=1e-10.
• maxK: maximum momentum of kgrid and qgrids. By default, maxK=10kF.
• minK: minimum interval of panel kgrid and qgrids. By default, minK=1e-7kF.
• Nk: number of grid points of panel kgrid and qgrids. By defalut, Nk=8.
• order: number of grid points of subgrid of kgrid and qgrids. By defalut, order=8.
• sigmatype: type of fermionic self-energy. (no self-energy :none, G0W0 approximation :g0w0)
• int_type: type of effective interaction. By default, int_type=:rpa.
• α: mixing parameter in the self-consistent iteration. By default, α=0.7.

Return:

• Inverse of low-energy linear response 1/R₀ ($R_0=R(\omega_0, k_F)$)
• Linear response R(ωₙ, k), which is calculated by the Bethe-Slapter-type equation

\[ R(\omega_n, k) = 1 - \frac{1}{\beta} \sum_m \int \frac{d^dp}{(2\pi)^d} \Gamma(\omega_n,k;\omega_m,p) G^{(2)}(\omega_m,p)R(\omega_m,p)\]

where $1$ is the default sourced term, $\Gamma(k,\omega_n;p,\omega_m)$ is the particle-particle four-point vertex with zero incoming momentum and frequency, and $G^{(2)}(p,\omega_m)$ is the product of two single-particle Green's function.