# Constraining the major merger history of z 3 - 9 galaxies using JADES: dominant in-situ star formation

Dávid Puskás, Sandro Tacchella, Charlotte Simmonds, Kevin Hainline, Francesco d'Eugenio, Stacey Alberts, Santiago Arribas, William M. Baker, Andrew J. Bunker, Stefano Carniani, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Dávid Puskás, Sandro Tacchella, Charlotte Simmonds, Kevin Hainline, Francesco d'Eugenio, et al. Constraining the major merger history of z 3 - 9 galaxies using JADES: dominant in-situ star formation. Monthly Notices of the Royal Astronomical Society, 2025, 9, p. 729-740. 10.1093/mn-ras/staf813. insu-05101653

HAL Id: insu-05101653

https://insu.hal.science/insu-05101653v1

Submitted on 4 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining the major merger history of $z \sim 3–9$ galaxies using JADES: dominant *in situ* star formation

Dávid Puskás<sup>1,2*</sup> Sandro Tacchella<sup>1,2</sup> Charlotte Simmonds<sup>1,2</sup> Kevin Hainline<sup>1,3</sup>

Francesco D'Eugenio<sup>1,2</sup> Stacey Alberts<sup>1,3</sup> Santiago Arribas<sup>1,4</sup> William M. Baker<sup>1,5</sup>

Andrew J. Bunker<sup>1,6</sup> Stefano Carniani<sup>1,7</sup> Stéphane Charlot<sup>1,8</sup> Qiao Duan<sup>1,2</sup> Daniel J. Eisenstein<sup>1,9</sup>

Zhiyuan Ji<sup>1,3</sup> Benjamin D. Johnson<sup>1,9</sup> Gareth C. Jones<sup>1,2</sup> Roberto Maiolino<sup>1,2,10</sup>

William McClymont<sup>1,12</sup> Marcia Rieke<sup>1,3</sup> Pierluigi Rinaldi<sup>1,3</sup> Brant Robertson<sup>11</sup> Hannah Übler<sup>1,12</sup>

Christina C. Williams<sup>1,13</sup> Christopher N. A. Willmer<sup>1,3</sup> Chris Willott<sup>1,14</sup> and Joris Witstok<sup>1,15,16</sup>

Accepted 2025 May 12. Received 2025 April 15; in original form 2025 February 3

## ABSTRACT

We present a comprehensive analysis of galaxy close-pair fractions and major merger rates to evaluate the importance of mergers in the hierarchical growth of galaxies over cosmic time. This study focuses on the previously poorly understood redshift range of $z \approx 3–9$ using JADES observations. Our mass-complete sample includes primary galaxies with stellar masses of $\log(M_*/M_\odot) = [8, 10]$, having major companions (mass ratio $\ge 1/4$) selected by 5–30 pkpc projected separation and redshift proximity criteria. Pair fractions are measured using a statistically robust method incorporating photometric redshift posteriors and available spectroscopic data. The pair fraction evolves with redshift and shows dependence on the stellar mass: at $\log(M_*/M_\odot) = [8.0, 8.5]$ there is an increase up to $z \sim 5–6$, followed by a turnover, while at higher stellar masses there is a flattening and weak decline with increasing redshift. Similarly, the derived galaxy major merger rate increases and flattens beyond $z \sim 6$ to $2–8 \text{ Gyr}^{-1}$ per galaxy, showing a weak scaling with stellar mass, driven by the evolution of the galaxy stellar mass function. A comparison between the cumulative mass accretion from major mergers and the mass assembled through star formation indicates that major mergers contribute approximately 3–13 per cent to the total mass growth over the studied redshift range, which is in agreement with the *ex situ* mass fraction estimated from our simple numerical model. These results highlight that major mergers contribute little to the direct stellar mass growth compared to *in situ* star formation but could still play an indirect role by driving star formation itself.

**Key words:** galaxies: formation – galaxies: high-redshift – galaxies: interactions.

## 1 INTRODUCTION

Galaxy mergers have long been predicted to play a significant role in the evolution and overall mass build-up of the galaxy population

through the hierarchical growth of dark matter haloes in the Lambda cold dark matter (ΛCDM) cosmological model (White & Rees 1978). A fundamental question from an evolutionary perspective is how massive galaxies have accumulated their stellar mass. It is believed that galaxies grow through two main channels: the formation of new stars via the accretion of cold gas from their surrounding intergalactic medium and mergers with nearby galaxies, resulting in a single,

<footnote n="1">¹Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="2">²Kavli Institute for Cosmology, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="3">³Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA</footnote>

<footnote n="4">⁴Centro de Astrobiología (CAB), CSIC-INTA, Cra. de Ajalvir Km. 4, E-28850, Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="5">⁵DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark</footnote>

<footnote n="6">⁶Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="7">⁷Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="8">⁸Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="9">⁹Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA</footnote>

<footnote n="10">¹⁰Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="11">¹¹Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 96054, USA</footnote>

<footnote n="12">¹²Max-Planck-Institut für extraterrestrische Physik (MPE), Gießenbachstraße 1, D-85748 Garching, Germany</footnote>

<footnote n="13">¹³NSF National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="14">¹⁴NRC Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="15">¹⁵Cosmic Dawn Center (DAWN), Jagtvej 128, DK-2200, Copenhagen, Denmark</footnote>

<footnote n="16">¹⁶Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark</footnote>

<footnote>*E-mail: dp670@cam.ac.uk</footnote>
more massive galaxy. The relative contribution of these channels to galaxy mass growth, particularly at early epochs, remains poorly constrained. While galaxy star formation rates (SFRs) have been reliably measured in the past, determining the galaxy merger rate over cosmic time and assessing the relative importance of these growth mechanisms is considerably more challenging. In this paper, we present the most extensive study to date about the importance of mergers in the mass growth of galaxies in the first 2 Byr of cosmic history.

In addition to understanding the mass build-up of galaxies, mergers were first proposed to explain observed trends in the physical properties of galaxies. In their seminal paper, Toomre & Toomre (1972) identified two spiral galaxies in the process of merging by observing morphological disturbances and tidal tails. Since then, mergers have been considered crucial in the structural evolution and morphological transformation of massive elliptical galaxies (Barnes & Hernquist 1996; Bell et al. 2006; Naab, Khochfar & Burkert 2006; Bournaud et al. 2011). Gas-rich mergers have been found to trigger starburst events (exceptionally high rates of star formation), and merging galaxies exhibit enhanced SFRs compared to isolated ones (Mihos & Hernquist 1994; Patton et al. 2011, 2013, 2020; Torrey et al. 2012; Lanz et al. 2013; Moreno et al. 2015, 2019; Pearson et al. 2019; Garduño et al. 2021; Ellison et al. 2022; Thorp et al. 2022; Montenegro-Taborda et al. 2023; Duan et al. 2024b; Reeves & Hudson 2024; Yuan et al. 2024). Mergers can also trigger active galactic nuclei (AGNs) activity (Silk & Rees 1998; Hopkins et al. 2008; Ellison et al. 2011; Satyapal et al. 2014; Gao et al. 2020; Bickley et al. 2023; Li et al. 2023; Duan et al. 2024b; Sharma et al. 2024), and the most luminous AGNs, ultra-luminous infrared galaxies, and extreme emission-line galaxies are often associated with mergers (Kartaltepe et al. 2010; Ellison et al. 2013; Gupta et al. 2023; Marshall et al. 2023; Perna et al. 2025). Moreover, Witten et al. (2024) recently discovered several Ly$\alpha$ emitters with close companions at $z > 7$, concluding that mergers may drive Ly$\alpha$ emission in these early systems and facilitate the escape of Ly$\alpha$ photons from ionized bubbles (see also Saxena et al. 2023; Witstok et al. 2024). This intense star formation activity in mergers is subsequently quenched (Hopkins et al. 2008; Toft et al. 2014; Ellison et al. 2022), leading to quiescent galaxies characterized by low or suppressed SFRs. This sequence of events ultimately contributes to bulge formation, resulting in the massive elliptical galaxies observed in the local Universe. Recently, it has been found that mergers can also induce bar formation in spiral galaxies at high redshifts (see e.g. Fragkoudi et al. 2025). Galaxy mergers are directly related to the mergers of supermassive black holes, which produce a low-frequency gravitational-wave background that has been measured by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) Collaboration (Agazie et al. 2023). Therefore, it is crucial to accurately measure and constrain the galaxy merger history throughout cosmic time to understand these physical processes better.

There are two main methods to empirically study the fraction of galaxies that are undergoing mergers: (i) counting the galaxies that are in close pairs on the projected plane of the sky and from this estimating the abundance of mergers – *close-pair method* (e.g. Zepf & Koo 1989; Le Fèvre et al. 2000; Patton et al. 2002; López-Sanjuan et al. 2015; Man, Zirm & Toft 2016; Mundy et al. 2017; Mantha et al. 2018; Duncan et al. 2019; Conselice et al. 2022; Duan et al. 2024a), (ii) finding systems that are at the late stages or have recently completed merging based on disturbed morphologies – *morphological method* – that can be further broken down to *quantitative* techniques, i.e. using morphological parameters (e.g. Conselice et al. 2003; Lotz et al. 2008b; Jogee et al. 2009; Desmons

et al. 2023; Rose et al. 2023; Dalmasso et al. 2024), and more recently, using deep learning to identify mergers (e.g. Pearson et al. 2019, 2022; Ferreira et al. 2020; Bickley et al. 2022; Margalef-Bentabol et al. 2024), or *qualitative* methods, for example, visual classifications (e.g. Kartaltepe et al. 2015), and studies that utilize the identification and presence of tidal signatures (e.g. Kado-Fong et al. 2018; Mantha et al. 2019). These two methods complement each other as they probe different phases of galaxy mergers with different observability time-scales. However, one of the main sources of uncertainty in the derived merger rates is the time-scale itself, which makes comparisons between the results obtained by the two different methods difficult and ambiguous. In addition, several other methods exist based on the internal kinematics from spectroscopy (e.g. velocity fields) to identify galaxies involved in mergers or recent merger remnants (e.g. Jesseit et al. 2007; Shapiro et al. 2008). Several JWST/NIRSpec IFU works, mostly as part of the Galaxy Assembly with NIRSpec IFS (GA-NIFS) survey, have found and studied mergers and multiple companion groups, both in star-forming galaxies (SFGs; e.g. Arribas et al. 2024; Jones et al. 2024; Lamperti et al. 2024; Marconcini et al. 2024; Rodríguez Del Pino et al. 2024; Scholtz et al. 2025) and in AGNs (e.g. Perna et al. 2023, 2025; Marshall et al. 2024). These are studies on individual systems with spatially resolved spectroscopy in 2D regions of typically $r < 10$ kpc, using a methodology rather different from this work. In the following text, we will only focus on the close-pair methodology to constrain the galaxy merger rate.

When determining merger rates by the close-pair method, the main galaxy and its companion(s) have to satisfy some selection criteria. The typical criteria are that galaxies have to be at the same (or similar) redshift, which translates to being at a close radial distance from each other and within some projected 2D physical separation (usually in the range 20–50 kpc) on the sky to be close-pairs. In observational studies, the challenge lies in constraining these potential pairs in the third spatial coordinate, i.e. requiring them to be close in redshift. In spectroscopic surveys this corresponds to a maximum velocity offset of ~500 km s$^{-1}$. However, few large spectroscopic surveys exist that probe deep enough to measure redshifts for a mass-complete sample at a meaningful scale, especially at high redshifts. Instead, most works resort to photometric redshifts from large deep photometric surveys and define a similar criterion in redshift space (although less restrictive due to the uncertainty in photometric redshifts), usually defined as $\delta z/(1+z) \simeq 0.01$ (e.g. Molino et al. 2014; Man et al. 2016), which translates into a velocity uncertainty of 3000 km s$^{-1}$. However, only using the peaks of the photometric redshift posterior distributions and their $1\sigma$ uncertainties (e.g. Man et al. 2016; Mantha et al. 2018) leads to information loss about the exact shape of the probability distribution function (PDF) and could lead to unreliable results. The method developed by López-Sanjuan et al. (2015) therefore incorporates the full posterior distribution of the photometric redshifts and propagates the associated uncertainties throughout the full close-pair analysis (see also Mundy et al. 2017; Duncan et al. 2019; Conselice et al. 2022). In this work we use the same method to obtain close-pair fractions.

The first observational studies of galaxy close-pair fractions relied on robust spectroscopic data sets but could only study bright and massive galaxies at low redshifts up to $z \sim 1.2$ (e.g. Patton et al. 2002; Lin et al. 2004; Kartaltepe et al. 2007; López-Sanjuan et al. 2012; Xu et al. 2012). All of these results agree that the major galaxy pair fraction increases with redshift, evolving as a power law. To probe merger fractions at higher redshifts, studies had to rely on large photometric surveys, either by identifying morphologically disturbed systems (using CAS and $M_{20}$/Gini parameters e.g. as in
Conselice, Rajgor & Myers 2008) or by looking for close-pairs using photometric redshifts derived from fluxes measured in multiple bands (e.g. Bluck et al. 2009; Williams, Quadri & Franx 2011; Man et al. 2016; Mantha et al. 2018). These latter works mainly utilize deep photometric data using *Hubble Space Telescope (HST)* from large surveys such as the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Koekemoer et al. 2011), and study close-pair fractions in the range of $z \sim 0-3$. Although at low redshifts, all studies agreed that the pair fraction rises, at these intermediate redshifts of $z \sim 2-3$, some find a further increase, a flattening evolution, or in some cases, even a turnover. This shows that there was already a significant uncertainty in close-pair fraction trends found by *HST* at $z \sim 3$. It is important to note that all these works only use the peak values of the derived photometric redshifts, which could be highly affected by uncertainties. Therefore, works propagating the entire probability distributions of photometric redshifts throughout their analysis by the method described in López-Sanjuan et al. (2015) are generally more robust. Although these studies provide more reliable results (Mundy et al. 2017; Conselice et al. 2022) and probe to higher redshifts (up to $z \sim 6$, see Duncan et al. 2019), they still disagree on the trends of the pair fraction and merger rate evolution. Other instruments and surveys also investigated the number of close-pairs at higher redshifts (up to $z \sim 6-7$), such as the Multi Unit Spectroscopic Explorer (MUSE) deep fields (Ventou et al. 2017, 2019), Subaru Hyper Suprime-Cam imaging surveys (Shibuya et al. 2022), or the ALMA [CII] surveys (e.g. ALPINE; Romano et al. 2021). Similar to previous studies, these results do not resolve the disagreement between a steep increase or turnover in pair fractions at high redshifts. It is also worth noting that most of these studies mentioned above use significantly different selection functions and employ different projected separation limits (see Table 2 for a list of a few relevant studies), which makes comparison difficult and might amplify the lack of agreement.

The advent of the *JWST* (Gardner et al. 2006, 2023) opens a new window upon the study of galaxy mergers at high redshifts using deep photometry by the NIRCam instrument (Rieke et al. 2023a), complemented by high-resolution spectra from NIRSpec (Jakobsen et al. 2022; Böker et al. 2023). A high number of recent studies present discoveries at high redshifts about early galaxy formation, and in particular related to galaxy mergers (e.g. Claeysens et al. 2023; Hashimoto et al. 2023; Hsiao et al. 2023, 2024; Jin et al. 2023; Suess et al. 2023; Tacchella et al. 2023; Treu et al. 2023; Alberts et al. 2024; Carnall et al. 2024; de Graaff et al. 2024; Decarli et al. 2024; Duan et al. 2024a, b). A recent study by Duan et al. (2024a) examines close-pair fractions and major merger rates in the redshift range $4.5 \le z \le 11.5$ using the probabilistic method mentioned previously, where they find rising merger rates up to $z \sim 6$ followed by a flattening. In the discussion section, we review the differences and improvements compared to their study.

The aim of this study is to investigate the redshift evolution of close-pair fractions and major galaxy merger rates at the ambiguous and previously underexplored redshift range of $z \sim 3-9$. We use deep photometric and spectroscopic observations of the GOODS-South and GOODS-North fields by the *JWST* Advanced Deep Extragalactic Survey (JADES) collaboration (Eisenstein et al. 2023a). We adopt the statistical method developed by López-Sanjuan et al. (2015), but instead of analysing a luminosity-selected galaxy sample, we perform a sample selection based on stellar masses as in Mundy et al. (2017). We find that the close-pair fraction rises and peaks at $z \approx 5-6$, followed by a turnover at higher redshifts that also scales with the stellar mass of the primary galaxies. In the case of merger rates, we find that there is an increase and a subsequent flattening

beyond $z \approx 6$. The resulting mass accretion rate is significantly (by a factor of $\sim 5$) lower than the star-forming main sequence (SFMS) SFR at these high redshifts. This suggests that mergers are not the dominant channel for the mass growth of galaxies and contribute to about 3–13 per cent of the *ex situ* mass fraction of an average galaxy.

The paper is structured as follows. Section 2 provides an overview of the catalogues and data products used in this study, starting with an overview of the JADES survey. We detail the methodology for obtaining photometric redshifts and compare them with existing spectroscopic surveys, as well as discussing the derived stellar masses and their completeness. Section 3 outlines the process for identifying close pairs, including selection criteria and sample selection based on redshift and stellar mass. We also describe the pair probability function and correct for selection effects. In Section 4, we present results on the pair fraction and its evolution with redshift. Section 5 covers the conversion of observational pair fractions into a physically meaningful major merger rate by assuming a merger observability time-scale. Section 6 discusses our findings, including evolutionary trends, comparisons with cosmological simulations, and the relationship between star formation and merger rates. Finally, Section 7 summarizes our work and highlights the main conclusions.

Throughout this paper, we adopt the AB magnitude system (Oke & Gunn 1983). We use a standard cosmology with $\Omega_m = 0.310$, $\Omega_\Lambda = 0.689$, and $H_0 = 67.66 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Planck Collaboration VI 2020), and we refrain from using the little *h* notation for the Hubble parameter (Croton 2013). Throughout our analysis, we use the *ASTROPY* python package (Astropy Collaboration 2022), and its subpackage *astrology*, where we assume a flat $\Lambda$CDM cosmology with parameters from Planck Collaboration VI (2020).

## 2 DATA

In this section, we present the various data sets utilized to measure the close-pair fractions of galaxies, which we will later convert into major merger rates. We start by discussing the data and footprints from the JADES survey, followed by a description of the estimated photometric redshifts and their quality assessment, as well as the available spectroscopic data. Finally, we explain how stellar masses are calculated for each galaxy using spectral energy distribution (SED) fitting, as well as the overall stellar mass completeness of the different survey regions.

### 2.1 JADES survey

The JADES (Eisenstein et al. 2023a) is the deepest and most extensive extragalactic survey to date, probing the areas of the Great Observatories Origins Deep Survey (GOODS; Giavalisco et al. 2004), covering the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006) in the GOODS-South region, and the GOODS-North region. JADES is a joint GTO program between the NIRCam and NIRSpec GTO teams that consists of NIRCam imaging, NIRSpec spectroscopy, and MIRI imaging. With DR3 and previous data releases (Eisenstein et al. 2023b; Rieke et al. 2023b; Bunker et al. 2024; D'Eugenio et al. 2025), complemented by JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023) and First Reionization Epoch Spectroscopically Complete Observations (FRESCO; Oesch et al. 2023) data, JADES has more than $\sim 175 \text{ arcmin}^2$ area covered by 8–10 photometric filters ($\sim 200 \text{ arcmin}^2$ in F444W) and over 6000 spectra. It probes down to an unprecedented photometric depth of 30.8 AB magnitude in F444W.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Footprint of the JADES GOODS-South and GOODS-North fields from the mosaics as RGB images (using the F090W, F200W, and F444W bands). Both the GOODS-South and GOODS-North field is further divided into deep and medium tiers based on the exposure time in the F444W band, with the deep tier indicated by the orange area with exposure time above the threshold of $T_{\text{exp}}^{\text{GS}} = 32.5$ ks and $T_{\text{exp}}^{\text{GN}} = 9.5$ ks, respectively. The GS-Deep tier corresponds to the JADES DR1 and JOF footprints.</figcaption></figure>

We divide the JADES photometric sample into four different tiers based on the area and exposure time. Both the GOODS-South and GOODS-North fields have different flux depths due to multiple observations and exposure times. Therefore, we divide the footprints by imposing an exposure threshold of $T_{\text{exp}}^{\text{GS}} = 32.5$ ks for GOODS-South, and $T_{\text{exp}}^{\text{GN}} = 9.5$ ks for GOODS-North. In the case of GOODS-South, the area that was observed with longer exposure time than $T_{\text{exp}}^{\text{GS}}$ is defined as GS-Deep and corresponds to the JADES footprint from DR1 (Rieke et al. 2023b) and the JADES Origins Field (JOF; Eisenstein et al. 2023b) and it is indicated on Fig. 1. The remaining area in the GOODS-South field is defined as GS-Medium, which corresponds to an exposure time lower than $T_{\text{exp}}^{\text{GS}}$ and was released as part of DR2 (Eisenstein et al. 2023b). Similarly, GOODS-North is divided into a deep and medium region based on the integration time.

Here, we list the coverage by different photometric bands that are later used in photometric redshift estimates and SED-fitting to obtain stellar masses (see Sections 2.2 and 2.4). We use JADES NIRCam photometry in bands F090W, F115W, F150W, F162M, F200W, F250M, F277W, F300M, F335M, F356W, F410M, and F444W, which are available in both GOODS-South and -North Deep and Medium tiers. This is supplemented by deep medium-band photometry in F162M, F250M, and F300M as part of the JOF observations (Eisenstein et al. 2023b), which is contained in our GS-Deep tier. In addition, we use JEMS (Williams et al. 2023) NIRCam photometry in the F182M, F210M, F430M, F460M, and F480M bands as part of GOODS-South. Both GOODS fields are further complemented by FRESCO (Oesch et al. 2023) and CONGRESS (GOODS-North) F444W NIRCam photometry, obtaining the most extensive coverage in this band. Furthermore, both GOODS fields have ancillary HST photometric observations using the ACS bands F435W, F606W, F775W, F814W, F850LP, and from the IR bands F105W, F125W, F140W, and F160W.

For the photometric catalogues, we are using internal versions v0.9.3 for GOODS-South and v0.9.1 for GOODS-North, most of which were released as part of DR2 (Eisenstein et al. 2023b) and DR3 (D'Eugenio et al. 2025) and is publicly available on Mikulski Archive for Space Telescopes (MAST; https://archive.stsci.edu/hlsp/jades).<fnref n="1" /> A detailed description about the exact catalogue construction can be found in Rieke et al. (2023b), and here we give a brief summary of the main steps as these could significantly affect the number of sources and satellites detected, and hence the pair fractions derived. A detection catalogue of blended sources is first generated using Photutils (Bradley et al. 2024), adopting a threshold of SNR ≥ 1.5. Deblending is performed on the segmentation map utilizing a logarithmically scaled F200W image. For larger segments, further deblending is applied with the `deblend.sources` function from `Photutils`. Satellite sources are identified by running `detect.sources` on the high-pass filtered outer light profiles of extended sources. Compact sources missed or excluded in previous steps are recovered by reapplying source detection to blank regions, using a threshold of SNR = 3.5. The final segmentation map serves as the basis for constructing the photometric catalogue. Object centroids are determined using the windowed position algorithm provided in `Photutils`, with implemented `Source Extractor` (Bertin & Arnouts 1996) methodology, applied to the NIRCam long-wavelength signal image. Elliptical sizes, orientations, and Kron radii (with $K = 2.5$) for each source are derived from the signal image. The photometric catalogues provide flux measurements in multiple apertures for each source, with forced photometry performed at object positions identified during detection. In this analysis, we use Kron apertures applied to images convolved (`KRON_CONV`) to a consistent resolution across photometric bands.

<footnote n="1">¹DOI: 10.17909/8tdj-8n28</footnote>
## 2.2 Photometric redshifts

### 2.2.1 EAZY

We obtain photometric redshifts for every galaxy in the GS and GN fields by using the photometric redshift code EAZY (Brammer, van Dokkum & Coppi 2008). EAZY is a template fitting code that combines galaxy spectra to fit the observed photometric data by searching for the best fit on a redshift grid. The templates and assumptions made for our fits are detailed in section 3.1 of Hainline et al. (2024b). We take $z_a$ as the photometric redshift, which corresponds to the minimum in the $\chi^2$ of the fit. The grid search is performed in the redshift range $z = 0.01–22$ with a step size of $\Delta z = 0.01$.

Here, we briefly summarize the templates used in the EAZY fits. We started with the modified version of the original seven templates from Brammer et al. (2008) to include line emission, and complemented them with a dusty and a high-equivalent width template (Erb et al. 2010). We further supplemented these by seven new JADES templates<fnref n="2" /> (Hainline et al. 2024b) that were optimized for better redshift estimates using mock-galaxy observations from the JAGUAR simulation (Williams et al. 2018). These templates better cover the colour space of JAGUAR galaxies, including blue UV-bright and red dusty galaxies, and were generated by using flexible stellar population synthesis (FSPS; Conroy & Gunn 2010). For a comprehensive description of the EAZY templates used for the fits, we refer the reader to Hainline et al. (2024b).

### 2.2.2 Photometric redshift probability distributions

We obtain the photometric redshift posterior distributions from EAZY by using the $\chi^2(z)$ outputs of the fits. Assuming a constant prior for the redshift, we calculate the posterior distribution by $P(z) = \exp[-\chi^2(z)/2]$ with a normalization of $\int P(z)dz = 1.0$. We note here that generally, $z_a$ does not coincide exactly with $z_{peak}$, which is internally calculated by EAZY as a probability-weighted average redshift. As described in Hainline et al. (2024b), we use $z_a$ as the best photometric redshift in the following analysis. We note here that recent studies found that photometric redshift estimation could be affected by the Eddington bias (see e.g. Serjeant & Bakx 2023; Donnan et al. 2024) that leads to higher estimated values than the true population, therefore we compare our sample with a robust spectroscopic redshift catalogue (see Section 2.3).

### 2.2.3 Odds quality parameter

The odds quality parameter $\mathcal{O}$ is a proxy for the reliability of the photometric redshift fit, and it is a useful measure to select a robust sample of galaxies with accurate photometric redshifts and a low rate of catastrophic outliers. The odds parameter is defined as the redshift PDF integrated over a small region $\pm K(1+z_a)$ (Benítez 2000; Molino et al. 2014) around the best photometric redshift $z_a$,

$$ \mathcal{O} = \int_{-K(1+z_a)}^{+K(1+z_a)} P(z-z_a)dz, \qquad (1) $$

where $K$ is a survey-specific constant. Previous studies have empirically chosen values for $K$ in the range 0.0125–0.05 (e.g. López-Sanjuan et al. 2014, 2015; Molino et al. 2014), which typically depend on the photometric redshift accuracy of the data. For example, Molino et al. (2014) adopt $K = 0.0125$ for medium-band filters,

<figure><img src="image_4.png" /><figcaption>**Figure 2.** One-to-one comparison between the estimated photometric redshifts from EAZY and all the available spectroscopic redshifts from JADES and other surveys as indicated by the colours. This comparison is for objects in the GOODS-South field, including only the highest-quality redshifts and the best coordinate matches between the photometric and spectroscopic sources. There is a good agreement between the two types of redshifts with an outlier fraction of 6.91 per cent. The cloud of outlier point at $z_{spec} \sim 3–4$ below the one-to-one line might be present due to EAZY confusing the Lyman-break with a Balmer-break, and hence estimating the photometric redshifts to be much lower.</figcaption></figure>

while Conselice et al. (2022) use a larger value of $K = 0.05$ in the case of broadband filters. In this work, we choose $K = 0.03$ since we have photometric redshifts that have relatively high accuracy obtained from multiple wide- and medium-band JWST filters and additional HST photometry.

In our initial sample selection, we require that all objects must have an odds parameter $\mathcal{O} \ge 0.3$. This choice of odds cut is consistent with previous studies (e.g. López-Sanjuan et al. 2015; Mundy et al. 2017; Conselice et al. 2022), and ensures that our initial sample has accurate and robust photometric redshifts. This odds quality cut significantly reduces the initial sample (discarding ~40 per cent of the sources in both GS and GN).

## 2.3 Spectroscopic redshift catalogue

To assess the accuracy of the estimated photometric redshifts from EAZY, we compare them to spectroscopic redshifts. We compile a catalogue containing all the available spectroscopic redshifts of objects falling in the JADES footprints. We include spectroscopic redshifts from the JADES NIRSpec observations (Bunker et al. 2024; D'Eugenio et al. 2025), the FRESCO and CONGRESS surveys (F. Sun, private communication; see also Meyer et al. 2024 for $H\alpha$ and Covelo-Paz et al. 2025 for $H\beta+[O III]$ catalogues, and survey paper by Oesch et al. 2023), the MUSE DR2 (Bacon et al. 2023), the MUSE-Wide DR1 (Urrutia et al. 2019), the ASPECS (ALMA SPECTroscopic Survey in the UDF) Large Program (Decarli et al. 2019), and a collection of all publicly available spectroscopic redshifts from legacy observations falling in the CANDELS fields (GDS collection on Fig. 2; N. Hathi, private communication). The catalogue is built by coordinate matching and correcting for any systematic coordinate offsets between surveys. The classification of reliability and quality of the spectroscopic redshifts from these different surveys is converted to a unified system, categorizing them as secure/best (1), less confident (2), and unreliable (3).

<footnote n="2">²Publicly available on Zenodo: https://zenodo.org/records/7996500</footnote>
Our final assembled catalogue of spectroscopic redshifts consists of 5382 sources in the GOODS-South and 2591 in the GOODS-North field, falling within the best category in quality and coordinate match. The assembled spectroscopic redshift catalogue is compared to the matched photometric redshift catalogue from EAZY. We plot the highest quality matched sources for GOODS-South in Fig. 2. We defined the outlier fraction as $|z_{spec} - z_{phot}|/(1 + z_{spec}) > 0.15$. Overall, there is a good agreement between the two redshifts, with an outlier fraction of 6.91 per cent for GS and 13.62 per cent for GN. The higher outlier fraction in case of GN is likely due to the shallower depth compared to GS and to the fewer number of medium photometric bands available, making the EAZY fits less constrained. The average difference between the spectroscopic and photometric redshifts is $\langle z_{spec} - z_{phot} \rangle = 0.026$ for GS. The scatter around the one-to-one relation is measured by the Normalized Median Absolute Deviation (NMAD), which is defined as

$$ \sigma_{\text{NMAD}} = 1.48 \times \text{median} \left( \left| \frac{\delta z - \text{median}(\delta z)}{1 + z_{\text{spec}}} \right| \right), \quad (2) $$

where $\delta z = z_{spec} - z_{phot}$. For the GS sample, this quantity is $\sigma_{NMAD} = 0.065$, which reflects a good agreement between the spectroscopic and photometric redshifts and the robustness of the catalogue. We note that at high redshift ($z > 8$), damped Lyman-alpha (DLA) absorption might bias photometric redshifts upwards (Fujimoto et al. 2023; Finkelstein et al. 2024; Willott et al. 2024; Helton et al. 2024a; Hainline et al. 2024b). However, our sample of interest ($3 \le z \le 9$) is barely affected by this bias.

Beyond this comparison being a useful and necessary check on the photometric redshifts estimated by EAZY, we also compile a catalogue of the best available redshifts, including the spectroscopic redshifts instead of the photometric ones if available. If multiple $z_{spec}$ exists for an individual object, we select the best match and highest quality for our catalogue. If multiple $z_{spec}$ of the same highest quality and match exist for a single source, we take the average of the spectroscopic redshifts if they are in good agreement ($\delta z_{spec} < 0.1$). In very few cases, there is a disagreement between the highest quality spectroscopic redshifts, and we do not include these values in our final catalogue of best available redshifts.

## 2.4 Stellar masses

For our close-pair selection to find major mergers, a crucial parameter for the galaxies considered is their stellar mass. We obtain stellar masses by two different codes, which allows us to perform two independent analyses and compare the results, providing a further check on the robustness of this study.

### 2.4.1 Prospector

To obtain a robust catalogue of stellar masses, we run Prospector (Leja et al. 2017; Johnson et al. 2021) for a subsample that is relevant to this study. Prospector is an SED-fitting code that uses Bayesian inference to estimate galaxy parameters using different priors, employing an FSPS framework (Conroy, Gunn & White 2009; Conroy & Gunn 2010) through python-fsps (Johnson et al. 2024), using Markov Chain Monte Carlo methods through emcee (Foreman-Mackey et al. 2013) to explore the parameter space and determine the probability distributions of various galaxy properties employing dynamic nested sampling by dynasty. The code incorporates models for stellar populations, star formation histories, including non-parametric SFH (Leja et al. 2019), dust attenuation, nebular emission (Byler et al. 2017), and other physical

processes, which allows for robust estimates and uncertainties for parameters such as stellar mass, SFR, metallicity, and dust content (see e.g. Tacchella et al. 2022; Robertson et al. 2023). On the other hand, it is computationally resource-intensive and significantly slower than eazy-py but provides more accurate and robust results.

Therefore, we only fit a subsample of galaxies using Prospector from our initial catalogue (57 863 fits in total for GS and GN), corresponding to the redshift range of $3 \le z \le 9$ above a signal-to-noise ratio (SNR) threshold in the flux of SNR > 3 measured in the F444W band, which significantly reduces the number of objects to be fitted. The SED fitting used in this work is detailed in Simmonds et al. (2024), and we briefly summarize the priors used and their allowed ranges. We assume a Chabrier (2003) initial mass function (IMF) with mass cut-off at 0.1 and 100 M$_{\odot}$, a top-hat distribution (between $\log(M_{\star}/M_{\odot}) = 6$ and 12) for the stellar mass prior, a non-parametric SFH prior (Leja et al. 2019) with eight time bins, a free ionization parameter with an upper limit of $\log\langle U \rangle_{\text{max}} = -1.0$, and we allow the stellar metallicity to be in the range $0.01 - 1 Z_{\odot}$. We assume a two-component dust model (Charlot & Fall 2000; Conroy et al. 2009), which accounts for the differential effect of dust on young stars (< 10 Myr) and nebular emission lines through different optical depths and a variable dust index (Kriek & Conroy 2013). Finally, we assume a constant photometric redshift prior based on the EAZY photo-z estimates, or fixing the redshift at the spec-z if available in our previously compiled catalogue.

We note here that the estimated stellar masses by Prospector are robust; nevertheless, they can be prone to uncertainties due to the different assumptions in the priors. One such source of uncertainty is the choice of the prior SFH (Tacchella et al. 2022; Whitler et al. 2023), for which we choose a non-parametric model (continuity SFH; Leja et al. 2019). In this model, the SFH is divided into eight different SFH bins, with the first and last bins fixed 5 Myr and $z = 20$, respectively, and the remaining time bins divided into equal intervals of $\log_{10}(t_{\text{lookback}})$. The ratios between adjacent time bins are determined by the bursty-continuity prior (Tacchella et al. 2022). In this model, we adopt the Student's t-distribution (with width $\sigma = 0.3$) for fitting for the $\Delta \log(SFR)$ between the neighbouring bins, which prevents sharp transitions (Leja et al. 2017). Furthermore, this model is also supported by Lower et al. (2020), who demonstrate that non-parametric SFHs outperform traditional parametric models and lead to significantly improved and robust stellar masses. Another source of uncertainty can be the choice of IMF that can introduce variations in the resulting stellar masses (Conroy 2013). These uncertainties caused by the different model assumptions would likely introduce systematic offsets in the estimated stellar masses, and since we estimate the completeness and define the mass bins self-consistently, the resulting pair fractions would likely be similar in the offset stellar mass bins.

### 2.4.2 eazy-py

We also obtain the stellar mass of each photometric source by eazy-py (Brammer et al. 2008; Brammer 2021), which is a set of python photometric redshift tools based on EAZY. It is a fast template-fitting SED code that outputs different galaxy parameters beyond only estimating the photometric redshift, such as the stellar mass, SFR, age, luminosity, and UVBJ colours. Although these parameter estimates are often crude and uncertain, the main advantage of this code is its speed compared to other more sophisticated SED-fitting codes, which makes it ideal for our large data set (over 300 000 photometric sources in total).
<table><caption>**Table 1.** Flux limits derived from the 5σ point-source depths measured in the different tiers for GOODS-North and GOODS-South fields, and the area of each footprint calculated from the non-zero pixels in the F444W band.</caption><thead><tr><th>Tier</th><th>5σ flux limit [nJy]</th><th>Flux limit [AB mag]</th><th>Area [arcmin<sup>2</sup>]</th></tr></thead><tbody><tr><td>GN-Deep</td><td>4.34</td><td>29.81</td><td>36.97</td></tr><tr><td>GN-Medium</td><td>6.99</td><td>29.29</td><td>68.29</td></tr><tr><td>GN</td><td>5.75</td><td>29.50</td><td>105.26</td></tr><tr><td>GS-Deep</td><td>1.97</td><td>30.66</td><td>33.59</td></tr><tr><td>GS-Medium</td><td>4.36</td><td>29.80</td><td>65.37</td></tr><tr><td>GS</td><td>3.37</td><td>30.08</td><td>98.96</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stellar mass distribution of our sample in the GOODS-South field, with grey dots representing masses estimated with Prospector. The stellar mass completeness limit is calculated by the method of Pozzetti et al. (2010). The red dots represent the $M_{lim}$ values of the faintest 20 per cent galaxies, while the red circles with black edges are the $M_{min}$ values in each redshift bin (width of $\Delta z = 0.5$), obtained by taking the 90th percentile of the $M_{lim}$ distribution. The black solid curve is the redshift-dependent completeness limit $M_{\star}^{\text{comp}}(z)$ obtained by fitting a logarithmic curve to the $M_{min}$ values. Throughout this work we focus on the redshift range of $z = 3-9$, and the stellar mass range of $\log_{10}(M_{\star}/M_{\odot}) = 8-10$. The red hatched box represents a fixed arbitrary stellar mass bin centred around $\log(M_{\star}/M_{\odot}) = 9.5 \pm 0.5$, which determines the primary galaxy sample (in the example above $N_{\text{primary}} = 308$) around which we search for major companions. The blue hatched box is the corresponding secondary sample (where $N_{\text{secondary}} = 1120$), which represents a larger pool of galaxies that contains potential satellites with up to 1/4 lower mass than the primary galaxies.</figcaption></figure>

These stellar masses are used for an independent analysis of pair fractions using a data set that is less robust, which is discussed in Appendix C.

## 2.5 Stellar mass completeness

It is important to accurately assess the completeness of our survey and different sub-tiers in order to perform a robust close-pair analysis at high redshift, where the incompleteness in stellar mass could potentially significantly affect the trends observed in the evolution of close-pair fractions and major merger rates.

### 2.5.1 5σ point-source flux depth estimates

The first step to assess the stellar mass completeness is to estimate the 5σ point-source flux depths of our different subfields. For this purpose, we use the full mosaics of the GS and GN fields

in the F444W band. We perform aperture photometry using the photutils python package (Bradley et al. 2024) for a fixed 0.3" diameter aperture. We use the segmentation map obtained from Source Extractor (Bertin & Arnouts 1996), perform a binary dilation (expands the regions flagged by the segmentation map to exclude any contamination caused by extended sources), and create a mask for any source pixels based on the segmentation map or any non-observed area. After sigma-clipping, we obtain the final mask to ensure that we only put apertures on the background, which are not affected by any source pixels. We randomly place apertures with a constant diameter of 0.3" where we measure the flux and use the 80 per cent enclosed energy radius to do aperture correction (resulting in an aperture correction factor of 1.52). Finally, we calculate the 5σ flux limit, which corresponds to the point source flux depth of the survey considered. Additionally, while we measure the depth for fixed 0.3 arcsec apertures, the Kron-apertures for the selected initial sample are similar on average, with a median value of ~0.23 arcsec, and therefore the overall effect is small and the choice of 0.3 arcsec is a good approximation for the depth estimates.

We perform this measurement for the four different tiers defined in Section 2.1. We obtain 5σ flux limits of 1.97 nJy (30.66 AB mag), 4.36 nJy (29.80 AB mag), 4.34 nJy (29.81 AB mag), and 6.99 nJy (29.29 AB mag) for the GS-Deep, GS-Medium, GN-Deep, and GN-Medium subfields respectively (see Table 1). These 5σ limits are lower (f fainter) than the values cited in previous JADES release papers (Eisenstein et al. 2023b; Rieke et al. 2023b; Hainline et al. 2024b; D'Eugenio et al. 2025), which is partly due to using different aperture sizes and using more data that was observed after the data releases which increases the flux depths.

### 2.5.2 Completeness analysis

We perform the stellar mass completeness analysis by adopting the method described in Pozzetti et al. (2010). We aim to find a redshift-dependent stellar mass completeness limit $M_{\star}^{\text{comp}}(z)$, above which we essentially have a complete sample in stellar mass and potentially can observe all types of galaxies. In our flux-limited sample, the limiting stellar mass we can observe depends on the redshift and the mass-to-light ratio $M_{\star}/L$. Therefore, we calculate a minimum mass $M_{min}$ for individual redshift bins of width $\Delta z = 0.5$, where we assume that all galaxies can be observed above a typical $M_{\star}/L$.

In each redshift bin, we calculate $M_{min}$ by first calculating the limiting stellar mass $M_{lim}$ of each galaxy. The limiting stellar mass is the mass a galaxy would have at its spectroscopic (or photometric) redshift if its apparent magnitude ($m_{AB}$ measured in F444W) were equal to the limiting magnitude of the survey or subfield considered, which we obtain from the 5σ flux limits. The limiting stellar mass can be obtained by

$$ \log(M_{\text{lim}}) = \log(M_{\star}) + 0.4(m_{\text{AB}} - m_{\text{AB}}^{\text{lim}}), \quad (3) $$

where $m_{AB}^{\text{lim}}$ is the limiting magnitude of the survey considered (measured in the F444W mosaic; see Section 2.5.1). We then look at the resulting distribution of $M_{lim}$ and only include the faintest 20 per cent galaxies in each redshift bin. To obtain the 90 per cent completeness limit $M_{min}$ in each redshift bin, we take the 90th percentile of the $M_{lim}$ distribution. Finally, we obtain the redshift-dependent stellar mass completeness limit $M_{\star}^{\text{comp}}(z)$ by fitting a logarithmic curve to the $M_{min}$ values at different redshifts, as shown on Fig. 3. These completeness limits are different for each subfield, but in general, are in the range $\log_{10}(M_{\star}/M_{\odot}) \approx 7-8$ depending on redshift, with GS-Deep having the most complete and GN-Medium being the least complete out of our four subfields.
# 3 CLOSE-PAIR METHODOLOGY

This section discusses our methodology for finding close pairs of major galaxy mergers at a wide redshift range. We adopt the close-pair selection method developed by López-Sanjuan et al. (2015), where we make use of the full photometric redshift PDF to propagate the associated uncertainties, as opposed to only using the peaks of these distributions. However, instead of selecting a sample based on the flux ratio of galaxies (as in López-Sanjuan et al. 2015), we perform a selection based on the stellar mass ratio, as it was done in Mundy et al. (2017). Although we mainly follow the selection method described in these two works with some minor modifications, we describe this method in full detail below. See Table 2 for examples of close-pair selection methods and criteria applied by other studies.

## 3.1 Close-pair selection criteria

All studies usually define three main selection criteria to find major galaxy mergers. As the name ‘close-pair’ suggests, galaxies have to be within some defined physical separation $r_{max}$ projected on the plane of the sky, which translates to an observed angular separation $\theta_{max}$. To avoid confusing star-forming clumps within the same galaxy with mergers and ensure clearly debленded sources, many studies define a minimum separation $r_{min}$ as well (and corresponding $\theta_{min}$). The second criterion for two galaxies to be in a close pair is close proximity in the radial direction. In spectroscopic studies, this is usually defined as $\Delta v_{max} \le 500$ km s$^{-1}$ in velocity space (e.g. Patton et al. 2000), while in photometric surveys this translates to close proximity in photometric redshift space (e.g. $\Delta z < 0.1$; Man et al. 2016). For finding major mergers, the third requirement is that the stellar mass ratio between the secondary and the primary galaxy is $\mu = M_2/M_1 \ge 0.25$ (e.g. Lotz et al. 2011). For minor mergers, this value is usually defined as $0.1 \le \mu \le 0.25$. Once the number of close pairs is determined, the pair fraction $f_{pair}$ can be easily computed

$$f_{pair} = \frac{N_{pair}}{N_1}, \qquad (4)$$

where $N_{pair}$ is the number of close-pairs, and $N_1$ is the number of primary galaxies within some pre-defined initial sample, e.g. a volume-limited sample of galaxies within a specified mass bin.

To summarize these three main selection criteria for finding close pairs of major galaxy mergers:

(i) Close projected physical separation $r_{min} \le r \le r_{max}$ (with our choice of $r_{min} = 5$ kpc and $r_{max} = 30$ kpc),<fnref n="3" /> translating to close angular separation on the sky $\theta_{min} \le \theta \le \theta_{max}$,

(ii) Close proximity in redshift or velocity space,

(iii) Stellar mass ratio above $\mu \ge 0.25$ to find major mergers.

In the following text, we describe in detail the motivation and implementation of these selection criteria for a large photometric survey, where we include the $P(z)$ distributions of the estimated photometric redshifts.

## 3.2 Initial sample selection

In this section, we describe how we selected the initial galaxy sample for our analysis. First, we limit our sample to the photometric redshift range of $3 \le z \le 9$, which is an underexplored redshift range in the study of galaxy mergers. We choose this lower and upper redshift limit, as the $z < 3$ range has been studied extensively in

previous research on pair fractions and merger rates, and at $z > 9$, our catalogues contain very few objects after pre-selection (less than 0.1 per cent of the $z > 3$ sample), making a meaningful analysis difficult.

We clean our initial sample by requiring a signal-to-noise ratio in the F444W_KRON_CONV photometry (Kron-aperture placed on images that have been convolved to the same resolution) of $SNR \ge 3$ to filter out background noise and false detections (reduction to 52 per cent of the parent photometric sample of ~110 000 objects). We perform an odds cut of $O \ge 0.3$ as described in Section 2.2.3 to only include galaxies with well-constrained photometric redshifts from EAZY (a further reduction to 63 per cent of the sample from the previous step). Furthermore, we remove any known interlopers, such as 9 brown dwarf candidates (Hainline et al. 2024a), and falsely identified objects belonging to diffraction spikes of bright stars (see Section 3.4.3), resulting in a further reduction by ~5 per cent of the remaining sample. Finally, we compute the stellar mass completeness limit for each tier (GS-Deep, GS-Medium, GN-Deep, GN-Medium) as described in Section 2.5.2, and require all primary (and also secondary) objects to be above $M_{\star}^{comp}(z)$.

### 3.2.1 Redshift and stellar mass bins

Following the initial sample selection, which is uniformly applied for the full analysis for each tier, we define different redshift and stellar mass bins, and we compute the pair fractions within each of these bins. We define the redshift bins equally between $3 < z < 9$ with widths of $\Delta z = 1$, i.e. we obtain results at half redshifts. We then define different constant stellar mass bins which do not depend on redshift. Due to the varying number density of objects having stellar masses above the completeness limit at different redshifts, we define bins with non-equal sizes, centred at $\log(M_{\star}/M_{\odot}) = 8.25 \pm 0.25, 8.75 \pm 0.25, 9.50 \pm 0.50$. These bins were chosen such that they are always above $M_{\star}^{comp}(z)$. We select primary galaxies from these stellar mass bins and will search around them for lower mass satellites in the following parts of the close-pair selection. We can immediately determine the secondary sample of galaxies based on the stellar mass from the primary sample, as the stellar mass ratio is fixed for major mergers at $\mu \ge 0.25$. The primary and secondary sample for a given fixed stellar mass bin is visualized in Fig. 3 (the number of primary galaxies within each bin can be found in Table 3).

## 3.3 Pair probability function

In this section, we discuss in detail how we select close pairs from our photometric survey in a probabilistic and statistically robust manner. We build the pair probability function (PPF) throughout this section, which mathematically codifies the selection criteria from Section 3.1.

### 3.3.1 Redshift probability function

As the first step and most crucial feature of this work, we assess the proximity of potential pairs in redshift space, i.e. closeness in radial distance. To perform this, we incorporate in our analysis the full posterior distributions of the estimated photometric redshifts output from EAZY. If there is an available spectroscopic redshift for a galaxy from our compiled $z_{best}$ catalogue described in Section 2.3, we use that instead. To include a $z_{spec}$ in this probabilistic framework, we fit a narrow Gaussian to this redshift assuming a standard deviation of $\sigma_z = 0.01$ and normalization of 1.

<footnote n="3">³For an analysis using different separation criteria see Appendix D.</footnote>
<table><caption>**Table 2.** Summary of the selection criteria of close-pair studies that are frequently referenced in this work, displaying the redshift and stellar mass ranges probed, the projected separation criteria, and the selection method, i.e. either spectroscopic, photometric (using the peak photometric redshifts), or probabilistic (using the full posterior distributions of the photometric redshifts).</caption><thead><tr><th>Study</th><th>Redshift range</th><th>Mass range, log(M<sub>*</sub>/M<sub>☉</sub>)</th><th>r<sub>sep</sub> [kpc]</th><th>Selection method</th></tr></thead><tbody><tr><td>López-Sanjuan et al. (2015)</td><td>[0.4, 1]</td><td>M<sub>B</sub> &le; -19, -19.5, -20<sup>a</sup></td><td>[10, 50] h<sup>-1</sup></td><td>Probabilistic, based on P(z)</td></tr><tr><td>Man et al. (2016)</td><td>[0.1, 3]</td><td>&gt; 10.8</td><td>[10, 30] h<sup>-1</sup></td><td>|z<sub>1</sub> - z<sub>2</sub>|/(1 + z<sub>1</sub>) &lt; 0.1<sup>b</sup></td></tr><tr><td>Mundy et al. (2017)</td><td>[0.005, 3.5]</td><td>&gt; 10, 11</td><td>[5, 30]</td><td>Probabilistic, based on P(z)</td></tr><tr><td>Ventou et al. (2017)</td><td>[0.2, 6]</td><td>&ge; 9.5, &lt; 9.5</td><td>[5, 25] h<sup>-1</sup></td><td>&Delta;v<sub>max</sub> &lt; 500 km s<sup>-1</sup></td></tr><tr><td>Mantha et al. (2018)</td><td>[0, 3]</td><td>&gt; 10.3</td><td>[5, 50]</td><td>(z<sub>1</sub> - z<sub>2</sub>)<sup>2</sup> &le; &sigma;<sub>z,1</sub><sup>2</sup> + &sigma;<sub>z,2</sub><sup>2</sup><sup>c</sup></td></tr><tr><td>Duncan et al. (2019)</td><td>[0.5, 6]</td><td>[9.7, 10.3], &gt; 10.3</td><td>[5, 30]</td><td>Probabilistic, based on P(z)</td></tr><tr><td>Conselice et al. (2022)</td><td>[0, 3]</td><td>&gt; 11</td><td>[5, 30]</td><td>Probabilistic, based on P(z)</td></tr><tr><td>Duan et al. (2024a)</td><td>[4.5, 11.5]</td><td>[8, 10]</td><td>[20, 50]</td><td>Probabilistic, based on P(z)</td></tr></tbody></table>

Notes. $^a$Sample selected by B-band luminosity.

$^b$For $z_1 < 1$, and $< 0.2$ for $z_1 > 1$.

$^c$This study used both photometric and spectroscopic redshifts, and $\sigma$ are the photometric redshift errors.

We define the redshift probability function $\mathcal{Z}(z)$ as the convolution of the two individual redshift posteriors $P_1(z)$ and $P_2(z)$ of a given projected close-pair as

$$ \mathcal{Z}(z) = \frac{2 \times P_1(z) \times P_2(z)}{P_1(z) + P_2(z)} = \frac{P_1(z) \times P_2(z)}{N(z)}, \quad (5) $$

where $N(z)$ is a normalization defined as

$$ N(z) = \frac{P_1(z) + P_2(z)}{2}, \quad (6) $$

which is implicitly constructed such that $\int_0^\infty N(z)dz = 1$. Therefore, the meaning of $\mathcal{Z}(z)$ is the number of fractional pairs for a given close-pair considered at redshift $z$. For a redshift bin, $\int_{z_{\min}}^{z_{\max}} \mathcal{Z}(z)dz$ gives the number of fractional close-pairs for a system within that redshift range considered. Therefore, the total number of fractional pairs for a given system of two galaxies can be computed as

$$ \mathcal{N}_z = \int_0^\infty \mathcal{Z}(z) dz \quad (7) $$

that ranges between 0 and 1. This quantity $\mathcal{N}_z$ essentially relates to the probability of the two galaxies within the system considered being at the same redshift. See Fig. 4 for examples.

The real power of this analysis is that we can propagate the uncertainty associated with the photo-$z$ estimates throughout our calculations by the redshift probability function $\mathcal{Z}(z)$.

### 3.3.2 Angular separation mask

The second step is to select galaxies from the secondary sample that are in close projected physical distance to galaxies from the primary sample. We start by defining physical separations between the primary galaxy and its companion, as this can be related to the dynamics of the merger. Many previous studies defined a constant physical separation ranging between $r = 5$ and 50 kpc $h^{-1}$ (e.g. López-Sanjuan et al. 2015; Man et al. 2016; Mundy et al. 2017; Duncan et al. 2019; Conselice et al. 2022; Duan et al. 2024a). The choice of maximum separation is crucial for the resulting pair fractions and the comparison to previous results, as we expect a higher number of companions within larger separation limits. Therefore, we choose a constant separation of $r_{max} = 30$ kpc for easier comparison, where we refrain from further usage of the dimensionless Hubble parameter $h$ (Croton 2013), as stated before. In the Appendix, we also discuss the physical motivation and effects of a redshift-dependent separation limit. As for the minimum separation, we choose $r_{min} = 5$ kpc, which is motivated by the higher resolving

power of JWST, taking care of deblending issues (i.e. single galaxy with clumps or two galaxies) even at $z \sim 9$. This lower separation criterion is further motivated by the typical sizes being ~0.5–2 kpc for galaxies at these redshifts (Allen et al. 2024; Miller et al. 2024; Morishita et al. 2024; McClymont et al. 2025b; Westcott et al. 2025). For a pair fraction analysis using different separation criteria, see Appendix D.

We convert the projected physical separation to an observable angular separation by using the angular diameter distance. Therefore, the angular separation will be redshift-dependent according to

$$ \theta_{\max}(z) = \frac{r_{\max}}{d_A(z)} \quad \text{and} \quad \theta_{\min}(z) = \frac{r_{\min}}{d_A(z)}, \quad (8) $$

where $d_A(z)$ is the angular diameter distance that depends on redshift and the assumed cosmology. To calculate $d_A(z)$ we use the angular_diameter_distance(z) function from as-tropy.cosmology.

Finally, using the notation of López-Sanjuan et al. (2015), we define an angular separation mask for each potential close-pair as

$$ \mathcal{M}^{\theta}(z) = \begin{cases} 1, & \text{if } \theta_{\min}(z) \le \theta(z) \le \theta_{\max}(z) \\ 0, & \text{otherwise,} \end{cases} \quad (9) $$

which is a binary mask having the same length as the $P(z)$ array from EAZY. For each redshift bin considered, this binary mask will have a window of allowed angular separations.

### 3.3.3 Pair selection mask

The third step is to ensure the correct stellar mass ratio for the pairs considered to find major mergers. We start by defining a mask to select the primary galaxies in a mathematical form corresponding to what is described in Section 3.2.1

$$ S(z) = \begin{cases} 1, & \text{if } M_{\star}^{\text{lim},1}(z) \le M_{\star,1}(z) \le M_{\star}^{\text{max}} \\ 0, & \text{otherwise} \end{cases}, \quad (10) $$

where $M_{\star,1}(z)$ is the redshift-dependent stellar mass of the primary galaxy and

$$ M_{\star}^{\text{lim},1}(z) = \max\{M_{\star}^{\text{min}}, M_{\star}^{\text{comp}}(z)\}, \quad (11) $$

where $M_{\star}^{\text{min}}$ and $M_{\star}^{\text{max}}$ are the lower and upper bounds of the stellar mass bin considered, and $M_{\star}^{\text{comp}}(z)$ is the redshift-dependent stellar mass completeness limit calculated in Section 2.5.2. We note here that in the computation of $M_{\star}^{\text{lim},1}(z)$, we use $M_{\star}^{\text{comp}}(z_{\text{max}})$ at the upper bound of the redshift bin considered $z_{\text{max}}$ since the stellar
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left panel: example of a candidate close-pair system at $z_1 = 4.22^{+0.08}_{-0.13}$ (primary galaxy), where the red curve indicates the photometric redshift posterior distribution $P_1(z)$ of the primary galaxy, and the blue curve represents $P_2(z)$ for the companion galaxy. Both are well-constrained posteriors with unique narrow peaks, where the dashed vertical lines indicate the location of the peaks. The black continuous curve is the cumulative integral of the redshift probability function $\mathcal{Z}(z)$, resulting in a high number of fractional pairs $N_z$ indicated on the figure, meaning that they are at the same redshift with a high probability. The two galaxies also satisfy the stellar mass ratio criterion, with values displayed in the upper left corner, where $\mu$ is the mass ratio. The total contribution of the two galaxies to the pair counts is given by $N_{\text{pair}} = \int w_2(z) \times \text{PPF}(z) dz = 1.09$, which is greater than unity due to a boosting by the correction weights. A cut-out of two galaxies from the mosaic is displayed on the inset, showing that the system is within a small angular separation. Right panel: another example of a close-pair candidate that only passed the close separation criterion by a narrow margin. Although the peak photo-$z$ values are close to each other ($z_a \approx 4.0$), the posteriors are less constrained, and the overlap between them is smaller, resulting in a lower $N_z$ (lower probability of being a real close-pair). The posteriors also have secondary peaks, suggesting that the estimated photometric redshifts are less reliable.</figcaption></figure>

mass completeness limit determined monotonically increases with redshift.

We further note that the stellar mass of the galaxy is expected to evolve with redshift. However, our fitting codes (EASY or Prospector) do not output the covariance of $M_*$ with $z$. Therefore, we simply assume that the stellar mass scales as

$$M_*(z) = M_*^{\text{peak}} \left( \frac{1+z}{1+z_{\text{peak}}} \right)^2, \qquad (12)$$

where $M_*^{\text{peak}}$ corresponds to the stellar mass output from easy-py or Prospector at the peak of the redshift posterior $z_{\text{peak}}$ (more specifically at $z_a$). This approximation is based on the redshift-luminosity relation by assuming a constant mass-to-light ratio.

Now that we have ensured the correct selection function for the primary galaxy, we have to define a pair selection mask for the secondary galaxy, satisfying the major merger criterion. We define the pair selection mask as

$$M^{\text{pair}}(z) = \begin{cases} 1, & \text{if } M_*^{\text{lim},1}(z) \le M_{*,1}(z) \le M_*^{\text{max}} \\ & \text{and } M_*^{\text{lim},2}(z) \le M_{*,2}(z) \le M_*^{\text{max}} \\ 0, & \text{otherwise,} \end{cases} \qquad (13)$$

where $M_{*,2}(z)$ is the redshift-dependent stellar mass of the secondary galaxy and

$$M_*^{\text{lim},2}(z) = \max\{\mu M_{*,1}(z), M_*^{\text{comp}}(z)\}, \qquad (14)$$

where the criterion for the mass ratio $0.25 \le \mu \le 1$ is imposed to search for major mergers.

### 3.3.4 Pair probability function

We are now equipped with the mathematical forms of the three selection criteria to combine them into a single function that can be used to select close pairs in a probabilistic way, which is the strength of this analysis. We define the PPF as

$$\text{PPF}(z) = \mathcal{Z}(z) \times \mathcal{M}^{\theta}(z) \times \mathcal{M}^{\text{pair}}(z). \qquad (15)$$

The integral of this function over the full redshift range gives the probability of a system of two galaxies being in a major merger close-pair, $\int_0^\infty \text{PPF}(z) dz$. Two examples of close-pair candidates at redshift $z \sim 4$ selected by the above method can be seen in Fig. 4.

## 3.4 Correction for selection effects

In this section, we address the different potential selection biases that affect the selection function we defined previously. We assign statistical weights to correct for each of these effects and include this in the final calculation to determine the close-pair fractions.

### 3.4.1 Mass incompleteness

In the case of primary galaxies that are close to the stellar mass completeness limit, the mass range to look for secondary galaxies will be more strongly affected by incompleteness. Mathematically, this happens when $\mu M_{*,1}(z) \le M_*^{\text{comp}} \le M_{*,1}(z)$. This is due to the requirement in equation (13) that the $M_{*,2}(z)$ has to be above the completeness limit. Therefore, there is a lower chance of finding satellites around a primary galaxy, and the results will be biased.

To account for this selection bias, we assign statistical weights to the secondary galaxies to boost their number counts, as we
might potentially miss some of them due to the incomplete stellar mass search range available. Similarly, we assign a weight to the primary galaxies close to the stellar mass completeness limit to reduce their number counts and, hence, their biasing effect. To compute these weights, we compare the expected and the observable number density of galaxies by using the galaxy stellar mass function (SMF), $\phi(M_*, z)$. We combine the most recent SMFs available from JWST observations (Navarro-Carrera et al. 2024; Harvey et al. 2025) and extrapolate it to the highest relevant redshifts, as described in Appendix B. Each secondary galaxy gets the following completeness weight assigned

$$w_2^{\text{comp}}(z) = \left[ \frac{\int_{M_{\star,1}(z)}^{M_{\star,1}(z)} \phi(M_{\star}, z) dM_{\star}}{\int_{\mu M_{\star,1}(z)}^{M_{\star,1}(z)} \phi(M_{\star}, z) dM_{\star}} \right]^{-1}, \quad (16)$$

which is the inverse of the ratio between the expected number density of observable galaxies (above the mass completeness limit) and the total number density of galaxies corresponding to the mass search range around the primary galaxy at a given redshift. Similarly, for each primary galaxy, we assign the weight

$$w_1^{\text{comp}}(z) = \frac{\int_{M_{\star,1}^{\text{lim},1}(z)}^{M_{\star,1}^{\text{max}}(z)} \phi(M_{\star}, z) dM_{\star}}{\int_{M_{\star,1}^{\text{min}}}^{M_{\star,1}^{\text{max}}} \phi(M_{\star}, z) dM_{\star}}, \quad (17)$$

which compares the expected number density of galaxies above the completeness limit to the number density for the full stellar mass bin.

As Fig. A1 shows, these mass completeness weights can drastically increase, approaching the mass completeness limit, with values blowing up at the limit due to the ratio in their definition. To avoid over-correcting our results, we clip the values of $w_2^{\text{comp}}$ at 10.

### 3.4.2 Photometric redshift quality

The next weighting is incorporated to account for the quality of the photometric redshifts. This is realized through the odds sampling rate (OSR), which was originally introduced in López-Sanjuan et al. (2015).

The odds quality parameter $\mathcal{O}$ (Benítez 2000; Molino et al. 2014) was already discussed in Section 2.2.3, which essentially encodes the probability of a galaxy being found within a narrow redshift interval centred on the peak of the photometric redshift posterior distribution. We also require an odds cut of $\mathcal{O} \ge 0.3$ in our pre-selection. The OSR is simply the ratio between the number of galaxies above this criterion and the total number of galaxies at magnitude $m$, defined as

$$\text{OSR}(m) = \frac{\sum N(\mathcal{O} \ge 0.3)}{\sum N(\mathcal{O} \ge 0)}. \quad (18)$$

We calculate the OSR in magnitude bins of width $\Delta m = 0.25$ and fit a sigmoid curve to the data points to find a continuous function for OSR($m$) at each value of $m$. For this purpose, we convert the fluxes of each galaxy measured in the F444W filter to AB magnitudes and get an OSR($m$) curve such as in Fig. A2.

Finally, to compute the photometric redshift quality weight corresponding to each galaxy, we take the inverse of the OSR value as defined below

$$w^{\text{OSR}} = \frac{1}{\text{OSR}(m)}, \quad (19)$$

where $m$ is the AB magnitude of the galaxy in the F444W filter.

### 3.4.3 Survey boundaries

Finally, we assign weights to correct for incomplete survey areas. Primary galaxies that lie close to the survey boundaries or masked regions (where artefacts were masked out during the data reduction) will have a reduced area around them to look for potential secondary galaxies. This effect does not significantly alter the final results, nevertheless, we correct for it.

Furthermore, we note here that additional masks were created for bright stars and their diffraction spikes. During the initial runs of the merger selection code, it was found that a large number of artificial sources along the diffraction spikes of bright stars (due to an erroneous ‘shredded’ source identification by the automatic source extraction) significantly boosted the close-pair fraction at certain redshifts. Therefore, we generously masked these artefacts, which further reduced the size of the search area and altered its shape.

We calculate the area weight by performing photometry on the annulus surrounding each primary galaxy at inner and outer radii corresponding to $r_{\text{min}}$ and $r_{\text{max}}$. In practice, since $\theta_{\text{min}}(z)$ and $\theta_{\text{max}}(z)$ are redshift-dependent through the angular diameter distance $d_{\text{A}}(z)$, the fraction of available area, $f_{\text{area}}(z)$ will be also redshift-dependent and has to be calculated at each redshift-step. We use the `photutils` python package (Bradley et al. 2024) to compute the sum over the area within a fixed annulus around a source (see for an example Fig. A3), containing pixels with non-zero values and the sum over the same area including pixels with zero flux which is equivalent to the masked area. The ratio of these two values gives the area fraction $f_{\text{area}}(z)$. To correct for the unavailable area and potentially missed secondary galaxies, we assign an area weight to each primary galaxy within our sample

$$w_{\text{area}}(z) = \frac{1}{f_{\text{area}}(z)}. \quad (20)$$

In practice, we find that the overall effect of these weights is small, on average being 1–2 per cent (see also e.g. Duncan et al. 2019). Therefore, to make our merger selection code more efficient, we only evaluate these area weights at the best available redshift (see Section 2.3) of each primary source. Hence, we assume that $w_{\text{area}}(z) = w_{\text{area}}$ for the rest of the analysis.

### 3.4.4 Combined weightings

Ultimately, we combine all the weightings mentioned above into single weights assigned to each primary and secondary galaxy. We then include these weights, accounting for different selection biases, in integrating the pair probability functions to calculate close-pair fractions.

In the case of each primary galaxy, two of the aforementioned weights are combined, and the total weight is given by

$$w_1(z) = w_1^{\text{comp}}(z) \times w_1^{\text{OSR}}. \quad (21)$$

For each secondary galaxy around a primary galaxy, we need to convolve the total weight of the primary with the weights corresponding to the secondary (as its selection is dependent on the former), getting a total combined pair weight of

$$w_2(z) = w_1^{\text{area}} \times w_1^{\text{comp}}(z) \times w_2^{\text{comp}}(z) \times w_1^{\text{OSR}} \times w_2^{\text{OSR}}, \quad (22)$$

where we include $w_1^{\text{area}}$ only for the secondary galaxy as it accounts for potentially missing satellites around a primary source, with the subscript ‘1’ further stressing this point.
We further note that the redshift-dependence of $w_1(z)$ and $w_2(z)$ arises due to the redshift-dependence of the mass completeness weights, which further depends on the previously calculated mass completeness limit. We calculate this on the same redshift grid as the posterior output of EAZY for consistency. The most significant contribution usually comes from $w_2^{\text{comp}}(z)$. Hence, a clipping is necessary to avoid over-boosting with this weight, as described previously.

These weights are combined with the pair probability function for each potential pair, allowing us to measure volume-limited pair fractions.

### 3.5 Final integrated pair fractions

Finally, using the results from all previous steps, the total integrated pair fractions can be calculated as follows. For each primary galaxy $i$ from the initial sample, the number of associated close-pairs in the redshift range $z_{\text{min}} \le z \le z_{\text{max}}$ (within each redshift bin) can be computed by

$$N_{\text{pair}}^{i} = \sum_{j} \int_{z_{\text{min}}}^{z_{\text{max}}} w_{2}^{j}(z) \times \text{PPF}_{ij}(z) dz, \quad (23)$$

where $j$ are the indices of each potential secondary galaxy around the primary that satisfies the pre-selection criteria, PPF$_{ij}$(z) are the pair probability functions for each pair considered given by equation (15) (see Section 3.3.4), and $w_2^j(z)$ are the pair weights associated with each secondary galaxy given by equation (22).

Similarly, the contribution to the number of primary galaxies for each system $i$ within a $z$-bin corresponding to $N_{\text{pair}}^i$, is given by<fnref n="4" />

$$N_{1}^{i} = \int_{z_{\text{min}}}^{z_{\text{max}}} w_{1}^{i}(z) \times P_{i}(z) \times S_{1}(z) dz, \quad (24)$$

where $P_i(z)$ is the normalized redshift posterior distribution function of the primary galaxy (see Section 2.2), $w_1^i(z)$ is the corresponding weight given by equation (21), and $S_1(z)$ is a selection function for the primary galaxy sample. This selection function incorporates the initial selection criteria described in Section 3.2, and its mathematical form is given by equation (10) (note that we added an extra ‘1’ index to emphasize that it is the primary galaxy selection function).

Finally, we have derived all the necessary quantities to redefine and calculate the galaxy close-pair fraction (given by equation (4)). The pair fraction for a redshift-bin of $z_{\text{min}} \le z \le z_{\text{max}}$, using the probabilistic methodology described above, has the following form

$$f_{\text{pair}} = \frac{\sum_{i} N_{\text{pair}}^{i}}{\sum_{i} N_{1}^{i}}, \quad (25)$$

and for the sake of completeness, the full form of the pair fraction is given by

$$f_{\text{pair}} = \frac{\sum_{i} \sum_{j} \int_{z_{\text{min}}}^{z_{\text{max}}} w_{2}^{j}(z) \times \text{PPF}_{ij}(z) dz}{\sum_{i} \int_{z_{\text{min}}}^{z_{\text{max}}} w_{1}^{i}(z) \times P_{i}(z) \times S_{1}(z) dz}. \quad (26)$$

In the case of different subfields (i.e. GS-Medium, GS-Deep, GN-Medium, GN-Deep) the total pair fraction can be calculated by

$$f_{\text{pair}} = \frac{\sum_{k} \sum_{i} N_{\text{pair}}^{i,k}}{\sum_{k} \sum_{i} N_{1}^{i,k}}, \quad (27)$$

where the index $k$ corresponds to the subfield, and the stellar mass completeness limit $M_*^{\text{comp}}(z)$ is calculated separately for each subfield.

### 3.6 Bootstrapping analysis

We estimate the uncertainties on the $f_{\text{pair}}$ values using the common bootstrapping technique (Efron 1979, 1981). The merger selection code estimates the pair fractions and their errors for each stellar mass–redshift bin by randomly drawing a primary sample of galaxies with replacement from the initial sample and performing for $N$ independent realizations. We set $N = 100$ for most of our analysis as we find that the uncertainties relatively quickly converge after a few resamples.

The standard error on the pair fraction in each bin from the bootstrapping analysis is given by

$$\sigma_{f_{\text{pair}}} = \sqrt{\frac{\sum_{i} (f_{\text{pair}}^{i} - \langle f_{\text{pair}} \rangle)^{2}}{N - 1}}, \quad (28)$$

where $\langle f_{\text{pair}} \rangle$ is the average pair fraction of the $N$ independent realizations given by

$$\langle f_{\text{pair}} \rangle = \frac{\sum_{i} f_{\text{pair}}^{i}}{N}. \quad (29)$$

All the above definitions and calculations are incorporated in the merger selection python code. We measure the pair fraction in this probabilistic approach for a large sample of ~110 000 galaxies (size of the original photometric catalogue before sample cleaning) for both GOODS-South and GOODS-North. For a few examples of close-pairs found by the code, see Fig. 5, which shows close-pairs with a high probability of being actual mergers ordered in a Toomre sequence, and Fig. 6 for examples of peculiar multimerger systems. Fig. 6 System B is an extended clumsy object where the failed deblending results in identifying a spurious close-pair system with a high multiplicity. A visual inspection of 500 randomly selected close-pairs from our final catalogue shows that approximately ~3 per cent of the sample may be affected by deblending issues. This percentage is likely even lower regarding the number of galaxies affected since these pairs often have a high multiplicity. However, these problems mainly occur at lower redshifts in extended objects, where the number counts are sufficiently high to avoid significant boosting by these spurious pairs.

## 4 PAIR FRACTION

In this section, we discuss the results of the probabilistic close-pair fraction analysis. We present our fiducial results for the previously defined data set and selection criteria and compare them to the existing literature. We also investigate the effect on the resulting pair fractions using different data sets and varying some parameters in the selection criteria. In the next section, we convert the pair fractions to the merger rate by assuming a merger time-scale.

### 4.1 Major merger pair fraction evolution with redshift

In our fiducial analysis, we use stellar masses obtained from Prospector and a redshift-independent separation criterion. See Appendix C for an alternative analysis using stellar masses from eazy-py, and Appendix D for a redshift-dependent separation selection criterion.

Pair fractions are calculated in the redshift range $z = 3–9$ and stellar mass range $\log(M_*/M_\odot) = 8.0–10.0$ and physical projected

<footnote n="4">⁴We note that many previous publications include a $\sum_i$ in front of this equation. We consider this to be a typo as it leads to mathematically erroneous results – in the final form of the pair fraction, there would be a double summation over $i$ for $N_1^i$.</footnote>
<figure><img src="image_7.png" /><figcaption>Merging sequence of close galaxy pairs at z~3</figcaption></figure>

**Figure 5.** Close galaxy pairs with a high probability of being actual mergers ($N'_z > 0.8$) at $3 < z < 4$. The RGB cutouts (with blue, green, and red assigned to the F090W, F200W, and F444W bands, respectively) are ordered by decreasing projected separation between the close pairs (from the top-left to the bottom-right) in an attempt to build a Toomre sequence (Toomre & Toomre 1972) of merging galaxies at $z > 3$. The primary galaxies are indicated by the longer half-crosshairs, and the secondary galaxies by the shorter crosshairs. The redshifts of the primary galaxies are shown on the top-right corner of each cutout (with the redshift of the secondary galaxies being at the same or similar values). The angular and proper physical scale (calculated at $z = 3.5$) is displayed on the top-left cutout, and all the other cutouts have the same scale (with some minor differences in the proper physical scale due to the deviation from $z = 3.5$). The order of these close pairs in the Toomre sequence is both determined by the projected separation and by visual inspection, where one might see tidal features and bridges forming on the bottom row of the cutouts.

separation between $r = 5–30$ kpc, as defined in Section 3.2.1. Results are plotted in Fig. 7 for the three stellar mass bins, and values can be found in Table 3. Pair fractions in the lowest stellar mass bin ($8.0 \le \log(M_*/M_\odot) < 8.5$) show the most significant evolution. There is an

increase up to $z \sim 6$, where 37 per cent of the galaxies have major companions, and it subsequently turns over to a declining trend. Pair fractions in the higher stellar mass bins show little evolution with redshift, with potentially a slight decrease. We note that data points
<figure><img src="image_8.png" /></figure>

**Figure 6.** Multimerger systems at $z > 3$ with a high probability ($N_z > 0.75$), showing an abundant variation of morphologies and multiplicities, reflecting the complexity of these major mergers at high redshift and the difficulty of finding them. The primary galaxy (the most massive galaxy considering all pairs in the system) is indicated by a white circle on each RGB cutout, with half-crosshairs pointing to the secondary galaxies. The redshifts of the primaries are displayed in the top-right corner, as well as the proper physical scale at that redshift for each cutout. *Systems A, C, and D* have multiple satellites around each primary galaxy. However, they likely have even more companions but with a lower probability than our selection threshold (and hence not displayed) due to less reliable photometric redshifts that spectroscopic observations could confirm in the future. *System C* falls within the FRESCO footprint, and therefore, all candidates have spectroscopic redshifts from grism observations. The bright galaxy at the lower left of the cutout is likely to be also a member, but it was not selected by the code as having a different spec-z from higher priority MUSE Wide observations (see Section 2.3 for the spec-z classification). *System B* is a peculiar case where we only display the primary galaxy, as the companions identified by the close-pair selection code are likely clumps associated with the same structure. This extended clumpy morphology causes deblending to fail, which in turn translates to likely false close-pairs. These few peculiar systems (likely only ~3 per cent of the sample) would require further morphological analysis for confirmation, which is beyond the scope of this paper.

<figure><img src="image_9.png" /></figure>

**Figure 7.** Redshift-evolution of the major galaxy close-pair fractions for three different stellar mass bins. Uncertainties are obtained by bootstrapping with 100 resamples. The lowest stellar mass bin shows the most clear evolution by increasing up to $z \sim 6$ and then turning over at higher redshifts. Pair fractions at higher stellar mass bins show little to slightly decreasing evolution with redshift overall. Each data point represents the pair fraction in a given redshift and stellar mass bin for the overall sample, including all subfields (e.g. GS and GN, deep and medium tier).

represent the total pair fractions in a given bin for the overall sample, which is obtained by equation (27) and not by averaging the values obtained by each subfield. Bins at high redshifts (above $z \sim 7–8$) can become unreliable as the number counts of primary (massive) galaxies become very low and are prone to variations, especially at the higher stellar mass bins (see Table 3). Therefore, we decided to discard the data points at the highest redshifts at $z \sim 8–9$ for the

$\log(M_*/M_\odot) = 8.5–9.0$ and $\log(M_*/M_\odot) = 9.0–10.0$ stellar mass bins as it can be seen in Fig. 7 (see Fig. 8 for measurements in the individual tiers for these bins).

We also look at the field-to-field variation of the pair fractions in each stellar mass bin, which can be seen in Fig. 8. The most significant variation can be found in the $\log(M_*/M_\odot) = 8.0–8.5$ stellar mass bin, with the GS-Medium results differing from the other three the most. We note that, for the $\log(M_*/M_\odot) = 8.0–8.5$ bin, the GS-Medium results drive the peak of $f_P$ to be at $z \sim 5–6$, and without GS-Medium the peak becomes less evident. This higher close-pair count could be driven by the presence of multiple overdensities at $z \sim 5–6$, discovered by Helton et al. (2024a, b). Since this lowest stellar mass bin is the closest to the mass-completeness limit, we expect that it is still somewhat affected by incompleteness, even after using the appropriate correction weights, defined in Section 3.4.1. While the statistical uncertainties estimated by bootstrapping do not reflect this, there could still be a systematic uncertainty affecting our completeness estimates due to a higher time-variability of star formation in low-mass and/or high-redshift galaxies (see e.g. Sun et al. 2023). In our stellar mass completeness estimates, based on the method of Pozzetti et al. (2010), we implicitly assume the same colours for the undetected galaxy population below the completeness limit, which might be strongly affected by the higher burstiness in star formation of low-mass galaxies at high $z$. Furthermore, based on these plots, we conclude that our samples in the different subfields are not significantly affected by cosmic variance.

## 4.2 Comparison to literature

In this section, we compare our results for the pair fractions with those of the existing literature. This is a difficult task as most other works use different stellar mass bins or separation limits or both (e.g. see Table 2). Nevertheless, we attempt to find results from the literature
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Redshift-evolution of the close-pair fractions for three different stellar mass bins, represented by the three subplots. Each subplot displays the field-to-field variation within each bin for the four subfields (e.g. GS and GN, deep and medium tier), with the lowest stellar mass bin showing the most variations. Overall, cosmic variance does not seem to significantly affect the resulting pair fractions. The averages here represent the overall pair fractions as calculated from the union of the four subfields (equation 27) instead of just calculating the mean. For the two highest stellar mass bins we do not show the averages at the highest redshift bin at $z \sim 8-9$ as the size of the initial sample is not sufficiently large to allow us measuring pair fractions robustly.</figcaption></figure>

<table><caption>**Table 3.** Number of primary galaxies ($N_{primary}$), and the measured major merger pair fraction ($f_p$) and merger rate ($\mathcal{R}_M$) values at each respective redshift and stellar mass bin, with their associated uncertainties. Uncertainties for pair fractions are estimated using bootstrapping (using 100 resamples) and carried forward to calculate uncertainties for merger rates.</caption><thead><tr><th>Redshift</th><th>[3.0, 4.0]</th><th>[4.0, 5.0]</th><th>[5.0, 6.0]</th><th>[6.0, 7.0]</th><th>[7.0, 8.0]</th><th>[8.0, 9.0]</th></tr></thead><tbody><tr><td>Mass range (log M<sub>*</sub>/M<sub>☉</sub>)</td><td colspan="6" style="text-align: center;">Number of primary galaxies, N<sub>primary</sub><sup>a</sup></td></tr><tr><td>[8.0, 8.5]</td><td>3655</td><td>2842</td><td>2170</td><td>1523</td><td>501</td><td>234</td></tr><tr><td>[8.5, 9.0]</td><td>2295</td><td>1397</td><td>662</td><td>350</td><td>95</td><td>37</td></tr><tr><td>[9.0, 10.0]</td><td>1613</td><td>760</td><td>312</td><td>104</td><td>34</td><td>10</td></tr><tr><td>Mass range (log M<sub>*</sub>/M<sub>☉</sub>)</td><td colspan="6" style="text-align: center;">Pair fraction, f<sub>P</sub></td></tr><tr><td>[8.0, 8.5]</td><td>0.19 ± 0.10</td><td>0.31 ± 0.15</td><td>0.37 ± 0.17</td><td>0.34 ± 0.16</td><td>0.19 ± 0.09</td><td>0.11 ± 0.05</td></tr><tr><td>[8.5, 9.0]</td><td>0.20 ± 0.10</td><td>0.21 ± 0.11</td><td>0.12 ± 0.06</td><td>0.14 ± 0.07</td><td>0.08 ± 0.04</td><td>0.02 ± 0.01<sup>b</sup></td></tr><tr><td>[9.0, 10.0]</td><td>0.14 ± 0.06</td><td>0.12 ± 0.07</td><td>0.09 ± 0.04</td><td>0.12 ± 0.05</td><td>0.06 ± 0.02</td><td>0.09 ± 0.05<sup>b</sup></td></tr><tr><td>Mass range (log M<sub>*</sub>/M<sub>☉</sub>)</td><td colspan="6" style="text-align: center;">Merger rate, R<sub>M</sub> (Gyr<sup>-1</sup>)</td></tr><tr><td>[8.0, 8.5]</td><td>1.60 ± 0.80</td><td>3.91 ± 1.92</td><td>6.46 ± 3.05</td><td>8.08 ± 3.66</td><td>5.64 ± 2.64</td><td>4.04 ± 1.88</td></tr><tr><td>[8.5, 9.0]</td><td>1.70 ± 0.85</td><td>2.67 ± 1.33</td><td>2.20 ± 1.02</td><td>3.29 ± 1.56</td><td>2.49 ± 1.30</td><td>0.57 ± 0.33<sup>b</sup></td></tr><tr><td>[9.0, 10.0]</td><td>1.18 ± 0.52</td><td>1.52 ± 0.85</td><td>1.62 ± 0.77</td><td>2.89 ± 1.14</td><td>1.75 ± 0.68</td><td>3.38 ± 2.01<sup>b</sup></td></tr></tbody></table>

Notes. <sup>a</sup>We note here that these values reflect the number of primary galaxies falling into each redshift-stellar mass bin based on their stellar mass and peak photometric (or spectroscopic) redshift ($z_a$) falling within the bin limits. For the purposes of the pair fraction measurement, the total numerical contribution of the primary galaxies in each bin is given by $\sum_i N_i^1$, where $N_i^1$ is given by equation (24).

<sup>b</sup>The values presented in the table above are included for the sake of completeness. However, we do not plot or use them in our further calculations, as the number counts from the initial sample selection in these bins are too low to perform a reliable analysis.

whose selection criteria are the closest to our choices and compare them to our results. In the Appendix, we discuss multiple variations of selection criteria for the fiducial analysis and their effects on the results.

For comparing results from the existing literature to our findings, we attempt to keep the distinction between the three stellar mass ranges. While this categorization is not always perfectly clear and well defined in the source papers, we nevertheless try to associate data points with one of the three stellar mass bins (e.g. if a data point is slightly above $log_{10}(M_*/M_{\odot}) = 10$, we will include it in the $9 \le log_{10}(M_*/M_{\odot}) < 10$ bin).

The final results for pair fractions and their comparison to values from the literature can be found in Fig. 9. The values associated with the three stellar mass bins are colour-coded, and uncertainty bands for results from this work are plotted but omitted in the case of literature results to avoid overcrowding the figure. The blue upward-pointing

triangles represent data from Conselice et al. (2003) for morphologically selected $M_* > 10^8$ M$_{\odot}$ major mergers at $z \le 3$, while the blue downward-pointing triangles data from Conselice et al. (2008) for galaxies selected by a similar method in the $10^8 - 10^9$ M$_{\odot}$ range. The blue circles represent results from a spectroscopic close-pair study by de Ravel et al. (2009) for a sample of $M_B \le -18$ galaxies at $z < 1$, that approximately corresponds to $> 10^8$ M$_{\odot}$. The blue, orange, and green diamonds correspond to measurements from Casteels et al. (2014) that match our stellar mass bins for $z < 0.2$ galaxies, selected by a morphological method. Orange and green squares represent measurements from deep MUSE observations by Ventou et al. (2017) extending up to $z \approx 6$. The green downward-pointing triangles show results from Mundy et al. (2017) for $M_* > 10^{10}$ M$_{\odot}$ galaxies with $r = 5-30$ kpc separations. The blue and orange rightward-pointing triangles show results from Ventou et al. (2019) of the MUSE deep fields that approximately match our stellar mass bins. The green
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The evolution of major galaxy pair fractions measured in this work and their associated uncertainty bands, colour-coded by stellar mass. Data points from the existing literature (Conselice et al. 2003, 2008, 2022; de Ravel et al. 2009; Casteels et al. 2014; Mundy et al. 2017; Ventou et al. 2017, 2019; Duncan et al. 2019; Dalmasso et al. 2024; Duan et al. 2024a) are plotted with uncertainties omitted to avoid overcrowding the figure, and their specific details are described in the text. Although most of these works use different selection criteria and stellar mass ranges, we attempted to select the closest, most relevant and colour-coded so that they approximately match our stellar mass bins. We find that the pair fractions weakly scale with stellar mass and that pair fractions at higher stellar masses peak at later epochs (i.e. at lower redshifts). We also observe a slight turnover in the pair fraction trends at the highest redshifts.</figcaption></figure>

upward-pointing triangles represent CANDELS data from Duncan et al. (2019) for $10^{9.7} < M_*/M_{\odot} < 10^{10.3}$, extending up to $z \approx 6$. The green circles show results from Conselice et al. (2022) for many different data sets for $M_* > 10^{11} M_{\odot}$. The open black circles show recent results from Duan et al. (2024a) using JWST data for the same overall mass range as this work (hence no colour coding) and similarly for results from Dalmasso et al. (2024) that use morphological parameters to select close-pairs.

We fit a power law + exponential function to all the data points from both our measurements and the literature, which is theoretically motivated by the Press–Schechter formalism for merging galaxies (Carlberg 1990; Conselice et al. 2008) and seems to better describe the pair fraction evolution in the literature. This function has the form

$$f_{\text{P}}(z, M_*) = f_0 \times (1+z)^m \times e^{\tau(1+z)}, \quad (30)$$

where $f_0$, $m$, and $\tau$ are parameters for which we perform weighted non-linear least squares fitting with uncertainties incorporated as weights. The fitted parameters and their uncertainties can be found in Table 4. We find that the fit is worse in the case of the individual stellar mass bins, and therefore, we omit these from Table 4 and do not plot them in Fig. 9.

There is a clear difference between the pair fraction evolution with redshift in the case of the three different stellar mass bins. Based on the measured values (displayed in Fig. 9), the peak of the pair fractions is at lower redshifts for increasing stellar masses of the primary galaxies. The pair fraction trend in the lowest stellar mass bin has the steepest evolution and is dominating at higher redshifts. The opposite can be said about pair fractions in the highest stellar mass bin. There is also a slight turnover in the trends of pair fraction evolution in all three cases.

<table><caption>**Table 4.** Fit parameters for the pair fraction ($f_p$) and merger rate ($R_M$) for three stellar mass bins for an assumed power law + exponential model of the form $C \times (1+z)^m \times e^{\tau(1+z)}$, where C is $f_0$ for the pair fraction, and $R_0$ for the merger rate. The fitting is done using the weighted non-linear least squares method with uncertainties incorporated as weights. We also provide overall fits for the full mass range of $log(M_*/M_{\odot}) = [8, 10]$ considered for all available data assuming the same functional form as before, for both the pair fractions and merger rates. Although this fit is somewhat forced and erases the scaling relation with stellar mass, it is useful and simplifies further calculations.</caption><thead><tr><th>Fit parameters</th><th>C &times; 10<sup>2</sup></th><th>m</th><th>&tau;</th></tr></thead><tbody><tr><td>Mass range (log M<sub>*</sub>/M<sub>&odot;</sub>)</td><td colspan="3">Pair fraction, f<sub>p</sub></td></tr><tr><td>All</td><td>3.71 &plusmn; 0.24</td><td>2.27 &plusmn; 0.20</td><td>-0.48 &plusmn; 0.07</td></tr><tr><td>Mass range (log M<sub>*</sub>/M<sub>&odot;</sub>)</td><td colspan="3">Merger rate, R<sub>M</sub> (Gyr<sup>-1</sup>)</td></tr><tr><td>[8.0, 8.5]</td><td>2.20 &plusmn; 0.20</td><td>4.18 &plusmn; 0.29</td><td>-0.42 &plusmn; 0.07</td></tr><tr><td>[8.5, 9.0]</td><td>2.25 &plusmn; 0.37</td><td>3.36 &plusmn; 0.40</td><td>-0.27 &plusmn; 0.11</td></tr><tr><td>[9.0, 10.0]</td><td>2.61 &plusmn; 0.17</td><td>2.90 &plusmn; 0.20</td><td>-0.20 &plusmn; 0.06</td></tr><tr><td>All</td><td>2.49 &plusmn; 0.16</td><td>2.85 &plusmn; 0.20</td><td>-0.16 &plusmn; 0.06</td></tr></tbody></table>

## 5 MAJOR MERGER RATE

While the galaxy close-pair fraction is a useful quantity to identify mergers, it is a purely observational one that depends on the exact methodology chosen, such as the selection criteria definitions and the stellar mass range being probed. Therefore, it is hard to compare directly to other studies that use slightly different selection parameters. Conversely, the merger rate is a physically more meaningful quantity that is universal (cf. SFR) and can be directly compared with results from different methodologies (such as derived from close-pair fractions or morphological identification). Converting the observed close-pair fractions to merger rates requires the assumption of a merger observability time-scale, which is discussed in the next section.
We define two merger rate quantities that are fundamental to the study of galaxies via mergers. The first is the merger rate per galaxy ($R_M$), which traces the number of mergers occurring per unit time per massive galaxy. The second is the merger rate density ($\Gamma_M$), which is the overall merger rate per comoving volume in cMpc$^3$.

## 5.1 Merger observability time-scale

Converting the inferred pair fractions to merger rates requires the assumption of a merger observability time-scale that corresponds to the duration when the merger can be selected as a close pair.

Typically, the merger rate is defined in the following form

$$R_M(z, M_*) = \frac{f_{\text{merg}}(z, M_*)}{\tau_{\text{merg}}(z)} = \frac{f_P(z, M_*) \times C_{\text{merg}}}{\tau_{\text{merg}}(z)}, \quad (31)$$

where $f_{\text{merg}}$ is the merger fraction, i.e. the fraction of galaxies within a sample that will actually merge, and $\tau_{\text{merg}}$ is the merger time-scale (dynamical time) that is redshift dependent. Here, we make a distinction between the pair fraction and merger fraction, where the former is the fraction of observed close pairs, and the latter is the fraction of galaxies that will actually merge. The conversion factor between the quantities is denoted by $C_{\text{merg}}$, which is necessary as two nearby galaxies will only have some probability to eventually merge over some time-scale and might orbit each other for a significantly longer time (e.g. due to a long dynamical friction time-scale; Duncan et al. 2019). A typical value of $C_{\text{merg}} = 0.6$ is assumed, which is derived from simulations for all possible merging scenarios (Lotz et al. 2011; Conselice 2014; Mundy et al. 2017). However, this value is approximate and other sources place it in the range $C_{\text{merg}} \sim 0.4-1.0$.

In this work, we define the merger rate as

$$R_M(z, M_*) = \frac{f_P(z, M_*)}{\tau_P(z)} (\text{Gyr}^{-1}), \quad (32)$$

where $f_P$ is the pair fraction that we measured before, and $\tau_P$ (or $\tau_{\text{obs}}$) is the merger observability time-scale, which factors in the effects associated with $C_{\text{merg}}$ in equation (31). The observability time-scale is the duration a merging pair can be identified in a galaxy catalogue.

The derived merger rate is sensitive to the right choice of the observability time-scale (hereafter just merger time-scale or time-scale for shortness) according to equation (32). Various previous studies established that this time-scale is roughly constant with time at lower redshifts (typically $z < 3$) ranging between $\tau_P(z) \sim 0.3-1.0$ Gyr (Conselice 2006; Kitzbichler & White 2008; Lotz et al. 2008a, 2010). However, more recent studies using the IllustrisTNG simulations (Genel et al. 2014; Vogelsberger et al. 2014) at higher redshifts discovered that the merger time-scale does have a redshift dependence and evolves as a power law (Snyder et al. 2017; Duncan et al. 2019)

$$\tau_P(z) = 2.4 \times (1+z)^{-2} \text{ Gyr.} \quad (33)$$

A further refined version of the equation above by Conselice et al. (2022) based on TNG300-1 is also dependent on the mass ratio such that

$$\begin{align}
\tau_P(z) &= a \times (1+z)^b \quad \text{with} \quad a = -0.65 \pm 0.08 \nonumber \\
&\quad \times \mu + 2.06 \pm 0.01 \nonumber \\
&b = -1.60 \pm 0.01, \tag{34}
\end{align}$$

We note these quantities do not differ greatly from the typical halo dynamical time-scale, $0.1t_H(z)$, where $t_H(z)$ is Hubble-time (i.e. age

<figure><img src="image_12.png" /><figcaption>Figure 10. Comparison between the merger observability time-scale derived by Snyder et al. (2017) that is redshift-dependent, and the time-scale of Conselice et al. (2022) that is also mass-ratio dependent, and 10 per cent of the Hubble time. They agree well, especially at higher redshifts, with the lowest being the time-scale by Snyder et al. (2017), which we choose as our fiducial time-scale for deriving the merger rates.</figcaption></figure>

of the Universe) at redshift $z$, as can be seen in Fig. 10. For our fiducial merger rate analysis, we advocate for the time-scale derived by Snyder et al. (2017), similarly to Duncan et al. (2019).

## 5.2 Merger rate per galaxy

Using all definitions and results from above, we are now able to calculate the galaxy major merger rate, given by equation (32). This quantity essentially measures the average number of mergers per massive (primary) galaxy. We note here that this is equivalent to setting $C_{\text{merg}} = 1$ in equation (31) and just using the appropriate merger time-scale as in Duncan et al. (2019), Conselice et al. (2022), and Duan et al. (2024a).

As the merger time-scale decreases with redshift, the merger rate seen in Fig. 11 increases at higher redshifts even though pair fractions are declining or constant in that range. Our results show that merger rates at different stellar mass bins steeply increase with redshift at intermediate redshifts ($z \sim 3-6$) and then become constant at high redshifts ($z \sim 6-9$) in the range of $2-8 \text{ Gyr}^{-1}$ per galaxy. At the lowest stellar masses, it seems that there might be a turnover at $z \sim 6-7$. However, this might be affected by the uncertainties in the stellar mass correction weights, especially at high redshifts and low stellar masses, or potentially due to the lower limit (i.e. 5 kpc) adopted for the distance to consider a close pair.

We show measurements from this work and results from the literature in Fig. 12. Our results are plotted with uncertainty bands and are colour-coded by stellar mass, as previously shown in the case of pair fractions. Values of merger rates and associated uncertainties can be found in Table 3. We also plot results from the existing literature, as in the case of pair fractions. In this case, however, we use values for merger rates as quoted in the respective papers if they exist, where these are often calculated from pair fractions by a different choice of merger time-scale that is more appropriate for the respective data set. Not all works discuss merger rates, and in these cases, we use the same merger time-scale as we used to convert our measured pair fractions (see Section 5.1).

We fit a power law + exponential model to our measured merger rates and results from the literature for each stellar mass bin, as in
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Evolution of the galaxy major merger rate (top panel) and merger rate density (bottom panel) with redshift for the three stellar mass bins between $z \approx 3-9$. Merger rates steeply increase at intermediate redshifts ($z \approx 3-6$), after which they flatten and become constant at high redshifts ($z \approx 6-9$). At lower stellar masses, potentially there could be a turnover beyond $z \approx 7$. The same can be observed in the case of the redshift evolution of the merger rate densities, which show a stronger scaling with stellar mass.</figcaption></figure>

the case of pair fractions, having the form

$$ \mathcal{R}_M(z, M_*) = R_0 \times (1+z)^m \times e^{\tau(1+z)}, \quad (35) $$

where $R_0$, $m$, and $\tau$ are parameters for which we perform weighted non-linear least squares fitting with uncertainties incorporated as weights. The fitted parameters and their uncertainties can be found in Table 4, and the fitted curves are plotted in Fig. 12 with the respective colour coding.

## 5.3 Merger rate density

Another useful quantity to measure the merger rate that is frequently used in the literature is the merger rate density or total merger rate. The merger rate density describes how many mergers are occurring per unit time per unit comoving volume at a given redshift. We define the comoving merger rate density as

$$ \Gamma_M(z, M_*) = \frac{f_P(z, M_*)n_c(z, M_*)}{\tau_P(z)}, \quad (36) $$

where $n_c(z, M_*)$ is the comoving number density of (primary) galaxies within a redshift and stellar mass bin. The comoving number density of galaxies is calculated by directly integrating the SMF (equation B1) for a given stellar mass and redshift bin. While the volume-averaged galaxy merger rate $\Gamma_M(z)$ is a useful measure, in

the following text, we will only use the number of mergers per massive galaxy, i.e. the merger rate $\mathcal{R}_M(z)$.

## 5.4 Mass accretion rates

We define the specific mass accretion rate (sMAR) as in Duncan et al. (2019) in the following way:

$$ \text{sMAR}(z) = \frac{\dot{M}}{M} = \mathcal{R}_M(z)\bar{\mu}, \quad (37) $$

where $\bar{\mu}$ is the median mass ratio (or merger ratio). This quantity is analogous to the sSFR and measures the rate of mass accumulated by mergers per unit mass. The merger ratio $\bar{\mu}$ can be computed in each stellar mass–redshift bin by finding the median of the cumulative distribution of pair fraction as a function of merger ratio (see Fig. E1 in Appendix E for an example). These median merger ratios vary with stellar mass and redshift, and their average is $\langle \bar{\mu} \rangle = 0.485$ (see Fig. E2 in Appendix E).

Another way of calculating the specific mass accretion rate is to simply take the mass accretion rate $\rho_{1/4}$ and divide it by the average stellar mass of primary galaxies (Duan et al. 2024a) as

$$ \text{sMAR}(z) = \frac{\dot{M}}{M} = \frac{\rho_{1/4}(z)}{\mathcal{M}_{*,1}(z)} = \mathcal{R}_M(z) \frac{\mathcal{M}_{*,2}(z)}{\mathcal{M}_{*,1}(z)}, \quad (38) $$

where the mass accretion rate (MAR; analogous to SFR), which is the average amount of mass added through mergers per galaxy per unit time, is defined as

$$ \rho_{1/4} = \mathcal{R}_M(z) \times \mathcal{M}_{*,2}(z), \quad (39) $$

where the subscript ‘1/4’ denotes that this quantity is for major mergers. The average stellar mass of the primary and secondary galaxies in each redshift bin, respectively are defined using the galaxy SMF $\phi(M_*, z)$ as

$$ \mathcal{M}_{*,1}(z) = \frac{\int_{M_*^{\min}}^{M_*^{\max}} \phi(M_*, z) M_* dM_*}{\int_{M_*^{\min}}^{M_*^{\max}} \phi(M_*, z) dM_*} \quad (40) $$

$$ \mathcal{M}_{*,2}(z) = \frac{\int_{\mu, \mathcal{M}_{*,1}}^{\mathcal{M}_{*,1}} \phi(M_*, z) \mathcal{M}_* d\mathcal{M}_*}{\int_{\mu, \mathcal{M}_{*,1}}^{\mathcal{M}_{*,1}} \phi(M_*, z) d\mathcal{M}_*} \quad (41) $$

Now the ratio of the average stellar mass of the secondary and primary galaxies will be essentially the same as $\langle \bar{\mu} \rangle$, and we find that $\langle \mathcal{M}_{*,2}/\mathcal{M}_{*,1} \rangle = 0.474$ (see Fig. E2). This average value is close to the median merger ratio found by the previous method, and although it has a stellar mass–redshift variation, we make the assumption that $\bar{\mu} = 0.48$ for the rest of the calculations.

Finally, we plot the specific mass accretion rate as a second y-axis on Fig. 12 since it is just a rescaled version of the merger rate by a factor of $\bar{\mu}$ that we calculated earlier. This allows us to directly compare the contribution to mass accretion by mergers via sMAR and star formation via sSFR on the same plot. We plot five different SFMS sSFRs for comparison from Speagle et al. (2014), Faisst et al. (2016), Popesso et al. (2023), McClymont et al. (2025a), and Simmonds et al. (in preparation). The model by Speagle et al. (2014) was empirically determined from a compilation of 25 studies from the literature, is valid in the $z \sim 0-6$ range, and is well described by a single time-dependent function. The study by Faisst et al. (2016) finds a SFMS sSFR with a stronger power-law evolution at $z < 2.2$
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution of the galaxy major merger rate with redshift for the three stellar mass bins between $z \approx 3-9$. Merger rates increase at intermediate redshifts ($z \approx 3-6$), after which they flatten and become constant at high redshifts ($z \approx 6-9$). By fitting a power law + exponential curve to our measured merger rates and the ones available from the literature calculated at different stellar mass bins (Conselice et al. 2003, 2008, 2022; de Ravel et al. 2009; Casteels et al. 2014; Mundy et al. 2017; Ventou et al. 2017, 2019; Duncan et al. 2019; Dalmasso et al. 2024; Duan et al. 2024a), it becomes evident that the evolutionary trends show a weak mass dependence. Overlaid is the specific SFR of the SFMS at $M_* = 10^9 M_\odot$ for different prescriptions from the literature (Speagle et al. 2014; Faisst et al. 2016; Popesso et al. 2023; McClymont et al. 2025a; Simmonds et al. in preparation) which can be directly compared to the specific mass accretion rate (sMAR) as plotted on the right-hand side axis (essentially a rescaling of the merger rate; see Section 5.4). The coloured bands show the range of the sSFRs for the stellar mass range considered in this study, where the hatched sections are extrapolations of each respective relation. The sSFR is comparable to or larger than the sMAR in all cases, indicating that in situ star formation dominates over stellar mass accretion via major mergers.</figcaption></figure>

(extrapolation for $z < 0.5$) and a weaker power-law evolution at $z > 2.2$ (extrapolation for $z > 6$) with redshift, using data from the literature and spectroscopic data from the COSMOS field (Scoville et al. 2007) with median stellar mass of $\sim 10^{9.8} M_\odot$. Popesso et al. (2023) use a large sample of SFGs at $0 \le z \le 6$ with stellar masses $10^{8.5} - 10^{11.5} M_\odot$ and find a flattening sSFR above $z \sim 3$. Since all these empirical SFMS prescriptions are extrapolations beyond $z \sim 6$, and they show flattening sSFRs at higher $z$ that could be affected by selection effects, we opt for the model of McClymont et al. (2025a). Their sSFR is derived from the THESAN-ZOOM simulation (Kannan et al. 2025) and increases at higher redshifts, being valid for the $z > 3$ range. This SFMS choice is further supported by a good agreement with an observational counterpart based on the JADES data set at $z \sim 3-9$ by Simmonds et al. (in preparation).

# 6 DISCUSSION

In this section, we discuss the potential origins and implications of our measurements of galaxy close-pair fractions and major merger rates in the broader context of galaxy formation and evolution, focusing on the mass assembly of galaxies through cosmic time. We give possible explanations for the evolutionary trends found and compare them to predictions and results from large-scale cosmological simulations. We compare the importance of galaxy mergers and star formation in terms of mass assembly, particularly focusing on the difference between *in situ* star formation versus *ex situ* stellar mass accretion.

## 6.1 Qualitative reasoning behind the evolutionary trends found

In this work, we utilize the deepest available data from JWST NIRCam photometry through JADES observations, combined with all available spectroscopic data, to measure galaxy close-pair fractions

and derive major merger rates in the poorly understood redshift range $z \sim 3-9$. We successfully incorporate the full posterior distributions of the estimated photometric redshifts and propagate them throughout our analysis to find robust results with reliable uncertainties.

The redshift evolution of the galaxy close-pair fraction shows a significant dependence on the stellar mass bin considered for selecting the primary galaxies. In the case of the $\log_{10}(M_*/M_\odot) = [8.0, 8.5]$ stellar mass bin, there is an increase in pair fractions at redshifts $z \sim 3-5$, peaking at $z \sim 5-6$, that is followed by a turnover (see Fig. 8). For $\log_{10}(M_*/M_\odot) = [8.5, 9.0]$ a similar, but weaker peak can be observed at $z \sim 4-5$, followed by a decreasing trend. For the third stellar mass bin at $\log_{10}(M_*/M_\odot) = [9, 10]$, a flat and slightly decreasing trend can be seen, and if there is a peak as in the previous cases, it is at $z \sim 3-4$ or at lower redshift. Therefore, we conclude that the redshift where the pair fraction turns over increases with decreasing stellar mass. The same behaviour was observed in the theoretical study of Huško et al. (2022), which is based on the Planck Millennium cosmological simulation. This can be explained by the stellar mass range of interest passing into the exponentially decreasing regime of the galaxy SMF at increasing redshifts. Based on recent studies (e.g. Navarro-Carrera et al. 2024; Weibel et al. 2024; Harvey et al. 2025) the characteristic stellar mass ($M^*$), where the SMF transitions from a power-law decline for massive galaxies to an exponential cut-off, has a weak dependence on redshift. This transition happens at different redshifts for the different stellar mass bins considered (the characteristic stellar mass being in the range $M^* \approx 10^{10-11} M_\odot$), and we simply expect the pair fractions to trace this sudden absence of massive galaxies at higher redshifts.

Although the $\sim 200$ arcmin² area studied in this work is small compared to the entire sky, we conclude that we are not significantly affected by cosmic variance since results from both GOODS-South
and GOODS-North agree well within the uncertainties (see Fig. 8). Based on our results, there is more variation between the Deep and Medium regions within each field. This suggests that there is more dependence on the depth of the regions considered and, subsequently, on stellar mass completeness.

We also note here that in the literature, two closely related quantities are reported in general, namely the close-pair fraction and the merger fraction, highlighted by Huško et al. (2022). While either of these can be converted to a merger rate by assuming an appropriate time-scale, there is a significant difference between them. The close-pair fraction is measured from dynamical and geometrical selection criteria, that is, close separation projected on the sky and proximity in redshift or line-of-sight velocity. On the other hand, merger fractions are determined by selecting close pairs by morphological parameters that show signs of interactions or an early merging phase. This latter quantity attempts to measure the true fraction of galaxies that will actually merge, whereas the former includes random projections as well, which will not merge. At high redshifts ($z \gtrsim 6$), the close-pair method is more reliable for large samples, as it is currently difficult to select mergers based on morphological parameters (although it has been attempted by e.g. Dalmasso et al. 2024).

We observe similar trends in the pair fractions and merger rates, although, for the latter, the redshift evolution is weaker. Our results suggest a flattening in merger rates beyond $z \sim 6$, which is also highlighted in the trends of the fitted power law + exponential curves for each stellar mass bin. Focusing only on the results from this work, one can also see a slight turnover at the highest redshift bins for the different stellar mass ranges. This can be again explained by entering the exponentially decreasing regime of the SMF at high redshifts, as in the case of close-pair fractions (Huško et al. 2022).

We also compare our results with recent studies of pair fractions and merger rates. Duan et al. (2024a) use a similar probabilistic approach to this work (although slightly different selection criteria, e.g. separation limits $r = 20–50$ kpc), using data from eight JWST Cycle-1 fields, adding up to an area slightly smaller (~190 armin$^2$) than in this work. They use Bagpipes (Carnall et al. 2018) for SED-fitting to obtain stellar masses adopting a log-normal SFH, compared to Prospector with a non-parametric continuity SFH used in this work. They look at the same overall stellar mass range of $log_{10}(M_*/M_\odot) = [8, 10]$, but they do not carry out a full analysis on smaller stellar mass bins within this mass range as compared to this current work. However, they find lower pair fractions in subsample selected at higher stellar masses ($log_{10}(M_*/M_\odot) = 9–10$) that declines beyond $z \sim 5$ compared to a subsample with lower stellar masses ($log_{10}(M_*/M_\odot) = 8–9$), which is qualitatively in agreement with our results. Overall, their data points match our results well, falling within our uncertainty ranges, as seen in Figs 9 and 12. Furthermore, Dalmasso et al. (2024) also investigate close-pair fractions at high redshifts in the same stellar mass range but using morphological parameters to identify mergers. They find a flat evolution with redshift, and we convert their pair fractions using the time-scale chosen for our analysis to calculate merger rates. Similarly, the merger rates converted from their pair fractions fall within our reported uncertainty ranges.

## 6.2 SFR versus merger rate

To further assess the importance of major mergers in a quantitative way, we directly compare them to the other main channel of mass growth-star formation. One direct comparison is based on calculating the mass accreted by mergers and by star formation.

To calculate the cumulative mass accretion by mergers, we simply integrate the mass accretion rate (MAR; equation 39) for the desired cosmic time or redshift range as

$$ \delta M(z) = \int_{z_0}^{z} \rho_{1/4}(z) dz, \qquad (42) $$

where $z_0$ is chosen at $z = 9$ (which is the upper redshift limit of this study) and $z$ is at each subsequent point, and $\rho_{1/4}(z)$ is the mass accretion rate defined previously. We calculate this cumulative mass accretion from mergers, where for simplicity, we assume an overall fit for the merger rate $R_M(z)$ (given in Table 4) for the calculation of $\rho_{1/4}(z)$.

To directly compare the mass accreted from major mergers over time, we calculated a similar quantity for the mass accreted by star formation. Here, we use the SFMS fit of McClymont et al. (2025a) to derive the sSFR at $log_{10}(M_*/M_\odot) = 9$, the central mass in our total stellar mass range of interest. This SFMS sSFR is derived from the THESAN-ZOOM simulation (Kannan et al. 2025) that shows a rising sSFR at high-$z$. We note here that this simulation-based SFMS is in good agreement with its observational counterpart, which is calculated using the JADES data set by Simmonds et al. (in preparation; see Fig. 12). In the next section, we will discuss the motivation and effect of this choice of sSFR compared to other prescriptions in more detail.

We integrate this SFR the same way as we did for the MAR of mergers to obtain the cumulative mass accretion by star formation. The ratio of these two quantities gives us an approximate way to measure the ex situ mass fraction (only due to major mergers) in an average galaxy lying on the SFMS. The resulting ex situ mass fraction ($f_{ex-situ}$) has a weak evolution in the redshift range of $z \approx 3–9$ between 1–8 per cent, with a median value of $\sim 5$ per cent (see the purple curve on Fig. 13). Although approximate, this is one direct calculation to quantitatively assess the importance of major mergers from our inferred quantities, which shows that in situ star formation is the dominant channel for the mass growth of galaxies compared to major mergers.

We note here that these results are in contrast with a similar calculation performed by Duan et al. (2024a), where they find that the contribution of major mergers to galaxy stellar mass is 71 (±25) per cent (equivalent to the contribution from star formation). As their merger rate values agree well with our results, this discrepancy might arise because they use sSFRs estimated by Bagpipes instead of using a SFMS prescription. To further look into this disagreement, we looked at the median sSFR values of the entire initial sample at 10, 50, and 100 Myr measured by Prospector and found a good agreement with the McClymont et al. (2025a) and Simmonds et al. (in preparation) SFMS sSFR trends. Therefore, we conclude that our ex situ fraction estimates are consistent between SED-derived and simulation-based measurements.

## 6.3 In situ versus ex situ star formation

To look in more detail at the difference between the in situ versus ex situ mass fraction of galaxies due to star formation and major mergers from another perspective, we build a simple model for galaxy mass growth. We would like to note that our numerical model only focuses on major mergers, and we discuss the potential effect of minor mergers at the end of this section. Our model includes both gas-poor mergers (also referred to as dry mergers in the literature) and gas-rich mergers (also known as wet mergers). However, in the latter case, it only considers the stellar component, i.e. it ignores the
<figure><img src="image_15.png" /><figcaption>**Figure 13.** The redshift evolution of the cumulative *ex situ* stellar mass fraction due to major mergers from our numerical model that simulates the mass growth of galaxies using the merger rate from our measurements (valid for the $10^8-10^{10}$ M$_\odot$ stellar mass range) and assuming an SFR which always lies on the SFMS of McClymont et al. (2025a). We plot the median and the corresponding 16$^{th}$ and 84$^{th}$ percentiles of the distribution of the results from our numerical model of 50 000 galaxies. The median reaches only up to about ~13 per cent of *ex situ* fraction even at low redshifts. In the inset figure, we plot the mass growth history of three example galaxies from our simulation that have undergone several major mergers throughout their evolution, reaching final stellar masses indicated under each respective curve. We also plot the $f_{ex-situ}$ found by weighting it by the SMF (purple curve), and the $f_{ex-situ}$ estimated by taking the ratio between the stellar mass accreted and by mergers by star formation (grey dotted line; Section 6.2). All three calculations give similar results ranging between 3–13 per cent, placing in *situ* star formation as the dominant channel for the mass growth of galaxies compared to major mergers.</figcaption></figure>

gas that would be accreted and fuel star formation (as well as the triggering of star formation itself).

In our numerical model, we start with a sample of 50 000 galaxies with initial stellar masses randomly sampled from a linear distribution of masses with slope $\alpha = -1.5$ between $10^{-1}$ M$_\odot$ and $10^6$ M$_\odot$ at $z = 16$. We then look at the stellar mass growth evolution of these galaxies over 10 000 time-steps until $z = 0.5$. Our choices for these values are reasonable but somewhat arbitrary; nevertheless, they do not significantly affect the outputs we are interested in since these initial conditions are quickly erased by the nearly exponential growth of galaxies at these early times. During the evolution, at each time-step, we increase the mass of each galaxy in our initial sample by the amount of mass assembled by star formation, where we use the SFMS from McClymont et al. (2025a) (see Fig. 12). We discuss the effects of different choices of sSFRs at the end of this section. We also require that star formation is shut down when the galaxy reaches the quenching threshold mass of $M_{\text{quench}} = 10^{10.8}$ M$_\odot$, which is the location of the knee of the SMF. We also introduce a burstiness parameter for the star formation, that is, $10^b$ with b being a random number drawn from a Gaussian distribution centred at 0 with $\sigma = 0.3$. We multiply the SFR with this burstiness parameter at each step, and it is independent of the previous step. To first order, this takes into account the dispersion around the SFMS (Tacchella, Forbes & Caplar 2020).

The probability that a merger takes place is given by the merger rate at each time-step. For the merger rate, we use a single mass-independent fit from this work for simplicity, that is parametrized in equation (35) with fitted parameters given in Table 4 (see also

Fig. 12). We assume a Poisson distribution for a merger occurring, as these events are independent and relatively rare over a short time period. If there is a merger in a given time-step, we calculate the mass accreted through the merger by multiplying the merger rate $\mathcal{R}_M$ at that time (redshift) by the time-step and a fraction of the main galaxy’s mass. This fraction is the mass ratio that is a random number in the range $\mu = [0.25, 1]$ for a major merger for the purpose of this simulation. In the subsequent time-step, there is an independent draw for a merger to take place or not (for examples of galaxy mass growth histories, see the inset in Fig. 13).

This model allows us to calculate the *ex situ* mass fraction of galaxies by using empirical values such as merger rates and SFMS SFRs, as well as some simple assumptions. We plot in Fig. 13 the cumulative *ex situ* mass fraction, $f_{ex-situ}$, that is, measuring the redshift-evolution of the mass accreted from mergers compared to the total mass of the galaxy. We plot the median of $f_{ex-situ}$ and the 16$^{th}$ and 84$^{th}$ percentiles of the distribution of the galaxy sample from our simulation. As we might have expected based on Fig. 12 by comparing the merger rate to the sSFR, the *ex situ* mass fraction is slightly lower at higher redshifts (~3 per cent at $z = 9$) and slowly increases with cosmic time (to ~13 per cent at $z = 3$). In other words, based on this calculation, major mergers contributed less to the mass build-up of galaxies at earlier times and subsequently, their relative contribution increased (see the blue curve in Fig. 13). However, even at low-z, the median of $f_{ex-situ}$ reaches about 13 per cent, and the 84$^{th}$ percentile goes slightly below ~40 per cent. This is in agreement with our previous estimate of the *ex situ* mass fraction in Section 6.2 and places major mergers as a significant but not primary contributor to the mass assembly of galaxies.

In a third calculation to estimate $f_{ex-situ}$, we integrate the SMF (Appendix B) between $M_{\text{min}} = 10^8$ M$_\odot$ and $M_{\text{max}} = 10^{12}$ M$_\odot$ to obtain the cosmic stellar mass density (cSMD), and perform the same integration for a SMF convolved with the median $f_{ex-situ}$ calculated by the numerical model. This latter calculation should, in principle, give the contribution of major mergers to the cSMD. Finally, we take the ratio between the two quantities at redshifts in the range $z = 3-9$ to obtain a SMF weighted $f_{ex-situ}$ (see the grey dotted curve in Fig. 13), which shows a similar redshift evolution as the previous fraction at slightly lower values. All three calculations give consistent values of $f_{ex-situ}$, ranging between 3–13 per cent, and showing a weak redshift evolution.

We also look at the *ex situ* mass fractions as a function of the stellar mass of the galaxy. We plot in Fig. 14 the median $f_{ex-situ}$ values of the distribution of the galaxies from our simulation at different redshifts. Notably, we observe a scaling of *ex situ* mass fractions with redshift, where $f_{ex-situ}$ is higher for lower redshifts, as expected from Fig. 13. There is a sudden steep increase in $f_{ex-situ}$ above a mass of $10^{10.8}$ M$_\odot$, which is not surprising, as this is our chosen threshold mass above which a galaxy quenches its star formation, and the only remaining way to accrete more mass is through major mergers in this model. This is in agreement with predictions from Behroozi, Wechsler & Conroy (2013) and Behroozi & Silk (2015) (c.f. fig. 2 in Behroozi & Silk 2015). These results are also in agreement with the findings of Rodriguez-Gomez et al. (2016), Tacchella et al. (2019), and Huško et al. (2022) based on the Illustris, IllustrisTNG, and GALFROM simulations at $z = 0$. We also observe the same scaling with redshift as the one reported in Tacchella et al. (2019). Furthermore, these results seem to be consistent with observations of quiescent galaxies at lower redshifts at $z \sim 2$ (Ji & Giavalisco 2022, 2024; Ji et al. 2024). Since these galaxies are already quenched by the time of observation, their star formation must have stopped in earlier times. Based on detailed analysis of morphology and kinematics of these galaxies,
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Ex *situ* mass fraction $f_{\text{ex-situ}}$ from our numerical model at different stellar masses and redshifts. We observe that the $f_{\text{ex-situ}}$ stays approximately constant in the 3–13 per cent range for stellar masses up to $\sim 10^{10.8} M_{\odot}$, after which it quickly increases due to our mass threshold for quenching star formation. After a galaxy completely quenches its star formation, it can only grow its mass by mergers in this model, which leads to the up-tick in $f_{\text{ex-situ}}$ above the quenching threshold. We can also observe a dependence on redshift, where $f_{\text{ex-situ}}$ values are lower for higher redshifts. For comparison, we plot the *ex situ* mass fraction trends predicted by the Illustris (Rodriguez-Gomez et al. 2016), IllustrisTNG (Tacchella et al. 2019), and GALFORM (Huško, Lacey & Baugh 2022) simulations at $z=0$, and find that they are in good agreement with our results. Ex *situ* mass fractions measured from observations of local galaxies at $z \approx 0$ (Seigar, Graham & Jerjen 2007; Spavone et al. 2017, 2018, 2020; Cattapan et al. 2019; Iodice et al. 2020) also seem to match well our model predictions.</figcaption></figure>

they conclude that quiescent galaxies with $\log(M_*/M_\odot) > 11$ should have a significant *ex situ* mass contribution (e.g. mergers) to alter their structure after quenching, while the *ex situ* contribution is much smaller for lower mass $\log(M_*/M_\odot) < 11$ quiescent galaxies.

We note here that our results for the *ex situ* stellar mass fractions depend on the choice of the SFMS prescription. We use the SFMS$_{10}$ fit of McClymont et al. (2025a), which is derived from the THESAN-ZOOM simulation (Kannan et al. 2025). This SFMS is in good agreement with observational measurements based on JADES data (Simmonds et al. in preparation; see Fig. 12), as well as Prospector-derived sSFR values of our sample. Other SFMS fits from the literature, which are based directly on observational measurements of stellar mass and SFR (e.g. Speagle et al. 2014; Faisst et al. 2016; Popesso et al. 2023) agree well with each other in the $z < 3$ range (see Fig. 12), however, they diverge at the higher redshifts (beyond $z \approx 6$) and lower masses ($\log(M_*/M_\odot) < 9$) which we consider in this work. Additionally, the flattening of the SFMS evolution at high-$z$, which is seen in observational SFMS prescriptions, can be mostly explained by selection effects (McClymont et al. 2025a). We therefore opt for a simulation-derived prescription, which is not subject to these effects and instead shows a rising sSFR at high-$z$. Nevertheless, we run our numerical model using an sSFR from Speagle et al. (2014), Faisst et al. (2016), and Popesso et al. (2023), and find *ex situ* mass fractions in a similar range (below $\sim 16$ per cent) and mostly even lower due to the higher sSFR values in the different prescriptions. However, in these cases, the weak redshift evolution of $f_{\text{ex-situ}}$ is reversed, and it is in general higher at higher redshifts and subsequently slowly decreases with decreasing redshift (as it is expected based on Fig. 12). As an additional check, we also run the numerical model for the observational SFMS by Simmonds et al.

(in preparation) and find similar but slightly higher $f_{\text{ex-situ}}$ values (5–12 per cent) than our fiducial calculation. Therefore, while the numerical values and the weak redshift scaling slightly change based on the assumed SFMS, our overall conclusion about dominant *in situ* star formation compared to the effect of major mergers still does hold overall, independent of the choice of the SFMS prescription.

We further note here that our findings regarding the *in situ* versus *ex situ* mass growth of galaxies are based only on major mergers and do not take into account any contributions from minor mergers, i.e. mergers with mass ratios of $\le 1/4$. According to the hierarchical growth model of galaxies, minor mergers are expected to be more common than major mergers (e.g. Fakhouri, Ma & Boylan-Kolchin 2010). Previous models predicted that quiescent galaxies mainly grow via gas-poor minor mergers rather than major mergers (e.g. Bezanson et al. 2009). Studies by HST found that the size growth of massive quiescent galaxies is caused by a significant number of minor mergers (see e.g. Newman et al. 2012; Belli, Newman & Ellis 2015). A recent study using JWST by Suess et al. (2023) found a large population of low-mass companions (mostly having mass ratios $< 1/10$) around massive quiescent galaxies ($\log(M_*/M_\odot) > 10.5$). Their results suggest that in the redshift range $0.5 \le z \le 3$, minor mergers are potentially the main drivers of quiescent galaxy size growth and can also introduce radially decreasing colour gradients. However, the study of minor mergers is beyond the scope of this paper and will require a further independent study to constrain their rate of occurrence.

## 6.4 Comparison to theoretical models

In this section, we compare our results to findings from large-scale cosmological simulations and semi-analytical models. The existing results from the literature can be divided approximately into two main categories, namely, studies that claim that merger rates flatten or turn over at high redshifts depending on the mass selection (e.g. Huško et al. 2022); and sources that advocate for the monotonic increase of merger rates with redshift (e.g. Hopkins et al. 2010a; Rodriguez-Gomez et al. 2015; O'Leary et al. 2021).

Our results agree well qualitatively with the predictions of Huško et al. (2022) using the Planck Millennium cosmological simulation and the GALFORM semi-analytical model of galaxy formation (Baugh et al. 2019). Interestingly, most models show an increase in merger rates at all redshifts, unlike GALFORM, which predicts a plateau and a turnover depending on the stellar mass of interest. Although, Huško et al. (2022) calculated the merger rate at $\log_{10}(M_*/M_\odot) = [9, 10.3, 11]$, all these results show similar trends to our findings, where they rise and then stabilize at intermediate redshifts, and subsequently turn over (see fig. 6 from the cited work). They explain this by the exponential drop-off in the galaxy SMF. They find similar behaviour in the close-pair fractions (see fig. 7 in their work), where the turnover at higher masses moves toward lower masses with increasing redshift, reflecting the pattern seen in the merger rate. Again, this shift occurs because fewer galaxies populate the steeply declining tail of the SMF. This trend is more prominent in our observational measurements, as can be seen in Figs 7 and 9. This is further supported by not only the decreasing ratio of massive and satellite galaxies (i.e. pair fractions) but also the absolute number counts of primary galaxies declining at high redshifts, as expected.

Notably, Huško et al. (2022) highlight an interesting difference between physical pairs (with actual 3D separations matching their projected distances) and projected pairs (with significantly larger line-of-sight separations than their projected distances, i.e. the galaxies will not merge), which cannot be disentangled in
observational studies as opposed to simulations. Their results reveal that physical pair fractions, analogous to merger fractions, remain nearly constant with stellar mass, while projected pair fractions decline with increasing mass, consistent with the decreasing trend of the SMF. Notably, even at high stellar masses, about half of all pairs are due to projection (reducing by a similar amount the merger rate). For these massive systems, projected pairs are predominantly found within the same dark matter halo.

Other theoretical models predict qualitatively similar trends for both pair fractions and merger rates. The empirical model EMERGE (O'Leary et al. 2021) predicts a flattening in pair fractions and a weak turn over at high redshifts, including a scaling with the stellar mass selection, as found by this work. Similarly, the merger rates calculated from their mock observations can be best described by a power-law exponential evolution. Moreover, Endsley et al. (2020) use the UniverseMachine semi-empirical model to produce mock observations for JWST and find similar results as from GALFORM and EMERGE. Findings from Illustris by Snyder et al. (2017) also agree with these results, especially at the lower redshift regime.

On the contrary, some simulation-based predictions claim an ever-
rising merger rate with redshift, even for the high stellar mass regime
(e.g. $M_* > 10^{11} M_\odot$). It has been well established in the past that
the dark matter halo-halo merger rate rises with redshift (e.g. Genel
et al. 2008, 2010; Fakhouri et al. 2010). Therefore, assuming a direct
correlation between the halo-halo and galaxy-galaxy merger rates,
it would be expected that the latter also increases with redshifts.
Using the Illustris simulation, Rodriguez-Gomez et al. (2015) claims
that the major merger rate increases as a power law with redshift, in
agreement with the semi-empirical models of Stewart et al. (2009)
and Hopkins et al. (2010a), having a similar redshift dependence
as the halo merger rate. As discussed in Huško et al. (2022), these
discrepancies could arise for several reasons, such as differences in
definitions and methodology. The Stewart et al. (2009) and Hopkins
et al. (2010a) models, like EMERGE and GALFORM, assign merger
times when a halo becomes a subhalo but do not track subhalo
evolution within the primary halo. Additionally, their results rely
on abundance matching up to $z = 2$, making predictions at higher
redshift extrapolations.

Our results relating to $f_{\text{ex-situ}}$ are consistent with the findings
of Baker et al. (2025), where they look at *in situ* versus *ex situ*
mass build-up for quiescent galaxies in the FLAMINGO simulations
(Schaye et al. 2023) and conclude that on average high-$z$ quies-
cent galaxies have not undergone multiple major mergers ($\sim$0.6
major+minor mergers per galaxy per Gyr using our definition) and
that their mass growth was dominated by *in situ* star formation.

**6.5 Caveats and further improvements**

Here, we discuss the main uncertainties related to our analysis of
merger rates. The three main sources that could bias our calcula-
tions are incompleteness, photometric redshift uncertainties, and the
assumed merger time-scale.

We conduct a detailed assessment of the stellar mass completeness,
as presented in Section 2.5.2, and define our primary mass bins to
be above the completeness limit. We assign weights to account for
cases where the mass approaches the completeness limit or when
the secondary mass bin partially enters the incomplete regime, as
described in Section 3.4.1. These weights depend on the assumed
SMF (see Appendix B) and might fully account for the incomplete-
ness. Therefore, we might still be affected by incompleteness at
high redshifts. Another systematic source of uncertainty simply lies

within the nature of photometric redshifts estimated by EAZY. We
perform an odds quality parameter cut (see Section 2.2.3) to keep
the sources with well-defined photo-z's and also assign a weight
accounting for less constrained redshifts in the case of fainter objects
(see Section 3.4.2). We also compare the photometric redshifts to
available spectroscopic redshifts and find a good agreement (see
Section 2.3). However, considering all this, our sample could still
be affected by uncertainties due to projection effects and increased
pair probability functions due to the larger overlap between wider
photometric redshift posteriors. Furthermore, another selection effect
could be related to the distance limits adopted in the definition of
close-pairs, and in particular, the lower limit of $r_{\min} = 5$ kpc. This
could affect the inferred merger rates in the higher $z$ samples, for
which the merger time-scale is shorter.

Finally, the assumption of a merger observability time-scale
has a significant impact on the conversion of pair fractions to
merger rates, which should be seriously considered. As we discuss
in Section 5.1, the choice of time-scale is crucial for estimating
merger rates, and various parametrizations exist in the literature. We
choose the observability time-scale given by Snyder et al. (2017),
which is obtained from the Illustris simulation (Genel et al. 2014;
Vogelsberger et al. 2014) in order to reconcile the merger rates with
the pair counts directly measured from the simulation. It has a redshift
dependence proportional to ~$(1+z)^{-2}$, which leads to shorter time-
scales at higher redshifts, which in turn raises the declining tail of
pair fractions. Merger time-scales unfortunately remain impossible
to measure directly at high redshifts, and the current theoretical
estimates differ by factors of 2-3 (Conselice 2006; Kitzbichler &
White 2008; Lotz et al. 2008a, 2010; Hopkins et al. 2010b; Snyder
et al. 2017; Conselice et al. 2022; Huško et al. 2022), therefore by
choosing one particular time-scale we are inevitably biased to that
assumed theoretical model.

7 SUMMARY AND CONCLUSIONS

This study presents a comprehensive analysis of the major merger
history of galaxies in the redshift range $z \approx 3-9$, leveraging the
unprecedented depth of the JADES data set. By combining advanced
probabilistic methodologies for incorporating photometric redshift
posteriors and robust corrections for selection effects, we provide
novel insights into the role of major mergers in galaxy evolution
during the early Universe.

We apply a probabilistic method that takes into account the
full photometric redshift posterior distributions and available spec-
troscopic redshifts to select major merger close-pairs with mass
ratios > 1/4 within projected separations of $r_p = 5-30$ kpc for
stellar masses in the range $\log_{10}(M_*/M_\odot) = [8, 10]$. We perform
the analysis in sub-tiers of the GOODS-South and GOODS-North
samples from JADES. We also correct for potential biases and
selection effects when compiling our close-pair sample. Our main
results can be listed as follows:

I. Close-pair fraction evolution: our analysis reveals that the evolution of galaxy pair fractions with redshift follows a distinct pattern depending on the stellar mass range considered. In the lowest stellar mass bin (8.0 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) < 8.5), the pair fractions increase up to z ~ 5–6 and exhibit a turnover at higher redshifts. In the case of the two higher mass bins, there is a flattening and weak decline with increasing redshift. This turnover suggests a reduced abundance of close pairs at the earliest epochs, possibly due to fewer galaxies populating the steeply declining tail of the SMF. At higher stellar masses, the pair fractions show less variation with redshift,
pointing to a more stable merger history for massive galaxies. Our
results agree well with those of previous close-pair studies.

II. **Major merger rate evolution:** we convert the pair fractions to major merger rates by assuming a merger time-scale, and we observe an apparent flattening in the merger rate evolution beyond $z \sim 6$, stabilizing in the range of $2-8$ Gyr$^{-1}$ per galaxy, depending on the stellar mass range considered. This behaviour reflects the interplay between the evolution pair fractions and the assumed merger observability time-scale and the influence of the observational selection criteria. The impact of stellar mass on inferred merger rates is weak, nevertheless evident across the studied redshift range. Our inferred merger rates match well with previous studies, and an exponential-power law functional form can well describe their overall redshift evolution.

III. **Merger contribution to mass growth:** a comparison between the cumulative stellar mass accretion from major mergers and the mass assembled through star formation through a numerical model and, similarly, integrating the SMF convolved with the *ex situ* fraction consistently reveal that major mergers (in the mass range $\log_{10}(M_*/M_{\odot}) = 8-10$) contribute approximately 3–13 per cent to the total mass growth over $z \approx 3-9$. This contribution is substantial but secondary to *in situ* star formation, which remains the primary driver of stellar mass build-up during this epoch in comparison with major mergers. These findings are consistent with predictions from cosmological simulations, which similarly estimate major mergers as a significant but not dominant channel for mass assembly in the early Universe. These results provide a critical constraint on the merger-driven growth of galaxies in the high-redshift Universe.

IV. **Methodological advances:** the probabilistic framework further refined and applied in this study ensures robust results by incorporating full photometric redshift probability distributions and applying correction weights for mass incompleteness, photometric redshift quality, and survey boundaries. This methodology overcomes key limitations of previous works and paves the way for future studies of galaxy interactions using *JWST* data. Our approach also highlights the importance of using mass-selected samples for deriving meaningful close-pair fractions and merger rates, particularly at high redshifts and low stellar masses probed.

This work establishes that major mergers, while an important aspect of hierarchical galaxy formation, play a supplementary role in galaxy mass growth during the early Universe. The trends in pair fractions and merger rates support the view that mergers are a secondary process, acting in concert with star formation to shape the evolution of galaxies. Nevertheless, mergers could still drive star formation within galaxies, and hence, the physical properties of these mergers should be studied next. By constraining the merger history of galaxies at $z \approx 3-9$, this study provides a benchmark for theoretical models and lays the groundwork for future observational efforts to explore galaxy assembly at even higher redshifts with *JWST* and beyond.

## ACKNOWLEDGEMENTS

We thank the referee for the helpful comments that improved the manuscript. DP acknowledges support by the Huo Family Foundation through a P.C. Ho PhD Studentship. DP and ST acknowledge support by the Royal Society Research Grant G125142. CS, FDE, GCJ, and RM acknowledge support by the Science and Technology Facilities Council (STFC), by the ERC through Advanced Grant 695671 'QUENCH', and by the UKRI Frontier Research grant

RISEandFALL. RM also acknowledges funding from a research pro-
fessorship from the Royal Society. SA acknowledges support from
the JWST Mid-Infrared Instrument (MIRI) Science Team Lead, grant
80NSSC18K0555, from NASA Goddard Space Flight Center to the
University of Arizona. SA acknowledges grant PID2021-127718NB-
100 funded by the Spanish Ministry of Science and Innovation/State
Agency of Research (MICIN/AEI/ 10.13039/501100011033). WMB
acknowledges support by a research grant (VIL54489) from Villum
Fonden. AJB acknowledges funding from the 'FirstGalaxies' Ad-
vanced Grant from the European Research Council (ERC) under the
European Union's Horizon 2020 research and innovation programme
(Grant Agreement No. 789056). SC acknowledges support by the
European Union's HE ERC Starting Grant No. 101040227-WINGS.
DJE, ZJ, BDJ, MR, BER, and CNAW acknowledge support by
a JWST/NIRCam contract to the University of Arizona, NAS5-
02015. DJE is also supported as a Simons Investigator. BER also
acknowledges support from the JWST Program 3215. WM thanks
the Science and Technology Facilities Council (STFC) Center for
Doctoral Training (CDT) in Data Intensive Science at the University
of Cambridge (STFC grant number 2742968) for a PhD studentship.
HÜ acknowledges funding by the European Union (ERC APEX,
101164796). Views and opinions expressed are however those of the
authors only and do not necessarily reflect those of the European
Union or the European Research Council Executive Agency. Neither
the European Union nor the granting authority can be held respon-
sible for them. The research of CCW is supported by NOIRLab,
which is managed by the Association of Universities for Research in
Astronomy (AURA) under a cooperative agreement with the National
Science Foundation. JW gratefully acknowledges support from the
Cosmic Dawn Center through the DAWN Fellowship. The Cosmic
Dawn Center (DAWN) is funded by the Danish National Research
Foundation under grant No. 140.

This work is based [in part] on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with PIDs 1180, 1181, 1210, 1286, 1895, 1963, and 3215. The authors acknowledge the teams led by PIs Daniel Eisenstein (PID 3215), Christina Williams, Sandro Tacchella, Michael Maseda (JEMS; PID 1963), and Pascal Oesch (FRESCO; PID 1895) for developing their observing program with a zero-exclusive-access period.

## DATA AVAILABILITY

The data underlying this article will be shared at a reason-
able request by the corresponding author. Fully reduced NIRCam
images and NIRSpec spectra are publicly available on MAST
(https://archive.stsci.edu/hlsp/jades), with doi:10.17909/8tdj-8n28,
doi:10.17909/z2gw-mk31, and doi:10.17909/fsc4-dt61 (Eisenstein
et al. 2023b; Rieke et al. 2023b; Bunker et al. 2024; D'Eugenio et al.
2025).

## REFERENCES

Agazie G. et al., 2023, ApJ, 952, L37

Alberts S. et al., 2024, ApJ, 975, 85

Allen N. et al., 2024, preprint (arXiv:2410.16354)

Arribas S. et al., 2024, A&A, 688, A146

Astropy Collaboration, 2022, ApJ, 935, 167

Bacon R. et al., 2023, A&A, 670, A4
Baker W. M. et al., 2025, MNRAS, 539, 557

Barnes J. E., Hernquist L., 1996, ApJ, 471, 115

Baugh C. M. et al., 2019, MNRAS, 483, 4922

Beckwith S. V. W. et al., 2006, AJ, 132, 1729

Behroozi P. S., Silk J., 2015, ApJ, 799, 32

Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143

Bell E. F. et al., 2006, ApJ, 640, 241

Belli S., Newman A. B., Ellis R. S., 2015, ApJ, 799, 206

Benítez N., 2000, ApJ, 536, 571

Bertin E., Arnouts S., 1996, A&AS, 117, 393

Bezanson R., van Dokkum P. G., Tal T., Marchesini D., Kriek M., Franx M., Coppi P., 2009, ApJ, 697, 1290

Bickley R. W., Ellison S. L., Patton D. R., Bottrell C., Gwyn S., Hudson M. J., 2022, MNRAS, 514, 3294

Bickley R. W., Ellison S. L., Patton D. R., Wilkinson S., 2023, MNRAS, 519, 6149

Bluck A. F. L., Conselice C. J., Bouwens R. J., Daddi E., Dickinson M., Papovich C., Yan H., 2009, MNRAS, 394, L51

Böker T. et al., 2023, PASP, 135, 038001

Bournaud F. et al., 2011, ApJ, 730, 4

Bradley L. et al., 2024, astropy/photutils: 1.12.0. Zenodo, available at: https://zenodo.org/records/10967176

Brammer G., 2021, eazy-py. Github, available at: https://github.com/gbramer/eazy-py

Brammer G. B., van Dokkum P. G., Coppi P., 2008, ApJ, 686, 1503

Bunker A. J. et al., 2024, A&A, 690, A288

Byler N., Dalcanton J. J., Conroy C., Johnson B. D., 2017, ApJ, 840, 44

Carlberg R. G., 1990, ApJ, 359, L1

Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379

Carnall A. C. et al., 2024, MNRAS, 534, 325

Casteels K. R. V. et al., 2014, MNRAS, 445, 1157

Cattapan A. et al., 2019, ApJ, 874, 130

Chabrier G., 2003, PASP, 115, 763

Charlot S., Fall S. M., 2000, ApJ, 539, 718

Claeysens A., Adamo A., Richard J., Mahler G., Messa M., Dessauges-Zavadsky M., 2023, MNRAS, 520, 2180

Conroy C., 2013, ARA&A, 51, 393

Conroy C., Gunn J. E., 2010, ApJ, 712, 833

Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486

Conselice C. J., 2006, ApJ, 638, 686

Conselice C. J., 2014, ARA&A, 52, 291

Conselice C. J., Bershady M. A., Dickinson M., Papovich C., 2003, AJ, 126, 1183

Conselice C. J., Rajgor S., Myers R., 2008, MNRAS, 386, 909

Conselice C. J., Mundy C. J., Ferreira L., Duncan K., 2022, ApJ, 940, 168

Covelo-Paz A. et al., 2025, A&A, 694, A178

Croton D. J., 2013, Publ. Astron. Soc. Aust., 30, e052

D'Eugenio F. et al., 2025, ApJS, 277, 4

Dalmasso N. et al., 2024, MNRAS, 533, 4472

de Graaff A. et al., 2024, A&A, 684, A87

de Ravel L. et al., 2009, A&A, 498, 379

Decarli R. et al., 2019, ApJ, 882, 138

Decarli R. et al., 2024, A&A, 689, A219

Desmons A., Brough S., Martínez-Lombilla C., De Propris R., Holwerda B., López-Sánchez Á. R., 2023, MNRAS, 523, 4381

Donnan C. T. et al., 2024, MNRAS, 533, 3222

Duan Q. et al., 2024a, preprint (arXiv:2407.09472)

Duan Q. et al., 2024b, preprint (arXiv:2411.04944)

Duncan K. et al., 2019, ApJ, 876, 110

Efron B., 1979, Ann. Stat., 7, 1

Efron B., 1981, Biometrika, 68, 589

Eisenstein D. J. et al., 2023a, preprint (arXiv:2306.02465)

Eisenstein D. J. et al., 2023b, preprint (arXiv:2310.12340)

Ellison S. L., Patton D. R., Mendel J. T., Scudder J. M., 2011, MNRAS, 418, 2043

Ellison S. L., Mendel J. T., Scudder J. M., Patton D. R., Palmer M. J. D., 2013, MNRAS, 430, 3128

Ellison S. L. et al., 2022, MNRAS, 517, L92

Endsley R., Behroozi P., Stark D. P., Williams C. C., Robertson B. E., Rieke M., Gottlöber S., Yepes G., 2020, MNRAS, 493, 1178

Erb D. K., Pettini M., Shapley A. E., Steidel C. C., Law D. R., Reddy N. A., 2010, ApJ, 719, 1168

Faisst A. L. et al., 2016, ApJ, 821, 122

Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267

Ferreira L., Conselice C. J., Duncan K., Cheng T.-Y., Griffiths A., Whitney A., 2020, ApJ, 895, 115

Finkelstein S. L. et al., 2024, ApJ, 969, L2

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Fragkoudi F., Grand R. J. J., Pakmor R., Gómez F., Marinacci F., Springel V., 2025, MNRAS, 538, 1587

Fujimoto S. et al., 2023, ApJ, 949, L25

Gao F. et al., 2020, A&A, 637, A94

Gardner J. P. et al., 2006, Space Sci. Rev., 123, 485

Gardner J. P. et al., 2023, PASP, 135, 068001

Garduño L. E., Lara-López M. A., López-Cruz O., Hopkins A. M., Owers M. S., Pimbblet K. A., Holwerda B. W., 2021, MNRAS, 501, 2969

Genel S. et al., 2008, ApJ, 688, 789

Genel S., Bouché N., Naab T., Sternberg A., Genzel R., 2010, ApJ, 719, 229

Genel S. et al., 2014, MNRAS, 445, 175

Giavalisco M. et al., 2004, ApJ, 600, L93

Gupta A. et al., 2023, ApJ, 957, L35

Hainline K. N. et al., 2024a, ApJ, 964, 66

Hainline K. N. et al., 2024b, ApJ, 964, 71

Harvey T. et al., 2025, ApJ, 978, 89

Hashimoto T. et al., 2023, ApJ, 955, L2

Helton J. M. et al., 2024a, ApJ, 962, 124

Helton J. M. et al., 2024b, ApJ, 974, L41

Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356

Hopkins P. F. et al., 2010a, ApJ, 715, 202

Hopkins P. F. et al., 2010b, ApJ, 724, 915

Hsiao T. Y.-Y. et al., 2023, ApJ, 949, L34

Hsiao T. Y.-Y. et al., 2024, ApJ, 973, 8

Huško F., Lacey C. G., Baugh C. M., 2022, MNRAS, 509, 5918

Iodice E. et al., 2020, A&A, 635, A3

Jakobsen P. et al., 2022, A&A, 661, A80

Jesseit R., Naab T., Peletier R. F., Burkert A., 2007, MNRAS, 376, 997

Ji Z., Giavalisco M., 2022, ApJ, 935, 120

Ji Z., Giavalisco M., 2024, preprint (arXiv:2410.00099)

Ji Z. et al., 2024, preprint (arXiv:2401.00934)

Jin S. et al., 2023, A&A, 670, L11

Jogee S. et al., 2009, ApJ, 697, 1971

Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021, ApJS, 254, 22

Johnson B. et al., 2024, dfm/python-fsps: v0.4.7. Zenodo, available at: https://zenodo.org/records/12447779

Jones G. C. et al., 2024, preprint (arXiv:2405.12955)

Kado-Fong E. et al., 2018, ApJ, 866, 103

Kannan R. et al., 2025, preprint (arXiv:2502.20437)

Kartaltepe J. S. et al., 2007, ApJS, 172, 320

Kartaltepe J. S. et al., 2010, ApJ, 721, 98

Kartaltepe J. S. et al., 2015, ApJS, 221, L11

Kitzbichler M. G., White S. D. M., 2008, MNRAS, 391, 1489

Koekemoer A. M. et al., 2011, ApJS, 197, 36

Kriek M., Conroy C., 2013, ApJ, 775, L16

Lamperti I. et al., 2024, A&A, 691, A153

Lanz L. et al., 2013, ApJ, 768, 90

Le Fèvre O. et al., 2000, MNRAS, 311, 565

Leja J., Johnson B. D., Conroy C., van Dokkum P. G., Byler N., 2017, ApJ, 837, 170

Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, 3

Li W. et al., 2023, ApJ, 944, 168
Lin L. et al., 2004, ApJ, 617, L9

López-Sanjuan C. et al., 2012, A&A, 548, A7

López-Sanjuan C. et al., 2014, A&A, 564, A127

López-Sanjuan C. et al., 2015, A&A, 576, A53

Lotz J. M., Jonsson P., Cox T. J., Primack J. R., 2008a, MNRAS, 391, 1137

Lotz J. M. et al., 2008b, ApJ, 672, 177

Lotz J. M., Jonsson P., Cox T. J., Primack J. R., 2010, MNRAS, 404, 575

Lotz J. M., Jonsson P., Cox T. J., Croton D., Primack J. R., Somerville R. S., Stewart K., 2011, ApJ, 742, 103

Lower S., Narayanan D., Leja J., Johnson B. D., Conroy C., Davé R., 2020, ApJ, 904, 33

Man A. W. S., Zirm A. W., Toft S., 2016, ApJ, 830, 89

Mantha K. B. et al., 2018, MNRAS, 475, 1549

Mantha K. B. et al., 2019, MNRAS, 486, 2643

Marconcini C. et al., 2024, MNRAS, 533, 2488

Margalef-Bentabol B. et al., 2024, A&A, 687, A24

Marshall M. A. et al., 2023, A&A, 678, A191

Marshall M. A. et al., 2024, preprint (arXiv:2410.11035)

McClymont W. et al., 2025a, preprint (arXiv:2503.00106)

McClymont W. et al., 2025b, preprint (arXiv:2503.04894)

Meyer R. A. et al., 2024, MNRAS, 535, 1067

Mihos J. C., Hernquist L., 1994, ApJ, 431, L9

Miller T. B. et al., 2024, preprint (arXiv:2412.06957)

Molino A. et al., 2014, MNRAS, 441, 2891

Montenegro-Taborda D., Rodriguez-Gomez V., Pillepich A., Avila-Reese V., Sales L. V., Rodríguez-Puebla A., Hernquist L., 2023, MNRAS, 521, 800

Moreno J., Torrey P., Ellison S. L., Patton D. R., Bluck A. F. L., Bansal G., Hernquist L., 2015, MNRAS, 448, 1107

Moreno J. et al., 2019, MNRAS, 485, 1320

Morishita T. et al., 2024, ApJ, 963, 9

Mundy C. J., Conselice C. J., Duncan K. J., Almaini O., Häußler B., Hartley W. G., 2017, MNRAS, 470, 3507

Naab T., Khochfar S., Burkert A., 2006, ApJ, 636, L81

Navarro-Carrera R., Rinaldi P., Caputi K. I., Iani E., Kokorev V., van Mierlo S. E., 2024, ApJ, 961, 207

Newman A. B., Ellis R. S., Bundy K., Treu T., 2012, ApJ, 746, 162

O'Leary J. A., Moster B. P., Naab T., Somerville R. S., 2021, MNRAS, 501, 3215

Oesch P. A. et al., 2023, MNRAS, 525, 2864

Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

Patton D. R., Carlberg R. G., Marzke R. O., Pritchet C. J., da Costa L. N., Pellegrini P. S., 2000, ApJ, 536, 153

Patton D. R. et al., 2002, ApJ, 565, 208

Patton D. R., Ellison S. L., Simard L., McConnachie A. W., Mendel J. T., 2011, MNRAS, 412, 591

Patton D. R., Torrey P., Ellison S. L., Mendel J. T., Scudder J. M., 2013, MNRAS, 433, L59

Patton D. R. et al., 2020, MNRAS, 494, 4969

Pearson W. J. et al., 2019, A&A, 631, A51

Pearson W. J. et al., 2022, A&A, 661, A52

Perna M. et al., 2023, A&A, 679, A89

Perna M. et al., 2025, A&A, 696, A59

Planck Collaboration VI, 2020, A&A, 641, A6

Popesso P. et al., 2023, MNRAS, 519, 1526

Pozzetti L. et al., 2010, A&A, 523, A13

Reeves A. M. M., Hudson M. J., 2024, MNRAS, 527, 2037

Rieke M. J. et al., 2023a, PASP, 135, 028001

Rieke M. J. et al., 2023b, ApJS, 269, 16

Robertson B. E. et al., 2023, Nat. Astron., 7, 611

Rodríguez Del Pino B. et al., 2024, A&A, 684, A187

Rodríguez-Gomez V. et al., 2015, MNRAS, 449, 49

Rodríguez-Gomez V. et al., 2016, MNRAS, 458, 2371

Romano M. et al., 2021, A&A, 653, A111

Rose C. et al., 2023, ApJ, 942, 54

Satyapal S., Ellison S. L., McAlpine W., Hickox R. C., Patton D. R., Mendel J. T., 2014, MNRAS, 441, 1297

Saxena A. et al., 2023, A&A, 678, A68

Schaye J. et al., 2023, MNRAS, 526, 4978

Schechter P., 1976, ApJ, 203, 297

Scholtz J. et al., 2025, MNRAS, 539, 2463

Scoville N. et al., 2007, ApJS, 172, 1

Seigar M. S., Graham A. W., Jerjen H., 2007, MNRAS, 378, 1575

Serjeant S., Bakx T. J. L. C., 2023, Nat. Astron., 7, 1143

Shapiro K. L. et al., 2008, ApJ, 682, 231

Sharma R. S. et al., 2024, MNRAS, 527, 9461

Shibuya T., Miura N., Iwadate K., Fujimoto S., Harikane Y., Toba Y., Umayahara T., Ito Y., 2022, PASJ, 74, 73

Silk J., Rees M. J., 1998, A&A, 331, L1

Simmonds C. et al., 2024, MNRAS, 535, 2998

Snyder G. F., Lotz J. M., Rodriguez-Gomez V., Guimarães R. d. S., Torrey P., Hernquist L., 2017, MNRAS, 468, 207

Spavone M. et al., 2017, A&A, 603, A38

Spavone M. et al., 2018, ApJ, 864, 149

Spavone M. et al., 2020, A&A, 639, A14

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Stewart K. R., Bullock J. S., Barton E. J., Wechsler R. H., 2009, ApJ, 702, 1005

Suess K. A. et al., 2023, ApJ, 956, L42

Sun G., Faucher-Giguère C.-A., Hayward C. C., Shen X., 2023, MNRAS, 526, 2665

Tacchella S. et al., 2019, MNRAS, 487, 5416

Tacchella S., Forbes J. C., Caplar N., 2020, MNRAS, 497, 698

Tacchella S. et al., 2022, ApJ, 927, 170

Tacchella S. et al., 2023, MNRAS, 522, 6236

Thorp M. D., Ellison S. L., Pan H.-A., Lin L., Patton D. R., Bluck A. F. L., Walters D., Scudder J. M., 2022, MNRAS, 516, 1462

Toft S. et al., 2014, ApJ, 782, 68

Toomre A., Toomre J., 1972, ApJ, 178, 623

Torrey P., Cox T. J., Kewley L., Hernquist L., 2012, ApJ, 746, 108

Treu T. et al., 2023, ApJ, 942, L28

Urrutia T. et al., 2019, A&A, 624, A141

Ventou E. et al., 2017, A&A, 608, A9

Ventou E. et al., 2019, A&A, 631, A87

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Weibel A. et al., 2024, MNRAS, 533, 1808

Westcott L. et al., 2025, ApJ, 983, 121

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Whitler L., Stark D. P., Endsley R., Leja J., Charlot S., Chevallard J., 2023, MNRAS, 519, 5859

Williams R. J., Quadri R. F., Franx M., 2011, ApJ, 738, L25

Williams C. C. et al., 2018, ApJS, 236, 33

Williams C. C. et al., 2023, ApJS, 268, 64

Willott C. J. et al., 2024, ApJ, 966, 74

Witstok J. et al., 2024, A&A, 682, A40

Witten C. et al., 2024, Nat. Astron., 8, 384

Xu C. K., Zhao Y., Scoville N., Capak P., Drory N., Gao Y., 2012, ApJ, 747, 85

Yuan F.-T., Zheng Z.-Y., Jiang C., Zhu S., Lin R., Cheng C., 2024, ApJ, 975, 53

Zepf S. E., Koo D. C., 1989, ApJ, 337, 34

## APPENDIX A: CORRECTION WEIGHTS

In this section, we discuss in more detail the various correction weights (introduced in Section 3.4) applied to the PPF for finding major galaxy close-pairs.

The first weight applied is a correction for stellar mass incompleteness (see Section 3.4.1). Fig. A1 shows the strength of the mass completeness weights as they approach the mass completeness limit $M_{\star}^{\text{comp}}(z)$. Each subfield (i.e. GS and GN, deep and medium tiers) will have its own corresponding mass completeness weights, as these depend on $M_{\star}^{\text{comp}}(z)$, which in turn directly depends on the $5\sigma$ point-source limit of the field. These weights are calculated
<figure><img src="image_17.png" /><figcaption>**Figure A1.** The magnitude of the mass completeness weights on the stellar mass–redshift plane, which are statistical weights incorporated in the pair probability function to account for potentially missing galaxies due to reduced search range in the proximity of the mass completeness limit. As the top panel for $w_1^{\text{comp}}$ and the bottom panel for $w_2^{\text{comp}}$ show, the strength of these weights exponentially increases, approaching the mass completeness limit.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure A2.** Dependence of the OSR on the AB magnitude of the galaxies (calculated from fluxes measured in F444W) for the GS-Medium field. As the objects get fainter, the OSR significantly drops, meaning they have poorly constrained photometric redshifts. We fit a sigmoid curve to these data points to obtain a continuous OSR($m$) relation.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure A3.** Example of missing search area around a galaxy at $z_{\text{phot}} = 2.73$ in the primary sample in GS-Medium. The inner circle (red) corresponds to $r_{\text{min}}$ and the outer circle to $r_{\text{max}}$. We are interested in the available search area in the annulus between the two circles. In this particular case, there is a missing area (represented by white) to the left of the source due to the survey boundary and to the top right due to a masked diffraction spike of a nearby bright star. Hence, we compute $f_{\text{area}} = 0.66$ for this primary galaxy.</figcaption></figure>

from the ratios of expected number densities of galaxies using the SMF for a given mass bin by equations (16) and (17). In the case of the primary galaxy’s mass completeness weights ($w_1^{\text{comp}}$; top panel of Fig. A1) there are just three mass bins, and they are not heavily affected by incompleteness (as they were chosen to be above $M_*^{\text{comp}}(z)$). On the other hand, weights corresponding to secondary galaxies have to be calculated for each possible stellar mass ($w_2^{\text{comp}}$; bottom panel of Fig. A1), and usually significantly affect the calculations (as mass bins of secondary galaxies can potentially enter the incomplete regime at higher redshifts). As these weights exponentially increase in magnitude when approaching the mass completeness limit and also depend on the assumed SMF, we balance their contributions by clipping $w_1^{\text{comp}}$ at the [0.1, 1.0] and $w_2^{\text{comp}}$ at the [1, 10] range.

The second weight applied to the PPF accounts for the quality of photometric redshifts for galaxies by their brightness. The odds quality parameter $\mathcal{O}$ effectively measures how well defined is the primary peak of photometric redshifts (as discussed in Sections 2.2.3 and 3.4.2). The odds sampling rate (OSR) is the ratio between the number of objects above an odds parameter threshold value ($\mathcal{O} \ge 0.3$) and the total number of objects within a magnitude bin of width $\Delta m = 0.25$, which translates to the quality of the photometric redshifts versus the brightness of the source, or how peaked are the photo-$z$ posteriors. Fig. A2 shows the OSR curve for the GS-Medium field. The deeper we probe objects at higher magnitudes, the worse the photometric redshift quality becomes; hence, the OSR decreases. According to this, we are potentially missing multiple primary and secondary galaxies due to their poorly constrained photo-$z$’s. Therefore, we account for these potentially missing close-pair candidates by multiplying the pair probability functions of faint objects by $w_{\text{OSR}}$ to boost their number counts. All magnitudes are calculated from F444W fluxes, as this filter has the most extensive coverage in our fields, and all objects should be detected in this band, even at high $z$.

The third weight we apply is to account for potentially missing secondary galaxies due to incomplete search annuli around primary
galaxies (see Section 3.4.3). This could be caused by two main reasons: a primary galaxy being close to the survey boundary (in most cases, less than 50 per cent of the search area missing) or due to masked regions (see Fig. A3 for an example). As mentioned earlier, we mask out bright stars and strong diffraction spikes, as they can introduce artificial objects and, hence, many close-pairs that are, in fact, artefacts. These area weights have the weakest effect on our results. Nevertheless, we implement them for completeness. We calculate the area weights for each field separately before using them for the close-pair finding algorithm.

## APPENDIX B: GALAXY STELLAR MASS FUNCTION

Here, we present the galaxy stellar mass function (SMF) assumed for this work based on two recent high-redshift JWST studies by Navarro-Carrera et al. (2024) and Harvey et al. (2025).

The first work, by Navarro-Carrera et al. (2024), uses JWST data from the HUDF and UKIDSS Ultra Deep Survey field, complemented with ancillary data from CANDELS to measure the galaxy SMF for redshifts at $3.5 < z < 8.5$ probing down to stellar masses of $\simeq 10^8 M_{\odot}$. Harvey et al. (2025) extend measurements to the redshift range of $6.5 < z < 13.5$ and uses deep JWST NIRCam observations from the PEARLS and other publicly available surveys (CEERS, GLASS, JADES GOODS-South, NGDEEP, and SMACS0723) to compile a sample of high-$z$ galaxy candidates and measure the SMF.

We adopt a Schechter function (Schechter 1976) as in both cited works and combine their results for the three fit parameters. This function has a power-law shape with an exponential fall-off and is characterized by three parameters: $\alpha$, $M^*$, and $\Phi^*$. The SMF has the following logarithmic form:

$$ \Phi(\log M) = \ln(10)\Phi^*e^{-10^{\log M-\log M^*}}(10^{\log M-\log M^*})^{\alpha+1}d\log M, \quad (B1) $$

with $\alpha$ being the slope of the low-mass end, $M^*$ is the characteristic stellar mass at which the function turns from a power law to an exponential fall-off, and $\Phi^*$ controls the overall scaling of the SMF.

For these three parameters, we take the results from the two papers above and fit lines using least-square fitting, including measurement uncertainties (see Fig. B1). In the case of the characteristic stellar mass, we fix it at a value of $\log_{10}(M^*/M_{\odot}) = 10.7$, which is

motivated mainly by previous studies and the seemingly weak dependence on redshift. For $\alpha$, we fit a line but cap it to a constant value of $\alpha = -2.2$ at $z > 9.5$. For $\log_{10}(\Phi^*)$, we simply fit a line. We use these combined parameters in our SMF to cover intermediate and high redshifts with deep JWST data.

## APPENDIX C: RESULTS WITH EAZY-PY

In this section, we discuss a variation of our fiducial analysis of major galaxy close-pairs, using stellar masses computed with eazy-py.

Running eazy-py takes significantly less computational resources than Prospector; however, it results in less reliable stellar masses and other estimated physical parameters. Therefore, we obtain stellar masses for each photometric source, resulting in a significantly larger overall sample. For consistency with the fiducial analysis, we choose the same redshift limit, SNR, and odds selection criteria, resulting in an initial sample with a size comparable to the Prospector analysis. Fig. C1 shows the stellar mass distribution of our sample against redshift. We perform the same completeness analysis as in Section 2.5.2 by Pozzetti et al. (2010), and we observe a significantly larger scatter in the 90th percentile completeness limits in the redshift bins, especially at $4 < z < 5$.

We perform the same close-pair selection as described in Section 3, with the exception of using a data set with stellar masses estimated by eazy-py instead of Prospector. We present our results in Fig. C2, where we find an almost constant evolution on pair fractions with redshifts in all three stellar mass bins. There is a weak decline after $z \approx 5$, which is more pronounced in the higher stellar mass bins. It is interesting to compare this result with Fig. 7, where the most striking difference is between the trends found in the $8.0 \le \log_{10}(M_*/M_{\odot}) < 8.5$ mass bin. In our fiducial model, there is a pronounced increase, a strong peak at $z \approx 6$, followed by a decline at high-$z$, whereas here, we find an almost flat evolution. Although we conduct the same analysis in this case, we are more confident in the results from the fiducial analysis, as we better trust the stellar masses estimated by Prospector than the ones from the fast template-fitting code eazy-py. Nevertheless, this additional analysis shows the effect of the adopted SED-fitting codes that can significantly alter the results. Ultimately, large spectroscopic surveys at high redshifts could give precise and accurate close-pair measurements that are currently not yet available for a large enough sample.

<figure><img src="image_20.png" /><figcaption>**Figure B1.** Redshift evolution of the parameters of the galaxy SMF used in this work based on recent JWST studies by Navarro-Carrera et al. (2024) and Harvey et al. (2025). For comparison, we also plot results from Weibel et al. (2024), and we find a good agreement overall. We fix the characteristic stellar mass at $\log_{10}(M^*/M_{\odot}) = 10.7$, assuming redshift independence. We fit a line for $\alpha$ and cap it at the lowest allowed value of $\alpha = -2.2$, while we fit a line to $\log_{10}(\Phi^*)$ using least-squares fitting with uncertainties.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure C1.** Stellar mass distribution and completeness calculation for the total GOODS-South sample using masses calculated by *eazy-py*. A significant scatter can be observed in the 90th percentile completeness limits between the different redshift bins, which is caused by the lower accuracy of *eazy-py* compared to *Prospector*. While the overall sample size is larger compared to the *Prospector* fits, as each photometric source will have a stellar mass estimate by *eazy-py*, the initial sample after pre-selection on which we perform our analysis has a similar size.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure C2.** Evolution of the major merger close-pair fractions with redshift stellar masses estimated by *eazy-py*, instead of *Prospector* as in the fiducial analysis. We observe an almost constant and slightly declining trend (beyond $z \approx 5$) with increasing redshifts, where the decline is more pronounced at the higher stellar mass bins. There is a significant contrast between the results from the fiducial analysis and this one, mainly in the lowest stellar mass bin, that might be caused by the lower accuracy of *eazy-py*.</figcaption></figure>

## APPENDIX D: PAIR FRACTIONS AT REDSHIFT-DEPENDENT SEPARATIONS

In this section, we discuss another variation of the fiducial close-pair analysis, where we use a redshift-dependent separation criterion instead of a constant one.

The separation criterion for this analysis is based on the virial radius of the dark matter halo hosting each primary galaxy at each redshift. To estimate this radius, we use the stellar mass-halo mass relation from Behroozi et al. (2019) using the UniverseMachine code. From this parametrization, we obtain the peak halo mass $M_{halo}$

<figure><img src="image_23.png" /><figcaption>**Figure D1.** Redshift dependence of the virial radii of dark matter haloes corresponding to the primary galaxies that we adopt as the redshift-dependent close-pair selection criterion. The square markers represent the values at each stellar mass-redshift bin. The constant separation criterion of $r_{sep} = [5, 30]$ kpc from the fiducial close-pair analysis is highlighted as a light blue rectangle for better comparison with the redshift evolution of the virial radii.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure D2.** Close-pair fraction evolution using stellar masses from *Prospector*, but adopting a redshift-dependent separation criterion instead of a constant one, based on the virial radii of the dark matter haloes hosting the primary galaxies. We find that the pair fractions are smoothly declining for all stellar mass bins, in contrast to the results found in the fiducial analysis.</figcaption></figure>

corresponding to the stellar mass $M_*$ in the centre of each stellar mass bin (at 8.25, 8.75, and 9.5 log $M_\odot$, respectively). In the case of $\log_{10}(M_*/M_\odot) = 9.5$ the parametrization of Behroozi et al. (2019) does not go beyond $z = 8.6$ and we extrapolate the trend up to $z = 9$, which gives reasonable results (see Fig. D1). We can compute the virial radius corresponding to the peak halo masses using the equation

$$R_{\text{vir}} = \left( \frac{2GM_{\text{halo}}}{\Delta_c H^2(z)} \right)^{1/3}, \qquad (\text{D1})$$

where $H(z)$ is the Hubble parameter, $G$ is the gravitational constant, and we assume $\Delta_c = 200$ for the density contrast. Using this equation, we get redshift-dependent virial radii for the DM haloes hosting the average primary galaxy in each redshift bin (see Fig. D1).
We adopt these redshift-dependent virial radii as the new separation criteria that we determine for each stellar mass–redshift bin (as indicated by the square markers on Fig. D1). Here, we make the assumption that if a secondary galaxy is closer to the primary than the corresponding virial radius of the DM halo hosting it, then the system forms a close pair and will eventually merge. This is a simple assumption that is physically motivated.

The resulting pair fractions are plotted in Fig. D2, which show somewhat different behaviour than the fiducial analysis and the results from eazy-py. For this analysis, we are using stellar masses from Prospector as in the fiducial analysis. We find that the pair fraction declines with redshift in all stellar mass bins considered. The decline is smoother and goes to lower values than in the fiducial analysis (c.f. Fig. 7). In comparison, this is simply caused by the smaller search radius around high-z primary galaxies according to the relation from Fig. D1. However, this might be physically more motivated than looking at a constant separation criterion, as the dark matter haloes at high-z were significantly smaller than the ones at lower z, and therefore gravitational interaction and mergers took place at smaller physical scales.

This analysis shows that choosing a different separation limit significantly affects the resulting close-pair measurements. We also note here that most previous works adopted a constant separation limit, whereas a redshift-dependent one is physically more motivated.

## APPENDIX E: MEDIAN MASS RATIO

Here, we discuss in more detail the analysis of the merger ratio distribution of our sample and the calculation of the average merger ratio. This step is essential to convert the merger rate $R_M(z)$ to the specific mass accretion rate (sMAR) (see equation 38).

First, we calculate the median mass ratio (merger ratio) $\bar{\mu}$ for our close-pair sample. This can be done by computing the pair fraction $f_P(> \mu)$ in each redshift–stellar mass bin above increasing values of $\mu$ (see Fig. E1 for an example). We then calculate the median value of the mass ratio from this cumulative pair fraction distribution by extrapolation to find the $\mu$ value at which the distribution reaches the 50th percentile. We also fit a curve to the average $f_P(> \mu)$ values that is described by a parametrization (Duncan et al. 2019) of the form

$$f_P(> \mu) = A \times \left( \frac{1}{\mu} - 1 \right)^B \qquad (E1)$$

We then calculate the median mass ratios in each stellar mass–redshift bin and plot them on Fig. E2. Interestingly, we find that there is a considerable variation of $\bar{\mu}$ values at different redshifts and stellar masses. This is mainly due to having fewer close pairs in the higher mass bins that are difficult to resolve into smooth $f_P(> \mu)$ distributions to infer $\bar{\mu}$. We simply take the average of all these median mass ratios to compute $\langle \bar{\mu} \rangle = 0.485$.

We compare the previously computed average median mass ratio $\langle \bar{\mu} \rangle$ to a predicted mass ratio based on the SMF. This is done by calculating the ratio between the predicted average mass in the secondary stellar mass bins $M_{*,2}(z)$, and the primary stellar mass bins $M_{*,1}(z)$ defined in equation (41). We calculate these ratios for the three stellar mass bins and find that they are similar and have a weak dependence on redshift (as seen in Fig. E2). We then take the

<figure><img src="image_25.png" /><figcaption>**Figure E1.** Example distribution of close-pair fractions with increasing mass ratio cuts for the stellar mass–redshift bin indicated on the plot for the four tiers and their averages. This is a cumulative distribution as pair fractions with lower $\mu$ values contain close-pairs with higher mass ratios than that $\mu$ limit (e.g. $f_P(> \mu = 0.3)$ also contains close-pairs included in $f_P(> \mu = 0.4)$). We calculate the median mass ratio in this bin by calculating the 50th percentile of this distribution by extrapolation and also fit a curve described by equation (E1).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure E2.** Empirically determined median merger ratios from the $f_P(> \mu)$ cumulative distributions in each stellar mass–redshift bin. We find that there is a significant scatter in these values a higher redshifts, which is due to the low number counts of close-pairs that results in a poor resolution of $f_P$ in terms of merger ratio, subsequently leading to less reliable $\bar{\mu}$. We also plot the redshift evolution of the predicted merger ratio of $\langle M_{*,2}(z)/M_{*,1}(z) \rangle$ using the SMF. The averages of these two quantities are displayed on the figure and agree well with each other.</figcaption></figure>

average of these ratios and find that $\langle M_{*,2}(z)/M_{*,1}(z) \rangle = 0.474$. This is in good agreement with our empirically determined value from the $f_P(> \mu)$ distributions.

This paper has been typeset from a TeX/Latex file prepared by the author.