Analysis of Fermi Large Area Telescope (Fermi-LAT) data has uncovered an extended gamma- ray source in the central region of the Milky Way. This Galactic Center Excess (GCE) has a spectral peak at a few GeV and appeared, initially, to have a spherically symmetric profile. These properties suggested that it may be evidence of self-annihilating weakly interacting massive particles (WIMPs) with a Navarro-Frenk-White (NFW) profile. However, the GCE gamma-ray spectrum is also similar to those of millisecond pulsars (MSPs) resolved in the Fermi-LAT data. An alternative possibility would then be that the GCE is produced by a population of MSPs in the Galactic Center too faint to be resolved individually as point sources at the present time. Providing further support for the MSP scenario, in recent years it has become apparent that the GCE may not be spherically symmetric, but may be spatially correlated with the distribution of stellar mass in the Galactic bulge.

In this thesis, we perform detailed modelling of the Galactic MSP population using data from the Fermi Large Area Telescope fourth source catalog data release 2 (4FGL-DR2) and the Australia Telescope National Facility (ATNF) pulsar catalog. Including in our model the spin down between formation and observation, we allow MSP luminosities to depend on intrinsic properties such as period P, magnetic field strength B and spectral energy cutoff Ecut. We find a model in which luminosity L ∝ Eaγ Bbγ E˙dγ provides the best fit to the data, where aγ = 1.2 ± 0.3, bγ = 0.1 ± 0.4 and dγ = 0.5 ± 0.1, and where E˙ ∝ B2/P4 is the spin-down power. This model is significantly better than one in which luminosity is independent of other properties of MSPs, with the luminosity distribution being fitted directly. The Milky Way disk is expected to be the source of the resolved MSPs, with the GCE potentially produced by MSPs in the Galactic bulge. The Galactic bulge is divided into two structures: the boxy bulge, a significant bar structure extending a few kpc from the Galactic Center; and the nuclear bulge, a less massive component with radius ∼ 230 pc. Due to differing star formation histories it is expected that the MSPs in the Galactic bulge are older and therefore dimmer than those in the Galactic disk. Additionally, correlations between the spectral parameters of the MSPs and the spin-down rate of the corresponding neutron stars have been observed. This implies that the bulge MSPs may be spectrally different from the disk MSPs. Although we confirm these correlations, we do not find they are sufficiently large to significantly differentiate the spectra of the bulge MSPs and disk MSPs when the uncertainties are accounted for. We find the age distributions of MSPs cannot be distinguished from a uniform birth rate, based on current data. Our results demonstrate that the population of MSPs that can explain the gamma-ray signal from the resolved MSPs in the Galactic disk and the unresolved MSPs in the boxy bulge and nuclear bulge can consistently be described as arising from a common evolutionary trajectory for some subset of astrophysical sources common to all these different environments. We do not require that there is anything systematically different about the inner Galaxy MSPs to explain the GCE. Additionally, we use a more accurate geometry for the distribution of bulge MSPs and incorporate dispersion measure estimates of the MSPs’ distances. We find that the elongated boxy bulge morphology means that some bulge MSPs are closer to us and so easier to resolve. We identify three resolved MSPs that may belong to the bulge population.

In the “recycling” channel of MSP formation the neutron star forms from a core collapse supernovae that undergoes a random “kick” due to the asymmetry of the explosion. This would imply a smoothing out of the spatial distribution of the MSPs. We use N-body simulations to model how the MSP spatial distribution changes. We estimate the probability distribution of natal kick velocities using the resolved gamma-ray MSP proper motions, where MSPs have velocities relative to circular motion of σv = 77 ± 6 km/s. The scale of these peculiar velocities are determined as part of our Galactic MSP population model. We find that, due to the natal kicks, there is an approximately 10% increase in each of the bulge MSP spatial distribution dimensions and also the bulge MSP distribution becomes less boxy. We estimate that natal kicks change the axis ratios of the MSP distribution in the Galactic boxy bulge from ∼ 1 : 0.43 : 0.40 to ∼ 1 : 0.46 : 0.44. Therefore, the bulge MSP distribution is still far from spherical.