Nevanlinna Analytic Continuation for Migdal-Eliashberg Theory
D. M. Khodachenko R. Lucrezi P. N. Ferreira M. Aichhorn C. Heil






Institute of Theoretical and Computational Physics, Graz University of Technology, NAWI Graz, 8010, Graz, Austria

Institute of Theoretical and Computational Physics, Graz University of Technology, NAWI Graz, 8010, Graz, Austria

Institute of Theoretical and Computational Physics, Graz University of Technology, NAWI Graz, 8010, Graz, Austria
Computational Materials Science Group (ComputEEL/MatSci), Universidade de São Paulo, Escola de Engenharia de Lorena, DEMAR, Lorena, Brazil

Institute of Theoretical and Computational Physics, Graz University of Technology, NAWI Graz, 8010, Graz, Austria

[Corresponding author: ]christoph.heil@tugraz.at
Institute of Theoretical and Computational Physics, Graz University of Technology, NAWI Graz, 8010, Graz, Austria

§ ABSTRACT 
In this work, we present a method to reconstruct real-frequency properties from analytically continued causal Green's functions within the framework of Migdal-Eliashberg (ME) theory for superconductivity. ME theory involves solving a set of coupled equations self-consistently in imaginary frequency space, but to obtain experimentally measurable properties like the spectral function and quasiparticle density of states, it is necessary to perform an analytic continuation to real frequency space. Traditionally, the ME Green's function is decomposed into three fundamental complex functions, which are analytically continued independently. However, these functions do not possess the causal properties of Green's functions, complicating or even preventing the application of standard methods such as Maximum Entropy. Our approach overcomes these challenges, enabling the use of various analytic continuation techniques that were previously impractical. We demonstrate the effectiveness of this method by combining it with Nevanlinna analytic continuation to achieve accurate real-frequency results for ME theory, which are directly comparable to experimental data, with applications highlighted for the superconductors MgB_2 and LaBeH_8.






Nevanlinna Analytic Continuation for Migdal-Eliashberg Theory
    C. Heil
    September 9, 2024
=============================================================





§ INTRODUCTION

Temperature-dependent quantum field theory via Green's functions is a state-of-the-art approach to theoretical condensed matter physics. Solving many-body Green's function equations directly on the real frequency axis is possible, but numerically very challenging, due to the existence of poles on the real axis, requiring the evaluation of principal value integrals. 
However, these equations simplify significantly when formulated in imaginary, i.e. Matsubara frequency space. 
By doing so, the problem becomes markedly simpler both analytically and numerically, as it leads to well-defined integrals and discrete Matsubara frequency sums <cit.>, reducing the numerical cost of calculations considerably compared to a direct evaluation in otherwise continuous real frequency space. 

The disadvantage that comes with this is that in order to be able to compare to real frequency properties obtained in experiment, such as angle-resolved photo emission spectra <cit.> or scanning tunneling spectroscopy <cit.>, an analytic continuation is required, which is an ill-conditioned problem.
As a result, 
a variety of analytic continuation methods have been developed, such as Nevanlinna analytic continuations (NAC) <cit.>, the Padé approximation <cit.>, Maximum Entropy (MaxEnt) formalism <cit.>, stochastic analytic continuation <cit.>, genetic algorithms, and machine learning <cit.> or causal projection methods <cit.>.

A highly topical and prominent use case for such a workflow is the state-of-the-art theoretical description of conventional superconductivity within Migdal-Eliashberg (ME) theory <cit.>, as for instance implemented in the EPW package <cit.> of the Quantum ESPRESSO software <cit.>, which allows for the ab initio computation of important properties of the superconducting phase, most notably the superconducting gap function Δ_n𝐤(iω_j), and thus , on the imaginary frequency axis. 

Many of the analytic continuation methods mentioned above require the causal Green's function in imaginary time or frequency. This is a problem in ME formalism, where the electron-phonon mediated self-energy Σ̂_n𝐤 and therefore the Green's function is split into three fundamental complex functions:
the already mentioned superconducting gap function Δ_n𝐤(iω_j), the mass renormalization function Z_n𝐤(iω_j) and the energy shift function χ_n𝐤(iω_j). It is important to have access to those in real frequency space, which is why current implementations of analytic continuation in the EPW package are restricted to methods that are able to continue these three functions directly. This is the case for the computationally lightweight Padé approximation <cit.> and a computationally very expensive iterative procedure <cit.> specific to ME theory. The latter is generally not feasible in most scenarios due to numerical cost, and the former is known to often have numerical artifacts which break the non-negative causality condition of spectral functions.

In this work, we will demonstrate how to reconstruct the complex functions Δ_n𝐤(ω), Z_n𝐤(ω) and χ_n𝐤(ω) in real frequency space by analytically continuing the components of the Nambu-Gor'kov Green's function, and we further introduce an efficient general analytic continuation workflow for ME theory.
Importantly, this workflow can be used with any analytic continuation method that requires causal Green's functions, therefore giving us access to many more methods without any loss of information for the complex functions Δ_n𝐤(ω), Z_n𝐤(ω) and χ_n𝐤(ω). Due to symmetry, only two analytic continuations are required, which is one less than previously needed for full-bandwidth calculations with χ_n𝐤(ω) ≠ 0 <cit.>. 

In this work, we employ the recently proposed NAC <cit.>, which is very well suited for continuing high-quality data and is able to reconstruct the spectral function with a high level of accuracy. We introduce a full implementation 
for both the isotropic and the anisotropic ME equations in the EPW package. We find that this method works exceptionally well for obtaining superconducting properties accurately, while keeping numerical costs comparable to the Padé approximation. NAC is numerically much more stable and, due to its analytic approach, always fulfills the causality conditions of normal spectral functions. 



§ THEORY




 §.§ Nevanlinna Analytic Continuation (NAC)


We start with a brief summary of the NAC as employed for our purposes. For a more detailed description of the method, we refer the reader to Refs. <cit.>. 

Our goal is to obtain the experimentally measurable spectral function A_n𝐤(ω), which is proportional to the imaginary part of the retarded Green's function G_n𝐤^ret,

    A_n𝐤(ω) = -1/π Im G_n𝐤^ret(ω+iη),

where n and 𝐤 are band and momentum indices, respectively. 
Due to the existence of poles on the real axis, the retarded Green's function is evaluated slightly above the real axis, described by the infinitesimal positive number η.
A causal spectral function such as the one in Eq. (<ref>) needs to fulfill two causality conditions,

    A_n𝐤(ω) ≥ 0     and  ∫_-∞^∞dω  A_n𝐤(ω) = 1.

We want to note at this point that NAC can still be used in cases where normalization is not fulfilled, as will be the case for the auxiliary Green's function introduced in section <ref>.

To obtain the retarded Green's function G_n𝐤^ret(ω+iη) from the imaginary frequency solution, analytic continuation in the form G_n𝐤(iω_j) → G_n𝐤^ret(ω+iη) is required. This problem, however, is ill-conditioned, i.e. small errors in the Matsubara Green's function G_n𝐤(iω_j) result in large errors in the retarded Green's function G_n𝐤^ret(ω+iη), regardless of the method one uses to analytically continue. The reason for that can be made apparent by considering
the kernel K(τ,ω) of the spectral representation of the Green's function,

    G_n𝐤(τ) = ∫dωe^-ωτ/1+e^-ωβ_K(τ,ω) A_n𝐤(ω),

which becomes exponentially small for high frequencies <cit.>. 
Determining the spectral function boils down to computing the inverse of the kernel, where for high frequencies divisions by very small numbers occur, resulting in small errors in the data having a huge impact on the result. Thus, the inverse of the kernel cannot be used to obtain the spectral function in practice.



The main interpolation procedure in NAC works according to the Schur algorithm <cit.>, performed on the unit circle. In particular, to construct Nevanlinna interpolants, the first step is to Möbius transform <cit.> the Green's function to a contractive function on the unit circle with 
θ(iω_j) = (-G_n𝐤(iω_j) - i) / (-G_n𝐤(iω_j) + i),
as can be appreciated in Fig. <ref>. We then use the following recursive formula from the Schur algorithm,

    θ(z)[z;θ_M(z)] = a(z)θ_M(z)+b(z)/c(z)θ_M(z)+d(z),

with

    [ a(z) b(z); c(z) d(z) ]
        = ∏_j = 1^M[                 z-iω_j/z+iω_j                   θ_j-1(iω_j); z-iω_j/z+iω_j [θ_j-1(iω_j)]^*                             1 ],

to obtain analytic continuation on the unit circle θ(z). 

By construction, θ_0(z)=θ(z) interpolates the Green's function for all Matsubara frequencies j=1,2,...,M, while each subsequent contractive function θ_1(z),...,θ_M(z) interpolates one less point than the previous one. As a result, the final function θ_M is an unconstrained arbitrary contractive function on the unit circle. For discrete spectral functions of Lorentzian shape, the method works well regardless of the choice of θ_M(z), which is, for example, the case for most isotropic ME solutions. 

The simplest approach is using a constant θ_M(z)=0, called the free solution.
However, in the case of smooth and rather featureless spectral functions, using a constant θ_M(z) can result in rapidly oscillating functions, 
in which case optimization of θ_M(z) can be required. We found that the free solution works very well in isotropic ME calculations, as discussed in section <ref>. However, in the case of anisotropic ME theory, individual spectral functions will usually oscillate for θ_M(z)=0. Optimization of θ_M(z) is possible <cit.>, but will not be considered in this work and is subject of future studies.

Finally, to obtain the desired G_n𝐤^ret(ω+iη), the function θ(z) in Eq. (<ref>) is evaluated at ω+iη (red dashed line) and then transformed back to the complex plane using the inverse Möbius transformation (see Fig. <ref>).

One big advantage of NAC is that we can use a generalization of the Pick criterion <cit.> to avoid the errors from higher frequency Matsubara points while also significantly improving numerical efficiency of the method. With the Möbius transformation h(z)=(z-i)/(z+i), the criterion is fulfilled if the matrix

    [ 1-θ(iω_j)θ(iω_j)^*/1-h(iω_i)h(iω_j^*)]_i,j   i,j = 1,2,...,M

is positive semi-definite, in our case verified by using Cholesky decomposition. When the Pick criterion is fulfilled, it is guaranteed that a solution to the Schur algorithm exists, which should always be the case for a noise-free imaginary solution. In practice, however, higher frequency points are much noisier, as has been discussed for the kernel in Eq. (<ref>), so they will generally not fulfill the Pick criterion. As a result, discarding frequencies that don't fulfill the criterion reduces noise and improves numerical speed of the method.



 §.§ Using NAC with Migdal-Eliashberg Theory



Contrary to a continued fraction approach like in the Padé approximation, NAC requires causal Green's functions to be employed. The Green's function we want to access in real frequency space using analytic continuation is the generalized Green's function Ĝ_n𝐤(τ) from the Nambu-Gor'kov formalism <cit.>,


    Ĝ_n𝐤(τ) = - 
        [    ⟨ T_τĉ_n𝐤↑(τ) ĉ_n𝐤↑^†(0) ⟩     ⟨ T_τĉ_n𝐤↑(τ) ĉ_n-𝐤↓(0) ⟩; ⟨ T_τĉ_n-𝐤↓^†(τ) ĉ_n𝐤↑^†(0) ⟩  ⟨ T_τĉ_n-𝐤↓^†(τ) ĉ_n-𝐤↓(0) ⟩ ],

where the main diagonal elements are the normal state Green's functions and the off-diagonal elements are the anomalous Green's functions, which describe Cooper pairs <cit.>.
Although the definition of the Nambu-Gor'kov matrix is different to a general matrix-valued Green's function, 
it is possible to transform one to the other using a spin-dependent particle-hole transformation such as ĉ_n-𝐤↓→ĉ_n-𝐤↓^†. As a result, analytic continuation methods for matrix-valued Green's functions do work for the Nambu-Gor'kov matrix. However, due to symmetry properties of both the main and off-diagonal elements of Ĝ_n𝐤, as well as its small size, an approach using auxiliary Green's functions as presented in this paper is much more efficient. As we will demonstrate later, only two analytic continuations are necessary to obtain all the components of the Nambu-Gor'kov matrix.

Due to the periodicity of the imaginary time Green's function, it can be expanded to imaginary Matsubara frequency space using a Fourier series of the form

    Ĝ_n𝐤(τ) = 1/β∑_n=-∞^∞e^-iω_jτĜ_n𝐤(iω_j),

and the Nambu-Gor'kov matrix in imaginary frequency space can be written as

    Ĝ_n𝐤(iω_j) =  
        [ G^11_n𝐤(iω_j)    F_n𝐤(iω_j);  F_n𝐤^*(iω_j) G^22_n𝐤(iω_j) ].


One can explicitly show that this is fulfilled for the negative of the main diagonal elements of the Nambu-Gor'kov matrix using the Lehmann representation

    G^11_n𝐤(iω_j) = 1/Z∑_mm'⟨ m'|ĉ_n𝐤↑|m⟩⟨ m|ĉ_n𝐤↑^†|m'⟩/iω_j + E_m' - E_m (e^-β E_m + e^-β E_m')

where Z is the partition function, and E_m and E_m' are the eigenenergies of the states |m⟩ and |m'⟩. Identifying the positive term K = Z^-1 |⟨ m|ĉ_n𝐤↑^†|m'⟩|^2 (e^-β E_m + e^-β E_m') ≥ 0 and using z=x+iy, with y>0, we find

    Im[-G_n𝐤(iω_j)] = ∑_mm'Ky/(x + E_m' - E_m)^2 + y^2≥ 0,

verifying the correct mapping 𝒩: ℂ^+ →ℂ^+. The prove for G^22_n𝐤(iω_j) follows the same lines. 

The situation for the anomalous Green's function F_n𝐤(iω_j) in the off-diagonal elements of Eq. (<ref>) however is different, as the product of the matrix elements ⟨ m'|ĉ_n𝐤↑|m⟩⟨ m|ĉ_n𝐤↑|m'⟩ cannot be rewritten as an absolute value. This means that the spectral function of the anomalous Green's function does not fulfill the causality condition of strictly non-negative spectral weight, 
and a direct analytic continuation using Nevanlinna functions is not possible.
This can also be seen from the anomalous spectral function A^an_n𝐤(ω) which is defined as usual from the imaginary part of F_n𝐤(ω). Due to the anti-commutation of the creation and annihilation operators with themselves, the normalization ∫_-∞^∞dω  A^an_n𝐤(ω) = 0 is fulfilled by definition.

Nevertheless, we can obtain the spectral functions of the anomalous Green's functions F_n𝐤(iω_j) by constructing auxiliary Green's functions that can be continued with NAC. 
In accordance with Refs. <cit.> we define a mixed operator of the form

    â_n𝐤(τ) = ĉ_n𝐤↑(τ) + ĉ_n-𝐤↓^†(τ).

The auxiliary Green's function G^aux_n𝐤(τ) = -⟨ T_τâ_n𝐤(τ) â_n𝐤^†(0) ⟩ is then guaranteed to have a spectral function with constant sign. Now, by rewriting the auxiliary Green's function in terms of the original creation and annihilation operators, one obtains

    G^aux_n𝐤(τ) =    - ⟨ T_τĉ_n𝐤↑(τ) ĉ_n𝐤↑^†(0) ⟩ - ⟨ T_τĉ_n𝐤↑(τ) ĉ_n-𝐤↓(0) ⟩
       - ⟨ T_τĉ_n-𝐤↓^†(τ) ĉ_n𝐤↑^†(0) ⟩ - ⟨ T_τĉ_n-𝐤↓^†(τ) ĉ_n-𝐤↓(0) ⟩,

or, in other words, the sum over all the components of the Nambu-Gor'kov matrix. In Matsubara frequency space, the auxiliary Green's function is then given by

    G^aux_n𝐤(iω_j) =     G^11_n𝐤(iω_j) + F_n𝐤(iω_j)   + 
       F^*_n𝐤(iω_j) + G^22_n𝐤(iω_j).


In ME theory, the components of the Nambu-Gor'kov matrix from Eq. (<ref>) have the form


    G^11_n𝐤(iω_j)    = -iω_j Z_n𝐤(iω_j) + ε_n𝐤 + χ_n𝐤(iω_j)/Θ_n𝐤(iω_j)
    
        F_n𝐤(iω_j)    = F^*_n𝐤(iω_j)=-Δ_n𝐤(iω_j) Z_n𝐤(iω_j)/Θ_n𝐤(iω_j)
    
        G^22_n𝐤(iω_j)    = -iω_j Z_n𝐤(iω_j) - ε_n𝐤 - χ_n𝐤(iω_j)/Θ_n𝐤(iω_j)

with

    Θ_n𝐤(iω_j) =    [ω_j Z_n𝐤(iω_j)]^2 + [ε_n𝐤 + χ_n𝐤(iω_j)]^2 + 
       [Δ_n𝐤(iω_j) Z_n𝐤(iω_j)]^2.

where we used the shorthand notation ε_n𝐤 = ϵ_n𝐤 - ϵ_F. Then the auxiliary Green's functions G^aux_n𝐤(iω_j) in Matsubara frequency space is given by

    G^aux_n𝐤(iω_j) = -2/Θ_n𝐤[ .    iω_j Z_n𝐤(iω_j) + Δ_n𝐤(iω_j) Z_n𝐤(iω_j) . ].



In contrast to the anomalous spectral function, the corresponding auxiliary spectral function A^aux_n𝐤(ω) is normalized to two, 
i.e. ∫_-∞^∞dω  A^aux_n𝐤(ω) = 2, and does fulfill the causality condition A^aux_n𝐤(ω) ≥ 0, allowing the application of NAC to G^aux_n𝐤(iω_j).

The symmetry G^11_n𝐤(-z)=-[G^22_n𝐤(z)]^*
is fulfilled in both real and imaginary frequency space, which we can use to avoid explicit analytic continuation of G^22_n𝐤. As a result, we only need two analytic continuations to obtain all the components of the Nambu-Gor'kov matrix in real frequency space, one for G^11_n𝐤 and one for G^aux_n𝐤. 
The analytically continued anomalous Green's functions F_n𝐤(ω) is then given by 


    F_n𝐤(ω) = 1/2[ . G^aux_n𝐤(ω) - G^11_n𝐤(ω) - G^22_n𝐤(ω) . ],

where we do not explicitly write the small imaginary part η to keep notations light.

We can now finally reconstruct all functions in ME formalism by using the Nambu-Gor'kov Green's function in real frequency space

    Ĝ_n𝐤(ω) = ω Z_n𝐤(ω)τ̂_0 + [ε_n𝐤 + χ_n𝐤(ω)]τ̂_3 + Δ_n𝐤(ω) Z_n𝐤(ω)τ̂_1/[ω Z_n𝐤(ω)]^2 - [ε_n𝐤 + χ_n𝐤(ω)]^2 - [Δ_n𝐤(ω)Z_n𝐤(ω)]^2.


Here, τ̂_j are the Pauli matrices. By solving Eq. (<ref>) like a nonlinear system of equations, we find for the gap function


    Δ_n𝐤(ω) = 2ω F_n𝐤(ω)/G^11_n𝐤(ω) + G^22_n𝐤(ω),

in agreement with Refs. <cit.>. For the remaining complex functions in ME formalism, we have the following two equations

    Z_n𝐤(ω) = - G^11_n𝐤(ω)+G^22_n𝐤(ω)/2ω[[F_n𝐤(ω)]^2-G^11_n𝐤(ω)G^22_n𝐤(ω)],


    χ_n𝐤(ω) = -2ε_n𝐤[F_n𝐤(ω)]^2+G^11_n𝐤(ω)[2ε_n𝐤G^22_n𝐤(ω)-1]+G^22_n𝐤(ω)/2[[F_n𝐤(ω)]^2-G^11_n𝐤(ω)G^22_n𝐤(ω)].

By obtaining the off-diagonal anomalous Green's function from the analytic continuation of the auxiliary Green's function in Eqs. (<ref>) and (<ref>), there is no longer any loss in information, regardless of the analytic continuation method. This is the key idea of the generalized analytic continuation workflow we will introduce in more detail in section <ref>.

Eqs. (<ref>-<ref>) represent the set of self-consistent equations for the fully anisotropic ME theory. To speed up the calculations, two approximations are usually employed. (i) The first one is the isotropic limit, in which the gap function and the pairing are averaged over the Brillouin zone. This can be justified by considering that the superconducting gaps are found to be isotropic for most materials (MgB_2 being a notable exception) and that, in real materials, defects are always present, which tend to average out anisotropic gaps. (ii) Secondly, the density of states is considered constant as only electronic states within an energy window of the phonon energy around the Fermi level will experience renormalization due to electron-phonon effects. Considering only assumption (i), commonly referred to as the isotropic full-bandwidth ME equations, our framework simplifies to 


    Δ(ω) = 2ω F(ω)/G^11(ω)+G^22(ω),


    Z(ω) = - G^11(ω)+G^22(ω)/2ω[[F(ω)]^2-G^11(ω)G^22(ω)],


    χ(ω) = G^22(ω)-G^11(ω)/2[[F(ω)]^2-G^11(ω)G^22(ω)].


Working with both assumptions, the energy shift function χ_n𝐤 vanishes in imaginary and real frequency space. As a result the main diagonal components o the Nambu-Gor'kov matrix become identical, i.e. G(ω)=G_n𝐤^11(ω)=G_n𝐤^22(ω). So, the ME equations in our framework are further simplified to


    Δ(ω) = ω F(ω)/G(ω),

and

    Z(ω) = - G(ω)/ω[[F(ω)]^2-[G(ω)]^2].


One particularly interesting quantity of the superconducting state, apart from Δ and  is the superconducting quasiparticle density of states, for the Fermi-surface restricted case given by

    N_S (ω)/N_F = -1/π∫_-∞^∞dϵ_n𝐤 Im  G^11_n𝐤(ω).


Although solving this integral numerically is possible from the G_ n𝐤^11 component, results will oscillate quite strongly beyond the initial gap for any method, as significantly finer 𝐤-grids would be necessary to average over them properly. 

In the isotropic BCS-limit, where it is assumed that χ_n𝐤=0 and Z_n𝐤=1 we have

    N_S (ω)/N_F = ∑_n ∫d𝐤/Ω_BZδ(ϵ_n𝐤 - ϵ_F)/N_FRe[ ω/√(ω^2 - Δ^2_n𝐤(ω))].

In other words, by using Eq. (<ref>), we also have access to the quasiparticle density of states in the BSC-limit, which generally results in more pronounced peaks and smoother curves.



 §.§ Analytic Continuation Workflow for Migdal-Eliashberg Theory using Green's functions


To perform the analytic continuation of the complex functions Δ_n𝐤, Z_n𝐤 and χ_n𝐤 from ME formalism using causal Green's functions, we introduce the following generalized workflow as shown in Fig. <ref>. The steps are as follows:



    
  * Start with the imaginary frequency solution from ME theory Δ_n𝐤(iω_j), Z_n𝐤(iω_j) and χ_n𝐤(iω_j).
    
  * Calculate the normal Green's function G^11_n𝐤(iω_j), as well as the auxiliary Green's function G^aux_n𝐤(iω_j) from Eq. (<ref>).
    
  * Perform the analytic continuation of these two Green's functions (a) to obtain them in real frequency (b).
    
  * Determine the remaining components of the Nambu-Gor'kov matrix in real frequency, by using the symmetry G^22_n𝐤(ω)=-[G^11_n𝐤(-ω)]^*, and Eq. (<ref>) for the anomalous Green's function F_n𝐤(ω).
    
  * Calculate any other real frequency properties from the obtained Nambu-Gor'kov matrix, as for instance the gap function Δ_n𝐤(ω) and the quasiparticle density of states in the BCS-limit N_S (ω)/N_F.


This workflow can be used with any analytic continuation method. Even in the case of the Padé approximation, which initially does not require causal Green's functions, it is possible to restrict the initial continuous fraction assumption to the analytic continuation of causal Green's functions <cit.>

    G_n𝐤(z) = a_1z+b_1-a_2z+b_2-a_2z+b_3-...,

where a_j > 0 and b_j∈ ℝ. 

Using this causal Green's function approach for the Padé approximation could make the method numerically more stable within ME theory, but is not the focus of this work.





§ RESULTS






We start by showing some results where NAC and the Padé approximation are in excellent agreement. First, for an isotropic full-bandwidth calculation of LaBeH_8, with a fine 𝐤-grid of 30 × 30 × 30 at a temperature of T=70 K. We used default parameters for the Padé approximation as implemented in EPW. For NAC, we used the Pick criterion to determine the optimal number of Matsubara frequencies and a broadening parameter of η=1 meV. The same broadening was added to the spectral function from the Padé approximation in post-processing to obtain Fig. <ref>, where we show a comparison between the two methods for both the normal and anomalous spectral functions A^11_n𝐤(ω) and A^12_n𝐤(ω), respectively. The pole position, and therefore the quasiparticle energy renormalized by superconducting pairing, is located around ω=± 35 meV. NAC is done without optimization, i.e. by using a constant free parameter of θ_M(z)=0. We found that this works very well for the isotropic approximation to ME theory, although small unphysical satellite peaks with a maximum of 1 eV^-1 to 2 eV^-1 can appear at higher frequencies. 

The most common outcome is this kind of agreement between NAC and the Padé approximation, which we found to be the case in about 95% of our testing.

NAC can also be used for anisotropic calculations, however spectral functions of individual n and 𝐤 index often have oscillations for θ_M(z)=0, meaning an optimization of the θ_M(z) parameter would be required. Work on this is ongoing and will be implemented in a future version. 



Nevertheless, we can use NAC to obtain excellent agreement in the quasiparticle density of states, demonstrated in Fig. <ref> for the well-studied two-gap superconductor MgB_2 <cit.>. Although slight oscillations persist for individual n and 𝐤 contributions when NAC is combined with anisotropic ME calculations, they effectively cancel out entirely in the quasiparticle density of states if the fine 𝐤-grid is sufficiently dense. As a result we have almost perfect agreement between the two methods.

From Fig. <ref>, we can, therefore, identify the gaps that originate from the different electronic π and σ orbitals. For example, in our calculation for T=15 K, we find a first gap of about 2Δ_π=5 meV and a second gap of about 2Δ_σ=16 meV. Finally, we obtained a critical temperature of = 42 K for our anisotropic calculation of MgB_2.



Beyond the most common 
cases with excellent agreement between NAC and the Padé approximation, we found two significant advantages of NAC: The first rather obvious one is the avoidance of non-causal negative A^11_n𝐤(ω) spectral functions as shown in Fig. <ref>, demonstrating a situation where the Padé approximation exhibits non-physical numerical artifacts at about ω=± 110 meV. 

By definition, the spectral functions from NAC is always positive, making the method preferable for analytic continuations of simple Lorentzian peak shaped spectral functions, which is almost always the case for isotropic ME calculations.

The second advantage is that NAC is able to analytically continue even in difficult cases where the Padé approximation fails entirely. Most common for isotropic calculations at very low temperatures, the computation of the Padé interpolants will run into division by zero issues. 
The reason for this is that at low temperatures, the ME calculation requires a lot of Matsubara points, which then enter the Padé approximation and significantly increase the probability of failure. 

This can be appreciated in Fig. <ref> for an isotropic calculation of MgB_2 and LaBeH_8 at a low temperature of T=0.05 K. NAC is numerically much more stable by introducing quad precision arrays in the Schur algorithms of the method. We tested a similar quad precision approach for the Padé approximation; however, the numerical cost becomes much larger by approximately a factor of ten times more than NAC. The reason for this is the efficiency of the Pick criterion, which allows NAC to require much less Matsubara frequencies. However, even when using the same number of Matsubara frequencies for both methods, i.e. by not using the Pick criterion, NAC has never run into division by zero issues for our testing, as numerical errors only appear in the form of the above-mentioned satellite peaks.




The numerical costs of NAC using quadruple precision are very similar to those of the Padé approximation using double precision. We reached this conclusion by comparing the runtimes of anisotropic calculations, as isotropic calculations are too brief to provide meaningful comparisons. However, it's important to note that the runtimes of anisotropic calculations depend highly on the chosen 𝐤-grid, so the following estimates are intended to give a general sense of the performance. In our calculations for MgB_2, we observed typical CPU runtimes of approximately 50 seconds per temperature for NAC with quadruple precision, comparable to the runtimes for Padé approximation calculations using double precision. In contrast, using quadruple precision for the Padé approximation significantly increases the computation time to around 500 seconds.



§ CONCLUSION

In this work, we demonstrate a method to obtain
real-frequency properties from the analytically continued Nambu-Gor'kov matrix within the framework of Migdal-Eliashberg (ME) theory for superconductivity.
By leveraging the symmetry properties of the Green's function, our method requires only two analytic continuations to obtain all components in real frequency space: one for the normal and one for the auxiliary Green's function. This significantly broadens the scope for applying a variety of analytic continuation techniques that were previously inaccessible in the Migdal-Eliashberg framework, retaining full information of the real frequency complex functions explicitely, such as Δ_n𝐤(ω), Z_n𝐤(ω), and χ_n𝐤(ω). Our generalized workflow opens up the possibility for employing advanced methods like MaxEnt, stochastic analytic continuation, genetic algorithms, and causal projection in the modeling of superconducting materials, thereby enhancing the applicability and accuracy of the Migdal-Eliashberg formalism.

In particular, as a very powerful method for analytic continuation, we implemented the NAC into this framework.

Compared to the standard Padé approximation, our method offers two key advantages: (i) it effectively avoids non-physical negative values in the spectral function, and (ii) it provides improved numerical stability, enabling accurate analytic continuations at significantly lower temperatures than previously achievable. Despite these enhancements, the NAC method maintains comparable CPU runtimes, thanks to the efficiency of the Pick criterion and the streamlined analytic continuation workflow we introduced.

Using NAC with an unoptimized θ_M(z) parameter performs well for isotropic Migdal-Eliashberg calculations at lower frequencies. However, minor satellite peaks can emerge at higher frequencies due to random oscillations, indicating the need for careful optimization of the θ_M(z) parameter. Likewise, optimization can be necessary to mitigate oscillatory behavior and ensure accurate results when dealing with spectral functions from anisotropic calculations.



§ DATA AVAILABILITY


The authors confirm that the data supporting the findings of this study are available within this article. The updated and newly developed code routines can be obtained upon request and will also be made available within a future EPW release.

We thank E. Gull, E. R. Margine, and H. Mori for their insightful discussions. Calculations have been performed on the dCluster and the lCluster of the Graz University of Technology and on the Vienna Scientific Cluster (VSC) under Project 71754. PNF acknowledges the São Paulo Research Foundation (FAPESP) under Grants 2020/08258-0 and 2021/13441-1.





§ COMPUTATIONAL DETAILS

DFT and DFPT calculations were performed using the Quantum ESPRESSO (QE) software, employing norm-conserving Vanderbilt (ONCV) pseudopotentials <cit.> with a PBE-GGA exchange-correlation potential.

For Fm3m-LaBeH_8 at 100 GPa, we used a 12×12×12 Monkhorst-Pack 𝐤-grid, a 6×6×6 Monkhorst-Pack 𝐪-grid, a kinetic energy cutoff of 80 Ry, a smearing of 0.01 Ry, an electronic convergence threshold of 10^-10 Ry, and a phonon self-consistency threshold of 10^-14 Ry. 

After that, we used the Wannier implementation of EPW to interpolate electron-phonon matrix elements onto dens 𝐤- and 𝐪-grids. We used the "selected columns of the density matrix" (SCDM) method <cit.> to automatically wannierize our electronic band structure, using a frozen window of about [-14 eV, +5 eV] around the Fermi energy E_F. For the SCDM entanglement, we use an error function with the two parameters σ=5 and μ_c at the Fermi energy E_F. Then for the solution of the isotropic full-bandwidth ME equations, we used a coarse 𝐤- and 𝐪-grids of 6×6×6 interpolated onto 30×30×30 fine grids. We set the lower boundary of phonon frequency to 10 cm^-1, the width of the Fermi surface window to 1 eV, and the imaginary convergence threshold to 10^-5 Ry, all for a Morel-Anderson pseudopotential of μ^*=0.16. The ME equations were solved while updating the chemical potential to keep the particle number constant, which requires a comparably high-frequency cut-off of 10 eV for convergence. We used sparse sampling and analytically continued with the Padé approximation and NAC simultaneously, starting at a temperature of 60 K increased by steps of 10 K.

For MgB_2, we used 80 Ry for the kinetic energy cutoff, 24^3 Monkhorst-Pack 𝐤-grid, and a Methfessel-Paxton first-order spreading <cit.> of 0.02 Ry for Brillouin zone integration. The dynamical matrices and the variations of the self-consistent potentials were obtained on a 6^3 Monkhorst-Pack 𝐪-grid with a phonon self-consistency threshold of 10^-16 Ry.

The electronic wave functions for the Wannier interpolation were calculated on a Γ-centered uniform 𝐤-grid of 6^3 points. We used 5 Wannier functions for describing the band structure of MgB_2 around the Fermi level: two B-p_z-like states and 3 sp^2 hybridized orbitals associated with the Mg atom.

To solve the Eliashberg equations for MgB_2, we adopted the Fermi-surface restricted approximation <cit.> to keep our calculations consistent with the approach adopted in Ref. <cit.>. We interpolated the electron-phonon matrix elements onto fine Γ-centered uniform grids containing 60^3-𝐤 and 30^3-𝐪 points. The Matsubara frequency cutoff is set to 1 eV, approximately 10 times the maximum phonon frequency. The smearing parameters in the Dirac δ functions for electrons and phonons are set to 100 meV and 0.5 meV, respectively. The Pick criterion was used to determine the optimal number of Matsubara points, and the real frequency broadening of NAC was set to zero.

For the Coulomb pseudopotential μ^* of MgB_2, we utilized the μ^* calculated in Ref. <cit.> according to the Morel-Anderson approximation at the GW level. To rescale the value of μ^* used for solving the Allen-Dynes equation in Ref. <cit.> to be used in the Eliashberg equation in the present work, we used the simple rule <cit.>

    1μ^*_Eliashberg = 1μ^*_AD + ln(ω_phω_c),

where ω_ph is the maximum of the phonon spectrum. By doing this, we obtained a μ^* of 0.11.





§ EPW INPUT FLAGS

We have implemented three new input flags that are all written into the input file of the EPW calculation. Although not necessary, it is good practice to perform analytic continuations using a restart calculation. We require the imaginary ME solution, so "limag=.true." has to be set. Then, to enable the NAC calculation, use the logical flag "lneva=.true." (default .false.), which can also be done in the same run as any other analytic continuation method. By default, the Pick criterion will be used to determine the optimal number of Matsubara points. However, if that is not desirable one can disable the criterion with "lpick=.false." (default .true.), where all Matsubara frequencies are used instead. Finally, we have implemented a flag for the real frequency broadening η (in eV) of NAC, called "neva_broad". By default, a small broadening of "neva_broad=0.001", or 1 meV, is used.


apsrev4-2


