We calculate the two-loop Bethe logarithm correction to atomic energy levels in hydrogen-like systems. The two-loop Bethe logarithm is a low-energy quantum electrodynamic (QED) effect involving multiple summations over virtual excited atomic states. Although much smaller in absolute magnitude than the well-known one-loop Bethe logarithm, the two-loop analog is quite significant when compared to the current experimental accuracy of the 1S-2S transition: it contributes -8.19 and -0.84 kHz for the 1S and the 2S state, respectively. The two-loop Bethe logarithm has been the largest unknown correction to the hydrogen Lamb shift to date. Together with the ongoing measurement of the proton charge radius at the Paul Scherrer Institute its calculation will bring theoretical and experimental accuracy for the Lamb shift in atomic hydrogen to the level of 10 −7 .PACS numbers: 12.20. Ds, 31.30.Jv, 06.20.Jr, In 1947 Hans Bethe explained the splitting of 2S 1/2 and 2P 1/2 levels in hydrogen by the presence of the electron self-interaction [1], and expressed it in terms of the "Bethe" logarithm. For S states this quantity may be represented as a matrix element involving the logarithm of the nonrelativistic Hamiltonian of the hydrogen atom.In natural units withh = c = ǫ 0 = 1 and m denoting the electron mass, it readsThis Bethe logarithm is due to the emission and subsequent absorption of a single soft virtual photon (it is independent of the nuclear charge number Z and depends only on the principal quantum number n and the orbital angular momentum which is zero for S states). Over the years, QED theory has been developed and refined [2], and various additional radiative, relativistic, and combined corrections have been obtained to face the increasing precision of the measurements of the hydrogen spectrum [3,4]. These include higher-order relativistic one-, two-, and three-loop corrections, nuclear recoil, finitesize corrections, and even the nuclear polarizability. The modern all-order calculation of the leading one-loop selfenergy was developed by Mohr in [5] and significantly improved recently using convergence acceleration techniques which led to a highly accurate evaluation of the fully relativistic Green function [6]. One of the conceptually most difficult and as well as interesting corrections involve nuclear recoil effects. The finite nuclear mass, although large as compared to the electron mass, prohibits the use of the one-body Dirac equation, and alternative approaches such as the Bethe-Salpeter equation or nonrelativistic QED [7] have been introduced. Although these methods are quite general, no compact formulas have been derived for relativistic recoil effects. A few years ago, Shabaev tackled the problem of recoil corrections to hydrogenic energy levels of first order in the mass ratio, deriving expressions which are nonperturbative in the nuclear charge (see a recent review in [8]), and this has led to the current highly accurate calculations of relativistic recoil corrections. Another class of effects, namely bin...