Numerical implementation of the Eliashberg equations on the imaginary axis at strong coupling

In the presence of electron-phonon coupling, phonon dynamics within a superconductor give rise to an electronic retardation effect, which in turn results in a frequency-dependent order parameter. In Eliashberg theory, the frequency-dependent pairing and renormalization functions are determined self-consistently by solving the Eliashberg equations, which take into account the dynamical phonon-mediated interaction. We review and extend the numerical implementation of the finite-temperature Eliashberg equations on the imaginary axis, highlighting the technical difficulties in calculating the pairing function away from the weak-coupling regime. We similarly review the numerical calculation of physical quantities of interest, such as the critical temperature, and make comparisons with applicable results in the literature.