Studying the Physical Properties of Tidal Features I. Extracting Morphological Substructure in CANDELS Observations and VELA Simulations
Kameswara Bharadwaj Mantha Daniel H. McIntosh Cody P. Ciaschi Rubyet Evan Henry C. Ferguson Logan B. Fries Yicheng Guo Anton M. Koekemoer Luther D. Landry Elizabeth J. McGrath Raymond C. Simons Gregory F. Snyder Scott E. Thompson Eric F. Bell Daniel Ceverino Nimish P. Hathi Camilla Pacifici Joel R. Primack Marc Rafelski Vicente Rodriguez-Gomez

	Non-circular Motions in the Outer Perseus Spiral Arm
    Thomas M. Dame
    December 30, 2023
====================================================


	
	
	
	
	
	
	
	
	
	
	
	
	
	
	The role of major mergers in galaxy evolution remains a key open question. Existing empirical merger identification methods use non-parametric and subjective visual classifications which can pose systematic challenges to constraining merger histories. As a first step towards overcoming these challenges, we develop and share publicly a new Python-based software tool that identifies and extracts the flux-wise and area-wise significant contiguous regions from the model-subtracted “residual” images produced by popular parametric light-profile fitting tools (e.g., GALFIT). Using Hubble Space Telescope (HST) H-band single-Sérsic residual images of 17 CANDELS galaxies, we demonstrate the tool’s ability to measure the surface brightness and improve the qualitative identification of a variety of common residual features (disk structures, spiral substructures, plausible tidal features, and strong gravitational arcs). We test our method on synthetic HST observations of a z∼ 1.5 major merger from the VELA hydrodynamic simulations. We extract H-band residual features corresponding to the birth, growth, and fading of tidal features during different stages and viewing orientations at CANDELS depths and resolution. We find that the extracted features at shallow depths have noisy visual appearance and are susceptible to viewing angle effects. For a VELA z∼3 major merger, we find that James Webb Space Telescope NIRCam observations can probe high-redshift tidal features with considerable advantage over existing HST capabilities. Further quantitative analysis of plausible tidal features extracted with our new software hold promise for the robust identification of hallmark merger signatures and corresponding improvements to merger rate constraints.
		
	
	
		Galaxies: evolution – Galaxies: statistics – Galaxies: high-redshift
	
	
	
	

§ MOTIVATION

	Merging of two similar-mass (stellar-mass ratio ≤4:1) galaxies is often referred to in the literature as major merging. Identifying and quantifying the rate at which galaxies experience such mergers (merger rates) over cosmic history is a key step to empirically quantify the role of major mergers in galaxy evolution. Conceptually, two broad observational approaches are employed to quantify major merging, namely: (1) close-pair methods  and (2) morphological methods. These approaches yield redshift-evolutionary trends that broadly agree with theoretical merger rate predictions that major mergers become more frequent with increasing redshift out to z<1.5. However, both methodologies suffer from unique, yet analogous systematic biases from selection effects and uncertain observability timescale assumptions, making it difficult to interpret the merger rates at z>1.5. <cit.> investigated the effects of such systematics in binary merger simulations using non-parametric identifiers of plausible merging (  Gini-M_20 and CAS metrics). In this study, we attempt to help overcome the subjective identification of `hallmark' merger signatures (  tidal arms, tails, bridges, and extended fans) and to robustly quantify their observability timescales  – the largest sources of uncertainty towards merger rate estimation in morphology-based methods – by developing a new public Python analysis tool to extract and quantify the morphological substructure of galaxies. When applied to both empirical and simulated data, this tool will facilitate future efforts to calibrate the observables associated with galaxy merging signatures and lead to improved merger rate constraints.
	
	Substantial efforts have attempted to constrain major merger rates as a function of cosmic time using empirical “merger fractions” based on close-pair statistics, typically identifying galaxy-galaxy pairs in close physical proximity <cit.> or by measuring clustering statistics <cit.>. This method leverages on the long-standing idea (supported by numerical expectations) that galaxies in close proximity become gravitationally bound and merge into a more-massive system <cit.>. Many previous close-pair-based works broadly agree that merger rates grow more frequent at earlier cosmic times during the redshift range 0≲ z≲1.5, albeit with a wide range of redshift dependencies ∝(1+z)^0.5-3 <cit.>. Such variance is largely caused by different study-to-study close-pair selection assumptions <cit.>. In <cit.>, we systematically quantified the effect of different close pair selection choices on the derived merger fractions.
	
	Recent efforts extend the close pair method to higher redshifts and find that merger rates may be un-evolving or diminishing with increasing redshift between 1.5≤ z ≤ 3 <cit.>. These measurements disagree with theoretical predictions of steadily increasing merger rates from z=0 to z=3 <cit.>. Leveraging recent simulation-based work by <cit.>, observational studies argue that the data-theory discrepancy at z>1.5 may be due to the assumption of simplistic (non-evolving) close-pair observability timescale employed in merger rate calculations <cit.>. What is needed is a thorough analysis of the observability timescale evolution for both close-pair and morphological methods. In a related study (Mantha et al. in prep), we are constraining the close-pair observability timescale as a function of different pair selection variables. The residual feature extraction tool that we introduce in this paper will aid new efforts to constrain the morphological feature observability timescales.
	
	Besides merger rates derived from close pairs, rates based on morphology have also been measured motivated by simulations showing that merging galaxies exhibit morphologically-disturbed appearance <cit.>. These merger constraints can be broadly categorized into rates based on visual identification of disturbed morphologies <cit.> and quantitative metrics of large-scale galaxy asymmetries or morphology such as CAS <cit.> and used by  <cit.>, and Gini-M_20 <cit.> and applied extensively by <cit.>. Early morphology-based studies find a broad range of merger rate evolutions R∝(1+z)^2-5 at z≲ 1.5 <cit.>, albeit with significant study-to-study scatter where some studies find no redshift dependence of merger rates <cit.>, and sometimes finding up to 25%-50% of their samples to be merging <cit.>. These studies suffer from low-number statistics due to small-volume pencil-beam surveys, redshift-dependent systematic biases induced by cosmological surface brightness dimming and strong morphological k-corrections when using rest-frame UV images, which can impact the purity and completeness of merger selection. Although the advent of large Hubble Space Telescope (HST) surveys like CANDELS <cit.> alleviate some issues, the observability timescales for morphology-based methods are highly uncertain and may be varying with redshift.
	
	
	A possible solution to overcome the systematic, purity, and timescale related challenges faced by the aforementioned morphology-based approaches is to employ more focused morphological signatures as merger identifiers. For example, <cit.> used the presence of multiple-nuclei separated by ≲8  kpc to select merging galaxies during a narrow window straddling the early and post-merger stages <cit.>, while aiming to minimize the contamination from spurious non-merging interlopers. Tidal features are a specific sub-class of morphological substructures that are predicted to be prevalent among merging galaxies <cit.>. A fruitful technique for improving the detection of faint and transient tidal features is through the analysis of residual images [Produced by subtracting the best-fit light profile of the host galaxy using tools such as GALFIT <cit.>.] either by visual inspection <cit.> or quantitative parameters <cit.>. In fact, galaxy light-profile fitting is a popular step employed in many large surveys to provide key insights into the structural (size, shape) evolution of galaxies <cit.>. In this context, our novel tool is designed to analyze the by-product residual images from galaxy image fitting routines and extract useful information about additional (often complex) substructures. 
	
	
	Detailed objective distinction of hallmark tidal features from other non-merging signatures (  lopsided disks, asymmetric or one-armed spiral structures) is yet to be achieved. Existing empirical merger-identification methods only provide a plausible indication of tidal signatures and do not attempt to quantitatively capture key information such as structure (length, width, and pitch), morphology (  tails and fans), color, and stellar mass of the features. Furthermore, constraining the role of observational effects (  depth, cosmological surface-brightness dimming) and observability timescales T_ obs for hallmark tidal features is an open question, and is a key hurdle in achieving robust tidal-feature-based merger rates.
	
	As a first step towards quantitatively identifying and extracting tidal signatures of galaxy merging, in this study, we introduce a new technique and its corresponding public Python software pipeline[The data products discussed in this paper along with the software are publicly available at <https://github.com/AgentM-GEG/residual_feature_extraction>] to extract and quantify  residual substructures and we demonstrate its application with a select sample of 16 CANDELS galaxies that exhibit a range of merging and non-merging morphological substructures. We also test our technique's ability to identify plausible interaction signatures by applying it to mock Hubble Space Telescope (HST) observations of a galaxy merger from the VELA zoom-in hydrodynamic simulations <cit.>. We structure this paper as follows: In <ref>, we introduce our CANDELS test galaxy sample and their GALFIT-derived data products, followed by a brief description of a galaxy merger from the VELA hydro-dynamic simulation and the generation of its mock HST observations.  In <ref>, we describe the step-by-step methodology of our residual substructure analysis pipeline, and we apply this method to the select CANDELS galaxies in <ref>. In <ref>, we apply our methodology to mock observations of a merger simulation and investigate the role of image depth and viewing angle on the feature extraction process. In <ref>, we present possible applications and limitations of our method. We present our conclusions in <ref>. Throughout this work, we adopt the <cit.> cosmological framework (H_0=67.7  km s^-1 Mpc^-1), and use the AB magnitude system <cit.>.
	
	
	
	
	
	

§ DATA AND TEST SAMPLES

	
	To develop and demonstrate our novel method for extracting and quantifying residual substructures, we use HST/WFC3 H-band (F160W) images and single-Sérsic GALFIT residual images from <cit.> for a select sample of galaxies from the CANDELS survey <cit.> that exhibit different residual substructures including plausible tidal signatures. To explore the application of our method for repeatable identification and extraction of bonafide tidal features, we also analyze mock images derived from a merging galaxy from the VELA zoom-in simulations <cit.>.
	

 §.§ CANDELS Data

	
	We make use of the H-band Sérsic-fitting data products from  (hereafter vdW12) for galaxies bright than H=24.5  mag in the CANDELS survey. Briefly, vdW12 used GALFIT <cit.> to generate an image-cube containing the original image used for fitting, the best-fit model image with the Sérsic model information stored in its header, and the residual (original-model) image. Following standard procedures outlined in <cit.>, vdW12 used <cit.> to extract sources from the CANDELS imaging <cit.>, estimated the background sky, and determined GALFIT parameter input values, and produced a rectangular image cutout of each galaxy based on its SExtractor parameters: x,y centroid, long side length of five (5) times the Kron radius, and orientated by the position angle θ. In practice, the postage stamp images typically have sides 20-30 times the galaxy's effective size (r_ e). To ensure optimal galaxy profile fitting with GALFIT, vdW12 masked neighboring objects 4 mag fainter than the central galaxy and simultaneously fit all remaining objects in the image cutout.
	
	The methodology we describe in <ref> starts with the vdW12 image cubes for a select sample of 16 galaxies from an ongoing effort (HST-AR 15040; PI: McIntosh) to visually characterize the H-band residual substructures hosted by CANDELS galaxies with 9.8 ≲(M_ stellar/M_⊙) ≲11 and 0.5≲ z ≲ 2.5. 
	As shown in Table <ref>, these galaxies sample different kinds of residual features: disk (4 galaxies), spiral arms (4 galaxies), plausible tidal signatures (7 galaxies), and gravitational strong-lensing arc (1 galaxy). We also note that our galaxies span two observational depth regimes – CANDELS/Wide (HST 2-orbit; 6 galaxies) and CANDELS/Deep (HST 10-orbit; 10 galaxies) depths as described in <cit.>. For each galaxy in our sample, we extract residual substructures using the images with best-available depth, unless otherwise specified. In Figure <ref>, we provide example postage-stamps of our test set galaxies and their respective residual features to be extracted and analyzed in this paper. For reference, we also show an example `clean' (i.e., visually no left-over) residual image.
	
	
	
	
	

	

 §.§ VELA Datasets 

	
	To demonstrate our new methodology on the extraction of residual substructure signatures hosted by a galaxy-galaxy major merger, we apply it to mock images from VELA#01[The VELA simulation-suite comprises of 35 intermediate-mass galaxies (log_10(M_ stellar/M_⊙) ∼ 10 at z∼ 2) and they are referred using an ID starting from 01 to 35.] of the VELA zoom-in hydro-dynamic simulations <cit.>. 
	Briefly, in this simulation baryonic physics[The VELA simulations include many physical prescriptions related to galaxy evolution – gas cooling, star formation, ISM metal enrichment, photo-ionization heating from the Ultra-Violet background radiation, stellar-mass loss, stellar feedback, tracing of smooth cosmological accretion of other galaxies and gas.] are solved in conjunction with gravitational N-body interactions using the Adaptive Refinement Tree (ART) code <cit.> with spatial resolution of 17.5  pc to 35  pc. The VELA simulation-suite has been extensively used in the context of galaxy morphology such as elongated galaxies <cit.>, giant clumps <cit.>, compaction <cit.>, galaxy spin and size <cit.>, H-alpha morphology <cit.>, and structural parameters <cit.>. Starting at a cosmological scale factor a= 0.125 (z∼7), the simulation information is recorded at equal increments of Δ a = 0.01, which approximately corresponds to a time frame Δ t=100  Myr at z=2. For the merging system VELA 01, we use the information during 15 time intervals spanning 0.37≤ a ≤0.47 (1.7≲ z ≤ 1.1), bracketing roughly 0.6 (0.8) Gyr before (after) the coalescence stage (occurring at a=0.41; z=1.44). We note that this system experiences a gas-rich, major merger (stellar mass ratio ∼1.33:1, gas mass ratio ∼ 1:1) resulting in post-merger stellar and gas masses of (M_ stellar/M_⊙) = 10.05 and (M_ gas/M_⊙)=9.55, respectively.
	
	At each time step, we use mock HST WFC3/F160W observations of the merging system generated using the methodology described in <cit.>. Briefly, they use the raw-simulation information at each time-step and assign stellar population synthesis model-informed <cit.> spectral energy distributions (SEDs) to the stellar particles using SUNRISE <cit.> with dust scattering and absorption assuming a dust-to-metals ratio of 0.4 <cit.> and a Milky Way-like dust grain size distribution <cit.>. At each time-step, they integrate this 3-dimensional data onto a 2D camera plane to produce an idealized image at  19 (camera) orientation angles <cit.>, which are then convolved with the point spread function (PSF) generated by Tiny Tim <cit.> to reach desired spatial resolution. These mock idealized images are available on MAST[<https://archive.stsci.edu/prepds/vela/>] <cit.>.
	
	
	To match the empirical CANDELS F160W 5σ limiting surface-brightness sensitivities at both Wide (25.25  mag arcsec^-2) and Deep (26.25  mag arcsec^-2) depths, we add Gaussian random noise to the F160W PSF convolved images of VELA 01. This process of generating mock HST CANDELS images from simulated data products is often referred to as CANDELization <cit.> and we call the resultant depth-matched images as CANDELized images.  We apply GALFIT to the VELA 01 CANDELized images to generate single-Sérsic residual images consistent with the vdW12 empirical data at each time step, viewing angle, and depth realization for a total of 380 unique residual images.
	
	
	
	
	
	
	
	
	

§ METHODOLOGY

	
	We demonstrate our residual substructure extraction method by analyzing the GALFIT-based residual images of a selective sample of galaxies hosting different left-over substructures (<ref>). 
	We organize the functionality of our method into three main objectives as follows. Starting with a H-band original image, we generate two masks to facilitate estimation of local-sky background (σ_ bkg) and residual feature extraction. Next, we estimate σ_ bkg in the smoothed single-Sérsic-based residual image and search (in proximity to the galaxy-of-interest) for contiguous-pixel regions that satisfy a sky-significance. Finally, we identify the flux-wise significant regions whose areas are larger than a pixel-contiguity area threshold (A_ thresh) informed using noise-only Monte-Carlo (MC) simulations to quantify their H-band surface brightness. We present a detailed flow-chart of our residual extraction process in Figure <ref> and label the key steps to be discussed in this section.
	
	We acknowledge that our methodology is closely dependent on the outputs generated by different software routines such as GALFIT and Source Extractor. Although these software are extensively used in the astronomy community, their specific functionality may change in the future updated versions. In this context, we emphasize that our modular algorithmic approach described in <ref> and illustrated in the flow chart (Figure <ref>) can be treated as a road map for future applications using the updated or new versions of the software analogous to GALFIT or Source Extractor.
	
	
	
	

 §.§ Masking
 
	
	We devise a dual masking scheme to facilitate the extraction and isolation of residual substructures at radii outside the each galaxy center. This requires robust estimation of the local-sky background, detection of faint sources not fit by vdW12, and a definition of galaxy center. The magnitude criterion imposed by vdW12 to mask nearby faint objects (<ref>) is a reasonable choice in the context of quantifying galaxy structural properties. However, such faint un-fit sources can have comparable surface brightness to our desired faint residual substructure (  tidal features).
	In addition, central galactic regions may house complex substructures (especially in young post-mergers) that may not be well represented by a Sérsic profile and may lead to significant central left-over light <cit.> that can obscure the detection of faint tidal features of interest. 
	
	

  §.§.§  (Step#1) Generate sky mask

	
	
	
	We start by performing faint source detection within the postage stamp of the galaxy 
	(original image) using a python-based implementation <cit.> of the standard source-extraction tool SExtractor <cit.>. 
	We find that the source extraction parameter values provided in Table <ref> yield optimal identification of small faint sources in the images, without deblending faint extended features (  plausible tidal signatures) into separate objects[A small value of DETECT_THRESH ensures that the resultant segmentation maps envelop the outer-most light of the sources An alternative to this step is to use a higher detection threshold value and grow the resultant segmentation regions to fully encompass the extended light.].
	
	Compared to the SExtractor configurations used to generate the CANDELS source catalogs <cit.>, we employ a smaller convolution filter size (five vs. nine) for better detection of smaller and fainter sources. 
	
	In the top panels of Figure <ref>, we demonstrate the thorough detection of all faint sources seen in the vdW12 residual image of an example galaxy (GDS 14637; center of image). Hereafter, we use this galaxy to demonstrate all of our analysis steps outlined in the flow-chart (Figure <ref>). In the bottom left panel, we provide the resultant background sky mask. Masked regions are shown in black (pixel value=0). The unmasked regions (white) represent the sky contribution in the vicinity of the galaxy, which use during our feature extraction in <ref>.
	
	
	
	
	
	
	
	

  §.§.§ (Step#2) Generate residual extraction mask

	We also generate a mask for extraction of residual substructures within a physically-motivated proximity of R_ proj < 30  kpc from the centroid of the galaxy-of-interest.
	
	This corresponds to nearly five (ten) times the effective size of late (early) type galaxies at z∼ 1 <cit.>. 
	We find that this extraction radius encompasses all faint residual substructure for the galaxies in our sample. As shown in Figure <ref> (bottom-right panel), we mask contributions
	from the central core regions of each galaxy residual so that we extract only outer features. We define the inner masking by an ellipse centered on the galaxy, with semi-major and semi-minor axes equal to 1.5 times those of the SExtractor detection ellipse from step#1, and orientation along the SExtractor position angle. Finally, we also mask any faint companions from step#1 that fall within the R_ proj < 30  kpc extraction region. In summary, 
	
	we use this mask during the extraction of residual substructures (<ref>).
	
	
	
	
	
	

 §.§ Final Extraction of Residual Features

	
	
	We start the feature extraction process by analyzing a Boxcar-smoothed version of the residual image (top-left panel; Figure <ref>) to facilitate the robust identification of faint features. 
	We utilize the masks described in <ref> to identify all unmasked pixels within the residual extraction mask region that are
	significantly above the background sky. Then, we identify contiguous groupings of such pixels to define unique residual features.
	
	

  §.§.§ (Step#3) Identify flux-wise significant pixels

	
	The first step in extracting residual features is to quantify the background sky. We compute the background σ_ bkg as the standard deviation of the best-fit Gaussian profile representing the normalized distribution of pixel values in the unmasked regions of the smoothed residual image after applying the sky mask (<ref>). In Figure <ref> (top-right panel), we show the sky pixel-value distribution (blue histogram) with its corresponding Gaussian curve. We denote the 2σ_ bkg with a vertical dashed line.
	
	Next, we apply the residual extraction mask and identify all unmasked pixels in the smoothed-residual image (red histogram in  Figure <ref>, top-right panel). In Figure <ref> (bottom-left panel), we plot all such pixels satisfying a fiducial flux-wise significance cut f_ pix>2σ_ bkg, where f_ pix is flux value of the pixel. Visually, we observe one large flux-wise significant residual feature and many small areas of >2σ_ bkg pixels. We note that a less stringent significance cut of f_ pix>1σ_ bkg increases both the number of small pixel groupings and the typical area per pixel-group, albeit with decrease in the overall signal-to-noise ratio. On the other hand, a more strict f_ pix>3σ_ bkg cut identifies fewer pixels (with high signal-to-noise), but with significant structural discontinuities (  an extended feature is broken into multiple smaller regions). As such, we elect to use f_ pix>2σ_ bkg as our fiducial criterion for flux-wise significant pixels in our analysis.
	
	
	
	

  §.§.§ (Step#4) Identify contiguous-pixel regions

	
	To isolate and separate larger, contiguous residual substructures from smaller groupings of flux-wise significant pixels, we use a python-based 2D pixel-connectivity algorithm skimage: label <cit.>.  This tool assigns unique identification number to each contiguous pixel region using a binary image in which all the flux-wise significant pixels are set to unity and the remainder are set to zero. In Figure <ref> (bottom-left panel), we show the contiguous flux-wise significant pixels (color coded regions).
	
	Next, we compute the area and enclosed flux of these contiguous regions using a companion module skimage: regionprops[The module computes cumulative flux per unique contiguous region in the residual image at their respective pixel coordinates.]  <cit.>. Re-affirming our visual interpretation, we find a contiguous substructure with large area (shown in green) in the immediate vicinity of our example galaxy and many regions with smaller areas (Figure <ref>; bottom-left panel). We note that correlated-noise fluctuations can cause several pixels to be flux-wise significant and contiguous by random chance. We quantify this effect and reject such spurious regions in the following section.
	
	
	

 §.§ Final Feature Selection

	
	
	We isolate the final, flux-wise and area-wise significant residual features following a MonteCarlo approach. We estimate an area threshold (A_ thresh) based on a noise-only pixel contiguity expectation and select those flux-wise significant contiguous features with areas larger than A_ thresh.
	
	

  §.§.§ (Step#5) Quantify area threshold A_ thresh

	As a first step to quantify an area threshold above which the noise contribution to pixel contiguity is minimum, we generate 200 Boxcar smoothed images with Gaussian noise matching the σ_ bkg value from step#3. In each random-noise image, we identify contiguous regions satisfying our fiducial flux-wise significance and quantify their area (as described in <ref>).  In Figure <ref>, we show an example smoothed random noise image (left panel) and its corresponding flux-wise significant contiguous regions (right panel; color coding). We notice that the noise-generated contiguous areas visually mimic the smaller area regions identified in step#4 (Figure <ref>; bottom-left panel). Combining the result of 200 MonteCarlo iterations for our example galaxy, we find an extraction area threshold A_ thresh = 20, where > 99% of the contiguous pixel regions in smoothed noise-only images have areas smaller than 20 pixels. We repeat this MonteCarlo process for all galaxies in our analysis that span different CANDELS fields and find a range of A_ thresh values between 17 to 23 pixels, owing to varying field-to-field noise properties.
	
	
	

  §.§.§ (Step#6) Identify and Quantify Final Residual Features

	To select final significant features, we identify the contiguous regions with areas larger than the threshold set by the noise properties (  A>A_ thresh). We then quantify the surface brightness of the final feature regions (μ_ feat) as:
	
    μ_ feat = -2.5(F) + zp + 2.5(A^') [ mag arcsec^-2] ,

	where F is the cumulative flux (in Jansky) of all the regions satisfying our criteria, zp = 8.9 is the zero-point magnitude, and A^' is the cumulative area of all these significant contiguous regions in units of arcsec^2. We also report the photometric error on the quantified surface brightness by adding the noise from σ_ bkg and Poisson noise from the feature in quadrature. In Figure <ref> (bottom-right panel; shading), we show the final residual features (an extended tidal-tail structure) and quantify its H-band surface brightness to be μ_ feat = 25.5±0.3  mag arcsec^-2.
	
	
	
	
	
	
	
	
	

§ APPLICATION TO CANDELS GALAXIES

	
	In this section, we demonstrate the application of our residual substructure extraction method on example CANDELS galaxies hosting a range of features – disk, spiral, and plausible interaction signatures (see Table <ref>). We also illustrate the role of image depth by comparing results at HST 2-orbit and 10-orbit depths for three galaxies (GDS 4608, GDS 14637, and GDS 14876). 
	
	
	

 §.§ Extracting disk and spiral features

	
	In Figure <ref> (top-two rows), we show the results of applying our feature extraction method on galaxies visually identified to be hosting residual disk structure. In the cases of  EGS 10518 and GDS 18959, we find similar looking, inclined disk-like residual structures, whereas the residuals of EGS 27018 and EGS 12782 exhibit circular features (likely owing to their face-on orientation). We also note hints of clump-like substructures embedded within the disk feature of EGS 12782. We find that these features have surface brightness measures ranging between 23.4 - 24.5  mag arcsec^-2.
	Re-affirming our visual interpretation, we find that our method extracts features that are consistent with a residual disk, likely owing to an underlying bulge+disk morphology being fit with a single-Sérsic model.
	
	In Figure <ref> (bottom-two rows), we show the results of our feature extraction on galaxies hosting residual spiral substructures (EGS 1690, GDS 18160, UDS 4165, and GDN 7220). We extract two-arm spiral patterns for all the four galaxies and find their surface brightness ranges between 23.6-24.9  mag arcsec^-2. We note that the spiral substructures exhibited by GDS 18160 and GDN 7220 are visually apparent (despite having fainter surface brightness) than EGS 1690 and UDS 4165, owing to the difference in their observational depths (deep vs. shallow; see Table <ref>). We also note that all except GDN 7220 appear to be hosting a central bar. Extracting such spiral substructures with our method will help future analyses quantify their physical parameters such as pitch angle, bar structure.
	
	
	
	
	
	

 §.§ Extracting plausible tidal features

	In Figure <ref>, we show the  results of our feature extraction method on galaxies exhibiting plausible tidal interaction signatures. These galaxies broadly fall into one of the following four categories – i) close-pair systems in which both galaxies exhibit plausible signatures of tidal interaction (GDN 14758, UDS 24437, and GDS 4608); ii) close-pair systems in which only one galaxy exhibits a plausible interaction feature (GDS 16440 and GDS 14637 ); one example galaxy experiencing a plausible minor interaction (GDS 14786);  and one plausible post-merger system (GDN 19863). The extracted residual substructures in these examples highlight our visual interpretations.
	
	First, we apply our feature extraction method to each galaxy in the three example close-pair systems with plausible dual interaction signatures. We note that all three pairs have |Δ z|/(1+z) ≤ 0.025 and residual substructure consistent with theoretical predictions for tidal features produced by major mergers <cit.>. We find a striking structural similarity between the features exhibited by UDS 24437 – UDS 23110 pair (z∼ 1.8) and the GDS 4608 – GDS 4529 pair (z∼ 1), such that each close pair hosts two tidal arms and an overlapping region between the host and companion galaxy. These broad and extended features are similar to those found in observations and simulations of major mergers between spheroid-dominated galaxies <cit.>. The residual substructure midway between each galaxy core in these two pairs may be a stellar bridge or may be an artifact of the simultaneous Sérsic fitting of two galaxies in very close proximity. In contrast, each galaxy in pair GDN 14758–GDN 14516 appears to have two tidal arms which are thin and clumpy in some places. The characteristics of the latter tidal arms are likely the result of interacting disk-dominated galaxies. 
	
	Next, we demonstrate the feature extraction on two pairs in which only one galaxy exhibits plausible interaction signatures. GDS 14637 has a single arm-like residual that is similar in appearance to the broad and extended tidal arms discussed above. In contrast, GDS 16440 shows two arm-like features, where one extends outward away from the host and is aligned with a nearby (compact) galaxy (at right) with a spectroscopic redshift difference of only |Δ cz|∼ 100  km s^-1 from GDS 16440. This picture is consistent with theoretical merger simulations <cit.> of a disk-dominated galaxy experiencing a major merger with a nearby compact companion. We note that GDS 16440 has a clear inner spiral pattern that extends to a radius of R∼5–10 kpc. Despite its point-like appearance suggesting a plausibly interaction-fueld AGN, the companion has no known X-ray detection to corroborate this suggestion. Finally, we acknowledge the possibility that the extracted arm-lke features of GDS 16440 may be extensions of the inner spiral structure, yet, their diffuse appearance and physical extent out to R>25 kpc is much more consistent with tidal interaction signatures. Admittedly, quantitatively distinguishing spiral versus interaction features is critical for robust identification of galaxy merging. While beyond the scope of this work, we believe that further quantification of extracted residual features will provide improved methods for distinguishing distinct physical processes.
	
	Finally, we demonstrate our residual feature extraction technique on two additional galaxies with plausible signs of tidal activity. First, GDS 14786 shows a minor companion galaxy with a narrow ∼ 10  kpc tail-like feature, consistent with numerical predictions of features caused by minor mergers  <cit.>.
	Second, GDN 19863 appears to have two tidal-tail residuals, each extending ∼ 10  kpc from the galaxy center. One is much fainter than the other. Multiple faint tidal tails are found briefly after coalescence in simulations of major mergers.
	
	

 §.§ Role of Depth during Feature Extraction

	
	
	To illustrate the role of imaging depth on the extraction and quantification of different plausible tidal features, we compare the results of residual substructures extracted from deep (CANDELS/Deep) and shallow (CANDELS/Wide) images for three examples in Figure <ref>.  We perform this comparison on the following examples from <ref>: GDS 4608 (close-pair system where both galaxies show plausible tidal arms), GDS 14637 (close-pair system where on galaxy shows plausible interaction feature), and GDS 14876 (galaxy experiencing a plausible minor interaction).
	
	We find that the morphologies of the extracted residuals at a shallow 2-orbit depth are qualitatively similar to their 10-orbit depth counterparts in all the three cases, albeit with noticeably noisy visual-appearance. We note that the features extracted from the shallow depth images tend to be confined to the brighter surface brightness (by 0.3-0.5  mag arcsec^-2) regions of the same features extracted from the deep images. The shallow residuals are also lower in signal-to-noise when compared to features extracted at a deeper depth. This is expected as an increase in noise per pixel reduces the structural contiguity, owing to the reduction in number of pixels meeting flux-wise significance cut, thereby, resulting in a less significant yet brighter surface brightness.  Although we re-calibrate the feature extraction employing a larger area threshold value for shallow depth images, we notice additional small faint sources in the nearby projected vicinity that are contaminating the extracted features. This is because such faint sources are not significantly detected during the source extraction and masking steps, but are extracted from the residuals by our method. This deep-shallow comparison highlights the role of  observational depth in identifying faint, high-redshift plausible tidal signatures.
	
	
	
	
	

§ EXTRACTING TIDAL FEATURES IN A MAJOR MERGER FROM COSMOLOGICAL SIMULATION

	
	In this section, we demonstrate the application of our residual substructure extraction method on synthetic HST images of a major merger from the VELA simulation <cit.>. We extract the residual features and quantify their areas during 10 time steps of the merging process spanning 1.7≲ z ≤ 1.1 (Δ t∼1.4  Gyr), each at 19 viewing orientation angles, and two observational depth configurations matching empirical CANDELS observations (Wide and Deep). We note that the first pericentric passage and the coalescence stages occur at z=1.63 and z=1.44, respectively. To visually comprehend the quantitative trends, we also discuss the qualitative appearance of example extracted features at one fixed viewing orientation.
	
	

 §.§ Quantitative Analysis of Extracted Features

	
	In Figure <ref>, we show the median extracted feature areas as a function of time before and after the coalescence stage (t=0  Gyr) for CANDELS Wide and Deep depths at all 19 viewing angles. We also quantify the detection probability as the percentage of viewing sightliness during which we detect features with their areas larger than the A_ thresh and μ_ feat brighter than the limiting surface brightness (Wide – 25.25  mag arcsec^-2;  Deep – 26.25  mag arcsec^-2). We find the deep and shallow extracted residuals follow qualitatively similar trends in terms of their changing prominence (i.e., feature area) during different pre- and post-merger stages. Owing to coarse time resolution (∼ 150  Myr) between each time step, we are only able to discuss broad trends and differences between the deep and shallow results. 
	
	At deeper image depth, we detect residual substructure at nearly all viewing angles and time steps. Prior to first pericentric passage (t=-0.46  Gyr), the feature area is minimal (A = 3± 1  A_ thresh). At first pericentric passage, we find a plausible increase in the feature prominence, albeit with a large spread between different viewing orientations. This increase clearly peaks at the next time step 150 Myr later, hinting at the formation of tidal interaction-triggered features. If true, these features abruptly diminish in significance (by a factor of two) at t=-0.15  Gyr, after which extracted feature areas steadily increase through the coalescence stage until they peak at 300 Myr later with a value 3–7 times greater than the initial minimum. This temporal behavior strongly implies the presences of tidal merger signatures. During the three time steps that follow (0.47≤ t ≤ 0.79), we note that the feature prominence follows a stochastic rising and falling trend with a range of area values between those of coalescence and the last maximum, or even slightly more at some viewing angles. As we discuss later, these stochastic fluctuations during different time steps correspond to noticeable changes in the visual appearance of the extracted residual features.
	
	
	In the case of shallow depth images, we typically find a similar temporal time evolution trend, albeit with some key distinctions. First, we note that the trends at shallow depths are shifted to systematically smaller areas when compared to the same features observed in deeper images. This is because of the decreased signal-to-noise per pixel causing less number of contiguous pixels to satisfy our feature extraction criteria. Next, contrary to what we observe in deeper images, we find no increase in the prominence of the extracted feature at the t=-0.31  Gyr time step following first pericentric passage. This is consistent with the idea that tidal interaction features are faint and their detection is highly sensitive to image depth. Indeed, we find the detection probability of any residual substructure in the shallow images is less than 50% until after first passage, and only reaches unity at the coalescence stage when the first noticeable increase in feature area occurs. During t>0  Gyr, we note that post-merger time steps with detection probability <100% correspond to smaller feature areas. Likewise, maximum feature areas correspond to high detection probability. These subtle trends highlight the role played by viewing lines-of-sight during the detection of certain merger features, especially in shallow depth images. Nevertheless, during the time steps with high detection probabilities (irrespective of the depth), we note that the viewing angles can induce a scatter of 20-40% in the measured feature prominences.
	
	
	
	
	
	
	
	

 §.§ Qualitative Interpretation of Extracted Features

	
	To illustrate the temporal quantitative trends shown in Figure <ref>, we discuss the qualitative appearances of the features extracted from one example viewing orientation. We show the extracted features at shallow and deep depths alongside their idealized images during 10 time steps in Figure <ref>. We discuss the extracted features from different image depths with respect to the distinct structures evident in their idealized (no noise added) images, and we compare the morphologies of the features extracted from the VELA mock data to those from CANDELS images of plausible merging systems shown in Figure <ref>.
	
	During the first two time steps (t≤ -0.46  Gyr), we extract small residual structures in both deep and shallow images that correspond to faint sources in proximity to the merging galaxies within their idealized images. At t=-0.31  Gyr (∼ 150  Myr after the first pericentric passage), we extract features at deeper depths that correspond to the interaction debris surrounding the merging galaxies  in the idealized images. However, we do not detect these features in shallow depth images, which aligns with the noticeable difference in feature prominence at t=-0.31  Gyr in Figure <ref>. Just prior to coalescence (t=-0.15  Gyr), we extract a small feature midway between the two merging galaxies at both depths, which corresponds to a plausible bridge-like structure in the idealized image. We find similar excess-light features in the extracted residuals for UDS 24437 and GDS 4608 in Figure <ref>. Yet, interestingly, the VELA merger does not exhibit tidal arms during its close-pair stage, which is in contrast to all three example plausible interacting systems (GDS 4608, UDS 24437, and GDN 14758). 
	
	
	Starting at coalescence (t=0  Gyr), we extract features at both depths that highlight the birth of a post-merger tidal-arm that grows in area over the next ∼ 300  Myr. This feature is similar in appearance to the single arm-like structure we extract from one of the two example CANDELS pairs in which only one galaxy has a plausible interaction signature (GDS 14637; Fig. <ref>). As such, GDS 14637 is also consistent with a post-merger in close proximity to another galaxy.
	At t=0.47  Gyr, we extract two narrow tail features with clumpy visual appearance in the deep image, which resemble the plausible post-merger signature we extract for GDN 19863. This picture is consistent with the theoretical expectations that major gas-rich mergers produce diffuse tidal features with clumpy knots <cit.>. We note that we only marginally detect a fraction of this tidal feature in the shallow depth image, which aligns with the abrupt decrease in feature prominence noted in Figure <ref>. Finally, during the last two time steps in both shallow and deep images we extract a large and extended plausible tidal loop or two-tailed feature (at t=0.63  Gyr) that becomes less prominent at t=0.79  Gyr. The changes in these features during the last two time steps illustrate the stochastic fluctuations in the extracted feature areas we find for all viewing angles in Figure <ref>. This exercise illustrates the interplay between interaction/merger feature extraction and observational depth.
	
	
	
	

§ DISCUSSION

	
	So far, we have presented the development of a residual substructure extraction tool and discussed its application to example galaxies from the CANDELS survey and mock CANDELized observations of a galaxy merger simulation. In this section, we discuss several additional possible applications of our residual extraction method.
	
	
	

 §.§ Possible Applications

	We present some possible applications of our residual feature extraction method towards improving popular structural-fitting routines, studying multi-wavelength properties of different single-component Sérsic residual features (  disk structures from host bulge+disk systems) and substructures (  spiral arms, plausible interaction signatures), and extracting strong gravitational arcs.
	
	

  §.§.§ Improving structural fitting
 
	An important step in many galaxy structural-fitting studies <cit.> is to choose the number of required parametric components during GALFIT modeling (or other similar tools), starting with a simple single-Sérsic model and advancing to more complex models (  bulge+disk) including even the inclusion of Fourier modes (  to model spiral features). For upcoming large-scale surveys, developing machinery to make such decisions is key for automated structural modeling for very large samples. One way to approach this is to use machine learning to inform the number of parametric components required to best represent a galaxy morphology <cit.>. Alternatively, we propose that our method can be used to inform additional residual components and repeat  GALFIT with updated number of components. 
	
	

  §.§.§ Multi-wavelength application
 
	Empirical studies often perform multi-wavelength decomposition of galaxy structure into their bulge and disk components to separately inform their physical parameters <cit.>, which can provide key insights into the physics of galaxy evolution. We propose that extending our feature extraction method to multiple wavelengths can inform the physical properties for disk, central bulge, spiral arms, and interaction signatures. Furthermore, recent studies are developing resolution-element-based spectral energy distribution fitting in conjunction with spectroscopic measurements to produce maps of physical parameters <cit.>. Using our method in conjunction with these semi-resolved maps can enable future investigations to study structurally-resolved properties (e.g., colors, star-formation, chemical composition, kinematics) of different galactic substructures <cit.>.
	
	
	

  §.§.§ Extracting Gravitational Arcs

	
	Yet another possible application of our method is the extraction of strong gravitational lensing arcs. Strong gravitational lensing provides an unique opportunity to observe the otherwise inaccessible high-redshift universe. Several previous studies have used strong lensing to study cosmology, the physical properties of lensing systems, and the properties of distant background lensed sources <cit.>. Identifying and extracting the gravitational lensing signatures is a key step in such works, which often use GALFIT to subtract the host galaxy light to unearth such lensing features <cit.>.
	
	In Figure <ref>, we show the result of our feature extraction method on an example CANDELS strong-lensing candidate (UDS 10713), which is identified using a machine learning based classification scheme by <cit.>. We use the galaxy's single-Sérsic residual image by vdW12 and extract significant contiguous substructure above a 3σ_ bkg threshold (instead of our fiducial 2σ_ bkg). We find arc-like features with surface brightness μ_ feat=24.3±0.2  mag arcsec^-2 in close vicinity to the galaxy, which qualitatively aligns with the blue-arc in its corresponding false-color RGB image (at roughly 10 o'clock). However, we note that the single-Sérsic residual of UDS 10713 shows excess galactic light and additional sources in its projected proximity, which are contributing to the surface brightness of our extracted features. Isolating such contaminants is beyond the scope of this current work and we recommend that future studies should employ more complex Sérsic models for better galactic light subtraction and use multi-wavelength information of the features for better automated lensing feature detection and extraction. Furthermore, future works can extend our residual extraction method to galaxy cluster environments where strong lensing can be more prevalent (usually with high magnifications) as a means to derive cluster mass measurements.
	
	
	
	
	

 §.§ Future Prospects with the James Webb Space Telescope

	
	In preparation to probe high-redshift galaxy merging using the James Webb Space Telescope (JWST), we analyze synthetic JWST observations of a different merger in the VELA cosmological simulation suite (VELA 21), which experiences a major merger (post-merger stellar mass log_10(M_ stellar/M_⊙) =10.5; stellar-mass ratio ∼ 3:1) at z∼ 3.2 (coalescence). In Figure <ref>, we compare and contrast the feature extraction process on its mock JWST and HST observations. 
	
	Motivated by our simulated observations of a major merger in <ref>,  we focus this exercise on one specific merger stage z=3 (∼120  Myr after coalescence), where we expect the galaxy to exhibit tidal features. We generate idealized mock HST WFC3/F160W and JWST NIRCam/F277W images using the process described in <ref>. We add Gaussian noise to them such that they match the CANDELS/Deep depth and the proposed Cosmic Evolution Early Release Science Survey (CEERS; PI: S. Finkelstein, JWST PID: 1345) depth [The proposed depth of the CEERS survey reaches greater than a 10σ point-source and 5σ extended-source sensitivity of 28  mag with a ∼ 2800 second exposure in NIRCam/F277W.], respectively. Finally, we generate the single-Sérsic residuals by applying GALFIT to the HST and JWST noise-added images following vdW12, and apply our feature extraction technique on these residuals (<ref>).
	
	While we detect no significant residual features at CANDELS/Deep HST depth, we extract extended tidal-fan features in the NIRCam/F277W JWST CEERS depth image, which mimic the structures exhibited by our example CANDELS close-pair system (GDS 4608; see Figure <ref>). Our demonstration asserts the capabilities of JWST to probe faint tidal features out to high-redshifts with considerably less exposure time and at a longer rest-frame wavelengths than the current extremes accessible to HST. This will enable future studies to probe high-redshift merging systems that are key towards understanding the role of galaxy merging in galaxy evolution. 
	
	
	
	
	
	
	
	
	
	
	
	
	

§ CONCLUSIONS

	
	As a first step towards quantitatively identifying `hallmark' signatures of galaxy merging (  tidal arms, tails, bridges, and extended fans), in this paper, we introduce a new analysis tool to extract and quantify residual substructures hosted by galaxies. Our method analyzes residual images produced by any standard light-profile fitting tools and extracts contiguous features that are significant above the sky background and a noise-only expectation of pixel-contiguity. When applied to both empirical and simulated data, this tool will facilitate future efforts to calibrate the observables associated with galaxy merging signatures and provide improved merger rate constraints.
	
	We apply our residual feature extraction method to single-Sérsic residual images of 16 select galaxies from the CANDELS survey <cit.> and extract different structures hosted by these galaxies  – disk structures (4 galaxies), spiral substructures (4 galaxies), and plausible interaction/merger-triggered tidal features (8 galaxies). The main conclusions of our empirical exercise are:
	
	∙  We extract residual disk features and two-armed spiral substructures for all example galaxies visually identified as hosting disk and spiral features, respectively.
	
	∙ We extract plausible structurally-similar, dual tidal interaction signatures hosted by three example close-pair systems, two close-pair systems where only one galaxy exhibits plausible interaction signatures, one system experiencing a plausible minor merger event, and one galaxy hosting plausible post-merger signature.
	
	∙ To illustrate the role of imaging depth during feature extraction, we compare the residual substructures extracted from deep (HST 10-orbit) and shallow (HST 2-orbit) images of three example galaxies hosting plausible tidal features. We find that the features extracted at shallow depths are qualitatively similar but noisier in appearance compared to their deeper counterparts owing to the limited sensitivity of shallow CANDELS data to only the higher-surface brightness substructure of the features.
	
	We also apply our tool to synthetic HST observations of a galaxy merger from the VELA hydro-dynamic simulations <cit.>. We extract the residual features and quantify their areas during 10 time steps of the merging process spanning 1.7≲ z ≤ 1.1 (Δ t∼1.4  Gyr), each at 19 viewing orientation angles, and two observational depth configurations matching empirical CANDELS observations (Wide and Deep). We also illustrate the qualitative appearance of these residual features at one fixed orientation. The main conclusions of our theoretical exercise are:
	
	∙ We find that the features typically rise in prominence (feature area) ∼ 150  Myr after the first pericentric passage (z=1.63; time since coalescence t=-0.46  Gyr), which visually correspond to close-pair interaction signatures. These features abruptly diminish by a factor of two in significance during t=-0.15  Gyr before coalescence.
	
	
	∙ Starting with the coalescence stage (z=1.44; t=0  Gyr), we extract post-merger features that grow in prominence until peaking at 3-7 times the initial minimum value after 300  Myr. These stages visually represent the birth and growing prevalence of tidal arm features. At later time steps 0.47≤ t ≤ 0.79  Gyr, the feature prevalence follows a stochastic rising and falling trend, where the extracted tidal features appear structurally fragmented and clumpy appearance. 
	
	∙ Imaging depth and viewing orientation play a key role in detecting certain merger features. At deep image depths, we extract all the features at nearly all viewing angles (∼ 100% detection probability). However,  the same features in shallow depth images are detected at ∼ 50-75% of the viewing sightliness during -0.46 ≤ t ≤ -0.15  Gyr prior to coalescence. They only reach 100% detection probability at t>0  Gyr, when the first noticeable increase in feature prominences occur. Additionally, the temporal trends in shallow depth images are shifted to systematically smaller areas than their deeper depth counterparts. At a fixed depth, the measured feature prominences between different viewing orientations can vary by 20-40%.
	
	While the application of feature extraction method provides important insights about the onset and evolution of tidal features hosted by an example galaxy major merger from the VELA hydro-simulations, we caution the extrapolation of our quantitative (e.g., observability timescales) and qualitative (e.g., tidal feature shapes) interpretations to a general sample of mergers. This is because the strength and appearance of tidal debris may depend on the intrinsic properties of merging galaxies such as: stellar mass-ratio, redshift (driven by the evolution gas-fractions), merger orbital configurations (e.g., prograde vs. retrograde approach). Moreover, it is also important to note that the tidal feature observability carried out in a cosmological context may differ from a similar exercise on idealized binary-merger simulations.
	
	Lastly, we discuss two additional possible future applications of our general purpose tool. First, we illustrate the extraction of strong gravitational arcs in a candidate strong lensing system identified by <cit.>.  Second, we demonstrate the future prospect of using our method to probe high-redshift tidal features by applying it to synthetic JWST observations of a z∼ 3 simulated major merger. We find that JWST can probe faint, high-redshift tidal features with a smaller exposure time investment and longer rest-frame wavelengths (2.7 μ m) than HST/WFC3 H-band 10-orbit (i.e., CANDELS/Deep) observations.
	
	

§ ACKNOWLEDGMENTS

	We thank Roberto Abraham for the constructive referee suggestions and feedback that improved this work. We thank Mark Brodwin, Gabriela Canalizo, Alexander De la Vega, Boris Häußler, Marc Huertas-Company, Kartheik Iyer, Marziye Jafariyazani, Erin Kado-Fong, Allison Kirkpatrick, David Koo, Peter Kurczynski, Bret Lehmer, Jennifer Lotz, Bahram Mobasher, Ripon Saha, and Xianzhong Zheng for their helpful suggestions. KBM, DHM, LDL, and RE acknowledge funding from the NASA Hubble Space Telescope Archival Research grant 15040. CPC, RE, LBF, LDL, and SET acknowledge support from the Missouri Consortium of NASA's National Space Grant College and Fellowship Program. KBM also acknowledges funding from the School of Graduate Studies (SGS) fellowship grant and the Ronald A. MacQuarrie graduate fellowship offered by the University of Missouri-Kansas City.  DC acknowledges funding by the ERC Advanced Grant, STARLIGHT: Formation of the First Stars (project number 339177). This work is based on observations taken by the CANDELS Multi-Cycle Treasury Program with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for Program number HST-GO-12060 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. The VELA simulations were performed at the National Energy Research Scientific Computing Center (NERSC) at Lawrence Berkeley National Laboratory, and at NASA Advanced Supercomputing (NAS) at NASA Ames Research Center. This publication also made use of NASA's Astrophysics Data System Bibliographic Services, TOPCAT <cit.>, the core python package for the astronomy community <cit.>, collage maker (<github.com/delimitry/collage_maker>).
	
	
	
	

§ AFFILIATIONS

	^1Department of Physics and Astronomy, University of Missouri-Kansas City, Kansas City, MO 64110, USA.

	^2Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA.

	^3Department of Physics and Astronomy, University of Missouri, Columbia, MO 65211-7010, USA.

	^4Department of Physics and Astronomy, Colby College, Waterville, ME 04961, USA.

	^5Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA.

	^6Department of Astronomy, University of Michigan, Ann Arbor, MI, 48109-1107, USA.

	^7Universität Heidelberg, Zentrum fur Astronomie, Institut für Theoretische Astrophysik, 69120, Heidelberg, Germany.

	^8Physics Department, University of California, Santa Cruz, CA, 95064, USA.

	^9Santa Cruz Institute for Particle Physics, University of California, Santa Cruz, CA 95064, USA.
	mnras 
	
