Electroweak second order shifts of muonium ($\mu^+e^-$ bound state) energy
levels are calculated for the first time. Calculation starts from on-shell
one-loop elastic $\mu^+ e^-$ scattering amplitudes in the center of mass frame,
proceed to renormalization and to derivation of muonium matrix elements by
using the momentum space wave functions. This is a reliable method unlike the
unjustified four-Fermi approximation in the literature. Corrections of order
$\alpha G_F$ (with $\alpha \sim 1/137$ the fine structure constant and $G_F$
the Fermi constant) and of order $\alpha G_F /(m_Z a_B)$ (with $m_Z$ the Z
boson mass and $a_B$ the Bohr radius) are derived from three classes of Feynman
diagrams, Z self-energy, vertex and box diagrams. The ground state muonium
hyperfine splitting is given in terms of the only experimentally unknown
parameter, the smallest neutrino mass. It is however found that the neutrino
mass dependence is very weak, making its detection difficult.