\section{Introduction}

Our understanding of the Universe is closely tied to our capabilities to observe it, and surface brightness limits create severe biases that must be taken into account in order to provide a correct representation of the low-mass galactic population \citep{1976Natur.263..573D, 1995AJ....110..573M}. While the new generation {of} deep optical surveys are providing a growing number of previously undetected low surface brightness galaxies \citep[e.g.,][]{2018ApJ...857..104G, 2019ApJS..240....1Z, 2020arXiv200604294T}, the surface brightnesses achieved by integrated photometry are far from the limits achieved by star counting in the local Universe, where galaxies of up to 32 mag arcsec$^{-2}$ have been discovered \citep{2019PASJ...71...94H}. There is, therefore, a huge population of currently undetectable low-mass galaxies beyond the distance at which we can resolve individual stars.

An example of the impact of surface brightness limitations is found in the search for the so-called dark galaxies, that is, dark matter halos that failed to convert their gas into stars. They are expected from a theoretical perspective because of the low star formation efficiency of low mass halos at low metallicity \citep{2010ApJ...714..287G, 2012ApJ...753...16K}. Objects compatible with dark galaxies have been detected at high redshifts in Ly$\alpha$ fluorescence by being illuminated by quasars \citep{2012MNRAS.425.1992C, 2018ApJ...859...53M, 2019ApJ...875..130L}. However, the systematic search for dark galaxies in the local Universe is strongly biased by the prevalence of H\,{\sc i} debris expelled from galactic interactions that can masquerade as gas within a dark matter halo that failed to form stars, {that is, a genuine dark galaxy} \citep[see][]{2005MNRAS.363L..21B, 2008ApJ...673..787D}. In fact, H\,{\sc i} clouds with no optical counterpart are common in high-density environments \citep[e.g.,][]{2012MNRAS.423..787T,2013A&A...555L...7O,2017MNRAS.464..530S,2020A&A...642L..10B}, but systematic searches in low-density environments (where the possibility of H\,{\sc i} tidal debris can be minimized) have not returned any clear candidates for true dark galaxies \citep{2015AJ....149...72C}, calling into question the existence of bona fide dark galaxies.

{Gas-rich} tidal interactions are commonly accompanied by bursts of star formation \citep{2001ApJ...550..204I, 2005A&A...430..443V, 2012MNRAS.426.2441D, 2014MNRAS.443.3601L} due to local instabilities of the gas \citep{2007MNRAS.375..805W}, forming the so-called tidal dwarf galaxies \citep[TDGs; ][]{1998A&A...333..813D}. Due to the difficulty in obtaining {a} high resolution surface density of the gas needed to establish star formation thresholds in these systems \citep{2004ApJ...609..667S}, the potential relation between tidal interactions and H\,{\sc i} clouds without an optical component is hard to study observationally. Additionally, effects such as fading luminosity due to the aging of stellar populations have to be properly taken into account in order to assert the presence or absence of a stellar component.

\begin{figure*}
    \centering
        \includegraphics[width=0.97\textwidth]{Panel.pdf}
    \caption{Illustration of the structures present in the northwestern region of HCG16. Upper left panel: Color-composed image using the \textit{g, r}, and \textit{z} optical bands from DECaLS. Lower left panel: Pseudo zenith view (R.A. vs. H\,{\sc i} velocity) of the H\,{\sc i} emission higher than 3$\sigma_\mathrm{rms}$. {Integrated moment-zero and first moment maps as well as detailed information about these data and maps for the whole group can be found in \cite{2019A&A...632A..78J}}. The dashed blue line in both left panels shows the trajectory of the H\,{\sc i} tail, the end of which coincides with the diffuse object detected in the optical. {We note that the progenitor of HCG 16-LSB1 is most likely HCG 16c and not the apparently closest galaxy (in projection), HCG 16a}. Right panels: Zoom-in of the object and its fitting using optical bands through a S\'ersic model.}
    \label{fig:Panel}
\end{figure*}

In this letter we study the properties of a peculiar object that was identified by \cite{2019A&A...632A..78J} as a very faint feature detected in optical data {and} located at the end of an H\,{\sc i} tail emerging from the complex interactions of Hickson Compact Group 16 (HCG 16), which \cite{2019A&A...632A..78J} suggest is a clear TDG candidate. For convenience, we have named this object HCG 16-LSB1. We aim to test its TDG nature, study its properties, and analyze its potential future evolution. We assume standard cosmology parameters, with a distance to the object of 54.7 Mpc and a spatial scale of 0.259 kpc arcsec$^{-1}$. We use the AB photometric magnitude system in this work.

\section{Properties of HCG 16-LSB1}

We used publicly available photometric data from the seventh {data} release of the Dark Energy Camera Legacy Survey  \citep[DECaLS; ][]{2016MNRAS.460.1270D} to obtain the photometric and structural properties of HCG 16-LSB1 in the optical. We performed a S\'ersic model fitting using \texttt{IMFIT} \citep{2015ApJ...799..226E} in the three available photometric bands (\textit{g}, \textit{r}, and \textit{z}). The parameters resulting from this modeling are shown in Table \ref{tab:phot}. These properties match with those of diffuse galaxies frequently found in groups of galaxies \citep[e.g.,~][]{2017MNRAS.468.4039R}, with $r_{e}$~=~1.5~kpc, $\mu(0)_{g}$~=~25.2~mag~arcsec$^{-2}$, and $n$~=~0.48. The colors of HCG 16-LSB1 are extremely blue, with \textit{g-r}~=~0.00~$\pm$~0.10 mag, suggesting very young stellar populations. In the right panels of Fig. \ref{fig:Panel} we plot the color composite model of the object together with its residuals. As can be seen, a S\'ersic model is a good representation of the object with a clean residual image.

We further explored other optical data sets. We found observations with the Canada-France-Hawaii Telescope (CFHT) in the CFHT Science Archive for the \textit{g}, \textit{r}, and \textit{i} optical bands. However, the presence of an internal reflection from a nearby star prevents photometry of the object in these images. We also found observations with the Hubble Space Telescope (HST) in the Hubble Legacy Archive (PropID = 10787) for the \textit{F814W}, \textit{F606W}, and \textit{F450W} bands. In this case, we ruled out the use of these images due to the presence of residuals in the reduction process that have a strong impact on low surface brightness objects in the HST images; these images would thus require complex data processing in order to obtain reliability \citep{2019A&A...621A.133B}. While photometric or structural measurements in these additional optical data are unfeasible, their high resolution allows us to explore in great detail the possible presence of HII regions. However, we do not find any such regions within the depth and resolution limits that these data provide: HCG 16-LSB1 appears to be a smooth and featureless galaxy.

\subsection{SED fitting}

We conducted a study of the stellar populations of HCG 16-LSB1 with the available photometric information in this region, aiming to constrain its age and metallicity. In addition to the aforementioned data in the \textit{g}, \textit{r}, and \textit{z} optical bands, we also added the ultraviolet data from the Galaxy Evolution Explorer \citep[GALEX;][]{2005ApJ...619L...1M} and infrared archival data from Spitzer \citep{2004ApJS..154....1W}. The photometry was carried out by fitting a S\'ersic model, in a similar way to what was done in the optical bands. The resulting photometry is given in Table \ref{tab:phot}. In the case of Spitzer data, HCG 16-LSB1 was not detected, so upper limits were estimated.

\begin{table}
\begin{center}
\begin{tabular}{cc}
Parameter & Value \\ \hline
R.A. & 32.3587\\
Dec. & -10.1194\\
r$_e$ & 1.5 $\pm$ 0.1 kpc\\
$<\mu _g>$ & 25.6 $\pm$ 0.1 mag/arcsec$^2$ \\
$\mu _g$(0) & 25.2 $\pm$ 0.1 mag/arcsec$^2$ \\
S\'ersic index & 0.48 $\pm$ 0.01\\
Axis ratio & 0.96 $\pm$ 0.01 \\
M$_{g}$ & $-$14.05 $\pm$ 0.04 mag\\
H\,{\sc i} velocity & 3552~$\pm$~2~km~s$^{-1}$\\
\\
Far-UV band & 20.24 $\pm$ 0.45 mag\\
Near-UV band & 20.30 $\pm$ 0.15 mag\\
g-band & 19.64 $\pm$ 0.04 mag\\
r-band & 19.64 $\pm$ 0.06 mag\\
z-band & 19.94 $\pm$ 0.20 mag\\
IRAC1$_{3.6 \mu m}$ band & < 20.19 mag\\
IRAC2$_{4.5 \mu m}$ band & < 20.19 mag\\
IRAC3$_{5.8 \mu m}$ band & < 18.68 mag\\
IRAC4$_{8.0 \mu m}$ band & < 18.68 mag\\
\newline
\end{tabular}
\caption{Summary of the properties of HCG 16-LSB1. All photometric values are corrected from extinction \citep{2011ApJ...737..103S}.}
\label{tab:phot}
\end{center}
\end{table}

\begin{figure}
    \centering
        \includegraphics[width=0.8\columnwidth]{SED.pdf}
    \caption{Spectral energy distribution of the galaxy. The filled circles are the observed magnitudes of the galaxy, while the inverted triangles correspond to the upper limits in those bands where the galaxy was not detected. The open circles are the expected flux after convolving the best fit model with the filter responses.}
    \label{fig:SED}
\end{figure}

As HCG 16-LSB1 shows no structure or "clumpiness" and has an expected low age, we assume that it formed in a single burst of star formation. Therefore, characterizing it as a single stellar population (SSP) with only a few parameters, such as metallicity, age, and initial mass function (IMF), is a good approximation. We compared the observed spectral energy distribution (SED) of HCG 16-LSB1 with SSP models \citep{2003MNRAS.344.1000B}, {using} a \citet{2003ApJ...586L.133C} IMF to obtain estimates of the age and metallicity of the galaxy. We used a $\chi^2$-minimization approach to obtain the best fit model to our data \citep[as described in][]{2014MNRAS.439..990M}. The uncertainties were estimated by marginalizing the 1D probability distribution functions. Figure \ref{fig:SED} shows the best fitting model to the SED of HCG 16-LSB1 with 57.9$^{+21.7}_{-9.4}$~Myr and [Fe/H]~=~$-$0.16$^{+0.43}_{-0.41}$ (solid blue line). All photometric points fit well with the model except for the one corresponding to the \textit{z} band. We think that this could be due to background sky fluctuations that are stronger in this band, producing a higher uncertainty due to the very low surface brightness of HCG 16-LSB1. The stellar-mass-to-light ratio is M/L$_{V-band}$~=~0.073$^{+0.097}_{-0.059}$ $\Upsilon_{\sun}$, giving a stellar mass of M$_\star$~=~4.2$^{+5.6}_{-3.4}\times$10$^6$~M$_\odot$ for HCG 16-LSB1.

The metallicity of HCG 16-LSB1 of [Fe/H]~=~$-$0.16$^{+0.43}_{-0.41}$ is much higher than what would be expected for a dwarf galaxy of its stellar mass, which is typically within the range [Fe/H]~=~[-2.0,~-1.5]  \citep{2003AJ....125.1926G, 2013ApJ...779..102K, 2017A&A...606A.115H}. This suggests that it formed from the pre-enriched gas of the galactic components of HCG 16, this high metallicity content being a fundamental characteristic defining TDGs. \citep[see a review by][]{2012ASSP...28..305D}.

\subsection{Dynamics and survival}
\label{sec:dynamic}

In the upper panels of Fig. \ref{fig:HI_spec} we show the H\,{\sc i} emission for different velocity channels corresponding to the tip of the H\,{\sc i} tail shown in Fig. \ref{fig:Panel}. {These data were obtained in the Very Large Array using configurations C and D and were processed by \cite{2019A&A...632A..78J}; more detailed information regarding the process is available in that work.} We plot the emission contours of the robust=0 (19.4\arcsec \ $\times$ 14.8\arcsec) and robust=2 (37.2\arcsec \ $\times$ 30.3\arcsec) cubes, with an available spectral resolution of $\approx$ 21 km s$^{-1}$.

\begin{figure}
    \centering
        \includegraphics[width=1\columnwidth]{Panel_HI.pdf}
    \caption{{H\,{\sc i} distribution of HCG 16-LSB1.} Upper panels: Velocity channel maps corresponding to high (light blue) and low (black) resolutions of the H\,{\sc i} emission plotted over a high contrast $g+r$ optical image. The beam size for both resolutions is placed in the upper left corner (see text). We note that the H\,{\sc i} emission in the plotted velocity channels is not associated with the bright galaxies (HCG 16a and HCG 16b) south of HCG 16-LSB1 (see Fig. \ref{fig:Panel}). Bottom panel: H\,{\sc i} spectrum obtained from the channels (high resolution data) presented in the upper panel (i.e., integrated over HCG 16-LSB1; see Sect. \ref{sec:dynamic} for details). The peak emission is fitted with a Gaussian function.}
    \label{fig:HI_spec}
\end{figure}

The velocity channels show spatial matching between the tip of the H\,{\sc i} tail and the diffuse emission in optical bands. Considering the young age and high metallicity of the stellar component, which can only be explained by a recent formation with pre-enriched material, and the fact that TDGs are preferentially formed at the tip of H\,{\sc i} tails, a chance superposition is extremely unlikely. This confirms the association between the H\,{\sc i} cloud located at the end of the H\,{\sc i} tail and the stellar component detected in optical bands. However, we find an offset between the peak emission in H\,{\sc i} and the center of the object in the optical. This offset is 5.5'' for the 3562 km s$^{-1}$ channel, which is approximately the value of the effective radius (1.5 kpc). While the association between the H\,{\sc i} tail and the stellar component seems clear, this potential shift could indicate a certain degree of gas stripping or a spatial asymmetry between both components that cannot be resolved with the resolution available in our data.

In the lower panel of Fig. \ref{fig:HI_spec} we show the integrated H\,{\sc i} spectrum. It was extracted using an aperture with dimensions twice  that of the synthesized beam of the robust=0 cube (i.e., four times the area). The robust=0 cube was favored over the robust=2 cube due to its better angular resolution; however, we note that this is at the cost of sensitivity, and thus we retain the initial H\,{\sc i} mass estimation from \cite{2019A&A...632A..78J} of M(HI)~=~3.9$\times$10$^8$~M$_\odot$, which used the robust=2 cube. We fitted the peak associated with the tip of the tail with a Gaussian function for the 3520, 3541, 3562, and 3583 km s$^{-1}$ channels, which have a spatial match with the optical counterpart. The residuals of the Gaussian fits are minimal up to 3600 km s$^{-1}$. From this velocity onward, an excess flux remains; this flux corresponds to the residual flux of the rest of the H\,{\sc i} tail, which shifts spatially to a southeastern direction toward HCG16-c (see Fig. \ref{fig:Panel}). The model parameters of the fitted Gaussian distribution are $\sigma$~=~24.2~$\pm$~2.8~km~s$^{-1}$ centered at 3552~$\pm$~2~km~s$^{-1}$.

In order to obtain a dynamic mass estimation, we used the approach from \cite{1996ApJS..105..269H}:

\begin{equation}
    M_{dyn} = 2.325\times10^5 \left( \frac{V_{rot}^2 + 3\sigma^2}{km^2\:s^{-2}}\right) \left(\frac{r}{kpc}\right)M_\odot.
\end{equation}

This equation is only valid in the case of dynamic equilibrium and self-gravitation. Some clues indicating that this could be the case here are the Gaussian shape of the H\,{\sc i} peak and the symmetry and sphericity of the optical counterpart, evidenced by the absence of residuals in the S\'ersic modeling, which contrasts with the complex morphology of similar galaxies in clear interaction \citep{2021MNRAS.tmp..419S}. However, a source of uncertainty is the low spectral resolution (21 km s$^{-1}$), which does not allow us to resolve dynamical substructures within the tip of the H\,{\sc i} tail. Taking these considerations into account, we proceeded to obtain this estimate.

For the radius we used a range between a minimum value of two times the effective radius, which is 3 kpc, and a maximum value corresponding to the selected aperture in which the H\,{\sc i} spectrum was extracted, 40 arcsec and approximately 10 kpc at the group distance. These values produced a range of possible dynamic masses derived from the low resolution detection of H\,{\sc i}. Due to the apparent absence of a rotation pattern within the resolution limits \citep[see also][]{ 2019A&A...632A..78J}, we assumed that the dynamic mass is pressure supported, producing a value interval of M$_{dyn}$~=~[1.2, 4.1]~$\times$~10$^9$~M$_\odot$; this value is in contrast with the baryonic content, which is mostly gas, of M(HI)~=~3.9~$\times$~10$^8$~M$_\odot$. This would imply a mass-to-light ratio in the interval M/L~=~[3,~10]. It is worth commenting that even assuming a purely rotating unresolved component and a very conservative inclination of $i$~=~45$\degree$, the derived mass-to-light ratio would be in the range M/L~=~[1.5,~5]. Therefore, although there are caveats due to the low spectral resolution of the H\,{\sc i}, our analysis suggests that HCG 16-LSB1 contains a significant amount of dark matter; this is in contrast to what is expected for TDGs, which are typically poor in dark matter content \citep[e.g.,][]{ 2012ASSP...28..305D}. It is also worth commenting that even exclusively considering the observed baryonic matter of M(HI)~=~3.9~$\times$~10$^8$~M$_\odot$, it would be higher than the minimum mass for a long-lived TDG, which is approximately M~=~10$^8$~M$_\odot$ \citep{2006A&A...456..481B}. This would make the system gravitationally bound in future evolution, even assuming a total absence of dark matter.

\begin{figure*}
    \centering
        \includegraphics[width=0.75\textwidth]{LSB_evo.pdf}
    \caption{Evolution of the surface brightness of HCG 16-LSB1 in time. The black diamond symbol marks the observed age and surface brightness. The tracks in different colors mark the theoretical temporal evolution of the surface brightness that would correspond to different metallicities. The dashed black line corresponds to the evolution calculated for the metallicity of HCG 16-LSB1 ([Fe/H] = -0.16). Approximate limits on the maximum surface brightness of detectable galaxies are marked with dotted lines for the Sloan Digital Sky Survey \citep[SDSS; ][]{2000AJ....120.1579Y}, DECaLS \citep{2016MNRAS.460.1270D}, the IAC Stripe82 Legacy Survey \citep{2016MNRAS.456.1359F,2018RNAAS...2..144R}, the HSC-SSP \citep{2018PASJ...70S...4A}, and the future LSST \citep{2009arXiv0912.0201L}.}
    \label{fig:LSB_evo}
\end{figure*}

Although dynamic stability by self-gravitation seems assured, interactions with massive galaxies could easily tear the galaxy apart, which is often the case in compact groups. However, tidal features tend to eject at high speeds, being able not only to perform a distant orbit but even being able to flee the system by exceeding the escape velocity. HCG 16-LSB1 seems to have the latter ability. The escape velocity in HCG16 was estimated by \cite{2019A&A...632A..78J} to be 400~km~s$^{-1}$, close to the radial velocity of HCG 16-LSB1 with respect to the average of the host group, 300~km~s$^{-1}$. We can make a rough estimate of the tangential component of the velocity. Assuming (i) a flight time equal to the age of HCG 16-LSB1 and (ii) a distance traveled equal to half the projected distance between HCG 16-LSB1 and HCG16c (see Fig. \ref{fig:Panel}) -- considering its formation not in the vicinity of HCG16c but once the H\,{\sc i} tidal tail has already been ejected, which has been observed in other TDGs \citep{2014MNRAS.443.3601L} -- results in a tangential component of ${26\,kpc}/{58\,Myr}$ $\approx$ 450~km~s$^{-1}$. This is of course an estimate with large uncertainties, but it is useful to show that the tangential component of the velocity is significant. Therefore, the 3D component of the velocity is expected to be higher than the escape velocity of 400~km~s$^{-1}$, and therefore HCG 16-LSB1 will likely escape from the gravitational influence of the group and its interactions that could destroy it.

\subsection{Evolution in time of the optical surface brightness}

Since HCG 16-LSB1 will probably survive, it is interesting to explore how its photometric properties in the optical will evolve in the future. For that we can consider it as an SSP that is passively evolving, which is a good approximation since no current star formation is detected, the H\,{\sc i} peak is diffuse enough to not trigger new star formation \citep[2.45~$\times$~10$^{20}$~cm$^{-2}$; see][]{2019A&A...632A..78J}, and the accretion of new gas to the system is unlikely. For this we used the stellar population models from \cite{2015MNRAS.449.1177V} with a Kroupa universal IMF \citep{2001MNRAS.322..231K}. In Fig. \ref{fig:LSB_evo} we plot the photometric tracks corresponding to a temporal evolution, starting at the observed age and surface brightness of HCG 16-LSB1. We also plot photometric tracks for a wide range of metallicities, from [Fe/H] = $-$2.3 to $+$0.2, for comparison. Typical surface brightness limits for detecting low surface brightness galaxies are marked for the different surveys. We used as a reference the work by \cite{2018ApJ...857..104G}, who used the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) with a detection limit for the mean surface brightness of approximately 28 mag arcsec$^{-2}$ in the \textit{g} band. We scaled this value using the limits derived from the Poisson noise through the relation given by \citet[][Appendix A]{2020A&A...644A..42R}, taking the various surveys' differing telescope apertures and exposure times into account. These calculated limits are an approximation since they not only depend on the depth of the data but also on the sizes of the sources and the method applied in the detection.

At the young age of HCG 16-LSB1, the stellar dimming due to the aging of the stellar populations is very strong. We can see that in just 400 Myr, the surface brightness will decrease enough to be undetectable in the optical data that we use in this work (DECaLS). In 1 Gyr of evolution, it would be undetectable in the deepest currently available survey, which is the HSC-SSP. And in 2 Gyr it would be undetectable in what will be the deepest optical survey in the near future, the Legacy Survey of Space and Time (LSST), reaching 29 mag arcsec$^{-2}$, which is the surface brightness of the faintest galaxy ever detected through integrated light \citep{2018ApJ...863L...7M}. After 10 Gyr of passive evolution, the stellar dimming is extreme, with a decrease in surface brightness of 5 mag arcsec$^{-2}$, making HCG 16-LSB1 completely undetectable, with an effective surface brightness above 30 mag arcsec$^{-2}$. It is interesting to note that the higher the metallicity, the stronger the dimming of the stellar populations. This means that in the case of TDGs, whose metallicity tends to be very high, this effect has a stronger impact on their detectability than that of ordinary dwarfs, whose metallicities are lower.

\section{Discussion}

The properties of HCG 16-LSB1 are peculiar. Its morphology and low surface brightness place it in the category of ultra diffuse galaxies {\citep[UDGs;][]{2015ApJ...798L..45V}}, that is, galaxies with large effective radii and low surface brightnesses {found in any environment \citep[e.g.,][]{2017MNRAS.468.4039R, 2017ApJ...842..133L, 2018MNRAS.481.4381M, 2019MNRAS.486..823R}. However,} its high metallicity and its location at the end of an H\,{\sc i} tail makes HCG 16-LSB1 compatible with a recently formed TDG. Strikingly, systematic detections of TDGs show compact-clumpy morphologies with bright surface brightnesses \citep[e.g.,][]{2020MNRAS.499.3399R}, contrary to that of HCG 16-LSB1. The dark matter content of HCG 16-LSB1 also appears anomalous for a TDG{, with a calculated mass-to-light ratio in the range M/L~=~[3, 10]. This is a standard for the general low-mass galactic population and would be in contradiction with the general characteristics of TDGs, which are generally poor in dark matter \citep[see e.g., ][]{2012ASSP...28..305D}. However, we warn that the spectral (and spatial) resolution that we have is not sufficient to obtain a reliable estimate of the dynamical mass of HCG 16-LSB1, dynamical calculations of such systems being extremely challenging in general \citep{2015A&A...584A.113L}. Future studies with better resolution will be necessary in order to confirm this significant presence of dark matter; however, if confirmed, this would make HCG 16-LSB1 an object of great interest, with {mixed} characteristics of both diffuse galaxies and TDGs.}

Due to the low surface brightness of HCG 16-LSB1, we argue that its peculiar characteristics could be due to an observational bias in which the surface brightness limits of typical optical observations are insufficient to systematically detect such extreme objects. In fact, very deep optical observations focused on the study of TDGs find objects as faint as or even fainter than HCG 16-LSB1 \citep{2014MNRAS.440.1458D}. Therefore, it is expected that the new era of deep optical surveys will reveal a systematic presence of objects of this type, preferentially in strongly interacting galactic associations such as the Hickson Compact Groups.

Indeed, a correct understanding of observational biases is key to providing a correct representation of the low-mass galactic population. Our analysis, exemplified by the case of HCG 16-LSB1, shows that dimming via the passive evolution of stellar populations is another important factor to consider when studying galaxies without optical counterparts \citep[see also][]{ 2017ApJ...836..191T}. The characteristics that give rise to TDGs would be optimal in producing fake dark galaxies. In the first place, the ejection of H\,{\sc i} tails in which TDGs are formed would allow them to escape from the potential well of the host system, avoiding interactions that could destroy them. Additionally, since they are formed by reprocessed gas from massive host galaxies, TDGs are high in metallicity, maximizing the fading of their stellar populations, as clearly shown in Fig. \ref{fig:LSB_evo}.

The connection between H\,{\sc i} debris and optically dark galaxies has been widely debated \citep[e.g.,][]{2005MNRAS.363L..21B, 2008ApJ...673..787D, 2017MNRAS.467.3648T}. While the fraction of H\,{\sc i} debris that actually contains a hidden stellar component is unknown, our work suggests that even {H\,{\sc i} debris galaxies} containing a detectable stellar component at the time of its formation, as is the case with TDGs, {are} able to reproduce the observational properties of optically dark galaxies observed in high-density environments \citep[e.g.,][]{2015AJ....149...72C, 2020A&A...642L..10B} after the passive evolution of their stellar content. This seems to be the case for HCG 16-LSB1. Our analysis shows that its most likely future involves its ejection out of the environment of HCG 16; it will appear as a gas cloud with a stellar population that will quickly fade beyond the optical detection limits. If the galaxy survives and is held together gravitationally -- even if it loses a significant fraction of its gas component, either by tidal interactions taking place in its formation or later on via photoevaporation -- its observational properties will be compatible with a dark matter halo in which only gas is detectable, with no optical counterpart, therefore mimicking the observational properties of dark galaxies without actually being one.

\begin{acknowledgements}
{We thank the anonymous referee for a careful review and useful comments}. We also thank Karen Lee-Waddell for helpful suggestions that improved this work. The authors acknowledge financial support from the grants AYA2015-65973-C3-1-R and RTI2018-096228-B-C31 (MINECO/FEDER, UE), as well as from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709). JR acknowledges support from the State Research Agency (AEI-MCINN) of the Spanish Ministry of Science and Innovation under the grant "The structure and evolution of galaxies and their central regions" with reference PID2019-105602GB-I00/10.13039/501100011033. This research is based in part on observations made with the Galaxy Evolution Explorer, obtained from the MAST data archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555, and is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.
\end{acknowledgements}


\section{Introduction} \label{sec:introduction}
    In a hierarchical model of structure formation, the structures observed today are a result of the merging of dark matter clumps at high redshift via gravitational collapse. As the clumped dark matter grows into larger objects to form haloes, they experience tidal torques from neighbouring objects. At the same time, baryons condense in the centre of such structures to form galaxies \citep{Peebles69,Doroshkevich70,White78,White84}. During their mutual formation both the dark matter and baryonic component experience the same tidal fields, and hence it is expected that they gain the same amount of specific angular momentum. Furthermore, considering that baryons evolve inside dark matter haloes, it is also expected that the kinematic evolution of galaxies is influenced by the halo in which they reside.

    However, galaxies and their haloes also feel the influence of the environment. Haloes falling towards and eventually orbiting within galaxy clusters are disrupted by a series of processes predominant in the cluster environment, e.g. ram-pressure stripping \citep{Gunn72,Abadi99,Bahe15,Arthur19,Mostoghiu21} that removes the gas in haloes and quenches the star formation of galaxies; galaxy harassment \citep{Moore96,Moore98,Smith10,Smith15},  mergers \citep{Dressler80,Hashimoto00,Behroozi14}, tidal torques \citep{Fujita98,Balogh00,Park07}, and interactions in general \citep{Knebe06a,Recchi14} that can disrupt the haloes' components; and dynamical friction \citep{Valtonen90,Jiang00,Fujii06,vanDenBosch17,Miller20} which slows down infalling haloes and, over time, causes them to fall to the centre of the cluster. While the stars residing deep inside the potential well of the halo are shielded from tidal effects, they might nevertheless feel and react to the change of its own halo caused by the aforementioned processes. Or put differently, if the kinematic evolution of galaxies is indeed tied to its halo, the cluster environment affecting the halo could eventually also disrupt the kinematics of galaxies at its centre. This question now lies at the heart of the present study.

    Previous numerical studies have shown how tidal interactions can disrupt infalling haloes in different environments. The seminal work of \citet{Hayashi03} studied how tides influence substructure in dark matter-only simulations. Their analysis showed that, although tides preferentially strip the outer regions of haloes, they also decrease the halo's central density after each pericentric passage. Subsequent studies improved the tidal disruption estimation by including a stellar component to their analysis  \citep[e.g.][]{Bullock05,Penarrubia08b}. The stellar component was found to be exceptionally resilient to tides, preserving its density profile shape even after losing a considerable amount of stars. Nevertheless, the stellar dynamics of such calculations were modelled by analytic profiles, which simplify processes such as the mass loss from dynamical friction, the halo phase-space evolution after merger interactions, or the stellar mass fractions determined by star formation. More recent studies overcame these limitation by introducing a stellar component modelled by full-physics hydrodynamics \citep[e.g.][]{Smith16,Lokas20,Errani20,Mazzarini20}. However, these studies primarily focus on mass-loss processes, hence the question of how environmental effects influence the kinematic properties of the stellar component of infalling haloes remains to be addressed.

    We approach these issues by analysing simulations from \textsc{The Three Hundred} project\footnote{\url{https://the300-project.org}}, i.e. a sample of over 300 galaxy clusters simulated with full-physics hydrodynamics \citep{Cui18}. These simulated clusters have been used for different studies, for example, environmental effects \citep{Wang18}, cluster profiles \citep{Mostoghiu19,Li20,Baxter21}, backsplash galaxies \citep{Arthur19,Haggar20,Knebe20}, cluster dynamical state \citep{Capalbo20, DeLuca20}, filament structures \citep{Kuchner20,Rost21,Kuchner21} and gravitational lensing \citep{VegaFerrero21}. In this work we extend the previous stellar angular momentum analysis to this set of massive simulated galaxy clusters to study the influence of the cluster environment on the stellar kinematics of infalling galaxies.

    This paper is organised as follows. In \Sec{sec:data} we present the data used for the analysis. \Sec{sec:data-subsec:simulations} briefly describes the simulations. In \Sec{sec:data-subsec:sample} to \Sec{sec:data-subsec:stellarcomponent} we define how we selected the sample of galaxies used for the analysis and we present the classification of our objects. \Sec{sec:results} describes our results: in \Sec{sec:results-subsec:kappa} we focus on the kinematic evolution of the sample, and in \Sec{sec:results-subsec:jM} we study their angular momentum-stellar mass relation. Finally, we conclude the study in \Sec{sec:conclusions}.

\section{The Data} \label{sec:data}

    \subsection{`The Three Hundred' Central Galaxy Clusters} \label{sec:data-subsec:simulations}

    \paragraph*{The Simulations}The simulated clusters in \textsc{The Three Hundred} dataset were created by extracting 324 spherical regions of $15 h^{-1}$ Mpc radius centred on each of the most massive haloes identified at $z=0$ within the dark-matter-only MDPL2 simulation \citep{Klypin16}\footnote{The MultiDark simulations -- incl. the MDPL2 used here -- are publicly available at \url{https://www.cosmosim.org}}. MDPL2 was simulated with a \textit{Planck} 2015 cosmology \citep{Planck15}, with $\Omega_{\rm M} = 0.307$, $\Omega_{\rm b} = 0.048$, $\Omega_\Lambda = 0.693$, $h=0.678$, $\sigma_8 = 0.823$, and $n_s=0.96$ and consists of a box of $1 h^{-1}$ Gpc side-length which contains $3840^3$ dark matter particles each of mass $1.5 \times 10^9 h^{-1}$ M$_{\odot}$. In order to model the galaxy clusters with all the relevant baryonic physics, those $15 h^{-1}$ Mpc regions were traced back to the initial conditions and there populated with gas particles by leading to a mass resolution of $m_{\rm DM}=1.27\times 10^{9} h^{-1}$${\ifmmode{{\rm {M_{\odot}}}}\else{${\rm{M_{\odot}}}$}\fi}$ and $m_{\rm gas}=2.36\times 10^{8} h^{-1}$${\ifmmode{{\rm {M_{\odot}}}}\else{${\rm{M_{\odot}}}$}\fi}$, respectively. Outside the re-simulated region, to reduce the computational cost of the original MDPL2 simulation, dark matter particles are degraded with lower mass resolution particles to maintain the same large scale tidal field. Using a Plummer equivalent softening of $6.5h^{-1}$ kpc for both the dark matter and baryonic component, the new initial conditions were now moved forward in time using \code{GADGET-X}  \citep{Beck16}.  \code{GADGET-X} is  a modified version of \code{GADGET3}  with a modern Smooth Particle Hydrodynamics (SPH) scheme which improves the treatment of gas particles \citep{Beck16, Sembolini16b}. Results of simulations of galaxy clusters based on \code{GADGET-X} have been presented in several previous papers \citep[e.g.][]{Rasia15,Planelles17} and in the \textit{nIFTy cluster comparison} project \citep{Sembolini16a,Elahi16,Cui16,Arthur17}. A total of 129 snapshots are saved from $z=16.98$ to $z=0$.

    \paragraph*{The Halo Finding}The halo analysis was done using the \code{AHF}\footnote{\url{http://popia.ft.uam.es/AHF}} halo finder \citep{Gill04a,Knollmann09}. \code{AHF} locates local overdensities in an adaptively smoothed density field as potential halo centres and automatically identifies haloes and substructure (subhaloes, subsubhaloes, etc.). The radius of a halo $R_{200}$ and the corresponding enclosed mass $M_{200}$ are calculated as the radius $r$ at which the cumulative density $\rho(<r)=M(<r)/(4\pi r^{3}/3)$ drops below $200\rho_{\rm crit}(z)$, where $\rho_{\rm crit}(z)$ is the critical density of the Universe at a given redshift $z$.

    \paragraph*{The Merger Trees}The progenitors of the haloes are tracked across the snapshots with \code{MergerTree}, a tool that comes with \code{AHF}. Each halo identified at redshift $z = 0$ is tracked backwards in time, identifying as the main progenitor at some previous redshift the halo that maximises the merit function $\mathcal{M} =N_{A\cap B}^2/(N_{A} N_{B})$, where $N_A$ and $N_B$ are the number of particles in haloes $H_A$ and $H_B$, respectively, and $N_{A\cap B}$ is the number of particles that are in both $H_A$ and $H_B$. The code further has the ability to skip snapshots, i.e. progenitors of haloes that are not found in the previous snapshot are still searched for in earlier snapshots, recovering an otherwise truncated branch of the merger tree \citep[see][]{Wang16}. However, to reduce errors during the tracking of the stellar kinematic evolution of our objects we consider an object `lost' if it cannot be found for 5 consecutive snapshots.\\

    In summary, we have 324 numerically modelled central galaxy clusters and all the haloes orbiting in and about them out to a distance of $15 h^{-1}$ Mpc available for our analysis. The details of this full data set are presented in \citet{Cui18}. Here we are though applying a few selection criteria to both the central clusters and the field haloes and subhaloes in the regions, to be described now.

    \subsection{Cluster Selection} \label{sec:data-subsec:sample}

    As we aim at tracing back all the objects in and about the central galaxy cluster, our analysis requires that we can always define a main progenitor for each of the 324 central haloes. As shown in \citet{Behroozi15}, major mergers during the formation of those central objects pose a challenge to this. Following the discussion in \citet{Haggar20}, we therefore identified and removed central cluster haloes whose main progenitor's position changes by more than half their radius $R_{200}(z)$ between two consecutive snapshots $[z,z+\Delta z]$. We further require the main branch to at least extend to redshift $z=2$. This reduces the number of regions entering the analysis from 324 down to 236. This is approximately 8 times more cluster regions than used in other state-of-the-art cluster simulation studies \citep[e.g.][]{ Bahe17,Barnes17}.

    \subsection{Halo Selection} \label{sec:data-subsec:haloselection}

    Each of our 236 selected central galaxy clusters is surrounded by a multitude of haloes out to the  $15 h^{-1}$ Mpc edge of the region that was modelled including all the relevant baryonic physics. But as we are interested in studying the angular momentum of the stellar component of these haloes, we are limiting our analysis to haloes with at least $M_{\rm {star}} > 5\times10^{10}{\ifmmode{h^{-1}{\rm {M_{\odot}}}}\else{$h^{-1}{\rm{M_{\odot}}}$}\fi}$. This corresponds to at least $1000$ star particles. Further, only those haloes that lie within $3R_{200}$ of the central galaxy cluster at redshift $z=0$ are traced backwards in time. For these haloes we then define their infall redshift $z_{\rm{inf}}$ as the redshift at which the halo crosses central galaxy cluster's $2R_{200}(z)$ for the first time: recent numerical studies show that around $\sim 1.5R_{200}-2R_{200}$, haloes experience a sharp cut-off in their gas content which could indicate the presence of an accretion shock \citep{Power18, Arthur19, Baxter21} and hence we decided to use $2R_{200}$ (instead of $R_{\rm 200}$) as our reference crossing radius. Haloes that cannot be assigned an infall redshift will be removed from the analysis.\\

    \subsection{Stellar Component} \label{sec:data-subsec:stellarcomponent}

    To select the stellar component of our objects -- that one might identify with their central galaxy -- we opted for a spherical region enclosing 15 per cent of the haloes' physical radius \citep[e.g.][]{Bailin05}. However, as tidal interactions with the central galaxy cluster impact on the size of subhaloes \citep[][]{Muldrew11,Onions12} we used the radius as found right before crossing $2R_{\rm 200}$ of the central galaxy cluster. This aperture is now interpreted as the `size of the galaxy' and kept fixed in physical coordinates across the snapshots. We performed a series of tests using different criteria (e.g. an aperture not depending on the halo's physical radius but rather fixing it to 30kpc) that showed that even though haloes suffer an overall mass loss while orbiting within their host halo, the central stellar mass in the aperture is mostly unaffected. This agrees with other simulations, which show that stellar stripping is rare, and happens only after the dark matter has suffered significant stripping \citep[$> 80$ per cent, e..g][]{Smith16,Bahe19}. It is important to point out again that even though the stellar mass of our galaxy does not change (something we also quantify below), we cannot rule out any back-reaction of the stars to the change in potential caused by tidal stripping of the halo. Finally, to take into account the spatial resolution of our simulations we restrict our objects to have a `galaxy size' of at least $2\epsilon = 13{\ifmmode{h^{-1}{\rm {kpc}}}\else{$h^{-1}{\rm{kpc}}$}\fi} \sim 20$ kpc. We remark that this additional condition removes only the smallest galaxies in the sample, which is less than 5 per cent of the objects selected so far. Moreover, practically all of the star particles reside within our aperture and hence we are not considerably reducing the number of star particles by cutting out the `galaxy' as defined here. The remaining 6509 objects, with halo masses from $9.2\times10^{10}{\ifmmode{h^{-1}{\rm {M_{\odot}}}}\else{$h^{-1}{\rm{M_{\odot}}}$}\fi}$ to $5.3\times10^{14}{\ifmmode{h^{-1}{\rm {M_{\odot}}}}\else{$h^{-1}{\rm{M_{\odot}}}$}\fi}$, constitute the analysis sample, built by combining all the objects that satisfy the aforementioned criteria from the 236 cluster regions considered here.

\section{Results} \label{sec:results}
    In what follows we are comparing the properties of our selected haloes orbiting in and about the central galaxy clusters at infall redshift $z_{\rm inf}$ (i.e. when crossing $2R_{\rm 200}$) and present day time $z=0$.

    \subsection{Mass Evolution}\label{sec:results-subsec:mass}

    \begin{figure}
    \includegraphics[width=\columnwidth]{plots/deltamcompare.pdf}
    \caption{Evolution of the different mass components of the haloes in the sample from their infall redshift $z_{\rm{inf}}$ to their $z=0$ values. Contour levels show the 75, 50, 25 and 5 per cent of the maximum counts of the distributions.  The horizontal lines at 0 and -1 correspond to no change and maximum mass loss, respectively. The position in the plane of objects beyond the 5 per cent contour are shown as square markers.}
    \label{fig:deltam}
    \end{figure}

    To understand the processes that shape the haloes after their first infall we show in \Fig{fig:deltam} the fractional change in mass since their infall redshift and as a function of infall mass. Contour levels show 75, 50, 25, and 5 per cent of the maximum counts of the distributions. The horizontal lines at 0 and -1 correspond to no change and maximum mass loss, respectively. Objects beyond the 5 per cent contour are shown as square markers. As expected, haloes crossing the denser regions of the central galaxy cluster lose a considerable amount of mass  \citep[e.g.][]{Klimentowski10}. Haloes that only entered the central halo once and did not leave, in general, conserve more of their initial mass than the ones that experienced multiple infalls into the denser regions of the galaxy cluster. In terms of the stellar mass in the haloes, the loss is significantly lower due to the fact that most of the stellar component resides in the central region of the halo mainly shielded from tidal forces. We nevertheless observe an increase in stellar mass loss at the 5 per cent level as we move to the massive end of the stellar mass distribution. Such haloes present substructure residing in the halo component (i.e. outside the central aperture) and thus they are prone to suffering mass loss from processes that strip the halo at its outskirts. The gas component, on the other hand, is mostly gone regardless of their mass, as gas is affected by an entirely different set of processes as they fall into the central halo, i.e. ram-pressure stripping \citep{Arthur19, Mostoghiu21}.  Studying the radial distribution of the gas inside infalling haloes (not shown here), we find that haloes that never reached the denser regions of the cluster halo are still considerably stripped of their gas component, in agreement with \citet{Power18, Arthur19, Mostoghiu21}.\\

    In summary, \Fig{fig:deltam} clearly shows that the stellar component of our haloes is the least affected by the environment of the central galaxy cluster. The question nevertheless remains if we will find changes in the kinematical properties. But the haloes of our `galaxies' (as defined by the star particles in the aperture, cf. \Sec{sec:data-subsec:stellarcomponent}) certainly undergo some changes as manifested by the mass loss. It is therefore not yet clear that even though the galactic stellar mass more or less remains constant that there will be no reaction of its internal dynamics to the varying influence of the central galaxy cluster.

    \subsection{Fraction of co-rotational energy}\label{sec:results-subsec:kappa}

    We study the influence of the central galaxy cluster on the internal dynamics of our galaxy sample after their first infall by following the evolution of their specific stellar angular momentum $j_{*}$.

    The galaxies in our sample can be classified by the fraction of stellar kinetic energy that is invested in co-rotation $\kappa_{\rm{corot}}$, as presented in  \citet{Sales10,Correa17}:
    %
    %
    \begin{equation}\label{eq:kappa_def}
        \kappa_{\rm{corot}} = \frac{E_{\rm{corot}}}{E_{\rm{kin}}}\textrm{   ,   }
        E_{\rm{corot}}=\sum_{j_{z,i}>0}{ \frac{1}{2}m_{i}\left(\frac{j_{z,i}}{r_{2d,i}}\right)^2}\textrm{   ,   }
    \end{equation}
    %
    %
    %
    where $E_{\rm{kin}}$ is the kinetic energy of the star particles in the aperture, $E_{\rm{corot}}$ is the rotational energy of the corotating star particles contributing to the rotation of the galaxy, $m_{i}$ the mass of the star particle;
    %
    %
    \begin{equation}
    \begin{split}
        j_{z,i}   &=\bm{j}_{i}\cdot\bm{\hat{J}_{\rm{tot}}}\textrm{ ,} \\
        \bm{j}_{i} &=(\bm{r}_i-\bm{r}_{\rm{halo}})\times(\bm{v}_i-\bm{v}_{\rm{halo}})\textrm{ , and } \\
        \bm{\hat{J}}_{\rm{tot}} &= \frac{1}{J_{\rm{tot}}} \frac{\sum{m_{i}\bm{j}_i}}{\sum{m_{i}}}\textrm{   ,   }
    \end{split}
    \end{equation}
    %
    %
    %
    the specific angular momentum along the direction of the total angular momentum, the specific angular momentum of a star particle in the rest frame of its halo, and the direction of the total specific angular momentum of the galaxy, respectively; and $r_{2d,i}=(\lvert\bm{r}_i-\bm{r}_{\rm{halo}}\rvert^2 -((\bm{r}_i-\bm{r}_{\rm{halo}}) \cdot \bm{\hat{J}}_{\rm{tot}})^2)^{1/2}$ the cylindrical radius of star particles. As we exclusively use the definition based on co-rotating star particles, we drop the subscript from the parameter $\kappa$ from now on.

    \subsubsection{Changes in $\kappa$ since infall}\label{sec:results-subsec:kappa-subsubsec:kappainf}

    \begin{figure}
    \includegraphics[width=\columnwidth]{plots/kappaCorotAp_compare.pdf}
    \caption{Fraction of ordered rotation $\kappa$ at $z=0$ and at their infall redshift $z_{\rm{inf}}$. The diagonal dashed line shows the $1:1$ relation.  Median values in bins of $\kappa (z=0)$ for the sample are represented by the solid cyan line. Bins are coloured by their corresponding normalised counts. The contour lines show where 75, 50, 25, and 5 per cent of the maximum counts lie.}
    \label{fig:kappaCorotAp_compare}
    \end{figure}

    We compared the fraction of ordered rotation at two different times in the evolution of the galaxies to find if the $z=0$ values are the result of the influence of the central galaxy cluster during their infall. In \Fig{fig:kappaCorotAp_compare} we show the $\kappa$ values of the sample at $z=0$ compared to their values at infall $z_{\rm{inf}}$, before experiencing the central cluster's influence. The black dashed diagonal line shows the $1:1$ relation. The median values in bins of $\kappa$ at $z=0$ are represented by the solid cyan line. Bins are coloured by their corresponding number of counts and contours show where 75, 50, 25, and 5 per cent of the maximum counts lie.

    We find that for most of the haloes in the sample the fraction of stellar ordered rotation is unaffected by their infall (Pearson coefficient of 0.56), with a median value at $z=0$ and at $z_{\rm{inf}}$ of $\kappa\sim 0.22$. However, beyond the 5 per cent contour we identify galaxies which experienced a considerable change in $\kappa$ since their infall, e.g. $\Delta\kappa \sim \pm0.2-0.4$.

    \begin{figure}
    \includegraphics[width=\columnwidth]{plots/deltamdeltakappa_compare.pdf}
    \caption{Evolution of the different mass components of the haloes in the sample and their fraction of ordered rotation $\kappa$ from their infall redshift $z_{\rm{inf}}$ to their $z=0$ values. Contour levels show the 75, 50, 25 and 5 per cent of the maximum counts of the distributions.  The horizontal lines at 0 and -1 correspond to no change and maximum mass loss, respectively. The position in the plane of objects beyond the 5 per cent contour are shown as square markers.}
    \label{fig:deltamdeltakappa}
    \end{figure}

    In \Sec{sec:results-subsec:mass} we found that the haloes in the sample are stripped of their material during their fall into the into the denser regions of the central galaxy cluster. Nevertheless, mass-loss does not appear to considerably affect the rotational properties of their galaxies, as seen from their $\kappa$ parameter. We quantify this change in \Fig{fig:deltamdeltakappa}, in which we show the same evolution of the mass component of the haloes in the sample presented in \Fig{fig:deltam} but this time as a function of the evolution in their fraction of ordered rotation.

    For the total halo and DM mass components of the haloes in the sample we identify objects that lost up to 90 per cent of their mass at infall, yet we do not find a correlation between the amount of mass-loss and their fraction of ordered rotation evolution trend. At the 5 per cent contour, regardless of the amount of stripped mass, their fraction of stellar ordered rotation changed within $\sim\pm60$ per cent since their infall. As we discussed in the previous section, the stellar component is less affected by mass-loss processes. Just so, we see that at the 5 per cent contour the objects experienced a change in stellar mass within $\sim\pm30$ per cent and that their change in $\kappa$ tends towards zero as their infall mass is disrupted. This is likely due to the fact that such stellar mass gain/loss tends to happen in the outer region of haloes, outside the apertures used to define the central galaxy and calculate $\kappa$. As for the gas component, even when haloes have been mostly stripped of their gas by the time they reach $z=0$, the stellar fraction of ordered rotation does not correlate with such mass-loss. Our galaxies can fall into galaxy cluster regions, losing most (if not all) of their gas during their infall, and still retain the stellar kinematic properties they had prior their infall. Similar results have been reported in \citet{Cortese19}: satellite galaxies go through significant changes in their specific star formation rate, but those are not necessarily accompanied by changes in their stellar spin parameter.

    Overall, we find no correlation between the amount of stripped mass and the change in their $\kappa$ parameter. Thus, we conclude that the fraction of ordered rotation in our simulations is hardly affected by the (potentially violent) stripping processes that disrupt the dark matter halo.

    \subsubsection{Temporal evolution of high and low $\kappa$ galaxies}\label{sec:results-subsec:kappa-subsubsec:kappaovertime}

    In \Fig{fig:kappaCorotAp_compare} we identified objects beyond the 5 per cent contour with a considerable change in their fraction of stellar co-rotational since their infall. To study which processes are responsible for such changes we selected from the $\kappa$ distribution at $z=0$ the galaxies within 5 per cent of the highest and lowest fractions of ordered rotation, i.e. the $\kappa > 95$-th percentile, the high $\kappa$ sub-sample; and $\kappa < 5$-th percentile, the low $\kappa$ sub-sample of the distribution. This corresponds to a high $\kappa$ threshold of 0.38 and a low $\kappa$ threshold of 0.15. Note that, our choice for the threshold values is purely motivated by kinematics: we simply aim to understand the origin of the substantial changes in the $\kappa$ distribution of a (relatively) small fraction of the total sample.  Adopting other threshold values \citep[which also take into account other non-kinematic properties such as star formation or colours, e.g. $\kappa = 0.4$ from][]{Correa17} for the stellar kinematic classification, only 4 per cent of our sample would be classified (in terms of their kinematics) as fast-rotating galaxies.

    Recent numerical results show that galaxies which continue to accrete gas and form stars are very efficient at spinning up \citep{Lagos17}. However, this is unlikely to be the cause here as gas accretion is expected to be hampered in clusters. \citet{Lagos18a} found that galaxies can be spin up or down by mergers depending on their orbital orientation and gas content \citep[see][for similar results]{Schulze18, Lagos18b}. To investigate the origin of the change in the fraction of co-rotation after infall found in \Fig{fig:kappaCorotAp_compare}, in \Fig{fig:kappa_peritime} we study the $\kappa$ evolution of the low and high $\kappa$ sub-samples found at $z=0$ as a function of the time since their pericentre passage, defined as the closest approach of an infalling halo to the central galaxy cluster \footnote{Our definition of `pericentre' does not necessarily imply they are within $R_{200}$ of the central galaxy cluster.}. Contours show 75, 50, 25, and 5 per cent of the maximum counts in the distribution. We show the median value at each time bin with a cyan solid line, and when the number of haloes in a bin is less than 50 per cent of the maximum count we use a dashed line instead. The low and high $\kappa$ threshold values are represented by the red and blue horizontal dash-dotted lines, respectively, and the pericentre time is marked with a vertical dashed line. In each panel we additionally show the $\kappa$ evolution of a galaxy from the respective sub-sample with green square markers, and we marked the allowed $\kappa$ values of each sub-sample with corresponding shaded regions.

    \begin{figure}
    \includegraphics[width=\columnwidth]{plots/kappa_periTime.pdf}
    \caption{Evolution of the fraction of co-rotational energy as a function of the time since the pericentre passage for low (left) and high (right) $\kappa$ sub-sample as classified by the $5-$th and $95-$th percentiles of the $z=0$ $\kappa$ distribution in \Fig{fig:kappaCorotAp_compare}. Contours show 75, 50, 25, and 5 per cent of the maximum counts of the 2D distribution. The cyan solid line shows the median value at each time bin. For bins in which less than 50 per cent of the sample entered the median calculation, we instead use a dashed line. The red and blue horizontal dash-dotted line shows the $\kappa$ threshold values for the low and high $\kappa$ sub-samples, respectively. The pericentre is marked with a vertical dashed line. In each panel we additionally show the evolution of a galaxy from the corresponding sub-sample using green square markers. The shaded red and blue regions show the allowed $\kappa$ values corresponding to each sub-sample.}
    \label{fig:kappa_peritime}
    \end{figure}

    As already anticipated by the $\kappa$ evolution in \Fig{fig:kappaCorotAp_compare}, most of the galaxies in each sub-sample retain their fraction of co-rotational energy $\kappa$ after their infall. However, over the course of 5 Gyrs after the pericentre passage we observe a decrease of $\sim0.2$ in $\kappa$ in both sub-samples ($\sim0.04/$Gyr) at the 5 per cent contour. Moreover, the $\kappa$ evolution of the low $\kappa$ sample galaxy shows that it had high stellar co-rotational energy at some point, but lost it after $\sim6$ Gyr. Such slow decrease can be attributed to the two-body heating of the stellar component: as galaxies infall into the central galaxy cluster, the background (more massive) dark matter particles residing inside the cluster that fly by these galaxies tend to increase the stellar mean interparticle distance due to the softening scale used for the stellar particles, heating up the stellar component and effectively puffing up the stellar distribution while spinning them down. As these galaxies fall further into the cluster, the density of dark matter particles increases and consequently the amount of fly-bys is boosted. Along with this slow decrease, we also identify processes which instead can spin-up galaxies on a faster scale. But these changes -- as observed for the high $\kappa$ sample galaxy -- are happening close to pericentre passage and are investigated in more detail now.

    \subsubsection{Spinning-up galaxies with low $\kappa$}\label{sec:results-subsec:kappa-subsubsec:spinup}

    As two-body heating slows down the stellar kinematics of every infalling galaxy in equal measure, i.e. $\kappa$ decreases at the same rate for both sub-samples, we focus on processes which are able to spin-up galaxies during their infall. To isolate the galaxies which cross the threshold $\kappa$ value to become high $\kappa$ galaxies according to our classification, from the ones that conserve their classification up to $z=0$, i.e. the ones that remained within the high $\kappa$ sub-sample after their infall into the cluster environment, we select galaxies from the high $\kappa$ sub-sample which experienced a rapid $\kappa$ increase (i.e. $\Delta\kappa > 0.2$) within 1 Gyr since their pericentre passage. These galaxies form only 12 per cent of the high $\kappa$ galaxies at $z=0$ (39 objects, 0.6 per cent of the total number of objects in the sample). Within this new sub-sample, we identify galaxies for which a sudden change in their star particle count was accompanied by a sudden spin-up of their specific angular momentum (and consequently an increase in their $\kappa$ parameter) near their pericentre, and galaxies for which the number of stars close to the pericentre does not seem to play a crucial role in their $\kappa$ evolution.

    Following the star particle distribution of these galaxies, we find that the fast increase in $\kappa$ and in the number of stars is a consequence of two processes: mergers with other haloes, and fly-bys of stellar remnants within the central aperture used to define the galaxies residing in each halo. In agreement with \citet{Lagos18a}, we find that the co-rotating (counter-rotating) infalling stars from mergers are able to spin-up (spin-down) our galaxies. On the other hand, acting on shorter timescales and without a significant stellar gain, we find that stellar remnants are able to temporarily disrupt the spin of our objects. These haloes are the residual cores of infalling haloes, where the original dark matter component of the halo has become lost and subsumed by the main halo. Considering that the stellar component of the remnant haloes is smaller than the central galaxy extension of the haloes in the sample ($\sim 1/3$ of the aperture size), and that the amount of remnant stellar cores increases close to the centre of the galaxy cluster, these objects contribute transiently to the co-rotational energy in the aperture and, as such, $\kappa$ increases near the pericentre of our objects.\\

    \subsubsection{Summary}\label{sec:results-subsec:kappa-subsubsec:summary}

    We conclude that for most of the galaxies in our sample entering massive galaxy clusters, their kinematic properties (as captured by the fraction of ordered rotation $\kappa$) do not change in a significant manner -- even though we have seen in \Sec{sec:results-subsec:mass} that their haloes undergo substantial changes. In general, low $\kappa$ galaxies at $z=0$ had low fractions of ordered rotation before entering the cluster halo, and high $\kappa$ galaxies at $z=0$ had such high fractions prior to entering the cluster environment. The two-body heating of the stellar component of infalling haloes induced by the more massive dark matter particles in the cluster environment affects equally both samples, slowly reducing their stellar specific angular momentum as they orbit the cluster region. As such, we find galaxies that at their infall time had high fractions of stellar co-rotational energy ($\kappa \gtrapprox 0.4$) that ended up as low $\kappa$ galaxies by the time they reached $z=0$, i.e. with a $\kappa<5-$th percentile of the $\kappa$ sample distribution at $z=0$. On the other hand, from the high $\kappa$ sample at $z=0$, i.e. the galaxies with $\kappa>95-$th percentile, we identify $\sim12$ per cent of them which previously had considerably lower fractions. We find that these two processes acting on different timescales, i.e. mergers with other infalling haloes and stellar remnants transiting the aperture used for defining the galaxies of the haloes in the sample, are able to spin-up the galaxies.

    \subsection{Stellar angular momentum and stellar mass relation}\label{sec:results-subsec:jM}

    \begin{figure}
    \includegraphics[width=0.9\columnwidth]{plots/jstarmstar_zcompare_all.pdf}
    \caption{Specific angular momentum and stellar mass relation of the sample of galaxies at $z=0$ and at their infall redshift $z_{\rm{inf}}$. The diagonal dashed lines show the $M_{\rm{star}}^{2/3}$ relation for zero-points $-3$ (top line) to $-6$ (bottom line). The contours show 75, 50, 25 and 5 per cent of the sample. The position in the plane of objects beyond the 5 per cent contour are shown as filled markers.}
    \label{fig:jstarmstar_zcompare}
    \end{figure}

    In the previous section we identified processes which are able to change specific angular momentum of infalling galaxies. Nevertheless, most of the objects remain unaffected by such interactions. In this section we study the influence of the central galaxy cluster on the internal dynamics of our galaxies by following their evolution in the specific angular momentum $j_{\rm{star}}$ -stellar mass M$_{\rm{star}}$ plane, an extensively studied fundamental kinematic property of galaxies which can be described a power-law of the form $j_{\rm{star}}\propto$ M$_{\rm{star}}^{\alpha}$, where $\alpha \approx 0.7$ \citep[e.g.][]{Fall83, Romanowsky12, Obreschkow14, Teklu15, Fall18}.

    As some of our galaxies experienced great changes in their fraction of ordered rotation during their evolution since their infall, we aim to find if these changes are also visible in the specific stellar angular momentum-stellar mass plane. We show their position in the $j-M$ plane after their infall redshift in \Fig{fig:jstarmstar_zcompare}. The contours represent 75, 50, 25, and 5 per cent of the sample. The position in the plane of galaxies beyond the 5 per cent contour are shown as filled markers. The dashed lines indicate the $j_{\rm{star}} = \lambda M_{\rm{star}}^{2/3}$ relation for different zero-points ($\log$ $\lambda$ = $-3$ to $-6$), as presented in \citet{Teklu15}. Note that, overall, the sample in this analysis does not include galaxies with stellar specific angular momenta typical of disk galaxies in the local Universe.

    Once again, we find that the distribution is mostly unchanged after entering the cluster halo. However, galaxies positioned furthest from the average position in the $j-M$ plane at their infall redshift move towards lower specific angular momentum and stellar mass values after infalling, as seen in the 5 per cent contours and points below that level.

    These results show that the internal dynamics of our infalling galaxies into massive galaxy clusters, as quantified by either the fraction of ordered rotation $\kappa$ or the evolution of their location in the specific angular momentum-stellar mass plane, are mostly unaffected by the processes in such environments even after their halos have been considerably stripped of their mass, and consequently, disrupted the halo potential in which these galaxies reside.

\section{Conclusions} \label{sec:conclusions}
    We analysed the central stellar component of haloes from \textsc{The Three Hundred} project, a suite of 324 galaxy cluster regions simulated with full-physics hydrodynamics, to study of the influence of galaxy clusters on the internal dynamics of infalling galaxies. We selected objects with M$_{\rm{star}}>5\times10^{10}$ ${\ifmmode{h^{-1}{\rm {M_{\odot}}}}\else{$h^{-1}{\rm{M_{\odot}}}$}\fi}$ within $3R_{200}$ from the galaxy cluster halo at $z=0$. After applying different selection criteria that ensure the correct tracking of the stellar component of each halo we obtain a sample of 6509 galaxies from 236 cluster regions, i.e. approximately 8 times more cluster regions than used in other state-of-the-art cluster simulation studies \citep[e.g.][]{Bahe17,Barnes17}. Here we summarise our main results:

    \begin{itemize}

        \item Using the fraction of stellar co-rotational energy $\kappa$ to track the internal dynamics of the infalling galaxies in our sample we find that it remains mostly unchanged for most of the galaxies after they enter the galaxy cluster environment, despite the (potentially violent) stripping processes that disrupt their haloes.

        \item We studied the time evolution of the $\kappa$ parameter. Even though most galaxies do not experience great $\kappa$ changes during their evolution, we identify an overall slow decrease in the fraction of ordered rotation close to their pericentre passage ($\sim0.04/$Gyr) which we attribute to numerical effects. Along with this effect we found that $\sim12$ per cent of the galaxies with high $\kappa$ at $z=0$ (i.e. $\kappa>95-$th percentile of the $z=0$ distribution), an already small population in our sample, had considerably lower fractions of co-rotational energy in the past and experienced a fast boost in their $\kappa$ parameter within the past few Gyrs ($\sim0.2/$Gyr) due to two different processes: stellar mergers of infalling haloes \citep[in agreement with ][]{Lagos18b,Lagos18a}, and transient encounters with the stellar remnants of haloes that have lost their dark matter component during their own passage through the cluster, entering within the aperture we used for defining the properties of our galaxies.

        \item We found that similar to the $\kappa$ evolution of our galaxies, the specific angular momentum-stellar mass relation for the galaxies in the sample showed no substantial change in their location on the $j-M$ plane from their infall redshift until $z=0$.
    \end{itemize}

    Previous numerical studies \citep[e.g.][]{Smith16,Lokas20} have shown that infalling and orbiting haloes within galaxy clusters are disrupted by a series of processes which dominate in such environments, such as ram-pressure stripping \citep[e.g.][]{Bahe15,Arthur19, Mostoghiu21}, mergers \citep[e.g.][]{Behroozi14,Lagos18a} or tidal torques \citep[e.g.][]{Park07}. Nevertheless, these numerical studies primarily analyse mass-loss processes in infalling haloes. Considering that cluster environments disrupt the different mass components of infalling haloes, such environments could also potentially disrupt their stellar kinematics.

    In this work we investigate such question by studying the internal central stellar dynamics of infalling objects towards the numerically simulated massive galaxy clusters in our sample. We found that, roughly, their $z=0$ classification holds even from before the time they entered the cluster environment. Moreover, we find that in terms of their position in the fundamental specific stellar angular momentum-stellar mass plane \citep[e.g.][]{Obreschkow14,Fall18} they remained mostly unaffected by such processes. For a small number of galaxies, we identified processes which are able to considerably boost the spin of infalling galaxies during their infall via mergers with other objects and via transient encounters with stellar remnants within the cluster halo. These two processes are able to offset the numerical slow-down induced by numerical effects, such that at $z=0$ they end up within the top 5 per cent galaxies with the highest amount of ordered rotation within the simulation.

    Due the mass resolution of our simulations, which favours the modelling of a large number of galaxy clusters and their environments, our sample contains a limited kinematic distribution of galaxies when compared with observations \citep[e.g.][]{Brough17}. Thus, generalising these environmental constrains to galaxies with higher stellar angular momentum is not possible, as mass-loss processes during the infall of such galaxies in cluster environments might prove to play a more crucial role in their stellar angular momentum evolution. Future simulations, where the stellar mass resolution is increased and a smaller softening scale is used to better resolve the stellar component, for similar cosmological volumes, would help us extend our analysis to a wider range of galaxies with higher fractions of ordered rotation and obtain environmental constraints with far greater statistical significance.

\begin{acknowledgements}
    This work has been made possible by the `The Three Hundred' collaboration.\footnote{\url{https://www.the300-project.org}} The project has received financial support from the European Union's Horizon 2020 Research and Innovation programme under the Marie Sklodowskaw-Curie grant agreement number 734374, i.e. the LACEGAL project.

    The authors would like to thank the anonymous referee for their constructive feedback, which improved the quality of this paper. Additionally, they would like to thank The Red Espa\~{n}ola de Supercomputaci\'{o}n for granting us computing time at the MareNostrum Supercomputer of the BSC-CNS where most of the cluster simulations have been performed. Part of the computations with \code{GADGET-X} have also been performed at the `Leibniz-Rechenzentrum' with CPU time assigned to the Project `pr83li'.

    RM, AK, and GY would like to thank MINECO/FEDER (Spain) for financial support under research grants AYA2015-63819-P and MICIU/FEDER for financial support under research grant PGC2018-094975-C21. RM further acknowledges support from the International Centre for Radio Astronomy (ICRAR) for funding this project. AK further acknowledges support from the Spanish Red Consolider MultiDark FPA2017- 90566-REDC and thanks Joy Division for unknown pleasures. CL has received funding from the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. WC acknowledges the supported by the European Research Council under grant number 670193. SB acknowledges financial support from PRIN-MIUR 2015W7KAWC, the INFN INDARK grant, and the EU H2020 Research and Innovation Programme under the ExaNeSt project (Grant Agreement No. 671553). KD acknowledges support through ORIGINS, founded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311.

    The authors contributed to this paper in the following ways: RM, AK, FRP and CP formed the core team. RM analysed the data, produced the plots and wrote the paper along with AK, FRP, CL, WC, SB, and KD.  SB, KD, GM, and GY supplied the simulation data. All authors had the opportunity to provide comments on this work.

    This work was created by making use of the following software: \code{Python}, \code{Matplotlib} \citep{Matplotlib-Hunter07}, \code{Numpy} \citep{Numpy-vanDerWalt11}, \code{scipy} \citep{Scipy-Virtanen19}, \code{astropy} \citep{AstropyI,AstropyII}, and \code{pynbody} \citep{pynbody}.
\end{acknowledgements}

\bibliographystyle{aa}