We compute the galaxy power spectrum at two-loop order in cosmological perturbation theory (effective field theory, EFT). We derive galaxy bias operators through the fifth order and obtain two-loop renormalization conditions for their bias coefficients. We also derive the relevant higher-derivative and stochastic contributions, and implement IR resummation using . We obtain complete dark-matter and galaxy renormalization conditions for the higher-derivative operators at the power spectrum, bispectrum, and trispectrum level. Having identified the redundant operators, we find that the two-loop galaxy-galaxy and galaxy-matter power spectra require 28 additional bias, higher derivative, and stochastic parameters per galaxy sample relative to the one-loop model. We compare our computation with the galaxy-galaxy and galaxy-matter power spectra from the PT Challenge N-body simulation at $z=0.61$ and find a per mille-level agreement up to $k=0.85h$Mpc${}^{-1}$. We show that even with conservative priors on all EFT parameters, the two-loop model produces an unbiased measurement of the mass fluctuation amplitude ${\sigma }_8$ with three times narrower error-bars than the linear theory model. The improvement over the one-loop model is $\simeq 40\%$. This suggests significant gains in the two-loop EFT analyses of galaxy clustering and galaxy–lensing two-point functions (``$2\times 2$ pt’’) from CMB lensing maps and imaging surveys like Euclid, LSST, and Roman. In addition, our two-loop computation offers a probe of new physics scenarios that modify the shape of the matter power spectrum at wavenumbers $(0.4-0.8)h$Mpc${}^{-1}$ such as the presence of ultra-light axion dark matter sub-components with masses $m_a\sim {10}^{-24}$ eV.