Spherical Stellar Systems Research Articles

Exploring the dynamics of collisionless spherical stellar systems using the matrix method: Insights from the dilation mode

Context. Analytical solutions to the perturbed equations that govern self-gravitating collisionless stellar systems are crucial for both code testing and theoretical insights. For spheres, a solution has been known for years that corresponds to the entire object’s shift from the origin. We recently introduced a new exact stationary solution, relevant for models with a single length parameter. This solution, referred to as the scale-invariant or dilation mode, has led to insights regarding the concept of perturbation energy within the linear theory framework. Aims. Our aim is to use Hénon’s isochrone model as an example to verify the ability of the standard matrix method to successfully predict the existence of a dilation mode, and to explore its potential application as a test disturbance. Methods. We used the standard matrix method for radial perturbations and applied Clutton-Brock potential-density pairs to determine the properties of the perturbations. Results. In this particular case of stationary radial perturbations, the typical relationship between the perturbations of the distribution function and the potential fails. This discrepancy poses a challenge when attempting to use the dilation mode as a test. When using Clutton-Brock pairs with the matrix method, a mass conservation equation as an additional equation to the ordinary set of linear equations is required. With this added equation, it’s possible to obtain the needed test: identical vanishing of the determinant of this modified set of equations with an increasing number of included basis functions.

SCALE-INVARIANT MODE IN COLLISIONLESS SPHERICAL STELLAR SYSTEMS

An analytical solution of the perturbed equations is obtained, which exists in all ergodic models of collisionless spherical stellar systems with a single length parameter. This solution corresponds to variations of this parameter, that is, the stretching/contraction of the sphere keeping the total mass constant. During the process, the system remains in an equilibrium state. The simplicity of the solution makes it possible to explicitly give expressions for the distribution function, potential, and density in all orders of perturbation theory. This, in turn, helps to clarify the concept of perturbation energy, which, being a second-order magnitude in amplitude, cannot be calculated in linear theory. The expression for the 2nd-order perturbation energy does not match the well-known quadratic form for perturbation energy derived from 1st-order perturbations in linear theory. However, both of these energies are integrals of motion and differ only by a constant. The obtained solution can be used to control the correctness of codes and the accuracy of calculations in the numerical study of collisionless stellar models.

Scale-Invariant Mode in Collisionless Spherical Stellar Systems

An analytical solution for the perturbed equations, applicable to all ergodic models of collisionless spherical stellar systems with a single length parameter, has been derived. This solution corresponds to variations in this parameter, i.e., the expansion or contraction of the sphere while conserving total mass. During this process, the system maintains an equilibrium state. The simplicity of the solution allows for the explicit expression of the distribution function, potential, and density across all orders of perturbation theory. This, in turn, aids in clarifying the concept of perturbation energy, which, being of second order in amplitude, cannot be calculated using linear theory. It is demonstrated that the correct expression for perturbation energy, accounting for second-order perturbations, does not align with the well-known expression for perturbation energy via a quadratic form, derived from first-order perturbations within linear theory. However, both these energies are integrals of motion and differ only by a constant. The derived solution can be utilized to verify the correctness of codes and the precision of calculations in the numerical study of collisionless stellar models.

New dipole instabilities in spherical stellar systems

ABSTRACT Spherical stellar systems have weakly damped response modes. The dipole modes are seiche modes. The quadrupole are zero pattern-speed prolate modes, the stable precursors to the radial orbit instability (ROI). We demonstrate that small wiggles in the distribution function (DF) can destabilize the dipole modes and describe the newly identified instabilities in NFW-like dark-matter (DM) haloes and other power-law spherical systems. The modes were identified in N-body simulations using multivariate singular spectrum analysis and corroborated using linear-response theory. The new mode peaks inside the half-mass radius but has a pattern speed typical of an outer-halo orbit. As it grows, the radial angle of the eccentric orbits that make up the mode correlates and loses angular momentum by a resonant couple to outer-halo orbits. This leads to an unsteady pattern with a density enhancement that swings from one side of the halo to another along a diameter, like the orbits that comprise the instability. In this way, the dipole mode is similar to the ROI. Since the DF found in Nature is unlikely to be smooth and isotropic with df(E)/dE < 0 necessary for Antonov stability, these modes may be ubiquitous albeit slowly growing. Haloes that are less extended than NFW, such as the Hernquist model, tend to be stable to this dipole instability. We present the critical stability exponents for one- and two-power models. These different critical outer power-law exponents illustrate that the gravitational coupling between the inner and outer DM halo depends on the global shape of density profile.

Stellar hardening of massive black hole binaries: the impact of the host rotation

ABSTRACT Massive black hole binaries (MBHBs) with masses of ∼104 to $\sim 10^{10} \, \mathrm{M}_{\odot {}}$ are one of the main targets for currently operating and forthcoming space-borne gravitational wave observatories. In this paper, we explore the effect of the stellar host rotation on the bound binary hardening efficiency, driven by three-body stellar interactions. As seen in previous studies, we find that the centre of mass (CoM) of a prograde MBHB embedded in a rotating environment starts moving on a nearly circular orbit about the centre of the system shortly after the MBHB binding. In our runs, the oscillation radius is ≈ 0.25 (≈ 0.1) times the binary influence radius for equal mass MBHBs (MBHBs with mass ratio 1:4). Conversely, retrograde binaries remain anchored about the centre of the host. The binary shrinking rate is twice as fast when the binary CoM exhibits a net orbital motion, owing to a more efficient loss cone repopulation even in our spherical stellar systems. We develop a model that captures the CoM oscillations of prograde binaries; we argue that the CoM angular momentum gain per time unit scales with the internal binary angular momentum, so that most of the displacement is induced by stellar interactions occurring around the time of MBHB binding, while the subsequent angular momentum enhancement gets eventually quashed by the effect of dynamical friction. The effect of the background rotation on the MBHB evolution may be relevant for LISA sources, that are expected to form in significantly rotating stellar systems.

Damped perturbations in stellar systems: genuine modes and Landau-damped waves

ABSTRACT This research was stimulated by the recent studies of damping solutions in dynamically stable spherical stellar systems. Using the simplest model of the homogeneous stellar medium, we discuss non-trivial features of stellar systems. Taking them into account will make it possible to correctly interpret the results obtained earlier and will help to set up decisive numerical experiments in the future. In particular, we compare the initial value problem versus the eigenvalue problem. It turns out that in the unstable regime, the Landau-damped waves can be represented as a superposition of van Kampen modes plus a discrete damped mode, usually ignored in the stability study. This mode is a solution complex conjugate to the unstable Jeans mode. In contrast, the Landau-damped waves are not genuine modes: in modes, eigenfunctions depend on time as exp (−iωt), while the waves do not have eigenfunctions on the real $v$-axis at all. However, ‘eigenfunctions’ on the complex $v$-contours do exist. Deviations from the Landau damping are common and can be due to singularities or cut-off of the initial perturbation above some fixed value in the velocity space.

Breaking beta: a comparison of mass modelling methods for spherical systems

ABSTRACT We apply four different mass modelling methods to a suite of publicly available mock data for spherical stellar systems. We focus on the recovery of the density and velocity anisotropy as a function of radius, either using line-of-sight velocity data only or adding proper motion data. All methods perform well on isotropic and tangentially anisotropic mock data, recovering the density and velocity anisotropy within their 95 per cent confidence intervals over the radial range 0.25 < R/R1/2 < 4, where R1/2 is the half-light radius. However, radially anisotropic mocks are more challenging. For line-of-sight data alone, only methods that use information about the shape of the velocity distribution function are able to break the degeneracy between the density profile and the velocity anisotropy, β, to obtain an unbiased estimate of both. This shape information can be obtained through directly fitting a global phase-space distribution function, by using higher order ‘virial shape parameters’ or by assuming a Gaussian velocity distribution function locally, but projecting it self-consistently along the line of sight. Including proper motion data yields further improvements, and in this case, all methods give a good recovery of both the radial density and velocity anisotropy profiles.

Relaxation of spherical stellar systems

ABSTRACTA total of 10 000 simulations of 1000 particle realizations of the same cluster are computed by direct force summation. After three crossing times self-gravity has amplified the original Poisson noise more than tenfold on large scales, and the amplification is still far from complete. The cluster’s fundamental dipole mode is strongly excited by Poisson noise, and this mode makes a major contribution to driving diffusion of stars in energy. The diffusive flow through action space is computed for the simulations and compared with the predictions of both Chandrasekhar’s local-scattering theory and the Balescu–Lenard (BL) equation. The predictions of local-scattering theory are qualitatively wrong because the latter neglects self-gravity. These results imply that local-scattering theory can account for only a fraction of a cluster’s relaxation. Future work on cluster evolution should employ either N-body simulation or the BL equation. However, significant code development will be required to make use of the BL equation practicable and the way forward may be to merge BL theory with local scattering theory so fluctuations of every scale are efficiently handled.

Resonant Relaxation in Globular Clusters

Abstract Resonant relaxation has been discussed as an efficient process that changes the angular momenta of stars orbiting around a central supermassive black hole due to the fluctuating gravitational field of the stellar cluster. Other spherical stellar systems, such as globular clusters, exhibit a restricted form of this effect where enhanced relaxation rate only occurs in the directions of the angular momentum vectors, but not in their magnitudes; this is called vector resonant relaxation (VRR). To explore this effect, we performed a large set of direct N-body simulations, with up to 512k particles (where k =1024) and ∼500 dynamical times. Contrasting these simulations, which naturally include the collective effects, with Spitzer-style Monte Carlo simulations, which by design only exhibit two-body relaxation, we show that the temporal behavior of the angular momentum vectors in N-body simulations cannot be explained by two-body relaxation alone. VRR operates efficiently in globular clusters with N > 104. The fact that VRR operates in globular clusters may open a way to use powerful tools in statistical physics for their description. In particular, since the distribution of orbital planes relaxes much more rapidly than the distribution of the magnitude of angular momentum and the radial action, the relaxation process reaches an internal statistical equilibrium in the corresponding part of phase space while the whole cluster is generally out of equilibrium, in a state of quenched disorder. We point out the need to include effects of VRR in Monte Carlo simulations of globular clusters.

Black Hole Disks in Galactic Nuclei.

Gravitational torques among objects orbiting a supermassive black hole drive the rapid reorientation of orbital planes in nuclear star clusters (NSCs), a process known as vector resonant relaxation. In this Letter, we determine the statistical equilibrium of systems with a distribution of masses, semimajor axes, and eccentricities. We average the interaction over the apsidal precession time and construct a MonteCarlo Markov chain method to sample the microcanonical ensemble of the NSC. We examine the case of NSCs formed by 16 episodes of star formation or globular cluster infall. We find that the massive stars and stellar mass black holes form a warped disk, while low mass stars resemble a spherical distribution with a possible net rotation. This explains the origin of the clockwise disk in the Galactic center and predicts a population of black holes (BHs) embedded within this structure. The rate of mergers among massive stars, tidal disruption events of massive stars by BHs, and BH-BH mergers are highly increased in such disks. The first two may explain the origin of the observed G1 and G2 clouds, the latter may be important for gravitational wave detections with LIGO and VIRGO. More generally, black holes are expected to settle in disks in all dense spherical stellar systems assembled by mergers of smaller systems including globular clusters.

How to break the density-anisotropy degeneracy in spherical stellar systems

We present a new non-parametric Jeans code, GravSphere, that recovers the density $\rho(r)$ and velocity anisotropy $\beta(r)$ of spherical stellar systems, assuming only that they are in a steady-state. Using a large suite of mock data, we confirm that with only line-of-sight velocity data, GravSphere provides a good estimate of the density at the projected stellar half mass radius, $\rho(R_{1/2})$, but is not able to measure $\rho(r)$ or $\beta(r)$, even with 10,000 tracer stars. We then test three popular methods for breaking this $\rho-\beta$ degeneracy: using multiple populations with different $R_{1/2}$; using higher order `Virial Shape Parameters' (VSPs); and including proper motion data. We find that two populations provide an excellent recovery of $\rho(r)$ in-between their respective $R_{1/2}$. However, even with a total of $\sim 7,000$ tracers, we are not able to well-constrain $\beta(r)$ for either population. By contrast, using 1000 tracers with higher order VSPs we are able to measure $\rho(r)$ over the range $0.5 < r/R_{1/2} < 2$ and broadly constrain $\beta(r)$. Including proper motion data for all stars gives an even better performance, with $\rho$ and $\beta$ well-measured over the range $0.25 < r/R_{1/2} < 4$. Finally, we test GravSphere on a triaxial mock galaxy that has axis ratios typical of a merger remnant, $[1:0.8:0.6]$. In this case, GravSphere can become slightly biased. However, we find that when this occurs the data are poorly fit, allowing us to detect when such departures from spherical symmetry become problematic.

Collisionless loss-cone refilling: there is no final parsec problem

Coalescing massive black hole binaries, formed during galaxy mergers, are expected to be a primary source of low frequency gravitational waves. Yet in isolated gas-free spherical stellar systems, the hardening of the binary stalls at parsec-scale separations owing to the inefficiency of relaxation-driven loss-cone refilling. Repopulation via collisionless orbit diffusion in triaxial systems is more efficient, but published simulation results are contradictory. While sustained hardening has been reported in simulations of galaxy mergers with $N \sim 10^6$ stars and in early simulations of rotating models, in isolated non-rotating triaxial models the hardening rate continues to fall with increasing N, a signature of spurious two-body relaxation. We present a novel approach for studying loss cone repopulation in galactic nuclei. Since loss cone repopulation in triaxial systems owes to orbit diffusion, it is a purely collisionless phenomenon and can be studied with an approximated force calculation technique, provided the force errors are well behaved and sufficiently small. We achieve this using an accurate fast multipole method and define a proxy for the hardening rate that depends only on stellar angular momenta. We find that the loss cone is efficiently replenished even in very mildly triaxial models (with axis ratios 1 : 0.9 : 0.8). Such triaxiality is unavoidable following galactic mergers and can drive binaries into the gravitational wave regime. We conclude that there is no final parsec problem.

Made-to-measure dark matter haloes, elliptical galaxies and dwarf galaxies in action coordinates

We provide a family of action-based distribution functions (DFs) for the double-power law family of densities often used to model galaxies. The DF itself is a double-power law in combinations of the actions, and reduces to the known limits in the case of a pure power-law at small and large radii. Our method enables the velocity anisotropy of the model to be tuned, and so the anisotropy in the inner and outer parts can be specified for the application in hand. We provide self-consistent DFs for the Hernquist and Jaffe models - both with everywhere isotropic velocity dispersions, and with kinematics that gradually becomes more radially anisotropic on moving outwards. We also carry out this exercise for a cored dark-matter model. These are tailored to represent dark haloes and elliptical galaxies respectively with kinematic properties inferred from simulations or observational data. Finally, we relax a cored luminous component within a dark matter halo to provide a self-consistent model of a dwarf spheroidal embedded in dark matter. The DFs provide us with non-rotating spherical stellar systems, but one of the virtues of working with actions is the relative ease with which such models can be converted into axisymmetry and triaxiality.

Radially anisotropic models of collisionless spherical stellar systems without central singularity

A new class of the simplest equilibrium two-parameter distribution functions for spherical stellar systems with a radially anisotropic stellar velocity distribution is investigated. The models under consideration are a less singular counterpart of the so-called generalized polytropes, which in the past were among the most popular models in works on the equilibrium and stability of gravitating systems. In contrast to the well-known generalized polytropes, the proposed models have finite density and potential at the center. The absence of a singularity is necessary for a proper consideration of the radial orbit instability, which is the most important instability of spherical stellar systems. The main observed parameters of the proposed models (potential, density, anisotropy) are compared with those in well-known equilibrium models.

Equilibrium models of radially anisotropic spherical stellar systems with softened central potentials

We study a new class of equilibrium two-parametric distribution functions of spherical stellar systems with radially anisotropic velocity distribution of stars. The models are less singular counterparts of the so called generalized polytropes, widely used in works on equilibrium and stability of gravitating systems in the past. The offered models, unlike the generalized polytropes, have finite density and potential in the center. The absence of the singularity is necessary for proper consideration of the radial orbit instability, which is the most important instability in spherical stellar systems. Comparison of the main observed parameters (potential, density, anisotropy) predicted by the present models and other popular equilibrium models is provided.

On the threshold anisotropy for unstable spherical stellar systems

The radial orbit instability generally arises in anisotropic collisionless stellar systems with the dominance of radial motions over transverse ones. Using the simplest anisotropic generalization of polytrope models for spherical clusters as an example, we show that the instability growth rates become exponentially small as the isotropic limit is approached. Given the finite lifetime of real astronomical objects, these systems should be assumed to become stable at some finite radial anisotropy.

Notes on the stability threshold for radially anisotropic polytropes

We discuss some contradictions found in the literature concerning the problem of stability of collisionless spherical stellar systems which are the simplest anisotropic generalization of the well-known polytrope models. Their distribution function $F(E,L)$ is a product of power-low functions of the energy $E$ and the angular momentum $L$, i.e. $F\propto L^{-s}(-E)^q$. On the one hand, calculation of the growth rates in the framework of linear stability theory and N-body simulations show that these systems become stable when the parameter $s$ characterizing the velocity anisotropy of the stellar distribution is lower than some finite threshold value, $s<s_\textrm{crit}$. On the other hand Palmer & Papaloizou (1987) showed that the instability remained up to the isotropic limit $s=0$. Using our method of determining the eigenmodes for stellar systems, we show that the growth rates in weakly radially-anisotropic systems are indeed positive, but decrease exponentially as the parameter $s$ approaches zero, i.e. $\gamma\propto \exp(-s_{\ast}/s)$. In fact, for the systems with finite lifetime this means stability.

Radial-orbit instability in modified Newtonian dynamics

The stability of radially anisotropic spherical stellar systems in modified Newtonian dynamics (MOND) is explored by means of numerical simulations performed with the N-body code N-MODY. We find that Osipkov-Merritt MOND models require for stability larger minimum anisotropy radius than equivalent Newtonian systems (ENSs) with dark matter, and also than purely baryonic Newtonian models with the same density profile. The maximum value for stability of the Fridman-Polyachenko-Shukhman parameter in MOND models is lower than in ENSs, but higher than in Newtonian models with no dark matter. We conclude that MOND systems are substantially more prone to radial-orbit instability than ENSs with dark matter, while they are able to support a larger amount of kinetic energy stored in radial orbits than purely baryonic Newtonian systems. An explanation of these results is attempted, and their relevance to the MOND interpretation of the observed kinematics of globular clusters, dwarf spheroidal and elliptical galaxies is briefly discussed.

Generalized Schwarzschild’s method

We describe a new finite element method (FEM) to construct continuous equilibrium distribution functions of stellar systems. The method is a generalization of Schwarzschild's orbit superposition method from the space of discrete functions to continuous ones. In contrast to Schwarzschild's method, FEM produces a continuous distribution function (DF) and satisfies the intra element continuity and Jeans equations. The method employs two finite-element meshes, one in configuration space and one in action space. The DF is represented by its values at the nodes of the action-space mesh and by interpolating functions inside the elements. The Galerkin projection of all equations that involve the DF leads to a linear system of equations, which can be solved for the nodal values of the DF using linear or quadratic programming, or other optimization methods. We illustrate the superior performance of FEM by constructing ergodic and anisotropic equilibrium DFs for spherical stellar systems (Hernquist models). We also show that explicitly constraining the DF by the Jeans equations leads to smoother and/or more accurate solutions with both Schwarzschild's method and FEM.

ALTERNATING DIRECTION IMPLICIT METHOD FOR TWO-DIMENSIONAL FOKKER-PLANCK EQUATION OF DENSE SPHERICAL STELLAR SYSTEMS

The Fokker-Planck (FP) model is one of the commonly used methods for studies of the dynamical evolution of dense spherical stellar systems such as globular clusters and galactic nuclei. The FP model is numerically stable in most cases, but we find that it encounters numerical difficulties rather often when the effects of tidal shocks are included in two-dimensional (energy and angular momentum space) version of the FP model or when the initial condition is extreme (e.g., a very large cluster mass and a small cluster radius). To avoid such a problem, we have developed a new integration scheme for a two-dimensional FP equation by adopting an Alternating Direction Implicit (ADI) method given in the Douglas-Rachford split form. We find that our ADI method reduces the computing time by a factor of <TEX>${\sim}2$</TEX> compared to the fully implicit method, and resolves problems of numerical instability.