We present a full analysis of galaxy major merger pair fractions merger rates and mass accretion rates thus uncovering the role of mergers in galaxy formation at the earliest previously unexplored epoch of 4.5ltzlt11.5. We target galaxies with masses log10mathrmMmathrmModot  8.0  10.0 utilizing data from eight JWST Cycle1 fields CEERS JADES GOODSS NEPTDF NGDEEP GLASS ElGordo SMACS0723 MACS0416 covering an unmasked area of 189.36 mathrmarcmin2. We develop a new probabilistic paircounting methodology that integrates full photometric redshift posteriors and corrects for detection incompleteness to quantify close pairs with physical projected separations between 20 and 50 kpc. Our analysis reveals an increase in pair fractions up to z  8 reaching 0.211 pm 0.065 followed by a statistically flat evolution to z  11.5. We find that the galaxy merger rate increases from the local Universe up to z  6 and then stabilizes at a value of sim 6 Gyr1 up to z  11.5. We fit both a powerlaw and a powerlaw  exponential model to our pair fraction and merger rate redshift evolution finding that the latter model describes the trends more accurately particularly at z  8.0  11.5. In addition we measure that the average galaxy increases its stellar mass due to mergers by a factor of 2.77 pm 0.99 from redshift z  10.5 to z  5.0. Lastly we investigate the impact of mergers on galaxy stellar mass growth revealing that mergers contribute 71 pm 25 as much to galaxy stellar mass increases as star formation from gas. This indicates that mergers drive about half of galaxy assembly at high redshift.
['Astrophysics - Astrophysics of Galaxies']
Galaxy Mergers in the Epoch of Reionization I A JWST Study of Pair Fractions Merger Rates and Stellar Mass Accretion Rates at z  4.511.5
{'Galaxy Mergers in the Epoch of Reionization I: A JWST Study of Pair Fractions, Merger Rates, and Stellar Mass Accretion Rates at 𝑧 = 4 . 5 -11 . 5': "Qiao Duan, 1 Christopher J. Conselice, 1 ★ Qiong Li, 1 † Duncan Austin, 1 ‡ Thomas Harvey, 1 Nathan J. Adams, 1 Kenneth J. Duncan, 2 James Trussler, 1 Leonardo Ferreira, 3 Lewi Westcott, 1 Honor Harris, 1 Rogier A. Windhorst, 4 Benne W. Holwerda, 5 Thomas J. Broadhurst, 6 , 7 , 8 Dan Coe, 9 , 10 , 11 Seth H. Cohen, 4 Xiaojing Du, 12 Simon P. Driver, 13 Brenda Frye, 14 Norman A. Grogin, 9 Nimish P. Hathi, 9 Rolf A. Jansen, 4 Anton M. Koekemoer, 9 Madeline A. Marshall, 15 , 16 Mario Nonino, 17 Rafael Ortiz III, 4 Nor Pirzkal, 9 Aaron Robotham, 13 Russell E. Ryan, Jr 9 Jake Summers, 4 Jordan C. J. D'Silva, 13 , 16 Christopher N. A. Willmer, 14 Haojing Yan 18  \n- 1 Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road, Manchester, UK\n- 2 Institute for Astronomy, University of Edinburgh Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK\n- 3 Department of Physics & Astronomy, University of Victoria, Finnerty Road, Victoria, British Columbia, V8P 1A1, Canada\n- 4 School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287-1404\n- 5 University of Louisville, Department of Physics and Astronomy, 102 Natural Science Building, 40292 KY Louisville, USA\n- 6 Department of Physics, University of the Basque Country UPV/EHU, E-48080 Bilbao, Spain\n- 7 DIPC, Basque Country UPV/EHU, E-48080 San Sebastian, Spain\n- 8 Ikerbasque, Basque Foundation for Science, E-48011 Bilbao, Spain\n- 9 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA\n- 10 Association of Universities for Research in Astronomy (AURA) for the European Space Agency (ESA), STScI, Baltimore, MD 21218, USA\n- 11 Center for Astrophysical Sciences, Department of Physics and Astronomy, The Johns Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA\n- 12 Department of Statistical Science, Wake Forest University, Winston-Salem, NC, USA\n- 13 International Centre for Radio Astronomy Research (ICRAR) and the International Space Centre (ISC), The University of Western Australia, M468, 35 Stirling Highway, Crawley, WA 6009, Australia\n- 14 Department of Astronomy/Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ, 85721-0009, USA\n- 15 National Research Council of Canada, Herzberg Astronomy & Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada\n- 16 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia\n- 17 INAF-Osservatorio Astronomico di Trieste, Via Bazzoni 2, 34124 Trieste, Italy\n- 18 Department of Physics and Astronomy, University of Missouri, Columbia, MO 65211, USA  \n27 November 2024", 'ABSTRACT': 'We present a full analysis of galaxy major merger pair fractions, merger rates, and mass accretion rates, thus uncovering the role of mergers in galaxy formation at the earliest previously unexplored epoch of 4 . 5 < 𝑧 < 11 . 5. We target galaxies with masses log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0, utilizing data from eight JWST Cycle-1 fields (CEERS, JADES GOODS-S, NEP-TDF, NGDEEP, GLASS, El-Gordo, SMACS-0723, MACS-0416), covering an unmasked area of 189.36 arcmin 2 . We develop a new probabilistic pair-counting methodology that integrates full photometric redshift posteriors and corrects for detection incompleteness to quantify close pairs with physical projected separations between 20 and 50 kpc. Our analysis reveals an increase in pair fractions up to 𝑧 = 8, reaching 0 . 211 ± 0 . 065, followed by a statistically flat evolution to 𝑧 = 11 . 5. We find that the galaxy merger rate increases from the local Universe up to 𝑧 = 6 and then stabilizes at a value of ∼ 6 Gyr -1 up to 𝑧 = 11 . 5. We fit both a power-law and a power-law + exponential model to our pair fraction and merger rate redshift evolution, finding that the latter model describes the trends more accurately, particularly at 𝑧 = 8 . 0 -11 . 5. In addition, we measure that the average galaxy increases its stellar mass due to mergers by a factor of 2 . 77 ± 0 . 99 from redshift 𝑧 = 10 . 5 to 𝑧 = 5 . 0. Lastly, we investigate the impact of mergers on galaxy stellar mass growth, revealing that mergers contribute 71 ± 25% as much to galaxy stellar mass increases as star formation from gas. This indicates that mergers drive about half of galaxy assembly at high redshift.  \nKey words: galaxies: merger -galaxies: high redshift - galaxies: formation  \n- ★ E-mail: conselice@manchester.ac.uk\n- † E-mail: qiong.li@manchester.ac.uk', '1 INTRODUCTION': "The formation and evolution of galaxies are intricately linked to their merger history. Galaxies are believed to form from the mergers of cold dark matter halos (Blumenthal et al. 1984), suggesting that contemporary galactic structures evolved from the mergers of smaller predecessors, a concept supported by the theory of dark matter halo mergers (Jenkins et al. 1998; Maller et al. 2006). Beyond the simple merging of galaxies, these activities profoundly impact various galaxy properties, affecting phenomena such as size growth (Trujillo et al. 2007; Buitrago et al. 2008; Naab et al. 2009; Bluck et al. 2012), star formation (Mihos & Hernquist 1995; Scudder et al. 2012; Patton et al. 2013; Pearson et al. 2019; Patton et al. 2020; Garduño et al. 2021; Ellison et al. 2022; Shah et al. 2022; He et al. 2023; Wolter et al. 2023; Montenegro-Taborda et al. 2023), AGN activity (Satyapal et al. 2014; Barrows et al. 2017; Gao et al. 2020; Li et al. 2023; Sharma et al. 2024), morphological transformations (Naab et al. 2006; Bournaud et al. 2011), as well as possible neutrino emissions (Bouri et al. 2024).  \nRecent advancements have also led to discoveries such as dual and triple AGN systems at redshifts around 𝑧 = [ 3 , 4 ] (Perna et al. 2023), as well as evidence that mergers may catalyze the formation of extreme emission line galaxies in the early Universe (Gupta et al. 2023). Mergers also provide possible insights into magnetic fields and how they play a pivotal role in the dynamics of disk galaxy mergers through magnetohydrodynamic (MHD) simulations (Whittingham et al. 2023). Galaxy mergers thus drive the formation of not only the dark matter and stellar masses of galaxies, but also in the construction of central massive black holes and triggering of star formation. Mergers are an essential process in galaxy formation, yet we do not have very much observational information on their influence across cosmic time.  \nUnderstanding pair fractions and merger rates across the Universe's history not only elucidates the merger process itself but also provides insights into broader cosmic environments and construction at different epochs. Two main methods are employed to study merger phenomena: morphology-based selection (e.g., Conselice et al. 2003; Conselice et al. 2008; Lotz et al. 2008a,b; Hwang & Chang 2009; Lotz et al. 2010a,b; Whitney et al. 2021; Pearson et al. 2022; Ren et al. 2023; Rose et al. 2023; Desmons et al. 2023; Wilkinson et al. 2024; Dalmasso et al. 2024) and close pair analysis (e.g., Zepf & Koo 1989; Burkey et al. 1994; Carlberg et al. 1994; Woods et al. 1995; Yee & Ellingson 1995; Patton et al. 1997; Le Fèvre et al. 2000; Patton et al. 2000; de Ravel et al. 2009; López-Sanjuan et al. 2013, 2015; Patton et al. 2016; Man et al. 2016; Mundy et al. 2017; Snyder et al. 2017; Ventou et al. 2017; Mantha et al. 2018; Duncan et al. 2019; Conselice et al. 2022; Inoue et al. 2024). The morphological selection method employs criteria such as the CAS system (Schade et al. 1995; Abraham et al. 1996; Bershady et al. 2000; Conselice 2003; Conselice et al. 2003), the Gini coefficient and 𝑀 20 (Lotz et al. 2004), to identify merging pairs. Close pair methodology, on the other hand, utilizes redshifts and projection separations to select galaxies likely to merge in the foreseeable future (López-Sanjuan et al. 2013, 2015; Man et al. 2016; Mundy et al. 2017; Ventou et al. 2017; Snyder et al. 2017; Mantha et al. 2018; Duncan et al. 2019; Conselice et al. 2022).  \nThe measurement of pair fractions and merger rates up to 𝑧 < 3 . 5 has been extensively explored (e.g., Conselice 2003; Conselice et al. 2008; Lotz et al. 2011; Bluck et al. 2012; Ownsworth et al. 2014; López-Sanjuan et al. 2015; Man et al. 2016; Mundy et al. 2017; Mantha et al. 2018; Duncan et al. 2019; Conselice et al. 2022; Patton et al. 2024). However, studies analyzing these metrics at redshifts of  \n3 . 0 ≤ 𝑧 ≤ 6are less frequent (Ventou et al. 2017; Duncan et al. 2019), primarily due to the limitations of earlier observational data. Previous observations lacked the necessary sensitivity and wavelength coverage to accurately determine the redshifts of high𝑧 galaxy samples, constraining our understanding of galaxy evolution during these epochs.  \nThehigh-redshift Universe is increasingly unveiled by JWST, highlighted by numerous recent studies on early galaxy discoveries (e.g., Harikane et al. 2022; Atek et al. 2022; Castellano et al. 2022; Donnan et al. 2022; Naidu et al. 2022; Curtis-Lake et al. 2022; Adams et al. 2023; Maiolino et al. 2023; Austin et al. 2023a; Donnan et al. 2023; Finkelstein et al. 2023; Arrabal Haro et al. 2023; Bouwens et al. 2023; McLeod et al. 2023; Hainline et al. 2023a; Franco et al. 2023; Casey et al. 2023; Maiolino et al. 2024; Adams et al. 2024; Chakraborty et al. 2024). These studies reveal a significant population of distant galaxies at 𝑧 > 6, previously underrepresented in HST data, providing a new window to explore merger evolution at even high redshifts.  \nGiven the challenges of identifying mergers through morphology at high redshifts, where galaxies are more compact and may appear as point-like sources, this paper pioneers the study of highredshift merger evolution using the close-pair method. We utilize data from JWST's Near-InfraRed Camera (NIRCam; Rieke et al. 2005, 2023c) in eight deep Cycle-1 blank imaging fields (CEERS, JADES GOODS-S,NEP-TDF,NGDEEP,GLASS,El-Gordo,SMACS-0723, MACS-0416), covering unmasked area of 189 . 36 arcmin 2 . We expand upon a novel close-pair methodology, building on foundational works (López-Sanjuan et al. 2015; Mundy et al. 2017; Duncan et al. 2019; Conselice et al. 2022) to compute pair fractions, merger rates, and mass accretions from 𝑧 = 4 . 5 to 𝑧 = 11 . 5 for the first time. We then use these to derive the history of galaxy merging and how mergers contribute to the mass assembly of galaxies at this very early cosmic times.  \nThe structure of this paper is outlined as follows. In Section 2, we detail the dataset sourced from 8 JWST fields. Our photometric redshift and stellar mass determinations are outlined in Sections 3 and 3.2. The close pair methodology is presented in Section 4. Our main findings and analysis are presnted in Section 5. A summary of our conclusions is provided in Section 7. Throughout this work, we adhere to the Planck 2018 Cosmology (Planck Collaboration et al. 2020) with 𝐻 0 = 67 . 4 ± 0 . 5 km s -1 Mpc -1 , Ω M = 0 . 315 ± 0 . 007, and Ω Λ = 0 . 685 ± 0 . 007 to facilitate comparison with other observational studies. All magnitudes reported are consistent with the AB magnitude system (Oke 1974; Oke & Gunn 1983).", '2 OBSERVATIONS AND DATA REDUCTION': "The launch of the James Webb Space Telescope (JWST) in December 2021 (Rigby et al. 2023) marked the beginning of a new era in the exploration of the high-redshift Universe. We utilized a consistent approach for reducing and analyzing publicly available deep JWST NIRCamdatafromthePEARLS,CEERS,GLASS,JADESGOODSS, NGDEEP, and SMACS0723 surveys. Specifically in this paper, for the PEARLS program, we included data from the NEP-TDF, MACS-0416, and El-Gordo fields, leading to a combined total of 189.36 arcmin 2 of unmasked sky when integrated with the other five surveys. From this comprehensive dataset, we have identified 3625 robust galaxy candidates at 𝑧 > 4 . 5, which includes 1276 candidates at 𝑧 > 6 . 5, forming our EPOCHS V1 sample. For the 𝑧 = 4 . 5 -6 . 5 range, only data from the CEERS, JADES GOODS-S, and NEP-TDF  \nsurveys are used due to the availability of HST ACS/WFC F606W and F814W band data.  \nFurthermore, for the NEP-TDF, MACS-0416, and El Gordo observations, we position one of the two NIRCam modules in each cluster pointing on the lensing cluster, while the other module is located approximately 3 arcminutes to the side, providing an effectively 'blank-field' view of the distant Universe. For this paper we do not include sources found in the cluster module to avoid sometimes uncertain corrections with mergers where there is more contamination from the cluster members and uncertainties due to the effects of lensing (e.g., Griffiths et al. 2018; Bhatawdekar & Conselice 2021). For these clusters it is only the NIRCam module which is not centred on the cluster is used.  \nIn this paper, we use our own EPOCHS reduction, with a detailed description given in Adams et al. (2024); Austin et al. (2024); Harvey et al. (2024), and Conselice et al. 2024 (in preparation). We do not repeat this information and process here, but many details are available in these papers. It is important to keep in mind that we re-reduce all these fields using the same methods and furthermore measure and detect our galaxies through a homogeneous method on each of these fields. We have uniformly reprocessed all lowerlevel JWST data products following our modified version of the official JWST pipeline. For detailed information on data reductions and catalogs, please refer to Adams et al. (2024); Austin et al. (2024); Harvey et al. (2024).  \nIn the following subsections, we describe the various fields studied and the relevant properties for this particular paper. The depth of our data, taken in the various fields and JWST filters as well as with ACS, is outlined in Table 1.", '2.1 CEERS': "The Cosmic Evolution Early Release Science Survey (CEERS; PI: S. Finkelstein, PID: 1345, Bagley et al. 2023b) is one of the 13 Director's Discretionary Early Release Science (DD-ERS) programs conducted during the first year of JWST's operation. CEERS features a ∼ 100arcmin 2 NIRCamimagingmosaiccoveringwavelengths from 1 to 5 microns, across 10 NIRCam pointings. Concurrently, CEERS includes imaging with the mid-infrared instrument (MIRI) over four pointings, covering 5 to 21 microns, and multi-slit spectroscopy using the near-infrared spectrograph (NIRSpec) across six pointings. NIRCam imaging data is used in this work. This includes image across seven distinct filters: F115W, F150W, F200W, F277W, F356W,F410M,andF444W,witha5 𝜎 depth of 28.6 AB magnitudes using 0.1 arcsec radius circular apertures. The dataset we use encompasses observations collected during June 2022, accounting for 40% of the total NIRCam area covered for CEERS with the remaining data taken in the latter half of the same year.", '2.2 JADES Deep GOODS-S': "The JWST Advanced Deep Extragalactic Survey (JADES; Rieke et al. 2023a; Eisenstein et al. 2023; Bunker et al. 2023; D'Eugenio et al. 2024; Hainline et al. 2023b) cover both the GOODS-S and GOODS-N fields. In this paper, we focus on the region covered by the JADES DR1 release which is in the GOODS-S field (PI: Eisenstein, N. Lützgendorf, ID:1180, 1210). The observations utilise nine filter bands: F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W, encompassing a spatial extent of 24.4 - 25.8 arcmin 2 . A minimum of six dither points was used for each observation, with exposure times spanning 14-60 ks. As shown in  \nHarvey et al. (2024) and Table 1 the depth of the JADES data within the JWST bands ranges from 29.58 to 30.21, with the deepest band being F277W. This field also includes deep F606W data from ACS.  \nAcross all filter bands, JADES ensures a high level of pixel diversity, thereby significantly reducing the impact of flat-field inaccuracies, cosmic ray interference, and other issues at the pixel level.", '2.3 PEARLS and NEP-TDF': "The North Ecliptic Pole Time-Domain Fields (NEP-TDF) is part of the JWST PEARLS observational program. The Prime Extragalactic Areas for Reionization and Lensing Science project (PEARLS; PIDs 1176, 2738, PI: Rogier Windhorst; Windhorst et al. 2023; Diego et al. 2023; Frye et al. 2023), is a JWST Guaranteed Time Observation (GTO) program. As a Cycle 1 GTO project, the PEARLS team was allocated an exposure time of 110 hours. The project's primary objective is to capture medium-deep NIRCam imaging of blank and cluster fields with an approximate depth of 28-29 AB magnitudes.  \nFour of these targets are located in and around gravitationally lensing galaxy clusters, while one is within a blank field - the NEP. The four clusters include MACS 0416, Clio, and El Gordo. The data acquired from PEARLS includes imaging in seven wide filter bands: F090W,F115W,F150W,F200W,F277W,F356W,andF444W,along with one medium band: F410M. In addition, we also include bluer F606W imaging from the HST Advanced Camera for Surveys (ACS) Wide Field Channel (WFC) in the NEP-TDF from the GO-15278 (PI: R. Jansen) and GO-16252/16793 (PIs: R. Jansen & N. Grogin, see O'Brien et al. 2024) HST programs.  \nFor the lensing fields of SMACS 0723, MACS 0416, and El Gordo, the observations include pointings with one NIRCam module centered on the lensing cluster, while the second module is offset by around 3 arcminutes in a 'blank' region. Although we process both modules in these fields, we do not include sources found in the module containing the lensing cluster in this study.", '2.4 NGDEEP': "The Next Generation Deep Extragalactic Exploratory Public Survey (NGDEEP; ID 2079, PIs: S. Finkelstein, Papovich, and Pirzkal) (Austin et al. 2023a; Bagley et al. 2023a; Leung et al. 2023) is the deepest JWST NIRCam data during the first year of operations. Initially scheduled for late January to early February 2023, NGDEEP's primary focus involved NIRISS Wide Field Slitless Spectroscopy of galaxies within the Hubble UltraDeep Field. However, due to an unforeseen suspension of NIRISS operations during the survey's observation window, only half of the planned observations were executed, with the remainder taken in early 2024. In this work we use exposures from the first epoch of observations. The parallel NIRCam observations from NGDEEP encompass six wide-band photometric filters (F115W, F150W, F200W, F277W, F356W, F444W). Some of our early work in this field was presented in Austin et al. (2023b).", '2.5 GLASS': 'The Grism Lens Amplified Survey from Space (GLASS) observational program (GLASS; PID 1324, PI: T. Treu; Treu et al. 2022) is centered on the Abell 2744 Hubble Frontier Field lensing cluster. GLASS primarily employs NIRISS and NIRSpec spectroscopy to study galaxies lensed by the Abell 2744 cluster. In parallel, NIRCam is used to observe two fields offset from the cluster center. Due to the minimal lensing magnification in these offset areas (Bergamini  \net al. 2023), they effectively serve as blank fields. The observations include seven wide filter bands: F090W, F115W, F150W, F200W, F277W, F356W, and F444W.', '2.6 SMACS J0723-7327 ERO': 'The SMACS-0723 lensing cluster is part of the very first ERO program (PID 2736, PI: Klaus Pontoppidan, Pontoppidan et al. (2022)). The program includes NIRCam, MIRI, NIRISS, and NIRSpec observations, and we use the NIRCam observations which contains six bands: F090W, F150W, F200W, F277W, F356W, and F444W. The exposure time of this observation is set to roughly reach the Frontier depth in F814W (which was reached in 100 ks rather than 4 ks), 1.5 times the F160W depth, and almost 10 times the Spitzer (Werner et al. 2004) IRAC 1+2 depth. Because of the inaccuracies of the photometric redshifts for galaxies within this field around 𝑧 ∼ 10, for the most part the pairs in this field are at the lower or higher redshifts of our range.', '3.1 Photometric Redshift Measures': 'Analyzing the quality and robustness of photometric redshifts is a key aspect of this study, as it directly influences the accuracy of galaxy close-pair selections. In this section, we describe how we measure the photometric redshifts of our samples.  \nWe use EAZY-PY (hereafter EAZY (Brammer et al. 2008)) as our primary photometric SED-fitting code. We employ the "tweak\\_fsps\\_QSF\\_12\\_v3" templates alongside Sets 1 and 4 of the SED templates developed by Larson et al. (2023). These additional templates are optimized for high-redshift galaxies, which exhibit bluer rest-frame UV colors (Topping et al. 2022; Cullen et al. 2023; Nanayakkara et al. 2023; Cullen et al. 2024; Austin et al. 2024) and stronger emission lines (e.g. Withers et al. 2023), characteristics frequently observed in such young and distant galaxies.  \nThe "tweak\\_fsps\\_QSF\\_12\\_v3" templates are generated using the flexible stellar population synthesis (FSPS) package (Conroy & Gunn 2010). Set 1 includes three pure stellar templates from Larson et al. (2023), which are based on the BPASS v2.2.1 stellar population synthesis (SPS) model (Eldridge et al. 2017; Stanway & Eldridge 2018; Byrne et al. 2022) for ages of 10 6 , 10 6 . 5 , and 10 7 years at a fixed metallicity 𝑍 ★ = 0 . 05 Z ⊙ , assuming a standard Chabrier (Chabrier 2003) IMF. Set 4 incorporates nebular continuum and line emission calculated using CLOUDY v17 (Ferland et al. 2017), with an ionization parameter log 𝑈 = -2, gas metallicity 𝑍 𝑔𝑎𝑠 = 𝑍 ★ , Lyman continuum escape fraction 𝑓 esc , LyC = 0, hydrogen density 𝑛 H = 300 cm -3 , and a spherical geometry, excluding Ly 𝛼 emission.  \nWe execute the EAZY-py fitting three times: first with a uniform redshift prior of 0 ≤ 𝑧 ≤ 25, and subsequently with a constrained upper redshift limit of 𝑧 ≤ 4 and 𝑧 ≤ 6. This triple approach allows us to compare the goodness of fit between high and low-redshift solutions, as indicated by the 𝜒 2 values, ensuring robust identification of high-redshift galaxies.  \nTo account for uncertainties in flux calibrations and aperture corrections, as well as potential mismatches between synthetic templates and observed galaxies, we impose a minimum flux uncertainty of 10% (Rieke et al. 2023b). Our selection criteria rely predominantly on specific cuts in computed quantities to minimize biases and incompleteness. Although we perform a visual review of cutouts and SED-fitting solutions, this step results in less than 5% rejection of our  \noriginal sample, a significantly lower rate compared to other studies (e.g., Hainline et al. 2024). We refer readers to §3 of Adams et al. (2024), §2.5 of Austin et al. (2024), §3.2 of Harvey et al. (2024), and Conselice et al. (2024, in prep) for details of our selection criteria for robust high-redshift galaxies.  \nThe photometric redshift quality of the EPOCHS sample has been investigated previously (Adams et al. 2024; Duan et al. 2024) by comparison with spectroscopically confirmed sources from JWSTs NearInfraRed Spectrograph (NIRSpec; Ferruit et al. 2022; Jakobsen et al. 2022; Böker et al. 2023). We give a brief summary of the quality of the data here. In Adams et al. (2024), we matched our photometric redshifts with spectroscopic redshifts for five fields-CEERS, JADES GOODS-S, SMACS-0723, MACS-0416, and El-Gordo-finding a ∼ 10% outlier fraction, with the largest discrepancies occurring around 2 < 𝑧 < 3. In Duan et al. (2024), we investigated the properties of 43 galaxies at 𝑧 > 7 with both NIRCam and NIRSpec PRISM observations in the CEERS and JADES GOODS-S fields. This investigation revealed an outlier fraction of 4.7%, defined as z phot > 1 . 15 ( z spec + 1 ) or z phot < 0 . 85 ∗ z spec + 0 . 15 (Duan et al. 2024), and a normalized median absolute deviation (NMAD) of 0.035. These results provide strong evidence of the high quality of our photometric redshifts, thereby ensuring the accuracy of closepair selections in this study.', '3.2 Stellar Mass Measures': 'To infer galaxy stellar masses from the photometry we use the Bagpipes (Carnall et al. 2018) bayesian SED fitting code, adopting redshift priors from the EAZY redshift probability distribution function for each galaxy. We use Log10 priors for dust, metallicity, and age. The reason for selecting Log10 priors is because we expect high redshift galaxies to be young, with lower metallicity and to contain little dust. We set prior limits for metallicity in the range of [ 5e-03 , 5 ] Z ⊙ , and we use a dust prior in the range of [ 0 . 0001 , 10 . 0 ] in 𝐴 V , ionization parameter Log 10 ( U ) in range of [-4 , -1 ] . The time assumed in our fitting models for star formation to start is at 0.001 Gyr, while the time assumed for star formation to stop is at 𝜏 U , with 𝜏 U denoting the age of the Universe at that redshift. In addition, Kroupa (2001) IMF, Bruzual & Charlot (2003) SPS model, and the Calzetti et al. (2000) dust attenuation model is implemented.  \nWe fit the galaxy data with six parametric SFH models - lognormal, delayed, constant, exponential, double delayed, and delayed burst-along with a non-parametric Continuity model (Leja et al. 2019). The integration of the SFH is the stellar mass. These models have been used in various works (e.g. Carnall et al. 2023; Endsley et al. 2023; Tacchella et al. 2023; Whitler et al. 2023). Different SFHs can have a large impact on inferred stellar mass and SFR (Furtak et al. 2021; Tacchella et al. 2022; Harvey et al. 2024; Wang et al. 2024). We have chosen to present results inferred from our fiducial model, the log-normal SFH, because we expect high-redshift galaxies to exhibit rising or bursting star formation. Please refer to Harvey et al. (2024) for a thorough discussion of the topic of priors in SED fitting.  \nAs we employ an integrated probability method that incorporates the entire photometric redshift probability distribution to compute the evolution of the merger history, we compute the stellar mass for each galaxy not only at the best-fit redshift but across all redshifts within the redshift posterior. For each galaxy, we take 400 independent draws fromthebest-fit posterior from the SED, each of which has a synthetic spectra for the galaxy, and taken from this the associated galaxy properties such as redshift and stellar mass. This lets us parameterise M ∗ ( 𝑧 ) , the best-fit galaxy stellar mass as a function of redshift. It is important to note that we give Bagpipes an informative redshift  \nTable 1. 5 𝜎 depths and unmasked areas of the eight JWST Cycle-1 fields used in this study. For CEERS, JADES GOODS-S, NEP-TDF, and NGDEEP, we also incorporate existing HST ACS/WFC observations. CEERSP9 is listed separately from the other 9 pointings in CEERS as it is significantly deeper than the other pointings.  \nFigure 1. Example of stellar mass M ∗ ( 𝑧 ) and redshifts 𝑃 ( 𝑧 ) distributions for an example galaxy from the NEP-TDF field, generated using Bagpipes and EAZY . The blue line represents the M ∗ ( 𝑧 ) distribution, while the pink line depicts the 𝑃 ( 𝑧 ) distribution. We compute the values of M ∗ ( 𝑧 ) only at redshifts where the value of 𝑃 ( 𝑧 ) is non-zero.  \n<!-- image -->  \nprior, taken from the highest likelihood draw from the EAZY posterior, which limits the possible range of redshifts Bagpipes will fit. An example of M ∗ ( 𝑧 ) together with its normalised redshift probability distribution for one of the NEP-TDF galaxies is shown in Figure 1.', '4 GALAXY PAIR FRACTION METHODOLOGY': 'The primary objective of this study is to accurately determine the fraction of galaxy close-pairs at high redshift, denoted as f p . Subsequently, we derive merger rates and mass accretion rates from this fraction. Accurately computing f p requires a statistically robust approach, particularly in light of the challenges presented by the photometric uncertainty, inherent incompleteness of galaxy samples, and the prevalent bias towards the brightest high-z galaxies. To address these issues, we have developed a galaxy pair selection methodology, building upon the foundational work of López-Sanjuan et al.  \n(2015), Mundy et al. (2017), Duncan et al. (2019), and Conselice et al. (2022). This method employs the integrated probability approach. Instead of using the best-fit redshifts of galaxies to search for mergers, we utilize their entire redshift probability distributions. The steps are described in Sections 4.1 - 4.3. Each step is carried out at every redshift, and the results are integrated to derive the final pair fractions as detailed in Section 4.4.', '4.1 Pair Probability Function': 'In this section, we illustrate the three close-pair selection criteria and describe how they combine to establish the final Pair Probability Function, PPF( 𝑧 ) , at each redshift.', '4.1.1 Projection Separation Selection': 'The initial selection of close pairs in our study is based on the projection separation. We consider the physical projection separation, 𝑟 𝑝 , between two galaxies, which is a function of their angular separation 𝜃 and the angular diameter distance 𝑑 𝐴 ( 𝑧 ) at redshift 𝑧 . The projected separation 𝑟 𝑝 is thus mathematically expressed as:  \n𝑟 𝑝 ( 𝑧 ) = 𝜃 × 𝑑 𝐴 ( 𝑧 ) . (1)  \nSubsequently, the projection separation mask is defined as follows:  \nM 𝑟 𝑝 ( 𝑧 ) = ( 1 if 20 kpc < 𝑟 𝑝 ( 𝑧 ) < 50 kpc 0 otherwise , (2)  \nwhere a value of 1 signifies that the pair falls within the desired range for projected separation, thus meeting our first close pair selection criteria.  \nThis selected range is strategically chosen. The lower limit of 20 kpc excludes very closely bound systems where galaxies are in advanced merger stages, often requiring additional complex analysis such as the CAS (Concentration, Asymmetry, Smoothness) system to identify (Abraham et al. 1994; Schade et al. 1995; Abraham et al. 1996; Bershady et al. 2000; Conselice 2003; Conselice et al. 2003). The upper limit of 50 kpc is set to omit pairs that are too widely separated to have a chance to merge soon at larger separations. The galaxy separation chosen in this work are different from those in the previous works of Mundy et al. (2017), Duncan et al. (2019),  \nand Conselice et al. (2022), and the logic behind this change will be explained. In Section 5.1.1, we will explore the effect of different 𝑟 𝑝 ranges on the derived close pair fractions. What we find is that the pair fractions are in anycase very similar when we use a lower limit of 5kpc (see Figure 5a).', '4.1.2 Line of Sight (Redshift) Selection': 'For each close pair that meets the projection separation criteria, we designate the more massive galaxy as the primary galaxy and the less massive one as the secondary galaxy. The combined redshift probability function, Z( 𝑧 ) , is calculated to assess the likelihood that the pair is at the same redshift and, therefore, represents a true merger pair. This function is defined as:  \nZ( 𝑧 ) = 2 𝑃 1 ( 𝑧 ) 𝑃 2 ( 𝑧 ) 𝑃 1 ( 𝑧 ) + 𝑃 2 ( 𝑧 ) = 𝑃 1 ( 𝑧 ) 𝑃 2 ( 𝑧 ) 𝑁 ( 𝑧 ) , (3)  \nwhere 𝑃 1 ( 𝑧 ) and 𝑃 2 ( 𝑧 ) are the normalised redshift probability distribution functions for the primary and secondary galaxies, respectively, and 𝑁 ( 𝑧 ) serves as a normalization factor such that Z( 𝑧 ) is properly scaled.  \nThe combined redshift probability function Z( 𝑧 ) incorporates the individual probability distributions of both galaxies, taking into account the uncertainty in their redshift determinations. This approach allows us to statistically account for the probability that both galaxies in a pair are at a similar redshift, which is a fundamental requirement for identifying potential mergers. However, we are careful to not double count pairs when examining the pair fraction.  \nAninteresting quantity that can be derived from Z( 𝑧 ) is N 𝑧 , which is calculated by integrating Z( 𝑧 ) across the entire redshift range:  \nN 𝑧 = ∫ ∞ 0 Z( 𝑧 ) 𝑑𝑧, (4)  \nwhere N 𝑧 values range between 0 and 1. This parameter quantifies the probability that a pair of galaxies in a potential merger system is a real pair. N 𝑧 can be very useful for assessing the proximity of two galaxies along the line of sight direction, even though it is not directly utilized in the derivation of PPF( 𝑧 ) .  \nAs each target galaxy can have more than one close companion, each potential galaxy pair is analysed separately and included in the total pair count. This ensures that our analysis accounts for potential multiple mergers for each galaxy.  \nWe demonstrate the impact of Z( 𝑧 ) and N( 𝑧 ) on the selection of close pairs in Figure 2a and Figure 2b. These figures present two close-pair systems from the NEP-TDF and JADES GOODS-S fields. In the example from NEP-TDF, two galaxies have a physical projection separation of 41.3 kpc and N( 𝑧 ) of 0.98, satisfying the projection separation criteria (20 kpc - 50 kpc) in this work, and have a very high probability of being a real pair. Conversely, in the JADES GOODS-S example, two galaxies have a physical projection separation of 6.8 kpc and N( 𝑧 ) of 0.02, which do not meet the projection separation criteria, and have a very low probability of being a real pair.', '4.1.3 Stellar Mass Selection': 'We have developed a modified version of the stellar mass selection criteria, building upon López-Sanjuan et al. (2015); Mundy et al. (2017); Duncan et al. (2019); Conselice et al. (2022). The primary distinction lies in the aspect of sample completeness. The aforementioned papers apply a completeness correction, employing a redshift-dependent mass-completeness limit influenced by the flux  \nFigure 2. Example of N( 𝑧 ) for two close pair systems from the NEP-TDF and JADES GOODS-S fields. The blue and red lines represent the normalized probability distributions of the galaxy redshifts within each pair. 𝑧 med is the best fit photometric redshift. The black line denotes the integrated Z( 𝑧 ) , with N( 𝑧 ) displayed in the bottom right corner of each panel. (a): Two galaxies with an angular separation of 7.32 arcsec and a physical projection separation of 41.3 kpc, with N( 𝑧 ) = 0 . 98, satisfy the projection separation criteria (20 kpc - 50 kpc) in this work and have a very high probability of being a real pair. (b): Two galaxies with an angular separation of 1.18 arcsec and a physical projection separation of 6.8 kpc, with N( 𝑧 ) = 0 . 02. This system does not meet the projection separation criteria and has a very low probability of being a real pair.  \n<!-- image -->  \nlimit set by the survey to accommodate varying survey depths across different fields. However, as can be seen from Table 1, there are no big differences in survey depths across the eight JWST fields we are using. What is crucial to consider is that we are likely missing faint high-redshift galaxies, which could reduce the observed number of close-pairs. Thus, we have developed our novel completeness correction methodology, which is outlined in Section 4.3. In this section, we apply a stellar mass selection to ensure the general completeness of the sample, which is log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0, as indicated in Figure 3, before proceeding to the completeness correction step.  \nThe stellar mass selection mask is defined as:  \nM M ( 𝑧 ) =         1 if M lim , 1 ∗ ( 𝑧 ) ≤ M ∗ , 1 ( 𝑧 ) ≤ M max ∗ and M lim , 2 ∗ ( 𝑧 ) ≤ M ∗ , 2 ( 𝑧 ) 0 otherwise , (5)  \nwhere M ∗ , 1 ( 𝑧 ) and M ∗ , 2 ( 𝑧 ) represent the stellar masses of the primary and secondary galaxies at a given redshift, respectively. The primary galaxy is the higher mass one in a potential pair. The definitions of M lim , 1 ∗ ( 𝑧 ) and M lim , 2 ∗ ( 𝑧 ) are given by  \nM lim , 1 ∗ ( 𝑧 ) = M min ∗ , (6)  \nand  \nM lim , 2 ∗ ( 𝑧 ) = 𝜇 M ∗ , 1 ( 𝑧 ) . (7)  \nThe parameter 𝜇 represents the mass ratio, where values greater than 𝜇 characterize major mergers. In this study, we have set 𝜇 to a value of 1 / 4 and defined major mergers as those with 𝜇 > 1 / 4. The stellar mass range of interest is defined by M max ∗ and M min ∗ . In this work, we use the range log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0. The last mask is the primary galaxy selection mask, which is the first equation of Equation 5, given by:  \nS( 𝑧 ) = ( 1 if M lim , 1 ∗ ( 𝑧 ) ≤ M ∗ , 1 ( 𝑧 ) ≤ M max ∗ 0 otherwise . (8)', '4.1.4 Pair Probability Function': 'The final pair probability function at a given redshift, PPF( 𝑧 ) , is derived by combining three selection criterias: projection separation, line of sight, and stellar mass. It is defined as:  \nPPF( 𝑧 ) = M 𝑟 𝑝 ( 𝑧 ) × Z( 𝑧 ) × M 𝑀 ( 𝑧 ) . (9)  \nIt is important to reiterate that the above processes are performed for each redshift within the entire photometric redshift fitting range. That is, at every redshift, we first generate an initial merger catalog based on projection separation criteria, then apply a mass selection mask, and compute Z( 𝑧 ) to produce the PPF( 𝑧 ) of each pair at that redshift. This is why we compute masses at all redshifts where the probability is non-zero for each galaxy (Section 3.2). In the following sections, we will discuss two additional parameters that, when combined with PPF( 𝑧 ) , will be integrated over the redshift bin of interest to compute pair fractions.', '4.2 Correction for Selection Effects': 'The PPF( 𝑧 ) function defined in Equation 9 is affected by two selection effects:  \n- (I) The incompleteness in the search volume near the boundaries of the images. (Section 4.2.1)\n- (II) Photometric redshift quality differences between survey regions. (Section 4.2.2)  \nWedetailthecorrections applied to address these issues in this section and discuss their implementation in the context of JWST data.', '4.2.1 Incompleteness Due to Image Borders and Masked Areas': "Primary galaxies located close to image boundaries or masked regions (e.g., bright stars) possess a reduced spatial search area for potential companions, consequently reducing the observed number of pairs. Since the search area is contingent upon a fixed physical  \nsearch radius, 𝑟 𝑝 , this correction factor varies with redshift and necessitates computation at each redshift of interest.  \nFor every primary galaxy, at each redshift, we delineate an annular region defined by an inner circle with a physical radius of 20 kpc and an outer circle with a radius of 50 kpc (projection separation criteria defined in 4.1.1), centered on the galaxy. We then calculate the fraction of the annular's area, 𝑓 area ( 𝑧 ) , that does not lie outside the image boundaries or within masked regions. Subsequently, each secondary galaxy is assigned a weight inversely proportional to the fraction of the area unmasked or inside the image. This weight is defined as:  \n𝑤 area ( 𝑧 ) = 1 𝑓 area ( 𝑧 ) . (10)  \nThis approach inherently compensates for the irregular geometries of survey areas and minor computational inaccuracies arising from finite pixel sizes. The impact of this correction to pair fractions is ∼ 7% -15%.Thisis analogous to the volume correction in luminosity functions.", '4.2.2 Photometric Redshift Quality Correction': "Theoddsparameter, O , conceptualized by Benítez (2000) and Molino et al. (2014), quantifies the confidence of a galaxy's photometric redshift falling within a certain interval around the best-fit value. This parameter has been effectively utilized in studies by LópezSanjuan et al. (2015) and Mundy et al. (2017).  \nAs analyzed in Section 3, our photometric redshifts exhibit exceptional quality, with only 4.7% outliers, and an NMAD value of 0.035 at 𝑧 > 7. Consequently, we do not use the O parameter in our analysis. Similarly, Duncan et al. (2019) did not employ the O parameter, attributing this to the comprehensive magnitude-dependent calibration of photometric redshift posteriors, which ensures that the redshift distributions are accurately calibrated across all magnitudes.", '4.3 Completeness Correction for the JWST High-z Sample': "TheincompletenessoftheJWSThigh-redshiftsampleisevidentfrom the mass distribution presented in Figure 3. A comprehensive sample would feature an increasing number of galaxies at lower masses, not the Gaussian-like distribution observed. More importantly, surface brightness dimming at 𝑧 = 4 . 5 -11 becomes significantly high. This results in a higher likelihood of missing faint companions at high redshifts. Consequently, the only objects for which pairs can be counted at 𝑧 > 4 . 5 are those with companions that are point sources with very high surface brightness. We aim to correct this effect and retrieve those missing faint companions. The details of our completeness correction methodology are outlined in this section.  \nAfter applying a stellar mass selection of log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0 to ensure general completeness, as described in Section 4.1.3, we proceed to correct for the incompleteness of each individual galaxy. The method we developed for correction involves simulating different incompleteness scenarios with low-redshift data ( 𝑧 < 6) and investigating its impact on the derived pair fractions. We employ the same pair fractions methodology and use the same HST CANDELS catalogs from Duncan et al. (2019) to calculate pair fractions under various scenarios of incompleteness. The CANDELS dataset encompasses five fields with a total number of 181,086 galaxies: GOODS-S (Guo et al. 2013), GOODS-N (Barro et al. 2023), COSMOS (Nayyeri et al. 2017), UDS (Galametz et al. 2013), and Extended Groth Strip (EGS) (Stefanon et al. 2017). This dataset is complete for log 10 ( M ∗ / M ⊙ ) > 10 . 3 and has undergone multiple quality checks  \nand flag applications to ensure its cleanliness and reliability, thus enabling us to trace pair fraction evolution up to 𝑧 = 6 . 0.  \nWe simulate incompleteness scenarios by systematically removing portions of the complete lower redshift sample, ranging from 10% to 90%, with a step size of 10%, thus generating catalogs with varying degrees of completeness from 0.9 to 0.1. For each level of completeness, we implement the pair fraction methodology defined above. This allows us to evaluate the impact of sample completeness on the pair fractions. To minimize random errors due to the random exclusion of galaxies, we repeat the entire procedure five times. The data from all fields are utilized for 𝑧 < 3 . 5, while for 3 . 5 < 𝑧 < 6 . 5, we exclusively use data from the HST GOODS-S field, noted for its geatest depth among the five.  \nUpon calculating the pair fraction, f p , across different redshift bins under nine scenarios of incompleteness, we plot the f ratio p against completeness for each redshift bin, as illustrated in Figure A2. Here, the f ratio p is defined as the pair fraction at a given completeness level divide by the complete sample case (completeness = 1). The resulting trend, indicative of a direct linear correlation between increasing completeness and f ratio p , is modeled with a linear fit (Equation 11) using the Bayesian Markov Chain Monte Carlo (MCMC) method and emcee package. Specifically, we employ 100,000 steps and 50 walkers to generate candidate gradients and y-intercept values. For both sets of values, we adopt the median as the representative value and use the 1 𝜎 deviation as the associated uncertainty, as the distributions follow a perfect Gaussian.  \nf ratio p = Gradient ( 𝑧 ) × Completeness ( 𝑧, M ∗ ) + Y-Intercept ( 𝑧 ) , (11)  \nSubsequently, we chart the evolution of the best-fit lines' gradients and y-intercepts over different redshift bins in Figure A3. This reveals an increasing trend in the gradient at higher redshifts and a decreasing trend in the y-intercept. These trends are further fitted linearly, with the best-fit parameters presented in the caption of Figure A3, using the MCMC method described above.", '4.3.1 The C( 𝑧, M ∗ ) correction factor': 'The completeness correction factor, C( 𝑧, M ∗ ) , depends on both mass and redshift. We compute the completeness level of JWST galaxies using the JAGUAR simulation (Williams et al. 2018), with details outlined in Appendix A and in Figure A1. Additionally, we obtain the gradient and y-intercept information for the best-fit line (Equation 11) at different redshifts based on Figure A3. By substituting these three parameters: gradient, completeness, and y-intercept to Equation 11, we derive the f ratio p , which ranges between 0 and 1. The completeness correction factor is then defined as:  \nC( 𝑧, M ∗ ) = 1 f ratio p ( Completeness ( 𝑧, M ∗ )) . (12)  \nWethenimplementthiscorrectioninto our pair fraction measurement as described in the next section.', '4.4 Final Integrated Pair Fractions': 'Integrating the PPF( 𝑧 ) , weights and C( 𝑧, M ∗ ) discussed in Sections 4.1, 4.2, and 4.3, we compute the final integrated pair fractions within each redshift interval of interest, denoted by [ 𝑧 min , 𝑧 max ] . Recall the more massive galaxy in each pair is designated as the primary, with the less massive serving as the secondary.  \nLetting 𝑖 index each primary galaxy and 𝑗 index every potential  \nsecondary galaxy surrounding the primaries, we express the number of pairs for each primary galaxy within each redshift interval as  \n𝑁 pair ,𝑖 = ∑︁ 𝑗 ∫ 𝑧 max 𝑧 min C( 𝑧, 𝑀 ∗ ) × PPF 𝑗 ( 𝑧 ) × 𝑤 area ( 𝑧 ) 𝑑𝑧 (13)  \nwhere PPF( 𝑧 ) is the pair probability function detailed in Section 4.1, and C( 𝑧, 𝑀 ∗ ) is the completeness correction factor introduced in Section 4.3. Note that usually the summation over 𝑗 is just for one single pair, with multiple pairs and higher 𝑗 values relatively rare.  \nEach primary galaxy\'s contribution to the overall distribution of primary galaxies within each redshift range is given by:  \n𝑁 1 ,𝑖 = ∫ 𝑧 max 𝑧 min P 1 ,𝑖 ( 𝑧 ) × S 1 ,𝑖 ( 𝑧 ) 𝑑𝑧 (14)  \nwhere P 1 ,𝑖 ( 𝑧 ) represents the normalized redshift distribution function for each primary galaxy, S 1 ,𝑖 ( 𝑧 ) denotes the selection mask for primary galaxies as defined in Section 4.1.3.  \nFinally, we define the pair fraction for each redshift interval as:  \nf p = " 𝑖 𝑁 pair ,𝑖 " 𝑖 𝑁 1 ,𝑖 (15)  \nThis is then the pair fraction which we use through the remainder of the paper and which we also use to calculate the merger rate and the merger mass accretion rate.', '4.4.1 Uncertainty in Pair Fractions': 'To compute the statistical uncertainties associated with the pair fraction, we employ the bootstrap resampling method (Efron 1979; Efron & Stein 1981). The standard deviation of the pair fraction, denoted as 𝜎 f p , is defined by the equation:  \n𝜎 f p = √︄ " 𝑁 𝑖 = 1 ( 𝑓 𝑚,𝑖 - ⟨ f p ⟩) 2 𝑁 -1 , (16)  \nwhere 𝑓 𝑚,𝑖 represents the pair fractions estimated from a randomly selected sample of galaxies, 𝑁 denotes the number of realizations, and ⟨ f p ⟩ is the average pair fraction, calculated as ⟨ f p ⟩ = 1 𝑁 " 𝑁 𝑖 = 1 f p, i . During the bootstrap process, resampling with replacement is accounted for. We use 100,000 iterations to ensure statistical stability. In fact, we find the pair fractions to be stable after 10,000 iterations.', '4.5 Simulation Data Pair Counts': 'Comparing observational results with simulations is imperative, as this can enhance our understanding of cosmic structures. We derive pair fractions using the Semi-analytic Forecasts for JWST (SC-SAM) simulation catalog. The "Semi-analytic Forecasts for the Universe" project is an ongoing project that offers physically grounded predictions regarding the characteristics and distributions of galaxies and quasars throughout cosmic history. These forecasts are generated using the well-established Santa Cruz semi-analytic model (SAM) for galaxy formation. This model stands out for its balance of physical rigor and computational efficiency, making it a valuable tool for understanding galactic evolution. Notably, the model has demonstrated its capability to reliably reproduce observational constraints up to a redshift of approximately 𝑧 ∼ 10 (Yung et al. 2019a).  \nThe project encompasses two primary data series: Semi-analytic Forecasts for JWST and Semi-analytic Forecasts for Roman . For our research, we utilize the Semi-analytic Forecasts for JWST data product. The methodology of the simulation, along with its comprehensive details, is thoroughly documented in a series of studies  \nFigure 3. Histogram displaying the distribution of best-fit stellar masses for all 1276 galaxies at 𝑧 > 6 . 5 and for 3452 galaxies at 𝑧 > 4 . 5 (an additional 2376 galaxies from the CEERS, JADES GOODS-S, and NEP-TDF fields). These samples are from eight JWST fields: CEERS, JADES GOODS-S, NEP-TDF, NGDEEP, GLASS, El-Gordo, SMACS-0723, and MACS-0416. The median stellar masses are found to be 7.89 [log 10 ( M ∗/ M ⊙) ] and 7.92 [log 10 ( M ∗/ M ⊙) ] for each group, respectively. At the same time, it is clear that there is a drop in completeness at log 10 ( M ∗/ M ⊙) < 8 which is also seen in image simulations of these fields (Harvey et al. 2024).  \n<!-- image -->  \n(Yung et al. 2019a,b, 2020a,b, 2021, 2022). This simulation comprises three products: wide , deep , and ultra-deep . We use the deep dataset, as it overlaps with the HUDF, where the JADES GOODS-S field is located.  \nThe methodology for calculating pair fractions from the SC-SAM catalogs differs from that used with actual observational data. In the simulation, precise redshift values, 𝑧 , are available, in contrast to the redshift probability distribution function, 𝑃 ( 𝑧 ) , used in observational analyses, thereby simplifying the pair counting process.  \nThe first step involves filtering the catalog. The simulation includes galaxies fainter than JWST\'s observational limit, which could potentially contaminate the pair fraction calculation. The median m F444W of our JWST sample is 27 . 8, with a 1 𝜎 spread of 1 . 05. Consequently, we impose a magnitude cutoff of 𝑚 F444W < 30 on the SC-SAM catalog, utilizing the catalog\'s NIRCam F444W dust column for this purpose. Additionally, we apply a mass selection criterion of log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0 to the SC-SAM catalog, ensuring consistency with the completeness selections of our highredshift JWST galaxies. The deep dataset contains eight different realizations. After applying magnitude and mass selections, a total of 522,249 galaxies remain across all eight realizations.  \nNext, we select galaxy pairs from the SC-SAM catalog based on the projection separation criteria: 20 kpc < 𝑟 𝑝 ( 𝑧 ) < 50 kpc. For redshift (line-of-sight) selection, we employ the velocity difference upper limit method. The rest-frame velocity difference between two galaxies is given by  \nΔ 𝑣 = 𝑐 × | 𝑧 1 -𝑧 2 | ( 1 + 𝑧 mean ) , (17)  \nwhere 𝑧 1 and 𝑧 2 are the redshifts of the two galaxies, respectively, and 𝑧 mean is their average redshift. While Ventou et al. (2017) use an upper limit of 500 km/s, higher limits have also been used, such as 1000 km/s (Endsley et al. 2020) and 10,000 km/s (Gupta et al. 2023). We calculate pair fractions at various velocity upper limits  \nscenarios, with the highest value set at 2000 km/s, to better match what we are likely obtaining with our photo-zs. We find that there are not big differences in the derived pair fractions using different velocity upper limits, consistent with the findings of Duncan et al. (2019).  \nFor each identified pair, we continue to apply the same definitions as in the observational data: the primary galaxy is the one with the greater mass in the pair, and the secondary galaxy is the less massive one. Only pairs associated with primary galaxies are counted to avoid over-counting identical pairs. We then calculate the pair fraction for each redshift bin using the formula:  \nf p = 𝑁 pair 𝑁 total , (18)  \nwhich reflects the ratio of the number of pairs to the total number of galaxies in each redshift bin. The calculation of the errors in the pair fraction is performed using Poisson statistics.  \nAs SC-SAM simulation deep dataset comprises eight different realizations, we have conducted analyses on all these realizations and averaged the pair fractions derived from these multiple datasets to compute the final pair fractions for every redshift bin.', '5 HIGH-REDSHIFT MERGER EVOLUTION AT Z > 6': 'In this section, we present the first close-pair analysis of the pair fraction, merger rate, and mass accretion rate evolutions at 𝑧 = 4 . 5 -11 . 5. Along with our high-redshift JWST data points, we also include lower redshift results from the literature to trace the evolution from the Early Universe to the Local Universe. We present catalogs and cutout images of high-redshift close-pair merger systems in Appendix B.', '5.1 Pair Fraction': 'Wecompute pair fractions across eight JWST fields and the SC-SAM simulation for 𝑧 = 4 . 5 -11 . 5, using the method detailed in Section 4. Considering the limited field area of each JWST field and the possibility of overdense regions to potentially elevate pair fractions (Pearson et al. 2024), we have combined the results from all eight fields to mitigate these effects.  \nIn certain redshift bins within some fields, the samples are almost all incomplete to some extent (Completeness < 0.2). If we blindly follow the completeness correction methodology outlined in Section 4.3, it will produce an extremely large completeness correction factor, C( 𝑧, M ∗ ) , which will greatly contaminate the derived pair fractions at very low completeness levels. As a result, these redshift bins are discarded and not included in the final combined f p derivation due to their being very incomplete and the corrections for these pair fractions is very high. These discarded redshift bins are shown in Table 2. Furthermore, we find a notably limited number of close pairs in SMCAS-0723, GLASS, MACS-0416, and El-Gordo fields due to small field area as shown in Table 1. Consequently, our calculation of f p is dominated by data from CEERS, JADES GOODS-S, NEP-TDF, and NGDEEP fields, where the more substantial number of pairs were detected. The completeness values are discussed in more detail in the appendix and are shown graphically in Figure A1. Most of our corrections for this completeness in the pair fraction are relatively modest. We only plot and analyze these relatively complete bins in redshift and stellar mass for our analysis.  \nThe evolution of the major merger pair fraction (f p ) with redshift is presented as yellow orange stars in Figure 4a. The values of pair  \nfractions for each redshift bin are shown in Table 3. In addition to our JWST high redshift points, we include lower redshift results derived from (Casteels et al. 2014; López-Sanjuan et al. 2015; Duncan et al. 2019; Conselice et al. 2022). It is also worthwhile to include many great lower redshift works (e.g., López-Sanjuan et al. 2013; Man et al. 2016; Snyder et al. 2017; Ventou et al. 2017; Mundy et al. 2017; Mantha et al. 2018; Kim et al. 2021; Inoue et al. 2024). However, considering the JWST data and methodology used in this work, it would be more appropriate to integrate with lower-redshift results that utilize similar datasets (e.g., HST-based) and methodologies to produce pair fraction evolution from 𝑧 ∼ 11 . 5 to the local Universes.  \nWe observe a gradual increase in the major merger pair fraction up to 𝑧 ∼ 8, reaching a value of 0 . 211 ± 0 . 065. Although the last two points at 𝑧 = 9 . 0 and 𝑧 = 10 . 5 are lower than the point at 𝑧 = 8 . 0, these two points have high uncertainty. Taking this uncertainty into account, we conclude that the overall evolution of pair fractions shows a flat trend for 𝑧 > 6. We treat the point at 𝑧 = 6 . 0 as an outlier, as it is noticeably high. The reason for this is that among the three largest fields that dominate the f p values - CEERS, JADES GOODSS, and NEP-TDF - the JADES GOODS-S and NEP-TDF galaxies exhibit a very high degree of incompleteness at 𝑧 = 5 . 5 -6 . 5 (Figure A1,Table 2), with the majority of galaxies having a completeness of less than 0.1, resulting in their exclusion. Thus, the f p at 𝑧 = 6 . 0 is dominated by the CEERS point, which is inherently high and could not be averaged down due to the exclusion of the two fields that could have balanced that value. The decline in f p at 𝑧 > 8 is perhaps because we are probing a smaller volume, which could reduce the pair fractions. The average f p between 𝑧 = [ 4 . 5 -11 . 5 ] is 0 . 17 ± 0 . 05.  \nThis indicates that f p at higher redshifts does not change as dramatically as observed at lower redshifts. Comparatively, the morphology method by Dalmasso et al. (2024) finds no evidence of f p evolutions between 𝑧 = [ 4 . 0 -9 . 0 ] , with an average f p = 0 . 11 ± 0 . 04. This value is approximately 0.06 lower than what we computed using the close-pair method, showing that the merger timescales of these two methods are different and/or that not all close-pair galaxies will ultimately merge into a single galaxy.  \nLiterature has established that fitting a Power-Law or a Power-Law + Exponential model to the evolution of pair fractions with redshift is appropriate (e.g., Patton et al. 2002; Conselice 2003; Bridge et al. 2007; Conselice et al. 2008; Man et al. 2016; Ventou et al. 2017; Mantha et al. 2018; Duncan et al. 2019; Conselice et al. 2022), expressed as:  \nf p = 𝑓 0 × ( 1 + 𝑧 ) 𝑚 , (19)  \nand  \nf p = 𝑓 0 × ( 1 + 𝑧 ) 𝑚 × 𝑒 𝜏 ( 1 + 𝑧 ) . (20)  \nWe apply both a power-law and a power-law + exponential model to fit our data. The fitting is performed using only our JWST data, supplemented by a zero-point data from Casteels et al. (2014), represented as a blue triangle. The fitted lines are shown as black and red lines in Figure 4a, and the fitted parameters are detailed in Table 4. The observation of a gradual increase in f p up to 𝑧 = 8, followed by a statistically flat trend towards 𝑧 = 11 . 5, indicates that a power-law + exponential model is better suited to describe the evolution of f p with redshift. This can be directly observed from the figure, as the power-law fitting results in a continuous increase in pair fractions even at very high redshifts ( 𝑧 > 8), where a flat trend is expected. Furthermore, we find that the fitted power-law + exponential line aligns very well with established literature, even though they are not included in the fitting. The higher f p values from VIDEO, GAMA, and UltraVISTA by Conselice et al. (2022) are because this work  \nTable 2. Redshift bins discarded due to the presence of a large number of incomplete samples, to prevent contamination of pair fractions.  \npresents total pair fractions, which are, by definition, higher than major pair fractions, as pair fractions include both major and minor pairs.  \nIt is also arguable that the stellar mass range of the literature results we used to compare and fit with our results are not consistent. We use lower redshift results with a generally higher stellar mass selection than our stellar mass range. However, this is simply because galaxies at high redshifts are less massive, and we are integrating our corrected mass-completed high-redshift sample with other mass-completed results at lower redshifts. We have done our best to normalize out the incompleteness with this, and as long as the pair fraction does not change significantly with mass at a given redshift without our range, then these comparisons are fine for a first examination of this merger evolution.  \nMinor mergers (1 / 10 < 𝜇 < 1 / 4) are observed to be more dominant than major mergers at low redshift (e.g., Duncan et al. 2019; Pearson et al. 2019; Conselice et al. 2022). However, minor mergers are harder to detect because they are less luminous and more compact in size, especially at our redshift range 𝑧 > 6 where galaxies are so distant. When using our method to compute minor merger pair fractions, we obtain values 3.5 times lower than those for major mergers, which is inaccurate due to the omission of many faint close-pair systems. Therefore, in this work, we do not present minor merger pair fractions and focus solely on major mergers.', '5.1.1 Pair Fractions with Varying Projection Separation Criteria': 'In Section 4.1.1, we discuss the physical projection separation criteria used to select initial merger samples, with the range being 𝑟 𝑝 = [ 20 , 50 ] kpc. However, the optimal range for this criterion is debatable and still uncertain; some studies set an upper limit of ∼ 30 kpc (Man et al. 2016; Mundy et al. 2017; Ventou et al. 2017; Duncan et al. 2019; Conselice et al. 2022), while others opt for 50 kpc (López-Sanjuan et al. 2015; Snyder et al. 2017; Gupta et al. 2023). Therefore, investigating the impact of different projection separation criteria on pair fractions at 𝑧 = 4 . 5 -11 . 5 is crucial.  \nIn Figure 5a, we present pair fractions derived from five different projection separation criteria- with inner, outer radii in kpc: [ 5 , 10 ] , [ 5 , 20 ] , [ 5 , 30 ] , [ 5 , 40 ] , [ 5 , 50 ] kpc-for 𝑧 = 4 . 5 -11 . 5, while maintaining consistency in all other methodologies. The result from our fiducial separation criteria, [ 20 , 50 ] kpc, is shown as the black dotted line. We observe that the f p derived from all projection separations exhibits a similar trend-a gentle increase in the first few data points followed by a statistically flat trend or gradual decline at higher redshifts. Additionally, as expected, smaller projection separation measures result in lower pair fractions. For every 10 kpc reduction in the projection separation search, f p decreases by a factor of 1 . 72 ± 0 . 47.  \nTable 3. Pair Fractions, Merger Rates, and Mass accretion rates and Specific Mass Accretion rates derived in this work across six different redshift bins.Table 4. Fitted Parameters for Power Law and Power Exponential Law Models applied to Pair Fractions, Merger Rates, and Specific Mass Accretion Rates, as described in Equation 19 and Equation 20. These fittings incorporate results from our high-redshift JWST data, supplemented by zero-point data from Casteels et al. (2014).', '5.1.2 Pair Fractions with Varying Stellar Mass Criteria': 'In this section, we investigate the evolution of pair fractions across different stellar mass bins, while maintaining consistency in all other methodologies. As outlined in Section 4.1.3, we initially apply a stellar mass selection of log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0 to create a general mass completed sample, before proceeding to correct for completeness using our novel methodology discussed in Section 4.3. The range log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0 is sufficiently broad to encompass all mass-completed samples, but it is also intriguing to investigate how pair fractions vary within this range. Therefore, in addition to computing the pair fractions at this mass range with a fiducial projection separation of 20 -50 kpc (shown as the black dotted line), we also calculate f p for narrower mass ranges: log 10 ( M ∗ / M ⊙ ) = 8 . 0 -9 . 0 and log 10 ( M ∗ / M ⊙ ) = 9 . 0 -10 . 0.  \nWepresent this comparison in Figure 5b, with f p computed for stellar mass ranges of log 10 ( M ∗ / M ⊙ ) = 8 . 0 -9 . 0 and log 10 ( M ∗ / M ⊙ ) = 9 . 0 -10 . 0 represented in red and blue diamonds. We observe that f p values for the log 10 ( M ∗ / M ⊙ ) = 8 . 0 -9 . 0 and log 10 ( M ∗ / M ⊙ ) = 8 . 0 -10 . 0 ranges are quite similar. However, the f p from the log 10 ( M ∗ / M ⊙ ) = 9 . 0 -10 . 0 range shows a markedly sharp declining trend with redshift. This distinction is mainly attributable to the stellar mass distribution of high-redshift galaxies, as outlined in Figure 3 and other recent literature (e.g., Duan et al. 2024; Harvey et al. 2024). High-redshift galaxies are generally less massive than low-redshift ones, and those within the log 10 ( M ∗ / M ⊙ ) = 9 . 0 -10 . 0 range are particularly rare, resulting in statistically insufficient sample sizes, thereby leading to distinct trends in f p .', '5.1.3 Simulation Pair fractions': "We calculate major merger pair fractions ( 𝜇 > 1 / 4) from the SCSAM simulations' deep dataset, employing the method described in Section 4.5. We consider various upper limits for rest-frame velocity differences between galaxies, with the maximum set at Δ 𝑣 = 2000 km/s. The evolution of major merger pair fractions with redshifts is illustrated in Figure 4b, represented by purple-yellow lines. Fitted  \npower-law and power-law + exponential curves based on observational data are depicted in black and red lines. Pair fractions from the simulation by Huško et al. (2022) are displayed in blue and pink lines, where blue represents log 10 ( M ∗ / M ⊙ ) > 9 . 0 and pink for log 10 ( M ∗ / M ⊙ ) > 8 . 0.  \nFrom SC-SAM simulation, we observe that major merger pair fractions reach a maximum value of approximately 0 . 3 at 𝑧 = 1 . 5 -3 . 5, before starting to decline at higher redshifts. Additionally, we find that the SC-SAM simulation inadequately describes the pair fractions at 𝑧 > 6, or at any redshift except near 𝑧 ∼ 4, due to the limited number of pairs seen in the simulation. A moderate agreement between the SC-SAM simulations and observational results is observed at 𝑧 < 4, where pair fractions can be adequately described by a powerlaw. However, at higher redshifts ( 𝑧 > 4), the observational results significantly deviate from the simulations. Our analysis reveals an increase in pair fractions up to 𝑧 = 8, reaching 0 . 211 ± 0 . 065, followed by a statistically flat evolution to 𝑧 = 11 . 5. This finding sharply contrasts with the fast decreasing trend predicted by the simulations after 𝑧 = 3, where the pair fractions decrease to around 0.01 at 𝑧 = [ 9 -11 . 5 ] .  \nThese are major differences between simulations and observations in the merger history, and this difference is now the most important conflict between simulations and observations for distant galaxies. There are several possibilities for explaining the difference we find, including observational as well as simulation input. For example, over-quenching of satellite galaxies may explain these distinctions. The Santa Cruz SAM, like many early models, excessively quenches satellite galaxies. In the current configuration, cosmological accretion halts once a galaxy becomes a satellite, rapidly depleting its gas supply and stifling star formation. This process is likely more extreme than what occurs in the real Universe, leading to simulated satellite galaxies that are less luminous and less massive than their observed counterparts. This would then imply that when these systems merge with most massive galaxies the merger is not seen as major, but minor. Updates to semi-analytic and empirical satellite models are necessary, and many researchers, such as Behroozi et al. and Pandya et al., are actively developing these enhancements.  \nIn addition to the SC-SAM simulation, we compare our results with major merger pair fractions derived by Huško et al. (2022), who employed the Planck Millennium cosmological dark matter simulation and the GALFORM semi-analytical model of galaxy formation. They computed pair fractions using a projection separation of 𝑟 𝑝 < 20 kpc and a velocity difference of Δ 𝑣 < 500 km/s. Their results are illustrated in Figure 4b with blue and pink lines, where blue lines represent f p for galaxies with log 10 ( M ∗ / M ⊙ ) > 8 . 0, and pink lines represent f p for galaxies with log 10 ( M ∗ / M ⊙ ) > 9 . 0. Comparing with SC-SAM f p , Huško et al. (2022) seems to agree better with our observational data points at 𝑧 > 6 . 5. However, at 𝑧 < 6 . 5, the f p from Huško et al. (2022) increases more steeply than observed in our data.  \nOverall, regardless of the simulations we compare with, we do not find strong agreement between our observations and simulations, as has been observed in other studies at lower redshift (e.g., Bertone & Conselice 2009; Jogee et al. 2009; Duncan et al. 2019). To explain this difference in models and observations, Man et al. (2016) proposes that the discrepancy arises from the different ways in which galaxies are selected. Observational studies typically rely on stellar mass, while simulations often use baryonic mass or dark matter mass. Since gas fractions in high-redshift galaxies are more abundant, some of the minor mergers ( 𝜇 < 1 / 4) selected based on stellar mass may actually be major mergers if we use baryonic mass. As this is the first observational work employing the close-pair method to study pair fractions at 𝑧 > 6 . 0, we lack information on the gas content of these systes to make this comparison. However, in the future it would be beneficial to include this aspect in the merger fractions for these galaxies.", '5.2 Merger Rates': "The merger rate is defined as the number of merger events per galaxy per unit time, denoted as R M ( 𝑧, M ∗ ) . Compared to pair fractions, merger rates are more fundamental in nature because the pair fraction is an observational metric, subject to variations due to differing methodologies across redshift bins. Such variations affect the sensitivity to the timescales being observed, thereby complicating direct comparisons across different studies. In contrast, comparing merger rates between literature becomes meaningful when the typical timescale for observing a galaxy merger is well-defined. The methodology for calculating R M ( 𝑧, M ∗ ) in this paper follows Duncan et al. (2019); Ferreira et al. (2020); Conselice et al. (2022), with specific details provided in following paragraphs.  \nThe merger rate per galaxy is expressed as  \nR M ( 𝑧, M ∗ ) = f p ( 𝑧, M ∗ ) × 𝐶 merg 𝜏 m ( 𝑧 ) , (21)  \nwhere f p ( 𝑧, M ∗ ) is the pair fraction at each redshift bin, 𝐶 merg is the fraction of pairs that will eventually result in a merger event, and 𝜏 m ( 𝑧 ) is the merger timescale at a certain redshift.  \nIn the literature, a value of 0 . 6 is commonly assigned to 𝐶 merg (Lotz et al. 2011; Conselice 2014; Mundy et al. 2017). However, this value is subject to significant uncertainty, as it depends on various factors including stellar mass, mass ratio, and redshift. Considering the high redshift scenario of our study, which further adds to the complexity of accurately defining 𝐶 merg , we opt to set this value to 1. This implies, under our working assumptions, that each pair we observe is anticipated to eventually merge.  \nThe merger timescale, 𝜏 m ( 𝑧 ) , refers to the duration over which a merging system remains in a specific observed physical state. Various methods have been developed to derive merger timescales, ranging from empirical measurements (Conselice 2009a) to simulations  \n(Conselice 2006; Kitzbichler & White 2008; Lotz et al. 2008b; Genel et al. 2014a; Vogelsberger et al. 2014a). Typical values haven been found to be around 𝜏 m ( 𝑧 ) = 0 . 3 -0 . 7 Gyr for pairs at 𝑧 < 3. Using IllustrisTNG simulation (Genel et al. 2014b; Vogelsberger et al. 2014b), Snyder et al. (2017); Duncan et al. (2019); Conselice et al. (2022) found this timescale evolve with redshifts and mass ratio 𝜇 . They determined that merger timescales scale as follows:  \n𝜏 m ( 𝑧 ) = 𝑎 × ( 1 + 𝑧 ) 𝑏 [Gyr] , (22)  \nwhere  \n        𝑎 = -0 . 65 ± 0 . 08 × 𝜇 + 2 . 06 ± 0 . 01 𝑏 = -1 . 60 ± 0 . 01 . (23)  \nIn Figure 6, we plot the evolution of merger timescales with redshifts and 𝜇 .  \nWe present major-merger rates in Figure 7 using this measure of the timescales. Merger rates derived in this work are represented by yellow-orange stars, with values for each redshift bin shown in Table 3. We include lower redshift results from Duncan et al. (2019); Conselice et al. (2022) in the plot and compute the merger rates from the pair fractions of Casteels et al. (2014) and López-Sanjuan et al. (2015) using our method. Combined with these previously published values from the literature, we observe that R M increases between 𝑧 = 0 . 0 -6 . 0 and stabilizes after 𝑧 = 6 . 0, with a mean number of merger events per galaxy of 5 . 78 ± 0 . 98 Gyr -1 . From 𝑧 = 8 -11 . 5, although the pair fraction decreases, the merger timescale also decreases by approximately the same ratio. According to Equation 21, this results in a stable R M within this redshift range.  \nPrevious studies have demonstrated a continuous power-law increase in R M from 𝑧 = 0 . 35 to 𝑧 = 6 . 0. In this study, we model the R M evolution using both a power-law and a power-law + exponential (as presented in Equation 19 and Equation 20). The fitting is performed using only our JWST data, supplemented by a zero-point data from Casteels et al. (2014), represented as a blue triangle. The fitted lines are shown as black and red lines in Figure 7, and the fitted parameters are detailed in Table 4. The discovery of saturation in merger rates at higher redshift suggests that the pattern is more accurately described by a power-law + exponential rather than a power law. The rapid merger rates observed during the Universe's first 1 Gyr indicate a rapid transformation of galaxy properties, which we will explore further in the remainder of this paper.", '5.3 Merger-Driven Stellar Mass Growth': 'In this section we elucidate the impact of merger activities on the stellar mass of galaxies, focusing specifically on the mass increase attributable to major mergers. This is important as it for the first time allows us to determine the role of mergers in galaxy formation - that is, the establishment of stellar mass within galaxies at 𝑧 > 6. We quantify this by the mass accretion rate (MAR), total mass accretion rate density (tMARD), and specific mass accretion rate (sMAR). These three quantities are examined and analyzed in this section.', '5.3.1 Mass Accretion Rates': "Themassaccretionrate, 𝜌 1 / 4 ( 𝑧 ) , is the average amount of mass added through mergers per unit time per galaxy. The subscript 1 / 4 signifies major mergers characterized by a mass ratio 𝜇 > 1 / 4. Recall that the primary galaxy is the one with the higher mass in a pair, while the secondary galaxy is the one with the lower mass. For each primary  \n(a) The evolution of our measured pair fractions with redshift. The new high-redshift major merger fraction (f p ) from our JWST analysis is depicted by the yellow-orange stars. The mass range used for our JWST data is log 10 ( M ∗/ M ⊙) = 8 . 0 -10 . 0. Additionally, we present lower redshifts results from four studies: green pentagons denote data from ALHAMBRA as reported by López-Sanjuan et al. (2015) for 𝑀 𝐵 < -20 . 0; rings represent CANDELS data by Duncan et al. (2019), with colored rings indicating each sub-field and black rings representing the combined results for masses log 10 ( M ∗/ M ⊙) > 10 . 3; diamonds symbolize data from the REFINE survey by Conselice et al. (2022) for masses log 10 ( M ∗/ M ⊙) > 11 . 0; and the blue triangle denotes GAMA data computed by Casteels et al. (2014). All literature except Conselice et al. (2022) present major merger f p ; Conselice et al. (2022) presents total f p . Power-law and power-law + exponential fits (Equation 19, Equation 20) are applied to our JWST data, with (Casteels et al. 2014) data used as a zero point. The fitted parameters are detailed in Table 4. We find an increase in pair fractions up to 𝑧 ∼ 8, reaching 0 . 211 ± 0 . 065, followed by a statistically flat evolution to 𝑧 = 11 . 5. The average f p between 𝑧 = [ 4 . 5 -11 . 5 ] is 0 . 17 ± 0 . 05.  \n<!-- image -->  \n(b) Evolution of pair fractions with redshift derived from the SC-SAM simulations and Huško et al. (2022). The fitted power and power exponential lines from observational data in Figure 4a are shown in black and red lines. Results from Huško et al. (2022) are displayed in blue and pink lines, where blue represents f p for galaxies with log 10 ( M ∗/ M ⊙) > 8 . 0, and pink for log 10 ( M ∗/ M ⊙) > 9 . 0. SC-SAM f p for the major merger ( 𝜇 > 1 / 4) is depicted by purple-orange lines, respectively. The color represents different upper rest-frame velocity limits between two galaxies in km/s, with each velocity maximum corresponding to a unique pair fraction curve. Major merger pair fractions from SC-SAM peak at 𝑧 = 2 -3 before declining. Observational data shows a gradual increase in the major merger pair fraction from 𝑧 = 4 . 5 to 𝑧 = 8 with a f p value of 0 . 211 ± 0 . 065, followed by a moderate decline to 0 . 117 ± 0 . 084 at 𝑧 = 10 . 5.  \n<!-- image -->  \nFigure 4. Two figures illustrate the evolution of pair fractions, f p , derived from observations (top panel) and the Santa Cruz-SAM simulations (bottom panel). A notable discrepancy in pair fractions between observations and simulations is observed at most redshifts.  \n(a) The pair fractions from 𝑧 = 4 . 5 to 𝑧 = 11 . 5, calculated using different projection separation criteria while maintaining consistent methodologies and a stellar mass range of log 10 ( M ∗/ M ⊙) = 8 . 0 -10 . 0. Results from our fiducial setting, [20,50] kpc, is shown by the black dotted line, while f p computed at [5,10], [5,20], [5,30], [5,40], and [5,50] kpc are shown by different colored circles. For every 10 kpc reduction in the projection separation search, f p decreases by a factor of 1 . 72 ± 0 . 47. Note that if we use a 5kpc inner limit we obtain very similar results as our default 20 kpc inner limit.  \n<!-- image -->  \nd  \nFigure 5. In these two figures, we present the pair fractions computed with varying (a) Projection Separation Criteria, and (b) Stellar Mass Selection Criteria, with all other methodologies remaining consistent.  \n<!-- image -->  \nRedshift  \n(b) Pair fractions from 𝑧 = 4 . 5 to 𝑧 = 11 . 5, calculated using different stellar mass selection criteria while maintaining consistent methodologies and a projection separation range of 20 -50kpc. Results from our fiducial setting, log 10 ( M ∗/ M ⊙) = 8 . 0 -10 . 0, are shown by the black dotted line, while f p computed at log 10 ( M ∗/ M ⊙) = 8 . 0 -9 . 0 and log 10 ( M ∗/ M ⊙) = 9 . 0 -10 . 0 are shown as red and blue diamonds. Pair fractions computed at log 10 ( M ∗/ M ⊙) = 8 . 0 -9 . 0 is very similar to our fiducial model. Pair fractions from log 10 ( M ∗/ M ⊙) = 9 . 0 -10 . 0 shows a sharp decline compared to the other two, which is because there are very few massive galaxies in this range at high redshift.  \ngalaxy within a given redshift bin, the additional stellar masses added through mergers is the average stellar mass of secondary galaxies, symbolized as M ∗ , 2 ( 𝑧 ) . Consequently, the mass accretion rate is expressed as:  \n𝜌 1 / 4 ( 𝑧 ) = R M ( 𝑧 ) × M ∗ , 2 ( 𝑧 ) , (24)  \nwhere M ∗ , 2 ( 𝑧 ) is calculated using the following integral:  \nM ∗ , 2 ( 𝑧 ) = ∫ M∗ , 1 𝜇 M∗ , 1 𝜙 ( 𝑧, M ∗ ) M ∗ 𝑑 M ∗ ∫ M∗ , 1 𝜇 M∗ , 1 𝜙 ( 𝑧, M ∗ ) 𝑑 M ∗ . (25)  \nHere, 𝜙 ( 𝑧, M ∗ ) is the Galaxy Stellar Mass Function (GSMF) at each redshift bin. We use stellar mass functions from Duncan et al. (2014)  \nFigure 6. Evolution of Merger Timescales as a Function of Redshift as well as mass ratio 𝜇 . The depicted data are derived using parameters reported by Conselice et al. (2022), based on analyses conducted with the IllustrisTNG 300-1 cosmological simulation. The numbers on the vertical lines correspond to the timescales (y-value) at corresponding redshifts, expressed in units of Gyr.  \n<!-- image -->  \nfor 𝑧 = 4 . 5 -6 . 5, and Harvey et al. (2024) for 𝑧 = 6 . 5 -11 . 5. The value M ∗ , 1 represents the average stellar mass of the primary galaxy in each redshift bin, computed as:  \nM ∗ , 1 ( 𝑧 ) = ∫ M max ∗ , 1 M min ∗ , 1 𝜙 ( 𝑧, M ∗ ) M ∗ 𝑑 M ∗ ∫ M max ∗ , 1 M min ∗ , 1 𝜙 ( 𝑧, M ∗ ) 𝑑 M ∗ , (26)  \nwhere M min ∗ , 1 and M max ∗ , 1 represent the stellar mass range we used, being 10 8 and 10 10 M ⊙ , respectively.  \nBy integrating the mass accretion rate over a specific time interval, we can retrieve the average amount of mass added to the primary galaxy's stellar mass during that time interval, denoted as  \n𝛿 M( 𝑡 ) = ∫ 𝑡 max 𝑡 min 𝜌 1 / 4 ( 𝑡 ) 𝑑𝑡. (27)  \nWe use the point at 𝑧 = 10 . 5 as the 𝑡 min value, corresponding to 441 Myr after the Big Bang. The redshift of each subsequent data point is used as 𝑡 max .  \nEvolution of 𝜌 1 / 4 ( 𝑧 ) and 𝛿 M( 𝑡 ) with redshifts are shown in Figure 8a and Figure 8b, respectively. Values of 𝜌 1 / 4 ( 𝑧 ) at each redshift bin are presented in Table 3. From Figure 8a, we observe that the galaxy mass accretion rates are generally higher at higher redshifts, but becomes stable after 𝑧 > 6 . 0. From 𝑧 = 4 . 5 to 𝑧 = 11 . 5, the average mass of secondary galaxies ( M ∗ , 2 ) decreases, but only by 2.5%. Thus, we conclude that the higher mass accretion rates at 𝑧 > 6 . 0 are primarily due to the observed higher merger rates.  \nIn Figure 8b, we present the increase in stellar mass resulting from major merger events. The value of each point indicates the integration of 𝜌 1 / 4 ( 𝑧 ) from 𝑧 = 10 . 5 to the redshift at that point. Consequently, this value is higher at lower redshifts, owing to a longer time interval allowing more mass to accumulate. On top of each point is the value of 𝛿 M( 𝑡 ) divided by the average primary galaxy stellar mass, M ∗ , 1 , in the corresponding redshift bin. We then fit the unlogged 𝛿 M( 𝑧 ) versus redshift using a linear line, depicted as the blue line in  \nthe plot. This linear relationship describes the trend very well, with uncertainties in the gradient and y-intercept being 4.24% and 3.12%, respectively.  \nExamining this figure in detail, we find that, on average, galaxy stellar masses have nearly tripled (2 . 77 ± 0 . 99) due to merger events from 𝑧 = 10 . 5 to 𝑧 = 5 . 0, over a period of ∼ 730 Myr. This trend is similar to findings at redshifts ( 𝑧 < 3), which suggest a doubling of stellar mass from mergers, according to Conselice et al. (2022). Therefore, we conclude merger events in the early Universe play a non-negligible role in growing galaxy stellar masses. One major question however is how important this process of stellar mass assembly is compared to the process of star formation, which we will investigate in next two sections.", '5.3.2 Total Mass Accretion Rate Density': 'In this section we investigate the contribution of mergers to galaxy stellar mass comparing with galaxy intrinsic star formation rates. To achieve this, we compare the Total Mass Accretion Rate Density (tMARD) with Total Star formation Rate Density (tSFRD), and Specific Mass Accretion Rate (sMAR) with Specific Star Formation Rate (sSFR). The later will be discussed in the next section.  \nWe compute the tMARD(z) by dividing the mass accretion rate per galaxy, 𝜌 1 / 4 ( 𝑧 ) , by the comoving volume, 𝑉 comoving ( 𝑧 ) , and then multiplying by the total number of galaxies in the corresponding redshift bin, before applying the area normalization to account for the entire sky. The mathematical way to express this is  \ntMARD ( 𝑧 ) = 𝜌 1 / 4 ( 𝑧 ) × 𝑁 ( 𝑧 ) galaxies 𝑉 comoving ( 𝑧 ) ∗ 𝐴 Total 𝐴 Survey , (28)  \nwith unit of M ⊙ yr -1 Mpc -3 . Note that this is a very narrow mass range for which the mass accretion is being measured, which we later compare with the star formation at the same mass range. We are not attempting here to obtain a \'total\' tMARD value over all stellar masses, which is often done to the star formation rate density measured as a function of redshift (e.g., Harikane et al. 2023; Adams et al. 2024).  \nThe tSFRD for each redshift bin is calculated by summing the SFRD of all galaxies within that bin, and then multiplying by the same area normalization factor. Note that we multiply the SFR by ( 1 -𝑅 ) , where 𝑅 is the return fraction that we calculated as 𝑅 = 0 . 423 from Madau & Dickinson (2014a), to account for the fraction of stellar mass that is returned to the interstellar medium through stellar processes such as supernovae and stellar winds. This adjustment ensures that we reflect only the mass that remains (the observed mass). Thus, the tSFRD at each redshift bin in this work can be written as  \ntSFRD ( z ) = " 𝑁 𝑖 SFR 𝑖 × ( 1 -𝑅 ) 𝑉 comoving ( 𝑧 ) × 𝐴 Total 𝐴 Survey . (29)  \nWe derive SFRs from three sources: Bagpipes , using both 10 Myr and 100 Myr averaged timescales based on a log-normal star formation history, and the SFRD inferred from the UV luminosity. As we directly utilize all galaxies from these eight JWST fields without any completeness corrections for galaxies fainter than the observational limit, the tSFRD reported here is generally significantly lower than tSFRD derived from UV luminosity functions (e.g., Harikane et al. 2023; Adams et al. 2024).  \nWe present the tMARD and tSFRD at 𝑧 = 4 . 5 -11 . 5 in Figure 9a. tMARD is represented by yellow-orange stars, while tSFRD is depicted with lines: red for Bagpipes 10Myraveragedtimescale, green  \nFigure 7. Evolution of the galaxy merger rate with redshift. High-redshift major merger rates from our JWST results are represented by yellow-orange stars. The mass range used for our JWST data is log 10 ( M ∗/ M ⊙) = 8 . 0 -10 . 0. Additionally, we present lower redshifts results from four studies: green pentagons denote data from ALHAMBRA as reported by López-Sanjuan et al. (2015) for 𝑀 𝐵 < -20 . 0; rings represent CANDELS data by Duncan et al. (2019), with colored rings indicating each sub-field and black rings representing the combined results for masses log 10 ( M ∗/ M ⊙) > 10 . 3; diamonds symbolize data from the REFINE survey by Conselice et al. (2022) for masses log 10 ( M ∗/ M ⊙) > 11 . 0; and the blue triangle denotes GAMA data computed by Casteels et al. (2014). All literature except Conselice et al. (2022) present major merger rate; Conselice et al. (2022) presents total merger rate. A power law and power-law + exponential models are fitted to the combined data, with parameters shown in Table 4. We observe a saturation in merger rates at 𝑧 > 6 . 0. Beyond this point, the merger rate becomes stable with a mean value of 5 . 78 ± 0 . 98 Gyr -1 .  \n<!-- image -->  \nfor 100 Myr averaged timescale, and purple for the UV luminosity based measurements. We find that the contribution of our measured major mergers to galaxy stellar mass assembly over this epoch is, on average, 45 (± 16 ) %, 69 (± 23 ) %, and 108 (± 44 ) %, equivalent to intrinsic star formation computed from Bagpipes 10 Myr averaged, UV, and Bagpipes 100 Myr averaged. Thus, a typical value is that major mergers are adding 74 (± 28 ) %tothe stellar masses of galaxies through this redshift range.', '5.3.3 Specific Mass Accretion rate': 'The specific mass accretion rate (sMAR) is defined as the mass accretion rate divided by the average stellar mass of primary galaxies:  \nsMAR ( 𝑧 ) = / M M ∗ = 𝜌 1 / 4 ( 𝑧 ) M ∗ , 1 ( 𝑧 ) = R M ( 𝑧 ) × M ∗ , 2 ( 𝑧 ) M ∗ , 1 ( 𝑧 ) . (30)  \nSimilar to the previous section, we compare the specific mass accretion rate (sMAR) from mergers with the specific star formation rate (sSFR) from intrinsic star formation in this section, to evaluate their respective contributions to the total galaxy stellar mass. The sSFR is computed from three sources: Bagpipes using 10 Myr and 100 Myr averaged, and UV luminosity. Note that we also multiply the SFR by ( 1 -𝑅 ) to reflect the observed SFR, as we did in the previous section.  \nWe display these two quantities in Figure 9b. The shaded red and green regions represent the sSFR Bagpipes averaged over 10 and  \n100 Myr respectively, and the purple region gives the sSFR UV , estimated directly from the rest UV photometry. The sMAR values from this study, alongside those derived using our methodology from literature-based pair fractions, are represented using distinct markers. Values of our JWST sMAR at each redshift bin are shown in Table 3. Wefitthe evolution of sMAR with a power-law and a power-law + exponential form, with the best-fitted parameters presented in Table 4. Discrepancies in sSFR between Bagpipes and the UV method arise from differences in star formation history assumptions. For the sSFR calculated with Bagpipes , we utilize a log-normal SFH, whereas for the UV-derived sSFR, we adopt the constant SFH model by Madau & Dickinson (2014b), incorporating the dust correction factor from Meurer et al. (1999).  \nFrom this comparison, we find results similar to we discussed in Section 5.3.2. The contribution of mergers to galaxy stellar mass is, on average, 45 (± 14 ) %, 51 (± 16 ) %, and 110 (± 34 ) %of the intrinsic star formation rates computed from Bagpipes averaged over 100 Myr, UV luminosity, and Bagpipes averaged over 10 Myr, respectively. In a broader context, combined with the findings from Section 5.3.2, we conclude that the contribution of mergers to galaxy stellar mass is 71 (± 25 ) %equivalent to the contribution from intrinsic star formation from gas. This implies that at least 42% ± 24% of the total stellar mass of galaxies arises from merger events. The reason for saying "at least" is that we are searching for close pairs, and some of the galaxies in these pairs may have already merged before we observe them (post-mergers). Identifying post-merger galaxies requires a more complex analysis, such as using deep learning con-  \n(a) Average Mass Accretion Rate per Galaxy vs. Redshift. The value on the yaxis is the amount of mass accretion onto an \'average\' primary galaxy within each redshift bin due to major mergers. We find that this value increases up to 𝑧 = 6, and then becomes stable. Values at each redshift are listedin Table 3.  \n<!-- image -->  \n(b) The integrated average stellar mass accreted per galaxy, attributed exclusively to major mergers. This figure is the integral of the top panel, such that the average total amount of accumulated stellar mass from mergers starting at 𝑧 = 10 . 5 is shown in the bottom panel. The text above each data point indicates the ratio of the total mass accreted due to mergers compared to the average stellar mass of the primary galaxies at the start of the integration. The fit shown with the blue line, applied to the unlogged mass scale, effectively captures this trend.  \n<!-- image -->  \nFigure 8. Two figures demonstrating: (a) Stellar mass accretion rates from major mergers, (b) Integrated stellar mass accreted from major mergers.  \nvolutional neural networks (Ferreira et al. 2020), which is beyond the scope of this paper. Thus, the intrinsic star formation rates from gas that we computed using SED and UV luminosity inevitably include contributions from mergers as well.', '6 DISCUSSION': "Utilizing the deepest and widest observational NIRCam imaging programs to date from JWST, we have conducted a new close-pair major merger analysis spanning redshifts from 𝑧 > 4 . 5 to 𝑧 = 11 . 5. This is a redshift range where the merger history is largely unknown observationally. Based on this imaging data and using new methods we have computed pair fractions, merger rates, and stellar mass accretion rates at 𝑧 > 6, highlighting the significant role of major mergers in galaxy formation.  \nNotably, we observed a roughly constant rapid merger rate at the highest redshifts of 5 . 78 ± 0 . 98 Gyr -1 mergers per galaxy, which is a rapid merger rate several orders of magnitude higher than the merger rate locally. We estimate that an average galaxy from 𝑧 = 10 . 5 to 𝑧 = 5 . 0 would become approximately 2 . 77 ± 0 . 99 times more massive by 𝑧 = 5 . 0. This is roughly similar to the amount of stellar mass added to these galaxies due to the star formation process. This also shows that mergers are a critical factor in galaxy formation early in the Universe.  \nAsaresult, mergers are most certainly a significant factor in galaxy evolution at early times and far more important than even at modest redshifts at 3 < 𝑧 < 6 (e.g., Duncan et al. 2019). In this paper, we undertake the first assessment of the contribution of mergers to stellar mass assembly across the redshift range 𝑧 = 4 . 5 -11 . 0 (see Section 5.3, Figure 8a, Figure 9a, and Figure 9b). We determine that mergers contribute 71 ± 25% of the mass formed through star formation, which implies that at least 42% ± 24% of an average galaxy's total stellar mass arises from mergers during this redshift range.  \nAforthcoming detailed analysis by Duan et al. (2024c, in preparation) will further explore how mergers influence star formation, stellar mass growth, and AGN activity. Additionally, our calculated pair fractions and merger rates are crucial for predicting detectable gravitational waves at 𝑧 > 6, providing essential forecasts for the LISA (Laser Interferometer Space Antenna), which is slated for launch in 2035.  \nDespite using nearly all available public JWST NIRCam Cycle1 programs, our combined masked area of 189.36 arcmin 2 covers only a minuscule portion of the entire sky, and our results could be skewed by cosmic variance. However, considering our field locations in both the Northern and Southern hemispheres and their random selection, the impact of cosmic variance is likely negligible, as discussed as well in the luminosity function based on these data (Adams et al. 2024). Our study's accuracy is primarily challenged by three factors: merger timescales, sample completeness, and photometric redshift uncertainties. Merger timescales are challenging to measure and have historically been assumed to be constant. In this study, we utilize merger timescale estimates from the IllustrisTNG simulation, which suggests a redshift evolution characterized by shorter timescales at higher redshifts, although these estimates are associated with significant uncertainties. In fact, it is the time-scale of these mergers which is the most uncertain part of our analysis or any analysis of galaxy merger histories, as there are few observational ways to determine this time-scale (c.f. Conselice 2009b).  \nFurthermore, while we utilize the JAGUAR simulation to compute sample completeness, and have developed an innovative method to correct for incompleteness, achieving an entirely complete model is unfeasible. Moreover, the absence of spectroscopic observations for high-redshift photometric samples necessitates reliance on photometric redshifts to conduct our merger analysis, which introduces potential uncertainties similar to those noted in lower redshift studies (López-Sanjuan et al. 2015; Mundy et al. 2017; Conselice et al. 2022). One issue is that the Ly 𝛼 damping wing, which could po-  \n(a) In this figure, we present the Total Mass Accretion Rate Density (tMARD) from major mergers, in units of M ⊙ yr -1 Mpc -3 , versus Redshift alongside the Total Star Formation Rate Density (tSFRD) versus Redshift. The values of tMARD are represented by yellow-orange stars, while we present three sources of tSFRD: red for Bagpipes 10 Myr averaged timescales (log-normal SFH), green for 100 Myr averaged timescales, and purple for the UV luminosity measurement. From this comparison, we find that the contribution of mergers to galaxy stellar mass assembly, is on average, is 45 (± 16 ) %, 69 (± 23 ) %, and 108 (± 44 ) %, equivalent to intrinsic star formation computed from Bagpipes 10 Myr averaged, UV, and Bagpipes 100 Myr averaged.  \n<!-- image -->  \nFigure 9. These two figures compare the contributions of stellar mass assembly from major mergers with the amount of mass added due to intrinsic star formation: (a) Shows the Total Mass Accretion Rate Density (tMARD) and the Total Star Formation Rate Density (tSFRD) across redshifts, and (b) Shows the Specific Mass Accretion Rates (sMAR) and Specific Star Formation Rates (sSFR) across redshifts.  \n<!-- image -->  \nRedshift  \n(b) Evolution of specific mass accretion rates (sMAR) from mergers and specific star formation rates (sSFR) with redshift. We present three sources of sSFR: Bagpipes with a 100 Myr average and 10 Myr average timescale, and UV luminosity, which are depicted by green, red, and purple shaded regions. These sSFR values are computed from the entire EPOCHS V1 high-redshift sample. High-redshift specific mass accretion rates derived in this work are represented by yellow-orange stars. Additionally, we use our methodology to derive the values of the sMAR from pair fractions from four lower redshifts results: green pentagons denote data from ALHAMBRA as reported by López-Sanjuan et al. (2015) for 𝑀 𝐵 < -20 . 0; rings represent CANDELS data by Duncan et al. (2019), with colored rings indicating each sub-field and black rings representing the combined results for masses log 10 ( M ∗/ M ⊙) > 10 . 3; diamonds symbolize data from the REFINE survey by Conselice et al. (2022) for masses log 10 ( M ∗/ M ⊙) > 11 . 0; and the blue triangle denotes GAMA data computed by Casteels et al. (2014). All literature except Conselice et al. (2022) present major merger result; Conselice et al. (2022) presents total merger result. From this comparison, we find the contribution of mergers to galaxy stellar mass is, on average, 45 (± 14 ) %, 51 (± 16 ) %, and 110 (± 34 ) % equivalent to the intrinsic star formation rates computed from Bagpipes averaged over 100 Myr, UV luminosity, and Bagpipes averaged over 10 Myr, respectively.  \ntentially affect photometric redshifts due to Ly 𝛼 photons passing through varying amounts of neutral IGM, which is also not taken into account. This omission may lead to under-predicting the number of close pairs, an effect that worsens at higher redshifts when there is more neutral IGM to traverse and at greater separations when the sightlines become more independent. This effect cannot be averaged out, but it does not play a role when the sightlines are exactly the same and is negligible when the sightlines are similar (as in close-pair systems). However, we have minimized these three factors (merger timescales, sample completeness, and photometric redshift uncertainties) using a statistical probabilistic pair counting methodology as in the previous literature.  \nCompared with a recent morphological merger study at 𝑧 = 4 -9 by Dalmasso et al. (2024), which reported pair fractions of f p = 0 . 11 ± 0 . 04 with no observed evolution with redshift, our findings show a higher average f p of 0 . 17 ± 0 . 05 in this redshift range. We find an increase in pair fractions up to 𝑧 = 8, reaching 0 . 211 ± 0 . 065. The fitted power-law + exponential model presents a flat trend at 𝑧 > 6. Our f p values at 𝑧 = 4 -9 consistently exceeds those reported in the morphological study, a discrepancy likely because either merger timescales are shorter, as we find, or that not all close-pair systems ultimately merge, resulting in lower morphology-based pair fractions. There needs to be a careful study of merging galaxies based on pairs and structure to resolve this issue, which sohuld inovlve the use of simulations.  \nAs discussed, we also compare our observational pair fractions with those predicted by simulations, specifically the SC-SAM and the Millennium cosmological dark matter simulation integrated with the GALFORM semi-analytical model. We find a significant disagreement between our observations and the simulations. Several reasons can explain this discrepancy. First it must be noted that these simulations are often calibrated on observed galaxies properties, such as the mass and luminosity function, but the merger history is not considered in this. Thus the merger history is an excellent test of the predictive power of simulations and models, and the disagreements should be considered in some detail.  \nSome of the possible explanations for the observed difference include: overly quenched satellite galaxies in the semi-analytical model, the different ways galaxies are selected between observational works and simulations (e.g., Bertone & Conselice 2009), and the use of stellar mass for observations versus baryonic mass for simulations, as discussed in Section 5.1.3, as well as in past comparison papers such as Duncan et al. (2019) and Bertone & Conselice (2009). More simulation work focused on the merger history is needed to reconcile observations with theory, as well as how to calibrate stellar and total mass as is done in simulations (Conselice et al. 2018).  \nOur results are revealing as galaxy formation occurs through a variety of mechanisms, with star formation from gas being largely the predominant process for building stellar mass. However, the literature indicates that the rate of galaxy mass growth often surpasses the contributions possible from the in situ gas within the galaxies themselves (Conselice 2012; Ownsworth et al. 2016; Padmanabhan &Loeb 2020; Walter et al. 2020). From our observations we are now able to trace this process back to the earliest times in which galaxies can be detected with JWST, and thus determine how galaxy assemble occurs rather than just when.", '7 CONCLUSIONS': "In this study, we conduct an early analysis of galaxy merger evolution up to 𝑧 = 12, utilizing the deepest and widest data from eight JWST  \nCycle-1 field programs, which collectively cover an unmasked area of 189.36 arcmin 2 . We have accurately measured major-merger pair fractions, merger rates, and stellar mass accretion rates for galaxies with stellar masses between 8 . 0 < log 10 ( M ∗ / M ⊙ ) < 10 . 0. Our methodology, while grounded in established literature, introduces a novel completeness correction approach to precisely measure these quantities in the epoch of reionization. Our key findings are summarized as follows:  \nI. We find an increase in the measured pair fractions up to 𝑧 ∼ 8, reaching a value of 0 . 211 ± 0 . 065, followed by a statistically flat evolution to 𝑧 = 11 . 5. We conclude that the power-law + exponential model, f p = ( 0 . 030 ± 0 . 016 )∗( 1 + 𝑧 ) 1 . 359 ± 0 . 725 ∗ 𝑒 (-0 . 138 ± 0 . 175 ) ( 1 + 𝑧 ) , is more appropriate than a traditional power law to describe the evolution of pair fractions with redshifts. This suggests that we have measured the peak merger fraction at 𝑧 ∼ 6 which stays constant at higher redshifts.  \nII. Using merger time-scales based on simulations we find that merger rates (number of mergers per unit time) increases from the local Universe to 𝑧 = 6 . 0, at which time they stabilize at a value of 5 . 78 ± 0 . 98 Gyr -1 . The rapid merger rates observed during the Universe's first 1 Gyr indicate a rapid transformation of galaxy stellar masses and thus also galaxy properties, showing that major mergers in the early Universe are of major importance.  \nIII. Based on our measured merger rates we calculate that an average galaxy during the epoch 𝑧 = 10 . 5 to 𝑧 = 5 . 0 grows by a factor of 2 . 77 ± 0 . 99 by 𝑧 = 5 . 0. It is important however to keep in mind that this is the average stellar mass addition, and that some galaxies will have more than this while others will have less.  \nIV. We analyze the stellar mass accretion into galaxies due to mergers and compare it with the intrinsic star formation process within galaxies themselves by examining the total mass accretion rate density (tMARD) from mergers versus the total star formation rate density (tSFRD) as well as the specific mass accretion rates (sMAR) from mergers with its star formation equivalent - the specific star formation rates (sSFR). We find that at 𝑧 = 4 . 5 -11 . 5, the stellar mass accreted from major mergers is 71 ± 25% of the stellar mass created from star formation.  \nBased on this approximately 42 ± 24% of a galaxy's total stellar mass accrues from major merger events at this time, while the remainder originates from star formation. This shows that galaxy mergers are at least as important as star formation for building the stellar masses of galaxies at early times. This does not include the star formation induced the mergers themselves, or minor mergers, which would only increase the importance of the merger process for adding stellar mass to galaxies.  \nV. We compare pair fractions from observational data with those derived from SC-SAM and Planck Millennium cosmological dark matter simulations. We find significant disagreements between our observations and simulations. Several reasons can explain this discrepancy, such as the overly quenched satellite galaxies in the semianalytical model, the different ways galaxies are selected between observational works and simulations, and the use of stellar mass for observations versus baryonic mass for simulations. More simulation work focused on the merger history is needed to reconcile observations with theory.  \nOverall, using an integrated probabilistic pair counting method, we have made the first attempt to extend the investigation of merger evolution beyond the 𝑧 = 6 limit to 𝑧 = 11 . 5. This study currently represents the highest redshift analysis of mergers and may continue to be one of the most advanced. A worthwhile future endeavor, when a substantial number of photometric and spectroscopic observations are available, would be to compare these results with more extensive  \ndatasets to trace the minor merger history as well as determine the properties of mergers and how they correlate with SF and AGN activity.", 'DATA AVAILABILITY': 'The data underlying this article is made available by Adams et al. (2023, 2024); Austin et al. (2023a); Austin et al. (2024); Harvey et al. (2024); Conselice et al. (2024). The catalogues of the sample discussed herein may be acquired by contacting the corresponding author.', 'ACKNOWLEDGEMENT': 'Weacknowledge support from the ERC Advanced Investigator Grant EPOCHS(788113),aswellasastudentship from STFC. LF acknowledges financial support from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) in the form of a PhD studentship. CCL acknowledges support from the Royal Society under grant RGF/EA/181016. CT acknowledges funding from the Science and Technology Facilities Council (STFC). We also extend our gratitude to Dr. Steven Willner for his constructive comments and suggestions on the manuscript. This work is based on observations made with the NASA/ESA Hubble Space Telescope (HST) and NASA/ESA/CSA James Webb Space Telescope (JWST) obtained from the Mikulski Archive for Space Telescopes ( MAST ) at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST, and NAS 5-26555 for HST. Some of the data products presented herein were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center, which is funded by the Danish National Research Foundation under grant No. 140. This research made use of the following Python libraries: N/u.pc/m.pc/p.pc/y.pc (Harris et al. 2020); S/c.pc/i.pc/p.pc/y.pc (Virtanen et al. 2020); M/a.pc/t.pc/p.pc/l.pc/o.pc/t.pc/l.pc/i.pc/b.pc (Hunter 2007); A/s.pc/t.pc/r.pc/o.pc/p.pc/y.pc (Astropy Collaboration et al. 2013, 2018, 2022); EAZY-PY (Brammer et al. 2008); L/e.pcP/h.pc/a.pc/r.pc/e.pc (Arnouts et al. 1999; Ilbert et al. 2006); B/a.pc/g.pc/p.pc/i.pc/p.pc/e.pc/s.pc (Carnall et al. 2018); /m.pc/p.pc/i.pc4/p.pc/y.pc (Dalcin & Fang 2021); P/i.pc/c.pc/k.pc/l.pc/e.pc (Van Rossum 2020).', 'REFERENCES': 'Rieke M., et al., 2023a, arXiv preprint arXiv:2306.02466  \nde Ravel L., et al., 2009, A&A, 498, 379', 'APPENDIX A: THE JWST COMPLETENESS LEVELS': "There are various incompletenesses that we have to consider within this merger work. Some of these issues are unique and differ from say incompletenesses that arise in detection of galaxies at various magnitudes (e.g., Harvey et al. 2024) that are needed to correct galaxy mass and luminosity functions (Adams et al. 2024). Namely, because we are incomplete in detections, we have to additionally consider how to deal with this in terms of galaxies which would be detected as pairs if the data was complete, but due to the depth of the imaging and the brightness of these galaxies remain undetected.  \nThe completeness level for our JWST high redshift galaxy sample at each redshift bin and mass range is computed utilizing the JAGUAR simulation (Williams et al. 2018). We consider this within each of our fields, which will have their own unique properties of incompleteness. The JAGUAR simulation produces an output catalog of simulated galaxies, each characterized by attributes such as redshift, stellar mass, photometry, etc. For a given field, we consider the known average depths (i.e., the average uncertainty in flux) and apply a Gaussian scatter to the galaxy's photometry accordingly. The fluxes, once scattered, are processed through our comprehensive EAZY SED fitting and selection pipeline to determine whether a galaxy is selected and if its redshift is accurately estimated. We apply the same Bagpipes fitting to obtain masses and compare these masses to known catalog masses. Applying this process to the entire JAGUAR catalog allows us to determine the fraction of true high redshift galaxies in our final sample (completeness) and the number of low-redshift interlopers. This varies depending on the field, as both the filters and depths differ across fields. We can parameterize the completeness in terms of known variables (stellar mass, M UV , apparent magnitude), by categorizing them into bins. A more detailed explanation is outlined in previous EPOCHS papers Adams et al. (2024); Austin et al. (2024); Harvey et al. (2024). We present the completeness level of all eight JWST fields we used in this work in Figure A1. In Figure A2 and Figure A3, we present the impact of incompleteness on the derived pair fractions using the lower redshift CANDELS catalog. This step is a crucial part of the completeness correction, which we have discussed in depth in Section 4.3.", 'APPENDIX B: MERGER CATALOG': "In this appendix, we present a catalog and cutout images of all merger close pairs detected with a probability ( N 𝑧 ) higher than 0.90 and physical projection separation 10 -50 kpc, across all eight JWST fields: CEERS, JADES GOODS-S, NEP-TDF, NGDEEP, GLASS, El-Gordo, SMACS-0723, and MACS-0416, in Table B1, Figure B1, and Figure B2. We also show the number of pairs selected with N 𝑧 greater than a certain threshold in Figure B3.  \nFor each close-pair cutout graph, a 0.32 arcsecond aperture is drawn in red on each galaxy, with the redshift labeled at the top left of each aperture. We employ the Nakajima template (Nakajima & Maiolino 2022) for photometric AGN detection, marking the identified AGN-host galaxies in orange. A notable example is the close pair system, JADES GOODS-S 32751 and 32927, which have been observed with NIRSpec as part of the JADES project and with the F444W grism by FRESCO (Oesch et al. 2023). This system provides compelling evidence of a dual AGN merger and is a strong Lyman𝛼 emitter. The detailed analysis of this system will be presented in a forthcoming paper by Li et al. (2024).  \nFigure A1. The level of completeness for our Eight JWST Fields computed using the JAGUAR Simulation. The grey areas indicate where no galaxies in those mass-redshift ranges within the simulation are located. As our analysis focuses on the mass range log 10 ( M ∗/ M ⊙) = 8 . 0 -10 . 0 and the redshift range 𝑧 = 4 . 5 -11 . 5, these gaps in the simulation data have a very minor impact on our results.  \n<!-- image -->  \nFigure A2. Evolution of the value f ratio p (Equation 11) across various levels of incompleteness in different redshift bins. f ratio p is the pair fraction at a given completeness level divide by the complete sample case (completeness = 1). To reduce random errors resulting from the arbitrary exclusion of samples, the procedure for introducing incompleteness into the pair fraction calculation is repeated five times, and the average f ratio p values are shown for each redshift bin. An observable trend shows a lower f ratio p corresponding to reduced completeness. Linear fits have been applied to each redshift bin, and the resulting gradients and y-intercepts of these best-fit lines are presented in Figure A3.  \n<!-- image -->  \nFigure A3. Evolution of the gradient and y-intercept for Best-Fit Lines Across Different Redshift Bins as shown in Figure A2. Straight lines are fitted to both the gradient and the y-intercept using the MCMC method, revealing an increasing trend in the gradient and a decreasing trend in the Y-intercept with redshift. The line of best fit for the gradient is 𝑦 = ( 0 . 028 ± 0 . 012 ) 𝑥 + ( 0 . 879 ± 0 . 044 ) , and for the Y-intercept, it is 𝑦 = (-0 . 013 ± 0 . 006 ) 𝑥 + ( 0 . 072 ± 0 . 021 ) .  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure B1. Cutout images for all close pair galaxies with N 𝑧 > 0 . 90 and Projection Separations within 10 -50 kpc.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure B2. Stacked images of all close-pair galaxies with N 𝑧 > 0 . 90 and projection separations ranging from 10 to 50 kpc. The images are generated by averaging data across multiple wavelength bands, with blue representing the ['F090W', 'F115W', 'F150W'] bands, green the ['F200W', 'F277W', 'F335M'], and red the ['F356W', 'F410M', 'F444W']. A 0.32 arcsecond aperture is drawn in red on each galaxy, with the redshift labeled at the top left of each aperture. Galaxies identified photometrically as AGN using the Nakajima template (Nakajima & Maiolino 2022) are highlighted in orange.  \nTable B1. Merger catalogs for galaxy pairs with N 𝑧 > 0 . 90 and projection separations 10 - 50 kpc.  \nFigure B3. Total Number of galaxy pairs detected by applying different N 𝑧 threshold in all eight JWST fields. The x-axis shows the threshold at which N 𝑧 is greater than a certain value, while the y-axis give the number of galaxies retrieved at that level. In general, the higher the value of N 𝑧 , the more likely that a galaxy pair is real.  \n<!-- image -->  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author."}