{'fraction_non_alphanumeric': 0.03302795031055901, 'fraction_numerical': 0.02020186335403727, 'mean_word_length': 4.649210526315789, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 9, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Machine Learning methods to estimate observational properties of galaxy clusters in large volume cosmological N-body simulations

In this paper we study the applicability of a set of supervised machine learning (ML) models specifically trained to infer observed related properties of the baryonic component (stars and gas) from a set of features of dark matter only cluster-size halos. The training set is built from THE THREE HUNDRED project which consists of a series of zoomed hydrodynamical simulations of cluster-size regions extracted from the 1 Gpc volume MultiDark dark-matter only simulation (MDPL2). We use as target variables a set of baryonic properties for the intra cluster gas and stars derived from the hydrodynamical simulations and correlate them with the properties of the dark matter halos from the MDPL2 N-body simulation. The different ML models are trained from this database and subsequently used to infer the same baryonic properties for the whole range of cluster-size halos identified in the MDPL2. We also test the robustness of the predictions of the models against mass resolution of the dark matter halos and conclude that their inferred baryonic properties are rather insensitive to their DM properties which are resolved with almost an order of magnitude smaller number of particles. We conclude that the ML models presented in this paper can be used as an accurate and computationally efficient tool for populating cluster-size halos with observational related baryonic properties in large volume N-body simulations making them more valuable for comparison with full sky galaxy cluster surveys at different wavelengths. We make the best ML trained model publicly available.


INTRODUCTION
Galaxy clusters are the largest gravitationally bound objects of the Universe and constitute one of the best cosmological probes to constrain cosmological parameters of the Universe. The main component of galaxy clusters is dark matter (DM), which accounts for 85 per cent of the total mass (for a full review see e.g. Allen et al., 2011;Kravtsov & Borgani, 2012). Although the existence of DM is now widely accepted by the scientific community and strongly supported by modern cosmological theories, it has not been directly detected yet. To study galaxy clusters, we have therefore to focus on their baryonic component, which represents the remaining 15 per cent of the mass. It is composed by the hot gas of the Intra Cluster Medium (ICM, around 10-15 per cent of the total cluster mass) and stars (less than 5 per cent of the mass).
Numerical simulations play a fundamental role to study the properties of galaxy clusters. In the simplest scenario, N-body simulations can easily describe the dark-matter component of clusters, which is governed only by gravity; nowadays it is computationally possible to perform very large cosmological simulations, up to a few Gpc 3 , e.g. M XXL (Angulo et al., 2012), MICE (Fosalba et al., 2015), M D (Klypin et al., 2016), D (Skillman et al., 2014), O R (Habib et al., 2016), FLAGSHIP (Potter et al., 2017), U (Ishiyama et al., 2021), BACCO (Angulo et al., 2021) and UNIT project (Chuang et al., 2019), which include thousands of galaxy clusters. Nevertheless, when aiming to describe the baryon component of clusters, due to the complex physics involved in the processes of cluster formation, radiative hydrodynamic numerical simulations have to be used. These simulations are computationally very expensive so this puts strong limitations to the size of the computational volumes. Examples of state-of-the art of such simulations are: I (Vogelsberger et al., 2014), E - (Schaye et al., 2015), H -AGN (Chisari et al., 2016), M (Dolag et al., 2016) or BAHAMAS (McCarthy et al., 2018). Hydrodynamical simulations are essential to calibrate mass proxies and to study the systematics affecting observational measurements. They are also essential to deeply understand the formation and evolution of clusters of galaxies and all their gas-dynamical effects. For this reason, numerical simulations have been a powerful tool to guide galaxy clusters observations for more than 20 years (Evrard et al., 1996;Bryan & Norman, 1998).
In an ideal scenario one would need to have a large sample of simulated galaxy clusters with enough numerical resolution, both in mass and in the gravity and pressure forces. This high resolution would allow to accurately resolve the internal substructures and to obtain a detailed modelling of the most relevant physical processes. The best way to achieve this would be by simulating large cosmological boxes containing up to tens of thousands of galaxy clusters. Unfortunately, due to the large computational effort demanded by these simulations, one needs to find a compromise between their three main components: volume size, mass resolution and physical processes included. A possible solution to the computational problems related with scalability of present-day hydrodynamical codes is to proceed with the so-called 'zoom' simulations, such us the MUSIC 1 simulation (Sembolini et al., 2013), the D clusters (Planelles et al., 2013), R -G (Wu et al., 2015), MACSIS (Barnes et al., 2016), C -EAGLE (Barnes et al., 2017), - (Bahé et al., 2017) clusters and T T H (The300) 2 simulation project (Cui et al., 2018). Zoom simulations are performed mimicking the observations, by creating a catalogue of resimulated galaxy clusters that are extracted from low-resolution N-body simulations. The regions containing clusters of galaxies are then resimulated at very high resolution, adding gas physics in the resimulated areas and keeping the rest of the box at low resolution in order to reproduce the same gravitational evolution.
An alternative approach to hydrodynamical simulations to describe the gas and stellar properties of galaxy clusters, is to use Semi-Analytic Models (SAMs), such us GALACTICUS (Benson, 2012), SAG (Cora et al., 2018), SAGE (Croton et al., 2016) and GALFORM (Lacey et al., 2016). In this approach, the numerous complex non-linear radiative physical processes associated to the 1 https://music.ft.uam.es 2 https://the300-project.org gas-star components are modelled using a combination of analytic approximations and empirical calibrations of many free parameters against a set of observational constrains (see e.g. Baugh (2006) for a review). Nevertheless, SAMs are also computationally expensive since most of them are based on the information provided by merger history of each individual dark matter halo. A complementary approach is the use of phenomenological models to derive physical properties of the ICM as in Zandanel et al. (2018) and Osato & Nagai (2022). Describing the gas physics in simulated galaxy clusters requires therefore a big computational effort and impose a compromise between numerical resolution and size of the cosmological volume to simulate.
The main goal of supervised Machine Learning (ML) is to generate models that can learn complex relationships between input and output variables from high-dimensional data that can later be used to make predictions on unseen data. In this scenario, ML could offer a powerful alternative to infer some fundamental information on the main properties (e.g. gas and star masses, gas temperature, etc) of the baryon component of galaxy clusters, without the large computational cost required by hydrodynamical simulations or SAMs. Applications of ML to find a mapping between hydrodynamical and N-body simulations have been already presented in previous works. Firstly, in Kamdar et al. (2016), a promising technique to study galaxy formation using numerical simulations and ML was presented; Jo & Kim (2019) estimated galactic baryonic properties mimicking the I TNG simulation (Nelson et al., 2019); Wadekar et al. (2021) generated neutral hydrogen from dark matter; Bernardini et al. (2022) predicted high resolution baryon fields from dark matter simulations; Moews et al. (2021) used hybrid analytic and machine learning model to paint dark matter galactic halos with hydrodynamical properties; Lovell et al. (2022) explore the halogalaxy relationship in the periodic EAGLE simulations, and zoom C-EAGLE simulations of galaxy clusters; and McGibbon & Khochfar (2022) consider a ML model that is built using the extremely randomised tree (ERT) algorithm and takes subhalo properties over a wide range of redshifts as its input features for galaxy scales. Recently, The CAMELS collaboration (Villaescusa-Navarro et al., 2022) has released results from almost ten thousands simulations (both hydrodynamical and N-body) with different cosmologies and baryon physical models that are an invaluable tool for training current and future Artificial Intelligence algorithms that will be very useful for galaxy formation studies. Unfortunately, given the box sizes, the number of cluster-size objects is poorly represented in these simulations.
The purpose of this study is to explore the applicability of ML techniques to generate baryon cluster properties from DM-only halo catalogues mimicking the results from T T H hydrodynamical simulations. More precisely, we use the properties of the cluster-sized halos extracted from parent dark matter only full box simulation MDPL2 as the features of our dataset. Then we collect several baryon properties of the objects that have been re-simulated with radiative processes and hydrodynamics as targets (the predicted variables) of the ML models. Our work differs from previous studies in that the baryon properties are extracted from 'zoom' MDPL2-based simulations and therefore, we have paired one to one the objects corresponding to the full N-body only simulations with their hydrodynamical counterparts. As explained below, The300 simulations corresponds to spherical regions centred on the 324 most massive clusters found in the MDPL2 box. But there are more cluster-size halos found within each region with lower masses. The masses of the cluster-size catalogue of hydrodynami-cal simulated objects we are using ranges from ∼ 10 13 ℎ −1 M up to ∼ 10 15 ℎ −1 M .
The article is structured as follows: In § 2, we describe how the training dataset is generated using The300 and the MDPL2 simulations. In § 3, we explain the different ML algorithms used in this work and the training setup. We also study the feature importance and selection of our feature space. In § 4, the main results for this work are shown, including an analysis of the performance of the ML models and their dependence on mass resolution of the simulations. In § 5, we study the scaling relations extracted from the new ML-generated catalogues and finally in § 6, we draw our main conclusions and propose possible future studies.

THE TRAINING DATASET
In order to create the database for training the ML models, we use the MDPL2 3 simulation, which has been run using the cosmological parameters measured by the Planck Collaboration (Planck Collaboration et al., 2016). The MDPL2 simulation consists of a periodic cube volume of comoving length 1 ℎ −1 Gpc containing 3840 3 dark-matter particles each with a mass 1.5 × 10 9 ℎ −1 M .
To build this training dataset, we need first to identify and extract from the MDPL2 simulation the same cluster objects that were used to run zoomed The300 hydrodynamical simulations. We then select the main properties of the dark matter clusters and associate them with the baryonic properties extracted from their The300 hydrodynamical counterparts.

MDPL2: Dark Matter input variables
In order to identify the dark matter halos and measure their internal properties in the MDPL2 N-body simulation we have used the R halo finder (Behroozi et al., 2012), complemented with additional information based on the halo mass accretion history from the C H M T analysis (Behroozi et al., 2013). We have extracted a total of 26 relevant physical R -+ C T variables 4 (masses at different radii, velocities, symmetry factors, properties related with mass accretion history, etc) to create our dark matter catalogue. In addition, we have also considered the scale factor ( ) of clusters as an input variable. Furthermore, we have introduced a cut-off in halo mass such that log( /( ℎ −1 M )) ≥ 13.5 and redshift ≤ 1.03.
In Fig. 1, we show the Spearman correlation matrix of the 26 R variables and the scale factor ( ). These variables are ordered using a hierarchical clustering algorithm based on Ward's linkage on a condensed distance matrix. We used the Python implementation of this algorithm from S P (Virtanen et al., 2020). We can easily identify 5 groups in the correlation matrix. The first group (variables 0 to 12) corresponds to masses and velocities at different radii. In a second group, different ellipticity shape factors (from 13 to 16) are included. Variables from 17 to 21 corresponds to the scale radius, the ratio between the kinetic and potential energy and the offsets between density peak and centre-of-mass, which are directly related to the dynamical state of the cluster halos. The next group of variables (22 and 23) correspond to the dimensionless spin parameters of the cluster. Finally, variables from 24 to 26 represent 3 www.cosmosim.org 4 More information regarding the selection of R variables can be found in Appendix A the scale factor (redshift) and the time evolution of mass accretion. As can be seen in the figure, feature variables inside the same block are strongly correlated among them and they are weakly, or not correlated to variables inside other blocks. This might imply that selecting more than one feature belonging to the same block could not add any new predictive information. This is studied in detail in section § 3. A more detailed description of the selected feature variables can be found in the Appendix A.

The300: baryonic output variables
Subsequently, for a subset of the MDPL2 cluster halos, we need to have their baryonic properties. For this purpose, we have used the results of The300 project, which has re-simulated spherical regions of radius 15ℎ −1 Mpc centred around the 324 most massive clusters found in the MDPL2 simulation at = 0. These regions were then mapped back to the initial conditions and their particles were split into gas and dark-matter, while the rest of the particles in the remaining box were re-sampled into different levels of lower resolution and larger masses. With this zoom-in technique, we ensure that the subsequent gravitational evolution will reproduce the same objects in the high resolution area while we minimise the effects of contamination of low resolution particles from external regions due to mass segregation. In any case, we checked that all the clusters used in this work are free from contamination of low mass resolution particles at least within their virial radii.
The300 project has produced different versions of hydrodynamical simulations from these initial zoomed conditions which include different baryonic physics modules: radiative cooling, star formation and Supernovae Feedback using the G -MUSIC SPH+TreePM code (Sembolini et al., 2013) and newer versions that include feedbacks from Super Massive Black Holes: G -X (Murante et al., 2010;Rasia et al., 2015), GIZMO-SIMBA (Davé et al., 2019;Cui et al., 2022).
However, in this work, we only make use of the G -X runs. The halos in these simulations are identified and analysed with the A H F (AHF) (Knollmann & Knebe, 2009), which is more suitable than R for simulations with multiple particles species (i.e. dark matter particles, gas, stellar particles and Black Holes). From the information contained in the AHF catalogues, we have collected the following baryon properties: • The total gas mass gas inside a spherical volume with a mass density that is 500 times larger than the critical density of the Universe at each redshift. The radius of this sphere is denoted as 500 .
• The Stellar mass star inside 500 .
• The gas temperature gas computed as the mass weighted temperature, inside 500 where and are respectively the temperature and mass of the gas particle. We additionally made a cut of > 0.3keV to exclude low-temperature gas particles.
• The X-ray Y-parameter X defined as gas × gas , which related with the total thermal energy of the gas and it has been shown that it is a good proxy of the total cluster mass . Note that this quantity can be derived from others. However, we prefer to treat it as an independent target, i.e the ML models are also trained to predict X as one of the target variables.
• The integrated Compton-y parameter SZ over 500 given by the Sunyaev-Zel'dovich (SZ) effect (Sunyaev & Zeldovich, 1972). Particularly, the integrated value SZ is computed from Compton-y parameter maps estimated as in the following: where T is the Thomson cross section, B is the Boltzmann constant, the speed of light, e the electron rest-mass, e the electron number density, e is the electron temperature and the integration is done along the observer's line of sight. Assuming = , Eq.
(2) is computed in our simulated data as in Sembolini et al. (2013) and Le Brun et al. (2015): Note that here we have used the number of electrons in the gas particles e given that e = e / . Moreover, ( , ℎ i ) is the same SPH smoothing kernel as in the hydrodynamical simulation with smoothing length ℎ i . The -maps are generated with the centre on the projected maximum density peak position of the halo. Each image has a fixed angular resolution of 5 that is extended to at least 200 in all the clusters. The clusters at = 0 are placed at = 0.05 to generate the mock images while the clusters at higher redshifts simply use its original value from the simulations. We then integrate the Compton-y map up to 500 using only the zplane projection. Since the dataset is large, the effect of projections is negligible. Note that this approach of estimating SZ gives us the cylindrical Compton-y parameter 500,cyl . In practice, 500,cyl is related to the corresponding spherically integrated value 500,sph such that 500,cyl / 500,sph 1.2 as given by the Arnaud et al. (2010) pressure profile and also compatible with results from numerical simulations (Sembolini et al., 2013).
Note that star forming gas particles with > 0.1 yr −1 are also excluded in the calculation of gas X and SZ . This is commonly adopted in simulations with a multiphase subgrid physics model in which gas particles are composed of hot and cold gas components (see e.g. Borgani et al. (2004)). These star forming gas particles are poorly modelled as they are mostly made by the neutral cold gas component which should not emit any X-rays or contribute to the total electron thermal pressure. 5

The Final Training Dataset
After defining our input and output variables, we finally match oneby-one the clusters between the two simulations that fulfil these two conditions for the relative shifts between the cluster centres and the halos mass differences: MDPL2,200 The300,200 − 1 < 0.1 .
Here, MDPL2 and The300 stand for the centre of mass of the clusters while MDPL2,200 and The300,200 stand for the mass inside a sphere of radius = 200 for each simulations (between DM 5 We have verified that changing to > 0 makes no difference in the results due to a very small number of star forming gas particles in galaxy clusters(see Li et al., 2020). only R catalogue and the AHF catalogue respectively). Due to both the baryon effect (see Cui et al., 2012Cui et al., , 2014 and to different algorithms used by the halo finders, it is not possible to determine with all certainty that all the halos are exactly matched. Notice that the centre difference can be as high as 0.4 200 . However, with this restrictive selection criteria, only the true/very close counterparts are selected. In this way, we finally provide the baryon properties for the matched MDPL2 clusters using the corresponding The300 objects.
After this procedure, our dataset is finally composed of 49540 different objects. Note that all the 33 halo catalogues available from = 0 to = 1.03 in the two simulations have been considered. Only 1264 objects correspond to clusters at = 0, the rest of them are the progenitors of the same objects at different redshifts. The number of objects as a function of their mass and redshift can respectively be found in Fig. 2 and Fig. 3. Our final dataset is composed of 27 DM input variables and 5 baryon output variables. These are the features and targets which are used for training and testing the ML algorithms described in the next section.

MACHINE LEARNING ALGORITHMS: DESCRIPTION AND TRAINING
In this section, we first describe the machine learning algorithms used in this work and the training setup. Then, we study the importance of our feature variables in order to reduce the dimensionality of our dataset.

Machine Learning Algorithms and Training Setup
In order to estimate the baryon properties of the dark matter only clusters, several effective supervised machine learning methods have been employed. We particularly focus on four methods: random forest (RF; Breiman, 2001), extreme gradient boosting (XGBoost; Chen & Guestrin, 2016), dense Neural Networks or Multilayer Perceptron (MLP; Schmidhuber, 2015) and Natural Gradient Boosting for Probabilistic Prediction (NGBoost; Duan et al., 2019). RF and XGBoost have shown to be among the best machine learning methods for tabular data (i.e. without a known grid-like topology, such as images) (Fernández-Delgado et al., 2014;Bentéjac et al., 2021;Zhang et al., 2017). Convolutional deep neural network models have shown spectacular performance for image-based and structured data in general (Schmidhuber, 2015). However, for tabular data, as is the case of this study, their performance is poor (Zhang et al., 2017). Notwithstanding, deep dense networks can perform well in these scenarios, so we will also consider these models. Random Forest and XGBoost are ensemble models composed of decision trees. During training, these algorithms build hundreds of decision trees from a single training dataset. The process for building these trees in random forest and XGBoost, is based on quite different ideas. Although, the objective in both cases is to build decision tree models that complement each other in order to obtain a classification/regression model better than any of its parts (Dietterich, 1998).
Random forest rely on stochastic techniques to generate many random solutions to the problem at hand. In order to generate each single tree, the random forest algorithm first generates a new dataset by extracting at random instances of the training data of size with replacement (i.e. bootstrap sample). This bootstrap sample is used to train a decision tree in which the best split at each node of the tree is selected from a random subsample of features of the data.  (1) M500c (2) M200c (3) Mpeak (4) mvir (5) Macc (6) Vpeak (7) Vmax\@Mpeak (8) vmax (9) Vacc(10) rvir(11) vrms(12) b_to_a(500c)(13) c_to_a (500c)(14) b_to_a (15) c_to_a (16) rs (17) Rs_Klypin (18) T/|U| (19) Xoff (20) Voff (21) Spin (22) Spin_Bullock (23) a(24) scale_of_last_MM (25) Halfmass_Scale (  Generally, the size of the random subset of features is of the order of √ or log 2 ( ), with the number of features of the problem. The final output of the random forest for a given instance is obtained as the mode or mean of all trees for classification and regression respectively. In addition, since the randomisation process to build  the trees is independent, the process of building a random forest can be easily parallelised.
On the other hand, XGBoost relies mainly in a gradient descend approach although it also incorporates stochastic techniques to further increase its performance. XGBoost is an additive model based on Gradient Boosting. The output of an additive model is the sum of the outputs of its components. In order, to create this ensemble, regression trees are trained sequentially to approximate the gradient of the loss function of the data in the previous iterations. Hence, each new tree learns the remainder of the concept not learned in previous steps. XGBoost also includes a penalisation term in the number of leaves of the trees to avoid over-fitting. In addition, XGBoost incorporates random feature selection, bootstrap sample and several other randomisation features.
In order to perform a fair comparison among algorithms and also to obtain good estimations of the performance of the different algorithms, we carried out the following experimental procedure based on K-fold cross-validation and grid-search. K-fold crossvalidation consist in splitting the data into K disjoint sets of approximately equal size and then to use iteratively − 1 sets for training the model and the remaining set for validation. The main experiment is performed using the same 10-fold cross-validation for the prediction of the five baryonic properties analysed in this study using the R halo catalogue from The300 hydro clusters. The steps for each of the 10 partitions of the cross-validation are: (i) Find the best hyper-parameters of each of the tested algorithms: RF, XGBoost and MLP. For that, a grid-search with 5-fold crossvalidation within the train dataset only was performed. The values for the grid of hyper-parameters are shown below; (ii) The best set of hyper-parameters for each method were used to train a single model using the whole training set; (iii) The models were validated using the test set; In order to generate dark matter only halo catalogues with hydrodynamic properties, the 10 trained models from each of the 10-folds of the cross-validation were used. The hydrodynamic features of each halo are then computed as the average of the inferred values from these 10 models.
For the grid search the set of values of the tested hyperparameters for each of the analysed methods are: Furthermore, MLP has been trained for 500 epochs or until the training loss is constant during 20 epochs. For more information of these hyper-parameters, we refer the reader to the Python libraries used throughout this work: for RF and MLP we have used scikit-learn 6 (Pedregosa et al., 2011) and for XGBoost its own library 7 . For the hyper-parameters not considered in the search grids, their defaults values were used.
In order to train these models the mean squared error of the logarithmic values of the targets (logarithmic MSE) was used as the loss function: where true, is the true value of the target extracted from the The300 simulation and <log pred, | > is the predicted target's value by our model. Note that, since the model is trained with log as targets, then the prediction of the model is directly the logarithm of the given target. In addition, corresponds to the number of objects in the dataset (e.g. train, validation and test set) where the MSE is computed. For numerical reasons, we have also used the logarithmic value of the features during the training process. Also note that these models infer only one prediction for every input value . These predictions are considered as the mean predictions without their statistical uncertainties. Moreover, Stiskalek et al. (2022) showed that by accounting for a proper modelling of uncertainties ML models can successfully mimic the statistics of the data, i.e. not only the mean but also the scatter. Moreover, modelling uncertainties using ML algorithms is a topic of recent studies (e.g Kodi Ramanah et al., 2020Ramanah et al., , 2021Ho et al., 2021;Eisert et al., 2022;de Andres et al., 2022) In order to address this issue we apply a generalisation of gradient boosting for probabilistic modelling: NGBoost. NGBoost is a gradient boosting algorithm as XGBoost, but that is based on assuming a particular parametric ( ) probability distribution ( | ). Thus, the loss function can be written in terms of these parameters as the negative log-likelihood: In our work we assume a Gaussian distribution with parameters mean =< log pred, | > and standard deviation . Therefore, the last equation can be written as Note that for every input example NGBoost predicts a mean value and its scatter which is given by the width of the Gaussian . NGBoost also implements the generalised natural gradient when minimising the loss function.
In order to tune NGBoost, we follow a procedure similar to the one proposed in the original paper. For every K-fold, we find the best value of the number of estimators by cross-validation. Moreover, for all the experiments the base learners are decision trees with a maximum depth of three levels and the learning rate is set to 0.01. In summary, the hyper-parameters for NGBoost are: • NGBoost -'n_estimators'= [best out of a maximum of 3000] -'max_depth'= 3 -'learning_rate'= 0.01 All other hyper-parameters were set to their default values. Note that NGBoost is the only model that can generalise the intrinsic scatter: where is the noise due to the scatter. Previous models (RF,XGBoost, MLP) are only designed to infer the mean value of < log | >. Moreover, NGBoost predictions on unseen data are computed by averaging over the 10 models' distributions corresponding to our 10 different K-folds:

Feature importance and selection
Although machine learning models can generalise complex functions, generally, it is not trivial to interpret their decisions. In fact, they are often referred as black box estimators (e.g. Barredo Arrieta et al., 2020). Therefore, it is of great value to be able to inspect what is the learnt relation between features and targets given a particular model. One such inspection technique is feature importance. Particularly, feature importance is a family of techniques that assigns a score ( , ) to each input features depending on how useful they are when it comes to predicting a particular target . Furthermore, feature selection is a family of techniques that aim at getting rid of non-informative variables from a model (e.g. Kuhn et al., 2013). In this section, we use a feature importance algorithm to determine what features are more relevant and therefore, reduce the dimensionality of our 27-dimensional input space. One commonly used algorithm to estimate feature importance for ensembles of decision trees (such as RF and XGBoost) is Permutation Importance (Breiman, 2001). In this algorithm, the importance of each feature is estimated as the decrease of the model score when the values of a feature are randomly shuffled. This technique, however, fails when correlated features are present in the dataset (Altmann et al., 2010). A second shortcoming of this algorithm is that it only considers the importance of individual features.
Other technique is the use of forest of trees to evaluate the importance of features computed as the mean and standard deviation of accumulation of impurity decrease within each tree (Breiman, 2001), which for regression is the variance reduction. In random forest, internal node features are selected with some criterion, or loss function. We can then measure how on average each feature decreases the criterion in the splits of the forests. Nevertheless, this technique also fails due to the fact that our features are highly correlated, and it is also known to be biased in favour of variables with many possible split points (e.g Nembrini et al., 2018).
Instead, we use the Greedy Search Feature Importance Algorithm (GSFIA, see for example Ferri et al. (1994). This technique considers the importance of the combination of features and not only the individual feature importance. It works iteratively by selecting and evaluating one variable at a time until all features are ordered from the most to the least relevant. The algorithm works with a list of selected variables, , initially empty, = [], and a pool of possible variables to be selected, , initially containing all variables of the problem, = ( 1 , . . . , ). Then, a procedure is repeated times in which, at each step, one variable from the pool is selected and moved to the list . In the − th step of the loop the procedure creates | | models trained on all of the features in plus one feature from . The model that minimises the MSE identifies the most important variable from in combination with the variables in . This variable is then removed from and appended to the list . At the end of the algorithm, all variables of the problem are sorted by importance in list together with the loss function associated with them. GSFIA it is depicted using pseudocode in Algorithm 1.
With this algorithm, we can define the feature importance score ( , ) as follows: • Run GSFIA to rank all features from the most to the least important variables and save the corresponding value of MSE.
• The score is then defined as the MSE of every iteration normalised to the corresponding value of the first iteration.
Note that the normalised MSE will be 1 for the first feature, and will decrease progressively as we consider more features until it converges to a minimum value. It could happen that after including several features, the normalised MSE increases as more features are included (see Fig. 4 for the case of star ). This indicates that the last features included do not improve or even degrade the performance of the model. The algorithm was run using random forest as model (line 9 of Algorithm 1). In addition, due to the randomness of the ML model, the inner loop of the algorithm was repeated 10 times in order to reduce the variability of the results. In Figure 4, the average of the normalised MSE and its standard deviation are shown for the different targets considered. In the horizontal axis, the final order for the feature variables is shown. Variables in red colour are the reduced set of features that will be considered for further analysis. These features are summarised in Table 1.
As shown in Table 1, we expect that the selected variables generally come from different correlation blocks, as shown in Fig. 1. This is so, since variables from the same block are correlated and  (2), Vpeak (7), scale_of_last_MM(25), Macc (6), a (24) star Vpeak (7), Halfmass_Scale (26), scale_of_last_MM(25), a (24) gas Vpeak (7), a(24), scale_of_last_MM(25) X Vpeak (7), M500c (2), scale_of_last_MM(25), vrms (12),Macc (6) SZ Vpeak (7), Mpeak(4), scale_of_last_MM(25), a(24), rs (17) once the algorithm chooses one feature, it skips using variables with the same information. However, this is not always the case (e.g. variables 2, 6 and 7 are selected for a couple of targets). This can be explained since the correlation between those variables is high, but it is not 1. Hence, for our case the marginal information that a second variable inside a correlated block gives, is higher than that given by other variables. As far as the meaning of selected variables is concerned, we can distinguish two different important blocks in the correlation matrix: The mass and velocity block (the first block from 0 to 12), and the time evolution block (from 24 to 26). The conclusion of this analysis is that the rest of the blocks are redundant or contribute little to the estimating baryon properties, i.e. the ellipticity block (from 13 to 16), the dynamical state block (from 17 to 21) and the spin block. Moreover, masses and velocities are the most important features for estimating baryon properties while the variables associated with the time evolution of the mass accretion into halos play a secondary role in the regression algorithms. The redundant role of the ellipticity variables can be explained by taking into consideration that we are estimating integrated quantities from the particles within spheres of radius = 500 , regardless of the shape of their 3D distributions. Note that, we combine data from different redshifts as our training and test samples. We do not think that the evolution of these baryon properties will affect our results because (1) as shown in Cui et al. (2022), these quantities in G -X simulations hardly depend on redshift, especially at 1 (see also Truong et al., 2018, for example); (2) we also include the scale factor as a feature variable in the training. If there were a clear redshift dependence on any target variable, the scale factor feature would show a higher contribution. However, as shown in Fig. 4, the scale factor contributes only weakly to the normalised MSE.
Furthermore, we have to highlight that although we have used Random Forest for the GSFIA, other Machine Learning algorithms might also be used. However, GSFIA is computationally expensive given the fact that its computing time increases with the number of features as ( 2 ). Therefore, we prefer to use RF because it is computationally more efficient and it does not have as many hyperparameters to tune in. Consequently, this choice might introduce a bias given the fact that a particular model is being used for the selection of the important variables. However, in the next section we will show that this particular selection of variables yields similar performance for the different ML algorithms considered throughout this work.

RESULTS
In this section, we first study what machine learning algorithm is of higher quality for our particular dataset and study the accuracy of our model predictions. Then, we populate the dark-matter-only MDPL2 simulation with baryon properties and determine whether  26 12 17 3 4 22 6 11 16 15 18 13 1 0 8 19 23 5 2 14 21 9 20 Table 1. we can also successfully use the trained machine learning model on dark-matter-only low resolution simulations.

Error analysis
In order to determine the accuracy of our ML models, we have trained our four models on the dataset composed of all features and on the dataset with the reduced set of features summarised in Table 1 using the experimental setup described in the previous section. The average performance of the models is shown in Fig. 5 Table 1.
As a general result, it can be observed from Fig. 5 that XG-Boost algorithm has the best performance for all targets. For RF, we find equivalent performances for both sets of features in gas , X and SZ ; a somewhat worse result for the reduced set on gas ; and better performance on star for the reduced set. For XGBoost, the trends are similar to those of Random Forest, although the difference in performance for XGBoost between both sets of features is negligible for gas and smaller for star . For the MLP model, all results using the reduced set of features are worse than those obtained when using all the features in the catalogues. These differences between the tree based approaches (RF and XGBoost) and MLP can be explained taking into consideration that the selection of important features was done using Random Forest. In any case, the performance of MLP is the worst for all targets even when all features are considered. For NGBoost, the MSE of mean predictions is similar in both set of features and worse than XGBoost and RF. After the previous analysis, we can conclude that XGBoost gives the most accurate model predictions. However, as shown later, NG-Boost is more accurate than XGBoost when it comes at mimicking the scatter of the true The300 data. Therefore, we will only consider XGBoost and NGBoost algorithms for the rest this work. A summary of the performance for all models can be found in Table 2 for the reduced set of features.
The scores shown in Fig. 5 summarise in a single value the performance of the models. However, they do not allow us to understand how the model performs in the different regions of the space of features and targets. In order to analyse this, we first define the relative difference in performance for a single target y as diff(y) = 10 <log pred, | > − true Note that in Eq.(12), we are not considering the logarithmic value of the targets, but the model aims at predicting the logarithmic values. One can interpret these differences as a probability distribution. This means that given a value of pred one might estimate the aleatory scatter to that particular predicted value. These differences are shown in Fig. 6 as a function of the predicted target pred (first column), the cluster mass 500 (second column) and the peak of the velocity profile along the mass accretion history, peak (third column) for all redshifts. In Fig. 6, instead of plotting the individual differences for all instances, the mean value (dashed black) and the 66% (red region) and 95% (blue region) confident intervals are represented for sliding windows (bins) containing roughly the same number of objects.
The main result that can be observed from Fig. 6 is that the predictions are unbiased with respect to the most important features ( 500 and peak ) and with respect to the predicted targets, since the mean is very close to 0 for all ranges. However, the scatter varies depending on the target as it is depicted in Fig. 5 and Fig. 6. Particularly, gas is the target most accurately predicted, with an average scatter of 7% (standard deviation of Eq. (12)) and X is the predicted variable with higher average scatter ( 16%). The numerical values corresponding to 1 of the distribution of the relative differences can be found in Table 2. In addition, we found a slight dependence of the scatter on the 500 , peak and on the inferred targets values (except for gas ). The scatter seems to decrease as these values increase. From a statistical point of view, the scatter of baryon properties for high mass halos is smaller compared to low mass halos. A possible physical explanation is that massive clusters are more self-similar than smaller groups that present a larger halo-to-halo variation due to the stronger impact of non-gravitational processes. The relative difference for the NGBoost mean predictions is similar to those of the other tested ML models (see Fig. 6). However, overall scatter is higher, as shown in Table 2.

XGBoost vs NGBoost
As stated in section § 3, NGBoost is a probabilistic model that can learn to infer not only the mean predictions, but also the scatter. In Table 2, we show that the best model is XGBoost only when taking into account its mean predicted values. However, when it comes at mimicking the complete behaviour of the data, the scatter of the predicted targets tends to be underestimated for deterministic ML models.
In order to test whether XGBoost and NGBoost predictions are similarly spread as in true The300 data, we bin our baryonic targets in mass ( 500 ) bins. This is done for The300 true targets, XGBoost and NGBoost predictions. Then, for every mass bin, the scatter (standard deviation) of the baryonic properties is computed. Note that this process is repeated ten times for our ten disjoint K-folds. The results are displayed in Fig. 7 where we show the standard deviation per mass bin of XGBoost and NGBoost predictions divided by the standard deviation of the true baryonic data. As a general result, NGBoost successfully mimics the true scatter for most targets. In contrast, XGBoost baryonic properties are ∼ 0.1dex less spread. Notwithstanding, for gas and SZ the scatter of XGBoost is closer to the The300 true data with pred / true = 0.95 on average. Also note that the scatter for star is not completely well predicted by either of the models. However, NGBoost prediction of the scatter is also more precise in this target. The reader should bear in mind that although NGBoost successfully mimics the behaviour of the data in terms of predicted scatter, the predicted baryonic values of the XGBoost models are always closer to The300 true data.   Figure 6. The relative difference (y-axis) defined in Eq.(12) as a function of the XGBoost-predicted target variable (first column), the cluster 3D dynamical mass 500 (2) (second column) and the peak of the maximum value of the radial circular velocity profile across the halo's mass accretion history peak (7) (third column). From top to bottom, each row corresponds to different baryon targets: gas , star , gas , X , SZ . Dashed black lines are the mean value of the relative difference, which are very close to diff=0. Red regions and blue regions represent the 1 and 2 scatter regions, respectively. Additionally, 0% and 20% relative difference lines are represented in dashed green colour. The data is binned using sliding windows that contain roughly the same number of clusters. In addition, we have also computed the covariance between different baryonic properties and checked that XGBoost and NG-Boost predictions have a similar covariance structure than that of The300 simulation. The interested reader can find these results in the Appendix B.

ML inference of Baryonic properties in Dark matter only datasets
We now proceed to apply the trained ML model to infer the different baryonic properties in the full set of MDPL2 halo catalogues. We will use the 10 different XGBoost models trained on the reduced set of features of The300 clusters. In order to create the catalogue, we first build a dataset with the reduced set of features (shown in will be taken into consideration. Where training is a feature corresponding to the training dataset and MLDP2 is the same feature for the full MDPL2 simulation. Only 397 clusters out of 1,306,185 are outside the hyper-cube defined by the most important features and therefore, they are not considered for the analysis. In order to evaluate if the generated catalogue presents properties that are coherent with the properties of fully simulated data, we will compare our baryon properties with the halo mass for The300 and the full MDPL2 generated catalogue. These results are shown in Fig. 8 for different redshift values (columns). In these plots, the values of the targets (rows) are plotted with respect to 500 . For the targets gas and star , the plots show the relative fractions: where can be either gas or star. Error bars represent the intrinsic scatter of our baryon ML estimates on the MDPL2 catalogue (1 ), orange/brown regions correspond to 1 region for The300 test set predictions (for all the k-folds) and blue regions are the equivalent but for The300 true targets. Moreover, in the last row of the figure the number of clusters per bin is represented as a function of 500 for both The300 and MDPL2 datasets. As a general result, the XGBoost-predicted values for MDPL2 objects (black error bars) are similar and also their distributions per mass bin are comparable with the true values (blue region), i.e. in agreement with Fig. 6. However, the scatter of the predictions is slightly smaller (around 10-20%) than the corresponding scatter using the true values of The300 data for gas and star . This issue can be solved by using probabilistic regression models such us NGBoost where the scatter is predicted more accurately. The mean predictions (green squares) and scatter (green dashed lines) for NGBoost is also shown in the same figure. Furthermore, a similar result to the ones shown in Fig. 8 are obtained when plotting as a function of peak instead of 500 . We need to point out that for massive clusters ( > 8 × 10 14 ℎ −1 M ), the number of objects is similar in the The300 and MDPL2 simulations. Particularly, the last two mass bins are mostly composed of the same objects and the difference lies in the baryon properties of the The300 simulation.

Dependence of ML model predictions on DM mass resolution
The ML models have been trained on a particular DM simulation with a fixed resolution in mass. Here we are interested to compare the predictions of the ML model when applied to halo catalogues from simulations with lower mass resolution. Since some of the features of the halos are expected to be affected by resolution, then the infer baryon quantities from the ML models could also be affect by that. Since our goal is to make our ML models as universal as possible so they can be applied to different DM-only simulations with larger volumes, it is important to test for these effects. In order to do that, we are going to apply the trained XGBoost and NGBoost models in two simulations run with identical initial conditions but with a difference of a factor 8 in particle mass. For this test, we are going to use also another completely different realisation than MDPL2, i.e. the UNIT project. The UNIT 8 N-body cosmological simulations (UNITSIM, Chuang et al., 2019) are designed to provide accurate predictions of the clustering properties of dark matter halos using the suppressed variance method proposed by Angulo & Pontzen (2016). We particularly focus on one of the UNIT simulations with the same box side length than MDPL2, (i.e. 1 ℎ −1 Gpc ) and similar number of particles (4096 3 ). Furthermore, this simulation has also been performed with 8 times less number of particles (2048 3 ). For simplicity we will refer to these two simulations as UNITSIM4096 and UNITSIM2048 for the high and low resolution versions respectively. Dark matter cluster-size halo catalogues from R + C T are then selected for UNITSIM4096 and UNIT-SIM2048 following the same procedure described in § 2. We then apply the trained XGBoost and NGBoost models to these catalogues to infer the target baryon properties for each DM halos in the two versions. These baryon properties present similar statistics (mean and scatter per mass bin) as those shown in Fig. 8. In order to make a more quantitative comparison of the results for the two UNIT simulations, we bin the data as in Fig. 8 according to 500 and compute the difference of the mean values and estimate an upper limit for its scatter as a bin. The particular values of¯and¯are shown for 3 different snapshots in Fig. 9. As can be seen in this figure,¯ 0 with a small value for the scatter for all mass bins. The scatter¯is within ∼ 2% for gas and ∼ 0.5% for star . For gas , the residuals amount to ∼ 0.1 dex and for X and SZ up to ∼ 0.2 dex. Therefore, we conclude that the baryonic properties predicted by the ML model for the same halos simulated with a factor of 8 difference in mass resolution are statistically equivalent. NGBoost predictions have slightly larger scatter overall (green shaded area) and the mean values (green squares) are similar to XGBoost predictions.

VALIDATION OF THE GAS SCALING RELATIONS
Scaling relations are generally power laws that relate properties in astrophysical systems, such us the Colour-Magnitude Relation or the Tully Fisher relation (Tully & Fisher, 1977) for galaxies. The applications of the scaling relations are manifold, such as inferring masses of galaxy clusters that are sensitive to cosmological parameters (e.g. Planck Collaboration et al., 2016). For a recent review of scaling relations for galaxy clusters we refer the reader to e.g. Lovisari & Maughan (2022). The temperature-mass relation can be written as where ( Here, and ( = , X, SZ) are the parameters that we are interested in obtaining by fitting the above equations to our data. Once we have generated baryon catalogues for different N-body simulations we apply a simple linear fitting function in logarithm space to fit the data to the equations listed above. However, selecting data from different snapshots gives us small variations of the and best fitting parameters with redshift. Therefore, we use the following parametrization to study the redshift dependence: where ,0 and ,0 are the values of the intercept and slope at z=0 and and describe their possible dependence with respect to the redshift. With this new parametrization we apply a non-linear least square fitting model to fit the function described by equations (16), (17) and (18) updated with equations (19) and (20). The best-fitting parameters are shown in Tables 3, 4, 5 and 6. Note that we have used the mass corresponding to the N-body simulation (the feature variable 500 (2) as the mass of the cluster). Moreover, in order to study the dependence of these parameters on the mass, we opt for a similar approach as Le Brun et al. (2017), i.e. a fixed broken power law. We split the data in different mass ranges and fit the above equations. The considered mass ranges are: 1) log  14 because all of the scaling relations appear to break at that point for the radiative models. Furthermore, in order to study how the data is spread around the scaling law, we also used the MSE defined in Eq.(6) and the relative error defined in Eq.(12). Note that here the mean predictions < pred | > will be given by the corresponding scaling law.
As a general result, the fitting parameters are in agreement among the three different N-body simulations and are slightly different from The300 hydrodynamical simulation. This deviation, though small, is caused by the fact of considering the full box of dark matter only simulations instead of the smaller volume of the 'zoom' simulation. The effect of resolution is negligible for galaxy clusters. There is also a small difference between The300 simulations true data (The300), and the fitting counterpart using the ML predicted data (The300*). This slight difference can be mainly appreciated in the intrinsic scatter of the linear fitting function, which is generally smaller in the case of The300*. It is important to note that the scatter of the scaling law for The300 simulations is generally larger when comparing it with the values shown in Table 2, where the scatter (standard deviation of the relative difference) is reduced by a factor of 0.5 for the gas temperature, 0.3 for X and 0.45 for SZ . Moreover, the most relevant variables for each gas properties presented in table Table 1 can be used for finding analytical expressions for scaling laws with a reduced MSE using genetic algorithms (Wadekar et al., 2022). The difference between NGBoost and XG-Boost is that the MSE with respect to the scaling law is generally bigger in the case of NGBoost-predicted baryonic properties. This means that probabilistic regression mimics the scatter of the true data more accurately than non-probabilistic models.
As far as the redshift dependence is concerned, it is negligible for X and SZ where the parameters and are of order 10 −3 . However, the parameter 0.3 cannot be ignored. This indicates that the evolution of gas is relevant as it can also be appreciated in

SUMMARY AND CONCLUSIONS
Numerical simulations are key to studying galaxy clusters. On the one hand, with the current technology it is possible to perform large volume N-body simulations that can be useful to describe the dark-matter component. However, big volume hydrodynamical simulations cannot be carried out due to their computational demands. We have therefore trained a set of machine learning models to populate high volume dark-matter-only simulations with baryonic properties. In particular, we have defined our feature space as the R variables of DM-only halos and our target variables are directly estimated from T T H hydrodynamical simulations: the mass of the gas gas , the mass of the stars star , the gas temperature gas , X-ray Y-parameter X and the integrated Compton-y parameter SZ . All these quantities are integrated quantities in spherical region of overdensity 500 times the critical density at their corresponding redshift.
Particularly, we have considered four different ML models, random forest (RF), extreme gradient boosting (XGBoost), MultiLayer Perceptron (MLP) and Natural Gradient Boosting for Probabilistic Prediction (NGBoost). We have determined that XGBoost is the algorithm that is more suitable to our dataset and whose predictions are closer to the true hydrodynamical targets, as shown in Table 2. However, as depicted in Fig. 7, probabilistic regression is needed for a proper modelling of the scatter and therefore, NGBoost is the best model in that regard. We have applied an algorithm -Greedy Search Feature Importance Algorithm (GSFIA)to identify the features that have more predictive information. By using GSFIA, we have managed to reduce the dimensionality of our feature space from 27 to approximately 5 variables depending on the target variable. We have demonstrated that masses and velocities have a higher amount of predictive information while time evolution variables play a secondary role in the prediction of our targets. What is more, ellipticity, dynamical state, and spin features are redundant. A possible explanation for this is that our baryon targets are integrated in spherical regions.
Then, we have applied our trained ML models to populate halo catalogues with baryonic properties from two full box N-body simulations: the M D simulation (MDPL2) and the UNIT Nbody cosmological simulations (UNITSIM). The MDPL2 predicted baryon properties are compatible to those of The300 simulations, as shown in Fig. 8. The application on two UNITSIM simulations with 1 ℎ −1 Gpc box size 2048 3 and 4046 3 particles has determined that our model can be successfully applied to boxes whose resolution is up to 1/8 of the corresponding simulation used for training. This suggests that this is a promising method to populate the UNITSIM large volume N-body halos with baryon properties up to 27 Gpc 3 (i.e a 3 ℎ −1 Gpc box size with 6144 3 particles). This will be an excellent tool to study the large scale distribution of galaxy clusters in an unprecedented way. For instance, we can estimate the cosmic variance in the number counts of X-ray detected clusters from the eROSITA satellite all-sky survey Liu et al. (2021) by extracting many different light-cones from this large computational volume. This will be the subject of a forthcoming paper.
Furthermore, the scaling relations are powerful massobservable proxies. We have check that the best-fitting parameters inferred using our three mock DM full-box baryon catalogues are compatible. They nevertheless differ slightly from those of the The300, partially because of the considerable smaller number of cluster objects in the hydrodynamical simulations used to get the best fit values. This would suggest that mass completeness have an small impact, thought not negligible, in the calibration of the mass-proxies.
Our ML models have been trained using the results from a simulation with fixed cosmological parameters and a particular selection of values for the parameters of the subgrid physics models implemented in GADGETX (Cui et al., 2018). This is clearly a limitation of the applicability of the ML models presented in this paper for simulations with different cosmologies. This limitation could be alleviated by running an ensemble of thousands of simulations varying cosmological and astrophysical parameters and training the ML models so they can marginalise over all the parameters (either explicitly as feature variables or implicitly inside other features such as mass, etc). A recent example of this method is the CAMELS project (Villaescusa- Navarro et al., 2022). Unfortunately, given the low density of cluster sized objects, the simulations needed to repeat the CAMELS technique for galaxy clusters would require of Gigaparsec volumes with multi-billion dark matter particles and then, using the zooming technique, generating many hydrodynamical simulations of the same region with different baryon physics models. This is certainly a way to go for the future, but it is well beyond the scope of this paper.
Moreover, our test concerning the dependence of ML models predictions on DM mass resolution in subsubsection 4.3.1 suggests that the joint distribution of halo properties ( , peak , etc) is not very sensitive of the numerical resolution.
To conclude, our work shows that ML models are very useful methods for finding a mapping between dark matter halo properties found in N-body and the complex hydrodynamical simulations. We checked that, on average, the generated catalogue for the 3 darkmatter-only simulations used throughout this paper have the same distributions to that of true training set and therefore, they can be used for painting dark matter halos with baryonic properties that are directly related with observed quantities, providing added value to large volume collisionless N-body simulations.

DATA AVAILABILITY
The trained models and data products for MDPL2, UNITSIM2048 and UNITSIM4096 are publicly available at https://github. com/The300th/DarkML. Table 5. NGBoost best fit parameters: The best fit parameters for the − , X − and SZ − relations for the different simulation sets. The MSE in Eq.(6) and average scatter of the relative difference in Eq.(12), (in parenthesis), are also shown. For The300 simulation, the true values of the baryon properties have been used while for The300* the predicted NGBoost values are used instead. The relative error in the estimated parameters , and is always ≤ 10 −3 .   Figure A1. The covariance matrices, Q, of predicted target values for XGBoost and NGBoost, normalised to the covariance of the ground truth data from The300 simulations. Table A1. The feature variables used in this text from the R catalogue. The first column represents the variable name and their enumeration in brackets.

Variable Units Description
M2500c (0) ℎ −1 M Mass inside a radius of a sphere where the matter density is 2500 times the critical density at the cluster's redshift num_prog (1) total number of progenitors of the cluster M500c (2) ℎ −1 M Mass inside a radius of a sphere where the matter density is 500 times the critical density at the cluster's redshift M200c (3) ℎ −1 M Mass inside a radius of a sphere where the matter density is 200 times the critical density at the cluster's redshift Mpeak (4) ℎ −1 M The peak value of the halo mass across its accretion history mvir (5) ℎ −1 M halo mass within the virial radius Macc (6) ℎ −1 M halo mass at accretion time. Vpeak (7) km/s Peak value of Vmax(9) across mass accretion history. Vmax\@Mpeak (8) km/s Vmax at the expansion time at which Mpeak was reach Vmax (9) km/s maximum value of the circular velocity. Vacc (10) km/s Vmax at accretion time rvir (11) ℎ −1 kpc halo radius at virial overdensity vrms (12) km/s root mean squared velocity dispersion b_to_a(500c) (13) ration between the second largest shape ellipsoid axis and largest shape ellipsoid axis, for particles within 500 c_to_a(500c) (14) ration between the third largest shape ellipsoid axis and largest shape ellipsoid axis, for particles within 500 b_to_a (15) ration between the second largest shape ellipsoid axis and largest shape ellipsoid axis determined by method in Allgood et al. (2006) c_to_a (16) ration between the third largest shape ellipsoid axis and largest shape ellipsoid axis determined by method in Allgood et al. (2006) rs (17) ℎ −1 kpc comoving scale radius from the fit to a NFW (Navarro et al., 1997) density profile Rs_Klypin (18) ℎ −1 kpc comoving scale radius determined using Vmax and Mvir formula (Klypin et al., 2011) T/|U| (19) the ratio between the total kinetic and potential energies of particles within virial radius. Xoff (20) ℎ −1 kpc Offset between comoving density peak and the particles center of mass position Voff (21) km/s Offset between halo core velocity and the center of mass velocity for particles within the virial radius Spin (22) Peebles's dimensionless Spin parameter of the halo (Peebles, 1969). Spin_Bullock (23) Bullock's dimensionless spin parameter (Bullock et al., 2001) a (24) Expansion scale factor of the corresponding simulation snapshot scale_of_last_MM (25) Expansion scale factor of the last major merger with a mass ratio greater than 0.3 Halfmass_Scale (26) Expansion scale factor when the most massive halo progenitor reached 0.5 × Mpeak (4) properties is similar to that of The300 simulation. In order to do that, we have computed the covariance matrix as Cov( , ) = [( − ( ))( − ( ))] (B1) and the element-wise quotient between the covariance matrix of ML predicted targets and the true targets, i.e. :
In this way, Q shows how similar the covariance matrix of the predicted targets is with respect to the corresponding one of true targets. Q is plotted in Fig. A1 for both models XGBoost and NG-Boost. As can bee seen in the figure, Q values are very close to 1 but always less than 1. This means that the predictions of both ML models are distributed around the mean true values, but they do not completely reproduce the tails of the distributions of the real data. As a general result, NGBoost baryonic properties show a covariance structure closer to the corresponding ground-truth values of The300 simulations.
This paper has been typeset from a T E X/L A T E X file prepared by the author.