© The Authors 2024

1 Introduction

Searching in the vicinity of the Milky Way, we find mostly spherically shaped stellar systems with an increasing density of stars toward their center. These are called globular clusters (GCs). They are commonly known for their very old age, which reaches ∼10 Gyr in most cases or even more.

GCs were thought to consist of a single stellar population with low metallicity. However, the detection of chemical and photometrical differences first suggested and then confirmed the presence of different generations of stars. A variety of questions were raised, establishing a renewed interest in GCs. Astronomers are studying the possible formation scenarios of such multiple stellar generations, as well as their initial properties, using the enormous amount of recently collected observational data and by developing advanced models (see the reviews by Bastian & Lardo 2018; Gratton et al. 2019; Milone & Marino 2022). Behind their intriguing formation and stellar generations, GCs have a highly articulated dynamical evolution. Their high-density environment makes them the only known collisional star-cluster systems in astronomy. The relaxation process driven by gravitational encounters between stars plays a fundamental role in determining the structure and evolution of GCs following their birth. The efficiency of relaxation brings these systems toward a degree of kinetic energy equipartition among stars (Spitzer 1987), altering their structural and kinematic features. Such internal dynamical processes affect the spatial distribution of stars, a phenomenon known as mass segregation. A system close to equipartition has massive stars with a lower velocity dispersion, which consequently sink toward the center. On the contrary, less massive stars are faster and migrate to the outer regions. Here, the gravitational pull of our Galaxy subtracts stars from the cluster, leading to a limited phase-space domain through what is known as evaporation (Ambartsumyan 1938; von Hoerner 1957; King 1958a,b; Spitzer 1987), and producing the tidal streams (Odenkirchen et al. 2001; Piatti & Carballo-Bello 2020). The joint effect of mass segregation and evaporation flattens the mass function, whose initial shape is largely unknown. Our comprehension of the internal dynamics of GCs can provide hints concerning their initial mass function (IMF). The mass loss from GCs following their formation also contributes to populating the Galactic environment. Depending on their orbits, they can suffer shocks and inhomogeneities in the still unknown early galactic potential (Webb et al. 2019). Regarding multiple stellar generations, these dynamical processes, coupled with stellar evolution, favor a strong and predominant loss of first-generation stars in the early evolutionary phases of clusters. The evaporation of such stars may explain the origin of galactic halo stars. Furthermore, the internal dynamics mixes these populations, mainly erasing their initial spatial and kinematic differences, further complicating attempts to understand the characteristics of the system. Studying GC dynamics can thus reveal fundamental insights concerning the early environment of our Galaxy, as well as their contribution to its evolution (see Chapter 8 of Gratton et al. 2019).

Although it is well known that massive stars are more centrally concentrated than less massive ones, a theoretical description of energy equipartition, mass segregation, and evaporation is still missing and astronomers are trying to quantify these processes. A satisfying explanation of this intricate stellar dynamics is required to extract information about the primordial environment and the formation scenario of GCs. Understanding the stellar dynamics may also help to solve the puzzle of the formation of multiple stellar generations, describing their dynamical evolution. However, even single stellar population dynamical phenomenology still needs to be fully understood.

Focusing on energy equipartition, several scientists are exploring GC evolution through N-body simulations to study the efficiency of this process and its relation to structural and internal properties, as well as initial conditions (Trenti & van der Marel 2013; Webb & Vesperini 2017; Pavlík & Vesperini 2021b,a, 2022). Among such works, the most important from our perspective is that of Bianchini et al. (2016). Here, with a set of N -body simulations, the authors look at the degree of energy equipartition through the velocity dispersion as a function of stellar mass in the inner regions of clusters (namely inside the half-mass radius). These authors find an analytic formula to fit the velocity dispersion–mass relationship that is $$\sigma (m)=\{\begin{array}{ll}{\sigma }_{0,\hspace{0.17em}\text{B}}\exp \left(-\frac{m}{2m_{\text{eq}}}\right)\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}m\le m_{\text{eq}},\hfill \\ {\sigma }_{\text{eq}}{\left(\frac{m}{m_{\text{eq}}}\right)}^{-1/2}\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}m>m_{\text{eq}}\hfill \end{array}$$(1)

where the equipartition mass m_eq is the parameter that quantifies the degree of energy equipartition, while σ_0,B is a normalization constant and σ_eq = σ_0,B exp (−1/2). In a subsequent paper, Bianchini et al. (2018) characterize the variation of m_eq during GC evolution and its dependence on the radial coordinate.

The provided fitting function includes the limit of complete equipartition for masses m > m_eq, which is σ(m) ∝ m^−1/2. Then, m_eq tells us which masses reach complete equipartition, with high-mass stars closer to equipartition than low-mass stars. Furthermore, as the degree of equipartition is expected to increase with the dynamical state of a cluster, the m_eq parameter appears to be related to the concentration parameter c and the numbers of relaxation times elapsed during the life of a cluster, namely N_relax = t_age/t_relax, where t_age is the age of a cluster and t_relax is the relaxation time, which is often evaluated in the core t_rc and at the half-mass radius, yielding the number of core and median relaxation timescales N_core and N_half respectively (Bianchini et al. 2016). The work of these authors has gained attention in the literature due to the simplicity of the analytical function. It was used in N -body simulations to quantify the degree of energy equipartition in GCs hosting multiple stellar populations (Vesperini et al. 2021) and to explore the effect of anisotropy, primordial binaries, tidal field, and black holes in equipartition (Pavlík & Vesperini 2021b,a, 2022; Aros & Vesperini 2023). The exponential fitting function in Eq. (1) is often compared with the classical power law σ(m) = σ_s(m/m_s)^−η with η = 0.5 for complete equipartition, where σ_s and m_s are scale values for the velocity dispersion and mass.

Thanks to Hubble Space Telescope (HST) proper motion measurements, observers can now both measure internal kinematics and estimate stellar masses, quantifying the equipartition level in a few GCs (Libralato et al. 2018, 2019, 2022). In Watkins et al. (2022, hereafter W22), the authors analyzed the degree of energy equipartition in nine GCs using the Bianchini et al. (2016) fitting function and the classical one, with the respective parameters. Most, if not all, of the stars of the studied clusters do not show complete equipartition.

Although the Bianchini function reproduces the observable well, it provides no information on the physics underlying the equipartition process. Our aim is to compare measured data and the Bianchini fitting function with theoretical results, which provide more fundamental information about cluster dynamics. Multi-mass King-like dynamical models have the advantage of being self-consistent, that is they are directly obtained from the physics of the system, namely the gravitational interaction between stars, the mass distribution, and the limited available phase space. As a consequence, the dependence of the velocity dispersion on mass is predicted by the model itself, and is not related to an interpolation function that fits simulation outputs and observations. In this paper, we carefully discuss the relation between the outcomes of our model and N-body-based results, namely the successful Bianchini fitting function. We note that N-body simulations are used to study the temporal evolution of stars in clusters, numerically integrating the Newtonian equation of motion after choosing suitable initial conditions. Current state-of-the-art simulations include a variety of phenomena, such as the treatment of stellar evolution, binary systems, compact object formation, and tidal forces, which make them a powerful numerical tool for simulating more realistic clusters. However, the definition of such recipes and the choice of initial conditions themselves require accurate modeling. Theoretical assumptions and prescriptions can significantly affect the simulation results. This is why most researchers explore the parameters of those models to determine the initial conditions and to describe the considered phenomena. As a result, N-body methods are a modeling tool needed to explore the temporal evolution, but they require theoretical references to consistently set up the simulations themselves. Despite large increases in computational resources in recent years, realistic simulations require large and advanced computer infrastructures as well as computational time, which grows with the square of the number of simulated stars.

On the contrary, exploring dynamical models for GCs is much less computationally expensive. In our approach, a single equilibrium configuration can be obtained in few seconds or minutes depending on the desired precision. Furthermore, these models can be used to predict radial and projected profiles for several observables, such as the surface density (both in mass and in the number of stars), the surface brightness (if a mass– luminosity relation is known), and the velocity dispersion. When using models including a mass distribution, our knowledge of the distribution function characterizes not only the phase-space but also the mass function. These models provide a description of the seven-dimensional space of masses, positions, and velocities. Such theoretical frameworks are clearly limited by the assumptions made, as in our case where isotropic velocities are considered. Furthermore, some theoretical quantities require observational constraints in order to get their numerical values, as is usually the case for physical models. The temporal evolution is assumed as described by equilibrium configurations with different structural parameters, whose values determine the dynamical state. The temporal dependence of such quantities can be addressed with numerical simulations, which offer a unique opportunity in that sense (see Wang et al. 2016, for a relation between the King single-mass model parameter W₀ and advanced simulations). Indeed, observing the historical properties of GCs would require observations of the satellites of distant galaxies (i.e., far enough to look back in time into the life of the GCs) under the strong assumption that cluster evolution is similar in every galactic environment, which may not be the case. Unfortunately, we are currently limited to GCs of the Milky Way and its satellites and nearby galaxies.

Numerical simulations and theoretical modeling can therefore provide a broader perspective as to the dynamics of GCs. However, we strongly emphasize that taking advantage of dynamical models is fundamental to understanding the underlying physics of the processes involved. Such models are an advanced interpretative and predictive tool for addressing the major open questions in the field.

2 Methods

Our main approach consists in comparing the velocity dispersion dependence on stellar mass predicted by our multi-mass dynamical model for GCs, with the Bianchini function in Eq. (1), fitting the observational data by W22. The fitting procedure leads to an estimate of model parameters, to be related with the Bianchini equipartition mass m_eq.

We also look at the surface brightness profiles (SBPs) of the clusters by fitting the measured data by Trager et al. (1995) with our theoretical prediction, which provides an estimate of parameters from a different observable. This allows us to see how confident our model is at reproducing both the observed internal kinematics and luminosity profile.

2.1 Dynamical model

Our dynamical model derives from the distribution function (DF) obtained by King (1965) as an approximated solution of the Fokker-Planck equation, which itself comes from the Boltzmann equation for collisional systems (Chandrasekhar 1943; Rosenbluth et al. 1957; Spitzer & Härm 1958; Chandrasekhar 1960). The King model (King 1966) was a single-mass one, with the DF having a cut-off in the phase-space that takes the Galactic tidal field into account. Such a model describes the equilibrium configuration of the system, once the DF is used to calculate the density profile and solve the Poisson equation for the gravitational potential. It was used for many years by astronomers to fit the SBPs of GCs and to obtain structural parameters (Peterson & King 1975; Trager et al. 1995; Harris 1996; Miocchi et al. 2013). A further step on was done first by Da Costa & Freeman (1976), introducing a discrete multi-mass model with ten mass classes, each having its own DF with an energy cut-off and a weight factor to be constrained by observations. They successfully fit M3 surface brightness, where the King model was failing to reproduce both the inner and outer profile.

Due to the general good agreement between the King model and observations, little effort was put in developing and exploring multi-mass models, although they are fundamental in the comprehension of several phenomena produced or altered by the mass distribution, like segregation, evaporation, and equipartition (Miocchi 2006; Gieles & Zocchi 2015; Peuten et al. 2017; Torniamenti et al. 2019; Dickson et al. 2023). Even if N-body simulations are offering great opportunities in understanding the dynamical evolution of GCs, we underline the importance of developing analytical models to describe physical processes, offering an important tool in the comprehension of such systems, like it was the King DF.

The formulation of a continuous multi-mass model was presented by Merafina (2019), as an extension of the Da Costa & Freeman one and recovering the King formalism. The model we present is similar, but with some improvements concerning the derivation and the relation with the global mass function. Our DF is an approximated solution of a generalized expression of the Fokker-Planck equation, valid for collisional systems with a mass distribution (see Eq. (A.1)). The DF keeps the property of being limited in the phase-space. Moreover, not only it brings information on the velocity distribution within the cluster, but it also concerns the mass distribution of the stars, meaning that the mass function is embedded in the DF. Its expression is $$g(r,v,m)=k(m){\text{e}}^{-m\left[\phi (r)-{\phi }_0\right]/\left(k_{\text{B}}\theta \right)}\left[{\text{e}}^{-\epsilon (v,m)/\left(k_{\text{B}}\theta \right)}-{\text{e}}^{-{\epsilon }_{\text{c}}(r,m)/\left(k_{\text{B}}\theta \right)}\right],$$(2)

where ε = mv²/2 is the kinetic energy of a star with mass m and ε_c = mv²/2 is its cut-off energy with v_e = v_e(r) the escape velocity, φ(r) is the gravitational potential (with φ₀ its value in the cluster center) and r is the radial coordinate. The variable θ is the thermodynamic temperature, a memory of the Boltzmann DF limit and constant all over the equilibrium configuration (Merafina 2017; Merafina 2018, 2019) and k_B is the Boltzmann constant. The multiplying factor k(m) weights the DF of each mass m, like it was in the Da Costa & Freeman (1976) model, although theoretically it gathers some mass-dependent functions resulting from the derivation, related with the mass function and the escape velocity at cluster center (see Appendix A for details).

The mass density radial profile ρ(r) is computed through an integration over the masses and velocities of the DF in Eq. (2), that is $$\rho (r)={\displaystyle {\int }_{m_{\min }}^{m_{\max }}\left[{\displaystyle {\int }_0^{v_{\text{e}}(r)}g}(r,v,m){\text{d}}^{3}v\right]}\text{d}m,$$(3)

with m_min and m_max the extremes of the mass range, where the mass function is valid. Equation (3) is used to solve the Poisson equation for the gravitational potential. Following the King formalism, one could introduce $W(r)={\epsilon }_{\text{c},1}(r)/\left(k_{\text{B}}\theta \right)=m_1{v}_{\text{e}}^2(r)/\left(2k_{\text{B}}\theta \right)=m_1\left[{\phi }_{\text{R}}-\phi (r)\right]/\left(k_{\text{B}}\theta \right)$ and solve the Poisson equation for W(r), obtaining a set of equilibrium configurations each identified by the initial condition W_0,1 = m₁ [φ_R − φ₀] /(k_Bθ), like shown in Merafina (2019). However, the parameter for determining the configuration should not depend on the scale mass (the models by Da Costa & Freeman 1976 and Merafina 2019 were using the greatest mass m₁ of the system). Using a generic mass unit m_u and solving for W_u(r) = m_u[φ_R − φ(r)]/(k_Bθ), the expression of the dimensionless Poisson equation keeps the same analytical expression, that is $$\frac{1}{x^2}\frac{\hspace{0.17em}\text{d}}{\hspace{0.17em}\text{d}x}\left(x^2\frac{\hspace{0.17em}\text{d}{W}_{\text{u}}}{\hspace{0.17em}\text{d}x}\right)=-9\frac{\rho }{{\rho }_0},$$(4)

with the dimensionless coordinate x = r/r_k,u in units of the King radius $r_{\text{k},\text{u}}=\sqrt{9k_{\text{B}}\theta /\left(4\pi G{m}_{\text{u}}{\rho }_0\right)}$. The initial conditions of Eq. (4) are $W_{\text{u}}^{\prime }(x=0)=0$ and W_ɪ(x = 0) = W_0,u = m_u(φ_R −φ₀)/(k_Bθ), which determines the solution. Here, we use a better parameter to identify the equilibrium configuration, that is Φ₀ = (φ_R − φ₀)/(k_Bθ), which measures the potential well without depending on a scale mass (as W_0,u do) and having the dimension of the inverse of a mass.

To solve Eq. (3), the quantities in the DF in Eq. (2) must be known. At the present state of development, the factor k(m) has to be constrained from the global mass function ξ(m) with a numerical procedure. Since the integration of g(r, v, m) over the radial coordinate r and the velocity v gives the mass function by definition, we can iteratively constrain k(m). In this work, we consider a single-power law mass function ξ(m) ∝ m^α , with the slope α taken from Baumgardt et al. (2023) or theoretically assumed. The code that numerically integrates Eq. (3), once given the slope α, the mass range and Φ₀, evaluates k(m) with a convergence procedure and draws radial profiles for each mass composing the system. These include the 3D velocity dispersion profile as function of mass $$\sigma (r,m)=\langle v^2(r,m)\rangle =\frac{{\displaystyle {\int }_0^{v_{\text{e}}(r)}g}(r,v,m)v^2\hspace{0.17em}{\text{d}}^{3}v}{{\displaystyle {\int }_0^{v_{\text{e}}(r)}g}(r,v,m){\text{d}}^{3}v},$$(5)

that we project in two dimensions to obtain σ(R, m), with R the projected distance from the cluster center. This profile is needed to fit the proper motion observations on σ(m).

The model also provides the density profile of each mass, which is used to predict the surface density profile for each mass Σ(R, m), by integrating along the line of sight. This allows us to estimate the surface brightness profile I(R), by introducing a mass-luminosity relation. This is a fundamental step when dealing with multi-mass systems in particular. We go through the details of the production of the theoretical SBPs in Sect. 2.2.2, where we discuss the fitting procedure on observational data.

2.2 Fitting procedure

2.2.1 Velocity dispersion–mass relationship

W22 provide the binned velocity dispersion of stars as function of their estimated mass for a sample of Galactic GCs (see their Table 2). Their dispersion measures are obtained from proper motion (then in the 2D plane of the sky) and inside a ring in the inner regions of the analyzed clusters.

In order to compare the model prediction with the data, we compute σ(m) from the projected profile σ(R, m) in the same region. When looking at a projected radial distance from cluster center R, the observer intercepts all stars having 3D radius r ∈ [R, R_t], where R_t = r_t is the tidal radius. Then, we average the 3D profile σ(r, m) in that range to obtain the 2D profile σ(R, m). Subsequently, the projected dispersion is averaged in the circular annulus covered by the data sample, drawn around the median radius (see W22, Table 1). We use the observed core radius r_c from the Harris (1996) catalog (2010 edition, hereafter indicated as the Harris catalog) and the theoretical core radius r_c,th to convert the dimensional coordinates into our dimensionless ones. Here, r_c,th is obtained from the surface density profile, predicted by our model. This coincides with r_c if the surface brightness profile is proportional to the surface density profile (with a factor only dependent on mass).

The overall average procedure returns the velocity dispersion as function of stellar mass, normalized to the dispersion of the lowest mass σ₀ . The shape of σ(m)/σ₀ depends on the model parameters. In particular, for an increasing Φ₀ that gives a more advanced dynamical state, the dispersion σ(m) is closer to the complete energy equipartition limit σ ∝ m^−1/2. Furthermore, more massive stars have a greater degree of equipartition than less massive ones for any value of Φ₀ (see Fig. B.5 available on Zenodo).

Concerning the choice of the mass function shape needed to draw theoretical profiles, we set m_min = 0.1 and m_max = 1.0, while for the slope α we select two approaches: the first takes it from Baumgardt et al. (2023), the second one (more general) explores a wide range α ∈ [−2.0, 0.0], actually adding the slope to the parameter space of the minimization procedure. As a consequence, the first approach gives the best-fit values of σ₀ and Φ₀, while the second also returns an estimate for α. For the fit, we numerically minimize the χ² test value between data and model prediction. We apply this procedure to the GCs in our sample, only removing NGC 2808 following W22. We also fit the velocity dispersion data with the Bianchini function in Eq. (1), to have an estimate of the χ² , to be compared with our model ones. Although such a fit is less detailed than in W22, we obtain the Bianchini equipartition mass m_eq and the normalization constant σ_0,B .

2.2.2 Surface brightness profiles

An extensive catalog of SBPs for Milky Way GCs was released by Trager et al. (1995). Although more accurate measurements are available from the work conducted by Noyola & Gebhardt (2006) with HST photometry, they are restricted to internal regions and to a few subsets of clusters. For these reasons, we choose the Trager et al. (1995) dataset. To analyze the profiles, we follow a procedure similar to that described by McLaughlin & van der Marel (2005) and Zocchi et al. (2012), which these authors used to estimate uncertainties on the data and fit SBPs with single-mass dynamical models. The Trager et al. (1995) SBPs consists of measurements for the logarithm of the projected radial coordinate R_i, the surface (apparent) magnitude in the V -band µ_V (R_i), as well as its best-fit with Chebyshev polynomials and a weight w_i for each point according to its by-eye reliability. The analysis performed by McLaughlin & van der Marel (2005) provides an estimate of the uncertainty for each surface brightness measure, namely δµ_V(R_i) = σ_µ/w_i, where σ_µ is a constant that varies from cluster-to-cluster. Following Zocchi et al. (2012), we remove data points with w_i < 0.15 and we correct each measure for the extinction as derived from the Harris catalog.

In the context of multi-mass dynamical models, the standard assumption of a constant mass-to-light ratio is too crude. To obtain the theoretical profile for the surface brightness, useful to fit the available data, we need a mass-magnitude relation in the V -band. Using theoretical isochrones we can establish such a relationship by assigning to the masses of the dynamical model a value for the corresponding absolute magnitude M_V . We use BaSTI isochrones (Hidalgo et al. 2018; Pietrinferni et al. 2021; Salaris et al. 2022; Pietrinferni et al. 2024), considering an age of 13 Gyrs, [α/Fe] = +0.4, the He mass fraction Y = 0.247 and a different metallicity [Fe/H ] for each cluster as a reference case. The metallicity value is taken from the Harris catalog, where we also get the distance modulus to convert the observed apparent magnitudes into absolute ones. Due to stellar mass loss, the final mass of each star is naturally different from its initial mass, depending on its evolutionary state. In particular, initially more massive stars are now much brighter than low-mass stars, even if their current mass may be similar. Furthermore, since GCs are very old, and excluding remnants, the most massive stars have around 0.8 M_⊙ and their mass loss occurred mainly in the last million years. As a consequence, they did not have the time to adjust their dynamical state, which is better described by their initial mass. We then use the initial mass as representative of the dynamical mass, in order to assign the corresponding M_V value. We note that when the stars evolve beyond the main sequence, a large variation in luminosity occurs for a very small variation in their mass. These evolved stars are also the more massive ones still alive in the cluster. For these reasons, they play an important role in shaping the surface brightness profile in the more central regions of clusters, because of the mass segregation process. To account for these effects, we choose a regular mass step, but substantially smaller than that used to determine the velocity dispersion. We note that the mass range, in particular the maximum mass of stars still alive, is provided by the isochrones and depends on the age of the cluster. While this maximum mass has a little effect on the prediction of σ(m), it will play an important role in shaping the surface brightness profile.

Finally, we follow Zocchi et al. (2012) regarding the minimization procedure of the χ², by converting both the observed and the theoretical SBPs into solar luminosities.

Concerning the theoretical prediction, the model gives the surface density profile for each mass Σ(R, m). The global SBP I(R) is defined as $$I(R)={\displaystyle {\int }_{m_{\text{min\hspace{0.17em}}}}^{m_{\text{max\hspace{0.17em}}}}\frac{L_{\text{V}}}{m}}\text{Σ}(R,m)\text{d}m,$$(6)

where L_V is the luminosity obtained by converting the theoretical absolute magnitudes M_V to solar luminosities, using M_V,⊙ = 4.83 (see the Sun Fact Sheet by NASA). Although M_V,⊙ may suffer from several non-negligible uncertainties, as noted by Torres (2010), its choice is irrelevant in our fits.

The computed theoretical profile is normalized to its central value I₀. The shape of I(R)/I₀ depends on the parameters of the model, and its values are given for the dimensionless radial coordinate. However, in our fitting procedure, we must be aware that the observed radial coordinate is in physical units. To deal with all these, we define a numerical auxiliary function h(R/R_c; Φ₀,α, m_min, m_max) = I(R)/I₀ where R_c is the core radius, α is the mass function slope that we take from Baumgardt et al. (2023), m_min and m_max are the mass extremes given by the isochrone. We numerically minimize the χ² and obtain an estimate of the best-fit Φ₀, I₀ and R_c . On the theoretical ground, since the 2D and 3D radial grids are the same, r_c = R_c and r_t = R_t. We apply this method to each cluster for which we have analyzed the velocity dispersion–mass relation. In this way, we can obtain two independent estimations of the Φ₀ parameter.

2.3 Error estimation

Our fitting procedure is based on the χ² minimization with respect to parameters. In our modeling for the velocity dispersion-mass relation, such parameters are Φ₀ and σ₀ for the first approach, while the second one also includes the slope α. The minimization of the χ² requires that дχ2/дa_j = 0 and д²χ²/дa_iдa_j < 0 where a_j represents the parameters. We obtain such conditions with a numerical procedure. For the uncertainties, we need to compute the matrix M_ij = (1/2) д²χ²/дa_iдa_j, whose inverse M⁻¹ is the covariance matrix, from which we obtain parameters errors $\delta {\text{Φ}}_0={\left(M^{-1}\right)}_{00}^{1/2}\text{\hspace{0.17em}and\hspace{0.17em}}\delta {\sigma }_0={\left(M^{-1}\right)}_{11}^{1/2}$ as well as error bands for the velocity dispersion. A similar procedure is done when the slope is considered as a free parameter.

Our theoretical prediction for the velocity dispersion as function of stellar mass depends on parameters. We can write that σ(m; Φ₀, σ₀) = σ₀f(m; Φ₀) where f(m; Φ₀) is the numerical relation obtained from the model (that generally also depends on α, fixed in the first approach). To obtain the matrix M, we must compute the derivative of the χ² with respect to parameters. Specifically, we have an analytical expression for the first and second derivatives with respect to σ₀. Concerning the derivatives with respect to Φ₀, we use the Finite Difference Methods for their approximation. For the first derivative, we use $$\frac{\partial \sigma }{\partial {\text{Φ}}_0}=\frac{1}{2}\frac{\sigma \left({\text{Φ}}_0+\text{Δ}{\text{Φ}}_0\right)-\sigma \left({\text{Φ}}_0-\text{Δ}{\text{Φ}}_0\right)}{\text{Δ}{\text{Φ}}_0}+\mathcal{O}\left(\text{Δ}{\text{Φ}}_0^2\right),$$(7)

where ΔΦ₀ is the numerical step in Φ₀ used to compute the function f . Regarding the second derivative, we take $$\frac{{\partial }^2\sigma }{\partial {\text{Φ}}_0^2}=\frac{\sigma \left({\text{Φ}}_0+\text{Δ}{\text{Φ}}_0\right)-2\sigma \left({\text{Φ}}_0\right)+\sigma \left({\text{Φ}}_0-\text{Δ}{\text{Φ}}_0\right)}{\text{Δ}{\text{Φ}}_0^2}+O\left(\text{Δ}{\text{Φ}}_0^2\right).$$(8)

Both approximations are considered at second order. When the slope is added to the parameter space, we have σ(m; Φ₀, α, σ₀) = σ₀ f (m; Φ₀, α). For the derivatives with respect to the slope, we consider the same Finite Difference Eqs. (7) and (8) for their approximation, using the numerical step Δα. The mixed derivatives are obtained similarly, applying the same method.

Concerning the fitting procedure on SBPs, we follow the same approach. We determine the covariance matrix by computing (and inverting) M_ij = (1/2) д²χ²/дa_iдa_j where the parameters are Φ₀, R_c and I₀. Except for the latter that has an analytical expression, the first, and second derivative of the χ² with respect to Φ₀ and R_c need the derivatives of I(R) with respect to the same parameters. These are approximated by using the Finite Difference Method, similarly to what already described, as well as for mixed derivatives. We obtain uncertainties δΦ₀ , δI₀ and δR_c , from which we can also compute the uncertainty on the tidal radius.

3 Results

We report here the main results concerning our fit of the projected velocity dispersion as function of stellar mass. The estimated model parameters are given and the relation between Φ₀ and the equipartition mass m_eq is shown and discussed. The relations with GCs structural parameters are also presented.

Finally, we show the estimated Φ₀ parameter by fitting the SBPs and compare them with the values obtained by fitting the σ(m) observable.

3.1 Fitting the velocity dispersion to constrain Φ₀

As a first approach, we evaluate Φ₀ by fitting the velocity dispersion data with the model prediction, considering a single power law mass function ξ(m) ∝ m^α with the slope α taken from Baumgardt et al. (2023), but with a mass range 0.1 < m/M_⊙ < 1.0, slightly broader than the usual one [0.2, 0.8] M_⊙. However, it can be shown that the effect of a different range of masses (when reasonable) is negligible with respect to variations in Φ₀ or α.

In Fig. 1 we plot σ(m) predicted by the model against dispersion data by W22 and the fit with the Bianchini function for NGC 6397.

The plot clearly underlines a great compatibility between the model prediction and the Bianchini fitting function on the same data. We confirm a partial degree of equipartition in all clusters and for all masses (see also Fig. B.1 for other clusters). We expect that a cluster in an advanced dynamical state is closer to energy equipartition and shows a smaller value of m_eq . In the model framework, dynamically old means a greater Φ₀ value, namely a deeper gravitational potential well. Then, there should be a relation between these parameters: a cluster showing a higher degree of equipartition should have a greater value of Φ₀ and a lower value of m_eq . This is what we obtain in Fig. 2, where we plot the relationship between m_eq by W22 and the estimated Φ₀, colored with the obtained χ² . In Table 1 we give the estimated parameters and the χ² test value for all analyzed clusters, also for our fit with the Bianchini fitting function. Our evaluation of m_eq , conducted on binned data and with a nonlinear least squares algorithm (available in the SciPy software package), is similar to that obtained by W22, who used unbinned data, set prior constraints on the parameter values, and used a Monte Carlo method to sample the parameter space. The χ² values of the two fits are close to each other, sustaining a similar confidence with the data. Furthermore, the scaling dispersion σ0 of our model is systematically lower than σ_0,B , a natural behavior due to the different meaning they have. While the first is the dispersion of the lowest model mass, the latter is the dispersion value in the limit m → 0.

The results of the fits made by using the classical power law σ(m) = σ_s(m/m_s)^−η, with m_s = 1.0 M_⊙, are shown in Table 2.

Our estimates of η, obtained with the same method as for m_eq , are compatible with those provided by W22, statistically equivalent to the other fitting functions, with small variations in the obtained χ² . Furthermore, η follows the expected trend with m_eq and Φ₀ . For a more advanced dynamical state, we have a larger degree of equipartition described by a higher η value, as well as a larger Φ₀ and a smaller m_eq , as expected. W22 already provided a discussion of the relationship between the classical power law and the Bianchini fitting function, as well as the obtained values for the equipartition mass and η, showing that they appear equivalent in measuring the equipartition degree. Here we focus on the link between the Bianchini function and the dynamical model prediction due to the closer statistical compatibility obtained in our fits.

Concerning the results of the fitting procedure and the relation between m_eq and Φ₀ , some considerations must be made. First, for NGC 5139 (ω-Cen) we have the lowest value of Φ₀ and the greatest value of m_eq, also coming with a very poor fit quality (largest χ² value). This cluster must be treated with caution due to its intricate history and structure. It shows evidence of multiple stellar populations that are dynamically interacting and are likely to have different degrees of equipartition. Furthermore, each has its own mass function and chemical pattern. It is also strongly suspected to be the survived core of a disrupted dwarf galaxy. The full picture behind this object requires further study with more sophisticated models, given its complexity. Although the overall trend for σ(m) looks promising, we emphasize that the uncertainties may be underestimated due to the fluctuations around the decreasing pattern (we come back to this point in the next lines). This suggests an additional source of statistical error that may increase the estimated uncertainty in the corresponding Φ₀ value. However, the long relaxation timescale of such a cluster (see Table 1 in W22 work) may indicate that our estimates are not far from reality, specifically a low energy equipartition level corresponding to a smaller Φ₀ and a larger m_eq . Furthermore, NGC 6397 and NGC 6752 are post-core collapsed objects according to the Harris catalog, since they have a very high value of the concentration, while NGC 6266 is suspected to be collapsed even if its concentration is not very high. Considering that only for three of the analyzed clusters, namely NGC 6266, NGC 6341 and NGC 6397, we have a high confidence level coming from the fitting procedure, excluding collapsed objects and poor fit quality results, we are left only with NGC 6341, which is also the cluster with the lowest m_eq and the highest Φ₀ . Moreover, this cluster has a low mass coverage in the velocity dispersion data, with a range of 0.59 < m/M_⊙ < 0.78. Looking at data points (see Fig. B.1), it looks evident a decreasing trend of σ(m) with the mass, but the oscillations around such tendency, which change from cluster to cluster, suggest that uncertainties can be underestimated. Indeed, we obtain a relation between our χ² test value and the average relative error on σ(m) data (see Fig. B.4), showing that low relative errors produce a greater χ² value in our fitting procedure (with ω-Cen showing the lower relative error and the larger χ²), while the opposite is vaguely true for greater relative errors. Furthermore, our χ² looks slightly asymmetric with respect to Φ₀ in the minimum, suggesting a more advanced uncertainty estimation method for such parameter, like evaluating the right and left errors. However, this does not affect the results that we present here.

The uncertainties given for Φ₀ reflect the observed ones and come from the minimization process. A possible source of error in our procedure is the choice of the core radius. It can affect the theoretical prediction, since it determines the radial shell where the projected velocity dispersion σ(R, m) is averaged (see Sect. 2.2.1). An underestimation of the core radius would map the observed radial shell into a more outer dimensionless radial shell. As a result, the velocity dispersion would be generally smaller and more dominated by the low-mass stars due to mass segregation. Then, σ(m) would be steeper and the corresponding best-fit Φ₀ would be higher. The opposite holds for an overestimate of the core radius. This introduces a systematic error into our procedure for determining model parameters, but we do not explore such an effect here. Actually, the expected relation between Φ₀ and r_c suggests that only a large relative uncertainty on the core radius may affect the estimate of Φ₀ . We quantify the effect of varying r_c on the Φ₀ determination in Sect. 3.4.

All these uncertainties motivate an independent examination of Φ₀ , as obtained by fitting the surface brightness profile (see Sect. 3.5).

Consequently, our work must be considered as a first guess in the evaluation of multi-mass dynamical model parameters from GCs kinematic observations concerning energy equipartition.

Fig. 1 Projected velocity dispersion as a function of stellar mass for NGC 6397. The error bars are the data from Watkins et al. (2022), the continuous green line is our best fit with the Bianchini et al. (2016) fitting function in Eq. (1) with its error band, the dotted brown line is the complete equipartition limit (σ ∝ m−1/2), and the dashed orange line is our model best-fit with its 68% confidence band.

Fig. 2 Relation between the estimated Φ0 and meq from Watkins et al. (2022). The χ2 test value is given by the color of the data point for each NGC according to the color bar to the right.

3.2 Adding the mass function slope to the parameter space

As an extension of the previous analysis, we add to the parameter space the mass function slope α. The fitting procedure is naturally more computationally expensive, because we have to iteratively solve equilibrium configurations with a different mass function. As previously done, we minimize the χ² for each cluster: we explore Φ₀ around a value close to the one obtained with the first approach and the slope α is taken in the range [−2.0, 0.0]. The adopted range for α accounts for the typical slopes values of GCs and is only a first guess, to be eventually modified to increase the precision. Unfortunately, we find a degeneration between Φ₀ and α, preventing a more precise evaluation of the slope.

In Fig. 3, we show the contour plot for NGC 6397 in the Φ₀–α plane, with some values around the χ² minimum, together with the mass function assumed in the first approach, with its uncertainty. Even if our minimum falls in the error band by Baumgardt et al. (2023) for NGC 6397, the same is not valid for other clusters (see Fig. B.2). Furthermore, the shape of the parameter space is preventing us to constrain the slope of the mass function, whose uncertainties are too big and overcome the considered range. In Table 3 we give the results for the estimated model parameters, with their errors.

Although the contour plots are not the same for all analyzed clusters, in some cases there is a small tendency of having a greater Φ₀ corresponding to a lower α, while in other ones the contour plot lines are almost horizontally oriented. Indeed, the velocity dispersion dependence on mass is steeper and closer to equipartition for both a larger Φ₀ or a greater slope α. Such effect is expected since a larger slope means a flatter mass function, with more massive stars than a steeper one. This results in a greater degree of energy equipartition, because massive stars are more efficient in reaching equipartition due to their higher mass, which cause them to suffer more gravitational encounters than less massive stars. Similarly, a greater Φ₀ implies a deeper potential well, namely a more dynamically evolved cluster, with a higher level of equipartition. As a result, a given degree of equipartition can be reached with more massive stars and a less deep gravitational well or by a stronger gravity and a steeper mass function. However, the effect of the gravitational potential through the Φ₀ parameter looks more important, with smaller changes affecting equipartition more than variations in the slope. This is also remarked by the estimated errors on Φ₀ , which are compatible with respect to the first approach, meaning that any reasonable assumption on the mass function slope would work well in setting up the model and constrain Φ₀ from the velocity dispersion σ(m).

Fig. 3 Contour plot for NGC 6397 in the Φ0−α plane, showing the χ2 levels. The blue band gives the slope from Baumgardt et al. (2023). Squares indicate the χ2 values around the minimum, which is identified by a large gray circle.

Fig. 4 Radial and projected profile of the equipartition mass meq, obtained by fitting the model prediction of the velocity dispersion with the Bianchini function at each radial coordinate, given in units of the core radius. The panel zooms onto the region r ≤ rc.

3.3 Shell selection and projection effects

As already outlined, the data we are fitting concerns projected quantities, and it is restricted to a circular ring close to the central regions of the clusters. We explore, by means of our model and the Bianchini fitting function, the effects of projection and shell selection in quantifying the degree of equipartition.

Our dynamical model predicts the shape of the function σ(r, m), which can be used to explore the degree of energy equipartition in different radial shells. Equipartition is more efficient in the core of clusters, and the maximum degree is expected to be measured there. Unfortunately, to date kinematic measurements in the core are difficult due to crowding effects and are limiting observers to look around it, as W22 do. As a result, the measured degree of equipartition is less than the core one. Furthermore, working with projected velocities also gives a smaller equipartition level than using 3D quantities. Indeed, at the projected radial distance R, the observer intercepts stars with a 3D radial distance r є [R, r_t]. Being the velocity dispersion σ(r, m) a decreasing function of r, its value is lowered for each mass, in average, with respect to the 3D one. Being massive stars more centrally concentrated than less massive ones, the projection decreases more σ(m) for low-mass stars and brings to a flatter profile, resulting in a lower equipartition degree with respect to the 3D case. Both these effects bias the estimation of m_eq, leading to higher values (i.e., lower equipartition). To visualize this, we plot in Fig. 4 the value of m_eq obtained by fitting with the Bianchini function the 3D and 2D theoretical profiles σ(r, m) and σ(R, m) respectively, for a particular choice of α and Φ₀. Such analysis is possible due to the very good similarity between our prediction for σ(m) and Eq. (1), which we check quantifying the uncertainty of the fitted m_eq and verifying its small value.

As expected, m_eq increases with the distance from the center of the cluster, due to the decreasing efficiency of the relaxation process. A measure of the maximum degree of equipartition should be done with the 3D equipartition mass at least in the core r < r_c, where an almost flat trend is seen, as shown in the zoom panel of Fig. 4. The projected profile of m_eq always gives bigger values, even in the core. Every analysis of the projected velocity dispersion as function of stellar mass should consider that the estimation of equipartition obtained fitting data with the Bian- chini function will give an overestimated m_eq (underestimating the degree of equipartition), with respect to the 3D value. This is not the case when the fit is done with our model, since the degree of equipartition depends on Φ₀ that uniquely determines the velocity dispersion value σ(r, m) in any position and for each mass, from which the projected one can be computed without any loss of information.

Furthermore, we give in Fig. 5 few 3D radial profiles of m_eq obtained changing Φ₀ and α, to show their effect. Here, for very low values of Φ₀, namely a dynamically young cluster, the profile m_eq(r) starts from ∼10 M_⊙ in the central regions and grows with radius, almost insensible to α variations. With a greater Φ₀ , the effect of a different slope is clearly visible in the outer regions. A flatter mass function will produce a more extended profile, as a greater Φ₀ do (due to the existent degeneration). Additionally, it also gives smaller values of m_eq at each radial coordinate than a steeper mass function, as already mentioned.

In the core, the value of m_eq is approximately constant and decreases with increasing Φ₀ , as expected from more advanced dynamical states. The zoom panel of Fig. 5 gives the trend for r ≤ r_c , showing the small differences between the considered slopes. Such approximately flatness of m_eq(r) in the central region, a common behavior of other observable profiles like the surface density, can be used to quantify the maximum degree of equipartition reached by stars in the cluster. Then, we build a relation between m_eq(r ≤ r_c) and Φ₀, which is highly useful in the perspective of quantifying the maximum degree of equipartition in terms of the equipartition mass. Indeed, the value of Φ₀ can be constrained from other observables more easily measurable than the velocity dispersion as function of mass in the center, which suffers crowding. Since the slope of the mass function has a negligible effect in the relation between m_eq and Φ₀ in the core, we choose α = −0.5 and plot such relationship in Fig. 6, comparing with the estimates of m_eq and Φ₀ already shown in Fig. 2 and given in Table 1. The plot shows how the parameters constrained from the projected σ(m) observations lie in the region above the curve. This underlines that the estimated m_eq by W22 do not measure the maximum degree of equipartition. They suffer the projection effect and the radial shell selection, as already said. Indeed, the radial coverage of the data sample has a median radius beyond the core radius for most clusters. Only for two clusters, namely NGC 5139 and NGC 6656, the selected shell lies around the core radius (0.4 ≲ r/r_c ≲ 1.2). Indeed, these clusters are also the closest points to the theoretical curve in Fig. 6.

Concerning the estimates of Φ₀ , we stress that this is not affected by projection or radial coverage, since it is a global parameter that identifies the equilibrium configuration with its profiles. As a result, once Φ₀ is determined or constrained from observations, even with a limited radial coverage, the equipartition degree is known in every region of the cluster. It can also be presented in terms of the Bianchini equipartition mass, as we have shown here.

Fig. 5 Theoretical radial profile of the equipartition mass obtained for different values of the model parameters Φ0 (same line style) and α(same colors and line width). The panel zooms onto the core region r ≤ rc.

Fig. 6 Equipartition mass meq in the core (r ≤ rc) as function of Φ0. The blue continuous line and the orange dashed line are the 3D and 2D theoretical predictions, respectively. The circles with error bars are the same as in Fig. 2.

3.4 Relation with structural parameters

The dynamical state of GCs is related to the amount of relaxation they experienced during their life. A dynamically old cluster has seen several relaxation timescales t_relax during its evolution, a time that is shorter in the central regions, being gravitational interactions more frequent. In literature, the relaxation process is often quantified with the relaxation timescale in the core t_rc (Djorgovski 1993) and at half-mass radius (Binney & Tremaine 2008), which have an analytical approximated formula for their evaluation. This leads W22 to compare their prediction of m_eq with the number of core and median relaxation timescales N_core and N_half, as shown in their Fig. 16.

From the King (1966) model, it is also expected that dynamically evolved clusters have a greater concentration c = log(r_t/r_k) value, with the King radius r_k often replaced with the core radius r_c . Indeed, a monotonically increasing relation exists between the concentration and the corresponding W₀ equilibrium parameter of the King model.

In the framework of our multi-mass dynamical model, the equilibrium configuration is identified by Φ₀ . However, the scaling radius of the model r_k,u loses, in principle, the relation with the core radius that the King model had. This leads us to define the equivalent King concentration c_k = log(r_t/ r_k,u) and the theoretical concentration c_th = log(r_t/r_c,th), with r_c,th estimated from the surface density profile, assuming a proportionality relationship between the mass and luminosity density profiles (as outlined before). These concentrations have a small relative difference and both monotonically increase with Φ₀, following a trend that depends on the slope α. In the perspective of comparing with observations, we consider only c_th.

About the relaxation timescale, the model derivation highlights a dependence on mass, consistent with what is observed. However, obtaining an analytical formula requires further investigation and development. We are then limited to take the values of N_core and N_half from literature.

We expect that Φ₀ increases with N_core, N_half and c_obs = log(r_t/r_c). For the core and median relaxation timescales numbers we take the same values considered by W22, while for the concentrations we consider c in the Harris catalog and its estimation obtained with two sources for the core radius, the one in the Harris catalog and the other by Baumgardt (2018–2023) online catalog (which used the Spitzer 1987 definition), while for the tidal radius we refer to Webb et al. (2013), reported in the same catalog.

In Fig. 7, the previously estimated values of Φ₀ are given with the concentration. We show the theoretical relation for some choices of the slope α, against the observational concentrations from different sources. The predicted concentration grows continuously with Φ₀, as already stated, but there are differences with the observational ones. From the model point of view, such discrepancy can be related to the fit quality and the value of the theoretical core radius. Low fit quality results can give an incorrect estimate of Φ₀, which affects the vertical positions of each marker in the pictures. Concerning the core radius, the adopted choice is reasonable, since it is based on the assumption that the luminosity density profile is proportional to the density profile, with a multiplying factor only dependent on mass. From the observational point of view, we first note that the Harris catalog concentrations come by fitting the surface brightness profile of clusters with a King single-mass model (King 1966), which has been found to be in disagreement with observations in the outer regions, being clusters more extended, like already pointed out so far by Da Costa & Freeman (1976). Indeed, the authors found that the outer regions of M3 can be described by a King profile with c = 1.98, while the inner ones using c = 1.29. This suggests that clusters are more extended and, consequently, have a higher concentration with respect to the single-mass model prediction. Multi-mass models can overcome such discrepancy, being more extended. In principle, a concentration value completely based on observations requires the knowledge of both the core radius and the tidal radius. While the first is easier to measure, the latter is quite tricky, since a clear distinction between stars that belong to the cluster or are leaving it is difficult to infer. This leads us to compute the concentration using the core and tidal radius from different sources (Baumgardt 2018–2023; Harris 1996; Webb et al. 2013). Indeed, Webb et al. (2013) found an empirical formula for the limiting radius of GCs through N-body simulations, which takes the orbital properties of the system in the galactic gravitational potential into account. This results in a greater tidal radius and consequently a greater concentration with respect to the King ones from the Harris catalog, as can be inferred from Fig. 7, comparing the orange squares and the green stars with the blue triangles. Such tendency slightly erases the differences with the theoretical prediction, with the only exception of NGC 5139 (ω-Cen), which however is a complex object, as already mentioned, whose deviation is not so surprising.

Our analysis of few clusters reveals that there are important inhomogeneities among observational-based concentrations. Indeed, we can state that the estimate of the concentration suffers several sources of errors that can differ between clusters. The difficulty of measuring the tidal radius and, consequently, the concentration suggests to recover the role of theoretical models in predicting its value, as it was for the King model. The concentration value could be estimated from the relationship with the gravitational potential well through Φ₀, to be constrained fitting observational profiles with model predictions. Anyway, an extensive comparison between the theoretically predicted concentrations and the observational ones must be addressed in future works, with a greater GCs sample and heterogeneous observational sources.

In Fig. 8, we show the relation between the gravitational potential and the number of core and median relaxation timescales. The plot shows a tendency of having greater Φ₀ values with a larger N_core and N_half , as expected for more dynamically evolved states. A similar trend was found by W22 with the m_eq parameter. Unfortunately, the great errors and the absence of a theoretical prediction on N_core and N_half do not allow us to further discuss the obtained pattern.

Another quantifier of the dynamical state of GCs is the so-called dynamical-clock, mainly measured through the A⁺ parameter that is the area between the cumulative distribution of Blue Straggler Stars (BSSs) and the reference stars (Ferraro et al. 2012, 2020; Alessandrini et al. 2016; Lanzoni et al. 2016). Only four of the analyzed GCs have a measure of A⁺, which also well correlates with the number of core relaxation timescales and the core radius, as shown by Lanzoni et al. (2016) and Ferraro et al. (2018). We give in Fig. 9 the estimated Φ₀ and the corresponding A⁺ and log(r_c), where we use again the core radius from the Harris catalog and Baumgardt (2018–2023) online catalog.

As for the number of core relaxation timescales, the Φ₀ parameter grows at increasing values of A⁺, both measuring the dynamical state of GCs, although the number of points is small and only a general increasing trend can be deduced. While Φ₀ measures the gravitational potential well, the A⁺ parameter tracks the segregation of BSSs with respect to the other stars, measuring the mass segregation process. They have a lower velocity dispersion value, corresponding to their expected higher mass with respect to reference stars, namely the main sequence turn-off mass (see Baldwin et al. 2016, and references therein). In principle, the mass of BSSs depends on the formation mechanisms of such objects, which is still unknown. The most spread idea is that they are binary stars or a result of stellar coalescence. To precisely measure the dynamical state of the cluster, the formation of such objects should be mostly coeval with the cluster itself. Any recently formed BSS will naturally be less segregated, meaning that the overall BSSs population segregation will decrease, resulting in a lower dynamical state for the cluster. Indeed, the dynamical processes that bring the collisional system toward mass segregation and energy equipartition are phenomena that strongly depend on the mass of the object and the relation between its age and the relaxation timescale. Our work, as well as several energy equipartition studies in literature, underlines that massive stars are faster in reaching equipartition than less massive ones. Consequently, the population of BSSs can have internal differences in their kinematic properties if a spread in mass exists. As we observe and predict a different degree of energy equipartition for reference stars, a similar phenomenology can state also for BSSs. One of the main strength of our modeling approach is the ability of predicting such kinematic differences for objects with whatever mass, once known their mass function and the structural parameter Φ₀. That would further broaden the role of BSSs in characterizing GCs dynamical state and understanding the complex internal dynamics, giving important insight in the energy equipartition and mass segregation topics.

Concerning the relation with the core radius, advanced dynamical states have a greater value of Φ₀ and a lower r_c, as shown in the lower panel of Fig. 9, which also affects the concentration that will be higher. Indeed, as the relaxation proceed, the cluster is driven toward the core collapse phase and the gravothermal catastrophe. Here, the structure of the system is altered and King-based models do not reproduce well the structural properties. However, the energy equipartition keeps going on in the survived core, where a higher degree of equipartition can be expected. As mentioned in Sect. 3.1, the estimates of Φ₀ from the observable σ(m) may suffer from a systematic effect related to the adopted value for the core radius, which is neglected in our analysis. With simple order-of-magnitude calculations and a linear interpolation of the Φ₀−log (r_c) relation, in the case of NGC 5904 and NGC 6397 a 10% variation in the core radius is requested to get a systematic shift on Φ₀ comparable to the estimated statistical error. For all the other clusters in our sample, the variation of rc should be much larger. Although the relative error for the core radius is typically smaller, as can be seen from Table 4, larger discrepancies are found when the Harris catalog values are compared to those of Baumgardt (2018–2023).

Fig. 7 Relation between Φ0 and the concentration c = log(rt/rc). The continuous, dashed, and dotted lines represent, respectively, the theoretical prediction with a mass function slope α = 0.0, α = −1.0, and α = −2.0, while the black circles show the constrained values of Φ0 and c for the analyzed clusters. The blue triangles show the King concentration from the Harris (1996) catalog (2010 edition), while the orange squares and green stars are computed using the core radius from the Harris (1996) and Baumgardt (2018–2023) catalogs respectively, and the tidal radius from Webb et al. (2013).

Fig. 8 Obtained values of Φ0 against the number of core relaxation timescales Ncore (upper panel) and the number of median relaxation timescales Nhalf (lower panel) from Watkins et al. (2022).

Fig. 9 Comparison of structural parameters that measure the dynamical state of globular clusters. Upper panel: relation between Φ0 and the area A+ between the cumulative distribution of blue straggler stars and the reference stars. Lower panel: Φ0 relation with the core radius, using values from the Harris (1996) catalog (blue triangles) and the Baumgardt (2018–2023) online catalog (orange squares).

3.5 Fitting surface brightness profiles

To provide an independent determination of the structural parameter Φ₀, we fit the SBPs data by Trager et al. (1995) with the prediction of the dynamical model, by using theoretical isochrones from BaSTI (Hidalgo et al. 2018; Pietrinferni et al. 2021; Salaris et al. 2022; Pietrinferni et al. 2024) for an age of 13 Gyrs, with [α/Fe] = +0.4, Y = 0.247 and a different metal- licity [Fe/H] for each cluster, taken from the Harris catalog, as described in Sect. 2.2.2.

In Fig. 10, we plot the obtained best-fit profile for NGC 6341. In Table 4, we give the estimated Φ₀ as well as the central surface brightness, the core and tidal radius and the normalized χ² value for the analyzed clusters.

The obtained best-fit curves reproduce well the observed data, with a very narrow error band and small uncertainties for the parameters. However, the obtained estimates neglect possible error sources, both statistical and systematic ones. Regarding the conversion of the observed magnitudes into luminosities, several systematic effects can occur, as due to the extinction values or the distance modulus value as well as in the assumed solar V- magnitude. However, all these parameters produce an offset in the magnitude, but they do not affect our minimization procedure and the obtained Φ₀. This is not the case for the parameters related to the isochrones, such as the age and the adopted values of [α/Fe], He abundance and [Fe/H]. As an example, we report in the last column of Table 4 the effect of a different age in the estimated value of Φ₀, using isochrones at 11 Gyr, keeping the previous adopted values for [α/Fe], Y and [Fe/H]. The overall impact of a reduced age is a smaller Φ₀ , resulting in a systematic effect. As already explained in Sect. 2.2.2, a smaller age implies a larger maximum mass. Since more massive stars are more segregated and affect the luminosity distribution in the innermost region of the cluster, a steeper brightness profile is expected. Therefore, in order to reproduce the observed SBP, the effect of an increase in the maximum mass should be compensated by a reduction of the Φ₀ value. We note that we cannot distinguish from the 11 Gyr and 13 Gyr cases, because the implied variations of the observable parameters (such as the central surface brightness, core radius and tidal radius) are smaller than their respective statistical errors. In addition, the isochrones vary from cluster to cluster also due to the different metallic- ity that we considered. The role of the metallicity, [α/Fe], He mass fraction Y, and cluster age in shaping GCs SBPs, especially in the framework of multi-mass King-like models, requires a broader exploration and a dedicated work. The same holds for the mass function shape, which can alter the predictions due to its structural role, being a fundamental physical ingredient in such dynamical models. In fact, the massive stars content is extremely important in shaping the SBPs and thus the slope α is expected to play a relevant effect.

Finally, we compare the independent estimates of Φ₀ obtained from internal kinematics, namely the velocity dispersion as function of stellar mass, with those from the SBPs. Concerning the latter, we first account for age variations by averaging the corresponding Φ₀ values and computing the associated uncertainty. In Fig. 11, we plot the relative variation ΔΦ₀/〈Φ₀〉, showing the 2σ confidence interval. The two estimates are compatible at 2σ level for most clusters, although some are compatible at 1σ confidence interval (namely NGC 5904, NGC 6266 and NGC 6656). The larger deviation is obtained for NGC 5139 (ω-Cen). As already mentioned, this cluster is a much more complex object. It is probably the remnant of a dwarf galaxy. It hosts multiple stellar generations, also showing metallicity variations, for which multiple isochrones should be used and potentially multiple dynamical states for each of the subsystems that make up such an object.

The plot also highlights how multi-mass King-like dynamical models can predict both internal kinematics and the surface brightness profiles. However, as already outlined, there are several possible sources of error that could further increase our uncertainties in the Φ₀ estimates. A full exploration of the assumptions made in the SBPs fitting procedure is needed to better constrain the Φ₀ value. Only afterward, its relation with the other structural quantities such as concentration, number of relaxation times, core radius and A⁺ can be addressed similarly to what done in Sect. 3.4.

Fig. 10 Surface brightness profile for NGC 6341. The black circles with error bars are the data from Trager et al. (1995) analyzed following the work by McLaughlin & van der Marel (2005) and Zocchi et al. (2012). The continuous blue line is our model best fit with its confidence band, which is obtained by assuming the Baumgardt et al. (2023) mass function slope and adopting the BaSTI isochrones (Hidalgo et al. 2018; Pietrinferni et al. 2021; Salaris et al. 2022; Pietrinferni et al. 2024) with 13 Gyr, [α/Fe] = +0.4, Y = 0.247, and metallicity [Fe/H]=−2.31, taken from the Harris (1996) catalog (2010 edition).

Fig. 11 Relative variation in the estimates of the Φ0 parameter for each of the analyzed clusters, obtained by fitting the velocity dispersion dependence on stellar mass σ(m) and the surface brightness profiles I(R). The error bars are the 2σ level uncertainties.

4 Conclusions

In this work, we explore the degree of energy equipartition in GCs by means of a multi-mass King-like dynamical model for collisional and phase-space-limited systems. This theoretical description offers an efficient way to predict several dynamical processes occurring in GCs, such as energy equipartition, mass segregation, and evaporation. These phenomena alter the properties of clusters, modifying their spatial and kinematic structure. Here we focus on the degree of equipartition, and determine the dynamical state of a few GCs by means of model parameters.

We fit the velocity dispersion dependence on the stellar mass – measured through HST proper motions by Watkins et al. (2022) – with our prediction, discussing the relation with the Bianchini et al. (2016) fitting function and the equipartition mass m_eq. We obtain a similar confidence level and a relation between the equipartition mass and the model parameter Φ₀ , which is a measure of the gravitational potential difference between the center and the edge of the cluster. An increasing value for this quantity describes a more dynamically old cluster, with a larger degree of energy equipartition, which is in turn associated with a lower value of m_eq.

We confirm that the energy equipartition in the analyzed GCs is only partial, even in the central regions. The equipartition degree appears similar in the core, while decreasing significantly in the outer regions. More massive stars are closer to equipartition than less massive ones, as already found and outlined in other works.

We extend our fitting procedure by adding the slope of the mass function α to the parameter space; α is a required input for the model, and we take its value from Baumgardt et al. (2023). This procedure provides similar results regarding the estimate of Φ₀ , but prevents any constraint on the mass function slope using the data. Although we obtain a degeneration between the two parameters, the effect of varying the gravitational potential well alters the observable σ(m) to a greater degree than the variations in the mass function when reasonable (i.e., α є [−2.0, 0]).

We carefully discuss the implications of quantifying the degree of energy equipartition through the Bianchini fitting function with respect to our dynamical model prediction, when working with a restricted radial shell and projected quantities. We find that the estimated m_eq from the projected dispersion is higher than the three-dimensional one. Moreover, for radial shells overcoming the core radius, the equipartition mass underestimates the maximum degree of equipartition reached in the core. On the contrary, the equilibrium parameter Φ₀ uniquely defines both projected and 3D radial theoretical profiles, without suffering shell selection or projection effects. We also find that variations in the slope α affect the degree of equipartition mainly in the outer regions, while their effect is relatively small in the core. Here, by fitting theoretical profiles with the Bianchini function, we obtain a mostly constant m_eq, suggesting that this value can be used as a measure of the maximum degree of energy equipartition in clusters. However, taking advantages of dynamical models like ours offers the opportunity to quantify the degree of equipartition in the cores of clusters by means of structural parameters. Here, observational data suffer several limitations and a theoretical tool for predicting observables can offer important support to the astronomical community.

We compare the estimates of Φ₀ with other structural properties of GCs that depend on the dynamical state of the system. The theoretical relation between Φ₀ and the concentration is presented, as is the relation between the estimated Φ₀ and observational values for the concentration from different sources (Harris 1996; Baumgardt 2018–2023; Webb et al. 2013). The spread among observations as well as the differences with respect to the theoretical prediction highlight the presence of several sources of error. However, the theoretical relation between Φ₀ and the concentration provides an important tool to constrain the latter on a more fundamental ground, once the former is properly estimated. Consequently, we advise to recover the predicting role of models concerning the determination of the concentration value for GCs. A positive trend is also seen between Φ₀ and the number of core and median relaxation timescales. We note a similar trend for the area A⁺ between the cumulative distributions of blue straggler stars and reference stars (Lanzoni et al. 2016; Ferraro et al. 2018), although the number of points is small (only four clusters). All such structural quantities, such as Φ₀ , increase for dynamically old systems. Conversely, the core radius decreases with the dynamical state, a trend we clearly see between Φ₀ and r_c. These outcomes strongly suggest that we can include Φ₀ among the properties that track the dynamical age of GCs. However, a statistically robust comparison between the structural parameters of Galactic GCs and those of our model must be carried out in future works, with a wider sample of clusters.

Finally, we successfully fit the surface brightness profiles by Trager et al. (1995) following a procedure similar to that described by Trenti & van der Marel (2013) and Zocchi et al. (2012). To make the theoretical prediction, we need the surface density profile and a mass–luminosity relation, which we take from theoretical isochrones (Hidalgo et al. 2018; Pietrinferni et al. 2021; Salaris et al. 2022; Pietrinferni et al. 2024). We obtain an estimate for Φ₀ , the central surface brightness, the core, and the tidal radius. We discuss the possible effects that can alter our predictions, mainly in regard to the mass–luminosity relation. A different age for the isochrone introduces a systematic error in the determination of Φ₀ that is larger than the statistical uncertainty inherent to the fitting procedure. We average these contributions to get Φ₀ estimates from the fitting procedure on each cluster SBP. We compare these values with the ones computed by fitting the velocity dispersion–mass relation, showing that they are compatible at the 2σ level for almost all clusters. The largest deviation is obtained for ω-Cen, whose complexity suggests that our fitting procedure needs to be extended to consider multiple stellar populations with different chemical contents and metallicities, and consequently different isochrones.

The obtained results strongly underline the central role of dynamical models in predicting the phenomenology of GCs, such as the highly debated energy equipartition process. With the increasing amount of information coming from internal kinematics observations, more effort must be put into advanced physical modeling. Such models are powerful tools for elucidating the mechanisms behind several dynamical phenomena, such as segregation and evaporation, as well as anisotropic velocity distributions and internal rotation. However, determinations of the dynamical state and predictions of the different observables are affected by several sources of error. These errors are related to, for example, the internal kinematics and the surface brightness profiles, where the role of the assumptions made regarding structural properties and the mass–luminosity relation must be further explored.

Data availability

Appendix B is available on Zenodo.

Acknowledgements

MT thanks Michele Bellanzini for important hints that expanded this work. The authors would thank the anonymous referee for the useful indications and comments that improved this paper.

Appendix A Model derivation

The distribution function (DF) that defines our model is given in Eq. (2). It comes from an analytical derivation, which replicates and adapt the analytical procedures that brought from the dynamical friction (Chandrasekhar 1943) to the King (1966) model. In King (1965), the author obtained a DF as an approximated solution of the Fokker-Planck equation, valid for collisional systems. Then, he generalized and used it in King (1966). Our DF keeps the memory of the King one, but takes the mass distribution into account, as was done by Da Costa & Freeman (1976). They presented a discrete model, where each mass class has a DF with an energy cut-off and a weight factor, to be constrained from observations. In principle, the assumption that each mass has that specific DF analytical shape must be proved similarly, solving the Fokker-Planck equation in a multi-mass system.

The procedure that lead to the King DF starts from the Boltzmann equation for collisional systems and follows the Chandrasekhar (1943) development: assuming local approximation and low energy exchanges, the equation is written in terms of the diffusion coefficients, which quantify the dynamical friction that a test star with mass m suffers due to field stars with mass m_a. Following Binney & Tremaine (2008), the Spitzer & Härm (1958) expression of the Fokker-Planck equation is obtained. King was able to solve such equation in 1965 with an approximation. He originally obtained a DF where both masses m and m_a appear, but then he takes them equal.

In principle, the same approach can be followed again, but taking care of considering the field stars’ mass function when evaluating the diffusion coefficients, as well as their velocity distribution, assumed Maxwell-Boltzmann like. A further needed step is the distinction between dynamical quantities and thermodynamic ones. In particular, when using a Maxwell-Boltzmann DF for field stars, which is $\propto \exp \left[-v^2/\left(2{\sigma }_{\text{a}}^2\right)\right]$, one gets the scaling factor σ_a which is related to the thermodynamic temperature θ (kept constant) and the field stars’ mass, namely ${\sigma }_a^2=k_{\text{B}}\theta /m_{\text{a}}$. In this framework σ is the 1D velocity dispersion of the Maxwell-Boltzmann DF, and it depends only on the mass, while there is another quantity that comes at the end of the procedure, namely the kinetic one. In the limit of a Maxwell-Boltzmann DF for the test star (that describes an isothermal sphere) the kinetic quantity converges to the thermodynamic one. This distinction allows us to explicit mass dependence in the obtained generalized expression of the Fokker-Planck equation, which finally states $$\frac{\text{d}f}{\hspace{0.17em}\text{d}t}=\frac{1}{t_{\text{R}}(m)}\frac{1}{x^2}\frac{\partial }{\partial x}\left[2x\mathcal{G}(x,m)\left(2xf+\frac{\partial f}{\partial x}\right)\right],$$(A.1)

where f = f (x, m; t) is the DF of the test star, x = v²/(2σ²) and $$\frac{\mathcal{G}(x,m)}{t_{\text{R}}(m)}={\displaystyle {\int }_{\text{Δ}{m}_{\text{a}}}\frac{1}{{\tau }_{\text{R},\text{a}}}}{\left(\frac{m}{m_{\text{a}}}\right)}^{3/2}\left[\frac{2}{\sqrt{\pi }{x}_{\text{a}}^2}{\displaystyle {\int }_0^{x_{\text{a}}}{y}^2}{\text{e}}^{-y^2}\hspace{0.17em}\text{d}y\right]\text{d}{m}_{\text{a}},$$(A.2)

with $x_{\text{a}}=v^2/\left(2{\sigma }_{\text{a}}^2\right)=m_{\text{a}}{v}^2/\left(2k_{\text{B}}\theta \right)\text{\hspace{0.17em}and\hspace{0.17em}}{\tau }_{\text{R},\text{a}}={\tau }_{\text{R}}\left(m,m_{\text{a}}\right)$ the relaxation time for binary encounters between m and m_a , with the latter distributed in the interval Δm_a. Here, t_R(m) is the relaxation time of the mass m, due to gravitational encounters with all the field stars.

Following King (1965), Eq. (A.1) can be solved similarly assuming f (x, m; t) = ց(x, m)exp[−λ(m) t/1_R(m)], with λ the evaporation rate of stars with mass m and ց(x, m) = A(m)ḡ(x, m), which brings the dependence on x and the dimensions in the multiplying factor A(m). With an expansion in power series of the evaporation rate, with boundary conditions ḡ(0) = 1, ḡ′(0) = 0 and ḡ(x_e) = 0 where x_e is the cut-off, we obtain the approximated solution for ց(x, m) in the central region, where the treatment is valid, that is $$g(x,m)=k(m)\left[\exp \left(-x^2\right)-\exp \left(-x_{\text{e}}^2\right)\right],$$(A.3)

with k(m) that gathers all the mass-dependent multiplying factors. From the reasonable assumption that field stars have a Maxwell-Boltzmann DF, to be recovered in the limit of an infinite escape velocity (i.e., the isothermal sphere), and since the DFs we are dealing with give the number of stars in the infinitesimal volume d³v d³ r dm, it results that $$k(m)={\left(\frac{m}{2\pi k_{\text{B}}\theta }\right)}^{3/2}\frac{\hspace{0.17em}\text{d}{n}_0(m)}{\text{d}m}{\left[1-\exp \left(-x_{\text{e}}^2\right)\right]}^{-1},$$(A.4)

where n₀(m) is the number density of stars with mass between m and m + dm in the central region and dn₀/dm = ξ₀(m)/ V with ξ₀(m) the mass function and V the volume of the central region.

The final step is the generalization to different radial regions. Applying the Jeans theorem leads to the DF in Eq. (2) that is given in terms of the kinetic energy and cut-off energy.

References

All Tables

All Figures

Fig. 1 Projected velocity dispersion as a function of stellar mass for NGC 6397. The error bars are the data from Watkins et al. (2022), the continuous green line is our best fit with the Bianchini et al. (2016) fitting function in Eq. (1) with its error band, the dotted brown line is the complete equipartition limit (σ ∝ m−1/2), and the dashed orange line is our model best-fit with its 68% confidence band.

In the text

	Fig. 2 Relation between the estimated Φ0 and meq from Watkins et al. (2022). The χ2 test value is given by the color of the data point for each NGC according to the color bar to the right.
In the text

	Fig. 3 Contour plot for NGC 6397 in the Φ0−α plane, showing the χ2 levels. The blue band gives the slope from Baumgardt et al. (2023). Squares indicate the χ2 values around the minimum, which is identified by a large gray circle.
In the text

	Fig. 4 Radial and projected profile of the equipartition mass meq, obtained by fitting the model prediction of the velocity dispersion with the Bianchini function at each radial coordinate, given in units of the core radius. The panel zooms onto the region r ≤ rc.
In the text

	Fig. 5 Theoretical radial profile of the equipartition mass obtained for different values of the model parameters Φ0 (same line style) and α(same colors and line width). The panel zooms onto the core region r ≤ rc.
In the text

	Fig. 6 Equipartition mass meq in the core (r ≤ rc) as function of Φ0. The blue continuous line and the orange dashed line are the 3D and 2D theoretical predictions, respectively. The circles with error bars are the same as in Fig. 2.
In the text

Fig. 7 Relation between Φ0 and the concentration c = log(rt/rc). The continuous, dashed, and dotted lines represent, respectively, the theoretical prediction with a mass function slope α = 0.0, α = −1.0, and α = −2.0, while the black circles show the constrained values of Φ0 and c for the analyzed clusters. The blue triangles show the King concentration from the Harris (1996) catalog (2010 edition), while the orange squares and green stars are computed using the core radius from the Harris (1996) and Baumgardt (2018–2023) catalogs respectively, and the tidal radius from Webb et al. (2013).

In the text

	Fig. 8 Obtained values of Φ0 against the number of core relaxation timescales Ncore (upper panel) and the number of median relaxation timescales Nhalf (lower panel) from Watkins et al. (2022).
In the text

Fig. 9 Comparison of structural parameters that measure the dynamical state of globular clusters. Upper panel: relation between Φ0 and the area A+ between the cumulative distribution of blue straggler stars and the reference stars. Lower panel: Φ0 relation with the core radius, using values from the Harris (1996) catalog (blue triangles) and the Baumgardt (2018–2023) online catalog (orange squares).

In the text

Fig. 10 Surface brightness profile for NGC 6341. The black circles with error bars are the data from Trager et al. (1995) analyzed following the work by McLaughlin & van der Marel (2005) and Zocchi et al. (2012). The continuous blue line is our model best fit with its confidence band, which is obtained by assuming the Baumgardt et al. (2023) mass function slope and adopting the BaSTI isochrones (Hidalgo et al. 2018; Pietrinferni et al. 2021; Salaris et al. 2022; Pietrinferni et al. 2024) with 13 Gyr, [α/Fe] = +0.4, Y = 0.247, and metallicity [Fe/H]=−2.31, taken from the Harris (1996) catalog (2010 edition).

In the text

	Fig. 11 Relative variation in the estimates of the Φ0 parameter for each of the analyzed clusters, obtained by fitting the velocity dispersion dependence on stellar mass σ(m) and the surface brightness profiles I(R). The error bars are the 2σ level uncertainties.
In the text