To obtain a state-of-the-art benchmark potential energy surface ͑PES͒ for the archetypal oxidative addition of the methane C-H bond to the palladium atom, we have explored this PES using a hierarchical series of ab initio methods "Hartree-Fock, second-order Møller-Plesset perturbation theory, fourth-order Møller-Plesset perturbation theory with single, double and quadruple excitations, coupled cluster theory with single and double excitations ͑CCSD͒, and with triple excitations treated perturbatively ͓CCSD͑T͔͒… and hybrid density functional theory using the B3LYP functional, in combination with a hierarchical series of ten Gaussian-type basis sets, up to g polarization. Relativistic effects are taken into account either through a relativistic effective core potential for palladium or through a full four-component all-electron approach. Counterpoise corrected relative energies of stationary points are converged to within 0.1-0.2 kcal/mol as a function of the basis-set size. Our best estimate of kinetic and thermodynamic parameters is Ϫ8.1 ͑Ϫ8.3͒ kcal/mol for the formation of the reactant complex, 5.8 ͑3.1͒ kcal/mol for the activation energy relative to the separate reactants, and 0.8 ͑Ϫ1.2͒ kcal/mol for the reaction energy ͑zero-point vibrational energy-corrected values in parentheses͒. This agrees well with available experimental data. Our work highlights the importance of sufficient higher angular momentum polarization functions, f and g, for correctly describing metal-d-electron correlation and, thus, for obtaining reliable relative energies. We show that standard basis sets, such as LANL2DZϩ1 f for palladium, are not sufficiently polarized for this purpose and lead to erroneous CCSD͑T͒ results. B3LYP is associated with smaller basis set superposition errors and shows faster convergence with basis-set size but yields relative energies ͑in particular, a reaction barrier͒ that are ca. 3.5 kcal/mol higher than the corresponding CCSD͑T͒ values.