Characterising Tidal Features Around Galaxies in Cosmological Simulations

Tidal features provide signatures of recent mergers and offer a unique insight into the assembly history of galaxies. The Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) will enable an unprecedentedly large survey of tidal features around millions of galaxies. To decipher the contributions of mergers to galaxy evolution it will be necessary to compare the observed tidal features with theoretical predictions. Therefore, we use cosmological hydrodynamical simulations NewHorizon, EAGLE, IllustrisTNG, and Magneticum to produce LSST-like mock images of $z\sim0$ galaxies ($z\sim0.2$ for NewHorizon) with $M_{\scriptstyle\star,\text{ 30 pkpc}}\geq10^{9.5}$ M$_{\scriptstyle\odot}$. We perform a visual classification to identify tidal features and classify their morphology. We find broadly good agreement between the simulations regarding their overall tidal feature fractions: $f_{\text{NewHorizon}}=0.40\pm0.06$, $f_{\text{EAGLE}}=0.37\pm0.01$, $f_{\text{TNG}}=0.32\pm0.01$ and $f_{\text{Magneticum}}=0.32\pm0.01$, and their specific tidal feature fractions. Furthermore, we find excellent agreement regarding the trends of tidal feature fraction with stellar and halo mass. All simulations agree in predicting that the majority of central galaxies of groups and clusters exhibit at least one tidal feature, while the satellite members rarely show such features. This agreement suggests that gravity is the primary driver of the occurrence of visually-identifiable tidal features in cosmological simulations, rather than subgrid physics or hydrodynamics. All predictions can be verified directly with LSST observations.


INTRODUCTION
In the 'hierarchical structure formation' model of the Universe (e.g.Press & Schechter 1974;Fall & Efstathiou 1980;Ryden & Gunn 1987;van den Bosch 2002;Agertz et al. 2011), mergers play an important role in the evolution of galaxies, transforming galaxies through a variety of non-secular processes e.g.merger-driven star formation, merger-enhanced active galactic nuclei activity, dynamical evolution and stellar mass accretion (e.g.Dubois et al. 2016;Martin et al. 2018;Davison et al. 2020;Martin et al. 2021;Remus & Forbes 2022;Cannarozzo et al. 2023).However, the intricacies of this picture, regarding the detailed merger statistics of galaxies, remain unresolved.There remain untested theoretical predictions from cosmological simulations.For example, Naab & Ostriker (2009) and Martin et al. (2018) predicted that low-mass mergers are the major contributor to galaxy morphological evolution since  ∼ 1.There are ★ E-mail: aman.khalid@unsw.edu.au also detailed predictions regarding the environmental dependence of galaxy interaction and merger rates from cosmological-N-body simulations (e.g.Gnedin 2003;Mihos 2003) and semi-analytic models (Jian et al. 2012), where the galaxy encounter rates increase with increasing halo mass to a peak at halo masses corresponding to galaxy groups before dropping off at halo masses above that, corresponding to large galaxy clusters.To test these predictions we need to better understand galaxy mergers.
There are several ways to measure galaxy mergers observationally.These include detecting visible signatures of ongoing and past mergers via the presence of tidal features and selecting close pairs of galaxies that are likely to merge in the near future.Tidal features are diffuse non-uniform regions of stars that extend out from a galaxy, signatures of proceeding and concluded mergers in the forms of 'tails', 'streams', 'asymmetric halos', 'double nuclei' and 'shells' (examples in  The galaxy marked by a red circle in the bottom right image is an example of one that we removed from the sample.While there is enough stellar mass within the 30 pkpc aperture for the object to be in our sample, the object itself is too dense and compact to be an accurate representation of a galaxy of  ★,30 pkpcs ≥ 10 9.5 M ⊙ .Each image is 124 pkpc × 124 pkpc (∼ 2400 pixels × 2400 pixels).
making tidal features crucial probes for a galaxy's recent merger history.Tidal feature detection requires very deep images to pick up signatures of ongoing and past mergers (e.g.Atkinson et al. 2013;Kado-Fong et al. 2018;Martin et al. 2022).Close pairs are detected most accurately using spectroscopic observations to measure robust 3D positions for each galaxy (e.g.Robotham et al. 2014;Banks et al. 2021).Close pair detection is important in understanding the role mergers play in driving galaxy evolution, however, this method does have limitations.It cannot detect the presence of a secondary galaxy if: that galaxy has been ripped apart by tidal forces; absorbed into a host galaxy; or the secondary galaxy is too low mass to detect via spectroscopy (e.g.Lotz et al. 2011;Desmons et al. 2023b).Studying the tidal features of a galaxy provides insight into these aspects of the merging process.
The formation pathways of tidal features have been studied using N-body and hydrodynamic simulations to obtain insight into the merger histories of observed galaxies that host tidal features.These simulations have shown that tidal features are tracers of galaxy interactions and mergers.They probe detailed properties of these interactions (e.g.dynamics, mass ratios and orbits).They have established that tail-like structures (top left panel Figure 1) are formed from high angular momentum passages between similar mass galaxies (common in major mergers), while streams form from almost circularly infalling lower mass satellite galaxies (minor mergers; e.g.Toomre & Toomre 1972;Hendel & Johnston 2015;Karademir et al. 2019).
Various formation scenarios have been proposed to explain the formation of shells, the scenario most supported by analytical and numerical models is formation predominantly through radial mergers (e.g.Hendel & Johnston 2015;Amorisco 2015;Pop et al. 2018;Karademir et al. 2019;Valenzuela & Remus 2023), resulting in overdensities of stripped stars accumulating at the apocentres of their orbits.Simulations have linked the properties of shells to the merger dynamics and mass ratios of the progenitor galaxies.N-body studies such as Hernquist & Spergel (1992) have found that shells formed from major mergers tend to be preferentially aligned with the major axis of the host galaxy.Pop et al. (2018) found that most of the  = 0 shell population of the hydrodynamic-cosmological simulation Illustris were formed predominantly through major mergers, suggesting that due to dynamical friction, pairs of more massive galaxies could probe a greater range of impact parameters (deviating further from a purely radial infall) and still form shells when compared to more minor mergers.These simulations provide strong evidence for shell formation occurring through mergers with more radial trajectories than those that form tails and streams.
The understanding of tidal feature formation from simulations is used to interpret observational studies of the characteristics of tidal features such as their colours, dynamics, appearance and the environment in which the tidal feature host resides.These are then used to constrain the properties of the interaction that caused them.For example Foster et al. (2014) and Martínez-Delgado et al. (2021) used spectroscopic and photometric data, respectively, in conjunction with numerical simulations, to study the formation of tidal features around the NGC 4651 (Foster et al. 2014) and M104 (Martínez-Delgado et al. 2021) galaxies and constrain their recent merger histories.Foster et al. (2014) were able to infer the recent passage of a satellite galaxy through the disk of NGC 4651 from their modelling of an observed stellar stream.Martínez-Delgado et al. (2021) were able to estimate the time around which the major wet merger forming M104 occurred.These particular examples illustrate the strength of tidal features in reconstructing merger histories of observed galaxies and hint at the potential of using tidal features to analyse larger statistical populations of galaxies to provide a detailed understanding of the merger process in our Universe.
One of the largest observational studies of tidal features to date, Kado-Fong et al. (2018), analysed the tidal features in 0.05 <  < 0.45 galaxies in the HSC-SSP (Hyper Suprime-Cam Subaru Strategic Program) Wide layer.They found evidence suggesting that shells are formed predominantly through minor mergers.
Cosmological hydrodynamical simulations enable the study of large statistically-significant samples of galaxies to predict the dominant pathways of merger signatures in the ΛCDM model of hierarchical structure formation.Observations provide us with snapshots of our Universe at various lookback times, giving us the ground truth to which we can compare our models.Studying both observations and simulations enables us to understand the merger histories that are inferred from observations.
With the upcoming Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST; Ivezić et al. 2019;Robertson et al. 2019;Brough et al. 2020) it will be possible to study tidal features around millions of galaxies (Martin et al. 2022), allowing for the most robust statistical survey of tidal features to date.To fully capitalise on this unprecedented wealth of data we will need robust predictions of the properties of these features from current hydrodynamiccosmological simulations.In this study, we will characterise the tidal features around galaxies in the four hydrodynamic-cosmological simulations (described in §2.1), NewHorizon (Dubois et al. 2021), Il-lustrisTNG (Springel et al. 2018;Nelson et al. 2018;Marinacci et al. 2018;Naiman et al. 2018), EAGLE (Schaye et al. 2015;Crain et al. 2015) and Magneticum (Teklu et al. 2015), in an observationallymotivated approach that facilitates future comparison to LSST.We do this by visually classifying the tidal features present in LSST-like mock images, produced following Martin et al. (2022).We study four simulations to probe whether the different subgrid physics models applied to each simulation result in different merger pathways and whether this presents differently in the tidal features of galaxies.We compare the limited observations made to date to test how well the simulations are able to replicate them.
Our sample selection and production of mock images are described in §2.2 and §2.3, respectively.Our methodology for visual classification, including the tidal feature morphologies considered, are given in §2.4 and we describe the methodology for calculating the errors on our tidal feature fractions in §2.5.We present a statistical analysis of our visual classifications, considering the occurrence rates of tidal features in each simulation in §3 and then present how the fraction of galaxies exhibiting tidal features changes with both stellar mass (in §3.1) and halo mass (in §3.2).We then discuss our results with respect to past analyses of simulations and observational results in §4.1 and §4.2 respectively.The implications of our visual classifica-tion method are addressed and discussed in §4.3.In §4.4 we consider the implications of our results on the occurrence and observability of tidal features across the cosmological hydrodynamical simulations that we study.We discuss the relationship between tidal feature occurrence in the hydrodynamic-cosmological simulations we study and the environment in which the galaxies reside in §4.5.In §4.6 we explore the predictions for LSST regarding tidal feature fractions and their trends with stellar and halo mass based on the hydrodynamiccosmological simulations studied here.We draw our conclusion in §5.
We use the native cosmology from each simulation for calculating the distances between particles and creating our mock images, these are given in Table 1.For distances we use a 'c' prefix to denote comoving coordinates and a 'p' prefix to denote proper coordinates, e.g.cMpc is comoving megaparsecs and pMpc is proper megaparsecs.

Simulations
To investigate the characteristics of tidal features in cosmological hydrodynamical simulations we explored four different simulations.These simulations evolve gas, stars and dark matter under gravity and hydrodynamics along with a range of models for subgrid processes to simulate the realistic hierarchical structure formation of galaxies.
For our study we selected four simulations, NewHorizon, EAGLE RefL0100N1504, IllustrisTNG L75n1820TNG, Magneticum Pathfinder Box4-uhr (from here on, NewHorizon, EA-GLE, TNG and Magneticum).These four simulations, summarised in Table 1, allow us to probe tidal features across a range of different resolutions, galaxy environments and subgrid physics models.As we approach the advent of large observational surveys of galaxies within the sensitivity to observe tidal features, it is important to characterise what each of these simulations predicts concerning these features.TNG and EAGLE probe a similar volume and therefore the same range of environments (isolated galaxies, groups and clusters).The Magneticum simulation probes an intermediate-sized box about one-third the volume of the two larger simulations and has a slightly lower stellar mass resolution.NewHorizon probes a much smaller volume with a much higher stellar mass resolution, allowing for more detailed tidal structures to be resolved by probing a more limited sample of galaxies.
In this section we introduce each of the simulations, describing their mass resolution, volume and how they were calibrated.We discuss how galaxy structure, stellar mass and halo mass are measured in §2.1.1 and the simulations' stellar mass-halo mass (SMHM) relations in §2.1.2.We focus on these elements of the simulations as they are the salient differences for this study, directly impacting the number of stellar particles available to form tidal features and the halo masses that are covered by the simulations.The simulation subgrid models for star formation, stellar feedback and black hole feedback, which could indirectly contribute to the factors above and therefore, the tidal feature occurrence, are described in Appendix A.
The NewHorizon (Dubois et al. 2021) simulation is a zoomin simulation from the parent Horizon-AGN simulation (Dubois et al. 2014), which adopts the WMAP-7 cosmology (Komatsu et al. 2011).NewHorizon employs a spherical volume of 16 3 cMpc 3 with varying dark matter resolution.The initial 10 cMpc radius has  DM = 1.2 × 10 6 M ⊙ , this high-resolution patch is embedded in spheres of decreasing mass-resolution of 10 7 , 8×10 7 and 6×10 8 M ⊙ corresponding to spheres of radius 10.6, 11.7, and 13.9 cMpc.The Table 1.Summary of the properties of the four cosmological hydrodynamical simulations.From left to right, the columns are the simulation name, matter density, dark energy density, Hubble constant, simulation box side length, stellar mass resolution, the redshift corresponding to the simulation timestep at which the particle data used was recorded, the number of galaxies in the simulations with  ★, 30 pkpc ≥ 10 9.5 M ⊙ , the number of galaxies classified, number of compact point-like galaxies discarded, and final sample size for analysis.Horizon-AGN is calibrated to match the black hole density and the scaling relations between black hole mass and galaxy properties (albeit at a lower spatial resolution to NewHorizon).The NewHorizon simulation also reproduces the observed galaxy surface brightness as a function of stellar mass (Dubois et al. 2021), without requiring any calibration for surface brightness.

Simulation
Table 1 shows that NewHorizon has stellar particles ∼2 orders of magnitude less massive than the other three simulations allowing it to resolve tidal features from higher mass ratio mergers.However, NewHorizon only covers field and group environments and does not include galaxy clusters of significant mass (Dubois et al. 2021).
The EAGLE (Evolution and Assembly of GaLaxies and their Environments) project is a large set of cosmological hydrodynamical simulations.EAGLE contains simulations of cubic volume 12 3 , 25 3 , 50 3 and 100 3 cMpc 3 .The data which will be used in this study is from the reference model (RefL0100N1504) which has a volume of 100 3 cMpc) 3 .The reference model is the simulation run with standard parameters and physics as described in Schaye et al. (2015) and We subsample 50 galaxies in each of 9 evenly spaced bins between 11 ≤ log 10 (  200, crit /M ⊙ ) ≤ 14.5.We also show the mean number of satellite galaxies per iteration using dashed lines and the mean number of centrals per iteration using circular points.We similarly depict the distributions of satellite and central galaxies for NewHorizon.The grey histograms show the parent distribution of log halo mass.The error bars give the standard error on the mean for each of the measurements obtained from the Monte Carlo iterations.Crain et al. (2015).The simulation uses the cosmological parameters advocated by the Planck 2013 results (Planck Collaboration 2014).The dark matter particle mass is  DM = 9.70 × 10 6 M ⊙ and the stellar particle mass is  ★ = 1.81 × 10 6 M ⊙ .
The simulation was calibrated to match the galaxy stellar mass function (GSMF) and the relationship between black hole mass and galaxy stellar mass at  ∼ 0 (details in Schaye et al. 2015).
TNG is calibrated to match the star formation rate density, the GSMF and the stellar-to-halo mass relation, the total gas mass contained within the virial radius of massive groups, the stellar massstellar size and the black hole -galaxy mass relations all at  = 0, in addition to the overall shape of the cosmic star formation rate density at  ≲ 10 ( Pillepich et al. 2018b).
Magneticum Pathfinder simulations are a suite of cosmological hydrodynamical simulations, ranging in box size from 25.6 3 to 3818 3 cMpc 3 .We use the Box4-uhr simulation, which models the physics of dark matter and baryons in a 68 3 cMpc 3 box.All simulations are performed with an updated version of the SPH code GADGET-3 (Springel et al. 2005b).We use the Box4-uhr run, which adopts a WMAP-7 fit cosmology (Komatsu et al. 2011).The mass resolution for stellar particles in this run is  ★ = 2.6 × 10 6 M ⊙ and  DM = 3.6 × 10 7 M ⊙ for dark matter particles.
The Magneticum simulations are calibrated to match the intracluster gas content of galaxy clusters, rather than matching any characteristics with individual galaxies.

Identifying structure and measuring stellar mass
To study the properties of galaxies and their tidal features the particles in the simulations must be grouped into central galaxies and their satellites.For TNG, EAGLE and Magneticum, the galaxy and halo structure is computed using the baryonic version of the SubFind code (Dolag et al. 2009, also Springel et al. 2001).The versions of structure finder used across these simulations are similar but not identical.In contrast, NewHorizon utilises AdaptaHOP for finding halos and subhalos (Tweed et al. 2009).It is important to note that the stellar masses shown in Figure 2 and Figure 3 are computed directly from each simulation's native structure finders.They are computed by summing the stellar masses of the particles assigned to each subhalo.
In general, SubFind and AdaptaHOP differ in how they assign particles to halos and subhalos: SubFind identifies substructure based on its local overdensity and gravitational boundedness.In its initial step SubFind uses a Friends of Friends (FOF) algorithm on dark matter particles to distribute them into halos, the baryons are then assigned to the halo of the nearest dark matter particle (if any).Substructures within each FOF halo are identified by searching for local density peaks, considering all particle types and resolution elements.Subhalos are identified using contours of isodensity that pass through limiting saddle points in the density field.Particles/resolution elements not gravitationally bound to the subhalo to which they have been assigned are removed by applying an iterative unbinding procedure.Finally, all particles and resolution elements that are unassigned to a subhalo are assigned to the central subhalo (the largest subhalo within a halo).One notable limitation of the Sub-Find method is that any resolution element/particle lying beyond the limiting isodensity contour is ignored, even if they are gravitationally bound to a subhalo (e.g.Muldrew et al. 2011;Cañas et al. 2019).
AdaptaHOP identifies substructure based on topology and does not carry out any unbinding.Particles are assigned to groups defined by peaks in the density field that are linked to other groups by saddle points in the density field.The substructure is then assigned by repeating this process with an increasing density threshold, resulting in smaller groups within groups.AdaptaHOP defines halos as a group of groups that exceeds 160 times the mean dark matter density.Groups within each halo are regrouped such that each subhalo has a mass smaller than the host halo.
The differences between the algorithms are non-trivial and will lead to differences in the stellar masses of galaxies.AdaptaHOP tends to identify fewer structures than FOF, i.e. several FOF halos are included in one AdaptaHOP halo, these differences can in-turn influence how SubFind assigns subhalos compared to AdaptaHOP (Tweed et al. 2009).We only use galaxy stellar masses derived from the simulation's native structure finding algorithm in Figure 2 and Figure 3, where we compare our visually-classified sample's stellar mass distribution to the parent stellar mass distribution in the simulation.In light of the differences between the algorithms we state that while the stellar mass distributions of both parent and child samples should not be compared trivially across the simulations, we feel comfortable in using the distributions to check how well our classified samples probe the stellar mass distributions of the simulations above a stellar mass of 10 9.5 M ⊙ .We also use the subhalo centres of potential as defined by SubFind and AdaptaHOP for the centres of our LSST-like mock images.Additionally, we rely on the halo assignment when loading particles to produce mock images as we only use particles assigned to the same halo.While for high halo masses (AdaptaHOP halo mass ≥ 10 13 M ⊙ ), FOF tends to assign ≳ 2 times as many halos to a given collection of particles as Adap-taHOP (Tweed et al. 2009).We do not probe such high halo masses with NewHorizon therefore the differences in halo assignment will not have as significant an impact.
From §3.1 onwards, when comparing the occurrence of tidal features as a function of stellar mass, to limit systematic biases between AdaptaHOP and SubFind derived masses we use the 30 pkpc spherical aperture,  ★, 30 pkpc , to measure the stellar masses.As discussed in Schaye et al. (2015), the 30 pkpc spherical aperture is comparable to the 2D Petrosian aperture often used in observational studies, making it a good choice for our aim to analyse the data in a manner directly comparable to observations.The spherical aperture measurement of stellar masses requires only the stellar particles assigned to a subhalo/galaxy within a 30 pkpc sphere centred on the subhalo's centre of potential to be included for the galaxy's stellar mass.
In §3.2, for our analysis of tidal feature fractions as a function of halo mass, we rely on the halo assignment from each simulation's native structure finding algorithms to define main halos, assign galaxies to main halos and the consequent measurement of  200, crit for the main halo. 200, crit is the total mass, including dark matter, gas and stars within the radius where the average density is 200 times the critical density of the universe.For TNG, EAGLE and Magneticum which rely on SubFind, the calculation of  200, crit only relies on FOF for the identification of halos, the calculation of their mass is independent of the halo assignment.In contrast, NewHorizon's ADAPTHOP measured  200, crit exclusively uses the resolution elements assigned to a particular halo.While each subhalo is assigned to a FOF identified host halo in SubFind, AdaptaHOP's use of different density thresholds to separate halo and subhalo structure and the application of minimum density of 160 times the mean dark matter density for a structure containing substructures to be identified as a group leads to some subhalos in NewHorizon not being assigned to any parent halos.

Stellar mass-halo mass relation
The number of stars involved in interactions between halos determines the number of stars available to form tidal features.Therefore, it is important to consider the efficiency of star formation at a given halo mass for each of the simulations.
In Figure 2 we show the SMHM distributions for each simulation, displaying all central galaxies within our stellar mass cut (sample sizes given by  centrals in Figure 2).For comparison, we provide the observationally-constrained empirical SMHM relations from Behroozi et al. (2013, constrained by stellar mass function, cosmic star formation rate and specific star formation rate) and Moster et al. (2013, constrained by the stellar mass function) and the weak lensing derived relation from Hudson et al. (2015) and give the median SMHM relation for each of the simulations.We note that the median relation for each simulation appears to flatten for  200, crit ≲ 10 11 M ⊙ .This is driven by our lower stellar mass cut meaning that in this regime we only sample high mass halos due to not including the scatter to lower stellar masses (e.g.Teklu et al. 2017).
NewHorizon qualitatively reproduces the shape of the observational relations for the halo mass range of the simulation.Relative to the empirical relations NewHorizon does deviate to higher stellar masses below the knee of the SMHM relation.This is at least partly due to the smaller simulation volume (Martin et al. 2022), resulting in under representation of group and cluster environments where star formation tends to be less efficient (e.g.Garrison-Kimmel et al. 2019;Samuel et al. 2022).Another contributing factor to the higher star formation efficiency of NewHorizon could be the subgrid physics implemented (described in Appendix A).We also note recent empirical studies of the low-mass end of the SMHM relation, constrained by the star formation rates of local group dwarf galaxies, have shown significant scatter to higher stellar masses for  200, crit ≲ 10 10 M ⊙ , which is qualitatively similar to what is observed for NewHorizon (O'Leary et al. 2023).
EAGLE, TNG and Magneticum fall within the scatter of the observationally-constrained relations for 10 11 M ⊙ ≲  200, crit ≲ 10 12 M ⊙ .For  200, crit ≳ 10 13.5 M ⊙ for EAGLE and  200, crit ≳ 10 12.5 M ⊙ for TNG and Magneticum we see that the centrals are sitting at systematically higher stellar masses than the observational relations.A factor contributing to this difference could be the differences between stellar mass measurements in simulations and observations.For example, Remus & Forbes (2022) found better agreement between the cosmological simulations EAGLE, Magneticum and TNG-300 and observational measurements of the SMHM relation (Kravtsov et al. 2018) when the light from the central galaxy and the intracluster light (ICL) around it was included in the observations.This deviation between simulations, observations, and models has been known for a while.Kravtsov et al. (2018) showed that the observed SMHM relation deviates from the model predictions if the brightest cluster galaxy is measured out to large radii and their ICL is included.This leads in those cases to larger observed stellar masses at a given halo mass than predicted by Behroozi et al. (2013) and Moster et al. (2013).This is found to be true for many simulations, where the SMHM relation from models and observations that exclude the ICL is usually reproduced if aperture cuts are applied to the simulations, but results in larger stellar masses if the stellar masses from the halo finders are used, that is the ICL and residuals from satellite galaxies are included in the stellar mass measurement of the brightest cluster galaxy (see for example Remus & Forbes (2022) for Magneticum, Pillepich et al. (2018b) for IllustrisTNG, Contreras-Santos et al. (2024) for The Three Hundred Simulations, and Schaye et al. (2023) for FLAMINGOS).And while the chosen apertures of 30-50 pkpc are somewhat artificial, it was recently demonstrated by Brough et al. (2024) by applying observational methods to separate brightest cluster galaxies from their ICL to mock images of simulations from Magneticum, EAGLE, IllustrisTNG and Horizon-AGN, that the best agreement is found for apertures of 50-100 pkpc, providing observational support for this previously artificial aperture cut.
There are small offsets between the SMHM relation for EAGLE, TNG and Magneticum, albeit smaller than the offset with NewHorizon.EAGLE is shifted to lower stellar masses compared to the other simulations and Magneticum shifted to marginally higher stellar masses than TNG at the low and high end of the relation.A detailed comparison between the models was carried out by Remus & Forbes (2022).They suggest that the likely cause of this is the star formation efficiencies of low-mass halos, where Magneticum is observed to convert more gas into galaxies.As a result of its SMHM relation, we expect NewHorizon to bring in more stars into mergers than the other three simulations for mergers containing halos of a given size.

Sample selection
We select all simulated galaxies with  ★, 30 pkpc ≥ 10 9.5 M ⊙ .We make no distinction between centrals and satellites in our sampling, allowing the underlying distributions of centrals/satellites to carry over to our samples.At a stellar mass of 10 9.5 M ⊙ Magneticum resolves these galaxies with ∼ 1200 stellar particles, EAGLE with ∼ 1750, TNG with ∼ 2250 and NewHorizon with ∼ 24300.With any lower stellar mass threshold, there are fewer particles representing galaxies and our mock images will struggle to resolve significant tidal distortions to the stellar components of the galaxies.We also have few observational results to compare to below this stellar mass.
We sampled our galaxies from similar redshift snapshots for TNG, EAGLE and Magneticum.However, for NewHorizon we used the two most recent available snapshots, taking care not to sample galaxies twice (Table 1).The number of galaxies with  ★, 30 pkpc ≥ 10 9.5 M ⊙ in each simulation are given in Table 1.
From the total sample of galaxies with  ★, 30 pkpc ≥ 10 9.5 M ⊙ , we randomly subsample 2000 galaxies to visually classify from each of Magneticum and EAGLE and 1907 galaxies from TNG.This is done so that the visual classification can occur over a reasonable timeframe.Within these galaxies, we note the presence of several galaxies which seemed to be very compact, point-like and have a stellar mass far higher than one would expect for such a compact object (see bottom right of Figure 1).We discarded these compact, point-like objects from our sample, the number discarded is given under  out in Table 1.Our final sample sizes are 62 for NewHorizon, 1978for EAGLE, 1826for TNG and 1990 for Magneticum.
To study how tidal features depend on stellar mass in each of the simulations it is necessary to match the stellar mass distributions between the classified simulation samples to avoid a mass-biased comparison.We subsample from the classified TNG and Magneticum samples to match the stellar mass distribution of the parent EAGLE sample, which was calibrated to match the local GSMF (Furlong et al. 2015).Matching to EAGLE's stellar mass distribution therefore makes our results more directly comparable with observational samples of tidal features.However, for NewHorizon, we do not subsample, leaving the sample as it is as we do not wish to reduce the sample size any further.
To construct our stellar mass-matched sample, we fit the probability densities for the simulation parent stellar mass distributions using Gaussian kernel density estimation (KDE).Due to the unavailability of  ★, 30 pkpc measurements for the parent samples of TNG and Magneticum we use SubFind-derived  ★ values to fit the KDE for the parent stellar mass distributions from EAGLE, TNG and Magneticum.Using these KDEs we infer the likelihoods of sampling galaxies with particular  ★, 30 pkpc from each of the simulations.We then use these likelihoods to construct weights with which we sample the classified TNG and Magneticum galaxies to match the parent distribution of EAGLE.Our samples are constructed using the following steps: (i) We draw, without replacement, N galaxies from EAGLE.(ii) We draw, without replacement, a sample of N galaxies from TNG and Magneticum, weighting the samples so that the likelihood of a galaxy being drawn is the same as it would be if we were sampling from the EAGLE parent stellar mass distribution.
(iii) We redraw N galaxies from TNG and Magneticum until the distributions reasonably match that of EAGLE, by requiring that a KS test (Hodges 1958) yield a p-value > 0.9.
(iv) The largest sample that could be subsampled from TNG and Magneticum, while satisfying condition (iii) is N=1300.
The exact weights we used for subsampling from our classified TNG and Magneticum galaxies are as follows: Where  TNG.( ★, 30 pkpc ) and  Mag.( ★, 30 pkpc ), are the weights assigned to a TNG and Magneticum galaxy respectively, based on their  ★, 30 pkpc . EAG.( ★, 30 pkpc ),  TNG.( ★, 30 pkpc ) and  Mag.( ★, 30 pkpc ) are the probability density functions from which we can calculate the likelihoods of galaxies with a particular  ★, 30 pkpc being sampled from the parent distribution of the respective simulations.
Figure 3 shows the stellar mass distributions of the parent samples for each of the simulations in grey histograms and the distributions of our selected sample in coloured histograms, black for NewHorizon, blue for EAGLE, red for TNG and cyan for Magneticum.TNG and Magneticum are mass-matched to EAGLE.The KS-test p-values for the matching of the  ★, 30 pkpc are 0.98 for EAGLE and TNG and 0.90 for EAGLE and Magneticum.
We also compare our classified samples with the parent distributions from which they were sampled.This is done by using a KS-test to compare their distributions in structure finder-defined stellar mass ( ★ ).The high p-value of 0.35 for our EAGLE sample is expected as we have not attempted to adjust the sampling with respect to its parent distribution, we see TNG (p=0.16) and Magneticum (p=0.001) are less like their parent distributions as they have been sampled to match the 1300 galaxies sampled from EAGLE.We note that the conclusions we draw are not qualitatively changed by stellar mass-matching the samples.
For our halo mass-matched comparison we construct a flat distribution of log halo masses, from 11 ≤ log 10 ( 200, crit /M ⊙ ) ≤ 14.5 for EAGLE, TNG and Magneticum.Due to the small sample size of NewHorizon we do not attempt to perform any subsampling for the simulation.There is greater degeneracy amongst the parent halo masses of the galaxies as many galaxies can belong to the same halo.This makes it more difficult to match the populations with the weighted sampling method described above, so we opt for a Monte Carlo-based approach.To perform this we subsample our classified sample 100 times.We choose a linear bin spacing in log 10 ( 200, crit /M ⊙ ), such that our tidal feature fraction with halo mass measurement would have a sufficient number of galaxies in each bin whilst having a large enough number of bins to resolve trends with halo mass.We found that a subsample of 450 galaxies per simulation, divided across 9 linearly spaced bins in log halo mass was appropriate for our purposes.
The leftmost panel in Figure 4 shows our NewHorizon sample with respect to halo mass with a solid line, we note that NewHorizon does not include higher mass halos, equivalent to galaxy group and cluster mass systems, due to the small volume of this simulation.Furthermore, four of the galaxies in our NewHorizon sample had no parent halo assigned, reducing our sample to 58 galaxies for this part of the analysis.We also consider the number of central and satellite galaxies in our haloes.The central galaxy is the subhalo consisting of the largest number of particles within a parent halo by the structure finder, the other galaxies in the halo are defined to be satellites.In Figure 4 the dashed line shows the number of satellite galaxies in each halo mass bin and the open circles show the number of central galaxies.The three rightmost panels in Figure 4 illustrate the flat sampling of galaxies with respect to their parent halo mass for EAGLE, TNG and Magneticum that we apply where analysing the halo mass and central and satellite galaxy properties of our sample.The solid lines in each of these panels show the 50 galaxies sampled in each of the halo mass bins.The mean number of satellites in each bin is shown using a dashed line and the mean number of centrals is shown using open circles.We see that the number of satellites in each bin increases with increasing halo mass whereas the number of centrals decreases with halo mass.This is expected given that field galaxies are classified as centrals while group and clusters contain many satellite galaxies orbiting around one central galaxy.
Examples of mock-observed images are shown in Figure 1.Described briefly the construction of mock images involves: (i) Extracting all the stellar particles associated with a group (FOF or AdaptaHOP identified halo, depending on the simulation) within a 1 Mpc cube of the simulated galaxy and calculating spectral energy distributions (SED) for each stellar particle using simple stellar population models (Bruzual & Charlot 2003) extrapolated to the age and metallicity of each stellar particle.Dust attenuation is accounted for in Martin et al. (2022)'s method but is not included in the mock images created here.The brightness for each of the stellar particles is determined by the redshift of the galaxy and the convolution between the SED and the transmission function for the relevant LSST filters (Olivier et al. 2008).Using this method we produce apparent magnitudes for the LSST ,  and  filters.
(ii) Smoothing the particle distribution in regions where there is a potential for undersampling in the final image.This smoothing is necessary to remove unrealistic variations in flux between pixels in the mock image.This is done by subdividing stellar particles into smaller mass particles with positions distributed normally around the original particle, with a standard deviation equal to the distance to the 5th nearest neighbour of the stellar particle being smoothed (Martin et al. 2022, see also Merritt et al. 2020).
(iii) Produce mock images by collapsing the cube along one of its axes with a grid size of 0.2 ′′ × 0.2 ′′ (the LSST pixel size; Ivezić et al. 2019).The apparent magnitudes of the light from each of the pixels are measured accounting for the dimming due to the redshift ( = 0.025).
(iv) To model the observational effects of seeing, the image is then convolved with the g-band Hyper-Suprime Cam point spread function (Montes et al. 2021) rescaled to match the LSST pixel size.
(v) To match observational background noise, Gaussian noise is added using the empirical relation between the standard deviation of the background noise and the limiting surface brightness, given in Román et al. (2020).
Where Ω is the size of one side of the box over which the surface brightness limit is computed, = 10 ′′ , pix is the pixel scale in arcseconds per pixel, = 0.2 arcseconds/pixel,  is the number of Gaussian standard deviations for the surface brightness limits, = 3 and  lim band is the limiting surface brightness for a particular photometric band.This model for noise intends to include all possible spatially invariant sources of noise, including instrumental noise, unresolved background sources, scattered light etc. without explicitly modelling them.However, we note that this model does not include spatially varying sources of noise, e.g.sky gradients, diffraction spikes, imperfect sky subtraction or galactic cirrus.Martin et al. (in prep) have shown using Sersic fits to mock images and galaxies in the HSC-SSP observed Cosmos field that this technique for producing mock images does not induce any measurable biases in galaxy structure measured from the mock images.They found that the variance between simulated and observed galaxies is mostly driven by the differences between the galaxy evolution models.

Tidal feature classification
We carry out our classification using single g, r and i band images along with the combined colour image, to give ourselves the best chance of detecting and correctly classifying these faint features.The images were viewed at a fixed range of contrast and brightness as shown in Figure 1.The classification is carried out by the lead author.To classify the tidal features we chose to follow a classification scheme similar to that used in Bílek et al. (2020) and Desmons et al. (2023b).The categories include: • Streams/Tails: Prominent, elongated structures orbiting or expelled from the host galaxy.These usually have similar colours to the host galaxy • Shells: Concentric radial arcs or ring-like structures around a galaxy.
• Plumes or Asymmetric Stellar Halos: Diffuse features in the outskirts of the host galaxy, lacking well-defined structures like stellar streams or tails.
1 Hint of tidal feature detected, classification difficult.
2 Even chance of correct classification of tidal feature presence and/or morphology.
3 High likelihood of the tidal feature being present and morphology being obvious.
Table 2. Descriptions of the classification confidence levels used for our analysis.
• Double nuclei: Two clearly separated galaxies within the mock image where merging is evident through the presence of tidal features.
In this scheme galaxies are allowed to host more than one feature type.A feature is allocated to the galaxy to which it appears to connect to, e.g. a tail/stream is allocated to the galaxy from which it stems, a shell is associated with the galaxy on which it is centred and a asymmetric halo is associated with the galaxy possessing the halo.In cases where features can be associated with multiple objects, e.g. a tail/stream between two galaxies (a tidal bridge) or an asymmetric halo enveloping a double nucleus, the tidal feature is counted for both objects.Following Desmons et al. (2023b) we have merged the categories of streams and tails as they are easily mistaken for one another at the surface brightness depth of LSST (Martin et al. 2022).We also use a system of confidence levels in our classification detailed in Table 2.This system is similar to that followed by Atkinson et al. (2013), but we group their confidence levels 3 and 4 which are represented as confidence level 3 here.Example features are shown in Figure 1.

Calculating the errors on our counts and fractions
We estimate the binomial errors on the tidal feature fractions presented here using 1 ≃ 0.683 confidence levels computed using a Bayesian beta distribution generator for binomial confidence intervals described in detail in Cameron (2011).This approach is robust against small and intermediate sample sizes, which is particularly important due to the small size of the NewHorizon sample.

RESULTS
We present the results of our classification in Figure 5 and Table 3. Figure 5 illustrates the fractions of galaxies exhibiting each kind of tidal feature for the four simulations: EAGLE in blue, TNG in red, NewHorizon in grey and Magneticum in cyan, for each of the confidence levels given in Table 2.The error bars show the 1 binomial errors on the fractions.We find excellent agreement between the four simulations at each confidence level.The differences between the tidal feature fractions are much smaller than the change in the fractions with feature morphology or with confidence level.The simulations all agree on the relative frequencies of tidal features, namely, asymmetries in the halo are the most common feature followed by double nuclei, streams/tails and lastly shells.
As expected, the fraction of galaxies exhibiting tidal features From left to right the panels show the results for each confidence threshold:   .≥ 1,   .≥ 2 and   .= 3 (described in Table 2).The fractions are given for 5 feature categories: total tidal feature fraction (Total), stream/tail (S/T), shell, asymmetric halo (Asym.)and double nucleus (DN).The error bars give the 1 binomial confidence levels.The occurrence of features relative to one another is similar across each confidence level, the fractions of galaxies exhibiting features drop as the confidence level increases.Generally, the simulations are in good agreement, with NewHorizon, tending to have a higher fraction tidal features that the other simulations.
decreases with an increasing classification confidence threshold (  .).The agreement between Magneticum, TNG and EAGLE gets stronger with the increasing confidence threshold, indicating that as we look at tidal features that are obvious/well resolved in the 3 simulations we have an increasing agreement.
We observe the inverse trend with NewHorizon, whilst still agreeing well with the frequencies of tidal features found in the other simulations.For   .≥ 2, NewHorizon systematically sits at higher fractions of tidal features than any of the other simulations.Table 3 shows that between   .≥ 1 and   .= 3, NewHorizon's tidal feature fraction drops by 0.1, whereas in EAGLE and TNG, it drops by 0.24 and in Magneticum the fraction drops by 0.15.The smaller change in tidal feature fractions when moving from the lowest to the highest confidence level shows that tidal features detected in NewHorizon tend to be higher confidence than the other simulations.
The qualitative threshold of confidence levels does not change the relative frequencies of tidal features within simulations, they act as a trade-off between the detectability of a tidal feature by visual classification and the accuracy of its classification.For this reason, we tread a middle path between detectability and accuracy for the remainder of the results and present the moderate confidence results (i.e.  .≥ 2), noting that the choice of confidence level does not qualitatively change the conclusions we draw.
We study the tidal feature fraction as a function of stellar and halo mass.In §3.1 we present the results of our visual classification in a mass-matched comparison of EAGLE, TNG and Magneticum and the entire sample of 62 galaxies from the higher resolution NewHorizon simulation.In §3.2 we study the trends with halo mass, considering this as a proxy for galaxy environment (e.g.Yang et al. 2005).

Tidal features and stellar mass
Table 4 shows the average stellar masses of the tidal feature hosts for every simulation at the moderate classification confidence level (  .≥ 2).These are in good agreement across each of the simulations.The detected shells have the highest overall mean stellar masses across the simulations,  ★, 30 pkpc ∼ 10 11.2 M ⊙ , not including NewHorizon which only has 1 shell in the classified sample.The other tidal features appear to have very similar mean stellar masses.
Figure 6 investigates the fraction of the galaxy population exhibiting tidal features above the LSST surface brightness limits as a function of stellar mass.We do this by constructing stellar mass bins with equal numbers of galaxies.For the stellar mass-matched samples of EAGLE, Magneticum and TNG we have 10 bins containing 130 galaxies each.For NewHorizon the 62 galaxies are distributed across 7 bins of ∼ 9 galaxies.We again see excellent agreement between the mass-matched simulations and broadly good agreement for NewHorizon, though we note that NewHorizon has systematically higher feature fractions, particularly for  ★ > 10 10 M ⊙ .
The top left panel of Figure 6 displays the stream/tail fractionstellar mass distribution of galaxies for each simulation.All simulations exhibit a trend of increasing stream/tail fraction with increasing stellar mass, with very strong agreement between TNG, EAGLE and Magneticum.NewHorizon exhibits the same trend but has systematically higher fractions than the 3 other simulations.However, it is important to keep in mind the much smaller sample size (and therefore larger uncertainties) and higher resolution of NewHorizon when considering the significance of these offsets.
The top right panel shows the shell fraction -stellar mass distributions.There is a notable increase in the fraction of galaxies exhibiting shells at  ★, 30 pkpc ≳ 10 11 M ⊙ for TNG, EAGLE and Magneticum.The singular shell in NewHorizon does not provide sufficient information to determine any trends.
Table 3.The tidal feature fractions for every feature type and classification confidence level (  .≥ 1,   .≥ 2 and   .= 3, described in Table 2).The 5 feature categories presented are: total tidal features, stream/tail, shell, asymmetric halo and double nucleus.The rightmost column gives the overall mean fraction across the four simulations.The uncertainties give the 1 binomial confidence intervals on each of the fractions.* NewHorizon is excluded from the calculation of the mean shell fraction as there is only one feature in that sample.The asymmetric halo fraction -stellar mass distribution for each simulation is plotted in the bottom left panel of Figure 6.Again we see excellent agreement between the simulations and a trend of increasing asymmetric halo fraction with increasing stellar mass.The rate of increase in tidal feature fraction with stellar mass is the largest here (from ∼ 0.2 to ∼ 0.6), indicating a strong dependence on stellar mass.
Lastly, the bottom right panel of Figure 6 shows the double nucleus fraction for each of the stellar mass bins.We see a sharply increasing double nucleus fraction with increasing stellar mass for TNG, EA-GLE and Magneticum, and excellent agreement between the three.For NewHorizon, there is a sharp increase with increasing stellar mass and an indication of a higher fraction in the highest mass bins.

Tidal features and halo mass
The left panel of Figure 7 shows our total tidal feature fractions as a function of halo mass for each of our simulations.We again find excellent agreement between our simulations, with EAGLE, TNG and Magneticum following the same trends and having consistent values for tidal feature fractions in each halo mass bin.For these simulations the fractions of galaxies hosting tidal features increase from halo masses of 10 11 to ∼ 10 12.7 M ⊙ but decreases with increasing halo mass for halo masses above this.NewHorizon follows a similar trend for the halo masses it covers, with systematically higher fractions.
By comparing the location of the peak in tidal feature fractions with halo mass in the left panel of Figure 7 to the proportions of central and satellite galaxies in our sample in Figure 4, we see the peak coincides with the transition from a central-dominated population to a satellite-dominated population.Therefore, we further analyse this trend by comparing the total tidal feature fractions for central and satellite galaxies as a function of the parent halo mass in the right panel of Figure 7.We see again similarly excellent agreement between the simulations, with EAGLE, TNG and Magneticum coinciding in their fractions for both central and satellite galaxies and NewHorizon following similar trends for centrals although at a systematically higher fraction.The fraction of central galaxies exhibiting tidal features increases with increasing halo mass.We note that there is substantial noise in the measurement in the ∼ 10 13.5 M ⊙ halo mass bin as it contains on average ≲ 10 galaxies.The fractions of central galaxies exhibiting tidal features increase to fractions of 0.86 +0.05 −0.2 , 1.0 −0.2 and 0.73 +0.1 −0.2 by  200, crit ≳ 10 13 M ⊙ for EAGLE, TNG and Magneticum.This indicates that at cluster-like halo masses a majority of central galaxies exhibit visually identifiable tidal features.
The satellite populations for EAGLE, TNG and Magneticum show a declining trend in tidal feature fraction with increasing halo mass.There are only 11 satellites in NewHorizon so it is difficult to comment on any trends with halo mass, but for the halo mass bin centred around 10 11.3 M ⊙ we find a significantly higher fraction than the closest halo mass bins for the other three simulations, in line with the higher tidal feature fractions we find for NewHorizon throughout.
Figure 8 shows our specific tidal feature fractions as a function of halo mass.The top left panel shows the stream/tail fraction as a function of halo mass.We find good agreement between the mean distributions for EAGLE, TNG and Magneticum.The fractions of galaxies with features appear to increase from a halo mass of 10 11 M ⊙ to ∼ 10 12.7 M ⊙ , from a fraction of ∼ 0 to ∼ 0.15 -0.2, before decreasing again down to a fraction of ∼ 0 -0.05.NewHorizon follows a similar shape for the range of halo masses it covers, however, with systematically higher fractions.
The top right panel of Figure 8 shows the shell fraction as a function of halo mass.While there are very few shells in the sample, there is a small peak in the shell fraction at a halo mass ∼ 10 12.7 M ⊙ .The single shell identified in the NewHorizon simulation does not allow for any meaningful comment on any potential trend.The peak is at a similar halo mass as for the other tidal features.
The fraction of asymmetric halos with halo mass is presented in the bottom right panel of Figure 8.The shape resembles that seen for streams/tails and shells, albeit the peak is more clearly resolved due to the larger number of asymmetric halos.The peak is at a similar The NewHorizon data is given for the sample of 58 classified galaxies with assigned parent halos and each point represents the feature fraction for a halo mass bin containing ≃ 12 galaxies.The points for EAGLE, TNG and Magneticum represent the mean of 100 Monte Carlo iterations, where each iteration subsampled 450 galaxies to the flat log halo mass distributions depicted in Figure 4.Each of the bins for EAGLE, TNG and Magneticum contain 90 galaxies, these galaxies are a varying mix of central and satellite galaxies as can be seen in Figure 4.The error bars show the mean 1 binomial errors.
halo mass (∼ 10 12.7 M ⊙ ) to the peak for streams/tails and shells, and the fraction is larger in size (maximum asymmetric halo fraction ∼ 0.5).NewHorizon appears to again have a systematically higher fraction for each halo mass, whilst following the same trend for the halo masses it covers.Lastly, the bottom right panel presents the double nucleus fraction as a function of halo mass.We again see a steady increase, peak and decrease in feature fraction, similar to the other tidal feature morphologies.The halo mass peak is consistent with the peak found for the other features.NewHorizon follows a similar trend for the halo masses it covers with a systematically higher double nuclei fraction.
Generally, we see a trend of increasing tidal feature fraction with halo mass to  200, crit ∼ 10 12.7 M ⊙ , after which the fraction appears to decrease.NewHorizon appears to resemble the distributions of the other simulations for the range of halo masses it covers, despite the lack of mass matching and much smaller sample size.
Figure 9 shows the specific tidal feature fractions for the central and satellite galaxy populations.For each of the feature types, we see that the tidal feature fractions of central galaxies tend to increase with increasing halo mass for each simulation, with the exception of the stream/tail fraction in the halo mass bin centred around ∼ 10 13.5 M ⊙ for EAGLE and Magneticum.This could be a consequence of the small number statistics of central galaxies in this bin, on average there are only 10.21 +0.04 −0.03 and 5.54 +0.04 −0.02 centrals in that halo mass bin for EAGLE and Magneticum respectively.We find that shells appear almost exclusively in central galaxies in all the simulations.
For tails/streams in Figure 9, satellite galaxies appear to follow a flat distribution with halo mass.For asymmetric halos we find that there is some evidence for a small peak in the distribution at  200,crit ∼ 10 12.7 M ⊙ .We see similar evidence for a peak in the distributions of double nuclei fractions in the satellite populations of TNG and Magneticum.However, the distribution of the fraction of satellites exhibiting double nuclei from EAGLE remains flat with respect to halo mass.

DISCUSSION
We have presented a visual classification of the tidal features around galaxies in four cosmological hydrodynamical simulations: NewHorizon, EAGLE, TNG and Magneticum.Our visual classification finds broad agreement between the four simulations regarding the tidal feature fractions and their behaviour with stellar mass and halo mass, while showing systematically higher tidal feature fractions for NewHorizon in the moderate and highest confidence levels.In §4.1 we compare our results to previous studies of tidal features in simulations and in §4.2 we compare them to observational results to ensure that they are consistent with previous work.In §4.3, we discuss the impact of our visual classification method, namely having the classifications done by an individual rather than a group.In §4.4,we discuss the possible implications of our results on the occurrence of tidal features in cosmological hydrodynamical simulations and in §4.5 we explore the potential of an environmental dependence for tidal feature occurrence in our results.We go on to give testable predictions for LSST in §4.6.

Comparison with other simulated analyses
We check whether our results are consistent with previous analyses of tidal features around galaxies in simulations.Martin et al. (2022) performed a study of visually-classified tidal features using NewHorizon across a similar range of stellar masses to our study.Their original sample consisted of 37 galaxies from the  = 0.2 The NewHorizon data is given for the sample of 58 classified galaxies with assigned parent halos and each point represents the feature fraction for a halo mass bin containing ≃ 12 galaxies.The points for EAGLE, TNG and Magneticum represent the mean of 100 Monte Carlo iterations, where each iteration subsampled 450 galaxies to the flat log halo mass distributions depicted in Figure 4.Each of the bins for EAGLE, TNG and Magneticum contain 90 galaxies.The error bars show the mean 1 binomial errors.We see excellent agreement between the four simulations plotted.  .Distribution of tidal feature fractions for satellite and central galaxies as a function of halo mass.From left to right, top to bottom, we have the: stream/tail, shell, asymmetric halo, and double nucleus fractions.The results for NewHorizon are given in black, EAGLE in blue, TNG in red, and Magneticum in cyan.The NewHorizon data is given for the sample of 58 classified galaxies with assigned parent halos.The dashed lines show the feature fractions for satellite galaxies and the circles show the feature fractions for centrals for each halo mass bin containing ≃ 10 galaxies.For the Monte Carlo sampled distributions the points for centrals and satellites show the means of 100 iterations, each bin contains 90 galaxies.Each bin contains a varying mix of central and satellite galaxies as can be seen in Figure 4.The error bars for NewHorizon shows the 1 binomial errors and the error bars for EAGLE, TNG and Magneticum data show the mean 1 binomial errors.
timestep and their progenitors at  = 0.4, 0.6 and 0.8, leading to a sample of 148 galaxies.They had 45 astronomers visually classify the images and varied galaxy orientations relative to the observer and the distances at which the object was placed relative to the observer in order to gather ∼ 8000 images, 600 of which were classified by up to 5 different classifiers.They applied an extended classification scheme, including bridges and merger remnants in addition to the classifications we use here.Comparing our results for the fractions of galaxies exhibiting at least one instance of a tidal feature morphology (Figure 6), we see qualitatively similar trends of the fraction Fraction of galaxies  Martin et al. (2022).For each simulation, each mass bin contains a similar number of galaxies,  bin,NewHorizon ≃ 9 and  bin,EAGLE =  bin,TNG =  bin,Magneticum = 130.Martin et al. (2022) bin their galaxies using stellar mass bins of equal width.
of galaxies exhibiting streams/tails, double nuclei and asymmetric halos all increasing with increasing stellar mass.
The most easily comparable measure between the two studies is the total tidal feature fraction.We compare to the Martin et al. (2022) results for all tidal features, excluding asymmetric halos to align with their approach, in Figure 10.The results are depicted similarly to Figure 6, with identical binning and including only classifications with   .≥ 2. Our NewHorizon results agree well for stellar masses ≥ 10 10 M ⊙ .For stellar masses below 10 10 M ⊙ , we report lower tidal feature fractions, with our fractions being lower by ∼ 0.4.For galaxies with 10 9.5 M ⊙ ≤  ★ ≤ 10 10 M ⊙ , Martin et al. (2022) sample 7 galaxies from the  = 0.2 timestep, excluding the remaining galaxies due to their halos being contaminated by lower resolution dark matter particles.This contamination is exclusively a concern for zoom-in simulations and therefore only impacts NewHorizon in our sample.We include these galaxies, sampling 32 galaxies in this range to maximise our number statistics for NewHorizon.This difference in sampling is the likely cause of the differences in the measured fraction.
We further check for consistency by comparing our to those from Valenzuela & Remus (2023).They visually classified the presence of streams, tails and shells in the same Magneticum simulation box we use for this study.Valenzuela & Remus (2023) performed their visual classification on the 3D stellar component of Magneticum galaxies with  ★ ≥ 2.4 × 10 10 M ⊙ and virial mass  vir ≥ 7.1 × 10 11 M ⊙ ( vir ≈  200, crit ).The most direct comparison to our results is using the total tidal feature fractions.In Figure 11, we show the fraction of galaxies exhibiting tidal features in each of our samples, for   .≥ 1 and   .= 3 classifications.The confidence levels increase with the opacity of the point.The left-hand panel compares the simulated tidal feature fractions.The stellar mass error bars for each of the fractions give the range of stellar masses studied (using the  ★, 30 pkpc masses for our results), while the point itself gives the fraction and the mean stellar mass where the measurement is available, and the midpoint of the stellar mass range where the mean is unavailable.Valenzuela & Remus (2023) measured a tidal feature fraction of 0.23 ± 0.02, which agrees well with our tidal feature fraction for Magneticum and other simulations at   .= 3. Limiting ourselves to just tails, streams and shells within the stellar mass range covered by Valenzuela & Remus (2023), we find a feature fraction of 0.17 ± 0.02 at   .≥ 1 and 0.09 +0.02 −0.01 for   .= 3.Therefore, our fractions sit a factor of ∼ 1.4 − 2.5 times below the fractions for the same Magneticum simulation box.Valenzuela & Remus (2023) also find an increasing fraction of galaxies exhibiting tidal features with stellar mass, qualitatively agreeing with our predictions.The differences in our fractions could be due to us visually detecting and classifying tidal features from mock images while Valenzuela & Remus (2023) classified based on unprojected 3D data.Mock images are subject to projection effects and the surface brightness limit we have applied which Pop et al. (2018); Kado-Fong et al. (2018) and Martin et al. (2022) suggest have a significant impact on the detectability of tidal features.Therefore, lower fractions from these projected mock images are expected.Pop et al. (2018) studied shell formation in an analogous run used from the previous generation Illustris simulation which has a similar resolution to that used in this study (Genel et al. 2014;Vogelsberger et al. 2014).They performed their analysis on 220 of the most massive galaxies in the Illustris simulation ( ★ ≳ 10 11 M ⊙ ) and visually detected the presence of shells using stellar surface density maps across 3 orthogonal projections.Their shell fraction of 0.18 +0.03 −0.02 is presented in the left panel of Figure 11.Furthermore, Pop et al. (2018) predict that in an observational survey (i.e. with only one projection available) they would measure a shell fraction of 0.14±0.03.For a direct comparison, we take our most optimistic shell fraction (  .≥ 1) for our stellar mass-matched TNG sample.This gives a substantially lower fraction of shells than Pop et al. (2018), 0.005 +0.003 −0.001 (Table 3).From Figure 11, we can see that our studies cover different mass ranges.Pop et al. (2018) probes galaxies that are much higher in stellar mass than ours.As shell fractions increase with stellar mass, we find it likely that the different stellar mass ranges of the samples contribute to the different tidal feature fractions.To test the impact of the different stellar mass ranges, we extrapolate our shell fraction to the midpoint of Pop et al. (2018)'s stellar mass range, 10 12.3 M ⊙ .While our   .≥ 2 results plotted in Figure 6 only have measured shell fractions in the highest stellar mass bin, if we include   .≥ 1 classifications we obtain measurements in the last two stellar mass bins of 0.008 +0.02 −0.002 and 0.05 +0.03 −0.01 .An approximate linear gradient of 0.1 ± 0.1 can be obtained using the mean galaxy stellar masses in these two bins.Using this as an estimate for the rate of shell fraction increase with stellar mass for our sample, we can estimate what our total shell fraction would be if our stellar masses are extrapolated from the mean mass in our highest stellar mass bin to the midpoint of Pop et al. (2018) stellar mass range.This gives a shell fraction 0.18 +0.17 −0.15 , which is consistent with the fraction measured by Pop et al. (2018), suggesting that the differences are mainly due to the stellar mass ranges of the respective samples.
Figure 8 shows the relationship between halo mass and tidal feature occurrence.It suggests that we may be able to test the rates of close encounters as a function of the environment predicted by N-body simulations (e.g.Gnedin 2003) using observations from LSST.The simulations predicted that the rate of close encounters per galaxy as a function of the relative velocities of the galaxies during the encounter has a peak at an encounter velocity ∼ 500 km s −1 .This corresponds to 2 galaxies moving at speeds ∼ 250 km s −1 relative to the centre of mass of the cluster.The halo mass-velocity dispersion relationship shows that  ≃ 250 km s −1 corresponds to a peak at  vir ∼ 10 13 M ⊙ (Elahi et al. 2018).We expect there to be a corresponding peak in the relationship we measure between tidal feature occurrence and halo mass.As tidal features are tracers of galaxy encounters, the fact that we see a peak at  200,crit ∼ 10 12.7 M ⊙ in Figure 8 is an encouraging prediction from our simulations.Jian et al. (2012) studied close pairs in the Millenium N-body simulation and a suite of semi-analytic models.They found a peak in the relationship between the fractions of close pairs undergoing mergers and halo mass between 10 12 -10 13 M ⊙ .This is in good agreement with our double nuclei panel in Figure 8.
We find good agreement between the tidal feature fractions in our simulations, suggesting that tidal features are a genuine probe of mass assembly and may not be sensitive to processes being modeled by subgrid physics.Canas et al. (2020) and Proctor et al. (2023) studied intra-halo light (IHL) in Horizon-AGN and EAGLE respectively and Brough et al. (2024) studied the IHL in Horizon-AGN, Hydrangea, Magneticum and IllustrisTNG.These works all found good agreement in IHL fractions despite the different subgrid physics models between simulations.As IHL is a similarly low surface brightness feature formed through mergers, the agreement between IHL fractions across these simulations suggests that gravitational physics is playing a dominant role in this probe of assembly history.
In general, we find an encouraging agreement with existing simulation results.The differences in the tidal feature fractions in our results and existing studies are likely explained by a combination of the differences in stellar mass ranges and the potential reduction of tidal feature detectability due to the surface brightness limits of the mock images.

Comparison with observational results
Observed tidal feature fractions measured by prior studies span a large range of values.Most of these differences are likely driven by the different surface brightness limits of the surveys, as tidal features become more visible with deeper imaging and at higher galaxy stellar masses (e.g.Kado-Fong et al. 2018;Martin et al. 2022).To ease comparison with the values measured here we have plotted the total tidal feature fractions measured in observational studies in the righthand panel of Figure 11.Atkinson et al. (2013) used observations from the wide component of the Canada-France-Hawaii Telescope Legacy Survey.They visually classified g'-r'-i' stacked images of 1781 galaxies with a redshift range of 0.02 <  < 0.4.The limiting surface brightness of their data was   ′ AB = 27.7 ± 0.5 mag arcsec −2 .The classification scheme used by Atkinson et al. (2013) identifies; streams, arms, linear features, miscellaneous structure, diffuse fans and shells and has significant overlap with our own classification types.Our streams/tails category is roughly analogous to their streams + arms + linear features, and our asymmetric halo category corresponds to their diffuse fans + miscellaneous structure.For their lowest confidence level Atkinson et al. (2013) measure a total tidal feature fraction of 0.37±0.01 and at their highest confidence level, they measure 0.18 ± 0.01.Their lowest confidence level sits below our   .≥ 1 fractions for NewHorizon, EAGLE and TNG: 0.45 ± 0.06, 0.48 ± 0.01 and 0.44 ± 0.01.However, their fraction is comparable to Magneticum: 0.38±0.01,which presents a tidal feature fraction slightly lower than the other simulations, likely due to its slightly lower resolution discussed in more detail in §4.4.Their highest confidence level tidal feature fraction sits significantly below our   .= 3 fractions for NewHorizon, EA-GLE, TNG and Magneticum: 0.35 ± 0.06, 0.24 ± 0.01, 0.20 ± 0.01 and 0.23 ± 0.01.This is expected given Atkinson et al. (2013) examined data at a limiting surface brightness ∼ 2.6 mag arcsec −2 brighter than our mock images.
There have been many studies of tidal features made using HSC-SSP observations, such as the work of Kado-Fong et al. (2018), Huang & Fan (2022) and Desmons et al. (2023b).Kado-Fong et al. (2018) studied 20,000 galaxies across a redshift range of 0.05 <  < 0.45 using the first data release of the HSC-Wide catalogue (Bosch et al. 2018).They used an automated method, creating high spatial frequency images to detect tidal features and then visually classifying the morphologies of these features.They found a much lower total tidal feature fraction of 0.056 ± 0.002.They cover a much larger range of stellar masses than our sample and detect tidal features to surface brightness depths of   ∼ 26.4 mag arcsec −2 , 3.3 mag arcsec −2 shallower than our mock images.With their stellar mass range probing down to 10 8 M ⊙ and redshift range reaching  ∼ 0.45, many tidal features will fall below their detection limit resulting in their lower fraction.
Huang & Fan (2022) used the coadded images from the Wide layer of HSC-SSP third data release, which reaches a surface brightness limit of   ∼ 28.5 mag arcsec −2 (∼ 1.2 mag arcsec −2 brighter than our mocks).Their final sample of galaxies consists of 2649 earlytype galaxies without companions (i.e. the sample has no double nuclei) with stellar masses ≥ 10 11 M ⊙ .They apply an automated procedure to model and remove the host galaxy and leave behind only the tidal features.Their total tidal feature fraction is 0.28 ± 0.01.Their fraction is slightly higher than our   .= 3 fractions for all of our simulations other than NewHorizon.Their fraction falls ∼ 0.2 below our   .= 1 fractions, this indicates that when we include all our galaxies that exhibit any trace of a tidal disturbance, even where it is difficult to determine its morphology, we visually detect far more than Huang & Fan (2022).This is likely due to the deeper surface brightness limits applied to our images.Desmons et al. (2023b) used the HSC-SSP second data release UltraDeep region images, with a limiting surface brightness of   ∼ 30.76,   ∼ 29.82 and   ∼ 29.41 mag arcsec −2 (3; 10 ′′ × 10 ′′ ).They limited their sample to galaxies in the redshift range 0.04 <  < 0.2 and applied stellar mass limits of 10 9.5 M ⊙ ≤  ★ ≤ 10 11 M ⊙ , constructing a volume-limited sample of 852 galaxies.Their classification scheme is identical to ours, making this a particularly relevant comparison.They measure a total tidal feature fraction of 0.23 ± 0.02.This is in good agreement with our   .= 3 results for EAGLE, TNG and Magneticum and ∼ 0.1 lower in fraction than NewHorizon, while we measure significantly higher fractions when including all possible tidal features (  .≥ 1).Given their lower stellar mass range, it could be that the differences in fractions are driven by the contribution of simulated galaxies with  ★, 30 pkpc ≥ 10 11 M ⊙ that tend to have higher tidal feature fractions than lower stellar mass galaxies.Furthermore, all the galaxies in the Desmons et al. (2023b) sample, are further from the observer than in our sample.Galaxies at  = 0.2 are further away than the galaxies in our mock images, more distant tidal features may be too dim to detect visually.Desmons et al. (2023b) found qualitatively similar trends with increasing tidal feature fractions as a function of stellar mass.They do detect a substantially higher fraction of galaxies exhibiting shells, 0.02 ± 0.01, higher than our mean   .≥ 1 shell fraction of 0.006 +0.001 −0.0007 (excluding NewHorizon as we only detect one shell in this simulation).This result is surprising as the Desmons et al. (2023b) sample is measured over lower stellar masses than ours.Given that shell fraction increases with increasing galaxy stellar mass we would expect our sample to return a higher shell fraction.Total fraction of galaxies exhibiting tidal features from various observations and simulations studies compared to the results presented here.On the left, we display our own simulation results and simulation results from the literature.On the right, we give the results of observational studies.The error bars along the stellar mass axis show the range of stellar masses for each measurement, with the marker situated on the mean where the sample mean stellar mass is provided (where unavailable the midpoint of the stellar mass range is used instead).Our   .≥ 1 results are presented with a transparent marker and our   .= 3 results are presented with the solid marker.We use a similar opacity to indicate the different confidence level derived fractions from observations.The dashed line links our TNG   .≥ 1 shell fraction with an estimated fraction when linearly extrapolated to higher stellar masses.We find our tidal feature fractions to be in the expected range for observational results and simulation results.
Constraining our sample to  ★, 30 pkpc ≤ 10 11 M ⊙ , we find EAGLE, TNG and Magneticum to have mean   .≥ 1 shell host masses of log 10 ( ★, 30 pkpc /M ⊙ ) = 10.8 ± 0.2, 10.7 ± 0.1 and 10.8 ± 0.1.These are consistent within 2 of the mean shell host stellar mass of log 10 ( ★ /M ⊙ ) = 10.56 ± 0.04, found by Desmons et al. (2023b).The higher fractions of galaxies hosting shells in observations could be a result of factors relating to the properties of the host galaxies, the characteristics of the shells themselves and the methods we use to produce the mock images.The detectability of shells was examined in Martin et al. (2022) and Bazkiaei (2023), who found that it depends significantly on the degree of contrast against the background galaxy.The variation in detectability is driven more by the physical parameters of a galaxy (stellar mass, Sersic index and galaxy size) than by the characteristics of the shells themselves (opening angle and shell width).Bazkiaei (2023) find that detection of shells is dependent on the stellar mass resolution of the simulations, which limits the amount of contrast a shell can achieve against the host galaxy in the inner radii close to the galaxy.Martin et al. (2022) found that within ∼ 4  eff., shells will not be detectable in NewHorizon mock images but could be detectable in LSST.We expect this region to extend to larger radii for the other simulations as they have lower stellar mass resolutions (Bazkiaei 2023).
It is possible that systematics introduced by the mock image construction (e.g.smoothing) could reduce the visibility of some shells, making shells more difficult to detect or potentially be misidentified as a different class (e.g.asymmetric halo).However, the agreement of our TNG shell fractions with the shell fractions measured directly from Illustris simulation data by Pop et al. (2018), suggests that the mock image creation method does not significantly impact the visibility of shells.Direct comparison with LSST data could help clarify this tension between simulation and observed shell fractions.

Discussion of visual classification methodology
The majority of the visual classification was performed by the lead author.There was consultation on a subset of 31 NewHorizon galaxies to ensure good agreement on the visually classified morphological types of features.Other studies have already used visual classifica- This illustrates that there are no differences in the relations followed by central or satellite galaxies across all simulations.The bottom left panel is identical to the right panel of Figure 7 and shows the tidal feature fraction as a function of halo mass for central and satellite galaxies.The bottom right panel plots the tidal feature fraction estimated from the tidal feature fraction-stellar mass relation shown in the top left panel, for the populations of centrals and satellites as a function of mean halo mass.For central galaxies we find reasonable agreement between the bottom two panels, indicating that the trends seen in the tidal feature fraction-halo mass relation are a result of the relationship of the tidal feature fractions with stellar mass.However, satellites in the bottom right panel have fractions that systematically increase with increasing halo mass compared to the bottom left panel, where they increase and then fall as a function of increasing halo mass.This could indicate that there are factors other than stellar mass that drive tidal feature fractions for satellite galaxies.The uncertainties in the tidal feature fraction are given by the 1 binomial errors for the relations with stellar mass and the mean 1 binomial errors for the relations with halo mass.The uncertainties on the mean halo mass are given by standard deviation of the mean halo mass over the 100 Monte Carlo iterations.
tions by groups of scientists and reported on the stability and convergence of such classification systems.Bílek et al. (2020) used a group of scientists of different experience levels to identify streams, shells and tails on the MATLAS Survey data.They found that there is a significant scatter between classifiers and that this scatter tends to decrease with the experience of the classifiers.Martin et al. (2022) also found that with faint tidal features, there was decreased concurrence between classifiers on the nature of the morphologies, increasing the scatter.The above suggests that using multiple classifiers helps to characterise the scatter in the morphology and detectability of the tidal feature and could increase the accuracy and reliability of visual detection and classification.However, we argue that having a single classifier makes for a similarly biased comparison across all of the simulations, as the biases of the single classifier are carried across the simulations as systematic errors, allowing for a fair consensus on whether the simulations agree with each other.

Comparison between the different simulations
The broad agreement between the four different simulations analysed here, and the behaviour of these fractions with respect to stellar mass and halo mass suggests that the different subgrid physics and the different hydrodynamics schemes applied are not playing a dominant role in the observability of tidal features.The impact of classification confidence on the feature fractions is presented in Table 2 and Figure 5. Importantly we see that the decrease in the fractions of galaxies exhibiting tidal features with an increasing confidence level is less significant in the NewHorizon simulation than in the other simulations.This implies that the classified tidal features are of a higher confidence and therefore easier to detect in NewHorizon than the other simulations.A contributing factor to the higher confidence tidal features in the NewHorizon simulation could be its ∼ 2 orders of magnitude higher stellar mass resolution compared to the other simulations.For a given number of particles, NewHorizon would be able to resolve tidal features from disrupted galaxies with ≳ 100 times lower stellar mass than the other simulations, which could result in the higher fractions (Appendix B).The slightly lower mass resolution for Magneticum when compared to the other simulations could also be contributing to its lower fraction when including   ≥ 1 tidal features that are likely to include fainter features where the mass resolution could significantly impact its visibility when near a bright host (Appendix B).
The galaxy SMHM relation (Figure 2), points to a different potential source of NewHorizon's higher tidal feature fractions as a function of stellar and halo mass.In this relation, we see that for a given parent halo mass NewHorizon has on average higher stellar mass galaxies populating the halo compared to the other simulations.This means that when galaxies undergo mergers, NewHorizon will tend to bring more stellar mass into the merger and therefore, is likely to produce more visible tidal features for a galaxy of a particular stellar mass than the other simulations.
These results suggest that the simulations are in reasonable agreement with one another regarding the occurrence of visuallydetectable tidal features.The higher resolution and higher galaxy stellar masses for a given halo mass of NewHorizon could be driving the higher fractions we measure.

Tidal feature fractions and environment
The relationships between tidal feature fractions and halo mass, plotted in the left hand panel of Figure 7 illustrate a peak in the occurrence of tidal features at  200, crit ∼ 10 12.7 M ⊙ .This peak coincides with the transition from a central-dominated population of galaxies to a satellite-dominated population.These observations together suggest that the the occurrence of tidal features could depend on environment both measured as a function of halo mass and as location in the halo.In Figure 12 we investigate whether the trends we measure with halo mass are indicative of an environmental influence on the occurrence of tidal features or whether the observed trends can be reproduced by accounting for the trends of tidal features with galaxy stellar mass.
We start in the top panel by investigating whether there is any separation in the tidal feature-stellar mass relation when dividing it into central and satellite galaxies.The top left panel is a version of Figure 10, including asymmetric halos and with the highest stellar mass bin in the EAGLE, TNG and Magneticum simulations, split into two bins linearly spaced in log-stellar mass to allow for the trends in the higher stellar mass ranges ( ★, 30 pkpc ≳ 5 × 10 10 M ⊙ ) to be explored.The top right panel shows tidal feature fractions with galaxy stellar mass for central and satellite galaxies in each of the bins in the top-left panel.We see for all the simulations that the tidal feature fraction with stellar mass does not depend significantly on whether the galaxy is a satellite or a central.
In the lower panels we explore whether the relationship observed with halo mass in Figure 7 could be a result of the observed relationship with stellar mass.We repeat Figure 7 in the bottom left panel for direct comparison.The bottom right panel shows the tidal feature fraction estimated using the mean stellar masses of the central and satellite galaxies in each of the bins in the bottom left panel.The trends seen for central galaxies in the bottom left panel are reproduced well in the bottom right panel, suggesting that their tidal feature fractions are mostly driven by their stellar mass.The satellite tidal feature fractions inferred from their stellar mass systematically increase with increasing halo mass reaching 0.3 -0.4 for  200, crit ≳ 10 13 M ⊙ .However, the true satellite tidal feature fractions increase with increasing halo mass to  200, crit ∼ 10 12.7 M ⊙ and decrease with increasing halo mass beyond this mass.While this is within 2 of the uncertainties on these fractions, the systematic nature of the difference suggests there might be factors other than stellar mass driving the relationship between tidal feature fractions and halo mass for satellite galaxies.This might indicate an enhancement of mergers in satellites residing in halos of  200, crit ∼ 10 12.7 M ⊙ and a reduction in mergers in satellites residing in halos of  200, crit ≳ 10 13 (e.g.Jian et al. 2012).The work of Omori et al. (2023a) provides another example of galaxy mergers rates exhibiting a dependence on environment.Omori et al. (2023a) used a deep learning model to detect mergers in HSC-SSP, finding that on scales of 0.5 to 8 Mpc h −1 , mergers tended to occur in lower density environments, showing that galaxy environment is a significant factor in the occurrence of mergers.The role of environment in the occurrence of galaxy mergers is further highlighted by Sureshkumar et al. (2024).Sureshkumar et al. (2024) studied the spatial clustering of merger and non-merger galaxies in the GAMA catalogue and found that galaxy mergers tend to occur in under-dense environments on spatial scales greater than 50 kpc h −1 .

Predictions for LSST
Given the broadly good agreement between the simulations regarding tidal feature fractions and their trends with stellar mass and halo mass, we can be confident in the simulations' predictive power for comparison with LSST.We note that NewHorizon's higher tidal feature fractions in some relations could potentially indicate that the higher mass resolution of this simulation allows it to resolve tidal features from lower stellar mass progenitors.We use a simple model to probe the limiting progenitor mass for tidal features that each simulation can detect in the mock images in Appendix B. In NewHorizon we can detect tidal features down to stellar masses of ∼ 10 6 M ⊙ , so long as they are not disrupted over an area ≳ 36 pkpc 2 and in the other simulations this limit is ∼ 10 8 M ⊙ with the tidal feature being disrupted over an area ∼ 2500 pkpc 2 so long as the light from the host galaxy is at a similar surface brightness as the feature or fainter within the same region as the tidal feature.This could suggest that the predictions here offer a lower bound on the potential detections in LSST.However, the higher stellar mass for a given halo mass in NewHorizon's SMHM relation as well as its smaller volume could also contribute to the higher fractions.The fact that the three simulations of a similar resolution are in agreement regarding the tidal feature fractions and their trends with stellar and halo mass, indicates that the simulations are well converged and offer robust predictions down to the lowest stellar mass tidally disrupted galaxy that can be both resolved in the simulation and visibly detected in LSST (e.g.Martin et al. 2022).
Given that the simulations offer robust predictions, complete to their lowest stellar mass tidally disrupted galaxy that can be resolved in our mock images, they offer at a minimum a lower bound on tidal feature fractions should LSST be able to probe tidally tidal features below simulation resolution limits which will depend on observational processing (Watkins et al. 2024).We inspect Figure 6 and Figure 7 to examine what the simulations may predict with respect to visually-detected tidal features in LSST.The increasing tidal feature fraction with stellar mass is commonly observed in both simulations and observations.This trend supports the observational evidence that the visual detectability of tidal features increases with the stellar mass of the host.A comparison of our trends of tidal feature fraction with stellar mass against future LSST observations will enable us to further test whether these predictions are consistent with the real Universe.We can also probe whether the tidal feature fraction relation with halo mass for satellite galaxies has any dependence on environment.
LSST will survey billions of galaxies.To take full advantage of this unprecedented data set we will need automated detection and classification of tidal features.There is ongoing work on automated detections of tidal features (e.g.Kado-Fong et al. 2018;Pearson et al. 2019;Desmons et al. 2023a;Omori et al. 2023b).Desmons et al. (2023a) have developed a self-supervised machine-learning method to detect tidal features in HSC-SSP galaxies.Further work to convert these detections into classifications of specific tidal feature morphologies would enable the measurement of trends of total and specific tidal feature fractions with stellar and halo mass over large samples for which visual classification would be prohibitively timeconsuming.

CONCLUSIONS
We have presented the results of our visual classifications of tidal features around galaxies with 10 9.5 M ⊙ ≤  ★ ≲ 10 12 M ⊙ in mock images made to predicted LSST 10-year surface brightness depths of   ∼ 30.3 mag arcsec −2 ,   ∼ 30.3 mag arcsec −2 and   ∼ 29.7 mag arcsec −2 from four cosmological hydrodynamical simulations.We compared the tidal feature fractions and their behaviour as a function of stellar and halo mass.We draw the following conclusions from our results: • Tidal feature fractions decrease with increasing classification confidence levels.This highlights the need to specify a confidence level of visual classification appropriate for the analysis.We present our results using   .≥ 2 (see Table 2).
• The total tidal feature fractions are consistent between the simulations:  NewHorizon = 0.40 ± 0.06,  EAGLE = 0.37 ± 0.01,  TNG = 0.32 ± 0.01 and  Magneticum = 0.32 ± 0.01 (Table 3).This shows that the occurrence of visually identified tidal features in the cosmological hydrodynamical simulations may not be sensitive to the different subgrid physics applied by each.
• The higher tidal feature fractions for NewHorizon as a function of stellar and halo mass may be driven by differences between its galaxy SMHM relation and that of the other simulations (Figure 2) as well as its higher stellar mass resolution (Appendix B).
• The impact of simulation stellar mass resolution is similar for all feature types, generally enhancing the likelihood of detecting a feature around a given galaxy.
• Tidal feature fractions increase with increasing halo mass to a peak mass of  200, crit ∼ 10 12.7 M ⊙ before declining with increasing halo mass (Figure 8).
• Central galaxies exhibit increasing tidal feature fractions with increasing halo mass, a trend that can be accounted for by the relationship between tidal feature fraction and galaxy stellar mass.We expect over half of the central galaxies in haloes with masses ≳  200,crit 10 12.7 M ⊙ to have evidence of at least 1 tidal feature.
• Satellite galaxies tend to exhibit a declining tidal feature fraction for  200,crit ≳ 10 13 M ⊙ , and this trend cannot be fully accounted for by the stellar masses of the galaxies alone, suggesting a potential additional environmental dependence (Figure 12).
• Comparison with observations indicate that our results are consistent with the tidal feature fractions that we would expect from LSST, complete to the tidal features with stellar masses ≳ 10 6 M ⊙ for NewHorizon that are disrupted to an area ≲ 36 pkpc 2 and ≳ 10 8 M ⊙ for EAGLE, TNG and Magneticum that are detectable even if disrupted to an area of ∼ 2500 pkpc 2 , provided the light from the host in the same region is of a similar surface brightness or fainter than the tidal feature (Appendix B).There is some indication from our comparisons with Desmons et al. (2023b) that shells might be detectable at lower host galaxy stellar masses in observations than in simulations.
We have produced predictions from cosmological simulations that are directly testable using LSST data, in particular exploration of trends of tidal feature fractions as a function of stellar and halo mass.
Of particular interest will be using observed tidal features as a proxy for close encounters and mergers between galaxies, and determining how this fraction depends on halo mass to probe the galaxy assembly histories and the impact of the environment on them.

A2 EAGLE
EAGLE runs on the Lagrangian Smooth Particle Hydrodynamic (SPH) code GADGET-3.Star formation is implemented stochastically, following the pressure law scheme of Schaye & Dalla Vecchia (2008), implementing the observed Kennicut-Schmidt (Schmidt 1959;Kennicutt 1998) relation into the simulation.Star formation is set to occur in regions with a metallicity-dependent hydrogen mass density threshold, formulated by Schaye (2004), to account for the more efficient radiative transfer in metal-rich gas clouds.Stellar particles are assumed to follow a Chabrier IMF (Chabrier 2003), with stars spanning from 0.1-100 M ⊙ .EAGLE accounts for stellar winds, radiation and SNe feedback by implementing their collective feedback through thermal heating, distributing the energy produced at each timestep by a stellar particle to the neighbouring particles using the stochastic thermal feedback scheme of Dalla Vecchia & Schaye (2012).
The radiative cooling and heating rates of gas resolution elements at a given density, temperature and redshift are computed with the software CLOUDY (Ferland et al. 1998).The gas is assumed to be optically thin, in an ionisation equilibrium and exposed to the cosmic microwave background and a spatially homogeneous UV/X-ray background that evolves over time (Haardt & Madau 2001;Wiersma et al. 2009).EAGLE implements one mode of AGN feedback through stochastic heating, following the prescriptions of Dalla Vecchia & Schaye (2008) and Booth & Schaye (2009).

A3 IllustrisTNG
TNG used the moving mesh code AREPO (Springel 2010), which differs from AMR code (e.g.RAMSES) and SPH code (e.g.GADGET-3) in that it uses an unstructured mesh defined by a Voronoi tessellation of a set of discrete points.This alternative methodology was formulated to address issues with AMR and SPH, which impact their accuracy in particular situations (e.g.suppression of fluid instabilities in SPH and lack of Galilean invariance and overmixing present in AMR codes).
TNG models star formation stochastically, treating the star formation and pressurisation of a multi-phase interstellar medium following Springel & Hernquist (2003).In the model by Springel & Hernquist (2003), cold gas above a density threshold of 0.1 H cm −3 forms star particles following the empirically defined Kennicutt-Schmidt relation (Schmidt 1959;Kennicutt 1998) and the Chabrier IMF (Chabrier 2003).Under the Springel & Hernquist (2003) prescription, SNe pressurises the gas and can lead to enhanced star formation.The stellar winds produced by star formation can carry and inject kinetic energy into the surrounding gas.They are modelled as described in Vogelsberger et al. (2013) and Pillepich et al. (2018a).
The radiative cooling of gas is modelled accounting for its metal enrichment (following Wiersma et al. 2009), and a time-evolving homogeneous UV background with self-shielding corrections in dense interstellar medium (following Katz et al. 1992;Faucher-Giguère et al. 2009).Following the models of Vogelsberger et al. (2013), radiative cooling is further influenced by the nearby radiation fields of AGN.
TNG implements subgrid AGN feedback with a radio and quasar mode.For accretion rates below 5% of the Eddington limit, the radio-mode feedback is active, injecting bursty thermal energy into a ∼ 50 pc bubble displaced away from the host galaxy (Sĳacki et al. 2007).For higher accretion rates the quasar mode is active, injecting thermal energy continuously into the adjacent gas (Springel et al. 2005a;Di Matteo et al. 2005) The details of all subgrid physics can be found in Weinberger et al. (2017); Pillepich et al. (2018a) and Nelson et al. (2019).

A4 Magneticum
Magneticum Pathfinder simulations are a suite of cosmological hydrodynamical simulations, ranging in box size from 25.6 3 to 3818 3 cMpc 3 .We use the Box4-uhr simulation, which models the physics of dark matter and baryons in a 68 3 cMpc 3 box.(2013), andBeck et al. (2016).Star formation and the stellar winddriven kinetic feedback are modelled following Springel & Hernquist (2003).Each star particle represents a stellar population following a Chabrier IMF (Chabrier 2003) and each gas particle can form up to 4 stars.Star formation and metal enrichment from supernova feedback and Asymptotic Red Giant Branch stars is modelled following the Springel & Hernquist (2003) as well as the local metallicity dependent processes (Wiersma et al. 2009;Dolag et al. 2017).The radiative cooling and heating rates of gas resolution elements at a given density, temperature and redshift are computed with the software CLOUDY (Ferland et al. 1998).The radiative cooling accounts for a UV/X-ray background following the prescriptions of Haardt & Madau (2001) and Wiersma et al. (2009).
The subgrid physics treatment of supermassive black holes and their AGN feedback is implemented as described by Fabjan et al. (2010) andHirschmann et al. (2014).The implemented black holes feedback scheme accounts for a transition from quasar to radio mode according to (Sĳacki et al. 2007).Note that the black holes in this simulation are not pinned to the potential minimum.Thermal conduction is implemented according to Dolag et al. (2004) but following Arth et al. (2017), with 1/20 of the classical Spitzer value (Spitzer 1962).The details of all subgrid physics are described by Teklu et al. (2015) and Hirschmann et al. (2014).

APPENDIX B: TIDAL FEATURE DETECTION LIMITS
The stellar mass resolutions of our simulations as well as the surface brightness limits of our mock images limits the minimum stellar masses of the tidal features we can resolve.Through some simple modelling we investigate the tidal feature detection limits for each of our simulations.Following Martin et al. (2022) we calculate the signal-to-noise ratio for a tidal feature against its background galaxy and take a SNR tidal > 5 as the range for a detectable tidal feature.
Here,  tidal is the number of particles that comprise a tidal feature,  host is the number of particles from the host in the same region as the tidal feature.Using this equation and the simulation mass resolutions (Table 1) we determined for each of our simulations, given a tidal feature of a specific stellar mass, what was the maximum amount of stellar mass from the host that could be in the same region without making the tidal feature undetectable.The results of this calculation are shown in Figure B1, where we show the maximum host stellar mass contamination in the same region as the tidal feature as a function of tidal feature stellar mass.We can see that each simulation has a drop off corresponding to where its mass resolution is limiting its ability to resolve tidal features with SNR> 5 without there being virtually no contamination from the host galaxy.From this calculation we can predict that NewHorizon is able to resolve tidal features that are ≳ 2 orders of magnitude less massive than in the other three simulations.NewHorizon should also be able to resolve tidal features with much greater contamination from the host galaxy when compared to the other three simulations.
To investigate the minimum mass galaxy that can be tidally disrupted and detected in our mock images, we must construct a model for the surface brightnesses of tidal features stemming from different stellar mass hosts.We build a simple model assuming the following: • the mass-to-light ratio for the stellar particles is 1, • the light from a tidal feature is homogeneously distributed across a square region.
While mass-to-light ratios vary as a function of galaxy morphology (e.g. de Zeeuw et al. 2002;Emsellem et al. 2004;Cappellari et al. 2007) and could vary significantly across a galaxy (e.g.Mehrgan et al. 2024), 1 is a physical mass-to-light ratio and eases computation substantially.The homogeneous distribution of light across a square area is a simplifying assumption that allows us to compute the surface brightness of tidal features without simulating features of varying light distributions and regions.If the mass-to-light ratio for the tidally disrupted galaxy is above 1 or the distribution of the light is inhomogenous resulting in some regions that are brighter than the average brightness of the tidal feature, the resultant tidal features would be detectable despite having lower stellar mass than our estimates.
We convert from the physical surface brightness of the stellar particles in L ⊙ pc −2 to LSST r-band surface brightness in mag arcsec −2 using the following relationship:   = −2.5 log 10 (SB(L ⊙ pc −2 )) +  ⊙,r + 5 log 10 648000 10 (B2) where we use the estimated absolute AB magnitude of the Sun in the LSST -band,  ⊙,r = 4.64 (Willmer 2018).
The left most panel in Figure B2 shows the surface brightness with tidal feature extent for tidal features of stellar masses 10 6 , 10 7 , 10 8 , 10 9 and 10 10 M ⊙ .The surface brightness profiles in remaining panels of Figure B2 indicate the maximum allowed host brightness within the same area to still allow the tidal feature to be resolved against the host background contamination.We compute these assuming a minimum SNR tidal > 5 and again assuming the host light is distributed evenly in the background.Our results highlight that even tidal features from galaxies that are 10 6 M ⊙ are above the mock image surface brightness limits provided their radius remains within ∼ 6 pkpc.However, only NewHorizon has the sufficient resolution to resolve them so long as the light from the host galaxy in the same region as the tidal feature is of a similar brightness or fainter.The surface brightness profile for the maximum allowed host contamination indicates that the these features are only detectable if there is a similar amount of light from the host in this region.From Figure B1 we can estimate that NewHorizon will resolve a 10 6 M ⊙ feature within the survey detection limits if there is ≲ 1.5 × 10 6 M ⊙ of stellar mass from the host within the same region, provided the feature is no more extended than ∼ 6 pkpc.For the larger cosmological simulations, we can detect tidal features of stellar mass ≳ 10 8 M ⊙ even if they have been disrupted over extents ≳ 50 pkpc, so long as the host stellar contamination is ≲ 1.5 × 10 8 M ⊙ within the same region.
As the tidal features become more massive than these lower limits, the maximum surface brightness of the host galaxy within the same region as the tidal feature increases for each of the simulations.

Figure 1 .
Figure1.Examples of gri-colour mock images from each simulation in our sample.Each of the tidal feature morphologies are illustrated from left to right, top to bottom: stream/tail, shell, asymmetric halo, double nuclei and featureless.The galaxy marked by a red circle in the bottom right image is an example of one that we removed from the sample.While there is enough stellar mass within the 30 pkpc aperture for the object to be in our sample, the object itself is too dense and compact to be an accurate representation of a galaxy of  ★,30 pkpcs ≥ 10 9.5 M ⊙ .Each image is 124 pkpc × 124 pkpc (∼ 2400 pixels × 2400 pixels).

Figure 2 .
Figure 2. The galaxy stellar mass-halo mass relation for all central galaxies within our stellar mass range for each of our simulations, the sample sizes are given in the bottom right corner.The relation is plotted using 25 equally spaced bins for 9.5 ≤ log 10 (  ★ /M ⊙ ) ≤ 12.6 and 10.3 ≤ log 10 (  200, crit /M ⊙ ) ≤ 14.7.The dashed and dotted lines show observationally-derived relations from the literature for redshifts comparable to the simulation snapshots, the dotted lines correspond to  = 0.26 (NewHorizon snapshot redshift) and the dashed lines correspond to  = 0.05 (EAGLE, TNG and Magneticum snapshot redshifts).The median relation for central galaxies in each simulation is illustrated with an amber line.On average we see that the relations for EAGLE, TNG and Magneticum are similar, however, NewHorizon is missing high mass halos and contains more stars for a given halo mass than the other simulations.Note for this Figure we present the AdaptaHOP (for NewHorizon) and SubFind (for EAGLE, TNG and Magneticum) based galaxy stellar masses ( ★ ) for the galaxies.The counts have been normalised such that the volume under the histogram integrates to 1.

Figure 3 .Figure 4 .
Figure3.Stellar mass distributions of parent (grey) and selected samples (NewHorizon -black, EAGLE -blue, TNG -red and Magneticum -cyan) for each simulation.The  ★ ≥ 10 9.5 M ⊙ range explored by our study is well sampled by the parent distribution for TNG, EAGLE, and Magneticum, but the smaller volume NewHorizon simulation has fewer galaxies at the high mass end.We include the results of the KS tests comparing the sample to the parent distributions.EAGLE roughly samples its parent distribution and TNG and Magneticum have been sampled to match EAGLE.The KS test results show that TNG remains similar to its parent, whereas Magneticum no longer resembles the parent distribution.Note for this Figure we present the AdaptaHOP (for NewHorizon) and SubFind (for EAGLE, TNG and Magneticum) based galaxy stellar masses ( ★ ) for the parent and subsampled galaxies.

Figure 5 .
Figure 5.The frequencies of the tidal feature morphologies for each simulation: NewHorizon in grey, EAGLE in blue, TNG in red and Magneticum in cyan.From left to right the panels show the results for each confidence threshold:   .≥ 1,   .≥ 2 and   .= 3 (described in Table2).The fractions are given for 5 feature categories: total tidal feature fraction (Total), stream/tail (S/T), shell, asymmetric halo (Asym.)and double nucleus (DN).The error bars give the 1 binomial confidence levels.The occurrence of features relative to one another is similar across each confidence level, the fractions of galaxies exhibiting features drop as the confidence level increases.Generally, the simulations are in good agreement, with NewHorizon, tending to have a higher fraction tidal features that the other simulations.

Figure 6 .
Figure6.Distributions of the fraction of galaxies exhibiting at least one instance of a tidal feature over the range of stellar masses.From left to right, top to bottom we have the fractions of galaxies hosting: streams/tails, shells, asymmetric halos, and a double nucleus.The fractions for NewHorizon are shown in black, EAGLE in blue, TNG in red, and Magneticum in cyan.The error bars indicate the 1 binomial errors.For each simulation, each mass bin contains a similar number of galaxies,  bin,TNG =  bin,EAGLE =  bin,Magneticum = 130 and  bin,NewHorizon ≃ 9.The distributions of tidal features as a function of stellar mass agree well across all simulations, and exceptionally well between the stellar mass-matched samples from TNG, EAGLE and Magneticum.

Figure 7 .
Figure7.Distributions of the total tidal feature fractions as a function of halo mass.The left panel shows the total tidal feature fractions for each of our simulations and the right panel shows the total tidal feature fractions for the central and satellite galaxy populations.The results for NewHorizon are given in black, EAGLE in blue, TNG in red, and Magneticum in cyan.The NewHorizon data is given for the sample of 58 classified galaxies with assigned parent halos and each point represents the feature fraction for a halo mass bin containing ≃ 12 galaxies.The points for EAGLE, TNG and Magneticum represent the mean of 100 Monte Carlo iterations, where each iteration subsampled 450 galaxies to the flat log halo mass distributions depicted in Figure4.Each of the bins for EAGLE, TNG and Magneticum contain 90 galaxies, these galaxies are a varying mix of central and satellite galaxies as can be seen in Figure4.The error bars show the mean 1 binomial errors.

Figure 8 .
Figure8.Distributions of tidal feature fractions as a function of halo mass.From left to right, top to bottom, we have the: stream/tail fraction, shell fraction, asymmetric halo fraction, and double nucleus fraction.The results for NewHorizon are given in black, EAGLE in blue, TNG in red, and Magneticum in cyan.The NewHorizon data is given for the sample of 58 classified galaxies with assigned parent halos and each point represents the feature fraction for a halo mass bin containing ≃ 12 galaxies.The points for EAGLE, TNG and Magneticum represent the mean of 100 Monte Carlo iterations, where each iteration subsampled 450 galaxies to the flat log halo mass distributions depicted in Figure4.Each of the bins for EAGLE, TNG and Magneticum contain 90 galaxies.The error bars show the mean 1 binomial errors.We see excellent agreement between the four simulations plotted.

Figure 9
Figure9.Distribution of tidal feature fractions for satellite and central galaxies as a function of halo mass.From left to right, top to bottom, we have the: stream/tail, shell, asymmetric halo, and double nucleus fractions.The results for NewHorizon are given in black, EAGLE in blue, TNG in red, and Magneticum in cyan.The NewHorizon data is given for the sample of 58 classified galaxies with assigned parent halos.The dashed lines show the feature fractions for satellite galaxies and the circles show the feature fractions for centrals for each halo mass bin containing ≃ 10 galaxies.For the Monte Carlo sampled distributions the points for centrals and satellites show the means of 100 iterations, each bin contains 90 galaxies.Each bin contains a varying mix of central and satellite galaxies as can be seen in Figure4.The error bars for NewHorizon shows the 1 binomial errors and the error bars for EAGLE, TNG and Magneticum data show the mean 1 binomial errors.
Figure 11.Total fraction of galaxies exhibiting tidal features from various observations and simulations studies compared to the results presented here.On the left, we display our own simulation results and simulation results from the literature.On the right, we give the results of observational studies.The error bars along the stellar mass axis show the range of stellar masses for each measurement, with the marker situated on the mean where the sample mean stellar mass is provided (where unavailable the midpoint of the stellar mass range is used instead).Our   .≥ 1 results are presented with a transparent marker and our   .= 3 results are presented with the solid marker.We use a similar opacity to indicate the different confidence level derived fractions from observations.The dashed line links our TNG   .≥ 1 shell fraction with an estimated fraction when linearly extrapolated to higher stellar masses.We find our tidal feature fractions to be in the expected range for observational results and simulation results.

Figure 12 .
Figure12.Investigating whether trends of tidal feature fraction with halo mass and as a function of central and satellite galaxy populations can be explained by the correlation between tidal feature fraction and galaxy stellar mass.The top left panel shows the tidal feature fraction-stellar mass relation shown in Figure10with the inclusion of asymmetric halos and with the final bin for EAGLE, TNG and Magneticum divided into two bins, linearly spaced in log-stellar mass.The top right panel uses identical binning to the top left, and plots the fractions of centrals and satellite galaxies exhibiting tidal features as a function of stellar mass.This illustrates that there are no differences in the relations followed by central or satellite galaxies across all simulations.The bottom left panel is identical to the right panel of Figure7and shows the tidal feature fraction as a function of halo mass for central and satellite galaxies.The bottom right panel plots the tidal feature fraction estimated from the tidal feature fraction-stellar mass relation shown in the top left panel, for the populations of centrals and satellites as a function of mean halo mass.For central galaxies we find reasonable agreement between the bottom two panels, indicating that the trends seen in the tidal feature fraction-halo mass relation are a result of the relationship of the tidal feature fractions with stellar mass.However, satellites in the bottom right panel have fractions that systematically increase with increasing halo mass compared to the bottom left panel, where they increase and then fall as a function of increasing halo mass.This could indicate that there are factors other than stellar mass that drive tidal feature fractions for satellite galaxies.The uncertainties in the tidal feature fraction are given by the 1 binomial errors for the relations with stellar mass and the mean 1 binomial errors for the relations with halo mass.The uncertainties on the mean halo mass are given by standard deviation of the mean halo mass over the 100 Monte Carlo iterations.

Figure B1 .
Figure B1.The maximum stellar mass from the host galaxy that in the same region as the tidal feature without making it undetectable as a function of the tidal feature stellar mass.For a tidal feature to be detectable it has SNR> 5.

Figure B2 .
Figure B2.The -band surface brightness of the model tidal features as a function of their radii and the corresponding maximum possible surface brightness of the light from the host galaxy within the tidal features area such that the tidal feature is still detectable.The tidal features panel shows the surface brightness of each tidal feature as a function of the radius of the circle the light is spread over.The remaining panels show the maximum surface brightness of the light from the host that can be within the same region as the tidal feature without making it undetectable for each simulation.The colours correspond to the stellar masses of the tidal features, 10 6 , 10 7 , 10 8 , 10 9 and 10 10 M ⊙ .The LSST 10 year -band, 3, 10 ′′ × 10 ′′ , limiting surface brightness is shown in Rubin Turquoise.

Table 4 .
The mean galaxy stellar masses for each of the tidal feature categories, from top to bottom: all tidal features, stream/tail, shell, asymmetric halo and double nucleus, for the four simulations, from left to right: NewHorizon, EAGLE, TNG and Magneticum.The final column shows the overall mean over the four simulations for the tidal feature host mass.The uncertainties show the standard deviations on the means.* NewHorizon is excluded from the calculation of the average shell host mass as there is only one feature in that sample. ★, 30 pkpc /M ⊙ )] Martin et al. (2022)ellar mass of the total tidal feature fraction (excluding asymmetric halos) presented with results fromMartin et al. (2022).Stellar mass is measured by  ★, 30 pkpc for our sample and Adap-taHOP stellar masses forMartin et al. (2022).NewHorizon is shown in black, EAGLE in blue, TNG in red and Magneticum in cyan.TheMartin et al. (2022)results are shown in orange.We compute the fractions using classifications with   ≥ 2. The error bars indicate the 1 binomial uncertainty for our results and the 1 uncertainty obtained from 100 000 bootstraps for