TESTING THE DARK MATTER CAUSTIC THEORY AGAINST OBSERVATIONS IN THE MILKY WAY

1.1. The Caustic Ring Dark Matter Halo

The caustic ring dark matter theory has been developed by Pierre Sikivie (Sikivie et al. 1995; Sikivie 2003, 2011; Natarajan & Sikivie 2007; Duffy & Sikivie 2008; Banik & Sikivie 2013) over the past 20 years. Although Sikivie started from the spherically symmetric self-similar models of Fillmore & Goldreich (1984) and Bertschinger (1985), it is important to note that the Sikivie model has evolved to include the physics of dark matter particles beyond a simple gravitational interaction. This physics results in dark matter flows that rotate with the stellar disk.

The caustic ring model (Duffy & Sikivie 2008) is consistent with a universe that expands with the Hubble flow and has power-law density fluctuations that create the seeds of individual galaxies. The dark matter starts out on these density fluctuations with temperature zero, so the velocity dispersion at any point in space is zero. As the dark matter falls into a galaxy it gains velocity, but since it initially occupies only three dimensions of six-dimensional phase space (three spatial dimensions and three velocity dimensions), it falls in on discrete flows that occupy only three dimensions of the six-dimensional phase space. The infalling flows, like the particles of which they are composed, fall in toward the galaxy center until they reach their closest approach, then fly out until all of the kinetic energy has been converted to gravitational potential energy, and then turn around and fall back in.

For special locations in the galaxy, the dark matter will form caustics, where the density is infinite in an infinitely small region of space (for a nonzero velocity dispersion, the caustic will not be as sharply defined). These locations are the turnaround points on the orbits of the infalling sheets of dark matter. The caustics are rings at the inner turnaround radius and spheres at the outer turnaround radius. In a self-similar model (Sikivie 2011), the positions of the caustics move outward with time, as initially more distant material collapses into the potential well of the galaxy. The number of caustics increases as more matter falls into the galaxy (with infall times of about 10⁸ years); each time the initial infalling particles reach perigalacticon, a new ring caustic forms near the center of the galaxy, at the original galactocentric distance of that earliest accreting matter. Thus, the number of caustics grows with time (about one additional caustic for each 10⁸ years); the most recently accreted matter is associated with the outermost ring. The locations of the ring caustics depend only on the slope of the evolved power spectrum of density perturbations on galaxy scales, the rotation speed of the galaxy, and the age of the universe. For the Milky Way (rotation speed of 220 km s⁻¹ at 8.5 kpc from the Galactic center) and the current standard cold dark matter (ΛCDM) model of the universe (age = 13.7 Gyr), the inner ring caustics are predicted to be at radius

$$\begin{eqnarray}&&{a}_{n}\simeq \displaystyle \frac{40\;\mathrm{kpc}}{n}\ \ \ (n=1,2,3,4,...).\end{eqnarray} \tag{ 1 }$$

The first four caustics in the Milky Way are rings at predicted distances of 40, 20, 13, and 10 kpc from the center of the Milky Way (Duffy & Sikivie 2008). In cross section, these rings have a tricusp shape, as shown in Figure 1 (top panel).

Zoom In Zoom Out Reset image size

1.2. Caustic Rings and Hierarchical Structure Formation

Although the caustic ring model (we will also refer to it as the caustic ring halo) does not include the effects of hierarchical merging, it is plausibly consistent with the merger history of the Milky Way. Large galaxy mergers will destroy dark matter caustics, similar to the way galaxy mergers disrupt disks of stars (Toth & Ostriker 1992), and any other coherent structure of collisionless particles. However, it is very important to note that any dark matter that falls in after the collision will remain intact on the three-dimensional (3D) sheets in six-dimensional (6D) phase space. A merger produces a local thickening in the phase space sheet and a velocity dispersion associated with the flow of dark matter in a caustic ring. Sikivie (2003) estimates the spread of the caustic ring radius, $\delta a$, due to the velocity dispersion, $\delta v$, of the flow of dark matter. Minor mergers will cause the caustics to spread and become fuzzy, but the caustic will become sharp again after the merger event, as material from the cold flow fills the caustic in again. Minor mergers like the Sagittarius dwarf galaxy (Ibata et al. 1994) that is currently being disrupted by the Milky Way's tidal forces would not affect the structure of the caustics.

The relationship between the redshift (z) at which the caustic ring forms and ring number (n) is approximately (P. Sikivie 2015, private communication)

$$\begin{eqnarray}&&z(n)={(3.8n-0.9)}^{0.4}-1.\end{eqnarray} \tag{ 2 }$$

Therefore, the last (n = 1 is the last ring to form) 14 caustic rings ($n=1,2,....14$) were formed after $z=3.9$, and the last two caustic rings ($n=1,2$) were formed after $z=1.1$. If there were no major mergers since $z=1.1$, then we would expect to see a ring caustic (n = 2) at 20 kpc from the Galactic center, just where we see the Monoceros Ring of stars.

In the Milky Way in particular, there is no evidence for a major merger in recent history. In particular, the boxy bulge of the Milky Way is consistent with formation from disk instabilities rather than from mergers (Shen et al. 2010; Ness et al. 2013). Even within the ΛCDM model, galaxies with the lowest ratios of bulge to total luminosity have no major mergers at $z\lt 2$ (Martig et al. 2012). Gilmore et al. (2002) and Wyse et al. (2006) have suggested that the thick disk may be the result of the last major merger of the Milky Way, 10–12 Gyr ago ($1.8\lt z\lt 4$). If this is true, the four ring caustics, at 40, 20, 13, and 10 kpc from the center of the Milky Way, should be undisturbed.

The details of the caustic ring model are still evolving; some of the effects of the Bose–Einstein condensate (BEC) model are still being incorporated in the model (Banik & Sikivie 2013). The axions in this special degenerate condensed state (the de Broglie wavelength, ${\rm{\hslash }}/p$, must be large compared to interparticle dimensions) set up a vorticity that is more permanent than the original model, when the axion self-interactions are included (the particles share angular momentum with each other). Just as water going down the drain of a sink has a curl (quasi-rigid rotation with the inner particles rotating faster than the outer particles), axions in a BEC also have a curl. In contrast, the velocity field for weakly interacting massive particles (WIMPs) has zero curl. Although the original density fluctuations are spherical by construction, the vorticity of the axions will result in an oblate dark matter halo. In addition, the baryons and the BEC dark matter may thermalize, thus explaining why disks have such a high angular momentum. Although very recent ΛCDM cosmological simulations can reproduce disks that are plausibly like those observed, it is unclear what process gives the disks rotation while leaving the dark matter halo without significant rotation (see Marinacci et al. 2014). In this paper, we do not include the modifications to the caustic ring model suggested by Banik & Sikivie (2013) because their implications on the Galactic potential have not yet been worked out.

1.3. Implications of a Caustic Ring Halo on the Nature of Dark Matter

1.4. Evidence for a Caustic Ring Dark Matter Halo

Note that there is no way to add the caustic rings to a standard halo density model, such as a Navarro–Frenk–White (NFW) profile (Navarro et al. 1996), because the acceleration calculated in the caustic ring model is not the acceleration due to an excess of matter on a ring; it is rather the acceleration due to one sheet of dark matter falling into the Milky Way. To derive the gravitational field for the entire halo, one must perform a vector sum of all of the accelerations due to all of the 20 $(n=1,2,3,4,...)$ sheets of matter that have fallen in within the last $2\times {10}^{9}$ Gyr (previous infalls are expected to be a small fraction of the total halo mass). All of the halo dark matter, including the ring caustics, is included in this final calculation of the acceleration. There is no freedom to adjust the overall shape of the halo; it simply falls out of the sum of the integrals.

Since the caustic ring model is axially symmetric, the equations for dark matter density and gravitational field in the Galaxy are conveniently expressed in cylindrical coordinates $(\rho ,z)$, centered on the Galactic center. There is one equation for the density at a position near a caustic ring and another for the density at a position far away from a caustic ring. The equation for the density near a caustic ring is given by Equation (2) in Natarajan & Sikivie (2007):

$$\begin{eqnarray}&&d{(\rho ,z)}_{\mathrm{near}}=\displaystyle \frac{1}{\rho }\displaystyle \sum _{j=1}^{N(\rho ,z)}\displaystyle \frac{{dM}}{d{\rm{\Omega }}{dt}}(\alpha ,\tau ){\left.\displaystyle \frac{\mathrm{cos}\alpha }{| D(\alpha ,\tau )| }\right|}_{({\alpha }_{j}(\rho ,z),{\tau }_{j}(\rho ,z))}.\end{eqnarray} \tag{ 3 }$$

Here the sum is over the number of flows of axions at $(\rho ,z)$, $\frac{{dM}}{d{\rm{\Omega }}{dt}}(\alpha ,\tau )$ is the rate at which mass (dM) falls into the position $(\rho ,z)$ per unit solid angle and unit time ($d{\rm{\Omega }}{dt}$), and $D(\alpha ,\tau )$ is the Jacobian determinant of the transformation between the coordinates $(\alpha ,\tau )$ and $(\rho ,z)$. Each of these quantities is a complex function of constants and coordinates in $(\alpha ,\tau )$, defined in Natarajan & Sikivie (2007). Here, α is an angle associated with the particle's last turnaround, and τ is the time when the particle crosses the z = 0 plane. The density of dark matter near a caustic ring is very large because there are many particles at turnaround points on their orbits around the Galaxy. The near equation shows this behavior since the density $d(\rho ,z)$ approaches infinity when $D(\alpha ,\tau )$ approaches 0. When $D(\alpha ,\tau )=0$, the transformation between the coordinates $(\alpha ,\tau )$ and $(\rho ,z)$ gives a set of points $(\rho ,z)$ outlined by a tricusp shape. The tricusp shape is the blue solid line in Figure 1 (top panel).

The equation for the density far from a caustic is given by Equation (3.45) in Duffy & Sikivie (2008):

$$\begin{eqnarray}&&d{(\rho ,z)}_{\mathrm{far}}=\displaystyle \frac{1}{v}\displaystyle \frac{{dM}}{d{\rm{\Omega }}{dt}}\displaystyle \frac{1}{\sqrt{{({r}^{2}-{a}^{2})}^{2}+4{a}^{2}{z}^{2}}},\end{eqnarray} \tag{ 4 }$$

where $r=\sqrt{{\rho }^{2}+{z}^{2}}$, v is the speed of the flow near the caustic, $\frac{{dM}}{d{\rm{\Omega }}{dt}}$ is the mass infall rate per unit solid angle and unit time, and a is the caustic ring radius from the Galactic center. The equations given above for density apply to one caustic ring and its associated infalling and outgoing axions. The predicted radii $({a}_{n})$, flow speeds $({v}_{n})$, and infall rates $(\frac{{d}^{2}M}{d{\rm{\Omega }}{dt}}{| }_{n})$ for the first 20 caustics in the disk of the Galaxy are listed in Table 3 of Duffy & Sikivie (2008). These parameters will now be used with the subscript n that specifies the flow number (n = 1–20).

The gravitational field for each caustic ring is similarly expressed by a near and a far equation. If the curved nature is neglected near the ring, it can be modeled as a straight tube. This leads to a gravitational field of the following form:

$$\begin{eqnarray}&&{\boldsymbol{g}}(\rho ,z)=-2G\int d\rho ^{\prime} {dz}^{\prime} d({\rho }^{\prime },z^{\prime} )\displaystyle \frac{(\rho -{\rho }^{\prime })\hat{\rho }+(z-z^{\prime} )\hat{z}}{{(\rho -{\rho }^{\prime })}^{2}+{(z-z^{\prime} )}^{2}},\end{eqnarray} \tag{ 5 }$$

(Natarajan & Sikivie 2007), where G is the gravitational constant and $d({\rho }^{\prime },z^{\prime} )$ is the dark matter density (Equation (3)). This equation describes the gravitational field from the flow of all axions flowing into and out of the ring, and not just the ring structure itself. Tam (2012) presents a method to calculate this rather difficult integral analytically. After transforming coordinates and substituting $d(\rho ,z)$ from Equation (3) into (5), Tam (2012) reduces ${\boldsymbol{g}}$ into the form

$$\begin{eqnarray}&&{\boldsymbol{g}}{(\rho ,z)}_{\mathrm{near}}=-\displaystyle \frac{8\pi G}{\rho {v}_{n}}{\displaystyle \frac{{dM}}{d{\rm{\Omega }}{dt}}| }_{n}[{I}_{\rho }(\rho ,z)\hat{\rho }+{I}_{z}(\rho ,z)\hat{z}],\end{eqnarray} \tag{ 6 }$$

where ${I}_{\rho }(\rho ,z)$ and ${I}_{z}(\rho ,z)$ are integrals defined in Tam (2012). Here we replace the b in Equation (16) from Tam (2012) with v_n. The parameter b is of the same order of magnitude as the speed of the particles in the caustic (Natarajan & Sikivie 2007). Figure 1 (top panel) shows the vector representation of ${\boldsymbol{g}}{(\rho ,z)}_{\mathrm{near}}$ for the n = 2 caustic ring, which lies at the radius a_n of 20.1 kpc. The field points toward the caustic in all directions. The tricusp is symmetric around a central point where the field goes to zero. For the n = 2 ring (Figure 1; top panel), the central point lies at $(\rho ,z)=(20.35,0)$ kpc.

The gravitational field at a position far from the caustic is

$$\begin{eqnarray}&&{\boldsymbol{g}}{(\rho ,z)}_{\mathrm{far}}=\displaystyle \frac{{A}_{n}}{s(2{a}_{n}^{2}+s)}[({r}^{2}-{a}_{n}^{2})\rho \hat{\rho }+({r}^{2}+{a}_{n}^{2})z\hat{z}],\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{A}_{n}={\left.\displaystyle \frac{4\pi G}{{v}_{n}}\displaystyle \frac{{dM}}{d{\rm{\Omega }}{dt}}\right|}_{n},\end{eqnarray} \tag{ 8 }$$

where $r=\sqrt{{\rho }^{2}+{z}^{2}}$ and $s=\sqrt{{({r}^{2}-{a}_{n}^{2})}^{2}+4{a}_{n}^{2}{z}^{2}}$ (P. Sikivie 2015, private communication).

The caustic ring halo acceleration versus Galactocentric distance in Figures 1 (bottom panel) and 2 was calculated as follows: if the distance (ρ, z) is inside the tricusp ring, use the field near the ring (Equation (6)); if the distance (ρ, z) is outside the ring, use the field far from the ring (Equation (7)). The gravitational field at a distance (ρ, z) is the vector sum from all caustic rings (n = 1–20):

$$\begin{eqnarray}&&{\boldsymbol{g}}{(\rho ,z)}_{\mathrm{halo}}=\displaystyle \sum _{n=1}^{20}{{\boldsymbol{g}}}_{n}(\rho ,z).\end{eqnarray} \tag{ 9 }$$

We spent considerable effort finding the spatial limits that made the near-caustic function match seamlessly to the far-caustic function and finding and deciphering the published model parameters for the Milky Way (Duffy & Sikivie 2008). There were essentially no tunable parameters in the model. Figure 1 (bottom panel) shows the gravitational field vectors for ${{\boldsymbol{g}}}_{\mathrm{halo}}$ as a function of position in the Milky Way. The blue lines show the positions of one cross section of the caustic rings in the $\rho z$ plane. The field generally points toward the Galactic center except near the innermost ring, where the field points away from the center. This is because there is no mass inside the n = 20 ring in our model, so objects not exactly at the galaxy center are pulled toward the innermost rings. The contribution to the halo mass from the additional caustics for $n\gt 20$ is assumed to be negligible. These additional caustics have likely been rearranged by early merger events. Although these figures show the cross section of the rings in the $\rho z$ plane, the tricusp is barely resolved on this scale. This is because the average horizontal size of the tricusp is on the order of ∼0.4 kpc. Note that the vertical size of the tricusp has been exaggerated in Figure 1 (blue lines in bottom panel); this average size is ∼0.5 kpc. The sizes of the tricusp are not predicted by the caustic theory. Duffy & Sikivie (2008) estimate these sizes from rises in the Milky Way rotation curves, which will be discussed later in this paper. The reasoning is that (at z = 0) there are discontinuities in the rotation curve at positions that lie on the boundary of the ring. For example, in Figure 1 the tricusp boundary for z = 0 lies at $\rho =20.1$ kpc and $\rho =21.1$ kpc. This implies there will be discontinuities in the rotation curve at $\rho =20.1$ kpc and $\rho =21.1$ kpc.

The source code for calculating the caustic ring halo acceleration given a position has been made available on the Milkyway@home GitHub repository (http://github.com/Milkyway-at-home/milkywayathome_client). The caustic ring halo has also been released as a gravitational potential module in the NEMO stellar dynamics toolbox (http://bima.astro.umd.edu/nemo/). Orbits and N-body simulations can be run on particles in a caustic ring halo potential using both NEMO and Milkyway@home.

The halo potential models (NFW and logarithmic) that will be compared to the caustic ring model throughout this paper are listed below. Also, X, Y, and Z are Galactocentric Cartesian coordinates for all models used in this paper.

The equation for the spherical NFW model (Navarro et al. 1996) was taken from Equation (2.67) in Binney & Tremaine (2008):

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{NFW}}=-4\pi G{\rho }_{0}{r}_{s}^{2}\displaystyle \frac{\mathrm{ln}(1+r/{r}_{s})}{r/{r}_{s}},\end{eqnarray} \tag{ 10 }$$

where $r=\sqrt{{X}^{2}+{Y}^{2}+{Z}^{2}}$, ${\rho }_{0}$ is a characteristic overdensity of the halo, and r_s is the scale radius. The parameters we used for the NFW model were adopted from Xue et al. (2008). In this paper, the observed radial velocities of Blue Horizon Branch (BHB) stars from the Sloan Digital Sky Survey (SDSS; York et al. 2000) were compared to those from simulations of the formation of galaxies like the Milky Way. In Xue et al. (2008), ${\rho }_{0}$ is defined as ${\rho }_{s}$ (Equation (14) in Xue et al. 2008), and ${r}_{s}={r}_{\mathrm{vir}}/c$. The characteristic density ${\rho }_{s}$ depends on the critical density of the universe, the matter contribution to the critical density, and the critical overdensity at virialization. The parameter r_vir is the virial radius, and c is a dimensionless concentration parameter. A Milky Way model including an NFW halo was fit to the rotation curve. Xue et al. (2008) found fit values for the virial radius and concentration. We use those fit parameters here for the NFW model. The values for the parameters in Equation (10) are listed in Table 1.

Table 1.  Halo Model Parameters

Model	Parameter	Value
NFW	Characteristic density	${\rho }_{o}=4.197\times {10}^{-3}\;{M}_{\odot }\;{\mathrm{pc}}^{-3}$
	Scale radius	${r}_{s}=22.25\;\mathrm{kpc}$
Logarithmic	Halo speed	${v}_{{\rm{l}},\mathrm{halo}}=114.0\;\mathrm{km}\;{{\rm{s}}}^{-1}$
	Z-direction flattening	${q}_{{\rm{l}}}=1.0$
Triaxial	Halo speed	${v}_{{\rm{t}},\mathrm{halo}}=121.9\;\mathrm{km}\;{{\rm{s}}}^{-1}$
	Z-direction flattening	${q}_{{\rm{t}}}=1.36$
	Scale length	$d=12.0\;\mathrm{kpc}$
	Constant	${C}_{1}=0.992947$
	Constant	${C}_{2}=0.532153$
	Constant	${C}_{3}=0.114889$

Download table as:  ASCIITypeset image

The parameters for a logarithmic (log) halo were taken from Law et al. (2005). Law et al. (2005) used a log halo in the Galaxy model to fit the orbit of the Sagittarius dwarf galaxy tidal stream. We discuss the Sagittarius stream in Section 5. Here, ${v}_{\mathrm{halo}}$ is the halo speed, q is the Z-direction flattening, and d is the halo scale length. The functional form of the logarithmic halo is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{log}}={v}_{{\rm{l}},\mathrm{halo}}^{2}\mathrm{ln}\left({X}^{2}+{Y}^{2}+\displaystyle \frac{{Z}^{2}}{{q}_{{\rm{l}}}^{2}}+{d}^{2}\right).\end{eqnarray} \tag{ 11 }$$

The parameters for a triaxial halo were taken from Law & Majewski (2010). The values ${v}_{\mathrm{halo}}$, q, and d are defined the same as they are for the logarithmic halo, and the functional form is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{triaxial}}={v}_{{\rm{t}},\mathrm{halo}}^{2}\mathrm{ln}\left({C}_{1}{X}^{2}+{C}_{2}{Y}^{2}+{C}_{3}{XY}+\displaystyle \frac{{Z}^{2}}{{q}_{{\rm{t}}}^{2}}+{d}^{2}\right).\end{eqnarray} \tag{ 12 }$$

The constants C₁, C₂, and C₃ are defined in Equations (4)–(6) of Law & Majewski (2010). These allow for rotation about the halo Z axis and flattening in the X and Y directions. These parameters were also fit by Law & Majewski (2010) to match the orbit of the Sagittarius tidal stream and are listed in Table 1.

The gravitational acceleration of the NFW and spherical (${q}_{{\rm{l}}}=1.0$) log halo models in Figure 2 was calculated using ${\boldsymbol{g}}=-{\boldsymbol{\nabla }}{\rm{\Phi }}$. It is important to note that unlike the NFW, log, and triaxial models, the caustic ring model parameters are predicted directly from theory. The caustic ring model parameters a_n, v_n, and ${\left.\frac{{dM}}{d{\rm{\Omega }}{dt}}\right|}_{n}$ (Equations (6)–(8)) follow from solutions to the equations of motion for axion flows falling into the Galaxy. The initial conditions are specified by the initial dark matter mass distribution, the rotation speed of the Galaxy, and the initial dark matter flow angular momentum. The "tunable" parameters are the angular momentum parameter j_max and the size p of the caustic ring (Duffy & Sikivie 2008). The angular momentum parameter j_max for the Milky Way was fit to the rotation curve by Sikivie et al. (1995). The model predicts the radii for the first 20 caustics (Duffy & Sikivie 2008). The contribution to the halo mass from the additional caustics for $n\gt 20$ is assumed to be negligible. The sizes p of the caustic rings are inferred from the analysis of bumps in rotation curves. Table 4 of Duffy & Sikivie (2008) lists the positions of observed rises in Milky Way rotation curves corresponding to the positions of the n = 3 and $n=5-14$ caustic rings. There should be an upward kink in the rotation at the position of one edge of the caustic and a downward kink at the position of the other edge of the caustic. In Figure 1 (top panel), this corresponds to the left edge at $\rho =20.1$ kpc and the right edge at $\rho =$ 20.4 kpc. For our calculations of the caustic ring model gravitational acceleration, we use the p values from Table 4 of Duffy & Sikivie (2008) for the n = 3 and $n=5-14$ rings. For n = 1–2, 4, and 15–20 we use ${p}_{n}/{a}_{n}=0.05$, which is the average value of the ratio ${p}_{n}/{a}_{n}$ for the other rings.

Zoom In Zoom Out Reset image size

The caustic ring halo acceleration is compared to that of an NFW halo and a spherical log halo in Figure 2. The acceleration for the NFW and log halos is given as a function of distance from the Galactic center. For the caustic ring halo, which is not spherically symmetric but is axially symmetric, we show acceleration in the plane (z = 0) as a function of cylindrical radius ρ, and along the axis of symmetry ($\rho =0$) as a function of z. In either case, the direction of the acceleration is toward the Galactic center. The obvious feature is that in the Galactic plane there are large jumps in the acceleration at the positions of the caustic rings. Because the NFW and log halos are spherically symmetric, there is no difference in these curves for z = 0 and $\rho =0$.

Given the complexity of calculating the acceleration, it is extremely gratifying and reassuring to discover that the resulting potential was plausibly similar to familiar Galactic potentials. Note that the caustic ring halo, overall, is reasonably similar to a log halo, but with jumps in the acceleration near the caustics. The shape of the dark matter profile differs tremendously from the NFW profile near the Galactic center; there is no strong central cusp in the caustic ring halo. In fact, it is just the opposite: there is essentially no dark matter in the center of a caustic ring halo. The sharp spike in Figure 2 near the center indicates an acceleration away from the Galactic center, toward the caustic rings. Figure 1 (bottom panel) shows the gravitational field vectors as a function of position in the Milky Way, showing the positive acceleration along the radius inside the innermost ring caustic. This effect disappears when bulge and disk potentials are included, as we shall see later. In the real galaxy, we can imagine that the caustic ring halo is slightly blurred due to the nonzero velocity dispersion of the dark matter. Also, past major mergers would disrupt the caustic pattern near the Galactic center.

The potential of the Milky Way Galaxy includes contributions from the dark matter halo, the disk of baryons, and the stellar bulge. In this paper, we use analytical disk and bulge potentials to model the baryons in stars and gas in the Milky Way. We use a Miyamoto & Nagai (1975) disk and Hernquist (1990) bulge:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{disk}}=-\displaystyle \frac{{{GM}}_{\mathrm{disk}}}{\sqrt{{X}^{2}+{Y}^{2}+{(a+\sqrt{{Z}^{2}+{b}^{2}})}^{2}}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{bulge}}=-\displaystyle \frac{{{GM}}_{\mathrm{bulge}}}{r+c},\end{eqnarray} \tag{ 14 }$$

where $r=\sqrt{{X}^{2}+{Y}^{2}+{Z}^{2}}$. The parameters ${M}_{\mathrm{disk}}$ and ${M}_{\mathrm{bulge}}$ are the masses of the disk and bulge; a and b are the disk scale length and height; c is the bulge scale radius; and X, Y, and Z are Galactocentric Cartesian coordinates for all models used in this paper. We will compare the results of three cases for the fixed analytical potential of the Milky Way:

• 1.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{caustic}\ \mathrm{halo}}\end{array}$
• 2.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{logarithmic}\ \mathrm{halo}}\end{array}$
• 3.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{triaxial}\ \mathrm{halo}}.\end{array}$

The three halo models used here are discussed in Section 3. The rotation curve was obtained as follows. For circular orbits, the centripetal acceleration at a distance r from the Galactic center is $g(r)=\frac{v{(r)}^{2}}{r}$, the square of the velocity at r divided by r. We get the acceleration from the Miyamoto–Nagai disk, Hernquist bulge, log halo, and NFW halo potential models using ${\boldsymbol{g}}(r)=-{\boldsymbol{\nabla }}{\rm{\Phi }}$. For the caustic ring halo, we already have the acceleration. To get the contribution from the gas in the disk, we used the velocity points of the gas disk model from Olling & Merrifield (2000). Our ${g}_{\mathrm{gas}}$ is a fifth-order polynomial fit to the HI gas model and a sixth-order polynomial fit to the ${{\rm{H}}}_{2}$ gas model. Because we have the acceleration for all of our models, we can solve for $v(r)$. We add all of the contributions to the Galaxy acceleration ${g}_{\mathrm{total}}={g}_{\mathrm{disk}}+{g}_{\mathrm{bulge}}+{g}_{\mathrm{gas}}+{g}_{\mathrm{halo}}$ and plot $v=\sqrt{{g}_{\mathrm{total}}\cdot r}$ for the rotation curve.

We use a reduced ${\chi }^{2}$ as a measure of fitness of the rotation curve data points to the models:

$$\begin{eqnarray}&&{\chi }_{v}^{2}=\displaystyle \frac{1}{N-n-1}\displaystyle \sum _{i}{\left(\displaystyle \frac{{v}_{\mathrm{model},i}-{v}_{\mathrm{data},i}}{{\sigma }_{v,i}}\right)}^{2},\end{eqnarray} \tag{ 15 }$$

where v is the velocity, ${\sigma }_{v,i}$ is the error in the measurement of the data point, N is the number of data points, and n is the number of parameters being fit. We performed the fit by finding the optimum ${\chi }^{2}$ value for rotation curves produced with a given set of parameters for the disk, bulge, and halo models. Reduced ${\chi }^{2}$ values around 1.0 represent a reasonable fit to the data, implying that most of the model point values are within the errors of the data point values. The parameters for the log (parameter values adopted from Law et al. 2005; see Section 3), NFW (parameter values adopted from Xue et al. 2008; see Section 3), and caustic ring halo were not fit.

In this case, the ${\chi }^{2}$ value now depends only on the disk and bulge parameters: ${\chi }^{2}({M}_{\mathrm{disk}},\ a,\ b,{M}_{\mathrm{bulge}},\ c)$. The fit parameters were found by parameter sweeps, varying one parameter while leaving the others constant and finding the chi-squared values. This consequentially finds the ${\chi }^{2}$ as a function of one parameter. For example, if we have an arbitrary parameter p, we can find ${\chi }^{2}(p)$. If the plot of ${\chi }^{2}(p)$ is close to a horizontal line, the fit does not depend on this parameter. Parameters that affect the fit produce ${\chi }^{2}(p)$ plots that are close to parabolas with a dip at the minimum chi-squared value. We found that the fit here depends on the disk mass $({M}_{\mathrm{disk}})$, disk scale length $(a)$, bulge mass $({M}_{\mathrm{bulge}})$, and bulge scale length $(c)$, but does not depend on the value for the disk scale height $(b)$. We have ${\chi }^{2}({M}_{\mathrm{disk}},\ a,\ {M}_{\mathrm{bulge}},\ c)$ with four free parameters. The number of data points is N = 121, and the number of parameters is n = 4 (Equation (15)). The fit values we found for the disk and bulge parameters are shown in Table 2.

Table 2.  Rotation Curve Fit Parameters

		${{\rm{\Phi }}}_{\mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{disk}}+$	${{\rm{\Phi }}}_{\mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{disk}}+$	${{\rm{\Phi }}}_{\mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{disk}}+$
Parameter		${{\rm{\Phi }}}_{\mathrm{caustic}\ \mathrm{ring}\ \mathrm{halo}}$	${{\rm{\Phi }}}_{\mathrm{log}\ \mathrm{halo}}$	${{\rm{\Phi }}}_{\mathrm{NFW}\ \mathrm{halo}}$
Disk mass	${M}_{\mathrm{disk}}\ ({M}_{\odot })$	$9.89\times {10}^{10}$	$9.11\times {10}^{10}$	$7.33\times {10}^{10}$
Disk scale length	a (kpc)	5.0	5.0	5.0
Disk scale height	b (kpc)	0.26	0.26	0.26
Bulge mass	${M}_{\mathrm{bulge}}\ ({M}_{\odot })$	$1.89\times {10}^{10}$	$1.89\times {10}^{10}$	$1.89\times {10}^{10}$
Bulge scale radius	c (kpc)	0.33	0.33	0.33
	${\chi }_{v}^{2}$	1.04	1.01	1.19

Download table as:  ASCIITypeset image

In Figure 3, the sum of a caustic ring halo (lower solid line), Hernquist bulge (dotted), Miyamoto–Nagai disk (dashed), HI gas layer (dash-dotted line), and ${{\rm{H}}}_{2}$ gas layer (dash-dotted line) produces a rotation curve (upper solid line) that matches observations. The blue data points are from Sofue (2012). The data from Sofue (2012) are a compilation of previous rotation curve data from the literature, with the velocities scaled to assume a rotation speed of $200\;\mathrm{km}\;{{\rm{s}}}^{-1}$ at 8.0 kpc. The dash-dotted line is a polynomial fit to the gas model (open circles for ${{\rm{H}}}_{2}$ and crosses for HI) from Olling & Merrifield (2000), as discussed above. The caustic ring halo rotation curve maintains kinks in the curve at the position of the caustics, similar to the acceleration plot (Figure 2). Although the data and the model are similar, the observations are not detailed enough, particularly for $r\gt 10$ kpc, to correlate bumps in the rotation curve with rises in the data points, as was suggested as evidence for the existence of caustic rings by Sikivie (2003).

Zoom In Zoom Out Reset image size

Figure 4 shows rotation curves generated with the same disk, bulge, and gas components as Figure 3, except with different halo models. The rotation curve fit for $0\lt r\lt 2$ depends primarily on the contribution of the bulge model. The bulge mass and scale radius determine the fit to the region near the central peak at $r\sim 0.3$. The fit for $2\lt r\lt 15$ depends mostly on the contributions of the disk and halo. The caustic ring model has a negligible amount of mass within 2 kpc of the Galactic center and also does not contribute very much to the rotation curve for the inner Galaxy. This is shown in Figure 4 (top panel), where the disk curve is much higher than the halo curve for $0\lt r\lt 15$. The same is true for the case of a log halo model (middle panel). The NFW halo conversely predicts a much larger amount of mass for the inner Galaxy. In Figure 4, the NFW halo curve is almost level with the disk model curve. This exposes a degeneracy in the choice of the disk or halo model to fit the rotation curve data. Here the choice for a caustic ring halo or log halo is consistent with the maximum disk hypothesis (discussed in Binney & Tremaine 2008). This is the idea that the disk should contribute most (50%–90% of the total radial force) to the rotation curve of the inner Galaxy, and it assumes that the mass-to-light ratio of the disk does not depend on the radius.

Zoom In Zoom Out Reset image size

The fit values for the disk and bulge parameters are shown in Table 2. The main difference in the fits for the three panels in Figure 4 is the value of the disk mass $({M}_{\mathrm{disk}})$: 9.89 $\times {10}^{10}$, 9.11 $\times {10}^{10}$, and 7.33 $\times {10}^{10}$ ${M}_{\odot }$. In Table 2, a, b, ${M}_{\mathrm{bulge}}$, and c are identical for each column. Estimates for the Milky Way disk and bulge masses vary depending on which particular model and technique are used to fit the observational data. Dehnen & Binney (1998) used a Galaxy model with a spheroidal bulge, exponential disk, and spheroidal halo to fit values of ${M}_{\mathrm{disk}}\approx 5.13\times {10}^{10}\ {M}_{\odot }$ and ${M}_{\mathrm{bulge}}\approx 0.518\times {10}^{10}\ {M}_{\odot }$ for the disk and bulge. Widrow & Dubinski (2005) fit values of ${M}_{\mathrm{disk}}\approx 4.58\times {10}^{10}\ {M}_{\odot }$ and ${M}_{\mathrm{bulge}}\approx 1.1\times {10}^{10}\ {M}_{\odot }$ for a Galaxy model with a central black hole, Hernquist bulge, exponential disk, and NFW halo. McMillan (2011) fit a value of ${M}_{\mathrm{disk}}\approx 5.76\times {10}^{10}\ {M}_{\odot }$ assuming a bulge mass of ${M}_{\mathrm{bulge}}=0.897\times {10}^{10}\ {M}_{\odot }$ for a Galaxy model with a Bissantz & Gerhard (2002) bulge, exponential disk, and NFW halo. Kafle et al. (2014) fit values of ${M}_{\mathrm{disk}}\approx {9.5}_{-3.0}^{+2.4}\times {10}^{10}\ {M}_{\odot }$ and ${M}_{\mathrm{bulge}}\approx {0.91}_{-0.38}^{+0.31}\times {10}^{10}\ {M}_{\odot }$ for a Galaxy model with a spheroidal bulge, Miyamoto–Nagai disk, stellar halo, and NFW dark halo. The disk masses we found here are larger than what was found for an exponential disk but within the errors of what was found for a Miyamoto–Nagai disk. Our bulge masses are sightly larger than what was found in the sources listed here.

In this context, tidal streams are collections of stars that have been gravitationally stripped from dwarf galaxies that in turn are orbiting much larger galaxies. As a dwarf galaxy is tidally disrupted, two streams or "tails" of stars are formed. The leading tail moves ahead of the orbit of the bound dwarf galaxy, and the trailing tail follows behind it. The Sagittarius (Sgr) dwarf galaxy (Ibata et al. 1994) tidal stream is the largest and best-studied tidal stream in the Milky Way Galaxy (Majewski et al. 2003, 2004; Johnston et al. 2005; Law et al. 2005; Vivas et al. 2005; Law & Majewski 2010).

The Sgr stream has been used to measure the Milky Way potential and thus constrain the shape of the dark matter halo. The data for Sgr tidal debris have, however, been difficult to fit to halo models. There have been claims for spherical, prolate, and oblate ellipsoid models. Since the Sgr debris is approximately confined to a plane, Majewski et al. (2003) suggested that the potential should be spherical. Vivas et al. (2005) compared radial velocities and distances from the leading tidal tails to numerical simulations for different halo models. They found that prolate and spherical models fit better but still fail to reproduce the entire angular span of the stream. The tilt of the tidal tails with respect to the Sgr orbital plane, however, suggests that the potential is oblate (Johnston et al. 2005). Law et al. (2005) conversely found that the prolate models best fit the leading tail velocities. More recently, Law & Majewski (2010) found a reasonable fit to both tidal tails if they assumed a triaxial halo, with both the major and minor axes in the plane of the Milky Way. Vera-Ciro & Helmi (2013) expanded upon this idea to fit a halo potential that varies as a function of Galactocentric radius. These authors obtain an oblate fit within 10 kpc of the Galactic center and a fit to the triaxial model for outer radii. Due to the controversy in modeling the tidal tails of the Sgr stream, it is worth asking whether substructures such as caustic rings could affect the shape and velocity structure of this tidal debris stream.

As a test of the viability of the caustic ring model, we determine whether Sgr stream observations are consistent with orbits and N-body simulations in a caustic ring halo potential. In Section 5.1, we discuss the orbit-fitting technique used and the results of comparing the model to the observational data. After the orbits are fit, N-body simulations are run with those orbits. We will compare the N-body models to the data in Section 5.2, using three cases for the fixed analytical potential of the Milky Way:

• 1.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{caustic}\ \mathrm{halo}}\end{array}$
• 2.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{NFW}\ \mathrm{halo}}\end{array}$
• 3.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{triaxial}\ \mathrm{halo}}.\end{array}$

The galaxy models used here are those that have been fit to the Milky Way rotation curve in Section 3. The disk, bulge, and halo parameters are from Tables 1 and 2. For the case of the triaxial halo potential, we use the disk and bulge parameters from the log halo fit in Tables 1 and 2.

5.1. Orbit Fitting

An orbit (line through phase space) of a test particle through a particular external potential is specified by initial conditions for position and velocity. We generate orbits using the mkorbit and orbint tools from the NEMO Stellar Dynamics Toolbox (Teuben 1995). The mkorbit tool generates an orbit file given a specified analytical external potential and starting values for position and velocity. The orbint tool integrates the orbit through a specified number of time steps. We use the orblist tool to list the position and velocity values calculated at each time step.

We can constrain the potential of the Milky Way by matching model orbits with observational data from stars in tidal streams. We used a gradient descent method to fit Sgr stream data to orbits using methods described in Section 6 of Willett et al. (2009). The NEMO code generally uses Galactocentric Cartesian coordinates for input of positions and velocities $(x,\ y,\ z,\ {v}_{x},\ {v}_{y},\ {v}_{z})$. For convenience, we will use the Sgr heliocentric coordinate system defined in Majewski et al. (2003) to designate positions $({d}_{\mathrm{sun}},\ {{\rm{\Lambda }}}_{\odot },\ {B}_{\odot })$. In this section and the next, we adopt a distance of ${R}_{\odot }=8.0$ kpc from the Sun to the Galactic center. Now the initial conditions for an orbit are specified by $({d}_{\mathrm{sun}},\ {{\rm{\Lambda }}}_{\odot },\ {B}_{\odot },\ {v}_{x},\ {v}_{y},\ {v}_{z})$.

We use the reduced ${\chi }^{2}$ as a measure of fitness of a model orbit to stream data points:

$$\begin{eqnarray}&&{\chi }_{{v}_{\mathrm{gsr}}}^{2}=\displaystyle \sum _{i}{\left(\displaystyle \frac{{v}_{\mathrm{gsr}\ \mathrm{model},i}-{v}_{\mathrm{gsr}\ \mathrm{data},i}}{{\sigma }_{{v}_{\mathrm{gsr}},i}}\right)}^{2},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\chi }_{B}^{2}=\displaystyle \sum _{i}{\left(\displaystyle \frac{{B}_{\ \mathrm{model},i}-{B}_{\mathrm{data},i}}{{\sigma }_{B,i}}\right)}^{2},\ \mathrm{and}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{N-n-1}({\chi }_{{v}_{\mathrm{gsr}}}^{2}+{\chi }_{B}^{2}),\end{eqnarray} \tag{ 18 }$$

where ${v}_{\mathrm{gsr}}$ is the Galactic standard of rest velocity, σ is the error in the measurement of the data point, N is the number of data points, and n is the number of parameters being fit. Since the orbit is a line through phase space, the best fit is one where the orbit line goes most closely through the set of stream data points. The fit is therefore determined by the starting positions and velocities and the choice of the external Galactic potential model. Here we can therefore compare orbit fits in different external potentials.

Since the current coordinates of the Sgr dwarf galaxy are well known, we fix the starting points for the distance and orbital longitude and fit the starting points for the orbital latitude and velocity $({d}_{\mathrm{sun}},\ {{\rm{\Lambda }}}_{\odot },\ {B}_{\odot },\ {v}_{x},\ {v}_{y},\ {v}_{z})$ $=\;(28.0\ \mathrm{kpc},\ 0\buildrel{\circ}\over{.} 126,\ {B}_{\odot },\ {v}_{x},\ {v}_{y},\ {v}_{z})$. We are fitting four parameters. The orbit was fit to Two Micron All Sky Survey (2MASS) M giant angular position/velocity data from Table 3 of Majewski et al. (2004). Only M giants with ${{\rm{d}}}_{\mathrm{sun}}\gt 13$ kpc and $-7.8\lt {B}_{\odot }\lt 7.8$ were selected for our fits. These points were split into $10^\circ$ bins in latitude ${{\rm{\Lambda }}}_{\odot }$, and points with velocities greater than $3\;\sigma$ of the bin average were removed from the data set. The error of each velocity point was taken to be ${\sigma }_{{v}_{\mathrm{gsr}},i}=6\;\mathrm{km}\;{{\rm{s}}}^{-1}$, as stated in Majewski et al. (2004). The error in latitude, ${\sigma }_{B,i}=3$, is the standard deviation of ${B}_{\odot }$ in our data set. The number of data points is N = 80, and the number of parameters is n = 4. We completed orbit fits for the three cases listed above for the external Galactic potential, which highlights the effects of varying the choice of the halo model.

The gradient descent method takes an initial set of parameters $({d}_{\mathrm{sun}},\ {{\rm{\Lambda }}}_{\odot },\ {B}_{\odot },\ {v}_{x1},\ {v}_{y1},\ {v}_{z1})$, generates an orbit, compares the velocities from this orbit to Sgr M giant data, and calculates a ${\chi }^{2}$ value, ${\chi }_{\mathrm{old}}^{2}$. Next the velocities are adjusted to new values $({d}_{\mathrm{sun}},\ {{\rm{\Lambda }}}_{\odot },\ {B}_{\odot },\ {v}_{x2},\ {v}_{y2},\ {v}_{z2})$ by moving in the direction of the ${\chi }^{2}$ gradient, a new orbit is generated, and ${\chi }_{\mathrm{new}}^{2}$ is calculated. If ${\chi }_{\mathrm{new}}^{2}$ is less than ${\chi }_{\mathrm{old}}^{2}$, the new set of parameters is kept. If ${\chi }_{\mathrm{new}}^{2}$ is greater than ${\chi }_{\mathrm{old}}^{2}$, the old set of parameters is kept. The process is repeated until ${\chi }_{\mathrm{new}}^{2}={\chi }_{\mathrm{old}}^{2}$, signaling that a minimum has been reached. The best-fit parameters, $({d}_{\mathrm{sun}},\ {{\rm{\Lambda }}}_{\odot },\ {B}_{\odot },\ {v}_{x},\ {v}_{y},\ {v}_{z})$ $=\;(28.0\ \mathrm{kpc},\ 0\buildrel{\circ}\over{.} 126,\$ ${B}_{\odot ,\ \mathrm{bestfit}},\ {v}_{x,\ \mathrm{bestfit}},\ {v}_{y,\ \mathrm{bestfit}},\ {v}_{z,\ \mathrm{bestfit}})$, correspond with a minimum ${\chi }^{2}$ value. The errors in the fit parameters are estimated using the method described in Cole et al. (2008). We assume that the shape of the ${\chi }^{2}$ surface near the minimum is Gaussian. The curvature of a surface is found using the Hessian or curvature matrix ${\boldsymbol{H}}$, a matrix containing all possible second derivatives of a function. We calculate the derivatives for the Hessian matrix numerically using Equation (23) from Cole et al. (2008). The inverse of the Hessian matrix gives the error with respect to each best-fit parameter:

$$\begin{eqnarray}&&{\boldsymbol{V}}=\displaystyle \frac{1}{N}{{\boldsymbol{H}}}^{-1}={\sigma }^{2}.\end{eqnarray} \tag{ 19 }$$

The diagonal elements of ${\boldsymbol{V}}$ are the variances, and the off-diagonal elements are the covariances. The variance is normalized by the number of data points N. The variable ${\sigma }_{i}=\sqrt{{V}_{{ii}}/N}$ gives the error estimate with respect to each best-fit parameter.

The starting velocities for the gradient descent were random values chosen from a range $({v}_{x},\ {v}_{y},\ {v}_{z})\;=$ (150–250, −40–20, 150–250). These values were chosen to ensure that the gradient descent algorithm begins close to the orbit of Sgr that is known from previous publications. The gradient descent, however, will not necessarily end with values in these ranges. A disadvantage to using the gradient descent method is that if the starting point is very far from the global minimum, the algorithm may become trapped in a local minimum. We start with these values to ensure the global minimum is found. Our best fits are chosen as the run with the lowest ${\chi }^{2}$ value after 10 separate runs. The results for the fits in the three cases listed above are $({B}_{\odot ,\ \mathrm{bestfit}},\ {v}_{x,\ \mathrm{bestfit}},\ {v}_{y,\ \mathrm{bestfit}},\ {v}_{z,\ \mathrm{bestfit}})\;=$

• 1.  
  $(0\buildrel{\circ}\over{.} 028\pm 0.7,\ 210.2\pm 1.1\;\mathrm{km}\;{{\rm{s}}}^{-1},\$ $-30.8\pm 1.1\;\mathrm{km}\;{{\rm{s}}}^{-1},\$ $212.7\pm 0.7\;\mathrm{km}\;{{\rm{s}}}^{-1})$, $\ {\chi }^{2}=52.2\$ (caustic),
• 2.  
  $(0\buildrel{\circ}\over{.} 622\pm 0.7,\ 222.9\pm 1.4\;\mathrm{km}\;{{\rm{s}}}^{-1},$ $\ -32.8\pm 1.4\;\mathrm{km}\;{{\rm{s}}}^{-1},\$ $208.7\pm 0.8\;\mathrm{km}\;{{\rm{s}}}^{-1})$, $\ {\chi }^{2}=45.8\$ (NFW),
• 3.  
  $(0\buildrel{\circ}\over{.} 441\pm 0.7,\ 224.6\pm 1.4\;\mathrm{km}\;{{\rm{s}}}^{-1},$ $\ -29.4\pm 1.4\;\mathrm{km}\;{{\rm{s}}}^{-1},\$ $198.7\pm 0.9\;\mathrm{km}\;{{\rm{s}}}^{-1})$, $\ {\chi }^{2}=53.7\$ (triaxial).

Figures 5 and 6 show the best-fit orbits to the M giant data (black points) for cases 1–3. The trailing data lie on the backward orbit (solid line), and the leading tail data lie on the forward orbit (dotted line). We have plotted the Galactic standard of rest velocity $({v}_{\mathrm{gsr}}=x\cdot {v}_{x}+y\cdot {v}_{y}+z\cdot {v}_{z})$ versus orbital longitude ${{\rm{\Lambda }}}_{\odot }$ in Figure 5. The angle along the tidal stream, ${{\rm{\Lambda }}}_{\odot }$, is zero at the position of the Sgr dwarf galaxy and is defined to increase in the direction of the trailing debris (Majewski et al. 2003). The ${\chi }^{2}$ values for each fit are on the order 50 because of the outliers in the data from ${{\rm{\Lambda }}}_{\odot }=0^\circ -50^\circ$ and the spread in range of the leading tail velocity points. The two data points on the edge of the leading tail are about $50-100\;\mathrm{km}\;{{\rm{s}}}^{-1}$ away from the orbit fit lines, yet from visual inspection of the plot the fits in the three panels are reasonably close to the data points. The trailing tail orbits in all three potentials are similar. In the leading tail, the orbits in the NFW and triaxial halos favor upper data points. Yet the caustic ring halo orbit moves through the center of the points. Since the ${\chi }^{2}$ values for all three fits are similar, perhaps more data are needed in the trailing arm to determine which fit is favored. We have also plotted the orbital latitude, ${B}_{\odot }$, versus orbital longitude, ${{\rm{\Lambda }}}_{\odot }$, in Figure 6. Here, ${B}_{\odot }\approx 0$ is in the plane of the Sgr dwarf galaxy orbit. The orbit fits in each panel are all well confined to the Sgr plane. The N-body simulations (red points) in Figures 5 and 6 are discussed in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We estimate the radial period of the best-fit orbits by calculating the time difference in going from the first backward orbit apocenter to the first forward orbit apocenter. We found the radial periods to be 1.09 Gyr (caustic), 1.18 Gyr (NFW), and 0.92 Gyr (triaxial). These values are similar to those found in Law et al. (2005): 0.85 Gyr, 0.88 Gyr, and 0.87 Gyr for oblate, spherical, and prolate potentials, respectively.

5.2. N-body Simulations of the Sagittarius Dwarf Tidal Stream

In this section, we will compare the results of N-body simulations of the tidal disruption of the Sgr dwarf galaxy. We use the gyrfalcON tool (Dehnen 2000, 2001, 2002) in the NEMO Stellar Dynamics Toolbox (Teuben 1995) to run the N-body simulations. Dwarf galaxies were modeled by a Plummer sphere (Plummer 1911):

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{P}}}=-\displaystyle \frac{{{GM}}_{{\rm{P}}}}{\sqrt{{r}^{2}+{r}_{o}^{2}}}\end{eqnarray} \tag{ 20 }$$

with initial mass ${M}_{{\rm{P}}}={M}_{\mathrm{Sgr},0}=1.8\times {10}^{9}\ {M}_{\odot }$ and scale length ${r}_{o}=1.0$ kpc with the number of bodies $N={10}^{4}$. We assume the initial mass is the sum of the light and dark matter components of the dwarf galaxy. The caustic ring model does not specifically address the dark matter distribution in dwarf galaxies, but it is reasonable to assume that, like the Milky Way, they would not include a central cusp as in an NFW profile. We use ${M}_{\mathrm{Sgr},0}\sim 1.8\times {10}^{9}\ {M}_{\odot }$, which is slightly larger than the value of ${M}_{\mathrm{Sgr},0}\sim 6.4\times {10}^{8}\ {M}_{\odot }$ used by Law & Majewski (2010) but smaller than ${10}^{10}-{10}^{11}\ {M}_{\odot }$ found by Niederste-Ostholt et al. (2010) and used by Purcell et al. (2011). The exact initial mass of the Sgr dwarf galaxy is not important to the results of our simulation.

The Plummer model is placed at a point along the backward orbit (in the past) and evolved to the position of the dwarf galaxy at the present day. The model is evolved in an external Galactic potential corresponding to the three cases listed in Section 5.1. The fixed parameters for the potential are also the same as listed above. The simulations in Figures 5–8 were run for a time equal to three orbital periods in the respective potentials: 3.28 Gyr (caustic), 3.54 Gyr (NFW), and 2.77 Gyr (triaxial).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures 5–8, the Plummer model (red points) is evolved along the best-fit orbits determined in the previous section. The panels show the results in the three halo potentials: caustic (top), NFW (middle), and triaxial (bottom) as a function of the angle along the tidal stream, ${{\rm{\Lambda }}}_{\odot }$. Figure 5 shows ${v}_{\mathrm{gsr}}$ versus ${{\rm{\Lambda }}}_{\odot }$. The caustic ring halo plot has much lower velocities at the end of the leading (dashed line) tidal tail, and one can even imagine that the two M giant stars are each from a different apparent branch of the N-body simulation. The N-body models in all three halos match the trailing tail data well. The caustic and NFW debris match the leading tail data well, whereas the triaxial debris matches five of the data points and misses the sixth one. More data points, however, would be needed to determine which fit is preferred in the leading arm. Figure 6 shows ${B}_{\odot }$ versus ${{\rm{\Lambda }}}_{\odot }$. The N-body models in all three potentials are well confined to the Sgr orbital plane.

In Figures 7 and 8 we compare the N-body simulations to a set of more recent Sgr data that extends farther along the leading and trailing tidal tails, in the range $150\lt {{\rm{\Lambda }}}_{\odot }\lt 250$. Figure 7 shows the simulations overplotted with recent velocity data for giant stars detected photometrically from SDSS by Belokurov et al. (2014). Here the caustic and NFW debris match the trailing data (blue points), but the triaxial debris does not. Conversely, the triaxial debris matches the leading tail data (green points), but the NFW and caustic debris do not for the points from ${{\rm{\Lambda }}}_{\odot }=200^\circ -250^\circ$. Figure 8 is plotted using the Sgr Galactocentric Cartesian coordinate system defined in Majewski et al. (2003). The center of the Sgr Galactocentric system is the point in the Sgr dwarf galaxy plane that is closest to the Galactic center. Because the orbit of the Sgr dwarf galaxy is nearly polar, $({X}_{\mathrm{sgr},\mathrm{gc}},{Y}_{\mathrm{sgr},\mathrm{gc}},{Z}_{\mathrm{sgr},\mathrm{gc}})$ are approximately in the directions (X, −Z, Y), where X points from the Sun to the Galactic center, Y is in the direction of the Sun's motion, and Z is in the direction of the north Galactic pole. The simulation in Figure 8 is overplotted with recent Sgr distance observations from Newberg et al. (2003) and Newby et al. (2013). Newberg et al. (2003) measured the distances of BHB stars detected photometrically from SDSS. Newby et al. (2013) measured Sgr stream center distances by fitting the spatial density of F turnoff stars in SDSS using the statistical photometric parallax technique. The caustic and NFW N-body models extend as far as 90 kpc from the Galactic center, much farther than the triaxial N-body model. This is plausible because the caustic ring halo potential well is shallow (see Figure 2), and the bodies can maintain enough kinetic energy to fly out to farther distances. The caustic and NFW debris, however, do not match the data for the leading arm positions. The triaxial debris is the best fit to the leading arm. We find here that none of the models is able to match both the leading and trailing tail data simultaneously.

The orbital precession of the stream is found by calculating the change in the angle ${{\rm{\Lambda }}}_{\odot }$ between the leading and trailing apocenters. Since the precession depends on gravitational potential, it is worth observing how it varies in the three potentials we compare in this section. We fit second-order polynomials to the N-body model debris and then found the angle between the apocenters. The precesssion values are $\delta {{\rm{\Lambda }}}_{\odot }=103\buildrel{\circ}\over{.} 6$ (caustic), $\delta {{\rm{\Lambda }}}_{\odot }=124\buildrel{\circ}\over{.} 2$ (NFW), and $\delta {{\rm{\Lambda }}}_{\odot }=132\buildrel{\circ}\over{.} 9$ (triaxial). In the Galactocentric coordinate system (Majewski et al. 2003), the values are $\delta {{\rm{\Lambda }}}_{\mathrm{GC}}=100\buildrel{\circ}\over{.} 8$ (caustic), $\delta {{\rm{\Lambda }}}_{\mathrm{GC}}=114\buildrel{\circ}\over{.} 2$ (NFW), and $\delta {{\rm{\Lambda }}}_{\mathrm{GC}}=123\buildrel{\circ}\over{.} 1$ (triaxial). Belokurov et al. (2014) found values of $\delta {{\rm{\Lambda }}}_{\odot }=99.3\pm 3\buildrel{\circ}\over{.} 5$ and $\delta {{\rm{\Lambda }}}_{\mathrm{GC}}=93\buildrel{\circ}\over{.} 2\pm 3\buildrel{\circ}\over{.} 5$ from detections of BHB stars from SDSS in the leading and trailing tails. The caustic halo values are closest to those found by Belokurov et al. (2014). Belokurov et al. (2014) show that for spherical potentials the precession angle increases as the rotation curve gets steeper. We can assume the caustic ring halo potential is approximately spherical since the curves along the radial and z directions have the same scale and general shape (see Figure 2). By inspection of Figure 4, it makes sense that the precession angle in the NFW halo is larger than that in the caustic ring halo since the NFW rotation curve is steeper.

In summary, we show that N-body simulations of the Sgr stream in a caustic ring halo match most of the recent data for stream velocities and distances except on the edge of the leading tail at ${{\rm{\Lambda }}}_{\odot }=200^\circ -250^\circ$. Neither of the caustic, NFW, or triaxial models is able to match velocity/distance data from the leading and trailing tails simultaneously. Lastly, the caustic halo N-body model can match Sgr distance observations in the trailing tidal tail, including debris out to 90 kpc.

We ran N-body simulations of dwarf galaxies represented by Plummer sphere models on orbits that pass close to caustic rings, to see whether the high-density areas would have a strong effect on tidal disruption. We found that although there is a large discontinuity in the acceleration near a caustic, the caustic ring gravitational field does not qualitatively affect the results of the simulation. There are, however, differences in the shapes of the orbits and resulting tidal tails in a caustic ring halo potential compared to the spherical log halo potential. In Figures 9–11, a Plummer model (Equation (20)) dwarf galaxy with ${M}_{P}={10}^{8}\ {M}_{\odot }$, ${r}_{o}=0.5$ kpc, and number of bodies $N={10}^{4}$ is evolved in two different Milky Way gravitational potentials:

• 1.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{caustic}\ \mathrm{halo}}\end{array}$
• 2.  
  $\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{Galaxy}} & = & {{\rm{\Phi }}}_{\mathrm{Hernquist}\ \mathrm{bulge}}+{{\rm{\Phi }}}_{\mathrm{Miyamoto}-\mathrm{Nagai}\ \mathrm{disk}}\\ & & +{{\rm{\Phi }}}_{\mathrm{log}\ \mathrm{halo}}.\end{array}$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The parameters for these models are listed in Tables 1 and 2.

In Figures 9 and 10, the Plummer sphere is evolved along different orbital paths that move near the n = 2 caustic ring at 20.1 kpc from the Galactic center (left boundary of tricusp in Figure 1). The simulations start along the backward orbit (solid line) and evolve to move along the caustic ring and then move along the forward orbit (dotted line). For Figure 9, the orbit is designed to move in the $+Y$ direction through the side of the caustic ring in the plane of the Galaxy at $(X,\ Y,\ Z,\ {v}_{x},\ {v}_{y},\ {v}_{z})$ $\;=\;(20.1\;\mathrm{kpc},\ 0\;\mathrm{kpc},\ 0\;\mathrm{kpc},$ $0\;\mathrm{km}\;{{\rm{s}}}^{-1},$ $200.0\;\mathrm{km}\;{{\rm{s}}}^{-1},\ 0\;\mathrm{km}\;{{\rm{s}}}^{-1})$. Here, X, Y, and Z are Galactocentric Cartesian coordinates with the position of the Sun taken as ${R}_{\odot }=8.0$ kpc or $(X,\ Y,\ Z)=(-8.0,\ 0,\ 0)$ kpc. The simulation time is 2 Gyr. The orbit in a caustic ring halo (top left panel) is approximately circular, whereas the log halo orbit (bottom left panel) has some precession. The dwarf galaxies also travel to different positions after 2.0 Gyr. When we look at histograms of the number of bodies along the Galactic longitude (l is convenient since we are looking in the Galactic plane), we find that the dwarf is less disrupted in a caustic halo. The central peak in a caustic halo has about 50% of the total number of bodies, and the log halo has about 30% of the total. We also ran simulations where the Plummer sphere ended in the same position for the two potentials; those are not shown here. In this case, we found the same behavior where the dwarf galaxy disrupts more slowly in a caustic potential. Here the dwarf galaxy may be on a stable orbit around the ring, and the bodies are stripped from the galaxy core more slowly.

In Figure 10, the orbit is along a polar path and moves in the $+Z$ direction through the side of the caustic at $(X,\ Y,\ Z,\ {v}_{x},\ {v}_{y},\ {v}_{z})\;=$ $(20.1\;\mathrm{kpc},\ 0\;\mathrm{kpc},\ 0\;\mathrm{kpc},$ $0\;\mathrm{km}\;{{\rm{s}}}^{-1},$ $\ 0\;\mathrm{km}\;{{\rm{s}}}^{-1},\ 200.0\;\mathrm{km}\;{{\rm{s}}}^{-1})$, normal to the Galactic plane. The simulation time is 1.1 Gyr. The dwarf galaxy moves to almost the same position after 1.1 Gyr in both the caustic ring and log potentials. In this case, the dwarf disrupts more slowly in the log halo; the central peak in a caustic halo has about 20% of the total number of bodies, and the log halo has about 40% of the total. The caustic ring may have the effect of stripping more bodies from the dwarf galaxy in its passage through the disk.

In Figure 11, the orbit is designed to move in the $+Y$ direction of the Galactic plane just inside and outside of the caustic at $(X,Y,Z,{v}_{x},{v}_{y},{v}_{z})=(19.1\mathrm{kpc},0\mathrm{kpc},$ $0\mathrm{kpc},0\mathrm{km}{{\rm{s}}}^{-1},200.0\mathrm{km}{{\rm{s}}}^{-1},0\mathrm{km}{{\rm{s}}}^{-1})$ and $(X,Y,Z,{v}_{x},{v}_{y},{v}_{z})=(21.1\mathrm{kpc},0\mathrm{kpc},0\mathrm{kpc},0\mathrm{km}{{\rm{s}}}^{-1},200.0\mathrm{km}{{\rm{s}}}^{-1},$ $0\mathrm{km}{{\rm{s}}}^{-1}).$ These simulations were run for 2 Gyr. Here an interesting result is found. If the orbit is started just inside the n = 2 caustic ring, most of the stars are disrupted interior to the rings with distances less than 20 kpc. If we start just outside the ring, the tidal debris is concentrated greater than 20 kpc just outside the ring.

In summary, we find that the caustic rings do have subtle influences on the formation of tidal tails, but the differences are not dramatic. The study in this section is motivated by the proposal that the Monoceros Ring, a ring of stars about 20 kpc from the Galactic center, is associated with a dark matter caustic ring at 20 kpc (Natarajan & Sikivie 2007, citing Newberg et al. 2002). Certainly we designed the orbits in Figures 9–11 to move along a caustic ring to study its effects, but realistically dwarf galaxies would not necessarily fall into the Milky Way with such prescribed orbits. Xu et al. (2015) found that the observed stellar ring structures in the disk are associated with radial waves with peaks above and below the Galactic plane. Xu et al. (2015) find an upward peak or overdensity of stars at 10 kpc, a downward peak at 12–14 kpc, then an upward peak at 16–18 kpc, and a downward peak at 20–24 kpc. These peaks are coincidentally near the positions of the predicted caustic rings at 10, 13, and 20 kpc.

In this paper we have presented the first study of the caustic ring model of the Milky Way halo (Duffy & Sikivie 2008), or "Sikivie halo," in relation to astronomical observations. We show that the caustic ring model is consistent with observations of stars in the Milky Way Galaxy.

We presented the formalism for calculating the gravitational acceleration of a caustic ring halo. Our open-source code for a caustic ring halo is available in the NEMO Stellar Dynamics Toolbox and the Milkyway@home client repository. We use this code to show that the caustic ring halo reproduces a roughly logarithmic halo, with large perturbations near the rings. The caustic ring halo has very little mass within 2 kpc of the Galactic center, resulting in a gravitational acceleration that points away from the Galactic center and toward the innermost rings for distances less than 2 kpc. This outward acceleration may not be realized in practice, due to hierarchical merging and the presence of baryons.

We show that the caustic ring halo can reasonably match the known Galactic rotation curve. We were not able to confirm or rule out an association between the positions of the caustic rings and oscillations in the observed rotation curve, due to insufficient rotation curve data. The parameters for a Galactic model that includes a Hernquist bulge, Miyamoto–Nagai disk, and caustic halo are listed in Table 2. The values for the disk mass, scale length, and scale height and bulge mass and scale length that we found are comparable to those found in other studies that fit Galaxy component models to observations (Dehnen & Binney 1998; Widrow & Dubinski 2005; McMillan 2011; Kafle et al. 2014).

We explore the effects of dark matter caustic rings on dwarf galaxy tidal disruption with N-body simulations. Simulations of the Sagittarius dwarf galaxy in a caustic ring halo potential, with disk and bulge parameters that are tuned to match the Galactic rotation curve, match observations of the Sgr trailing tidal tails as far as 90 kpc from the Galactic center. Like the NFW halo, they are, however, unable to match the leading tidal tail. None of the caustic, NFW, or triaxial logarithmic halos are able to simultaneously match observations of the leading and trailing arms of the Sagittarius stream.

We further show that simulations of dwarf galaxies that move through caustic rings are qualitatively similar to those moving in a logarithmic halo. The dwarf galaxy disrupts more or less slowly depending on whether it moves on a planar or polar orbit in the caustic and logarithmic halo potentials. In simulations of dwarf galaxies that move just inside or outside a caustic ring, the tidal debris is shown to be disrupted mostly to the interior or exterior of the ring, respectively.

We gratefully acknowledge Pierre Sikivie for helpful discussions and clarifications about the caustic ring model. We acknowledge Peter Teuben for assisting with the integration of our code into the NEMO Stellar Dynamics Toolbox. We thank Heywood Tam for sending us his program used to calculate the gravitational field near a caustic ring. We also thank Larry Widrow for insightful discussions.

J.D. acknowledges support from NSF grant AST 10-09670, the NASA-NY Space Grant, the American Fellowship from AAUW, and generous gifts from the Marvin Clan, Babette Josephs, Manit Limlamai, and the MilkyWay@home volunteers.