The Three Hundred: Cluster Dynamical States and Relaxation Time Scale

Using the galaxy clusters from The Three Hundred Project, we define a new parameter: $\lambda_{DS}$ to describe the dynamical state of clusters, which assumes a double-Gaussian distribution in logarithm scale for our mass-complete cluster sample at $z=0$ from the dark-matter-only (DMO) run. Therefore, the threshold for distinguishing relaxed and unrelaxed clusters is naturally determined by the crossing point of the double-Gaussian fitting which has a value of $\lambda_{DS} = 3.424$. By applying $\lambda_{DS}$ with the same parameters from the DMO run to the hydro-dynamically simulated clusters (Gadget-X run and GIZMO-SIMBA run), we investigate the effect of baryons on the cluster dynamical state. We find a weak baryon-model dependence for the $\lambda_{DS}$ parameter. Finally, we study the evolution of $\lambda_{DS}$ along with clusters mass accretion history. We notice an upper limit of halo mass change $\frac{\Delta M_{200}}{M_{200}} \sim 0.12$ that don't alter the cluster dynamical state, i.e. from relaxed to unrelaxed. We define relaxation period (from the most relaxed state to disturb and relaxed again) which reflects how long the dynamical state of a cluster restores its relaxation state, and propose a correlation between this relaxation period and the strength of halo mass change $\frac{\Delta M_{200}}{M_{200}}$. With the proposed fitting to such correlation, we verify the relaxation period can be estimated from $\frac{\Delta M_{200}}{M_{200}}$ (including multi mass change peaks) with considerably small error.


INTRODUCTION
Galaxy clusters are the largest gravitationally bounded structures in the Universe.Masses of clusters, ranging from ∼ 10 14 M to over 10 15 M , are dominated by dark matter (∼80% -85%) (e.g.White et al. 1993;Fukugita et al. 1998).The gravitational processes exhibiting by the dynamical properties of dark matter halo can be viewed from the baryonic compositions in galaxy clusters, such as the galaxies and intra-cluster medium (ICM), which emit photons at different wavelengths and thus can be observed through telescopes.As the "brightest" objects in the sky, especially in Xray band, the galaxy clusters, such as the Coma cluster, used to be the main targets of astronomical observation.Therefore, numerous studies have focus on the cluster properties.Cluster dynamical state -one unique feature to describe the cluster virialized state, connects to many other important cluster quantities such as its mass (e.g.Nelson et al. 2012;Biffi et al. 2016;Gianfagna et al. 2021), formation time (e.g.Mostoghiu et al. 2019), and concentration (e.g.Neto et al. 2007), etc. Galaxy clusters mergers, the most powerful events which E-mail: bz287@cam.ac.uk (BZ) † E-mail: weiguang.cui@ed.ac.uk (WC) violate the cluster equilibrium, provide unique conditions to study a range of physics (e.g.Poole et al. 2006;Zenteno et al. 2020).From the perspective of cosmology, disturbed clusters provide test bed for ΛCDM model (e.g.Thompson et al. 2015;Kim et al. 2017;Sereno et al. 2018), and their enhanced strong lensing efficiency provide powerful tools to investigate the universe at high redshift (see Baldi et al. 2013;Acebron et al. 2019, for example).
From the perspective of galaxy evolution, cluster dynamical state is also valuable.For example, Morell et al. (2020) found that galaxies evolve in the same way in Gaussian and Non-Gaussian systems, but their formation histories leads to different mixtures of galactic types and infall patterns.Furthermore, disturbed clusters can be used to examine the effect of the ram pressure from the cluster ICM thanks to these jellyfish galaxies inside (McPartland et al. 2015).Lastly, we found that the halo formation time affect the central BCG properties (Cui et al. in prep, see Cui et al. (2021) for a similar result at lower halo mass), thus this can be also linked to the cluster dynamical state.
At the late stage of the hierarchical structure formation process, the recent formation history of massive structures such as galaxy clusters, correlates strongly with the degree of its dynamical equilibrium (e.g.Wong & Taylor 2012).This is easy to understand as a recent major merger will result the cluster in dynamical disequilibrium for a certain time, while earlier major merger means a following decreased mass accretion.As a result of those highly complicated dynamical processes of merger events, galaxy clusters can have a wide range dynamical states, which can be reduced to two categories: relaxed(or virialized) and non-relaxed(or non-virialized).Relaxed clusters are expected to have nearly spherical shape and Gaussian velocity distribution (e.g.Faltenbacher & Diemand 2006), while non-relaxed clusters can show elongated shapes (Gouin et al. 2021), non-Gaussian velocity distribution (Hou et al. 2009), the presence of massive substructures (Lopes et al. 2018) and irregular morphology properties (e.g.De Luca et al. 2021).
With the purpose of better utilization of galaxy clusters as probes of cosmology and testbeds for galaxy formation and evolution, it is of great importance to interpret the dynamical processes in galaxy clusters.Different methods for classifying cluster dynamical states and assessing the relaxation degree of clusters are found in the literature (see Cui et al. 2017, and references therein for theoretical approaches), (see De Luca et al. 2021, and references therein for observationallike approaches).In observations, cluster dynamical states is classified by and related to the ICM hydrostatic equilibrium: clusters that are undergoing, or have undergone a merger process, leave the ICM in turbulent.Wen & Han (2013) developed method to diagnose the substructure and dynamical state of galaxy clusters by using photometric data of Sloan Digital Sky Survey.Capalbo et al. (2021) and De Luca et al. (2021) investigated the correlation between cluster dynamical states and cluster morphology, which is directly measured through images of the surface brightness in the X-ray band and of the y parameter from the thermal Sunyaev-Zel'dovich effect in the millimetre range.In theory, there are a vague of ways with which halo dynamical state can be evaluated.Using DMO simulations, Bett et al. (2007) used integrated virial ratio 2T /|W |+1 to classify dynamical states, and suggested 2T /|W |+1 < 1.5 to select halos in quasi-equilibrium states.Neto et al. (2007) expanded the criteria by including substructure mass fraction and centre-of-mass offset, which contain the information of the constituents in cluster and the shape of cluster respectively.Shaw et al. (2006) additionally took the surface pressure energy Es into account in virial ratio calculation (see also Cui et al. 2017, for detailed calculation for hydrodynamic simulations).Davis et al. (2011) found the effect of the potential energy from particles outside of halos is negligible.Nevertheless, all of these methods require manually set thresholds to split the cluster dynamical state into relaxed and un-relaxed (disturbed) in both theory and observation, which doesn't generally showing up in their distributions (in either single or combined parameters De Luca et al. 2021;Haggar et al. 2020).Therefore, our first task in this paper is to remove these thresholds and present a new parameter-free cluster dynamical state classification.
One key question which is not well thoroughly studied in the literature is the relaxation time scale of the cluster dynamical state, i.e. how long does it need from getting disturbed to relaxed in hydrostatic equilibrium.This will help us to understand the cluster thermalisation (Sereno et al. 2021).Furthermore, it is also interesting to see the effect of baryons on the cluster dynamical state and their evolution.
The layout of this paper is as following: we introduce the The Three Hundred project in section §2.The new parameter-free cluster dynamical classification method and separation of relaxed and un-relaxed clusters are presented in §3.Our main results on the cluster dynamical state is shown in §4.We finally conclude and discuss our study on cluster dynamical state in §5.

THE THREE HUNDRED PROJECT
The The Three Hundred1 consists of 324 re-simulated clusters and 4 field regions extracted from the MultiDark Planck simulation, MDPL2 (Klypin et al. 2016).The MDPL2 simulation has cosmological parameters of ΩM = 0.307, ΩB = 0.048, ΩΛ = 0.693, h = 0.678, σ8 = 0.823.All the clusters and fields have been simulated using the full-physics hydrodynamic codes Gadget-X (GX in short, Rasia et al. 2015;Steinborn et al. 2015;Beck et al. 2016), Gadget-MUSIC (Sembolini et al. 2013), which are updated versions of Gadget2 (Springel 2005) and GIZMO-SIMBA (GIZMO in short, Davé et al. 2019, andsee Cui et al. 2021 in prep.for the details of the The Three Hundredcluster run) which is developed from MUFASA (Davé et al. 2016) using the GIZMO code (Hopkins 2015).In the re-simulation region, the mass of a dark matter particle is 12.7 × 10 8 h −1 M and the mass of a gas particle is 2.36×10 8 h −1 M .Each cluster resimulation consists of a spherical region of radius 15 h −1 Mpc at z = 0 centred on one of the 324 largest objects within the host MDPL2 simulation box, which is 1 h −1 Gpc on a side.The halo masses of central galaxy clusters range from ∼ 6.4 × 10 14 h −1 M to 2.63 × 10 15 h −1 M .
A more detailed introduction of The Three Hundred can be found in Cui et al. (2018).Besides these studies on the cluster dynamical state which has been mentioned in the introduction, these simulated galaxy clusters have been used for different proposes: the filaments around the clusters (Rost et al. 2021;Kuchner et al. 2020Kuchner et al. , 2021;;Kotecha et al. 2021); the backsplash galaxies (Haggar et al. 2020;Knebe et al. 2020) and shock radius (Baxter et al. 2021;Anbajagane et al. 2021b).The advanced baryon models in hydrodynamic simulations allow us to perform a detailed investigation on the cluster properties, such as profiles (Mostoghiu et al. 2019;Li et al. 2020), substructure and its baryonic content (Arthur et al. 2019;Haggar et al. 2021;Mostoghiu et al. 2021b,a), the cluster (non-)thermalization (Sayers et al. 2021;Sereno et al. 2021), the fundamental plane (Díaz-García et al. 2021), and the cluster mass bias (Ansarifard et al. 2020;Li et al. 2021;Anbajagane et al. 2021a).Additional runs allow us to investigate more things: such as the effect of environment by comparing to void/field regions (Wang et al. 2018); constraining the dark matter cross-section with the self-interacting dark matter run (Vega-Ferrero et al. 2021); examining the chameleon gravity (Tamosiunas et al. 2021).
In this paper, we only use the halos identified by the Amiga's Halo Finder (AHF Knollmann & Knebe 2009) with a spherical overdensity of 200 ×ρcrit.The progenitors of these halos are tracked and identified using the Mergertree that is part of the AHF package.We only focus on the main progenitors of the cluster, which is defined as the highest matched halo in previous snapshot, for tracking their mass accretion history.

THE PARAMETER-FREE CLASSIFICATION OF CLUSTER DYNAMICAL STATE
3.1 Dynamical parameters and previous work on classifying cluster dynamical states In the literature, for example, Cui et al. (2017), different parameters are used to describe the dynamical states of clusters.
The commonly used parameters are: • The Virial Ratio, η.
The exact expression for the virial theorem is 1 2 where I is the moment of inertia, T and W are kinetic energy and potential energy respectively, and Es is the energy from surface pressure P .
If the cluster system is in dynamical equilibrium, the Equation 1 will reduce to which can be rewritten as Therefore, the virial ratio is defined as and a relaxed cluster is expected to have η ≈ 1.
• Subhalo Mass Fraction, fs.fs represents the fraction of total mass of cluster contained in subhalos, which are identified by AHF.However, this fraction doesn't include the most massive substructure since it only includes the bounded components in main halos.For the most relaxed clusters, the subhalo mass fraction should be very small.
• Center of Mass Offset, ∆r.
The offset of the center of mass of cluster is defined as where Rvir is virial radius, within which virial theorem applies for a bounded system.Rc is the position of the peak of the density field of cluster, and Rcm is the position of the center of mass, which can be calculated by where mi and ri are the mass and position of the i th particle, M is the virial mass of the halo and n is the total number of particles within Rvir.
Empirically, a gravitationally bounded system in equilibrium has symmetric mass distribution, which requires small distance between the center of mass and the peak of density.Therefore, small ∆r is expected for a relaxed cluster.
We emphasise here that all the three parameters are only phenomenal descriptions2 of the cluster dynamical state due to the lack of a physically defined quantity for it.Furthermore, it is not clear which parameter contribute more to or describe better the cluster dynamical state.Therefore, different criteria are applied to classify a cluster as relaxed.For example, Cui et al. (2017) concluded that a relaxed cluster should satisfy three criteria: ∆r < 0.04, fs < 0.1 and 0.85 < η < 1.15 .With these criteria, Haggar et al. (2020) combined these three parameters which are normalised to their thresholds but with equal weight, to a continuous, nonbinary measure of cluster dynamical states, which is defined as the "relaxation" parameter of cluster, χDS: χDS = 3 For a cluster to be dynamically relaxed, it requires ∆r and fs to be minimised, and η ≈ 1.Therefore, the most relaxed clusters are expected to have large χDS (χDS > 1).However, there is no clear separation between the relaxed and unrelaxed clusters.Actually, the distribution of the χDS exhibits a single peak Gaussian curve.Therefore, we need to manually set a threshold to break the cluster dynamical state into relaxed and unrelaxed for other studies.

The threshold-free λDS function
As discussed before, the common issue in all previous works of classifying cluster dynamical states with either single or multiple dynamical parameters is that the thresholds for these parameters are chosen arbitrarily.In order to overcome such issue, we assume that a mass-complete cluster sample at z = 0 can be roughly separated into dynamically relaxed and unrelaxed from the DMO simulations.The DMO instead of hydro simulations are chosen because the DMO simulations are very robust, different codes give very small difference (e.g.Sembolini et al. 2016a,b) -unlike hydro-simulations of which the internal structures can be significantly altered due to different baryon models (e.g.Cui et al. 2016;Elahi et al. 2016).
Here, we introduce a new relaxation parameter, λDS here, which is a modification version of Equation 7: Instead of using prefactors of ∆r, fs and |1−η| terms representing their thresholds of cluster dynamical states, our new λDS completely remove these threshold with its prefactors a and b fitting-determined from our assumption -λDS has a Double-Gaussian distribution over the DMO clusters.Since we only care about the distribution pattern of λDS, instead of their absolute values, the only important thing in Equation 8 is the relative contributions from ∆r, fs and |1 − η|.Therefore, we can arbitrarily set one prefactor from one term in the denominator to be 1, i.e here we choose |1-η| term to have unit prefactor.The prefactors a and b can then be determined by fitting the distributions of λDS with a double-Gaussian function, and finding out the pair of a and b which can give the "best" Double-Gaussian distribution.After the two families of clusters are fitted, it is naturally to use the crossing point of the two Gaussian curves as the threshold for separating the cluster into relaxed and unrelaxed.

Determine coefficients a and b
To determine the a and b parameters, a two dimensional array is made.Each element in this array is a pair of trial a and b, with the separation of 0.01 between neighbouring elements.For each pair of a and b, the distribution of λDS, estimated by Equation 8, is used to fit a double-Gaussian function with the free parameters c1, c2, µ1, µ2 and σ1, σ2 to be determined by Several constraining criteria are made to select a and b values.
Firstly, the list of λDS must not pass the Shapiro-Wilk test (SHAPIRO & WILK 1965), which aims at testing normal distribution with a single peak.shapiro() function can be directly imported from scipy.stats, and it will return an indicator called P value when acting on a list-like object.The distribution is normal if its P value is greater than 0.05.Because we are looking for a double-Gaussian distribution, which shouldn't be normal.Therefore, we require P value < 0.05 for the list of λDS.
These four constraints are summarised below, as: Lastly, with all qualified pairs of a and b selected, the best candidate is the one which has the least square fitting error E. The square fitting error can be calculated from the real distribution of λDS with the predicted values from the fitting result, which is defined as In this case, yn is the height of the n th bin, which represents the number of clusters with λDS values within the range of the bin centring at xn. N is the total number of bins and f (xn) is the fitting value of the Double-Gaussian function at xn.
The best fitting parameters for the DMO run are a = 7.30 and b = 0.30.The λDS distribution with the fitting results is plotted in Figure 1.The best fit parameters for the Double-Gaussian function are listed in Table 1.The value of a is much larger than b, which indicates that the λDS function place far more emphasis on centre-of-mass offset than both the sub-halo mass fraction and the virial ratio.This implies that centre-of-mass offset ∆r plays an dominated role in λDS.We also notice that there is a linear correlation between ∆r and fs, so the small contribution from sub-halo mass fraction provides additional information going beyond the linear correlation.
We specially note here that although the method is reliable, the fitting parameters a and b, thus the threshold for separating relaxed and unrelaxed clusters, can be sample dependent, i.e. reducing or increase the minimum cluster mass in the mass-complete sample, a and b can be changed slightly.However, we are limited to our sample in this study, and as long as we are consistent, our results won't quantitatively change.This is because these key quantities: fs, ∆ and η are all physical, which should not change along mass and redshift.Therefore, the same λDS classified as relaxed at z = 0 or for cluster with higher mass, should be equally relaxed at high z or a lower mass.Further investigation regarding the changes of a and b at with different samples requires much large simulation, we leave it for a later study.

The threshold for separating relaxed and unrelaxed clusters
The Double-Gaussian distribution of λDS naturally avoid the arbitrary choice of the threshold for separating dynamically relaxed and unrelaxed clusters.The threshold is defined as the x coordinate of the crossing point of two Single-Gaussian functions, see Figure 1.The threshold value is λDS = 3.424 in normal scale.
In different hydro-dynamical simulations, different baryonic models are used, which can result in different best-fitted Double-Gaussian functions.For the simplicity in this investigation, and in order to highlight the changes due to different baryon models, we apply the same fitting results from DMO fitting as the baseline, i.e. with a = 7.3, b = 0.3, to calculate λDS in GX and GIZMO runs.The distributions of λDS in GX and GIZMO runs are plotted in Figure 2, along with the fitting results.Note that the distribution of λDS from GIZMO can not be fitted to double-Gaussian and we never expect that as the baryon models will change cluster dynamical state.Although the distribution of λDS from GX can be fitted to double-Gaussian, there is a shift of the thresh- old value compared to the DMO result.Nevertheless, applying the same parameters with threshold, we can examine the effects of baryons.For example, with the same threshold, λDS = 3.424, applied to GX and GIZMO run, we find that 151/171/170 clusters are classified as relaxed clusters in DMO/GX/GIZMO run.It looks that hydro-simulations with baryon model tend to increase the number of relaxed clusters.More details will be presented in §subsection 4.1.

The relationship between λDS and χDS
This new relaxation parameter, λDS, is compared with the corresponding old relaxation parameter, χDS, in Equation 7     Table 2.The median numbers of the differences of (from left to right) η, ∆r, fs and λ DS between GX run (the first row)/GZIMO run (the second row) and DMO run (numbers in brackets are standard errors), The third row displays the median numbers of each parameter in DMO run.

The baryon effect on the cluster dynamical state
We further investigate the baryon effect on dynamical parameters η, fs, ∆r and λDS.We matched the corresponding clusters from different runs at z = 0.For each cluster, the differences of these parameters between hydro-dynamical simulation (GX or GIZMO) and DMO simulation are shown in Figure 4. Distributions of these differences are plotted in histograms.For each parameter, the median numbers of these differences are used to quantify the baryon effect.Those median numbers and standard errors are marked as vertical lines in Figure 4 and listed in Table 2. Main results of baryon effects are discussed below: • η from the GX run and the GIZMO run are reduced by about 2% compared to the median value from DMO run.The differences distribution between the two hydro-dynamical simulations is very small, which gives the insight that the effect on η depends weakly on baryon models.Cui et al. (2017) showed a similar result on the weakly-model-dependent effect of the decrease in η, but with a larger difference, about 10%, for CSF run and AGN run.Here CSF run referred to a hydro-dynamical simulation ignoring the AGN feedback, and AGN run includes AGN feedback.They also concluded that the ratio between η from hydro-dynamical run and η from DMO run shows no dependence on cluster mass.
• Standard deviations of the differences between ∆r from GX/GIZMO run and from DMO run are comparable to the scale of the median number of ∆r from DMO run, which shows the scattering distribution of ∆r in hydro-dynamical simulation, in agreement with the result in Cui et al. (2017).This is the consequence of that the position of substructure can be largely affected by baryons.However, the average amount of change for all clusters is small, which behaves as a small decrease about 5% compared to DMO run.This could be mainly caused by the central galaxy formation which deeps the potential and increases the halo concentration.Thus, more weights are contributed from the central region.
• Comparing to the DMO value, fs increases by 17% in GIZMO run.This is in agreement with the result on Cui et al. (2017): their CSF run increases the fs by 40% for higher cluster mass, and by 20% for clusters with lower masses.
However, the median change of fs in GX run is negligible, smaller than 5%.The difference between fs from GX run and GIZMO run should come from the feedback models that control the galaxy formation in these less massive substructures.Through the comparison of satellite stellar mass function in Cui et al. 2021 (in prep.), it is clear that the satellite stellar mass function from the GIZMO run agrees better with the observation results at lower galaxy mass than GX which is about 5 times lower.
• λDS in GX (GIZMO) run is 9 (6) per cent higher than in DMO run, which indicates a weak the baryon-model dependence of λDS.This is not surprising as the baryon models have a weak influence on the individual parameters.

Dynamical state and cluster mass accretion history
It is clear that the cluster dynamical state changes are caused by the accretion of mass, especially in the case of major merger events.This is also revealed from the three key items in Equation 8.However, it is unclear how significant the cluster dynamical state can be altered and how long the cluster will return to relaxed state after a merger event.In this section, we will try to quantify the relationship between cluster dynamical states and the mass changes, and investigate the relaxation time scale -from the beginning of a disturbance to the final relaxed state (see more details in the following section).In Figure 5, we illustrate the evolution tracks of λDS and ∆M 200 M 200 over time for one arbitrary example cluster, where the variation ∆M200is estimated by M200 in the snapshot i minus M200 in snapshot i-1.The original behavior of the evolution track of λDS is highly jagged because of the frequent mergers along with difficulties in correctly tracking the progenitors in the simulation, so the function savagol_f ilter from SciP y.signal is applied to smooth the evolution curve (see scipy.org).In this function, the length of the filter window is set to be 11 data points, and the order of the polynomial to fit the sample is set to be 3.

The cluster relaxation time scale
In order to investigate the evolution of cluster dynamical states, we define the "relaxation period" to describe the time taken by a cluster to evolve from relaxed state to unrelaxed state and then return to relaxed states.As shown in Figure 6, one relaxation period starts with the local maximum of λDS above the threshold before decreasing, after which the λDS of cluster continues decreasing until reach some local minimum below the threshold.Then, the relaxation period ends with the first crossing point between λDS evolutionary track and the threshold, through which the cluster return to a relaxed state again.Note that we exclude the evolution track in the very beginning 4 Gyrs.This is because the halo is still very small and its dynamical state can be dramatically changed due to very frequent merging events.Our definition of this relaxation time scale is very similar to the merger time which is defined in Contreras-Santos et al. 2021.We share the same initial point to mark the start of relaxation time scale.However, Contreras-Santos et al. 2021 require the cluster returns to a following peak of the dynamical relaxation parameter for the end instead of the crossing of the threshold (our case).Besides that, they used χDS parameter to quantify the cluster dynamical state which is very similar to our λDS as shown in Figure 3. Therefore, we expect a similar scale between their merger time and our relaxation time.It worth noting that their studies focus on major merger events (∆M /M ≥ 0.5), while we will provide a more statistical view of the relaxation time scale.
As shown in Figure 6, one cluster can have more than one relaxation periods during its evolution process.The distributions for relaxation periods for clusters in samples are shown in Figure 7.The relaxation time scale is quantified as the median number of relaxation periods, which are 1.9 (1.8) Gyr, 1.6 (1.6) Gyr and 1.4 (1.6) Gyr for DMO run, GX run and GIZMO run respectively, the numbers inside brackets are standard deviations.The red brackets labels two relaxation periods identified in the evolutionary process of the cluster.The red horizontal line represent the threshold of cluster dynamical states, above and below which are relaxed and unrelaxed states respectively.

Connection to the halo mass changes
Relaxation time scale provides useful information about the evolution of cluster dynamical states, especially when connecting with the merger events.As mentioned above, it is intuitively correlated to the mass accretion history of cluster.
In order to investigate such correlation, we first select out those relaxation periods each with only one ∆M 200 M 200 peak inside, and make the scatter plots of relaxation time period vs the maxima of fractional halo mass change, ∆M 200 M 200 .However, no direct correlation can be found in this scattering plot.Then, we find that the relaxation period, after normalized with a dynamical time scale,t dyn , t relax /t dyn then shows a moderate linear correlation with ∆M 200 M 200 peak.The dynamical time scale t dyn is define as where Rvir and Mvir are virial radius and virial mass respectively.Here, we simply adopt R200 as the virial radius and M200 as the virial mass.
From Equation 15, it is easy to show that the dynamical time scale, t dyn , is only a function of critical density, ρc, which solely depends on redshift z.Hence, the dynamical time scale can be determined only with a given redshift.In this study, t dyn is determined from the redshifts at which the relaxation periods start.
The scatter plots of relaxation period/dynamical time scale vs ∆M 200 M 200 are shown in Figure 8.The Pearson correlation coefficient is 0.57/0.55/0.51 for DMO/GX/GIZMO run, which indicates a moderate correlation.Including more halo properties may give a better correlation, we retain that for a future study.Note that, as have discussed in subsection 4.3, it is not surprise to see a similar distribution of t relax/t dyn in Figure 8 compared to the Fig. 3  Then, we fitted these scatter plots with a linear function: and we obtain k = 12.047/10.706/9.794and h = 1.169/1.111/1.183from DMO/GX/GIZMO run, respectively.For the two hydro-dynamical simulations, we calculate their mean square fitting errors by Equation 14 and compare them with the mean square errors from DMO fitting function, i.e. data are from hydro-dynamical run, but predictions are made with Equation 16with parameters k and h yielded from DMO fitting.For GX run, the mean square error from DMO fitting (k = 12.047, h = 1.169) is 0.36, and that from GX fitting (k = 10.706,h = 1.111) is 0.34.For GIZMO run, the mean square error from DMO fitting is 0.34, and that from GIZMO fitting (k = 9.794, h = 1.183) is 0.33.The differences in mean square fitting errors from DMO fitting and hydro-dynamical fitting are small in both cases.Therefore, we simply use the values of k and h from the DMO fitting for all three simulations.
The distributions of fitting errors for relaxation periods with single ∆M 200 M 200 peak inside are showed in Figure 9.Most errors between predicted and real relaxation periods (89%/91%/89% for DMO/GX/GIZMO run) are less than ∼0.5 Gyr, which are considerably less than the median length of those relaxation periods with single ∆M 200 M 200 peak, which is 1.847/1.577/1.413Giga years for DMO/GX/GIZMO run.The median fitting errors of three distributions are close to 0, slightly deviate towards a positive direction.
Given the linear correlation shown in Figure 8, it is not surprising to see such a relatively small fitting error.To verify this fitting function, we adopt it for making predictions of the relaxation periods with more than one ∆M 200 M 200 peaks by simply linear summation of contributions from all ∆M 200 M 200 peaks: where t dyn is calculated by the redshift at which the relaxation period starts, and n represents the total number of ∆M 200 M 200 peaks that happen within the relaxation period.Note that a different t dyn for each peak may give a better prediction.The distributions of fitting errors for those relaxation periods with multiple peaks inside are showed in Figure 10.Most errors (82%/88%/88% for DMO/GX/GIZMO run) are less than ∼2 Gyrs.However, the median numbers of these distributions deviate towards the positive direction ( < ∼ 0.5 Gyr), which means that Equation 17 slightly overestimate the length of relaxation period.
The fractional fitting error distributions for DMO, GX and GIZMO run are plotted altogether in Figure 11.The histograms are not normalized, the total number of relaxation periods in two hydro-dynamical runs are significantly larger than that in DMO run, which implies an increased merger events by the baryon effect.81.3%/74.1%/71.2% of fractional errors in DMO/GX/GIZMO run are less than 0.6.In agreement with the behaviors in absolute error distributions, all fractional error distributions deviate towards positive direction.The deviation of the median number of fractional fitting error is strongest in GIZMO run, and the median number in GX run also has a larger deviation than that in DMO run.

CONCLUSIONS AND DISCUSSIONS
In this work, we use the mass-complete cluster sample from The Three Hundred to study the cluster dynamical states and proposed a new parameter λDS to classify the clusters into dynamical relaxed and unrelaxed without a manually set threshold.Benefited from the different runs (DMO, GX and GIZMO) within this project, we are also able to investigate the baryon effect on the cluster dynamical state.Furthermore, we define a relaxation time scale and connect it to the halo mass changes.Main findings are summarised below: • Based on the relaxation parameter χDS in Haggar et al. (2020), a new threshold-free function of λDS is proposed to classify cluster dynamical states, which is λDS = 3 (7.30× ∆r) 2 + (0.30 × fs) 2 + |1 − η| 2 (18) The threshold distinguishing relaxed and unrelaxed state is naturally set by the double-Gaussian fitting of the λDS distribution.At redshift z = 0, 151/171/170 clusters of all 324 clusters are classified to be dynamically relaxed in DMO/GX/GIZMO run.The λDS parameter is linearly correlated to χDS parameter, and it preserves the classification results based on χDS.
• Including baryons to simulations can slightly reduce the virial ratio η, which is 2% lower in GX and GIZMO run compared to DMO run.
Baryonic effect results in the scattering distribution of the center of mass offset, ∆r, the standard deviation of the difference between ∆r from GX/GIZMO run and DMO run is large (more than 50%) compared to the scale of ∆r from DMO run.
Sub-halo mass fraction fs is 17% higher in the GIZMO run than in the DMO run, while the GX run is about 2 per cent lower.
In combination, the λDS in GIZMO run is 3% lower than in GX run which has about 10 per cent higher value than the DMO run.Therefore, more relaxed clusters are presented in the hydrodynamic simulations.Nevertheless, the baryons play a weak role in altering the cluster dynamical state.
• The median number of relaxation periods (the time taken by a cluster to evolve from most relaxed state to unrelaxed state and then return to relaxed state) also regarded as relaxation time scale, is 1.913/1.610/1.419for DMO/GX/GIZMO run, respectively.
• The relaxation period is correlated to cluster mass accretion history.For relaxation periods with single ∆M 200 M 200 peak inside, a moderate linear correlation is observed, which is described as In general case, the length of relaxation periods can be predicted from the heights of ∆M 200 M 200 peaks with with a considerable small error, basically less than 2 Gyrs.
As shown in Figure 3, the new proposed λDS is basically linear correlated with the χDS.So it can be correlated with these observational measured quantities, such as M (Cialone et al. 2018;De Luca et al. 2021) and C (Capalbo et al. 2021) parameters.By applying its threshold from a double-Gaussian fitting, the clusters can be naturally separated into relaxed and unrelaxed.With this single and non-arbitrary classification, it is straightforward to define some time scale to describe the transition rate of dynamical states of a cluster, and such time scale can be determined completely from the features of the evolution track of λDS (see Figure 6), which makes it applicable to be analyzed statistically for large number of clusters, thus evaluate its overall correlation with other observables (e.g fractional mass change of cluster).
In this work, we only impose two constraints on the λDS parameter: having well-behaved Double-Gaussian distribution over clusters and preserving classification results with χDS.Meanwhile, the observed linear correlation between sub-halo mass fraction fs and center of mass offset ∆r is likely to introduce additional degrees of freedom in Double-Gaussian fitting.Therefore, we acknowledge there may be some other values of a and b in Equation 8, or even a different form of function to combine dynamical parameters together which can make λDS satisfy our requirements.In future work, it will be worthy to investigate the potential improvement of the formalism of λDS with some advanced statistical methods.
Although the baryons can affect the cluster properties in different aspects (see Cui et al. 2016, for example), the cluster dynamical state seems to be less influenced by the baryons.That is understandable as baryons have the strongest effect at very small scale, while the dynamical state describes the whole dynamical information of the cluster.This is similar to baryon effects on the total cluster mass (see Cui et al. 2012Cui et al. , 2014, for example), for example).Agreed with Zhang et al. ( 2016), the baryon effect does shirk the cluster's relaxation period, which results in slightly more relaxed clusters in the hydro-dynamical runs.However, different to their ideal case study which doesn't include any baryon processes in twohalos merger event, the hydro-simulated clusters from The Three Hundred project do not show significantly change in the relaxation time scale.This can be explained as in reality the merger speed and the gas content are relatively low, which is in agreement with their results -∼ 70 per cent reduction in the merger time scale.
Note that our definition of the cluster relaxation time scale is slightly different to the merger time which is widely used in the Semi-Analytical models (for example Boylan-Kolchin et al. 2008;Jiang et al. 2008;Jiang & van den Bosch 2014).Our definition focus on the overall cluster dynamical state, while the merger time mainly interest in the dynamical friction, for example, a satellite galaxy moving in a dark matter halo.The two time scales are very similar when a major merger happens.Moreover, by using the relation between the cluster dynamical state relaxation time scale with the cluster mass changes in this study, one can roughly predict how long will the cluster will get back to a relaxed state.
As the merger events can lead to the cluster/galaxy property changes.Contreras-Santos et al. 2021 (in review) using the cluster dynamical changes (similar to our relaxation time scale definition) to define pre-and post-merger phases, found that stellar content of BCGs grows significantly during mergers: the main growth mechanism is the accretion of older stars; there is a burst in star formation induced by the merger.Furthermore, the evolution of the hydrodynamic equillium bias can be also tightly connect to the major mergers (Gianfagna et al. 2021 in prep.).Therefore, though the observed accretion in mass, we can predict the cluster relaxation time scale which can be used to predict the changes of these quantities.
gramme under the Marie Sklodowskaw-Curie grant agreement number 734374, i.e. the LACEGAL project.
The simulations used in this paper have been performed in the MareNostrum Supercomputer at the Barcelona Supercomputing Center, thanks to CPU time granted by the Red Española de Supercomputación.The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP).The MultiDark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064.
WC is supported by the STFC AGP Grant ST/V000594/1.He further acknowledges the science research grants from the China Manned Space Project with NO.CMS-CSST-2021-A01 and CMS-CSST-2021-B01.For R500 data, a halo mass cut,M500 = 4.6e14 is applied to exclude low mass clusters.Then the same method is applied to the left 246 clusters, and the free coefficients for λ in eqn 8 are determined to be a = 15.85 and b = 1.04.The distribution of λDS for R500 is showed in Figure A1.
The threshold is λDS = 2.61, as the X coordinate of the crossing point of two Single-Gaussian functions.With this threshold, 101 in 246 clusters are classified as relaxed.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Distributions of the relaxation parameter, log 10 λ DS , for the mass-complete cluster sample from the DMO run, at redshift z = 0.The best fitting parameters are a = 7.30 and b = 0.30.Red line represent the fitted Double-Gaussian distribution.The two single-Gaussian functions are represented by orange and green line.
. Based on the value of χDS,Haggar et al. (2020) split the sample as relaxed clusters (χDS > 1.030), unrelaxed clusters (χDS < 0.619) and intermediate with 0.619 < χDS < 1.030.The two thresholds fromHaggar et al. (2020) are represented by two red vertical lines in Figure3.In our work, a single threshold, λDS = 3.424 is determined from a systematical way, which is represented by the green horizontal line.The correlation between λDS and χDS is almost linear, and the most clusters classified as relaxed or unrelaxed byHaggar et al. (2020) have a similar classification with our parameter.This means although our new relaxation parameter λDS adjusts the relative contributions between dynamical parameters η, fs and ∆r to re-scale and re-distribute the old relaxation parameter χDS, it is still monotonic correlated with old one and almost doesn't qualitatively change the results from previous work on classifying cluster dynamical states.Whilst, λDS can provide a single and non-arbitrary threshold from its special Double-Gaussian distribution, which makes the classification of cluster dynamical states straightforward and clear.

Figure 2 .
Figure 2. Distributions of relaxation parameter, λ DS , in log10 scale (the same as Figure 1) for 324 clusters at redshift z = 0, in GX run (a) and GIZMO run (b).Red lines represent Double-Gaussian fit, single-Gaussian fits are represented by orange and green lines in each plot.

Figure 3 .
Figure 3. λ DS vs χ DS in logarithm scale for 324 clusters, at z = 0, in GX run, which Haggar et al. (2020) used.Two red vertical lines represent the threshold on χ DS for relaxed clusters and unrelaxed clusters, which are χ DS = 1.030 and χ DS = 0.619 respectively.The green vertical line represents the threshold on λ DS , which is λ DS = 3.424.

Figure 4 .
Figure 4. Distributions of the relative differences for (from top to bottom) η, ∆r, fs, λ DS between GX (blue histograms) or GIZMO (orange histograms) run and DMO run.Green (red) vertical lines represent the median number of differences for GX (GIZMO) run.

Figure 5 .
Figure 5. Evolution tracks of λ DS (orange) and ∆M 200 M 200 (green) over time for the 5th cluster, in (a) DMO run; (b) GX run; (c) GIZMO run.The red horizontal line represents the threshold λ DS = 3.424.The region above this line represents the cluster in a dynamically relaxed state.

Figure 6 .
Figure 6.λ DS evolutionary track for the 5th cluster in DMO run.The red brackets labels two relaxation periods identified in the evolutionary process of the cluster.The red horizontal line represent the threshold of cluster dynamical states, above and below which are relaxed and unrelaxed states respectively.

Figure 7 .
Figure 7. Relaxation period distribution for all clusters, in (a) DMO run; (b) GX run; (c) GIZMO run.Blue bins represent all relaxation periods, orange bins represent relaxation periods with only one ∆M 200 M 200 peak inside.
Figure 8. t relax t dyn vs ∆M 200 M 200 , for relaxation periods with single

Figure 9 .
Figure 9. Fitting errors of relaxation periods with single ∆M 200 M 200 peak inside the duration, for (a) DMO run; (b) GX run; (c) GIZMO run.Red vertical line represents the median number of fitting errors.

Figure 10 .
Figure 10.Prediction errors of relaxation periods with multiple ∆M 200 M 200 peaks inside the duration, for (a) DMO run; (b) GX run; (c) GIZMO run.Red vertical lines represent the median prediction errors.

Figure 11 .
Figure 11.Distributions of the relative prediction errors, (t relax,predict -t relax,real )/t relax,real , for clusters in DMO(filled blue bins), GX(unfilled orange line) and GIZMO(unfilled green line) run.Blue/Orange/Green vertical line represents the median number of prediction error from DMO/GX/GIZMO run.

Figure A1 .
Figure A1.Distributions of the relaxation parameter, log 10 λ DS , for the mass-complete cluster sample from the DMO run, R500, at redshift z = 0.The best fitting parameters are a = 15.85 and b = 1.04.Red line represent the fitted Double-Gaussian distribution.The two single-Gaussian functions are represented by orange and green line.

Table 1 .
P value , coefficients c, µ, σ, square error E for the fitting result with the best parameters: a = 7.30 and b = 0.30.