We calculate ionization energies and fundamental vibrational transitions for H + 2 , D + 2 , and HD + molecular ions. The NRQED expansion for the energy in terms of the fine structure constant α is used. Previous calculations of orders mα 6 and mα 7 are improved by including second-order contributions due to the vibrational motion of nuclei. Furthermore, we evaluate the largest corrections at the order mα 8 . That allows to reduce the fractional uncertainty to the level of 7.6 × 10 −12 for fundamental transitions and to 4.5 × 10 −12 for the ionization energies.The hydrogen molecular ions (HMI) play an essential role in testing molecular quantum mechanics [1,2]. From the theoretical point of view the HMI is one of the simplest nonintegrable quantum system, which still allows very accurate numerical treatment. As was pointed out already some time ago [3], and recently discussed more extensively [4], if the theory would be sufficiently precise, the spectroscopy of HMI may be used for determining of the fundamental physical constants such as the protonto-electron mass ratio. The ionization energy of HMI is also of high importance for the determination of ionization and dissociation energies of the hydrogen molecule from spectroscopic studies [5][6][7] as well as for the determination of atomic masses of light nuclei [8][9][10].On the experimental side there are many new projects started, which are now oriented towards Doppler-free spectroscopy with accuracy targeted to 1 ppt (one part per trillion) or better [4,[11][12][13]. These perspectives bring strong motivation for theory.The aim of this Letter is to improve the theoretical precision of spin-averaged energies and ro-vibrational transition frequencies in HMI. To this end we consider the largest QED contributions which had not been evaluated in our previous works [14][15][16], namely, corrections to orders mα 6 and mα 7 due to vibrational motion of nuclei and the leading contributions to order mα 8 . As it was shown recently [17], taking into account the vibrational motion of nuclei is essential for accurate theoretical description. It has allowed to resolve the longstanding discrepancy between theory and experiment in the hyperfine structure of H + 2 ion. These new achievements reduce the relative uncertainty in the fundamental vibrational transitions of HMI to the level of 7.6 × 10 −12 and allow to obtain the most precise theoretical values for the ionization energies of H + 2 , D + 2 , and HD + molecular ions. In conclusion we discuss how these new results may have impact on fundamental physical constants such as the Rydberg constant, proton-to-electron mass ratio, and proton charge radius.We use atomic units throughout this paper.The terms of mα 6 and higher orders are calculated in the adiabatic approximation. For this purpose we use the Born-Oppenheimer formalism. In this approach the states of the molecule are taken in the formThe electronic wave function obeys the clamped nuclei Schrödinger equation for a bound electronwhere H el = p 2 /(2m) + V + Z 1 Z ...