Monte-Carlo Stellar Dynamics near Massive Black Holes Two-dimensional Fokker-Planck solutions of multiple mass components

,


INTRODUCTION
It is widely believed that massive black holes (MBHs) reside in the center of most galaxies, embedded by dense clusters consisting of stars and compact objects (the socalled nucleus star clusters).The deep potential of the MBH significantly changes the stellar dynamics, and the strong interplay between the stars and the MBH leads to a number of violent phenomena in galactic nuclei.
Many of these phenomena are tightly connected to the evolution of MBHs and the star formation history of galactic nuclei, e.g., the tidal disruption of stars (e.g., Rees 1988), ejection of hyper-velocity stars (e.g., Hills 1988), and stellar collisions (e.g., Duncan & Shapiro 1983;Quinlan & Shapiro 1990).
In order to precisely reveal the characteristics of these events, it is of fundamental importance to understand the stellar dynamics that drive the evolution of the or-bits of stars and compact remnants in the vicinity of MBHs.
The stellar dynamics around the MBHs are driven by a number of different relaxation processes, e.g., the twobody relaxation, resonant relaxation (RR) (Rauch & Tremaine 1996;Hopman & Alexander 2006), and other dynamical effects, e.g., the tidal field of the MBH, relativistic orbital precession, gravitational wave orbital decay.
In the case of multiple mass components, mass segregation can drive the lighter components outward and the heavier ones migrate in, which is important for a number of investigations of the properties of the cluster, such as density profiles (e.g., Bahcall & Wolf 1977;Alexander & Hopman 2009; Preto & Amaro-Seoane 2010) and rates of EMRIs (e.g., Hopman & Alexander 2006;Amaro-Seoane & Preto 2011).If the particles are binaries, the dynamics can be even more complex, e.g., the ionization of binaries, exchanges of binary components, Kozai-Lidov effects (e.g., Hopman 2009;Zhang et al. 2019).These complicated effects, which cover from small to large scales of the cluster, impose difficulty in obtaining accurate pictures of the evolution and outcomes of particles of different types and masses near the MBH.
The two-body relaxation process, which is essential for stellar dynamics near the MBH, has already been significantly investigated over the past several decades.Most of them adopt the one-dimensional Fokker-Planck (FP) analytical method (hereafter 1D-FP method), which mainly considers only the evolution in energy under the assumption of an isotropic distribution of angular momentum (Lightman & Shapiro 1977;Bahcall & Wolf 1976, 1977;Bar-Or et al. 2013;Hopman & Alexander 2005;Alexander & Hopman 2009;Keshet et al. 2009).
As the relaxation on angular momentum is usually much faster than that in energies (unless in a nearcircular orbit), its evolution can be considered approximately independent, and the effects of the loss cone can be approximately included in the evolution of energy (Lightman & Shapiro 1977;Bahcall & Wolf 1977).Usually, the orbits of particles can be considered approximately Keplerian near the MBH, and the effects of the stellar potential are important only at the outer parts of the cluster (i.e., outside the influence radius of the MBH) (Marchant & Shapiro 1979).
The numerical methods that consider both the evolution of energy and angular momentum in the FP equations have also been developed, either by Monte-Carlo methods (Shapiro & Marchant 1978;Bar-Or & Alexander 2016), finite differential methods (Cohn & Kulsrud 1978;Takahashi 1995), or finite element methods (Takahashi 1995).However, these previous studies are limited to only one mass component (assuming equal masses for all particles) and thus cannot handle the process of mass segregation if there are multiple mass components in the cluster.Takahashi (1997) extends the finite differential/element methods to include multiple mass components, although we notice that it is difficult to expand further on the method by including other dynamical effects, such as stellar collisions and those of the binaries.
Another approach is to use direct N-body numerical simulations (e.g., Antonini 2014;Merritt et al. 2010;Preto et al. 2004;Panamarev et al. 2019;Baumgardt et al. 2004), most of which adopt the N-body code series (Aarseth 1999).However, it is usually expensive to integrate the orbits of stars for more than several relaxation timescales, and sometimes it is necessary to use specialized hardware such as GRAPE (Preto et al. 2004;Baumgardt et al. 2004;Antonini 2014).Limited by the resolution and the number of particles in the simulation, it is challenging to simulate scales below 10 −3 ∼ 10 −2 of the gravitational influence radius of the MBH (Panamarev et al. 2019).Also, the N-body simulations are difficult to include various processes that happen at small scales as mentioned earlier.
Another important approach is the Monte-Carlo method based on the H'enon scheme (Hénon 1961;Freitag & Benz 2001).The relaxation in such a method can be considered by shell-like particles, and the density of particles is updated after each time step.The method is flexible to include various dynamical effects and multiple mass components, which has already been applied for simulations of globular clusters (e.g., Joshi et al. 2000) and galactic nuclei (Freitag & Benz 2001, 2002;Freitag et al. 2006).
A number of interesting dynamical phenomena, e.g., tidal disruption events, stellar collisions, and different kinds of gravitational wave sources such as EMRIs, are results of particles evolving in both energy and angular momentum and embedded in a cluster consisting of multiple mass components.Thus, it is necessary to develop a Monte-Carlo method to simulate the evolution of particles in such an environment with acceptable accuracies and numerical costs.Here we develop a different Monte-Carlo approach (GNC)1 based on Shapiro & Marchant (1978), which is extended such that, for the first time in the literature, it can obtain the dynamical evolutions of objects with multiple masses based on twodimensional (both energy and angular momentum) FP equations.The method is also an overhaul of the pre-vious version of the code in (Zhang et al. 2019(Zhang et al. , 2021)).One of the advantages of GNC is that it is very flexible in including various kinds of stellar objects and different dynamical effects mentioned above.GNC can obtain statistically accurate results for rare objects by varying the weightings of particles.Our code can be applied to a number of important dynamical phenomena in galactic nuclei in the future, providing a more consistent picture of dynamics for particles of various types and across a spectrum of masses.
This paper is organized as follows.In Section 2, we introduce the core method of GNC, including the Monte-Carlo method of the two-body relaxation, the boundary conditions, the weighting of particles, and the general Monte-Carlo routines in obtaining steady-state solutions.In Section 3, we compare our simulation results with those of previous studies for models consisting of one or multiple mass components.We test the relaxation processes, the density profiles, and the consumption rates of particles in the loss cone.Discussion and conclusions of the work will be presented in Section 4.

THE METHOD
GNC is a Monte Carlo code that can calculate the evolution of particles in energy and angular momentum spaces when a central massive black hole (MBH) is present.In our previous works Zhang et al. (2019) and Zhang et al. (2021), we have considered the evolution of binary black holes in a fixed background of single stars or black holes.In this work, we have made many improvements to GNC to transform it into a general tool for obtaining steady-state solutions of multiple mass components.It is important to note that we only provide a description of the basic version of GNC in this context, and we do not incorporate other dynamical effects (e.g., gravitational wave orbital dissipation) apart from the two-body relaxation and Newtonian loss cone.Further details of the method are presented in the following sections.

Two body FP evolution of multiple mass components
Per definition, the MBH dominates the dynamics of the cluster within the gravitational influence radius r h .According to the M • − σ relation from Kormendy & Ho (2013), where σ h is the velocity dispersion of the galaxy, the gravitational influence radius r h is given by, where r 0 = 3.1 pc is the influence radius for the MBH in Milky Way.This radius is roughly consistent with the breaking radius measured about ∼ 3pc in our Galactic center (Schödel et al. 2018).Then the density at r h , i.e., n h , is given by
If the cluster contains multiple mass components, the two body relaxation on energy E and angular momentum J of each of the components follows the FP differential equations that is coupled with those of other components.Here we define E = GM • /(2a 2 ) and 2 ), where a 2 and e 2 are the semimajor axis (sma) and eccentricity of the outer orbit of the particle circling MBH.For the α-th component, where N m is the total number of mass components, it is destribed by ) where N α (E, J) is the number distributions of particles of the α-th mass component in E − J spaces.
The diffusion coefficients D E , D J describe the drift, and D EE and D JJ describe the scatterings of E and J. D EJ describe the correlation of scatterings between E and J.These diffusion coefficients are coupled with the distribution functions in other mass bins.Specifically, for the α-th component with mass m α , the drift terms are given by the diffuse and cross terms are given by where β runs over all mass components.j = J/J c is the dimensionless angular momentum and is the angular momentum for circular orbits.Γ β functions can be found in Bar-Or & Alexander (2016) or Zhang et al. (2019), which depends on the dimensionless function ḡβ (E) of the β-th component, and it will be explained later.
The relation between the number density N α (E, J) and the distribution function f α (E, J) in phase spaces is given by where P (E) = 2πM • /(2E) 3/2 is the orbital period around a MBH with mass M • , V and v is the 3dimensional position and velocity, respectively.When discussing the dynamics and distributions of particles, it is more convenient to adopt the following dimensionless form: where σ ⋆ = σ h , n ⋆ = n h is the velocity dispersion and density of the reference star at r h , respectively.The dimensionless number distribution N α (x, j) in the α-th mass bin is then given by N According to Equation 8, we have We can now define ḡα (x), the j-averaged dimensionless parameter by ḡα (x) = 2 1 0 g α (x, j)jdj. (10) Then ḡα (x) can then be used for calculations of the Γ α functions in the corresponding mass bin.

Boundary conditions
Particles will be destroyed if they approach too close to the MBH.Therefore, it is natural to impose a vanishing boundary condition in the inner part of the cluster.In this case, we set x max = 10 5 as the inner boundary, which corresponds to a distance of approximately ∼ 3.2 AU for an MBH with M • = 4 × 10 6 M ⊙ .
On the other hand, if a particle moves too far away from the influence radius, it may become unbound from the MBH.Ideally, the boundary should be set at x = 0.However, approaching this value numerically is not feasible, so we set the outer boundary at x = x B = 0.05, which is smaller than the value used by Shapiro & Marchant (1978) (they set x B = 0.2).Particles with energies larger than this value are considered unbound from the MBH.The unbound stellar populations should be asymptotically similar to those of the bulge stars.Therefore, it is reasonable to set a fixed given number density (or number ratio) for different particle types, which corresponds to the Dirichlet boundary condition at x B for Equation 4.
Suppose there are multiple mass components for x < 0 (the unbound population), denoted by m α = 1, • • • , N m , and each of them consists of β = 1, • • • , N tα different types of particles.For example, we can set β = 1 for stars, 2 for stellar-mass black holes (SBHs), 3 for neutron stars (NSs), 4 for white dwarfs (WDs), and 5 for brown dwarfs (BDs).Let s β (m α ) represent the number density ratio of the α-th mass component of type β relative to the total number of stars at the outer boundary.For stars, we have s β (m ⋆ )dm β = 1.If we assume that all stars have equal masses, then we simply have s β (m ⋆ ) = 1.
Assume violent relaxation boundary conditions, the velocity dispersion of all components is the same as those of the reference star (Alexander & Hopman 2009).For the objects with type β in the α-th mass bin, the distribution function g αβ (x, j), x ≤ x B , at the boundary is given by As we will discuss later (in Section 3.3), when the mass of the MBH is around M • = 10 4 − 10 8 M ⊙ , the particles, especially compact objects, slightly inside the outer boundary, follow a form closer to the full loss cone (See also Equation D20).However, if the mass of the MBH is larger than 10 8 M ⊙ , the distribution gradually transitions to the empty loss cone condition (See also Equation D21).Since we primarily focus on galaxies containing M • < 10 8 M ⊙ , which are more commonly found in the nearby universe, we simply assume a full loss cone condition for all simulations conducted in our study.Therefore, we set2 G(j) = 1.
For the isotropic boundary condition, we follow the approach of Shapiro & Marchant (1978) and assume that the unbound population (x < 0) consists of bulge stars with isothermal distributions and a Maxwellian velocity distribution.Thus, we have g αβ (x < x B , j) = s β (m α )e x for x < 0. This isotropic outer boundary is consistent with observations of our Galactic center (Feldmeier-Krause et al. 2017).
With these boundary conditions, GNC will generate a continuous flow of particles into the cluster.Further details about the Monte Carlo methods can be found in Appendix A. Figure 1 provides an illustration of the boundary methods.

Particle Weights and Particle Splitting Method
We acknowledge that there is a difficulty in the Monte Carlo simulation due to the significant variation in the number of particles in different mass bins.For instance, in a two-component model, the number ratio of single black holes (∼ 10 − 40M ⊙ ) to reference stars (∼ 1M ⊙ ) is on the order of 10 −3 (Hopman & Alexander 2006) or even smaller.Therefore, if each particle has the same weight, the number of stars will dominate over the number of single black holes.To increase the efficiency of the simulation and improve the statistics of rare particles, it is more convenient to assign different weights to particles in different mass bins.We denote this weighting factor as w α for the α-th mass bin.For example, in the case of M • = 4 × 10 6 M ⊙ , we can set w 1 = 400 for stars (or w 2 = 10 for single black holes), meaning that one particle of stars (or single black holes) in the simulation represents 400 (or 10) particles.In this way, we find that in the steady state of the simulation, the total number of single black holes is comparable to the number of stars.For particles in other mass bins and for different MBH masses, we can adjust the value of w α to ensure an adequate number of particles for accurate statistics.
Previously, Shapiro & Marchant (1978) used a particle cloning scheme to increase the statistics of particles in the inner regions.The previous version of GNC also adopted this scheme (Zhang et al. 2019).However, this scheme requires labeling particles as "original" or "clones," which increases the complexity for future expansions of the code.For example, when a "clone" star merges with an "original" star after a stellar collision, or when a "clone" single black hole becomes a component of an "original" binary system after a three-body exchange, it becomes challenging to determine whether the merged/exchanged product should be classified as an "original" particle or a "clone" particle.
To address this issue, we have improved the cloning scheme to a more general splitting method.When a particle moves inward and crosses the energy x/x 0 = 10 i , the particle splits into Π clones (or creates Π−1 clones), where i = 0, • • • , 4 and x 0 is an arbitrarily dimensionless energy above which the clone scheme is implimented.Conversely, when any particle moves upward and crosses back to those energies (regardless of whether it is a clone or an original particle), there is a chance of 1 − 1/Π that it will be deleted from the simulation.To obtain correct statistical results, we need to correct the weight of the split particle by another factor, which is given by w c = 1/Π L , where L = max[Int(log 10 (x/x 0 ) + 1), 0].The particle splitting method is also illustrated in Figure 1.We note that if Π is the same for particles in different mass bins, there will be an overabundance of massive particles concentrated in the inner parts of the cluster due to mass segregation.To address this issue, we can assign different values of Π for each mass bin and decrease Π to a lower value for mass bins with larger masses.In the models presented in this work, typically Π varies from 30 to 4 − 10 as the masses increase from ≲ 1M ⊙ to > ∼ 10 − 40M ⊙ .Suppose the "uncorrected" dimensionless distribution function at the outer boundary for an arbitrary mass bin α is given by ḡu α (x B ).It needs to be corrected such that: This correction determines a constant w n = s(m α )/ḡ u α (x B ), which can be used for the correction to obtain the true weights of all particles.It is important to note that we can select an arbitrary mass bin for this correction, as the number ratio of particles beyond the outer boundary, s β (m α ), remains unchanged during the simulation.Therefore, the value of w n obtained from the α mass bin is the same for other mass bins as well.
Finally, the true weights of any particle are determined as: Normalization is usually necessary before the 3rd or 4th iteration, after which the weighting w n remains constant as the distribution near x B converges.It is important to note that, according to the normalization method and Equation 9, the true weight of particles corresponds to the objects in reality, as well as the associated results (e.g., event rates of tidal disruptions for stars), will be proportional to n 0 r3 0 , where r 0 and n 0 are given by Equations 2 and 3, respectively.This should be kept in mind when comparing our results with those from other authors.

Steady-state solution of distribution, density and anisotropy
The Monte Carlo simulation begins with an initial guess of distributions for different mass components, denoted as N α (E, J).All particles undergo a simulation duration of δτ = 0.005 ∼ 0.01T rlx (r h ), where T rlx (r h ) is the two-body relaxation time at r h .In the case of multiple components, it is given by (Binney & Tremaine 1987): where approximately Λ ≃ M • /m ⋆ .After that, the distribution of particles ḡ(x) is updated to a new one, which is then used to calculate the diffusion coefficients for the next iteration.Repeat this process for a sufficient number of times, and the distribution will gradually converge to a steady-state solution.More details are shown in Appendix C. At any snapshot in time, the number density of particles in the mass bin α at a distance of r is then given by 3 We can also estimate the anisotropy at position r by the parameter where and is the velocity dispersion in tangential and radial directions, respectively.These values can be estimated according to

SIMULATIONS AND RESULTS
To test the results from GNC, we need to compare them to previous studies, particularly those that have used one-dimensional Fokker-Planck (1D-FP) methods.The 1D-FP method assumes an isotropic distribution of angular momentum and considers only one-dimensional energy relaxations.This method has been extensively tested in the past and has shown general consistency with results obtained from other methods.Therefore, our main focus here is to compare our results with those obtained using the 1D-FP method, as it serves as a simplified version of our approach.We provide a review of the relaxation process and the loss cone effect of the 1D-FP method in Appendix D.

The Relaxation and density profile of nuclear star cluster
To test whether GNC accurately simulates the dynamical evolution of particles around the central supermassive black hole (MBH), we first examine the relaxation and density profile of the nuclear star cluster using various test models with one or multiple mass components.The boundary conditions for different test models, ranging from 1 to 5 components, are shown in Table 1.We also discuss and analyze the results for models with a spectrum of masses in this section.We begin with a MBH mass of M • = 4 × 10 6 M ⊙ , which corresponds to the mass of the Milky Way black hole (Gillessen et al. 2009), but we vary the MBH mass later in this section.

Single-mass component
We first examine the results from model M1, which consists of stars with equal masses (m ⋆ = 1M ⊙ ).The left and right panels of Figure 2 show the evolution of the density profile and the energy distribution at different times of the simulation, respectively.We observe that the convergence time is around ∼ 0.4T rlx (r h ) for particles inside a radius of 0.1r h , where T rlx (r h ) is the two-body relaxation time at r h given by Equation 14.This timescale is consistent with the results of Preto et al. (2004), who also used N-body simulations.It is evident that the relaxation of particles becomes faster as they get closer to the central MBH.
Figure 3 compare the results of M1 from GNC and those from 1D-FP methods.The two methods are well consistent, for the energy distribution ḡ(x), the density profile n(r), and the cumulative number distribution N (< r).
The discrepancies of results between these two methods at the inner parts of the cluster are slightly larger if additionally the loss cone is considered.The difference is mostly apparent in ḡ(x), as it will amplify the difference in density profile by a factor of ∝ x 3/2 .
As we assume isotropic boundary conditions, the anisotropy of the particles is small (A ∼ 0) at the outer parts of the cluster.However, it becomes slightly tangentially anisotropic (A ∼ −0.2) in the middle of the cluster, even without the presence of a loss cone (top right panel of Figure 3).If the particles are additionally depleted by the loss cone, the tangential anisotropy is increased to A ≲ −0.6 for r < 10 −3 r h (bottom right panel of Figure 3).This small anisotropy developed in the cluster may explain the slight differences between our method and the 1D-FP methods, as the latter assumes isotropic distributions.
The tangential anisotropy observed in the inner parts of the cluster is consistent with the findings of Cohn & Kulsrud (1978), who used finite differential methods in two-dimensional E − J spaces.Observations of the Galactic center also suggest that the nuclear star cluster exhibits isotropy at outer parts but becomes tangentially anisotropic at inner parts (Feldmeier-Krause et al. 2017).
The slope index of the density profile (γ = d ln n(r)/d ln r) of stars in model M1 can be seen in the third panel of Figure 3.We observe that without the loss cone, γ ≃ −1.75 for 10 −4 ≲ /r h ≲ 0.1r h , which is the expected value for a star cusp around a MBH (Bahcall & Wolf 1976;Lightman & Shapiro 1977).The slope of stars in the inner parts of the cluster (r < 0.01r h ) becomes slightly shallower (γ = −1.6 ∼ −1.7) when the loss cone effects are included, and the dimensionless distribution function ḡ(x) drops quickly as the inner boundary is approached.
The left panel of Figure 4 shows examples of star trajectories in the simulation of model M1 when the loss cone is considered.It is evident that in most parts of the cluster, stars evolve along elongated trajectories, allowing them to efficiently reach high eccentricities.This is primarily because the relaxation in angular momentum is usually much faster than the relaxation in energy, except for nearly circular orbits.
The dynamics can be divided into two regions: the full loss cone region with x ≪ x c (or a 2 ≫ a c = M • G/(2x c )), as defined in Equation D23, and the empty loss cone region with x ≫ x c (or a 2 ≪ a c ).For this specific case with M • = 4 × 10 6 M ⊙ , we find that x c ≃ 0.9, corresponding to a c ≃ 0.6r h ≃ 1.7 pc. Figure 4 demonstrates that near the full loss cone region, orbits of stars can still evolve without being immediately disrupted by    tidal forces.They are only consumed if they remain within the loss cone for a sufficient duration to pass through the pericenter.Some of these stars can even move out of the loss cone before being depleted.Stars that remain inside the empty loss cone are highly likely to be disrupted since, within an orbital period, their angular momentum changes are too small to remove them from the loss cone.In the inner part of the cluster, where the relaxation timescale is much shorter, some low-eccentricity stars can sink towards the inner boundary (x > 10 5 , or a 2 ≲ 10 −2 mpc for M • = 4 × 10 6 M ⊙ ) without being disrupted.
In the empty loss cone region (a 2 ≪ 0.6r h for model M1), the distribution of angular momentum follows the expectation of Equation D21.The right panel of Figure 4 shows that for particles with a 2 = 0.01r h or 10 −3 r h , the distribution N (j) obtained from GNC is in good agreement with the expected distribution.For particles with semi-major axes a 2 = r h , since a 2 ∼ 2a c rather than a 2 ≫ a c , the distribution N (j) follows the form of the full loss cone region (N (j) ∝ j) more closely than that of the empty loss cone region.
All of these results, including the relaxation timescale, the energy and density profile evolution, and the E − J evolution of particles, suggest that GNC yields results that are consistent with expectations from the literature, particularly in the case of single-mass particles in galactic nuclei.

Multiple mass components
When there are multiple mass components in the cluster, the dynamics become more complicated due to mass segregation, which drives a flux of heavy components into the central regions while pushing the light ones out.In this study, we use GNC to investigate the steady-state distribution for systems with multiple (up to 12) mass components.
We start with a two-component model, M2, which consists of stars and SBHs with an asymptotic relative number ratio of 10 −3 at the outer boundary (Table 1).The evolution of the density profile for stars and SBHs in M2 from GNC at different times is shown in Figure 5. Comparing these results to those of model M1 in Figure 2, it is apparent that the relaxation of stars is accelerated due to the presence of SBHs.The relaxation timescale for both stars and SBHs, which both have cavities inside ∼ 0.1r h , is approximately 0.15T rlx .These convergence timescales are generally consistent with the regrowth timescale derived by Amaro-Seoane & Preto (2011).
The steady-state results of M2 from GNC are presented in Figure 6.We can observe that, if the loss cone is ignored, the slope index of the stars in M2 from GNC is well consistent with γ = −1.5 across the cluster.When the loss cone effect is included, the density profile of the stars becomes shallower in the inner regions (γ = −1.3∼ −1.4).The SBHs in M2 concentrate near the center, forming a steeper density profile as expected from mass segregation.The density profile of SBHs varies from −2.3 ∼ −1.7, but can be approximately considered as ∼ −2.0 in general, for particles between 10 −3 r h < r < 0.1r h .
In model M2, the differences in the density profiles of both stars and SBHs between GNC and the 1D-FP method are quite small, except in the very inner part of the cluster.When the loss cone is considered, the density profiles of the stars from GNC are usually slightly shallower than those in the 1D-FP method.This discrepancy is partly due to the small anisotropy developed in the cluster.The right panels of Figure 6 show the results of the anisotropy for different components obtained from GNC.The maximum tangential anisotropy for all types of objects is approximately −0.3 (or ≲ −0.6) when the loss cone is not considered (when the loss cone is considered).For SBHs, the anisotropy is very similar to that of stars.
The effects of mass segregation in the case of two component system can be characterized by a ∆ parameter introduced by (Alexander & Hopman 2009) where m H (N H ) is the mass (asymptotic number density) of the heavy component and m L and N L is for the lighter component.We vary the asymptotic number s β at the outer boundary for model M2.The results of the slope index of the two components are shown in Figure 7.The slope index γ varies from −1.75 to −2.25 (−1.5 to −1.6) for the heavy component (light component) when the delta parameter varies from 30 to 10 −3 .Thus, the smaller the ∆ parameter, the stronger the effect of mass segregation.The trend is generally consistent with those of the 1D-FP method calculated in this work, and those from Alexander & Hopman (2009); Preto & Amaro-Seoane (2010).However, we notice that the density slope of the lighter components from GNC is slightly shallower than those of the 1D-FP method.
We further consider models with more mass components.Figure 8 shows the simulation results of model M3, which includes stars, SBHs, and NSs, and model M5, which additionally includes WDs and BDs.We find that in both of these models, the density profiles of the lighter components (e.g., WDs, BDs, NSs) have a slope index of −1.3 ∼ −1.5, which is very similar to those of  1).
the stars.The density profiles of the heavy component, i.e., the SBHs, in these two models are also around −2.0, similar to those in model M2.

Components with a spectrum of masses
We also tested the profiles for systems that consist of components with a broad spectrum of masses in GNC.First, we tested a simple system with a power-law mass function for stars, which was designed for testing pur-poses.In the α-th mass bin m α , the asymptotic number ratio is given by where We use 12 mass bins range from 0.1 − 100M ⊙ , each of which centered at 0.133, 0.237, 0.422, 0.750, 1.33, 2.37,   1.In all panels, the loss cone is considered, symbols are results of GNC, and orange lines are results of 1D-FP method.Note-The fitting results of the slope index as a function of mass from GNC or 1D-FP method to the function γ = −γ0 − ηm = −γ0 − p0/mmaxm.∆ parameter is estimated by Equation 23. , where γ0 and p0 are free parameters and mmax is the maximum mass of the particles.The best-fit values for different models can be found in Table 3.In all panels, the reference lines are as follows: γ = −1.4,which represents the observed density profile of stars in the Galactic center (Genzel et al. 2003;Schödel et al. 2018); γ = −1.5, which represents the expected index of profiles for light components from 1D-FP methods in a two-component model (Bahcall & Wolf 1977); γ = −1.75, which represents the slope index for equal-mass star cusps (Bahcall & Wolf 1976); γ = −2.0,which represents the slope index for stellar-mass black holes in a two-component model (Alexander & Hopman 2009; Amaro-Seoane & Preto 2011).The reference line γ = −2.75represents the slope index expected from the dynamical friction limit (Alexander & Hopman 2009).Note that the error bar for each symbol represents the variations of γ in the range 0.001r h < r < 0.1r h ..22, 7.50, 13.3, 23.7, 42.2, and 75.0 (M ⊙ ).Each mass bin is composed by stars only, for simplicity.We test three models with γ ⋆ = −1.3,−2.3 and −3.3, which is named model M12A, M12B and M12C, respectively.

4
In the case of systems with a broad spectrum of masses, it has already been found that the steadystate number density of particles in each mass bin follows its own power-law profiles n(r) ∝ r γm (Bahcall & Wolf 1977;Keshet et al. 2009;Freitag et al. 2006).1D-FP methods similar to those in Appendix D found that (Bahcall & Wolf 1977;O'Leary et al. 2009) where γ 0 = 3/2, p 0 is a parameter mainly determined by mass segregation and m max is the maximum mass among all bins.The ∆ parameter defined by Equation 19now becomes (Alexander & Hopman 2009) where the components are separated into high and low mass groups if specifying a given boundary of mass.
For all models we adopt 5M ⊙ as the boundary separating the low and high mass groups.The ∆ parameters for different models are shown in Table 3.We fit simultanously Equation 22to the slope index results of GNC, and the best fitting value of γ 0 and 0 parameter for different models are shown in Table 3.
Figure 10 show the γ m as a function of m in different models from GNC.In most models the slope index of light components (≲ 5M ⊙ ) varies very slowly around −1.3 ∼ −1.4.We notice that the slope index of stars is consistent with those of the observered stars at the Galactic center (∼ −1.4,Genzel et al. 2003;Schödel et al. 2018).
For heavy components ( > ∼ 5M ⊙ ) the density profile is steeper for more massive particles.In most models (expect M12C) the slope index of the most massive bin is between −1.75 and −2.0.For weak mass segregation (∆ > 0.5, e.g., in model M12A), p 0 ∼ 0.35, which is slightly larger than ∼ 0.3 expected from the 1D-FP method Bahcall & Wolf (1977).When the mass segregation is strong (∆ ≲ 0.5), p 0 can be 0.5 ∼ 0.9 and apparently increased with the parameter ∆ and is related to the the maximum mass of particles.
Model M12C have the strongest mass segregation effects among all models, where the most heavy com-ponent (75M ⊙ ) has a density slope of −2.2 ∼ −2.3, which is close to the expected limit of dynamical friction (= −11/4) (Alexander & Hopman 2009).
Table 3 and Figure 10 suggest that, even in the case of a broad spectrum of masses, both results from GNC and 1D-FP are in general consistent with each other, although we find that the 1D-FP usually over predict the density of particles in the inner parts of the cluster, compared with those from GNC.
In reality, however, the mass function of particles, including stars and compact objects, is a slightly complex function and does not strictly follow a continuous power-law profile.For a more realistic treatment of the boundary condition, here we use MOBSE (Giacobbo et al. 2018;Giacobbo & Mapelli 2018) 4 to generate stellar populations after a continuous star formation last for 10Gyr.We use Kroupa initial mass function (Kroupa 2001) for stars between 0.01 and 150M ⊙ .The results are shown in Figure 9 and the asymptotic number ratio of each mass bin is listed in Table 2.We test two different models, named M7 and M10, of which adopt metalicity Z = 0.002 and Z = 0.02, respectively.
The slope index as a function of masses can be also found in Figure 10.The ∆ parameter of M7 and M10 is given by 0.16 and 0.57 (see Table 3), respectively, suggesting that the mass segregation in M7 is relatively stronger than those in M10.The power-law index for SBHs is around −2.0 for M7 (at ∼ 10M ⊙ ) and around ∼ −1.8 for M10 (at ∼ 50M ⊙ ).These results suggest that in real cluster system the ∆ parameter is ∼ 0.1−0.6, and that the mass segregation is modest such that the slope index of the most heavy single component is about −1.8 ∼ −2.0.

Consumption rates of particles in the loss cone
Stellar objects will be destroyed if they are wandering too close to the MBH.In our simulations this happens when particles move inside the loss cone, and that the particles must have passed through the pericenter of the orbit.In steady state there is a constant flux of particles that vanish inside the loss cone.In this section we examine the flux of different kinds of particles vanished in the loss cone by both GNC and the 1D-FP method.
Left panel of Figure 11 show the results of the dimensionless flux xI(x) (per log x) for model M1 with M • = 4 × 10 6 M ⊙ .The flux at x ≪ x c are expected to be different with those at x ≫ x c , where x c ≃ 0.9 (≃ 0.56r h ) (Equation D23), of which the former follows Equation D16, and the later Equation D17.We can see that there is a general consistency between the results of GNC and those from the 1D-FP method.The difference are within a factor of ∼ 2, and slightly larger at the inner and outer parts of the cluster.We notice that 1D-FP method assumes that ḡ(x) follows a power-law profile, and that the angular momentum distribution of particles is isotropic, which are both not entirely true in GNC.The discrepancy at inner regions is possibly because 1D-FP method over predict the number of particles at that region (See the left panel in Figure 3).The loss cone at the outer part of the cluster (x ∼ 1) is not either full (q ≫ 1) nor empty (q ≪ 1) (for M • = 4 × 10 6 M ⊙ ), thus the predicted flux rates from 1D-FP method is less accurate.
We also investigate the consistencies of these two methods for MBHs with different masses.For a given type of stellar object, the size of the loss cone, and also the consumption rate in the loss cone vary with the mass of MBH.The transition dimensionless energy x c decrease with the mass of MBHs.As shown in Figure 11, we find general consistency between GNC and the 1D-FP method for MBHs with masses range from 10 4 M ⊙ to 10 8 M ⊙ .We notice that the dimensionless flux of empty loss cone appears in a self-similar way above the transition energy, although the mass of MBH is different.
We also examine the consumption rates for models consist with multiple mass components and different stellar types.Figure 12 show the results of the flux of stars, NSs and SBHs in model M3.We find that 1D-FP method also have general consistencies with GNC for all mass components.The transition energy and the amount of consumption rates of stars are weakly affected by the presence of other components.As the loss cones of NSs and SBHs are much smaller than those of stars, the position of transition is deeper, which is located at x c ∼ 1.7 (a c ∼ 0.3r h ).It seems that the 1D-FP method slightly lower estimates the consumption rates of SBHs at the empty loss cone region by about a factor of < 2 ∼ 3, compared to those from GNC.
The above estimations are for stellar populations bound to the cluster.The total consumption rate can then be obtained by summing both the bound and unbound populations (Equation D25).Note that the consumption rates of bound population will be larger than those of unbound one when x c > ∼ 1.6.For stars, that happens when the mass of MBH is ≤ 10 6 M ⊙ .

DISCUSSION AND CONCLUSIONS
In this work, we have developed a Monte Carlo method (GNC) to study the dynamical evolution of multiple mass components in a star cluster containing a massive black hole at its center.The basic version of the method presented here incorporates two-body relaxation and the effects of the loss cone.For the first time, we calculate the two-body relaxation of multiple mass components based on a two-dimensional (energy and angular momentum) Fokker-Planck (FP) Monte Carlo method, allowing us to obtain steady-state solutions for multiple mass components under given boundary conditions.
We find that GNC yields results that are consistent with those obtained from one-dimensional FP (1D-FP) methods.This consistency holds for various models consisting of equal or multiple mass components and for black holes spanning different mass scales.We also investigate the anisotropy of the cluster and observe the development of tangential anisotropy in the inner regions, while the outer parts can remain isotropic.This anisotropy is the primary cause of the discrepancies between our method and 1D-FP methods, as the latter assume isotropy of cluster.
Nevertheless, it is important to emphasize that GNC is a Monte Carlo method that offers flexibility and advantages over both 1D-FP methods and N-body simulations.Some of these advantages include: The ability to incorporate complex dynamics in a straightforward manner, such as resonant relaxations, stellar collisions, tidal dissipations, gravitational wave dissipations of orbits, complex dynamics of binaries, and more, while simultaneously considering the twodimensional evolution of particles in energy and angular momentum.
Also, the ability to increase the number of particles in the simulation arbitrarily using weighting methods.This is particularly useful for rare objects such as stellar black holes or stellar binary black holes, allowing us to obtain better statistical results for these particles.
The capability to achieve well-converged density profiles and accurate particle evolution down to distance scales of 10 −5 r h from the black hole using particle splitting and weighting methods.This level of accuracy is challenging to achieve with N-body simulations, which typically approach 0.01−0.001rh .Accurate particle evolution in these regions is crucial for studying phenomena such as tidal disruption of stars and extreme-mass-ratio inspirals (EMRIs).
However, it is important to acknowledge some limitations of our current version of the method that need to be addressed in future studies.The most significant limitation is the absence of the stellar object potentials in estimating particle energy and angular momentum.This omission introduces inaccuracies in the density profile and flux into the loss cone from our method, particularly in regions beyond r h .As a result, the results in the outer parts of the cluster are only approximate.Additionally, due to this limitation, we can currently only consider the steady state of the cluster by implementing outer boundary conditions.In the future, when we include the effects of stellar potentials, we will be able to obtain the time-dependent evolution of the cluster and the mass of the black hole.
Nevertheless, our method remains accurate below regions of 0.1r h , where the potential of the black hole dominates over that of the stars.Therefore, we can still apply our method to study various interesting events in galactic nuclei driven by dynamical processes occurring in the inner regions of the cluster.Examples of such events include partial or full tidal disruption of stars of different types and masses, as well as the properties of different types of gravitational wave sources, both first-generation and multiple-generation sources.These applications are particularly relevant for future multimessenger observations in galactic nuclei.4. Thus, it is equivalent to a scheme where every particle that has changed its energy and angular momentum (E and J) values is immediately replaced by a new particle with the same E and J values, ensuring that the distribution of particles on the boundary remains unchanged.However, in the case where a particle successfully moves into the cluster by crossing the boundary (E < E B → E > E B ) at a specific time t = t 0 , to compensate for the void created by the injection of the particle, a new particle with the exact same previous values of E and J should be created at that time.By repeating this process for all particles, there will be a constant flow of particles injected into the cluster from the outer boundary.
On the other hand, if a particle is located inside the boundary E B and attempts to move outside of it (E > E B → E < E B ) in the next step, its simulation should be stopped.In other words, the particle is eliminated at that boundary.The injection of particles from the outer boundary and the elimination of particles that flow out of it will gradually balance each other, ultimately creating a smooth distribution of particles near the outer boundary.
According to the described schemes, the Monte-Carlo simulations of particles in GNC proceed as follows, one-by-one in sequence: If the particle is outside the boundary of the cluster, with E = E 0 < E B and J = J 0 , then it follows this process: 1. Obtain the changes δE and δJ, as well as δt using the functions provided in Appendix B. Update t by setting t → t + δt. 5. Perform the Monte-Carlo steps for particles inside the cluster, as described below.
6.If the particle ends up with E < E B , set t = t 0 , E = E 0 , and J = J 0 , then return to Step 1. Otherwise, proceed to the next step.
7. Create a new particle at time t = t 0 , E = E 0 , and J = J 0 .
8. Terminate the simulation for the current particle.
If a particle is within the cluster, the simulation runs in the following way: 1. Calculate the time step, δE, δJ, and the change in mean anomaly δM using the functions provided in Appendix B. Update t by setting t → t + δt.
2. If J < J lc , where J lc is the loss cone angular momentum, test whether the mean anomaly of the orbit has passed through the pericenter.If M < π and M + δM > π, or if δt > P , the particle is considered destroyed by the loss cone, and the simulation for that particle stops.If not, proceed to the next step.

Save
h , the simulation for the particle stops because it moves outside the outer boundary.5.If E > E max = x max σ 2 h , the simulation for the particle stops because it moves outside the inner boundary.6.If E i /(x 0 σ 2 h ) < 10 i and E f /(x 0 σ 2 h ) > 10 i , create Π − 1 clone particles.7. If E i /(x 0 σ 2 h ) > 10 i and E f /(x 0 σ 2 h ) < 10 i , and if RND(0, Π) > 1, the simulation for the particle stops.Otherwise, proceed to the next step.8.If the simulation time has reached the limit, the simulation stops.Otherwise, return to Step 1.The details of the particle simulation in these two cases are also illustrated in Figure 13.
In our simulation, the dimensionless angular momentum of particles can range from j min to 1.If j moves inside j min (or outside 1), we apply j → 2j min − j (or j → 2 − j), as D EE diverges near j → 0. The value of j min should be much smaller than the minimum size of the loss cone for all types of particles in the simulation.In most cases, it is sufficient to set j min to a value between 10 −4 and 10 −3 .

B. TIME STEPS AND EVOLUTION IN TWO-BODY RELAXATION
The time step δt of simulation for a particle required by the two-body relaxation is given by (Shapiro & Marchant 1978) If the loss cone effect is additionally considered, and the size of the loss cone is J LC , then the time steps need to be adjusted according to the following requirement: If the particle is within the loss cone and has not been previously destroyed by it in a previous step, we additionally require that the next time step cannot be larger than the period of the orbit: ) where P (E) is the orbital period of the particle.
Given the time step, we evolve the E, J and M according to two-body relaxation where y 1 and y 2 are two unit normal random numbers with correlation ρ = Note that in this work, we only consider the effects of two-body relaxation on the orbits to test the basic version of GNC.It is straightforward to add the effects of resonant relaxation, gravitational wave dissipation, and other factors in each time step.These aspects will be discussed in our future works.

C. MONTE-CARLO SIMULATIONS AND STEADY STATE SOLUTIONS
We take the following steps to get the steady-state distributions of particles: 1. Initially, we generate a number of particles of type β with mass m α , where w α is the weighting factor of the particles in mass bin α (See Section 2.3), and s β (m α ) is the number ratio of the particle beyond the outer boundary (See Section 2.2).
The initial energy distribution E of all particles follows a guessed distribution N ) between 0.03 < x < 100x B , where α E,ini ranges from 0 to 0.25.The distribution of particles in the boundary regions (0.03 < x < x B = 0.05) is determined by the value of α E,ini .However, as long as the thickness of the boundary is sufficiently small, the model results are insensitive to the specific value of α E,ini .Changing it to other values usually does not lead to noticeable differences in the model results.The angular momentum j of each particle at the boundary follows the distribution given by G(j) in Equation 11.Initially, a random mean anomaly between 0 and 2π is assigned to the outer orbit M for each particle.
2. We calculate the number distribution in the mass bin α, N α (E, J), and the dimensionless functions g α (x, j) and ḡα (x) according to Equation 9 and Equation 10, respectively.We then calculate the weighting constant w n by normalizing ḡα (x B ) (usually selected from the first mass bin, α = 1) in Equation 12.This normalization is typically necessary before the third or fourth iterations.After that, the weighting w n remains constant as the distributions of particles near x B converge.
3. We calculate the tables of diffusion coefficients shown in Equation 5 and 6 with size of 96 × 96 (or 120 × 120) in the space of log x − log j.
4. Run the Monte-Carlo subroutines for each particle as described in Appendix A. Run the simulation for a duration that is about a few fractions of the two-body relaxation time at r h , typically between 0.005 and 0.01 times T rlx (r h ).
During the simulation, obtain the diffusion coefficients at any given values of x and j through interpolation or finding the nearest value from the table prepared in the previous step.Keep track of particles that exit normally (exit due to reaching the end of the simulation time).To reduce Monte-Carlo errors, refresh the angular momentum of each particle outside the boundary before the next iteration by replacing j with a new value that follows the distribution G(j).Repeat the process from step 2 until the profiles of ḡα (x) in all mass bins have converged.
Similar to Cohn & Kulsrud (1978), to increase the computational efficiency we can generate the following table of an auxiliary function C(s, e) that can be used for all simulations: where κ = 16π 2 G ln Λ.Our computational cost is approximately proportional to the number of particles and slightly increases with the number of mass bins.The most computationally expensive case in our simulation is model M12A, which takes around 30 hours to simulate 1 million particles with 12 mass bins (as discussed in Section 3.1.2)for a duration of approximately one relaxation timescale (T rlx (r h )).This computation is performed on an Intel Xeon CPU with a clock speed of 2.20 GHz, utilizing 24 threads.For the other models in this work, the simulation times typically range from one to several hours.We generally run the simulations for a duration of at least 1.5 times T rlx (r h ), although the profiles usually converge after approximately 0.3-0.4times T rlx (r h ).

D. REVIEW OF THE ONE-DIMENSIONAL FP METHOD IN CALCULATING THE RELAXATION-AND LOSS-CONE EFFECTS
In one-dimensional FP methods, the evolution of particles primarily focuses on the evolution of energy and assumes that the angular momentum distribution of any particle in the α-th mass bin is always isotropic, i.e., f α (E, J) = fα (E) (Bahcall & Wolf 1976).Consequently, the j-averaged diffusion coefficients of the α-th bin reduce to (Bahcall & Wolf 1977): where A = 2 √ 2πG ln ΛE 5/2 = 2 √ 2πG ln Λx 5/2 σ 5 h , Λ = M • /m α .Under the assumption of isotropic angular momentum distribution, we substitute Equation D8 in Equation 4, and then substitute f (E) by the dimensionless distribution ḡ(x), finally the equation reduces to (Bahcall & Wolf 1977)  where τ 0 is about half of the relaxation time at r h .The loss cone effect can be incorporated into Equation D9 by adding an additional term F lc,α (x) that describes the consumption rates of objects in the loss cone.The value of F lc,α (x) depends on a dimensionless parameter q = D JJ,0 P (E)/J 2 lc , where D JJ,0 = D JJ (J → 0) ≃ 2JD J .This parameter distinguishes between the empty (q ≪ 1) and full (q ≫ 1) loss cone regions.The term F lc,α (x) is related to the physical flux F lc,α (E) (per unit time t and unit E) by the following equation: where we have defined I(x) = F lc,α (x)x −5/2 , (D13) which is a dimensionless flux and F 0 is given by According to the physical flux given by Lightman & Shapiro (1977), we have F lc,α (E) = q 4π 2 J 2 lc fα (E) ln J c /J lc , q ≪ 1 = 1.442π 2 J 2 lc fα (E), q ≫ 1. (D15) Substituting the above equation into Equation D12, we obtain that, when q ≪ 1, which is exactly the Equation 6 of Hopman & Alexander (2006).If q ≫ 1, similarly, we have F lc,α (x) = 0.54 π 3/2 x 5/2 ḡα (x) x 5/2 ḡα (x), q ≫ 1 (D17) where Note that when deriving Equation D16 and D17, it is necessary to make the assumption that ḡβ (x) (for all mass bins β) scales with some power of x (Lightman & Shapiro 1977;Bahcall & Wolf 1977), even though the actual form

Figure 1 .
Figure 1.This illustration shows the treatment of boundary and the evolution of multiple types of particles.Particles outside the outer boundary (the left panel) flow into the cluster following the procedure described in Section 2.2 and Appendix A. Inside the outer boundary particles move in the two dimensional E − J space (mid and right panel).Particles pass through x = 10 i , i = 1, • • • , 4 from the left will be split to multiple clone particles, of which have smaller weights to ensure the conservation of particle number (the right panel).The boundary zone ranges from x = 0.03 to xB, where xB = 0.05 and the simulation zone ranges from x = xB to xmax = 10 5 .Symbols filled with light colors show the possible positions where the simulation of a particle stops.

Figure 2 .
Figure 2. The evolution of distributions of equal mass star cluster (m⋆ = 1M⊙, model M1) as a function of time.Loss cone effects are ignored.Left panel: Density profile evolution of a cluster as a function of time in unit of T rlx , where T rlx is the two-body relaxation time at r h .The mass of the MBH is M• = 4 × 10 6 M⊙.Right panel: Same as the left one but for the dimensionless distribution function ḡ(x).

AFigure 3 .
Figure3.Comparison of our method and those from 1D-FP methods for single mass of stars (model M1, See Table1).Top panels and bottom panels are results ignoring and including the loss cone effect, respectively.Panels from the left to the right are: The dimensionless energy distribution function ḡ(x); the number density n(r); The slope index of the number density γ = d ln n(r)/d ln r; and the aniostropy.The colored lines in the first three panels are results from the 1D-FP method.In the top left panel, the black solid line are the results fromBahcall & Wolf (1976), which adopt xmax = 10 4 .The magenta dashed line are the results from our method adopting xmax = 10 4 . .

Figure 5 .AFigure 6 .
Figure 5. Evolution of the distributions of stars (left panel) and SBHs (right panel) in model M2 as a function of time by GNC.Loss cone effects are ignored.Different lines are results at different time of simulation, in unit of T rlx (r h ), where T rlx (r h ) is the two-body relaxation time at r h .The mass of the MBH is M• = 4 × 10 6 M⊙.

Figure 7 .
Figure 7.The slope index γ = d ln n(r)/d ln r of density (averaged between 0.001r h < r < 0.1r h ) as a function of the ∆ parameter in model M2 which consists of 1M⊙ stars and 10M⊙ SBHs.Those colored lines are results from 1D-FP methods.The black (or colored) filled circle/cross mark the results from GNC adopting M• = 10 5 M⊙ (or M• = 4 × 10 6 M⊙).

AFigure 8 .
Figure8.Top panels: Similar to Figure3, but for the model M3, which consists of stars, neutron stars and SBHs; Bottom panels: Similar to the top ones but for the model M5, consisting of five components.For details of model M3 and M5 see Table1.In all panels, the loss cone is considered, symbols are results of GNC, and orange lines are results of 1D-FP method.

Figure 9 .
Figure9.Relative abundance of stellar populations from MOBSE after 10 Gyrs of continuous star formation with Kroupa initial mass function(Kroupa 2001).Left and right panel shows the results given metalicity of z = 0.02 and z = 0.002, respectively.The asymptotic number ratio of different stellar objects at different mass bins are shown in Table2.

Figure 10 .
Figure10.Top panel: The symbols represent the slope index γ = d n(r)/d ln r of the density (averaged between 0.001r h < r < 0.1r h ) as a function of the mass of components obtained from GNC (including the loss cone).The dashed lines with the same color as the symbols represent the results from the 1D-FP method.Bottom panel: Similar to the top panel, but the dashed lines with the same color as the symbols represent the fitting results of γ = −γ0 − p0/(mmaxm), where γ0 and p0 are free parameters and mmax is the maximum mass of the particles.The best-fit values for different models can be found in Table3.In all panels, the reference lines are as follows: γ = −1.4,which represents the observed density profile of stars in the Galactic center(Genzel et al. 2003;Schödel et al. 2018); γ = −1.5, which represents the expected index of profiles for light components from 1D-FP methods in a two-component model(Bahcall & Wolf 1977); γ = −1.75, which represents the slope index for equal-mass star cusps(Bahcall & Wolf 1976); γ = −2.0,which represents the slope index for stellar-mass black holes in a two-component model(Alexander & Hopman 2009; Amaro-Seoane & Preto 2011).The reference line γ = −2.75represents the slope index expected from the dynamical friction limit(Alexander & Hopman 2009).Note that the error bar for each symbol represents the variations of γ in the range 0.001r h < r < 0.1r h .

Figure 11 .Figure 12 .
Figure 11.The dimensionless flux xI(x) (per log x, by Equation D13) of stars to the loss cone for model (M1) which consists of equal mass of stars.The mass of MBH is M• = 4 × 10 6 M⊙ (left panel), 10 4 M⊙ (middle panel) and 10 8 M⊙ (right panel).Symbols are results of GNC.The dotted (or dashed) orange lines show the analytical results from 1D-FP method in the case of full (or empty) loss cone, according to Equation D16 (or Equation D17).The arrow marks the transition position (xc) according to Equation D23.The black stars in the middle and right panel are the results of GNC given M• = 4 × 10 6 M⊙ (the stars in the left panel), which is used for comparison.
supported in part byNational Natural Science Foundation of China under grant No. 12273006, 11603083, U1731104, the Natural Science Foundation of Guangdong Province under grant No. 2021A1515012373.This work was also supported in part by the Key Project of the National Natural Science Foundation of China under grant No. 11733010,12133004.The simulations in this work are performed partly in the TianHe II National Supercomputer Center in Guangzhou.We acknowledge the funds from the "European Union NextGenerationEU/PRTR", Programa de Planes Complementarios I+D+I (ref.AS-FAE/2022/014).APPENDIX A. MONTE-CARLO SIMULATION OF A PARTICLE In Monte-Carlo simulations, particles outside the boundary E < E B (E B = x B σ 2 h ), which do not move into the cluster in the next timestep (E < E B → E < E B ), should be considered as fixed and not move in E − J spaces.This is because the outer boundary, by definition, represents the Dirichlet boundary condition of Equation

2
. If E + δE > E B , proceed to Step 3. Otherwise, save t 0 = t and continue to Step 4. 3.If the simulation time has reached the limit, terminate the simulation.Otherwise, return to Step 1.4.Update E by setting E → E + δE, J → J + δJ, and M → M + δM .

Figure 13 .
Figure 13.Flow charts depicting the of particles inside the boundary (left panel) and those outside the boundary (right panel).Π is the clone factor for particles crossing the dimensionless energy position x/x0 = 10 i , i = 1, • • • , 4, where x = E/σ 2 h , and x0 is an arbitrary energy above which the clone scheme is implimented.

Table 3 .
Fitting of the slope index of density profiles for multiple mass component models