© The Authors 2026

1 Introduction

The long-term evolution of collisional star clusters is driven by two-body relaxation, which operates on timescales that are much longer than the dynamical time. For single-mass systems, this evolution may be approximated by self-similar solutions in which all explicit time dependence is absorbed into a single scale. Such homologous models were first identified by Hénon (1961) and further developed by Hénon (1965) and Lynden-Bell & Eggleton (1980), and they form the basis of the standard theoretical description of gravothermal evolution (see also Binney & Tremaine 2008; Ito 2021).

In systems with a spectrum of stellar masses, relaxation drives the system towards energy equipartition, but full equipartition is generally unattainable (see Spitzer 1969; Inagaki & Saslaw 1985, and subsequent works). These studies have established that sufficiently massive components decouple dynamically from the lighter background and segregate on accelerated time scales - a phenomenon commonly referred to as the Spitzer instability.

Despite this, many numerical studies (including orbit-averaged Fokker-Planck, Monte Carlo, gaseous, and N-body models) have shown that multi-mass star clusters can exhibit approximately self-similar evolution globally, particularly in the post-core-collapse phase (e.g. Cohn 1980; Giersz & Heggie 1994; Takahashi 1995; Giersz & Heggie 1996; Giersz & Spurzem 2000; Pavlik & Subr 2018). They demonstrate long-lived, near-homologous density profiles accompanied by persistent mass segregation. These results suggest that self-similarity is not destroyed by the presence of multiple masses, but rather modified systematically. However, to date, this behaviour has largely been described phenomenologically. In particular, the structural stability of the classical single-scale self-similar solution of Lynden-Bell & Eggleton (1980) under mass-dependent relaxation has not been established, nor its connection to the Spitzer instability and the near-homologous evolution seen in numerical experiments.

In this paper, we analyse the stability of self-similar evolution in collisional star clusters containing multiple mass components. By ‘instability’, we mean the loss of a single, mass-independent similarity scale, rather than a disruption of the overall homologous structure. Therefore, we focus on whether different mass components can consistently share an identical self-similar scaling in a collisional system. Using a gaseous-model formulation, we study the response of a self-similar background to massdependent perturbations driven by energy exchange. We first consider the isotropic case and derive the corresponding instability growth rates and then relax this assumption to examine the role of velocity anisotropy. Finally, we discuss the implications of the analysis for multi-mass self-similarity and its thermodynamic interpretation.

2 Theoretical framework

2.1 Physical assumptions and definitions

We adopted the standard assumptions underlying gaseous models of collisional star clusters:

1. Spherical symmetry;

2. Collisional regime dominated by two-body relaxation, with the characteristic timescale t_rel;

3. Slow evolution, that is, dynamical interactions happen on a much shorter time scale than relaxation and than the global structural evolution of the system (t_dyn ≪ t_rel ≪ t_evol);

4. Locally near-Maxwellian velocity distributions;

5. Single-mass background population;

6. Isotropic velocities (later also relaxed to the anisotropic regime);

7. Continuum limit (many stars per volume element).

These assumptions allow us to work with moment equations of the Boltzmann equation and to define a homologous background solution.

We use ρ(r, t) to denote the mass density, u(r, t) the bulk radial velocity, σ²(r, t) the one-dimensional velocity dispersion, P = ρσ² the pressure, and Φ(r, t) the gravitational potential. The continuity equation then is $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2\rho u\right)=0.$$(1)

The momentum equation follows from the first velocity moment $$\frac{\mathrm{\partial }u}{\mathrm{\partial }u}+u\frac{\mathrm{\partial }u}{\mathrm{\partial }r}=-\frac{1}{\rho }\frac{\mathrm{\partial }P}{\mathrm{\partial }r}-\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }r},$$(2)

with ∂_rΦ = GM(r)/r². The isotropic energy equation (second moment) is $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\frac{3}{2}\rho {\sigma }^2\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left[r^2(\frac{3}{2}\rho {\sigma }^{2}u+F)\right]=-P\frac{\mathrm{\partial }u}{\mathrm{\partial }r},$$(3)

where F = −κ ∂_rσ² is the conductive heat flux (with conductivity set by the local relaxation time as κ ∝ P/t_rel).

2.2 Self-similar (homologous) solution

A solution is self-similar when all explicit time dependence can be absorbed into global scale factors. We can therefore write the background (single-scale) homologous solution as $$\rho (r,t)={\rho }_0(t)f(x),$$(4) $${\sigma }^2(r,t)={\sigma }_0^2(t)s(x),$$(5) $$u(r,t)={\dot{r}}_0(t)x,\text{where}x\equiv r/r_0(t).$$(6)

The assumed linear radial velocity field is the standard consequence of homology and follows from dimensional balance in the momentum equation (see, e.g. Hénon 1961, 1965; Lynden-Bell & Eggleton 1980). Substitution into the continuity equation, in Eq. (1), yields $$\frac{{\dot{\rho }}_0}{{\rho }_0}f-\frac{{\dot{r}}_0}{r_0}x{f}^{\prime }+\frac{{\dot{r}}_0}{r_0}\frac{1}{x^2}(x^{3}f{)}^{\prime }=0$$(7)

(with the dot being ∂_t and the prime ∂_x). The separation of variables (valid for arbitrary x) requires $$\frac{{\dot{\rho }}_0}{{\rho }_0}=-\alpha \frac{{\dot{r}}_0}{r_0},\text{hence}{\rho }_0\propto r_0^{-\alpha }.$$(8)

The exponent α is the usual similarity exponent of the inner halo. Numerical N-body and Fokker-Planck studies give a value of α ≃ 2.2 in core-collapse models (see, e.g. Cohn 1980; Takahashi 1995; Pavlik & Subr 2018). With this condition, the continuity equation is satisfied for arbitrary f (x) and the remaining moment equations reduce to ordinary differential equations for the similarity profiles f (x) and s(x), together with ordinary temporal evolution of the scales r₀(t) and σ²₀(t).

We emphasise two points that are used repeatedly below: (i) self-similarity is a global statement about the separation of time and space variables and not a local approximation; and (ii) the single-scale ansatz fixes a unique relation among the scale functions, which component-dependent relaxation effects cannot change unless those effects are non-negligible and do not share the same scale-free structure.

2.3 Mass perturbation: Formulation and linearisation

We introduce a dilute tracer population of particles with mass m, density ρ_m ≪ ρ (i.e. the equations for ρ_m are linear), isotropic velocity dispersion σ²_m(r, t), and pressure $P_m={\rho }_m{\sigma }_m^2$. The tracer does not perturb the background hydrostatic fields, so the background solution remains essentially self-similar. The tracer obeys the same continuity equation as background matter, advected by the homologous velocity field, expressed as $$\frac{\mathrm{\partial }{\rho }_m}{\mathrm{\partial }t}+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_{m}u\right)=0,$$(9)

while its internal (kinetic) energy evolves through encounters with the background. Using the standard kinetic result for two-body equipartition (Spitzer 1969), we write the local, orbit-averaged energy-exchange $$\frac{\mathrm{d}{\sigma }_m^2}{\mathrm{d}t}=-\frac{1}{t_{\mathrm{e}\mathrm{q}}(m)}({\sigma }_m^2-\mu {\sigma }^2),$$(10)

where d/dt = ∂_t + u∂_r is the convective derivative, μ ≡ 〈m〉/m, with 〈m〉 being the mass of the background stars, and t_eq(m) is the equipartition time for the tracer - to leading order t_eq(m) ~ μ t_rel. The tracer energy equation (second moment) is then expressed as $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\frac{3}{2}{\rho }_m{\sigma }_m^2\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left[r^2(\frac{3}{2}{\rho }_m{\sigma }_m^{2}u+F_m)\right]=-P_m\frac{\mathrm{\partial }u}{\mathrm{\partial }r}+Q,$$(11)

where $Q=(3/2){\rho }_m{\mathrm{\partial }}_t{\sigma }_m^2$ is the volumetric energy-exchange to the background and Fm is the (generally small) conductive flux of the tracer. We linearise about the local equipartition, $${\sigma }_m^2=\mu {\sigma }^2+\delta {\sigma }_m^2,$$(12)

where δ denotes small perturbations. After neglecting the small tracer flux and using the background isotropic energy equation, all zeroth-order terms cancel out identically. Consequently, Eq. (11) reduces to $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m\delta {\sigma }_m^2\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m\delta {\sigma }_m^{2}u\right)=-\frac{{\rho }_m}{t_{\mathrm{e}\mathrm{q}}}\delta {\sigma }_m^2,$$(13)

where all terms are linear in the perturbation (see Appendix B for the derivation).

For the pressure perturbation, we write to linear order $$\delta P_m=(\mu {\sigma }^2+\delta {\sigma }_m^2)\delta {\rho }_m+{\rho }_m\delta {\sigma }_m^2\approx \mu {\sigma }^2\delta {\rho }_m+{\rho }_m\delta {\sigma }_m^2.$$(14)

Under the quasi-hydrostatic assumption - i.e. t_dyn ≪ t_self, where t_dyn is the local dynamical (e.g. crossing) time and t_self is the timescale over which the background self-similar solution evolves (e.g. core collapse or cluster expansion) - pressure readjusts rapidly. To leading order in the perturbation, the local pressure perturbations vanish (δP_m ≈ 0), with the readjustment controlled by the equipartition timescale, t_eq¹. Consequently, $$\frac{\delta {\rho }_m}{{\rho }_m}=-\frac{1}{\mu }\frac{\delta {\sigma }_m^2}{{\sigma }^2}$$(15)

and we obtain the linearised evolution for the density perturbation, $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\delta {\rho }_m\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2\delta {\rho }_{m}u\right)=\frac{\delta {\rho }_m}{t_{\mathrm{e}\mathrm{q}}}.$$(16)

This compact equation contains the two competing processes: advection by the similarity flow on the left-hand side and the source coming from the local equipartition on the right-hand side. Equation (16) is the basis for the instability analysis below.

2.4 Instability growth rate

For the tracer density, we define $${\rho }_m(r,t)={\rho }_0(t)g(x,t)$$(17)

where the shape function g may, in principle, depend on time. We also introduce the dimensionless (logarithmic) similarity time variable $$\tau \equiv \ln r_0(t)(+\text{const.}),$$(18)

with the following identity: $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}=\frac{\mathrm{\partial }\tau }{\mathrm{\partial }t}\frac{\mathrm{\partial }}{\mathrm{\partial }\tau }=\frac{{\dot{r}}_0}{r_0}\frac{\mathrm{\partial }}{\mathrm{\partial }\tau }=\frac{1}{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}\frac{\mathrm{\partial }}{\mathrm{\partial }\tau }.$$(19)

Expressing Eq. (16) in (x,τ) coordinates and using the leadingorder advection result from the continuity equation, in Eq. (9), yields the relation $$\frac{\mathrm{\partial }g}{\mathrm{\partial }\tau }+(3-\alpha )g=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}}g,$$(20)

with α from Eq. (8). We again note that the source term on the right-hand side arises directly from energy exchange leading to equipartition. For fixed x, Eq. (20) reduces to an ordinary differential equation in τ, with the solution $$g(\tau )=g_0(x)e^{\lambda (m)\tau },\text{where}\lambda (m)\equiv \frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}(m)}-(3-\alpha ).$$(21)

The quantity λ(m) is, therefore, the linear growth (or decay) rate of the tracer in similarity time. Converting back to physical variables, exponential growth in τ corresponds to a power-law separation in r₀ and a faster-than-homologous concentration of the mass component.

The instability condition is simply λ> 0, that is, $$\frac{m}{⟨m⟩}>(3-\alpha )\frac{t_{\mathrm{r}\mathrm{e}\mathrm{l}}}{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}.$$(22)

Physically, this compares the equipartition (or mass segregation) time of a given component with the homologous rescaling time. When two-body encounters relax the mass component faster than the similarity flow can rescale it, the mass component departs exponentially from the single-scale self-similarity². However, this does not mean that the entire cluster loses its scale-free evolution, nor that the density profiles cannot appear self-similar. What it says is that there is no single exponent λ(m) in the solution of e^λ(m)τ that could make $\lambda (m)=0,\mathrm{\forall }m$, and consequently, different mass components cannot share the same similarity scaling.

We note that the growth rate λ(m) depends only on the ratio t_self/t_eq(m) and on the structural exponent α. It does not depend sensitively on the detailed form of the conductive flux. Thus, any conductivity of the form κ ∝ P/t_rel would lead to the same scaling of the growth rate, so the breaking of singlescale self-similarity is robust with respect to the heat-transport model. More to the point, although we derived this result within the gaseous model, the present instability and its growth rates also arise in orbit-averaged Fokker-Planck treatments through the same energy-exchange terms (see, e.g. Spitzer 1987; Inagaki & Wiyanto 1984). Therefore, the result reflects the underlying collisional kinetics rather than the specific fluid approximation.

2.5 Velocity anisotropy

Velocity anisotropy is introduced using the Osipkov-Merritt model (Osipkov 1979; Merritt 1985), in which the anisotropy profile is controlled by a single radius r_a and given by $$\beta (r)\equiv 1-\frac{{{\sigma }_{\mathrm{t}\mathrm{a}\mathrm{n}}}^2(r)}{{{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}}^2(r)}=\frac{r^2}{r^2+{r_{\mathrm{a}}}^2},$$(23)

where σ²_rad denotes the radial velocity dispersion and σ²_tan is the tangential velocity dispersion. The value β ≈ 0 corresponds to an isotropic velocity distribution, while β → 1 is fully radial, and β < 0 is a tangentially biased velocity distribution.

For the tracer population, we define the density ρ_m(r, t), the mean radial velocity u(r, t) (as the background), and velocity dispersions ${{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2$ and ${{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2$, together with the tracer anisotropy parameter β_m, defined analogously to the expression in Eq. (23). As before, the tracer is dynamically passive and does not modify the background potential or flow.

The anisotropic momentum equation for the background follows directly from the Jeans equation (e.g. Binney & Tremaine 2008): $$\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }u}{\mathrm{\partial }r}=-\frac{1}{\rho }\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\rho {{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}}^2\right)-\frac{2\beta {{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}}^2}{r}-\frac{GM(r)}{r^2}.$$(24)

Anisotropy affects only the pressure terms; therefore, the tracer continuity equation is the same as in the isotropic case in Eq. (9).

In contrast, the tracer energy equation must be split into radial and tangential components. The corresponding second-moment equations are $$\begin{array}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m{{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m{{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^{2}u\right)+2\frac{{\rho }_m}{r}({{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2-{{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2)\\ =-\frac{{\rho }_m}{t_{\mathrm{e}\mathrm{q}}[,rad]}({{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2-\mu {{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}}^2)\end{array}$$(25)

and $$\begin{array}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m{{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m{{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^{2}u\right)-\frac{{\rho }_m}{r}({{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2-{{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2)\\ =-\frac{{\rho }_m}{t_{\mathrm{e}\mathrm{q}}[,tan]}({{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2-\mu {{\sigma }_{\mathrm{t}\mathrm{a}\mathrm{n}}}^2).\end{array}$$(26)

Here the additional terms proportional to $({{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2-{{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2)/r$ arise from the divergence of the anisotropic pressure tensor and represent a centrifugal coupling between the radial and tangential motions.

As in the isotropic case, we linearise about local equipartition by writing $${{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2=\mu {{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}}^2+\delta {{\sigma }_{m,\mathrm{r}\mathrm{a}\mathrm{d}}}^2,{{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2=\mu {{\sigma }_{\mathrm{t}\mathrm{a}\mathrm{n}}}^2+\delta {{\sigma }_{m,\mathrm{t}\mathrm{a}\mathrm{n}}}^2.$$(27)

The tracer density is similarly decomposed into a homologously advected background and a small perturbation. Substitution into Eqs. (25) and (26) shows that, as in the isotropic case, the leading-order advection terms cancel identically when the background is self-similar. The remaining terms govern the evolution of the perturbations.

Introducing similarity variables (x, τ) and decomposing the density perturbation into amplitudes associated with the radial and tangential degrees of freedom, we define the dimensionless functions g_rad(x, τ) and g_tan(x, τ) through $$\delta {\rho }_m\propto {\rho }_0(t)g_{\mathrm{r}\mathrm{a}\mathrm{d}}(x,\tau )\text{and}\delta {\rho }_m\propto {\rho }_0(t)g_{\mathrm{t}\mathrm{a}\mathrm{n}}(x,\tau )$$

when projected onto the radial and tangential moment equations, respectively. We then obtain the coupled system of first-order differential equations: $$\begin{array}{r}\frac{\mathrm{\partial }{g}_{\mathrm{r}\mathrm{a}\mathrm{d}}}{\mathrm{\partial }\tau }+(3-\alpha )g_{\mathrm{r}\mathrm{a}\mathrm{d}}+2\beta \frac{r_0}{r}(g_{\mathrm{r}\mathrm{a}\mathrm{d}}-g_{\mathrm{t}\mathrm{a}\mathrm{n}})=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,rad]}{g}_{\mathrm{r}\mathrm{a}\mathrm{d}},\end{array}$$(28) $$\begin{array}{r}\frac{\mathrm{\partial }{g}_{\mathrm{t}\mathrm{a}\mathrm{n}}}{\mathrm{\partial }\tau }+(3-\alpha )g_{\mathrm{t}\mathrm{a}\mathrm{n}}-\beta \frac{r_0}{r}(g_{\mathrm{r}\mathrm{a}\mathrm{d}}-g_{\mathrm{t}\mathrm{a}\mathrm{n}})=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,tan]}{g}_{\mathrm{t}\mathrm{a}\mathrm{n}}.\end{array}$$(29)

Anisotropy, therefore, affects the dynamics solely by introducing a linear coupling between radial and tangential perturbations through the anisotropic pressure terms.

The self-similar solution implies a flat central core (r ≪ r_a), in which the Osipkov-Merritt anisotropy satisfies β ≪ 1. In this limit, the coupling terms in Eqs. (28) and (29) are negligible and the equations decouple $$\frac{\mathrm{\partial }{g}_{\mathrm{r}\mathrm{a}\mathrm{d}}}{\mathrm{\partial }\tau }+(3-\alpha )g_{\mathrm{r}\mathrm{a}\mathrm{d}}=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,rad]}{g}_{\mathrm{r}\mathrm{a}\mathrm{d}},\frac{\mathrm{\partial }{g}_{\mathrm{t}\mathrm{a}\mathrm{n}}}{\mathrm{\partial }\tau }+(3-\alpha )g_{\mathrm{t}\mathrm{a}\mathrm{n}}=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,tan]}{g}_{\mathrm{t}\mathrm{a}\mathrm{n}}.$$(30)

The solutions are $$g_{\mathrm{r}\mathrm{a}\mathrm{d}}(\tau )\propto e^{{\lambda }_{\mathrm{r}\mathrm{a}\mathrm{d}}\tau },g_{\mathrm{t}\mathrm{a}\mathrm{n}}(\tau )\propto e^{{\lambda }_{\mathrm{t}\mathrm{a}\mathrm{n}}\tau },$$

with the growth rates $${\lambda }_{\mathrm{r}\mathrm{a}\mathrm{d}}=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,rad]}-(3-\alpha ),{\lambda }_{\mathrm{t}\mathrm{a}\mathrm{n}}=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,tan]}-(3-\alpha ).$$(31)

The form is identical to the isotropic case, but the radial and tangential modes are no longer degenerate when t_eq,rad ≠ t_eq,tan.

Away from the centre (β ≠ 0), the centrifugal terms modify the local growth rates. In the regime relevant for mass segregation, we find to leading order $${\lambda }_{\mathrm{r}\mathrm{a}\mathrm{d}}=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,rad]}-(3-\alpha )-2\beta ,$$(32) $${\lambda }_{\mathrm{t}\mathrm{a}\mathrm{n}}=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,tan]}-(3-\alpha )-\beta ,$$(33)

where λ_rad > 0 or λ_tan > 0 lead to radial or tangential instability, respectively. The full coupled eigenvalue problem, including higher-order terms, is solved in Appendix C. The different coefficients in front of β in Eqs. (32) and (33) show that anisotropy both shifts the overall instability strength and lifts the degeneracy between predominantly radial and predominantly tangential modes. This behaviour is in agreement with physical intuition and numerical results.

Specifically, for radially anisotropic systems (β > 0), both eigenvalues are reduced relative to the isotropic case, with the predominantly radial mode receiving the larger correction. Physically, this reflects the enhanced radial kinetic support and the associated radial-orbit heating seen in N-body simulations (e.g. Pavlik & Vesperini 2021, 2022a,b; Aros & Vesperini 2023; Pavlik et al. 2024), which delays the core contraction and slows central relaxation.

For tangentially biased systems (β < 0), both growth rates increase relative to the isotropic case, corresponding to reduced radial heating and more rapid central contraction. This behaviour is consistent with the numerical results of Pavlik et al. (2024), where tangentially anisotropic models exhibit shorter corecollapse times and more rapid inner mass segregation than isotropic and radially anisotropic models.

Moreover, since t_eq,rad/tan ~ μ t_rel, massive tracers violate single-scale self-similarity first, typically in the core. Velocity anisotropy modifies the rate at which this departure occurs, but the ordering of instability among different mass components remains governed primarily by the equipartition timescale, t_eq(m), rather than by anisotropy itself.

2.6 Multi-mass self-similarity

We first recall the isotropic, single-mass self-similar solution from Eq. (4) - i.e. ρ(r, t) = ρ₀(t) f (x), with x ≡ r/r₀(t) - in which all explicit time dependence is absorbed into the single scale r₀(t). We now generalise this framework to a system composed of discrete mass components m_i, with the densities ρ_i(r, t). The total density and gravitational potential satisfy $$\rho (r,t)=\sum _i{\rho }_i(r,t),{\mathrm{\nabla }}^2\mathrm{\Phi }=4\pi G\sum _i{\rho }_i.$$(34)

Each component obeys its own continuity and energy equations, while all components share the same potential.

Suppose all components share a common similarity scale, that is, $${\rho }_i(r,t)={\rho }_{0,i}(t)f_i(x),x=r/r_0(t).$$(35)

Using the linearised tracer evolution in Eq. (20), we find $$\frac{\mathrm{d}}{\mathrm{d}\tau }\left(\frac{{\rho }_i}{\rho }\right)=[\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}(m_i)}-(3-\alpha )]\left(\frac{{\rho }_i}{\rho }\right),\text{where\hspace{0.17em}}\tau =\ln r_0(t).$$(36)

Because the equipartition time t_eq(m_i) depends on mass, the right-hand side is different for different components. Thus, no single choice of r₀(t) can make the evolution of all mass components stationary. Single-scale self-similarity is thus structurally unstable in a multi-mass system, which contradicts Eq. (35).

We must, therefore, relax the assumption of a common scale and allow each mass component to evolve self-similarly with its own characteristic radius, that is, $${\rho }_i(r,t)={\rho }_{0,i}(t)f_i(x_i),x_i\equiv r/r_{0,i}(t),$$(37)

as well as its own power-law normalisation, and dimensionless time (i.e. also similarity time t_self,i), $${\rho }_{0,i}\propto r_{0,i}^{-{\alpha }_i},{\tau }_i\equiv \ln r_{0,i}(t).$$(38)

In general, α_i ≠ α_j for two different masses m_i and m_j.

To remain consistent with a slowly evolving gravitational potential, expressed via Eq. (34), we require that $\sum _i{\rho }_{0,i}{f}_i(x_i)$ remains approximately scale-free over any given radial range. This is possible only if different mass components dominate different regions, which implies a radial stratification by mass. Since heavier components have shorter equipartition times, their instability growth rates are larger, and their scale radii shrink more rapidly. Consequently, we have $$r_{0,i}(t)<r_{0,j}(t)\text{for}{m}_i>m_j.$$(39)

Hence, the system evolves towards a multi-scale, mass-segregated configuration rather than a single-scale similarity solution.

This structural behaviour is seen in numerical simulations (e.g. Inagaki & Lynden-Bell 1983; Giersz & Heggie 1996), which show long-lived, near-homologous post-collapse states with persistent mass segregation. In these models, the global density profile evolves approximately self-similarly, while different mass components remain radially stratified. The numerical results therefore support the interpretation that multi-mass clusters evolve toward a configuration that is effectively a superposition of component-wise self-similar distributions, rather than a strictly single-scale homologous solution.

3 Conclusions

We analysed the stability of self-similar evolution in collisional multi-mass star clusters. Using a gaseous-model formulation, we show that mass-dependent relaxation breaks single-scale self-similarity and drives the system towards a multi-scale configuration in which different mass components may evolve self-similarly on distinct characteristic scales. We further show that velocity anisotropy modifies this behaviour by lifting the degeneracy between radial and tangential modes, while the overall structure of the evolution remains unchanged. The predicted ordering of growth rates and characteristic similarity scales is consistent with trends reported in orbit-averaged Fokker-Planck, Monte Carlo, and N-body studies of multi-mass star clusters. Specifically, this helps explain the results seen in numerical simulations by Pavlik et al. (2024).

When applied to a system containing multiple stellar mass components, our analysis shows that mass-dependent relaxation renders the classical single-scale self-similar solution structurally unstable. The resulting instability naturally leads to mass segregation and to a configuration in which different mass components evolve with distinct characteristic similarity scales. This does not imply a breakdown of global homology; rather, the cluster can remain approximately self-similar in its overall structure while individual mass components follow their own scale-dependent evolution. Such multi-scale, near-homologous, mass-segregated states are seen in numerical models (e.g. Giersz & Heggie 1996). Taken together, these results provide a framework that links the classical self-similar theory to more realistic star clusters.

Acknowledgements

I am grateful to Steve Shore for many stimulating conversations that inspired this work, to Douglas Heggie and Luca Ciotti for helpful insights into the multi-mass description of self-similar star clusters, and to Enrico Vesperini for valuable guidance on the kinematics of star clusters. I have received funding from the European Union’s Horizon Europe and the Central Bohemian Region under the Marie Sklodowska-Curie Actions - COFUND, Grant agreement ID:101081195 (“MERIT”). I also acknowledge support from the project RVO:67985815 at the Czech Academy of Sciences. I thank the referee for a careful and constructive report that helped improve the clarity and presentation of this paper.

References

Appendix A Pressure perturbations and thermalisation

This appendix provides a brief clarification of the assumptions underlying the velocity dispersion (or pressure) perturbation δσ²_m used in Section 2.

The analysis in this paper is based on a moment description of a collisional stellar system, in which the velocity distribution is assumed to be locally relaxed and well approximated by a Maxwellian distribution. Departures from equilibrium are treated at the level of low-order moments. In particular, the perturbation variable $\delta {\sigma }_m^2$ represents a small deviation of the second velocity moment from local equipartition, and may be interpreted as a pressure perturbation of the tracer population.

This approach is related in spirit to the Chapman-Enskog expansion in kinetic theory (see Chapman 1916; Enskog 1917; Chapman & Cowling 1970), where the distribution function is written as $$\mathcal{F}={\mathcal{F}}^{(0)}+ϵ{\mathcal{F}}^{(1)}+\dots ,$$

where F⁽⁰⁾ is a local Maxwellian and the correction F⁽¹⁾ is introduced to account for non-equilibrium fluxes. Although the correction F⁽¹⁾ is often described as a perturbation of the distribution function itself, it only serves to generate corresponding corrections to macroscopic quantities such as the stress tensor and the heat flux, and has no independent dynamical significance beyond its contribution to the moments.

In the present work, we formulate the problem directly at the level of moments. Writing $${\sigma }_m^2=\mu {\sigma }^2+\delta {\sigma }_m^2,$$

as in Eq. (12), does not imply a non-thermal or non-Maxwellian state. Instead, it expresses that the massive component carries slightly more or slightly less kinetic energy than implied by local equipartition at fixed density. Such pressure perturbations are physically realizable and are continuously generated and damped by two-body relaxation. Their decay is governed by the equipartition time t_eq, as expressed by the local energy-exchange term $$\frac{\mathrm{d}{\sigma }_m^2}{\mathrm{d}t}\propto -\frac{1}{t_{\mathrm{e}\mathrm{q}}}({\sigma }_m^2-\mu {\sigma }^2),$$

given in Eq. (10).

Pressure perturbations evolve on the same timescale as mass segregation and structural evolution. By contrast, genuinely non-thermal perturbations of the distribution function (such as beams, phase-space substructure, or higher-order velocity moments) are erased on much shorter timescales and do not persist long enough to influence secular evolution of the star cluster. For this reason, pressure perturbations constitute the relevant slow degrees of freedom in a collisional, self-similar system.

The instability analysed in this paper is, therefore, not a kinetic instability of the distribution function, but a structural instability of a thermalised background when subject to massdependent relaxation. Its origin lies in the competition between the rescaling imposed by the self-similar flow and the local tendency toward equipartition. When the latter acts more rapidly, small pressure perturbations grow in similarity time, leading to a separation of characteristic scales among different mass components.

The same considerations apply in the presence of velocity anisotropy, which is introduced at the level of second moments and does not invalidate the assumption of local thermalisation. It modifies the coupling between radial and tangential components of the pressure tensor and leads to distinct growth rates for radial and tangential perturbations, but does not introduce oscillatory or collisionless behaviour. In the physically relevant regime, the eigenvalues governing the evolution remain real, which reflects the strongly dissipative, entropy-producing nature of collisional stellar dynamics.

Appendix B Tracer energy equation

To show the derivation of the tracer energy equation in terms of the velocity dispersion perturbation, we start from Eq. (11), that is, $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\frac{3}{2}{\rho }_m{\sigma }_m^2\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left[r^2(\frac{3}{2}{\rho }_m{\sigma }_m^{2}u+F_m)\right]=-P_m\frac{\mathrm{\partial }u}{\mathrm{\partial }r}+Q.$$

Assuming F_m → 0 and taking $$Q=-\frac{3}{2}\frac{{\rho }_m}{t_{\mathrm{e}\mathrm{q}}}\delta {\sigma }_m^2,$$(B.1)

which comes from the linearisation about equipartition in Eq. (12), then Eq. (11) becomes $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m{\sigma }_m^2\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m{\sigma }_m^{2}u\right)=-\frac{2}{3}{P}_m\frac{\mathrm{\partial }u}{\mathrm{\partial }r}-\frac{{\rho }_m}{t_{\mathrm{e}\mathrm{q}}}\delta {\sigma }_m^2.$$

Using the linearisation $${\sigma }_m^2=\mu {\sigma }^2+\delta {\sigma }_m^2,$$

we decompose all terms into zeroth-order (background) and first-order (perturbation) contributions. The time-derivative term becomes $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m{\sigma }_m^2\right)=\mu ({\sigma }^2\frac{\mathrm{\partial }{\rho }_m}{\mathrm{\partial }t}+{\rho }_m\frac{\mathrm{\partial }{\sigma }^2}{\mathrm{\partial }t})+\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m\delta {\sigma }_m^2\right),$$

the advective flux term is $$\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m{\sigma }_m^{2}u\right)=\mu \frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m{\sigma }^{2}u\right)+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m\delta {\sigma }_m^{2}u\right),$$

and the pressure term on the right-hand side is $$-\frac{2}{3}{P}_m\frac{\mathrm{\partial }u}{\mathrm{\partial }r}=-\frac{2}{3}{\rho }_m\mu {\sigma }^2\frac{\mathrm{\partial }u}{\mathrm{\partial }r}-\frac{2}{3}{\rho }_m\delta {\sigma }_m^2\frac{\mathrm{\partial }u}{\mathrm{\partial }r}.$$

Collecting all zeroth-order terms (proportional to μσ²), we obtain $$\mu [{\sigma }^2\frac{\mathrm{\partial }{\rho }_m}{\mathrm{\partial }t}+{\rho }_m\frac{\mathrm{\partial }{\sigma }^2}{\mathrm{\partial }t}+\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\rho }_m{\sigma }^{2}u\right)+\frac{2}{3}{\rho }_m{\sigma }^2\frac{\mathrm{\partial }u}{\mathrm{\partial }r}]=0,$$

which vanishes thanks to the background isotropic energy equation (3) applied to the passive tracer density ρ_m. Therefore, only the first-order terms proportional to δσ²_m survive, yielding Eq. (13).

Appendix C Anisotropy eigenvalue problem

We analyse the coupled evolution of radial and tangential density perturbations in the presence of velocity anisotropy Starting from Eqs. (28) and (29), we write $$\begin{array}{r}\frac{\mathrm{\partial }{g}_{\mathrm{r}\mathrm{a}\mathrm{d}}}{\mathrm{\partial }\tau }+(3-\alpha )g_{\mathrm{r}\mathrm{a}\mathrm{d}}+2\beta \frac{r_0}{r}(g_{\mathrm{r}\mathrm{a}\mathrm{d}}-g_{\mathrm{t}\mathrm{a}\mathrm{n}})=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,rad]}{g}_{\mathrm{r}\mathrm{a}\mathrm{d}},\\ \frac{\mathrm{\partial }{g}_{\mathrm{t}\mathrm{a}\mathrm{n}}}{\mathrm{\partial }\tau }+(3-\alpha )g_{\mathrm{t}\mathrm{a}\mathrm{n}}-\beta \frac{r_0}{r}(g_{\mathrm{r}\mathrm{a}\mathrm{d}}-g_{\mathrm{t}\mathrm{a}\mathrm{n}})=\frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,tan]}{g}_{\mathrm{t}\mathrm{a}\mathrm{n}}.\end{array}$$

For an Osipkov-Merritt profile with $$\beta (r)=1-\frac{{{\sigma }_{\mathrm{t}\mathrm{a}\mathrm{n}}}^2(r)}{{{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}}^2(r)}=\frac{r^2}{r^2+{r_{\mathrm{a}}}^2},$$

see also Eq. (23), the velocity-dispersion difference scales as $({{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}}^2-{{\sigma }_{\mathrm{t}\mathrm{a}\mathrm{n}}}^2)\sim \beta {\sigma }^2$. The anisotropy coupling term, therefore, scales as β r₀/r. In the central region (r ≪ r_a), where β(r) ∝ r², the product β r₀/r vanishes and no divergence arises. In the transition region (r ~ r_a), both β and r₀/r are order-unity quantities varying slowly with radius. Retaining only this leading-order scaling in β, and deferring the precise normalisation by orderunity coefficient for the next section, the coupled system can be written in matrix form as $$\frac{\mathrm{\partial }}{\mathrm{\partial }\tau }\left(\begin{array}{c}{g}_{\mathrm{r}\mathrm{a}\mathrm{d}}\\ g_{\mathrm{t}\mathrm{a}\mathrm{n}}\end{array}\right)=\left(\begin{array}{cc}{A}_{\mathrm{r}\mathrm{a}\mathrm{d}}-2\beta & 2\beta \\ \beta & A_{\mathrm{t}\mathrm{a}\mathrm{n}}-\beta \end{array}\right)\left(\begin{array}{c}{g}_{\mathrm{r}\mathrm{a}\mathrm{d}}\\ g_{\mathrm{t}\mathrm{a}\mathrm{n}}\end{array}\right)$$(C.1)

where, for compactness, we substituted $$A_{\mathrm{r}\mathrm{a}\mathrm{d}}\equiv \frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,rad]}-(3-\alpha ),A_{\mathrm{t}\mathrm{a}\mathrm{n}}\equiv \frac{t_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}{t_{\mathrm{e}\mathrm{q}}[,tan]}-(3-\alpha ).$$(C.2)

C.1 Solution and limiting cases

The eigenvalues λ satisfy $$\begin{array}{r}0=det\left(\begin{array}{cc}{A}_{\mathrm{r}\mathrm{a}\mathrm{d}}-2\beta -\lambda & 2\beta \\ \beta & A_{\mathrm{t}\mathrm{a}\mathrm{n}}-\beta -\lambda \end{array}\right)\\ =(A_{\mathrm{r}\mathrm{a}\mathrm{d}}-2\beta -\lambda )(A_{\mathrm{t}\mathrm{a}\mathrm{n}}-\beta -\lambda )-2{\beta }^2,\end{array}$$(C.3)

which yields $${\lambda }_{\pm }=\frac{1}{2}[(A_{\mathrm{r}\mathrm{a}\mathrm{d}}+A_{\mathrm{t}\mathrm{a}\mathrm{n}}-3\beta )\pm \sqrt{(A_{\mathrm{r}\mathrm{a}\mathrm{d}}-A_{\mathrm{t}\mathrm{a}\mathrm{n}}-\beta {)}^2+8{\beta }^2}].$$(C.4)

In the isotropic limit β → 0, the two modes decouple and $${\lambda }_{\pm }\to A_{\mathrm{r}\mathrm{a}\mathrm{d}},A_{\mathrm{t}\mathrm{a}\mathrm{n}},$$(C.5)

where we are recovering the isotropic self-similar result. We note that for A_rad ≡ A_tan, this corresponds to the classical Spitzer instability (Spitzer 1969).

In the regime relevant for mass segregation, A_rad, A_tan ≫ β, the off-diagonal terms in Eq. (C.1) are smaller than the diagonal growth terms. The eigenvectors align asymptotically with the radial and tangential directions. A first-order expansion then gives $${\lambda }_{\mathrm{r}\mathrm{a}\mathrm{d}}\simeq A_{\mathrm{r}\mathrm{a}\mathrm{d}}-2\beta ,{\lambda }_{\mathrm{t}\mathrm{a}\mathrm{n}}\simeq A_{\mathrm{t}\mathrm{a}\mathrm{n}}-\beta ,$$(C.6)

which are the expressions quoted in the main text in Eq. (33).

C.2 Dependence on the anisotropy coupling

The gaseous model constrains only the leading-order scaling of the anisotropy coupling and does not fix its precise numerical normalisation. In particular, retaining only the leading-order dependence on β in Eq. (C.1) leaves an undetermined order-unity factor associated with the local value of r₀/r in the anisotropy-dominated region. We therefore introduce an arbitrary orderunity normalisation coefficient C to assess the sensitivity of the instability to the strength of the anisotropic coupling. Equation (C.1) then generalises to $$\frac{\mathrm{\partial }}{\mathrm{\partial }\tau }\left(\begin{array}{c}{g}_{\mathrm{r}\mathrm{a}\mathrm{d}}\\ g_{\mathrm{t}\mathrm{a}\mathrm{n}}\end{array}\right)=\left(\begin{array}{cc}{A}_{\mathrm{r}\mathrm{a}\mathrm{d}}-2C\beta & 2C\beta \\ C\beta & A_{\mathrm{t}\mathrm{a}\mathrm{n}}-C\beta \end{array}\right)\left(\begin{array}{c}{g}_{\mathrm{r}\mathrm{a}\mathrm{d}}\\ g_{\mathrm{t}\mathrm{a}\mathrm{n}}\end{array}\right).$$(C.7)

In the isotropic limit (β → 0), the eigenvalues again reduce to λ_± = A_rad, A_tan, independently of C. For β ≠ 0, the eigenvalue separation is $${\lambda }_{+}-{\lambda }_{-}=\sqrt{(A_{\mathrm{r}\mathrm{a}\mathrm{d}}-A_{\mathrm{t}\mathrm{a}\mathrm{n}}-C\beta {)}^2+8C^2{\beta }^2}.$$(C.8)

Thus, the ordering of the eigenvalues cannot be altered by any finite order-unity rescaling of the anisotropy coupling. We also note that although tangential anisotropy (β < 0) modifies the growth rates, the two modes remain non-degenerate, and their ordering cannot be reversed solely by changing the sign of β.

In the mass-segregation regime A_rad, A_tan ≫ Cβ, we find $${\lambda }_{\mathrm{r}\mathrm{a}\mathrm{d}}=A_{\mathrm{r}\mathrm{a}\mathrm{d}}-2C\beta +\mathcal{O}({\beta }^2),{\lambda }_{\mathrm{t}\mathrm{a}\mathrm{n}}=A_{\mathrm{t}\mathrm{a}\mathrm{n}}-C\beta +\mathcal{O}({\beta }^2).$$(C.9)

The coefficient C affects only the magnitude of the anisotropic correction, not its sign or the instability criterion.

We therefore reach the following overall conclusion: the anisotropic instability, its growth-rate ordering, and its isotropic limit are robust to order-unity uncertainties in the normalisation of the anisotropy coupling.