LAUSR

Migdal-Eliashberg theory describes the properties of the normal and superconducting states of electron-phonon-mediated superconductors based on a perturbative treatment of the electron-phonon interactions. It is necessary to include both electron and phonon self-energies self-consistently in Migdal-Eliashberg theory in order to match numerically exact results from determinantal quantum Monte Carlo in the adiabatic limit. In this work we provide a method to obtain the real-axis solutions of the Migdal-Eliashberg equations with electron and phonon self-energies calculated self-consistently. Our method avoids the typical challenge of computing cumbersome singular integrals on the real axis and is numerically stable and exhibits fast convergence. Analyzing the resulting real-frequency spectra and self-energies of the two-dimensional Holstein model, we find that self-consistently including the lowest-order correction to the phonon self-energy significantly affects the solution of the Migdal-Eliashberg equations. The calculation captures the broadness of the spectral function, renormalization of the phonon dispersion, enhanced effective electron-phonon coupling strength, minimal increase in the electron effective mass, and the enhancement of superconductivity which manifests as a superconducting ground state despite strong competition with charge-density-wave order. We discuss surprising differences in two common definitions of the electron-phonon coupling strength derived from the electron mass and the density of states, quantities which are accessible through experiments such as angle-resolved photoemission spectroscopy and electron tunneling. An approximate upper bound on $2\mathrm{\ensuremath{\Delta}}/{T}_{c}$ for conventional superconductors mediated by retarded electron-phonon interactions is proposed.