A Simulation-driven Deep Learning Approach for Separating Mergers and Star-forming Galaxies: The Formation Histories of Clumpy Galaxies in All of the CANDELS Fields

IllustrisTNG is a suite of cosmological, gravo-magneto-hydrodynamical simulation runs with a diverse set of particle resolutions. From highest to lowest resolution, simulations were realized in three comoving simulation boxes of 50, 100, and 300 Mpc h⁻¹ length size, aptly named TNG50, TNG100, and TNG300 (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019; Pillepich et al. 2018a, 2019; Springel et al. 2018).

Our analysis makes use of the TNG100-1 simulation, which has proven to be a good compromise between resolution and volume. ⁴ TNG100 has already been used extensively in studies that analyze galactic morphologies and structures, including the comparison between simulations and observations (Huertas-Company et al. 2019; Blumenthal et al. 2020), and the use for deep learning (Wang et al. 2020; Bickley et al. 2021; Bottrell et al. 2022). Specifically, Zanisi et al. (2021) showed that TNG100 galaxies reproduce observed objects well, especially disk-dominated sources. While there are some deviations in the small-scale structure of highly concentrated spheroidal systems, this is a minor issue in our analysis since they only make up a small fraction of our sample. In addition, our galaxies are resolution limited at the current redshift of interest, meaning that tiny details of structure are not relevant in this analysis.

To counterbalance any limitations from resolution, we limit our analysis to galaxies with M_* > 10^9.5 M_⊙. Above this limit, and at our explored redshift range z > 0.5, galaxies are represented by thousands of stellar particles. This enable sampling the simulated galaxies into resolutions comparable to that of the observed CANDELS data. Specifically, the gravitational softening length of the simulation, , is not a limitation when compared to the HST Advanced Camera for Surveys (ACS) and WFC3 cameras resolution.

This approach is a noticeable refinement to our previous treatment in Ferreira et al. (2020), where the research question did not demand the resolution of fine morphological features like clumpy regions and tidal features, which are required for the present analysis.

To select appropriate galaxies from TNG100-1, we isolate galaxies with M_* > 10^9.5 M_⊙, in the redshift range 0.5 < z < 3. To limit contamination in our sample, we use a minimum dark matter to total mass ratio of

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{DM}}}{{M}_{\mathrm{total}}}\gt 0.1,\end{eqnarray} \tag{ 1 }$$

as a way to avoid subhalos created as a result of disk fragmentation. This means that at least 10% of the subhalo's mass needs to be in the form of dark matter. We acknowledge that this could also inadvertently remove galaxies that had their dark matter stripped; however, this number is small and does not impact the final sample. This criteria removes ≈2% of galaxies from the overall pool of available sources (subhalos) in TNG100.

We also remove objects that are smaller than the ACS point-spread function (PSF) size from the selection. To identify these objects, we first convert the half-mass–radius R_1/2 provided in the simulation group catalogs in kiloparsecs, to a pixel scale based on the cosmological model adopted by IllustrisTNG and the ACS pixel scale

$$\begin{eqnarray}&&{R}_{1/2\mathrm{Mass},\mathrm{pix}}=\displaystyle \frac{a(z)}{h}{R}_{1/2\mathrm{Mass},\mathrm{kpc}},\end{eqnarray} \tag{ 2 }$$

where a(z) is the angular size at z, and h is the Hubble constant/100. Any galaxy with R_1/2Mass,pix < 3 pix was then filtered out from our selections. This step removes ≈3% of galaxies from the total pool of available sources.

Our goal is to separate star-forming galaxies from post-mergers at intermediate to high redshifts based on their morphology. We define post-mergers as galaxies with at least one major-merger event with a mass ratio

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{M}_{2}}{{M}_{1}},\quad \mu \geqslant 0.25,\end{eqnarray} \tag{ 3 }$$

where M₁ and M₂ are the stellar masses of the galaxy pair involved in the merging event, ranked by their stellar mass, respectively, with M₁ > M₂. Galaxies are considered post-mergers if they have coalesced into a single galaxy in the past 500 Myr, where a single galaxy is represented by a subhalo in the simulation as identified by friends-of-friends algorithms (Rodriguez-Gomez et al. 2015). This selection window timescale is motivated by the observability timescales of disrupted structures caused by mergers identified by structure measurements in IllustrisTNG (Whitney et al. 2021), and is higher than what was previously used in Ferreira et al. (2020). We allow post-mergers to have low sSFRs. Their asymmetric features likely arise from the merging process rather than from star-forming clumps. In contrast, noninteracting star-forming galaxies are defined here as galaxies that have sSFRs above the following threshold:

$$\begin{eqnarray}&&\mathrm{sSFR}\gt {10}^{-9.5}\ {\mathrm{yr}}^{-1}\end{eqnarray} \tag{ 4 }$$

and are not interacting with other galaxies. To isolate noninteracting cases, we exclude any galaxy from the simulation that had major or minor merger events (μ > 0.1) around ±1 Gyr of its current redshift. Minor mergers are excluded completely from both definitions, and any conclusions presented in this paper should be considered with this in mind. Importantly, this selection is not intended to limit the noninteracting cases to extreme starbursting episodes alone, but to select noninteracting galaxies with sufficiently high sSFR to produce clumpy and asymmetric features that could be mistaken for merging signatures.

In summary, this selection results in a sample of ∼6000 post-mergers and ∼110,000 noninteracting star-forming galaxies. While this may be a realistic representation of actual fractions (only ∼5% of the sample are post-mergers), training the network requires a balanced data set. We thus use the post-merger sample as the baseline and separate it in bins of redshift, stellar mass, and size, randomly sampling the same number of noninteracting galaxies within each bin. We remove bins without adequate matched numbers of star-forming galaxies. This becomes noticeable in the higher mass bins where post-mergers dominate and very few star-forming galaxies are present.

After matching the samples, we count ∼4000 galaxies in each class as our final sample. A summary of this sample separated by class and redshifts is available in Table 1. The distribution of redshifts, star-forming rates, stellar masses, and stellar half-mass–radius is shown in Figure 1 for post-mergers in red, and star-forming galaxies in blue. Both classes have very similar physical properties, with a small excess of large, passive and massive post-mergers in comparison to the star-forming galaxies. Additionally, the top-right and bottom-right panels of Figure 1 show the time since the last major-merger event, τ, and the mass ratio, respectively, for post-mergers. The nature of the distribution for τ arises from the average time between snapshots in the simulation of around ∼0.15 Gyr. This timescale represents one to three snapshots after the coalescence of stellar masses.

Zoom In Zoom Out Reset image size

One of the main goals of this work is to predict star-forming and post-merger galaxies in the observed CANDELS imaging data (Grogin et al. 2011; Koekemoer et al. 2011), which comprises high-quality HST observations from COSMOS, UDS, EGS, GOODS-South, and GOODS-North (Grogin et al. 2011; Koekemoer et al. 2011). CANDELS data has been used extensively for galaxy merger studies, with estimated merger rates up to z ∼ 6 (e.g., Mantha et al. 2018; Duncan et al. 2019; Whitney et al. 2021). Importantly, CANDELS also provides visually classified morphologies (Kartaltepe et al. 2015), as well as photometric redshifts, star formation rates, and stellar mass estimates (Conselice et al. 2007; Duncan et al. 2014, 2018a, 2018b, 2019) from spectral energy distribution (SED) fitting, essential for creating matched samples. Stellar masses are estimated within a 0.2 dex systematic uncertainty, and the outlier fraction of photometric redshifts is smaller than 5%. We do not account for these uncertainties directly in our methods, but we refer the reader to the detailed discussion on the uncertainties associated with these physical measurements in the aforementioned references. Furthermore, the depth of the wide field data is comparable to the detection limit of galaxies in the IllustrisTNG simulations in the stellar mass range studied here.

To select CANDELS galaxies, we first remove all problematic objects according to their quality flags as recorded in the photometric catalog and in the Kartaltepe et al. (2015) catalogs to avoid edges, artifacts, and stars. Following Huertas-Company et al. (2016) and Kartaltepe et al. (2015), we then select galaxies with H-band magnitudes H < 24.5 mag. Because this cut can bias our sample against extended sources, we also include a signal-to-noise ratio (S/N) lower limit of S/N > 50 to exclude any compact source with only a few bright pixels. This magnitude cut removes 1074 sources, while the S/N cut further removes 430 sources. Then we proceed with the same cuts we used to select IllustrisTNG galaxies, using 0.5 < z < 3 and M_* > 10^9.5 M_⊙. We apply a final cut using the asymmetry (A > 0.1; Section 3.3) to remove regular unambiguous galaxies with no apparent disturbed or asymmetric features. This ultimately results in a sample of 23,494 galaxies from all of the CANDELS fields combined.

Finally, we produce cutouts for I₈₁₄, J₁₂₅, and H_160W bands centering on each selected CANDELS galaxy, each with a field of view of 50 kpc × 50 kpc, using photometric redshifts from Duncan et al. (2019), preserving relative sizes between galaxies. Importantly, this selection does not rely on size measurements that could easily be spurious in interacting or merging galaxies. We do not find any bias in our classifications that could be attributed to small changes of the field of view caused by the photometric redshift uncertainties.

In order to guarantee realistic representations of CANDELS galaxies in the simulated sample, we must include instrumental and cosmological effects to the images of the IllustrisTNG galaxies. An overview of the steps is shown in Figure 2 and is detailed in this section. IllustrisTNG data holds information on the stellar, gas, and dark matter particles for each source. Each particle represents a large physical region that can be described by rich stellar populations, which vary depending on age, mass, and metallicity. The resampling of the star-forming regions is particularly important to avoid problems with the coarse representations (Camps et al. 2016; Trayford et al. 2017).

Zoom In Zoom Out Reset image size

To create mock broadband images, we thus process each stellar particle with a population synthesis model following the recipes from Trayford et al. (2017) and Vogelsberger et al. (2020). This entails post-processing the simulation data with the Monte Carlo dusty radiative transfer code SKIRT (Camps & Baes 2015, 2020).

Each stellar particle in the simulation is considered as a single stellar population with GALAXEV (Bruzual & Charlot 2003) or MAPPINGSIII (Allen et al. 2008) SEDs based on its stellar mass, absolute metallicity, and age. We choose to adopt these particular templates because, first, they are implemented in SKIRT and, second, they had been tested previously in similar pipelines to generate mock observations from cosmological simulations (Trayford et al. 2017; Rodriguez-Gomez et al. 2019). Finally, the templates of Bruzual & Charlot (2003) are also those used to derive stellar masses and star formation rates for all of the CANDELS fields in Duncan et al. (2019) that are used in this study.

To account for the fact that each stellar particle represents an extended area (rather than treating them as a point source), we model the particles with a smoothing length of a truncated Gaussian emissivity profile equal to the distance to its 64th neighbor particle (Trayford et al. 2017). We then define a grid of wavelengths covering all spectral features we want to probe within the HST filter response functions, similar to the grid used in Trayford et al. (2017). For each wavelength bin of this grid, we launch 10⁶ photon packets, assuming isotropic emission until they reach the virtual detector.

This process produces IFU data cubes over the SKIRT wavelength grid, which we then reduce to broadband images with the same properties as the CANDELS HST images. SKIRT's reference frame used to generate the data cubes is located at a distance of 10 Mpc (initial redshift z₀) of the sources. We must therefore shift the IFU data to each target's redshift, z_t , by (1 + z_t ) while dimming its flux by

$$\begin{eqnarray}&&\displaystyle \frac{{f}_{t}}{{f}_{o}}={(1+{z}_{t})}^{-1}{\left(\displaystyle \frac{{D}_{{Lo}}}{{D}_{{Lt}}}\right)}^{2},\end{eqnarray} \tag{ 5 }$$

due to cosmological dimming (Hogg 1999; Equation (15)). Next, we convolve the IFU data with the broadband filter's response functions for I₈₁₄, J₁₂₅, and H_160W . The results are clean, noiseless images from the simulation galaxies at 30 mas pix⁻¹ (matching the ACS pixel scale) before adding any PSF effects (Figure 2, images (a)–(c)). We rebin the J and H bands from 30 mas pix⁻¹ to the WFC3 images pixel scales of 60 mas pix⁻¹. Examples for stamps where the background was added can be seen in Figure 2, images (e) and (f).

Figure 3 gives randomly selected examples of galaxies in our sample before any contamination from the CANDELS sky is included, separated by their class.

Zoom In Zoom Out Reset image size

The data-driven paradigm of deep learning methods imposes high requirements on the amount of data necessary to train a model that is capable of generalizing the training data well. In practice, this means that for the majority of models, a successful approach requires tens, hundreds, or even millions of examples. We are far away from these numbers in cosmological simulations. Our initial selection results in a balanced set of ∼4000 examples of each class (Section 2.1). Fortunately, in the case of galaxy images, there are ways to increase the initial data set by exploiting aspects of the final image that do not depend directly on the simulated galaxy. In our case, we apply data augmentation to our data set in three ways outlined below. An example of this approach is shown in Figure 4, following the same galaxy in each possible combination of orientation/field.

Zoom In Zoom Out Reset image size

First, since IllustrisTNG provides the 3D distribution of all particles associated with a galaxy, we generate each galaxy with different line-of-sight projections, treating each new representation as a new galaxy. We select four different projections, three aligned with the axis of the simulation, XY, XZ, and YZ, respectively, and a fourth line of sight aligned with one octant of the simulation cube.

Second, each CANDELS field has unique observational properties (e.g., different noise levels, and depth). We exploit this aspect and reproduce each of the different orientations from the previous step on top of a random patch of sky of each CANDELS field, taking care to use appropriate noise levels for the simulated galaxy. To find empty patches of sky, we randomly sample the R.A. and decl. within each field, and make a large cutout of the area that is four times larger than the final size of the cutout. Using positions given in the CANDELS catalogs, we then identify all sources within this cutout and reselect a new R.A. and decl. location within the cutout that does not centrally overlap with another source. We allow some degree of overlapping source, but require a unique central position. We do this interactively until a patch of sky that matches all of above criteria is found. This, combined with all of the orientations, augments our data set 20 times. In addition, this also helps the network to generalize the impact of contamination from neighboring sources, as the same galaxy in one field might be isolated in its cutout, but in a denser environment in another.

Finally, we apply random flips, rotations, and small zoomed-in/zoomed-out regions around the central source on the fly during training as a regularization technique. This does not increase the overall size of the sample, but at each training epoch, the network sees different realizations of the same sample.

Overall, our sample increases from ∼8000 examples to ∼160,000. However, having multiples of similar galaxies in our data set can result in overfitting. To reduce this risk, we do not allow different realizations of the same galaxy to fall in both the training sample and the test sample. This ensures that testing and validating are performed on unique data sets.

Providing realistic levels of contamination and the inclusion of neighboring sources are some of the most important requirements for a good generalization between samples of simulated galaxies and real observations (Bottrell et al. 2019). In an update to what was done in Ferreira et al. (2020), we included realistic contamination in our IllustrisTNG sample, as described in Section 2.4. By comparing clean galaxy realizations to their respective background-added images, we can thus test how our methods behave when faced with a variety of contamination levels, drawing direct conclusions for real world applications. We quantify the degree of contamination in each image using two measurements, which are also listed in Figure 4.

First, we define how much of the galaxy is covered by a background source. We call this the overlapping percentage, Θ. For this, we measure segmentation maps both for the central source and all background sources of each image stamp. Θ is the percentage of the segmentation map of the central galaxy that is covered by segmentation map(s) of background sources and ranges from 0%, for no overlap, to 100%, where the central galaxy is completely covered by another galaxy in the field.

Second, we estimate the average flux (per pixel) of all background sources, BG_flux, by averaging the flux of the sources within the segmentation map over its area. BG_flux values are given in units of e s⁻¹ pix⁻¹. This ranges from BG_flux ∼ 0, where there is no apparent or very faint source in the background, to values that are comparable to or even higher than the flux of the central source. Very high values may be due to bright neighboring sources that outshine the central galaxy. Stars can also be identified by this method.

We use the overlapping percentage, Θ, and the flux of background sources, BG_flux, to define galaxy images with low contamination. Figure 5 shows the parameter space formed by these two measurements for the entire sample of ∼160,000 simulated and candelized images. The blue box framed by the dashed line defines a region of galaxies with low contamination,

$$\begin{eqnarray}&&{\rm{\Theta }}\lt 15 \% ,\quad {\mathrm{BG}}_{\mathrm{flux}}\lt \times {10}^{-1.5}\ {\rm{e}}\ {{\rm{s}}}^{-1}\ {\mathrm{pix}}^{-1},\end{eqnarray} \tag{ 6 }$$

which can be considered as a conservative choice. We find that ∼90% of our sample is located in this region. We do not remove the remaining 10% of the galaxies from our sample, because such highly contaminated cases will also be present in observations. We use these contamination estimates to understand how our methods are impacted by it.

Zoom In Zoom Out Reset image size

These two properties form a simple and powerful way to characterize the contamination of our sample, as they control different contributions to contamination. Because these are challenging to measure directly in real CANDELS observations, we trained a deep learning model to predict the same values in real images. We describe this exercise in the Appendix. By inference, any discussion based on contamination measurements in our simulation sample is also valid for the CANDELS observations.

We employ neural networks to forward model the simulations into the observational domain. The neural network takes galaxy images as input and outputs a probability associated with its classification, in this case whether it is a post-merger or a star-forming galaxy.

Neural networks are known for being able to approximate complex functions where no analytical approach is feasible, based on the universal approximation theorem (Lu et al. 2017). Deep neural nets combine several layers of nodes (neurons) in a feed-forward fashion, mapping inputs to outputs using nonlinear activation functions. As a data-driven method, the underlying rules are not explicitly programmed into the network but learned from pattern recognition on the relationship between inputs and outputs of data. These rules are found by minimizing a loss function between the true outputs and the predicted outputs. It is optimized by adjusting the weights and biases of the network so that the loss function reaches a minimum.

A CNN is an end-to-end method, where the most meaningful spatial features are also learned from the data itself through convolution operations. These features are then combined for a classification task, producing the desired outcome based on the input.

In this work, we use an improved version of the CNN architecture described in Ferreira et al. (2020). This consists of a feed-forward network with an input image size of 128 × 128 pixels, where the number of convolutional blocks, convolutional layers, fully connected layers, number of filters, and kernel sizes, are all defined by the following hyperparameters:

• 1.  
  number_conv_blocks define the number of convolutional blocks, each will probe features of different scales;
• 2.  
  number_conv_per_block describe how many convolutions each block will have;
• 3.  
  initial_number_filters define the starting number of filters, which are then doubled after each convolutional block;
• 4.  
  initial_kernel_size is the initial size of the convolutional kernel, which is then reduced by two after each block, down to a minimum of three;
• 5.  
  n_fc_layers and size_fc_layers define the number of hidden layers and their respective sizes, respectively;
• 6.  
  l2_regularization and dropout are the degrees for each regularization technique used, respectively. l2 regularization is applied to all convolutional layers, where dropout is applied only after the hidden layers.

The approach of variable depth and width for neural networks is similar to the family of networks described in Tan & Le (2019). However, in our case, the networks are smaller due to the smaller image size used.

We modify the methods from Ferreira et al. (2020) to improve generalization of our models. First, instead of using two binary classification networks and combining their predictions to construct a multiclass classification, we now only use one network for the binary classification of post-mergers and star-forming galaxies.

Second, we treat the learning rate differently. In Ferreira et al. (2020), we monitored the learning rate decays during training as a hyperparameter. Here, we use cosine annealing, a type of learning rate scheduling (e.g., Loshchilov & Hutter 2016, for an explanation) combined with a regular stochastic gradient descent optimizer (Zhou et al. 2020). This approach probes several different learning rate regimes during training and uses cyclic resets that serve as a way to avoid unstable local minima, improving generalization of the solutions.

All hyperparameters are determined by a Bayesian optimization process (The GPyOpt 2016), and the values for the best model used here are summarized in Table 2. These values can be directly used in conjunction with our public keras implementation.

To avoid overfitting pitfalls from using our CANDELS background augmentation pipeline (Section 2.4), we train a suite of models, one for each CANDELS field. Because we have included areas of all of the CANDELS fields as background in our training set, the network could potentially memorize these and use them for predictions, impairing the results. To ensure this is not the case, each CANDELS field has two models—one at low redshift, 0.5 < z < 1.5, and one at high redshift, 1.5 < z < 3.0—trained only with images augmented with regions of the other four fields. All data sets (training, validation, and test) are restricted in this way, guaranteeing that any overfitting of the CANDELS background will have no impact on the final application of our models.

An example of this process is outlined in Figure 6, for the models that will be used for predictions in the GOODS-North (GDN) field. The training set contains galaxies augmented with the COSMOS (COS), GOODS-South (GDS), Extended Groth Strip (EGS), and The UltraDeep Survey (UDS) fields while the validation and test sets only contain galaxies from GDN.

Zoom In Zoom Out Reset image size

This ensures that each model is tailored to one CANDELS field and that no source from that particular field is used during training, i.e., the network never sees any of its data. We further apply a regularization method that makes use of random rotations and image flips on the fly during the training time.

Nonparametric structure measurements of galaxies are a traditional way to select galaxy mergers (Conselice et al. 2003; Lotz et al. 2004, 2008; Snyder et al. 2017). To measure structures for our sample, we fit Sérsic profiles to all galaxies, using the software Morfometryka (Ferrari et al. 2015; Albernaz Ferreira & Ferrari 2018; Lucatelli & Ferrari 2019). Morfometryka measures asymmetry (A), concentrations (C), the Gini coefficient (G), moment of light of the brightest pixels (M₂₀), normalized information entropy (H), and others. It also measures several structural parameters and fits 1D and 2D Sérsic profiles. For our purpose, we are particularly interested in the asymmetry of the galaxies (A), as well as their smoothness (S) since, together, they define a common criterion for finding galaxy mergers:

$$\begin{eqnarray*}&&(A\gt 0.35),\quad (A\gt S).\end{eqnarray*}$$

The asymmetry is defined as the pixelwise normalized difference between the original image and the same image rotated by 180°,

$$\begin{eqnarray*}&&A=\displaystyle \frac{{\rm{\Sigma }}| I-{I}_{180}| }{{\rm{\Sigma }}| I| }-{A}_{\mathrm{bg}},\end{eqnarray*}$$

where I is the image, I₁₈₀ is the rotated image, and A_bg is an asymmetry term associated with the background (e.g., Conselice 2014). We measure A_bg in each cell of a meshgrid overlaid onto the image, omitting the area occupied by the segmentation map of the central galaxy. We then use the median of these values as A_bg. This ensures a robust modeling of the impact of the background in the resulting asymmetry of the image (e.g., Tohill et al. 2021).

Finally, as we are especially interested in investigating the nature of the peculiar/irregular cases, we follow the hybrid method proposed by Bickley et al. (2021). We first filter out regular symmetric galaxies from the sample using the asymmetry (A). Instead of using the widely used cut for selecting mergers (A > 0.35), we choose a conservative selection of galaxies with,

$$\begin{eqnarray*}&&(A\gt 0.1).\end{eqnarray*}$$

This will remove cases that are irrelevant for our research question. These are galaxies without any disturbances that would classify them as peculiar or irregular.

In Figure 7 we show the distribution of asymmetries A measured with Morfometryka for star-forming galaxies (in blue) and post-mergers (in red) for the simulated galaxies. The distributions largely overlap, though asymmetries for post-mergers are generally slightly higher. The difference between both distributions is small enough that using solely the asymmetry (A > 0.35) will produce samples with low completeness and purity, and given that the fraction of merging galaxies is lower than regular star-forming galaxies, it is likely that this approach produces very contaminated samples.

Zoom In Zoom Out Reset image size

We measure the performance of our trained models in our prepared test sets. This is done by training the network with two realizations of the test data sets, one with full HST-matched properties including a CANDELS background patch of the sky (Sections 2.4, 3.2, which we call realistic mocks) and one with clean mocks with no sky noise and contamination included (which we call pristine mocks). For simplicity, in cases where we only mention the realistic mocks without specifying which CANDELS fields it was augmented with, we consider the average of all 20 models described in Section 3.2.

To compare between models and realizations of these data sets, we use traditional performance metrics common for evaluating machine-learning model performance. These consist of receiver operating characteristic (ROC) curves and precision-completeness diagrams (Powers 2011), as well as confusion matrices and their individual indices. Here, we are dealing with a single binary classification task, such that the probabilities of both classes respect the condition P(NMSF) + P(PM) = 1. Figure 8 displays the overall performance for each network.

Zoom In Zoom Out Reset image size

The left panel shows four different realizations of the network for comparison purposes. The network is trained twice to generate two different types of models: one labeled base that consists of a network trained with the realistic mocks, and a second labeled clean, which is trained with the pristine mocks. Then, each model is applied to both data sets. We do this to measure the best-case scenario within the simulations, in the absence of any contamination or impact from observational effects. Models trained with the realistic mocks data set are plotted as solid lines, while models trained with the pristine mocks data set are shown by dashed lines. Furthermore, the color conveys the data set in which the model was applied to, red and blue for pristine mocks and realistic mocks, respectively. In addition to these, a single parameter classifier based on the asymmetry (A) is also evaluated and displayed as the green dotted line. The area under the curve for each case can be found in the legend of the left panel.

The different realizations of our network (base and clean) cross-correlated with the realistic mocks and pristine mocks data sets confirm the importance of realistic observational modeling of the mocks (discussed in detail by Bottrell et al. 2019). This is especially important when crossing domains from cosmological simulations to real observations. Figure 8 shows that the base network performs just as well as the clean network when applied to the pristine mocks, resulting in similar performance metrics, as can be seen by the overlapping red curves. However, the base network outperforms by ∼10% the clean network when applied to the realistic mocks data set, as displayed by the difference between the blue curves in Figure 8. This demonstrates that correctly modeled observational features increase the generalization capabilities of the resulting models. A network that is only trained on pristine images will perform poorly in the real observations domain.

Importantly, all cases outperform the asymmetry by 20%–30%. To some extent, this is expected because asymmetries of post-mergers are lower than asymmetries of galaxies that are just in the beginning of their merging event, including cases of closely interacting galaxies. Evidently, the asymmetry function is a much more general morphological descriptor while the network is very specialized for the particular task of dividing post-mergers from star-forming galaxies.

We compare the performance of asymmetry (A) and CNN predictions further and show completeness-purity diagrams in the right panel of Figure 8. It displays outcomes for our ensemble of CNN models in the inferno color map, and for the classic asymmetry parameter in the viridis color map. The commonly used asymmetry value to classify galaxy mergers is generally higher than (A > 0.35), which is shown in the figure by the red patch over the curve. However, here we compare an asymmetry classifier with our neural network to exemplify how one can use the classification threshold of the network as a way to control the trade-off between precision and completeness. This is a useful feature when dealing with unbalanced data sets, like the case for galaxy mergers.

The precision and completeness of the asymmetry behave in unpredictable ways. First, the precision of the selection increases slowly, then it decreases again around (A ∼ 0.2), and spikes above 0.6 precision for (A > 0.8), but with very low completeness. We do not seek to redefine its use, but merely contrast it with our deep learning approach, and show in broad terms when it might fail when dealing with ambiguous morphologies.

Our network is able to correctly identify post-mergers and star-forming galaxies from the IllustrisTNG simulation in ∼80% of the cases. Figure 9 shows the confusion matrix for the realistic mocks data set identified within each individual CANDELS field, as well as for the pristine mocks sample, where accuracy reaches ∼90%. All classifications are done with the model trained with the realistic mocks. We show true positives (TPs) and true negatives (TNs) in blue, and false positives (FPs) and false negatives (FNs) in pink. The CLEAN case represents the best-case scenario, where our current method and data set achieves an even higher performance of ∼91% TPs. A histogram of the redshift distribution for each cell helps to visualize any possible biases in redshift for the misclassification cases. This demonstrates that the models are more likely to correctly classify low-redshift galaxies, as they represent the majority of the samples.

Zoom In Zoom Out Reset image size

4.1.1. Impact of Redshift

With the goal of applying our models to a wide range of redshifts, we explore how our performance metrics are impacted by increasing redshifts. Following the angular size–distance-relation, galaxies at increasingly larger distances from low to intermediate redshifts will be greatly impacted by decreasing resolution, which means that morphological features are less well sampled. The right panel of Figure 10 shows this effect on the performance of our models, where the scores of the metrics gradually decrease with increasing redshift, going from 85% accuracy at z = 0.5 to around 80% at z = 2. The error bars—sampled from bootstrapping our testing samples—follow accordingly.

Zoom In Zoom Out Reset image size

4.1.2. Contamination Impact on Classification

We use the contamination estimates measured in Section 2.5 to find the contamination failure threshold of our classifier, comparing performance metrics for subsets of the test set selected in bins of both the overlapping percentage, Θ, and the average background flux per pixel, BG_flux, as shown in Figure 10. The horizontal black dashed line at 0.8 shows the accuracy of the model when evaluated in the complete test set (80%). The metrics outperform this baseline in subsamples of images with low contamination, decreasing as we increase each of the contamination factors.

As described in Section 2.5, we select the point where the average mean values for each metric fall below the dashed line, which is our contamination cutoff, i.e.,

$$\begin{eqnarray*}&&{\rm{\Theta }}\sim 15 \% ,\quad {\mathrm{BG}}_{\mathrm{flux}}\sim {10}^{-1.5}\ {\rm{e}}\ {{\rm{s}}}^{-1}\ {\mathrm{pix}}^{-1}.\end{eqnarray*}$$

Since it is not possible to directly measure the contamination parameters in the real observations, we refer the reader to our deep learning model trained to measure the contamination in the Appendix.

We use our network to carry out predictions in all real CANDELS galaxies at 0.5 < z < 3, M_* ≥ 10^9.5 M_⊙, S/N > 50, and H_MAG < 24.5. We filter out regular galaxies using a conservative asymmetry cut of A > 0.1 as we are interested only in asymmetric, irregular/peculiar systems. This selection results in a sample of 23,494 galaxies, for which 14,410 have visual classifications from Kartaltepe et al. (2015). Based on the classifications from our networks, we separate these galaxies in post-mergers and noninteracting star-forming galaxies using a threshold probability of 60%. Galaxies with probabilities 50% < P(PM) ∧ P(SF) < 60% are not considered in any class. These represent 2125 galaxies (≈15%) of the sample with visual classifications. Figure 11 showcases some examples of galaxies in the CANDELS fields separated by the classification of our models. Post-mergers are shown in the left panel, and star-forming galaxies are shown in the right panel.

Zoom In Zoom Out Reset image size

To investigate how the relative number of post-mergers and star-forming galaxies changes over cosmic time, we divide the CANDELS sample in bins of Δz = 0.25. Figure 12 shows the change of class fractions change with redshift. We do this analysis in two mass regimes: low-mass galaxies with $9.5\lt \mathrm{log}({M}_{* }/{M}_{\odot })\lt 10.0$ (left panel) and high-mass systems with $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10.0$ (right panel). Noninteracting star-forming galaxies are shown by the blue circles while post-mergers are shown in by the red squares. The upper dashed line displays the fraction of galaxies that are not mergers, including the star-forming galaxies and other low-probability cases not included in any class.

Zoom In Zoom Out Reset image size

For the lower-mass post-mergers, we see an upward trend from ∼15% at z = 0.5 to ∼35% at z = 2, then a slight decrease beyond z = 2. This is still consistent with a ∼35% fraction within the error bars. The star-forming galaxies behave in the opposite way, decreasing from ∼70% at z = 0.5 to around ∼55% at z = 2. This suggests that among asymmetric galaxies of this mass range, there is an exchange between the classes as we go to higher redshifts up to z = 2. This once again emphasizes that classifications of local galaxies that are purely based on the asymmetry (A) are highly contaminated with noninteracting star-forming galaxies. However, this is mitigated at higher redshifts where we find more post-mergers. Nevertheless, samples selected based on A are still dominated by star-forming galaxies, albeit to a lesser extent.

Trends for higher-mass galaxies are substantially different (right panel of Figure 12). While the post-mergers exhibit a similar but steeper upward trend from ∼20% at z = 0.5 fraction to ∼50% at z = 2, the relative fraction of star-forming galaxies shows a constant value of ∼50% at 0.5 < z < 3, while the fraction of the rest of the sample (dashed line) goes from ∼75% to ∼50%. This supports the idea that for more massive systems, post-mergers at higher redshifts will eventually become massive passive galaxies with no significant asymmetric features.

At the highest redshift bins, the error bars are large, an indication for our networks to perform less accurately above z = 2. The fraction of post-mergers changes from 30% to around 50% at z = 2. We therefore attribute the downtrend in post-mergers beyond z = 2 to the poor performance of our models at high redshifts and do not take this to imply a real evolutionary effect.

We know that mergers are more common in the past (e.g., Mundy et al. 2017; Duncan et al. 2019; Ferreira et al. 2020; Whitney et al. 2021), and here we find further evidence that this is also the case for peculiar galaxies, indicating that the nature behind these disturbed morphologies at earlier times can be attributed to merging. To further investigate this, we select all galaxies from Kartaltepe et al. (2015) that are classified as an irregular/peculiar with f_Irr > 0.75, i.e., cases where more than 50% of the visual classifiers agree on the classification, and check how our networks perform on this subset. We observe similar trends with redshift, with the fraction of post-mergers increasing by ∼20% from z = 0.5 (∼30%) to z = 2 (∼50%), which agrees with the results for the complete sample. Furthermore, our methods classify ∼50% of the galaxies visually classified as potential mergers (f_merger > 0.75) in Kartaltepe et al. (2015) as post-mergers. This is higher than random, but does show the difficulty of obtaining exact matches between mergers determined visually compared with a quantitative process.

4.2.1. Visual Representation of the Classification

As a way to visualize how our networks organize the features extracted from the images to produce the final classification, we generate a 2D representation of the final dense layer of the network corresponding to 128 neurons (128 dimensions) using a UMAP (Figure 13). The color code of the points expresses their respective labels, red for post-mergers and blue for star-forming galaxies. Then, we overplot the positions assigned by the network for unlabeled CANDELS galaxies. We also include some examples of images of CANDELS galaxies close to their original position in this manifold as a way to visualize how the morphologies change with its position. Each region in the parameter space of this diagram is directly related to a probability. The maximum probability is found in the extreme regions farther away from the center, which represents how different these objects are for the network. Images of galaxies that the network struggles to identify are mixed in the bottom middle, representing the region where both probabilities are similar P(PM) ∼ P(SF).

Zoom In Zoom Out Reset image size

The influence of merging on the structure of peculiar/irregular galaxies at intermediate redshifts (0.5 < z < 3.0) is directly related to the question of whether merging galaxies can induce more starbursting episodes than galaxies evolving secularly. Enhanced star formation can then lead to more clumpy and asymmetric structures, and thus can impact the morphological appearance of galaxies greatly. By examining the SFMS of galaxies, one can investigate the nature of galaxies with unusually high SFRs and the formation path that resulted in this physical effect.

In order to investigate this, we select only galaxies in our CANDELS fields sample that lie above the SFMS as parameterized by Schreiber et al. (2015). We separate these sources by stellar masses, redshifts, and their post-merger/star-forming classification, measuring the mean distance to the SFMS (Δ_MS), as:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{MS}}=\mathrm{log}(\mathrm{SFR})-\mathrm{log}({\mathrm{SFR}}_{\mathrm{MS}}),\end{eqnarray} \tag{ 7 }$$

where $\mathrm{log}(\mathrm{SFR})$ is the log star formation rate of a particular galaxy, and $\mathrm{log}({\mathrm{SFR}}_{\mathrm{MS}})$ is the parameterization from Schreiber et al. (2015). The SFRs and stellar masses used here for CANDELS galaxies were compiled by Duncan et al. (2019) through SED fitting. We refer the reader to this publication for further details. In Figure 14 we show the mean value of each stellar mass bin, for four redshift ranges (one in each panel), separated into star-forming and post-merger galaxies by our classifications. For the 0.5 < z < 2.5 redshift range (panels A, B, and C), we do not find any impactful difference between the classes and Δ_MS, with all offsets well within the error bars. However, for redshifts 2.5 < z < 3.0, post-mergers with $\mathrm{log}({M}_{* }/{M}_{\odot })\lt 10.0$ are on average ∼0.1 dex higher than star-forming galaxies of the same mass. The opposite is found for $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10.0$; however, uncertainty is higher here. Additionally, Δ_MS increases with redshift in all cases, which describes a larger scatter above the SFMS. However, given the performance metrics of our models at high redshift (Figure 10), we can not claim that this is a real effect. We stress that in Figure 14 we only select galaxies above the SFMS, which is why the distance is always positive.

Zoom In Zoom Out Reset image size

In summary, locations of post-mergers and noninteracting galaxies in the SFMS diagram are comparable, with the possible exception at the highest redshifts. This suggests one of the following: within our sample of CANDELS galaxies, major-merging is not playing a major role in enhancing starbursting episodes; or the timescale probed by our method is too large and the SFR enhancement from the captured post-mergers is short lived.

A relevant result was discussed in Hani et al. (2020), who investigated TNG300-1 post-mergers at 0.0 < z < 1.0. They showed that post-mergers have enhanced sSFRs by a factor of ∼2, but that this effect decays on timescales of ∼0.5 Gyr, which can be driven in part by minor mergers. Although we do not find evidence for an enhancement in starbursts due to major mergers, we do not rule out the importance of minor mergers to this effect. We trained our models without the presence of minor mergers, but we can not be sure that the star-forming galaxies classified by our models are not in some cases triggered by minor mergers.

Our deep learning classifications relate to two different formation pathways. These formation scenarios could result in structures that differ for post-mergers and star-forming galaxies. To verify if in fact their structures are diverse from one another, we investigate light profile fitting by using Sérsic profiles measured by Morfometryka.

Figure 15 shows the distribution of Sérsic indices for post-mergers in red, and star-forming galaxies in blue. In general, each class presents very distinct distributions: the post-mergers have a mean Sérsic index $n\sim {1.8}_{-0.6}^{+0.7}$ roughly representative of a transition from disks to spheroids; star-forming galaxies have systematically lower Sérsic indexes with $n\sim {1.1}_{-0.5}^{+0.5}$, which is more consistent with disk-dominated galaxies. This offset of ∼1 dex increases for classification thresholds at higher values. The average Sérsic profile (n) of post-mergers increases while the distribution for star-forming galaxies continues with a similar shape. This is quantitative evidence that (1) post-mergers with higher light concentrations are more easily separable from noninteracting star-forming galaxies, and (2) these types of galaxies are intrinsically different from each other.

Zoom In Zoom Out Reset image size

By using the new classifications from this work, we can update classifications from Ferreira et al. (2020) for cases where the previous method had ambiguous probabilities for some major mergers and nonmergers.

Our new data set accounts for the effects of dust; it is not limited by orientation and probes the rest-frame optical. Thus we can check if any major-merger classifications in the previous work can be attributed to noninteracting star-forming galaxies or if any nonmergers can be reclassified as post-mergers. This is done by comparing the probabilities for major mergers and nonmergers, P(MM) and P(NM), respectively, from Ferreira et al. (2020) to the new probabilities P(PM) and P(SF). We update a nonmerger classification to post-merger if

$$\begin{eqnarray*}&&P(\mathrm{PM})\gt P(\mathrm{NM}),\end{eqnarray*}$$

and update the major-merger classifications to nonmerger if

$$\begin{eqnarray*}&&P(\mathrm{SF})\gt P(\mathrm{MM}).\end{eqnarray*}$$

In other words, we reclassify galaxies from the previous sample where our new method is more certain about its classification than the previous one. This leads to ∼5% of major mergers reclassified as star-forming nonmergers, which lowers the overall merger fractions at lower redshifts and keeps it similar at higher redshifts. In Figure 16 we compare the new merger fraction measurements, in green, to the results from Ferreira et al. (2020), in gray.

Zoom In Zoom Out Reset image size

The updated fit of the cosmic evolution of the merger fraction, f_m (z)

$$\begin{eqnarray}&&{f}_{m}{(z)=0.011\pm 0.002\times (1+z)}^{2.71\pm 0.31},\end{eqnarray} \tag{ 8 }$$

with errors estimated with bootstrapping, agrees with the previous measurement in Ferreira et al. (2020) within errors. To measure the galaxy major-merger rate (${ \mathcal R }$), we combine the timescale (τ_m = 0.5 Gyr) used in our selection (Section 2.1) with this merger fraction through

$$\begin{eqnarray}&&{ \mathcal R }=\displaystyle \frac{{f}_{m}}{{\tau }_{m}}.\end{eqnarray} \tag{ 9 }$$

The updated galaxy major-merger rate is

$$\begin{eqnarray}&&{ \mathcal R }=0.022\pm 0.006\times {(1+z)}^{2.71\pm 0.31}.\end{eqnarray} \tag{ 10 }$$

We emphasize that this correction is a minor adjustment to the galaxy major-merger rates presented in Ferreira et al. (2020), which remain broadly consistent with each other.

We now investigate the possible contamination in merger samples that are selected through our method. This approach is fairly direct and based on Bayesian statistics, and relies on some understanding of the true intrinsic merger fraction and how it evolves with time. It also requires that we have a good understanding of the fraction of contamination in merger samples (Bickley et al. 2021). The basic Bayesian formula to understand this is given by the following:

$$\begin{eqnarray}&&P({\boldsymbol{M}}| {\boldsymbol{S}})=P({\boldsymbol{S}}| {\boldsymbol{M}})\times \displaystyle \frac{P\left({\boldsymbol{M}}\right)}{P\left({\boldsymbol{S}}\right)},\end{eqnarray} \tag{ 11 }$$

where P( M ∣ S ) is the probability of a merger, given that a method used to select mergers, ( S ), identifies it as such. The value of P( M ) is the probability that an object is a merger before a selection of merger is made. P( S ) is the probability that a galaxy is selected as a merger, whether a real merger or a false positive. Because of the results of this paper, we know that this last number is very likely not equal to unity. It in fact can depend on various factors and methods of finding mergers. We can write the probability P( S ) as:

$$\begin{eqnarray}&&P\left({\boldsymbol{S}}\right)=P({\boldsymbol{S}}| {\boldsymbol{M}})\times P({\boldsymbol{M}})+P({\boldsymbol{S}}| {\boldsymbol{NM}})\times P(\mathrm{NM})\end{eqnarray} \tag{ 12 }$$

where NM standard for nonmergers, where P( S ∣NM) is the probability of identifying correctly a nonmerger, and the value of P(NM) is the probability that the galaxy is not a merger. We can simplify this if we know, a priori, what the merger fraction is based on previous work. If we denote the merger fraction as f_m, and the machine-learning probability of finding a merger/nonmerger as p_m and p_mn, respectively, then we can rewrite Equation (11), as

$$\begin{eqnarray}&&P({\boldsymbol{M}}| {\boldsymbol{S}})=\displaystyle \frac{{p}_{{\rm{m}}}\times {f}_{{\rm{m}}}}{({p}_{{\rm{m}}}{f}_{{\rm{m}}}+(1-{p}_{\mathrm{nm}})(1-{f}_{m}))}.\end{eqnarray} \tag{ 13 }$$

Thus, for example, if the accuracy of a machine-learning method for finding a merger is 0.9 and the accuracy for finding a nonmerger is 0.9, and the merger fraction f_m = 0.1, then the probability that a galaxy identified as a merger is actually a merger is P( M ∣ S ) = 0.5. This implies that even when the accuracy of finding mergers and nonmergers is 90%, at the lowest redshifts, where the merger fraction is low ∼10%, there is still a 60% chance that an identified merger is identified incorrectly as such. At higher redshifts, where the intrinsic merger fraction is higher, the probability of finding a merger correctly increases to ∼80% when the merger fraction is as high as f_m ∼ 0.3.

We can generalize the equation for P( M ∣ S ), as a function of z, by considering how the merger fraction f_m evolves with redshift, such that:

$$\begin{eqnarray*}&&P({\boldsymbol{M}}| {\boldsymbol{S}},z)=\displaystyle \frac{{p}_{{\rm{m}}}(z)\times {f}_{m}(z)}{{p}_{{\rm{m}}}(z){f}_{m}(z)+(1-{p}_{\mathrm{nm}}(z))(1-{f}_{m}(z))},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{f}_{m}{(z)\equiv {f}_{0}(1+z)}^{m},\end{eqnarray*}$$

which gives us a tool to understand how our classifications might be contaminated by sample unbalance effects with respect to redshift.

From this we can conclude that a significant fraction of individual galaxies within the CANDELS imaging may be incorrectly identified as either mergers or nonmergers. From our results here, our method effectiveness for correctly classifying mergers increases from ∼40% at z ∼ 0.5 to ∼70% at z ∼ 3. This is likely what can account for some of our misidentified galaxies as discussed in Section 4 when discussing the success of our method of separating star-forming systems from those that are undergoing mergers.

These are conservative estimates that do not include the fact that we pre-select CANDELS galaxies based on their asymmetry. This should increase P( M ∣ S , z) further since f_m is higher among galaxies with A > 0.1.

There is a growing concern on the applicability of simulation-trained deep learning models when applied to an intrinsically different domain. For us, this is the case with going from cosmological simulations to real observations. When transferring from one domain to another, deep learning models might fail due to relying too much on domain-specific features. Several techniques were developed to address this problem, focused on forcing neural networks to learn domain invariant features, leading to more robust models. Ćiprijanović et al. (2021) showed that adopting techniques for domain adaptation could increase model performance when applying to the target domain by 20%.

In our case, the source domain comprises the IllustrisTNG galaxies and the target domain comprises the CANDELS observations. To check if we need to apply domain adaptation techniques to this particular problem, we used uniform manifold approximation and projections (UMAPs; as described in Section 4.2.1; McInnes et al. 2018) to reduce the high dimensional space generated by the features extracted by our network to a 2D-space that is easy to visualize. ⁵ Then, for each of our trained models, we compare whether the features extracted by the network show similar distributions for Illustris and CANDELS galaxies. In Figure 17 we show UMAPs for each of the CANDELS field models, for low redshift (left) and high redshift (right), color-coded by their class in the case of Illustris and in black for real CANDELS galaxy images. As can be seen, these distributions of simulated galaxies and real observations are clustered together, with very few outliers not following the main cluster. Additionally, we can see that each class forms its own cluster, with overlapping regions, showing that features between classes are distinct and in general not domain specific. ⁶

Zoom In Zoom Out Reset image size

We attribute the generalization success of our models to our mock data pipeline, which is tailored to mimic each individual CANDELS field with maximum fidelity—with their instrumental and observational features. Also augmentations with patches of the sky from CANDELS introduce real observations into our source domain, which not only make our training sets big (∼140.000 images) but also help with domain confusion within the network. Thus, we do not include any domain adaptation process in our pipeline.

Based on the contamination measurements described in Section 2.5, we devise a new neural network with the goal of predicting the overlapping percentage (Θ) and the background flux (BG_flux) measurements from real observations. In this way, contamination thresholds can be applied to real observational samples in a similar way to what is done in the simulations.

The contamination quantification depends on our ability to separate the background patch of the sky from the central source, a feature that is only available when we are post-processing simulated galaxies. In the case of real CANDELS observations, directly measuring these properties is difficult, because it is not straightforward to de-blend background/foreground sources if they are projected on top of one another or are close enough to be a potential interaction.

We use all of the contamination information from our data pipeline (Section 2.4) to train a neural network to predict these values from the final image, without separating source and background. We use the same network architecture described in this work, but we replace the final sigmoid layer with a linear activation function, and change the loss function as well. The result is a model that can be directly applied to real observations, where the image is the input and the outputs are values for Θ and BG_flux.

Figure 18 displays the performance of these predictions based on the original measurements, together with Pearson and Spearman correlation indices. In general, the performance of the model is in good agreement with the original measurements, with rms errors on the order of ∼ 10⁻³ for BG_flux and ∼5% for Θ. These limits are well within the region of the parameter space formed by these indices that we defined as a low-contamination region. Apart from the small bias making the predictions undervalue the truth values, the performance is good enough to separate high-contamination cases from the rest of the sample, which is ultimately our goal. In Figure 19 we show examples of different combinations of Θ and BG_flux.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even though this network is designed to be used within the context of this work as a way to reproduce contamination quantification in the same manner as what was done with the simulated images, we recognize that this can be useful for a wider application. For example, this can be used as a fast selection tool that can remove catastrophically bad cases from big samples in just a couple of seconds; thus, it can be a powerful tool for quick exploration. In this regard, we release this model independent of the classification models presented in Section 3.1.