Resolved Stellar Streams around NGC 4631 from a Subaru/Hyper Suprime-Cam Survey

In this study, we use the HSC⁵ imager on the 8.2 m Subaru Telescope. HSC consists of 104 2048 × 4096 science CCDs with a scale of $0\buildrel{\prime\prime}\over{.} 17$ per pixel and covers a total field of view (FoV) of 15 in diameter (Miyazaki et al. 2012). In addition, the HSC FoV corresponds to about 190 kpc at the distance of NGC 4631. The HSC observations of the current target were taken during two nights in 2015 March (S15A-046; PI: Tanaka). Our HSC field includes both NGC 4631 and NGC 4656, as shown in Figure 1. The HSC observations were made with HSC-g and HSC-i filters, with total exposure times of 6.0 hr and 11.8 hr, respectively. The weather conditions were excellent, with stable seeings of $0\buildrel{\prime\prime}\over{.} 48$ in HSC-g and $0\buildrel{\prime\prime}\over{.} 60$ in HSC-i, which are measured on the final stacked images. Table 1 summarizes the full details of the observations.

Zoom In Zoom Out Reset image size

Table 1.  Observation Status for a Target Field

Coordinates	Filter	${t}_{\exp }$	${N}_{\exp }$	Airmass	Seeing
(${\alpha }_{2000}/{\delta }_{2000}$)		(s)			(arcsec)
${12}^{{\rm{h}}}{42}^{{\rm{m}}}{59}^{{\rm{s}}}0$	HSC-g	21660	92	1.02–1.60	048
+32°22'170	HSC-i	42601	180	1.02–1.92	060

Download table as:  ASCIITypeset image

To reduce the raw HSC data, we basically adopt the HSC pipeline^{6 ,7} 3.8.5, which is based on an earlier version of the LSST pipeline (Axelrod et al. 2010), calibrated against SDSS Data Release 7 astrometry. Since the sky subtraction process in this version of the pipeline did not work well for the case in which CCD chips suffer from apparently large objects such as nearby galaxies, we develop and apply our own sky subtraction method in such a case, as summarized in the Appendix. Since the data were taken from different airmasses, as shown in Table 1, flux scaling of each CCD chip is carried out before combining the science images (Bosch et al. 2017, in preparation). First, zero-point magnitudes of each CCD chip were estimated based on the stars with SDSS photometry imaged in each chip, Then, the zero-point magnitude of each exposure and the spatial variation over the focal plane, which is expressed as the 5th-order Chebyshev polynomial, were fitted by comparing the fluxes of the same star imaged in different exposures. By adopting this flux scaling, the spatial variation of flat-fielding and the extinction by airmass are corrected. As a result, the flux is calibrated within a photometric accuracy of 0.034 mag in the g-band, and 0.028 mag in the i-band. Finally, we obtained a $3\sigma$-clipped mean-stacked image for each band.

We conducted point-spread function (PSF)-fitting photometry using PyRAF in a similar manner as Tanaka et al. (2010). The PSF models were statistically chosen based on a reduced chi-squared computed by the PyRAF/DAOPHOT PSF task. When conducting detection and photometry, we divided a HSC field into several subfields and photometry was carried out for each subfield separately, which is shown by the red dashed grids of Figure 1 because spatial variations of PSF were caused by a huge FoV of HSC. We detect 2,301,816 (1,594,340) objects in the g-band (i-band). Finally, we merged two independent g-band and i-band catalogs into a combined catalog (including 1,218,065 objects) using a 3-pixel matching radius.

We calculated the color terms and photometric zero-points per second of both bands by comparing stellar objects from SDSS Data Release 12. The standard objects are selected under the restrictions of SDSS parameters ($\mathrm{psfMagErr}\lt 0.1$), with those from our photometric catalog selected under the restrictions of PyRAF parameters ($\mathrm{merr}\lt 0.01$, $| \mathrm{sharpness}| \lt 0.1$, and $\mathrm{chi}\lt 2$). Figure 2 shows the best-fitted relations between color terms and photometric zero-points. For further inspection, we checked to what extent photometric zero-points vary with position in our HSC field. The differences between the best-fitted zero-points and the ones in the subfields are 0.03 mag in the g-band and 0.02 mag in the i-band, which are $3\sigma$-clipped standard deviations. This result is consistent with the photometric accuracy of the flux calibrations conducted in the reduction process.

Zoom In Zoom Out Reset image size

Reddening correction is applied to each object in this HSC field based on the extinction maps of Schlegel et al. (1998) and its spatial variation is negligible. In addition, we adopt standard conversion laws, A_g = 3.793E(B – V) and A_i = 2.086E(B – V), in the SDSS's filter system (Schlafly & Finkbeiner 2011).

In order to evaluate incompleteness due to a low signal-to-noise ratio and crowding, we performed the artificial star experiments basically in the similar manner as Tanaka et al. (2010). We have added a sufficient number of artificial stars to each original image in each subfield and each magnitude using the PyRAF addstar task, and the ranges of the magnitudes are $24.0\leqslant g\leqslant 29.5$ and $22.0\leqslant g\leqslant 28.0$, with a binning step of about 0.25–1.0 mag. To prevent these stars from interfering with one another, we divided the HSC frames into grids of cells of 100 pixels of width and randomly added one artificial star to each cell for each run. In addition, we constrained each artificial star to have 10 pixels from the edges of the cell (see also Tanaka et al. 2010). About 5400 stars were configured for each run, and the total number is about 6,274,800 in each image. After detection and photometry were performed in the same manner as in the previous section, we regarded an output object as the same one as an input object based on the standard deviation of the magnitude difference between the input and output objects within the matching radius. Then, the 50% and 80% completeness limits were determined by fitting the analytical model described by Fleming et al. (1995). An example for gird ID = (5, 4) is shown in Figure 3.

Zoom In Zoom Out Reset image size

Figure 4 shows the spatial distributions of magnitudes with 80% and 50% completeness limits in each band. Near the edge of the FoV, the completeness limits are somewhat shallow due to a low signal-to-noise ratio, and close to the two main galaxies in this galaxy group, the limits are somewhat shallow due to crowding. On the other hand, the completeness limits of other fields such as NGC 4631's halo are deep enough to reach the tip of the red giant branch (TRGB) of the galaxy system.

Zoom In Zoom Out Reset image size

A wide-field, deep, and high-resolution HSC image enables us to confirm already-known satellite galaxies and search for still-unknown ones in the galaxy group consisting of NGC 4631 and NGC 4656. Figure 5 shows three color composites of dwarf candidates, and Table 2 lists detailed information on their coordinates, stellar populations, and corresponding references. Through visual inspection, we detect eight new candidate satellite galaxies around NGC 4631 (HSC-5 to 12), and we confirm three previously reported ones (HSC-1, 2, and 4). In the composite image at the location of DGSAT-3, a low surface brightness galaxy near NGC 4631 reported by Javanmardi et al. (2016), we do not find a galaxy, and suggest that this was a misclassified object due to blending of foreground stars and background galaxies (see our object HSC-3).

Zoom In Zoom Out Reset image size

Table 2.  Candidates for Satellite Galaxies

ID	R.A.(J2000)	Decl.(J2000)	Stellar Population	References
HSC-1	${12}^{{\rm{h}}}{42}^{{\rm{m}}}53\buildrel{\rm{s}}\over{.} 1$	+32°27'190	old? (unresolved)	2, 3, 4
HSC-2	${12}^{{\rm{h}}}{42}^{{\rm{m}}}06\buildrel{\rm{s}}\over{.} 1$	+32°37'148	RGB	2, 3, 4
HSC-3	${12}^{{\rm{h}}}{41}^{{\rm{m}}}08\buildrel{\rm{s}}\over{.} 0$	+32°26'504	not satellite?	4
HSC-4	${12}^{{\rm{h}}}{41}^{{\rm{m}}}50\buildrel{\rm{s}}\over{.} 3$	+32°31'019	MS	1
HSC-5	${12}^{{\rm{h}}}{43}^{{\rm{m}}}44\buildrel{\rm{s}}\over{.} 8$	+32°32'029	RGB	5
HSC-6	${12}^{{\rm{h}}}{43}^{{\rm{m}}}24\buildrel{\rm{s}}\over{.} 8$	+32°28'553	AGB/RGB	5
HSC-7	${12}^{{\rm{h}}}{40}^{{\rm{m}}}10\buildrel{\rm{s}}\over{.} 0$	+32°39'314	MS/RSG/AGB/RGB	5
HSC-8	${12}^{{\rm{h}}}{41}^{{\rm{m}}}47\buildrel{\rm{s}}\over{.} 2$	+32°51'242	MS/RSG/AGB/RGB	5
HSC-9	${12}^{{\rm{h}}}{40}^{{\rm{m}}}53\buildrel{\rm{s}}\over{.} 0$	+32°16'558	AGB/RGB	5
HSC-10	${12}^{{\rm{h}}}{42}^{{\rm{m}}}31\buildrel{\rm{s}}\over{.} 4$	+31°58'093	AGB/RGB	5
HSC-11	${12}^{{\rm{h}}}{42}^{{\rm{m}}}01\buildrel{\rm{s}}\over{.} 7$	+32°24'066	young? (unresolved)	5
HSC-12	${12}^{{\rm{h}}}{43}^{{\rm{m}}}07\buildrel{\rm{s}}\over{.} 1$	+32°29'273	MS/RSG/AGB/RGB	5

Note. The number of the ID column corresponds to the IDs in Figures 5 and 8. The coordinates are the central positions of each panel in Figure 5. The stellar population column shows the RGB/AGB/RSG/MS maps in which each candidate is detected. HSC-1 and HSC-11 are unresolved in our HSC image, while HSC-3, which was previously reported as a low surface brightness galaxy by a small amateur telescope (Javanmardi et al. 2016), is apparently a blending object in our HSC image.

References. (1) Seth et al. (2005a), (2) Karachentsev et al. (2014), (3) Martínez-Delgado et al. (2015), (4) Javanmardi et al. (2016), (5) this work.

Download table as:  ASCIITypeset image

HSC-1 and HSC-11 show somewhat clumpy features located in each of these galaxy systems, but are unresolved in our HSC image. The former has a relatively red diffuse structure in Figure 5, suggesting that it predominantly consists of an old stellar population, whereas a relatively blue image in the latter implies a young stellar object such as a dwarf irregular galaxy. The other candidates are clearly resolved in this HSC image; for instance, a diffuse low surface brightness galaxy, HSC-2, which was also reported by previous unresolved studies based on small amateur telescopes and HSC-1 (Karachentsev et al. 2014; Martínez-Delgado et al. 2015; Javanmardi et al. 2016), has an extended structure consisting of red, old stars, as can be seen in the composite image. On the other hand, the HSC-4 reported by the HST observation of Seth et al. (2005a) is located in the thick disk of NGC 4631 (see also Figure 8), and the existence of the resolved blue stars supports their conclusion that it is a young dwarf galaxy with an age of ∼30 Myr, estimated from the comparison of their resolved CMD with theoretical isochrones.

Both HSC-7 and HSC-8 look like dwarf irregular galaxies with blue, young stellar populations, while HSC-12, which is located near a saturated bright foreground star, looks like a compact dwarf irregular galaxy with blue stars. On the other hand, HSC-5, 6, 9, and 10 are classified as dwarf spheroidal galaxies dominated by red, old stellar populations in appearance. Further details of these candidate dwarf galaxies will be discussed in the forthcoming paper.

The left panel of Figure 6 shows the log-scaled CMD over the entire HSC field derived by plotting the number density of stars within $0.05\times 0.05$ mag boxes. It is compared with PARSEC isochrones (Bressan et al. 2012) for old populations with fixed ages (10 Gyr) and different metallicties ([M/H] = −2.28, −1.68, −1.28, −0.68, and −0.38) (orange solid lines) and young populations with different ages (4.0, 10.0, 17.8, 31.6, 56.2, 100, 177.8 Myr) and fixed metallicity ([M/H] = −0.38) (orange dashed lines), assuming that the distance modulus of NGC 4631 is 7.4 Mpc (Radburn-Smith et al. 2011). In order to make the CMD, we select stellar-like objects under the restrictions of PyRAF parameters ($\mathrm{merr}\lt 0.5$, $| \mathrm{sharpness}| \lt 0.5$ and $\mathrm{chi}\lt 5$). Apparently, the CMD indicates the existence of the broad red giant branch (RGB) at ${(g-i)}_{0}\sim 1$ and ${i}_{0}\sim 26$ and a young main-sequence (MS) at ${(g-i)}_{0}\sim -1$.

Zoom In Zoom Out Reset image size

To examine an effect of Galactic foreground contaminations, we produce a synthetic CMD in the same Galactic coordinate and the FoV as this study, based on the Besançon models (Robin et al. 2003), convolving photometric errors and the observational incompleteness (see the right panel of Figure 6). The synthetic CMD is consistent with the observational one, and both of them clearly indicate foreground dwarf stars in the MW halo at ${(g-i)}_{0}\sim 0.3$, as well as those in the MW disk at ${(g-i)}_{0}\sim 2.5$ in a bright magnitude range of ${i}_{0}\lesssim 24$. Because NGC 4631 is located at a high Galactic latitude ($b\sim 84^\circ$), foreground contaminations of the MW disk are relatively small. In Figure 6, the colored boxes show the same criteria of stellar populations as the observed CMD, indicating that a region of MS is little contaminated by foreground stars. In fact, the ratio of the predicted number of foreground stars to the total number of detected objects is about 0.5%. Although regions of asymptotic giant branch (AGB) and RGB seem largely contaminated, these ratios are actually 1.4% and 1.0%, respectively. This ratio for a region of red super giants (RSGs) is also small, 2.3%. Therefore, dominant contaminations in a faint magnitude range of ${i}_{0}\gtrsim 24.5$ of the observed CMD are high-redshift background galaxies, which are unresolved point sources in ground-based telescopes. In spite of the fact that analysis based on selection boxes of stellar populations on CMDs is expected to largely reduce effects from background galaxies, it is not easy to completely distinguish local signals of substructures from heavy contaminations due to local density peaks originating from cosmic variance. We thus confine ourselves to consider the local overdensities, which are evident from the visual inspection of the color composite image (Figure 6) in this study or those that were already suggested in some previous studies. On the other hand, regarding statistical estimation of the physical properties of spatially spread structures in Sections 3.4 and 3.5, we note that the usage of control fields shown in the next section is effective, provided that background galaxies are uniformly distributed in these fields.

We divide our photometric objects into four groups, which are characterized by MS, AGB, RSG, and RGB, based on the color cuts shown in the CMD. However, in selections of MS and RGB we limit the groups to objects brighter than 50% of completeness limits. In Figure 7, we display the spatial distribution for the number density of photometric objects per 0.01 × 0.01 square degrees corresponding to 74 × 74 kpc², for each of the CMD areas enclosed with different color lines in Figure 6. These maps indicate many interesting substructures characterized by previously known/unknown dwarf companion galaxies and stellar streams around NGC 4631. In particular, some of the dwarf candidates discussed in Section 3.1 are identified as overdensity regions in these maps. For example, the spatial density map of MS stars shows a remarkable overdensity region, HSC-8, which is apparently a star-forming dwarf galaxy, at the north part of NGC 4631, as well as the main bodies of NGC 4631 and NGC 4656. The spatial density map of RSG stars clearly shows a core of NGC 4627, a dwarf Elliptical (dE) satellite of NGC 4631, in addition to the objects presented in the MS map. In the AGB map,  several other dwarf satellites with intermediate age populations such as HSC-6 exist between NGC 4631 and NGC 4656.

Zoom In Zoom Out Reset image size

The RGB map clearly indicates the two resolved stellar streams at the northwest side (Stream NW) and the southeast side (Stream SE) of NGC 4631 that were previously reported by Martínez-Delgado et al. (2015), although there are apparently artificial patchy structures due to faint magnitudes near the completeness limits and crowding. Both streams mainly consist of likely old (possibly 10 Gyr or even older) stellar populations, because we cannot detect any counterparts to the two streams in the other maps, although AGB stars are slightly concentrated in the core of Stream NW. Furthermore, the lack of young stellar populations in Streams SE and NW is consistent with the lack of spatial correlation between these streams and the ${{\rm{H}}}_{{\rm{I}}}$ distribution discussed in Martínez-Delgado et al. (2015). In the subsequent sections, we investigate the photometric properties of these two stellar streams based on the resolved stellar population.

The star count maps whose noises are reduced also exhibit some as-yet-unknown candidates of dwarf galaxies around NGC 4631, as well as the stellar streams and vertically extended disk of the galaxy. Figure 8 shows the zoomed-in spatial density distribution of RGB stars (displayed in the bottom right panel of Figure 7) to highlight substructures around NGC 4631. The density is estimated by counting the number of photometric objects within a bin of 0.005 × 0.005 square degrees corresponding to 37 × 37 kpc², and overdensities with a signal over an S/N of 3 are marked in the map, assuming that the global background noise are distributed in Poisson statistics. The two red, dashed ellipses in the map indicate Streams SE and NW, while the two red dashed rectangles illustrate control fields for statistically estimating background and foreground contaminations in these stream fields (see also the next section). Therefore, these control fields are located at almost the same projected vertical distance from NGC 4631's disk as the stream fields, and do not have any remarkable substructures of high S/N within these fields. The four orange circles present the candidates of dwarf galaxies reported by the previous studies (Seth et al. 2005a; Karachentsev et al. 2014; Martínez-Delgado et al. 2015; Javanmardi et al. 2016). The eight blue ones show candidate new dwarf galaxies identified through the visual inspection in this study (see Section 3.1). We note that one of these candidates is located in the control field for Stream SE, but our analysis based on the resolved stellar populations is not affected because this is unresolved. Additionally, the other dwarf candidates without high S/Ns in the figure indicate unresolved objects such as HSC-1 and HSC-11. The reason that the signal of HSC-12 is apparently deflected to the right side from the center is due to a saturated bright star overlapping substantially with the dwarf galaxy. In addition, since HSC-4 located at the old thick disk is a young dwarf galaxy detected by the HST observation (Seth et al. 2005a), it mainly consists of MS stars. The resolved stellar populations based on which of the new dwarf candidates are visually identified are summarized in Table 2.

Zoom In Zoom Out Reset image size

The TRGB, which is driven by core helium ignition, is a useful indicator to estimate the distance to resolved old stellar systems such as nearby galaxies as well as globular clusters (e.g., Salaris & Cassisi 2005; Tanaka et al. 2010, 2011).

The i-band magnitude of the TRGB in metal-poor populations of $[\mathrm{Fe}/{\rm{H}}]\lesssim -0.7$ dex changes by less than 0.1 mag (Lee et al. 1993). Therefore, the TRGB was historically detected as a sharp cutoff of the i-band luminosity function (LF) with the application of  edge-detection algorithms such as the Sobel filter (e.g., Madore & Freedman 1995; Sakai et al. 1996). However, the edge-detection algorithm was confronted with the high levels of Poisson noise that abound in the poorly populated structures of galaxy halos. On the other hand, the maximum-likelihood TRGB detection method, in which a predefined model LF is fitted to the observed distribution of the stars, is robust against the strong Poisson noise (e.g., Mèndez et al. 2002; Makarov et al. 2006). In a more recent variation on the maximum-likelihood method, a Bayesian inferential approach for the TRGB detection was developed due to dramatic computational improvement (e.g., Conn et al. 2011, 2012; Tollerud et al. 2016). The traditional maximum-likelihood approach is based on the uncertain assumption that an estimate is normally distributed, whereby a Bayesian approach can estimate a full picture of a complex probability distribution even in a poorly populated structure.

Therefore, we adopt the Bayesian method as described in Conn et al. (2011), to derive the distance to the two streams around NGC 4631. The plots for Streams SE and NW are summarized in Figures 9 and 10, respectively. The top left and top right panels of these figures show log-scaled CMDs of the streams and their control fields, respectively. Using all detected objects, within the red rectangles of these CMDs, we construct completeness-corrected LFs smoothed with a Gaussian kernel, assuming that each object normally distributes at the measured i-band magnitude, ${ \mathcal N }(m,\sigma )$, where m and σ are the i-band magnitude and its photometric error calculated from artificial star experiments in Section 2.3, respectively (Sakai et al. 1996). These LFs are shown in the bottom right panels of Figures 9 and 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The background contaminations are not negligible in the magnitude range of faint RGBs in this study. Therefore, we reduce the effect based on the matched filter method adopted in Conn et al. (2012). To do so, we introduce a kernel density estimation (KDE), ${ \mathcal W }(\alpha ,\delta )$, which is the probability distribution function of individual object positions transformed into a smooth surface density function with a Gaussian kernel and a bandwidth h that is tuned to each distribution with the prescription $h=1.06\sigma {N}_{* }^{-0.2}$ (Silverman 1986, where σ is the standard deviation of the samples and N_* is the total number of samples.), instead of the radial density profiles of dwarf galaxies of Conn et al. (2012). Then, the LF, ${\rm{\Phi }}(m)$, within a particular magnitude bin, m, is given as

$$\begin{eqnarray}&&{\rm{\Phi }}(m)=\displaystyle \sum _{k=1}^{\mathrm{ndata}}{ \mathcal N }({m}_{k},{\sigma }_{k})(m){{ \mathcal C }}_{k}{ \mathcal W }({\alpha }_{k},{\delta }_{k}),\end{eqnarray} \tag{ 1 }$$

where ${\sigma }_{k}$, ${{ \mathcal C }}_{k}$, and ${ \mathcal W }({\alpha }_{k},{\delta }_{k})$ are the photometric error, the inverse of completeness and the KDE weight of kth object, respectively, and ndata is the total number of objects within the red parallelograms of the CMDs. Using resampling data of extinction-corrected i-band magnitude, $\{{m}_{1},\,\cdots ,\,{m}_{N}\}$, derived from Equation (1), we calculate the following likelihood function,

$$\begin{eqnarray}&&{ \mathcal L }(\{{m}_{1},\cdots ,{m}_{N}\}| {m}_{\mathrm{TRGB}},a)=\displaystyle \prod _{n=1}^{N}\phi ({m}_{n}| {m}_{\mathrm{TRGB}},a),\end{eqnarray} \tag{ 2 }$$

where the model LF, $\phi ({m}_{n}| {m}_{\mathrm{TRGB}},a)$, is described by

$$\begin{eqnarray}\left\{\begin{array}{ll}{\phi }_{\mathrm{RGB}}({m}_{n})+{\phi }_{\mathrm{BG}}({m}_{n}), & ({m}_{n}\geqslant {m}_{\mathrm{TRGB}})\\ {\phi }_{\mathrm{BG}}({m}_{n}), & ({m}_{n}\lt {m}_{\mathrm{TRGB}})\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where ${\phi }_{\mathrm{RGB}}({m}_{n})$ is ${10}^{a(m-{m}_{\mathrm{TRGB}})}$ (Makarov et al. 2006) and ${\phi }_{\mathrm{BG}}({m}_{n})$ is a 6th-order polynomial. Then, we manually determine the fraction of background objects, $f={D}_{\mathrm{BG}}/{D}_{\mathrm{SIGNAL}}$, by calculating the average density of objects in the control field ${D}_{\mathrm{BG}}$ and in the stream field ${D}_{\mathrm{SIGNAL}}$ in order to match the observational LF from Equation (1) with the model LF from Equation (3), namely,

$$\begin{eqnarray}&&{\int }_{{m}_{\mathrm{TRGB}}}^{{m}_{2}}{\phi }_{\mathrm{RGB}}({m}_{n}){dm}=1-f,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\int }_{{m}_{1}}^{{m}_{2}}{\phi }_{\mathrm{BG}}({m}_{n}){dm}=f,\end{eqnarray} \tag{ 5 }$$

based on Conn et al. (2011). Eventually, a posterior distribution function is described by

$$\begin{eqnarray}&&p({m}_{\mathrm{TRGB}},a| \{{m}_{1},\cdots ,\,{m}_{N}\})\\ &&\quad \propto { \mathcal L }(\{{m}_{1},\cdots ,\,{m}_{N}\}| {m}_{\mathrm{TRGB}},a)p({m}_{\mathrm{TRGB}})p(a),\end{eqnarray} \tag{ 6 }$$

where $p({m}_{\mathrm{TRGB}})$ and p(a) are the simplest non-informative prior distributions, that is, uniform, due to the principle of insufficient reason. The parameters, ${m}_{\mathrm{TRGB}}$ and a, are currently chosen for the model by a Markov Chain Monte Carlo (MCMC) algorithm. For sampling of the parameters and marginal likelihood estimation, we apply the standard Metropolis-Hasting algorithm with the symmetric proposal distribution, that is, Gaussian. Our algorithm based on 32 MCMC chains and 101,000 iterations, for which the first 1000 steps are discarded as burn-in, each chain adequately converges beyond burn-in, and in total has over 3 million effective random numbers in each stream.

Figure 11 shows posterior distributions and contour maps for the extinction-corrected TRGB magnitude and the LF slope a of Streams SE and NW based on the effective random numbers. The finally estimated distances to the two streams are summarized in Table 3, assuming that the absolute TRGB magnitude is ${M}_{i,\mathrm{SDSS}}^{\mathrm{TRGB}}=-3.44\pm 0.10$ (Bellazzini 2008), suggesting that Stream NW is relatively more distant from us than Stream SE. The best-fit LFs based on expected a posteriori (EAP) estimation are shown with red solid lines in the bottom right panels of Figures 9 and 10.

Zoom In Zoom Out Reset image size

Table 3.  Distance to the Streams

$(\mathrm{Mpc})$	MAPa	EAPb	5%	25%	50%	75%	95%
Stream SE	6.90	7.10	6.22	6.84	6.97	7.09	7.29
Stream NW	7.48	7.91	6.44	7.41	7.59	7.87	7.97

Notes. The 90% (50%) credible intervals for the distance to Streams SE and NW are [6.22, 7.29] ([6.84, 7.09]) and [6.44, 7.97] ([7.41, 7.87]) Mpc, respectively.

^aMaximum a posteriori. ^bExpected a posteriori.

Download table as:  ASCIITypeset image

In order to construct the metallicity distributions (MDFs) for the two stellar streams, we first assume that all of stars in the CMDs are old (t = 10 Gyr), so that the location of the RGB stars in the CMD depends almost solely on metallicity. The photometric metallicity of each star is derived from comparison with theoretical PARSEC isochrones (Bressan et al. 2012) on the CMD. Figure 12 shows the interpolated metallicity map on the CMD. The interpolation is conducted based on the radial basis function (Rbf) of the Python/SciPy package. Then, we select the target RGB stars brighter than ${i}_{0}=27.0$ mag, considering the uncertainty of incompleteness. Also, since metallicity conjecture outside the criteria of the orange solid lines in Figure 12 is quite uncertain, we select only the data inside the orange ones.

Zoom In Zoom Out Reset image size

The stream fields are highly contaminated by the background galaxies and stars of NGC 4631's halo. Therefore, we subtract the MDF of the control field, which is a nearly flat distribution, from the MDF of the stream field. Figure 13 shows the resultant contamination-subtracted, normalized MDFs for the two streams. The vertical error bars denote a nominal uncertainty in each metallicity bin as derived from the Poisson errors. Both MDFs have a broad distribution ranging from the most metal-poor part of the isochrones to the most metal-rich part and there is a clear high-metallicity peak at $[M/H]\gt -1$. The comprehensive shapes of both MDFs are quite similar as reproduced by a Gaussian mixture model (GMM), which is a weighted sum of two-component Gaussian densities.

Zoom In Zoom Out Reset image size

According to the N-body simulation of Martínez-Delgado et al. (2015), both streams originated from a tidal interaction between NGC 4631 and a single dwarf satellite. If it is true, Streams SE and NW have the same stellar population. Therefore, it is worth comparing the MDs in order to examine the validity of their model. Our approach is based on Bayesian model selection, and to do so, we prepare two models, H₀ and H₁, as follows.

In model H₀, we assume that both MDFs are reproduced by the same GMM, such as the following equation:

$$\begin{eqnarray}{x}_{i}^{\mathrm{SE}} & \sim & \displaystyle \sum _{k=1}^{2}{w}_{k}{ \mathcal N }({\mu }_{k},{\sigma }_{k}),\\ {x}_{j}^{\mathrm{NW}} & \sim & \displaystyle \sum _{k=1}^{2}{w}_{k}{ \mathcal N }({\mu }_{k},{\sigma }_{k}),\end{eqnarray} \tag{ 7 }$$

where ${x}_{i}^{\mathrm{SE}}$ and ${x}_{j}^{\mathrm{NW}}$ are the data of Streams SE and NW, respectively, which are re-sampled based on a probability distribution function, that is, a contamination-subtracted, normalized MDF, and the mixture weights satisfy the constraint that ${\sum }_{k=1}^{2}{w}_{k}=1$. The likelihood then takes the form

$$\begin{eqnarray}{ \mathcal L }({\boldsymbol{x}}| {{\boldsymbol{\theta }}}_{{H}_{0}}) & = & { \mathcal L }({{\boldsymbol{x}}}_{\mathrm{SE}},{{\boldsymbol{x}}}_{\mathrm{NW}}| {{\boldsymbol{\theta }}}_{{H}_{0}})\\ & = & \displaystyle \prod _{i=1}^{{N}_{\mathrm{SE}}}f({x}_{i}^{\mathrm{SE}}| {{\boldsymbol{\theta }}}_{{H}_{0}})\displaystyle \prod _{j=1}^{{N}_{\mathrm{NW}}}f({x}_{j}^{\mathrm{NW}}| {{\boldsymbol{\theta }}}_{{H}_{0}}),\end{eqnarray} \tag{ 8 }$$

where ${{\boldsymbol{\theta }}}_{{H}_{0}}=(w,{\mu }_{1},{\mu }_{2},{\sigma }_{1},{\sigma }_{2})$. By Bayes' theorem, we have the posterior probability distribution for the parameters, $f({{\boldsymbol{\theta }}}_{{H}_{0}}| {\boldsymbol{x}})\propto { \mathcal L }({\boldsymbol{x}}| {{\boldsymbol{\theta }}}_{{H}_{0}})f({{\boldsymbol{\theta }}}_{{H}_{0}})$, where each prior, $f({{\boldsymbol{\theta }}}_{{H}_{0}})$, obeys visually restricted normal distributions.

On the other hand, in model H₁ we assume that both MDFs are respectively reproduced by different GMM such as follows:

$$\begin{eqnarray}{x}_{i}^{\mathrm{SE}} & \sim & \displaystyle \sum _{k=1}^{2}{w}_{k}^{\mathrm{SE}}{ \mathcal N }({\mu }_{k}^{\mathrm{SE}},{\sigma }_{k}^{\mathrm{SE}}),\\ {x}_{j}^{\mathrm{NW}} & \sim & \displaystyle \sum _{k=1}^{2}{w}_{k}^{\mathrm{NW}}{ \mathcal N }({\mu }_{k}^{\mathrm{NW}},{\sigma }_{k}^{\mathrm{NW}}).\end{eqnarray} \tag{ 9 }$$

The likelihood function is

$$\begin{eqnarray}&&{ \mathcal L }({\boldsymbol{x}}| {{\boldsymbol{\theta }}}_{{H}_{1}})=\displaystyle \prod _{i=1}^{{N}_{\mathrm{SE}}}f({x}_{i}^{\mathrm{SE}}| {{\boldsymbol{\theta }}}_{{H}_{1}}^{\mathrm{SE}})\displaystyle \prod _{j=1}^{{N}_{\mathrm{NW}}}f({x}_{j}^{\mathrm{NW}}| {{\boldsymbol{\theta }}}_{{H}_{1}}^{\mathrm{NW}}),\end{eqnarray} \tag{ 10 }$$

where ${{\boldsymbol{\theta }}}_{{H}_{1}}=({{\boldsymbol{\theta }}}_{{H}_{1}}^{\mathrm{SE}},{{\boldsymbol{\theta }}}_{{H}_{1}}^{\mathrm{NW}})$,

$$\begin{eqnarray*}{{\boldsymbol{\theta }}}_{{H}_{1}}^{\mathrm{SE}} & = & ({w}^{\mathrm{SE}},{\mu }_{1}^{\mathrm{SE}},{\mu }_{2}^{\mathrm{SE}},{\sigma }_{1}^{\mathrm{SE}},{\sigma }_{2}^{\mathrm{SE}}),\\ {{\boldsymbol{\theta }}}_{{H}_{1}}^{\mathrm{NW}} & = & ({w}^{\mathrm{NW}},{\mu }_{1}^{\mathrm{NW}},{\mu }_{2}^{\mathrm{NW}},{\sigma }_{1}^{\mathrm{NW}},{\sigma }_{2}^{\mathrm{NW}}).\end{eqnarray*}$$

From Bayes' theorem, we have the posterior probability distribution for the parameters, $f({{\boldsymbol{\theta }}}_{{H}_{1}}| {\boldsymbol{x}})\propto { \mathcal L }({\boldsymbol{x}}| {{\boldsymbol{\theta }}}_{{H}_{1}})f({{\boldsymbol{\theta }}}_{{H}_{1}})$, where each prior, $f({{\boldsymbol{\theta }}}_{{H}_{1}})$, also obeys visually restricted normal distributions.

In order to estimate the posterior distributions, we apply an approximate Bayesian inference: automatic differentiation variational inference (ADVI) of the PyStan package, which is known as a faster algorithm than MCMC sampling (Kucukelbir et al. 2015). However, the GMM inference strongly depends on initial values, because the model has many local optimum and almost non-informative priors. Therefore, we explore a maximum of posterior distribution by updating initial values with higher likelihood through 1000 iterations. Then, 20,000 samples per iteration are made by the ADVI algorithm. Best-fit GMMs of both streams based on a maximum of posterior are represented by the solid lines in Figure 13. In addition, we summarize posterior predictive distributions for the best-fit GMM parameters for each model in Table 4. According to the estimations from model H₁, Stream SE seems to be systematically more metal-poor than model H₀, while Stream NW is systematically more metal-rich.

Table 4.  Posterior Predictive Distributions of the Best-fit GMM Parameters

Model	${\boldsymbol{\theta }}$	MAP	EAP	5%	25%	50%	75%	95%
	w	0.175	0.189	0.166	0.179	0.189	0.199	0.214
	${\mu }_{1}$	−1.247	−1.249	−1.257	−1.252	−1.249	−1.246	−1.241
H0	${\mu }_{2}$	−0.861	−0.855	−0.869	−0.861	−0.855	−0.850	−0.841
	${\sigma }_{1}$	0.323	0.333	0.308	0.322	0.332	0.343	0.358
	${\sigma }_{2}$	0.205	0.205	0.197	0.202	0.205	0.208	0.213
	${w}^{\mathrm{SE}}$	0.177	0.181	0.148	0.166	0.180	0.194	0.217
	${\mu }_{1}^{\mathrm{SE}}$	−1.309	−1.310	−1.322	−1.315	−1.310	−1.305	−1.298
	${\mu }_{2}^{\mathrm{SE}}$	−0.915	−0.905	−0.923	−0.912	−0.905	−0.897	−0.886
	${\sigma }_{1}^{\mathrm{SE}}$	0.341	0.326	0.286	0.308	0.325	0.342	0.368
H1	${\sigma }_{2}^{\mathrm{SE}}$	0.190	0.189	0.179	0.185	0.189	0.192	0.198
	${w}^{\mathrm{NW}}$	0.248	0.261	0.223	0.244	0.260	0.276	0.300
	${\mu }_{1}^{\mathrm{NW}}$	−1.159	−1.155	−1.170	−1.161	−1.155	−1.148	−1.139
	${\mu }_{2}^{\mathrm{NW}}$	−0.789	−0.785	−0.809	−0.795	−0.785	−0.776	−0.761
	${\sigma }_{1}^{\mathrm{NW}}$	0.317	0.321	0.293	0.309	0.321	0.333	0.351
	${\sigma }_{2}^{\mathrm{NW}}$	0.200	0.203	0.193	0.199	0.203	0.208	0.214

Download table as:  ASCIITypeset image

Our goal in this section is to test our hypotheses, H₀ and H₁. Here, we introduce the Watanabe–Akaike (or Widely Applicable) information criterion (WAIC; Watanabe 2010) for the model selection. It can be viewed as an improvement on the deviance information criteria for Bayesian models (Spiegelhalter et al. 2002), and one of a family of criteria that estimate the predictive power of a model, that is, how a model can anticipate new data. In addition, the better known Akaike information criterion (AIC; Akaike 1974) uses the maximum-likelihood estimate, while the WAIC averages over the posterior distribution of the parameters. Namely, the WAIC is more effective than the traditional AIC in the Bayesian approach. WAIC_H for a model, H, is defined as

$$\begin{eqnarray}&&{\mathrm{WAIC}}_{H}=-2\mathrm{Mean}(\mathrm{ln}({{ \mathcal L }}^{* }({\boldsymbol{x}}| {{\boldsymbol{\theta }}}_{H})))+2\mathrm{Var}(\mathrm{ln}({{ \mathcal L }}^{* }({\boldsymbol{x}}| {{\boldsymbol{\theta }}}_{H}))),\end{eqnarray} \tag{ 11 }$$

where ${{ \mathcal L }}^{* }({\boldsymbol{x}}| {{\boldsymbol{\theta }}}_{H})$ is a posterior predictive distribution for a model, H, that is, H₀ or H₁. As a result, the WAICs of our models, H₀ and H₁, are 7271.9 and 9000.0, respectively. Therefore, we can automatically adopt the model H₀, and we conclude that both MDFs are reproduced by the same GMM in this study, suggesting that Streams SE and NW have the same stellar population originating from a single dwarf satellite. Based on the model of the same stellar population, we can immediately estimate the metallicity probability distribution of the progenitor from a posterior predictive distribution as summarized in Table 5.