Context. Computational astronomy has reached the stage where running a gravitational Nbody simulation of a stellar system such as a Milky Way star cluster is computationally feasible but a major limiting factor that remains is the ability to set up physically realistic initial conditions. Aims. We aim to obtain realistic initial conditions for Nbody simulations by taking advantage of machine learning with emphasis on reproducing smallscale interstellar distance distributions. Methods. The computational bottleneck for obtaining such distance distributions is the hydrodynamics of star formation which ultimately determine the features of the stars including positions velocities and masses. To mitigate this issue we introduce a new method for sampling physically realistic initial conditions from a limited set of simulations using Gaussian processes. Results. We evaluated the resulting sets of initial conditions based on whether they meet tests for physical realism. We find that direct sampling based on the learned distribution of the star features fails to reproduce binary systems. Consequently we show that physicsinformed sampling algorithms solve this issue as they are capable of generating realisations closer to reality.
['gravitation', 'hydrodynamics', 'methods: numerical', 'open clusters and associations: general', 'Astrophysics - Astrophysics of Galaxies']
A machine learning framework to generate star cluster realisations
{'A machine learning framework to generate star cluster realisations': "George P. Prodan , ⋆ 1 , 8 , Mario Pasquato 2 , 3 , Giuliano Iorio 1 , 4 , 5 , Alessandro Ballone 1 , 4 , 5 , Stefano Torniamenti 1 , 4 , 6 , Ugo Niccolò Di Carlo 7 , and Michela Mapelli 1 , 4 , 6  \n- 1 Dipartimento di Fisica e Astronomia, Università di Padova, Vicolo dell'Osservatorio 3, 35122, Padova, Italy\n- 2 IASF Milano, via Alfonso Corti 12, Milano, Italy\n- 3 Ciela - Montreal Institute for Astrophysical Data Analysis and Machine Learning, Montréal, Canada\n- 4 INFN-Padova, Via Marzolo 8, 35131, Padova, Italy\n- 5 INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, Padova, Italy\n- 6 Institut für Theoretische Astrophysik, ZAH, Universität Heidelberg, Albert-Ueberle-Str. 2, 69120, Heidelberg\n- 7 SISSA - Scuola Internazionale Superiore di Studi Avanzati, via Bonomea 365, I-34136 Trieste, Italy\n- 8 Faculty of Sciences, University of Craiova, A.I. Cuza 13, 200585 Craiova, Romania  \nReceived June 05, 2024 ; accepted September 04, 2024", 'ABSTRACT': 'Context. Computational astronomy has reached the stage where running a gravitational N -body simulation of a stellar system, such as a Milky Way star cluster, is computationally feasible, but a major limiting factor that remains is the ability to set up physically realistic initial conditions.  \nAims. We aim to obtain realistic initial conditions for N -body simulations by taking advantage of machine learning, with emphasis on reproducing small-scale interstellar distance distributions.  \nMethods. The computational bottleneck for obtaining such distance distributions is the hydrodynamics of star formation, which ultimately determine the features of the stars, including positions, velocities, and masses. To mitigate this issue, we introduce a new method for sampling physically realistic initial conditions from a limited set of simulations using Gaussian processes.  \nResults. We evaluated the resulting sets of initial conditions based on whether they meet tests for physical realism. We find that direct sampling based on the learned distribution of the star features fails to reproduce binary systems. Consequently, we show that physics-informed sampling algorithms solve this issue, as they are capable of generating realisations closer to reality.  \nKey words. gravitation - hydrodynamics- (Galaxy:) open clusters and associations: general- methods: numerical', '1. Introduction': "The target of humanity's first deliberate attempt at interstellar radio communication was the star cluster M13 (Sta ff at the National Astronomy & Ionosphere Center 1975), and this act is a testament to the importance of star clusters to astronomy, which can hardly be overstated. For instance, star clusters are the likely birthplace of most stars (Lada & Lada 2003).  \nAstronomers use gravitational N -body simulations with the purpose of studying the evolution of star clusters (Aarseth 2003; Spurzem & Kamlah 2023). The goals, among others, are to predict gravitational wave emission from black holes and neutron stars formed by cluster stars (Rastello et al. 2021; Dall'Amico et al. 2021; Iorio et al. 2023), constrain the place of origin of the Solar System (Pichardo et al. 2012; Pfalzner et al. 2015), determine the stability and habitability of exoplanets (Spurzem et al. 2009; Malmberg et al. 2011; Parker & Quanz 2012), and understand how star clusters contribute to the overall evolution of the Milky Way as part of the field called Galactic archaeology (Chung et al. 2019; Pang et al. 2021, 2022).  \nStar clusters are understood to form within molecular clouds. Molecular clouds are high-density, low-temperature regions in the interstellar medium predominantly composed of molecular hydrogen. Star formation in molecular clouds occurs in small regions where the density becomes high enough (e.g. by tur-  \nbulent colliding flows) to start a fast runaway gravitational collapse (Evans 1999). Gas in molecular clouds is often animated by turbulent motions, which are ultimately responsible for the irregular and complex spatial distribution of young cluster stars (Krumholz et al. 2019; Krause et al. 2020; Ballone et al. 2020). The scientific usefulness of gravitational N -body simulations depends on the realism of initial conditions (primordial positions, velocities, and masses of stars), which hinges on the ability to properly model star formation in molecular clouds. Realistic initial conditions may be obtained by running hydrodynamical simulations. These are computationally expensive, with one simulation taking around 10 5 CPU hours (Ballone et al. 2021).  \nMachine learning can be leveraged to obtain samples of valid initial conditions at a fraction of the computational price by learning a generative model from the outputs of a limited number of simulations. The first attempt in this direction was presented by Torniamenti et al. (2022). Their work consists of a bespoke manipulation of hierarchical clustering models that is simple and e ff ective but lacks theoretical guarantees. Here we propose a different approach based on Gaussian processes (GPs).", '2. Gaussian processes': 'Gaussian processes serve as probabilistic models capable of predicting function values based on noisy observations (Wang 2022). Especially valuable when data is limited, GPs present a  \nmore fitting alternative to deep neural networks given the constraints of our application (Gri ffi ths 2023). For star cluster realisations, the distribution of stars in a cluster can be represented as a function of key physical parameters, such as star masses, coordinates, and velocities. Leveraging GPs, we can learn the star distribution in clusters in the space defined by these parameters, enabling the generation of synthetic clusters.  \nThe essence of GP modelling lies in Bayesian inference, where model beliefs are updated upon receiving new observations. A GP model is characterised by a mean function, µ ( x ) and a kernel function, K ( x , x \' ), with the kernel function reflecting the similarity degree between input points and influencing predictions (Rasmussen 2004). In our case, the kernel function aims to represent the spatial relations among stars in a cluster by capturing the relationships among the stars set at star formation and arising through dynamics via, for example, the e ff ects of gravitational attraction between the individual stars.  \nThe kernel function hyperparameters of GP models include the signal variance, σ y ; the length scale, l ; and the noise level, σ n . The signal variance measures signal amplitudes, the length scale indicates covariance decay distance, and the noise caters to errors. All of these parameters are trainable. Usually, the radial basis function kernel (RBF) is employed to compute the amplitude of the covariance function,  \nK ( x , x \' ) = σ 2 y exp " -( x -x \' ) 2 2 l 2 # . (1)  \nThe hyperparameters of the mean function can be defined depending on the choice of the prior mean function. For instance, a constant mean function µ ( x ) = c leads to one additional trainable hyperparameter.  \nTraining aims to minimize the log marginal likelihood (Rasmussen 2004) with respect to the above-mentioned hyperparatemers θ = { σ y , σ n , l , ... } , namely  \nlog P ( y | X , θ ) = -1 2 ( y -µ ) T ( K + σ 2 n I ) -1 ( y -µ ) -1 2 log | K + σ 2 n I | -n 2 log 2 π. (2)  \nThe first term ensures data fitting, the second one is responsible for regularisation, and the last one is a constant. Regularisation is an important part of model training that prevents overfitting (where the model learns irrelevant patterns in the data instead of the underlying relationships). By incorporating regularisation terms, the model generalizes better on unseen data, thus improving its predictive performance. Training with noise is an essential part of this work. By inserting noise, we obtained new cluster realisations that at the same time preserve some of the features of the training samples. Our aim for the GP fitting is to create a model that can reproduce these features.  \nWe drew the predictions from the posterior GP,  \ny ∼ GP ( µ, K + σ 2 n I ) , (3)  \nwhere µ , K are the posterior mean and kernel function and I is the identity matrix. One can draw two kinds of predictions. The first is the mean prediction, which generates the expected values that are similar to those determined for the training samples. This prediction is deterministic and represents the noiseless prediction of the GP. The second involves sampling from the GP\'s posterior distribution, incorporating noise. Here, the predictions are more diverse because they consider the inherent uncertainty  \nin the model. This can lead to a broader range of potential outcomes and possibly to the formation of di ff erent clusters. A complete mathematical overview on GP models is presented in Rasmussen & Williams (2006).', '3. Learning framework': "We have introduced a learning framework that utilizes GP modelling to learn the feature space distributions of stellar clusters, followed by sampling from these distributions to produce new clusters. This is an inverse problem, as the framework obtains new cluster realisations via simulation-based inference (Cranmer et al. 2020; Lueckmann et al. 2021). While training, the inputs of the GP model are the parameters of each star, θ i ( train ) with i = 1 , N for a cluster of N stars. The parameters are physical quantities of the stars, with each star having a specific mass, M ; location, r = ( x , y , z ); and velocity, v = ( vx , vy , vz ), with respect to the cluster's centre of mass and all of the parameters M , r , and v are not normalised. The GP learns a distribution over the probability density function of the feature space, f ( θ ).  \nOnce the GP model input features were determined, our next step involved computing the target distribution, y = f ( θ ). For this purpose, k -nearest neighbours density estimators (Mack & Rosenblatt 1979) were implemented to ascertain the probability density values, yi , which are normalised accordingly,  \ny ( train ) i = Nyi Σ N j = 1 yj . (4)  \nFurther, we implemented Markov chain Monte Carlo (MCMC) approaches to simulate new cluster realisations based on the Metropolis-Hastings algorithm (Chib & Greenberg 1995; Eckhardt 1987; Robert & Casella 2011). The states proposed in the Markov Chain include a set of N features at a given iteration step. We considered an iteration step n corresponding to a state Sn = { θ ( n ) 1 , θ ( n ) 2 , ..., θ ( n ) N } . The next candidate state, Sn + 1, was obtained by a jump with probability T ( Sn → Sn + 1). The acceptance probability of this jump is given by  \nα = min 1 , π ( Sn + 1) T ( Sn → Sn + 1) π ( Sn ) T ( Sn + 1 → Sn ) ! , (5)  \nwhere π is the stationary distribution of the chain. The probability of the new state can be written as the joint probability of the candidate features,  \nπ ( Sn ) = N Y i = 1 f ( θ i ) , (6)  \nwhere f is a probability density function drawn from the GP model that is already trained on the features θ ( train ) and the corresponding target values y ( train ) of the density function.  \nThis sequential method first crafts a GP statistical model based on the probability density function of the features space for the cluster stars. While sampling from a standard probability density function is generally straightforward, our specific problem requires meticulous adjustments due to the need for maintaining physical constraints and the presence of complex correlations. This ensures that the generated samples accurately reflect the underlying distribution and maintains the validity of the initial conditions. We employed traditional MCMC methods based on the Metropolis algorithm to directly sample in a sevendimensional space (DMCMC). Moreover, we propose a physicsinformed sampling approach based on learning the energy space distribution of nearest-neighbour star pairs (EMCMC).", '3.1. Direct sampling': "The direct approach is to define a seven-dimensional feature space that includes all the physical quantities of the stars in the clusters, meaning that for every star at the iteration of state n , we have θ ( n ) i = { M ( n ) i , x ( n ) i , y ( n ) i , z ( n ) i , v ( n ) xi , v ( n ) yi , v ( n ) zi } with i in 1 , N . Considering the large number of stars in the cluster, achieving convergence becomes challenging if all parameters are subject to change each time we propose a new candidate cluster. However, Eq. 5 remains valid also if we propose only a new candidate star θ ( n + 1) i such that the new state becomes Sn + 1 with  \nθ ( n + 1) j = (1 -δ i j ) θ ( n ) i + δ i j θ ( n + 1) i (7)  \nfor j in 1 , N , where δ i j is the Kronecker delta, and it is equal to one only if i = j ; otherwise it is zero. The acceptance rate will be  \nα = min         1 , f GLYPH<16> θ ( n + 1) i GLYPH<17> f GLYPH<16> θ ( n ) i GLYPH<17>         . (8)  \nTherefore, at each iteration a new subset of seven features θ ( n + 1) i is proposed with i chosen randomly. We drew new candidates using a normal distribution such that θ ( n + 1) ik ∼ N ( θ ( n ) ik , σ ) for every feature k of θ i , where σ is called the step size, and it is the amplitude of the perturbation applied on the previous state of the Markov chain and, mathematically, the standard deviation of the normal distribution centred on the previous sampled feature from which we sampled the new candidate.  \nThe new candidate is accepted or rejected based on the rate established by Eq. 8 and employing the same exact density function f drawn from the GP model, meaning that the GP noise seed must be kept unchanged during the sampling of one cluster realisation. The last condition is that it is essential not to break the detailed balance of the Metropolis-Hastings algorithm.  \nOne can notice that there is no conditioning on the global properties of the cluster and the sampler 'builds' the cluster by sampling stars individually. Therefore, this direct approach works similar to a black box that relies solely on the correlations captured by the GP model.", '3.2. Physics-informed sampling': "Sampling star systems, especially binary systems in close interaction, by generating new samples in a seven-dimensional space is challenging due to the low probability of proposing candidate stars in proximity within the physical space. This may result in erroneous modelling of star clusters at small scales, even if the sampling algorithm successfully replicates the density function given by the GP. Our solution to this problem is based on incorporating physical laws to the sampling process by focusing on binary stars creating chains of pairs based on the nearest neighbours and learning the distribution over the potential and kinetic energies of these pairs. Thus, the cluster is reduced to a chain of this kind, and from this chain we can reconstruct other clusters based on similar chains.  \nThis approach substitutes learning the full distribution of the mutual potential energies that results from N ( N -1) / 2 pairs of stars. This can be categorised as feature dimensionality reduction through feature selection (Jia et al. 2022) based on the heuristic that most of the binding energy is concentrated in the binary systems of the cluster (Torniamenti et al. 2021).  \nWe defined for each pair a feature set based on the mutual potential energy of the stars and their kinetic energies,  \nθ i = ( Ui , i + 1 , Ki , Ki + 1) with i in 1 , N -1. In this way, following the nearest neighbours chain, we get N -1 sets of features in the energy space. We considered that these sets of energy values fully define the energy state of the nearest neighbours chain. New energy states are sampled making use of a GP model trained on the probability density function of the energy space.  \nThe sampling scheme therefore starts by proposing new energy states. The first step of the solution is the EMCMC algorithm that samples the energy state of the nearest neighbours chain. The energy state at the n -th iteration, Sn = { θ ( n ) 1 , θ ( n ) 2 , ..., θ ( n ) N } , implies the mutual potentials of the nearest neighbours, U ( n ) 12 , U ( n ) 23 , ..., U ( n ) N -1 , N , and all the kinetic energies, K ( n ) 1 , K ( n ) 2 , ..., K ( n ) N . As before, we sampled for a random i a new candidate θ ( n + 1) i = ( U ( n + 1) i , i + 1 , K ( n + 1) i , K ( n + 1) i + 1 ) from a normal distribution based on the values of the previous iteration. This also implies a change for θ ( n ) i -1 and θ ( n ) i + 1 . In this case, applying Eq. 5 and Eq. 6 leads to the following acceptance rate:  \nα = min         1 , f GLYPH<16> θ ( n + 1) i -1 GLYPH<17> f GLYPH<16> θ ( n + 1) i GLYPH<17> f GLYPH<16> θ ( n + 1) i + 1 GLYPH<17> f GLYPH<16> θ ( n ) i -1 GLYPH<17> f GLYPH<16> θ ( n ) i GLYPH<17> f GLYPH<16> θ ( n ) i + 1 GLYPH<17>         , (9)  \nwhere f is the probability density function drawn from the GP model trained on the energy space. We note that at the boundaries, only two sets of features are subject to change.  \nUsing the new sampled chain, one needs to 'reconstruct' the cluster by defining its stars as entities with a position, velocity, and mass. The magnitudes of the relative distances between the stars, ri , i + 1 = ri + 1 -ri , and the velocities are derived according to the physical laws of Newtonian gravity,  \n| ri , i + 1 | = -GMiMi + 1 Ui , i + 1 (10)  \n| vi + 1 | = r 2 Ki + 1 Mi + 1 , (11)  \nwhere G is the universal gravitational constant and Mi and Mi + 1 are the masses.  \nFollowing the chain sequentially implies that the properties of the star i are known when deriving those of the next star in the chain (the star i + 1). The mass of the next star, Mi + 1, is sampled independently with the help of another GP model trained only on mass distributions. One can then determine the magnitudes of ri , i + 1 and vi + 1 . The next step is to explore the full seven-dimensional feature space of the stars to determine the optimal directions of ri , i + 1 and vi + 1 . For this, we used a GP model trained on the phase-space, GP ( r , v ). The exploration consists of proposing several Nc candidate directions for each vector. The candidates were chosen randomly from a uniform distribution. Using the GP model predictions, we chose the most probable position in the phase-space to determine the next star i + 1. The procedure was repeated until the end of the chain.  \nIn Algorithm 1, we provide the pseudocode of the proposed algorithm that works based on three GP models: GP ( θ ), GP ( M ), and GP ( r , v ). The function 'randomDirection' is employed to generate a random direction based on two random generated numbers, R θ and R ϕ , which are used to generate pairs of angles ( θ, ϕ ) that correspond to directions sampled uniformly on a three dimensional sphere. These directions are used to generate new candidates.  \nWe define two hyper-parameters, maxJump and Nc . The first parameter is used to control the maximum distance between the  \npairs of stars. In this way, the mass choice is conditioned on the spatial distribution of the cluster. Without this conditioning, the chain can be easily broken into multiple parts. This would lead to generating several smaller clusters. The other parameter, Nc , is the number of star candidates we generate randomly based on the values of | r i j | and | v j | . We use GP ( r , v ) afterwards to find the best candidate, that is, the star with the most probable position in the phase-space. As Nc is large enough, the resulting star features always point towards the densest region of the phase-space that is attainable at a given moment.  \nAlgorithm 1 EMCMC algorithm in pseudocode to find and select new star candidates for the generated cluster.  \nRequire: GP ( θ ), GP ( M ) , GP ( r , v ) , randomDirection , N , Nc , maxJump  \ncluster ← empty array to store the parameters of N stars cluster [0] ← initialize parameters of the first star  \nfor i from 0 to N-1  \nUi , i + 1 , Ki , Ki + 1 ← MCMC ( GP ( θ ))  \nMi , r i , v i ← cluster [ i ]  \n| r i , i + 1 | ← ∞  \nwhile | r i , i + 1 | > maxJump  \nMi + 1 ← MCMC ( GP ( M ))  \n| r i , i + 1 | ← Mimi + 1 Ui , i + 1  \nend while  \n| v i + 1 | ← q 2 Ki + 1 mi + 1  \ncandidates ← empty array to store Nc candidates  \nfor each candidate in candidates  \nλ r ← randomDirection()  \nλ v ← randomDirection()  \nr i + 1 = r i + | r i , i + 1 | λ r  \nv i + 1 = | v i + 1 | λ v  \ncandidate ← Mi + 1 , r i + 1 , v i + 1  \nend for  \nprobs ←GP ( r , v )( candidates )  \nbestCandidateIdx ← argmax ( probs )  \ncluster [ i + 1] ← candidates [ bestCandidateIdx ]  \nend for", '4.1. Dataset': 'We utilized a dataset consisting of ten clusters extracted from hydrodynamical simulations conducted by Ballone et al. (2020). These simulations were designed to accurately reproduce not only the clumpiness but also the fractal nature observed in starforming regions. Each cluster in the dataset has approximately 2500 to 4000 stars, with a cumulative mass ranging between 4000 and 40000 M ⊙ (solar masses). For our analysis, we divided this dataset into a training subset, consisting of seven clusters, and a validation subset made up of the remaining three clusters.', '4.2. Model training': "We used the gpytorch library (Gardner et al. 2018) to implement our models. This library o ff ers a modular and e ffi cient framework for GPs and supports graphics processing unit (GPU) acceleration. For all of our experiments, we employed one NVIDIA GeForce RTX 3060 GPU. Our GP model was trained using the RBF kernel, employing the exact marginal log likelihood as the loss function and using the Gaussian likelihood to  \nFig. 1. Training parameters (loss and noise level) with respect to the number of clusters used in the training. One validation cluster was used in the case of two or three training clusters; otherwise, we used two validation clusters. The training parameters correspond to the model's version at early stopping.  \n<!-- image -->  \ncalculate the posterior distribution. All models were initialised with a zero-prior mean. The necessary features for model training were derived from simulation data. Our training approach for each model integrates cross-validation (Hastie et al. 2001) and early stopping (Prechelt 1996), the latter having a patience threshold set at ten epochs. The properties of each trained model are summarised in Table 1.  \nThe models were trained on ten clusters, among which three are for validation. The training itself can be done using more or fewer clusters, depending on the context. We observed that the loss slightly increases when decreasing the number of clusters used in training (see Fig. 1), suggesting a slight decline of the model's performance. Also, we noticed a sudden increase of the loss when using eight clusters for training (and the other two for validation). This could be caused by underfitting, as suggested by the larger noise level. Therefore, we continued our experiments using a training validation ratio of seven to three.", '4.3. Sampling': "The sampling is based on a probability density function drawn from GP ( M , r , v ) for DMCMC or GP ( Uij , Ki , Kj ) for EMCMC. The initial state, in both cases, is initialised from a normal distribution N ( µ, σ ), with µ and σ computed over the feature distributions retrieved from the simulation data.  \nThe cluster size, that is, the number of stars in the cluster, N , is established when defining the initial state of MCMC proposals. The new candidate stars are proposed such that the state perturbations are uniformly distributed along the stars' features. This means that we do not define a maximum number of iterations, but several accepted states, NA . For instance, if NA = mN , new features for each star will be accepted exactly m times. In this way, one can recover m cluster realisations or use the first m -1 iterations as a burn-in strategy (Roy 2020).", '4.4. Generated clusters': 'The newly generated clusters were evaluated based on the distributions of the inter-particle distance, velocity, and mass. In addition, we computed global physical quantities such as the virial ratio, total mass, and the number of stars that are part of a bound system. The evaluation was based on comparing our generated clusters to the hydrodynamically simulated ones retrieved from Ballone et al. (2020). These clusters are our baseline, and our aim was to reproduce similar features. A two-dimensional spatial projection of such clusters obtained with DMCMC and EMCMC on the Oxy plane is shown in Fig. 2.  \nTable 1. Training results.  \nNotes. Columns: (1) Model; (2) best epoch; (3) training time per epoch; (4) scale length, l ; (5) noise level, σ n .', 'Hydrodynamical Simulations': "Fig. 2. Two-dimensional representation of one low-mass and one high-mass cluster generated with DMCMC and EMCMC, respectively. For comparison, the corresponding projections are illustrated also for the hydrodynamically simulated clusters m 1 . e 4 and m 9 . e 4 from the work of Ballone et al. (2020). The x and y coordinates are measured in parsecs. For EMCMC, we show the value of the hyper-parameter J (maxJump) and the total mass expressed in M 10 3 ⊙ = 1000 solar masses.  \n<!-- image -->  \nWe performed a set of experiments to sample a low-mass cluster (4000 M ⊙ ) of N = 3000 stars. The results are shown in Fig. 3. Comparing to the simulation data, for the inter-particle distance distribution of DMCMC, we noticed that there is a lack of stars at low inter-particle distances smaller than 10 -2 pc. This suggests a lack of close interactions between stars (i.e. binary systems). However, when analysing the total energy of the star pairs, we observed that there is a significant number of stars that are part of bound systems. For instance, the low-mass DMCMCgenerated cluster that corresponds to Figure 2 exhibits a number of 276 stars that are bound to other stars, whereas in the baseline training clusters there are between 400 and 800 binaries. According to the inter-particle distance distribution of the DMCMC clusters (Fig. 3), the orbital separations of this binary systems are large enough that gravitational encounters will likely unbind them during a cluster's evolution. This number is thus an upper limit to stable binaries, and the actual number of stable binaries could be significantly less, given the lack of close  \nbinaries. This is not the case for the EMCMC clusters, as there are also star pairs with low inter-particle distance that are able to keep the bound system intact.  \nOn the other hand, the physics-informed algorithm excels at sampling clusters with a balanced spectrum of inter-particle distances, which direct sampling in a seven-dimensional feature space fails to achieve. We also note that the energy spectrum of DMCMC exhibits a linear trend that di ff ers from those of the EMCMC-generated clusters or hydrodynamical simulations of clusters (see Fig. 5).  \nOther relevant physical quantities are presented in Fig. 4. In the left subfigure, it is shown how the number of stars varies with respect to the cluster's radius. We observed that there is not any specific trend in the distribution derived for the reference clusters. We noticed sudden jumps, meaning that there are regions of higher density, but there is one exception for one of the DMCMC clusters, where the number of stars increases smoothly. On the right side, we have plotted the histogram illustrating the  \n<!-- image -->  \n<!-- image -->  \nFig. 3. Distributions of inter-particle distances, velocity, and mass for low-mass clusters with 3000 stars generated using DMCMC and EMCMC methods. The mass distribution has a lower bound of 0 . 1. A simulated cluster from Ballone et al. (2020) is represented for comparison.  \n<!-- image -->  \n<!-- image -->  \nFig. 4. The distributions of key physical paramters for the new sampled clusters. The left subfigure shows the number of stars as a function of the distance to the cluster's centre of mass for a set of simulated clusters, DMCMC and EMCMC (three clusters each). On the right subfigure is plotted a histogram of the angular momentum distribution of the stars for the same sets of clusters as in the left subfigure. The angular momentum is measured in the units of the international system.  \n<!-- image -->  \ndistribution of angular momenta resulting after computing the magnitude of angular momentum of each star as  \n| L | = q ( ryvz -rzvy ) 2 + ( rzvx -rxvz ) 2 + ( rxvy -ryvx ) 2 . (12)  \nHere, we note that DMCMC clusters lack stars with high angular momentum. Probably, this happens because of the lower sampling rate of stars with high velocity magnitudes (see Fig. 3).  \nWe present a summary of additional experiments in Table 2. We show the global properties (virial ratio, estimated number of binary systems, and total mass) of the six realisations obtained via DMCMC and EMCMC. We performed the experiments in order to reproduce clusters with a low, intermediate, and high total mass.  \nWhen examining virial ratios, that is, the total kinetic energy divided by half of the total binding energy, we observed that the Metropolis algorithm delivers varied dynamical states. The EMCMC approach consistently produces realistic clusters; 60% of the clusters we sampled fall within a virial ratio between one and two, which is consistent with the training samples.  \nAll the experiments were aimed at generating clusters containing 3000 stars, which would allow us to provide computational performance comparisons. The number NA = 3000 was chosen arbitrarily; one can try di ff erent cluster sizes. However, the number of stars in the clusters should be the same order of magnitude as those used in the training (2000 - 5000 stars). More massive clusters could be sampled, but the framework performance has been tested primarily for lower-mass clusters, so  \napplying it to significantly larger clusters might necessitate retraining or fine-tuning the model to ensure that it handles the increased star count e ff ectively. Sampling with DMCMC takes around 15 minutes for NA = 3000 accepted samples, while EMCMC samples 3000 new accepted states during the same time, but it needs additional time, around 10 minutes, to reconstruct the cluster from the chain. The running time will scale linearly with the number of stars, as the sampling is done sequentially. One would expect EMCMC sampling to take longer, as we are estimating the probability density for three stars at each iteration. However, the acceptance rate is also three times higher for EMCMC. The experiments were done using an NVIDIA GeForce RTX 3060 GPU and a 12th Gen Intel Core i7-12700H CPU at 4.70 GHz. The training, sampling, and evaluation routines are available online on GitHub, 1 and there one can find several results that we have generated and the models we trained to obtain these results.", '5. Conclusions': "We have introduced a learning framework to generate new realisations of positions, velocities, and masses of star clusters to be used as initial conditions in N -body simulations. This framework allowed us to bypass the computational bottleneck represented by hydrodynamical simulations of star formation in molecular clouds, when provided an adequate training set. Based on well  \nTable 2. Monte Carlo sampling results using EMCMC and DMCMC algorithms.  \nNotes. The generated clusters contain 3000 stars. Columns: (1) #, cluster index; (2) method; (3) the lower bound for the mass distribution, Mlower , in solar masses ( M ⊙ ) ; (4) the step size in the Metropolis algorithm; (5) the acceptance rate; (6) the virial ratio; (7) the number of binaries; (8) the total mass of the sampled cluster.  \nFig. 5. Histogram with the energies of each pair of stars in the cluster. The energy for each pair was calculated as the sum of kinetic energies of the stars and their mutual potential energy. Negative energies correspond to gravitationally bound pairs. The resulting energy spectra for the clusters generated with EMCMC and DMCMC are compared to the clusters obtained via hydrodynamical simulations (Ballone et al. 2020).  \n<!-- image -->  \nunderstood GPs, our method yields predictable results and can be readily integrated into a simulation pipeline.  \nThe learning framework consists of two steps. First, a statistical model of the data is constructed through the GP, and second, new clusters are generated by sampling from the target distribution. The probability density function of the star features was employed to train the GP model. We tackled the problem in two ways: by sampling directly in the seven-dimensional feature space of the stars (DMCMC) and by sampling in a specific energy space (EMCMC). That was followed by reconstructing the cluster using the proposed algorithm.  \nIn the pursuit of sampling clusters with realistic spatial distributions and dynamical states, the EMCMC algorithm has proven to be particularly e ff ective, especially in addressing the crucial distribution of inter-particle distances, whereas DMCMC led to undesired distance distributions. The latter method and other conventional samplers work similar to a 'black box' in this case, relying only on the correlations learned by the GP model. The main di ff erence between the two approaches is that EMCMC incorporates physical principles. This physics-informed algorithm ensures that the generated samples are not only mathematically plausible but that they also possess physical significance.  \nThe DMCMC experiments showed that knowing the probability density function of the mass, position, and velocity is not su ffi cient to fully reproduce a cluster realisation. These properties are a consequence of physical interactions, and trying to reproduce them by sampling an estimated probability density function proved to be challenging due to a lack of binary stars. In this  \ncontext, the learned distribution can inform us about the possible location of a star in the feature space, but it cannot determine whether that star is part of a binary system, and it cannot reproduce its companion star. In contrast, with EMCMC, we did not change the sampling algorithm, but we sampled from a different target distribution that helps in answering both questions, namely, where the star is and whether it is a binary. The EMCMCgenerates clusters that not only match statistical properties but also reflect physical realities, such as binary star formation, making it a more robust and significant method for studies of stellar evolution and cluster dynamics.  \nAcknowledgements. This work acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017 (PI: Mapelli). MM and ST also acknowledge financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) STRUCTURES.", 'References': "Aarseth, S. J. 2003, Gravitational N-Body Simulations  \nBallone, A., Mapelli, M., Carlo, U. N. D., et al. 2020, Monthly Notices of the Royal Astronomical Society, 496, 49  \nBallone, A., Mapelli, M., Di Carlo, U. N., et al. 2020, MNRAS, 496, 49  \n- Ballone, A., Torniamenti, S., Mapelli, M., et al. 2021, MNRAS, 501, 2920\n- Chib, S. & Greenberg, E. 1995, The American Statistician, 49, 327\n- Chung, C., Pasquato, M., Lee, S.-Y., et al. 2019, ApJ, 883, L31\n- Cranmer, K., Brehmer, J., & Louppe, G. 2020, Proceedings of the National Academy of Sciences, 117, 30055  \nDall'Amico, M., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 508, 3045 Eckhardt, R. 1987, Los Alamos Science, 15, 131  \nEvans, Neal J., I. 1999, ARA&A, 37, 311  \nGardner, J. R., Pleiss, G., Bindel, D., Weinberger, K. Q., & Wilson, A. G. 2018,  \nGPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU  \nAcceleration  \nGri  \nffi  \nths, R.-R. 2023  \nHastie, T., Tibshirani, R., & Friedman, J. 2001, The Elements of Statistical Learning, Springer Series in Statistics (New York, NY, USA: Springer New York Inc.)  \nIorio, G., Mapelli, M., Costa, G., et al. 2023, MNRAS, 524, 426 Jia, W., Sun, M., Lian, J., & Hou, S. 2022, Complex & Intelligent Systems, 8 Krause, M. G. H., O ff ner, S. S. R., Charbonnel, C., et al. 2020, Space Science Reviews, 216  \nKrumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, Annual Review of Astronomy and Astrophysics, 57, 227  \nLada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57 Lueckmann, J.-M., Boelts, J., Greenberg, D. S., Gonçalves, P. J., & Macke, J. H.  \n2021, Benchmarking Simulation-Based Inference Mack, Y. P. & Rosenblatt, M. 1979, Journal of Multivariate Analysis, 9, 1 Malmberg, D., Davies, M. B., & Heggie, D. C. 2011, MNRAS, 411, 859 Pang, X., Tang, S.-Y., Li, Y., et al. 2022, ApJ, 931, 156 Pang, X., Yu, Z., Tang, S.-Y., et al. 2021, ApJ, 923, 20 Parker, R. J. & Quanz, S. P. 2012, MNRAS, 419, 2448 Pfalzner, S., Davies, M. B., Gounelle, M., et al. 2015, Phys. Scr, 90, 068001 Pichardo, B., Moreno, E., Allen, C., et al. 2012, AJ, 143, 73  \nPrechelt, L. 1996, in Neural Networks  \nRasmussen, C. E. 2004, Gaussian Processes in Machine Learning, ed. O. Bousquet, U. von Luxburg, & G. Rätsch (Berlin, Heidelberg: Springer Berlin Heidelberg), 63-71  \n- Rasmussen, C. E. & Williams, C. K. I. 2006, Gaussian processes for machine learning., Adaptive computation and machine learning (MIT Press), I-XVIII, 1-248  \nRastello, S., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 507, 3612 Robert, C. & Casella, G. 2011, Statistical Science, 26  \nRoy, V. 2020, Annual Review of Statistics and Its Application, 7, 387 Spurzem, R., Giersz, M., Heggie, D. C., & Lin, D. N. C. 2009, ApJ, 697, 458 Spurzem, R. & Kamlah, A. 2023, Living Reviews in Computational Astrophysics, 9, 3  \nSta ff at the National Astronomy & Ionosphere Center. 1975, Icarus, 26, 462 Torniamenti, S., Ballone, A., Mapelli, M., et al. 2021, MNRAS, 507, 2253 Torniamenti, S., Pasquato, M., Di Cintio, P., et al. 2022, MNRAS, 510, 2097 Wang, J. 2022, An Intuitive Tutorial to Gaussian Processes Regression"}