# A Simulation Driven Deep Learning Approach for Separating Mergers and Star
Forming Galaxies: The Formation Histories of Clumpy Galaxies in all the
CANDELS Fields

Leonardo Ferreira Centre for Astronomy and Particle Physics, University of
Nottingham, Nottingham, UK  Christopher J. Conselice  Centre for Astronomy and
Particle Physics, University of Nottingham, Nottingham, UK  Jodrell Bank
Centre for Astrophysics, University of Manchester, Oxford Road, Manchester UK
Alan Turing Fellow, Alan Turing Institute, 96 Euston Rd, London NW1 2DB, UK
Ulrike Kuchner  Centre for Astronomy and Particle Physics, University of
Nottingham, Nottingham, UK  Clár-Bríd Tohill  Centre for Astronomy and
Particle Physics, University of Nottingham, Nottingham, UK

(Received January 31, 2022; Revised March 20, 2022)

###### Abstract

Being able to distinguish between galaxies that have recently undergone major
merger events, or are experiencing intense star formation, is crucial for
making progress in our understanding of the formation and evolution of
galaxies. As such, we have developed a machine learning framework based on a
convolutional neural network (CNN) to separate star forming galaxies from
post-mergers using a dataset of 160,000 simulated images from IllustrisTNG100
that resemble observed deep imaging of galaxies with Hubble. We improve upon
previous methods of machine learning with imaging by developing a new approach
to deal with the complexities of contamination from neighbouring sources in
crowded fields and define a quality control limit based on overlapping sources
and background flux. Our pipeline successfully separates post-mergers from
star forming galaxies in IllustrisTNG 80% of the time, which is an improvement
by at least 25% in comparison to a classification using the asymmetry (A) of
the galaxy. Compared with measured Sérsic profiles, we show that star forming
galaxies in the CANDELS fields are predominantly disc-dominated systems while
post-mergers show distributions of transitioning discs to bulge-dominated
galaxies. With these new measurements, we trace the rate of post-mergers among
asymmetric galaxies in the universe finding an increase from 20% at z=0.5 to
50% at z=2. Additionally, we do not find strong evidence that the scattering
above the Star Forming Main Sequence (SFMS) can be attributed to major post-
mergers. Finally, we use our new approach to update our previous measurements
of galaxy merger rates $\mathcal{R}=0.022\pm 0.006\times(1+z)^{2.71\pm 0.31}$

Galaxy mergers (608) — Convolutional neural networks (1938)

††journal: ApJ

##  1 Introduction

The first deep Hubble Space Telescope (HST) images of the distant universe
revealed that many distant and faint galaxies are in fact irregular/peculiar
in appearance (e.g., Williams et al., 1996). Because the first cameras on HST,
WFPC1/WFPC2 were sensitive in optical wavelengths only, probing distant
galaxies was limited to their rest-frame ultra-violet light, due to the
effects of redshift. It was unclear whether the peculiar appearances were the
result of observational limitations or real. The question thus remained
whether the observed irregularities were in fact just the star forming areas
of these galaxies, while the older stars remained below detection. When the
NICMOS camera was launched in 1998 on HST it became clear that the
morphologies of distant galaxies were peculiar in their rest-frame optical
wavelengths as well, implying that the bulk stellar mass in these galaxies was
indeed out of equilibrium (e.g., Dickinson et al., 2000; Papovich et al.,
2005; Conselice et al., 2005; Mortlock et al., 2013; Whitney et al., 2021).
The common census was that distant galaxies are indeed intrinsically peculiar.
However, it remained unclear why and how this finding relates to the various
possible modes that could be responsible for producing these irregularities in
galaxy structures at high redshifts. The peculiar appearance is likely linked
to the formation process of the galaxies, but details of the origin of the
observed irregular structures have proven difficult to fully understand.

Since then, it has become clear that, overall, galaxies gradually transition
from peculiar galaxies at higher redshifts to ellipticals and disc systems at
lower redshifts (e.g., Conselice et al., 2003; Lotz et al., 2004; Mortlock et
al., 2013; Huertas-Company et al., 2015). This conclusion was made possible
with the advent of the WFC3 camera on HST which allowed astronomers to trace
the morphological evolution of galaxies over large areas of the sky. Galaxies
are therefore undergoing a transformation, and their irregular origins reveal
clues about the processes which drive galaxy formation. One popular and well
explored hypothesis is that these systems are in fact undergoing hierarchical
mergers to form larger systems. The basic idea is that two galaxies in the
early universe smash together to form a larger galaxy, a process which is
predicted to be a critical element in the cosmological context of galaxy
formation within a Cold Dark Matter (CDM) universe, with well defined
predictions of this process (e.g., Jogee et al., 2009; Bertone & Conselice,
2009; Mundy et al., 2017).

To make progress in understanding the evolution of galaxies, it is crucial to
identify merging galaxies correctly. In order to separate galaxies into
mergers and non-mergers – initially focusing on the nearby universe –
quantitative morphology tools were developed that use parameters such as the
asymmetry index ($A$) (e.g., Conselice et al., 2000; Conselice, 2003). Merging
galaxies are often identified through a combination of these morphology
measurements such as the CAS parameters, which links measurements of the
Concentration, Asymmetry, and Smoothness (e.g., Conselice et al., 2003; Lotz
et al., 2004). However, mergers are not uniquely identifiable in this
parameter space and some do not fall into the selection criteria at all (e.g.,
Conselice, 2006; Lotz et al., 2008). Therefore, care has to be taken to
calibrate their usage.

Today, the merger rate can be accurately measured to high redshifts ($z\sim
3$) using galaxy structure (e.g., Conselice et al., 2003, 2008; Man et al.,
2016; Mantha et al., 2018; Ferreira et al., 2020; Whitney et al., 2021). Using
e.g., CAS parameters, the measurements show that the merger rate increases at
higher redshifts up to $z\sim 3$, such that f${}_{\rm merger}\sim(1+z)^{2-3}$
(e.g., Conselice, 2014), an evolution which scales similarly to the density of
the universe, which evolves as $\sim(1+z)^{3}$. This implies that with
identifications of mergers at both high and low redshifts, we are able to
trace the galaxy merger history and investigate the role of mergers within the
formation of galaxies over time (e.g. Conselice, 2006; Mundy et al., 2017).

In addition to high merger rates, distant galaxies have much higher star
formation rates than today, peaking at $z\sim 2$ (e.g., Madau & Dickinson,
2014). We further know that galaxy structure is highly dependent on the star
formation rate in the sense that intensely star forming galaxies generally
appear more clumpy and irregular than quiescent galaxies at all redshifts
(e.g., Windhorst et al., 2002; Mager et al., 2018; Gao et al., 2018; Sazonova
et al., 2021). In fact, these two different types of galaxies – mergers and
non interacting intensely star forming galaxies – can look very similar by
eye, which complicates visual classifications. Even kinematically it can be
challenging to distinguish mergers from rotating galaxies with high
dispersions (e.g., Simons et al., 2019; Bottrell et al., 2021). In addition,
classifications and selections of galaxies after a merger event (post-mergers)
are highly contaminated by misclassified isolated galaxies with high specific
star formation rates (sSFR). This is because their star forming regions and
dusty inter stellar medium (ISM) can generate asymmetric features reminiscent
of (post-) merger features. It is therefore currently unknown if and how we
can correctly distinguish whether a galaxy is undergoing intense star
formation, or some type of merger using galaxy structures and morphologies.

One way to approach this question is through novel techniques using machine
learning. Recently, tremendous progress has been made in applying supervised
deep learning methods to investigate galaxy morphology (e.g. Huertas-Company
et al., 2018; Reiman & Göhre, 2019; Huertas-Company et al., 2019a; Cheng et
al., 2019; Martin et al., 2020; Walmsley et al., 2020, 2021). These end-to-end
techniques are also very promising for investigating galaxy mergers
specifically (Ackermann et al., 2018; Pearson et al., 2019a, b; Bottrell et
al., 2019; Wang et al., 2020; Ferreira et al., 2020; Bickley et al., 2021).
Additionally, one can also leverage information not only from visual
classifications and observations, but also by forward modeling cosmological
simulations to the observational domain (Ćiprijanović et al., 2020, 2021).

We have recently started a machine learning exercise to determine the merger
history of galaxies using cosmological simulation runs from IllustrisTNG
(Vogelsberger et al., 2014; Pillepich et al., 2018a; Nelson et al., 2019). In
Ferreira et al. (2020) we were able to separate mergers from other types of
galaxies in IllustrisTNG to a success rate of 90% up to $z\sim 3$. The present
paper is a followup to our first paper, in which we now investigate whether it
is possible to distinguish merging galaxies from intensely star forming
galaxies (above the SFR-stellar mass main sequence). These galaxies have the
lowest success rates in classifications from Ferreira et al. (2020). Our task
in this paper is to correctly distinguish mergers from star forming galaxies
by only using their morphology and structure.

This paper is organized as follows: in §2 we describe the data sets we
constructed for this task, from IllustrisTNG (simulations) and CANDELS
(observations). A description of the methods we used to train a Deep Learning
model and how we measure the structure of the galaxies in our samples is given
in §3. We present our results in §4 while a discussion on the implications is
laid out in §5. Finally, we summarize and conclude our findings in §6.

##  2 Data

To test our new Deep Learning approach, we use simulated galaxies from
cosmological simulations post-processed with the SKIRT (Camps & Baes, 2015,
2020) dusty radiative transfer code. The simulations are based on IllustrisTNG
(Sec. §2.1) which are used for the construction of the training sample for a
Convolutional Neural Network (CNN) that is subsequently applied to observed
galaxies from the CANDELS fields (Sec. §2.3). Our sample definitions for post-
mergers and star forming galaxies are given in Sec. §2.2. We discuss the
pipeline used to generate CANDELIZED mock images from IllustrisTNG in Sec.
§2.4.

###  2.1 IllustrisTNG

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,300\rm\ Mpc\ h^{-1}$ length size, aptly named TNG50, TNG100
and TNG300 (Pillepich et al., 2018b; Naiman et al., 2018; Nelson et al., 2018;
Springel et al., 2018; Marinacci et al., 2018; Nelson et al., 2019; Pillepich
et al., 2019).

Our analysis makes use of the TNG100-1 simulation, which has proven to be a
good compromise between resolution and volume111We tested our selections on
TNG50-1 but resulting samples are too small for Deep Learning training..
TNG100 has already been used extensively in studies that analyze galaxy
morphologies and structures, including the comparison between simulations and
observations (Huertas-Company et al., 2019b; 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) show that TNG100 galaxies
reproduce observed objects well, especially disc-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_{\odot}$. 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, $\epsilon$, is not a
limitation when compared to the HST 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_{\odot}$, 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

$\frac{M_{\rm DM}}{M_{\rm total}}>0.1,$ (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 $\approx 2\%$ of galaxies from
the overall pool of available sources (subhalos) in TNG100.

We also remove objects that are smaller than the ACS 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 kpc, to a pixel scale
based on the cosmological model adopted by IllustrisTNG and the ACS pixel
scale

$R_{1/2Mass,pix}=\frac{a(z)}{h}R_{1/2Mass,kpc},$ (2)

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

###  2.2 Sample Definitions

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

$\mu=\frac{M_{2}}{M_{1}},\quad\mu\geq 0.25,$ (3)

where $M_{1}$, $M_{2}$ are the stellar masses of the galaxy pair involved in
the merging event, ranked by their stellar mass respectively, with
$M_{1}>M_{2}$. Galaxies are considered post-mergers if they have coalesced
into a single galaxy in the past $500$ Myrs, 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 are 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, non-interacting star forming galaxies are defined here as galaxies
that have sSFRs above the following threshold,

$sSFR>10^{-9.5}\rm\ yr^{-1}$ (4)

and are not interacting with other galaxies. To isolate non-interacting cases,
we exclude any galaxy from the simulation that had major or minor merger
events ($\mu>0.1$) around $\pm 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 non-interacting cases to extreme
starbursting episodes alone, but to select non-interacting 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 $\sim 6,000$ post-mergers
and $\sim 110,000$ non-interacting star forming galaxies. While this may be a
realistic representation of actual fractions (only $\sim 5\%$ of the sample
are post-mergers), training the network requires a balanced dataset. 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 non-
interacting 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 $\sim 4,000$ 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 are 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, $\tau$, and the mass ratio, respectively, for post-mergers. The
nature of the distribution for $\tau$ arises from the average time between
snapshots in the simulation of around $\sim 0.15\ \rm Gyr$. This timescale
represents 1 to 3 snapshots after the coalescence of stellar masses.

Redshift | Post-Mergers | Star Forming | Total
---|---|---|---
$0.5\leq z<1.0$ | 1214 | 1167 | 2381
$1.0\leq z<1.5$ | 1082 | 1140 | 2222
$1.5\leq z<2.0$ | 847 | 881 | 1728
$2.0\leq z<2.5$ | 589 | 556 | 1145
$2.5\leq z<3.0$ | 333 | 321 | 645
Table 1: Summary of the initial IllustrisTNG sample. The numbers in this table
present the sample before each galaxy was post-processed with SKIRT and
CANDELIZER (see text for detail), during which each image was augmented by 20
for 4 orientations and 5 different fields. Figure 1: Physical properties of
the 8,000 IllustrisTNG TNG100-1 selected simulated galaxies. For both types of
galaxies we show distributions for redshifts (top left), star formation rates
(top middle), stellar masses (bottom left) and stellar half-mass radius
(bottom middle) in red for post-mergers, and blue for star forming galaxies.
Distributions agree in general, with a small excess of stellar mass and size
for the post-mergers. The time since the last major merging event and the mass
ratio, $\mu$, – properties unique to the post-mergers – are shown in top right
and bottom right, respectively.

###  2.3 CANDELS Fields

One of the main goals of this work is to predict star forming and post-mergers
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, 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\sim 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, 2018b, 2018a, 2019) from SED fitting, essential for
creating matched samples. Stellar masses are estimated within a $0.2\rm\ 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 catalogue and in the
Kartaltepe et al. (2015) catalogues 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\rm\ mag$. Because this cut can
bias our sample against extended sources, we also include a signal-to-noise
($SNR$) lower limit of $\rm SNR>50$ to exclude any compact source with only a
few bright pixels. This magnitude cut removes 1074 sources, while the SNR 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_{\odot}$.
We apply a final cut using the asymmetry ($A>0.1$, Sec. §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 the
CANDELS fields combined.

Finally, we produce cutouts for $I_{814}$, $J_{125}$ and $H_{160W}$ bands
centering on each selected CANDELS galaxy, each with a field of view of $50\
\rm kpc\times 50\rm\ 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.

###  2.4 Pipeline to produce CANDELIZED Mocks

Figure 2: Example of the processing steps of our mock pipeline. a) Noiseless
F814W broadband image generated from the simulated galaxy datacube with 0.03
${}^{\prime\prime}/pix$ pixel scale. b) The same image after rebinning from
z=0.5 to z=0.6. c) Image convolved by the HST F814W PSF. d) Image with
Gaussian noise added. e) Image added on top of a random patch of the sky
within a CANDELS field with no neighbouring sources. f) Image added randomly
to a patch of sky with other sources in the field of view. As this patch of
the sky is randomly selected, all final images have varying levels of
contamination from nearby sources. We quantify this by the total flux in the
sky patch before adding the simulated source to it.

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 are shown in
Fig. 2 and are 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, that 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).

To create mock broadband images, we thus process each stellar particle with a
population synthesis model following the recipes from Trayford et al. (2017) &
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 (SSP) 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, firstly, they are
implemented in SKIRT and, secondly, they had been tested previously in similar
pipelines to generate mock observations from cosmological simulations
(Trayford et al., 2017; Rodriguez-Gomez et al., 2018). Finally, Bruzual &
Charlot (2003) are also the templates used to derive stellar masses and star
formation rates for all the CANDELS fields in Duncan et al. (2019) which 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^{6}$ photon
packets, assuming isotropic emission until they reach the virtual detector.

Figure 3: A random selection of IllustrisTNG simulated galaxies in our test
sample; post-mergers (left) and star forming galaxies (right), with their
redshifts, SFRs, and stellar masses printed in each stamp. Images are ordered
from left to right in redshift, and top to bottom in SFR. For post-mergers we
also display the time since merger, $T_{m}$, and the mass ratio $\mu$. All
stamps use a square-root normalization. Figure 4: Demonstration of the
augmentation pipeline for one random galaxy from TNG100-1 (ID=192802, z=0.55,
at different orientations). We increase our sample by augmenting the dataset,
reproducing it in four orientations (rows) in each of the CANDELS fields
(columns). The simulated galaxy is placed in a random patch of the sky in the
CANDELS fields and thus can have other sources in the final cutout. The amount
of contamination from neighbouring sources varies widely due to the random
sampling of the background described in §2.4. This contamination is quantified
by the overlapping percentage, $\Theta$, and the average flux of the
background patch, $BG_{flux}$

This process produces IFU datacubes 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 datacubes is located at a
distance of 10 Mpc (initial redshift $z_{0}$) 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

$\frac{f_{t}}{f_{o}}=(1+z_{t})^{-1}\left(\frac{D_{Lo}}{D_{Lt}}\right)^{2},$
(5)

due to cosmological dimming (Hogg, 1999, eq. 15). Next, we convolve the IFU
data with the broadband filters response functions for $I_{814}$, $J_{125}$
and $H_{160W}$. The results are clean, noiseless images from the simulation
galaxies at $30{\rm\ mas/pix}$ (matching the ACS pixel scale) before adding
any PSF effects (Fig. 2, images a) to c). We rebin $J$ and $H$ bands from $30\
{\rm mas/pix}$ to the WFC3 images pixel scales of $60\ {\rm mas/pix}$.
Examples for stamps where the background was added can be seen in Fig. 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.

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 $\sim 4000$ examples of each
class (Sec. 2.1). Fortunately, in the case of galaxy images there are ways to
increase the initial dataset 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 dataset in three ways outlined below. An example of this
approach is shown in Fig. 4, following the same galaxy in each possible
combination of orientation/field.

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.

Secondly, each CANDELS field has unique observational properties (e.g.
different noise levels, 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 RA and
DEC within each field, and make a large cutout of the area that is 4 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
RA and DEC 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 above criteria is found. This, combined with all the orientations,
augments our data set 20 times. In addition, this also helps the network to
generalize the impact of contamination from neighbouring 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 zoom-in/zoom-outs 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 $\sim 8,000$ examples to $\sim 160,000$.
However, having multiples of similar galaxies in our dataset 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 datasets.

Figure 5: Contamination characterization for $162,000$ IllustrisTNG simulated
images in our sample. We show the logarithm of the average flux per pixel of
the background measured in each cutout, $\log(\rm BG_{flux})$ vs. the
overlapping percentage, $\Theta$, which indicates how much the central galaxy
segmentation map is covered by the segmentation map of the sources in the
background. We define a conservative region of low contamination shown by the
dashed line and blue area, which contains $90\%$ of the whole sample. Every
point represents at least one image.

###  2.5 Contamination Quantification

Providing realistic levels of contamination and the inclusion of neighbouring
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 §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 Fig. 4.

First, we define how much of the galaxy is covered by a background source. We
call this the overlapping percentage, $\Theta$. For this we measure
segmentation maps both for the central source and all background sources of
each image stamp. $\Theta$ 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,
$\rm BG_{flux}$, by averaging the flux of the sources within the segmentation
map over its area. $\rm BG_{flux}$ values are given in units of $e\ s^{-1}\
pix^{-1}$. This ranges from $\rm BG_{flux}\sim 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 neighbouring sources that outshine the central galaxy. Stars can
also be identified by this method.

We use the overlapping percentage, $\Theta$, and the flux of background
sources, $\rm 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 $\sim 160,000$ simulated and candelized images. The blue box
framed by the dashed line defines a region of galaxies with low contamination,

$\Theta<15\%,\quad\rm BG_{flux}<\times 10^{-1.5}\ e\ s^{-1}\ pix^{-1},$ (6)

which can be considered as a conservative choice. We find that $\sim 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.

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 Appendix A. By inference,
any discussion based on contamination measurements in our simulation sample is
also valid for the CANDELS observations.

##  3 Methods

We use a deep learning framework with a Convolutional Neural Network (CNN)
based on Ferreira et al. (2020) but with significant updates related to the
improved and more robust data pipeline that was discussed in §2.4. In this
section we describe our Deep Learning analysis (§3.1), where we also highlight
the improvements to Ferreira et al. (2020). In Section §3.2, we discuss how to
avoid overfitting due to the augmentation of the TNG sample, which was part of
our sample pipeline. We further wish to compare the resulting classifications
to “traditional” classifications. We thus measure non-parametric morphology
indices, structural parameters and Sérsic profiles for both the TNG sample and
the CANDELS sample with Morfometryka (Ferrari et al., 2015; Albernaz Ferreira
& Ferrari, 2018; Lucatelli & Ferrari, 2018), for which we provide a brief
overview in Section §3.3.

###  3.1 Deep Learning Classifications

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 non-
linear 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 Convolutional Neural Network (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 128x128 pixels, where the number of convolutional blocks,
convolutional layers, fully connected layers, number of filters, and kernel
sizes, are all defined by the following hyper-parameters:

  * •

number_conv_blocks define the number of convolutional blocks, each will probe
features of different scales;

  * •

number_conv_per_block describe how many convolutions each block will have;

  * •

initial_number_filters define the starting number of filters, that are then
doubled after each convolutional block;

  * •

initial_kernel_size is the initial size of the convolutional kernel, that is
then reduced by 2 after each block, down to a minimum of 3;

  * •

n_fc_layers and size_fc_layers define the number of hidden layers and their
respective size, respectively;

  * •

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 multi-class 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 hyper-parameter. 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 (SGD) 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 hyper-parameters 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.

Hyperparameter | Best Model
---|---
batch_size | 128
number_conv_blocks | 3
number_conv_per_block | 2
initial_number_filters | 32
initial_kernel_size | 11
number_fc_layers | 2
size_fc_layers | 128
l2_regularization | 0.1
dropout | 0.5
Table 2: The best hyper-parameters of our architecture found through Bayesian
Optimization (The GPyOpt, 2016). These define the depth, width, and number of
trainable parameters of our architecture. This process is done using our set-
aside validation samples. The same model is used for all the CANDELS datasets.

###  3.2 Augmentations and Overfitting Avoidance

To avoid overfitting pitfalls from using our CANDELS background augmentation
pipeline (§2.4), we train a suite of models, one for each CANDELS field.
Because we have included areas of all 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 datasets (training, validation and test) are
restricted in this way, guaranteeing that any overfitting of the CANDELS
background will have no impact in 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.

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.

Figure 6: Schematics of the training pipeline leveraging multiple fields for
augmentation. Each pair of models, at low redshift (LZ), and high redshift
(HZ), is trained only with data that is augmented with the CANDELS fields that
are not the target for the model. In this example we show a model designed for
predictions on GOODS North (GDN), trained on data augmented with
characteristics of all the remaining four fields (GDS, COS, EGS, UDS). This
model is also tuned and evaluated in validation and test sets that have only
of target CANDELS field augmentations, ensuring that no overfitting of
neighbouring sources is part of the predictive process.

###  3.3 Galaxy Structure and Morphology

Non-parameteric 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, 2018). Morfometryka measures
asymmetry ($A$) concentrations ($C$), the Gini coefficient ($G$), moment of
light of the brightest pixels ($M_{20}$), 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:

$(A>0.35),\quad(A>S).$

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

$A=\frac{\Sigma|I-I_{180}|}{\Sigma|I|}-A_{bg},$

where $I$ is the image, $I_{180}$ 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 overlayed 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,

$(A>0.1).$

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

Figure 7: Distribution of asymmetries $A$ measured with Morfometryka for our
TNG100-1 sample of galaxies. Star forming non-mergers and post-mergers are
shown in blue and red, respectively. The dashed vertical line illustrate the
typical threshold ($A>0.35$) used to classify galaxies as mergers.

##  4 Results

Here we discuss what our trained models reveal, first from the test dataset of
IllustrisTNG selected galaxies (§4.1), and then applied to the CANDELS fields
(§4.2).

###  4.1 Predictions within IllustrisTNG

Figure 8: Performance metrics for our four trained models and comparison with
the classical Asymmetry index $A$ for the simulated images. Left: ROC curves
for both the network trained with the pristine mocks dataset (dashed lines)
and with the realistic mocks dataset (solid lines) applied to both datasets,
color coded in red (pristine) and blue (realistic). The green dotted line
indicates the ROC curve for a classifier using only the asymmetry $A$. The
area under each curve can be read in the label. Right: Precision-Completeness
diagrams for the baseline network trained with the CANDELS matched mocks with
asymmetry $A$, color-coded by classification threshold levels for CNN
(inferno) and asymmetry (viridis). A small region in red is printed over the
asymmetry curve to point out the region where the classification threshold is
$A>0.35$.

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
datasets, one with full HST-matched properties including a CANDELS background
patch of the sky (§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 §3.2.

To compare between models and realizations of these datasets, 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 $\rm P(NMSF)+P(PM)=1$. Fig. 8 displays the overall
performance for each network.

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 datasets. 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
dataset are plotted as solid lines, while models trained with the pristine
mocks dataset are shown by dashed lines. Furthermore, the color conveys the
dataset 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 (AUC) 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 datasets confirms 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 $\sim
10\%$ the clean network when applied to the realistic mocks dataset, 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 inferno colormap, and for the
classic asymmetry parameter in viridis corlormap. 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 datasets, like the case for galaxy mergers.

The precision and completeness of the asymmetry behaves in unpredictable ways.
First, the precision of the selection increases slowly, then it decreases
again around ($A\sim 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 $\sim 80\%$ of the cases. Figure
9 shows the confusion matrix for the realistic mocks dataset identified within
each individual CANDELS field, as well as for the pristine mocks sample, where
accuracy reaches $\sim 90\%$. All classifications are done with the model
trained with the realistic mocks. We show True Positives (TP) and True
Negatives (TN) in blue, and False Positives (FP) and False Negatives in pink.
The CLEAN case represents the best case scenario, where our current method and
dataset achieves an even higher performance of $\sim 91\%$ True Positives. A
histogram of the redshift distribution for each cell helps to visualize any
possible biases in redshift for the miss-classification cases. This
demonstrates that the models are more likely to correctly classify low
redshift galaxies, as they represent the majority of the samples.

Figure 9: Confusion matrix for all the samples matched to CANDELS fields as
well as the the pristine sample (highlighted by gray shading in the bottom
right). These confusion matrices were evaluated with the ensemble of models
trained with the CANDELS matched mocks. We show True Negative (TN) and True
Positives (TP) highlighted in blue while the False Negatives (FN) and False
Positives (FP) are shown in pink. The colors are based on the rate percentage,
which is also printed in each cell. All the CANDELS fields have TP and FN
rates of around $\sim 80\%$. For the pristine case performance can reach as
high as $\sim 90\%$, marking the intrinsic limit of our method based on the
data available. The histograms show the redshift distribution for the galaxies
in each category, which demonstrate that it is easier to recover correct
classifications at lower redshifts.

####  4.1.1 Impact of Redshift

With the goal to apply 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 errorbars
– sampled from boostraping our testing samples – follow accordingly.

Figure 10: Impact of contamination and redshifts on the performance of our
models. The left and central panels show how accuracy, precision and recall
(blue squares, orange circles and green hexagons, respectively) behave for
increasing percentages of overlap ($\Theta$) and for increased background flux
($BG_{flux}$). In the right panel we show how the accuracy, precision and
recall of our methods change in bins of $\Delta z=0.25$ redshift. Errorbars
are sampled from bootstraping the test sample. The performance gradually
decreases with $z$, decreasing below $80\%$ beyond $z=2$. There is a slight
uptick at $z=2.5$, but with large errorbars. The cutoff at $z>2.5$ is the
result of a combination of small sample size and redshift effects. The black
dashed line at $score=0.8$ indicates the overall accuracy of the model in the
complete test set.

####  4.1.2 Contamination impact on classification

We use the contamination estimates measured in §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, $\Theta$,
and the average background flux per pixel, $BG_{flux}$, as shown in Fig.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 sub-samples 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 falls below the dashed line, which is our contamination
cutoff, i.e.,

$\Theta\sim 15\%,\quad BG_{flux}\sim 10^{-1.5}\ e\ s^{-1}\ \rm pix^{-1}.$

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 Appendix §A.

###  4.2 Classifications on CANDELS

We use our network to carry out predictions in all real CANDELS galaxies at
$0.5<z<3$, $M_{*}\geq 10^{9.5}M_{\odot}$, $S/N>50$ and $H_{\rm 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 non-interacting star forming galaxies using a threshold probability of
$60\%$. Galaxies with probabilities $50\%<P(PM)\wedge P(SF)<60\%$ are not
considered in any class. These represent 2125 galaxies ($\approx 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 in the left panel and star forming galaxies in the right panel.

Figure 11: Examples of CANDELS galaxies with $A>0.1$ classified by our models
into Post-mergers (left) and star forming galaxies (right), with their
redshifts, SFRs, and stellar masses. Images are ranked from left to right with
increasing redshift and top to bottom with increasing SFR. All stamps use a
square-root normalization.

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
$\Delta 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<\log(M_{*}/M_{\odot})<10.0$ (left panel) and high mass systems with
$\log(M_{*}/M_{\odot})>10.0$ (right panel). Non-interacting star forming
galaxies are shown in blue circles while post-mergers are shown in 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.

Figure 12: Relative class fractions for post-mergers and star forming galaxies
vs. redshift for real galaxies in the CANDELS fields. The fraction of post-
mergers increases from 30% at $z\sim 0.75$ to 50% by $z\sim 2$. Errorbars are
drawn from bootstrapping the samples and applying the underlying uncertainty
associated with the performance of our models, which decreases with redshift.

For the lower mass post-mergers we see an upward trend from $\sim 15\%$ at
$z=0.5$ to $\sim 35\%$ at $z=2$, then a slight decrease beyond $z=2$. This is
still consistent with a $\sim 35\%$ fraction within the error-bars. The star
forming galaxies behave in the opposite way, decreasing from $\sim 70\%$ at
$z=0.5$ to around $\sim 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 non-interacting 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 more steep upward
trend from $\sim 20\%$ at $z=0.5$ fraction to $\sim 50\%$ at $z=2$, the
relative fraction of star forming galaxies show a constant value of $\sim
50\%$ at $0.5<z<3$, while the fraction of the rest of the sample (dashed line)
goes from $\sim 75\%$ to $\sim 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, 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 $\sim 20\%$ from $z=0.5$ ($\sim 30\%$) to $z=2$
($\sim 50\%$), which agree with the results for the complete sample.
Furthermore, our methods classify $\sim 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 an UMAP (Fig. 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 further away from the
center, which represents how different these objects are for the network.
Images of galaxies the network struggles to identify are mixed in the bottom
middle, representing the region where both probabilities are similar
$P(PM)\sim P(SF)$.

Figure 13: UMAP representation of the output from the last dense layer of the
network. This representation shows the parameter space used for the network to
generate the final probability. Probabilities are highest in the extremes at
the top, and uncertainty increases due to increased contamination as we go
along this structure towards the middle. Same random examples of CANDELS
galaxies are placed close to their points in this manifold. Small regions
identified by circles show the clustering of non-galactic detections in this
parameter space, located close to the region of uncertain classifications at
the bottom of the UMAP. ”Stars” are stars in the center and ”Stars in FOV”
correspond to stars at the edge of the stamps.

##  5 Implications

Figure 14: Mean distance to the star formation main sequence ($\Delta_{MS}$)
vs. log stellar mass in bins of redshift for CANDELS galaxies above the star
formation main sequence. Post-mergers are plotted as red squares, star forming
galaxies as blue circles. Errorbars are estimated using bootstrapping and show
$\pm 1\ \sigma$. The classes are indistinguishable between $0.5<z<2$, both
increase similarly in $\Delta_{MS}$ as we increase redshift. This represents
the increase in scatter above the main sequence. All galaxies included in this
diagram lie above the SFMS, as we are only interested in exploring the
scattering above the SFMS. For the last redshift bin ($2<z<2.5$), there is a
significant difference between the two classes both at the low mass end and at
the high mass end. At low masses ($\log(M_{*}/M_{\odot})<10.0$), post-mergers
scatter higher than star forming galaxies with a difference of
$\Delta_{MS}\sim 0.1\ \rm dex$. At high masses ( $\log(M_{*}/M_{\odot})>10.0$)
the trend reverses and star forming galaxies scatter higher with a
$\Delta_{MS}$ difference of $\sim 0.1\ \rm dex$. However, care needs to be
taken in the interpretation of this trend as it could be spurious or
insignificant given the errorbars and the performance metrics of our models at
high redshift (Fig. 10).

Making use of the classifications from our deep learning models, we first
explore the impact of major-mergers on classifications above the Star Forming
Main Sequence as parametrized by Schreiber et al. (2015) (§5.1). We then
discuss the structure of the two galaxy classes using Sérsic profile
measurements (§5.2). In §5.3, we update classifications from Ferreira et al.
(2020) with our new specialized model, thus increasing certainty for
previously undefined classifications. We then add to the discussion proposed
by Bickley et al. (2021) regarding the Bayesian limitations of classifying
post-mergers by considering an evolving merger fraction. We finish with §5.5,
in which we compare extracted features from real CANDELS galaxies to features
extracted from IllustrisTNG galaxies, as a way to address the challenges of
transferring the model from simulations to real observations.

###  5.1 Classifications above the star forming main sequence

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 star formation main sequence 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 star forming main sequence as parametrized 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 star formation main sequence ($\Delta_{MS}$), as:

$\Delta_{MS}=\log(\rm SFR)-\log(\rm SFR_{MS}),$ (7)

where $\log(\rm SFR)$ is the log star formation rate of a particular galaxy
and $\log(\rm SFR_{MS})$ is the parametrization 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 $\Delta_{MS}$, with all offsets well within
the errorbars. However, for redshifts $2.5<z<3.0$ post-mergers with
$\log(M_{*}/M_{\odot})<10.0$ are on average $\sim 0.1\ \rm dex$ higher than
star forming galaxies of the same mass. The opposite is found for
$\log(M_{*}/M_{\odot})>10.0$, however uncertainty is higher here.
Additionally, $\Delta_{MS}$ increases with redshift in all cases, which
describes a larger scatter above the star forming main sequence. However,
given the performance metrics of our models at high redshift (Figure 10), we
cannot 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.

In summary, locations of post-mergers and non-interacting galaxies in the star
forming main sequence 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 specific star formation rates by a factor of $\sim 2$, but that this
effect decays in timescales of $\sim 0.5\ \rm 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.

###  5.2 Structure and light profiles

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 $\sim 1\ \rm 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 non-interacting star
forming galaxies, and 2) these types of galaxies are intrinsically different
from each other.

Figure 15: Sersic index distribution for post-mergers and star forming
galaxies, in red and blue respectively. Post-mergers display more concentrated
light distributions with $n\sim 1.8_{-0.6}^{+0.7}$ while the star forming
galaxies have $n\sim 1.1_{-0.5}^{+0.5}$ consistent with disk dominated
galaxies.

###  5.3 Merger Fractions and Rates

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

Our new dataset 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 non-
interacting star forming galaxies or if any non-mergers can be re-classified
as post-mergers. This is done by comparing the probabilities for major-mergers
and non-mergers, $P(MM)$ and $P(NM)$, respectively, from Ferreira et al.
(2020) to the new probabilities $P(PM)$ and $P(SF)$. We update a non-merger
classification to post-merger if

$P(PM)>P(NM),$

and update the major-merger classifications to non-merger if

$P(SF)>P(MM).$

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 $\sim 5\%$ of major-mergers reclassified as star forming non mergers,
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.

The updated fit of the cosmic evolution of the merger fraction, $f_{m}(z)$

$f_{m}(z)=0.011\pm 0.002\times(1+z)^{2.71\pm 0.31},$ (8)

with errors estimated with boostraping, 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 ($\tau_{m}=0.5\ \rm Gyr$) used
in our selection (§2.1) with this merger fraction through

$\mathcal{R}=\frac{f_{m}}{\tau_{m}}.$ (9)

The updated galaxy major merger rate is

$\mathcal{R}=0.022\pm 0.006\times(1+z)^{2.71\pm 0.31}.$ (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.

Figure 16: Major-merger fractions as a function of redshift. We show corrected
merger fractions from Ferreira et al. (2020) by re-classifying galaxies with
our new method in mergers and non-mergers, shown in green. The original
estimates are shown in gray.

###  5.4 Bayesian Analysis of Mergers

Figure 17: Extracted features by our networks in a UMAP 2D representation. For
each model in our ensamble we generate a UMAP from the extracted features from
the last convolutional layer of the trained networks, both when applied to the
Illustris galaxies, color-coded by the class, and to the unlabeled CANDELS
galaxies shown in black dots. Both Illustris and CANDELS extracted features
populate the same region of this representation, showing that the features
used by the network to then perform the classification task are in general
domain invariant. Additionally, both classes – post-mergers and star forming
galaxies – form separated clusters with some overlapping. Classification could
be done in this representation alone, but it is then better organized by the
fully connected layers that combine these features to produce the final output
probability.

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.

$P(\rm\mathbf{M|S})=P(\mathbf{S|M})\times\frac{P\left(\rm\mathbf{M}\right)}{P\left(\rm\mathbf{S}\right)},$
(11)

where $P(\rm\mathbf{M|S})$ is the probability of a merger, given that a method
used to select mergers, ($\rm\mathbf{S}$), identifies it as such. The value of
P($\rm\mathbf{M}$) is the probability that an object is a merger before a
selection of merger is made. P($\rm\mathbf{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($\rm\mathbf{S}$)
as:

$P\left(\rm\mathbf{S}\right)=P(\rm\mathbf{S|M})\times
P(\rm\mathbf{M})+P(\rm\mathbf{S|NM})\times P(\rm\mathbf{NM})$ (12)

where $\rm\mathbf{NM}$ standard for non-mergers, where P$(\rm\mathbf{S|NM})$
is the probability of identifying correctly a non-mergers, and the value of
P($\rm\mathbf{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_{\rm m}$, and the
machine learning probability of finding a merger/non-merger as $p_{\rm m}$ and
$p_{\rm mn}$ then we can rewrite equation (11), as:

$P(\rm\mathbf{M|S})=\frac{p_{\rm m}\times f_{\rm m}}{(p_{\rm m}f_{\rm
m}+(1-p_{\rm nm})(1-f_{m}))}.$ (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 non-merger is 0.9, and the merger
fraction $f_{\rm m}$ = 0.1, then the probability that a galaxy identified as a
merger is actually a merger is $P(\rm\mathbf{M|S})=0.5$. This implies that
even when the accuracy of finding mergers and non-mergers is $90\%$, at the
lowest redshifts, where the merger fraction is low $\sim 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 $\sim 80\%$ when the
merger fraction is as high as $f_{\rm m}\sim 0.3$.

We can generalise the equation for $P(\rm\mathbf{M|S})$, as a function of $z$,
by considering how the merger fraction $f_{\rm m}$ evolves with redshift, such
that:

$P(\rm\mathbf{M|S},z)=\frac{p_{\rm m}(z)\times f_{m}(z)}{p_{\rm
m}(z)f_{m}(z)+(1-p_{\rm nm}(z))(1-f_{m}(z))},$ $f_{m}(z)\equiv
f_{0}(1+z)^{m},$

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
non-mergers. From our results here, our method effectiveness for correctly
classifying mergers increases from $\sim 40\%$ at $z\sim 0.5$ to $\sim 70\%$
at $z\sim 3$. This is likely what can account for some of our misidentified
galaxies as discussed in §4 when discussing the success of our method of
separating star forming systems from those that are undergoing mergers.

These are conservative estimates which do not include the fact that we pre-
select CANDELS galaxies based on their asymmetry. This should increase
$P(\rm\mathbf{M|S},z)$ further since $f_{\rm m}$ is higher among galaxies with
$A>0.1$.

###  5.5 On Domain Adaptation Issues

There is a growing concern on the applicability of simulation trained deep
learning models when applied to a 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) show that adopting techniques for domain adaptation could increase
model performance when applying to the target domain by $20\%$.

In our case, the source domain are the IllustrisTNG galaxies and the target
domain the CANDELS observations. To check if we need to apply domain
adaptation techniques to this particular problem, we used UMAPs (Uniform
Manifold Approximation and Projection, as described in Sec. 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 visualize222We also
tested with t-SNEs with similar results.. 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 fields 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.333We also
tested generating random noise images to check their position in this
parameter space. As expected, they cluster away from the image regions,
forming its own outlier region which is far from the main locus where galaxies
are found.

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 ($\sim
140.000$ images) but also help with domain confusion within the network. Thus,
we do not include any domain adaptation process in our pipeline.

##  6 Summary

To shed light on the nature of peculiar/irregular objects at intermediate to
high redshifts, we have constructed a framework based on forward-modelling of
cosmological simulations with deep learning algorithms, that allows
classifications with physically motivated labels based on the formation
history of galaxies.

We used data from the IllustrisTNG TNG100-1 simulation to create realistic
mocks of galaxies with CANDELS-like properties, including a full radiative
transfer treatment with SKIRT for two specific classes of galaxies: post-
mergers and non-merging star forming galaxies. These are selected so that
their main difference is their formation history.

We produced a dataset of $\sim 160,000$ images of simulated IllustrisTNG
galaxies with realistic visual properties that mimick CANDELS observations in
the redshift range $0.5<z<3.0$. The images are used to train Deep
Convolutional Neural Networks to distinguish between formation histories of
post-mergers and star forming galaxies. The main conclusions drawn from this
work are summarized as follow:

  * •

The classifier network combined with our new dataset produces classification
models with a balanced performance of $\sim 80\%$ accuracy, precision, and
completeness when applied to a single-band imaging dataset, outperforming the
asymmetry ($A$) by at least 25% within the simulated data. Additionally, for
pristine images without any contamination and observational effects, the
theoretical limit of our model is $\sim 91\%$ accuracy. This is evidence that
using the asymmetry ($A$) alone for ambiguous morphological cases might
generate highly contaminated samples.

  * •

We define two new contamination indicators, the overlapping percentage,
$\Theta$, and the average flux of the background sources, $BG_{flux}$, by
leveraging how simulated galaxies are combined with true CANDELS background
sky patches. $\Theta$ controls how sources overlap and are projected in the
same stamp, while the $BG_{flux}$ value probes the effect of the brightness of
external sources on the classification of the central object. These allow us
to explore in detail how deep learning classifications are impacted by
contamination. We show that both crowded environments and projections and the
relative brightness of external sources to the central galaxy negatively
impacted our results. Based on this, we define quality control limits to our
approach within the CANDELS fields as $\Theta\sim 10\%$ and
$BG_{flux}<10^{-3}$ e s-1 pix-1. Although not universal, these limits provide
guidelines for sample selection when applying our models to data.

  * •

By applying our model to real CANDELS observations of galaxies with high
asymmetries, we show that the relative fraction of post-mergers to star
forming galaxies increases with higher redshift for two mass regimes. For low
mass sources ($9.5<\log(M_{*}/M_{\odot})<10.0$), the post-merger fraction
increases by $\sim 20\%$ within $0.5<z<2.0$, while the fraction of star
forming galaxies decreases by $\sim 15\%$ in the same redshift range. In the
high mass case ($\log(M_{*}/M_{\odot})>10.0$), the post-merger fraction
increases by $\sim 25\%$ at $0.5<z<2.0$, while the fraction of star forming
galaxies stays broadly constant.

  * •

We explore the impact of major mergers on galaxies located above the star
formation main sequence (SFMS) as parametrized by Schreiber et al. (2015). We
separate CANDELS galaxies above the SFMS in the classes provided by our model
and in bins of stellar mass. At $0.5<z<2.0$ we do not find any clear signs
that major mergers play a critical role on the scattering above SFMS, with
similar trends for post-mergers and star forming galaxies. However, in the
highest redshift bin with good sample statistics ($2.0<z<2.5$) we see a post-
merger driven SFR enhancement at lower masses of about $\sim 0.1\rm\ dex$.

  * •

We show that the light distribution parametrized through Sérsic profiles of
the CANDELS galaxies classified by our models as post-mergers are
intrinsically distinct from those classified as star forming galaxies. The
star forming galaxies sample is dominated by disk-like objects with an average
Sérsic index of $n=1.1_{-0.5}^{+0.5}$ while the post-mergers have more
concentrated light profiles corresponding to higher central concentration with
$n=1.8_{-0.6}^{+0.7}$, with a long tail at higher Sérsic indices. Moreover,
when we increase the probability threshold of our classifications to improve
the purity of our selections, only the post-merger distribution display higher
Sérsic indices. Evidently, our model predicts that post-mergers are more
likely to be bulge-dominated galaxies.

  * •

By using our updated data pipeline and models specifically tailored to
distinguish between post-mergers and star forming galaxies, we revisit the
merger fractions and merger rates from Ferreira et al. (2020) by correcting
ambiguous cases. This leads to updated galaxy merger rates that are slightly
lower, but consistent with previously reported rates: $\mathcal{R}=0.022\pm
0.006\times(1+z)^{2.71\pm 0.31}$.

  * •

We show that our models use similar features to classify IllustrisTNG and real
CANDELS galaxies, with no clear discrepancy between the two domains. Using the
features extracted by the convolutional layers of our network, we generate
UMAPs, which visualize the complex parameter space in 2D. Features of
IllustrisTNG galaxies and CANDELS galaxies overlap for all the CANDELS fields.
Although the CANDELS galaxies do not span the entire feature space of the
IllustrisTNG galaxies used here, they are contained within that feature space.

Our machine learning driven approach provides a new way to investigate the
formation history of galaxies with models that are informed by cosmological
simulations. This includes the use of the models themselves, and the
application of these models within accurate observing conditions.

Nevertheless, currently we are still limited to high-mass major merger cases
due to resolution limitations from the simulations and mass completeness from
the observations. In the upcoming years, combining the next generation of high
resolution, small box simulations (e.g. TNG50-1, New Horizons) with
observational data from the James Webb Space Telescope (JWST) and Euclid
Telescope will open a new window to incorporate the effect of minor mergers
and lower mass systems. Together, this will represent a major step towards
uncovering unresolved questions of galaxy evolution.

## acknowledgments

The authors thank the anonymous referee for the detailed review that improved
the paper greatly. The authors thank Centre for Astronomy and Particle Theory
of University of Nottingham for providing all computational infrastructure
necessary to run the training steps to produce the model described here. This
study was financed in part by the Coordernação de Aperfeiçoamento de Pessoal
de Nível Superior - Brazil (CAPES). CJC acknowledges support from the European
Research Council (ERC) Advanced Investigator Grant EPOCHS (788113). UK
acknowledges support from the Science and Technology Facilities Council
through grant number RA27PN. The IllustrisTNG simulations were undertaken with
compute time awarded by the Gauss Centre for Supercomputing (GCS) under GCS
Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the
supercomputer Hazel Hen at the High Performance Computing Center Stuttgart
(HLRS), as well as on the machines of the Max Planck Computing and Data
Facility (MPCDF) in Garching, Germany.

## Data Availability

Ready to use models are publicly available for anyone to
download444https://github.com/astroferreira/FERREIRA2022, with accompanying
code.

The post-processed IllustrisTNG data used for training is stored in large
TFRecords binary files and are available upon request.

## References

  * Ackermann et al. (2018) Ackermann, S., Schawinski, K., Zhang, C., Weigel, A. K., & Dennis Turp, M. 2018, Monthly Notices of the Royal Astronomical Society, 479, 415
  * Albernaz Ferreira & Ferrari (2018) Albernaz Ferreira, L., & Ferrari, F. 2018, Monthly Notices of the Royal Astronomical Society, 473, 2701
  * Allen et al. (2008) Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20
  * Bertone & Conselice (2009) Bertone, S., & Conselice, C. J. 2009, MNRAS, 396, 2345
  * Bickley et al. (2021) Bickley, R. W., Bottrell, C., Hani, M. H., et al. 2021, MNRAS, 504, 372
  * Blumenthal et al. (2020) Blumenthal, K. A., Moreno, J., Barnes, J. E., et al. 2020, MNRAS, 492, 2075
  * Bottrell et al. (2021) Bottrell, C., Hani, M. H., Teimoorinia, H., Patton, D. R., & Ellison, S. L. 2021, MNRAS, arXiv:2201.03579
  * Bottrell et al. (2022) —. 2022, MNRAS, 511, 100
  * Bottrell et al. (2019) Bottrell, C., Hani, M. H., Teimoorinia, H., et al. 2019, MNRAS, 490, 5390
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000
  * Camps & Baes (2015) Camps, P., & Baes, M. 2015, Astronomy and Computing, 9, 20
  * Camps & Baes (2020) —. 2020, Astronomy and Computing, 31, doi:10.1016/j.ascom.2020.100381
  * Camps et al. (2016) Camps, P., Trayford, J. W., Baes, M., et al. 2016, MNRAS, 462, 1057
  * Cheng et al. (2019) Cheng, T.-y., Conselice, C. J., Ara, A., & Li, N. 2019
  * Ćiprijanović et al. (2020) Ćiprijanović, A., Snyder, G. F., Nord, B., & Peek, J. E. G. 2020, Astronomy and Computing, 32, 100390
  * Ćiprijanović et al. (2021) Ćiprijanović, A., Kafkes, D., Downey, K., et al. 2021, MNRAS, 506, 677
  * Conselice (2003) Conselice, C. J. 2003, ApJS, 147, 1
  * Conselice (2006) —. 2006, Monthly Notices of the Royal Astronomical Society, 373, 1389
  * Conselice (2014) —. 2014, Annual Review of Astronomy and Astrophysics, 52, 291
  * Conselice et al. (2003) Conselice, C. J., Bershady, M. A., Dickinson, M., & Papovich, C. 2003, The Astronomical Journal, 126, 1183
  * Conselice et al. (2000) Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, ApJ, 529, 886
  * Conselice et al. (2005) Conselice, C. J., Blackburne, J. A., & Papovich, C. 2005, ApJ, 620, 564
  * Conselice et al. (2008) Conselice, C. J., Rajgor, S., & Myers, R. 2008, MNRAS, 386, 909
  * Conselice et al. (2007) Conselice, C. J., Bundy, K., Trujillo, I., et al. 2007, MNRAS, 381, 962
  * Dickinson et al. (2000) Dickinson, M., Hanley, C., Elston, R., et al. 2000, ApJ, 531, 624
  * Duncan et al. (2014) Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960
  * Duncan et al. (2019) Duncan, K., Conselice, C. J., Mundy, C., et al. 2019, ApJ, 876, 110
  * Duncan et al. (2018a) Duncan, K. J., Jarvis, M. J., Brown, M. J. I., & Röttgering, H. J. A. 2018a, MNRAS, 477, 5177
  * Duncan et al. (2018b) Duncan, K. J., Brown, M. J. I., Williams, W. L., et al. 2018b, MNRAS, 473, 2655
  * Ferrari et al. (2015) Ferrari, F., de Carvalho, R. R., & Trevisan, M. 2015, ApJ, 814, 55
  * Ferreira et al. (2020) Ferreira, L., Conselice, C. J., Duncan, K., et al. 2020, ApJ, 895, 115
  * Gao et al. (2018) Gao, Y., Bao, M., Yuan, Q., et al. 2018, ApJ, 869, 15
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, Astrophysical Journal, Supplement Series, 197, doi:10.1088/0067-0049/197/2/35
  * Hani et al. (2020) Hani, M. H., Gosain, H., Ellison, S. L., Patton, D. R., & Torrey, P. 2020, MNRAS, 493, 3716
  * Hogg (1999) Hogg, D. W. 1999, arXiv.org, astro-ph, 5116.  http://arxiv.org/abs/astro-ph/9905116
  * Huertas-Company et al. (2015) Huertas-Company, M., Gravet, R., Cabrera-Vives, G., et al. 2015, ApJS, 221, 8
  * Huertas-Company et al. (2016) Huertas-Company, M., Bernardi, M., Pérez-González, P. G., et al. 2016, Monthly Notices of the Royal Astronomical Society, 462, 4495
  * Huertas-Company et al. (2018) Huertas-Company, M., Primack, J. R., Dekel, A., et al. 2018, ApJ, 858, 114
  * Huertas-Company et al. (2019a) Huertas-Company, M., Rodriguez-Gomez, V., Nelson, D., et al. 2019a, MNRAS, 489, 1859
  * Huertas-Company et al. (2019b) —. 2019b, MNRAS, 489, 1859
  * Jogee et al. (2009) Jogee, S., Miller, S. H., Penner, K., et al. 2009, ApJ, 697, 1971
  * Kartaltepe et al. (2015) Kartaltepe, J. S., Mozena, M., Kocevski, D., et al. 2015, ApJS, 221, 11
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, Astrophysical Journal, Supplement Series, 197, doi:10.1088/0067-0049/197/2/36
  * Loshchilov & Hutter (2016) Loshchilov, I., & Hutter, F. 2016, arXiv e-prints, arXiv:1608.03983
  * Lotz et al. (2008) Lotz, J. M., Jonsson, P., Cox, T. J., & Primack, J. R. 2008, Monthly Notices of the Royal Astronomical Society, 391, 1137
  * Lotz et al. (2004) Lotz, J. M., Primack, J., & Madau, P. 2004, The Astronomical Journal, 128, 163
  * Lu et al. (2017) Lu, Z., Pu, H., Wang, F., Hu, Z., & Wang, L. 2017, arXiv e-prints, arXiv:1709.02540
  * Lucatelli & Ferrari (2018) Lucatelli, G., & Ferrari, F. 2018, 000
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, Annual Review of Astronomy and Astrophysics, 52, 415
  * Mager et al. (2018) Mager, V. A., Conselice, C. J., Seibert, M., et al. 2018, ApJ, 864, 123
  * Man et al. (2016) Man, A. W. S., Zirm, A. W., & Toft, S. 2016, ApJ, 830, 89
  * Mantha et al. (2018) Mantha, K. B., McIntosh, D. H., Brennan, R., et al. 2018, MNRAS, 475, 1549
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113
  * Martin et al. (2020) Martin, G., Kaviraj, S., Hocking, A., Read, S. C., & Geach, J. E. 2020, MNRAS, 491, 1408
  * McInnes et al. (2018) McInnes, L., Healy, J., & Melville, J. 2018, arXiv e-prints, arXiv:1802.03426
  * Mortlock et al. (2013) Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2013, MNRAS, 433, 1185
  * Mortlock et al. (2013) Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2013, MNRAS, 433, 1185
  * Mundy et al. (2017) Mundy, C. J., Conselice, C. J., Duncan, K. J., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 3507
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, Monthly Notices of the Royal Astronomical Society, 475, 624
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2
  * Papovich et al. (2005) Papovich, C., Dickinson, M., Giavalisco, M., Conselice, C. J., & Ferguson, H. C. 2005, ApJ, 631, 101
  * Pearson et al. (2019a) Pearson, W. J., Wang, L., Trayford, J. W., Petrillo, C. E., & van der Tak, F. F. S. 2019a, Astronomy & Astrophysics, 626, A49
  * Pearson et al. (2019b) Pearson, W. J., Wang, L., Alpaslan, M., et al. 2019b, 1
  * Pillepich et al. (2018a) Pillepich, A., Springel, V., Nelson, D., et al. 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018b, MNRAS, 475, 648
  * Pillepich et al. (2019) Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196
  * Powers (2011) Powers, D. M. W. 2011, Journal of Machine Learning Technology, 2, 37
  * Reiman & Göhre (2019) Reiman, D. M., & Göhre, B. E. 2019, Monthly Notices of the Royal Astronomical Society, 485, 2617
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, Monthly Notices of the Royal Astronomical Society, 449, 49
  * Rodriguez-Gomez et al. (2018) Rodriguez-Gomez, V., Snyder, G. F., Lotz, J. M., et al. 2018, 000.  http://arxiv.org/abs/1809.08239
  * Sazonova et al. (2021) Sazonova, E., Alatalo, K., Rowlands, K., et al. 2021, ApJ, 919, 134
  * Schreiber et al. (2015) Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74
  * Simons et al. (2019) Simons, R. C., Kassin, S. A., Snyder, G. F., et al. 2019, The Astrophysical Journal, 874, 59
  * Snyder et al. (2017) Snyder, G. F., Lotz, J. M., Rodriguez-Gomez, V., et al. 2017, Monthly Notices of the Royal Astronomical Society, 468, 207
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676
  * Tan & Le (2019) Tan, M., & Le, Q. V. 2019, arXiv e-prints, arXiv:1905.11946
  * The GPyOpt (2016) The GPyOpt, A. 2016, GPyOpt: A Bayesian Optimization framework in python, \url{http://github.com/SheffieldML/GPyOpt}, ,
  * Tohill et al. (2021) Tohill, C., Ferreira, L., Conselice, C. J., Bamford, S. P., & Ferrari, F. 2021, ApJ, 916, 4
  * Trayford et al. (2017) Trayford, J. W., Camps, P., Theuns, T., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 771
  * Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Monthly Notices of the Royal Astronomical Society, 444, 1518
  * Vogelsberger et al. (2020) Vogelsberger, M., Nelson, D., Pillepich, A., et al. 2020, Monthly Notices of the Royal Astronomical Society, 492, 5167
  * Walmsley et al. (2020) Walmsley, M., Smith, L., Lintott, C., et al. 2020, MNRAS, 491, 1554
  * Walmsley et al. (2021) Walmsley, M., Lintott, C., Géron, T., et al. 2021, MNRAS, arXiv:2102.08414
  * Wang et al. (2020) Wang, L., Pearson, W. J., & Rodriguez-Gomez, V. 2020, A&A, 644, A87
  * Whitney et al. (2021) Whitney, A., Ferreira, L., Conselice, C. J., & Duncan, K. 2021, ApJ, 919, 139
  * Williams et al. (1996) Williams, R. E., Blacker, B., Dickinson, M., et al. 1996, \aj, 112, 1335
  * Windhorst et al. (2002) Windhorst, R. A., Taylor, V. A., Jansen, R. A., et al. 2002, ApJS, 143, 113
  * Zanisi et al. (2021) Zanisi, L., Huertas-Company, M., Lanusse, F., et al. 2021, MNRAS, 501, 4359
  * Zhou et al. (2020) Zhou, P., Feng, J., Ma, C., et al. 2020, arXiv e-prints, arXiv:2010.05627

##  Appendix A Contamination Network

Based on the contamination measurements described in §2.5, we devise a new
neural network with the goal to predict the overlapping percentage ($\Theta$)
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 the contamination information from our data pipeline (§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 replacing the final sigmoid layer with a linear
activation function, changing 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 $\Theta$ 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 root mean square errors in the order of $\sim 10^{-3}$ for
$BG_{flux}$ and $\sim 5\%$ for $\Theta$. 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 $\Theta$ and
$BG_{\rm flux}$.

Figure 18: Performance of the contamination quantification network. (top) The
relationship between true and predicted values for $BG_{\rm flux}$ and
(bottom) relationship between truth and predicted values for $\Theta$. Pearson
and Spearman correlation indices are displayed for each case, as well as the
root mean square error.

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

Figure 19: Panel with cutouts of IllustrisTNG galaxies demonstrating four
different selections on the $O_{f}$ and $<f_{BG}>$ parameter space. Isolated
galaxies with almost no noticeable contamination have low overlap and low
background flux (top left). Low overlap and high background flux show cases
where the central galaxy is overshadowed by a bright companion, but with no
overlapping (top right). High overlap and low background flux show galaxies
overlapping with similar brightness, cases where projection effects can be
misinterpreted as a major merger (bottom left). High overlap and high
background flux show central galaxies with very large and bright companions
that extend over its segmentation map (bottom right). This illustrates how
useful these two measurements can be for proper selections.