# The effect of kick velocities on the spatial distribution of millisecond
pulsars and implications for the Galactic center excess

Harrison Ploeg and Chris Gordon

###### Abstract

Recently it has become apparent that the Galactic center excess (GCE) is
spatially correlated with the stellar distribution in the Galactic bulge. This
has given extra motivation for the unresolved population of millisecond
pulsars (MSPs) explanation for the GCE. However, in the “recycling” channel
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 random velocities relative to the circular motion
with a scale parameter of $77\pm 6$ km/s. 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
but is still far from being spherical.

##  1 Introduction

The Galactic Center Excess (GCE) is an extended source of gamma radiation in
the central region of the Galaxy found in the Fermi Large Area Telescope
(Fermi-LAT) data. When first discovered, its apparently spherically symmetric
profile and spectrum peaking at a few GeV suggested that it may be evidence of
self-annihilating weakly interacting massive particles (WIMPs) distributed
according to a Navarro-Frenk-White (NFW) profile [1, 2, 3, 4]. An alternative
scenario that was proposed was one in which the GCE was produced by a
population of unresolved millisecond pulsars (MSPs) [5]. More recently,
evidence has suggested that in fact the GCE is not spherically symmetric but
actually is correlated with the stellar mass in the Galactic bulge and so has
a more boxy morphology [6, 7, 8, 9, 10, 11]. Although one recent study, using
different methods, still argues for a spherically symmetric GCE [12]. If the
GCE does trace the stellar mass of the bulge this would favor the MSP (or
other dim unresolved astrophysical point source) explanation. Additionally,
several studies have claimed to find a non-Poissonian component to the GCE
[13, 14] which may be further evidence for the MSP explanation. However, there
is some controversy regarding the level of systematics in this approach [e.g.,
15, 16, 17, 18, 19, 11].

In the “recycling” model of MSP formation a neutron star is spun up to
millisecond periods through the transfer of mass from a binary companion. This
requires that the binary system survives the kick produced by any asymmetry in
the core collapse supernova explosion [20]. However, an alternative to the
recycling channel is accretion induced collapse of white dwarfs into neutron
stars which may produce more than half of all observed MSPs [21, 22, 23, 24].
In this case the system does not receive a significant natal kick [25, 26].
This would imply that the MSPs have much smaller peculiar velocities in
comparison to the recycling model [27, 28, 29].

Ploeg et al. [30] (hereafter referred to as P20) modelled Fermi-LAT detected
MSPs as having a Maxwell distributed peculiar velocity with the scale
parameter $\sigma_{v}$ found to be $77\pm 6$ km s-1 where we quote error bars
at the 68% confidence interval throughout this article. This velocity applies
for disk MSPs and we assume it will not be significantly different for bulge
MSPs. Although the star formation histories are very different in the bulge
and disk, as can be seen from Fig. 6 of P20 the probability distribution of
luminosities in the bulge and disk only differ by a few percent. Also, as can
be seen from Fig. 11 of P20, the bulge and disk have a ratio of number of MSPs
formed per solar mass which is within one error bar of each other. Based on
this, we assume that the bulge and disk have the same mix of MSP formation
channels and thus the same probability distribution of natal kick velocities.

If the source of the GCE is a population of unresolved MSPs, then the spatial
distribution may be smoothed to some degree relative to the stellar mass in
the bulge. Eckner et al. [31] used the virial theorem to estimate the
“smoothing length” of MSPs as $700$ – $900$ pc for kicks $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 3.0pt\hbox{$\sim$}}\hss}\hbox{\raise
2.0pt\hbox{$<$}}}}70$ km s-1. However, they assumed a spherically symmetric
spatial distribution for the MSPs.

In this article we use $N$-body simulations to estimate what are the effects
of MSP kicks for a boxy bulge like distribution. In Section 2 we explain our
method. Our results are given in Section 3, and our conclusions in Section 4.

##  2 Method

For this work we use the code of Bedorf et al. [32] to run $N$-body
simulations in order to model the Milky Way.111Available at:
https://github.com/treecode/Bonsai We use parameters corresponding to initial
conditions labelled MWa, MWb and MWc0.8 as denoted by Fujii et al. [33] as
they were the initial conditions that led to the best fitting simulations to
Galaxy observations that Fujii et al. found. The initial conditions consist of
a parametric description of the dark matter halo, Galactic disk and bulge.
Comparing to bulge kinematics, bar length, and pattern speed observations
Fujii et al. found MWa was the initial conditions which gave the best fit to
the data, followed by MWb, and then MWc0.8. For each initial condition model
we generated a total of $30$ million disk, bulge and dark matter halo
particles. These initial populations are generated using the methods of
Kuijken and Dubinski [34], Widrow and Dubinski [35], and Widrow et al.
[36].222We used the implementation at:
https://github.com/treecode/Galactics.parallel As in Fujii et al. [33], we use
time-steps of $\sim 0.6$ Myr, an opening angle of $0.4$ and ran the simulation
for $10$ Gyr. However, we use a softening length of $30$ pc. Also, our dark-
matter halo particles have a mass $8$ times larger than the disk and bulge
particles. Taking into account the masses of the various components, this
implies that, out of the 30 million particles, approximately $10$ million
particles represent stellar mass and the remainder represent dark matter.

In order to model the density of MSPs, we additionally include massless (so
they do not affect the simulation) disk and bulge particles that are given a
normally distributed perturbation to each component of their velocity vector
with mean zero and standard deviation $\sigma_{k}$. The kick velocity
magnitude is therefore Maxwell distributed. The probability density function
of a Maxwell distribution can be written as:

$p(x)=\sqrt{\frac{2}{\pi}}\frac{x^{2}\exp\left(-x^{2}/2\sigma^{2}\right)}{\sigma^{3}}$
(2.1)

where $x$ is the magnitude of a three dimensional vector with components
sampled from the normal distribution $\mathcal{N}(0,\sigma^{2})$. For each
model we try a case where the kicks occurred at the beginning of the $N$-body
simulations and a case where the kicks occur randomly with a uniform rate over
the course of the $10$ Gyr. These cases approximate scenarios in which MSPs
are largely born early or at a relatively steady rate. From these two
extremes, we can estimate the sensitivity of our results to the MSP age
distribution.

The first step is to estimate the kick velocity scale required to produce a
peculiar velocity distribution consistent with P20 where for the best model
the peculiar velocity scale parameter was $\sigma_{v}=77\pm 6$ km s-1. We do
this by running each initial condition model with $41$ populations of $10^{5}$
kicked particles with $\sigma_{k}$ between $70$ and $110$ km s-1. We separate
the velocity of each particle into two components:

$\bm{v}=\bm{v}_{c}+\bm{v}_{p}$ (2.2)

where $\bm{v}_{c}$ is the velocity of a particle on a circular orbit around
the center of the galaxy and $\bm{v}_{p}$ is the peculiar velocity. The
magnitude of $\bm{v}_{c}$ for a particle with coordinates $x$, $y$, and $z$
can be evaluated using the centripetal force:

$\lVert\bm{v}_{c}\rVert=\sqrt{\lVert\bm{a}_{c}(x,y,0)\rVert R}$ (2.3)

where $\bm{a}_{c}(x,y,z)$ is the acceleration toward the center of the galaxy
which can be obtained from a small modification to the the $N$-body code and
$R^{2}=x^{2}+y^{2}$. For $R$ outside the bulge region, and for small peculiar
velocity, we are therefore assuming that particles are rotating with the disk,
with $\bm{v}_{c}$ the rotation velocity of the disk at $R$ [37].

We use a maximum likelihood estimate of the final $\sigma_{v}$ for each
initial $\sigma_{k}$. For a set of $N$ particles with peculiar velocities
$v_{1}$, …, $v_{N}$, the log-likelihood is obtained by assuming velocities
have a Maxwell distribution:

$\log(L)=\frac{N}{2}\log(\frac{2}{\pi})-3N\log(\sigma_{v})+\sum_{i=1}^{N}2\log(v_{i})-\frac{v_{i}^{2}}{2\sigma_{v}^{2}}$
(2.4)

and therefore

$\frac{\differential\log(L)}{\differential\sigma_{v}}=-\frac{3N}{\sigma_{v}}+\sum_{i=1}^{N}\frac{v_{i}^{2}}{\sigma_{v}^{3}}\,.$
(2.5)

Then solving for $\sigma_{v}$ where
$\differential\log(L)/\differential\sigma_{v}=0$, we find the maximum
likelihood estimate for $\sigma_{v}$ is:

$\hat{\sigma_{v}}=\sqrt{\frac{\sum_{i=1}^{N}v_{i}^{2}}{3N}}$ (2.6)

This is done for particles where $4\textrm{ kpc}\leq R\leq 12\textrm{ kpc}$
and $\lvert z\rvert\leq 2\textrm{ kpc}$, ensuring we are estimating the
peculiar velocity distribution scale parameter for particles in the disk
region from which gamma-ray MSPs are most likely to be resolved and where
$\bm{v}_{c}$ approximately represents disk rotation.

Once we have a best fitting $\sigma_{k}$ we use that $\sigma_{k}$ to generate
an additional population of kicked particles. We rerun the $N$-body
simulations for each set of initial conditions (MWa, MWb, MWc0.8) including in
each initial condition case two new populations of massless particles where
each new population has $2\times 10^{6}$ particles. These new populations are
made massless so as not to influence the gravitational potential of the
$N$-body simulation. The first massless populations is sampled from the
initial condition model. For the second population we duplicate the first
population and then add a Maxwell distributed kick with scale parameter
$\sigma_{k}$ at the beginning of the simulation to each particle in the second
population. The simulations are then run for 10 Gy to see what the final state
of the two massless populations is for each set of initial conditions. We then
repeat this process but instead of the second population having a Maxwell
distributed kick added at the beginning of the simulation we add it randomly
to members of the second population at a uniform rate throughout the
simulation. So in the end for each initial condition we have two massless
populations of unkicked particles, one massless population of particles kicked
at the beginning and one for particles kicked at a uniform rate throughout the
simulation. We used two populations of massless unkicked particles to make the
kicked at beginning cases and uniform kick rate cases as independent as
possible from each other.

In order to understand the effect of MSP kicks on the final spatial
distribution of MSPs we found it clearer to have a parametric description of
the spatial distribution. Then by comparing the parameters for the kick and
no-kick cases we can have a succinct way of characterising the effects of the
kicks. We use Markov Chain Monte Carlo (MCMC) to fit the parametric model to
both final particle distributions. The form of our parametric model for the
final spatial distribution of MSPs was chosen based on previous studies and
also such that there were not excessive residuals between the model fits and
the simulations. Our model of the final distributions consists of four
components: a spherically symmetric bulge, a potentially non-spherical bulge
also known as a bar, a long bar and a disk. The spherically symmetric bulge
component uses the Hernquist model [38]:

$\rho_{\rm
Hernquist}(r)\propto\frac{1}{\left(r/a_{b}\right)\left(1+r/a_{b}\right)^{3}}$
(2.7)

where $r^{2}=x^{2}+y^{2}+z^{2}$ and $a_{b}$ is a free parameter. Our initial
conditions models, based on those of Fujii et al. [33], also include a
component distributed according the Hernquist model which would naturally have
different parameter values to the one we fitted after 10 Gy of evolution. Our
bar model is distributed as:

$\rho_{\rm bar}(R_{s})\propto K_{0}(R_{s})\times\begin{cases}1&R\leq R_{\rm
end}\\\ \exp(-(R-R_{\rm end})^{2}/h^{2}_{\rm end})&R>R_{\rm end}\\\
\end{cases}$ (2.8)

where $K_{0}$ is the modified Bessel function of the second kind and where:

$R_{s}=\left(R_{\perp}^{C_{\parallel}}+\left(\frac{\absolutevalue{z}}{z_{b}}\right)^{C_{\parallel}}\right)^{1/C_{\parallel}}$
(2.9)
$R_{\perp}=\left(\left(\frac{\absolutevalue{x}}{x_{b}}\right)^{C_{\perp}}+\left(\frac{\absolutevalue{y}}{y_{b}}\right)^{C_{\perp}}\right)^{1/C_{\perp}}\,.$
(2.10)

The free parameters are $C_{\parallel}$, $C_{\perp}$, $x_{b}$, $y_{b}$,
$z_{b}$ and $R_{\rm end}$ with $h_{\rm end}$ fixed at $\sqrt{\frac{1}{2}}$
kpc. The effective radius is $R_{s}$; the scale lengths are $x_{b}$, $y_{b}$,
and $z_{b}$; and $C_{\perp}$ and $C_{\parallel}$ are the face-on and edge-on
shape parameters. The bar shape is elliptical in the corresponding direction
when $C_{\perp},C_{\parallel}=2$, diamond-shaped when
$C_{\perp},C_{\parallel}<2$, and boxy when $C_{\perp},C_{\parallel}>2$. The
Gaussian function with scale length $h_{\rm end}$ in Eq. (2.8) truncates the
bar at radius $R_{\rm end}$. The modified Bessel function was also used in Cao
et al. [39] to model the distribution of red clump giants, but with no cutoff
and with $C_{\parallel}=4$ and $C_{\perp}=2$. For the long bar we use [40]:

$\rho_{\rm
long~{}bar}(x,y,z)\propto\exp\left(-\left(\left(\frac{\absolutevalue{x}}{x_{\rm
lb}}\right)^{C_{\perp,\rm lb}}+\left(\frac{\absolutevalue{y}}{y_{\rm
lb}}\right)^{C_{\perp,\rm lb}}\right)^{1/C_{\perp,\rm
lb}}\right)\exp\left(-\frac{\absolutevalue{z}}{z_{\rm
lb}}\right)\textrm{Cut}\left(\frac{R-R_{\rm out}}{\sigma_{\rm
out}}\right)\textrm{Cut}\left(\frac{R_{\rm in}-R}{\sigma_{\rm in}}\right)$
(2.11)

where $C_{\perp,\rm lb}$, $x_{\rm lb}$, $y_{\rm lb}$, $z_{\rm lb}$, $R_{\rm
out}$ and $R_{\rm in}$ are free parameters, $\sigma_{\rm out}=\sigma_{\rm
in}=\sqrt{\frac{1}{2}}$ kpc and:

$\textrm{Cut}\left(x\right)=\begin{cases}\exp(-x^{2})&x>0\\\ 1&x\leq 0\\\
\end{cases}$ (2.12)

Note that Fujii et. al [33] did not need to include any bar components in
their initial conditions as the bar components evolve naturally over the 10 Gy
from their Hernquist bulge, disk and halo model initial conditions.

Finally, we have a disk with a central hole:

$\rho_{\rm
disk}(x,y,z)\propto\exp(-R^{2}/2\sigma_{r}^{2})\exp(-\absolutevalue{z}/z_{0})H(x,y)$
(2.13)

where $\sigma_{r}$ and $z_{0}$ are free parameters and for the hole we use the
form adopted by Freudenreich [41]:

$H(x,y)=1-\exp\left(-\left(R_{H}/O_{R}\right)^{O_{N}}\right)$ (2.14)

with:

$R_{H}^{2}=(x)^{2}+(\epsilon y)^{2}$ (2.15)

where $\epsilon$, $O_{R}$ and $O_{N}$ are also free parameters. Note that
Fujii et. al [33] did not have a hole in their initial disk model. This hole
evolves naturally through the process of the formation of the bar over their
10 Gy simulation.

The total number of particles in our simulations are fixed. So we do not have
to include the number of particles as part of our likelihood. Therefore the
probability of having an $N$-body particle at position ${x,y,z}$ will be
proportional to the density of our model ($\rho$) at ${x,y,z}$. We have for
each component of the model a parameter giving the probability a particle is
from that component. We treat the probability of an $N$-body particle being
from a component of the density distribution as parameters. These parameters,
$P({\rm Disk})$, $P({\rm Bar})$, $P({\rm Hernquist})$ and $P({\rm
Long~{}Bar})$, have a Dirichlet prior [42]. This prior constrains

$P({\rm Disk})+P({\rm Bar})+P({\rm Hernquist})+P({\rm Long~{}Bar})=1$

and is uniformly distributed over any values of these parameter satisfying
this condition. The log-likelihood is then:

$\log(L)=\sum_{i}^{N}\log(\rho(x_{i},y_{i},z_{i}))$ (2.16)

where $x_{i}$, $y_{i}$ and $z_{i}$ are the coordinates of a particle, $N$ is
the number of particles, and $\rho$ is the density of the model:

$\rho(x,y,z)=P({\rm Disk})\rho_{\rm disk}(x,y,z)+P({\rm Bar})\rho_{\rm
bar}(x,y,z)+P({\rm Hernquist})\rho_{\rm Hernquist}(x,y,z)+P({\rm
Long~{}Bar})\rho_{\rm long~{}bar}(x,y,z)$ (2.17)

All scale parameters are given a prior so they are uniform in $\log(\theta)$
where $\theta\in\\{a_{b},x_{b},y_{b},z_{b},x_{\rm lb},y_{\rm lb},z_{\rm
lb},\sigma_{r},z_{0},O_{R}\\}$ and this implies $p(\theta)\propto 1/\theta$.
In calculating the likelihood, we do not include particles for which $R>12$
kpc or $\absolutevalue{z}>3$. It can be seen in Fujii et al. [33] that the
scale height of the disk may start to decline between $10\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 3.0pt\hbox{$\sim$}}\hss}\hbox{\raise
2.0pt\hbox{$<$}}}}R\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$<$}}}}15$ kpc. We also don’t
want the fit to be affected by particles that may have been kicked well out of
the galaxy. The likelihood ($L$) is insensitive to being multiplied by a
constant but that constant has to be the same for all parameters of our
combined model. To accommodate this we normalize each density component such
that

$\int_{R\leq 12~{}{\rm kpc},\absolutevalue{z}\leq 3~{}{\rm
kpc}}\rho_{i}(x,y,z)\,{\rm d}x\,{\rm d}y\,{\rm d}z=1$

where $i\in$ {disk, bar, Hernquist, long bar}. This integral is estimated with
importance sampling. We use a set of random numbers which are transformed into
the points at which we evaluate the density models in order to estimate the
normalization constant. In order to stabilize the estimation of the likelihood
function these numbers are always the same every time we perform the
importance sampling within a particular chain.

After running the $N$-body simulations, we shift the coordinates of the
particles so that the center of mass is at the origin, then rotate so the bar
is along the $x$-axis. The bar angle is estimated using the method described
in Fujii et al. [33]. However, we add four parameters that we expect to be
near zero to allow a further shift in the center and clockwise rotation of the
model. These are $\alpha$, $x_{\rm center}$, $y_{\rm center}$ and $z_{\rm
center}$, with the latter three parameters in parsecs, so:

$\displaystyle x_{\rm data}$ $\displaystyle=\cos(\alpha)x+\sin(\alpha)y+x_{\rm
center}/1000$ (2.18) $\displaystyle y_{\rm data}$
$\displaystyle=-\sin(\alpha)x+\cos(\alpha)y+y_{\rm center}/1000$
$\displaystyle z_{\rm data}$ $\displaystyle=z+z_{\rm center}/1000$

where $x_{\rm data}$, $y_{\rm data}$ and $z_{\rm data}$ are coordinates in the
coordinate system of the $N$-body simulation. In estimating the peculiar
velocity distribution scale parameter, $\hat{\sigma_{v}}$, above, we assumed
$x\approx x_{\rm data}$, $y\approx y_{\rm data}$ and $z\approx z_{\rm data}$.

To compare our model fits to the final $N$-body simulations we sample from the
model directly, then we bin the samples. We used a simple MCMC to sample from
the bar and long bar. Since the disk (excluding the hole) can be turned into
an invertible cumulative distributions in $R$ and $z$, and the Hernquist bulge
in $r$, we used inverse transform sampling for those. The disk hole was
handled by accepting the disk samples with probability $H(x,y)$, i.e. we used
rejection sampling for the disk.

Our MCMC algorithm used to fit the models is similar to that of Foreman-Mackay
et al. [43] with a mixture of the Differential Evolution [44] and snooker
updates [45]. We use a simple annealing method in which we divide the log-
likelihood by a temperature $T$ which is gradually reduced to $1$. This occurs
during the first half of each Markov chain, which we discard. See appendix B
for more details.

##  3 Results

In Table 1 we present the kick velocities $\sigma_{k}$ that produce peculiar
velocity distributions close to $\sigma_{v}=77\pm 6$ km s-1 as estimated by
P20. We display the rotation curves at $t=10$ Gyr for the three initial
condition models (MWa, MWb and MWc0.8) in Fig. 1. The central values in Table
1 were used to run $N$-body simulations with a larger number of particles to
which we fitted a parametric model. The fitted parameters are shown in Tables
2, 3 and 4. There are three potential sources of uncertainty in the model
parameters: the posterior density function, variation in the likelihood
between chains as a result of the importance sampling method used to estimate
the normalization constant for each model component, and the possibility that
Markov chains may get stuck in different local likelihood maxima. We found
that for most parameters the posterior distributions overlapped significantly
or, in many cases, were indistinguishable. For a few parameters we had outlier
chains; these can be seen in the tables as parameters with large, highly
asymmetric uncertainties. For the kicked distributions, the disk model hole
parameters were very uncertain with different Markov chains settling on a wide
range of different values. As this hole tended to be smaller and/or less
sharp, we simply removed it, in the kicked cases, by setting $H(x,y)=1$. We
found the long bar model often acts like a second disk component when fitting
to the kicked particles, this disk-like long bar was still allowed a hole
through the parameter $R_{\rm in}$.

We ran the $N$-body simulations for the two kick rate scenarios separately,
but we combine the two sets of Markov chains for the “No Kick” columns in the
parameter tables. Typically, the posterior distributions were very similar for
Markov chains generated using these two sets of data, slightly expanding our
$68\%$ intervals. However, we found a difference of order 10 pc in the
location of the center between simulations. For this reason, we only report
the change in parameters $\alpha$, $x_{\rm center}$, $y_{\rm center}$ and
$z_{\rm center}$ for the fits to the kicked distributions by subtracting off
the median of the corresponding Markov chains. For $\alpha$ the fitted $68\%$
intervals for no kick were $-2^{\circ}$ to $-1^{\circ}$ for all models,
$\absolutevalue{x_{\rm center}}$ and $\absolutevalue{y_{\rm center}}$ were
always $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$<$}}}}40$ pc, and
$\absolutevalue{z_{\rm center}}$ was $<2$ pc. We are not concerned about this
difference as it is negligible compared to the bar extent.

We show the final density maps generated from the MWa initial conditions in
Fig. 2. Those for MWb and MWc0.8 were similar and are given in the appendix as
Figs. A.1 and A.2 respectively. These figures include both the $N$-body
particle data as well as the fitted models. These density maps were produced
by binning particles (either the $N$-body particles or particles drawn from
the fitted model) within $0.25$ kpc of the $x$-$y$, $x$-$z$ and $y$-$z$
planes. The density maps for the fitted models were generated by taking the
average in each bin for simulated data generated using 500 random parameter
sets from our Markov chains in each case. In Fig. 3 we show the $N$-body
simulation data and simulated particle distributions along the $x$, $y$ and
$z$ axes for the MWa case with the kicks occurring at the beginning. For the
model we show the predicted mean number of particles in each bin and also the
standard deviation of that prediction. This can be used to test the model fit
as it provides the 68% posterior predictive intervals [46]. The corresponding
uniform kick rate, MWb, and MWc0.8 cases are displayed in the appendix as
Figs. A.3, A.4, A.5, A.6 and A.7. In these figures, we bin all particles
within $0.25$ kpc in the two perpendicular axes. We display in Fig. 4 the
integrated flux along lines of sight in the central $50^{\circ}\times
50^{\circ}$ of the galaxy. The Sun is placed at a distance of $7.9$ kpc, at an
angle relative to the bar of $20^{\circ}$ and at a height of $15$ pc [10].
This figure was generated by binning particles in Galactic latitude and
longitude with weights of $1/d^{2}$ where $d$ is the distance of a particle
from the Sun. We excluded particles within $1$ kpc of the Sun to reduce noise.
The corresponding figures for MWb and MWc0.8 are display in the appendix as
Figs. A.8 and A.9 respectively.

| Kick At Beginning | Uniform Kick Rate
---|---|---
MWa | $93\pm 10$ | $85\pm 9$
MWb | $97\pm 10$ | $84\pm 10$
MWc0.8 | $95\pm 10$ | $84\pm 9$
Table 1:  Kick velocity Maxwell distribution (Eq. (2.1)) parameters ($\sigma_{k}$ in km/s) that produce a peculiar velocity distribution where $\sigma_{v}=77\pm 6$ km/s.  Parameter | No Kick | Kick At Beginning | Uniform Kick Rate
---|---|---|---
$P({\rm Disk})$ | $0.499\begin{subarray}{c}+0.006\\\ -0.005\end{subarray}$ | $0.418\begin{subarray}{c}+0.006\\\ -0.006\end{subarray}$ | $0.365\begin{subarray}{c}+0.006\\\ -0.004\end{subarray}$
$P({\rm Bar})$ | $0.354\begin{subarray}{c}+0.004\\\ -0.005\end{subarray}$ | $0.3220\begin{subarray}{c}+0.0019\\\ -0.0021\end{subarray}$ | $0.364\begin{subarray}{c}+0.005\\\ -0.007\end{subarray}$
$P({\rm Hernquist})$ | $0.0145\begin{subarray}{c}+0.0013\\\ -0.0018\end{subarray}$ | $0.059\begin{subarray}{c}+0.003\\\ -0.002\end{subarray}$ | $0.0647\begin{subarray}{c}+0.0022\\\ -0.0024\end{subarray}$
$P({\rm Long~{}Bar})$ | $0.1323\begin{subarray}{c}+0.0025\\\ -0.0025\end{subarray}$ | $0.201\begin{subarray}{c}+0.007\\\ -0.007\end{subarray}$ | $0.206\begin{subarray}{c}+0.007\\\ -0.010\end{subarray}$
$\sigma_{r}$ (kpc) | $4.92\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$ | $5.66\begin{subarray}{c}+0.05\\\ -0.04\end{subarray}$ | $6.07\begin{subarray}{c}+0.05\\\ -0.08\end{subarray}$
$z_{0}$ (kpc) | $0.2238\begin{subarray}{c}+0.0008\\\ -0.0027\end{subarray}$ | $1.203\begin{subarray}{c}+0.023\\\ -0.013\end{subarray}$ | $1.222\begin{subarray}{c}+0.018\\\ -0.004\end{subarray}$
$O_{R}$ (kpc) | $2.82\begin{subarray}{c}+0.08\\\ -0.11\end{subarray}$ | $-$ | $-$
$O_{N}$ | $4.3\begin{subarray}{c}+0.3\\\ -0.2\end{subarray}$ | $-$ | $-$
$\epsilon$ | $0.780\begin{subarray}{c}+0.020\\\ -0.013\end{subarray}$ | $-$ | $-$
$a_{b}$ (kpc) | $0.205\begin{subarray}{c}+0.009\\\ -0.010\end{subarray}$ | $0.509\begin{subarray}{c}+0.019\\\ -0.010\end{subarray}$ | $0.555\begin{subarray}{c}+0.022\\\ -0.021\end{subarray}$
$C_{\perp}$ | $1.84\begin{subarray}{c}+0.02\\\ -0.05\end{subarray}$ | $1.831\begin{subarray}{c}+0.013\\\ -0.015\end{subarray}$ | $1.821\begin{subarray}{c}+0.017\\\ -0.015\end{subarray}$
$C_{\parallel}$ | $3.08\begin{subarray}{c}+0.11\\\ -0.11\end{subarray}$ | $2.49\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$ | $2.542\begin{subarray}{c}+0.019\\\ -0.019\end{subarray}$
$x_{b}$ (kpc) | $0.557\begin{subarray}{c}+0.008\\\ -0.008\end{subarray}$ | $0.685\begin{subarray}{c}+0.005\\\ -0.005\end{subarray}$ | $0.673\begin{subarray}{c}+0.007\\\ -0.011\end{subarray}$
$y_{b}$ (kpc) | $0.367\begin{subarray}{c}+0.006\\\ -0.003\end{subarray}$ | $0.455\begin{subarray}{c}+0.003\\\ -0.003\end{subarray}$ | $0.445\begin{subarray}{c}+0.006\\\ -0.011\end{subarray}$
$z_{b}$ (kpc) | $0.2557\begin{subarray}{c}+0.0014\\\ -0.0010\end{subarray}$ | $0.3328\begin{subarray}{c}+0.0014\\\ -0.0015\end{subarray}$ | $0.3155\begin{subarray}{c}+0.0026\\\ -0.0029\end{subarray}$
$R_{\rm end}$ (kpc) | $2.00\begin{subarray}{c}+0.05\\\ -0.06\end{subarray}$ | $4.6\begin{subarray}{c}+0.4\\\ -0.3\end{subarray}$ | $4.9\begin{subarray}{c}+0.6\\\ -0.6\end{subarray}$
$\Delta\alpha$ ($\deg$) | $-$ | $-0.34\begin{subarray}{c}+0.15\\\ -0.14\end{subarray}$ | $-0.07\begin{subarray}{c}+0.13\\\ -0.13\end{subarray}$
$\Delta x_{\rm center}$ (pc) | $-$ | $0.8\begin{subarray}{c}+1.0\\\ -1.1\end{subarray}$ | $-0.7\begin{subarray}{c}+1.0\\\ -1.0\end{subarray}$
$\Delta y_{\rm center}$ (pc) | $-$ | $-1.8\begin{subarray}{c}+0.9\\\ -0.9\end{subarray}$ | $-0.3\begin{subarray}{c}+0.7\\\ -0.7\end{subarray}$
$\Delta z_{\rm center}$ (pc) | $-$ | $-0.9\begin{subarray}{c}+0.6\\\ -0.5\end{subarray}$ | $2.7\begin{subarray}{c}+0.4\\\ -0.5\end{subarray}$
$x_{\rm lb}$ (kpc) | $5\begin{subarray}{c}+5\\\ -1\end{subarray}$ | $2.78\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$ | $2.46\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$
$y_{\rm lb}$ (kpc) | $1.43\begin{subarray}{c}+0.16\\\ -0.07\end{subarray}$ | $2.78\begin{subarray}{c}+0.04\\\ -0.04\end{subarray}$ | $2.49\begin{subarray}{c}+0.05\\\ -0.04\end{subarray}$
$z_{\rm lb}$ (kpc) | $0.315\begin{subarray}{c}+0.029\\\ -0.010\end{subarray}$ | $0.541\begin{subarray}{c}+0.006\\\ -0.011\end{subarray}$ | $0.495\begin{subarray}{c}+0.010\\\ -0.014\end{subarray}$
$C_{\perp,\rm lb}$ | $0.88\begin{subarray}{c}+0.04\\\ -0.11\end{subarray}$ | $1.764\begin{subarray}{c}+0.026\\\ -0.025\end{subarray}$ | $1.983\begin{subarray}{c}+0.028\\\ -0.026\end{subarray}$
$R_{\rm out}$ (kpc) | $2.86\begin{subarray}{c}+0.04\\\ -0.05\end{subarray}$ | $7.46\begin{subarray}{c}+0.06\\\ -0.04\end{subarray}$ | $8.37\begin{subarray}{c}+0.06\\\ -0.06\end{subarray}$
$R_{\rm in}$ (kpc) | $1.676\begin{subarray}{c}+0.018\\\ -0.016\end{subarray}$ | $2.09\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$ | $2.03\begin{subarray}{c}+0.07\\\ -0.12\end{subarray}$
Table 2:  Fitted parameters for model MWa. We have used the median of the MCMC chains for the central value and also included 68% confidence intervals.  Parameter | No Kick | Kick At Beginning | Uniform Kick Rate
---|---|---|---
$P({\rm Disk})$ | $0.506\begin{subarray}{c}+0.010\\\ -0.005\end{subarray}$ | $0.6144\begin{subarray}{c}+0.0024\\\ -0.0026\end{subarray}$ | $0.36\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$
$P({\rm Bar})$ | $0.322\begin{subarray}{c}+0.005\\\ -0.003\end{subarray}$ | $0.244\begin{subarray}{c}+0.007\\\ -0.003\end{subarray}$ | $0.317\begin{subarray}{c}+0.006\\\ -0.003\end{subarray}$
$P({\rm Hernquist})$ | $0.0140\begin{subarray}{c}+0.0012\\\ -0.0009\end{subarray}$ | $0.0470\begin{subarray}{c}+0.0022\\\ -0.0021\end{subarray}$ | $0.0819\begin{subarray}{c}+0.0017\\\ -0.0016\end{subarray}$
$P({\rm Long~{}Bar})$ | $0.154\begin{subarray}{c}+0.012\\\ -0.010\end{subarray}$ | $0.094\begin{subarray}{c}+0.004\\\ -0.006\end{subarray}$ | $0.24\begin{subarray}{c}+0.04\\\ -0.04\end{subarray}$
$\sigma_{r}$ (kpc) | $5.64\begin{subarray}{c}+0.06\\\ -0.07\end{subarray}$ | $5.36\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$ | $6.9\begin{subarray}{c}+0.7\\\ -0.3\end{subarray}$
$z_{0}$ (kpc) | $0.2249\begin{subarray}{c}+0.0003\\\ -0.0003\end{subarray}$ | $1.116\begin{subarray}{c}+0.004\\\ -0.014\end{subarray}$ | $1.30\begin{subarray}{c}+0.05\\\ -0.06\end{subarray}$
$O_{R}$ (kpc) | $3.19\begin{subarray}{c}+0.05\\\ -0.09\end{subarray}$ | $-$ | $-$
$O_{N}$ | $4.25\begin{subarray}{c}+0.09\\\ -0.17\end{subarray}$ | $-$ | $-$
$\epsilon$ | $0.762\begin{subarray}{c}+0.019\\\ -0.019\end{subarray}$ | $-$ | $-$
$a_{b}$ (kpc) | $0.202\begin{subarray}{c}+0.004\\\ -0.005\end{subarray}$ | $0.585\begin{subarray}{c}+0.023\\\ -0.025\end{subarray}$ | $0.614\begin{subarray}{c}+0.014\\\ -0.011\end{subarray}$
$C_{\perp}$ | $2.04\begin{subarray}{c}+0.08\\\ -0.10\end{subarray}$ | $2.02\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$ | $1.873\begin{subarray}{c}+0.016\\\ -0.017\end{subarray}$
$C_{\parallel}$ | $3.85\begin{subarray}{c}+0.14\\\ -0.08\end{subarray}$ | $3.06\begin{subarray}{c}+0.05\\\ -0.05\end{subarray}$ | $2.90\begin{subarray}{c}+0.04\\\ -0.04\end{subarray}$
$x_{b}$ (kpc) | $0.548\begin{subarray}{c}+0.029\\\ -0.009\end{subarray}$ | $0.587\begin{subarray}{c}+0.026\\\ -0.009\end{subarray}$ | $0.710\begin{subarray}{c}+0.010\\\ -0.006\end{subarray}$
$y_{b}$ (kpc) | $0.339\begin{subarray}{c}+0.006\\\ -0.005\end{subarray}$ | $0.406\begin{subarray}{c}+0.013\\\ -0.005\end{subarray}$ | $0.411\begin{subarray}{c}+0.010\\\ -0.005\end{subarray}$
$z_{b}$ (kpc) | $0.243\begin{subarray}{c}+0.004\\\ -0.002\end{subarray}$ | $0.306\begin{subarray}{c}+0.007\\\ -0.003\end{subarray}$ | $0.303\begin{subarray}{c}+0.004\\\ -0.002\end{subarray}$
$R_{\rm end}$ (kpc) | $2.53\begin{subarray}{c}+0.10\\\ -0.15\end{subarray}$ | $3.37\begin{subarray}{c}+0.10\\\ -0.08\end{subarray}$ | $5.2\begin{subarray}{c}+0.6\\\ -0.4\end{subarray}$
$\Delta\alpha$ ($\deg$) | $-$ | $-0.04\begin{subarray}{c}+0.13\\\ -0.13\end{subarray}$ | $0.30\begin{subarray}{c}+0.11\\\ -0.12\end{subarray}$
$\Delta x_{\rm center}$ (pc) | $-$ | $4.4\begin{subarray}{c}+1.2\\\ -1.2\end{subarray}$ | $2.5\begin{subarray}{c}+1.2\\\ -1.3\end{subarray}$
$\Delta y_{\rm center}$ (pc) | $-$ | $1.8\begin{subarray}{c}+0.9\\\ -0.8\end{subarray}$ | $-1.2\begin{subarray}{c}+0.9\\\ -0.8\end{subarray}$
$\Delta z_{\rm center}$ (pc) | $-$ | $-0.6\begin{subarray}{c}+0.7\\\ -0.8\end{subarray}$ | $1.4\begin{subarray}{c}+0.5\\\ -0.5\end{subarray}$
$x_{\rm lb}$ (kpc) | $2.7\begin{subarray}{c}+0.8\\\ -0.3\end{subarray}$ | $2.1\begin{subarray}{c}+0.3\\\ -0.1\end{subarray}$ | $2.78\begin{subarray}{c}+0.04\\\ -0.02\end{subarray}$
$y_{\rm lb}$ (kpc) | $1.14\begin{subarray}{c}+0.17\\\ -0.05\end{subarray}$ | $0.93\begin{subarray}{c}+0.08\\\ -0.03\end{subarray}$ | $2.65\begin{subarray}{c}+0.07\\\ -0.04\end{subarray}$
$z_{\rm lb}$ (kpc) | $0.414\begin{subarray}{c}+0.006\\\ -0.016\end{subarray}$ | $0.656\begin{subarray}{c}+0.011\\\ -0.020\end{subarray}$ | $0.56\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$
$C_{\perp,\rm lb}$ | $0.97\begin{subarray}{c}+0.03\\\ -0.11\end{subarray}$ | $1.61\begin{subarray}{c}+0.04\\\ -0.04\end{subarray}$ | $1.86\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$
$R_{\rm out}$ (kpc) | $3.66\begin{subarray}{c}+0.22\\\ -0.11\end{subarray}$ | $4.67\begin{subarray}{c}+0.05\\\ -0.05\end{subarray}$ | $9.32\begin{subarray}{c}+0.07\\\ -0.25\end{subarray}$
$R_{\rm in}$ (kpc) | $1.69\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$ | $1.78\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$ | $1.71\begin{subarray}{c}+0.10\\\ -0.06\end{subarray}$
Table 3:  Fitted parameters for model MWb.  Parameter | No Kick | Kick At Beginning | Uniform Kick Rate
---|---|---|---
$P({\rm Disk})$ | $0.555\begin{subarray}{c}+0.003\\\ -0.009\end{subarray}$ | $0.467\begin{subarray}{c}+0.012\\\ -0.010\end{subarray}$ | $0.357\begin{subarray}{c}+0.013\\\ -0.003\end{subarray}$
$P({\rm Bar})$ | $0.318\begin{subarray}{c}+0.004\\\ -0.005\end{subarray}$ | $0.262\begin{subarray}{c}+0.003\\\ -0.005\end{subarray}$ | $0.305\begin{subarray}{c}+0.003\\\ -0.002\end{subarray}$
$P({\rm Hernquist})$ | $0.0099\begin{subarray}{c}+0.0006\\\ -0.0005\end{subarray}$ | $0.0603\begin{subarray}{c}+0.0021\\\ -0.0025\end{subarray}$ | $0.0619\begin{subarray}{c}+0.0019\\\ -0.0020\end{subarray}$
$P({\rm Long~{}Bar})$ | $0.116\begin{subarray}{c}+0.014\\\ -0.004\end{subarray}$ | $0.210\begin{subarray}{c}+0.014\\\ -0.010\end{subarray}$ | $0.277\begin{subarray}{c}+0.005\\\ -0.017\end{subarray}$
$\sigma_{r}$ (kpc) | $5.49\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$ | $6.11\begin{subarray}{c}+0.07\\\ -0.09\end{subarray}$ | $6.99\begin{subarray}{c}+0.06\\\ -0.19\end{subarray}$
$z_{0}$ (kpc) | $0.2246\begin{subarray}{c}+0.0007\\\ -0.0022\end{subarray}$ | $1.266\begin{subarray}{c}+0.016\\\ -0.028\end{subarray}$ | $1.262\begin{subarray}{c}+0.026\\\ -0.022\end{subarray}$
$O_{R}$ (kpc) | $2.5\begin{subarray}{c}+0.4\\\ -0.16\end{subarray}$ | $-$ | $-$
$O_{N}$ | $4.9\begin{subarray}{c}+2.6\\\ -1.1\end{subarray}$ | $-$ | $-$
$\epsilon$ | $0.82\begin{subarray}{c}+0.12\\\ -0.12\end{subarray}$ | $-$ | $-$
$a_{b}$ (kpc) | $0.220\begin{subarray}{c}+0.006\\\ -0.006\end{subarray}$ | $0.72\begin{subarray}{c}+0.01\\\ -0.04\end{subarray}$ | $0.74\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$
$C_{\perp}$ | $1.88\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$ | $1.861\begin{subarray}{c}+0.020\\\ -0.020\end{subarray}$ | $1.909\begin{subarray}{c}+0.024\\\ -0.022\end{subarray}$
$C_{\parallel}$ | $3.37\begin{subarray}{c}+0.05\\\ -0.05\end{subarray}$ | $2.69\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$ | $2.69\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$
$x_{b}$ (kpc) | $0.551\begin{subarray}{c}+0.012\\\ -0.014\end{subarray}$ | $0.657\begin{subarray}{c}+0.006\\\ -0.007\end{subarray}$ | $0.620\begin{subarray}{c}+0.006\\\ -0.005\end{subarray}$
$y_{b}$ (kpc) | $0.342\begin{subarray}{c}+0.007\\\ -0.008\end{subarray}$ | $0.434\begin{subarray}{c}+0.005\\\ -0.007\end{subarray}$ | $0.378\begin{subarray}{c}+0.005\\\ -0.003\end{subarray}$
$z_{b}$ (kpc) | $0.2395\begin{subarray}{c}+0.0015\\\ -0.0018\end{subarray}$ | $0.3217\begin{subarray}{c}+0.0023\\\ -0.0023\end{subarray}$ | $0.2845\begin{subarray}{c}+0.0025\\\ -0.0018\end{subarray}$
$R_{\rm end}$ (kpc) | $2.03\begin{subarray}{c}+0.07\\\ -0.10\end{subarray}$ | $4.81\begin{subarray}{c}+0.29\\\ -0.25\end{subarray}$ | $5.5\begin{subarray}{c}+2.1\\\ -0.9\end{subarray}$
$\Delta\alpha$ ($\deg$) | $-$ | $-0.79\begin{subarray}{c}+0.18\\\ -0.18\end{subarray}$ | $0.32\begin{subarray}{c}+0.13\\\ -0.13\end{subarray}$
$\Delta x_{\rm center}$ (pc) | $-$ | $1.4\begin{subarray}{c}+1.4\\\ -1.3\end{subarray}$ | $1.0\begin{subarray}{c}+1.2\\\ -1.1\end{subarray}$
$\Delta y_{\rm center}$ (pc) | $-$ | $-0.2\begin{subarray}{c}+0.9\\\ -0.9\end{subarray}$ | $0.5\begin{subarray}{c}+0.8\\\ -0.8\end{subarray}$
$\Delta z_{\rm center}$ (pc) | $-$ | $2.4\begin{subarray}{c}+0.7\\\ -0.6\end{subarray}$ | $-0.5\begin{subarray}{c}+0.5\\\ -0.5\end{subarray}$
$x_{\rm lb}$ (kpc) | $3\begin{subarray}{c}+7\\\ -1\end{subarray}$ | $2.533\begin{subarray}{c}+0.023\\\ -0.022\end{subarray}$ | $2.533\begin{subarray}{c}+0.020\\\ -0.021\end{subarray}$
$y_{\rm lb}$ (kpc) | $1.5\begin{subarray}{c}+0.4\\\ -0.2\end{subarray}$ | $2.584\begin{subarray}{c}+0.027\\\ -0.025\end{subarray}$ | $2.417\begin{subarray}{c}+0.021\\\ -0.022\end{subarray}$
$z_{\rm lb}$ (kpc) | $0.335\begin{subarray}{c}+0.006\\\ -0.005\end{subarray}$ | $0.564\begin{subarray}{c}+0.019\\\ -0.006\end{subarray}$ | $0.530\begin{subarray}{c}+0.010\\\ -0.011\end{subarray}$
$C_{\perp,\rm lb}$ | $0.95\begin{subarray}{c}+0.07\\\ -0.11\end{subarray}$ | $1.964\begin{subarray}{c}+0.029\\\ -0.029\end{subarray}$ | $1.954\begin{subarray}{c}+0.022\\\ -0.022\end{subarray}$
$R_{\rm out}$ (kpc) | $3.0\begin{subarray}{c}+0.1\\\ -0.4\end{subarray}$ | $7.88\begin{subarray}{c}+0.07\\\ -0.09\end{subarray}$ | $9.13\begin{subarray}{c}+0.13\\\ -0.08\end{subarray}$
$R_{\rm in}$ (kpc) | $1.61\begin{subarray}{c}+0.03\\\ -0.06\end{subarray}$ | $1.91\begin{subarray}{c}+0.05\\\ -0.07\end{subarray}$ | $1.60\begin{subarray}{c}+0.05\\\ -0.03\end{subarray}$
Table 4:  Fitted parameters for model MWc0.8.  Figure 1: The rotation curves
at $t=10$ Gyr for the three initial condition models MWa, MWb and MWc0.8. They
are consistent with the local circular velocity of the Sun which is measured
to be $238\pm 15$ km s-1 [47]. The distance between the Sun and the Galactic
Center is approximately $8$ kpc.  Figure 2: Final density map of particles
with no kick, a kick at the beginning and a uniform kick rate generated from
the MWa initial conditions. Every second row shows the fitted model.  Figure
3: Profile along $x$, $y$ and $z$ axes with kicks occurring at the beginning
for the case generated from the MWa initial conditions. We show both the final
$N$-body simulation data and data simulated using the final fitted model. For
the final fitted model we show the mean number of particles in each bin and
the standard deviation.  Figure 4: Final flux distribution in Galactic
coordinates generated from the MWa initial conditions. The contours for each
distribution are at $1$, $2$, $4$, $8$ and $16$ times the mean in this region.
The Sun is placed at a distance of $7.9$ kpc, at an angle relative to the bar
of $20^{\circ}$ and at a height of $15$ pc.

In Table 5 we show the change in $-2\log(L)$ when replacing $\rho_{\rm
bar}(R_{s})$ with a range of different models from Freudenreich [41] and Cao
et al. [39]. As there is variation in the likelihood between chains, we also
show the standard deviation in $-2\log(L)$. Our choice of $\rho_{\rm
bar}(R_{s})\propto K_{0}(R_{s})$ is clearly preferred over the others. The
worst form, where $\rho_{\rm bar}(R_{s})\propto\exp(R_{s}^{-n})$, was entirely
removed with $P({\rm Bar})=0$ and the long bar component taking over the fit
in the central region.

$\rho_{\rm bar}(R_{s})$ | Mean $-2(\Delta\log(L))$ | Standard Deviation
---|---|---
$K_{0}(R_{s})$ | $0$ | $739$
$\exp(-R_{s})$ | $3581$ | $1606$
$\sech^{2}(R_{s})$ | $12795$ | $1985$
$\exp(-0.5R_{s}^{2})$ | $19369$ | $391$
$(1+R_{s}^{n})^{-1}$ | $20403$ | $953$
$\exp(R_{s}^{-n})$ | $57114$ | $3152$
Table 5:  Change in mean $-2\log(L)$ using different bar models $\rho_{\rm
bar}(R_{s})$. The mean was taken over all the samples in the MCMC chains.

In order to estimate the kick effects on the actual Milky Way Galactic bar we
used a linear fit to our simulation results of the form

$\theta_{{\rm kicked},i}=\alpha\theta_{i}+\beta$ (3.1)

where $\theta_{{\rm kicked},i}$ are the bulge parameters, $x_{b}$, $y_{b}$,
and $z_{b}$ for the kicked distribution and $\theta_{i}$ are the corresponding
parameters in the non-kicked case. The $\alpha$ and $\beta$ were found by
performing a least squared fit for the values given in Tables 2, 3, and 4. We
also did a similar fit for $C_{\perp}$ and $C_{\parallel}$. The results are
shown in Table 6 and Fig. 5. A prediction for the Milky Way bar parameters
found in ref. [39] are shown in Table 7. The predicted line of sight contours
for the kicked and unkicked Milky Way bar are shown in Fig. 6.

Parameter | $\alpha_{i}$ | $\beta_{i}$
---|---|---
$x_{b},y_{b},z_{b}$ | 1.13$\pm$0.05 | 0.03$\pm$0.02
$C_{\perp},C_{\parallel}$ | 0.56$\pm$0.02 | 0.81$\pm$0.05
Table 6:  Least square fit values with 68% confidence intervals for Eq. 3.1 fitted to the points shown in Figure 5.  Parameter |  $x_{b}$ (kpc) |  $y_{b}$ (kpc) |  $z_{b}$ (kpc) | $C_{\perp}$ | $C_{\parallel}$
---|---|---|---|---|---
Not kicked | 0.67 | 0.29 | 0.27 | 2 | 4
Kicked | 0.79$\pm$0.02 | 0.36$\pm$0.01 | 0.35 $\pm$ 0.01 | 1.93$\pm$0.02 | 3.05$\pm$0.03
Table 7:  Predictions with 68% confidence intervals for the kicked spatial
distribution for the Milky Way bar model found in ref. [39] using Eq. 3.1 and
the parameter values given in Table 6.

Figure 5: Simulation parameters with 68% confidence interval bands for
straight line model fits. The closed symbol values are obtained from values
given in Tables 2, 3, and 4. The predictions for the ref. [39] model of the
Milky Way Galaxy are given as open symbols.  Figure 6: Line of sight contours
for model of the Milky Way bar given by Cao et al. [39]. Both the model and
its kicked version, obtained from the parameters in Table 7, are shown. As in
ref. [39] the Galactic Centre is taken to be 8.13 kpc away and the angle of
the bulge to be 29∘. Contours are from 1 to 32 times the mean of the
corresponding case given in steps of factors of 2.

##  4 Discussion and Conclusions

Our goal was to investigate the effect of neutron star birth kicks on the
distribution of MSPs in the Galactic center. We began by running $N$-body
simulations with small populations of particles kicked with a range of scales
in order to estimate the required Maxwellian kick to produce a peculiar
velocity distribution similar to that of resolved gamma-ray MSPs. We then
reran the simulations with a larger number of particles at the required kick
velocity scale and used MCMC to fit the data with a model.

We used three initial condition models intended to approximate the Milky Way,
these were the MWa, MWb and MWc0.8 initial condition models of Fujii et al.
[33]. Our results were consistent with theirs as can be seen, for example, by
comparing our Fig. 1 to the top left hand panels of their Figs. 1, 2, and 3
which are also consistent with the observed circular rotation velocity at the
Sun’s location. In Cao et al. [39] the bar scale lengths for a modified Bessel
function of the second kind model fitted to red clump giant data are $0.67$,
$0.29$ and $0.27$ for the $x$, $y$ and $z$ axes respectively, with the
parameters $C_{\parallel}$ and $C_{\perp}$ fixed at $4$ and $2$. Our fits to
$N$-body models without kicks find $(x_{b},y_{b},z_{b})$ of $(0.56,0.37,0.26)$
for MWa, $(0.55,0.34,0.24)$ for MWb and $(0.55,0.34,0.24)$ for MWc0.8. We also
have $C_{\parallel}$ between $3$ and $4$, producing a boxy structure in
$x$-$z$ and $y$-$z$, this is visible in Fig. 2 for the MWa case and Figs. A.1
and A.2 for the MWb and MWc0.8 cases. The other shape parameter $C_{\perp}$
was relatively close to $2$ in all cases, resulting in a more elliptical shape
in $x$-$y$. We found it was necessary to extend the bar structure using the
long bar component given in Eq. 2.11. Without it the bar scale parameters
would be larger, but the bar would not be long enough to explain the structure
seen for $2\mbox{ kpc}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise
2.0pt\hbox{$<$}}}}\absolutevalue{x}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$<$}}}}4$ kpc in Fig. 3 in the
MWa, kicked at the beginning case, and Figs. A.3, A.4, A.5, A.6 and A.7 in the
uniform kick rate, MWb and MWc0.8 cases. We find the disk scale height to be
$0.22$ kpc for the MWa, MWb, and MWc0.8 unkicked cases. This is at the lower
end of the range of $220$ to $450$ pc given in Bland-Hawthorn and Gerhard
[47]. There is also a small spherically symmetric Hernquist component $\sim
1\%$ of the particles in the region of interest for the unkicked cases.

Eckner et al. [31] argued using the virial theorem that kicks $\langle
v^{2}\rangle\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$<$}}}}(70{\rm
km~{}s^{-1}})^{2}$ would lead to a “smoothing” of the distribution of
$700$–$900$ pc. We show the effect of a $400$ pc and an $800$ pc Gaussian
smoothing on the fitted bulge (bar plus Hernquist bulge) distribution in Figs.
7 and 8 for the kicked at the beginning case. The uniform kick rate case is
shown in Fig. A.10. Those figures also show the bulge component of the models
with and without kicks for comparison. It is clear that a Gaussian smoothing
kernel will remove the peak that survives in the $N$-body simulations of
kicked distributions and, particularly for the $800$ pc case, will produce an
apparently spherically symmetric bulge. From the peculiar velocity data we
inferred kicks that are larger than assumed by Eckner et al. [31] with
$\sigma_{k}$ at around $80$–$100$ km s-1 ($\langle
v^{2}\rangle=3\sigma_{k}^{2}$ for a Maxwell distribution where the angular
brackets signify the mean value) so the smoothing effect of the Gaussian would
be even more severe. We show an even smaller Gaussian smoothing of $200$ pc in
Figs. 9 and 10 for the kicked at beginning case and in Fig. A.11 for the
uniform kick rate case. We also show in Fig. 11 the profile for particles with
smaller kick scales between $0$ km s-1 and $80$ km s-1 for the kicked at the
beginning case. The corresponding uniform kick rate case is shown in Fig.
A.12. In these two figures, each kick scale has only $4\times 10^{5}$
particles; therefore, to reduce noise, the bins in each of the other two
dimensions are twice as big as in previous single dimensional plots, and
particles within $0.5$ kpc (previously $0.25$ kpc) of the axis are included.
The profiles along the $z$ axis in particular show that there is a reduction
in the slope as the kick velocities increase, along the other two axes the
general increase in scaleheight is seen as a reduction in density. These
results demonstrates that Gaussian smoothing is not a good way of modelling a
kicked version of a boxy bulge/bar template.

In every case, the bar fitted to the kicked data is both broader, with larger
scale parameters $x_{b}$, $y_{b}$ and $z_{b}$, and less boxy, with smaller
$C_{\parallel}$. For example, for the MWa initial conditions
$(x_{b},y_{b},z_{b})$ increases from $(0.56,0.37,0.26)$ to $(0.69,0.46,0.33)$
and $(0.67,0.45,0.32)$ for the kick at beginning case and the uniform kick
rate case respectively, while $C_{\parallel}$ declines from $3.08$ to $2.49$
and $2.54$. The spherically symmetric Hernquist bulge increases from $\sim
1\%$ of the particles to $6\%$ for MWa. Like the bar, it becomes broader with
$a_{b}$ increasing from around $0.2$ kpc to $0.51$ kpc and $0.56$ kpc. For the
other two models similar changes occur, $P({\rm Long~{}Bar})$ and $a_{b}$ both
increase significantly. In MSP model A1 of P20, the disk parameters were
$\sigma_{r}=4.5\begin{subarray}{c}+0.5\\\ -0.4\end{subarray}$ kpc and
$z_{0}=0.71\begin{subarray}{c}+0.11\\\ -0.09\end{subarray}$. In the current
article, after being kicked, the disk scale heights $z_{0}$ of all models
increase from $0.22$ kpc to $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$>$}}}}{1}$ kpc, while
$\sigma_{r}$ is in the range $5$–$7$ kpc. However, we find that in all kicked
cases, except for MWb with kicks occurring at the beginning, the long bar
behaves like a relatively thin disk component. We have $x_{\rm lb}\approx
y_{\rm lb}$, $C_{\perp,\rm lb}\sim 2$ and $R_{\rm out}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 3.0pt\hbox{$\sim$}}\hss}\hbox{\raise
2.0pt\hbox{$>$}}}}{7}$, resulting in a density $\sim\exp(-R/R_{0})$ in $R$ for
scalelength $R_{0}$. These scalelengths would then be between approximately
$2.4$ kpc and $2.8$ kpc. For comparison, in Bland-Hawthorn and Gerhard [47]
the Milky Way disk scalelength is reported as $2.6\pm 0.5$ kpc. The
exponential scaleheights of these “long bars” range between about $0.5$–$0.6$
kpc. In Fig. 3, for the MWa kicked at the beginning case and in Figs. A.3,
A.5, A.6 and A.7 for the uniform kick rate and MWb, and MWc0.8 case, there may
be, to varying degrees, an excess of kicked particles over the model in the
region of $2\mbox{ kpc}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise
2.0pt\hbox{$<$}}}}\absolutevalue{x}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$<$}}}}4$ kpc. This can be
seen as a correlated run of roughly two standard deviation difference between
the prediction and the data in this region.

Although, a more complicated parametric model could remedy these fit defects,
our main aim in this article was to estimate the effect of the MSP kicks on
their distribution in the Milky Way Galactic bar. As can be seen for the MWa
case from Fig. 4, the line of sight integral of the models are a good match to
the simulations. In particular the noisy simulation contours scatter in an
unbiased way around the smooth model contours. Similar results can be seen for
the MWb and MWc0.8 cases in Figs. A.8 and A.9 respectively.

We used the same bar parametric model as Cao et al. [39] who fit the best fit
parameters to the red clump luminosity density distribution of the Galactic
bar measured by the Optical Gravitational Lensing Experiment (OGLE) III
survey. Comparing to their fit, our final simulation fits of the unkicked
particles had somewhat different bulge parameters. But, there appears to be a
linear relationship between the unkicked scale parameters $x_{b}$, $y_{b}$,
and $z_{b}$ and their kicked counterparts. Similarly, there appears to be a
linear relationship between $C_{\perp}$, $C_{\parallel}$, and their kicked
counterparts. Therefore, we were able to estimate the Milky Way Galactic bar
kicked parameters as shown in Fig. 5 and Table 7. The residuals in this least
squares fit are partially due to the systematic error of our model misfit and
also the Cao et al. [39] model misfit to the red clump data. As can be seen,
there is more scatter in the $x_{b}$ parameter. This is not unexpected due to
the already mentioned degeneracy with the long bar. Also, as can be seen,
estimating the Milky Way bar kicked $x_{b}$ parameter did involve a reasonable
amount of extrapolation and so future simulations which have a larger $x_{b}$
will be needed to check it. A made-to-measure [48, 49] approach may be needed.
This would also be advantageous as it could take into account the X/peanut
shaped morphology of the bar [50, 51, 40] as done in ref. [52]. There is some
preliminary evidence that the X-shape may improve the fit to the Fermi-LAT
gamma-ray data [10].

In conclusion, we used $N$-body simulations to explore the effect of a Maxwell
distributed kick on the distribution of MSPs in the Galactic center. We find
that while a $700$–$900$ pc Gaussian smoothing of the stellar mass would be
too aggressive, the bulge distribution of the kicked particles is slightly
broader and less boxy. From these results, we expect the GCE to deviate by a
small amount from the stellar mass spatial distribution in the Galactic
center. Also, as can seen from Table 7, we would not expect the GCE to appear
spherically symmetric due to the MSP kicks as that would require
$x_{b}=y_{b}=z_{b}$ and $C_{\perp}=C_{\parallel}=2$ which are far from our
inferred points relative to their error bars.

The amount of spatial smoothing of the bulge MSPs will depend on the
proportion of MSPs in the bulge that are made from the recycling channel and
the proportion that are made from the accretion induced collapse channel.
Motivated by similarities between the bulge and disk population seen by P20 we
have assumed this mixture has the same proportions as the disk MSPs. If the
GCE is due to bulge MSPs, its morphology could be used to check our smoothing
prediction by comparing if there are any deviations between the GCE morphology
and the stellar spatial distribution. A complication to this approach would be
the possibility of some smearing of the GCE due to cosmic ray electron
diffusion [53, 54]. An additional complication is that if the MSP is spun up
by a captured star then the MSP spatial distribution would be proportional to
the stellar density squared [31, 8]. We have been assuming that the MSPs
formed in a binary system and so have a spatial distribution proportional to
the stellar density. Eventually, once the bulge MSPs are resolved [30, 55],
comparing their spatial distribution to the stellar distribution should
provide independent information to more robustly estimate the natal kick
distribution.

Figure 7: Final profile along $x$, $y$ and $z$ axes with kicks occurring at
the beginning generated from the MWa initial conditions. Here we show the
fitted bulge components, which consist of the bar plus Hernquist bulge, as
well as the no kick bulge smoothed with $400$ pc and $800$ pc Gaussians. We
show both $N$-body simulation data and data simulated using the fitted model.
For the fitted model we show the mean number of particles in each bin and the
standard deviation.  Figure 8: Final bulge, consisting of bar plus Hernquist
bulge, flux distribution in Galactic coordinates generated from the MWa
initial conditions. We also show the no kick bulge smoothed with $400$ pc and
$800$ pc Gaussians. The contours for each distribution are at $1$, $2$, $4$,
$8$ and $16$ times the mean in this region. The Sun is placed at a distance of
$7.9$ kpc, at an angle relative to the bar of $20^{\circ}$ and at a height of
$15$ pc.  Figure 9: Final profile along $x$, $y$ and $z$ axes with kicks
occurring at the beginning generated from the MWa initial conditions. Here we
show the fitted bulge components as well as the no kick bulge smoothed with a
$200$ pc Gaussian. We show both $N$-body simulation data and data simulated
using the fitted model. For the fitted model we show the mean number of
particles in each bin and the standard deviation.  Figure 10: Final bulge flux
distribution in Galactic coordinates generated from the MWa initial
conditions. We also show the no kick bulge smoothed with a $200$ pc Gaussian.
The contours for each distribution are at $1$, $2$, $4$, $8$ and $16$ times
the mean in this region. The Sun is placed at a distance of $7.9$ kpc, at an
angle relative to the bar of $20^{\circ}$ and at a height of $15$ pc.  Figure
11: Final profile along $x$, $y$ and $z$ axes with kicks occurring at the
beginning for kicks between $0$ km s-1 and $80$ km s-1 in $20$ km s-1
increments. Generated from the MWa initial conditions.

## Acknowledgments

We thank Roland Crocker and Oscar Macias for helpful comments. HP was
supported by a University of Canterbury Doctoral Scholarship.

## References

  * [1] L. Goodenough and D. Hooper, _Possible Evidence For Dark Matter Annihilation In The Inner Milky Way From The Fermi Gamma Ray Space Telescope_ , _arXiv e-prints_ (Oct., 2009) arXiv:0910.2998, [0910.2998].
  * [2] D. Hooper and T. Linden, _Origin of the gamma rays from the galactic center_ , _Physical Review D_ 84 (Dec, 2011) .
  * [3] K. N. Abazajian and M. Kaplinghat, _Detection of a Gamma-Ray Source in the Galactic Center Consistent with Extended Emission from Dark Matter Annihilation and Concentrated Astrophysical Emission_ , _Phys. Rev._ D86 (2012) 083511, [1207.6047].
  * [4] C. Gordon and O. Macias, _Dark Matter and Pulsar Model Constraints from Galactic Center Fermi-LAT Gamma Ray Observations_ , _Phys. Rev._ D88 (2013) 083521, [1306.5725].
  * [5] K. N. Abazajian, _The Consistency of Fermi-LAT Observations of the Galactic Center with a Millisecond Pulsar Population in the Central Stellar Cluster_ , _JCAP_ 1103 (2011) 010, [1011.4275].
  * [6] O. Macias, C. Gordon, R. M. Crocker, B. Coleman, D. Paterson, S. Horiuchi et al., _Galactic bulge preferred over dark matter for the galactic centre gamma-ray excess_ , _Nature Astronomy_ 2 (Mar, 2018) 387–392.
  * [7] R. Bartels, E. Storm, C. Weniger and F. Calore, _The Fermi-LAT GeV excess as a tracer of stellar mass in the Galactic bulge_ , _Nature Astronomy_ (Aug., 2018) , [1711.04778].
  * [8] O. Macias, S. Horiuchi, M. Kaplinghat, C. Gordon, R. M. Crocker and D. M. Nataf, _Strong evidence that the galactic bulge is shining in gamma rays_ , _JCAP_ 2019 (Sept., 2019) 042, [1901.03822].
  * [9] K. N. Abazajian, S. Horiuchi, M. Kaplinghat, R. E. Keeley and O. Macias, _Strong constraints on thermal relic dark matter from Fermi-LAT observations of the Galactic Center_ , _Phys. Rev. D_ 102 (2020) 043012, [2003.10416].
  * [10] B. Coleman, D. Paterson, C. Gordon, O. Macias and H. Ploeg, _Maximum entropy estimation of the Galactic bulge morphology via the VVV Red Clump_ , _MNRAS_ 495 (May, 2020) 3350–3372, [1911.04714].
  * [11] F. Calore, F. Donato and S. Manconi, _Dissecting the inner Galaxy with $\gamma$-ray pixel count statistics_, _arXiv e-prints_ (Feb., 2021) arXiv:2102.12497, [2102.12497].
  * [12] M. Di Mauro, _Characteristics of the Galactic Center excess measured with 11 years of F e r m i -LAT data_ , _Phys. Rev. D_ 103 (Mar., 2021) 063029, [2101.04694].
  * [13] R. Bartels, S. Krishnamurthy and C. Weniger, _Strong support for the millisecond pulsar origin of the Galactic center GeV excess_ , _Phys. Rev. Lett._ 116 (2016) 051102, [1506.05104].
  * [14] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer and W. Xue, _Evidence for Unresolved $\gamma$-Ray Point Sources in the Inner Galaxy_, _Phys. Rev. Lett._ 116 (2016) 051103, [1506.05124].
  * [15] R. K. Leane and T. R. Slatyer, _Dark Matter Strikes Back at the Galactic Center_ , _arXiv e-prints_ (Apr., 2019) arXiv:1904.08430, [1904.08430].
  * [16] R. K. Leane and T. R. Slatyer, _The enigmatic Galactic Center excess: Spurious point sources and signal mismodeling_ , _Phys. Rev. D_ 102 (Sept., 2020) 063019, [2002.12371].
  * [17] R. K. Leane and T. R. Slatyer, _Spurious Point Source Signals in the Galactic Center Excess_ , _Phys. Rev. Lett.,_ 125 (Sept., 2020) 121105, [2002.12370].
  * [18] M. Buschmann, N. L. Rodd, B. R. Safdi, L. J. Chang, S. Mishra-Sharma, M. Lisanti et al., _Foreground mismodeling and the point source explanation of the Fermi Galactic Center excess_ , _Phys. Rev. D_ 102 (July, 2020) 023023, [2002.12373].
  * [19] L. J. Chang, S. Mishra-Sharma, M. Lisanti, M. Buschmann, N. L. Rodd and B. R. Safdi, _Characterizing the nature of the unresolved point sources in the Galactic Center: An assessment of systematic uncertainties_ , _Phys. Rev. D_ 101 (Jan., 2020) 023014, [1908.10874].
  * [20] D. Bhattacharya and E. P. J. van den Heuvel, _Formation and evolution of binary and millisecond radio pulsars_ , _Physics Reports_ 203 (Jan., 1991) 1–124.
  * [21] L. Ferrario and D. Wickramasinghe, _The birth properties of Galactic millisecond radio pulsars_ , _Mon. Not. Roy. Astron. Soc._ 375 (01, 2007) 1009–1016.
  * [22] J. R. Hurley, C. A. Tout, D. T. Wickramasinghe, L. Ferrario and P. D. Kiel, _Formation of binary millisecond pulsars by accretion-induced collapse of white dwarfs_ , _Mon. Not. Roy. Astron. Soc._ 402 (02, 2010) 1437–1448.
  * [23] A. J. Ruiter, L. Ferrario, K. Belczynski, I. R. Seitenzahl, R. M. Crocker and A. I. Karakas, _On the formation of neutron stars via accretion-induced collapse in binaries_ , _Mon. Not. Roy. Astron. Soc._ 484 (Mar., 2019) 698–711, [1802.02437].
  * [24] A. Gautam, R. M. Crocker, L. Ferrario, A. J. Ruiter, H. Ploeg, C. Gordon et al., _Millisecond Pulsars from Accretion Induced Collapse naturally explain the Galactic Center Gamma-ray Excess_ , _arXiv e-prints_ (June, 2021) arXiv:2106.00222, [2106.00222].
  * [25] C. Fryer, W. Benz, M. Herant and S. A. Colgate, _What can the accretion-induced collapse of white dwarfs really explain?_ , _Astrophys. J._ 516 (may, 1999) 892–899.
  * [26] F. S. Kitaura, H.-T. Janka and W. Hillebrandt, _Explosions of o-ne-mg cores, the crab supernova, and subluminous type ii-p supernovae_ , _Astronomy & Astrophysics_ 450 (Apr, 2006) 345–350.
  * [27] A. G. Lyne and D. R. Lorimer, _High birth velocities of radio pulsars_ , _Nature_ 369 (May, 1994) 127–129.
  * [28] Wongwathanarat, A., Janka, H.-Th. and Müller, E., _Three-dimensional neutrino-driven supernovae: Neutron star kicks, spins, and asymmetric ejection of nucleosynthesis products_ , _Astron. and Astrophys._ 552 (2013) A126.
  * [29] E. Bear and N. Soker, _Neutron star natal kick and jets in core collapse supernovae_ , _Astrophys. J._ 855 (mar, 2018) 82.
  * [30] H. Ploeg, C. Gordon, R. Crocker and O. Macias, _Comparing the galactic bulge and galactic disk millisecond pulsars_ , _JCAP_ 2020 (Dec., 2020) 035, [2008.10821].
  * [31] C. Eckner, X. Hou, P. D. Serpico, M. Winter, G. Zaharijas, P. Martin et al., _Millisecond Pulsar Origin of the Galactic Center Excess and Extended Gamma-Ray Emission from Andromeda: A Closer Look_ , _ApJ_ 862 (July, 2018) 79, [1711.05127].
  * [32] J. Bédorf, E. Gaburov and S. Portegies Zwart, _A sparse octree gravitational N-body code that runs entirely on the GPU processor_ , _Journal of Computational Physics_ 231 (Apr., 2012) 2825–2839, [1106.1900].
  * [33] M. S. Fujii, J. Bédorf, J. Baba and S. Portegies Zwart, _Modelling the Milky Way as a dry Galaxy_ , _MNRAS_ 482 (Jan., 2019) 1983–2015, [1807.10019].
  * [34] K. Kuijken and J. Dubinski, _Nearly Self-Consistent Disc / Bulge / Halo Models for Galaxies_ , _MNRAS_ 277 (Dec., 1995) 1341.
  * [35] L. M. Widrow and J. Dubinski, _Equilibrium Disk-Bulge-Halo Models for the Milky Way and Andromeda Galaxies_ , _ApJ_ 631 (Oct., 2005) 838–855, [astro-ph/0506177].
  * [36] L. M. Widrow, B. Pym and J. Dubinski, _Dynamical Blueprints for Galaxies_ , _ApJ_ 679 (June, 2008) 1239–1259, [0801.3414].
  * [37] F. Verbunt, A. Igoshev and E. Cator, _The observed velocity distribution of young pulsars_ , _A &A _ 608 (Dec., 2017) A57, [1708.08281].
  * [38] L. Hernquist, _An Analytical Model for Spherical Galaxies and Bulges_ , _ApJ_ 356 (June, 1990) 359.
  * [39] L. Cao, S. Mao, D. Nataf, N. J. Rattenbury and A. Gould, _A New Photometric Model of the Galactic Bar using Red Clump Giants_ , _Mon. Not. Roy. Astron. Soc._ 434 (2013) 595–605, [1303.6430].
  * [40] C. Wegg, O. Gerhard and M. Portail, _The structure of the Milky Way’s bar outside the bulge_ , _MNRAS_ 450 (July, 2015) 4050–4069, [1504.01401].
  * [41] H. T. Freudenreich, _A COBE Model of the Galactic Bar and Disk_ , _ApJ_ 492 (Jan., 1998) 495–510, [astro-ph/9707340].
  * [42] M. Betancourt, _Cruising the simplex: Hamiltonian Monte Carlo and the Dirichlet distribution_ , in _American Institute of Physics Conference Series_ (P. Goyal, A. Giffin, K. H. Knuth and E. Vrscay, eds.), vol. 1443 of _American Institute of Physics Conference Series_ , pp. 157–164, May, 2012.  1010.3436.  DOI.
  * [43] D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, _emcee: The MCMC Hammer_ , _Publications of the Astronomical Society of the Pacific_ 125 (Mar., 2013) 306, [1202.3665].
  * [44] C. ter Braak, _A markov chain monte carlo version of the genetic algorithm differential evolution: Easy bayesian computing for real parameter spaces_ , _Statistics and Computing_ 16 (09, 2006) 239–249.
  * [45] C. ter Braak and J. Vrugt, _Differential evolution markov chain with snooker updater and fewer chains_ , _Statistics and Computing_ 18 (12, 2008) 435–446.
  * [46] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari and D. B. Rubin, _Bayesian Data Analysis_.  CRC Press, 3 ed., 2013.
  * [47] J. Bland-Hawthorn and O. Gerhard, _The Galaxy in Context: Structural, Kinematic, and Integrated Properties_ , _Annual Review of Astronomy and Astrophysics_ 54 (Sept., 2016) 529–596, [1602.07702].
  * [48] D. Syer and S. Tremaine, _Made-to-measure N-body systems_ , _MNRAS_ 282 (Sept., 1996) 223–233, [astro-ph/9605061].
  * [49] F. de Lorenzi, V. P. Debattista, O. Gerhard and N. Sambhus, _NMAGIC: a fast parallel implementation of a $\chi$2-made-to-measure algorithm for modelling observational data_, _MNRAS_ 376 (Mar., 2007) 71–88, [astro-ph/0701582].
  * [50] D. M. Nataf, A. Udalski, A. Gould, P. Fouqué and K. Z. Stanek, _The Split Red Clump of the Galactic Bulge from OGLE-III_ , _Astrophys. J. lett._ 721 (Sept., 2010) L28–L32, [1007.5065].
  * [51] A. McWilliam and M. Zoccali, _Two Red Clumps and the X-shaped Milky Way Bulge_ , _ApJ_ 724 (Dec., 2010) 1491–1502, [1008.0519].
  * [52] M. Portail, C. Wegg, O. Gerhard and I. Martinez-Valpuesta, _Made-to-measure models of the Galactic box/peanut bulge: stellar and total mass in the bulge region_ , _MNRAS_ 448 (Mar., 2015) 713–731, [1502.00633].
  * [53] D. Song, O. Macias and S. Horiuchi, _Inverse Compton emission from millisecond pulsars in the Galactic bulge_ , _Phys. Rev. D_ 99 (June, 2019) 123020, [1901.07025].
  * [54] O. Macias, H. van Leijen, D. Song, S. Ando, S. Horiuchi and R. M. Crocker, _Cherenkov telescope array sensitivity to the putative millisecond pulsar population responsible for the galactic center excess_ , _Monthly Notices of the Royal Astronomical Society_ (May, 2021) .
  * [55] F. Calore, M. Di Mauro, F. Donato, J. W. T. Hessels and C. Weniger, _Radio Detection Prospects for a Bulge Population of Millisecond Pulsars as Suggested by Fermi-LAT Observations of the Inner Galaxy_ , _ApJ_ 827 (Aug., 2016) 143, [1512.06825].
  * [56] H. Haario, E. Saksman and J. Tamminen, _An adaptive metropolis algorithm_ , _Bernoulli_ 7 (2001) 223–242.

##  Appendix A Uniform kick rate, MWb, and MWc0.8 figures

As the uniform kick rate, MWb, and MWc0.8 gave similar results to the MWa
kicked at beginning case we have moved many of their figures to this appendix.

Figure A.1: Final density map of particles with no kick, a kick at the
beginning and a uniform kick rate generated from the MWb initial conditions.
Every second row shows the fitted model.  Figure A.2: Final density map of
particles with no kick, a kick at the beginning and a uniform kick rate
generated from the MWc0.8 initial conditions. Every second row shows the
fitted model.  Figure A.3: Final profile along $x$, $y$ and $z$ axes with a
uniform kick rate generated from the MWa initial conditions. We show both
$N$-body simulation data and data simulated using the fitted model. For the
fitted model we show the mean number of particles in each bin and the standard
deviation.  Figure A.4: Final profile along $x$, $y$ and $z$ axes with kicks
occurring at the beginning generated from the MWb initial conditions. We show
both $N$-body simulation data and data simulated using the fitted model. For
the fitted model we show the mean number of particles in each bin and the
standard deviation.  Figure A.5: Final profile along $x$, $y$ and $z$ axes
with a uniform kick rate generated from the MWb initial conditions. We show
both $N$-body simulation data and data simulated using the fitted model. For
the fitted model we show the mean number of particles in each bin and the
standard deviation.  Figure A.6: Final profile along $x$, $y$ and $z$ axes
with kicks occurring at the beginning generated from the MWc0.8 initial
conditions. We show both $N$-body simulation data and data simulated using the
fitted model. For the fitted model we show the mean number of particles in
each bin and the standard deviation.  Figure A.7: Final profile along $x$, $y$
and $z$ axes with a uniform kick rate generated from the MWc0.8 initial
conditions. We show both $N$-body simulation data and data simulated using the
fitted model. For the fitted model we show the mean number of particles in
each bin and the standard deviation.  Figure A.8: Final flux distribution in
Galactic coordinates generated from the MWb initial conditions. The contours
for each distribution are at $1$, $2$, $4$, $8$ and $16$ times the mean in
this region. The Sun is placed at a distance of $7.9$ kpc, at an angle
relative to the bar of $20^{\circ}$ and at a height of $15$ pc.  Figure A.9:
Final flux distribution in Galactic coordinates generated from the MWc0.8
initial conditions. The contours for each distribution are at $1$, $2$, $4$,
$8$ and $16$ times the mean in this region. The Sun is placed at a distance of
$7.9$ kpc, at an angle relative to the bar of $20^{\circ}$ and at a height of
$15$ pc.  Figure A.10: Final profile along $x$, $y$ and $z$ axes with a
uniform kick rate generated from the MWa initial conditions. Here we show the
fitted bulge components, consisting of the bar plus Hernquist bulge, as well
as the no kick bulge smoothed with $400$ pc and $800$ pc Gaussians. We show
both $N$-body simulation data and data simulated using the fitted model. For
the fitted model we show the mean number of particles in each bin and the
standard deviation.  Figure A.11: Final profile along $x$, $y$ and $z$ axes
with a uniform kick rate generated for the MWa initial conditions. Here we
show the fitted bulge components as well as the no kick bulge smoothed with a
$200$ pc Gaussian. We show both $N$-body simulation data and data simulated
using the fitted model. For the fitted model we show the mean number of
particles in each bin and the standard deviation.  Figure A.12: Final profile
along $x$, $y$ and $z$ axes with a uniform kick rate for kicks between $0$ km
s-1 and $80$ km s-1 in $20$ km s-1 increments. Generated from the MWa initial
conditions.

##  Appendix B Markov Chain Monte Carlo Method

In previous work, such as ref. [30], we used for our MCMC the adaptive
Metropolis algorithm of Haario et al. [56]. In this article, however, we found
it was necessary to replace adaptive Metropolis algorithm with an alternative
algorithm to ensure rapid convergence to the peak likelihood region of the
parameter space. The MCMC algorithm used in this article is similar to that of
Foreman-Mackay et al. [43] with a mixture of the Differential Evolution [44]
and snooker updates [45]. Instead of performing a single random walk through
the parameter space where proposed moves are accepted with a probability based
on the likelihood of the current point of the chain and the proposed next
point, we use an ensemble of $K$ “walkers” where a proposed update for a
walker $j$ depends on the distribution of the other walkers.

A single step of the stretch move update suggested in Foreman-Mackay et al.
[43] involves updating the $K$ walkers sequentially. Let $\boldsymbol{x}_{j}$
be the state of walker $j$, an update for walker $\boldsymbol{x}_{j}$ is
performed as follows:

  1. 1.

Draw $k$ from $1,2,...,K$ where $k\neq j$

  2. 2.

Draw $z$ from the probability density function with parameter $a$ (Foreman-
Mackay et al. [43] suggest $a=2$):

$g(z)\propto\begin{cases}\frac{1}{\sqrt{z}}&z\in\left[\frac{1}{a},a\right]\\\
0&{\rm otherwise}\\\ \end{cases}$ (B.1)

  3. 3.

Calculate proposal
$\boldsymbol{y}=\boldsymbol{x}_{k}+z\left(\boldsymbol{x}_{j}-\boldsymbol{x}_{k}\right)$

  4. 4.

Calculate acceptance probability $r$:

$r=z^{d-1}\frac{p(\boldsymbol{y})}{p(\boldsymbol{x}_{j})}$ (B.2)

where $d$ is the number of dimensions

  5. 5.

Set $\boldsymbol{x}_{j}=\boldsymbol{y}$ with probability $\min(1,r)$

We found better results using a mixture of two alternative updates: $80\%$ the
Differential Evolution update of ter Braak [44] and $20\%$ the snooker update
of ter Braak and Vrugt [45].

To update walker $\boldsymbol{x}_{j}$ using the Differential Evolution update:

  1. 1.

Draw $k$ and $l$ from $1,2,...,K$ where $k\neq j$, $l\neq j$ and $l\neq k$

  2. 2.

Propose
$\boldsymbol{y}=\boldsymbol{x}_{j}+\gamma\left(\boldsymbol{x}_{k}-\boldsymbol{x}_{l}\right)+\boldsymbol{e}$
where $\gamma$ is a parameter and where $\boldsymbol{e}$ is drawn from a small
$d$ dimensional symmetric probability distribution

  3. 3.

Calculate acceptance probability $r$:

$r=\frac{p(\boldsymbol{y})}{p(\boldsymbol{x}_{j})}$ (B.3)

  4. 4.

Set $\boldsymbol{x}_{j}=\boldsymbol{y}$ with probability $\min(1,r)$

We drew $\boldsymbol{e}$ from a $d$ dimensional Gaussian with standard
deviation $10^{-5}$ in each dimension. In case the likelihood distribution had
multiple modes, we used $\gamma=1$ with probability $0.1$ as suggested by ter
Braak [44], otherwise we used the default value of $\gamma=2.38/\sqrt{2d}$.

Using the snooker update, we update $\boldsymbol{x}_{j}$ as follows:

  1. 1.

Draw $k$, $l$ and $m$ from $1,2,...,K$ with no index repeated or equal to $j$

  2. 2.

Calculate the orthogonal projections of $\boldsymbol{x}_{l}$ and
$\boldsymbol{x}_{m}$ onto the line $\boldsymbol{x}_{j}-\boldsymbol{x}_{k}$,
$\mathrm{proj}_{\boldsymbol{x}_{j}-\boldsymbol{x}_{k}}(\boldsymbol{x}_{l})$
and
$\mathrm{proj}_{\boldsymbol{x}_{j}-\boldsymbol{x}_{k}}(\boldsymbol{x}_{m})$,
where:

$\mathrm{proj}_{\boldsymbol{u}}(\boldsymbol{v})=\frac{\boldsymbol{v}\cdot\boldsymbol{u}}{\boldsymbol{u}\cdot\boldsymbol{u}}\boldsymbol{u}$
(B.4)

  3. 3.

Propose
$\boldsymbol{y}=\boldsymbol{x}_{j}+\gamma_{s}\left(\mathrm{proj}_{\boldsymbol{x}_{j}-\boldsymbol{x}_{k}}(\boldsymbol{x}_{l})-\mathrm{proj}_{\boldsymbol{x}_{j}-\boldsymbol{x}_{k}}(\boldsymbol{x}_{m})\right)$
where $\gamma_{s}$ is a parameter

  4. 4.

Calculate acceptance probability $r$:

$r=\frac{p(\boldsymbol{y})\absolutevalue{\boldsymbol{y}-\boldsymbol{x}_{k}}^{d-1}}{p(\boldsymbol{x}_{j})\absolutevalue{\boldsymbol{x}_{j}-\boldsymbol{x}_{k}}^{d-1}}$
(B.5)

  5. 5.

Set $\boldsymbol{x}_{j}=\boldsymbol{y}$ with probability $\min(1,r)$

We use $\gamma_{s}=2.38/\sqrt{2}$ as suggested by ter Braak and Vrugt [45].

We also used a simple annealing method in which we divide the log-likelihood
by a temperature $T$ which is gradually reduced to $1$. The posterior
probability density for a parameter set $\theta$ at MCMC iteration $t$ is:

$p(\theta|\textrm{$N$-body data})\propto p(\theta)L^{1/T_{t}}$ (B.6)

where $p(\theta)$ is the prior, $L$ is the likelihood, and $T_{t}$ is the
temperature. We used a linearly decreasing $\log(T)$ from $\log(1000)$ to
$\log(1)$ during the first half of each Markov chain, which we discard. This
method allows the algorithm to explore a broad region in the parameter space,
while slowly converging to the desired posterior distribution where $T=1$.
This appeared to help the Markov chains avoid getting stuck in local
likelihood maxima.