  There is growing evidence that the Galactic Center Excess identified in the
$\textit{Fermi}$-LAT gamma-ray data arises from a population of faint
astrophysical sources. We provide compelling supporting evidence by showing
that the morphology of the excess traces the stellar over-density of the
Galactic bulge. By adopting a template of the bulge stars obtained from a
triaxial 3D fit to the diffuse near-infrared emission, we show that it is
detected at high significance. The significance deteriorates when either the
position or the orientation of the template is artificially shifted, supporting
the correlation of the gamma-ray data with the Galactic bulge. In deriving
these results, we have used more sophisticated templates at low-latitudes for
the $\textit{Fermi}$ bubbles compared to previous work and the
three-dimensional Inverse Compton (IC) maps recently released by the ${\tt
GALPROP}$ team. Our results provide strong constraints on Millisecond Pulsar
(MSP) formation scenarios proposed to explain the excess. We find that an
$\textit{admixture formation}$ scenario, in which some of the relevant binaries
are $\textit{primordial}$ and the rest are formed $\textit{dynamically}$, is
preferred over a primordial-only formation scenario at $7.6\sigma$ confidence
level. Our detailed morphological analysis also disfavors models of the
disrupted globular clusters scenario that predict a spherically symmetric
distribution of MSPs in the Galactic bulge. For the first time, we report
evidence of a high energy tail in the nuclear bulge spectrum that could be the
result of IC emission from electrons and positrons injected by a population of
MSPs and star formation activity from the same site.
