nProFit: a tool for dynamical models fitting

The surface brightness profiles (SBPs) of star clusters hold invaluable information on the dynamical state of clusters. The observed SBPs of star clusters, especially that of globular clusters, are in good agreement with the SBPs expected for isothermal spheres containing stars of reduced kinetic energies. However, the SBPs of configurations that satisfy these theoretical criteria cannot be uniquely expressed by analytical formulae, which had hindered the analysis of dynamical state of observed clusters in external galaxies. To counter this shortcoming, it has become a practice to use empirical fitting formulae that best represent the core and halo characteristics of theoretical models. We here present a general purpose code, named nProFit, that allows fitting of the surface brightness profiles of extragalactic star clusters to theoretical star clusters, defined by dynamical models of King (1966) and Wilson (1975). In addition, we also incorporated theoretical models that result in power-law surface brightness profiles represented by Elson et al. 1987. The code returns the basic size parameters such as core radius, half-light radius and tidal radius, as well as dynamically relevant parameters, such as the volume and surface density profiles, velocity dispersion profile, total mass and the binding energy for a user-fixed mass-to-light ratio. The usefulness of the code in the dynamical study of extragalactic clusters has been already illustrated in Cuevas-Otahola et al. 2020. The code, which is python-based at the user end, but makes calls to advanced routines in Pyraf and Fortran, is now available for public use. We provide example scripts and mock clusters in the installation package as guide to users.


INTRODUCTION
Structural parameters of star clusters serve as proxies to give insights on their dynamical evolutionary state.These parameters, namely, core radius, half-mass and half-light radius, tidal radius and concentration index, can be obtained by fitting theoretical intensity profiles to the surface brightness profiles (SBP) of clusters.For instance, the King (1962) is one of the most widelyused profiles to obtain the structural parameters of old clusters, such as globular clusters (GCs).In the pioneering work, King (1966) demonstrated that the observed form of the profiles of GCs belong to the family of the surface mass density profiles corresponding to self-gravitating isothermal spheres of lower kinetic energies.Wilson (1975) proposed a dynamical profile resembling the structure of a King (1966) profile, with a larger halo, in order to fit the SBP of elliptical galaxies.Years later, Elson et al. (1987) found that the SBPs of clusters in the Large Magellanic Cloud (LMC) do not have a noticeable break corresponding to tidal radius of King profiles.They found their profiles are better fit by power-law functions rather than King profiles.In these power-law profiles, the shapes of the extended haloes are characterised by γ, with γ =2 corresponding to the profiles of infinite mass isothermal models.Profiles with γ >2 are steeper and have finite masses.
Over the last two decades, the observed SBPs of extragalactic star clusters have been analysed in several studies to obtain structural parameters using theoretical profiles.The most frequently used tools in these studies are ishape (Larsen 1999) and galfit (Peng et al. 2010), both of which are available for public use.Both these tools fit two-dimensional empirical profiles to the sky-subtracted observed 2D images of the objects under study.Geometrical parameters such as position angle and ellipticities are fitted as well.These codes are mainly used to obtain the size parameters such as core radius, half-light radius (R h ) and tidal radius for assumed shape of the profile.This is the best one can hope to do in images where star clusters are only marginally resolved and the background subtraction errors do not permit an analysis of the shapes of profiles in their external parts.
One of the defining parameters of star clusters is their mass, which is determined either using photometric techniques, or using dynamical models.The photometric mass is routinely determined using the observed luminosity along with a value for the mass-to-light ratio appropriate to the population of stars in the clus-ter.The mass-to-light ratio is calculated in population synthesis models, and is a function of age and metallicity (Bruzual & Charlot 2003).On the other hand, the dynamical mass is based on the determination of motions of stars under the influence of the collective force of all its stars, and is defined as the Virial mass, M vir = ησ 2 p R h /G, where σ p is the velocity dispersion projected along the line of sight, and η is the Virial factor which depends on the shape of the profile (Gieles et al. 2010).The photometric and dynamical masses do not always agree for clusters for which both the measurements are available (see e.g.McLaughlin et al. 2008;Gieles et al. 2010).Almost a factor of 10 uncertainty in η is one of the sources of the disagreements between the photometric and dynamical masses.This uncertainty can be avoided by theoretically calculating the radial profiles of dispersion velocities for the model that fits the observed SBP.Such calculations can be carried out for profiles that have an underlying physical model as illustrated by Barmby et al. (2007) and McLaughlin et al. (2008), who characterized the SBPs of globular clusters in M31 and NGC5128, respectively.The recent availability of high resolution multi-object or integral field unit-fed spectrographs on large telescopes (e.g.Gil de Paz et al. 2018) make it possible the determination of σ p of large samples of star clusters in nearby galaxies, which calls for the recovery of the σ p corresponding to the best-fit models.
The Hubble Space Telescope (HST) images of nearby galaxies (distance 5 Mpc) contain star clusters whose profiles are good enough for a characterisation of the shape of the outer halo, using physical models.However, the absence of a publicly available code is a handicap to analyse the profiles of star clusters in these images.The purpose of the present work is to develop a user-friendly code that can analyse the SBPs of star clusters on the HST images to obtain simultaneously the core and halo parameters, in addition to σ p .In order to achieve this, we follow the procedure outlined by Elson et al. (1987), as well as the prescription to fit dynamical models (King and Wilson) adopted by McLaughlin (2000), McLaughlin & van der Marel (2005) and Sollima et al. (2015).
We here introduce nProFit (Profile Fitting tool of nobjects) 1 for fitting dynamical models to 1D SBPs for a user-given list of star clusters in a single image.The dynamical models considered are King (King 1966), Wilson (Wilson 1975), and EFF (Elson et al. 1987).For the former two models, SBPs are generated for dynam-ically stable isothermal clusters of reduced kinetic energies, whereas for the latter profile Jeans' equation and the subsequent Poisson's equations are solved to obtain dynamical parameters that are consistent with the observed SBPs.From the structural parameters obtained by nProFit (scale radius r d for EFF or r 0 for King and Wilson, and shape parameters γ for EFF and W 0 for King and Wilson), our proposed tool computes dynamically relevant parameters such as mass, surface and volume mass densities, velocity dispersion and binding energy, as well as tidal radius and core radius.
In §2, we describe the nProFit underlying algorithm, starting with the program initialization, moving subsequently to the background subtraction, SBP extraction, and subsequently describing the fitted models and the corresponding parameters derivation, following the χ 2 minimization technique.In §4 we introduce the simulation tool mksample, and use it to create a mock clusters sample to illustrate the operation of nProFit.Finally, in §5 we show our conclusions and the future directions for nProFit.

THE ALGORITHM
nProFit stands for n-Profile Fitting Tool.This tool was developed to extract and fit the observed SBPs of n objects in an image.This task is carried out by following a series of steps implemented in Fortran, Python and PyRAF routines.In this section, we describe the structure of our code.

The scope
nProFit is developed to obtain the structural parameters of star clusters on the fits format science images of nearby galaxies taken with the HST, such as those in the ACS Nearby Galaxy Survey Treasury (ANGST, Dalcanton et al. 2009), which has more than 60 galaxies at distances 4 Mpc.There are two kinds of star clusters that can be easily detected on the HST images of nearby galaxies -GCs and Super Star Clusters (SSCs).These clusters typically have a half-light radius less than 10 pc, which allows them to be distinguished from stars on the HST images up to distances ∼5 Mpc.The cores of star clusters are supported against the gravity by the pressure exerted by the random motions of stars and hence are well modelled as isothermal spheres of finite kinetic energies.We hence structured our code to fit the observed SBPs with theoretical profiles for families of stable clusters.Theoretical profiles are defined in mass surface density, which is related to the observed SBPs through the mass-to-light ratio, which is assumed to be be 1 M /L and independent of radius in this work.For theoretical clusters containing stars of equal mass and without internal dust such as those defined by King (1966) models, this is a good assumption.However, real clusters have stars of a range of masses, with a tendency for the most massive stars to sink to the center as the cluster dynamically evolves, which produces a color gradient that is bluer towards the center (see e.g.Djorgovski & Piotto 1993).Given that the bulk of the mass of a cluster is in low-mass red stars for Kroupa (2001) and other initial mass functions, SBP in a red filter is expected to trace the mass-density profile better than that in a blue filter.Filters at longer wavelengths also are less affected by possible presence of differential reddening due to patchy internal dust (Cardelli et al. 1989).Among the commonly used HST/ACS filters, F814W is the reddest filter, and hence we recommend the use of F814W images for obtaining SBPs.
Morphological analysis of extragalactic star clusters suffers from two problems: (1) background variationstar clusters are often encountered in zones with varying background values such as in the spiral arms of galaxies, which requires a measurement of a local background value for each object, and (2) object crowding -star clusters are hardly isolated.Both these facts affect the surface brightness of the profile in its external parts.We have built-in algorithms in the code to address both these issues.

Initialization and data preparation
The analysis of SBPs of objects depends on a number of data-dependent parameters.As a first step, the code reads an input ASCII file containing specific analysis parameters.An example of program input can be found in Appendix A. One of the parameters in the input file is the name of a list containing coordinates of n objects.The code obtains 2D sub-images centered on each object in the list.Each sub-image is obtained by trimming the original image of user-defined sizes and centred on the coordinates in the object list, either world coordinates (WCS) or in pixels.These individual images centred at each object coordinates allow us to study the objects separately, since datasets are typically constituted by several objects.

Background estimation and subtraction
Estimating and subtracting the background value accurately is crucial in the determination of structural parameters.An erroneous estimation of the background value may result in an overestimation or underestimation of the derived structural parameters.We implemented two techniques to estimate the background values locally in each sub-image image, which are explained below.
Statistical determination in the corners: This strategy is based on the computation of median values in four corners of each sub-image over box sizes of 10% the size of each sub-image.The use of median, rather than the mean, ensures that each estimated background value is not much affected by any contaminating object in the corners.The availability of four measurements allows us to check the uniformity of the background.The minimum of the four median and root mean square (rms) values are taken as the optimum background and rms values for the object under analysis.
k×σ clipped images: This method uses the whole subimage to obtain an optimum background value for the object.A background image is obtained by iteratively rejecting pixels with values above and below k×σ around the median value in user-defined small boxes.The median and rms values of this background image are taken as the optimum background and rms values for the object under analysis.The obtained values are found to stabilise after ∼ten clipping iterations with k = 3 (i.e. the procedure is performed rejecting values above 3σ to estimate the background value).For this method to obtain reliable background value, the box size should be at least twice the size of the analyzed object.
Both these methods ignore variations in the background value over scales of the sub-image.In principle, gradient in the background values over scales of the cluster size can be obtained by including the background function while fitting the 1D SBPs (Mackey & Gilmore 2003).However, we found that the background value determined by this method does not easily converge to stable values, like the two methods described above.We have used these two techniques to analyse SSCs in M82, and obtained reliable SBPs and therefore accurate structural parameters (Cuevas-Otahola et al. 2020, 2021).The obtained background level values are subtracted to their corresponding 2-D individual images, resulting in background-subtracted sub-images.

Masking contaminants
In some cases, the distance between objects are too small, hindering an accurate profile extraction.Analysis of such clusters requires setting specific constraints on the fitting procedure by either masking the contaminants or setting smaller fitting radii.
nProFit uses the masking capability that the ellipse task provides to obtain 1D SBPs of sources without contribution from close contaminant sources.The mask input file is an ASCII file containing the number of circular masks, along with their spatial coordinates and radii in pixels.
We draw particular attention to cases where a contaminant source is located in the corners of the figure, causing noticeable variations in the background level.If a contaminant source fits the latter scenario, it will be masked by nProFit prior the background level estimation.

Surface brightness profile extraction
Star clusters on the HST images are not resolved enough to obtain SBPs using the star count method that is normally used to analyse star clusters in the Milky Way and the LMC (Brandl et al. 1996;Mackey & Gilmore 2003).We used the ellipse task in the IRAF/STSDAS package (Jedrzejewski 1987), one of the most widely used tools to extract the SBPs from 2D images, for analysing extragalactic star clusters.The ellipse task obtains SBPs by azimuthally averaging intensities in elliptical annular zones.The task allows for variation of center, ellipticity ( = 1 − b/a, where a and b are major and minor axes of the ellipse) and position angle (PA) of the major axis of the successive ellipses.Linear increments in the semi-major axis of the successive ellipses were used to calculate the intensity profiles.We defined concentric ellipses with their centers fixed at the user-supplied object coordinates.The and PA were also fixed either to the user-provided values or to their asymptotic values.In the latter case, nProFit computes the ellipticity value by running the ellipse task without fixing the ellipticity values in a first iteration, and obtains the SBPs in a second run of ellipse after fixing the values of and PA.In this scenario, nProFit analyzes the radial profiles of ellipticities from the first run, setting the ellipticity and PA values as the ones corresponding to the radius at which the ellipticity values are nearly constant.In the majority of the cases, the value where the ellipticity stabilizes matches the radius containing half the total cumulative flux, i.e. the effective radius.During the ellipticity estimation, nProFit avoids the most inner radii that does not provide reliable measurements.The ellipse task obtains cumulative intensities in addition to the SBPs for every fitted ellipse.The k −sigma clipped rms fluctuations in the azimuthal intensities for each fitted ellipse are taken as the error on the measured intensities.All masked pixels are excluded from the analysis during the ellipse fitting.

Fitting radius
Considering that in the majority of datasets, clusters and extended objects in general are not isolated, it is necessary to remove the contribution of contaminant sources (even after masking them) in the fitting procedure.With this aim, we define the fitting radius.For isolated clusters, the background-subtracted SBPs are Illustration of the choice of the fitting radius.The plots show isophotal intensity of two simulated clusters as a function of semi-major axes of the ellipses that best fit the isophotes.The top and bottom panels show clusters without and with a nearby contaminating source, respectively.The axis units correspond to natural units for the processed HST/ACS images, which is in electron/s for intensities and pixels of 0.05 arcsec for the semi-major axis (SMA).
expected to monotonically decrease until the intensity values reach one σ bg , with σ bg being the dispersion of the measured background value intensity.We limit the fitting up to a radius at which the SBP has a value of 3σ bg .We refer to such a radius as R 3σ .For objects located in relatively crowded regions, the SBPs show a bump instead of monotonically decreasing well before the intensity reaches the 3σ bg level.In such cases, we de-fine R ip , the inflection point such that at R ip = d 2 I dR 2 = 0.In Fig 1, we illustrate the fitting radius selection in the case of a cluster with a contaminant source nearby as well as for an isolated cluster.For each object under analysis, nProFit computes both radii and sets as the fitting radius the minimum of R 3σ and R ip .

Dynamical model fitting to observed SBPs
The main goal of our code is fitting the observed SBPs of star clusters to derive the basic structural parameters as well as dynamically relevant parameters.The structural parameters are obtained by fitting the observed SBPs with the theoretical SBPs for selfgravitating static models supported by pressure exerted by the stellar motions.We describe the theoretical SBPs implemented in nProFit below.

THEORETICAL SBPS
The theoretical SBPs are related to the surface stellar density distributions (SDPs), through the mass-to-light ratio of stars.Clusters are self-gravitating and hence the SDPs determine the potential of the cluster.SDPs are obtained based on self-consistent potential-density pairs, given by phase-space distribution functions f (x, v, t), depending on the positions and velocities of stars at a given time (Binney & Tremaine 1987).Distribution functions (DF) allows us to derive the cluster's spatial density profile and, the mean stellar velocity is defined by ν is directly related to the luminosity profile through the mass-to-light ratio associated to the objects' age (Bruzual & Charlot 2003), and is constrained by the collisionless Boltzmann equation with x i , v i and a i , the positions, velocities and accelerations of the stars.It can be noticed that the Boltzmann equation depends on 7 variables, which hinders obtaining a direct solution of the equation.In order to simplify the problem, the Jeans Equation is used.
with β the degree of anisotropy of the velocity distribution, v 2 r the radial velocity, φ the cluster potential.
For the sake of this work, we assume isotropic velocity distributions, for which β = 0, reducing the previous equation to two terms.The Jeans equation is obtained by taking moments of the Boltzmann Equation and integrating over all the velocities.

King models
King models are based on a modified isothermal sphere.The density of the dynamical King models is given by a distribution function (King 1966).The King distribution is described as follows with ε the relative energy of the system (defined in terms of the central gravitational potential and the total system energy as ε = −E+φ 0 , (or equivalently ε = Ψ− 1 2 v 2 ) per unit mass, σ o the model-associated velocity dispersion described as follows with ρ 0 the central density and r 0 the scale factor, referred to as the King Radius in the literature.
From the Poisson's equation (Eq.7) and Eq. 5, setting W = ψ σ 2 , the density function is obtained, in terms of W , the dimensionless potential with erf(x) the error function erf (x) = 2 √ π x o e −t 2 dt.W is obtained by solving ∇ 2 W = 4πGρ(W ), with boundary conditions W (0) = W 0 , with W 0 the central potential, and W (r t ) = 0, with r t the tidal radius.
In Fig. 2 (left panel), we show the density profiles for 4 King models, obtained from the previously described procedure.We show concentrated (W 0 = 2 and W 0 = 6) as well as more extended profiles (W 0 = 10 and W 0 = 15).

Wilson models
Wilson models are based on isothermal sphere models as the King models.These models were originally proposed by Wilson (1975) to fit the observed surface brightness profiles of elliptical galaxies, having larger haloes than the King models, produced by an extra term in the energy distribution function.The distribution function of these models is as follows with ε the relative energy of the system, and σ o the model-associated velocity dispersion described in Eq. 6.
The corresponding mass density function is obtained proceeding analogously as in the King model, and follows The theoretical surface brightness for King and Wilson models can be found from the mass density function by means of the following integral with ζ the mass-to-light ratio.For the sake of simplicity, we assume ζ = 1 M /L throughout this work.In Fig. 2 (middle panel), we show the density profiles for 4 Wilson models, obtained from the previously described procedure.We show concentrated (W 0 = 2 and W 0 = 6) as well as more extended profiles (W 0 = 10 and W 0 = 15).

Moffat-EFF empirical profiles
In the pioneering work by Elson et al. (1987), King models were used to fit the observed profiles of 10 intermediate-age star clusters in the Large Magellanic Cloud, that have masses and densities similar to that of old Globular Clusters (Portegies Zwart et al. 2010).They found King models did not provide good fits to the SBPs because these clusters displayed more extended haloes instead of truncated outer parts.Elson et al. (1987) proposed an empirical profile, based on Moffat (Moffat 1969) profile (henceforth Moffat-EFF where EFF stands for Elson, Freeman and Fall), given by with R the observed profile projected semi-major axis, r d the characteristic radius or model scale radius, L tot the total luminosity, and γ the Moffat-EFF index, which provides information on the shape of the halo.Moffat-EFF profile does not have an implicit distribution function.However, its 3-D luminosity density profile can be calculated using the expression

Volume density and dispersion velocity profiles for isothermal spheres
The King and Wilson models compute the ρ(r) as a function of a numerically defined potential W for isothermal spheres, given by equation 8 and 10, respectively.The ρ(r) and W define the core radius, dynamical mass, surface and volume mass densities, binding energy, bound mass and central velocity dispersion, and hence the dynamically useful parameters corresponding to the best-fit model are known a priori.
In order to obtain ρ(W ), it is necessary to solve the following expression, obtained from the Laplacian oper-ator in spherical coordinates and in terms of the dimensionless potential W and since W depends only on r, and following the prescription by King (1966) the Poisson's equation follows with r = r/r 0 .With the aim of solving these differential equations, we proceed to re-write the equation with W becoming the independent variable, to be consistent with Eqs. 8 and 10, following King (1966); Heggie & Aarseth (1992); Küpper et al. ( 2011) with X = r 2 .Following the dependence of ρ on W , the first step is to use the 4th order Runge-Kutta method (Runge 1895;Kutta 1901) to obtain W in terms of X and, thus ρ in terms of W for each obtained value from each iteration, following where, Re-writing Eq. 16 as and setting y = X, we have y = dX dW , and y = d 2 X dW 2 .The latter expression can be expressed by substituting Eq. 19 Finally we proceed to apply Runge-Kutta two times to obtain the values of W and For King and Wilson models, nProFit calculates the corresponding profiles following the prescription in Binney & Tremaine (1987) where resulting in the following equations for King and Wilson profiles, respectively (25) In Fig. 3, we show the velocity dispersion profiles for King and Wilson models with W 0 = 7.1 and W 0 = 5.7, with the same tidal radius r t for illustration purposes.

Potential and velocity dispersion profiles for Moffat-EFF models
On the other hand, some initially empirical models, such as the Moffat-EFF model, can be physically moti-vated.For instance, by means of Eq. 12, the Moffat-EFF volume density yields where with, I 0 the central surface brightness in units of L /pc 2 and Γ the usual gamma function.We bear in mind that in general, mass profiles do not strictly follow light profiles over all radii due to the effects of mass segregation (Shanahan & Gieles 2015;Baumgardt 2017), resulting in mass functions flatter than the luminosity functions.However, for the sake of simplicity, we assume that mass profiles follow light profiles, and convert from such quantities by means of the mass-to-light ratios.In Fig. 2 (right panel), we show the density profiles for 4 Moffat-EFF models, obtained from the previously described procedure.We show concentrated (γ = 4 and γ = 8) as well as more extended profiles (γ = 2 and γ = 3).
nProFit also computes the theoretical velocity dispersion profile for Moffat-EFF model following the Eq.16 in the prescription by Elson et al. (1987), under the assumption of a spherical cluster under hydrostatic equilibrium and Ω and κ the circular and epicyclic frequencies of the galaxy at the pericenter of the orbit of the cluster, and M (x) the mass enclosed at radius x.
In Fig. 3, we show the velocity dispersion profile for a Moffat-EFF profile with γ = 2.7, compared with King and Wilson profiles with the same r t .
The corresponding potential φ for Moffat-EFF model is computed by means of the Poisson Equation in spherical coordinates resulting in the following equation where (30) with 2 F 1 the hypergeometric function and Γ the usual gamma function.

Derived parameters
We now use our results for best-fit models to extract the most commonly used structural parameter, namely core radius and half-light radius.The latter quantity depends on the concentration parameter for the King and Wilson models, and on the gamma for the Moffat-EFF profile.We also provide three additional parameters namely projected central velocity dispersion, total mass and binding energy of the clusters, for an assumed mass-to-light ratio of 1.0

Concentration index c
For King (1966) and Wilson (1975) models, the central potential W 0 is related to the concentration parameter c = log(r t /r 0 ) obtained by fitting empirical King (1962) formula.In Figure 4, we show this relation, which is obtained from the solutions of the models described in the previous section.

Core radius Rc
The scale size of the isothermal spheres (r 0 ), and Moffat-EFF profiles (r d ) are related to the core radius R c .This quantity is computed recalling that R c is the radius at which the luminosity density reaches half its peak value.
For Moffat-EFF profiles, R c is given in terms of the characteristic radius or model scale radius (r d ) as follows The corresponding values for King and Wilson models are obtained by nProFit by interpolating over the profile density values.In the left-most panel of Fig. 5, we show an auxiliar dimensionless value used by nProFit to compute R c for King and Wilson models.For the sake of completeness, we show the same dimensionless function in terms of γ for Moffat-EFF models.

Half-light radius R h
For a given core size of isothermal spheres, the halfmass radius is related to the concentration index.We show such a relation in the second panel from left to right in Fig. 5, with Wilson models having larger values at the same R h value due to the more extended profiles of the Wilson models.McLaughlin (2000) characterized such a relation by fitting a 9th order polynomial for King models.We have also performed a polynomial fitting, in this case with a 5th order polynomial, for both King (Eq.32) and Wilson (Eqs.33 and 34, for values below or equal and above c = 3.25, respectively) models.log (34) On the other hand, for the Moffat-EFF profile, the R h is analytically related to the fitted structural parameters r d and γ.
We show the dimensionless function R h /r 0 for illustration purposes in the second panel from left to right in Fig. 5, along with the corresponding R h /r d relatiion for Moffat-EFF models, for the sake of completeness.

Central velocity dispersion σ0
nProFit computes de central velocity dispersion values from the previously computed velocity profiles described in Sec.3.4 for King and Wilson models and 3.5 for Moffat-EFF empirical profiles.The central velocity dispersion profile projected into the plane of the sky σ p,0 is computed for isothermal models as well as for Moffat-EFF models following with v 2 r the quadratic velocity dispersion profile computed in Sec.3.4 and 3.5, and I(R) the SBP.

Tidal Radius rt
From the computed velocity dispersion profiles corresponding to Moffat-EFF models, we can compute the tidal radius, by finding the radius r at which σ p (r) = 0. On the other, hand, we find such quantity for King and Wilson models, from the concentration parameter and scale radius, following the expression for the concentration index r t = 10 c r 0 .

Total mass Mtot
Another critical parameter these models provide is the total mass.For Moffat-EFF models, nProFit derives the model mass from the total luminosity L tot , using Eq. 12 and assuming a mass-to-light ratio ζ.The mass corresponding to King and Wilson models is computed from Eqs. 38 and 40 in King (1966) instead, resulting in with j = j j0 (expressed as a normalised quantity, as obtained in the previous sections), r t = rt r0 and r = r r0 .In the right-most panel in Fig. 5, we show auxiliar dimensionless function used by nProFit to compute M tot for King and Wilson models, along with Moffat-EFF profiles for completeness..

Binding energy E b
Following the prescription by McLaughlin (2000), we have also implemented in nProFit the computation of the binding energy, (E b ) using their Eq. 1 nProFit computes E b for Moffat-EFF models, substituting the potential described in Eq. 30 along with the density ρ in Eq. 26 into the following equation On the other hand, for King and Wilson models, nProFit uses the dimensionless function in the rightmost panel in Fig. 6 to compute the integral in the following equation 3.6.8.Central surface magnitude µ0 The central surface brightness I 0 in L pc −2 is converted to observational units (µ F,0 ) using with M F, the absolute magnitude of the Sun in the filter F, provided by the user (see Willmer ( 2018)).

Library of dynamical models
For King and Wilson models, nProFit solves the previously described set of differential equations with boundary conditions for W 0 ∈[2,15] in steps of 0.1 in terms of r r0 .The latter results in a library constituted by files containing r r0 , W , ρ, and Σ.These libraries are the results of the most computer-intensive module of nProFit.They are pre-evaluated to speed up computation times and the corresponding results are used to compute the previously described derived parameters.nProFit uses this library for any new fit, unless the user specifically asks to solve the equations.

Convolution with the Point Spread Function (PSF)
The theoretical profiles are convolved with the userprovided PSF 2 by nProFit using a Fortran routine.We recall that in general, for finite datasets as the models we described, the convolution operator is given by 2 PSF for the HST images can be obtained using the Tinytim tool (see https://www.stsci.edu/hst/instrumentation/focus-andpointing/focus/tiny-tim-hst-psf-modeling).Alternatively, it can be defined for each frame using the tasks for that purpose in daophot package or Sextractor command PSFex.In the present work, we assumed a Gaussian profile of FWHM=2.1 pixels, which represents well the point sources on the HST/ACS images.with g the signal and f the response function, which in our case are the dynamical models and the PSF, respectively.Computing Eq. 41 requires N 2 operations to perform the convolution.In order to reduce the number of operations, nProFit uses the fast fourier transform (FFT), resulting in N log 2 N operations, reducing the number of operations in 65%.nProFit computes the FFT using the numerical recipes routine convlv (Press et al. 1992).In general, the FFT splits the convolution expression in terms of odd and even indices in the sum.Hence, FFT routines require the input sampled in number of points being a power of 2. It is also required for the convlv routine, that the PSF or response function is wrapped around and filled with zeros to obtain a vector with a length of N-1 (odd size).N is the size of the vector containing the dynamical model, which is a power of two.The model needs to be mirrored and wrapped around.In Fig 7, we illustrate the convolution procedure along with the required zero-padding and wrapping.

Selection of the best-fit model
From the previously determined fitting radius, nProFit proceeds to compute which model provides the most accurate and reliable representation of the observed SBPs.To this aim, we use the non-parametric statistical test χ 2 (Bevington et al. 1993).The χ 2 test determines the goodness of a fit to data, suitable for observed data with gaussian errors, which by virtue of the central limit theorem (Alexander McFarlane Mood 1974) are good approximations of Poisson distributions.Poisson statistics are crucial in the determination of the data noise, having different variances, requiring to use a χ 2 test weighted by errors (Wall & Jenkins 2003) with N pts the number of points, the azimuthallyaveraged profile intensities I obsi , σ i the corresponding errors computed during the isophotal fitting, and Ĩmodeli the PSF-convolved model intensities at i, varying from 1 to N pts.nProFit performs the χ 2 minimization technique to find the best-fit model in the parameter space for each one of the available theoretical models in nProFit.In order to determine which of these models represents most accurately the observed data, we implemented the prescription by McLaughlin & van der Marel (2005) to compare the obtained best-fit models where χ 2 ref and χ 2 alt are the χ 2 min values of the reference model and the model to be compared, respectively.
This procedure allows to determine in which cases two models are equally good |∆χ 2 | ≤0.2 or whether the alternative model provides a better fit than the reference model.

Determination of errors
nProFit computes the errors on the obtained parameters by considering 1-σ significance regions.Considering that each fit depends basically on three free parameters, the intervals for a 1 − σ level of confidence are defined by means of the following equation (Wall & Jenkins 2003) with χ 2 min the χ 2 value of the best fit model, and χ 2 the value corresponding to the rest of the models.These intervals are computed for Moffat-EFF, King and Wilson models, separately.The errors on the derived parameters described in Sec.3.6 are computed by nProFit by propagating the errors on the basic parameters for each models, namely, r d and γ for Moffat-EFF models and r 0 and W 0 for King and Wilson models.

CODE ILLUSTRATION WITH SIMULATED CLUSTERS
In Cuevas-Otahola et al. (2020, 2021), we have applied the techniques described in this algorithm, and compared our results with those obtained using the publicly available tools Galfit and Ishape.Our results for the sample of super star clusters in M82 were in agreement with the results obtained by Galfit for Moffat-EFF models.More specifically in Cuevas-Otahola et al. ( 2021), we computed the derived parameters of the sample for Moffat-EFF models.In order to illustrate these techniques for King and Wilson models as well, we simulated a sample of 105 clusters following King, Wilson and, for the sake of completeness, we also simulated clusters following Moffat-EFF model profiles.With the aim of testing the code in a realistic scenario, considering that typically, star clusters are embedded in crowded regions, we used an archive image as the background of our mock clusters.We used an image in F555W filter extracted from the HST Legacy Survey, provided by the Hubble Heritage Team, of the prototype starburst galaxy M82, which is a complex study case, due to its high crowding, high inclination angle and background gradient.In Fig 8, we show a simulated image containing 105 clusters along with the zoom on sub-images around 5 of the clusters in the mock sample (the background of the mock sample has considerable gradient).The mock sample is based on the parameters set in Tab. 1, simulated several times for different central surface brightness.

Sample simulation
We have designed and implemented the subroutine mksample to generate a mock sample, from user given coordinates and model type.
Mock sample of clusters can be generated using our module mksample from user given coordinates and model type.mksample generates 2D images following the projected profiles of the available models.The synthetic profiles data are stored in a 2D matrix, which is generated following a similar procedure as that in the IRAF task mkobjects.We draw particular attention to the geometrical features of the models.As in the case of mkobjects, we consider the profiles of models with spherical symmetry.However, to reproduce the axysymmetry of some observed objects, we introduce an artificial ellipticity in the x, y model coordinates as Table 1.Set of initial parameters used to build the mock data with Moffat-EFF (M), King (K), and Wilson (W) 45) with PA and AR (AR=b/a, with a and b, the semimajor and semi-minor axis, respectively), the position angle and axis ratio, respectively.x i and y j are the (i, j) coordinates in the x and y direction of the image, and x c , and y c the centers of each object.We draw attention to the mock elliptic clusters, whose parameters are well recovered using spherical models, since as we have shown in Appendix A in Cuevas-Otahola et al. (2020), for the most elongated cluster in M82, M82-F, the ellipticity does not considerably affects the integrated profile obtained from the isophotal fitting.Each profile is simulated as follows sim(i, j) = S 0 mod mod(r, par1, par2) + bg (46) with r = x 2 mod + y 2 mod , mod the selected model (Moffat-EFF, Wilson or King) S 0 mod the central surface brightness used for the simulation, par1, par2, the model parameters (r d and γ for Moffat-EFF models, and r 0 and W 0 for King and Wilson models), bg is the simulated background value for each coordinate.For the sake of testing nProFit in a realistic background scenario, we used a real archive image of a galaxy instead of the usual background image drawn from a Gaussian distribution with mean and sigma values provided by the user.We used an image from the HST Legacy Survey made publicly available by the Hubble Heritage Team (Mutchler et al. 2007) corresponding to the prototype starburst galaxy M82, located at a distance of 3.63 Mpc (Freedman et al. 1994).Such a galaxy represents an interest study case, ideal to determine the extents of the code, considering that it has a strong background gradient, a large number of clusters along its disk and nucleus (around 600) (Mayya et al. 2008) and a high inclination degree 77 • (Mayya et al. 2005).In addition to the realistic background conditions, the mock sample spans a wide range of central surface brightness values (spanning around 4.5 mag/arcsec 2 , from 18.5 mag/arcsec 2 to 14 mag/arcsec 2 , in 5 bins), initially set to test the accuracy of the code and its dependence on the surface brightness profiles of the clusters. In

Structural parameters of the mock sample
We run the nProFit code over our synthetic data, and obtained in the first place the sub-images centered in each object.Background subtracted images are subsequently obtained.
In Fig. 10 we show examples of SBPs on our simulated image for two clusters.The intensities are given in units of raw counts per second (cps) on the background subtracted images.These are converted to mag arcsec −2 units using the zeropoint and image scale given by the user.The profile semi major axis is given in units of pixels, which are converted to parsecs, using the image scale given by the user, to match the theoretical models units for fitting purposes.
In Fig. 11, we show an example of an extracted surface brightness profile along with the corresponding best-fit model.The corresponding residuals are also shown.In order to compare in a more homogeneous way the results obtained by nProFit, we computed one of the most relevant quantities, the half-light radius (R h ), and compared it with their initial values.In order to ensure that the obtained results are reliable, regarding the fitted model, we determined, the best fitting model in each case, following Eq.43, as we show in Fig. 12.In panel (a), we show the comparison between the ref-   We fit the three models to all the mock sample data, regardless of the models they were drawn from.
In order to test the accuracy of nProFit, in Fig. 13 (a), we show the difference between the initial R h values and the values obtained by nProFit, weighted by the initial R h values, and compared with them, as a function of the difference between the central surface brightness and the measured background value (∆µ).As expected, we observe a better recovery for larger ∆µ values, even for small clusters, with larger errors, considering that the latter ones, have sizes close to the limit image resolution.
For the sake of this analysis, we removed extremely faint clusters in very crowded areas, resulting in fitting radius values shorter than 8 pixels, resulting in a subsample of 74 clusters.On average, we notice that the recovered values for clusters with ∆µ < 5 mag/arcsec 2 are within 20% of the initial values, whereas, for clusters with ∆µ > 5 mag/arcsec 2 , the recovered values are within 10% of the initial values.
In order to validate our results and compare them with the corresponding ones obtained by other publicly available tools, we carried out the structural parameters fitting using Galfit (Peng et al. 2010) and Ishape (Larsen 1999).In Fig. 13 (b), we compare the results obtained by Galfit (version 3.0.5),with the initial simulated values.We obtained converging fits for 55 clusters, which are on average among the brightest ones.We notice that Galfit provides fits comparable to nProFit in all size ranges.On the other hand, in Fig. 13 (c), we compare the results obtained by Ishape, with the initial simulated values.For clusters with large contrast (∆µ ≥5 mag/arcsec 2 ), Ishape provides equally good results for 70% of the clusters.However, the R h values values obtained by Ishape have larger dispersion and errors for ∆µ <5 mag/arcsec 2 .Cuevas-Otahola et al. (2020) had demonstrated that the core radii for real M82 clusters are well reproduced by Ishape.Hence, the large error on R h is most likely due to poor recovery of halo parameters in Ishape for clusters located in high background regions.
We have also empirically found the R h values by integrating the observed surface brightness profiles of each cluster (and corrected them by the PSF radius) and we compare them against the initial R h values in Fig. 13 (d).The recovery is especially poor for R h <4 pc, with the recovered values systematically larger.

CONCLUSIONS
In this work, we present the numerical code nProFit, devoted to obtain the best-fit structural parameters of star clusters in the HST images of nearby (distance < 5 Mpc) galaxies.The code is Python-based at the user end, but uses modules of Pyraf and Fortran.nProFit extracts sub-images centered in each analyzed object coordinates.Subsequently, a local background estimation is carried out by nProFit determining the median values in the corners or by scanning the whole images and estimating the background value by a σclipping procedure.The estimated value is subsequently subtracted from the extracted SBPs by nProFit from isophotal fittings.PSF-convolved Moffat-EFF, King and Wilson models are then fitted to background subtracted azymuthally averaged surface brightness profiles.nProFit uses a χ 2 -minimization technique to fit the models.As a result, the tool provides the set of basic structural parameters, scale parameters (r d for Moffat-EFF and r 0 for King and Wilson models) and shape parameters (γ for Moffat-EFF and W 0 for King and Wilson models).Since nProFit fits dynamical models, it offers a valuable opportunity to derive physically-relevant parameters.Among these parameters are central volume and luminosity densities (ρ 0 and j 0 ), total masses and luminosities (M and L), central velocity dispersions (σ 0 ), core radius (R c ), half-light radius (R h ), tidal radius (R t ) and binding energy (E b ).
We have tested nProFit on simulated clusters superposed on real HST images.For the simulated clusters, the surface brightness difference between the cluster maximum and the local background varies between 3 to 8 mag/arcsec 2 .We demonstrate that the input values are recovered within the 1-σ errors for majority of the simulated clusters for all the three theoretical models we have explored.The R h values are recovered within 10 percent for clusters with ∆µ >5 mag/arcsec 2 and 20 percent for clusters with ∆µ <5 mag/arcsec 2 .The accuracy of our recovery is comparable to that of Galfit, whereas it is clearly better than that with Ishape, especially for clusters ∆µ <5 mag/arcsec 2 .We illustrate that nProFit is a tool suitable to fit the structural parameters in samples with considerable crowding, such as the M82 disk, providing reliable values for clusters for which a neighbouring cluster does not contribute significantly within a distance of 8 HST/ACS pixels.
As a final note, we clarify that in this work, we considered that mass profiles follow light profiles over all clusters' radii.This is a simplification, since, mass segregation influences the mass density profiles of clusters, as well as their central velocity dispersions, causing the mass profiles to depart from the corresponding light profiles.For this reason, we will include mass segregation prescriptions in the upcoming version of nProFit .

Figure 2 .
Figure 2. Normalized volume densities corresponding to King (left panel), Wilson (middle panel) and Moffat-EFF (right panel).The densities corresponding to isothermal-based models (King and Wilson), were computed by nProFit.

Figure 3 .
Figure 3. Velocity dispersion profiles for King, Wilson and Moffat-EFF models, with the same tidal radius rt, and volume density ρ = 10 4.2 M /pc 3 .The King and Wilson profiles have W0 = 7.1 and W0 = 5.7, respectively, and the corresponding γ index for the Moffat-EFF profile is 2.7

Figure 4 .
Figure 4. Relation between the concentration index c and the dimensionless potential W0, computed by nProFit following the procedure described in Sec.3.4 for King and Wilson models.

Figure 5 .
Figure 5. Normalised functions used to compute Rc (left-most panel), R h (second panel from left to right), and auxiliar functions for I0/jor0 (third panel) and Mtot (right-most panel).These functions are used by nProFit, to speed up calculations.

Figure 6 .
Figure 6.Auxiliar function to compute the integral in Eq. 39 for the binding energy.

Figure 7 .
Figure 7. Convolution scheme performed by nProFit.The model is wrapped around and the PSF is mirrored and zeropadding is performed to meet the required size criterion for the FFT.The plotted example corresponds to 2 7 = 128 pixels.

Figure 8 .
Figure 8. Top panel: Synthetic clusters (105 green circles) superposed on a real HST image, which corresponds to the F555W image.Bottom panels: zoom in of five clusters of the mock sample, in images of 101x101 pixels, centered in clusters S24, S43, S49, S93 and S94 (from left to right, respectively).We show in the images bars indicating 1 kpc and 10 pixels, in the top and bottom panels, respectively.Throughout this work we use an image scale of 0.88 pc pixel −1 , which corresponds to the physical size of the HST/ACS pixels at the distance of M82 (3.63 Mpc).

Figure 9 .
Figure 9.Initial conditions of the mock simulation, along with the values obtained by nProFit, with the corresponding 1−σ confidence intervals for Moffat-EFF (left-most panels), King (middle panels), Wilson (right-most panels) models for faint (top panels) and bright (bottom panels) clusters.

Figure 10 .
Figure 10.Azimuthally averaged surface brightness profiles computed by nProFit through the ellipse IRAF task for the objects S93 and S28, in the left and right panels, respectively.

Figure 11 .
Figure 11.Dynamical models (Moffat-EFF, King and Wilson) fitting performed by nProFit to the surface brightness profile of a cluster in the mock sample generated with mksample (top panel).Fitting residuals (bottom panel).The fitting radius is shown with the vertical dashed line, the observed SBP with the empty dots and the fitted models with solid lines.

Figure 12 .
Figure12.∆χ 2 values to determine the best fitting model, setting as a model for comparisonMoffat-EFF (a),King (b)   andWilson (c)  versus percentage error for R h fits performed by nProFit, for the clusters initially fitted with the models shown in the figure legend.The green horizontal gaps represent a difference of 20%, between the compared models, where the compared fits provide equally good results.

Figure 13 .
Figure 13.Percentage error for R h fits performed by nProFit (a), Galfit (b), Ishape (c), and from an empirical estimate (d) versus the initial half-light radius, with the point sizes coded as a function of the difference between central surface brigthness and local background values for each cluster.The dotted horizontal line represents errors equal to zero.
King Radius for King and Wilson models, respectively.Each pair of parameters was simulated for five different values of central surface brightness (mag/arcsec 2 ): 18.5, 17.0, 16.3, 15.3, and 14.The varying galaxy background reaches values from 22.2 to 18.6 mag/arcsec 2 , with a median value of 21 mag/arcsec 2 , resulting in a heterogeneous sample with differences between the central surface brightness and the background value ∆µ from 2 to 7.5 mag/arcsec 2 , with a median value of 4.7 mag/arcsec 2 .

Table 2 .
Set of initial and fitted parameters obtained by nProFit for the mock sample , and with background values between 18.6 mag/arcsec 2 and 22.2 mag/arcsec 2 .We performed the fits using nProFit, and find that, the obtained parameters, are well within the 1 − σ confidence intervals centred at the initial values, in the majority of cases, showing fitted values closer to the initial values for brighter clusters.Our mock sample is constituted by these 21 combinations of structural parameters (7 per each model, summarized in Tab.1), along, with the 5 bins in surface brightness, resulting in a mock sample of 105 clusters, with each combination of structural parameters (γ and r d for Moffat-EFF, and W 0 and r 0 for King and Wilson models), simulated for each one of the initial central surface brightness profiles.The mock sample clusters positions were drawn from a uniform distribution limited by the galaxy geometry.
King andWilson models are well-fitted either by King or Wilson models.We computed the R h values, from the fitted r 0 and W 0 , for clusters best fitted either by King or Wilson models, and from r d and γ for those best represented by Moffat-EFF models, following the prescription in Sec.3.6.3.