Hierarchical Bayesian Inference of Globular Cluster Properties

We present a hierarchical Bayesian inference approach to estimating the structural properties and the phase space center of a globular cluster (GC) given the spatial and kinematic information of its stars based on lowered isothermal cluster models. As a first step towards more realistic modelling of GCs, we built a differentiable, accurate emulator of the lowered isothermal distribution function using interpolation. The reliable gradient information provided by the emulator allows the use of Hamiltonian Monte Carlo methods to sample large Bayesian models with hundreds of parameters, thereby enabling inference on hierarchical models. We explore the use of hierarchical Bayesian modelling to address several issues encountered in observations of GC including an unknown GC center, incomplete data, and measurement errors. Our approach not only avoids the common technique of radial binning but also incorporates the aforementioned uncertainties in a robust and statistically consistent way. Through demonstrating the reliability of our hierarchical Bayesian model on simulations, our work lays out the foundation for more realistic and complex modelling of real GC data.


INTRODUCTION
Globular clusters (GCs) are massive collections of stars found in the outskirts of every type of galaxy.The dominant mechanisms that govern a cluster's evolution -stellar evolution, two-body relaxation, and the external tidal field of their host galaxy -lead to clusters becoming mass segregated over time as they evolve towards a state of partial energy equipartition while playing host to stellar collisions and mergers (Spitzer 1987;Heggie & Hut 2003;Trenti & van der Marel 2013).These internal and external processes can also lead to clusters developing radial anisotropy (Watkins et al. 2015;Tiongco et al. 2016;Jindal et al. 2019).To gain a deeper understanding of how these various dynamical processes shape cluster evolution, we need to accurately measure the current spatial and kinematic distribution of stars within a given cluster, which can then be used to constrain the conditions under which GCs form.Ultimately, knowing the formation histories of a population of GCs provides insight into the formation and evolution of their host galaxy (West et al. 2004).
A large number of distribution functions (DFs) have been proposed to model the observed distribution of stellar positions and velocities in GCs (Woolley 1954;Michie 1963;King 1966;Wilson 1975;Gieles & Zocchi 2015;Claydon et al. 2019).These DFs share the similarity that both density and velocity dispersion profiles decrease to zero out to a truncation radius.How the DF drops to zero varies from model to model, with additional treatments necessary to address radial anisotropy (Michie 1963), globular cluster rotation ★ E-mail: ywen@caltech.edu(Varri & Bertin 2012), or mass segregation (Da Costa & Freeman 1976;Sollima et al. 2015).
Historically, one can fit GCs with the aforementioned models by comparing the observed and theoretical surface brightness profiles or density profiles and sometimes kinematic profiles to determine the model parameters.Several different DF-based models have been fit to Galactic (McLaughlin & van der Marel 2005;Miocchi et al. 2013;de Boer et al. 2019;Cohen et al. 2021;Cheng & Jiang 2023) and extragalactic GCs (Woodley & Gómez 2010;Webb et al. 2013;Usher et al. 2013).Beyond determining structural parameters, these DF-based models have also been used to probe more complex features of GCs such as the existence of intermediate-mass black holes (IMBHs) (Zocchi et al. 2019) and their stellar initial mass functions (IMFs) (Dickson et al. 2023).In addition to the traditionally-used surface brightness, number density, and velocity dispersion profiles, the comparisons to DF-based models have increasingly considered other observables such as velocity anisotropy, stellar mass functions, and pulsar timings (Gieles et al. 2018;Hénault-Brunet et al. 2020).
The typical approach to fitting observational data with models is to first radially bin the observed stars, then estimate the uncertainties on the counts in each bin, and finally minimize the  2 between the observed profile and the corresponding theoretical model profile.However, binning data and minimizing  2 are undesirable for several reasons: (i) Information about each individual star is lost.(ii) Systematic errors are introduced due to the exact binning scheme used.
(iii) Error estimation for the counts of each bin can be complex, depending on the completeness of the dataset, contamination from non-cluster stars, and measurement errors (which lead to inter-bin covariances).
(iv) When trying to simultaneously fit multiple profiles, one must assume how to weight the importance of each fit.It must be decided whether the total  2 is simply the sum of the individual  2 values calculated for the density, kinematic and other profile fits or if they should be weighted differently, which is not a statistically robust procedure.
To avoid radial binning and to make the inference more statistically robust, Eadie et al. (2022) (hereafter EWR22) has demonstrated the use of Bayesian inference to estimate the model parameters, cumulative mass profile, and mean-square velocity profile of a GC based on the positions and velocities of individual stars and assuming the lowered isothermal DF proposed in Gieles & Zocchi (2015) (hereafter GZ15).Focusing on the ideal cases with simulated data, EWR22 assumed that the positions and velocities of individual stars are completely known without any measurement errors, and their work investigated the catastrophic effects of selection bias on parameter fits if selection bias is not properly modelled.
The purpose of this work is to follow up on EWR22 to address the issues of incomplete data and measurement errors through hierarchical Bayesain inference (HBI).We build our hierarchical Bayesian model and demonstrate the use of HBI on simulated GC data to show the reliability and robustness of HBI to recover a cluster's structural parameters in the presence of missing data and measurement uncertainties.Furthermore, we move away from the Cartesian coordinate used in EWR22 to the projected space on the plane of the sky (i.e., the reference frame in which actual data are measured), and we include the inference for the positions and velocities of the cluster center into our HBI model as well.The code to model GCs and perform inference with our hierarchical model can be found at https://github.com/y52wen/hbmlimepy.
The paper is organized as follows.In Sec. 2, we review the lowered isothermal DF in GZ15 and introduce the simulated dataset generated from their limepy python package that we use to validate our statistical method.We describe our hierarchical model in Sec.3.1 and specify our model priors in Sec.3.2, followed by a review of Hamiltonian Monte Carlo (HMC) in Sec.3.3 and a description of our approach to emulating the lowered isothermal DF through linear interpolation in Sec.3.4, both of which are necessary for performing inference on our hierarchical model.In Sec. 4, we report our inference results on different GC simulations to demonstrate the reliability and robustness of our approach to model GCs.We summarize our findings and discuss future works for modelling real GC data in Sec. 5.

The Lowered Isothermal DF
In this work, we work with simulated data generated from the limepy code of lowered isothermal models for GCs introduced in GZ15.The limepy code has been successfully used to fit GCs from both realworld observation and simulated data.For example, Zocchi et al. (2016), together with Baumgardt & Hilker (2018) and Hénault-Brunet et al. (2019), successfully demonstrated that direct N-body simulations of star clusters could be well-fit by the lowered isothermal DF, while de Boer et al. (2019) used this class of DF to determine the structure of galactic GCs using Gaia data (Gaia Collaboration 2018b).The parameters for the lowered isothermal DF are  lp = (, Φ 0 ,  tot ,  h ), and the physical descriptions of these parameters are summarized in Table .1.In general,  and Φ 0 impact the shape and concentration of the GC profile, while  tot and  h are scaling parameters that determine the mass and size of the GC.
In this paper, we restrict our attention to isotropic GCs (i.e., the anisotropic radius   → ∞ under the default in limepy) with only a single mass component (i.e.we assume stars in a GC have the same mass).In reality, GCs often develop some radial anisotropy due to various dynamical processes and GC stars have a wide range of masses, so it is better to treat a GC system as the combination of several single mass models (Da Costa & Freeman 1976;Hénault-Brunet et al. 2019).However, anisotropic, multi-mass models pose significant challenges for hierarchical Bayesian inference, which we discuss in Sec. 5, and so we leave its implementation to future work.
In the case of isotropic GCs, a value of  = 0 in the lowered isothermal model is equivalent to the Woolley (1954) model, while the  = 1 and the  = 2 cases reduce to the King (1966) model and the Wilson (1975) non-rotating model respectively.Therefore, the lowered isothermal DF can be considered as the generalization of all of the above familiar DFs for GCs through varying , the truncation parameter that decides how fast the DF falls off to zero at large radii.
For the single-mass, isotropic models used in this work, the DF is given in GZ15: The DF in Eq. 1 depends on the total energy  in Eq. 2, where  is the velocity and () is the potential at distance  from the centre.In Eq. 1, the energy  is lowered by the potential at the truncation radius  ( t ).The function   (, ) is defined as where (, ) ≡ (, )/Γ() is the regularized lower incomplete gamma function (see Appendix D of GZ15 for properties of   (, ) and (, )).The potential function () in Eq. 1 and 2 needs to be solved consistently through the following non-linear Poisson's equation 1 where DF  ( (, )) given in Eq. 1 also depends non-linearly on ().To solve Eq. 4 numerically, the equation is often nondimensionalized as in King (1966) and GZ15.
In the dimensionless case, the model is completely specified by two parameters: the central potential Φ 0 , which is a required boundary condition for solving Poisson's equation and defines how concentrated the model is; and the truncation parameter , which controls the sharpness of the truncation of the model.The physical units of the model are defined by two scales: the velocity scale  and the phasespace normalization constant , which in turn determine the total mass  total and the GC radius scale  h .These two scale parameters  and  in Eq. 1 need to be solved consistently with the Poisson's Equation of Eq. 4 to ensure proper normalization of the DF to 1, so  and  are in fact functions of  lp .The limepy code fully implements and solves the above Eq.1-4 and is able to give the DF  (, ) for any  lp .We consider all possible combinations of these four globular cluster parameters.
The parameter combination ( = 2.0, Φ 0 = 8.0) is excluded due to its closeness to the cutoff boundary (shown in Fig. 2) that distinguishes realistic GC models from unrealistic ones (see Sec. 3.2 for a more detailed explanation).This leaves us with 32 sets of different parameterss.We generate 10 simulations for each of the parameter sets and test our model on these simulations.

Description Possible Values
number of stars measured in cluster 100, 500, 1000  c distance of cluster center to earth (kpc) 1, 2, 5, 10  measurement uncertainties (mas, mas/year) 0.02, 0.1, 0.5 Table 2. Hyperparameter values used to simulate stars in GCs that reflect realistic observation conditions.We fix the structural GC parameters at (, Φ 0 ,  tot ,  h ) = (2.0,5.0, 10 5 , 3) while changing the above hyperparameters of the simulations individually.In the third row,  denotes the measurement errors for angular positions (mas), parallax (mas), and proper motions (mas/year) for an individual star, where we assume the errors for all five components are the same.We also assume that every star in GC shares the same measurement uncertainties.

Simulation Parameters
We develop and test our method for GC parameter inference with simulated spatial and kinematic data of GC stars in the Heliocentric equatorial coordinate system.The limepy code can directly sample stars from the lowered isothermal DF in Cartesian coordinates (x GC , v GC ) within a GC-centric reference frame, and we translate the stars' coordinates to the Heliocentric Cartesian coordinate (x, v) with the GC center parameters   (the default value given in Eq. 5).
We then convert the Heliocentric Cartesian coordinates (x, v) to the Heliocentric equatorial coordinates (, ), for which the transformations (, ) =  −1 (x, v) are explicitly given in Eq.A7-A12 of Appendix A.
To demonstrate the reliability and the robustness of our method, we test and validate our HBI model on different structural parameters of GCs covering a wide range of GC morphology.For the values given in Table 1; we consider all possible combinations of these four globular cluster parameters (excluding parameters with ( = 2.0, Φ 0 = 8.0) due to the cutoff boundary shown in Fig. 2), which result in 32 sets of parameters with different values.We generate 10 simulations, with each simulation being a different realization of the lowered isothermal DF, for every parameter set.We will then obtain constraints though our HBI model on these 320 simulations.By default, we draw  = 1000 stars from the lowered isothermal DF in each simulation.In the Gaia DR2 data, the number of stars with available proper motion in a Galactic GC can range from tens to thousands (Vasiliev 2019), with 1000 being a representative figure that balances the amount of information available and the complexity of the model to demonstrate the performance of our inference procedure.
The default position parameters for the centers of globular clusters are chosen to be  (Vasiliev & Baumgardt 2021).The choice of  c impacts the accuracy of the constraints, since GCs located further away have more uncertain distance estimates.We assume the cluster is observed at a distance of  c = 1 kpc, which corresponds to the parallax value   = 1 mas.In reality, GCs can have distances of a few kpcs or tens of kpcs away, and we investigate the impacts of changing  c values in Sec.4.2.
We also add measurement errors to the simulated data  = (, ) for each star.The position and proper motion measurements  obs for each star are assumed to follow a normal distribution with a certain standard deviation  (which we refer to as the measurement error).Based on Gaia Collaboration (2018a), we assume the measurement errors for the declination angle , the right ascension , and the parallax  for all stars are   =   =   = 0.1 mas, and the measurement errors for the proper motions are assumed to be    * =    = 0.1 mas/year.The above assumed astrometric measurement errors correspond to assuming our samples stars with G-magnitude at around 17 (Gaia Collaboration 2018a).By default, we assume that the radial velocities for all stars are unknown, which is a reasonable assumption for cluster stars measured through the Gaia data.
In each case, we assume that there are 1000 stars with measured properties in each cluster.However, we also explore the effects of different numbers of star as given in the first row of Table 2 (see also Sec. 4.2) on the inference results.We also explore the impact of different measurement uncertainties (as given in the third row of Table 2) in the angular positions and angular velocities on the parameter inference.When we run the 320 simulations with different lowered isothermal model structural parameters (Φ 0 , ,  tot ,  h ), we keep all the hyperparameters in Table 2 at their default values ( c ,   R , , ) = (1 kpc, 0, 0.1, 1000).When we explore the impacts of different hyperparameter values, we keep the structural parameters fixed at (Φ 0 , ,  tot ,  h ) = (5, 2, 10 5 , 3).

A Hierarchical Model
Using the simulated spatial and kinematic data of stars from each GC mentioned in Sec. 2, we take a hierarchical Bayesian approach to infer the model parameters.With the maturation of inference techniques for large, complex Bayesian models over the last two decades, hierarchical modelling has gained momentum and been applied to a wide range of scientific problems across many disciplines.In the context of astronomy, HBI has been used to estimate photometric redshift in galaxy surveys (Leistedt et al. 2016(Leistedt et al. , 2023)), and estimate the mass of the Milky Way in (Eadie et al. 2018;Eadie & Jurić 2019;Shen et al. 2022), to name a few examples.A graphical representation for our hierarchical model to determine GC properties is shown in Fig. 1.
From Bayes' theorem, the posterior probability of a vector of model parameters  given data , is For each star Graphical representation of the hierarchical Bayesian model built in this work to analyze simulated GC data with the lowered isothermal distribution function.All six phase-space measurements in the heliocentric sky coordinates are given a multilevel treatment; the observed variables are assumed to be drawn from Gaussian distributions centered around true (latent) phase-space parameters, with standard deviations equal to the measurement errors.For the positions and proper motions, the measurement uncertainties from Gaia are incorporated into the model.The heliocentric latent phasespace parameters are then transformed into the Cartesian coordinate and used to estimate the four structural parameters of the lowered isothermal distribution function and the six phase-space parameters that describe the GC center.
For the parallax and the radial velocity of the individual star, we sample in ) -that is, the difference between the star's radius (radial velocity) and those of the GC center.We use these differences as the sampling parameter to reduce the correlations between parameters and improve the performance of HMC.
where ( |) is the probability of the data conditional on the model parameters (also known as the likelihood L (; ), a function of model parameters for fixed data), () is the prior distribution on the model parameters, and ( ) is the "evidence" or prior predictive density.The latter is a constant, leaving us with a target distribution proportional to the posterior distribution ( | ), which we estimate through Markov Chain Monte Carlo (MCMC) sampling (see Sec.3.3).We devote the rest of this section to specifying our model, which is summarized with the directed acyclic graph (DAG) shown in Fig. 1.
Our simulated data products, described in Sec. 2, are the six observed Heliocentric equatorial phase-space components of each individual star ), which are shown in the base level of the graphical representation for our Bayesian model in Fig. 1.
Each component of a star's observed equatorial phase-space coordinates  obs1 is assumed to follow the normal distribution where  is the individual component of the true heliocentric equatorial phase-space coordinates (, )), and   is the observed measurement error for the particular phase-space component of the star.
We have assumed the measurements for each phase-space component are independent, which is often not the case for Gaia measurements since the measurements for parallax, angular position, and proper motions are often correlated with non-negligible covariance.Here we assumed independent Gaussian distributions for simplicity, and we can potentially change to multivariate Gaussian distributions to incorporate covariance.Any missing data can be incorporated as a variable with a flat (uniform) prior.Our model therefore maximizes the available phase-space information.We are able to make parameter estimates using any combination of sources, whether they are missing positions, velocities, neither, or even both.
The true heliocentric equatorial phase-space components (, ) of each star are then transformed into the Heliocentric Cartesian coordinates (x, v) =  (, ), where the coordinate transformation  is given in Eq.A1-A6 of Appendix A. We next transform the Heliocentric coordinates into the GC-centric coordinates for each star: where the center of GC in Heliocentric Cartesian coordinates is ).For Eq. 8, we assume there is no rotation present for stars around GC centers.In the future, we could include rotation for more complex GC modelling (see Sollima et al. (2019) for an example).
As seen in Fig. 1, we treat the the coordinates of the GC center  c as parameters and directly infer them from our hierarchical model.We therefore obtain estimates for GC position and velocity centers in a Bayesian framework, which is in contrast to the usual approach of setting the mean of GC stars' positions as the GC center (Gaia Collaboration 2018b; de Boer et al. 2019) and other more sophisticated approaches to determine centers such as the pre-slice method and the density contour method (Goldsbury et al. 2010).The angular positions of GCs can be treated as fixed; each star's angular position has been measured with exquisite accuracy through Gaia.However, for the distance and kinematic centers of GCs, our Bayesian method can robustly include measurement errors and missing data for stellar position and kinematic data, with the potential for providing better and more statistically consistent constraints on GC phase-space centers.Our inference results on the GC centers for our simulations are discussed in detail in Sec. 4.
Given the GC structural parameters  lp , the likelihood for observed GC stars having true positions and velocities (ì where the lowered isothermal DF  lp is specified in Eq. 1,  is the number of stars, and the DF is normalized by the total mass to become a probability distribution function.By taking the products of DF evaluated at each star, Eq. 9 assumes that the sampling of observed stars from all GC stars are uniform and each observed star is independently selected, which is an overtly idealistic assumption for real surveys.We will try to incorporate more realistic modelling for selection functions in future works.
To write the DF in Eq. 9 in terms of the true Heliocentric equatorial coordinates, we use the coordinate transformations in Eq. 8 and obtain where the second term is the Jacobian for the coordinate transformations from the Heliocentric equatorial coordinates (ì p, ì q) to the GC-centric Cartesian coordinates (ì x GC , ì v GC ).The translation by the GC center does not change the Jacobian.Ignoring the Jacobian terms in Eq. 10 will bias the inference results, especially for the distance of the cluster center  c .
For the Heliocentric equatorial coordinates of each star, the parallax (distance) is known with the least accuracy among the position coordinates, while the radial velocity can be unknown.As a result,   and  c (and  , and  , ) are highly correlated with each other, since a shift in the cluster center  c (or  , ) can be compensated by a reverse shift in   (or  , ) for all stars.To reduce the correlations among parameters for better inference performance, we define the new parameters Δ  =   −  c and Δ , =  , −  ,c .In these new definitions, Δ  and Δ , are centered at zero when no radial distances or radial velocities of the stars are known, and these variables for individual stars will no longer shift along with the values of the cluster centers.We will infer Δ  and Δ , instead of   and  , .For convenience, we write the coordinate parameters in which we perform inference as which are circled in blue in Fig. 1.The relationship between the sampling coordinates and the Heliocentric equatorial coordinates is where  = 1 kpc/1 mas is the conversion factor between distance and parallax.
For the error modelling of the individual star given in Eq. 7, since ( ì  obs , ì  obs ) only depends on the GC parameters ( lp ,  c ) through the sampling coordinates (ì , ì ).The lowered isothermal DF 2 Here the bold vector notation ì x GC is used to represent the Cartesian position coordinates of all stars in contrast to the coordinates of one star x GC ≡ x GC,i .written in terms of the Heliocentric equatorial coordinates is given in Eq. 10, which models the relationship between each individual star and the entire GC.Using Bayes' theorem, we can now write down the posterior of our whole model, which is illustrated in Fig. 1: where Eq. 14 is obtained using Eq. 13 and the definition of the conditional probability, and Eq. 15 is obtained by substitution with Eq. 10.
The Jacobian term in the second line of Eq. 15 is explicitly specified by Eq.A14 of Appendix A. Eq. 15 reflects the hierarchical structure of our model: the first line of Eq. 15 gives the individual-level (starlevel) modelling, while the second line specifies the population-level (GC-level) modelling, and they together constitute our model likelihood.The third line   lp ,  c is our model prior, which we will discuss in the following section.

Prior Distributions
Two advantages of Bayesian inference are the ability to incorporate meaningful prior information and the necessity to state this explicitly.Our goal is to choose priors that are broad and generally nonrestrictive so that we can apply the same priors to models with a wide parameter range, while excluding parameters regions that are physically impossible or irrelevant.We now clarify our choice of the model prior   lp ,  c , which has the following structure: where the priors for all the GC parameters are independent except for Φ 0 and .For the structural parameters  lp , we have some prior knowledge on the total mass  total and half-mass radius  h based on observations of Milky Way GCs, while having considerably less prior information on the values of  and Φ 0 , aside from the physically allowable, positive values.We choose the priors for  and Φ 0 as the following: Our prior for Φ 0 is the same as the one in EWR22, while we modify the prior on the truncation parameter  in EWR22 by imposing an upper boundary   (Φ 0 ) that depends on Φ 0 to exclude models that are inapplicable to GCs.As discussed in Gomez-Leyton & Velazquez (2014) and GZ15, the density profile for the lowered isothermal model is infinite when  ≥ 3.5.There is also a class of isotropic models that are finite in extent but showing a sharp upturn in the density profile at large radii, which makes these models inapplicable to star clusters.To separate these inapplicable models from those realistic models for GCs, we follow GZ15 and use limepy to compute the ratio of the dimensionless virial radius r = − M2 /(2 Û) over rℎ for a grid of models with 1 ≤ Φ 0 ≤ 20 and 0.001 ≤  ≤ 3.49 (Fig. 2).One can see that for any given Φ 0 value, once  increases beyond a certain threshold, the change in r / rℎ becomes abrupt, which signals the separation between two classes of models.We use r / rℎ > 0.64 as the criterion for identifying the class of models relevant to modelling GCs.The upper boundary   (Φ 0 ) function for  is chosen to be a 10 th -degree polynomial fitted to the boundary corresponding to r / rℎ ≈ 0.64 (the blue line in Fig. 2), shifted downward by 0.2.  (Φ 0 ) is plotted as the black dashed line in Fig. 2. The 0.2 downward shift is chosen to exclude parameter spaces that are too close to the r / rℎ ≈ 0.64 boundary.The black dashedline represents the   (Φ 0 ) that we used as the upper boundary in the uniform prior of  for our Bayesian model.
For the mass and radius of GCs, we sample the parameters and set the priors in log scale where log 10  total and log 10  h are assumed to follow normal distributions (in units of pc and  ⊙ respectively): where log 10  h is truncated by the lower bound  and the upper bound  to ensure the GC radius is both positive and not overly large ( h ⪅ 30 pc).Our prior for  total is the same as the one chosen in EWR22.In general, GC masses span about an order of magnitude and setting the log-normal prior on  total is supported by the near universal GC mass function (Harris 2010).For  h , EWR22 chose a more restrictive prior that also depends on the value of  h , while we choose a broader prior that reflects the distributions for the half-mass radii of 168 Milky Way GCs included in the Baumgardt's catalogue 3 , of which the structural parameters are determined in Baumgardt & Hilker (2018).In general, our inference results are not sensitive to the exact choice of priors for  total and  h .
For the GC center parameters  c , we assume the prior for each 3 https://people.smp.uq.edu.au/HolgerBaumgardt/globular/phase-space component is independent as seen in Eq. 18.For each component  c (where  = , , ,   * ,   ,  R ), we adopt an empirical Bayes prior and assume such that the component of the center is bounded by the minimum and maximum of the corresponding components of all observed stars in GCs.The above prior ensures the basic physical validity of the model while remaining generally uninformative.One can possibly put a tighter prior on the GC center compared to Eq. 23.
For the coordinate parameters of each individual star (ì , ì ) (Eq.11), which we also need to sample in MCMC, they are entirely determined by the conditional probability  ì , ì | lp ,  c , which is considered as part of the likelihood in Eq. 15.From a sampling perspective, (ì , ì ) can be considered as having a flat (non-informative) prior.For practical purposes, flat priors are hard to sample, so we recast the measurement error distribution ( obs |) in Eq. 7 (the individual-level modelling) as the prior for , ,   * ,   and remove these distributions from the likelihood part of Eq. 15.For Δ and Δ R , the priors are still chosen as flat and ( obs |) and ( obs R | R ) remain as part of the likelihood due to our change of sampling coordinates from  and  R to Δ and Δ R .In general, the measurement error modelling distributions  ì  obs , ì  obs | ì , ì  in Eq. 14 has the flexibility to be either treated as part of the likelihood or the prior for the true phase-space coordinates ( ì , ì ) in MCMC.

Inference with HMC
The standard approach of Bayesian inference is to use Markov Chain Monte Carlo (MCMC) algorithms, where the posterior ( | ) in Eq. 6 is directly sampled based on the likelihood and the priors without the need to perform integration to determine the evidence ( ).At each step, the MCMC algorithm makes proposals to new positions in parameter space and checks whether the posterior at the proposed position is favorable relative to the current position.
The likelihood and priors of our hierarchical model have been specified in Sec.3.1 and 3.2.In total, there are 6 + 10 parameters (Fig. 1 and Eq.15), where  is the number of observed stars in the GC.Traditional MCMC algorithms such as Metropolis-Hastings and Gibbs samplers (Metropolis et al. 2004;Hastings 1970;Casella & George 1992) are unable to perform inference for high-dimensional problems with hundreds or thousands of parameters, which is the case for our model.Other more modern MCMC samplers such as affine-invariant ensemble sampler and slice sampler also struggle in high-dimensional problems (Huĳser et al. 2015;Neal 2003).For hierarchical models with analytical probability distributions in the exponential family, it is often possible to marginalize over certain lowerlevel parameters (or nuisance parameters) to significantly reduce the dimensions of problems (Eadie et al. 2015;Gelman et al. 2004).However, for our model likelihood in Eq. 15, the DF is based on a numerical ordinary differential equation (ODE) solver implemented in GZ15 without analytical solutions.Furthermore, the coordinate transformations from ( ì , ì ) to (ì  GC , ì  GC ) are highly non-trivial, which prevents the marginalization over the phase space coordinates ( ì , ì ) of each individual star.Therefore, our inference problem remains intrinsically high-dimensional and intractable with traditional MCMC methods.
Hamiltonian Monte Carlo (HMC) is a gradient-based MCMC algorithm that avoids random-walk behavior by making better proposals using Hamiltonian dynamics (Duane et al. 1987;Brooks et al. 2011).By using the gradient of the log-posterior to inform the movement of the particle, HMC avoids staying in the same position for long periods due to bad proposals, and ensures that subsequent proposals are distant in parameter space.HMC provides numerous benefits over traditional MCMC algorithms, particularly in rending high-dimensional inference problems feasible due to its fast, scalable algorithm.For a more detailed introduction to HMC and a discussion of its benefits, we refer the reader to Gelman et al. (2004) and Betancourt (2017).
For the practical use of HMC into inference problems, there are two main challenges.First of all, HMC requires the gradient of the likelihood with respect to all inference parameters, which is computationally expensive and not always available in the analytical form, although automatic differentiation (Baydin et al. 2015) makes this tractable.We will discuss our approach to rendering our model likelihood automatically differentiable in Sec.3.4.The second drawback of HMC is that the gradient, through the geometry of the target distribution depends on the specific parameterization of the problem.This dependence means that the step size, the Euclidean metric (which accounts for linear correlations in the posterior), and the number of steps in each HMC iteration need to be tuned manually to ensure good and efficient sampling of the parameter space.The tuning of these hyper-parameters is highly non-trivial and also model dependent.In the worst case of a poorly tuned HMC, the particle may repeatedly loop back in a way that subsequent samples are very close together, which is known as the "U-Turn problem" (Hoffman et al. 2014).
The No-U-Turn sampler (NUTS) is an adaptive extension of HMC with several self-tuning strategies that remove the need to manually select HMC parameters (Hoffman et al. 2014).As a result, the end user normally only needs to worry about the design of their model and not about computational inefficiencies or the detailed implementation of the NUTS itself.It is the advent of NUTS and its efficient implementation in probabilistic programming packages like PyMC (Salvatier et al. 2016) and Stan (Carpenter et al. 2017; Stan Development Team 2023) that make the inference of our hierarchical model and other recent works (e.g., Shen et al. 2022) a tractable problem.
In this work, we make use of PyMC, which is an open-source probabilistic programming package written in Python for Bayesian inference.It implements many inference algorithms including the aforementioned Metropolis-Hasting sampling, slice sampling Neal (2003), and NUTS.Treating the NUTS implemented in PyMC as a black-box can be insufficient for certain models, since the selftuning strategies adopted in NUTS are not always guaranteed to work.The specific parametrization of the problem and the details in the numerical implementation can still impact or even determine the quality and efficiency of the sampling process.As discussed in Sec.3.1, we choose Δ and Δ R as the sampling parameters instead of  and  R to decorrelate the radial coordinates of individual stars from those of the GC center.If we sample in  and  R , the sampled posteriors of these highly-correlated radial coordinates behave badly with additional nonphysical oscillatory features, though the mean and variance of the estimated posteriors remain sensible, which motivates our use of Δ and Δ R .In addition, there are certain cases when the inference of our model fails with the NUTS implemented in PyMC, which we will discuss in Sec.4.1 and 5.

Emulation of the lowered isothermal DF
As discussed in the previous section, HMC requires the gradient of the likelihood with respect to all model parameters, which can be difficult and expensive to obtain for functions without simple analytical derivatives.For our model likelihood in Eq. 15, it is difficult to obtain derivatives with respect to terms involving the lowered isothermal DF, which is obtained through the limepy code by numerically solving the Poisson Equation.For every different set of  lp , we have to resolve the Poisson Equation, which is a time-consuming process for performing inference using MCMC.Obtaining reliable numerical derivatives on limepy DF through simple methods like finite difference is challenging for multi-parameter numerical functions.Luckily, automatic differentiation (AD) allows one to surpass this challenge and obtain reliable derivatives for complicated numerical programs.
In general, AD techniques can transform a program that calculates numerical values of a function into a program that calculates numerical values for derivatives of that function, with about the same accuracy and efficiency as the function values themselves (Bartholomew-Biggs et al. 2000;Baydin et al. 2015).Using AD requires one to write the numerical program in specialized packages and frameworks, and in this work, we use JAX (Frostig et al. 2018), a Python framework capable of automatically differentiating native Python and NumPy functions.In addition, JAX allows just-in-time (JIT) compilation (Sabne 2020) of Python codes, which can significantly speed up the evaluation of Python programs.In our case, JIT allows quicker evaluation of our model likelihood in Eq. 15, thereby speeding-up the HMC inference process.
To render the limepy DF differentiable, we would therefore need to rewrite the limepy code in JAX, which is a highly non-trivial task requiring major code restructuring.The latter would be needed to eliminate branching, and significant efforts would be required to ensure the reliability of gradients under numerical ODE solvers.This is beyond the scope of this work.However, for the single-mass, isotropic limepy models, we can instead build a simple emulator of the lowered isothermal DF through linear interpolation, which we describe next.
In the DF for the single-mass, isotropic model (Eq.1), ,  t , , and  need to be determined by solving Poisson's equation (Eq.4).Note that ,  t , and  are functions of  lp , while  is a function of both  lp and .These functions only depend on the size  h and mass  total of GCs through the following scaling relations: Therefore, to bypass the ODE solver, we (i) linearly interpolate Â(Φ 0 , ), rt (Φ 0 , ), ŝ(Φ 0 , ), and φ(Φ 0 , , ) evaluated at  h,ref = 1 pc and  total,ref = 10 5  ⊙ on a pre-computed rectangular grid with size of 200 × 200 × 1250 forΦ 0 ∈ [1.4,16],  ∈ [0.001, 3.49], and  ∈ [0, 69]4 variables respectively using the limepy code, (ii) obtain the true values of ,  t , , and  at certain  h and  total with the scaling relations in Eq. 24, and (iii) substitute the interpolated and re-scaled values of ,  t , , and  into the analytical functions in Eq. 1-3 to calculate the limepy DF  lp (, | lp ).
Our emulator is now entirely based on analytical functions, and the gradient of the emulator with respect to any parameters can then be reliably obtained with AD in JAX.For functions with only two or three variables, linear interpolation is both accurate and practical with well-understood error properties compared to more complicated approaches that might rely on, e.g., neural networks.
To show the performance and accuracy of our interpolation-based emulator, we compare evaluations of the logarithmic lowered isothermal DF log  (, ), which was the quantity used in MCMC, from both the original limepy code and our emulator.As seen in Fig. 3, our emulator agrees well with the original function and has errors on the order of or smaller than 0.1% across different parameter combinations and ranges.In particular, our emulator reproduces accurate DFs even when  and Φ 0 are close to the boundary values used in the interpolation, as seen in Fig. 3a.By using JIT in JAX and using interpolation to avoid solving ODEs, we achieve a speed-up of more than 50 times for the likelihood evaluation step.Evaluating the DF from limepy code with a given set of (Φ 0 , ,  tot ,  h , ì , ì ) parameters takes about 0.06 s, while the same evaluation from our emulated code takes only about 0.001 s.This speed-up enables our HMC sampling to be completed in about one or two hours instead of days on a single CPU.
To show that our emulator is sufficiently accurate for performing Bayesian inference, we use both the limepy code and our emulator to perform the non-hierarchical Bayesian inference of GCs introduced in EWR22.Assuming complete and accurate knowledge for both the phase-space information of all  stars in the GC and the phasespace center of the GC, we want to infer  lp , the four structural parameters of the GC.The likelihood of this simpler model is given in Eq. 9, while we choose Eq. 19-22 as the model prior, which is the same prior for our hierarchical model.To also test the performance of HMC, we use slice sampling (Neal 2003), a popular alternative to the Metroplis-Hasting algorithm for non-hierarchical Bayesian inference, when we run the limepy code for likelihood evaluation; we rely on HMC sampling when our emulator is used to calculate the lowered isothermal DF.The results of MCMC runs under these two settings for the two perfect simulations are shown in Fig. 4. We see that the inference results using the emulator+HMC method agree with the baseline results from the limepy code+slice sampling method in the two cases shown in Fig. 4 (the overlapping blue and red contours), demonstrating the effectiveness and accuracy of our emulator for the limepy code.The limepy code+slice sampling method takes about an hour to sample 4000 points in a single chain for inferring the four GC structural parameters, while running HMC with our emulated code takes only about two to three minutes for sampling the same number of points.In addition, we also wish to highlight the partial degeneracy between Φ 0 and  based on the inference results shown in Fig. 4: a decrease in  can be partially compensated by an increase in Φ 0 to preserve the GC morphology.

RESULTS
We report the inference results on the simulations described in Sec. 2 using our hierarchical model.Given that the likelihood is defined by the DF that was used to generate the data, we expect to obtain reasonable parameter estimates through inference made from the posterior distribution using MCMC.However, given the complexity of our  Figure 4. Bayesian inferences of (Φ 0 , , log 10  tot , log 10  h ) GC structural parameters for two different GC simulations of 1000 stars generated using the limepy code.For each simulation, we assume the 3D radius and velocity of each star from the GC center are completely and accurately known without any measurement errors.The red contours (showing 50%, 75%, and 95% credible regions) are the inference results obtained with HMC sampling, and the likelihood is evaluated using our interpolation emulator of the lowered isothermal DF.The blue contours are obtained with Slice Sampling, and the likelihood is evaluated using the original limepy code.The dashed line gives the underlying correct parameter value used for generating the simulations.We see that the two runs with different DF evaluation methods and sampling methods agree with each other for both simulations, which suggests that our emulator is sufficiently accurate for performing Bayesian inference.models (∼6010 parameters, numerical DFs, and non-trivial coordinate transformations) and the use of techniques including emulation, auto-differentiation, and HMC, it is critical to validate our model performance on simulations, in which case the true parameters are known.We are also going to impose prior distributions that are at least weakly informative, and so it is good practice to test whether the posterior can still be used to reliably infer the model parameters.Furthermore, we can explore the impacts of different observational settings on inference results through changing the hyperparameters (e.g.number of stars observed , distance of cluster centers  c ) described in Table 2, which further validates the performance of our model.

GC Structural Parameters
Here we analyze and summarize the inference results of the 32 × 10 simulations for the 32 sets of structural parameter values, which were given in Table 1 of Sec. 2. With the same underlying structural parameters, the 10 simulations are different realizations of the same distribution function achieved through random sampling of stars.Using one CPU per chain, it takes about one to two hours for HMC to sample 4000 points in each chain under our hierarchical model, which is significantly more time-consuming than sampling the nonhierarchical model described in Sec.3.4.Further speed-up can be achieved by using multiple CPUs per chain or GPUs for DF evaluations and sampling.By default, we run four chains for each HMC run and discard the first 50% of the samples for burn-in and adaptive tuning.
We notice that one or two out of the four chains can get stuck in some local minimum (four chains do not agree), which occur roughly one fourth of the time out of all HMC runs.One particular noticeable local minimum exists around Φ ≈ 14, where the sampling often gets stuck.Rerunning HMC one or two times can usually fix the problem and produce converging inference results for all four chains.However, there are 10 simulations (out of 320) where one chain is persistently trapped in some local minimum (disagreeing with the other three chains) even after multiple reruns of HMC.For these 10 simulations, we simply discard the chain stuck in the local minimum and keep the other three chains (which are consistent with each other) as the inference results.
In addition, the NUTS algorithm in PyMC fails to produce any results for 3 simulations (out of 320).In all three cases, the adaptive tuning of the mass matrix fails and the diagonal of the mass matrix contains zero value, which halts the sampling process.Since this problem occurs relatively infrequently, we simply replace these three simulations with new simulations generated under different random seeds.We will address the local minimum problem and the adaptive tuning failure of NUTS in future works.
We also notice that when the simulations were generated with Φ 0 = 2 (not as part of the 32 × 10 simulations analyzed in this section), which is close to the lower boundary 1.5 for the prior of Φ 0 , the sampled posterior of Φ 0 can become overly non-Gaussian with features such as small peaks or oscillations in the shape of the 1D probability density function, despite the accuracy of our emulated DF at Φ 0 = 1.5 as shown in Fig. 3a.To address this sampling issue for Φ 0 values near the boundaries, we might need to replace the uniform priors in Eq. 19 and 20 with smooth distributions that do not possess sharp boundary cutoffs.
In general, our HBI gives reasonable parameter estimates for all 320 simulations outlined in Table 1.We show the 95% credible intervals of the GC parameters for two different sets of structural  lp , as examples, in Fig. 5.We summarize the inference results of all 320 simulations in Fig. 6 by plotting the parameter recovery rate on a blue scale.This is an estimation of the frequentist coverage probability of our Bayesian credible intervals; it represents the number of times out of 10 that the 95% credible interval overlaps the true  lp values for each parameter set.The 95% credible interval recovers the true parameter value 9 or 10 times out of the 10 simulations in most cases.
Overall, we take these results to mean that our emulator for limepy is working well, and that except for perhaps some extreme cases, the credible regions are reliable.In general, Fig. 5 and 6 demonstrate that our hierarchical model can give robust and reliable estimates of GC properties across different values for GC structural parameters.

Hyperparameters
Here we explore the impact of different observational conditions on the inference results to further showcase the performance of and validate our hierarchical model.The different observational conditions are reflected by the hyperparameters listed in Table 2.Note that we use the same structural parameters for  lp and the same random seed to generate all the simulations analyzed in these section.
Increasing the number of observable stars in a GC significantly improves the constraints for all GC parameters (Figure 7).At the same time, we note that meaningful and relatively accurate estimates of the GC properties can still be obtained with only 200 uniformly selected stars due to the use of both spatial and kinematic information in our model.
In Fig. 8, we show the inference results for simulations placed at different radial distances  c .We see that the posterior distributions on Φ 0 , , and  h are relatively insensitive to the change of  c , while increasing distance noticeably enlarges the uncertainty of the mass estimate (the first column of Fig. 8a).This behavior is expected; the morphology and shape of GCs are relatively robust to the measurement error on stars' positions and velocities in the Cartesian coordinates when we have hundreds of stars, while  ∝  2 is sensitive to the error on , which depends on the observer's distance to GC.With the same uncertainties for the angular proper motions of individual stars, a larger radial distance  c,true leads to a larger error on  2 and therefore a larger error on .In addition, there is a noticeable degeneracy between the GC radial distance and the GC mass (seen in the lower left corner of Fig. 8a): a lower radial distance estimate leads to a lower mass estimate, which suggests that any bias in estimating the radial distance to the cluster center can lead to a bias in estimating the GC mass.The measurement of the half-mass radius  h is also somewhat degraded, especially when  c,true increases from 5kpc to 10kpc.This suggests the importance of including the radial distance parameter  c into the model, and our hierarchical model provides the clear benefit of robustly including the uncertainty of the parallax of each individual star, thereby inferring the GC mass, GC radius, and the radial distance to the GC center at the same time and properly accounting for the degeneracy among these parameters.
In Fig. 8b, the uncertainty of the angular GC center parameter decreases as the distance increases, which is expected since a GC with a fixed physical size occupies a smaller area of the sky as its distance increases.However, there is no correlation between the radial distance estimate and the estimates of the angular coordinates and proper motions, which is evident from the last row of Fig. 8b.As shown in the 1D posterior plot (the lower right corner of both subplots in Fig. 8) for the ( c −  c,true )/ c,true parameter, which indicates the relative measurement error of the radial distance to the GC center, our HBI procedure can reliably recover the distance of the cluster center regardless of distance, but that the relative precision of the cluster distance estimate decreases as the true distance increases.This behaviour reflects the fact that the same parallax measurement errors become less informative as we increase the radial distance (decrease the parallax value).
We next show the impact of different , measurement errors of angular positions, parallax, and proper motions, in Fig. 9.For the structural parameters presented in Fig. 9a, an increase of measurement errors from 0.02 to 0.1 (in units of mas or mas/year) make little difference to the parameter constraints, while the uncertainty of the mass estimate noticeably increases when the measurement errors increase from 0.1 to 0.5.The constraints of GC angular positions and proper motions are relatively insensitive to the measurement errors (Fig. 9b), while the uncertainty of the radial distance estimate increases proportionally with respect to the measurement errors, mostly due to the dependence of radial distance estimates on the measurement errors of parallax.
By utilizing the posterior distribution estimation of a single GC, we can not only estimate parameters but also functions such as the density profile and the cumulative mass profile.The procedure we use to estimate the density profile is the same as the one outlined in Eadie & Jurić (2019) for estimating the Milky Way's cumulative mass profile.We apply this procedure to each set of GC structural parameters sampled through HMC to calculate the density profile () using the limepy code.As we have thousands of samples in our Markov chain, we generate thousands of density profile estimates.These estimates offer both a visual and quantitative representation, enabling us to compute Bayesian credible regions and to compare the estimates directly to the GC's true density profile.Another quantity of interest that we can estimate using the posterior distribution for a single GC is its velocity dispersion.
We plot the density and velocity dispersion profile constraints of the simulated GCs for different number of stars  in Fig. 10, different distances  c in Fig. 11, and different measurement errors  in Fig. 12.For these selected clusters, we recover both the density and velocity dispersion profiles very well through our HMC inference.The accuracy of density and velocity dispersion profiles reflect the accuracy of our constraints on  lp : in Fig. 10, we see that the constraints of both density and velocity dispersion profiles are significantly improved when the number of stars increases, which increases the amount of available information about the GC and thereby improves the constraints on  lp .In Fig. 11 and 12, the constraints of () are comparable for different  c and , while the uncertainties for the velocity dispersion profiles  2 noticeably increase as  c and  increase.Since  ∝  2 , the amplitude of the velocity dispersion profiles  2 indicates the value of  total .The uncertainties in  2 () therefore reflect the degradation of constraints on  total , illustrated in Fig. 8a and 9a.
In general, our hierarchical model performs well for simulations with different hyperparameters reflecting different observational conditions.The samples approximating the posterior distribution both (1) recover the underlying true parameters, and (2) give reasonable uncertainty estimates that are consistent with the assumed observation conditions.Together, these two achievements demonstrate the power of HBI to incorporate measurement errors and missing data into parameter estimation.

CONCLUSIONS
Using the lowered isothermal DF proposed by GZ15, we build a hierarchical Bayesian model for GC inference that can directly take data in the heliocentric equatorial coordinates (the reference frame in which we obtain measurements) and can incorporate incomplete data and measurement errors.We extend the basic Bayesian modelling of GCs in EWR22 into the hierarchical version and further include the inference for the positions and velocities of the GC center into the hierarchical model.To render the inference process computationally feasible, we build an interpolation-based emulator for the singlemass, lowered isothermal DF, written in JAX, which allows autodifferentiation, and we run Hamiltonian Monte-Carlo to estimate the posterior.We test our hierarchical model on hundreds of GC simulations with different parameters.We show that our hierarchical model can give accurate and reliable estimates of GC properties across different GC morphologies and distances, and our approach gives robust uncertainty estimates that incorporate measurement errors and are consistent with the assumed observation conditions.There are still several technical problems related to running HMC that need to be addressed in order to allow efficient and accurate hierarchical Bayesian inference on the entire parameter range.In addition to addressing the local minimum problem and updating the uniform priors with smoother priors to avoid sharp boundary cutoffs, which we have already discussed in Sec.4.1, we also face difficulty with using HMC to sample posteriors for simulations with  ≲ 0.6: the chains fail to converge due to issues with adaptive tuning and produce unstable inference results.Even though most of the Galactic GCs have been fitted with  ≥ 0.8 where our code does produce reliable results, there is still a sizable fraction of GCs fitted with  ≲ 0.6 (de Boer et al. 2019;Dickson et al. 2023).It is therefore important to continue our investigation into the parameter specification of our model and the auto-tuning procedures adopted in NUTS so that we can perform HBI in the entire prior range of (Φ 0 , ).
Another technical area for further efforts is to improve the speed of the HMC sampling.Currently, it takes one to a few hours to sample hierarchical models with 6010 parameters for 1000 stars, using one CPU per chain during HMC.This speed is sufficient for including stars with kinematic information, which range from  (10) −  (10 3 ) number of stars for GCs in the Gaia data (Vasiliev 2019).However, in a real GC there can be many more stars with known angular positions but without any kinematic information.To include all these stars, one In the left subplot, we see that the uncertainty of radial distance has a noticeable impact on the uncertainty of the mass estimate and, to a less extent, the half-mass radius estimate for the simulated GC.A lower radial distance measurement correlates with a lower mass and a lower radius estimates for GCs with a radial distance at 5 and 10 kpc.
The uncertainty of the angular positions and proper motions significantly increases with a decreasing distance, but there is no obvious correlation between the radial distance and the angular coordinates, as seen in the right subplot.The constraints of the density  and velocity dispersion  2 profiles from HMC runs on simulations with different radial distance  c .The solid black curves show the actual velocity dispersion profile evaluated at the true structural parameters that are used to generate these simulations.We set the y axis with the same range for all the panels in each row to allow easy comparison.The constraints for ( ) are similar for different  c , while the uncertainties for  2 increase as  c increases.
would need to infer very large models with  (10 4 ) numbers of stars, which renders the speed-up of HMC necessary.Such speed-up can be achieved by just-in-time compilation, GPU-acceleration, and further parallelization of the HMC sampler.
For the likelihood of our model in Eq. 9, the independence and uniformity assumption does not take into account of the selection bias occurred in real surveys such as Gaia and HST due to sourcecrowding, magnitude limits, the contamination of field stars, and other observational effects.EWR22 has emphasized the importance of selection function modelling for real survey data on GCs in order to avoid any severe inference bias.Such selection function modelling is naturally addressed in a hierarchical framework, and we plan to extend our model for observed data in future work.This is especially relevant for cases where biased sampling also requires combining observations from different datasets.However, our HBI model can already be used to fit real observational data when the sampling of stars across the clusters are uniform, and we aim to apply our HBI model to GC data compiled in de Boer et al. (2019) in future work.
With the ability to properly address selection functions and combine different datasets, the natural next step to extend and test our hierarchical model will be to include radial velocity measurements, which are available for a sizable number of stars in a large fraction of Milky Way GCs from ground-based spectroscopy (Abdurro'uf et al. 2022;Yan et al. 2022).It is also important to compare the method presented here to traditional methods in the literature that (1) use the projected distances of stars to estimate density and mass profiles with radial binning, and (2) combine data sets from different telescopes to use stars at all radii.We have so far only been able to incorporate the single mass, isotropic, lowered isothermal DF into our hierarchical model, which limits the applicability and flexibility of the model and can lead to significant biases on inferred GC properties such as mass and radius due to ignoring dynamical effects such as mass segregation (Sollima et al. 2015;Hénault-Brunet et al. 2019).In order to explore interesting GC features such as dark remnants, binaries, and mass segregation and to achieve more realistic GC modelling, future work should focus on building emulators for the multi-mass, anisotropic, lowered isothermal GC model, therefore allowing the hierarchical inference to utilize the full power of the limepy code.However, building emulators of the anisotropic, multi-mass model is more difficult due to a substantial increase in the complexity of the model and the number of parameters, which makes the linear interpolation process less efficient, if not impossible.One can perhaps build upon the simple interpolation approach laid out in Sec.3.4 for the anisotropic, multi-mass model or use emulators based on machine learning.Alternatively, one could rewrite the full limepy code to make it auto-differentiable so that we can obtain reliable gradients for HBI.
Astronomers are interested not only in the inherent characteristics of globular clusters (GCs), but also in comparing and selecting GC models.The latter is crucial for comprehending the internal dynamics of GCs and the broader narrative of their evolution as they move through the Galactic potential.Being able to associate an observed GC with a particular dynamical model is a critical step in unraveling the cluster's present properties and its evolutionary past.Understanding the distribution function of stars within a cluster enables more intricate features to be explored in greater detail, like its dark remnant population, binary population, degree of mass segregation, and tidal history.By utilizing a model that integrates all of these components and also enhancing the statistical framework to account for observational sampling bias, astronomers can gain a better understanding of the dynamical state of GCs.Knowing a cluster's dynamical state also places constraints on the cluster's properties at birth and how it has evolved over time.By demonstrating the robustness and practicality of HBI, our work lays out the foundation for more realistic and complex GC modelling that has the potential to maximize the cluster's utility as a tool to study the universe around it.

Figure 2 .
Figure 2. Ratio of the dimensionless virial radius to the half-mass radius, r rℎfor models with different Φ 0 and .The blue line shows the boundary corresponding to r / rℎ ⪆ 0.64, which distinguishes the lowered isothermal models (in the red region) that are relevant to modeling GCs.The black dashed-line represents the 10 th -degree polynomial   (Φ 0 ) that we used as the upper boundary in the uniform prior of  for our Bayesian model.Notice that our prior function   (Φ 0 ) cuts off some physically relevant (Φ 0 , ) parameters close to the blue boundary.

Figure 3 .
Figure 3.The evaluations of the lowered isothermal DF  ( , ) in Eq. 1 under four different parameter combinations (subplots a, b, c, and d) by both the original limepy code (red solid line) and our interpolation-based emulator (blue dashed line) introduced in Sec.3.4.The four parameter combinations are indicated by the titles of the subplots and are chosen to represent a broad range of values for  lp .For visualization, we fix  at a certain value in the left panels to show the variation of the DF with respect to , and we fix  in the right panels to show the change of DF with respect to .We plot the values of DF in the log scale, which is the scale used for computation in MCMC, in the upper panels.The relative errors of log  , which are defined as Δ log  /log  ≡ (log  emulator − log  limepy )/log  limepy , are shown in the lower panels.In general, we see excellent agreement between the two methods in all four parameter combinations.Our emulator has an error smaller than 0.1% in most of the parameter regions for log  .The relative error only exceeds 0.1% when the true log  values are close to 0.

Figure 7 .
Figure 7.The impact of different  , the number of observable stars in GC, on inference results.For this and all the following figures, we use the dark and the ligher color to indicate the 75% and 95% credible regions of the parameter posteriors respectively.The difference between the sampled angular position parameters and their true values,  c −  c,true and  c −  c,true , are shown in the unit of arcsec.Increasing the number of stars clearly improve the inference results for all GC parameters.

Figure 8 .
Figure8.The impact of different  c,true , the true radial distance to GC, on the inference results.We use the difference between the observed and the true distances normalized by the true distance, ( c −  c,true )/ c,true , to indicate the uncertainty on the radial distance estimate.In the left subplot, we see that the uncertainty of radial distance has a noticeable impact on the uncertainty of the mass estimate and, to a less extent, the half-mass radius estimate for the simulated GC.A lower radial distance measurement correlates with a lower mass and a lower radius estimates for GCs with a radial distance at 5 and 10 kpc.The uncertainty of the angular positions and proper motions significantly increases with a decreasing distance, but there is no obvious correlation between the radial distance and the angular coordinates, as seen in the right subplot.

Figure 9 .Figure 10 .
Figure9.The impact of different , measurement errors of angular positions, parallax, and proper motions, on inference results.The uncertainty of the mass estimate noticeably increases when the measurement errors increase from 0.1 to 0.5 (in units of mas or mas/year).The uncertainty of the radial distance estimate increases proportionally with respect to the measurement errors.

Figure 12 .
Figure 12.The constraints of the density  and velocity dispersion  2 profiles from HMC runs on simulations with different measurement uncertainties .

Table 1 .
Limepy model parameter values used to simulate stars in GCs.