Estimate of the Mass and Radial Profile of the Orphan–Chenab Stream's Dwarf-galaxy Progenitor Using MilkyWay@home

From the SDSS DR16 release (Ahumada et al. 2020), we selected all stars with R.A. (α) between and 12375 and 1720 and decl. (δ) between −250 and 600. We define our on field to be the remaining stars within 2° of B_corr = 0° and our off field to be the stars with 20 < ∣B_corr∣ < 40. The difference between these two fields in a given Λ bin is defined to be the excess in that bin. To avoid overlap between our two data sets, we also apply a Λ cut to the SDSS data, keeping only stars with Λ < 21°.

We select F-turnoff stars with a ${\left(g-i\right)}_{0}$ color between 0.12 and 0.47. We select this color range to maximize the signal-to-noise ratio between the on-field and off-field stars (see Figure 1). We also filter stars such that 17.0 < g₀ < 22.5. To remove stars in front of or behind the OCS, we enforced a distance cut of 20.7 < g_corr < 21.7, which is the same cut applied in Newberg et al. (2010).

Zoom In Zoom Out Reset image size

The errors in the counts are assumed to follow a Poisson distribution. Thus, the error of the counts in the ith Λ bin N_i , for both the on field and the off field, is simply given by

$$\begin{eqnarray}&&{\sigma }_{{N}_{i}}=\sqrt{{N}_{i}}.\end{eqnarray} \tag{ 3 }$$

To account for the incompleteness in our fields at Λ_OCS < −21°, we must calculate what fraction of each Λ bin is filled. Because the incompleteness results from the g₀ = 22.5 magnitude limit, we can use Newberg et al. (2010) to exactly determine the boundary g_b (Λ_OCS):

$$\begin{eqnarray}&&{g}_{b}({{\rm{\Lambda }}}_{\mathrm{OCS}})={a}_{g}{{\rm{\Lambda }}}_{\mathrm{OCS}}^{2}+{b}_{g}{{\rm{\Lambda }}}_{\mathrm{OCS}}+{c}_{g},\end{eqnarray} \tag{ 4 }$$

where a_g = −0.00022, b_g = 0.034, and c_g = 22.5. As can be seen in Figure 2, there are more F-turnoff stars at higher magnitudes in both the on field and the off field. We therefore correct for the missing stars by fitting a linear model to the star counts as a function of g_corr, as described in Appendix A.

Zoom In Zoom Out Reset image size

The DEC data was taken from Grillmair et al. (2015). Because the DEC uses the same g-band and i-band filters as the SDSS, we applied the same color and magnitude cuts as in the SDSS data, keeping stars with 21° ≤ Λ ≤ 48°. We define our on field and off field in the same way as we did for the SDSS data. However, as can be seen in Figure 3, the DEC data do not completely fill our on and off fields.

Zoom In Zoom Out Reset image size

To correct our fields for the lack of data, we apply Monte Carlo approximations to "fill in" the missing patches in the sky. For the ith Λ bin in the on field, we randomly populate the bin with 8192 (2¹³) test stars. For each test star, we check whether the test star is within 018 of a real star in that bin. We selected 018 as our threshold because this was slightly larger than the maximum nearest-neighbor angular distance between two stars in the DEC data. We count the number of test stars within 018 of a real star, p_i , and divide it by the total number of test stars, M = 8192, to get the "filled" fraction of the ith bin, k_i :

$$\begin{eqnarray}&&{k}_{i}=\displaystyle \frac{{p}_{i}}{M}.\end{eqnarray} \tag{ 5 }$$

Assuming a Poisson distribution for p_i , we derive the following expression for the error in k_i :

$$\begin{eqnarray}&&{\sigma }_{{k}_{i}}=\displaystyle \frac{\sqrt{{p}_{i}}}{M}.\end{eqnarray} \tag{ 6 }$$

By dividing the number of counts in an on-field bin N_i by the ratio k_i , we can approximate the true number of stars in that bin:

$$\begin{eqnarray}&&N{{\prime} }_{i}=\displaystyle \frac{{N}_{i}}{{k}_{i}}.\end{eqnarray} \tag{ 7 }$$

This means that the errors in counts in the Λ bin are given by

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{N{{\prime} }_{i}} & = & \displaystyle \frac{1}{{k}_{i}}\sqrt{{\sigma }_{{N}_{i}}^{2}+{\left(\displaystyle \frac{{N}_{i}}{{k}_{i}}\right)}^{2}{\sigma }_{{k}_{i}}^{2}}\\ & = & \displaystyle \frac{1}{{k}_{i}}\sqrt{{N}_{i}\left(1+\displaystyle \frac{{N}_{i}}{{{Mk}}_{i}}\right)}.\end{array}\end{eqnarray} \tag{ 8 }$$

We apply the same process to each bin in the off field as well.

After combining the stellar data from the SDSS with those from the DEC, we look at each bin and calculate the number of stars within the on field and off field. Using the off field as a background, we subtract it from the on field and determine the excess within the stream. The excess of the ith bin (E_i ) and its error (${\sigma }_{{E}_{i}}$) are calculated using the following formulae:

$$\begin{eqnarray}{E}_{i}=\left\{\begin{array}{cl}{N}_{\mathrm{on},i}-{N}_{\mathrm{off},i} & {N}_{\mathrm{on},i}\geqslant {N}_{\mathrm{off},i}\\ 0 & {N}_{\mathrm{on},i}\lt {N}_{\mathrm{off},i}\end{array}\right.\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\sigma }_{{E}_{i}}=\sqrt{{\sigma }_{{N}_{\mathrm{on},i}}^{2}+{\sigma }_{{N}_{\mathrm{off},i}}^{2}}.\end{eqnarray} \tag{ 10 }$$

The excess and errors of each bin are listed in Table 1 along with the number of stars in each on-field and off-field bin. We also present a histogram representing this data in Figure 4. Using our cuts, we find an excess of 5631 ± 356 F-turnoff stars in the OCS within the range of −33° ≤ Λ ≤ 48°. For the purposes of our likelihood calculation, which we will describe later, we also calculate the normalized excess star count (e_i ) and its error $\left({\sigma }_{{e}_{i}}\right)$:

$$\begin{eqnarray}&&{e}_{i}=\displaystyle \frac{{E}_{i}}{{\sum }_{j}{E}_{j}}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\sigma }_{{e}_{i}}=\displaystyle \frac{1}{{\sum }_{j}{E}_{j}}\sqrt{(1-2{e}_{i}){\sigma }_{{E}_{i}}^{2}+{e}_{i}^{2}\sum _{j}{\sigma }_{{E}_{j}}^{2}}\end{eqnarray} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

The simulations include bodies that represent both baryons (stars) and dark matter. The baryonic simulation bodies represent all of the stellar mass, not just the mass of the turnoff stars. Therefore, to constrain the mass of the OCS, we must relate the number of F-turnoff stars in the sky to an amount of baryonic mass. To accomplish this, we search for a globular cluster whose CMD properties most closely match those of the OCS. Once we find a globular cluster with a turnoff star color that is close enough to that of the OCS, we use that globular cluster to estimate the stellar mass per F-turnoff star in the OCS. Simply multiplying this ratio by the number of F-turnoff stars we detect in the OCS will give us an approximate stellar mass. For our analysis, we look at 11 candidate globular clusters using the globular cluster data from An et al. (2008).

It should be noted that because An et al. (2008) use the dust map from Schlegel et al. (1998) to calculate their unextincted magnitudes and our stellar data calculates extinctions using Schlafly & Finkbeiner (2011), we need to correct the magnitudes from our globular clusters to adjust for overreddening. Before we perform our globular cluster analysis, we add the extinction calculated from Schlegel et al. (1998) to each star's magnitude in An et al. (2008) and then subtract the extinction derived from Schlafly & Finkbeiner (2011).

When selecting globular cluster stars for our F-turnoff color calculation, we must search over a defined magnitude range for each cluster. To select this range, we plot a color–magnitude diagram (CMD) in ${\left(g-i\right)}_{0}$ and g₀. We set the minimum (brightest) turnoff magnitude at the intersection of the subgiants with the red giant branch. We record the ${\left(g-i\right)}_{0}$ color of the intersection and determine the brightest magnitude for main-sequence stars of that color, defining that g₀ as the maximum turnoff magnitude. The ranges are listed Table 2, and the CMDs of each candidate can be seen in Figure 5.

Zoom In Zoom Out Reset image size

Table 2. Color Peaks and the Magnitude Range for Each of Our 11 Globular Clusters

Globular Cluster	g0 F-turnoff Range	${\left(g-r\right)}_{0}$	${\left(g-i\right)}_{0}$
M2	19.0–21.0	0.351 ± 0.015	0.457 ± 0.028
M3	18.5–20.5	0.309 ± 0.022	0.382 ± 0.047
M5	18.0–20.0	0.326 ± 0.024	0.456 ± 0.051
M13	18.0–20.0	0.318 ± 0.035	0.425 ± 0.056
M15	18.75–20.75	0.338 ± 0.020	0.483 ± 0.033
M53	19.75–21.5	0.264 ± 0.034	0.351 ± 0.042
M92	18.0–20.0	0.262 ± 0.027	0.359 ± 0.043
NGC 4147	20.0–21.5	0.267 ± 0.038	0.353 ± 0.050
NGC 5053	19.5–21.5	0.241 ± 0.025	0.326 ± 0.041
NGC 5466	19.25–21.25	0.283 ± 0.024	0.351 ± 0.029
Palomar 5	20.5–22.25	0.361 ± 0.040	0.505 ± 0.053

Download table as:  ASCIITypeset image

Using these stars, we generate two one-dimensional histograms with bins in ${\left(g-r\right)}_{0}$ and ${\left(g-i\right)}_{0}$ color, respectively. We then fit the histogram to a skew-normal distribution with a linear background. The model has the form

$$\begin{eqnarray}&&f(x)={mx}+b+{{Ae}}^{-\tfrac{{\left(x-\xi \right)}^{2}}{2{\omega }^{2}}}\left(1+\mathrm{erf}\left(\gamma \left(\displaystyle \frac{x-\xi }{\omega \sqrt{2}}\right)\right)\right),\end{eqnarray} \tag{ 13 }$$

where m is the linear slope, b is the linear y intercept, A is the amplitude, ξ is the "unskewed" mean, ω is the "unskewed" standard deviation, and γ is the skew parameter. In our fits, we assume that γ ≥ 0 because we expect the number of stars to drop off more slowly toward the redder end of the distribution. We fit these six parameters using a χ² best-fit method and a differential evolution algorithm. We run the algorithm 10 independent times and keep the best fit to ensure our χ² fit does not fall into a local minimum. After fitting the six parameters, we calculate the Hessian of the natural logarithm of the χ² surface around that minimum. By inverting the Hessian, we derive the covariance matrix of the fitted parameters:

$$\begin{eqnarray}\begin{array}{rcl}{\left(H\mathrm{ln}{\chi }^{2}\right)}^{-1} & = & {\left(\begin{array}{ccc}\displaystyle \frac{{\partial }^{2}\mathrm{ln}{\chi }^{2}}{\partial {{x}_{1}}^{2}} & \displaystyle \frac{{\partial }^{2}\mathrm{ln}{\chi }^{2}}{\partial {x}_{1}\partial {x}_{2}} & \cdots \\ \displaystyle \frac{{\partial }^{2}\mathrm{ln}{\chi }^{2}}{\partial {x}_{2}\partial {x}_{1}} & \displaystyle \frac{{\partial }^{2}\mathrm{ln}{\chi }^{2}}{\partial {{x}_{2}}^{2}} & \cdots \\ \vdots & \vdots & \ddots \end{array}\right)}^{-1}\\ & = & \left(\begin{array}{ccc}{{\sigma }_{{x}_{1}}}^{2} & {\sigma }_{{x}_{1}{x}_{2}}^{2} & \cdots \\ {\sigma }_{{x}_{2}{x}_{1}}^{2} & {{\sigma }_{{x}_{2}}}^{2} & \cdots \\ \vdots & \vdots & \ddots \end{array}\right)\end{array}\end{eqnarray} \tag{ 14 }$$

Each second derivative is calculated numerically using a finite step size, where each ith parameter is assigned a step size h_i . We approximate each second derivative using the following formula:

$$\begin{eqnarray}\displaystyle \frac{{\partial }^{2}f}{\partial {x}_{i}\partial {x}_{j}}\simeq \left\{\begin{array}{cl}\displaystyle \frac{f({x}_{i}-{h}_{i},{x}_{j}-{h}_{j})-f({x}_{i}-{h}_{i},{x}_{j}+{h}_{j})-f({x}_{i}+{h}_{i},{x}_{j}-{h}_{j})+f({x}_{i}+{h}_{i},{x}_{j}+{h}_{j})}{4{h}_{i}{h}_{j}} & i\ne j\\ \displaystyle \frac{f({x}_{i}-{h}_{i})-2f({x}_{i})+f({x}_{i}+{h}_{i})}{{{h}_{i}}^{2}} & i=j.\end{array}\right.\end{eqnarray} \tag{ 15 }$$

After calculating the errors of each parameter, we take these errors and set them as the new step size for each of their respective parameters. We then calculate the Hessian again, taking those errors and feeding them back into the Hessian calculation over and over until the errors converge within a significantly small enough threshold (0.0001).

Taking our six fitted parameters and their errors, we calculate the mode of the fitted distribution. While skewed normal distributions do not have an analytical solution for the mode Mo, there does exist a numerical approximation for it while γ ≥ 0:

$$\begin{eqnarray}&&{Mo}\simeq \xi +\omega \left(\zeta \sqrt{\displaystyle \frac{2}{\pi }}-\left(\displaystyle \frac{4-\pi }{4}\right)\displaystyle \frac{{\left(\zeta \sqrt{\tfrac{2}{\pi }}\right)}^{3}}{1-\tfrac{2}{\pi }{\zeta }^{2}}-\displaystyle \frac{1}{2}{e}^{-\tfrac{2\pi }{\gamma }}\right),\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&\zeta =\displaystyle \frac{\gamma }{\sqrt{1+{\gamma }^{2}}}.\end{eqnarray} \tag{ 17 }$$

Running our SDSS OCS data through this pipeline, we find that the excess of the OCS has a peak ${\left(g-r\right)}_{0}$ color of 0.214 ± 0.066 and a ${\left(g-i\right)}_{0}$ color of 0.286 ± 0.083. We similarly run our 11 candidate globular clusters through our algorithm to determine which has the closest color peak. Table 2 lists the color peaks for each of these globular clusters.

The globular cluster whose peak F-turnoff star color is most similar to that of the OCS is NGC 5053.

After determining which globular cluster to use, we count the total number of F-turnoff stars we see in NGC 5053 from the photometric data provided by An et al. (2008), applying the same color cuts we employed in our stream separation and selecting stars with $0.12\lt {\left(g-i\right)}_{0}\lt 0.47$ and 19.5 < g₀ <21.5 as our turnoff stars. However, because globular clusters have a high stellar density in the sky, it is difficult to resolve all the stars within the center of the cluster. Although An et al. (2008) use a DAOPHOT pipeline that is specifically designed for crowded-field stellar photometry (Stetson 1987) to parse the SDSS data, we find that there are still stars that are not resolved. To correct for the missing stars, we fit the radial distribution of NGC 5053 to that of a Plummer sphere. The 2D surface number density of a Plummer sphere is given by

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=\displaystyle \frac{{N}_{\mathrm{total}}{a}^{2}}{\pi {\left({R}^{2}+{a}^{2}\right)}^{2}},\end{eqnarray} \tag{ 18 }$$

where the number of stars ΔN within ΔR of radius R from the center of the cluster is represented by the formula:

$$\begin{eqnarray}&&{\rm{\Delta }}N=2\pi R{\rm{\Sigma }}(R){\rm{\Delta }}R.\end{eqnarray} \tag{ 19 }$$

We fit the total number of F-turnoff stars (N_total) and the angular scale radius (a) to NGC 5053 using the method of least squares, setting the radius R for a given star as the angular distance between that star and the center of the cluster, which we find to be at (α, δ) = (199107, 176927). To ensure our fit is not heavily impacted by the overcrowding near the center of the globular cluster, we exclude all stars within the central core by requiring R > θ_1/2 ∼ 0035, where θ_1/2 is the angular distance from the globular cluster's center where the stellar density drops to half of its maximum value. Figure 6 shows the radial distribution of F-turnoff stars in NGC 5053 and the fitted Plummer curve.

Zoom In Zoom Out Reset image size

Correcting for missing stars near the center of the cluster, we find a total of 3932 ± 110 F-turnoff stars within NGC 5053. Taking the mass of NGC 5053 to be (5.37 ± 1.32) × 10⁴ M_⊙ (Baumgardt 2017), we calculate that each turnoff star represents a stellar mass of roughly (13.7 ± 3.4)M_⊙. A calculation of the mass per turnoff star using isochrones, which finds a similar value, is done in Appendix B. Multiplying this number by the 5631 ± 356 F-turnoff stars we calculate in the OCS gives us a total baryonic mass of (7.7 ± 2.0) × 10⁴ M_⊙ within the range −330 < Λ < 480.

2.4.1. Systematic Error from Using a Redder Globular Cluster

While the turnoff color of NGC 5053 matches that of the OCS within errors, the globular cluster does have a somewhat redder measure of turnoff color. To illustrate how a redder color would impact our estimate of the mass per turnoff star, we recalculate this quantity using all of our other globular clusters as the reference. Table 3 below shows the mass per F-turnoff star we find if we use each of these redder globular clusters as a reference.

Table 3. Mass per F-Turnoff Star Using Different Globular Clusters as a Reference

Globular Cluster	${\left(g-r\right)}_{0}$	${\left(g-i\right)}_{0}$	Total Mass (M⊙)	Num. F Stars	Mass/Turnoff Star (M⊙)
M2	0.351 ± 0.015	0.457 ± 0.028	(5.75 ± 2.65) × 105	16,701 ± 1687	34.4 ± 15.9
M3	0.309 ± 0.022	0.382 ± 0.047	(4.68 ± 1.40) × 105	15,808 ± 261	29.6 ± 8.9
M5	0.326 ± 0.024	0.456 ± 0.051	(3.89 ± 1.16) × 105	12,061 ± 376	32.3 ± 9.7
M13	0.318 ± 0.024	0.425 ± 0.056	(6.05 ± 0.33) × 105	17,464 ± 554	34.6 ± 2.2
M15	0.338 ± 0.020	0.483 ± 0.033	(5.13 ± 1.06) × 105	6543 ± 1843	78.4 ± 27.4
M53	0.264 ± 0.034	0.351 ± 0.042	(3.48 ± 1.04) × 105	10,913 ± 328	31.9 ± 9.6
M92	0.262 ± 0.027	0.359 ± 0.043	(2.75 ± 0.57) × 105	11,382 ± 1107	24.2 ± 5.5
NGC 4147	0.267 ± 0.038	0.353 ± 0.050	(3.72 ± 2.82) × 104	1564 ± 87	23.8 ± 18.1
NGC 5053	0.241 ± 0.025	0.326 ± 0.041	(5.37 ± 1.32) × 104	3932 ± 110	13.7 ± 3.4
NGC 5466	0.283 ± 0.024	0.351 ± 0.029	(3.80 ± 2.45) × 104	4425 ± 126	8.6 ± 5.5
Palomar 5	0.361 ± 0.040	0.505 ± 0.053	(5.2 ± 0.7) × 103	528 ± 37	9.8 ± 1.5

Note. The total mass of each globular cluster, with the exception of M13, NGC 5053, and Palomar 5, came from Kimmig et al. (2015). M13's total mass was pulled from Leonard et al. (1992) and Palomar 5's mass came from Odenkirchen et al. (2002).

Download table as:  ASCIITypeset image

Plotting the mass per turnoff star as a function of the peak color (Figure 7), we see that redder globular clusters typically produce a larger calculated mass per turnoff star. Fitting these points to an exponential curve of the form $\mathrm{ln}(M/{M}_{\odot })=r{\left(g-i\right)}_{0}+\mathrm{ln}(A)$, we find r = 17.7 ± 4.6 and $\mathrm{ln}(A)=-3.96\pm 1.89$ (χ² = 14.8). This implies that using a globular cluster for the mass estimate of a stellar system with a bluer peak color overestimates its mass. Because we select a fixed (blue) color range to select our turnoff stars $(0.12\lt {\left(g-i\right)}_{0}\lt 0.47)$, it makes sense that using a redder globular cluster would result in fewer stars falling within this color range and thus give us a larger stellar mass per turnoff star. If we extrapolate the curve to find the mass per turnoff star using a globular cluster with a color similar to the OCS, we find 3.0 ± 6.9M_⊙ per turnoff star, roughly one-quarter of the mass per turnoff star estimated using NGC 5053. It should be noted, however, that we do not use this extrapolated mass per turnoff for analyzing the OCS due to the extrapolation's poor fit and large errors. We instead use the mass per turnoff derived from NGC 5053, keeping in mind that any masses we fit using this data are overestimates.

Zoom In Zoom Out Reset image size

To fit the apparent width of the stream as a function of Λ, we measure the Beta dispersion of the OCS by splitting each Λ bin into 20 B bins, covering both the on and off fields. In each Λ bin, we fit the distribution of B values to a Gaussian distribution with a linear background, taking the fitted σ value as the Beta dispersion:

$$\begin{eqnarray}&&f(x)={mx}+b+{{Ae}}^{-\tfrac{{\left(x-\mu \right)}^{2}}{2{\sigma }^{2}}},\end{eqnarray} \tag{ 20 }$$

where m is the slope of the background, b is the y intercept, A is the amplitude, μ is the mean, and σ is the standard deviation of the Gaussian distribution.

To ensure a larger sample size, we group each consecutive set of three Λ bins into a single bin for the purposes of the Beta dispersion calculation. We employ the same Monte Carlo method of bin completion used in Section 2.2 to determine the number of stars in each Beta bin. However, we exclude any bins whose data coverage is less than half of the bin area from our fit. We use the "curve fit" algorithm from SciPy (Virtanen et al. 2020) to fit the beta dispersions and their respective errors. Figure 8 shows the fit in each set of Λ bins while Table 4 lists the numerical values of the Beta dispersion and their errors.

Zoom In Zoom Out Reset image size

We represent our simulated dwarf galaxy as a collection of 40,000 bodies: half of the bodies represent visible baryonic matter, and the other half represent dark matter. Regardless of the baryonic and dark matter masses, we employ the same number of bodies, changing the mass per particle to represent the full mass of a component. We choose 40,000 bodies because it is the minimum number of particles needed to consistently generate the same stream over different random seeds (see Figure 9).

Zoom In Zoom Out Reset image size

To model our progenitor dwarf galaxy, we arrange our 40,000 bodies into two nested, concentric Plummer spheres (Plummer 1911). The radial profile of a Plummer sphere is given by

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{3M}{4\pi {a}^{3}}{\left(1+\displaystyle \frac{{r}^{2}}{{a}^{2}}\right)}^{-5/2},\end{eqnarray} \tag{ 21 }$$

where a is a scale radius and M is the total mass of the Plummer sphere. Each progenitor can thus be described by four parameters, the baryonic mass (M_B ), the dark matter mass (M_D ), the baryon scale radius (a_B ), and the dark matter scale radius (a_D ). However, because the properties of our dwarf galaxies have a degree of scale invariance, it is more efficient to optimize over the total baryonic component (M_B and a_B ) and two ratios (ξ_R and ξ_M ). These ratios are defined as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\xi }_{R} & = & \displaystyle \frac{{a}_{B}}{{a}_{B}+{a}_{D}}\\ {\xi }_{M} & = & \displaystyle \frac{{M}_{B}}{{M}_{B}+{M}_{D}}.\end{array}\end{eqnarray} \tag{ 22 }$$

The method we use to randomly generate these progenitor dwarf galaxies is detailed in Shelton et al. (2021).

We determine the starting location of the progenitor by first specifying its present-day position and velocity, which we calculated using an orbit for the OCS determined from Newberg et al. (2010). Using its present-day position and velocity, we integrate a single-body orbit backward in time up to the evolution time τ_evolve to find the progenitor's starting position and velocity, centering the 40,000 body simulated progenitor at that point. Although the orbit remains fixed across simulations, we determined it from information different from what we will use to measure the dwarf-galaxy parameters; the orbit was fit using sky position, line-of-sight velocity, and heliocentric distance.

Our progenitor dwarf galaxy is thus completely characterized with five parameters: the evolution time τ_evolve, the baryonic Plummer radius a_B , the radius ratio ξ_R , the baryonic Plummer mass M_B , and the baryonic mass ratio ξ_M . As such, these are the parameters we wish to optimize over using MilkyWay@home.

Our Milky Way potential consists of three components: a Miyamoto‒Nagai disk, a Hernquist bulge, and a logarithmic halo (Miyamoto & Nagai 1975; Hernquist 1990; Sakamoto et al. 2003). Our spherical Hernquist bulge has a scale radius of 0.7 kpc and a mass of 9.9 ×10¹⁰ M_⊙ (Law et al. 2005). The Miyamoto‒Nagai disk in our simulations has a scale radius of 6.5 kpc, a scale height of 0.26 kpc, and a mass of 3.4 ×10¹⁰ M_⊙ (Law et al. 2005). We implement a logarithmic halo with a circular velocity of 73 km s⁻¹ (74.61 kpc Gyr⁻¹) (Newberg et al. 2010) and a scale radius of 12 kpc (Law et al. 2005) to model the dark matter halo. We translate these potentials into accelerations by calculating the negative gradient of each potential ( a = − ∇ϕ). Within our simulations, we represent distances in kiloparsecs (kpc), time in gigayears (Gyr), and mass in structure masses (SM), where 1 SM = 222,288.47 M_⊙. Using these units, G = 1.

Our algorithm uses classical Newtonian gravitation to calculate the acceleration between two bodies. We employ a Barnes‒Hut tree algorithm (Barnes & Hut 1986), which has a complexity of $O(N\mathrm{log}N)$ to speed up calculations. Our N-body simulations employ a Velocity Verlet method (Verlet 1967) to calculate the new positions and velocities of each particle at the end of each time step. We use this Velocity Verlet algorithm because it is a symplectic integrator, meaning the area of its phase space is conserved as the system evolves through time. This property is especially attractive for our integrator because it implies that energy is conserved as the dwarf galaxy falls into the Milky Way, an important constraint that must be maintained throughout our N-body simulations. The calculation of our accelerations, time step, and softening length is further detailed in Shelton et al. (2021).

After establishing the Milky Way potential and creating the progenitor dwarf, we determine the starting point of our simulation by defining approximately where we expect the progenitor's core to lie in the stream at the present time and thus where we would like the simulated progenitor core to end up at the end of the simulation. For the OCS, we place this starting point at (l, b, R) = (2580, 458, 21.5 kpc) with a velocity (v_x , v_y , v_z ) = (−185.5, 54.7, 147.4) kpc Gyr⁻¹. This velocity is given in right-handed Galactic Cartesian coordinates, in which the X-axis is the direction from the Sun to the Galactic Center and the Y-axis is in the direction of the Sun's motion. This starting point defines the orbit of the OCS and was calculated using the orbit fit from Newberg et al. (2010). We place a single body at the starting point and evolve that body backward in time by our evolution time, τ_evolve. We replace the body with the simulated dwarf-galaxy progenitor, and the N-body integration begins.

After running a simulation to completion, we need to translate the resulting tidal stream into a histogram to compare with the data histogram. We transform our bodies into a Lambda−Beta (Λ, B) system that follows the OCS along its B = 0 great circle. We use the same coordinate system as derived in Newberg et al. (2010). We split the stream into a histogram that consists of 28 Λ bins subtending the range (−33°, 48°) and one Beta bin covering the range (−15°, 15°). Within each bin, we record the number of bodies and normalize the counts with respect to the total number of bodies within the histogram's range. We also calculate the Beta dispersion at each Λ where we have a measured stream width.

To compare the similarity between two histograms, we employ a metric called the likelihood score. This value represents the natural logarithm of the probability that two histograms match. As such, the likelihood score ranges from negative infinity (worst case set to −9999999.9) to 0, where a higher likelihood is indicative of a closer match. Our likelihood score consists of three components: the Earth Mover Distance (EMD) Component, the Cost Component, and the Beta Dispersion Component.

The EMD Component (Rubner et al. 2000) measures how well the shapes of two normalized histograms match each other. Given two histograms, one representing the number of stream stars as a function of Λ and the other representing the number of simulation bodies as a function of Λ, the code calculates the minimum amount of "work" necessary to deform one histogram into the other and translates that work into an EMD score. If the EMD score is greater than some maximum EMD score we define (${\mathrm{EMD}}_{\max }=50$), the code outputs the worst-case likelihood score. Otherwise, we calculate the EMD component to be

$$\begin{eqnarray}&&\mathrm{ln}({{ \mathcal L }}_{\mathrm{EMD}})=300\mathrm{ln}\left(1-\frac{\mathrm{EMD}}{{\mathrm{EMD}}_{\max }}\right).\end{eqnarray} \tag{ 23 }$$

The Cost Component compares the total mass within each histogram and adds a penalty commensurate to the mass difference. We utilize this component because the EMD score can only be calculated if both histograms are normalized, so the information about their masses is lost. The Cost Component is given by

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}({{ \mathcal L }}_{\mathrm{Cost}})\\ \ =\ -\displaystyle \frac{{\left({M}_{\mathrm{sim}}{N}_{\mathrm{sim}}-{M}_{\mathrm{data}}{N}_{\mathrm{data}}\right)}^{2}}{2\left({M}_{\mathrm{data}}^{2}{N}_{\mathrm{data}}+{M}_{\mathrm{sim}}^{2}{N}_{\mathrm{sim}}\left(\tfrac{{N}_{\mathrm{sim}}}{{N}_{\mathrm{total}}}\right)\left(1-\tfrac{{N}_{\mathrm{sim}}}{{N}_{\mathrm{total}}}\right)\right)},\end{array}\end{eqnarray} \tag{ 24 }$$

where M_sim is the mass per body in the simulation, M_data is the mass per turnoff star we calculated in Section 2.4, N_sim is the number of baryonic bodies within the histogram's range, N_data is the number of turnoff stars within the histogram's range (5631 stars), and N_total is the total number of baryonic bodies we ran in the simulation (20,000).

The Beta Dispersion Component compares the width of the streams as a function of Λ. For each Λ bin in the simulation, we calculate the Beta dispersion ${\sigma }_{\beta ,\mathrm{sim},i}$ and its error ${\delta }_{{\sigma }_{\beta ,\mathrm{sim},i}}$, and then compare them to the values determined from the stellar data, σ_β,data,i and ${\delta }_{{\sigma }_{\beta ,\mathrm{data},i}}$. The formula for the Beta Dispersion Component is given by

$$\begin{eqnarray}&&\mathrm{ln}({{ \mathcal L }}_{\mathrm{Disp}})=-\displaystyle \frac{1}{2}\sum _{i}^{{N}_{\mathrm{bins}}}\displaystyle \frac{{\left({\sigma }_{\beta ,\mathrm{sim},i}-{\sigma }_{\beta ,\mathrm{data},i}\right)}^{2}}{{\delta }_{{\sigma }_{\beta ,\mathrm{sim},i}}^{2}+{\delta }_{{\sigma }_{\beta ,\mathrm{data},i}}^{2}}.\end{eqnarray} \tag{ 25 }$$

The total likelihood score is the sum of each of these components:

$$\begin{eqnarray}&&\mathrm{ln}({ \mathcal L })=\mathrm{ln}({{ \mathcal L }}_{\mathrm{EMD}})+\mathrm{ln}({{ \mathcal L }}_{\mathrm{Cost}})+\mathrm{ln}({{ \mathcal L }}_{\mathrm{Disp}})\end{eqnarray} \tag{ 26 }$$

The details of how each component is calculated can be found in Shelton et al. (2021). Only the locations of the baryonic bodies contribute to the calculation of the likelihood score because we are only able to observe the stars in the tidal stream.

Due to the chaotic behavior of N-body systems, the slightest deviation in the progenitor parameters (or the random seed) can shift the core of our progenitor further ahead or behind in the stream, drastically affecting the likelihood. To mitigate this chaos, we check the likelihood score over a range of evolution times to "match" the cores of the simulated stream to those of the data. Once the simulation reaches 98% of the evolution time, the code calculates the likelihood at each time step and keeps the snapshot with the best likelihood. This comparison continues until the simulation time reaches 102% of the evolution time. The best likelihood in this range of forward evolution times is returned as the likelihood for the given parameter set.

MilkyWay@home starts an optimization by first generating a random population of 50 parameter sets, each containing the five parameters described in Section 3.1, on the main server. For each parameter set in the population, MilkyWay@home sends out a package containing an N-body executable, a Lua file containing important settings for the simulation, the aforementioned parameter set, and a comparison histogram containing the real stellar data to several of the 26,000 volunteered computers. The volunteered computer runs the N-body executable with the chosen parameters and generates a simulated tidal stream, calculates the density of bodies across the stream, and logs the information into a histogram. The newly generated histogram is compared to the data histogram, and the computer generates a likelihood score based on how well the two histograms match. This score and its associated parameter set are sent back to the main server.

After calculating the likelihood scores of the initial population, we employ a differential evolution algorithm to generate new parameter sets to test. For each member of our population x, we select three other random distinct members (a, b, and c) to act as "genetic donors" to the new parameter set y. We then generate a random integer p between 1 and n, where n is the number of elements in our parameter set (in this case 5). Also, we generate a random number r_i between 0 and 1 for each element x_i . Given these numbers, the ith element y_i of the new parameter set is thus defined as

$$\begin{eqnarray}{y}_{i}=\left\{\begin{array}{cl}{a}_{i}+F({b}_{i}-{c}_{i}) & {r}_{i}\lt {CR}\vee i=p\\ {x}_{i} & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 27 }$$

where CR is the crossover rate (set to 0.9) and F is the differential weight (set to 0.8). These newly created sets are then sent to our volunteers to calculate a likelihood score. Once a parameter set returns with its likelihood score, it is compared against the rest of the population. Only the parameter sets with the highest likelihood scores are kept to act as "genetic donors for new parameter sets.

Ultimately, the population of parameter sets converges to the phase-space point with the highest likelihood and, thus, to a set of parameters that most accurately matches the stellar data. We say an optimization has converged when all members of the population are identical. MilkyWay@home can take between two to four weeks to converge.

Although Shelton et al. (2021) give us a realistic physical basis for our choice of time step, there are still regions of our parameter space that could lead to unwieldy run times. In particular, these occur when our generated dwarf progenitor is extraordinarily dense. MilkyWay@home searches over a five-dimensional parameter space for the parameters with the highest likelihood (search range defined in Table 5). However, it would not be wise to explore the full 5D box as the corners of our parameter space with high mass and low radius would produce exceedingly dense progenitors that could take up to 80 days to run. Regions of parameter space that would require more than 1.5 million time steps in the simulation are excluded from our optimization. Figure 10 shows the excised region of parameter space.

Zoom In Zoom Out Reset image size

Table 5. The Search Ranges for Our Optimizations Against Stellar Data

τevolve (Gyr)	RB (kpc)	Radius Ratio $\left({\xi }_{R}=\tfrac{{R}_{B}}{{R}_{B}+{R}_{D}}\right)$	MB (SM)	Mass Ratio $\left({\xi }_{M}=\tfrac{{M}_{B}}{{M}_{B}+{M}_{D}}\right)$
[2.0–6.0]	[0.01–0.50]	[0.1–0.6]	[0.1–100.0]	[0.001–0.95]

Download table as:  ASCIITypeset image

Placing a nonlinear cut on our phase space like this does remove the ability for MilkyWay@home to fit the densest ultracompact dwarf (UCD) galaxies. For example, simulating a dwarf galaxy like M60-UCD1, one of the densest dwarf galaxies ever measured (Strader et al. 2013), would require roughly 17 million time steps to evolve the progenitor through our gravitational well for 6.0 Gyr. Fortunately, our fits to the OCS do not suggest a UCD progenitor, but rather a dwarf galaxy with a similar radius but over 100 times less mass.

We use the Hessian method as described in Section 2.4 to calculate the errors in each of our fitted parameters. However, to accurately calculate these errors, we need our step size to already be comparable to the errors we wish to calculate. To accomplish this, we perform a marginalized one-dimensional parameter sweep along each dwarf parameter to visualize an adequate step size. In this section, we describe the methods by which we quickly and efficiently calculate the marginalized likelihood scores for these parameter sweeps.

Calculating a marginalized likelihood score from N-body simulations is difficult for two reasons. First, the marginalization computation must be completed within the likelihood space, rather than the log(likelihood) space we represent in our version of the likelihood score. This requires us to numerically compute integrals over regions of our parameter space where the likelihood is several orders of magnitude lower than machine precision. Second, any marginalization method we implement over our likelihood space is inherently computationally expensive, requiring us to calculate a complicated four-dimensional integral over a rough likelihood surface. Calculating one likelihood score from a set of parameters requires an N-body simulation. If we were to use a crude three-point Gaussian quadrature method to calculate each integral, that would require us to run 81 N-body simulations to calculate the marginalized likelihood for one point of our parameter sweep. Assuming we use about 30 points for our parameter sweep, that brings the total number of simulations to 2430, requiring one to three months for each parameter sweep.

To address the first issue, we define a special convention for adding together likelihood scores. Given two likelihood scores $\mathrm{ln}({{ \mathcal L }}_{\gt })$ and $\mathrm{ln}({{ \mathcal L }}_{\lt })$ (where ${{ \mathcal L }}_{\gt }\gt {{ \mathcal L }}_{\lt }$), we define their sum $\mathrm{ln}({{ \mathcal L }}_{\mathrm{sum}})$ to follow the following equation:

$$\begin{eqnarray}&&{e}^{\mathrm{ln}({{ \mathcal L }}_{\mathrm{sum}})}={e}^{\mathrm{ln}({{ \mathcal L }}_{\gt })}+{e}^{\mathrm{ln}({{ \mathcal L }}_{\lt })}.\end{eqnarray} \tag{ 28 }$$

Doing some simple math, we can rewrite this equation as

$$\begin{eqnarray}&&\mathrm{ln}({{ \mathcal L }}_{\mathrm{sum}})=\mathrm{ln}({{ \mathcal L }}_{\gt })+\mathrm{ln}\left(1+{e}^{\mathrm{ln}({{ \mathcal L }}_{\lt })-\mathrm{ln}({{ \mathcal L }}_{\gt })}\right).\end{eqnarray} \tag{ 29 }$$

Overall, larger likelihood scores dominate smaller ones. In a double-precision computation, if the two likelihood scores differ by more than 64, the smaller likelihood score is completely ignored in the calculation. Therefore, given a list of likelihoods to add, it is most prudent to add together the smaller likelihood scores first.

To calculate the marginalized likelihood score $\mathrm{ln}(\widetilde{{ \mathcal L }})$, we use a Monte Carlo integration method using N random points in our parameter space of volume V. To determine the minimum number of random points we need for our calculation, we need to know how N affects the error in $\mathrm{ln}(\widetilde{{ \mathcal L }})$. Given a list of likelihoods $\{{ \mathcal L }\}$, the marginalized likelihood $\widetilde{{ \mathcal L }}$ calculated using Monte Carlo is given by the following equation:

$$\begin{eqnarray}&&\widetilde{{ \mathcal L }}=\displaystyle \frac{V}{N}\sum _{i}{{ \mathcal L }}_{i}=V{\mu }_{{ \mathcal L }},\end{eqnarray} \tag{ 30 }$$

where ${\mu }_{{ \mathcal L }}$ is the average of ${ \mathcal L }$. From this, it is clear that the main source of error in this equation comes from the calculation of ${\mu }_{{ \mathcal L }}$, which is simply the error in the mean:

$$\begin{eqnarray}&&{\delta }_{\widetilde{{ \mathcal L }}}=V\displaystyle \frac{{\sigma }_{{ \mathcal L }}}{\sqrt{N}},\end{eqnarray} \tag{ 31 }$$

where ${\sigma }_{{ \mathcal L }}$ is the standard deviation of the population where ${ \mathcal L }$ was pulled from. Performing basic error propagation, we find the error in the marginalized likelihood score to be

$$\begin{eqnarray}&&{\delta }_{\mathrm{ln}(\widetilde{{ \mathcal L }})}=\displaystyle \frac{{\delta }_{\widetilde{{ \mathcal L }}}}{\widetilde{{ \mathcal L }}}=\displaystyle \frac{{\sigma }_{{ \mathcal L }}}{{\mu }_{{ \mathcal L }}\sqrt{N}}.\end{eqnarray} \tag{ 32 }$$

It can then be shown that

$$\begin{eqnarray}&&{\delta }_{\mathrm{ln}(\widetilde{{ \mathcal L }})}=\displaystyle \frac{1}{\sqrt{N-1}}\sqrt{N\displaystyle \frac{\sum {{ \mathcal L }}_{i}^{2}}{{\left(\sum {{ \mathcal L }}_{i}\right)}^{2}}-1}.\end{eqnarray} \tag{ 33 }$$

Equation (33) is what we use in Figure 11 to determine the errors in our parameter sweeps (shaded in blue). However, given this equation, we can generate a rough estimate for what this error looks like on average. We can reorder our list $\{{ \mathcal L }\}$ in descending order, where ${{ \mathcal L }}_{0}$ is the highest likelihood and ${{ \mathcal L }}_{1}$ is the second-highest likelihood. Because we anticipate that each pair of adjacent entries in the reordered list differs by tens of orders of magnitude, we can make the following approximations:

$$\begin{eqnarray}&&\sum _{i}{{ \mathcal L }}_{i}\approx {{ \mathcal L }}_{0}+{{ \mathcal L }}_{1}={{ \mathcal L }}_{0}\left(1+k\right),0\lt k\leqslant 1,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\sum _{i}{{ \mathcal L }}_{i}^{2}\approx {{ \mathcal L }}_{0}^{2}+{{ \mathcal L }}_{1}^{2}={{ \mathcal L }}_{0}^{2}\left(1+{k}^{2}\right),0\lt k\leqslant 1,\end{eqnarray} \tag{ 35 }$$

where k = p₁/p₀. We can make this approximation because we expect only the largest few likelihoods in our list to impact the summations. Substituting these quantities into Equation (33) gives us

$$\begin{eqnarray}\begin{array}{rcl}{\delta }_{\widetilde{{ \mathcal L }}} & \approx & \displaystyle \frac{1}{\sqrt{N-1}}\sqrt{N\displaystyle \frac{\left(1+{k}^{2}\right)}{{\left(1+k\right)}^{2}}-1}\\ & = & \sqrt{\displaystyle \frac{1+{k}^{2}-2k/(N-1)}{{\left(1+k\right)}^{2}}}.\end{array}\end{eqnarray} \tag{ 36 }$$

Taking the limit where k approaches 0, we find that ${\delta }_{\widetilde{{ \mathcal L }}}$ takes the form

$$\begin{eqnarray}&&{\delta }_{\widetilde{{ \mathcal L }}}\approx \sqrt{1-\displaystyle \frac{2k}{N-1}}\approx 1.\end{eqnarray} \tag{ 37 }$$

Surprisingly, we find that the effect of N on ${\delta }_{\widetilde{{ \mathcal L }}}$ is incredibly muted in our likelihood space. In fact, looking at the parameter sweeps in Figure 11, we find that the errors associated with each point are extremely close to 1. It is therefore unreasonable to use a bound in ${\delta }_{\widetilde{{ \mathcal L }}}$ to select a minimum value for N.

Zoom In Zoom Out Reset image size

Instead, we should set N to be the minimum number of measurements by which the unbiased value of σ_p can be most efficiently estimated to a set percentage. To determine σ_p to within 1%, the unbiased variance ${\sigma }_{p}^{2}$ must be determined to within 2%. Also, the relationship between ${\sigma }_{p}^{2}$ and the sample variance ${s}_{p}^{2}$ can be estimated using Bessel's correction:

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }_{p}^{2}}{{s}_{p}^{2}}=\displaystyle \frac{N}{N-1}.\end{eqnarray} \tag{ 38 }$$

Therefore, to know ${\sigma }_{p}^{2}$ to within 2%, we measure the likelihoods from 50 random points for each Monte Carlo approximation to calculate our parameter sweeps.

The likelihood surfaces we observe are not smooth, but we still see an overall dip near the optimal parameters. A Gaussian peak in the likelihood surface translates to a parabola in the log(likelihood) parameter sweeps. Using the apparent width of these peaks, we determine the appropriate step size to use when calculating the Hessian for the purposes of error analysis.

After finding the converged values for each of our optimizations, we reran the simulations using each best-fit parameter set to check that the generated simulated tidal stream matched the data, as shown in Figures 12 and 13. In all optimizations, we see a general trend where the B dispersion is relatively well fit for Λ < 10° but poorly fit for Λ > 10°, especially near the core. This is likely due to the larger errors in the B dispersion associated with the region near the core in the DEC data at higher Λ_OCS. Most of the optimizations (Runs 2, 3, and 5) converged to similar tidal debris that closely resembles the OCS. Runs 1, 4, and 6 had a poorer likelihood score compared to the other optimizations, and we see that the distribution of bodies near its core is wider across Λ_OCS than the data. Due to the generally poor fit of these runs, we exclude them from our analysis of the OCS progenitor.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It should be noted that none of these optimizations reproduce the gap (σ = 1.06) in the tidal stream at Λ = 105. This is not only because our data can barely resolve this gap, but also because such a gap would only appear as a result of substructures within the Milky Way's dark matter halo (Koposov et al. 2019), which is not accounted for in our current halo potential.

Taking the average of the three best-fitted parameters (Runs 2, 3, and 5), we estimate the following physical properties of the OCS progenitor dwarf galaxy:

$$\begin{eqnarray*}\begin{array}{rcl}{\tau }_{\mathrm{evolve}} & = & 3.6337\pm {0.0004}^{(\mathrm{stdev}.)}\pm {0.0004}^{(\mathrm{Hess}.)}\mathrm{Gyr}\\ {R}_{B} & = & 0.1994\pm {0.0167}^{(\mathrm{stdev}.)}\pm {0.0006}^{(\mathrm{Hess}.)}\mathrm{kpc}\\ {R}_{D} & = & 0.773\pm {0.046}^{(\mathrm{stdev}.)}\pm {0.005}^{(\mathrm{Hess}.)}\mathrm{kpc}\\ {M}_{B} & = & \left(2.68\pm {0.07}^{(\mathrm{stdev}.)}\pm {0.03}^{(\mathrm{Hess}.)}\right)\times {10}^{5}{M}_{\odot }\\ {M}_{D} & = & \left(1.94\pm {0.28}^{(\mathrm{stdev}.)}\pm {0.04}^{(\mathrm{Hess}.)}\right)\times {10}^{7}{M}_{\odot }\\ {M}_{\mathrm{total}} & = & \left(1.97\pm {0.28}^{(\mathrm{stdev}.)}\pm {0.04}^{(\mathrm{Hess}.)}\right)\times {10}^{7}{M}_{\odot }\end{array}\end{eqnarray*}$$

where the first error is the standard deviation of the mean and the second error is the propagated error from each individual optimization's Hessian error. We see that our largest source of error comes from the optimization process rather than the curvature of the likelihood surface. So, we only propagate those errors. From these parameters, we calculate a fitted mass-to-light ratio of γ = 73.5 ± 10.6.

Note that these results and errors are derived under the assumption of a perfect Milky Way potential, orbit, and functional form of the dwarf galaxy and do not include systematic errors due to these assumptions.

Our estimated mass for the OCS's progenitor is not consistent with the masses of UFDs measured from the dispersion of the line-of-sight velocities, which assume spherical symmetry and virial equilibrium. To illustrate this point more clearly, let us use these parameters to calculate the total mass within 300 pc of the progenitor's center. Within this radius, we calculate (1.6 ± 0.1) ×10⁵ M_⊙ of baryonic matter and (9.2 ± 2.2) × 10⁵ M_⊙ of dark matter. In total, we find a mass of (1.1 ± 0.2) × 10⁶ M_⊙ within 300 pc of the progenitor's center, roughly one order of magnitude smaller than what is claimed possible in Strigari et al. (2008). However, these results are only a factor of ∼2 smaller than the mass estimated in Newberg et al. (2010), who estimated the total OCS progenitor mass from N-body simulations to be closer to 2.5 × 10⁶ M_⊙. We note that this mass calculation was done under the assumption that the dark matter followed a Plummer profile. However, we note that 300 pc is very close to the half-light radius of the best-fit dwarf progenitor, and Shelton et al. (2021) showed that our algorithm was most sensitive to the mass within the half-light radius.

We use our fitted simulations to estimate the current location of the OCS's progenitor. Note that there is a peak in the star counts of both the data and the stream simulation around Λ ∼ 35° in Figure 12. We know this point of high stellar density must be the location of the progenitor remnant because there is no other reason for stars to pile up at this position, so far from apogalacticon. In this section, we will determine the position of the progenitor remnant and show that it is no longer gravitationally bound.

For each optimization, we translate the positions of all the bodies into spherical coordinates (r, θ, ϕ), where the origin is the Galactic Center and the z-axis is perpendicular to the Galactic Plane. We bin the baryonic (stellar) bodies in the polar (θ) and azimuthal (ϕ) angles and select the bin with the highest density of bodies. Using the bodies in the selected bin, we calculate the average Galactocentric distance (r) and its standard deviation. To improve the calculation, we remove all bodies with distances more than 2.5 standard deviations away from the mean and recalculate the average, a process known as sigma clipping. We perform this sigma clipping 50 times for the sake of being thorough. After calculating the 3D point in space with the highest density of baryonic bodies, we translate it into Galactic and equatorial coordinates. The progenitor position of each run is shown in Table 9.

Our simulations predict that the remnant progenitor core is expected to be at around (l, b) = ((264.9 ± 2.9)°, ((43.6 ±2.8)°) or (α, δ) = ((166.0 ± 0.9)°, (−11.1 ± 2.5)°). This is close to the sky position where Grillmair et al. (2015) detected the OCS's possible missing progenitor, (α, δ) ≈ (167°,−14°).

We next check whether or not the simulated remnant progenitor is still gravitationally bound. For each optimization, we calculate the total kinetic energy of the remnant and compare it to its gravitational potential energy. We select bodies that are within 0.5 kpc of the progenitor's center because that distance is roughly the cylindrical radius of our simulated stream (see Figure 14). We transform our bodies into the center-of-momentum reference frame and add together the kinetic energies $({E}_{K}=\tfrac{1}{2}{\sum }_{i}^{N}{m}_{i}{{v}_{i}}^{2})$ of all the bodies. Similarly, we add together the gravitational energies from each pair of bodies $\left({E}_{P}={\sum }_{i}^{N}{\sum }_{i\lt j}^{N}\tfrac{{{Gm}}_{i}{m}_{j}}{| {r}_{i}-{r}_{j}| }\right)$. We record the energies we calculate for each optimization in Table 10.

Zoom In Zoom Out Reset image size

We find that the simulated remnant progenitor has a kinetic energy of (6.4 ± 0.6) SM kpc² Gyr⁻² and a gravitational potential energy of −(0.094 ± 0.001) SM kpc² Gyr⁻² and is thus gravitationally unbound. Looking at the stream directly in Figure 14 makes this readily apparent.

The tidal debris we observe in our simulations reveals some interesting aspects of the distribution of baryonic and dark matter in the stellar stream. We analyze the distribution of matter within the tails of the simulated OCS. Figure 15 shows that the tails of the OCS are thicker than the central region nearer to the progenitor remnant. There are two possible contributions to this thickening: either the tails are naturally fanning out at apogalacticon due to the symplectic (phase-space area-preserving) property of gravitational systems, or material with a higher velocity dispersion in the progenitor is being tidally stripped first, populating the regions of the stream farther from the core.

Zoom In Zoom Out Reset image size

We can test this by looking at the simulated tidal stream generated using Run 3, our best run, when the core is at apogalacticon. If a high-velocity dispersion has a stronger impact on the thickness of the tails than the fanning at apogalacticon, then the tails should maintain their thickness even when they are not at apogalacticon. However, as can be clearly seen in Figure 16, while the core is at apogalacticon, it bears a similar thickness to the tails in the fully evolved stream. Also, the tails of the stream in the core‒apogalacticon state share a similar thickness to the core in the final state. This indicates that the thickening of the tails in the final state is mostly caused by the natural fanning that occurs at apogalacticon.

Zoom In Zoom Out Reset image size

In addition to being thicker than the core, the tails also possess a higher mass-to-light ratio. We see this in Figure 17, which maps out the density of baryons and dark matter along Λ_OCS. While the core of the OCS has a relatively high dark matter density, it is still not as high as in its tails, where the mass-to-light ratio jumps to as high as 142 in Run 3. Due to its high dark matter concentration and low baryonic contamination, the tails of the OCS might serve as a better candidate for indirect dark matter detection than the central core.

Zoom In Zoom Out Reset image size

When creating our data histogram, we excluded the southern tail of the OCS because our current N-body simulations do not account for the gravitational influence of the LMC. While the OCS's northern tail is not greatly influenced by the presence of the LMC, its southern tail exhibits peculiar perturbations in its stellar velocities, which cannot be explained by a static Milky Way potential or a low-mass LMC (Erkal et al. 2019). While this paper does not fit data from the southern tail, we do plan on including the tail and the gravitational effects of the LMC in a future MilkyWay@home paper.

While we measured the OCS progenitor properties without the assumption of dynamical equilibrium, our measurement still relies on a few key presuppositions. First, we assume the Milky Way has a static, axisymmetric Galactic potential for the entirety of the simulation. There currently exists no mechanism in our N-body simulations for our progenitor to perturb the Milky Way to give rise to second-order effects. We do not include any time dependence in our Milky Way potential. We do not take into account other gravitating systems, which could interact with the OCS, such as the LMC or other Milky Way satellites. Adding the LMC to our simulations would not only further perturb the stream but also induce a reflex motion in the entirety of the Milky Way, forcing the dwarf-galaxy progenitor to fall into an accelerating potential. It has also been recently discovered that most of the mass in the inner disk of the Milky Way is located in its fast rotating bar (Portail et al. 2017), which is also not included in our Galactic model.

Second, we fix the orbit of our simulated OCS throughout each of our optimizations to the values determined in Newberg et al. (2010). Because the orbit is fit to the northern (leading) tail of the OCS, which should be at lower energy than the orbit of the dwarf galaxy itself, we expect the orbit is not exact. We have not explored how the orbit affects the fit properties of the dwarf-galaxy progenitor. However, we are only using the stellar distributions along and across the stream, which are not thought to be strongly affected by small differences in the orbit, to constrain the stream. Additionally, the stream created in our best-fit simulation is similar to the stream observed, which gives us some confidence that our orbit and progenitor properties are reasonable.

Third, we assume that the OCS progenitor has the form of a two-component Plummer sphere. We used a Plummer profile for our progenitor because it is a simple model with an analytic energy distribution function. However, it is probably not the best (and certainly not a perfect) model, particularly for the dark matter halo. We are unsure whether using a different profile, such as a Hernquist or Navarro–Frenk–White profile, would significantly affect the optimization results.

Finally, there is the measurement of baryonic mass in the OCS. As discussed in Section 2.4.1, using the redder NGC 5053 to estimate the mass of the bluer OCS overestimates the mass per turnoff star. From our analysis, we find that the stellar mass could be reduced by a factor of 4 due to inaccuracies in our calculation of the stellar mass per turnoff. This would affect the mass-to-light ratio but would not significantly change the total mass within 300 pc of the center of the progenitor.