Resolved stellar streams around NGC4631 from a Subaru/Hyper Suprime-Cam survey

TL;DR

This study uses Subaru/HSC imaging to identify and analyze stellar streams around NGC4631, revealing their properties, distances, and likely common origin from a dwarf galaxy interaction.

Contribution

First detailed analysis of stellar streams around NGC4631 using resolved stellar populations and Bayesian methods to determine their properties and origin.

Findings

Streams are located in front of and behind NGC4631.

Streams likely originated from a single dwarf galaxy interaction.

Progenitor's stellar mass estimated at ~3.7 x 10^8 Msun.

Abstract

We present the first results of the Subaru/Hyper Suprime-Cam (HSC) survey of the interacting galaxy system, NGC4631 and NGC4656. From the maps of resolved stellar populations, we identify 11 dwarf galaxies (including already-known dwarfs) in the outer region of NGC4631 and the two tidal stellar streams around NGC4631, named Stream SE and Stream NW, respectively. This paper describes the fundamental properties of these tidal streams. Based on the tip of red giant branch method and the Bayesian statistics, we find that StreamSE (7.10 Mpc in Expected a posteriori, EAP, with the 90% credible intervals of [6.22, 7.29] Mpc) and StreamNW (7.91 Mpc in EAP with the 90% credible intervals of [6.44, 7.97] Mpc) are located in front of and behind NGC4631, respectively. We also calculate the metallicity distribution of stellar streams by comparing the member stars with theoretical isochrones on the color-magnitude diagram. We find that both streams have the same stellar population based on the Bayesian model selection method, suggesting that they originated from a tidal interaction between NGC4631 and a single dwarf satellite. The expected progenitor has a positively skewed metallicity distribution function with [M/H]_EAP=-0.92 with the 90% credible intervals of [-1.46, -0.51]. The stellar mass of the progenitor is estimated as 3.7 x 10e+8 Msun with the 90% credible intervals of [5.8 x 10e+6, 8.6 x 10e+9] Msun based on the mass-metallicity relation for Local group dwarf galaxies. This is in good agreement with an initial stellar mass of the progenitor presumed in the previous N-body simulation.

Tables

Table 1. Table 1 Observation Status for a Target Field

Coordinates	Filter	$t_{\exp }$	$N_{\exp }$	Airmass	Seeing
(${\alpha }_{2000}/{\delta }_{2000}$)		(s)			(arcsec)
${12}^h{42}^{m}59\stackrel{s}{.}0$	HSC-$g$	21660	92	$1.02-1.60$	$0\stackrel{\prime \prime }{.}48$
$+{32}^{\circ }22\prime 17\stackrel{\prime \prime }{.}0$	HSC-$i$	42601	180	$1.02-1.92$	$0\stackrel{\prime \prime }{.}60$

Table 2. Table 2 Candidates for Satellite Galaxies

ID	RA(J2000)	Dec(J2000)	Stellar Population	Ref.††footnotemark:
HSC-1	${12}^h{42}^{m}53\stackrel{s}{.}1$	$+{32}^{\circ }27\prime 19\stackrel{\prime \prime }{.}0$	old? (unresolved)	2,3,4
HSC-2	${12}^h{42}^{m}06\stackrel{s}{.}1$	$+{32}^{\circ }37\prime 14\stackrel{\prime \prime }{.}8$	RGB	2,3,4
HSC-3	${12}^h{41}^{m}08\stackrel{s}{.}0$	$+{32}^{\circ }26\prime 50\stackrel{\prime \prime }{.}4$	not satellite?	4
HSC-4	${12}^h{41}^{m}50\stackrel{s}{.}3$	$+{32}^{\circ }31\prime 01\stackrel{\prime \prime }{.}9$	MS	1
HSC-5	${12}^h{43}^{m}44\stackrel{s}{.}8$	$+{32}^{\circ }32\prime 02\stackrel{\prime \prime }{.}9$	RGB	5
HSC-6	${12}^h{43}^{m}24\stackrel{s}{.}8$	$+{32}^{\circ }28\prime 55\stackrel{\prime \prime }{.}3$	AGB/RGB	5
HSC-7	${12}^h{40}^{m}10\stackrel{s}{.}0$	$+{32}^{\circ }39\prime 31\stackrel{\prime \prime }{.}4$	MS/RSG/AGB/RGB	5
HSC-8	${12}^h{41}^{m}47\stackrel{s}{.}2$	$+{32}^{\circ }51\prime 24\stackrel{\prime \prime }{.}2$	MS/RSG/AGB/RGB	5
HSC-9	${12}^h{40}^{m}53\stackrel{s}{.}0$	$+{32}^{\circ }16\prime 55\stackrel{\prime \prime }{.}8$	AGB/RGB	5
HSC-10	${12}^h{42}^{m}31\stackrel{s}{.}4$	$+{31}^{\circ }58\prime 09\stackrel{\prime \prime }{.}3$	AGB/RGB	5
HSC-11	${12}^h{42}^{m}01\stackrel{s}{.}7$	$+{32}^{\circ }24\prime 06\stackrel{\prime \prime }{.}6$	young? (unresolved)	5
HSC-12	${12}^h{43}^{m}07\stackrel{s}{.}1$	$+{32}^{\circ }29\prime 27\stackrel{\prime \prime }{.}3$	MS/RSG/AGB/RGB	5

Table 3. Table 3 Distance to the Streams

$\left(Mpc\right)$	MAPaaMaximum a posteriori	EAPbbExpected a posteriori	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

Table 4. Table 4 Posterior Predictive Distributions of the Best Fit GMM Parameters

Model	$𝜽$	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$
$H_0$	${\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$
$H_1$	${\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$

Table 5. Table 5 Metallicity Probability Distribution of the Progenitor

	MAP	EAP	5%	25%	50%	75%	95%
[M/H]	$-0.85$	$-0.92$	$-1.46$	$-1.07$	$-0.89$	$-0.73$	$-0.51$

Equations

$$\Phi (m)=\sum _{k=1}^{ndata}\mathcal{N}(m_k,{\sigma }_k)(m){\mathcal{C}}_k\mathcal{W}({\alpha }_k,{\delta }_k),$$

$$\mathcal{L}(\{m_1,\cdots ,m_N\}|m_{\mathrm{TRGB}},a)=\prod _{n=1}^N\varphi (m_n|m_{\mathrm{TRGB}},a),$$

$$\displaystyle\left\{\begin{array}[]{ll}\phi_{\rm RGB}(m_{n})+\phi_{\rm BG}(m_{n}),\ (m_{n}\geq m_{\rm TRGB})\\ \phi_{\rm BG}(m_{n}),\ (m_{n}<m_{\rm TRGB})\end{array}\right.$$

$${\int }_{m_{\mathrm{TRGB}}}^{m_2}{\varphi }_{\mathrm{RGB}}(m_n)dm$$

$${\int }_{m_1}^{m_2}{\varphi }_{\mathrm{BG}}(m_n)dm$$

$$p$$

$$x_i^{\mathrm{SE}}$$

$$x_j^{\mathrm{NW}}$$

$$\mathcal{L}(\mathbf{x}|{\mathbf{\theta }}_{H_0})$$

$$x_i^{\mathrm{SE}}$$

$$x_j^{\mathrm{NW}}$$

$$\mathcal{L}(\mathbf{x}|{\mathbf{\theta }}_{H_1})=\prod _{i=1}^{N_{\mathrm{SE}}}f(x_i^{\mathrm{SE}}|{\mathbf{\theta }}_{H_1}^{\mathrm{SE}})\prod _{j=1}^{N_{\mathrm{NW}}}f(x_j^{\mathrm{NW}}|{\mathbf{\theta }}_{H_1}^{\mathrm{NW}}),$$

$${\mathbf{\theta }}_{H_1}^{\mathrm{SE}}$$

$${\mathbf{\theta }}_{H_1}^{\mathrm{NW}}$$

$$WAI{C}_H=-2Mean(\ln ({\mathcal{L}}^{\ast }(\mathbf{x}|{\mathbf{\theta }}_H)))+2Var(\ln ({\mathcal{L}}^{\ast }(\mathbf{x}|{\mathbf{\theta }}_H)))$$

$$g-V=0.371(g-i)+\mathrm{0.068.}$$

Full text

1. Introduction

The currently favored cosmology based on $\Lambda$-dominated cold dark matter ($\Lambda$CDM) predicts the hierarchical formation of large galaxies like the Milky Way and M31. Indeed, several $N$-body simulations based on this cosmology (e.g., Bullock & Johnston, 2005; Cooper et al., 2010) succeeded to reproduce diverse tidal stellar streams with various surface brightnesses, metallicities, velocities and morphologies depending on mass accretion history of galaxies (e.g., Johnston et al., 2008). Therefore, stellar streams as observed in nearby galaxies provide an important opportunity to assess the outputs of the $\Lambda$CDM model.

The Milky Way and M31 play a special role in these studies because they offer us very detailed information on stellar streams through their resolved stars. For instance, prominent stellar streams, the Sagittarius stream in the Milky Way (e.g., Ibata et al., 1994; Majewski et al., 2003) and Giant Southern Stream in M31’s halo (e.g., Ibata et al., 2001, 2007; Tanaka et al., 2010; Ibata et al., 2014) are the most spectacular examples, thereby enabling us to derive their dynamical histories in light of the $\Lambda$CDM model. Recently, large ground-based telescopes such as the 8.2m Subaru telescope made it possible to resolve stellar halos of other nearby galaxies beyond the Local Group, which reveal the diversity of stellar structures of galaxy outskirts with various morphologies (e.g., Mouhcine et al., 2010; Tanaka et al., 2011; Greggio et al., 2014; Okamoto et al., 2015; Crnojević et al., 2016). Furthermore, the GHOSTS survey with the Hubble Space Telescope (HST), which has observed various field along the minor and major axes of 16 nearby disk galaxies (Monachesi et al., 2013), also reveals the diversity in halo color profiles and stellar halo masses, supporting the galaxy-to-galaxy scatter in halo stellar properties, a consequence of the stochasticity inherent in the assembling history of galaxies (Monachesi et al., 2016; Harmsen et al., 2016). However, the pencil beam surveys with the HST such as GHOSTS have a drawback that they are less able to identify new streams and overdensities because of their narrow fields.

On the other hand, Martínez-Delgado et al. (2010) exploited small reflecting telescopes of no more than about 50 cm in order to study stellar halos and tidal features in more distant galaxies through their integrated light. This is because integrated surface brightness is independent of distance. Subsequently, other groups inspired by Martínez-Delgado et al. (2010) have adopted this observational methodology to perform a direct test of the $\Lambda$CDM cosmology (e.g., van Dokkum et al., 2014; Javanmardi et al., 2016; Merritt et al., 2016). These authors also found the varieties in the stellar halos in more distant nearby spiral galaxies beyond the Local Volume, implying stochasticity in the accretion histories of galaxies. These integrated light surveys are complementary to the GHOSTS and surveys with large ground-based telescopes based on resolved stellar populations of galaxies.

In this work, we select NGC 4631 as another test-bed aiming for a better understanding of the varieties in stellar halos. NGC 4631 is an edge-on spiral galaxy (Sc) at a distance of about $7$ Mpc (e.g., Tikhonov et al., 2006; Radburn-Smith et al., 2011). It appears that its companion galaxy, NGC 4656, is interacting with NGC 4631. In fact, the outskirt of NGC 4631 is interesting as an interacting system, because the presence of its extraplanar components is reported by a variety of observations, from X-rays to radio wavebands. In particular, Rand (1994) shows five HI spurs surrounding NGC 4631. Schechtman-Rook & Hess (2012) discovered a UV-bright, tidal dwarf galaxy candidate in the NGC 4631/4656 galaxy group, which may be originated from extremely disturbed HI, suggesting that these objects are young with age of 200-300 Myr. However, the relation between these peculiar gas structures and resolved stellar populations in this galaxy system remains unclear. Thus, the outskirts of this galaxy probed by its resolved stars will provide a unique laboratory to study a nearby system which is violently perturbed by neighbors. In particular, we expect that such an interacting galaxy is surrounded by some stellar streams or substructures with/without HI gas. In fact, Seth et al. (2005b) found somewhat thickened disk structures by HST/ACS observations, and Martínez-Delgado et al. (2015) discovered stellar streams at the north-west side (Stream NW) and the south-east side (Stream SE) of NGC 4631 based on the integrated surface light analysis. However, they found that these streams do not overlap with any HI components.

In this paper, we report on the resolved stellar populations of the two tidally-spreading streams around NGC 4631 using the Subaru telescope. The layout of this paper is as follows. In Section 2, we show our observations and detailed procedures of reduction and photometry for our data. In Section 3, we present the discovery of new dwarf galaxies around NGC 4631, the color-magnitude diagrams (CMDs) and spatial distribution of various stellar populations in our observing field. By focusing on the two stellar streams, we show more about the derivation of their distance and metallicity distributions, and discuss the origin of the streams later in the section. In Section 4, we draw discussion and conclusions.

2. Observation and Data Reduction

2.1. HSC Observations

In this study, we use the HSC111http://www.subarutelescope.org/Observing/Instruments/HSC/index.html imager on the 8.2-m Subaru Telescope. HSC consists of 104 $2048\times 4096$ science CCDs with a scale of $0\farcs 17$ per pixel and covers a total field-of-view (FoV) of 1.5 degree 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 March 2015 (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 the total exposure time of 6.0h and 11.8h, respectively. The weather conditions were excellent with stable seeing of $0\farcs 48$ in HSC-$g$ and $0\farcs 60$ in HSC-$i$ which are measured on the final stacked images. Table 1 summarizes full details of the observations.

To reduce the raw HSC data, we basically adopt the HSC pipeline 3.8.5222http://hsca.ipmu.jp333http://hsc.mtk.nao.ac.jp/pipedoc_e/index.html, 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 that 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 in different airmass 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, are fitted by comparing the fluxes of 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 the photometric accuracy of 0.034 mag in $g$-band, 0.028 mag in $i$-band. Finally, we obtained a $3\sigma$-clipped mean-stacked image for each band.

2.2. Photometry and Calibration

We conducted point-spread function (PSF)-fitting photometry using PyRAF in the similar manner to Tanaka et al. (2010). The PSF models were statistically chosen based on reduced chi-squared computed by the PyRAF/DAOPHOT PSF task. When conducting detection and photometry, we divide a HSC field into several subfields and photometry is carried out for each subfield separately 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 $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 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 ($psfMagErr<0.1$) with those from our photometric catalog selected under the restrictions of PyRAF parameters ($merr<0.01$, $|sharpness|<0.1$ and $chi<2$). Figure 2 shows best-fitted relations between color terms and photometric zero-points. For further inspection, we check to what extent photometric zero-points vary with position in our HSC field. The differences between the best-fitted zero-points and ones in the subfields are 0.03 mag in $g$-band and 0.02 mag in $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.

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.793$E(B-V) and $A_{i}=2.086$E(B-V), in the SDSS’s filter system (Schlafly & Finkbeiner, 2011).

2.3. Completeness

In order to evaluate incompleteness due to low signal-to-noise ratio and crowding, we have performed the artificial star experiments basically in the similar manner to Tanaka et al. (2010). We have added a sufficient number of artificial stars to each original image in each subfield and each magnitude by using the PyRAF $addstar$ task, and the range of the magnitudes is $24.0\leq g\leq 29.5$ and $22.0\leq g\leq 28.0$ with a binning step of about 0.25-1.0 mag. To prevent these stars from interfering with one another, we have divided the HSC frames into grids of cells of 100 pixels width and have randomly added one artificial star to each cell for each run. In addition, we constrain each artificial star to have 10 pixels from the edges of the cell (see also Tanaka et al. (2010)). About 5,400 stars were configured for each run, and the total number is about 6,274,800 in each image. After detection/photometry were performed in the same manner as in the previous section, we regard an output object as the same one as an input object based on standard deviation of the magnitude difference between the input and output objects within the matching radius. Then, the 50 and 80 percent 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.

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

3. Results

3.1. Satellite Galaxies around NGC 4631

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 of their coordinates, stellar populations and corresponding references. Through the visual inspection, we detect 8 new candidate satellite galaxies around NGC 4631 (HSC-5 to 12), and we confirm 3 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).

HSC-1 and HSC-11 show somewhat clumpy features located at each of these galaxy systems, but unresolved in our HSC image. The former has a relatively red diffuse structure in Figure 5, suggesting that it dominantly consists of an old stellar population, whereas a relatively blue image of 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 is 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, 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 of an age of $\sim 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 also 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.

3.2. Color-Magnitude Diagrams

The left panel of Figure 6 shows the log-scaled CMD over the entire HSC field derived by plotting number density of stars within $0.05\times 0.05$ mag boxes. It is compared with PARSEC isochrones (Bressan et al., 2012) for old population 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 population with different age (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 ($merr<0.5$, $|sharpness|<0.5$ and $chi<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$.

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-hand 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}\hskip 3.00003pt\raisebox{1.72218pt}{<}\raisebox{-3.01385pt}{\sim}\hskip 3.00003pt24. Because NGC 4631 is located at the 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 (RSG) is also small, 2.3%. Therefore, dominant contaminations in a faint magnitude range of i_{0}\hskip 3.00003pt\raisebox{1.72218pt}{>}\raisebox{-3.01385pt}{\sim}\hskip 3.00003pt24.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 originated 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 already suggested in some previous studies. On the other hand, regarding statistical estimation of physical properties of spatially-spread structures in Section 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.

3.3. Spatial Density Distributions

We divide our photometric objects into the four groups, which are characterized by MS, AGB, RSG and RGB, based on the color cuts shown in the CMD. However, we limit to brighter objects than 50% completeness limits in selections of MS and RGB. In Figure 7, we display the spatial distribution for the number density of photometric objects per 0.01 $\times$ 0.01 square degrees corresponding to 74 $\times$ 74 kpc2, 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 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 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, other several dwarf satellites with intermediate age populations such as HSC-6 exist between NGC 4631 and NGC 4656.

The RGB map clearly indicates the two resolved stellar streams at the north-west side (Stream NW) and the south-east side (Stream SE) of NGC 4631 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 Stream SE and NW is consistent with spatial consistent with lack of spatial correlation between these streams and the HI 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 the vertically extended disk of the galaxy. Figure 8 shows the zoomed 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 $\times$ 0.005 square degrees corresponding to 37 $\times$ 37 kpc2, 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 Stream 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 does 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. Besides, the other dwarf candidates without high $S/N$ 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 the new dwarf candidates are visually identified are summarized in Table 2.

3.4. TRGB distance to the streams

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 [Fe/H]\hskip 3.00003pt\raisebox{1.72218pt}{<}\raisebox{-3.01385pt}{\sim}\hskip 3.00003pt-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 an edge-detection algorithm such as 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, 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 has been 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 Stream SE and NW are summarized in Figure 9 and 10, respectively. The top-left and top-right panels of these figures show log-scaled CMDs of the stream and its control field, 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 $\sigma$ 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 Figure 9 and 10.

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}_{\ast}$ (Silverman, 1986, where $\sigma$ is the standard deviation of the samples and $N_{\ast}$ is the total number of samples.), instead of radial density profiles of dwarf galaxies of Conn et al. (2012). Then, the LF, $\Phi(m)$, within a particular magnitude bin, $m$, is given as

[TABLE]

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 $k$th object, respectively, and $ndata$ is 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,

[TABLE]

where the model LF, $\phi(m_{n}|m_{\rm TRGB},a)$, is described by

[TABLE]

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

[TABLE]

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

[TABLE]

where $p(m_{\rm TRGB})$ and $p(a)$ are the simplest non-informative prior distributions, that is, $uniform$, due to the principle of insufficient reason. The parameters, $m_{\rm 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, which the first 1,000 steps are discorded as burn-in, each chain adequately converge beyond burn-in, and has totally 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 Stream SE and NW based on the effective random numbers. The finally-estimated distance to the two streams are summarized in Table 3, assuming that the absolute TRGB magnitude is $M^{\rm TRGB}_{i,{\rm SDSS}}=-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 Figure 9 and 10.

3.5. Metallicity Distributions of the Streams

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 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.

The stream fields are highly contaminated by 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 contaminations-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]$>-1$. The comprehensive shapes of both MDFs are quite similar as reproduced by Gaussian mixture model (GMM) which is a weighted sum of 2 component Gaussian densities.

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, Stream 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_{0}$ and $H_{1}$ as follows.

In model $H_{0}$, we assume that both MDFs are reproduced by the same GMM such as the following equation,

[TABLE]

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

[TABLE]

where $\bm{\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(\bm{\theta}_{H_{0}}|\bm{x})\propto\mathcal{L}(\bm{x}|\bm{\theta}_{H_{0}})f(\bm{\theta}_{H_{0}})$, where each prior, $f(\bm{\theta}_{H_{0}})$, obeys visually-restricted normal distributions.

On the other hand, in model $H_{1}$, we assume that both MDFs are respectively reproduced by different GMM such as the following equation,

[TABLE]

The likelihood function is

[TABLE]

where $\bm{\theta}_{H_{1}}=(\bm{\theta}^{\rm SE}_{H_{1}},\bm{\theta}^{\rm NW}_{H_{1}})$,

[TABLE]

By Bayes’ theorem, we have the posterior probability distribution for the parameters, $f(\bm{\theta}_{H_{1}}|\bm{x})\propto\mathcal{L}(\bm{x}|\bm{\theta}_{H_{1}})f(\bm{\theta}_{H_{1}})$, where each prior, $f(\bm{\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 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 1,000 iterations. Then, 20,000 samples per an iteration are made by the ADVI algorithm. Eventually, best fit GMMs of both streams based on a maximum of posterior are described by solid lines of 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_{1}$, Stream SE seems to be systematically more metal-poor than model $H_{0}$, while Stream NW is systematically more metal-rich.

Our goal in this section is to test our hypotheses, $H_{0}$ and $H_{1}$. 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 (DIC) 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. WAICH for a model, $H$, is defined as

[TABLE]

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

4. Discussion and Concluding Remarks

In this paper, we report the first results of the Subaru/HSC survey of the interacting galaxy system consisting of NGC 4631 and NGC 4656. We have conducted the careful analysis for our HSC data based on the HSC pipeline and our own sky-subtraction procedure. Eventually, we have detected 11 dwarf galaxies (including already-known dwarfs) in the outer region of NGC 4631 and the two tidal stellar streams around NGC 4631 in the stellar density maps divided into each resolved stellar population of MS, RSG, AGB and RGB. To discuss individual properties of these dwarf galaxies is beyond the scope of this paper, hence it will be described in the forthcoming paper.

In this paper, we particularly focus on the fundamental properties of the two tidal streams. Based on the TRGB method and the Bayesian statistics, we find that Stream SE (at the heliocentric distance of 7.10 Mpc in EAP with the 90% credible intervals of [6.22, 7.29] Mpc) is located in front of NGC 4631 along the line-of-sight, while Stream NW (at 7.91 Mpc in EAP with the 90% credible intervals of [6.44, 7.97] Mpc) is nestled behind NGC 4631. On the other hand, we calculate metallicity distribution of each stream by comparing each star with theoretical isochrones on the color-magnitude diagram (CMD), and we find the possibility that both streams have the same stellar population based on the Bayesian model selection method, suggesting that they originated from a tidal interaction between NGC 4631 and a single dwarf satellite.

The expected progenitor has a positively skewed MDF with $\rm[M/H]_{\rm EAP}=-0.92$ with the 90% credible intervals of [$-1.46,-0.51$]. Provided that $\rm[M/H]=[Fe/H]+log(0.694\times[\alpha/Fe]+0.306)$ (Cassisi & Salaris, 2013) and $\rm[\alpha/Fe]=0.3$ ($\alpha$-enhancement of old stellar populations of the Milky Way), we estimate the mass of the progenitor as $3.7\times 10^{8}M_{\sun}$ with the 90% credible intervals of [$5.8\times 10^{6},8.6\times 10^{9}$] $M_{\sun}$ based on the mass-metallicity relation for Local group dwarf galaxies (Kirby et al., 2013). This is in good agreement with a total initial stellar mass ($5.4\times 10^{8}M_{\sun}$) of the progenitor presumed in the $N$-body simulation of Martínez-Delgado et al. (2015). However, red/metal-rich ends of these MDFs suffer from incompleteness even if we carefully construct them by considering completeness function derived in Section 2.3. Therefore we cannot exclude that the true MDFs of the streams may be more metal-rich if such red stars are present in these fields. This also suggests that the mass of the progenitor galaxy estimated in this study may be a lower limit.

To compare with other resolved stellar streams in galaxies of the local universe, we calculate the surface brightness of Stream SE and NW. First, we extract the secure RGB stars for both the two streams and each control field within orange lines shown in Figure 12. In this procedure, we adopt the lower limit of $i$-band magnitude to 1.5 mag fainter than each TRGB magnitude estimated in Section 3.4. Second, we convert the summed-up flux counts of selected stars to the surface brightness in mag arcsec*-2*, and subtract the surface brightness of the control field, for which remaining foreground and background contaminations are removed based on the statistical method as discussed in previous sections. To convert the surface brightness in HSC filter systems to the one in standard Johnson $V$-band, we adopt the following formula (Komiyama et al. 2017, in preparation):

[TABLE]

Finally, we measure the mean surface brightnesses within the red-dashed circles of Figure 8 of $\langle\mu_{V}\rangle=31.0\pm 0.02$ mag arcsec*-2* for Stream SE and $\langle\mu_{V}\rangle=32.2\pm 0.02$ mag arcsec*-2* for Stream NW. The errors are estimated from the Poisson statistics from the finite number of observed stars and subtracted contaminations. Although Stream SE and NW are the brightest stellar streams in NGC 4631’s system, they may be relatively faint in comparison with other resolved stellar streams detected in galaxies of the local universe according to observational evidence of the relation between metallicity and surface brightness of stellar substructures of M31 and NGC 55 (Gilbert et al., 2009; Tanaka et al., 2011).

We examine the spatial structures of surface brightness of the two streams and find that the brightest regions of Stream SE and NW are located at $(\alpha,\delta)\sim(190.75\arcdeg,32.38\arcdeg)$ and $(\alpha,\delta)\sim(190.24\arcdeg,32.79\arcdeg)$, and their surface brightnesses within 0.01 deg2 bin corresponding to 1.4 kpc2 and 1.2 kpc2 are $\mu_{\rm V}=30.1$ mag arcsec*-2* and $\mu_{\rm V}=29.5$ mag arcsec*-2*, respectively. Although the errors estimated from Poisson noise are less than 0.01 mag arcsec*-2*, the values of these surface brightness are probably lower limits due to the uncertainties of incompleteness corrections of detection and blending. Their positions are consistent with density peaks of the KDE maps of Figure 9 and 10. Therefore, Stream NW has a brighter and more concentrated core than Stream SE, implying that the main body of the progenitor is probably associated with Stream NW rather than Stream SE. In that case, it is reasonably thought that Stream SE formed through a tidal interaction between a dwarf satellite embedded in Stream NW and NGC 4631 several Gyrs ago. This interpretation is consistent with the prediction by the $N$-body simulation of (Martínez-Delgado et al., 2015). On the other hand, a continuous, bridge-like structure between Stream SE and NW that the model predicts are not detected in this study. If such a faint structure exists, its surface brightness is less than background, that is, $\mu_{\rm V}=32.8\pm 0.3$ mag arcsec*-2* in the control field for Stream SE and $\mu_{\rm V}=33.1\pm 0.3$ mag arcsec*-2* in the control field for Stream NW. These errors indicate standard deviations of the surface brightness within the total bins in each control field. Furthermore, the projected length of that structure is $\sim 100$ kpc, as measured between the cores of the two stream components, which is consistent with that simulation. Meanwhile, that structure might have the large extent more than several hundred kpc along the line-of-sight.

Currently, there is observational evidence that stellar halos may become less common at lower stellar masses than Milky Way mass spiral galaxies (e.g., Tanaka et al., 2011; Streich et al., 2016) although there is a variation in the masses of stellar halos of spiral galaxies with stellar masses similar to that of the Milky Way (Merritt et al., 2016; Harmsen et al., 2016). NGC 4631 is interpreted as a large Magellanic-type galaxy for many years, and it is classified as a relatively late-type spiral galaxy in the Local Volume (de Vaucouleurs & de Vaucouleurs, 1963). Therefore, we can infer that NGC 4631 is less massive galaxy than the Milky Way. Notwithstanding, NGC 4631 has a relatively active accretion history that continues to influence growth of its stellar halo characterized by the two large tidal stellar streams. It is expected that the further analysis of its stellar halo constrains formation mechanism of low mass stellar halos.

We thank the staff at the Subaru telescope for their excellent support during our observing runs. We wish to recognize and acknowledge the very significant cultural role and reverence that the sacred summit of Mauna a Wakea, which is the dwelling place of the goddess Poli’ahu, has always had within the endemic Hawaiian community. We are the most fortunate to have the opportunity to conduct observations from this sacred mountain. We are grateful to the referee for helpful comments and suggestions. This work was supported by JSPS KAKENHI Grant Number JP25800098 for MT and JP15K05037 for YK. This work was supported in part by MEXT Grant-in-Aid for Scientific Research on Innovative Areas (No.JP15H05889, JP16H01086 for MC).

Appendix A technique of sky subtraction

The raw data were basically reduced with HSC pipeline 3.8.5. We used SDSS Data Release 7 as an astrometric catalog in the pipeline. However, the sky-subtraction process of the pipeline didn’t work well in regard to CCD chips suffered from an apparently large object such as nearby galaxies (see also the left panel of Figure 15). Therefore, we applied our own sky-subtraction method to such contaminated CCD chips. The top panel of Figure 14 shows a color map of unmasked flat-fielded data consisting 10 ccd chips surrounding NGC 4631’s main body. First, we masked large objects, that is, NGC 4631 and NGC 4656, in flat-fielded data considering their coordinates, major-axis and minor-axis sizes based on NED as well as saturated/bad pixels based on mask information made from the HSC pipeline and ones with higher signal-to-noise ratio than $2\sigma$. Second, we conducted binning with 128 pixels for the data due to the limitation of our computational power. Third, assuming that the masked flat-fielded data is corresponding to natural sky gradients, we reconstructed sky-gradient model (see the bottom panel of Figure 14) using the fourth order Chebyshev polynomial.

The right hand panel of Figure 15 shows the sky-subtracted image based on our sky-gradient model, suggesting that our sky-subtraction process works better than that of the HSC pipeline. In order to statistically compare the sky-subtracted data by the HSC pipeline with the one by our procedure, we investigate how pixel counts (median values stacked along the $x$-axis) along the $y$-axis of the 2k $\times$ 4k CCD chips of Figure 15 change. Figure 16 shows that the data reduced by our procedure is reaching to no background gradient around zero with increasing distance from NGC 4631’s main body, while the data reduced by the HSC pipeline is strongly swelling with highly-deviated background. On the other hand, sky of unsuffered CCD chips was subtracted using the sky models estimated in each CCD chip based on our sky-subtraction procedure.