A method for constructing an analytical model for fitting globular cluster projection surface density profile

Abstract

The application discloses an analytic model construction method of fitting a globular cluster projection plane density profile, relates to the technical field of cluster radial profile analysis in astronomical observation, and introduces a turning radius parameter to divide the profile into a power-law exponent inner region and a classic King model outer region, wherein the power-law exponent directly quantifies the central collapse degree, the turning radius calibrates the spatial scale of abnormal peaks, and the classic King model describes the tidal truncation characteristics; when the data itself does not support the central peak, the influence of the power-law segment in the inner region can be ignored, overfitting of ordinary clusters without the central peak is avoided, a mixed prior constraint mechanism is designed in the sampling of independent parameters, namely, in the Bayesian sampling process, a uniform prior is adopted for the parameters in the boundary, and a wide Gaussian soft constraint is introduced for the parameters outside the boundary, so that Markov chain stagnation and boundary distortion caused by the hard boundary are avoided, and finally the globular cluster projection plane density profile is fitted with high precision.

Description

Technical Field

[0001] This invention relates to the field of radial profile analysis technology for star clusters in astronomical observations, and particularly to a method, apparatus, equipment, and medium for constructing an analytical model for fitting the projection surface density profile of a globular star cluster. Background Technology

[0002] As one of the oldest celestial systems in the Milky Way, globular clusters are crucial observables for studying stellar dynamics, mass stratification, and the effects of galactic tidal fields, particularly the radial profile of their projected areal density. To quantitatively extract cluster structural parameters from observational data, analytical models are widely used in astronomy to fit the areal density profile. The classic and most widely applied model is the King model, proposed by King in 1962. This model, by introducing an energy cutoff mechanism, successfully describes the finite radius (tidal radius) formed by the tidal stripping of the cluster's periphery by galactic forces, overcoming the limitation of infinite mass in earlier isothermal spherical models. The King model includes three parameters: central areal density amplitude, core radius, and tidal radius. Its physical picture is clear, and it can effectively describe clusters with flat cores and tidal cutoff characteristics.

[0003] However, with the application of high-resolution observation instruments such as the Hubble Space Telescope, researchers have discovered that the central region of some globular clusters (especially those that have undergone core collapse) is not a flat "core," but rather has a significant density spike, meaning that the surface density increases sharply in a power-law manner as the radius decreases, until it reaches the limit of observational resolution. For such clusters, there are obvious technical defects when using a single King model to fit the entire radial region: the model underestimates the true density in the central region and produces a systematic bias at the core radius, resulting in a non-random distribution of the fitting residuals, which cannot accurately capture the dynamic state inside the cluster.

[0004] To overcome these shortcomings, researchers have attempted various improvement schemes, such as designing generalized King models and multi-mass King models, and combining them with nonparametric Bayesian methods for final fitting. However, when fitting the entire radial region, such King models require the introduction of additional physical assumptions that are difficult to directly observe. Introducing additional parameters can lead to overfitting of ordinary star clusters without a central peak, making it difficult to uniformly describe the central density peak and the tidal truncation of globular clusters to accurately capture the dynamic state inside the cluster. Furthermore, during the fitting process, the Markov chain Monte Carlo method is used for Bayesian parameter estimation. However, in parameter estimation, a "hard boundary" prior constraint is generally applied to the parameters. But when the sampling chain attempts to cross the preset boundary, the acceptance probability of the "hard boundary" prior constraint drops to zero instantly, causing the Markov chain to be repeatedly rejected near the boundary and fall into "stagnation". This severely reduces the traversal efficiency of the parameter space and even causes artifact bias in the boundary estimation of the posterior distribution, leading to chain stagnation and boundary distortion. Ultimately, this makes it difficult to achieve high-precision fitting of the projected surface density profile of globular clusters. Summary of the Invention

[0005] This invention provides a method for constructing an analytical model for fitting the projection surface density profile of globular clusters, which can solve the problems existing in the prior art.

[0006] This invention provides a method for constructing an analytical model for fitting the projection surface density profile of globular clusters, comprising the following steps: We acquire observational data of globular clusters and predetermine a turning radius parameter to divide the radial profile of the globular clusters within the observational data into two independently analytically describable regions: an inner region and an outer region. In the inner region, a power-law descent function is used to describe the density peak in the central region of the cluster, while in the outer region, the classical King projection analytical model is used to describe the tidal truncation characteristics of the cluster. This is to construct a piecewise analytical function to describe the variation of the projection surface density of the globular cluster with radial distance, and to obtain multiple independent parameters of the piecewise analytical function. Based on the observation data, the maximum likelihood estimation method is used to obtain the initial estimates of each independent parameter; Based on the initial estimates of each independent parameter, Markov chain Monte Carlo sampling is performed on each independent parameter under the Bayesian inference framework. A hybrid prior constraint mechanism is adopted in the log probability function during sampling: for independent parameters within a preset uniform prior range, the sum of the uniform prior log probability and the log likelihood of the independent parameter is directly obtained as the prior log probability; for independent parameters exceeding the preset range, the log probability density of each independent parameter based on a preset width Gaussian distribution is obtained, and the log probability density of each independent parameter is summed as the prior log probability. Statistically analyze the posterior samples obtained by sampling to obtain the median estimated values and confidence intervals of each parameter, so as to construct a piecewise analytical model for fitting the projected surface density profile of globular clusters and generate the projected surface density profile of globular clusters after fitting.

[0007] Preferably, the construction of the piecewise analytical function includes: Preset a turning radius parameter z to divide the radial profile of the globular cluster into an inner region and an outer region, and set the projected surface density of the globular cluster as a function Σ(r) of the radial distance r, where r is the angular distance from the center of the cluster. Then the piecewise analytical function is expressed as: When r < z, the piecewise analytical function is expressed as: ; When r ≥ z, the piecewise analytical function is expressed as: ; Where: a Represents the central surface density amplitude parameter; Represents the core radius; Represents the tidal radius; k represents the inner region power-law index; z Represents the turning radius; Represents the King projected analytical model.

[0008] Preferably, the acquisition of the initial estimated values of each independent parameter includes: Assume that the observed data contains N radial intervals, the radial distance of the i-th interval is r i , the measured surface density is D i , the number of star counts in this interval is N i , assuming that the surface density errors of each point in the observed data follow a Gaussian distribution, construct the log-likelihood function as: ; Where: ; Represents the predicted value of the piecewise analytical function at the corresponding radius of the i-th interval; Based on the numerical optimization algorithm, maximize the constructed log-likelihood function to obtain a set of initial estimated values θ MLE .

[0009] Preferably, the sampling process of performing Markov chain Monte Carlo sampling on each independent parameter respectively includes: Based on a set of initial estimated values θ MLE , perform Markov chain Monte Carlo sampling on each independent parameter respectively under the Bayesian inference framework, and in the calculation of the log-probability function of the Markov chain Monte Carlo sampling, judge whether each independent parameter is within the uniformly prior range set independently for each; When the independent parameters are within the set uniform prior range [ Within this range, the sum of the uniform prior log probability and the log-likelihood of the independent parameters is directly obtained as the prior log probability. ; When the independent parameters are within the set uniform prior range [ Outside of this range, obtain the normal log probability density of each independent parameter with the midpoint of the preset range as the mean and one-five-thousandth of the range width as the standard deviation, and sum the log probabilities of each independent parameter as the prior log probability. .

[0010] Preferably, generating the fitted globular cluster projection surface density profile includes: Statistical analysis is performed on the posterior samples obtained from the sampling to obtain the sample set of each independent parameter. The preset first percentile is used as the median estimate, and the preset second and third percentiles are used as the lower and upper error bounds, respectively, to quantify the confidence interval of the fitting result; Based on the selected median estimate and confidence interval, the fitted segmented surface density profile curves of globular clusters and their confidence interval filling bands are generated using the sample sets of each independent parameter.

[0011] Preferably, the power-law exponent k is used to quantify the steepness of the density peak at the center of the globular cluster; When k>0, it indicates that there is a density peak at the center of the globular cluster; the larger the value of k, the steeper the peak. When k approaches 0, the piecewise analytical model degenerates into the truncation of the classical King projection analytical model within the turning radius z.

[0012] Preferably, the Markov chain Monte Carlo sampling employs an affine invariant stretching-movement algorithm and is configured with multiple parallel chains. Its sampling parameters include the number of chains, the number of sampling steps, the number of combustion phase steps, and the refinement interval.

[0013] This invention also provides an analytical model construction device for fitting the projection surface density profile of globular clusters, comprising: The piecewise analytical module is used to acquire observational data of globular clusters and pre-set a turning radius parameter to divide the radial contour of the globular cluster within the observational data into two independently analytically described regions: an inner region and an outer region. In the inner region, a power-law descent function is used to describe the density peak in the central region of the cluster, while in the outer region, the classical King projection analytical model is used to describe the tidal truncation characteristics of the cluster. This is to construct a piecewise analytical function to describe the variation of the density of the projected surface of the globular cluster with radial distance and to acquire multiple independent parameters of the piecewise analytical function. The analytical model building and fitting module is used to obtain initial estimates of each independent parameter based on the observed data using the maximum likelihood estimation method. Based on the initial estimates of each independent parameter, Markov chain Monte Carlo sampling is performed on each independent parameter under the Bayesian inference framework. A hybrid prior constraint mechanism is adopted in the log probability function during sampling: for independent parameters within a preset uniform prior range, the sum of the uniform prior log probability and the log likelihood of the independent parameter is directly obtained as the prior log probability; for independent parameters exceeding the preset range, the log probability density of each independent parameter based on a preset width Gaussian distribution is obtained, and the log probability density of each independent parameter is summed as the prior log probability. The posterior samples obtained from the sampling are statistically analyzed to obtain the median estimate and confidence interval of each parameter, so as to construct a piecewise analytical model for fitting the projection surface density profile of the globular cluster and generate the fitted projection surface density profile of the globular cluster.

[0014] This invention also provides an electronic device, including a memory and a processor; The memory is used to store computer programs; When the processor executes the computer program stored in the memory, it implements the steps of the analytical model construction method for fitting the projection surface density profile of a globular cluster as described above.

[0015] This invention also provides a computer-readable storage medium for storing a computer program, which, when executed by a processor, implements the steps of an analytical model construction method for fitting the projection surface density profile of a globular cluster as described above.

[0016] This invention provides a method for constructing an analytical model for fitting the projection surface density profile of globular clusters. Compared with existing technologies, its advantages are as follows: This invention creatively divides the contour into a power-law inner region and a classical King model outer region by introducing a turning radius parameter. This allows the model structure itself to describe two distinctly different physical regions. The power-law exponent directly quantifies the degree of central collapse, the turning radius calibrates the spatial scale of the anomalous peak, and the classical King model describes tidal truncation characteristics. This piecewise continuity mathematical construction means that when the data itself does not support a central peak, the power-law exponent in the fitting result will naturally approach zero, or the turning radius will shrink to near the minimum observation radius. In this case, the influence of the inner power-law segment can be ignored, avoiding overfitting to ordinary star clusters without a central peak. First, a unified description of the density peak at the center of the globular cluster and the tidal cutoff at the periphery are established, and independent parameters describing the internal dynamic state of the cluster are accurately quantified. Then, a hybrid prior constraint mechanism is designed in the Markov chain Monte Carlo sampling of the independent parameters. That is, in the Bayesian sampling process, a uniform prior is used for the parameters inside the boundary, and a wide Gaussian soft constraint is introduced for the parameters outside the boundary. This transforms the boundary from an infinitely deep "cliff" into a "slope" with a finite probability, so as to avoid Markov chain stagnation and boundary distortion caused by hard boundaries. This ensures the statistical rigor of the estimation of the uncertainty of independent parameters and the calculation of confidence intervals. Finally, the projection surface density profile of the globular cluster is fitted with high accuracy. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the overall process provided for an embodiment of the present invention. Detailed Implementation

[0018] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0019] Globular clusters are among the oldest celestial systems in the Milky Way, formed by tens of thousands to millions of stars bound together by gravity. As a natural laboratory for studying stellar dynamics, mass stratification, and the gravitational potential constraint of the Milky Way, the internal structure of globular clusters has always been a core research topic in astrophysics. Among these, the radial distribution profile of the projected surface density—that is, the variation of the stellar number density per unit angular area with distance from the cluster center—is the most directly observable structural characteristic, containing rich information about cluster formation, evolution, and interaction with the Milky Way's tidal field.

[0020] From an observational perspective, the surface density profile of globular clusters typically exhibits the highest density in the central region, gradually decreasing with increasing radial distance, and approaching background levels at a certain radius. Early star counting studies based on photographic plates revealed this basic morphological characteristic. With the advent of the charge-coupled device (CCD) era and the introduction of space telescopes, observational resolution and photometric depth have been greatly improved, allowing researchers to detect fainter stars and resolve regions closer to the cluster center. High-resolution observations have revealed that the central regions of some globular clusters do not present a flat "core" structure as early theories predicted, but rather exhibit significant density spikes, meaning that the surface density continuously and sharply increases with decreasing radius until reaching the observational resolution limit. This phenomenon is particularly prominent in core-collapsed clusters.

[0021] From a theoretical perspective, the density profile of a globular cluster reflects the result of the system's long-term evolution under the combined influence of its own gravity and external tidal forces. Classical theory suggests that during the relaxation process, the core region of the cluster tends to contract, and the outer stars are tidally stripped away, gradually forming a high-density core and a truncated periphery. For the core collapse process, theory predicts that the central density should exhibit a power-law peak shape, and the power-law exponent is related to the stellar mass spectrum and physical processes such as binary star heating.

[0022] King (1962) proposed an analytical model to describe the density distribution of globular clusters. This model was developed to address the limitations of earlier isothermal spherical models. The isothermal spherical models assumed that the velocity distribution of stars followed a Maxwell-Boltzmann distribution, resulting in a non-zero stellar density at infinity, which contradicted the observed definite boundaries of globular clusters. King's model, by introducing an energy cutoff mechanism, successfully constructed a dynamical model with finite mass and finite radius, enabling a more realistic description of the distribution of stellar systems. The King model includes three free parameters: the central surface density amplitude *a*, the core radius, and the... and tidal radius The core radius defines the radial distance at which the surface density drops to half its central value, while the tidal radius defines the outer boundary of the star cluster constrained by the Milky Way's tidal forces, beyond which the surface density drops to zero.

[0023] The King model, due to its clear physical picture—the core region is approximately isothermal and the periphery is tidal-truncated—and its relatively simple mathematical form, has been widely used to fit the areal density profile of globular clusters since its inception, becoming a standard tool in the field. The King model fitting parameters of a large number of clusters have been compiled into star catalogs, providing basic data for statistical research.

[0024] However, the King model is not applicable to all globular clusters. For clusters that have undergone core collapse, the density profile in the central region deviates from the isothermal core and shows a steeper power-law decline. In this case, using a single King model to fit the entire radial range will result in under-fitting in the central region - the model underestimates the central density and overestimates the core radius, and the fitting residuals show systematic deviations in the central region. To address this problem, researchers have proposed various improvement schemes, including the generalized King model, multi-mass King model, non-parametric Bayesian method, etc. These methods either introduce additional physical assumptions or have a high computational complexity.

[0025] Aiming at the problems existing in the current generalized King model, multi-mass King model, non-parametric Bayesian method, etc., the present invention proposes a piecewise analytical model for fitting the projected surface density profile of globular clusters, as Figure 1 shown. Within the turning radius, the surface density uses a power-law form to describe the possible density spike in the center. Outside the turning radius, the classical King model is retained to utilize its mature physical framework for describing tidal truncation. The two-segment function is forced to be continuous at the turning point, reducing the redundant degrees of freedom between parameters. The model contains a total of five free parameters, adding two parameters for describing the power-law behavior in the inner region compared to the classical King model. In terms of parameter estimation method, Markov chain Monte Carlo sampling under the Bayesian framework is adopted to obtain the posterior probability distribution of the parameters, rather than just giving a point estimate. This choice enables the uncertainty of the fitting result to be fully quantified, and the covariance structure between parameters can also be intuitively presented through the joint posterior distribution. To improve the MCMC sampling efficiency, a hybrid prior constraint mechanism is designed to avoid chain stagnation caused by hard boundaries while maintaining parameter constraints. At the program implementation level, a centralized parameter configuration architecture is adopted, encapsulating all user-adjustable parameters in an independent configuration area at the top of the code. This design allows researchers who are not professional in programming to apply the program to different cluster data by only modifying a small number of variables, greatly reducing the usage threshold of the method. Specifically, it includes: 1. Piecewise power-law-King model.

[0026] 1. Model definition and mathematical expression.

[0027] Let the projected surface density of the globular cluster be a function of the radial distance r, where r is the angular distance from the cluster center, usually in arcseconds. The piecewise model proposed by the present invention is defined as follows: When r < z: (1) When r ≥ z: (2) Where: aThis represents the center surface density amplitude parameter, corresponding to the theoretical value of the King model at r=0 (the actual King projection model takes a finite value at r=0, and a controls the scaling of the overall amplitude). It represents the core radius and defines the transition scale of the areal density profile as it descends from the flat core outwards; Represents the tidal radius, defining the outer boundary of the star cluster, as r→ The surface density tends to zero; k represents the power law exponent of the inner region, describing the steepness of the surface density increase in the central region as the radius decreases. k>0 indicates that there is a density peak in the center, and the larger the value of k, the steeper the peak. z This represents the turning radius, the boundary radius that separates the inner power-law behavior from the outer King behavior.

[0028] function For the King projection analytical model, the astropy.modeling.models.KingProjectedAnalytic1D class is called in the program implementation of this invention for calculation. The mathematical definition of this class follows King's (1962) projection surface density formula, which is monotonically decreasing in the range r≥0 and decreasing in the range r = The value is zero at that point.

[0029] 2. Ensuring segmented continuity.

[0030] A key design feature of this model is the continuity of the two function segments at the inflection point r = z; as shown in the above equation, when r = z: (3) This value is exactly equal to the value of equation (2) at r = z; therefore, the entire function Σ(r) is continuous at r = z. This continuity constraint has a dual significance: firstly, from a physical intuition perspective, the surface density distribution of a star cluster should be a continuous function, without any density jumps; secondly, from the perspective of parameter estimation, continuity reduces one degree of freedom—if the two segments are defined independently, six parameters are required, while this model reduces the inner region amplitude from Σ through the continuity condition. King (z) is determined, the total number of parameters is reduced to five, and the redundancy of the solution space is reduced.

[0031] 3. The physical motivation for the power-law behavior in the inner region.

[0032] Internal power-law form The introduction of this principle has a clear dynamic motivation. Theoretical studies have shown that the central density profile of a self-gravitating system undergoing core collapse tends to follow a power-law distribution. For globular clusters, if there is an intermediate-mass black hole or a large amount of dark matter in the central region, the distribution of stars may also deviate from the isothermal core and exhibit a spike. Observationally, high-resolution imaging of multiple globular clusters by the Hubble Space Telescope has indeed revealed the phenomenon of the surface density in the central region continuously increasing until the resolution limit, which is consistent with power-law behavior.

[0033] In this model, the power-law exponent k is a parameter that can be freely fitted. If there is no significant peak at the center of a star cluster, k will approach 0 in the fitting result. At this time, the surface density in the inner region is approximately constant, and the model automatically degenerates into a single King model with a cutoff within the turning radius. Similarly, if the fitted value of z approaches the minimum observation radius, it indicates that there is no data constraint in the power-law segment of the inner region, and the model is essentially equivalent to a single King model in the entire radial region. Therefore, the piecewise model can be regarded as an adaptive generalization of the classic King model—when the data supports the central peak, the model is characterized by non-zero k and reasonable z; when the data does not support it, the model is automatically simplified and will not cause overfitting.

[0034] 4. Interpretability of model parameters.

[0035] The physical / geometric meanings of the five parameters in this model are clear: k: Central density peak intensity, k = 0 indicates no peak, k≈0.5-1.0 corresponds to moderate to significant peak, k>1.5 corresponds to extremely steep peak.

[0036] a: Overall density scale, related to the total luminosity / total mass of the star cluster.

[0037] Core radius: The characteristic scale of the transition of the marker outline from the inner area to the outer area.

[0038] Tidal radius: marks the outer boundary of a star cluster's gravitational bond.

[0039] z: The ending radius of the inner power-law segment, which can be interpreted as the "peak scale", that is, the spatial range of the central density anomaly extension.

[0040] The estimated values ​​of these parameters can be directly correlated with the physical state of the star cluster, facilitating subsequent statistical analysis and physical interpretation.

[0041] II. Parameter estimation methods.

[0042] 1. Initial values ​​are obtained through maximum likelihood estimation.

[0043] MCMC sampling requires a reasonable starting position to ensure that the chain can quickly enter the high-probability region. This invention uses maximum likelihood estimation as a method to obtain initial values.

[0044] Suppose the observation data contains N radial intervals, and the radial distance of the i-th interval is... The measured areal density is The number of stars in this interval is The measurement error of surface density is mainly determined by Poisson counting statistics, and can be approximated by a Gaussian distribution when the star count is sufficiently large; this invention adopts the following error model: (4) This formula assumes that the surface density error is proportional to the surface density itself, with the proportionality constant being the square root of the reciprocal of the star count. This approximation... A larger value is reasonable.

[0045] Assuming that the errors at each data point are independent and follow a Gaussian distribution, the log-likelihood function is: (5) in: This represents the predicted value of the piecewise model at the radius corresponding to the i-th interval.

[0046] To maximize The objective is equivalent to minimizing the negative log-likelihood. This invention calls the `scipy.optimize.minimize` function and uses the L-BFGS-B algorithm for numerical optimization. The optimization starting point is formed by a user-provided initial value guess plus a small-amplitude random perturbation. The purpose of introducing random perturbation is to avoid the optimization getting trapped in local extrema. After the optimization converges, the maximum likelihood estimate θ is obtained. MLE This estimated value will serve as the initial center for MCMC sampling.

[0047] 2. Bayesian framework and posterior sampling.

[0048] The core of Bayesian inference is calculating the posterior probability distribution of the parameters. According to Bayes' theorem: (6) in: That is, the likelihood function. The prior distribution is the distribution of the parameters; the posterior distribution integrates data information and prior knowledge, and is a complete description of parameter estimation and uncertainty quantification.

[0049] For multi-parameter nonlinear models, the posterior distribution usually has no analytical form; the Markov chain Monte Carlo method constructs a chain in the parameter space... For a stationary Markov chain, various statistics of the posterior distribution are approximately obtained from the samples of the chain. This invention uses an affine invariant stretching algorithm implemented in the emcee software package (Foreman-Mackey et al. 2013). This algorithm has the following advantages: it is adaptive to the linear correlation between parameters; it requires only a small number of parameter adjustments; and it is applicable to low-dimensional to medium-dimensional parameter spaces.

[0050] 3. Hybrid prior constraint mechanism.

[0051] Prior distribution The choice of parameters has a significant impact on Bayesian inference. This invention sets independent uniform prior ranges for each of the five parameters. Within the prior scope, It is a constant; if a hard boundary uniform prior is directly adopted, the acceptance probability will be zero when a chain tries to step out of the preset range, which may cause the chain to be repeatedly rejected near the boundary, reducing the sampling efficiency.

[0052] To address this issue, this invention designs a hybrid prior constraint mechanism. In the implementation of the logarithmic probability function, it first checks whether the parameters are within the range. If they are within the range, it directly returns the sum of the uniform prior logarithmic probability and the log-likelihood. If they are outside the range, the program does not return -∞ to reject the step, but instead executes an alternative branch: it calculates the normal logarithmic probability density of each parameter with the midpoint of the preset range as the mean and one-five-thousandth of the range width as the standard deviation, and sums the logarithmic probabilities of each parameter as the prior logarithmic probability.

[0053] The effects of this hybrid mechanism are as follows: within the preset range, the prior is a uniform distribution, without introducing additional parameter preferences; outside the preset range, the prior is an extremely wide Gaussian distribution, allowing the chain to occasionally step out of the range and have a chance to return; since the probability of the tail of the Gaussian distribution is much lower than the uniform probability within the range, the vast majority of samples are still concentrated within the preset range in actual operation, and the preset range still plays a dominant constraining role. This "soft boundary" design effectively avoids the sampling efficiency problem caused by hard boundaries.

[0054] 4. MCMC sampling setup and post-processing.

[0055] The MCMC sampling parameters are set as follows: the number of chains is set to 64, the number of sampling steps is set to 10000, the combustion period steps are set to 1000, and the refinement interval is set to 15. The initial position is generated by adding random perturbations near the maximum likelihood estimate, so that each chain starts from a different position in the high probability region, thus accelerating the mixing of chains.

[0056] In this design, 64 chains start simultaneously from different random perturbation positions near the maximum likelihood estimate, which is equivalent to deploying 64 independent detectors in the high-probability region of the parameter space. Compared with the cascaded single-chain scheme, this design greatly reduces the risk of a chain getting stuck in a local extremum due to poor initial position. The sampling number of 10,000 steps is much greater than the general autocorrelation length requirement of the five-dimensional parameter space. This ensures that each chain can fully traverse its local posterior region after the burn-out period, effectively eliminating sequence correlation between samples and making the final flattened sample set have the statistical characteristics of independent and identically distributed sampling. The first 1,000 steps are explicitly marked as the burn-out period and discarded, ensuring that the statistical analysis objects are strictly composed of steady-state posterior samples through simple operation. At the same time, by taking a sample every 15 steps, the autocorrelation coefficient of the final sample set is significantly reduced while retaining enough effective samples for accurate percentile calculation, making the subsequent confidence interval calculation and distribution visualization more reliable.

[0057] After sampling was completed and the combustion period was discarded, and the data was refined by interval, a flattened posterior sample array was obtained. The 50th percentile of each parameter was calculated as the median estimate, and the 16th and 84th percentiles were used as the lower and upper error bounds, respectively, corresponding to the 68% confidence interval.

[0058] The interval formed by the 16th and 84th percentiles accurately captures the asymmetric tail characteristics of this distribution, and the upper and lower error bounds provided can be numerically different, thus truly reflecting the directional differences in parameter uncertainty. At the same time, the selection of the 68% confidence interval corresponds to the probability of the interval of one standard deviation of the normal distribution, which has an intuitive statistical interpretation. The 16th and 84th percentiles precisely enclose the boundary of the highest density region of 68% probability quality of the posterior sample center, allowing researchers to clearly state that "based on the current data and model, the true value of the parameter has a 68% probability of falling within this interval."

[0059] 5. Visual evaluation of goodness of fit.

[0060] To visually assess the model fit quality, the program generates two combined plots: a density profile master plot (plotting observed data points and error bars on logarithmic scales, the best-fit piecewise model curve, the 68% confidence interval filled band, and core radius markers) and a residual plot (logarithmic residual Δ = log...). 10 Σ model -log 10 Σ obs The residuals are randomly distributed near the zero line without a clear trend, indicating a good model fit. In addition to the single-parameter marginal distribution, the corner package was used to plot the two-dimensional joint distribution corner plots between each pair of parameters.

[0061] III. Program Implementation and Architecture Design.

[0062] 1. Dependencies and runtime environment.

[0063] This method is implemented using the Python 3 programming language and relies on the following open-source scientific computing libraries: numpy (array operations), scipy (numerical optimization and statistical distribution), matplotlib (plotting), pandas (data reading), emcee (MCMC sampler), astropy.modeling.models (King model), and corner (corner plotting). The program can run on both Unix-like and Windows environments.

[0064] 2. Centralized parameter configuration architecture.

[0065] To improve the usability and portability of the program, this invention designs a centralized parameter configuration architecture; all user-adjustable parameters are defined in an independent configuration area at the top of the source code, including the following categories: (1) Basic information of the star cluster: CLUSTER_NAME (star cluster name), DATA_FILE (data file path), FIND_RA and FIND_DEC (center coordinates, for recording only).

[0066] (2) Initial value guesses for the model: k_true, a_true, rc_true, rt_true, z_true.

[0067] (3) Parameter perturbation step size: steps list.

[0068] (4) Prior ranges: k_range, a_range, rc_range, rt_range, z_range.

[0069] (5) MCMC sampling parameters: N_WALKERS (number of chains), N_STEPS (number of sampling steps), DISCARD (number of combustion phase steps), THIN (refinement interval).

[0070] (6) Output settings: OUTPUT_DIR (output directory), OUTPUT_PREFIX (file name prefix).

[0071] Users only need to modify the above variable values ​​according to the star cluster to be analyzed; no other code needs to be changed.

[0072] 3. Program flow and output.

[0073] The main program flow includes: importing dependent libraries and global plotting settings, reading configuration parameters, checking data files, reading data and calculating errors, defining piecewise model functions, solving maximum likelihood estimation, mixed prior Bayesian sampling, posterior sample processing and result output, drawing density profile master plot and residual plot, and drawing parameter posterior joint distribution angle plot.

[0074] The program outputs the following after running: the maximum likelihood estimate, prior range, and MCMC fitting results are printed to the console; density profile combination plot (PNG format); and parametric posterior angle plot (PNG format).

[0075] IV. Scope of application and parameter settings of the method of the present invention.

[0076] 1. Data requirements.

[0077] The method of this invention has the following basic requirements for the input areal density data: (1) Data format: The input file should be a comma-separated text format, containing at least three columns of data: radial distance (angular distance from the center of the star cluster), surface density value (star count per unit area or surface brightness), and star count within the radial interval (for error estimation).

[0078] (2) Radial coverage: Data points should cover a sufficient radial range from the central region of the star cluster to the outer region. Ideally, it should span at least an order of magnitude of radius to ensure that there are sufficient data constraints in both the inner and outer regions of the model.

[0079] (3) Number of data points: It is recommended that there be no less than 10 effective data points. Too few data points may lead to insufficient constraints of the five-parameter model and excessive broadening of the posterior distribution.

[0080] (4) Error estimation: The surface density error can be approximated by Poisson statistics, i.e. ,in D For surface density, N Count the stars within this interval; given a sufficiently large number of stars (usually N>20), the Gaussian error is approximately reasonable.

[0081] 2. Suggestions for initial parameter settings.

[0082] The initial values ​​of the parameters in the user configuration area are used to construct the starting point for the maximum likelihood estimation. Appropriate initial values ​​can accelerate optimization convergence and avoid getting trapped in local optima. The following are general suggestions for setting the initial values ​​of each parameter: (1) Power law exponent k_true: If there is no prior information, it is recommended to set it to a value between 0.5 and 0.8. k=0 corresponds to no central peak, and k>1 corresponds to a strong peak.

[0083] (2) Center amplitude a_true: can be set according to the area density of the data center area, usually set to 1-2 times the density of the data center.

[0084] (3) Core radius rc_true: This can be determined by visually inspecting the density profile, specifically the radial distance where the areal density drops to approximately half of the center value. Initial value guess.

[0085] (4) Tidal radius rt_true: can be set to 2-5 times the maximum radial distance of the data, or refer to the typical value of the same type of star cluster.

[0086] (5) Turning radius z_true: can be set to 1-3 times the minimum radial distance of the data, or 0.1-0.3 times rc_true.

[0087] 3. Principles for setting prior scope.

[0088] The prior range defines the allowed range of values ​​for the parameter, and the setting principles are as follows: (1) k_range: The power law exponent is usually set in the range of (0, 2). The lower limit is usually set to 0, and the upper limit is set to 2 to cover the peak intensity that most actual star clusters may have.

[0089] (2) a_range: The range of the amplitude parameter should cover the order of magnitude of the data surface density. The lower limit can be set to 0.1 times the minimum surface density, and the upper limit can be set to 10 times the maximum surface density.

[0090] (3) rc_range: The range of the core radius should include the radial interval where the contour transitions from flat to sloping. Typically, the lower limit is not less than 0.5 times the minimum data radius, and the upper limit is not greater than the maximum data radius.

[0091] (4) rt_range: The tidal radius is usually greater than the maximum data radius. The lower limit can be set to 0.5 times the maximum data radius, and the upper limit can be set to 3-5 times the maximum data radius.

[0092] (5) z_range: The turning radius should be between the minimum data radius and Between. The lower limit can be set to 0 or 0.1 times the minimum data radius, and the upper limit can be set to 1-2 times rc_true.

[0093] If prior knowledge of a certain parameter is lacking, a wider range can be set, allowing the data to dominate the posterior constraints.

[0094] 4. Adjust MCMC sampling parameters.

[0095] (1) N_WALKERS (number of chains): Usually set to 10-20 times the parameter dimension. For a five-parameter model, 32-64 chains are a reasonable choice.

[0096] (2) N_STEPS (number of sampling steps): For a five-dimensional unimodal posterior, 5000-10000 steps are usually sufficient for the chain to reach a steady state. The convergence can be determined by checking the chain trajectory diagram.

[0097] (3) DISCARD (burning period steps): Generally, it is 10%-20% of the total number of steps.

[0098] (4) THIN (refinement interval): usually 10-20, to control storage requirements while retaining a sufficient number of effective samples.

[0099] 5. Interpretation of the output results.

[0100] After the program runs, it outputs the median values ​​of each parameter and the 68% confidence interval; in the density profile plot, the red dashed line represents the inner power-law segment, the blue dashed line represents the outer King segment, the dark blue filled band represents the 68% confidence interval, and the vertical dashed line marks... Median; Residual plots are used to assess goodness of fit—residuals that are symmetrically distributed near the zero line and have no obvious trend indicate a good fit, and corner plots can be used to assess the correlation between parameters and the posterior distribution shape.

[0101] V. Advantages of the present invention.

[0102] 1. The physical meaning and applicable scenarios of the segmented model.

[0103] The piecewise power-law King model proposed in this invention is designed for the surface density profile of star clusters with a central density excess characteristic. The inner power-law segment can describe the steep density rise in the central region, a morphology commonly seen in globular clusters undergoing core collapse, clusters with intermediate-mass black holes at the center, or clusters exhibiting significant mass stratification effects. The outer King segment retains the mature framework of classical models describing tidal truncation, where r ≥ r t The surface density naturally drops to zero.

[0104] For star clusters with no density excess at the center, the piecewise model can automatically degenerate into a single King model through parameter adaptation: when the fit yields k ≈ 0 or z ≈ r min When the inner power-law segment tends to be flat or has a very small influence range, the model is equivalent to the King model in the entire radial region. This characteristic makes the method of the present invention widely applicable, without the need to pre-determine the cluster type.

[0105] 2. Technical advantages of hybrid prior mechanisms.

[0106] When using a hard-boundary uniform prior in traditional MCMC sampling, the acceptance probability immediately drops to zero once the parameters step out of the preset range, causing the chain to be repeatedly rejected near the boundary, reducing sampling efficiency, and even causing posterior estimation bias at the boundary.

[0107] The hybrid prior mechanism adopted in this invention has the following advantages by introducing a wide Gaussian soft constraint instead of a hard boundary rejection: (1) it maintains a uniform prior within a preset range without introducing additional parameter preferences; (2) when the boundary is exceeded, the parameters can still obtain a limited log probability, allowing the chain to have a chance to return from outside the boundary; (3) the standard deviation of the Gaussian distribution is taken as 1 / 5000 of the range width, which is wide enough to not significantly affect the uniformity within the range, while being narrow enough to maintain effective constraints on the parameters; this design improves the sampling behavior near the boundary and enhances the hybrid efficiency of the chain.

[0108] 3. The engineering value of centralized parameter configuration architecture.

[0109] In traditional scientific code, model parameters, file paths, sampling settings, etc., are scattered throughout the code, requiring users to search and modify them line by line, which is error-prone and inefficient. This method adopts a centralized parameter configuration architecture, encapsulating all user-adjustable parameters in an independent configuration area at the top of the source code, thus decoupling program logic from input parameters. This design allows non-programming astronomy researchers to apply the program to different star cluster data without needing to understand the core algorithm, improving code portability and reusability, and facilitating batch processing and version management.

[0110] This invention defines a segmented analytical model at the inflection radius z—the inner region uses a power-law descent form to describe the possible density peak at the center, while the outer region uses the classic King model to describe the tidal truncation; the two segments are continuous at the inflection point. The model contains five free parameters: the power-law exponent k, the central amplitude a, the core radius, and the... tidal radius And the turning radius z; This invention adopts a two-stage strategy. The first stage obtains the initial values ​​of the parameters through maximum likelihood estimation; the second stage obtains the posterior distribution of the parameters through MCMC Bayesian sampling, using the median as the best estimate and the 16th and 84th percentiles as the 68% confidence interval; This invention designs a hybrid prior mechanism that uses a uniform prior within a preset range and introduces a wide Gaussian soft constraint when the range is exceeded, effectively avoiding the MCMC sampling efficiency problem caused by hard boundaries; This invention adopts a centralized parameter configuration architecture, encapsulating all user-adjustable parameters in an independent configuration area at the top of the code. Users only need to modify the variables in the configuration area to apply the program to different star clusters, which has strong portability and ease of use.

[0111] This invention proposes an inner-region power-law-outer-region King segmented combination model, which divides the radial surface density profile of a star cluster into two analytical representation regions by using a turning radius parameter, solving the technical problem that a single King model cannot simultaneously describe the central density peak and the outer-region tidal truncation. This invention designs a hybrid prior constraint mechanism, employing a uniform prior for parameters within the boundary during Bayesian sampling and introducing a wide Gaussian soft constraint for parameters outside the boundary, avoiding Markov chain stagnation caused by hard boundaries and improving parameter space traversal efficiency. This invention adopts a centralized parameter configuration architecture, encapsulating all variable parameters in an independent configuration area at the top of the source code, achieving decoupling between the program and data, allowing non-programming astronomy researchers to reuse the code without modifying the core algorithm. This invention quantifies the median estimate and confidence interval of the fitted parameters through full-parameter posterior sampling and visualization, providing a visual corner plot of the joint distribution between parameters, overcoming the limitations of traditional methods that only provide point estimates and linear error approximations.

[0112] The newly added power-law exponent k and turning radius z in this invention provide direct, cross-sample comparable quantitative indicators of the dynamic state of the central region of star clusters. Researchers can use these new parameters to conduct large-sample statistics, such as exploring the correlation between the power-law exponent and the age of star clusters, the orbital parameters of the Milky Way, and the existence of X-ray sources. This transforms the problem, which could only be qualitatively described ("there is" or "there is no" spike), into a quantitative scientific problem in a continuous parameter space, significantly improving the depth of scientific data mining.

[0113] In the hybrid prior mechanism of this invention, the standard deviation of the wide Gaussian is set to one-thousandth of the boundary width. This refined design has a profound balancing effect: within the allowable physical range, the uniform prior dominates, ensuring the objectivity of data-driven inference; while outside the boundary, the extremely wide Gaussian distribution provides an extremely low background probability, which can prevent parameters from escaping to completely non-physical numerical regions without exerting any substantial pressure on the posterior form within the boundary. This achieves a "strong guidance, weak constraint" prior strategy, which is more robust and intelligent than simple hard boundaries or purely weak information wide priors.

[0114] The centralized parameter configuration architecture of this invention converges all user interaction interfaces into an independent block at the top of the code. It also completely separates the roles of "model user" and "algorithm maintainer," achieving a user-friendly "black box" encapsulation for the transformation of scientific research results. Users do not need to understand the internal principles of the MCMC sampler, the details of maximum likelihood optimization, or the logic of piecewise function construction. They only need to modify the variables in the configuration area as if filling out a form to reuse the same high-precision analysis process for hundreds or thousands of different star cluster data. This represents a qualitative improvement in portability, reproducibility, and batch processing efficiency compared to existing technologies.

[0115] The embodiments described above are merely illustrative of several implementations of the present invention, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these all fall within the protection scope of the present invention. Therefore, the protection scope of this invention patent should be determined by the appended claims.

Claims

1. A method for constructing an analytical model to fit the density profile of a globular cluster's projected surface, characterized in that, Comprising the following steps: Obtain the observational data of the globular cluster, and preset a turning radius parameter to divide the radial profile of the globular cluster in the observational data into two independently analytically described regions, namely the inner region and the outer region. In the inner region, use a power-law decline function form to describe the density spike in the central region of the cluster, and in the outer region, use the classical King projection analytical model to describe the tidal truncation characteristics of the cluster, so as to construct a piecewise analytical function for describing the variation of the projected surface density of the globular cluster with the radial distance, and obtain multiple independent parameters of the piecewise analytical function; Based on the observational data, use the maximum likelihood estimation method to obtain the initial estimated values of each independent parameter; Based on the initial estimated values of each independent parameter, perform Markov chain Monte Carlo sampling on each independent parameter respectively within the Bayesian inference framework, and adopt a mixed prior constraint mechanism in the logarithmic probability function during sampling: for the independent parameters within the preset uniform prior range, directly obtain the sum of the uniform prior logarithmic probability and the log-likelihood of the independent parameter as the prior logarithmic probability, and for the independent parameters beyond the preset range, obtain the logarithmic probability density of each independent parameter based on a Gaussian distribution with a preset width, and sum the logarithmic probability densities of each independent parameter as the prior logarithmic probability; Statistically analyze the posterior samples obtained by sampling to obtain the median estimated values and confidence intervals of each parameter, so as to construct a piecewise analytical model for fitting the projected surface density profile of the globular cluster, and generate the projected surface density profile of the globular cluster after fitting.

2. The analytical model construction method for fitting the projected surface density profile of a globular cluster according to claim 1, characterized in that, The construction of the piecewise analytical function includes: Preset a turning radius parameter z to divide the radial profile of the globular cluster into an inner region and an outer region, and set the projected surface density of the globular cluster as a function Σ(r) of the radial distance r, where r is the angular distance from the center of the cluster, then the piecewise analytical function is expressed as: When r < z, the piecewise analytical function is expressed as: ； When r ≥ z, the piecewise analytical function is expressed as: ； in: a This represents the center surface density amplitude parameter; Indicates the core radius; The tidal radius is represented by k; the power-law exponent of the inner zone is represented by k. z Indicates the radius of inflection; This represents the King projection analytical model.

3. The analytical model construction method for fitting the projection surface density profile of a globular cluster according to claim 2, characterized in that, The obtaining of the initial estimated values of each independent parameter includes: Suppose the observation data contains N radial intervals, and the radial distance of the i-th interval is... The measured areal density is The number of stars in this interval is Assuming the areal density error at each point in the observed data follows a Gaussian distribution, the log-likelihood function is constructed as follows: ； in: ; This represents the predicted value of the piecewise analytic function at the radius corresponding to the i-th interval; By maximizing the constructed log-likelihood function using a numerical optimization algorithm, initial estimates θ of a set of parameters are obtained. MLE .

4. The analytical model construction method for fitting the projected surface density profile of a globular cluster according to claim 3, characterized in that, The sampling process of performing Markov chain Monte Carlo sampling on each independent parameter respectively includes: With an initial estimate of a set of parameters θ MLE Based on this, Markov chain Monte Carlo sampling is performed on each independent parameter within the Bayesian inference framework. In the calculation of the log probability function of the Markov chain Monte Carlo sampling, it is determined whether each independent parameter lies within its independently defined uniform prior range. ]Inside; When the independent parameters are within the set uniform prior range [ Within this range, the sum of the uniform prior log probability and the log-likelihood of the independent parameters is directly obtained as the prior log probability. ; When the independent parameters are within the set uniform prior range [ Outside of this range, obtain the normal log probability density of each independent parameter with the midpoint of the preset range as the mean and one-five-thousandth of the range width as the standard deviation, and sum the log probabilities of each independent parameter as the prior log probability. .

5. The analytical model construction method for fitting the projection surface density profile of a globular cluster according to claim 4, characterized in that, The generation of the projected surface density profile of the globular cluster after fitting includes: Statistical analysis is performed on the posterior samples obtained from the sampling to obtain the sample set of each independent parameter. The preset first percentile is used as the median estimate, and the preset second and third percentiles are used as the lower and upper error bounds, respectively, to quantify the confidence interval of the fitting result; Based on the selected median estimated values and confidence intervals, use the set of independent parameter samples to generate the piecewise surface density profile curve of the globular cluster after fitting and its confidence interval filling band.

6. The analytical model construction method for fitting the projected surface density profile of a globular cluster according to claim 5, characterized in that, The power-law exponent k is used to quantify the steepness of the density spike in the center of the globular cluster; When k > 0, it indicates that there is a density spike in the center of the globular cluster, and the larger the k value, the steeper the spike; When k approaches 0, the piecewise analytical model degenerates into the truncation of the classical King projection analytical model within the turning radius z.

7. The analytical model construction method for fitting the projection surface density profile of a globular cluster according to claim 1, characterized in that, The Markov chain Monte Carlo sampling adopts an affine invariant stretching movement algorithm and configures multiple parallel chains, and its sampling parameters include the number of chains, the number of sampling steps, the number of burn-in steps, and the thinning interval.

8. An analytical model construction device for fitting the projection surface density profile of a globular cluster, characterized in that, Including: The piecewise analytical module is used to acquire observational data of globular clusters and pre-set a turning radius parameter to divide the radial contour of the globular cluster within the observational data into two independently analytically described regions: an inner region and an outer region. In the inner region, a power-law descent function is used to describe the density peak in the central region of the cluster, while in the outer region, the classical King projection analytical model is used to describe the tidal truncation characteristics of the cluster. This is to construct a piecewise analytical function to describe the variation of the density of the projected surface of the globular cluster with radial distance and to acquire multiple independent parameters of the piecewise analytical function. The analytical model building and fitting module is used to obtain initial estimates of each independent parameter based on the observed data using the maximum likelihood estimation method. Based on the initial estimates of each independent parameter, Markov chain Monte Carlo sampling is performed on each independent parameter under the Bayesian inference framework. A hybrid prior constraint mechanism is adopted in the log probability function during sampling: for independent parameters within a preset uniform prior range, the sum of the uniform prior log probability and the log likelihood of the independent parameter is directly obtained as the prior log probability; for independent parameters exceeding the preset range, the log probability density of each independent parameter based on a preset width Gaussian distribution is obtained, and the log probability density of each independent parameter is summed as the prior log probability. The posterior samples obtained from the sampling are statistically analyzed to obtain the median estimate and confidence interval of each parameter, so as to construct a piecewise analytical model for fitting the projection surface density profile of the globular cluster and generate the fitted projection surface density profile of the globular cluster.

9. An electronic device, characterized in that, include: Memory and processor; The memory is used to store computer programs; When the processor executes the computer program stored in the memory, it implements the steps of the analytical model construction method for fitting the projection surface density profile of a globular cluster as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, Used to store a computer program, which, when executed by a processor, implements the steps of an analytical model construction method for fitting the projection surface density profile of a globular cluster as described in any one of claims 1 to 7.

Images