Predicting Steady States of One-dimensional Collisionless Gravitating Systems

TL;DR

This paper introduces a novel Hermite-Legendre analysis technique to determine steady states in one-dimensional collisionless gravitating systems, aiding understanding of their relaxation and implications for dark matter structures.

Contribution

It presents a new method for identifying steady states and analyzing their stability, emphasizing the role of initial conditions in system evolution.

Findings

Physically motivated parameter ranges for steady states identified.

Any initial condition in the linear regime can be decomposed into unique modes.

Steady states can be described by combinations of time-independent modes.

Abstract

Building on the development of a Hermite-Legendre analysis of one-dimensional gravitating collisionless systems, we present a technique for determining the steady states of such systems. This provides an important component for understanding the physics involved in the relaxation of these kinds of systems. As the dark matter structures in the universe should have traits in common with these systems, insight into this relaxation can provide clues to larger astrophysical questions. For large perturbation strengths, we determine physically motivated parameter ranges for the simplest families of steady states as well as their stability. We also demonstrate that any set of initial conditions in the linear regime can be resolved into unique time-independent and time-dependent modes. Combinations of time-independent modes then describe the steady state of any system linearly perturbed from equilibrium. These results highlight the importance of initial conditions over relaxation mechanisms in the evolution of these systems.

Equations

$$\frac{\mathrm{d}f}{\mathrm{d}t}=\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}+a(x)\frac{\partial f}{\partial v}=0,$$

$$a(x)=-g{\int }_{-\infty }^x\lambda (s)\mathrm{d}s+g{\int }_x^{\infty }\lambda (s)\mathrm{d}s,$$

$$\lambda (x)=M{\int }_{-\infty }^{\infty }f(x,v)\mathrm{d}v.$$

$$ϵ=\frac{v^2}{2}+\varphi (x),$$

$$f_0(x,v)=A{\mathrm{sech}}^2(\frac{\beta g{M}_{\mathrm{total}}}{2}x)e^{-\frac{\beta v^2}{2}},$$

$${\varphi }_0(x)=\ln (2\cosh \frac{\beta g{M}_{\mathrm{total}}}{2}x).$$

$$f_0(ϵ)=A{e}^{-\beta ϵ}.$$

$$\chi =\frac{\beta g{M}_{\mathrm{total}}}{2}x,\varpi =\sqrt{\frac{\beta }{2}}v,\text{mbox}and\tau =\sqrt{\frac{\beta }{2}}g{M}_{\mathrm{total}}t.$$

$${\tilde{f}}_0(\chi ,\varpi )=\frac{2}{g{M}_{\mathrm{total}}}\sqrt{\frac{{\beta }^3}{2}}{f}_0=\frac{1}{2\sqrt{\pi }}{\mathrm{sech}}^2\chi e^{-{\varpi }^2}.$$

$$\frac{\partial \tilde{f}}{\partial \tau }+\varpi \frac{\partial \tilde{f}}{\partial \chi }+\alpha (\chi )\frac{\partial \tilde{f}}{\partial \varpi }=0,$$

$$\alpha (\chi )=-{\int }_{-\infty }^{\chi }\Lambda ({\chi }^{\prime })\mathrm{d}{\chi }^{\prime }+{\int }_{\chi }^{\infty }\Lambda ({\chi }^{\prime })\mathrm{d}{\chi }^{\prime },$$

$$\Lambda (\chi )={\int }_{-\infty }^{\infty }\tilde{f}(\chi ,\varpi )\mathrm{d}\varpi$$

$$\tilde{f}={\tilde{f}}_0+\Delta {\tilde{f}}_1,$$

$$\frac{\partial {\tilde{f}}_1}{\partial \tau }+\varpi \frac{\partial {\tilde{f}}_1}{\partial \chi }+{\alpha }_0(\chi )\frac{\partial {\tilde{f}}_1}{\partial \varpi }+{\alpha }_1(\chi )\frac{\partial {\tilde{f}}_0}{\partial \varpi }=0.$$

$${\alpha }_0(\chi )=-{\int }_{-\chi }^{\chi }{\int }_{-\infty }^{\infty }{\tilde{f}}_0(\chi ,\varpi )\mathrm{d}\varpi \mathrm{d}{\chi }^{\prime }=-\tanh \chi .$$

$$\frac{\partial {\tilde{f}}_0}{\partial \varpi }=-2\varpi {\tilde{f}}_0.$$

$${\tilde{f}}_1=\sum _{m,n}{c}_{m,n}\mu \nu H_m(\varpi )P_n(\tanh \chi ){\mathrm{sech}}^2\chi e^{-{\varpi }^2},$$

$${\alpha }_1(\chi )$$

$$\int P_n(s)\mathrm{d}s=\frac{P_{n+1}(s)-P_{n-1}(s)}{2n+1},$$

$$\frac{\partial {\tilde{f}}_1}{\partial \tau }+\varpi (1-u^2)\frac{\partial {\tilde{f}}_1}{\partial u}-u\frac{\partial {\tilde{f}}_1}{\partial \varpi }-$$

$${\dot{c}}_{m,n}$$

$$L_{m,n}^{m-1,n-1}$$

$$L_{m,n}^{m-1,n+1}$$

$$L_{m,n}^{m+1,n-1}$$

$$L_{m,n}^{m+1,n+1}$$

$$\tilde{f}=\sum _{m,n}{A}_{m,n}\mu \nu H_m(\varpi )P_n(\tanh \chi ){\mathrm{sech}}^2\chi e^{-{\varpi }^2}.$$

$$P_j{P}_k=\sum _{\genfrac{}{}{0ex}{}{s=0}{\mathrm{even}}}^{j+k}{Q}_s^{(j,k)}{P}_{j+k-s},$$

$$Q_s^{(j,k)}=\frac{2j+2k-2s+1}{2j+2k-s+1}\frac{{\lambda }_{s/2}{\lambda }_{j-s/2}{\lambda }_{k-s/2}}{{\lambda }_{j+k-s/2}},$$

$${\lambda }_B=\frac{(2B)!}{2^B(B!{)}^2},$$

$${\dot{A}}_{m,n}$$

$$S_1=2\sqrt{\sqrt{\pi }m(2n+1)}\times$$

$$S_2=2\sqrt{\sqrt{\pi }m(2n+1)}\times$$

Full text

1 Introduction

The current paradigm surrounding the formation of large-scale structure in the universe relies on the behavior of collisionless dark matter (e.g., Spergel et al., 2003; Springel et al., 2005). Investigations of three-dimensional systems, such as individual galactic-scale dark matter haloes, involve a wide range of evolutionary processes that contribute to the relaxation from initial conditions to a final equilibrium state (e.g., Navarro, Frenk, & White, Navarro *et al.*1996; Moore et al., 1998). The radial orbit instability (Merritt & Aguilar, 1985) along with evaporation and ejection (Binney & Tremaine, 1987) are commonly discussed examples of these processes. The sheer variety of processes occurring during the relaxation of a three-dimensional object significantly complicates any attempt to disentangle essential behaviors. Our overall goal is to illuminate the roles of phase mixing and violent relaxation in self-gravitating collisionless evolution.

To do this, we concern ourselves only with a one-dimensional version of a self-gravitating collisionless system. Dropping to a one-dimensional system has several advantages. Immediately, the evolution of the fine-grained distribution function $f$, which describes how many particles exist in infinitesimal volumes of phase space, is completely defined by the following collisionless Boltzmann (or Vlasov) equation,

[TABLE]

where $a(x)$ is the acceleration. This acceleration is a battle between mass to the left and mass to the right of any location,

[TABLE]

where $g$ is the gravitational coupling constant and $\lambda$ is the density distribution. For a system of mass $M$,

[TABLE]

The relative simplicity of this version of the Boltzmann equation provides a manageable starting point for analytical treatment. A two-dimensional phase-space structure is straightforward to visualize, and removes the need for surfaces of section or other techniques for analyzing higher dimensional spaces. At the same time, Section 1.1 presents the details of a one-dimensional, collisionless equilibrium with a separable form. Additionally, the simplicity of the phase space for these systems allows one to take advantage of highly accurate and efficient $N$-body simulation schemes. Unlike three-dimensional situations where issues such as softening lengths and potential-calculation parameters can blur insight regarding relaxation processes, simulations of one-dimensional systems have no free parameters and rely on simple kinematics with constant acceleration to evolve. The simplicity of the analysis of one-dimensional systems has led to decades of work. Barnes & Ragan (2014) has a brief discussion of these investigations. Attacks on systems far from equilibrium (e.g., Joyce & Worrakitpoonpon, 2011) and near equilibrium (e.g., Reidel & Miller, 1987) have also been undertaken. Decompositions using action-angle variables (Weinberg, 1991; Barré, Olivetti, & Yamaguchi, 2011) have also led to insights into the dynamics of these types of systems.

While the previously highlighted differences are positive for this work, we freely admit that what follows will be necessarily unrealistic. We work under the assumption that the generic characteristics of phase mixing and violent relaxation are independent of the dimensionality of the system. For the results presented here to have any usefulness in a wider context, these processes must be simply linked to the collisionless self-gravitating natures of the systems. Another shortcoming of this work is that predictive capabilities are confined to relatively small perturbations from equilibria. Initial conditions such as those that would more closely resemble cosmological conditions evolve non-linearly and lead to families of time-independent solutions. Without linearity, such initial states cannot be uniquely decomposed into these solutions, nor is it possible to exclude certain families based on initial conditions.

With these qualifications in mind, what follows is a discussion of a method for finding time-independent solutions to Equation 1. The remainder of this introduction is devoted to reviewing the basics of the Hermite-Legendre expansion that underlies our analysis.

1.1 Separable Solution Equilibrium

Based on the structure of Equation 1, it is well known that any function of the specific energy,

[TABLE]

is a solution. We are specifically interested in the separable solution to Equation 1, which is commonly written as,

[TABLE]

where $\beta$ is an inverse energy, $g$ is the gravitational coupling constant, $M_{\rm total}$ is the total mass of the system, and $A=(gM_{\rm total}/4)\sqrt{\beta^{3}/2\pi}$ is the normalization constant. The energy scale $\beta$ is related to the total energy of the equilibrium, $E_{0}=3M_{\rm total}/2\beta$. For brevity, we will refer to this solution as the separable equilibrium from here on. With this distribution function, it can be shown that the corresponding potential is,

[TABLE]

The importance of this fact lies in its ability to transform Equation 5 into the Boltzmann distribution function,

[TABLE]

The Boltzmann nature of the one-dimensional self-gravitating equilibrium is a vital difference from the three-dimensional case. This simple form is key to the mathematical approach for dealing with perturbations to this equilibrium. Additionally, this form guarantees that the kinetic temperature of such an equilibrium is uniform and allows one to view $\beta$ as analogous to $kT$ in a collisional system’s distribution function.

For simplicity, we transform to dimensionless coordinates using the definitions,

[TABLE]

This leaves us to write the scaled equilibrium distribution function as,

[TABLE]

Equation 1 transforms to,

[TABLE]

where $\alpha(\chi)$ is the dimensionless acceleration function. From Equations 2 and 3, this acceleration is given by

[TABLE]

where

[TABLE]

is the dimensionless density.

1.2 Linear Perturbations

Our goal is to investigate the relaxation of systems initially not in equilibrium. A useful first step in this direction is to deal with linear perturbations to equilibrium,

[TABLE]

where $\tilde{f}_{1}$ is the perturbing function and $\Delta$ controls the perturbation strength. For linear perturbations, we will consider $\Delta\ll 1$. The remainder of this section is a brief review based on work in Barnes & Ragan (2014).

Using this perturbed $\tilde{f}$ in Equation 1 produces a modified Boltzmann equation for the perturbing function (in terms of the previously defined dimensionless quantities),

[TABLE]

Using Equation 10, it is straightforward to find that,

[TABLE]

It is also useful to note that,

[TABLE]

We continue by expressing the perturbing distribution function in terms of Hermite and Legendre functions,

[TABLE]

where $\mu=1/\sqrt{2^{m}\sqrt{\pi}m!}$ and $\nu=\sqrt{(2n+1)/2}$ are related to the normalization functions of Hermite and Legendre functions, respectively. Also, note that $c_{0,0}=0$ since the equilibrium has already been removed. With this identification, the perturbing acceleration can be written as

[TABLE]

where the $\tanh\chi$ argument of the Legendre polynomials has been omitted for simplicity.

The integrals in the perturbed acceleration can be performed if we take advantage of the following substitutions; $u=\tanh\chi$, $\operatorname{sech}^{2}\chi=1-u^{2}$, $\mathrm{d}u=(1-u^{2})\,\mathrm{d}\chi$. The integrals involving the Legendre function become,

[TABLE]

which reduces Equation 12 to

[TABLE]

Expanding the remaining $\tilde{f}_{1}$ functions in Equation 15 with Equation 13 and using the orthogonality of the Hermite and Legendre functions produces the following recursion relation version of the collisionless Boltzmann equation,

[TABLE]

The factors $L_{m,n}^{i,j}$, which can be arranged as matrix elements, are given by

[TABLE]

where $m,n,i,j\geq 0$. In a situation that does not involve the gravitational pull of the perturbation on equilibrium, one can imagine test particles moving in a perturbed potential. For such a case, the recursion relation is obtained by omitting the Kronecker $\delta_{1,m}$ terms.

1.3 Non-linear Perturbations

For large amplitude perturbations, the approach is similar to that for linear perturbations. We again decompose a distribution function into Hermite and Legendre polynomials,

[TABLE]

The significant difference from the linear case is the $\alpha\partial\tilde{f}/\partial\varpi$ term in Equation 9. In the non-linear case, a perturbation will act on itself as well. The result of this self-interaction is that a non-trivial triple product of Legendre polynomials appears. In the linear case, one of the Legendre functions in the triple product is $P_{1}(u)=u$ and the product can be handled more simply. We take advantage of the fact that,

[TABLE]

to reduce any triple product to a product that can be simplified using the Legendre orthonormality relationship. The $Q$ functions are defined by (Dougall, 1953),

[TABLE]

where

[TABLE]

if $B\geq 0$ and is zero otherwise.

With this complication, the recursion relation version of Equation 9 expands to,

[TABLE]

where

[TABLE]

and

[TABLE]

The matrix elements $R_{m,n}^{i,j}$ are just the test-particle versions of the $L_{m,n}^{i,j}$ for the linear case;

[TABLE]

The Kronecker delta terms of the linear case are simply single terms from the $S_{1}$ and $S_{2}$ sums when $s$ has its maximum value.

2 Linear Perturbation Time-independent Solutions

Based on the results of Section 1.2, we next discuss routes to steady states of the linearized collisionless Boltzmann equation in test-particle and self-gravitating regimes using coefficient recursion relations. For test particles, phase-mixing will be the only process active in the phase-space evolution of the system. As a result, any evolution will occur on the time-scale of phase-mixing. In the self-gravitating case, violent relaxation will occur as well. The last term on the left-hand side of Equation 12 makes the system self-gravitating, and the perturbed acceleration is the only term that can relate to this kind of relaxation.

The collisionless Boltzmann equation possesses an infinite set of steady-state solutions, as does its linearized version (Binney & Tremaine, 1987). In the following, we describe a procedure to obtain the general solution of the steady-state linear problem for perturbations of the separable equilibrium. The analysis yields a set of time-independent modes which form a complete orthonormal basis that span a sub-space of all possible configurations. These modes can then be used to construct any steady state. More importantly, the projection of an arbitrary small-amplitude initial perturbation onto this basis produces the steady state that would result from evolving the system according to Equation 12.

2.1 Recursion Relation Procedure

The form of the time-independent modes of the linear problem is suggested by the solutions of the test-particle case, which have the following form,

[TABLE]

where $\phi_{0}=\log(2\cosh\chi)$ is the external potential. For small deviations from equilibrium, the distribution function can be expanded in a power series in $\beta\epsilon$ times a Boltzmann kernel,

[TABLE]

where the $a_{k}$ are time-independent coefficients.

The set of linearly independent functions $\{\beta\epsilon(\chi,\varpi)^{k}\}$ can be rendered into an orthonormal basis $\{F^{(k)}(\chi,\varpi)\}$ via a Gram-Schmidt process, where the test-particle $\{F^{(k)}(\chi,\varpi)\}$ are $k$th order polynomials in $\varpi^{2}$ and $\phi_{0}(\chi)$. Likewise, the solutions of the self-gravitating Boltzmann equation can be written as

[TABLE]

The $b_{k}$ coefficients define the relative strengths of the various time-independent mode contributions to the perturbation distribution function.

We can further break the orthonormal $F^{(k)}$ functions into combinations of velocity and position polynomials,

[TABLE]

where the $G_{k,m}(\chi)$ can be written in terms of Legendre polynomials,

[TABLE]

For a given value of $k$, the $d^{(k)}_{m,n}$ coefficients obey the recursion relation in Equation 16. The difference between test-particle and self-gravitating coefficient recursion relations is what distinguishes the $F^{(k)}$ functions for the two cases. Since we are now looking at time-independent solutions, the time derivative term in Equation 16 must be set to zero. As an example, a test-particle simulation will have $d$-coefficient values that must obey this recursion relation,

[TABLE]

The problem of solving the time-independent Boltzmann equation is transformed into solving recursion relations on the $(m,n)$ grid in the region $0\leq m\leq k$, $n\geq 0$, where the $m=k$ row is set to zero, except $d^{(k)}_{k,0}$ which is left as a free parameter. Any coefficient with an odd $m$ or $n$ index must be zero as those coefficients give rise to non-zero center-of-mass position and/or velocity values that are incompatible with a time-independent state. Figure 1 illustrates an example layout of a $d^{(k)}$-coefficient grid. The coefficients on the left-most column $d^{(k)}_{m,0}$ are likewise left as free parameters for the subsequent Gram-Schmidt orthogonalization procedure. For a given $k$, one starts at the upper left-hand corner ($m=k$) and uses the appropriate recursion relation to find the coefficients for $(m=k-2,n>0)$ up to some $n_{max}$, working left to right. The cut-off $n_{\rm max}$ is chosen so that the series for $G_{k,m}$ is well-behaved (see § 2.2). Once the $m=k-2$ row has been completed, the process can be repeated for all $k-2>m>0$, working downward. Once the coefficients for each $k$ are determined, the free parameters are used to construct an orthonormal basis via a Gram-Schmidt process.

This is essentially the path we follow, with a few important details to be added. Thinking of the dynamics problem in general, we re-cast the coefficient Boltzmann equation (Equation 16) as a matrix equation. First, arrange the $(0\leq m\leq k,0\leq n\leq n_{\rm max})$ $\dot{d}^{(k)}_{m,n}$ and $d^{(k)}_{m,n}$ terms involved in Equation 16 as vectors. The $L^{i,j}_{m,n}$ factors can then be organized into a two-dimensional matrix. This results in the following relationship,

[TABLE]

Assuming that the time-dependence of a $\mn@boldsymbol{d}^{(k)}$ is given by $\exp{(\lambda^{(k)}t)}$, then Equation 28 transforms into an eigenvalue equation,

[TABLE]

Solving this equation is straightforward, but a complication arises that impacts any subsequent Gram-Schmidt process. The matrix L is not symmetric. This means that there are so-called left- and right-handed eigenfunctions (sets of $d^{(k)}_{m,n}$), but left/right eigenfunctions are only orthogonal to their opposite-handed counterparts. As the Gram-Schmidt process relies on orthogonal functions, we need to find both right- and left-handed eigenfunctions for $\lambda^{(k)}=0$. This means we must also find eigenfunctions of $\textbf{{L}}^{\rm T}$. With both left- and right-handed eigenfunctions, the Gram-Schmidt process can proceed as long as left/right pairs of functions are used for orthonormality. We note that the same approach could be taken to determine time-dependent solutions ($\lambda^{(k)}\neq 0$). A full discussion of such solutions will be postponed to maintain focus on the process of determining steady states, however the general behavior of any $(m,n)$ perturbation is to couple to higher-index coefficients (Barnes & Ragan, 2014). Essentially, any initial condition contains an infinite number of time-dependent modes. These modes then quickly phase mix, leaving only the time-independent modes as the visible remnant of the initial conditions.

As examples of the procedure, we consider the determination of the first two time-independent modes. The first mode is built from an initially unknown normalization constant $d^{(2)}_{2,0}$. The right- and left-handed coefficient sets are found by applying the L and $\textbf{{L}}^{\rm T}$ recursion relation operations, respectively. Recall that in our scheme, all coefficients with $(m\geq 2,n>0)$ are zero at this stage – only the coefficients with $m=0$ and $n\geq 2$ ($n$ even) will be non-zero. Equation 27 with $m=1,n=1$ gives the link between the normalization constant and the first $m=0$ coefficient. For $n>1$, the recursion relation of Equation 27,

[TABLE]

links the remaining coefficient values, which are then known in terms of the normalization constant. The value of this constant is determined by the orthogonality relationship between the left- and right-handed coefficients,

[TABLE]

where the $L$ and $R$ superscripts on the coefficients indicate their handedness. At this point, we have a time-independent mode (set of coefficients) that is orthonormal to equilibrium.

To continue, we allow for two initially undetermined constants, $d^{(4)}_{4,0}$ and $d^{(4)}_{2,0}$. We again apply the L and $\textbf{{L}}^{\rm T}$ recursion relations to determine right- and left-handed coefficient sets in terms of the uknown constants. At this stage all coefficients with $(m\geq 4,n>0)$ are zero. The orthogonality relation between the first and second mode coefficient sets,

[TABLE]

allows us to determine the $d^{(4)}_{2,0}$ value. Finally, the normalization condition provides us with the condition to find $d^{(4)}_{4,0}$. This procedure continues similarly for higher-order time-independent mode calculations. Normalization provides one of the free parameters, while orthogonality with the previous functions sets the remainder.

The importance of these time-independent modes lies in the ability to predict a steady state based on initial conditions. If a system is gently perturbed from the $f_{0}$ equilibrium, the perturbations evolve according to the linearized Boltzmann equation by dephasing in phase-space and, if self-gravity is present, by transferring particles into and out of equilibrium. At any point in time in an evolution, one can imagine the system being composed of a steady-state component and a decaying, fluctuating component. In the linear problem, the final state can be predicted by projecting the initial conditions into the time-invariant sub-space. A time-independent coefficient $b_{k}$ is found simply by taking the inner product of the initial conditions with $F^{(k)}$,

[TABLE]

For $N$-body distribution functions composed of delta functions, the coefficients can be calculated from the average value of $F^{(k)}$,

[TABLE]

Note that the $F^{(k)}$ functions involved in Equations 33 and 34 can be either left- or right-handed, as they depend on the likewise handed $d^{(k)}_{m,n}$ values through Equations 25 and 26. The handedness chosen for the $F^{(k)}$ in these expressions must be opposite to that assumed for the distribution function in Equation 24. For concreteness in what follows, we have assumed that the perturbation distribution function is composed of right-handed $F^{(k)}$ functions, which then demands that the orthonormal left-handed $F^{(k)}$ be used to calculate time-independent mode coefficients.

2.2 Limitations

In order for our Gram-Schmidt approach to be implemented, we have to truncate infinite series, i.e., recursion relations need to be solved over finite regions of coefficient space. Another example where truncation plays a role, the $G_{k,m}(\chi)$ functions of Equation 26 must be approximated using a finite number of terms. If we focus on only $m=0$ terms, we are looking at the essential behavior of densities associated with the various time-independent modes. The actual densities of the modes involve multiplying by a $\operatorname{sech}^{2}\chi$ term, so the behavior at large $\chi$ is effectively killed. Figure 2 shows the impact of changing mode and maximum $n$ value on these curves. Through trial and error, we have settled on $n_{\rm max}=400$ as an acceptable limit for this work. We have also set $k_{\rm max}=16$ as the range of time-independent modes created. In the end, these limits have been adopted because of the success the scheme has had in describing the results of $N$-body simulations (see Section 2.6).

2.3 Analytical Comparison

As a check on the scheme described above, we compare our set of $d^{(2)}_{m,n}$ with coefficients found through a different route. We imagine changing the temperature of an equilibrium, related to $\beta$ from Section 1, by a small amount. This should produce another equilibrium, which should be time-independent. Thinking of this as a series expansion,

[TABLE]

to first order in the temperature change. The derivative is

[TABLE]

The term in square brackets is the perturbation distribution function. It is straightforward to show that this distribution function is time-independent, as intended. We have calculated the $(m=0,n)$ coefficient values that correspond to this perturbation. They are not normalized in the same way as the $d^{(2)}_{m,n}$ values, but their successive values have the same ratios ($m=0,n=2$ over $m=0,n=4$, for example) as for the $d^{(2)}_{0,n}$. In other words, our recursion relation approach reproduces a known time-independent mode.

2.4 Energy Characteristics

With the linear time-independent modes identified, we next report on their energy properties. Here, we focus on self-gravitating situations, as test particle systems do not have interesting energy behaviors. The kinetic energy content of any mode is determined solely by the value of the $d^{(k)}_{2,0}$ coefficient. This simple form results from the fact that calculating the second velocity moment of any distribution function that is expanded as in Equation 13 is non-zero only when $m=2$ and $n=0$. The potential energy content of any mode involves only coefficients with $m=0$ and even $n>0$, the same as acceleration (Barnes & Ragan, 2014).

With the $d^{(k)}_{m,n}$ coefficients from above, we find that only the $F^{(2)}$ mode contains energy. This is reasonable, as this mode corresponds to changing the temperature of the system. All other $F^{(k)}$ with $k\geq 4$ have kinetic and potential energies that are equal in magnitude and opposite in sign. The values also indicate that all $F^{(k)}$ modes are in virial equilibrium. Again, this is unsurprising as it must be time-independent. The connection with these functions and energy lead us to refer to these functions as $E$ modes.

For general linear perturbations, there will be some energy-bearing component and some non-energetic components. It is important to note that the non-energetic components can still affect the spatial and velocity density of a system. As a result, the structure of a steady state composed of a combination of time-independent modes is uniquely determined by its initial conditions, and not by some general principle such as entropy maximization.

2.5 Alternative Approaches

The scheme we have laid out here is not unique. It is convenient because the coefficients involved in each mode have simple links to quantities like kinetic and potential energy. However, one could choose to form different time-independent modes. For example, if all coefficients with $m>2,n>2$ are set equal to zero and $d_{0,2}$ is left as an undetermined constant, then the recursion relations can be used to calculate $d_{m,0}$ values, for $m\leq m_{\rm max}$. In essence, Figure 1 could be flipped about the $m=n$ diagonal.

The recursion relations change for this approach, but the analogue to L remains non-symmetric. As a result, both left- and right-handed eigenfunctions must be determined as before, and the Gram-Schmidt technique involving both functions must be employed. To contrast with the $E$ modes described above, we refer to these alternative functions as $B$ modes. As distinction from the $b_{k}$ values for $E$ modes, the time-independent coefficients related to $B$ modes are labeled as $y_{k}$.

2.6 Simulations

Test particle evolutions in the equilibrium potential use an adaptive time step, Runge-Kutta scheme to track particles. Particle accelerations are determined by the equilibrium potential only. Tolerances are chosen so that the total energies of test particle systems experience fractional variations on the order of $10^{-11}$.

When allowing the perturbation to self-consistently evolve, substantially more care must be taken with a simulation. Fortunately, the distance-independence of the gravitational force in the one-dimensional problem allows one to take advantage of a key simplification. During an evolution, all particles move with constant acceleration between crossings. As a result, kinematic equations precisely predict the motions of particles (only numerical round-off errors degrade the process) and total system energies vary by approximately $10^{-9}$ during thousand-particle evolutions over a thousand dynamical times.

Rather than setting a fixed time step for numerical evolution, the conditions of the simulation determine when particle positions and velocities are updated. Initially, the time until the next collision of nearest neighbor pairs is calculated for every pair. The pair with the shortest interval sets the time step and the two particles that exist at a common location switch their constant acceleration values as they pass one another. By keeping track of a particle’s last and next crossing times, only a few particle must be updated after a time step. The bookkeeping is made easier when crossing times are stored in a heap structure that can be quickly re-sorted (Noullez et al., 2003; Joyce & Worrakitpoonpon, 2011).

Initial conditions for simulations are created by exciting specific modes to perturb the separable equilibrium. In practice, this process requires some caution during implementation. Perturbing modes can involve negative distribution function values, at least for some values of $\chi$ and $\varpi$. Since we do not have a simple way of incorporating negative masses into our simulations, care must be taken with the amplitudes of any such modes. If one cavalierly assigns a single mode amplitude, other, unintended modes can appear in the following manner. Any generated $N$-body initial distribution function is max[0,$\tilde{f}(\chi,\varpi)$], not $\tilde{f}(\chi,\varpi)$, which can be negative. The unintended modes are those that are needed to make $\tilde{f}\geq 0$. Unless otherwise noted, we have fixed perturbations strengths at values that render this problem negligible.

Our simulations consist of ensembles of 100 distinct realizations of a given initial distribution function, each with $N=1024$ particles. Each realization is evolved independently and ensemble-averaged quantities are then created. Typically, evolutions end at $\tau=5T$, where $T$ is the crossing time-scale for a constant-density system with mass $M_{\rm total}$. That this time range is adequate to guarantee that simulations reach steady states will be made clear in the following discussion. Our evolution code tracks quantities like energies (kinetic, potential, total), coefficient values, and entropy. Entropy in these $N$-body simulations is calculated using a particle counting scheme,

[TABLE]

where $n$ is the particle count and $i$ enumerates different areas of phase space (all of size $\Delta\chi\Delta\varpi$). Unlike in quantum situations where $\Delta\chi$ and $\Delta\varpi$ can be related to Planck’s constant, we have simply used trial and error to set sizes of the phase-space boxes. After investigating a wide range, we have found that values near the adopted $\Delta\chi=\Delta\varpi=2\times 10^{-2}$ produce entropy values that show the most obvious changes during evolution. Smaller values result in almost no particles falling into the boxes, while larger values produce boxes so large that variation is basically absent. In either case, resulting entropy changes are small.

To investigate how well our time-independent modes describe steady states, we have run several suites of simulations. Initial conditions consist of simple perturbations to equilibrium due to single coefficients; $c_{2,0}$, $c_{0,2}$, and $c_{2,2}$. Note that these coefficients are distinct from the time-independent coefficients we denote by $d^{(k)}_{m,n}$. For each perturbation, we have varied the strength, $0.05\leq\Delta\leq 0.30$. In this way, we map the range of steady states that are well-described by the time-independent modes.

Figure 3 compares predictions from our time-independent $E$ mode coefficient sets to the outcome of self-gravitating simulations with $c_{0,2}$ perturbations. To simplify this discussion, unless otherwise specified, time-independent coefficients discussed will be those for $E$ modes. Solid lines show how the $b_{k}$ coefficients should vary with $\Delta$ for linear perturbations. These examples extend only to $k=8$ for brevity, but similar plots up to $k=16$ show similar levels of agreement. Note that we have multiplied the perturbation strengths by 10 for the horizontal axes and the coefficient values by 100 for the vertical axes. The open circles are centered at the initial values of the coefficients while the crosses indicate final values. In general, there is good agreement between the predictions and simulated results. Unsurprisingly, as the perturbation strength grows, the separation between the initial and final values grows. Figure 4 is the $B$ mode analogue to Figure 3. As with the $E$ modes, the $B$ mode coefficients derived from simulations match predictions well.

Figure 5 is a focused version of panel b from Figure 3. In this figure, the thick error bars show the error-in-the-mean range. The thin error bars represent the full range of coefficient values for an ensemble. The changes in coefficient values during an evolution and the overall ranges in the coefficient values in an ensemble both grow with perturbation strength. The linear assumption underlying our prediction line is breaking down at the highest perturbation strengths investigated here. Figure 6 shows coefficient behaviors as a function of time for a perturbation strength $\Delta=0.25$. The changes in coefficient values early in the evolution highlight the non-linearity of this situation.

For comparison, the results of self-gravitating simulations with initial $c_{2,0}$ perturbations are shown in Figure 7. The same basic agreement between predictions and simulations is evident, and the discrepancies set in around the same perturbation strength as previously noted.

Looking in more detail at the $b_{2}$ values resulting from a $c_{2,0}$ perturbation with different strengths shows how simulation non-linearities impact the coefficient values. Figure 8 shows time evolutions of ensemble-averaged $b_{2}$ coefficients over the range of perturbations shown in Figure 7. The thin lines bounding the various evolutions indicate the ensemble error-in-the-mean ranges for each set of simulations. As the perturbation strength increases, the non-linearity of the simulations increases, but stays roughly within the statistical uncertainty of the coefficients. Higher time-independent-mode coefficient evolutions can be noisier than those for $b_{2}$, but overall any non-linearity due to the $N$-body nature of the simulations can be considered small.

As noted earlier, increasing perturbation strength can lead to unintended modes being populated in a simulation. A good example of this occurs in self-gravitating simulations with initial $c_{2,2}$ perturbations. In these systems, there should be no possibility of having the first time-independent mode (with coefficient $b_{2}$). Figure 9 shows that for the lowest perturbation strengths, this is reasonably achieved. However, for even modest strengths ($\Delta=0.15$) we see this first time-independent mode appearing in Figure 9a. This is a consequence of higher $m,n$ perturbations causing negative distribution functions at lower perturbation strengths. Subsequently, our simulations leave the linear regime for smaller $\Delta$ compared to those for $m=0,n=2$ and $m=2,n=0$ cases.

All of the previous discussion has involved self-gravitating simulations. Analyses of test-particle simulations yield very similar results. Figures 10 and 11 are the test-particle analogues to Figures 7 and 8, respectively. Note that even the small amount of non-linearity in the self-gravitating case is absent.

3 Non-Linear Time-independent Solutions

As mentioned previously, we need to discuss families of time-independent solutions for the non-linear perturbation case. The specific modes identified for the linear case result because time-independent coefficients values ($d^{(k)}_{2,0}$, $d^{(k)}_{4,0}$, etc.) can be fixed via orthonormality.

Due to the non-linear nature of Equation 19, the procedure for calculating time-independent $A$ coefficients changes. First, the highest $m$ row for the solution cannot be truncated after the $n=0$ term. Second, an iterative approach needs to be taken. This is analogous to a relaxation approach to solve Poisson’s equation in two-dimensions (Press et al., 1994). Families of solutions are determined by choosing the maximum $m$ that will be allowed (denoted by $k$ in analogy to the linear case), and family members are distinguished by the value of $A_{k,0}$. As with the linear perturbation results, we will indicate the solution family with a superscripted index, $A^{(k)}$.

As before, we use the specific example of the $A^{(2)}$ family to illustrate the process. Initially, $A_{2,0}^{(2)}$ is the only non-zero coefficient. All $A_{m>2,n}^{(2)}$ are zero and will remain so. Additionally, any odd-parity coefficients will be zero and fixed as well. We use the fact that $\dot{A}_{1,1}^{(2)}=0$ and $\dot{A}_{3,1}^{(2)}=0$ to start finding time-independent coefficients. Equation 19 with the conditions given results in two equations that only involve $A_{0,2}^{(2)}$, $A_{2,0}^{(2)}$, and $A_{2,2}^{(2)}$. Again, any higher $n$ terms are assumed to be zero at this point. The two equations can be solved simultaneously, giving first estimates of $A_{0,2}^{(2)}$ and $A_{2,2}^{(2)}$. Next, we use $\dot{A}_{1,3}^{(2)}=0$ and $\dot{A}_{3,3}^{(2)}=0$. With our $n=2$ coefficient estimates, we can solve the resulting equations for $A_{0,4}^{(2)}$ and $A_{2,4}^{(2)}$. In this manner, estimates for the time-independent coefficients can be found up to some $n_{\rm max}$. Once the $n_{\rm max}$ passes are completed, the next iteration begins again with $\dot{A}_{1,1}^{(2)}=0$ and $\dot{A}_{3,1}^{(2)}=0$. The now non-zero coefficients at higher $n$ values enter non-linearly and affect the new estimates of $A_{0,2}^{(2)}$ and $A_{2,2}^{(2)}$. Marching back out to $n_{\rm max}$ likewise updates all other coefficients. Repeating iterations, the true time-independent coefficient values are approached.

In practice, we have also implemented the same kind of numerical dissipation that one would use in a relaxation Poisson solver. At the end of an iteration, coefficient values are reset to the average of the current and previous sets. We have found that this technique reduces required iterations by at least a factor of two, for a given level of convergence. After 10 iterations, we find that the average change in coefficient values is less than one percent for the $m=k$ set and 1-2 orders of magnitude smaller for the lower $m$ sets. Likewise, $\dot{A}_{m,n}^{(2)}=\mathcal{O}(10^{-6})$ for all $m$ and $n$. Finally, the virial ratio for our solutions, $2K/U$, is typically $1+\varepsilon$, where $|\varepsilon|$ is $\mathcal{O}(10^{-4})$, but this does increase with increasing $A_{2,0}^{(2)}$.

3.1 Solution Behavior

Typically, $n_{\rm max}=64$, but values up to 256 have been used successfully. The major impact of increasing $n_{\rm max}$ is to smooth the central regions of the distribution function. The density curves in Figure 12 highlight the impact of $n_{\rm max}$ and perturbation strength $\Delta=A_{2,0}^{(2)}/A_{0,0}$. The top panel compares the equilibrium density to a time-independent density distribution with $\Delta=0.5$ and $n_{\rm max}=64$. The small variations seen near $\chi=0$ become magnified in the middle panel as $\Delta$ increases to 1.0. By increasing $n_{\rm max}$ to 128 in the bottom panel, the central oscillations are made smaller.

We have extended this scheme to also create higher-order families. As an example, we describe how the scheme changes by examining the $A^{(4)}$ family. This is a two parameter family described by $A_{2,0}^{(4)}$ and $A_{4,0}^{(4)}$. With this family, three simultaneous equations need to be solved; $\dot{A}_{1,1}^{(4)}=0$, $\dot{A}_{3,1}^{(4)}=0$, and $\dot{A}_{5,1}^{(4)}=0$. While we have not investigated higher-order families in the same detail as the $A^{(2)}$ and $A^{(4)}$, we have successfully found solutions for $k\geq 6$ by solving $k/2+1$ equations simultaneously and following the general procedure. As with the $A^{(2)}$ family, we show a few representative density profiles in Figure 13. Unlike the linear perturbation case, the non-linear coupling in a large amplitude perturbation makes it impossible to determine from initial conditions which family the steady state will belong to. However, the total energy of the initial system could be used to select compatible family members. For example, an initial perturbation that makes the energy differ from the equilibrium value must evolve to a steady state with the same energy. Calculating energies for $A^{(2)}$ and $A^{(4)}$ solutions with positive kinetic energies and non-negative density distributions (see § 3.2) reveals that there are one-to-one correspondences between energy values and $\Delta$ values. For the $A^{(2)}$ family, the energy follows $\epsilon\approx 1.5+2.1\Delta$ from the equilibrium value $\epsilon=1.5$. With the $A^{(4)}$ family, this becomes $\epsilon\approx 1.5+2.1\Delta_{2}-0.01\Delta_{4}$. This weak dependence on $\Delta_{4}$ reflects that while the velocity distributions of these solutions vary substantially, their potential energies are nearly the same. For a fixed energy, one can determine the relationship between $\Delta_{2}$ and $\Delta_{4}$ that the steady state must have.

As previous work has found that the Lynden-Bell distribution function has some success in describing steady states (Joyce & Worrakitpoonpon, 2011), we note that none of the $A^{(2)}$ family members investigated here closely resemble the Lynden-Bell form as $A^{(2)}$ densities have more extended core structure. As a quick summary of the Lynden-Bell distribution function, it serves as a distinguishable particle counterpart to the Fermi-Dirac distribution function (Lynden-Bell, 1967),

[TABLE]

where $\eta$ is a normalization constant and $\mu$ is the chemical potential which parameterizes the distribution. With this, it is straightforward to find the density as a function of potential, but to get $\Lambda(\chi)$, Poisson’s equation must be solved. For a $\beta\mu$ value, we numerically determine Lynden-Bell potential and density distributions.

We have decomposed the difference between a normalized Lynden-Bell distribution function and the equilibrium distribution function into coefficients. By removing equilibrium, we focus on just the perturbation that the Lynden-Bell function represents. We find that the $n=0$ coefficients decrease in magnitude with increasing $m$. As a result, we expect that higher-order time-independent functions with $\Delta_{m}$ values that decrease as $m$ increases should be able to provide more accurate matches to a Lynden-Bell distribution function. Without doing an optimized search, we have found a decent approximation for a Lynden-Bell density, with $\beta\mu=1$, in an $A^{(4)}$ solution with $\Delta_{2}=0.6$ and $\Delta_{4}=0.15$. Figure 14 shows comparisons between the $A^{(2)}$ ($\Delta=0.5$), the $A^{(4)}$ mentioned above, a higher-order $A^{(8)}$ solution, and the Lynden-Bell densities. The progression in the figure supports the conjecture that Lynden-Bell equilibria are examples of these time-independent solutions.

3.2 Family Boundaries

The $A^{(2)}$ family of solutions extends from $\Delta=-1/\sqrt{2}$ to $\Delta\gg 1$. The lower limit is set by the point at which the kinetic energy of the system is zero. As pointed out in Barnes & Ragan (2014), the kinetic energy is simply related to the $A_{2,0}$ coefficient,

[TABLE]

The upper limit to $\Delta$ is undetermined from our explorations. As $\Delta$ is increased from 1 to 2, we find that the minimum density (which occurs at $\chi=0$) appears to slowly, possibly exponentially, approach zero.

As the $A^{(4)}$ family has two parameters, we have investigated how $A_{2,0}^{(4)}$ and $A_{4,0}^{(4)}$ interact. To isolate the impact of the $A_{4,0}^{(4)}$ value, we initially set $A_{2,0}^{(4)}$ to zero. In this situation, any negative value of $A_{4,0}^{(4)}$ results in negative density values. Similarly, $\Delta_{4}\ga 0.7$ also produces negative densities. As $\Delta_{2}$ is changed from zero, the upper limiting value of $\Delta_{4}$ also changes. Figure 15 shows the approximate boundary for non-zero density distributions for the $A^{(4)}$ family. The points mark locations where a rough grid search of parameter space result in negative densities. The thick, solid lines result from linear fits to the two sets of dots, while the hatched area denotes $\Delta_{2}$ and $\Delta_{4}$ values that produce positive densities at any location. The kinetic energy still only depends on the $A_{2,0}^{(4)}$ term, so there remains the same limit of $\Delta_{2}=-1/\sqrt{2}$ for positive kinetic energy. However, we have not found any system with continuously positive density (for all $\chi$) when $\Delta_{2}\la-0.3$. As mentioned previously, numerical simulations can settle down into states that are nearly, but not quite, Lynden-Bell equilibria. From the point of view of this work, this behavior reflects that there is a parameter sub-space in Figure 15 that not only guarantees non-negative densities, but also produces systems similar to Lynden-Bell models ($\Delta_{4}<\Delta_{2}$).

3.3 Stability Analysis

In addition to identifying these equilibria, we have also investigated their stability. Rather than producing a linearized collisionless Boltzmann equation about the $f_{0}$ equilibrium, we have linearized about an non-linear steady state. This produces a recursion relation for the coefficients of these linear perturbations $d_{m,n}$. It is very similar in structure to Equation 19,

[TABLE]

The $R$ matrix elements maintain the same form as before, but in this case the $V_{1}$ and $V_{2}$ terms remain linear in the unknown $d_{m,n}$ coefficients,

[TABLE]

and

[TABLE]

Similar to the tactic described in Section 2.1, a matrix equation can be developed for these recursion relations. The eigenvalues of the matrix that represents the right-hand side of Equation 38 then provide information about the stability of the equilibrium described by the $A_{m,n}$.

We have tracked the maximum real eigenvalue $\lambda_{\rm max}$ for several members of the $A^{(2)}$ family, $0.1\leq\Delta\leq 1.0$ with $n_{\rm max}=32$, 64, and 128. Figure 16 shows the general increase in $\lambda_{\rm max}$ with increasing $\Delta$ across the $n_{\rm max}$ values. For this figure, we have set the size of the perturbation coefficient system based on $n_{\rm max}$. The R+V matrix that represents the non-linear recursion relations is square, with optimal dimension $J_{\rm max}=(n_{\rm max}/2)^{2}$. By restricting the size of the R+V matrix for a fixed $n_{\rm max}$, we have found that the differences in the Figure 16 curves is due to the size of the R+V matrix. For example, using an $n_{\rm max}=64$ solution with a R+V matrix with size $(16)^{2}$, not $(32)^{2}$, produces a stability curve that looks the same as the $n_{\rm max}=32$ curve in Figure 16. As a result, it appears that in the large $J_{\rm max}$ (and $n_{\rm max}$) limit, these solutions become more and more stable.

These results have been compared with $N$-body simulations like those described earlier. Initial conditions are drawn from $A^{(2)}$ solutions, and the systems are allowed to evolve for 100 $T$. Figure 17 shows the entropy evolutions of several $N$-body ensembles based on $n_{\rm max}=64$ solutions. The $N$-body curves shown represent time-independent solutions with $\Delta$ values around 0.5. The dichotomy of $N$-body behaviors about $\Delta=0.5$ is more reminiscent of the behavior seen for the $J_{\rm max}=(16)^{2}$ curve in Figure 16. We suggest that the impact of finite particle numbers degrades the equivalent modal resolution. A point in support of this is the fact that initial $N$-body density distributions based on $n_{\rm max}=32$, 64, and 128 solutions are essentially indistinguishable. Figure 18a shows the ensemble-average initial spatial distribution of particles when an $A^{(2)}$ solution with $\Delta=0.5$ and $n_{\rm max}=128$ is used. Figure 18b shows differences in the number of particles per bin when $n_{\rm max}=32$ or $n_{\rm max}=64$ solutions are used. To focus on variations due only to $n_{\rm max}$, the random number sequence used for each ensemble has been held fixed. The thin lines in Figure 18 indicate the error in the mean for the ensemble average at each bin (which are also the error bars in Figure 18a). Essentially, our $N$-body simulations cannot resolve the underlying distribution function well enough to allow the stability differences to appear. Further support of this idea is given in Figure 19. For an $A^{(2)}$ solution with $\Delta=1.0$ and $n_{\rm max}=128$, we have run additional $N$-body ensembles with $N=2048$ and $N=4096$. The slower rise in entropy with larger particle number is in line with our hypothesis that more particles result in a higher fidelity time-independent solution representation, which results in less instability. We have also created a $N=2048$ ensemble with $n_{\rm max}=64$ to confirm that this behavior occurs generally.

The stability analysis can be extended to members of the $A^{(4)}$ family as well. While we have not done an exhaustive search over the $(\Delta_{2},\Delta_{4})$ plane, we find that systems that would lie within the hatched region of Figure 15 have maximum positive eigenvalues that are approximately an order of magnitude smaller than $\lambda_{\rm max}$ for a system that lies outside the hatched region. Specifically, the Lynden-Bell-like system with $\Delta_{2}=0.6$ and $\Delta_{4}=0.15$ and a mildly perturbed system with $\Delta_{2}=\Delta_{4}=0.2$ have $\lambda_{\rm max}\approx 1\times 10^{-4}$ and produce stable $N$-body systems. On the other hand, the $\Delta_{2}=0.2$ and $\Delta_{4}=-0.4$ model is not stable in an $N$-body simulation and has $\lambda_{\rm max}\approx 1\times 10^{-3}$. From the values in Figure 16, a factor of 10 indicates a significant difference in stability.

Non-linear solutions are approximated up to some $n_{\rm max}$. Our results indicate that those solutions are stable in the limit $n_{\rm max}\rightarrow\infty$. $N$-body realizations of these non-linear steady states appear to have a perturbation amplitude dependence (Figure 17), but we argue that this is a result of finite particle numbers. Modest numbers of particles do not allow simulated systems to capture the small, but important, differences between models with small and large $n_{\rm max}$. Simulations with increasing $N$ result in a clear trend towards increasing stability.

4 Summary

We have presented a procedure for calculating time-independent solutions of arbitrary perturbations of one-dimensional gravitating systems. Both test-particle and self-gravitating systems can be analyzed with this approach. Sets of coefficients that describe Hermite-Legendre polynomial products form the time-independent solutions. In the case of linear perturbations from equilibrium, the solutions are independent modes. Suites of highly efficient and accurate $N$-body simulations have been created to test predictions based on steady-state solutions.

Starting with linear perturbations, we find that there are two routes to determining these coefficient sets. For what we term $E$ modes, coefficient sets are limited in their Hermite index. These modes are directly related to energies of a system. As expected with time independence, all $E$ modes are in virial equilibrium. The first $E$ mode has non-zero energy and represents changing the temperature of the separable equilibrium state. The alternative $B$ modes are formed by limiting the Legendre index, but the lack of an obvious, related physical quantity makes them less appealing than their $E$-mode analogues. However, steady states can be predicted just as well using either set of modes.

For large amplitude perturbations, we follow only the $E$ mode approach. Non-linear terms in the collisionless Boltzmann equation necessitate an iterative solution approach to solving the coupled coefficient recursion relations. There are boundaries on the parameters of these solutions based on maintaining positive kinetic energy and continuously positive density values. We have found that a subset of these solutions are similar to Lynden-Bell models. Analyzing Lynden-Bell distribution functions suggests that increasing the maximum Hermite order of solutions allows for better approximations to the Lynden-Bell form. Analyzing the stability of non-linear perturbation solutions via coefficient dynamics and $N$-body simulations indicates that physically relevant (positive kinetic energy and density) steady states are stable. However, that stability can be upset by insufficient particle numbers in $N$-body simulations.

Arguably the most important use of these time-independent modes is the prediction of steady states from linear perturbation initial conditions. We find that for modest strength perturbations in self-gravitating systems, any non-linearities present are of the same order as statistical uncertainties and that the time-independent modes accurately predict the simulated steady states. Unfortunately, such an approach is not possible for non-linear initial conditions. At best, one might be able to use coefficient values from initial conditions to determine what families of time-independent solutions may be present.