EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 671 (March 2023) A&A, 671 (2023) A141 Full HTML
Open Access
Issue
A&A
Volume 671, March 2023
Article Number A141
Number of page(s) 24
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202245011
Published online 16 March 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Within the hierarchical framework for galaxy formation, the stellar halos of massive galaxies are expected to form predominantly through a succession of mergers with lower-mass satellite galaxies. Cosmological numerical simulations of galaxy assembly in the Lambda-cold dark matter (ΛCDM) paradigm predict that satellite disruption occurs throughout the lifetime of all massive galaxies (e.g., Bullock & Johnston 2005; De Lucia & Helmi 2008; Cooper et al. 2010, 2013; Pillepich et al. 2015; Rodriguez-Gomez et al. 2016). As a consequence, the stellar halos of massive galaxies at the current epoch should contain a wide variety of diffuse remnants of tidally disrupted dwarf satellites. The most spectacular examples of such debris are long, dynamically cold stellar streams1 that wrap around the host galaxy and roughly trace the orbit of the progenitor satellite (e.g., Sanders & Binney 2013). Although these fossils of galactic cannibalism disperse into amorphous clouds of stars (through phase-mixing) in a few gigayears, ΛCDM simulations predict that significant numbers of coherent tidal features may be detectable with sufficiently deep observations in the outskirts of the majority of nearby massive galaxies (Font et al. 2011; Cooper et al. 2013; Pillepich et al. 2015).

The detection and characterisation of the faint remnants of dwarf galaxy accretion is a vital test of the hierarchical nature of galaxy formation (Johnston et al. 2008) which has not yet been fully exploited, mainly because their extremely low surface brightness makes them challenging to observe. Although some of the most luminous examples of tidal tails and shells around massive elliptical galaxies have been known for many decades (e.g., Arp 1966; Malin & Carter 1980), recent studies have shown that fainter analogues of these structures are common around spiral galaxies in the local Universe, including the Milky Way and Andromeda. The first known galactic tidal stream around the Milky Way (the Sagittarius stream) was discovered almost three decades ago (Ibata et al. 1994; Mateo et al. 1998; Martínez-Delgado et al. 2001; Majewski et al. 2003). In recent years, a new generation of wide-field digital imaging surveys (including the Sloan Digital Sky Survey, Pan-STARRS, and the Dark Energy Survey) have revealed a multitude of fainter streams (Belokurov et al. 2006; Shipp et al. 2018) and other amorphous sub-structures in our Galaxy’s stellar halo (e.g., Newberg et al. 2002; Jurić et al. 2008; Bernard et al. 2016). The Pan-Andromeda Archaeological Survey (PAndAS; Ibata et al. 2007; McConnachie et al. 2009) has provided an equally revolutionary panoramic view of similar features in the stellar halo of M31. These observations provide sound empirical support for the ΛCDM prediction that dwarf galaxy disruption is responsible for the streams and diffuse stellar halos around massive late-type galaxies.

A search for comparable fossil records in a much larger sample of nearby galaxies is required to allow for the possibility that the recent merging histories of the two Local Group spirals may not be ‘typical’ (e.g., Mutch et al. 2011). Improvements in ΛCDM numerical simulations have helped to guide the quest to discover and characterise stellar streams (e.g., Bullock & Johnston 2005; Cooper et al. 2010). Recent simulations have demonstrated that the characteristics of sub-structures similar to those observed are sensitive to the recent merger histories of their host galaxies. Models predict that a survey of ∼100 parent galaxies reaching a surface brightness limit of ∼29 mag arcsec−2 would reveal many tens of tidal features, with approximately one detectable feature per galaxy (Johnston et al. 2001). One of the core motivations for the survey we describe here is that the statistical comparison of these detailed simulations with observations has not yet been possible, because no suitably deep data are available for a representative sample of galaxies.

The primary reason for the highly incomplete observational portrait of dwarf satellite accretion events is the inherent difficulty of detecting the low surface brightness features they create, even in the local Universe. Only a few instances of extragalactic stellar tidal streams had been reported around other nearby galaxies before the beginning of this century. Using contrast enhancement techniques on deep photographic plates, Malin & Hadley (1997) highlighted two possible tidal streams surrounding the galaxies M83 and M104. In the last decade, the way forward has been demonstrated by a programme of deep, wide-field images of nearby Milky-Way analogue galaxies in the local volume taken with small telescopes (the Stellar Tidal Stream Survey, STSS; Martínez-Delgado 2019). This survey has revealed, in many cases for the first time, an assortment of large-scale tidal structures around nearby spiral galaxies, with striking morphological characteristics that are qualitatively consistent with the predictions of N-body simulations (Johnston et al. 2008; Martínez-Delgado et al. 2010, see their Fig. 2). Further deep imaging surveys of the outskirts of local galaxies have been completed over the past decade (Tal et al. 2009; Ludwig et al. 2012; Duc et al. 2015; Laine et al. 2014; Rich et al. 2019) or are currently ongoing (Kado-Fong et al. 2018). Resolved stellar population studies of stellar halo substructure (similar to PAndAS, McConnachie et al. 2009, 2018) are only feasible for a select number of nearby galaxies within a distance of a few Megaparsecs using very large telescopes (Mouhcine & Ibata 2009; Crnojević 2017) or from space (Mouhcine et al. 2010; Crnojević et al. 2016; Carlin et al. 2019).

Many of the studies to date have targeted spiral galaxies with known features visible in shallower amateur imaging or the Sloan Digital Sky Survey (SDSS, Abazajian et al. 2009). However, subjective target selection criteria, such as the morphological type of the host galaxy, make it challenging to compare the extremely heterogeneous samples currently available with one another, or with predictions from cosmological simulations. Therefore, a meaningful comparison between data and models demands deep imaging data for a large number of galaxies with a simple, quantitative selection function, based on fundamental quantities such as the stellar mass of the host galaxy. Although several studies listed above have taken important steps towards this goal, there are significant differences between their results, likely due to different sample selections, detection criteria and surface brightness limits (see Hood et al. 2018, for a recent discussion).

In this paper we present the new Stellar Stream Legacy Survey. This survey focuses on stellar tidal features formed during the accretion of low-mass dwarf galaxies orbiting host galaxies of luminosity comparable to the Milky Way. For brevity, we refer to all such features as stellar tidal streams, or simply streams, regardless of their detailed morphology. This definition excludes tidal tails and other structures arising from mergers with a characteristic mass ratio greater than 1:10, in which the gravitational potential of the host is strongly perturbed. In contrast to the study of major merger remnants, the systematic study of ‘mini-merger’ events (with a characteristic mass ratios of ∼1:50–1:100) is a relatively new area of exploration.

Due to their extremely faint surface brightness, clear examples of stellar streams are rarely detected even in ultra-deep imaging surveys. In a recent systematic search combining SDSS DR9 images in multiple bands, Morales et al. (2018) detected streams around only ∼10% of a mass-selected sample of isolated Milky Way analogues (see Sect. 2.2). Moreover, the majority of previously known streams were only barely visible in the deep images constructed by Morales et al. (2018), making it infeasible to analyse their photometric and structural properties from SDSS data alone. The exploration of a few thousand galaxies with deeper data is needed to obtain the uniform, statistically representative samples required for meaningful comparisons to theory. Undertaking this type of survey with small robotic telescopes would be very expensive, given the substantial fraction of non-detections even in the ultra-deep imaging (∼29.0 mag arcsec−2) obtained with such facilities (see Sect. 3.4).

The promising results of previous forays into this topic argue strongly that a comprehensive study of tidal streams in the nearby universe provide a new perspective on galaxy formation. With the Stellar Stream Legacy Survey we aim to collect a systematic deep, broad-band images of stellar streams around a much larger sample of galaxies than any previous survey of such features. This is now possible thanks to the latest generation of large-area imaging surveys. The quality of these data offers a unique opportunity to reveal as-yet hidden details of galaxy assembly in the local Universe and to make meaningful connections with theoretical predictions from the ΛCDM paradigm.

In Sect. 2 we describe the scientific motivation and design of our Stellar Stream Legacy Survey, starting with an overview of the underlying data from the DESI Legacy Imaging Surveys (Dey et al. 2019) in Sect. 2.1. In Sect. 2.2 we describe how we extend the reduction pipeline of the Legacy Surveys to produce image cutouts and sky background subtraction suitable for the discovery and photometry of very low surface brightness features. In Sect. 2.3 we explain how we select a parent sample of target galaxies in the Local Universe. To validate our approach to discovery and photometry, in Sect. 3 we present a photometric analysis of 24 stream candidates spanning the range of distance, surface brightness and morphology in our parent sample. This section introduces our baseline approach to aperture measurements of surface brightness and colour. Even our small proof-of-concept sample adds significantly to the number of robust colour measurements for low-mass streams in the literature; we briefly discuss our findings in the context of the colour distribution of dwarf satellite galaxies. Sections 3.2 and 3.3 discuss sources of confusion associated with stream detections, with which any future statistical analysis will have to contend. Having estimated our surface brightness limit, we compare our survey with others in Sec. 3.4. In Sect. 4 we look ahead to improvements in the theoretical techniques that are necessary to exploit our full sample. We pay particular attention to requirements for meaningful comparisons with cosmological simulations, and show two examples of how state-of-the-art simulations can be used to produce mock images tailored to our survey. We conclude in Sect. 5. All magnitudes are in the AB system (Oke & Gunn 1983) unless otherwise noted.

2. The Stellar Stream Legacy Survey

2.1. The DESI Legacy Imaging Surveys

The Stellar Stream Legacy Survey is based on a custom re-reduction of public data from the DESI Legacy Imaging Surveys project (Dey et al. 2019)2, which combines images from several new wide area photometric surveys, by using a single analysis pipeline. This provides high quality multi-colour imaging with significantly greater depth and resolution than previous surveys such as SDSS or Pan-STARRS. The diffuse light feature detection limit of co-added images from the DESI Legacy Imaging Surveys is similar to that previously obtained with amateur telescopes (Martínez-Delgado et al. 2010, 2015), but additionally provides broadband photometry in standard filters3 (see Sect. 3.4) and sub-arcsec seeing, which allows us to trace fine details of streams at larger distances. The deep images shown in this paper are based on data release 8 (DR8) of the DESI Legacy Imaging Surveys, although future work will exploit the slightly greater coverage, depth and improved reduction of the more recent DR9.

Figure 1 shows in blue the sky coverage of the DESI Legacy Imaging Surveys. The combined footprint of its constituent surveys covers approximately 20 000 deg2 of the extra-galactic sky, approximately bounded by −60° < δ < +84° in celestial coordinates and |b|> 18° in Galactic coordinates.

thumbnail Fig. 1.

The Stellar Stream Legacy Survey sky coverage. Top: The DESI Legacy Imaging Surveys footprint (blue), the source of optical data for our survey. The 2MASS Redshift Survey catalogue (Huchra et al. 2012) beyond the local group and below VLG < 7000 km s−1 is plotted over the whole sky with a marker size scaled to the K-band magnitude. Those galaxies matching our selection criteria (2000 < VLG < 7000 km s−1, K-band absolute magnitude MK < −19.6) are displayed in bright red colour. The Stellar Stream Legacy Survey sub-samples comprise approximately 3000 of these galaxies in bright red within the blue area. Low galactic latitude regions of the survey affected by the higher stellar density of the Galaxy are perceived with this Gaia overlay. The galactic referential is plotted in green. Bottom: E(B − V) distribution across the DESI Legacy Imaging Surveys footprint (Planck Collaboration XI 2014). Areas of greater dust density (dark blue) are also contaminated with the Galactic light diffused by the dust grains (galactic cirrus) causing a slightly higher background and occasional overlapping structures at scales similar to the faint streams motivating this study. The same 2MRS galaxies are represented on this second map with a different colour scheme.

The DESI Legacy Imaging Surveys comprise data in three optical bands (g, r and z) coupled with all-sky infrared imaging from Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010; Meisner et al. 2019). The optical data were obtained by separate imaging projects on three different telescopes: The DECam Legacy Survey (DECaLS), the Beijing-Arizona Sky Survey (BASS) and the Mayall z-band Legacy Survey (MzLS; Zou et al. 2019; Dey et al. 2019). The DECaLS and MzLS surveys, and the data reduction efforts, were undertaken to provide targets for the Dark Energy Spectroscopic Instrument (DESI) survey (DESI Collaboration 2016a,b). The DESI Legacy Imaging Surveys data releases also include re-reduced public DECam data from other projects, including the Dark Energy Survey (DES; Abbott et al. 2018).

DECaLS provides complete g, r and z-band DECam imaging with excellent seeing over 6200 deg2 of the SDSS/BOSS/eBOSS footprint (Smee et al. 2013; Dawson et al. 2016), covering both the North Galactic Cap region at δ ≤ +32° and the South Galactic Cap region at δ ≤ +34°. Both DECaLS and DES data were obtained with the Dark Energy Camera (DECam) mounted on the Blanco 4-m telescope, located at the Cerro Torrolo Inter-American Observatory (Flaugher et al. 2015).

In the Northern Hemisphere, the DESI Legacy Imaging Surveys comprise images from the BASS and MzLS surveys. These surveys have imaged the Dec. ≥ + 32° footprint of the upcoming Dark Energy Spectroscopic Instrument (DESI) survey, as a complementary program to DECaLS. BASS provides g and r coverage using the 90Prime camera at the prime focus of the Bok 2.3 m telescope located on Kitt Peak, adjacent to the 4-m Mayall Telescope. The 90Prime instrument is a prime focus 8k × 8k CCD imager, yielding a 1.12 deg field of view, with a scale of 0 . ″ $ \overset{\prime \prime }{.} $45 pixel−1 (Williams et al. 2004). The survey has g- and r-band photometry over 5000 deg2. The BASS survey tiled the sky in three passes, similar to the DECaLS survey strategy. At least one pass was observed in both photometric conditions and seeing conditions better than 1.7″. The Mayall z-band Legacy Survey (MzLS) imaged the same sky area with a similar observing strategy, using the MOSAIC-3 camera (Dey et al. 2016) at the prime focus of the 4 m Mayall telescope at Kitt Peak National Observatory, in photometric conditions and with seeing better than 1.3 arcsec.

2.2. Image cutouts and sky background subtraction

The unique low surface brightness demands of our tidal stream survey requires bespoke sky-subtraction and co-addition of the calibrated survey images in the vicinity of each target galaxy. Here, we briefly describe our custom analysis, which utilises large portions of the legacypipe4 software infrastructure, the dedicated pipeline developed by the DESI Legacy Imaging Surveys team for source detection and model fitting of this survey imaging. One advantage of utilising legacypipe is that our custom analysis produces a standard and well-documented set of imaging and catalogue data products. Relative to the standard Legacy Survey co-added image products, these bespoke images reduce sky subtraction systematics around brighter host galaxies and so provide a more robust and uniform starting point for the visual detection of faint features, automated classification and quantitative photometry.

We begin from the photometrically and astrometrically calibrated CCD images produced by the Legacy Surveys imaging team (see Dey et al. 2019, for details). Next, for the predefined custom footprint of each galaxy in our sample, we read the set of CCD exposures contributing to that field. Typically there are many separate partially overlapping CCD frames for a given target; an example is shown in Fig. 2. Then, we subtract the sky background from each CCD using a custom algorithm, described below, which preserves the low surface-brightness galactic features of interest. Finally, we use Lanczos3 (sinc) resampling to project each CCD onto a common (predefined) tangent plane and pixel scale ( 0 . ″ 262 $ 0{{\overset{\prime\prime}{.}}}262 $ pixel−1), and combine all the data using inverse variance weights.

thumbnail Fig. 2.

Example of the deep mosaics we construct from the individual DESI Legacy Imaging Surveys exposures overlapping our target galaxies. From left to right, panels show g, r, and z bands. The white rectangle shows the region we process with our custom pipeline around a given target (here NGC 2460), coloured rectangles show each individual CCD frame in the stack. The overlapping tiling pattern is visible. Since this galaxy is in the northern footprint of the DESI Legacy Imaging Surveys, the g and r bands are from BASS, where the CCD pixel scale is larger than in the z band, which is from MzLS. The right-most panel shows the final mosaic.

Our custom background subtraction algorithm consists of the following steps: We first bring all the input imaging in a given bandpass to a common (additive) pedestal by solving the linear equation Ax = b for the full set of N images and M overlapping images, where A is an M × N matrix of the summed variances of overlapping pixels, x is an N-length vector of the additive offsets we must apply to each image (in calibrated flux units), and b is an M-length vector of the summed inverse-variance weighted difference between the pixels of overlapping images. We solve this equation using standard linear least-squares and then add the derived constant offsets to each CCD image. Next, we resample and coadd the data onto a common tangent plane and measure the median sky background (after aggressively masking sources) in an annulus defined to be 2 − 4 times the masking radius of each galaxy. We define an unique masking radius for each galaxy based on its approximate projected angular size (in arcseconds) on the sky, although our results are not sensitive to the exact choice of masking radius since our background annulus is well outside the outer, low surface-brightness envelope of the galaxy.

2.3. Survey design and sample selection

To facilitate comparisons with simulations, we define a straightforward sample of isolated galaxies based only on luminosity and recessional velocity given by the HyperLeda database5 (Makarov et al. 2014). We consider galaxies with line-of-sight velocities in the Local Group rest frame (Karachentsev & Makarov 1996) in the range 2000 < VLG < 7000 km s−1 (30 ≲ D ≲ 100 Mpc), with a K-band absolute magnitude MK < −19.6 and within the footprints of the surveys described in Sect. 2.1.

The majority of K-band magnitudes are taken from the 2MASS survey (Skrutskie et al. 2006) and corrected for foreground extinction using reddening values from Schlafly & Finkbeiner (2011), assuming RV = 3.1. No colour or morphological selections are imposed on our sample. We exclude targets that lie close to the Zone of Avoidance (Galactic latitude |b|< 20 deg) to avoid the area severely contaminated by Galactic dust and high stellar densities.

To select ‘isolated’ galaxies we use a simple criterion which requires that there is no neighbour brighter than Kgal = 2.5 mag with |ΔV|< 250 km s−1 within a projected radius of 1 Mpc around each target. This criterion includes the brightest galaxies of ‘fossil groups’ but excludes strongly interacting groups of galaxies with comparable masses. The criterion is based on a typical turn-around radius of 1 Mpc for the Local Group and other nearby groups such as M81 and Centaurus (Karachentsev et al. 2009). The typical velocity dispersion of such groups is 70–80 km s−1 (Karachentsev 2005); we use a limit of 250 km s−1 as a 3σ threshold. Applying these criteria to the HyperLeda database yields a total sample of ∼3100 target galaxies.

The luminosity range of our sample ensures that the target galaxies have a broad range of morphologies, star formation rates and stellar masses. This selection includes ∼940 Milky-Way analogue systems, which are defined by their K-band luminosity (−24.6 < MK < −23.0) and the local environment criteria given in Geha et al. (2017). Our lower distance bound of ∼30 Mpc ensures that the virial radius of a typical galaxy in the sample (∼250–300 kpc) can be covered with a field of view (FOV) of 30′×30′ centred on a target (see Sect. 2.2). Furthermore, the image quality from the DESI Legacy Imaging Surveys is sufficiently high to discover stellar streams even at the upper distance bound of our survey, ∼100 Mpc. Galaxies at that distance require only 10′×10′ cutouts. This outer boundary of the survey volume was chosen to obtain a sample containing a statistically significant number of tidal features based on the detection frequency in shallower surveys (Morales et al. 2018).

3. Data analysis and photometry

To date, we have visually inspected a total of 689 cutouts from all the targets included in the DES sky footprint, yielding a sample of 89 candidate streams. Some additional conspicuous stellar streams have also found in an incomplete visual inspection of the target cutouts from the DECaLs and MsLS sky footprints. In this paper, as a proof of concept, we present a quantitative study of the detection significance and photometric properties for a representative sub-sample of 24 of those candidates. Images of this sub-sample are shown in Figs. 3 and 4, and the properties of their host galaxies are listed in Table 1. These features are reported for the first time here. The sub-sample was chosen to span a representative range of morphological types, surface brightnesses and distances in the DES and DECaLs sky footprints. Our aim in the following analysis is to asses the detection limits of our survey and the viability of our baseline photometric approach. A catalogue of candidates for the complete sample of hosts are published in a forthcoming paper (Miró-Carretero et al., in prep.).

thumbnail Fig. 3.

Stacked g and r image cutouts of 12 stellar streams discovered during the proof-of-concept phase of the Stellar Stream Legacy Survey introduced in this paper. For illustrative purposes, colour insets of the central region of the host galaxies have been added to the negative version of the images.

thumbnail Fig. 4.

As Fig. 3, but showing an additional 12 new stellar streams from our Stellar Stream Legacy Survey.

Table 1.

List of galaxies with stellar streams discovered during the proof-of-concept phase of our new survey.

3.1. Detection and photometry of stellar tidal streams

Our strategy for searching for stellar streams is mainly based on visual inspection of DESI LS imaging of the outskirts of the target galaxies. Streams are inherently non-symmetric and elongated, making them extremely hard to distinguish from the noise with classical detection methods. Therefore, we use NoiseChisel (Akhlaghi & Ichikawa 2015), a component of the GNU Astronomy Utilities 0.14 (Gnuastro6), to confirm detections and estimate the significance of each detected structure in our survey.7

NoiseChisel uses thresholds that are below the sky level to detect extremely faint signals without making assumptions about the shape or surface brightness profile of the sources. It then uses simple mathematical morphology operators (erosion and dilation) to extract a contiguous signal at low surface brightness. For example Akhlaghi (2020) describe a successful application to detect the outer wings of M51 to a surface brightness of 25.97 mag arcsec−2 (r-band) in a single ∼1 min SDSS exposure.

Figure 5a shows an example of NGC 4981 in the z-band from our dataset8. All pixels that do not belong to NGC 4981 and its stream have been masked in panel (b), as well as all foreground and/or background sources overlapping the stream. Line-of-sight contamination was identified as clumps over NoiseChisel’s detections, using Segment9 (Akhlaghi 2020, also part of Gnuastro) and masked in Fig. 5b. Figure 5c shows the host and the stream with all masked clumps interpolated.

thumbnail Fig. 5.

Detecting and measuring streams. (a) The input image: NGC 4981 in z band. The red aperture has a radius of 7 arcsec and lies on the thinnest and most diffuse end of the stream. (b) All pixels not belonging to the galaxy (detected by NoiseChisel) or its stream are masked, including the pixels of objects that are in the background or foreground of the stream (‘clumps’ in Segment). (c) All masked pixels interpolated by the median of ten neighbors. (d) The flux in the red aperture (dotted line), compared to a distribution of 10 000 similarly sized apertures, randomly placed over undetected pixels of the image (4.32σ). (e) The input dataset warped to a pixel grid where each pixel is 8 × 8 larger than the input’s pixels, to demonstrate the detection success on the input resolution.

The 7 arcsec radius aperture shown in (a) is positioned along the thinnest and faintest end of the detected stream at a point with no background or foreground source. The flux within this aperture has a surface brightness of 26.81 mag arcsec−2, which has a significance of 4.32σ compared to 10 000 random apertures of the same size placed over undetected regions (see panel (d)). As the area increases, the significance of the detection improves, enabling us to extract meaningful colour measurements. However, we are not limited to measuring colours in simple aperture geometries. Because of Gnuastro’s modularity (which separates image labelling, through detection or segmentation, from measurements and catalogue production), magnitudes and colours can be measured on any feature with arbitrary morphology, like the stellar streams (Akhlaghi 2019).

Another factor to consider is that the detectability of these features is closely related to the surface brightness limit of the image cutouts obtained in Sect. 2.2. This limit depends on the region of sky in which the galaxy is located, and in general on the exposure time and observing conditions. In this survey, the surface brightness limit of the images was calculated for the three bands following the standard method in Appendix A of Román et al. (2020), as implemented in Gnuastro’s MakeCatalog10, that is the measured 3σ value in an area of 100 arcsec2.

A high degree of accuracy in data reduction and post-processing is required to measure surface brightness, colours and integrated luminosities for extremely faint and diffuse features. The procedures detailed in previous sections provide a necessary improvement in the stability of the sky background estimate. However, other systematic uncertainties still need to be addressed. For example, some tidal features overlap with their host galaxy along the line of sight. Deblending these features to measure directly their total luminosity requires a precise model of the host galaxy. This is a challenging problem which we leave for future work.

We carried out a photometric analysis on the set of streams identified in our proof-of-concept sample by using Gnuastro’s MakeCatalog on the sky-subtracted image generated by NoiseChisel (Akhlaghi & Ichikawa 2015). The surface brightness of each stream was measured in the r, g and z bands, together with the colours (g − r) and (r − z). The measurements were performed using circular apertures placed on the detected track of the stream, after all foreground and background sources were removed. Regions where the stream surface brightness could be contaminated by the host galaxy were avoided. This method is illustrated in Fig. 6, which shows an example application to a tidal stream around NGC 3451. The average surface brightness of this stream is only 1.24 mag arcsec−2 brighter than the surface brightness limit of the image.

thumbnail Fig. 6.

Example of our aperture photometry approach applied to the stellar stream detected around NGC 3451. Top panel: input image in g-band, with subtracted sky, showing the stream around it. Middle panel: the detection contour of the galaxy and its stream are highlighted (black line) on a pixel grid 8 × 8 courser than the input image to highlight the low surface brightness signal. Bottom panel: the apertures used to measured the surface brightness, magnitude and colours are placed along the stream (red circles). The foreground stars and background sources are masked (white areas).

The resulting average surface brightnesses and colours of the stellar streams from our proof-of-concept sample are given in Table 2 along with the corresponding errors, computed by Gnuastro’s MakeCatalog11. This table also includes the surface brightness limit for the r band, which is representative of the depth of the corresponding images.

Table 2.

Photometry of stellar streams discovered in our proof-of-concept analysis for a sub-sample of Stellar Stream Legacy Survey targets, listed in Table 1.

In order to quantify the significance of individual stream detections, we define the Detection Index, DIstream = (Fstream − Fblank)/σ. We note that Fstream is the flux measured in a single aperture on the stream, Fblank is the median flux measured in 10 000 apertures of the same size randomly placed in regions of the sky-subtracted image with no detection12, and σ is the standard deviation of the same 10 000 background flux measurements. Columns 3 and 4 in Table 2 report the maximum and mean value of DIstream, respectively, taken over all the apertures placed on each stream.

Figure 7 shows distributions of the surface brightness limit of our images and the surface brightnesses and colours of the stellar streams identified around the galaxies listed in Table 2. Taking the colour distributions in Fig. 7 at face value, it appears that the g − r colours of streams are systematically redder than those of typical dwarf satellite galaxies. For example, the satellite galaxies around Milky Way analogues observed in the SAGA survey13 have a range 0.2 < g − r < 0.7 (Mao et al. 2021), whereas we find stream colours 0.5 < g − r < 0.8 around visually similar hosts. Under the simplistic assumption that the colours measured in SAGA are representative of satellite dwarf galaxies up to the point which their star formation is truncated by tidal disruption, after which they redden, our result may imply that the stream population is, on average, dynamically older than surviving luminous satellites. Redder colours may also imply that progenitors of the streams we detect in our survey are more massive than typical surviving satellites. Both interpretations would be consistent with theoretical expectations (e.g., Cooper et al. 2010). However, a robust statistical analysis requires both the larger sample of our complete survey and a careful exploration of the systematic uncertainties on our colour measurements, in addition to a more detailed study of model predictions.

thumbnail Fig. 7.

Histograms showing the distribution of the depth of the image cutouts and of the photometry measurements of the stellar streams around the galaxies listed in Table 2. From top to bottom: Surface brightness limit for the r band, surface brightness for the r band, and colours (g − r) and (r − z).

3.2. Distinguishing stellar streams from other features

As stated in Sect. 1, our aim is to study dwarf galaxy accretion around nearby galaxies that have not been strongly perturbed by major merger events in the recent past. Therefore we need to develop a strategy to identify and exclude from our sample features associated with the partially relaxed debris of major mergers and perturbations to a central stellar disk that are visually similar to stellar streams.

3.2.1. Old major-merger remnant systems

Our isolation criterion should exclude the majority of massive galaxies and coalescence phases of major mergers. It should also remove from our sample very bright galaxies undergoing tidal disruption within massive galaxy clusters. However, galaxies in the late stages of major mergers may still enter our parent sample. Such ‘post-merger’ systems may have a well-relaxed core surrounded by a large number of faint remnant features generated by past tidal interactions between the progenitors14.

It may be possible to disentangle these cases by exploiting the morphology of their tidal tails, in particular their width. The width of tidal features scales with the progenitor mass, m, as m1/3 (e.g., Johnston et al. 2001; Erkal et al. 2016). Older, more massive mergers have significantly broader features. A similar distinction was used in Johnston et al. (2001) to classify debris.

We demonstrate how this classification would work by running simulations of the tidal debris from different mass progenitors on a range of orbits. These simulations were carried out with the modified Lagrange Cloud stripping technique from Gibbons et al. (2014) as implemented in Erkal et al. (2017). This technique faithfully reproduces the debris of an N-body disruption by ejecting particles from the Lagrange points of the progenitor, allowing us to model the disruption in a few CPU-seconds instead of a few CPU-hours. The progenitors are modelled as Plummer spheres following the size–luminosity relation for observed local group dwarfs (Brasseur et al. 2011). The host galaxy potential was modelled using the MWPotential2014 from Bovy (2015) which provides a good match to the Milky Way potential. For computational simplicity, we replace the bulge with a Hernquist profile (Hernquist 1990) with a mass of 5 × 109 M⊙ and a scale radius of 0.5 kpc.

Examples from our test simulations are shown in the right panel of Fig. 8. It is clear that more massive progenitors create broader tidal debris streams. In addition, as the orbit is made more radial, the debris transitions from being stream-like to being shell-like (e.g., Hendel & Johnston 2015). The highest mass considered here corresponds to a major merger with a total mass ratio of ∼1 : 3. We have used abundance matching for the stellar mass-halo mass relation (e.g., Behroozi et al. 2013; Moster et al. 2013), and evolved the tidal debris for 2 Gyr.

thumbnail Fig. 8.

Left panels: examples of likely major merger remnants in our proof-of-concept sample (see Table 1): ESO474-26, NGC 655, NGC 3303 and NGC 7288. These were identified as massive merger remnants based on the similarity to the outcome of idealised simulations: see right panels. Right panels: simulation snapshots of galaxy mergers evolved for 2 Gyr (see Sect. 4.1). Each corresponds to different initial conditions: rows show different orbital eccentricities given by their circularity L/Lc as defined by Hendel & Johnston (2015); columns show different progenitor stellar masses. The 5 × 109 M⊙ progenitor corresponds to a major merger with a mass ratio ∼1:3; the other two models are for minor mergers. Progenitors were initialised on the size–luminosity relation of Local Group satellites (McConnachie 2012). The colour bar shows the surface brightness of the tidal debris in mag arcsec−2.

The simulations suggest a correlation between stream width, halo mass and dynamical age. We find the thinnest and most coherent structures are generated on the early phases of the interaction, when the central galaxy shows the strongest morphological distortions. This phase is also expected to generate starbursts and AGN activity that clearly distinguish merging systems (Springel & Hernquist 2005; Lotz et al. 2008; Scott & Kaviraj 2014; Pearson et al. 2019a). The debris is broader, and broadens more rapidly, than in a minor merger. In minor mergers, by definition, the central potential is not significantly perturbed, and the central galaxy remains almost quiescent with a small star formation enhancement and no AGN activation (Gordon et al. 2019). For comparison, the left-hand panel of Fig. 8 shows likely major merger remnants identified in our proof-of-concept sample. Quantitative comparison of a larger set of simulations and observations of partly-relaxed mergers is useful to develop more robust estimators of progenitor mass.

3.2.2. Galactic feathers

Another possible source of confusion when classifying streams could come from mis-classification of perturbed stellar disc material in the form of tidal tails (Laporte et al. 2019) that can be excited through interactions with low-mass companions (Toomre & Toomre 1972).

Minor mergers in late-type galaxies can excite various structures such as rings and tidal tails resembling structures seen at the disc-halo interface in the Milky Way (Kazantzidis et al. 2009; Slater et al. 2014; Gómez et al. 2016; Bauer et al. 2018). Among such structures showing stream-like morphology on the sky, we note the eastern banded structure (EBS) and the anticenter stream (ACS; Grillmair 2006), which have long been suspected to be tidal debris of accreted satellites or massive globular clusters (e.g., Peñarrubia et al. 2005; Amorisco 2015). However, recent N-body models of the interaction of a Sagittarius-like (Sgr) dwarf with the Milky Way have presented an alternative picture (Laporte et al. 2018). In this model, structures like the ACS and EBS were re-interpreted as tidal tails in the outer disc resulting from the interaction with Sgr which, with later passages, were excited to larger heights, giving them the appearance of thin stream-like structures of disc material or “feathers” (Laporte et al. 2019). A recent analysis of the chemical and age properties of the ACS supports this re-interpretation (Laporte et al. 2020). Due to the long orbital timescales in the outer disc (torb ≥ 1 Gyr), these structures can survive and remain coherent for long timescales of the order ∼ 6 Gyr or more and get gradually excited to larger heights above the Galactic midplane (z ∼ 30 kpc, making these potentially relevant sources of confusion as tidal streams in the stellar halo (Laporte et al. 2019).

In some instances, ongoing mergers leave behind detectable signatures, such as an easily identifiable remnant core if the central galaxy is viewed at a favourable angle. However, in other merging configuration, it may be difficult to pinpoint whether a tidal stream could be a feather or the remains of a shredded satellite. Stars from the host can be distinguished by their ages and metallicities, inferred through their resolved stellar populations (for example, from the ratio between M-giants and RR Lyraes, Price-Whelan et al. 2015), although diagnostics that do not rely on spectroscopy or resolved star counts are very limited and have not been explored. In the case of the Milky Way, the ACS has a median metallicity [Fe/H] ∼ − 0.7 and magnesium abundance [Mg/Fe], placing it at the α-poor end of the distribution for thin disc stars with ages of τ > 8 Gyr (Laporte et al. 2020), consistent with this structure originating from the disc. However, provided one has a constraint on the host’s stellar halo mass, (M⋆ ∼ 109 M⊙ in the Milky Way (see Deason et al. 2019), metallicity alone can suffice to rule out an ex-situ origin for substructures through the mass-metallicity relation of dwarf galaxies (Kirby et al. 2013). Thus, in the case of an external galaxy, provided colour information is available to derive robust photometric metallicities, it may thus be possible to distinguish streams from feathers using colour information alone in tandem with information on the total stellar halo mass of the system through deep exposures like those of this survey.

In addition to feathers, other diffuse overdensities may arise at the stellar halo-disc interface. An example is structures left behind by bending waves created during satellite encounters, which dissipate and flare the outer disc (Kazantzidis et al. 2009; Gómez et al. 2013; Chequers et al. 2018; Laporte et al 2018). Known cases in the Milky Way include the Monoceros Ring (Newberg et al. 2002; Slater et al. 2014) or TriAnd (Rocha-Pinto et al. 2004). Given their complex and amorphous shapes, such structures would be difficult to disentangle from minor merger events with photometry alone, and may ultimately require bespoke dynamical models of the outer discs of individual galaxies. Given the negligible self-gravity of the outer discs, it may be possible to use cost-effective modelling in the same spirit as Lagrangian Cloud Stripping methods (e.g., Gibbons et al. 2014) adapted to galactic discs.

In Fig. 9, we present a gallery of simulated feathers produced in Milky Way-like galaxies subjected to an ongoing merger with a 1:10 companion as presented in Laporte et al. (2018), viewed from different projection angles. It can be seen that, depending on the viewing angle, feathers can exhibit a wide range of complex morphologies which may lead to confusion with the features generated by dwarf galaxy disruption. In Fig. 10, we showcase possible candidates of ‘in situ’ feathers in our proof-of-concept sample.

thumbnail Fig. 9.

Gallery of simulated galactic feathers shown at different stages of interaction (different columns) and different viewing angles (rows). Top: edge-on view; middle: face-on view; and bottom: inclined view. Many of these structures show morphologies reminiscent of tidal streams from disrupting satellites although in this case this material is purely in-situ from the disc. A colour characterisation or metallicity measurement could potentially disentangle the two possibilities.

thumbnail Fig. 10.

Image cutouts of galactic feather candidates with diverse morphologies identified in the proof-of-concept phase of our survey (see Table 1): NGC 4981, NGC 5492, NGC 5635 and ESO245-010.

Although it remains to be seen if colour information alone is sufficient to reliable distinguish feathers from streams, our survey provides a useful test bed to systematically study in-situ structures at the disc-halo interface in external galaxies, and an opportunity to devise tests to disentangle them from genuine streams of accreted substructures (e.g., through residuals from direct image modelling).

3.3. Disentangling tidal features from Galactic cirri

In this section we describe a potentially useful test to identify false detections due to Galactic dust, which often mimics the morphology and surface brightness of tidal features in deep images. Our selection criteria reject galaxies in the zone of avoidance (Galactic latitudes |b|< 20deg). However, low-surface brightness, optically-thin dust clouds at relatively high Galactic latitudes are expected to be detected in our deep data.

Since both cirrus and diffuse features have a wide range of properties, it is helpful to use multiple diagnostics to distinguish them. First of all, visual analysis of the images is often sufficient. Most tidal feature morphologies are unambiguous and, together with the absence of structures of similar surface brightness in the field surrounding the host galaxy, allow us to rule out false positives. Also, detections of dust features in far infrared (FIR) data, such as those from Herschel (or in its absence, IRAS 100 μm or Planck 857 Ghz, Miville-Deschênes & Lagache 2005) can be used. Recent work by Román et al. (2020) provides another diagnostic to distinguish tidal streams from Galactic cirrus clouds, which makes use of the optical colours of the cirri in the g, r, i and z bands. Román et al. (2020) show that the colours of the cirrus clouds are significantly different from those of any galactic stellar population, being characterised by a bluer r − i colour for a given g − r colour. This method has the great advantage of being based on the optical observations themselves, independent of complementary IR data which may not be available in all cases. This diagnostic also benefits from the high spatial resolution in our optical images.

In Fig. 11, we plot the colours of the features described in Table 2 against the colours of Galactic cirri reported in Román et al. (2020). The figure shows that the colours of cirri are mostly clustered in a relatively compact region of g − r vs. r − z space. In contrast, we see that the colours of the tidal streams on this map occupy a much larger area. Although the two colour distributions partly overlap (relative to Román et al. 2020), the grz bands are not optimal for this separation), the separation is nevertheless significant enough to provide a useful estimate of the probability that a given feature is due to cirrus.

thumbnail Fig. 11.

g − r vs. r − z colours for the tidal streams presented in Table 2 (we exclude tidal features with colour errors > 0.25 mag for clarity) and the Galactic cirri regions characterised by Román et al. (2020). The black dashed line, (r − z) = 0.20 × (g − r)+0.20, limits the regions where the colours are compatible/incompatible (below/above of the line) with Galactic dust.

Colours (r − z) > 0.20 × (g − r)+0.20 (see Fig. 11) are highly likely to indicate diffuse extragalactic stellar populations rather than Galactic dust. We caution that, although potentially useful, this method is by no means definitive. It can be applied as a last resort in the absence of additional data, but it is most powerful in combination with other more complex approaches.

3.4. Comparison with existing and future ultra-deep imaging of nearby galaxies

The analysis of low surface brightness observations is hindered by the lack of consensus on the most appropriate method for measuring the surface brightness limit of an astronomical image (Mihos 2019). One of the most common methods is to measure the width of the background noise distribution over pixels of a fixed size (e.g., Trujillo & Fliri 2016). We follow this approach throughout, by quoting surface brightness limits (i.e. residual sky backgrounds, depths of our mosaics) as 3σ of the surface brightness distribution measured over an area of 10″ × 10″. This method only takes into account random statistical noise, not systematic errors. Images suffering from severe over-subtraction of the sky background have the same surface brightness limit as well-calibrated images according to this definition (Borlaff et al. 2019).

In Fig. 12 we compare the surface brightness magnitude limit (calculated as above, and averaged over all the available filters for each field or survey) as a function of the observed area for several current and planned deep imaging surveys. Extremely deep surveys like the Hubble Ultra Deep Field (Illingworth et al. 2013; Borlaff et al. 2019) have a much smaller area (a few arcmin2) compared to extremely wide but shallower surveys like SDSS (York et al. 2000). For the time being, depths below 30 mag arcsec−2 are incredibly hard to achieve with ground-based telescopes (Trujillo & Fliri 2016) and prohibitive for a large number of galaxies. Figure 13 shows that there are no striking differences in limiting depth between images taken with amateur telescopes from the Stellar Tidal Stream Survey (Martínez-Delgado 2019) and the results of our custom pipeline for the new Stellar Stream Legacy Survey applied to data from the DESI Legacy Imaging Surveys. Ten minutes of computing time to reprocess images from the latter can achieve the same depth as ∼10 h of dark observing time with amateur telescopes for each galaxy. A survey using small telescopes equivalent to compiling a sample of ∼3000 galaxies out to 100 Mpc would need ∼20 000 h of dark time. These results support our strategy of using archival DESI Legacy Imaging Surveys data, which also have superior seeing.

thumbnail Fig. 12.

Comparison of the surface brightness limit (3σ, 10 × 10 arcsec2) as a function of the observed area for a selection of deep optical and NIR surveys. Depths are measured in the r-band (noted when not applicable). The Stellar Stream Legacy Survey presented in this paper is indicated with a blue star. The other surveys comprise: Stellar Tidal Stream Survey (Martínez-Delgado 2019, STSS), Euclid VIS (Euclid Collaboration 2022a), SDSS (York et al. 2000), IAC Stripe 82 (Fliri & Trujillo 2016), the MATLAS deep imaging Survey (g-band magnitudes Duc et al. 2015), Hyper Suprime-Cam DR2 (Aihara et al. 2018), NGC 4565 Dragonfly observation (Gilhuly et al. 2020), HST WFC3 ABYSS HUDF (NIR bands Borlaff et al. 2019), XDF (Illingworth et al. 2013), UGC00180 10.4m GTC exploration (Trujillo & Fliri 2016), the Burrell Schmidt Deep Virgo Survey (Mihos et al. 2017), the CFIS/UNIONS survey (Sola et al. 2022), the VEGAS/VLT survey (Ragusa et al. 2021), and the Rubin Telescope Legacy Survey of Space and Time (10 yr full survey integration, Ivezić et al. 2019).

thumbnail Fig. 13.

Comparison of the image cutouts obtained by applying the methodology described in this paper to data from the DESI Legacy Imaging Surveys (top panels) versus long exposure times (6–10 h) images obtained with amateur robotic telescopes from the Stellar Tidal Stream Survey (Martínez-Delgado 2019, bottom panels): our new Stellar Stream Legacy Survey is well positioned to uncover hundreds of low surface brightness features in the local Universe.

A number of other deep imaging surveys are expected to approach 30 mag arcsec−2 in the next decade, including the space-based 1.2-m visual to NIR Euclid space telescope, the ground-based visual wavelength 8.4-m Vera C. Rubin Observatory (Ivezić et al. 2019, hereafter Rubin) and the space-based 2.4-m, NIR Nancy Grace Roman Space Telescope. The expected sensitivities and other basic parameters are compared to the Legacy Survey in Table 3. Points corresponding to the depths of the Rubin and Euclid surveys (according to our standard definition) appear in the low right corner of Fig. 12.

Table 3.

Comparison between the Legacy Survey and future surveys.

Although these future surveys do not cover a significantly larger area than our survey, they will eventually provide a significant improvement in depth and spatial resolution. Greater depth allows the detection of fainter stellar streams arising from ultra-faint satellites and globular clusters (see, e.g., Grand et al. 2017). Higher spatial resolution permits more accurate masking of distant background galaxies, and hence aids the recovery of thin stellar streams from globular clusters (see Pearson et al. 2019b). Improved masking of superposed sources could also increase the accuracy of colour (hence age and metallicity) measurements for streams. Deep NIR images improve stellar mass estimates and population constraints for streams, by sampling the peak of their SEDs.

4. Outlook and comparison with cosmological simulations

In the previous sections we have described the ideas and techniques of our new Stellar Stream Legacy Survey and we have shown proof-of-concept results for 24 new stellar streams detected in the local Universe. In this section, we discuss how the depth, coverage and uniformity of the full survey dataset can be exploited to address two key methodological and scientific problems related to stellar streams that we plan to undertake in future: automatic detection of stellar streams, morphological classification and comparison to cosmological simulations.

4.1. Automatic detection of stellar streams with Deep Learning algorithms

We complement our visual identification of tidal streams by employing convolutional neural networks (CNN; Simonyan & Zisserman 2014) to detect stellar tidal streams. CNNs, a form of deep learning (e.g., Lecun et al. 2015), take multi-dimensional data as an input, such as an image, and perform a series of non-linear convolutions before outputting a classification. This technique has been increasingly used in astronomy in recent years to identify rare objects in imaging surveys, such as quasars (e.g., Pasquet-Itam & Pasquet 2018), strong lenses (e.g., Petrillo et al. 2017; Davies et al. 2019) and galaxy mergers (e.g., Ackermann et al. 2018; Pearson et al. 2019d). Classification of mergers, in particular, relies on the identification of diffuse features (Pearson et al. 2019c), demonstrating that this technique is ideally suited to identifying stellar streams. CNNs also have the advantage of being reproducible: the same classification (stellar stream or no stellar stream, in this case) is reported every time the image is passed through the network.

Walmsley et al. (2019) presented the first attempt to identify faint tidal features around galaxies from the CFHTLS-Wide Survey using CNNs trained on a relatively small number of streams (305). The models achieve 76% completeness and 20% contamination. This relatively poor performance, compared to other classification tasks for galaxy images obtained with CNNs, may be due to the combination of the faintness of the tidal features combined with the lack of a large training sample.

To address the limitation imposed by the small size of the sample of known tidal streams, our survey trains the CNN using mock images generated from particle-spray models (see Sect. 4.2). Accurate classification of galaxy mergers with CNNs has been demonstrated using a few thousand merger examples (e.g., Ackermann et al. 2018; Pearson et al. 2019d; Bickley et al. 2021). As our objective is arguably more challenging than merger classification, we generate of the order of ten thousand mock tidal streams. These mock streams are combined with an identically sized sample of mock galaxies without streams, to provide a class-balanced training sample. As tidal streams are rare, a sample of galaxies that is matched to the true ratio of galaxies with and without tidal streams runs the risk that the CNN never learn to identify the features of interest. Performance metrics can be improved by assigning all objects to the majority class.

The mock images used for training are single channel to prevent any incorrect or non-physical colour information from the particle-spray model influencing the detection of streams, choosing the deepest band within each survey to maximise our sample size. Alternatively, it is possible to normalise each band independently, as done in Domínguez Sánchez et al. (2018) and Vega-Ferrero et al. (2021), which again remove any colour dependence of the classification. Once the CNNs are trained, they are able to identify the stellar streams in a fraction of the time it would take a human. Moreover, as demonstrated in Vega-Ferrero et al. (2021), machines are able to recover features hidden to the human eye when trained on a sample with correct labels (e.g., with labels coming from simulations). This technique allows the detection of faint tidal stream structures which may otherwise be missed by visual inspection.

Due to the differing angular resolutions, we train one CNN per survey. Training multiple CNNs has the benefit that it is possible to generate a sub-sample of very robust stellar streams in the regions where the surveys overlap. This is done by selecting galaxies that are identified as having stellar streams by the CNNs for two, or more, different surveys.

The performance of deep learning models trained with simulations and applied to real galaxies can be affected if the simulations do not capture the full complexity of the real images (e.g., noise, PSF, seeing and other observational effects). The mock images reproduce the processing steps applied to the real images as closely as possible. This includes applying a realistic noise model and convolving with the PSF of the survey in each band. Such techniques have also been demonstrated in the context of galaxy merger detection (Pearson et al. 2019c; Wang et al. 2020; Bickley et al. 2021). Domain adaption is also be used, a technique in which a CNN is paired with a Generative Adversarial Network to force the CNN to give identical, or at least similar, results for the simulation and real images (e.g., Wang & Deng 2018). Combining simulations and real images in this way help ensure the reliability of results obtained by applying the CNN to the real survey data. The sample of tidal streams identified by the Stellar Stream Legacy Survey is, in turn, useful for training of future CNNs for the mophological analysis of galaxies.

4.2. Morphological classification of stellar streams

Johnston et al. (2008) showed that estimates of the frequency, sky coverage and fraction of stars residing in substructure in stellar halos can be used as a diagnostic of recent merger events. In particular, they found that the surface brightness and morphology of substructure can constrain the masses and orbital distributions of recently accreted satellites. Hendel & Johnston (2015) found that the observed ratio between the number of great circle streams (formed by disrupted satellites on near-circular orbits) and the number of umbrellas/shells (formed by satellites on more radial orbits) can distinguish between the infall distributions predicted by different cosmological simulations, and that meaningful constraints could be obtained by an unbiased, uniform survey containing a few hundred stellar streams. Our survey yields even better statistics for substructure morphologies in the local Universe and hence place tighter constraints on the predictions of ΛCDM simulations. Figure 14 shows examples from our Stellar Stream Legacy Survey images compared to the theoretical morphological classification of tidal debris by Johnston et al. (2008).

thumbnail Fig. 14.

Bottom panels: external perspective snapshots for the three morphological types classified by Johnston et al. (2008): “Great circle” streams (left) originated from satellites accreted 6–10 Gyr ago on nearly circular orbits; “shells” or “umbrellas” (middle) arise from accretion events within the past ∼8 Gyr with radial orbits; and “mixed”-type tidal remnants (right) formed from ancient (more than 10 Gyr ago) accretion events that have fully mix along their orbits. Top panels: image cutouts obtained from the DESI Legacy Imaging Surveys data for each morphological case of tidal stream: a great circle stream in NGC 577 (left), shell-like features around NGC 681 (middle); and a “mixed” remnant around NGC 5866 (right). Figure adapted from Carlin et al. (2016).

However, we first need to characterise hundreds of debris structures in our survey. Hendel et al. (2019) developed a new machine vision technique which can automatically differentiate between shell-like and stream-like features through identification of density ‘ridges’ that define substructure morphology, which we can apply to our image data. More generally, our survey is well placed to make progress on techniques for the automatic morphological classification of streams. We therefore use two different approaches to characterise debris structure and streams in our sample:

Visual classification. Diffuse structures in stellar halos are typically classified into three morphological types: streams, umbrella/shells and mixed (see e.g., Johnston et al. 2008). There are currently few methods that can surpass visual inspection for the identification of spurious detections. It is important to note that surveys are often not sensitive to the most diffuse parts of the debris structure. The two panels of Fig. 15 illustrate this by showing additional structure revealed after an accurate sky background subtraction. Additionally, streams on specific orbits can lead to ‘fans’ of low-density tidal debris (Pearson et al. 2015; Price-Whelan et al. 2016; Yavetz et al. 2021), which are potentially missed. In our survey, we therefore use visual classification to exclude amorphous, diffuse structures in our images that are likely to correspond to more ancient accretion events or map out complex orbits in their host galaxy.

thumbnail Fig. 15.

Evidence of very faint, extended outer loops around the edge-on spirals UGC08717 (top panels) and 2MASXJ12284541-0838329 (bottom panels) detected as a result of the accurate sky background subtraction described in Sect. 2.2. These examples illustrate the effect of the surface brightness cut-off of the data on the interpretation of the overall morphology of the streams.

‘Particle spray’ library and fits. Bespoke N-body simulations can be used to generate streams that provide a qualitative resemblance to observed features (Martínez-Delgado et al. 2008). However, a full exploration of a wide parameter space using N-body simulations is usually too expensive. In our survey, we instead use a ‘particle spray’ technique to rapidly generate an extensive library of mock streams (see e.g., Bonaca et al. 2014; Küpper et al. 2015; Gibbons et al. 2014; Fardal et al. 2015). In the ‘particle spray’ technique, the stars, which escape a progenitor to form leading and trailing stellar arms, are massless test particles. Test particles are released uniformly in time from the progenitor’s Lagrange points, whose locations are calculated based on the progenitor’s initial mass. The progenitor mass can be updated throughout the simulation to mimic its mass-loss. Once the test particles are released from the Lagrange points with a dispersion in the velocity with respect to the progenitor, the positions and velocities of the test particles are evolved in the combined gravitational potential of the host and the progenitor. Gibbons et al. (2014) showed that including the self-gravity of the progenitor is essential for the ‘particle spray’ technique to reproduce the length of streams generated with N-body simulations. The fact that the test particles do not mutually interact with each other makes this technique extremely fast compared to direct N-body simulations. The morphology the resulting mock streams compare remarkably well to those generated with direct N-body simulations (see e.g., Gibbons et al. 2014).

We plan to use the modified Lagrange Cloud stripping variation of the ‘particle spray’ technique presented in Gibbons et al. (2014), an example of which is shown in Fig. 8 (right), to generate a large library of mock tidal disruptions with orbits and masses based on distributions obtained from the cosmological simulations described below. ‘Particle spray’ simulations are carried out in a range of galactic potentials which span the types of galaxies in our survey (i.e. various mass profiles and dark matter halo shapes). Mock debris structures in this library are compared visually to our observations of tidal debris. In particular, we can inject mock debris into real images and test how well our visual inspection works. For example, if a shell-like feature is sufficiently faint, it is possible that only the rim of the shell would be visible above the background (e.g., the bottom, middle panel of Fig. 3). These tests indicate whether such features can be distinguished from a thinner tidal stream on a near-circular orbit.

In addition to constraining the build-up of stellar halos through comparisons to cosmological models, we identify a subset of cases in our survey with prominent streams and use these to fit the potential of individual host galaxies. While many Galactic streams have been used to constrain the potential of the Milky Way (e.g., Koposov et al. 2010; Law & Majewski 2010; Küpper et al. 2015; Bovy et al. 2016; Malhan & Ibata 2019), this has only been attempted for two external galaxies (e.g., Fardal et al. 2013; Amorisco et al. 2015). The ‘particle spray’ technique can also be used to determine which mass distributions best reproduce the properties of observed streams in the Legacy Survey (Gibbons et al. 2014). This technique has already been applied to map the mass distribution in the Milky Way (e.g., Erkal et al. 2019). Constraints on the dark matter potentials of survey galaxies from this technique can be contrasted with those obtained through satellite kinematics and used to constrain relationships between galaxy properties and halo mass, which are fundamental to the galaxy formation models described in the following section.

4.3. Comparison with cosmological simulations

In the previous sections we have demonstrated that all-sky imaging surveys now reach the depth required to assemble statistically representative samples of diffuse tidal features in the Local Universe, and to quantify their frequency of occurrence, mass, morphology, spatial scale and stellar populations. Cosmological simulations are the only practical way to interpret the observed distributions of these properties. In turn, such observations can constrain the ‘subgrid’ models on which the simulations are built.

To exploit the statistical power of our survey, cosmological simulations suitable for comparison must have comparable volume and dynamic range from which to draw representative samples of galaxies. The physical models underlying these simulations are typically calibrated against statistical observables of the galaxy population, such as the stellar mass function, which allows the same selection function to be applied to both observed and simulated galaxies.

For a thorough comparison with our survey data that takes detectability and projection effects into account, it is necessary to derive realistic mock images (e.g., Torrey et al. 2015; Snyder et al. 2015) from a range of different simulations. These are stringent requirements, but nevertheless well-matched to the capabilities of state-of-the-art simulations from multiple groups. Below we describe the potential for comparison between our Stellar Stream Legacy Survey data and two state-of-the-art cosmological simulations, the CoCo and TNG50 simulations, and illustrate the techniques for generating mock images.

Cosmological simulations can be used to quantify intrinsic galaxy-to-galaxy scatter in the properties of streams, for comparison to observed distributions. In addition, we seek to quantify variance in the predicted structure of stellar streams and the disruption times of satellites between different galaxy formation models and simulation codes. For this purpose, we use the cosmological runs of the AGORA project15. AGORA is a comparison project exploring code-to-code (or group-to-group) differences in galaxy simulations (Kim et al. 2014, 2016). In its third phase, the AGORA community is comparing zoom-in cosmological simulations of the same Mvir ∼ 1012 M⊙ halo at redshift zero, with gravitational softening length ϵ = 80 pc, and particle mass resolutions of MDM = 2.8 × 105 M⊙ and M⋆|gas = 5.65 × 104 M⊙ (Roca-Fàbrega et. al in preparation). Using these simulations, we explore the numerical and resolution dependence of predictions for observables relevant to our tidal stream survey, with the ultimate goal of finding new ways in which the survey data can be used to constrain subgrid models.

4.3.1. The Copernicus Complexio models

The Copernicus Complexio (CoCo) cosmological N-body simulations (Hellwing et al. 2016) provides both high mass resolution (1.6 × 105 M⊙ per particle) and a representative analogue of the Local Volume (a high-resolution spherical region of radius ∼25 Mpc embedded in a lower-resolution box of side length 100 Mpc), well suited to comparison with the Stellar Stream Legacy Survey.

We have applied the technique known as particle tagging to CoCo, following the STINGS methodology of Cooper et al. (2017), building on previous applications to the Aquarius and Millennium II simulations (Cooper et al. 2010, 2013). By coupling collisionless simulations to semi-analytic models of galaxy formation, STINGS creates dynamically self-consistent cosmological N-body simulations of galactic stellar halos and their associated tidal streams, at a resolution comparable to current hydrodynamical simulations and without many of the simplifications required by a earlier generation of ‘fast but approximate’ galaxy models (e.g., Bullock & Johnston 2005). In this example we have used the Galform semi-analytic model (Lacey et al. 2016) to predict the evolution of stellar mass, size, and chemical abundances in every CoCo dark matter halo, constrained by comparisons to the large-scale statistics of the cosmological galaxy population (for example, optical and infrared luminosity functions). Individual dark matter particles in CoCo were then “painted” with single-age stellar populations according to a dynamical prescription that specifies the binding energy distribution of stars at the time of their formation (this prescription differs between particle tagging implementations by different groups; the details for STINGS and a comparison to other approaches are given in Cooper et al. 2017). Using the tagged particles, the individual star formation history of each satellite (and hence properties such as stellar mass, luminosity and metallicity) can be studied alongside the full phase-space evolution of its stars, before and after tidal disruption. Particle tagging can be thought of as an extension to the semi-analytic models that allows individual stellar populations to be tracked in phase space, while retaining the computational efficiency and flexibility of the semi-analytic galaxy formation model.

These simulations have sufficient dynamic range both to draw representative samples of galaxies and to study debris features from dwarf galaxies comparable to the classical satellites of the Milky Way. Importantly, CoCo adequately resolves the cores of massive satellite halos in these systems and the fine structure of the debris from their tidal disruption. Their softening scale of 330 pc (comparable to the width of the Milky Way’s Orphan stream) and particle mass resolution of 1.6 × 105 M⊙ are similar to the Medium-Resolution suite of the Apostle Milky Way Zoom simulations (comprising 11 galaxies). CoCo, however, provides a representative cosmological volume containing ∼50–70 Milky Way mass dark matter halos. Assuming a total-mass-to-light ratio of 100 in the region of phase space corresponding to the stellar body of a dwarf like the Fornax dSph, its debris would be resolved with several thousand tagged particles in these simulations.

Figure 16 shows mock images of three approximately L⋆ galaxies from CoCo, embedded in a simulated field constructed using DESI Legacy Imaging Surveys model fits to all sources in a random patch of sky. These examples adopt a fiducial distance of 50 Mpc. From left to right in Fig. 16, we vary the surface brightness limit to simulate, respectively, an SDSS-like observation (≈27 mag arcsec−2), a Stellar Stream Legacy Survey observation (≈29 mag arcsec−2) and an even deeper image comparable to what may be achievable with HSC or Rubin (≈31 mag arcsec−2). In the first two rows, new features appear in the Stellar Stream Legacy Survey-like images that are not visible at the SDSS depth, and these are further enhanced (hence easier to interpret) at HSC-like depth.

thumbnail Fig. 16.

Example r-band mock images of three ∼L⋆ galaxies from CoCo at a fiducial distance of 50 Mpc. From left to right, the images have surface brightness limits of μlim = 27, 29 and 31 magnitudes per square arcsecond, respectively. The dynamic range of each image runs from μlim + 1 (white) to μlim − 9 (black). We have embedded these models in the same fiducial field from the DESI Legacy Imaging surveys. Further details of how we generate these mock images are given in the text. From top to bottom, the stellar masses of the galaxies are log10M⋆/M⊙ = (10.1, 10.6, 10.1), and their virial masses are log10Mvir/M⊙ = (11.8, 12.2, 12.2).

To create these images, we first smooth the projected N-body tagged particle distribution in a box of side length 500 kpc around each galaxy, using a conventional two-dimensional ‘SPH’ cubic spline kernel. The smoothing length of each tagged particle is proportional to the RMS distance to its 16th nearest neighbour16. This smoothed image is further convolved with a simple analytic approximation to the Legacy Surveys PSF. Gaussian noise is added to simulate a sky-subtracted image of the required depth. We use the same definition of limiting surface brightness, μlim, as elsewhere in this paper, using 3σ of the flux distribution in 10″ × 10″ regions. Finally, we superpose the Tractor models of all sources from the DESI Legacy Imaging Surveys catalogue in a random patch of blank sky.

4.3.2. The TNG50 hydrodynamical model

It is reasonable to expect that other cosmological simulations, such as Illustris (Vogelsberger et al. 2014; Genel et al. 2014; Nelson et. al 2014) and IllustrisTNG (Nelson et al. 2018, 2019b; Naiman et al. 2018; Pillepich et al. 2018b; Springel et al. 2018; Marinacci et al. 2018), may make different predictions for the properties of tidal stellar streams and the correlation of such properties with other galactic observables and intrinsic properties. Crucially, these differences arise not only from the different galaxy formation prescriptions used in various models, but also from fundamentally different choices in their underlying dynamical assumptions.

For example, the CoCo simulations assume pure N-body dynamics, with star formation treated semi-analytically17. IllustrisTNG, on the other hand, includes not only gravity but also a self-consistent description of the gas hydrodynamics together with prescriptions for gas cooling and heating, formation and evolution of stellar particles, metal enrichment, galactic scale gas outflows due to star formation, and energy injections from super massive black holes (Weinberger et al. 2017; Pillepich et. al 2018a). Namely, the formation and evolution of massive thin disks of stars and gas are emerging phenomena in these models, and with them so are for example the star formation quenching of low-mass galaxies in dense environments and the assembly of galactic stellar haloes, whether of ex-situ or in-situ origin. The ‘higher order’ processes that are ignored by semi-analytic models but which arise naturally in a hydrodynamical simulation are likely to be important in the detection and interpretation of low surface brightness structures. For example, satellites can interact with galactic disks, which may significantly affect the satellite disruption rate (Errani et al. 2017; Garrison-Kimmel et al. 2017). Those interactions in turn may drive disk heating, warping and buckling, as well as asymmetric features, such as tidal arms.

Large-volume galaxy simulations like Illustris and the TNG100 and TNG300 runs of the IllustrisTNG project18 have already been queried to characterise the imprint of different accretion and merger histories on the stellar components of galaxies, from stellar disks to stellar haloes and intra-cluster light (e.g., Pillepich et al. 2014; Rodriguez-Gomez et al. 2016), including the formation and properties of shell galaxies (Pop et al. 2017, 2018). However, because of the somewhat limited numerical resolution, such studies have been limited to either global properties of stellar haloes or to the fine features in stellar haloes, such as streams and shells, of galaxies a few times more massive than the Milky Way.

Among the currently available cosmological hydrodynamical simulations of Milky Way analogues, the more recent TNG50 simulation (Pillepich et al. 2019; Nelson et. al 2019a) offers a unique combination of statistical sampling and numerical resolution. TNG50 is the third and most ambitious realisation of the IllustrisTNG project. It simulates a 52 Mpc/side volume at numerical resolution similar to zoom-in simulations like Auriga (i.e. about 8 × 104 M⊙/stellar particle). The resulting volume contains ∼200 among Milky Way and Andromeda like galaxies at z = 0, in addition to a few hundred more massive galaxies and thousands of dwarf galaxies, in both isolation or as satellites of more massive systems. The high mass resolution in this volume enables a detailed investigation of the present-day structural properties of galaxies, such as the thickness of gaseous and stellar disks, perturbations in those disks due to interactions with satellites, and low surface brightness features in stellar halos. Figure 17 provides a snapshot of the diversity of stellar halo features of TNG50 simulated galaxies chosen among hundreds with stellar or halo mass similar to those of the Milky Way. In particular, represented is a subset of Milky Way-like galaxies (i.e. with disk-like stellar structures) featuring prominent stellar streams, either linear or curved.

thumbnail Fig. 17.

Example galaxies from the TNG50 simulation. Images are in projected stellar column density above 104 M⊙ kpc−2 and represent a subset of Milky Way-like galaxies featuring prominent stellar streams, either linear or curved. The stellar masses of the galaxies are in the range log10M⋆/M⊙ = 10.3 − 10.8, as indicated. Notice that two of the galaxies are repeated in different projections.

The comparison we anticipate here between the results of the Stellar Stream Legacy Survey and the outcome of the TNG50 simulation can offer great insight on both the galaxy formation model underlying the simulation as well as on the assembly histories of the observed galaxies.

5. Conclusions

In this paper we have introduced the Stellar Stream Legacy Survey, a new survey of low surface brightness features associated with the tidal disruption of dwarf galaxies around isolated massive galaxies in the local Universe. Our survey is designed to constrain the rate of occurrence of those features and to provide robust distributions of their mass, morphology and colour. It is based on an optimised re-reduction of deep g, r and z-band imaging data from the DESI Legacy Imaging Surveys, which are publicly available and cover about 20 000 deg2 (Sect. 2.1). Our parent sample of isolated galaxies comprises ∼3100 systems with K-band absolute magnitude MK < −19.6, at distances of 30 < d < 100 Mpc (Sect. 2.3). The data we obtain allow us to make direct, quantitative comparisons of tidal stream statistics to predictions from state-of-the-art ΛCDM cosmological simulations, and hence to constrain cosmological and sub-grid parameters controlling the assembly rates of massive galaxies and the stellar populations of dwarf satellites.

As a proof of concept of the survey, we have reported here the detection of 24 new stellar streams (Figs. 3 and 4). An initial list of these discoveries and previously known features is presented in Table 1, which we will expand upon in future work. In this paper we have introduced the techniques that are used for the entire survey and have applied them to a subset of galaxies in our sample. We have demonstrated that:

  • Our custom reprocessing of data from the DESI Legacy Imaging Surveys has a r-band limiting depth of ∼28.5 ± 0.4 mag arcsec−2 (3σ of the sky background variation measured in 10 × 10″ pixels). This is comparable to other ultra-deep observations of diffuse structures in the literature and sufficient for the detection of a significant sample of minor merger events associated with our target galaxies. In our proof of concept sample, these features are typically 1 to 3 mag arcsec−2 brighter than our detection limit.

  • The re-reduced Legacy Survey images provide sufficient signal-to-noise to obtain reliable colours for streams using our baseline method for aperture photometry. In our proof-of-concept sample, we find a range of colours 0.5 ≲ g − r ≲ 0.8, somewhat redder than the range of massive satellite galaxies in the SAGA survey. We extend and refine this analysis in subsequent work describing the full sample.

  • Methods to distinguish features arising from dwarf galaxy accretion from those generated by major mergers and perturbations to thin galactic discs can be tested at the resolution of our survey (Sect. 3.2).

  • Our approach offers a much larger sample size, is competitive with the Euclid Wide Survey and only marginally shallower than the 10-year stacked data expected from the Vera C. Rubin Telescope Legacy Survey of Space and Time (Sect. 3.4).

  • The ‘particle spray’ technique can rapidly generate a large library of bespoke N-body realisations of individual streams, with which to test both visual and automated classification methods. We use this technique in combination with convolutional neural networks to quantify stream morphologies in our sample for comparison with cosmological simulations. This approach can constrain galactic accretion histories and satellite orbital distributions. We also use it to identify a large ‘gold standard’ sample of thin streams, in order to constrain the potentials of their host galaxies and provide a new test of the galaxy–halo connection (Sect. 4).

  • The current generation of cosmological volume simulations is well matched to the spatial scale, resolution and sample size of our data. Realistic mock images generated from these models enable quantitative comparisons with the results of our survey. (Sect. 4.3).

The full dataset of the Stellar Stream Legacy Survey is presented in a forthcoming paper. Our analysis of its images acts as a pathfinder for future surveys with higher spatial resolution and greater depth, including those with Euclid, the Vera C. Rubin Observatory and the Nancy Grace Roman Space Telescope.

1

Features arising from tidal disruption are often classified by their visual morphology, for example tidal tails or arms from mergers of the host galaxy with a gas rich, disk-like companion (Byrd & Howard 1992); fans and plumes from ‘dry’ major mergers (Feldmann et al. 2008) and shells around massive early-type galaxies (Schweizer & Seitzer 1992; Cooper et al. 2011). Here we follow the terminology of Duc et al. (2015).

3

The deep images with the robotic telescopes from the Stellar Tidal Stream Survey (Martínez-Delgado 2019) were obtained with a non-standard luminance filter under mediocre seeing conditions; see e.g., Fig. 1 in Martínez-Delgado et al. (2015).

7

Kado-Fong et al. (2018) proposed a method to automatically detect which galaxies have streams (or any kind of non-smooth faint structure) based on successive convolution of the images with the same kernel. Although this approach is robust, their suggested kernel and parameters are specific to their study of a more distant galaxy sample at 0.05 < z < 0.45 with Subaru Hyper Suprime-Cam imaging. We have tried to apply such a solution for the relatively massive, closer galaxies in our sample, without success. We find that this approach is not able to remove the host galaxy as clearly as they showed in their images, mainly because our target galaxies have a much larger sizes (in pixels) than the small kernel suggested by Kado-Fong et al. (2018). Further investigation into a suitable form of the kernel are needed to optimise this approach for samples of nearby galaxies such as ours.

8

Because the focus here is on the large and extended streams, NoiseChisel was given a Gaussian kernel of FWHM = 3 pixels (default 2 pixels), truncated at 4 times the FWHM. As well as --tilesize=40 and --detgrowmaxholesize=10 000, the default values were used for all other parameters.

9

To better separate the extended emission of background galaxies or foreground stars as clumps (to be masked), Segment was run with a Gaussian kernel of FWHM = 5 (default 1.5 pixels), the default values were used for all other parameters.

12

This is the quantity referred to as the ’Upper Limit’ in Gnuastro, see https://www.gnu.org/software/gnuastro/manual/html_node/Upper-limit-magnitude-of-each-detection.html

13

The SAGA colour-magnitude selection excludes galaxies with g − r ≳ 0.7 at r ≳ 17.

14

Although they are not our primary target, features associated with low mass ratio mergers may also provide useful constraints on galaxy assembly. Our survey is likely to discover a significant number of such features.

16

Each stellar population associated with a single N-body particle is smoothed by the same amount. This is for cosmetic purposes and the optimal choice of smoothing length is somewhat arbitrary; (see e.g., Lowing et al. 2011) for a more rigorous approach to smoothing in mock observations of a sparsely sampled stellar phase space.

17

The gravitational dynamics of baryons are included in the semi-analytic cooling, star formation and feedback models, albeit approximately.

18

http://www.illustris-project.org/ and http://www.tng-project.org/

Acknowledgments

We thank Eva K. Grebel, Facundo A. Gómez, Aaron Romanowski and Emilio Galvez for useful comments. We thank Mireia Montes for her insight and expertise regarding the capabilities of Vera C. Rubin Telescope Survey. DMD acknowledges financial support from the Talentia Senior Program (through the incentive ASE-136) from Secretaría General de Universidades, Investigación y Tecnología, de la Junta de Andalucía. DMD acknowledges funding from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709) and AEI project (PDI2020-114581GB-C21/AEI/10.13039/501100011033). APC is supported by the Taiwan Ministry of Education (MoE) Yushan Fellowship and Taiwan National Science and Technology Council (NSTC) grant 109-2112-M-007-011-MY3. This work made use of HPC facilities at the National Tsing Hua University Center for Informatics and Computation in Astronomy (CICA), supported by grants from MoE and NSTC. This work has also been supported by the Spanish Ministry of Economy and Competitiveness (MINECO) under grant AYA2014-56795-P. This work was supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr). AP and DMD acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subprojects A1 and A2). JR acknowledges financial support from the grants AYA2015-65973-C3-1-R and RTI2018-096228- BC31 (MINECO/FEDER, UE) and support from the State Research Agency (AEI-MCINN) of the Spanish Ministry of Science and Innovation under the grant “The structure and evolution of galaxies and their central regions” with reference PID2019-105602GB-I00/10.13039/501100011033. MA acknowledges support from the Spanish Ministry of Economy and Competi-tiveness (MINECO) under grant number AYA2016-76219-P and EU’s Horizon 2020 research and innovation programme under Marie Sklodowska-Curie grant agreement No 721463 to the SUNDIAL ITN. W.J.P. has been supported by the Polish National Science Center project UMO-2020/37/B/ST9/00466. CSF acknowledges support by the European Research Council (ERC) through Advanced Investigator grant DMIDAS (GA 786910) and by the Science and Technology Facilities Council (STFC) consolidated grant ST/P000541/1. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (https://www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. A.B. was supported by an appointment to the NASA Postdoctoral Program at the NASA Ames Research Center, administered by Universities Space Research Association under contract with NASA. This project uses data from observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. We acknowledge support from the Spanish Ministry for Science, Innovation and Universities and FEDER funds through grant AYA2016-81065-C2-2. We also used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovacão, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Enérgeticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. Support for this work was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51466.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This work was partly done using GNU Astronomy Utilities (Gnuastro, ascl.net/1801.009) version 0.14. Work on Gnuastro has been funded by the Japanese MEXT scholarship and its Grant-in-Aid for Scientific Research (21244012, 24253003), the European Research Council (ERC) advanced grant 339659-MUSICOS, and from the Spanish Ministry of Economy and Competitiveness (MINECO) under grant number AYA2016-76219-P.

References

  1. Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543 [Google Scholar]
  2. Abbott, T. M. C., Abdalla, F. B., Allam, S., et al. 2018, ApJS, 239, 18 [Google Scholar]
  3. Ackermann, S., Schawinski, K., Zhang, C., Weigel, A. K., & Turp, M. D. 2018, MNRAS, 479, 415 [NASA ADS] [CrossRef] [Google Scholar]
  4. Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8 [NASA ADS] [Google Scholar]
  5. Akhlaghi, M. 2019, ASPC, 521, 299A [Google Scholar]
  6. Akhlaghi, M. 2020, IAU Symposium 355, in press, [arXiv:1909.11230] [Google Scholar]
  7. Akhlaghi, M., & Ichikawa, T. 2015, ApJS, 220, 1 [Google Scholar]
  8. Amorisco, N. C. 2015, MNRAS, 450, 575 [NASA ADS] [CrossRef] [Google Scholar]
  9. Amorisco, N. C., Martinez-Delgado, D., & Schedler, J. 2015, ArXiv e-prints [arXiv:1504.03697] [Google Scholar]
  10. Arp, H. 1966, ApJS, 14, 1 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bauer, J. S., Widrow, L. M., & Erkal, D. 2018, MNRAS, 476, 198 [NASA ADS] [CrossRef] [Google Scholar]
  12. Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]
  13. Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJ, 642, L137 [Google Scholar]
  14. Bernard, E. J., Ferguson, A. M. N., Schlafly, E. F., et al. 2016, MNRAS, 463, 1759 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bickley, R. W., Bottrell, C., Hani, M. H., et al. 2021, MNRAS, 504, 372 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bonaca, A., Geha, M., Küpper, A. H. W., et al. 2014, ApJ, 795, 94 [NASA ADS] [CrossRef] [Google Scholar]
  17. Borlaff, A., Trujillo, I., Román, J., et al. 2019, A&A, 621, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bovy, J., Bahmanyar, A., Fritz, T. K., & Kallivayalil, N. 2016, ApJ, 833, 31 [NASA ADS] [CrossRef] [Google Scholar]
  20. Brasseur, C. M., Martin, N. F., Macciò, A. V., et al. 2011, ApJ, 743, 179 [NASA ADS] [CrossRef] [Google Scholar]
  21. Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931 [Google Scholar]
  22. Byrd, G. G., & Howard, S. 1992, AJ, 103, 1089 [NASA ADS] [CrossRef] [Google Scholar]
  23. Carlin, J. L., Beaton, R. L., Martínez-Delgado, D., et al. 2016, Tidal Streams in the Local Group and Beyond, 219 [NASA ADS] [CrossRef] [Google Scholar]
  24. Carlin, J. L., Garling, C. T., Peter, A. H. G., et al. 2019, ApJ, 886, 109 [NASA ADS] [CrossRef] [Google Scholar]
  25. Chequers, M. H., Widrow, L. M., & Darling, K. 2018, MNRAS, 480, 4244 [CrossRef] [Google Scholar]
  26. Cooper, A. P., Cole, S., Frenk, C. S., et al. 2010, MNRAS, 406, 744 [Google Scholar]
  27. Cooper, A. P., Martínez-Delgado, D., Helly, J., et al. 2011, ApJ, 743, L21 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cooper, A. P., D’Souza, R., Kauffmann, G., et al. 2013, MNRAS, 434, 3348 [Google Scholar]
  29. Cooper, A. P., Cole, S., Frenk, C. S., et al. 2017, MNRAS, 469, 1691 [NASA ADS] [CrossRef] [Google Scholar]
  30. Crnojević, D. 2017, Galaxies, 5, 40 [CrossRef] [Google Scholar]
  31. Crnojević, D., Sand, D. J., Spekkens, K., et al. 2016, ApJ, 823, 19 [Google Scholar]
  32. Davies, A., Serjeant, S., & Bromley, J. M. 2019, MNRAS, 487, 5263 [NASA ADS] [CrossRef] [Google Scholar]
  33. Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44 [Google Scholar]
  34. Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, MNRAS, 490, 3426 [NASA ADS] [CrossRef] [Google Scholar]
  35. De Lucia, G., & Helmi, A. 2008, MNRAS, 391, 14 [NASA ADS] [CrossRef] [Google Scholar]
  36. DESI Collaboration (Aghamousa, A., et al.) 2016a, ArXiv e-prints [arXiv:1611.00036] [Google Scholar]
  37. DESI Collaboration (Aghamousa, A., et al.) 2016b, ArXiv e-prints [arXiv:1611.00037] [Google Scholar]
  38. Dey, A., Rabinowitz, D., Karcher, A., et al. 2016, Proc. SPIE, 9908, 99082C [NASA ADS] [Google Scholar]
  39. Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168 [Google Scholar]
  40. Domínguez Sánchez, H., Huertas-Company, M., Bernardi, M., et al. 2018, MNRAS, 476, 3661 [CrossRef] [Google Scholar]
  41. Duc, P.-A., Cuillandre, J.-C., Karabal, E., et al. 2015, MNRAS, 446, 120 [Google Scholar]
  42. Erkal, D., Sanders, J. L., & Belokurov, V. 2016, MNRAS, 461, 1590 [NASA ADS] [CrossRef] [Google Scholar]
  43. Erkal, D., Koposov, S. E., & Belokurov, V. 2017, MNRAS, 470, 60 [NASA ADS] [CrossRef] [Google Scholar]
  44. Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, MNRAS, 487, 2685 [Google Scholar]
  45. Errani, R., Penarrubia, J., Laporte, C. F. P., et al. 2017, MNRAS, 465, 59 [Google Scholar]
  46. Euclid Collaboration (Borlaff, A. S., et al.) 2022a, A&A, 657, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Euclid Collaboration (Scaramella, R., et al.) 2022b, A&A, 662, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Fardal, M. A., Weinberg, M. D., Babul, A., et al. 2013, MNRAS, 434, 2779 [Google Scholar]
  49. Fardal, M. A., Huang, S., Weinberg, M. D., et al. 2015, MNRAS, 452, 301 [NASA ADS] [CrossRef] [Google Scholar]
  50. Feldmann, R., Mayer, L., & Carollo, C. M. 2008, ApJ, 684, 1062 [NASA ADS] [CrossRef] [Google Scholar]
  51. Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, AJ, 150, 150 [Google Scholar]
  52. Fliri, J., & Trujillo, I. 2016, MNRAS, 456, 1359 [Google Scholar]
  53. Font, A. S., Benson, A. J., Bower, R. G., et al. 2011, MNRAS, 417, 1260 [NASA ADS] [CrossRef] [Google Scholar]
  54. Garrison-Kimmel, S., Wetzel, A., Bullock, J. S., et al. 2017, MNRAS, 471, 1709 [NASA ADS] [CrossRef] [Google Scholar]
  55. Geha, M., Wechsler, R. H., Mao, Y.-Y., et al. 2017, ApJ, 847, 4 [NASA ADS] [CrossRef] [Google Scholar]
  56. Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175 [Google Scholar]
  57. Gibbons, S. L. J., Belokurov, V., & Evans, N. W. 2014, MNRAS, 445, 3788 [NASA ADS] [CrossRef] [Google Scholar]
  58. Gilhuly, C., Hendel, D., Merritt, A., et al. 2020, ApJ, 897, 108 [NASA ADS] [CrossRef] [Google Scholar]
  59. Gordon, Y. A., Pimbblet, K. A., Kaviraj, S., et al. 2019, ApJ, 878, 88 [NASA ADS] [CrossRef] [Google Scholar]
  60. Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [Google Scholar]
  61. Gómez, F. A., White, S. D. M., Marinacci, F., et al. 2016, MNRAS, 456, 2779 [Google Scholar]
  62. Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179 [NASA ADS] [Google Scholar]
  63. Grillmair, C. J. 2006, ApJ, 651, L29 [Google Scholar]
  64. Hellwing, W. A., Frenk, C. S., Cautun, M., et al. 2016, MNRAS, 457, 3492 [NASA ADS] [CrossRef] [Google Scholar]
  65. Hendel, D., & Johnston, K. V. 2015, MNRAS, 454, 2472 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hendel, D., Johnston, K. V., Patra, R. K., et al. 2019, MNRAS, 486, 3604 [NASA ADS] [CrossRef] [Google Scholar]
  67. Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
  68. Hood, C. E., Kannappan, S. J., Stark, D. V., et al. 2018, ApJ, 857, 144 [NASA ADS] [CrossRef] [Google Scholar]
  69. Huchra, J. P., Macri, L. M., Masters, K. L., et al. 2012, ApJS, 199, 26 [Google Scholar]
  70. Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194 [Google Scholar]
  71. Ibata, R., Martin, N. F., Irwin, M., et al. 2007, ApJ, 671, 1591 [CrossRef] [Google Scholar]
  72. Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6 [Google Scholar]
  73. Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111 [Google Scholar]
  74. Johnston, K. V., Sackett, P. D., & Bullock, J. S. 2001, ApJ, 557, 137 [NASA ADS] [CrossRef] [Google Scholar]
  75. Johnston, K. V., Bullock, J. S., Sharma, S., et al. 2008, ApJ, 689, 936 [Google Scholar]
  76. Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [Google Scholar]
  77. Kado-Fong, E., Greene, J. E., Hendel, D., et al. 2018, ApJ, 866, 103 [NASA ADS] [CrossRef] [Google Scholar]
  78. Karachentsev, I. D. 2005, AJ, 129, 178 [NASA ADS] [CrossRef] [Google Scholar]
  79. Karachentsev, I. D., & Makarov, D. A. 1996, AJ, 111, 794 [NASA ADS] [CrossRef] [Google Scholar]
  80. Karachentsev, I. D., Kashibadze, O. G., Makarov, D. I., & Tully, R. B. 2009, MNRAS, 393, 1265 [NASA ADS] [CrossRef] [Google Scholar]
  81. Kazantzidis, S., Zentner, A. R., Kravtsov, A. V., et al. 2009, ApJ, 700, 1896 [NASA ADS] [CrossRef] [Google Scholar]
  82. Kim, J., Abel, T., Agertz, O., et al. 2014, ApJS, 210, 14 [Google Scholar]
  83. Kim, J., Agertz, O., Teyssier, R., et al. 2016, ApJ, 833, 202 [NASA ADS] [CrossRef] [Google Scholar]
  84. Kirby, E., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102 [NASA ADS] [CrossRef] [Google Scholar]
  85. Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260 [Google Scholar]
  86. Küpper, A. H. W., Balbinot, E., Bonaca, A., et al. 2015, ApJ, 803, 80 [Google Scholar]
  87. Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2016, MNRAS, 462, 3854 [Google Scholar]
  88. Laine, S., Knapen, J. H., Muñoz-Mateos, J.-C., et al. 2014, MNRAS, 444, 3015 [NASA ADS] [CrossRef] [Google Scholar]
  89. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., et al. 2018, MNRAS, 481, 286 [NASA ADS] [CrossRef] [Google Scholar]
  90. Laporte, C. F. P., Johnston, K. V., & Tzanidakis, A. 2019, MNRAS, 483, 1427 [Google Scholar]
  91. Laporte, C. F. P., Belokurov, V., Koposov, S. E., et al. 2020, MNRAS, 492, L61 [NASA ADS] [CrossRef] [Google Scholar]
  92. Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229 [Google Scholar]
  93. Lecun, Y., Bengio, Y., & Hinton, G. 2015, Nature, 521, 436 [CrossRef] [PubMed] [Google Scholar]
  94. Lotz, J. M., Jonsson, P., Cox, T. J., & Primack, J. R. 2008, MNRAS, 391, 1137 [Google Scholar]
  95. Lowing, B., Jenkins, A., Eke, V., et al. 2011, MNRAS, 416, 2697 [NASA ADS] [CrossRef] [Google Scholar]
  96. Ludwig, J., Pasquali, A., Grebel, E. K., et al. 2012, AJ, 144, 190 [NASA ADS] [CrossRef] [Google Scholar]
  97. Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 3 [Google Scholar]
  98. Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., et al. 2003, ApJ, 599, 1082 [Google Scholar]
  99. Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Malin, D. F., & Carter, D. 1980, Nature, 285, 643 [CrossRef] [Google Scholar]
  101. Malin, D., & Hadley, B. 1997, PASA, 14, 52 [NASA ADS] [CrossRef] [Google Scholar]
  102. Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 [NASA ADS] [Google Scholar]
  103. Martínez-Delgado, D., 2019, Highlights on Spanish Astrophysics X, 146 [Google Scholar]
  104. Martínez-Delgado, D., Aparicio, A., Gómez-Flechoso, M. Á., et al. 2001, ApJ, 549, L199 [CrossRef] [Google Scholar]
  105. Martínez-Delgado, D., Peñarrubia, J., Gabany, R. J., et al. 2008, ApJ, 689, 184 [CrossRef] [Google Scholar]
  106. Martínez-Delgado, D., Gabany, R. J., Crawford, K., et al. 2010, AJ, 140, 962 [Google Scholar]
  107. Martínez-Delgado, D., D’Onghia, E., Chonis, T. S., et al. 2015, AJ, 150, 116 [CrossRef] [Google Scholar]
  108. Mateo, M., Olszewski, E. W., & Morrison, H. L. 1998, ApJ, 508, L55 [NASA ADS] [CrossRef] [Google Scholar]
  109. McConnachie, A. W. 2012, AJ, 144, 4 [Google Scholar]
  110. McConnachie, A. W., Irwin, M. J., Ibata, R. A., et al. 2009, Nature, 461, 66 [Google Scholar]
  111. McConnachie, A. W., Ibata, R., Martin, N., et al. 2018, ApJ, 868, 55 [Google Scholar]
  112. Mao, Y.-Y., Geha, M., Wechsler, R. H., et al. 2021, ApJ, 907, 85 [NASA ADS] [CrossRef] [Google Scholar]
  113. Meisner, A. M., Lang, D., Schlafly, E. F., et al. 2019, PASP, 131, 124504 [CrossRef] [Google Scholar]
  114. Mihos, J. C. 2019, IAU Symposium, 355, [arXiv:1909.09456] [Google Scholar]
  115. Mihos, J. C., Harding, P., Feldmeier, J. J., et al. 2017, ApJ, 834, 16 [Google Scholar]
  116. Miville-Deschênes, M.-A., & Lagache, G. 2005, ApJS, 157, 302 [Google Scholar]
  117. Morales, G., Martínez-Delgado, D., Grebel, E. K., et al. 2018, A&A, 614, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]
  119. Mouhcine, M., & Ibata, R. 2009, MNRAS, 399, 737 [CrossRef] [Google Scholar]
  120. Mouhcine, M., Ibata, R., & Rejkuba, M. 2010, ApJ, 714, L12 [NASA ADS] [CrossRef] [Google Scholar]
  121. Mutch, S. J., Croton, D. J., & Poole, G. B. 2011, ApJ, 736, 84 [NASA ADS] [CrossRef] [Google Scholar]
  122. Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 [Google Scholar]
  123. Nelson, D., Pillepich, A., Genel, S., et al. 2014, A&C, 13, 12 [NASA ADS] [Google Scholar]
  124. Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624 [Google Scholar]
  125. Nelson, D., Springel, V., Pillepich, A., et al. 2019a, ComAC, 6, 2 [NASA ADS] [Google Scholar]
  126. Nelson, D., Pillepich, A., Springel, V., et al. 2019b, MNRAS, 490, 3234 [Google Scholar]
  127. Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245 [Google Scholar]
  128. Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713 [NASA ADS] [CrossRef] [Google Scholar]
  129. Pasquet-Itam, J., & Pasquet, J. 2018, A&A, 611, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. Pearson, S., Küpper, A. H. W., Johnston, K. V., et al. 2015, ApJ, 799, 28 [NASA ADS] [CrossRef] [Google Scholar]
  131. Pearson, W. J., Wang, L., Alpaslan, M., et al. 2019a, A&A, 631, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  132. Pearson, S., Starkenburg, T. K., Johnston, K. V., et al. 2019b, ApJ, 883, 87 [NASA ADS] [CrossRef] [Google Scholar]
  133. Pearson, W. J., Wang, L., Trayford, J. W., Petrillo, C. E., & van der Tak, F. F. S. 2019c, A&A, 626, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  134. Pearson, W. J., Wang, L., Alpaslan, M., et al. 2019d, A&A, 631, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  135. Peñarrubia, J., Martínez-Delgado, D., Rix, H. W., et al. 2005, ApJ, 626, 128 [Google Scholar]
  136. Petrillo, C. E., Tortora, C., Chatterjee, S., et al. 2017, MNRAS, 472, 1129 [Google Scholar]
  137. Pillepich, A., Vogelsberger, M., Deason, A., et al. 2014, MNRAS, 444, 237 [Google Scholar]
  138. Pillepich, A., Madau, P., & Mayer, L. 2015, ApJ, 799, 184 [Google Scholar]
  139. Pillepich, A., Springel, V., Nelson, D., et al. 2018a, MNRAS, 473, 4077 [Google Scholar]
  140. Pillepich, A., Nelson, D., Hernquist, L., et al. 2018b, MNRAS, 475, 648 [Google Scholar]
  141. Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196 [Google Scholar]
  142. Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  143. Pop, A., Pillepich, A., Amorisco, N. C., et al. 2017, Galaxies, 5, 34 [NASA ADS] [CrossRef] [Google Scholar]
  144. Pop, A., Pillepich, A., Amorisco, N. C., et al. 2018, MNRAS, 480, 1715 [NASA ADS] [CrossRef] [Google Scholar]
  145. Price-Whelan, A. M., Johnston, K. V., Sheffield, A. A., Laporte, C. F. P., & Sesar, B. 2015, MNRAS, 452, 676 [Google Scholar]
  146. Price-Whelan, A. M., Johnston, K. V., Valluri, M., et al. 2016, MNRAS, 455, 1079 [NASA ADS] [CrossRef] [Google Scholar]
  147. Ragusa, R., Spavone, M., Iodice, E., et al. 2021, A&A, 651, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  148. Rich, R. M., Mosenkov, A., Lee-Saunders, H., et al. 2019, MNRAS, 490, 1539 [Google Scholar]
  149. Rocha-Pinto, H. J., Majewski, S. R., Skrutskie, M. F., Crane, J. D., & Patterson, R. J. 2004, ApJ, 615, 732 [Google Scholar]
  150. Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371 [Google Scholar]
  151. Román, J., Trujillo, I., & Montes, M. 2020, A&A, 644, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  152. Sanders, J. L., & Binney, J. 2013, MNRAS, 433, 1826 [CrossRef] [Google Scholar]
  153. Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [Google Scholar]
  154. Schweizer, F., & Seitzer, P. 1992, AJ, 104, 1039 [Google Scholar]
  155. Scott, C., & Kaviraj, S. 2014, MNRAS, 437, 2137 [NASA ADS] [CrossRef] [Google Scholar]
  156. Shipp, N., Drlica-Wagner, A., Balbinot, E., et al. 2018, ApJ, 862, 114 [Google Scholar]
  157. Simonyan, K., & Zisserman, A. 2014, ArXiv e-prints [arXiv:1409.1556] [Google Scholar]
  158. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
  159. Slater, C. T., Bell, E. F., Schlafly, E. F., et al. 2014, ApJ, 791, 9 [Google Scholar]
  160. Smee, S. A., Gunn, J. E., Uomoto, A., et al. 2013, AJ, 146, 32 [Google Scholar]
  161. Snyder, G. F., Torrey, P., Lotz, J. M., et al. 2015, MNRAS, 454, 1886 [NASA ADS] [CrossRef] [Google Scholar]
  162. Springel, V., & Hernquist, L. 2005, ApJ, 622, L9 [Google Scholar]
  163. Sola, E., Duc, P.-A., Richards, F., et al. 2022, A&A, 662, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  164. Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676 [Google Scholar]
  165. Tal, T., van Dokkum, P. G., Nelan, J., et al. 2009, AJ, 138, 1417 [NASA ADS] [CrossRef] [Google Scholar]
  166. Toomre, A., & Toomre, J. 1972, ApJ, 178, 623 [Google Scholar]
  167. Torrey, P., Snyder, G. F., Vogelsberger, M., et al. 2015, MNRAS, 447, 2753 [NASA ADS] [CrossRef] [Google Scholar]
  168. Trujillo, I., & Fliri, J. 2016, ApJ, 823, 123 [Google Scholar]
  169. Vega-Ferrero, J., Domínguez Sánchez, H., Bernardi, M., et al. 2021, MNRAS, 506, 1927 [NASA ADS] [CrossRef] [Google Scholar]
  170. Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [Google Scholar]
  171. Walmsley, M., Ferguson, A. M. N., Mann, R. G., et al. 2019, MNRAS, 483, 2968 [NASA ADS] [CrossRef] [Google Scholar]
  172. Wang, M., & Deng, W. 2018, Neurocomputing, 312, 135 [CrossRef] [Google Scholar]
  173. Wang, L., Pearson, W. J., & Rodriguez-Gomez, V. 2020, A&A, 644, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  174. Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291 [Google Scholar]
  175. Williams, G. G., Olszewski, E., Lesser, M. P., et al. 2004, Proc. SPIE, 5492, 787 [NASA ADS] [CrossRef] [Google Scholar]
  176. Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
  177. Yavetz, T. D., Johnston, K. V., Pearson, S., et al. 2021, MNRAS, 501, 1791 [Google Scholar]
  178. York, D. G., Adelman, J., Anderson, J. E., et al. 2000, AJ, 120, 1579 [Google Scholar]
  179. Zou, H., Zhou, X., Fan, X., et al. 2019, ApJS, 245, 4 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

List of galaxies with stellar streams discovered during the proof-of-concept phase of our new survey.

In the text
Table 2.

Photometry of stellar streams discovered in our proof-of-concept analysis for a sub-sample of Stellar Stream Legacy Survey targets, listed in Table 1.

In the text
Table 3.

Comparison between the Legacy Survey and future surveys.

In the text

All Figures

thumbnail Fig. 1.

The Stellar Stream Legacy Survey sky coverage. Top: The DESI Legacy Imaging Surveys footprint (blue), the source of optical data for our survey. The 2MASS Redshift Survey catalogue (Huchra et al. 2012) beyond the local group and below VLG < 7000 km s−1 is plotted over the whole sky with a marker size scaled to the K-band magnitude. Those galaxies matching our selection criteria (2000 < VLG < 7000 km s−1, K-band absolute magnitude MK < −19.6) are displayed in bright red colour. The Stellar Stream Legacy Survey sub-samples comprise approximately 3000 of these galaxies in bright red within the blue area. Low galactic latitude regions of the survey affected by the higher stellar density of the Galaxy are perceived with this Gaia overlay. The galactic referential is plotted in green. Bottom: E(B − V) distribution across the DESI Legacy Imaging Surveys footprint (Planck Collaboration XI 2014). Areas of greater dust density (dark blue) are also contaminated with the Galactic light diffused by the dust grains (galactic cirrus) causing a slightly higher background and occasional overlapping structures at scales similar to the faint streams motivating this study. The same 2MRS galaxies are represented on this second map with a different colour scheme.
In the text
thumbnail Fig. 2.

Example of the deep mosaics we construct from the individual DESI Legacy Imaging Surveys exposures overlapping our target galaxies. From left to right, panels show g, r, and z bands. The white rectangle shows the region we process with our custom pipeline around a given target (here NGC 2460), coloured rectangles show each individual CCD frame in the stack. The overlapping tiling pattern is visible. Since this galaxy is in the northern footprint of the DESI Legacy Imaging Surveys, the g and r bands are from BASS, where the CCD pixel scale is larger than in the z band, which is from MzLS. The right-most panel shows the final mosaic.
In the text
thumbnail Fig. 3.

Stacked g and r image cutouts of 12 stellar streams discovered during the proof-of-concept phase of the Stellar Stream Legacy Survey introduced in this paper. For illustrative purposes, colour insets of the central region of the host galaxies have been added to the negative version of the images.
In the text
thumbnail Fig. 4.

As Fig. 3, but showing an additional 12 new stellar streams from our Stellar Stream Legacy Survey.
In the text
thumbnail Fig. 5.

Detecting and measuring streams. (a) The input image: NGC 4981 in z band. The red aperture has a radius of 7 arcsec and lies on the thinnest and most diffuse end of the stream. (b) All pixels not belonging to the galaxy (detected by NoiseChisel) or its stream are masked, including the pixels of objects that are in the background or foreground of the stream (‘clumps’ in Segment). (c) All masked pixels interpolated by the median of ten neighbors. (d) The flux in the red aperture (dotted line), compared to a distribution of 10 000 similarly sized apertures, randomly placed over undetected pixels of the image (4.32σ). (e) The input dataset warped to a pixel grid where each pixel is 8 × 8 larger than the input’s pixels, to demonstrate the detection success on the input resolution.
In the text
thumbnail Fig. 6.

Example of our aperture photometry approach applied to the stellar stream detected around NGC 3451. Top panel: input image in g-band, with subtracted sky, showing the stream around it. Middle panel: the detection contour of the galaxy and its stream are highlighted (black line) on a pixel grid 8 × 8 courser than the input image to highlight the low surface brightness signal. Bottom panel: the apertures used to measured the surface brightness, magnitude and colours are placed along the stream (red circles). The foreground stars and background sources are masked (white areas).
In the text
thumbnail Fig. 7.

Histograms showing the distribution of the depth of the image cutouts and of the photometry measurements of the stellar streams around the galaxies listed in Table 2. From top to bottom: Surface brightness limit for the r band, surface brightness for the r band, and colours (g − r) and (r − z).
In the text
thumbnail Fig. 8.

Left panels: examples of likely major merger remnants in our proof-of-concept sample (see Table 1): ESO474-26, NGC 655, NGC 3303 and NGC 7288. These were identified as massive merger remnants based on the similarity to the outcome of idealised simulations: see right panels. Right panels: simulation snapshots of galaxy mergers evolved for 2 Gyr (see Sect. 4.1). Each corresponds to different initial conditions: rows show different orbital eccentricities given by their circularity L/Lc as defined by Hendel & Johnston (2015); columns show different progenitor stellar masses. The 5 × 109 M⊙ progenitor corresponds to a major merger with a mass ratio ∼1:3; the other two models are for minor mergers. Progenitors were initialised on the size–luminosity relation of Local Group satellites (McConnachie 2012). The colour bar shows the surface brightness of the tidal debris in mag arcsec−2.
In the text
thumbnail Fig. 9.

Gallery of simulated galactic feathers shown at different stages of interaction (different columns) and different viewing angles (rows). Top: edge-on view; middle: face-on view; and bottom: inclined view. Many of these structures show morphologies reminiscent of tidal streams from disrupting satellites although in this case this material is purely in-situ from the disc. A colour characterisation or metallicity measurement could potentially disentangle the two possibilities.
In the text
thumbnail Fig. 10.

Image cutouts of galactic feather candidates with diverse morphologies identified in the proof-of-concept phase of our survey (see Table 1): NGC 4981, NGC 5492, NGC 5635 and ESO245-010.
In the text
thumbnail Fig. 11.

g − r vs. r − z colours for the tidal streams presented in Table 2 (we exclude tidal features with colour errors > 0.25 mag for clarity) and the Galactic cirri regions characterised by Román et al. (2020). The black dashed line, (r − z) = 0.20 × (g − r)+0.20, limits the regions where the colours are compatible/incompatible (below/above of the line) with Galactic dust.
In the text
thumbnail Fig. 12.

Comparison of the surface brightness limit (3σ, 10 × 10 arcsec2) as a function of the observed area for a selection of deep optical and NIR surveys. Depths are measured in the r-band (noted when not applicable). The Stellar Stream Legacy Survey presented in this paper is indicated with a blue star. The other surveys comprise: Stellar Tidal Stream Survey (Martínez-Delgado 2019, STSS), Euclid VIS (Euclid Collaboration 2022a), SDSS (York et al. 2000), IAC Stripe 82 (Fliri & Trujillo 2016), the MATLAS deep imaging Survey (g-band magnitudes Duc et al. 2015), Hyper Suprime-Cam DR2 (Aihara et al. 2018), NGC 4565 Dragonfly observation (Gilhuly et al. 2020), HST WFC3 ABYSS HUDF (NIR bands Borlaff et al. 2019), XDF (Illingworth et al. 2013), UGC00180 10.4m GTC exploration (Trujillo & Fliri 2016), the Burrell Schmidt Deep Virgo Survey (Mihos et al. 2017), the CFIS/UNIONS survey (Sola et al. 2022), the VEGAS/VLT survey (Ragusa et al. 2021), and the Rubin Telescope Legacy Survey of Space and Time (10 yr full survey integration, Ivezić et al. 2019).
In the text
thumbnail Fig. 13.

Comparison of the image cutouts obtained by applying the methodology described in this paper to data from the DESI Legacy Imaging Surveys (top panels) versus long exposure times (6–10 h) images obtained with amateur robotic telescopes from the Stellar Tidal Stream Survey (Martínez-Delgado 2019, bottom panels): our new Stellar Stream Legacy Survey is well positioned to uncover hundreds of low surface brightness features in the local Universe.
In the text
thumbnail Fig. 14.

Bottom panels: external perspective snapshots for the three morphological types classified by Johnston et al. (2008): “Great circle” streams (left) originated from satellites accreted 6–10 Gyr ago on nearly circular orbits; “shells” or “umbrellas” (middle) arise from accretion events within the past ∼8 Gyr with radial orbits; and “mixed”-type tidal remnants (right) formed from ancient (more than 10 Gyr ago) accretion events that have fully mix along their orbits. Top panels: image cutouts obtained from the DESI Legacy Imaging Surveys data for each morphological case of tidal stream: a great circle stream in NGC 577 (left), shell-like features around NGC 681 (middle); and a “mixed” remnant around NGC 5866 (right). Figure adapted from Carlin et al. (2016).
In the text
thumbnail Fig. 15.

Evidence of very faint, extended outer loops around the edge-on spirals UGC08717 (top panels) and 2MASXJ12284541-0838329 (bottom panels) detected as a result of the accurate sky background subtraction described in Sect. 2.2. These examples illustrate the effect of the surface brightness cut-off of the data on the interpretation of the overall morphology of the streams.
In the text
thumbnail Fig. 16.

Example r-band mock images of three ∼L⋆ galaxies from CoCo at a fiducial distance of 50 Mpc. From left to right, the images have surface brightness limits of μlim = 27, 29 and 31 magnitudes per square arcsecond, respectively. The dynamic range of each image runs from μlim + 1 (white) to μlim − 9 (black). We have embedded these models in the same fiducial field from the DESI Legacy Imaging surveys. Further details of how we generate these mock images are given in the text. From top to bottom, the stellar masses of the galaxies are log10M⋆/M⊙ = (10.1, 10.6, 10.1), and their virial masses are log10Mvir/M⊙ = (11.8, 12.2, 12.2).
In the text
thumbnail Fig. 17.

Example galaxies from the TNG50 simulation. Images are in projected stellar column density above 104 M⊙ kpc−2 and represent a subset of Milky Way-like galaxies featuring prominent stellar streams, either linear or curved. The stellar masses of the galaxies are in the range log10M⋆/M⊙ = 10.3 − 10.8, as indicated. Notice that two of the galaxies are repeated in different projections.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.