Accurate Identification of Galaxy Mergers with Stellar Kinematics

As in N19, we utilize GADGET-3/SUNRISE simulations of merging galaxies. GADGET-3 (Springel & Hernquist 2003; Springel 2005) is a smoothed particle hydrodynamics (SPH) and N-body code, modeling processes such as radiative heating, radiative cooling, star formation, supernova feedback, and the multi-phase interstellar medium (ISM), using sub-resolution prescriptions. GADGET-3 also includes SMBH accretion, as well as AGN feedback (this is achieved by coupling 5% of the accreted luminosity to the gas as thermal energy). GADGET has been used for many different astrophysical applications, including widespread use in studies of merging galaxies (e.g., Di Matteo et al. 2005; Cox et al. 2006; Hopkins et al. 2006, 2008, 2013a; Robertson et al. 2006; Blecha et al. 2011, 2013b; Hopkins et al. 2013b).

We present the five galaxy merger simulations and their matched isolated simulations in Table 1. The framework for these simulations is established in Blecha et al. (2018), and the simulations themselves are presented in N19. Three of the simulations are major mergers (where the mass ratio, q, of the progenitors is greater than q = 0.25, or 1:4; Rodriguez-Gomez et al. 2015; Nevin et al. 2019) and two of the simulations are minor mergers. The major mergers have mass ratios of 1:2, 1:3, and 1:3. We define the gas fraction of these simulations as f_gas = M_gas,disk/(M_gas,disk + M_*,disk). The 1:2 and 1:3 mass ratio major mergers have a relatively high gas fraction of 0.3, and one of the 1:3 mass ratio major mergers has a relatively low gas fraction of 0.1.

We verify that the different gas fractions of the simulations (0.1 and 0.3) cover the full range of the gas fractions in the MaNGA galaxies. The mean gas fraction in MaNGA is defined by Barrera-Ballesteros et al. (2018) as

where σ_gas is the gas mass density, and σ_* is the stellar mass density. The stellar mass density is derived from the stellar template fitting and the gas mass density is derived from CO measurements in the EDGE-CALIFA survey (Bolatto et al. 2017; Barrera-Ballesteros et al. 2020).

Barrera-Ballesteros et al. (2018) find that the mean gas fraction for the MaNGA sample ranges from 0.16 to 0.32. A f_gas value of 0.1 is therefore below the mean for MaNGA galaxies, and as such is relatively gas poor, while f_gas = 0.3 is at the top of the range of mean values for the sample, and relatively gas-rich.

These simulations are named for their mass ratio and gas fraction; for instance, the gas-rich 1:2 mass ratio major merger is q0.5_fg0.3, the gas-rich 1:3 mass ratio major merger is q0.333_fg0.3, and the gas-poor 1:3 mass ratio major merger is q0.333_fg0.1. All of the major merger progenitors have a bulge-to-total (B/T) mass ratio of 0, meaning that they are a pure disk initially. Both of the progenitor galaxies of the minor mergers have a B/T ratio of 0.2, and both are gas-rich. These simulations are q0.2_fg0.3_BT0.2, which is the 1:5 mass ratio minor merger, and q0.1_fg0.3_BT0.2, the 1:10 mass ratio minor merger.

We build the matched sample of isolated galaxy simulations for each merger simulation based on two sources. Firstly, we use a standalone sample of isolated galaxies, matched for mass and gas fraction to each of the simulations. Some simulations have more than one matched isolated galaxy, but for a case where there is only one isolated galaxy, it is matched to the mass of the larger merging galaxy from the corresponding merger simulation. Secondly, we define snapshots of each simulated merger falling before the first pericentric passage, or >0.5 Gyr after final coalescence, as isolated galaxies. We refer to those isolated galaxies from snapshots occurring prior to the first pericentric passage as “pre-merger” isolated galaxies, and those snapshots occurring > 0.5 Gyr after final coalescence as “post-merger” isolated galaxies. This distinction is useful, as the properties of these two populations differ.

We couple GADGET-3 with SUNRISE in order to directly compare the simulated galaxies with observations. SUNRISE is a 3D polychromatic Monte Carlo dust radiative transfer (RT) code (Jonsson 2006; Jonsson et al. 2010), which is used to model a wide range of isolated and merging galaxies (e.g., Narayanan et al. 2010; Hayward et al. 2011, 2014; Blecha et al. 2013a; Snyder et al. 2013). The full details of the SUNRISE prescription are presented in Blecha et al. (2013b, 2018) and N19. Briefly, SUNRISE performs Monte Carlo radiative transfer on a 3D adaptively refined grid to compute the emission from stars, H II regions, and AGN. SUNRISE uses the STARBURST99 stellar population synthesis models (Leitherer et al. 1999) to calculate the age- and metallicity-dependent spectral energy distributions for each star particle. The treatment for dust includes dust self-absorption and thermal re-emission, as well as polycyclic aromatic hydrocarbon absorption and emission. In addition, we include kinematic (Doppler) effects, which requires very high spectral resolution. Ultimately, SUNRISE calculates the emergent, attenuated, resolved UV-to-IR spectra (3300–6990 Å, Δλ = 0.3 Å) for seven isotropically positioned viewing angles.

We utilize a data cube of SUNRISE optical synthetic spectra for the seven isotropically positioned viewpoints from each merger snapshot to produce the mock data cubes. In N19, a “snapshot” is the SUNRISE image; in this work, we use the term “snapshot” to refer to the full data cube from a specific point in time. These snapshots occur at 50–100 Myr intervals during each merger, and we refer to them as early-stage, late-stage, and post-coalescence stage snapshots. We define these stages using the r-band images from N19. The early-stage mergers occur after the first pericentric passage, and have (viewpoint) average stellar bulge separations of Δx ≥ 10 kpc, late-stage mergers have separations of 1 kpc < Δx < 10 kpc, and post-coalescence mergers are no longer resolvable with separations Δx ≤ 1 kpc until 0.5 Gyr after final coalescence. With a 50–100 Myr cadence for snapshots, we therefore have 5–10 snapshots for each of these stages. In total, there are ∼20 snapshots per simulation, and seven viewpoints per snapshot, which amounts to 100–200 observations per merger simulation.

We further discuss the importance of running RT, and incorporating dust attenuation and scattering into the merger snapshots in Appendix A; briefly, the stellar kinematic maps are affected by both the presence of dust, and dust scattering. The implication is that for this type of kinematic analysis, it is important to use velocities derived directly from the RT product (synthetic spectra), as opposed to the original SPH particle velocities.

To produce stellar kinematics for our sample of simulated galaxies, we use the specifications of MaNGA to create a data cube of spectra, and then we mimic the MaNGA Data Analysis Pipeline (DAP) to extract stellar kinematics for use in our kinematic classification. Examples of finalized “MaNGA-ized” stellar velocity, and stellar velocity dispersion maps, are presented in Figure 1. In this section, we describe how we mimic the specifications of MaNGA. This involves reducing the spatial and spectral resolution of the simulations to create a MaNGA-ized data cube, placing an appropriately sized fiber bundle over each galaxy, and fitting each spaxel with ppxf (a penalized pixel fitting method from Cappellari & Emsellem 2004; Cappellari 2017) to obtain the velocity and velocity dispersion of the stars at each spatial position.

Zoom In Zoom Out Reset image size

SDSS-IV/MaNGA is an IFS survey, targeting a sample of >10,000 nearby galaxies, and selected to span a wide range of environments and stellar masses (Gunn et al. 2006; Smee et al. 2013; Bundy et al. 2015; Drory et al. 2015; Law et al. 2015; Blanton et al. 2017). Spectra are obtained using the BOSS spectrograph on the 2.5 m telescope at Apache Point Observatory (Gunn et al. 2006), which has a spectral resolution of R ∼ 2000. Fibers are bundled into integral field units (IFUs); MaNGA has five different fiber bundles, equipped with 19, 37, 61, 91, and 127 fibers (the largest fiber bundle is known as the “frankenbundle”); each individual fiber has a 20 diameter with 25 spacing between fibers (Smee et al. 2013; Drory et al. 2015; Yan et al. 2016a; Wake et al. 2017). These fiber bundles range from 125 to 325 in diameter.

MaNGA has a median point-spread function (PSF) of 25, which roughly corresponds to a spatial resolution of 1–2 kpc. The primary sample of galaxies (which is 2/3 of the full sample) has coverage out to 1.5 times the r-band effective radius (R_e ), and the secondary sample has coverage to 2.5 R_e (Yan et al. 2016a; Wake et al. 2017). The redshift range of the MaNGA survey is 0.01 ≤ z ≤ 0.15.

The most recent internal MaNGA product launch (MPL-9) includes 8000 unique galaxies, observed and reduced by the Data Reduction Pipeline (DRP; Law et al. 2016). The publicly available version is released as DR-15, which includes 4621 unique galaxies (Aguado et al. 2019). The derived properties (including the stellar kinematics) are produced by the DAP (Belfiore et al. 2019; Westfall et al. 2019) in the format of a single data cube per galaxy (Yan et al. 2016b). The MaNGA team also creates and maintains Marvin, which is a useful tool for visualizing MaNGA data (Cherinka et al. 2019).

To create the mock data cubes, we begin with the SUNRISE synthetic spectra, which we extract at the median redshift of the MaNGA survey (z = 0.03). We also use the SUNRISE SDSS r- and g-band images to construct the mock data cubes, as these are essential for certain steps in the process. We follow this procedure (which mirrors the MaNGA DAP whenever possible):

1. 1.  
   Convolve the data cube with the 25 MaNGA PSF.Here we use a 25 Gaussian kernel, which is a good approximation of the effective PSF of the MaNGA data cubes. A model of the effective PSF is automatically computed for each data cube as part of the MaNGA DRP, which convolves a simulated point source with the fiber footprint of a given set of observations, incorporating as-observed details of seeing, transparency, differential atmospheric refraction, dithering, and other instrumental effects (Law et al. 2016). We briefly investigate the difference between using our simplified Gaussian kernel and the effective PSF model provided by the DRP, and find that there are small differences in the maps when the reconstructed PSF is used, in particular a slight increase in the spread of values in the velocity dispersion maps. This is to be expected, given that the reconstructed PSF does not have a perfectly Gaussian shape. However, the differences in the kinematic maps are minimal, and do not cause significant differences in classification. We define the classification “significantly different” in Section 3.8. We therefore conclude that using the 25 Gaussian PSF is adequate for this work.
2. 2.  
   Rebin to match the spatial (05 spaxels) and spectral sampling (R ∼ 2000) of MaNGA. The spectral sampling varies as a function of wavelength.
3. 3.  
   Use a mock g-band image, rebinned to the 05 spatial scale, to mask all spaxels in the data cube that fall below a g-band signal-to-noise ratio (S/N) cutoff value of 1.Follow the procedure from N19 to convolve, rebin, and introduce noise characteristic of SDSS imaging to the mock g-band images, to match the 05 spatial binning of the MaNGA cubes. Next, find the average g-band S/N per spaxel, and mask all spaxels that fall below an S/N cutoff value of 1. This procedure directly follows the MaNGA DAP (Westfall et al. 2019), which masks all spaxels using the same g-band S/N cutoff.
4. 4.  
   Use the MaNGA procedure to select which size of fiber bundle to use for each mock data cube, and mask all spaxels that are external to this hexagonal footprint.We use statmorph (Rodriguez-Gomez et al. 2019) to measure the effective radius of the mock r-band images from N19. We then determine the smallest fiber bundle needed to cover each galaxy to 1.5 R_e (this is how MaNGA’s primary sample is defined). We select the smallest fiber bundle if the total angular extent of the galaxy (2 × 1.5 R_e ) is smaller than 125, and the largest fiber bundle if the angular size exceeds 325.
5. 5.  
   We introduce noise to each spaxel to produce a data cube with noise, and a sqrt(variance) data cube (from here on, “error data cube”).We first produce a typical noise spectrum to demonstrate how the noise trends with wavelength for MaNGA observations. We then normalize this noise spectrum, using the g-band S/N value for each spaxel. The end result is a sqrt(variance), or error spectrum, which we use to introduce random noise to each spaxel in the data cube. The noisy spectra, and the accompanying error spectra, are the inputs to ppxf. Further details of this process can be found in Appendix B.To verify that the S/N of the simulated spectra are representative of the MaNGA sample, we use the peak g-band S/N as a comparison statistic, defined as the maximum value of the g-band S/N (per pixel) from a single galaxy observation. The peak g-band S/N ratio of a sample of MaNGA galaxies spanning the full range of sizes, surface brightnesses, and stellar masses of the MaNGA sample ranges from 10–60, with a median of 25. The same statistic for the simulation suite ranges from 10–100, with a median of 30. In Section 3.9, we experiment with changing the S/N of simulated spectra, and investigate how this affects the classification.Since the MaNGA data cubes oversample the effective PSF, they also contain significance covariance in the errors between adjacent spaxels, such that the S/N ratio of binned spectra does not increase as . This covariance is irrelevant for the fitting of individual spectra, but we account for it in our Voronoi binning by following the analytic approximation given by Law et al. (2016), as discussed below.
6. 6.  
   Having completed the masking steps, we further exclude regions that are background dominated.At this stage, we note that the data cubes have “patchy” outskirts, or regions of low S/N data, surrounded by masked regions. The MaNGA data cubes do not have this feature; instead, they exclude regions that could be characterized as “background dominated.” This patchiness does not affect the results of the classification; instead, we choose to correct it for cosmetic purposes. To do this, we mask spaxels where the g-band signal is less than 3σ above the background value, where σ is the standard deviation of the noise given above. This produces the desired effect, whereby the mask has a sharper cutoff, matching the appearance of the MaNGA cubes.
7. 7.  
   Rebin spatially, using a Voronoi binning scheme with a g-band S/N of 10 (Cappellari & Copin 2003).We create spatial bins with a g-band S/N of 10, reproducing the procedure described in Westfall et al. (2019). When a Voronoi bin contains more than one spaxel, the new spectrum is the masked average of all constituent spectra, while the error spectrum for that bin is determined by co-adding the error spectra. It is important to account for covariance between neighboring spaxels in our Voronoi bin calculation. In order to avoid the computational cost of calculating the covariance matrix for all spaxels, we instead use the correction given in Law et al. (2016). This correction is an analytic function of the number of spaxels in a bin (N_bins):

   

   where n_measured is the corrected noise level after the correction is applied to the co-added error, where covariance is not considered (n_{no covar}), and N_bins is the number of spaxels in a bin.

The final step in the creation of mock kinematic maps is to pass the Voronoi binned spectra through ppxf (Cappellari & Emsellem 2004; Cappellari 2017). ppxf is a penalized pixel fitting method, which assumes that a galaxy spectrum is a combination of stellar templates, convolved with line-of-sight velocity distribution (LOSVD).

To run ppxf, we follow these steps from the DAP:

1. 1.  
   Normalize the flux data so that the mean over all templates is unity.
2. 2.  
   Mask the spectra to match the wavelength range of the MILES-HC library (3600–7400 Å).
3. 3.  
   Mask the emission lines using the DAP module StellarContinuumBitMask().
4. 4.  
   Use the 42 template MILES-HC spectral library to globally fit each data cube.These templates are first convolved to the spectral resolution of MaNGA. ¹³
5. 5.  
   We use the “NZT,” or non-zero template iteration mode to fit all bins with ppxf.In this mode (which is also used in the DAP), we first fit the masked average of all spectra in the data cube, and use this global fit to isolate the subset of templates that are allocated non-zero weights. This template subset is then used to individually fit each bin.
6. 6.  
   Each fit iteration of ppxf uses an additive eight-order Legendre polynomial, and a Gaussian line of sight velocity dispersion (LOSVD) with two moments. As in the DAP, due to limited spectral resolution, we do not solve for the higher order moments, h₃ and h₄ (Westfall et al. 2019).

The final product of our MaNGA-izing procedure is the first two moments of the LOSVD, or a stellar velocity map and a stellar dispersion map, both with associated error maps from the fit to the stellar continuum.

Here, we define and describe the predictors extracted from the stellar kinematic maps. The goal is to create a set of kinematic predictors to adequately describe the different types of merger-induced kinematics in the velocity and velocity dispersion maps.

To develop this kinematic identification tool, we use stellar kinematics, rather than warm ionized gas kinematics (henceforth, “gas kinematics”). The stellar and gas kinematics trace different physical regions and processes in the merging galaxies. We select the stellar kinematics because they directly trace the assembly history of a galaxy’s stellar population. On the other hand, gas kinematics can be subject to a number of non-gravitational forces. Stellar and gas kinematics diverge in the presence of shocks, inflows, and/or outflows, all of which are processes not limited to merging galaxies. An analysis built on gas kinematics is a compelling direction for future work, but falls beyond the scope of this paper (i.e., see Khim et al. 2020 or Feng et al. 2020). For an analysis comparing the gas kinematics and the stellar kinematics of non-merging and merging galaxies in CALIFA, see Barrera-Ballesteros et al. (2014) and Barrera-Ballesteros et al. (2015), respectively. ¹⁴

The kinematic predictors are based on previous work to identify merging galaxies from the stellar kinematics of observed and simulated galaxies. All of these predictors are sensitive to the different orientations, merger stages, mass ratios, and/or gas fractions of merging galaxies. Our goal is to combine them into one LDA classification in order to best identify a variety of different types and epochs of merging galaxies. In total, we extract the following predictors (which are all introduced in Table 2): A, A₂, ΔPA, v_asym, σ_asym, resid, , , Δx_V , Δx_σ , μ_1,V , μ_1,σ , μ_2,V , μ_2,σ , ∣μ_3,V ∣, ∣μ_3,σ ∣, μ_4,V , and μ_4,σ . We include a brief definition of all predictors in Table 2, but focus the remainder of this section on those kinematic predictors selected by the random forest term selection technique described in Section 3.3: A₂, ΔPA, resid, , μ_1,V , μ_1,σ , μ_2,V , μ_2,σ , ∣μ_3,V ∣, ∣μ_3,σ ∣, μ_4,V , and μ_4,σ . These terms are the most informative for identifying the merging galaxies, and we discuss them throughout the remainder of this paper. We further describe those kinematic predictors that were not selected in Appendix D.

Table 2. Synthesis of All Kinematic Predictors Measured in This Paper

Predictor Name	Description	Derivation
A	The weighted asymmetry in the position angle,
	calculated from the Radon profile	where is the best-fit kinematic position angle
A2	The error-weighted asymmetry in the position angle
ΔPA	The difference between the global kinematic and	ΔPA = ∣PAkin − PAimg∣
	photometric position angles, measured
	based on kinemetry and the g-band image
σasym	Describes the degree of smoothness of the
	velocity dispersion map	which is the sum of the higher
		order coefficients from kinemetry
vasym	The deviation of the velocity dispersion map
	from ordered rotation	as above, but excluding the k1,v term
resid	The residual between the best-fit kinemetry
	model and the velocity map
	The approximate spin parameter
	Galaxy ellipticity	Measured using statmorph from the r-band imaging
ΔxV	The spatial distance between the center of the	The imaging center is measured from the r-band image, and
	velocity map and the imaging center in kpc	the kinematic center is based on the Radon Transform
Δxσ	Same as above, but for the velocity dispersion map	The center of the velocity dispersion map is determined
		using a low-pass filter
μ1,V and μ1,σ	The mean distribution of the	The distributions for each map are created
	velocity and velocity dispersion maps	by collecting the values of all spaxels
μ2,V and μ2,σ	The variance of the distributions
∣μ3,V ∣ and ∣μ3,σ ∣	The skewness of the distributions
μ4,V and μ4,σ	The kurtosis of the distributions

Note. We bold those predictors designated as important. We include a brief description and derivation for each predictor. For further details, see Section 2.3 for those predictors designated as important, and Appendix D for predictors not used in the classification.

Download table as:  ASCIITypeset image

To define the asymmetry in the kinematic position angle (A₂), we utilize the Radon Transform from Stark et al. (2018). We transform the velocity maps into circular coordinates (ρ,θ), where ρ is the distance from the spaxel to the center of the velocity map, which is the kinematic center (defined below), and θ is the angle between the positive x-axis and the line segment from the kinematic center to the spaxel. The angle θ ranges from 0 to 180 in the CCW direction. Positive values of ρ are regions of the velocity map above the x-axis, and negative values of ρ are below the positive x-axis.

The Radon Transform is defined as

where the velocity is summed along line integrals centered on the point (ρ, θ), and perpendicular to the kinematic center of the galaxy. The Radon Transform is a 2D array, calculated at all values of ρ and θ.

We then calculate the bounded Absolute Radon Transform, R_AB , which is integrated over a distance R_e , and denotes the absolute value of the difference between the velocity at each point, and the mean value along the line segment.

We present the bounded Absolute Radon Transform and the Radon profile in Figure 2. The Radon profile is computed by determining the minimum value of θ (, where the hat operator denotes an estimated value) for each value of ρ from the bounded Absolute Radon Transform. The value of traces the direction of maximal rotation in the stellar velocity maps at each radial position.

Zoom In Zoom Out Reset image size

We follow the procedure of Stark et al. (2018) to determine the galaxy’s kinematic center, which we describe in more detail in Appendix D.

We quantify the asymmetry of the Radon profiles using the kinematic predictor A₂, defined in Stark et al. (2018):

where is the absolute magnitude of the difference between the values of θ_i on one side of the Radon profile and the other (same ρ, different sign), is the uncertainty on , and the expression is summed over the i values of .

The A₂ predictor incorporates the absolute magnitude of the difference between the measured kinematic PA on one side of the galaxy from the other. We therefore expect that A₂ will be enhanced for merging galaxies, since mergers can cause warps in the stars of a galaxy (e.g., Shapiro et al. 2008).

We use kinemetry to measure both ΔPA and residuals (resid) from the LOSVD (Krajnovic et al. 2006). Functionally, kinemetry measures the kinematic asymmetry from the line of sight velocity maps by dividing them into a set of nested elliptical rings. The best-fit model at each radius is determined using a ring defined by the kinematic PA, and the flattening factor q_f = 1-e, where e is the ellipticity of the ring in the plane of the sky. These models use a decomposition of the moment maps into harmonic Fourier coefficients in polar coordinates. For instance, a velocity map, K(r, ψ), can be expanded into a finite number of harmonic frequencies:

where r is the semimajor axis of the ellipse, ψ is the azimuthal angle, A₀(r) is the systemic velocity, N is the number of ellipses fit, and A_n and B_n are the coefficients of the harmonic expansion. The best-fitting ellipses are obtained by minimizing χ² for the linear combination of the A_n and B_n coefficients.

An ideal rotating disk can be described using only the B₁ term, which represents the cosine term for the circular velocity of a galaxy’s rotating disk:

where r is the radius in the plane of the galaxy, ψ is the azimuthal angle, V_c (r) is the circular velocity, and i is the inclination of the galaxy disk.

To determine the best-fit Fourier coefficients, we run kinemetry multiple times. We first allow the best-fit kinematic PA and value of q_f to vary for each radius. We define the kinematic position angle (PA_kin) to be the median value of the best-fit kinematic PAs. We then allow the value of q_f to vary and determine the median value. After determining the global values for kinematic PA and q_f , we do a final run to determine the values of the higher order kinematic moments and therefore the best-fit disk model. We then compare PA_kin to the imaging major axis (PA_img, measured using statmorph from the r-band imaging) to create the predictor, ΔPA. Since ΔPA traces the recent global misalignments of stars, it should be elevated for those merging galaxies with misaligned stellar disks.

We use the global kinematic position angle from kinemetry to measure ΔPA, rather than the median of the kinematic position angles from the Radon Transform. The main motivation for this choice is that kinemetry uses an adaptive binning scheme; at each step outwards, the ellipses are larger, which gives less weight to the kinematic confusion at the outskirts of the galaxy. The Radon Transform, however, is equally sampled in ρ (see Figure 2), and can therefore be more influenced by the measurement of the kinematic PA at the outskirts of the galaxy. In most cases, the two measurements agree within error, but in cases where the kinematic maps are disturbed, the global kinematic PA from kinemetry more closely matches our by-eye assessment of the kinematic PA.

We also extract “resid,” or the kinemetry residuals between the best-fit rotating disk and the velocity map. This predictor is defined as

where V_* is the observed velocity map, V_model(r, Ψ) the circular velocity model from kinemetry, and N is the number of fitted spaxels. We include this normalization factor in order to penalize those fits that converge to a very inclined galaxy. For these galaxies, the fit is attempting to avoid fitting disordered kinematics in the exterior regions of the galaxy by fitting a smaller region. We show an example of a simulated galaxy snapshot from the q0.5_fg0.3 simulation fit with kinemetry and its velocity residuals in Figure 3.

Zoom In Zoom Out Reset image size

We measure , the approximate spin parameter, from the stellar velocity and velocity dispersion maps, as defined by Emsellem et al. (2007):

where F_n is the (r-band) flux of a spaxel, R_n is the distance from the kinematic center, V_n is the stellar velocity, and σ_n is the stellar velocity dispersion. We measure to the r-band effective radius. Since the fiber bundles are designed to provide coverage of each galaxy to 1.5R_e , if a secondary nuclei falls toward the outside edge of the hexagonal FOV, it is excluded from the measurement of . This effect is more relevant to minor mergers, where the secondary component covers a smaller effective area of the hexagonal FOV.

We measure the ellipticity of a galaxy, , based on r-band photometry, using statmorph. This is distinct from the ellipticity parameter used by kinemetry to fit rotation curves. We do not use as a kinematic predictor. Instead, we use it to construct the – diagnostic diagram in Section 4.1, where the division between fast and slow rotators is as defined by Cappellari (2016):

where slow rotators fall below this line.

In the – diagram, is the more predictive of the two axes; it decreases dramatically for the “slow-rotating” population of galaxies, which are dynamically disordered and dispersion-dominated. We predict that will decrease for merging galaxies, since mergers are kinematically disordered, and can contribute to bulge growth, which is associated with enhanced velocity dispersion.

In addition to the kinematic predictors utilized in previous works, we define a new set of predictors, based on the distributions of values in the velocity and velocity dispersion maps. These predictors include μ_1,V /μ_1,σ , μ_2,V /μ_2,σ , ∣μ_3,V ∣/∣μ_3,σ ∣, and μ_4,V /μ_4,σ , which are the standardized moments of the stellar velocity/velocity dispersion maps. These predictors are similar to the formulation given by Sweet et al. (2020), which calculates the moments of PDF(s), where s is the normalized specific angular momentum.

To determine the values of these predictors, we measure the four standardized moments of the distribution; mean (μ₁), variance (standard deviation; μ₂), skewness (μ₃), and kurtosis (μ₄). This produces eight different predictors (four each from the velocity and velocity dispersion distributions). These quantities are different from the higher order moments, h₃ and h₄, which are typically measured by ppxf. We show an example of these predictors, measured from a velocity dispersion map, in Figure 4.

Zoom In Zoom Out Reset image size

We expect to see an offset in the mean velocity (μ_1,V ) from the systemic for merging systems, and an enhanced mean velocity dispersion (μ_1,σ ). The spread in the velocity distribution (μ_2,V ) and the dispersion of the velocity dispersion distribution (μ_2,σ ) could identify superpositions of dynamically distinct stellar components. This could include a secondary merging galaxy, or features such as a stellar bulge.

The higher order moments could be useful for identifying subtler features of mergers, beyond bulk shifts in μ_1,V , for example. The skewness of a distribution is sensitive to the tails; we take the absolute value to treat positive and negative skew identically. A skewed velocity or velocity dispersion distribution (∣μ_3,V ∣ and ∣μ_3,σ ∣) could have a faint secondary source in the field of view, where the distribution is actually a combination of two galaxy rotation curves. Kurtosis measures how peaked a distribution is, relative to normal; a flatter distribution has a negative kurtosis, and a more peaked distribution has a higher peak. A smoothly rotating velocity dispersion field has a normally shaped distribution, whereas a disturbed field may have a negative (flatter) kurtosis (μ_4,σ ). On the other hand, post-coalescence mergers with recent bulge growth could have a positive kurtosis in the velocity dispersion distribution.

To summarize, we extract the following kinematic predictors: A, A₂, ΔPA, v_asym, σ_asym, resid, , , Δx_V , Δx_σ , μ_1,V , μ_1,σ , μ_2,V , μ_2,σ , ∣μ_3,V ∣, ∣μ_3,σ ∣, μ_4,V , and μ_4,σ . We then use the techniques described in the following sections (Sections 3.1–3.3) to select the most informative of these predictors. We ultimately use the following predictors in the LDA classification: A₂, ΔPA, resid, , μ_1,V , μ_1,σ , μ_2,V , μ_2,σ , ∣μ_3,V ∣, ∣μ_3,σ ∣, μ_4,V , and μ_4,σ .

The classification in this work relies on an LDA technique that separates nonmerging galaxies from merging galaxies, based on a combination of input predictors (for a review of LDA, see James et al. 2013). This approach was first presented in N19 for imaging predictors; here, we use this approach for kinematic predictors.

LDA is one of many statistical learning tools for performing classification tasks. Using pre-defined features (predictors) as inputs, LDA solves for the hyperplane in multi-dimensional predictor space that maximizes the separation between different classes of objects (i.e., mergers and nonmergers). The solution is a linear combination of the input predictors; the classification is therefore relatively easy to interpret, since its complexity is low. Recent works have employed other techniques to identify merging galaxies, such as random forest regressors (e.g., Goulding et al. 2018; Snyder et al. 2019), and convolution neural networks (CNNs; e.g., Bottrell et al. 2019). These techniques have various advantages and disadvantages, based on the data set at hand, and the goals of the work. Since we aim to optimize the interpretability of the method, we select LDA over, e.g., a CNN. CNNs may increase the number of correct classifications, but they achieve this using complex nonlinear features, which are not easily interpreted.

In this work, we have made several important changes to the technique employed in N19. As such, we review the relevant details from the LDA in N19, prior to discussing these changes. Relevant details of the LDA technique from N19 include:

1. 1.  
   All predictors are linearly standardized prior to the LDA technique, meaning that predictors with large numerical values (such as A₂) do not have an outsized effect on the analysis.
2. 2.  
   We utilize priors on the relative fraction of merging and nonmerging galaxies in nature versus that in the simulations. This accounts for the presence of a greater number of merging galaxy snapshots (relative to nonmerging snapshots) for each simulation. We use the same priors as given in N19; f_merg = 0.1 for major mergers, and f_merg = 0.3 for minor mergers. These priors are based on the fraction of nonmerging and merging galaxies, taken from both observation and simulations (e.g., Conselice et al. 2009; López-Sanjuan et al. 2009; Shi et al. 2009; Lotz et al. 2011; Rodriguez-Gomez et al. 2015).
3. 3.  
   We include interaction terms to explore correlations between predictors.
4. 4.  
   In order to select which coefficients are necessary for the classification, we use a forward stepwise selection technique, which orders and includes only the most important terms. This technique adds additional terms to LD1 only if they improve the F1 statistic, which is defined in Section 3.5. It also protects against the unnecessary addition of terms by finding the minimum number of terms that produce an F1 statistic, consistent with the maximum (within 1σ errors, where σ is the standard deviation on the F1 statistic, measured from each k-fold cross-validation set).
5. 5.  
   We use k-fold cross-validation to obtain 1σ errors on the predictor coefficients. At each step of the forward stepwise selection process, we divide the sample into k subsets. We then train the LDA on the first k − 1 subsamples, and test on the remaining subsample, which functions as the “cross-validation” set. We repeat these steps k times for all combinations of subsamples, and the variation in predictor coefficient values from the cross-validation subsamples is the 1σ error.

For complete details, including the full mathematical formulation for LDA, see N19.

We make several changes to the above technique, motivated by the additional challenges arising from the kinematic data:

1. 1.  
   Due to the number of predictor terms in this work, instead of including all of the kinematic predictors in the final classification, we first utilize the RFR technique as a selection technique to eliminate uninformative predictors from our analysis (see Section 3.3).
2. 2.  
   We adjust the model optimization statistic. In N19, we minimized the number of misclassifications, in order to both select predictors, and determine their coefficients. Here, we utilize the F1 statistic defined in Section 3.5 instead; this does a better job of balancing the number of false negatives and false positives in each classification.
3. 3.  
   We also adjust the k-fold cross-validation, Rather than using k = 10, we find that a value of k = 5 improves the performance of the LDA, creating a training set, constituting 80% of the sample, and a cross-validation set, constituting 20% of the sample (as opposed to the 90%/10% divide in N19). We then train the LDA on nine of the subsamples, and test on the tenth sample. We repeat this procedure ten times, and the mean number of misclassifications across all ten test samples allows us to determine which set of input predictors to select

We use the LDA both as a term selector, and to determine the coefficients and standard errors for each selected predictor. In order to directly compare imaging classification to kinematic classification, we utilize the same snapshots from all simulations, and we rerun the imaging analysis using all the same methods as for the kinematic classification.

Prior to term selection and classification, we examine the distributions of predictor values. We screen for outliers and examine the data in the context of the assumptions made by the LDA. The goal is to gain an understanding of the properties of the data, prior to classification.

Firstly, we remove outliers by transforming the distribution of each predictor into log space. We define data points falling more than 5σ above or below the mean of the distribution for each predictor as outliers. This combination of log transformation and 5σ cutoff allows us to identify outliers caused by errors in the creation of the mock maps, and not simply related to very disturbed kinematics. There are ∼4 outliers per simulation, amounting to 4 out of 100 or 200 data points.

Secondly, we check the input data for significant violations of the LDA’s assumptions. LDA operates under the assumptions that the predictors are normally distributed (multivariate normality), the covariance among the merging and nonmerging classes is equal (homoscedasticity), and the predictors are not strongly correlated with one another (multicollinearity).

Here, we test these three assumptions by closely examining the data. We carry out the same statistical tests described in N19 to test for normality, homoscedasticity, and multicollinearity. We find that the data violates all three assumptions. In addition, we plan to introduce interaction terms into the LDA classification, which will further increase the multicollinearity.

As discussed in N19, LDA has proven to be robust to violations of multivariate normality and homoscedasticity (Duda et al. 2001; Li et al. 2006). To ensure that the LDA technique used in this work is robust to these violations, we directly compare the LDA results to those of a logistic regression. The logistic regression and the LDA produce similar results, both in terms of the relative importance of the predictors for each simulation, and in the performance of the method. This is an indication that the LDA is converging, even though it nominally violates several assumptions.

In N19, we used a forward stepwise selection technique within the LDA to select informative predictors. Here, motivated by both an increase in the number of initial terms, ¹⁵ and a decrease in the predictive power of these terms, we modify the term selection procedure. We introduce a random forest regressor (RFR) into the methodology to select a subset of predictors for each simulation, which will then be presented to the LDA classifier.

An RFR (Ho 1995) is an ensemble learning technique that aggregates the result of many individual decision trees, run in parallel. We specifically utilize the scikit-learn implementation of RFR (Pedregosa et al. 2011). In an RFR, the number of features that can be used to split at each node of the decision tree is limited to a percentage of the total number of features, ensuring that the ensemble model does not rely too heavily on any one feature. This means that the RFR is able to combine all potentially predictive variables in a fair way. It is also able to incorporate nonlinear features to capture some higher order interaction terms. In practice, we find that the RFR is an efficient method of initially identifying useful features in the data set from an extensive list of kinematic predictors. ¹⁶

In order to select the informative terms from the RFR, we include an additional predictor. This predictor is assigned a random number for each galaxy snapshot, and therefore shows no significant difference between nonmerging and merging galaxies. We use this technique to eliminate all terms with a feature importance lower than the random term for all simulations. In this step, we eliminate the v_asym, σ_asym, A, Δx_V , and Δx_σ predictors. Next, for each individual simulation, we additionally eliminate predictors with an importance lower than the random value prior to initiating the LDA classification. The terms eliminated in this step vary from simulation to simulation.

Having utilized RFR term selection to narrow the number of kinematic terms down to 11 (ΔPA, resid, , A₂, μ_1,V , μ_2,V , μ_2,σ , ∣μ_3,V ∣, ∣μ_3,σ ∣, μ_4,V , and μ_4,σ ), we then run the LDA classification for each individual simulation. We also combine the three major mergers into a combined major merger classification, and the two minor mergers into a combined minor merger classification. We run the LDA with interaction terms; the result is a linear combination of selected predictors and coefficients which is unique for each simulation. We present the term coefficients and standard errors for the four most important terms and the intercept term in Table 3. Finally, we briefly discuss the main results of the LDA classification for each simulation, which we will examine in more detail in Section 4.

Table 3. Final LD1 Predictor Coefficients () with 1σ Confidence Intervals after Term Selection for the First Four Most Important Terms, and the Intersect () for All Simulations

Simulation
	1	2	3	4
All Major	−6.76 ± 0.45	4.99 ± 0.6 ∣μ3,σ ∣	4.54 ± 0.36 μ1,σ ∗	−4.44 ± 0.51 μ1,σ ∗ ∣μ3,σ ∣	−1.21 ± 0.07
All Minor	−4.99 ± 0.74 μ2,σ	−4.97 ± 0.59 μ2,σ ∗ μ4,V	3.47 ± 0.62 μ4,V ∗ μ4,σ	2.44 ± 0.38 μ4,V	−0.76 ± 0.04
q0.5_fg0.3	−7.15 ± 0.78 μ1,σ ∗ ∣μ3,σ ∣	−6.7 ± 0.63 μ1,σ ∗ μ2,σ	6.65 ± 0.53 μ2,σ	5.75 ± 0.2 μ1,σ	−2.57 ± 0.05
q0.333_fg0.3	8.27 ± 0.35 μ1,σ	−7.84 ± 0.71 μ1,σ ∗ μ2,σ	5.92 ± 0.52 μ2,σ	5.21 ± 0.73 ∣μ3,σ ∣	−0.77 ± 0.18
q0.333_fg0.1	−7.78 ± 0.91 μ1,σ ∗ μ2,σ	7.09 ± 0.59 μ2,σ	5.97 ± 0.61 μ1,σ	⋯	−0.26 ± 0.28
q0.2_fg0.3_BT0.2	−6.51 ± 1.09 μ1,σ	−6.2 ± 0.93 μ2,σ ∗	−5.75 ± 1.65 A2	5.5 ± 0.67 μ1,σ ∗	−0.79 ± 0.05
q0.1_fg0.3_BT0.2	25.06 ± 5.11 μ1,σ ∗ μ4,V	16.02 ± 3.19 μ1,σ	−12.88 ± 2.51 μ4,V	6.8 ± 0.99 μ4,σ	−1.06 ± 0.07

Download table as:  ASCIITypeset image

We first introduce the mechanics of the classification. LD1, which is the first linear discriminant axis, is formed from the linear combination of coefficients, multiplied by the standardized predictors and an intercept term:

where C is the matrix of coefficients, X is the standardized values of the selected predictors, and B is the intercept term.

LD1 is the hyperplane that best separates the populations of merging and nonmerging galaxies for each simulation. We use the result from the major merger classification as an illustrative example of how to interpret the LDA results. The LD1 for the major merger combined run (truncated after seven terms) is

where LD1 is a linear combination of all selected terms, which comprise a coefficient (positive or negative), followed by the standardized value of a predictor. The last term is the intercept term.

The higher the value of LD1, the more likely the galaxy is to be classified as merging. We calculate the values of LD1 for all snapshots from the major and minor combined simulations in Figure 5. The horizontal line at an LD1 value of zero is the decision boundary, corresponding to a p_merg value of 0.5; all galaxies with an LD1 value greater than zero would be classified as merging, using a threshold value of 0.5. ¹⁷ We find that the classification is better able to separate the merging and nonmerging classes for the major merger simulations; this is reflected in Figure 5.

Zoom In Zoom Out Reset image size

There are important nuances to the interpretation of the selected predictors and their coefficients in Equation (8), because the interaction terms complicate the analysis. For instance, in Equation (8), the first selected predictor is , which has a negative coefficient. Ignoring the rest of the equation, this means that if the value is large, then the probability that a given galaxy is merging will decrease. However, there are other terms in the equation that are coupled with other predictors in interaction terms. This means that tweaking the value of will not linearly change the value of LD1. While the interaction terms complicate the analysis, they are an integral part of the classification. Many of the most important terms for LD1 in Table 3 are interaction terms, and including them significantly improves the performance of the LDA. As we discuss in more detail in Section 4.4, these terms are able to capture the non-monotonic movement of mergers through predictor parameter space.

While it may be difficult to untangle many of the contributing terms, we can use Table 2 to determine which predictors are most prevalent, and therefore informative, for each simulation. For instance, the μ_2,σ and μ_1,σ predictors are selected as either primary or interaction terms for all simulations. They are therefore universally useful kinematic predictors (for a full discussion of why these terms are important, see Section 4.2.2). The selected predictors from the q0.333_fg0.1 simulation are similar to those of the q0.333_fg0.3 simulation. These two simulations are matched for mass ratio, but not for gas fraction. However, the difference between a gas fraction of 0.1 and 0.3 is insignificant, so we hesitate to form any conclusions regarding the impact of gas fraction on stellar kinematics. On the other hand, the minor merger simulations differ from the major merger simulations in the selected predictors. We find that and ∣μ_3,σ ∣ are important for the major mergers, whereas the minor mergers are more reliant on certain higher order terms, such as μ_4,V and μ_4,σ . We explore the implications of these findings, in terms of the physical nature of the kinematics of mergers, in Section 4.

Here, we define and measure the accuracy, precision, recall, and F1 statistic of the simulations (for a review, see Fawcett 2006). We present these results using a confusion matrix for the major and minor combined simulations in Figure 6, which shows the relative fraction of known mergers and nonmergers in the cross-validation samples, classified by the LDA as merging and nonmerging. These quantities are derived by taking the mean of the performance statistics measured on each of the cross-validation samples. We quantify the accuracy, precision, recall, and F1 score for all simulations in Table 4.

Zoom In Zoom Out Reset image size

The accuracy for a given simulation is defined as the number of correct classifications of mergers as mergers (true positives), and the number of correct classifications of nonmergers as nonmergers (true negatives), divided by the number of total classifications:

where FP is the number of false positives, or nonmerging galaxies classified as mergers, and FN is the number of false negatives, or mergers classified as nonmerging. A classifier has a higher accuracy when it is able to increase the number of true classifications relative to false classifications.

Precision is defined as the number of true positive classifications over the total number of positive classifications:

A precise classifier maximizes the fraction of true positive classifications relative to false positives. Precision is also known as the “positive predictive value.” In this work, we seek to eliminate false positives from the sample, i.e., nonmerging galaxies, incorrectly classified as mergers.

Recall is defined as the number of true positive classifications over the total number of known mergers:

A classifier with high recall is also known as “complete,” because it correctly identifies the majority of mergers as such.

Finally, we measure the F1 score, or the F1 statistic, which is the harmonic mean of recall and precision:

F1 ranges in value from 0 to 1, and is strongly penalized if either precision (P) or recall (R) is small. We maximize the F1 statistic within the LDA during cross-validation, in order to select the predictor terms used in the classification.

Figure 6 presents the number of true negatives, false positives, false negatives, and true positives (left to right, top to bottom) for the combined major and minor merger simulations. It also quantifies the accuracy, precision, recall, and F1 score. The major merger classification performs better, with accuracy/precision/recall values of 0.81/0.95/0.76, respectively, whereas the minor mergers have values of 0.69/0.87/0.51, respectively. The imbalance between precision and recall is due to the priors utilized in the classification (we use the priors from N19, where f_merg = 0.1 for the major mergers, and 0.3 for the minor mergers). We have designed the classification with these strong priors so that when it is applied to galaxy surveys (where there are fewer mergers), the classifier will be more balanced. As a result, the classifier produces more false negatives than false positives when tested on the training set.

We experiment with adjusting the performance statistics of the classification, a process also known as “hyperparameter tuning.” It is possible to increase the number of false positives while decreasing the number of false negatives, either by adjusting the decision boundary (i.e., the threshold of p_merg), or by changing the priors. This could be a direction to pursue in future work, if we find that we are no longer tolerant of false negatives, or if we wish to adjust the priors on f_merg. As a test, we adjust the priors so that f_merg = 0.5, finding that it produces a similar classification, lower precision, and higher recall for the major merger classifications, and results in slightly different selected predictors and higher performance statistics for the minor mergers. While the classification with adjusted priors performs better on the training and cross-validation data sets, we find that it is not a fair representation of the fraction of merging galaxies in nature, so we present the original classification, where f_merg = 0.1 and 0.3, in this work.

Overall, we find that the kinematic classifications generally score lower for all performance statistics than the imaging classifications. This is true for all viewpoints (see Section 3.8 for a discussion of the effect of viewpoint on kinematic classification). For instance, the accuracy/precision/recall/F1 values for the combined major merger run with the imaging predictors are 0.88/0.98/0.84/0.90. For the combined minor merger run with imaging predictors, the values are 0.80/0.89/0.72/0.80. Another result is that kinematic minor merger classifications generally score ∼10% lower across all statistics than major mergers. We discuss the implications of the performance of kinematic predictors in comparison to imaging predictors, together with the performance of major versus minor merger classifications, in Section 4.8.

The LDA observability timescale is defined as the sum of all consecutive snapshots where the viewpoint-averaged LD1 value for a given snapshot is greater than zero. We present the observability timescales for all simulations in Table 5, along with the total merger duration for each simulation, and the fraction observability, or the fraction of a merger’s duration observable via the LDA technique. We exclude the combined major and minor mergers from this table, since they are built from mergers that progress at different rates.

All of the simulations have a relatively long timescale of observability (2–6 Gyr). The exception is q0.5_fg0.3, where the observability timescale is 0.9 Gyr, due to a decline in the LD1 values for a handful of snapshots in the late stages of the merger. We present a visualization of how the mean values of LD1 change throughout the lifetime of each merger in Figure 7. Here, and throughout the remainder of this paper, we show the q0.5_fg0.3 and q0.2_fg0.3_BT0.2 simulations as examples of major and minor mergers, respectively. This shows the viewpoint-averaged value of LD1 for each snapshot, as well as the 1σ confidence interval on this value, and for the total range. We also plot the decision boundary for each simulation, which falls at an LD1 value of zero (horizontal line). The minor mergers do not fall significantly above this line; even though the viewpoint-averaged LD1 values for the q0.2_fg0.3_BT0.2 simulation fall above the decision boundary, they overlap the decision boundary to 1σ confidence at almost all points in time. This means that not all viewpoints are significantly above this boundary. On the other hand, the major merger simulations are significantly above this boundary for the majority of their duration. While this is not shown in Figure 7, the q0.5_fg0.3 simulation is an outlier with regard to LDA observability time.

Zoom In Zoom Out Reset image size

Figure 7 also demonstrates the sensitivity of the LDA classification to the merger stage. For instance, there are a number of false negatives at the early stages of mergers, where the galaxies are more disk-like. Another key finding is that the sensitivity of the technique gradually decays during the post-coalescence and post-merger stages. We discuss the implications of the different observability timescales, and variations with time, in more depth in Sections 4.5, 4.6, and 4.7.

Here, we summarize the factors most likely to lead to false classifications (false positives and false negatives) for the different simulations. Our goal is to identify the primary failure modes of the classification system, and assess whether or not it is making reasonable choices. In other words, we should be concerned if our by-eye classification disagrees with the majority of the false classifications.

We present a visual version of a confusion matrix for the q0.5_fg0.3 and q0.2_fg0.3_BT0.2 classifications in Figures 8 and 9, respectively. These simulations are representative of the results for major and minor mergers, respectively. We generate example velocity and velocity dispersion maps for each classification category (in rows, top to bottom: TP, TN, FN, FP) by combining the results of each iteration of the k-fold cross-validation, then randomly selecting example snapshots from each category.

Zoom In Zoom Out Reset image size

In Figure 8, following a by-eye examination, it makes sense that many of the false negatives and false positives are incorrectly classified in the q0.5_fg0.3 simulation. The false negatives (third row) exhibit orderly rotation with relatively low velocity dispersions. These resemble the pre-merger isolated population shown in the true negatives row. The majority of the false positives (fourth row) are post-merger snapshots with kinematic disturbances.

Incorrect classifications of major mergers can mainly be attributed to two factors. Firstly, the false negatives are due to a lack of disturbed features, meaning that it is difficult to correctly classify many of these snapshots as mergers. Despite these limitations, the classification correctly identifies the majority of early-stage snapshots as mergers, meaning that it is out-performing the by-eye assessment in many cases. Secondly, kinematic disturbances induced by the merger persist into the post-merger stages, producing a number of false positives. These kinematic features are very similar to the features in the post-coalescence stages, so it makes sense that these are commonly classified as false positives. As we discuss in more detail in Section 4.6, our definition of the “end” of the merger (the dividing line between post-coalescence and post-merger stages) is somewhat arbitrarily defined, and results in a number of misclassifications from these two stages.

In Figure 9 we find that it is more challenging to correctly classify the nonmerging and merging galaxies in the q0.2_fg0.3_BT0.2 simulation using a by-eye assessment. For instance, the false negatives are very visually similar to the true negatives in terms of their kinematic features. The same is true for the false positives, which are similar to the example true positives. The exception is a number of obvious merger snapshots (we show an example of one such snapshot in the upper middle pair of panels) where the kinematics are dramatically affected. However, these disturbances are short-lived, so that the majority of the merging snapshots appear as per the example in the upper right corner of the diagram. This diagram illustrates the crux of the problem for minor mergers. While the LDA is able to pick up on a number of subtle features (such as stellar bulge enhancements), it ultimately struggles with a number of limitations, related to the lack of identifiability of all stages of a merger as such. These challenges contribute to minor merger classification having lower performance statistics as compared to major merger classification.

Zoom In Zoom Out Reset image size

Overall, the LDA is not misclassifying obvious (by-eye) mergers or nonmergers. The lack of identifiability of mergers/nonmergers given their kinematic maps is therefore the largest challenge for this technique. Other work highlights this same challenge with kinematic predictors; Hung et al. (2016) find that a significant fraction of merging galaxy kinematics remain indistinguishable by-eye relative to the nonmerging kinematics. This indicates something fundamental about galaxy kinematics: namely, that we are not missing obvious features and instead that merging galaxies are often indistinguishable from nonmergers.

It is well known that many kinematic predictors, such as , are correlated with galaxy inclination (e.g., Cappellari et al. 2007; Emsellem et al. 2011; Harborne et al. 2019). In this section, we examine how the viewing angle, which is a proxy for inclination, affects kinematic predictors, and ultimately, the LDA.

As described in Section 2.1, there are seven isotropically distributed viewpoints (0–6) for each snapshot. Critically, the inclinations are not an exact match between standalone isolated galaxies and merging galaxies. For instance, viewpoints 0 and 4 are the most face-on viewpoints for the merging galaxies, but viewpoints 4, 5, and 6 are the most face-on for the standalone isolated galaxies.

We first explore how inclination affects the predictor in Figure 10, where increases as the galaxy inclination increases. For instance, viewing angles 0 and 4 are the most face-on, and also have the lowest values of . When the 1σ error bars are taken into consideration, the difference in values is marginally significant. This is fully consistent with the results from Emsellem et al. (2011), who predict that the measured and values of an axisymmetric rotating oblate spheroid vary, relative to the viewing angle (see Figure 3 of Emsellem et al. 2011). These error bars represent the standard deviation in values for all of the different moments in time of the merger. We observe a larger difference in the values as a function of time (Section 4.2.1) than as a function of viewpoint.

Zoom In Zoom Out Reset image size

We next investigate how the LDA classification changes as a function of viewpoint. To visualize this, we plot the distribution of LD1 values in Figure 11. We include the histograms of the LD1 value for nonmerging (blue) and merging (red and orange) galaxies in both the q0.5_fg0.3 simulation (left panel), and the q0.2_fg0.3_BT0.2 simulation (right panel). We then overplot the mean and standard deviation of the LD1 values for all snapshots of each specific viewpoint. Focusing on the mean LD1 values for the merging sample, we can determine whether the LDA shows variation as a function of viewpoint.

Zoom In Zoom Out Reset image size

Focusing first on the left panel of this figure, which refers to the q0.5_fg0.3 simulation, the means for both nonmerging and merging galaxies are not significantly different when considered as a function of viewing angle. In fact, the mean LD1 values are more similar than the variation we observe in LD1 as a function of time in Figure 7. The implication is that the major mergers are LDAs that are fairly robust to changes in viewing angle. For the q0.2_fg0.3_BT0.2 simulation, we observe slightly more variation in the LD1 distribution as a function of viewpoint, and the most face-on viewpoints (0 and 4) have lower LD1 values, which result in more false-negative detections at these viewpoints.

To further quantify changes in the LDA as a function of viewpoint, we iteratively drop the merging galaxies at each viewpoint from the analysis, and rerun the classification for the q0.5_fg0.3 and q0.2_fg0.3_BT0.2 simulations. If the classification changes, this could indicate that a given viewing angle and/or inclination is significantly more or less accurate than the other viewpoints, which would point to inclination itself being the primary driver of this difference.

From here on, we determine that the LDA is “significantly different” from the fiducial run if either of the following criteria are met: firstly, the majority of the selected predictors in the top four selected terms must change, or secondly, the performance statistics in Table 4 must change by more than 10% on average. This quantification of a significantly different classification applies to this section, where we explore the role of different viewing angles, and also to Section 3.9, where we experiment with changes in the data reduction (i.e., changing the S/N or redshift of the simulated galaxies).

When we rerun the LDA classification for the q0.5_fg0.3 simulation iteratively without each viewpoint, the LDAs are not significantly different from the fiducial run. This confirms our findings in Figure 11. For the q0.2_fg0.3_BT0.2 classification, we find that when viewpoints 2, 5, and 6 are absent, the classification is significantly different, with lower performance statistics, and different selected predictors. Our interpretation is that the minor mergers are best identified when the secondary nuclei are within the field of view, which occurs most often in viewpoints 2, 5, and 6. Therefore, the significant changes to the classification as a function of viewpoint both in Figure 11, and in the rerun of the LDA without these viewpoints, can be attributed to the chance positioning of the secondary galaxy as a function of viewpoint. This means that inclination-related effects on the intrinsic kinematic properties of the primary galaxy are not primarily responsible for the differences in the q0.2_fg0.3_BT0.2 LDA as a function of viewpoint. As a final note, in Appendix E in our discussion of between-class biases, we introduce inclination itself as a predictor in the LDA. We ultimately determine that , which we use as a proxy for inclination, is not an important predictor. This further supports the finding that changes in the kinematic predictors purely due to inclination effects do not bias the LDA classification itself.

To conclude, we have determined that while the kinematic predictors themselves can vary as a function of viewing angle and/or galaxy inclination, the LDA classification is only sensitive to viewpoint in the sense of the relative positioning of the secondary nucleus relative to the line-of-sight vector.

As with the imaging identification technique in N19, the kinematic technique is sensitive to both S/N and resolution, meaning that as the S/N decreases in the spectra and/or as the redshift of the galaxy increases, the technique will undergo significant changes. We test this sensitivity by decreasing the S/N, and by moving the mock galaxies to a higher redshift.

To test the sensitivity of the classification to decreased S/N, we decrease the average S/N of q0.5_fg0.3 simulation by factors of 1.5 and 2, respectively. In Figure 12, we compare a snapshot with S/N that has been decreased by a factor of 2 to the same snapshot from the fiducial run. When the S/N is decreased by a factor of 2, the classification is significantly different. While many of the predictors remain the same, the performance statistics decrease overall, and there is an increase in the number of false negatives during the early stages of the merger. When the S/N is decreased, the Voronoi bins increase in size in the exterior regions of the galaxy. This obscures the large-scale kinematic features, which lowers the performance of the classification. We therefore predict that MaNGA galaxies with low S/N are more likely to be misclassified.

Zoom In Zoom Out Reset image size

We plan to explore two possible approaches to classifying MaNGA galaxies with different S/N ratios. One option is to implement an S/N cut when we apply the fiducial classification to the MaNGA galaxies. However, MaNGA galaxies with lower S/N need not be excluded from classification. Instead, another option is to use the classification and completeness correction from the lower S/N ratio LDA run, so as to classify these galaxies separately. We predict that we should be able to classify the majority of the MaNGA sample, since the fiducial (non-decreased) S/N of the simulation suite is representative of the MaNGA sample.

For the purposes of comparison, in Figure 13 we present a sample of MaNGA galaxies across a range of surface brightness, redshift, and stellar mass. The approximate stellar mass is taken from the NSA catalog, and is estimated using the kcorrect code (Blanton & Roweis 2007). The simulated galaxies, whose stellar masses are log M_* ∼ 10.6, are intermediate-mass galaxies, relative to the full MaNGA sample.

Zoom In Zoom Out Reset image size

We also experiment with increasing the S/N by a factor of 2. The classification undergoes minimal changes, with a slight increase in the performance statistics. When the S/N is higher, we predict that the current classification will be better able to determine whether a galaxy is merging. Since the selected predictors do not significantly change, we can conclude that the fiducial run is successfully identifying all relevant kinematic features.

The mock galaxies are placed at a redshift of z = 0.03, which is the median redshift of galaxies observed by MaNGA. In order to understand the limitations of identification over the full range of redshifts for the MaNGA survey (0.01 < z < 0.15), we experiment with increasing the redshift of the mock galaxies. To do this, we increase the spaxel size from 05 to 10, and 15, and we increase the PSF size to 50 and 75. This mimics the effects of moving the simulated galaxies to a redshift of z = 0.07 and z = 0.1, respectively. When we artificially redshift the galaxies, we do not introduce cosmological dimming, meaning that the galaxies have the same S/N as the sample at z = 0.03. This is because we want to understand the effects of the apparent size of galaxies independently of S/N effects.

The classification does not change significantly when the galaxies are placed at a redshift of z = 0.07. At z = 0.1, the classification is significantly different; the number of misclassifications increases, and the selected terms are different. However, these results are based on a galaxy–galaxy merger, where each galaxy has a stellar mass on the order of 10¹⁰ M_⊙. While it is a valid conclusion that the technique will struggle with an intermediate-mass galaxy at z = 0.1, the MaNGA sample does not tend to include this type of galaxy. Instead, MaNGA is designed to maintain roughly uniform coverage in log M_* and radial coverage, meaning that higher mass galaxies (>10¹¹ M_⊙), which are more luminous and have larger angular sizes (i.e., the fourth row of Figure 13), are observed primarily at higher redshift, somewhat alleviating this concern (Bundy et al. 2015; Wake et al. 2017).

The simulation suite of merging galaxies is limited in stellar mass and B/T ratio. The mergers can all be characterized as intermediate-mass disk-dominated galaxies, spanning a range of (3.9–4.7) × 10¹⁰ M_⊙ in stellar mass, and 0–0.2 in initial B/T ratio. These limitations are of particular importance, given that many of the leading kinematic predictors (, μ_1,V , μ_1,σ , μ_2,V , and μ_2,σ ) are related to the intrinsic kinematic properties of galaxies. For instance, μ_1,σ is a proxy for stellar mass, so we are skeptical if this classification can reliably be applied to galaxies that differ in properties, i.e., bulge-dominated, elliptical galaxies.

One possible approach to circumvent these concerns is to remove these predictors from the classification. We rerun the LDA for all simulations, without the , μ_1,V , μ_1,σ , μ_2,V , and μ_2,σ predictors, and find that performance significantly decreases. Specifically, the accuracy, recall, and F1 score of the major mergers decrease by 20%–50%. This is unsurprising, given that the leading predictors presented in Table 2 include all of the predictors that are tied to the intrinsic kinematic properties of galaxies (, μ_1,V , μ_1,σ , μ_2,V , and μ_2,σ ). Interestingly, when the intrinsic kinematic predictors are excluded, the performance of the minor merger simulations is not significantly affected. In fact, the performance of the major merger simulations is comparable to, or worse than, that of the minor merger simulations. This highlights the fact that major mergers undergo a more dramatic global transformation during the merging process; this is reflected in the intrinsic kinematic properties of the remnant galaxy.

Since removing these predictors significantly decreases the performance of the classification, we choose to include all predictors in this work, and to attach the following caveat to this paper: Since this analysis focuses on kinematic predictors that are sensitive to intrinsic galaxy properties, we advise against applying this classification to all galaxy types in MaNGA. In Section 4.9, we discuss possible strategies for carefully applying the classification to MaNGA galaxies.

Some of the kinematic predictors traditionally utilized to identify merging galaxies are uninformative in this analysis. In this case, “uninformative” means that the predictor is discarded during the RFR term selection steps, or that it has a small relative coefficient value in the LDA. The uninformative predictors are ΔPA, v_asym, σ_asym, Δx_V , and Δx_σ .

4.1.1. The Misalignment between Kinematic PA and Imaging PA (ΔPA) Is Most Sensitive to Galaxy Inclination

4.1.2. The Asymmetry in the Velocity and Velocity Dispersion Maps (v_asym and σ_asym) Is Only Sensitive to the Most Disturbed Times in the Major Mergers

Previous works based on the gas kinematics of simulated and observed mergers have found that merging galaxies have enhanced values for both v_asym and σ_asym (e.g., Shapiro et al. 2008; Bellocchi et al. 2012; Hung et al. 2016; Bloom et al. 2017). These studies define a threshold value in K_asym to identify merging galaxies, where . For instance, Bellocchi et al. (2012) studied local luminous infrared galaxies (LIRGs), finding a threshold value of K_asym > 0.135. Hung et al. (2016) defined a K_asym threshold of 0.15 for the galaxies in their work, calculated from the velocities of gas particles for simulated SUNRISE mergers.

The σ_asym and v_asym predictors are ultimately unimportant in our analysis. We present the viewpoint-averaged values of σ_asym-v_asym in Figure 14 for the q0.2_fg0.3_BT0.2, and the q0.5_fg0.3 simulations. We include the K_asym diagnostic lines used to identify the mergers from Bellocchi et al. (2012) and Hung et al. (2016). The top panel of the plot demonstrates that there is minimal time evolution in the predictor values for minor mergers. The predictor values for major mergers are only slightly enhanced, falling above the diagnostic line of Hung et al. (2016) for a few snapshots during the late stages of a merger. While Figure 14 presents the predictor values in log space, the standardized predictor values of σ_asym and v_asym used to construct the classification also have minimal separation between the merging and nonmerging populations. The insensitivity of any of the mergers to an enhancement in the values of these predictors results in their exclusion during the RFR selection step.

Zoom In Zoom Out Reset image size

Ultimately, the v_asym and σ_asym predictors are unimportant in our work, because v_asym is only elevated at specific points in time during the late phases of merging. They are more useful in studies such as that of Bellocchi et al. (2012), who focused on LIRGs, which are prototypical gas-rich major mergers, or that of Hung et al. (2016), where the simulated galaxies include gas-rich major mergers. It is also of interest to note that Hung et al. (2016) found that the mergers in their sample exceeded the K_asym value only during the “strong interaction,” or late stage of merging.

4.1.3. The Offset between the Center of the Velocity and Velocity Dispersion Maps and the Imaging Center (Δx_V and Δx_σ ) Are Not Very Sensitive to Mergers

Here, we examine the most important predictors in the LDA for all simulations, and make connections between these predictors and the dynamical evolution of stars during a merger. We split the discussion by predictor, and for brevity, we focus only on the leading predictors presented in Table 3. The top predictors include , μ_1,σ , μ_2,σ , ∣μ_3,σ ∣, and μ_4,V .

4.2.1. The Approximate Spin Parameter Tracks a Global “Slowdown” in the Velocity Maps of the Major Mergers

The approximate spin parameter () is a key predictor for all merger simulations, and is especially important for the major mergers. The angular momentum of merging galaxies is therefore significantly different from that of the nonmerging population; this effect is more apparent in major mergers. In this section, we examine how changes with time for the various simulations, and how this compares to previous works.

Firstly, we directly examine the pre-standardized values of and for the q0.2_fg0.3_BT0.2 and q0.5_fg0.3 simulations in Figure 15. In this diagram, we indicate the “slow rotator” territory, located in the lower left corner. The – diagram is often used to kinematically distinguish the slow rotator population of early-type galaxies from the fast-rotating population. This fast-slow rotator distinction probes the evolutionary histories of galaxies, based on disk assembly (see Cappellari 2016 for a review).

Zoom In Zoom Out Reset image size

Much recent work has focused both on examining the observed populations of fast and slow rotators, and on making predictions as to how merging galaxies move through this territory. For instance, Naab et al. (2014) utilized cosmological merger tree simulations to show that major mergers significantly affect the angular momentum content of a galaxy; they can either spin up or spin down the remnant. In our case, all of the simulated galaxies begin with a value of ∼0.7. For the major mergers, the value dramatically decreases, to the boundary of the slow-rotator region. This confirms that major mergers can dramatically affect the kinematic properties of the remnant, kinematically transforming the galaxy from one that can be described as disk-dominated to one that is still rotating, but is dispersion-dominated.

4.2.2. The Mean and Dispersion of the Velocity Dispersion Distribution (μ_1,σ and μ_2,σ ) Track the Growth of a Stellar Bulge Component

The overall importance of the μ_1,σ and μ_2,σ predictors reflects the fact that the velocity dispersion map is informative for identifying mergers. We examine how these two predictors evolve with time during a merger in Figure 16, for the q0.5_fg0.3 and q0.2_fg0.3_BT0.2 mergers. Here, we present the average value for all of the viewpoints of a given snapshot. We also include representative velocity dispersion maps for a handful of informative snapshots.

Zoom In Zoom Out Reset image size

For all merger simulations, the μ_1,σ value increases throughout a merger, tracing the assembly of a stellar bulge component. This increase is more dramatic for the major mergers, which increase to a μ_1,σ value of ∼200 km s⁻¹. Even for minor mergers, the merger incites growth in the central velocity dispersion with time. This enhancement is still present 0.5 Gyr after coalescence, such that isolated post-coalescence stages are mixed with the merger snapshots along the μ_1,σ axis in Figure 16. The signatures of bulge growth are therefore dynamically long-lived, as opposed to imaging features that fade quickly with time following a merger (i.e., in N19, the imaging predictors fade within 0.5 Gyr of final coalescence).

The μ_2,σ predictor serves different roles in major versus minor mergers, as reflected in the differing evolution of μ_2,σ values with time. An increase of μ_2,σ with time for major mergers traces the presence of two kinematic components by capturing the “bridge” of higher velocity dispersion values between two merging galaxies. This is formed by two overlapping counter-rotating features. Moreover, post-coalescence, major mergers exhibit more significant bulge growth, which is reflected both in the enhancement in μ_1,σ , and in an increase in μ_2,σ , since the entire distribution is broadened in this process.

Minor mergers show less change in the value of μ_2,σ with time. While μ_2,σ is still informative (since it increases during the late stages of the merger to track the bridge of higher velocity dispersion), it does not continue to increase into the post-coalescence stages. This could indicate that a smaller fraction of stars are involved in the buildup of the stellar bulge in the case of minor mergers.

4.2.3. The Skewness of the Velocity Dispersion Distribution (∣μ_3,σ ∣) Identifies Secondary Kinematic Components

4.2.4. The Kurtosis of the Velocity Distribution (μ_4,V ) Identifies the Superposition of Two Merging Galaxies

Past studies have investigated how the properties of simulated mergers affect the kinematic predictors. Hung et al. (2016) investigated this phenomenon in relation to a set of simulated merging galaxies with mass ratios of 1:1 and 1:4. They found that the merger signatures in kinematics are most apparent for the 1:1 major merger, where they may be visible for up to twice as long as for the 1:4 major merger. While there are significant differences between the work in this paper and that of Hung et al. (2016) (i.e., we perform a full RT and create mock IFS maps, whereas Hung et al. (2016) use the GADGET-3 particle velocities to create velocity maps), we also find that the classification differs significantly with mass ratio.

The major and minor merger classifications are different in several ways. The minor merger classifications have performance statistics that are ∼10%–30% lower. The minor mergers are also more unstable in the terms that are selected, meaning that the coefficient values fluctuate slightly when the classification is rerun. The LDA for the minor mergers therefore uses a small sample of key snapshots to create the classification. If these snapshots are excluded from the training set, and instead fall in the CV set, then the classification is slightly different. The overall effect is that the minor mergers are less stable, and many of the selected terms have similar coefficient values, making it difficult to assess which are the most important.

Another difference between the major and minor merger classifications is that they are composed of different predictors. While some predictors, such as μ_1,σ and μ_2,σ are important for all classifications, is more important for major mergers and μ_4,V is more important for minor mergers. As we discuss in Section 4.2, both major and minor mergers demonstrate bulge growth, which leads to an enhancement in μ_1,σ , but the change is more apparent in major mergers. The global kinematic properties of major mergers are more significantly impacted; this includes the predictor, which traces a global slowdown in the velocity maps of the major mergers. On the other hand, the minor mergers are most sensitive to smaller-scale changes in the kinematic maps; this can be traced by predictors such as μ_4,V , which track the superposition of the secondary stellar nuclei.

Many of the kinematic predictors in this analysis evolve with time throughout the merger. In most cases, this evolution is non-monotonic, meaning that merging galaxies evolve back and forth in predictor space as a function of merger time. The LDA technique accounts for this behavior using interaction terms.

An example of an interaction term in action is the μ_σ ∗ σ_σ term, which has a negative coefficient for several of the major merger simulations. This means that if the value of μ_1,σ is relatively large, then μ_2,σ must be relatively small in order for the merger probability to increase. The opposite is also true: if μ_1,σ is relatively small, then μ_2,σ must increase. To be clear, “relatively large” and “relatively small” refer to the standardized values of these predictors, which are measured relative to the distribution of values for the entire merger. As such, if a term is relatively small, the predictor value will be negative. For instance, if μ_1,σ is large and μ_2,σ is small, then the term becomes: (− coefficient) ∗ (positive) ∗ (negative) = positive value = increase in LD1.

Consider the μ_σ − σ_σ diagram for the q0.5_fg0.3 merger, presented in the left panel of Figure 16. The μ_1,σ and μ_2,σ predictors are correlated, so they occupy the diagonal space of this diagram. This correlation and the interaction term function so that the center of the μ_σ − σ_σ diagram is the “merger territory.” This picture is somewhat complicated by the fact that there are other terms in the LDA, but if the μ_σ ∗ σ_σ interaction term has a large coefficient, then this interpretation generally applies.

The interaction terms account for the non-monotonic evolution of the predictor values with time. When we create a classifier with only the primary predictors, it is fundamentally different from the LDA with interaction terms. The linear classifier is inaccurate, classifying many of the isolated post-merger snapshots as mergers. When the LDA is forced to generalize to one direction of movement across these diagrams in the monotonic case, it loses key information.

When the interaction terms are included, the LDA becomes sensitive to the values of the other predictors, and is therefore able to create “if-then” cases for the predictor space. For example, if a galaxy is in the later stages of merging (characterized by significant bulge growth and a large value of μ_1,σ ), then to avoid ambiguity with post-merger isolated galaxies, the classification would like to see a relatively negative standardized value of μ_2,σ in order to classify the galaxy as merging. On the other hand, if the galaxy has a relatively small value of μ_1,σ , then μ_2,σ should be large in order for the galaxy to be classified as merging. An example of this type of galaxy can be seen in the early stages of a merger, where the bulge growth has not yet begun, so μ_1,σ is small. However, the galaxy has undergone its first pericentric passage, and is experiencing an enhancement in velocity dispersion, which leads to an increase in μ_2,σ . An example of this type of galaxy is shown in the left panel of Figure 16.

We conclude that it is critical to include the interaction terms in the LDA classification. Not only do they improve the performance of the technique, but they are physically motivated by the non-monotonic evolution in kinematic predictors over the course of a merger lifetime.

In Section 3.6, we present the observability timescales for the various simulated mergers. We conclude that the kinematic LDA technique lengthens the observability timescale of the simulated mergers as compared to that measured from the individual predictors. Here, we focus specifically on how the observability of a merger changes with the merger stage, and we refer the reader to Figure 7 for a useful visualization of the mean value of LD1 with time in all of the simulated mergers.

When we examine the observability of the mergers in terms of the predefined early-, late-, and post-coalescence stages, we notice several stage-related differences. Firstly, during the early stages of merging, some simulations show a larger standard deviation in the LD1 values. For all simulations, the value of LD1 tends to fall below the decision boundary during these early stages. Hung et al. (2015) found that using kinematic predictors to identify mergers results in a significant fraction of false negatives from epochs where the merger is indistinguishable from a rotating disk. We also find that a significant fraction of false negatives occur during the early stages, where the rotation is indeed disk-like.

During the late stage of the merger, the minor mergers do not show much variation from one snapshot to the next; the LD1 values are relatively flat. On the other hand, the major mergers, in particular q0.5_fg0.3, show variation between the late-stage snapshots. For q0.5_fg0.3, this factor contributes to the relatively short observability timescale. These changes in LD1 values are significantly greater than the variance due to different viewpoints. The kinematic features are therefore seen to change significantly and rapidly with time for some of the major mergers during their late stages. We also find that most simulations have relatively high LD1 values during the late stage of a merger. The late stage of a merger is therefore characterized by short-lived dramatic kinematic features. This is consistent with Hung et al. (2015), who found that kinematic tracers of mergers tend to be most informative during the late stage of the merger, which is when the imaging predictors are also most useful.

As the merger progresses into the post-coalescence and post-merger isolated stages, we find that the LD1 value is more stable, which is characteristic of longer-lived kinematic features. The LD1 value does not significantly decline during the post-merger stages. We focus on the kinematics of the post-coalescence and post-merger stages in Sections 4.6 and 4.7.

In N19, we define the end of the merger as 0.5 Gyr after final coalescence. This cutoff was selected so that the imaging predictors, and therefore the value of LD1, would decay smoothly until the end of the merger. The imaging technique is therefore very accurate and precise during the transition from a post-coalescence merger to an isolated post-merger galaxy. In contrast, in this work, the kinematic predictors, and therefore the LD1 value, remain elevated during the period of post-merger isolated snapshots. There are visually apparent warps in the velocity maps, and the velocity dispersion maps have elevated values of μ_1,σ and μ_2,σ for post-merger isolated galaxies. Hung et al. (2016) also found a persistence of kinematic merger signatures up to ∼Gyr following coalescence. We find that the kinematic disturbances fade 2–2.5 Gyr after coalescence, meaning that in order to improve the classification, we would need to significantly extend the post-coalescence phase.

Rather than changing the definition of the end of a merger, a more relevant task might be to define the merger more specifically by stage. For instance, if it is a priority to distinguish post-coalescence from post-merger isolated galaxies (post-coalescence snapshots occur immediately following coalescence until 0.5 Gyr after coalescence, and post-merger snapshots occur >0.5 Gyr after coalescence), we should rely on the imaging predictors, which do a better job in this specific case. In our future work we plan to combine the imaging and kinematic tools, and address this question more directly. One potential path forward could be to create separate classifications to target different stages of a merger. This could provide a more flexible definition of merging galaxies, and allow other users of the technique to target stages of interest in a merger.

Here, we focus on the kinematics of the post-coalescence and post-merger stages. During these stages we observe the buildup of a central component in the stellar velocity dispersion maps. The change is more dramatic for major mergers, and can best be explained by tracing the growth of a stellar bulge component. Hopkins et al. (2009) investigated the effect of mass ratio on the merger remnant, finding that the fraction of the primary stellar disk that is relaxed into the bulge is directly proportional to the mass ratio of the merger. This supports the hypothesis that a stellar bulge is built in the post-merger stages, since we would expect major mergers to contribute a larger fraction of stars to the bulge component.

The q0.5_fg0.3 merger has unique kinematic features in its stellar velocity map during the post-merger stages, so we focus the remainder of our discussion on this merger. We present several post-coalescence and post-merger snapshots from the q0.5_fg0.3 merger in Figure 17. The stellar velocity maps are particularly intriguing, as they have a distinct central component that appears in the post-coalescence phase, and persists into the post-merger phase. This central feature in the stellar velocity maps is spatially coincident with the bulge-like feature in the stellar velocity dispersion maps. It is not fully counter-rotating, but is misaligned from the main stellar disk. We therefore hypothesize that we have discovered a decoupled kinematic component.

Zoom In Zoom Out Reset image size

Previous theoretical works have predicted that major mergers with mass ratios of 1:1 or 1:2 may produce this type of intriguing kinematic component (e.g., Bendo & Barnes 2000; Jesseit et al. 2007; Crocker et al. 2009). For instance, Bendo & Barnes (2000) and Jesseit et al. (2007) found that equal mass simulated mergers display a far wider range of kinematic features, including counter-rotating cores and global misalignments, whereas unequal mass mergers tend to have disk-like kinematics. The merger remnants with complex kinematics are intriguing, because these decoupled central components have also been discovered in observational studies. For example, the ATLAS^3D survey found that a significant fraction of slow-rotating ETGs have decoupled kinematic components (e.g., Emsellem et al. 2011).

Our finding of a decoupled kinematic core in the remnant of the 1:2 (q = 0.5) mass ratio merger supports these past findings, i.e., that mergers with a large mass ratio can produce dramatic central features in the kinematic maps. Furthermore, this result suggests that selecting galaxies with kinematically decoupled cores would produce a sample of post-coalescence major mergers with a mass ratio q ≳ 0.5.

The kinematic LDA technique gives rise to a significant number of false negatives, which drive down the accuracy and recall of the technique. This can be partially ascribed to the chosen priors, which skew the classification toward minimizing false positives (see Section 3.5 for a full discussion). As we discuss in Section 3.7, this is also due to the lack of identifiability of certain snapshots as mergers, meaning that they are indistinguishable from nonmerging disks and/or post-merger remnants, both visually, and based on their kinematic predictors.

The imaging LDA performs better for all runs, with improvements in all performance statistics. This means that fewer galaxies in total will be correctly classified by the kinematic technique, due to the shortcomings listed above. While the major merger combined run is slightly improved when performed using imaging predictors (this improves the accuracy/recall/F1 score by 9%/11%/7%), the minor merger combined run experiences significant improvement (16%/41%/25%). The kinematic minor merger combined classification therefore scores lower across all performance statistics, relative to both the imaging minor merger classification and the kinematic major merger classification. This reflects the specific inability of the kinematic classification to identify minor mergers.

The imaging LDA also has longer observability times on average, ranging from observability times of 2.2–2.8 Gyr for major mergers, and 3.5–9.2 Gyr for minor mergers. The kinematic LDA has observability times of 0.9–2.4 Gyr for the major mergers, and 3.0–6.6 Gyr for the minor mergers. As such, while the observability times for the kinematic LDA significantly improve upon the observability times for individual kinematic predictors, the mergers are still observable for slightly longer timescales with the imaging LDA.

This has important implications for the relative capabilities of imaging versus kinematic predictors in identifying merging galaxies. The stellar kinematics of mergers take longer to exhibit disturbance, and remain disturbed for a longer time following a merger. The imaging predictors are better contained to the duration of the merger, and better able to identify all of the different merger stages. Hung et al. (2016) recommended using kinematic predictors in combination with imaging predictors, due to the frequency of false negatives in their investigation of kinematically identified merging galaxies. We support this conclusion, finding that the kinematic predictors have failure modes that may be ameliorated by incorporating imaging predictors.

However, there are some advantages of kinematic predictors relative to imaging predictors. In Sections 4.6 and 4.7, we find that the kinematic predictors are particularly useful for identifying post-coalescence and post-merger stages, because the kinematic disturbances persist long after the imaging predictors fade. The inference is that while imaging predictors are more informative overall, kinematic predictors are powerful in certain domains, and provide additional information. If forced to select between imaging and kinematic classification methods, the imaging approach is better. However, the best overall approach is to combine the two techniques into one imaging + kinematic classification. We plan to discuss this topic further in future work.

In Section 3.10, we discuss the implications of creating a kinematic classification from a suite of simulations with a narrow range in stellar mass (3.9 × 10¹⁰ < M_⊙ < 4.7 × 10¹⁰), and in an initial B/T ratio (0–0.2). Many of the kinematic predictors, specifically the most important predictors for the major merger simulations, also probe the intrinsic properties of galaxies. In other words, while these predictors are useful in identifying major mergers, they also change as a function of stellar mass and morphology. The implication is that a classification created from disk-dominated intermediate-mass mergers may not extrapolate well to the MaNGA sample, which spans a wider range in stellar mass (10⁸ < M_⊙ < 10¹¹) and morphology.

At present, it is unclear whether this will be a concern only for extreme cases, i.e., the most massive bulge-dominated ETGs, or if it will also cause concern in systems with a mix of rotation and dispersion support in their kinematics. We have preliminarily investigated the distribution of bulge- versus disk-dominated galaxies in MaNGA; while Wang et al. (2020) have found that MaNGA galaxies are predominantly disk-dominated, Graham et al. (2018) found that a significant fraction of MaNGA galaxies (across all masses) are bulge-dominated.

Since the MaNGA sample includes a diversity of different galaxy types, we plan to tread carefully when we apply the classification. One option may be to select those MaNGA galaxies with high values of to include in the classification; in this way, we could exclude bulge-dominated galaxies. Another option could be to de-emphasize certain domains of predictor space in the classification, or to remove those kinematic predictors that are most sensitive to intrinsic galaxy properties altogether. The details of this approach will be developed in future work, as they fall beyond the scope of this paper, which focuses mostly on the creation of the merger identification technique.