{'fraction_non_alphanumeric': 0.022690992018244014, 'fraction_numerical': 0.013721018624097301, 'mean_word_length': 4.6393741292465975, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantifying Departures from Equilibrium with the Spherical Jeans Equation

Proper motions of collisionless pointlike objects in a spherically symmetric system, for example stars in a galaxy, can be used to test whether that system is in equilibrium, with no assumptions regarding isotropy. In particular, the fourth order spherical Jeans equation yields expressions for two observable quantities characterizing the departure from equilibrium, both of which can be expressed in terms of time derivatives of first and third moments of the velocities. As illustrations, we compute these quantities for tracer distributions drawn from an exact equilibrium configuration and also from near equilibrium configurations generated using the N-body code GALIC.


INTRODUCTION
The ΛCDM standard cosmological model successfully reproduces a wealth of observations at large scales (Komatsu et al. 2011;. At the kpc scale the confrontation of ΛCDM predictions with observations is more challenging. For example, ΛCDM makes the robust prediction that pure dark matter halos will have cuspy profiles (Navarro et al. 1997). However in Nature no pure dark matter halo has yet been observed, and the dark matter halos hosting galaxies may be deformed by baryonic physics.
Various studies disagree on just how much baryonic matter is necessary for such a deformation, although many studies (Governato et al. 2012;Oñorbe et al. 2015), but not all (Read et al. 2016), find that ultrafaint dwarf (UFD) galaxies have little enough baryonic matter so that pure dark matter simulations should reliably reproduce their profiles. However the dark matter profiles of UFDs are as of yet completely unknown.
In the case of the larger classical dwarf spheroidal (dSph) galaxies, various analyses of observations suggest cored (Goerdt et al. 2006;Walker and Penarrubia 2011) or cusped (Richardson and Fairbairn 2014) density profiles. Unfortunately ΛCDM simulations also differ in their predictions of whether these systems should have cusped or cored profiles, with the recent FIRE simulations (Oñorbe et al. 2015) suggesting that baryonic physics generically leads to cored halos, while Governato et al. (2012) suggests that only the heaviest are cored. On the other hand Illustris simulations (Taylor et al. 2016) find that even the Milky Way is well fit by a cusped NFW profile, at least at the radii considered, and Eagle simulations (Schaller et al. 2015) find that baryonic physics leads to an even cuspier profiles in the case of massive galaxies. Therefore even if observations some day are able to clearly distinguish a core from a cusp in classical dSphs, they will nonetheless not provide a foolproof test of ΛCDM.
In the next two decades, 30 meter class telescopes will revolutionize our understanding of the shapes of the dark matter halos of the Milky Way's satellite dSphs and UFDs by providing proper motions for hundreds to tens of thousands of members in each (Skidmore et al. 2015). These telescopes will be able to cleanly distinguish an inner core from a cusp (Evslin 2015) and may well be able to confront the ellipticities of the halos with simulations. Observation programs are being developed for these projects now, and for their planning one needs to know just which parameters will be determined.
Critical to all of these analyses is an understanding of the assumptions which go into the determination of the halo parameters. To date, only line of sight velocities have been obtained and analyzed for these systems. Three strong assumptions have been used so far in such analyses. First, it has often been assumed that the systems are spherically symmetric (Walker and Penarrubia 2011;Richardson and Fairbairn 2014) or at least axisymmetric (Jardel & Gebhardt 2012;Hayashi & Chiba 2015) despite the fact the CDM simulations predict triaxial halos (Frenk et al. 1988) and observations demonstrate that the stellar distribution is aspherical (Irwin and Hatzidimitriou 1995). Second, often assumptions are made concerning the isotropy of the velocities. In particular one often assumes that the anisotropy is independent of the radius (Walker et al. 2009) or that the radial dependence is of a given functional form (Merrifield and Kent 1990;Battaglia et al. 2011;Bonnivard et al. 2015a,b). Finally, it has been assumed that these systems are in equilibrium (Merrifield and Kent 1990;Richardson and Fairbairn 2014). Eventually all three of these assumptions will need to be relaxed. As a small step in this direction, in the present note we present a method which allows the third assumption to be tested given the proper motions of sufficiently many tracers. This test is entirely independent of the second assumption, however it still depends on the first. In the future, we intend to use the results presented below as a starting point for a treatment of the first assumption.
The goal of this program is to provide the necessary tools to resolve the cusp/core debate using the precise astrometry which will be available on future extremely large telescopes. This resolution can be achieved only if one is able to determine which systems are in equilibrium and also if and at what radii the equilibrium assumption breaks down. More precisely the establishment of equilibrium is necessary to turn these observations into a subhalo density profile, which in turn can teach one about the formation histories of the systems and about the nature of dark matter itself. Given a dispersion-supported equilibrium configuration in a gravitational well, many formalisms exist to determine the underlying density profile which sources the gravitational potential. Below we will consider the Jeans equations which work with the moments of the velocities. However there are also orbit-based methods (Schwarzschild 1979;Rix et al. 1997) which incorporate the fact that the objects sourcing the potential must trace orbital solutions given the potential.
With proper motion data it will be possible to distinguish members of galaxies from contaminants at high radii, where foreground stars and background galaxies dominate. There have been numerous reports of stellar excesses beyond the putative tidal radius (Komiyama et al. 2007;Vivas & Mateo 2013), sometimes with expected nonequilbrium features (Roderick et al. 2015;Battaglia et al. 2012) but sometimes without (Hayashi et al. 2003). If it can be shown that they are in equilibrium, these distant stars can be used to probe a dark matter subhalo at large radii, providing a unique test of the 1/r 3 density profile common to both cold and warm dark matter. One may also test to see whether there is a gravitationally self-consistent solution to the subhalo profile given the tidal interactions with the host which may be inferred from the satellite galaxy proper motions. If either of these tests is negative, it would provide evidence for a long range, nongravitational force.
The simplest way to test the equilibrium of a spherically symmetric system is to measure the odd moments of the stellar velocities. In equilibrium these must vanish. This provides a necessary but not sufficient condition for equilibrium. In this paper we introduce a novel, complimentary necessary condition for equilibrium. We will use the second and fourth order spherical Jeans equations to determine the time derivatives of certain odd moments of stellar velocities in terms of the even moments. A nonvanishing time derivative of any moment implies that the system is not in equilibrium, although the converse is not necessarily true. We will make no assumptions regarding either equilibrium or the isotropy of the velocities, nor regarding the form of the gravitational potential, however we will make the rather poor assumption of spherical symmetry. We will find that, even in this ideal setting, the second order Jeans equation in no way restricts the time derivatives of moments as a result of a degeneracy between the gravitational potential and the time derivative of the mean radial motion. However the fourth order Jeans equation fixes two distinct combinations of time derivatives of velocity moments, and explicit formulae for these moments will be given below.
In Sec. 2 we will use the second and fourth order Jeans equations to derive the formulas for the time derivatives as functions of velocity moments in spherical coordinates. Then in Sec. 3 we will provide an explicit map which yields the spherical coordinate moments given the observed projected 3-dimensional velocity moments. Next in Sec. 4 we will use the N -body code GALIC (Yurin and Springel 2014) to generate an example which both illustrates the results of Secs. 2 and 3 and also provides an estimate of the precision with which these time derivates may be determined. In Sec. 5 we describe future observations which may allow such an analysis. Concluding remarks appear in Sec. 6.

JEANS ANALYSIS
Let f (x, v) be the phase-space distribution function of stars, in other words the phase-space density at a position x and velocity v. By Liouville's theorem (Boltzmann 1868;Gibbs 1884), in a collisionless system f (x, v) is constant along stellar trajectories df dt = 0. (1) One obtains Boltzmann's equation (Boltzmann 1896) by expanding the total differential in terms of partial differentials via the chain rule where Φ(r) is the gravitational potential. The gravitational potential and the distribution f (x, v) are both assumed to be spherically symmetric in the spherical coordinate system (r, θ, φ) centered at some point in space. The first term on the right hand side of Eq. (2) vanishes for a system in equilibrium. Our analysis differs from many in the literature in that we do not assume that this term is zero, it is in no way constrained in our analysis. Following Merrifield and Kent (1990) but keeping this extra term, we obtain the second order Jeans equation (Jeans 1915) by multiplying by v r and integrating over the full 3-dimensional space of velocities where the stellar (position) density is The first term in Eq. (3) vanishes for a system in equilibrium. It characterizes the radial acceleration of the system. Note that observations can only determine v 2 r and v 2 θ , not Φ, and so observationally it is not possible to separate the nonequilibrium and gravitational terms in (3), at best one may determine their sum As a result, it is not possible to test the equilibrium hypothesis using the second order spherical Jeans equation alone. Similarly, the two fourth order spherical Jeans equations are obtained by multiplying the Boltzmann equation by v 3 r and v r v 2 θ respectively and integrating over velocities (Merrifield and Kent 1990) and which are consequences of spherical symmetry. Again each equation contains both a nonequilibrium term and a gravitational potential term, and so once the velocity moments have been measured at a fixed time the fourth order equations consist of three unknowns but only two equations. Including the second order equation, there is an additional unknown, the time derivative of v r , and so three equations and four unknowns. However, it is possible to determine linear combinations of the unknown time derivatives, which vanish in an equilibrium system and so can be used to test the equilibrium hypothesis.
One simple way to combine the time derivatives into an expression that can be determined observationally is simply to use the second order Jeans equation to eliminate dΦ/dr in the fourth order equations. The fourth order equations then become and Now the left hand side of each equation characterizes the departure from equilibrium, and is equal to zero for an equilibrium configuration. The right hand side consists of moments of instantaneous velocities, and so can in principle be extracted from observations of the 3-dimensional velocities, if projection effects can be removed. In simulations (see Sec. 4) we have found that the imprecision in this method is dominated by statistical fluctuations in the velocities. These statistical fluctuations can be reduced if instead of v θ one uses the tangential velocity v T , defined to be the square root sum in quadrature of v θ and v φ . The reason for this is that, due to spherical symmetry, v θ and v φ are distributed identically. As a result the total uncertainty on the sum v 2 √ 2 times the uncertainty on v 2 θ or v 2 φ . Therefore the expected fractional fluctuations expected for the quantity v 2 T are suppressed with respect to those of v 2 θ or v 2 φ by a factor of √ 2. A similar argument shows that the fluctuations of the fourth moments are suppressed even more strongly when using v T .
In this case the above equations can be written and 3. DEPROJECTION Due to the projection of stars onto the celestial sphere, one cannot determine the spherical coordinates (r, θ, φ) at a point, rather one measures the cylindrical coordinates (Z, R, Ψ) where Z is the line of sight direction, R is the orthogonal to the line of sight and Ψ is the azimuthal angle. In practice Z cannot be constrained with a precision which would be useful for our analysis. Velocities in this cylindrical coordinate system are related to those in spherical coordinates via Integrating over the line of sight one recovers the standard formulas for the observed second velocity moments where η is the observed two dimensional density. The 3-dimensional density ν at r can be obtained from the observed 2-dimensional density η via an inverse Abel transform Similarly an inverse Abel transform of the last equation in (14) yields the second moment Transforming the sum of the first two equations in (14) one finds Subtracting (16) one finds the last second moment Summarizing so far, Eqs. (15), (16) and (18) provide the 3-dimensional density and second moments, the expressions appearing in the second order Jeans equation (3), from the observed projected density and second moments. The above equations express three unknowns v 2 r , v 2 θ and ν in terms of four observables v 2 R , v 2 Z , v 2 Ψ and η using only three equations. The fact that there is one more observable than equations implies that this system contains one exact degeneracy. More precisely, there is a one-parameter family of possible observations for which these equations yield the same set of unknowns. As a consequence it would be possible to modify the set of three equations by shifting the observables along the degenerate direction without affecting the obtained estimates of the unknowns. These equations are therefore not unique. Such a shift of the equations along the degenerate direction would not affect their validity or physical content. However the precisions with which the observables are measured do depend on the shift, and so such a shift would affect the resulting uncertainties in the estimates of the unknowns. In summary, while the equations above are correct, they do not offer optimal precision in the determination of the unknowns. This could be achieved, for example, with a χ 2 fit of the observables to the line of sight projections of the unknowns.
In summary, the above set of equations is not unique and does not minimize the uncertainties achieved on the parameters. We will nonetheless use them because they provide a fast estimation of the deprojected parameters, and computational speed will be important given the large number of simulations which we will analyze. However, the deprojection of real data can, at the expense of greater computation time, be performed more precisely using standard methods such as a χ 2 fit.
Similarly, the projected fourth moments are integrals along the line of sight of the spherical coordinate fourth moments As in the case of second moments, one may obtain a combination in a form which is easily Abel-transformed by noting and using the three combinations of moments which are products of these two invariants Linear combinations of the inverse Abel transforms then yield the spherical coordinate moments in terms of the observed projected moments The above equations express three unknowns v 4 and η using only three equations. As there are more observables than equations, again the obtained unknowns are invariant under transformations of the observables. Again this means that the equations may be transformed into a form which reduces the resulting uncertainties, and so they are suboptimal.
In summary, one first observes the projected, r-dependent 2-dimensional density η, and the second and fourth order moments. These can be transformed into spherical coordinates via the inverse Abel transforms (15), (16), (18) and (22). Finally, these spherical-coordinate fixed-time moments can be substituted into Eqs. (11) and (12) to obtain an estimate of the time-derivatives of the odd order moments. These time derivatives are nonzero only if the system is not in equilibrium. Thus, given spherical symmetry and sufficiently many proper motions 1 , one can test the equilibrium of the system without knowledge of the potential or any assumptions regarding the isotropy of the orbits. This is our main result.

The Setup
In this section we will present a simple example of Eqs. (11) and (12) at work. We will use the N -body code GALIC (Yurin and Springel 2014) to produce systems of 10 5 particles of mass 1.86 × 10 4 M which interact gravitationally. As the only interaction is gravitational, at this point each particle could be either dark or stellar matter. We will consider both spherical and axisymmetric mass distributions, so that we may estimate the systematic errors created in our approach as a result of our false assumption of spherical symmetry.
These correspond to Hernquist dark matter halos (Hernquist 1990), with no other components, described by the GALIC standard models H1 and two variations of the model H6, each rescaled by a factor of ten to very roughly yield typical dwarf spheroidal galaxy dimensions. In particular, the Virial radius r 200 of each halo is 20 kpc and the Hernquist parameter is a = 3.45 kpc.
For now we will assume that the true radial coordinates of the particles are known. In reality, this requires parallax measurements with a precision unlikely to be obtained for stars in dwarf spheroidal galaxies in the foreseeable future. As a result, we only test Sec. 2, not Sec. 3. In Subsec. 4.7 we will use the projected observables together with the deprojection in Sec. 3 to provide a comprehensive test of our results.
These made-to-measure simulations fix the density distribution of the particles and adjust the velocities towards an equilibrium configuration one at time, beginning with Gaussian density profiles that approximate those of Hernquist profiles. Unfortunately such initial configurations are already quite close to equilibrium, and so this approach will not allow us to test our formalism far from equilibrium 2 . Each time that a given particle is adjusted, the system is driven closer to, or at least not further from, equilibrium. Thus by measuring the system after different numbers of iterations per particle, one can obtain velocity profiles for various distances from equilibrium.
The arguments in Sec. 2 imply that, as the systems approach equilibrium, the right hand sides of Eqs. (11) and (12) should approach zero for a spherically symmetric density profile, whereas one expects that a breakdown of spherical symmetry will lead to a nonzero limiting value. This is because spherical symmetry was used in the derivation of Eqs. (3,6,7). In the absence of spherical symmetry, in general the left hand side of each equation will be an arbitrary function of the position. For a spherically symmetric equilibrium system, the time derivatives in these equations vanish and the right hand sides of Eqs. (11) and (12) are linear combinations of the left hand sides of Eqs. (3,6,7), and so are also arbitrary functions of position.
To test this prediction, we let α i and β i denote the four terms of the right hand sides of Eqs. (11) and (12) respectively, so that these equations imply that for an equilibrium configuration We divide the particles into bins corresponding to spherical shells of thickness 200 pc each, and for the jth bin we calculate α ij and β ij . Here i labels the four terms in Eqs. (11) and (12) and j labels the bin, so that α ij and β ij are just given by α i and β i evaluated using the results of GALIC for the jth bin. The parameters α ij and β ij are proportional to the particle density ν and so are smaller at high radii, where the particle number falls. Here the derivatives in Eqs. (11) and (12) are evaluated using the discrete difference between neighboring radial bins.
To obtain a measure of the deviation from equilibrium which is independent of such normalizations, we define the dimensionless index We will often omit the index j. For a generic configuration, one would expect all of the α ij and also all of the β ij to be uncorrelated. In this case the cross terms in the numerator of X would be of the same order as the squared terms and so X would be the sum of two terms each of order unity. Therefore one would expect X itself to be of order unity. In Subsec. 4.6 we will see that this is indeed the case for a system far from equilibrium. On the other hand, a lower value of X implies that the terms on the right hand sides of Eqs. (11) and (12) nearly cancel each other. If Sec. 2 is correct, then X should approach 0 as a system tends to equilibrium (where Eq. (23) holds) and so as the number of steps in the spherically symmetric simulations increases. On the other hand, for an aspherical equilibrium system we have argued that the right hand sides of Eqs. (11) and (12) are arbitrary functions of position. These right hand sides are integrated into the radial bins i α ij and i β ij , and so these will be arbitrary functions of the bin index j. As a result, X j will generically be nonzero even for equilibrium aspherical configurations. Therefore, as an aspherical system approaches equilibrium one expects X to tend to a nonzero function corresponding to the systematic error caused by the unwarranted assumption of spherical symmetry.

Results
Each simulation produces a near equilibrium configuration of 10 5 gravitationally interacting particles. These particles can be stellar or dark matter. To compare with observations, one needs to decide which particles are dark matter and which are stellar. This must be done such that each population is separately near equilibrium. We will define this segregation by simply imposing that each particle with specific energy less than −1000 (km/sec) 2 is a stellar tracer. As specific energy is a conserved quantity, such subpopulations of an equilibrium system will each be in equilibrium given the total gravitational potential. With this tracer condition, each simulation produces about 1.8 × 10 4 stars. As will be explained in Sec. 5, this corresponds roughly to the number of stars per classical dSph whose proper motions will be measured by upcoming extremely large telescopes. The radial distribution of the tracers is also deep inside of the dark matter halo, as is expected for a typical dSph.
The tracer condition is of course not realistic, most obviously because the stars are far too massive. However it is computationally cheap and we believe that this crude approximation does not affect our conclusions regarding the effectiveness of the parameter X as a test of equilibrium.
In Fig. 1 we plot our results for the value of X in all 200 pc bins up to 4.5 kpc, beyond which the tracer condition cannot be satisfied. In the inner 800 pc and also beyond 3 kpc, the number of particles per bins is small, leading to considerable shot noise which is compounded by the large bin size. The large bin size is an important source of error because the derivatives are calculated using differences in neighbouring bins. The top, middle and bottom panels represent a spherically symmetric, isotropic density distribution and two prolate axisymmetric distributions with eccentricities of 0.85 and 0.8. The number of tracers per bin in each case is plotted in the panels on the right. The lowest values of X are achieved at radii of 1.5 to 3 kpc, where each bin contains a large number of particles and also the derivatives are relatively small and so well approximated by the differences between neighboring bins. The error caused by finite binning is estimated in Subsec. 4.6.
The error bars represent only statistical fluctuations, not the finite binning uncertainties. The statistical fluctuations are estimated by partitioning the tracers and calculating X for each partition. Then the uncertainty is estimated as the square root of the variance of the X values of these subsamples divided by the square root of the number of partitions. We have checked that this value is independent of the number of partitions so long as there are few enough partitions that the fluctuations are smaller than the mean value. As the finite binning uncertainties are not included, we expect these errors to be grossly underestimated at small radii and also near 4 kpc.
The main result of this section is the fact that X indeed is quite close to 0, as compared with the expected nonequilibrium value of order unity, when these errors are small. In general it is of order 10 −2 in this case of roughly 10 4 tracers. Interestingly, this remains true even in the strongly prolate configurations. Thus the breakdown of spherical symmetry does not appear to significantly compromise the ability of X identify equilibrium configurations.
There is little evolution in X between the 0th and 500th steps of the relaxation, which indicates that this tracer population already begins quite close to equilibrium. We will see in Subsec. 4.5 that this is not true for all choices of tracer. This is because higher energy tracers, which were removed by the energy cut in this subsection, more often sample the nongaussian part of the Hernquist distribution, whereas GALIC configurations begin with purely Gaussian velocity distributions. Thus the higher energy particles in GALIC begin further from equilibrium. In a real galaxy one does not know a priori which subpopulation best approximates equilibrium, however one may use the methods of this paper to measure X for each. The subpopulation with the lowest value of X may be the closest to equilibrium and so the most reliable for the construction of a mass model.
The individual terms in X are displayed in Fig. 2 in the case of the spherically symmetric halo. We can see that while individual terms increase or decrease their relative significance as a function of radius, nonetheless the equilibrium condition always conspires to make them nearly cancel.

Aspherical Halos
The formulae presented in this paper were derived assuming a spherically symmetric halo, and so one only expects X to approach zero in this case. The purpose of the simulation of the prolate halos is to determine whether the resulting increase in X is sufficiently dramatic to render X insensitive to the distance between the dark matter configuration and equilibrium. On the contrary, in Fig. 1 one can observe that the eccentricity has only a modest effect on X, with the highly prolate bottom panel displaying an X which is larger than that in the spherical top panel by only about one standard deviation in each bin. As a result we are unable to quantify the effect of asphericity on the X value of the tracer population at this time. In Subsec. 4.5 we will present an example in which the eccentricity has a larger effect on X.
What about realistic dark matter halos? The GALIC code can only generate ellipsoidal halos. However in ΛCDM simulations halos are triaxial (Frenk et al. 1988) with radius-dependent axis ratios (Hayashi et al. 2007). The triaxiality and radial-dependence may provide serious complications to certain quantitative tasks, such as determining the shape of the dark matter halo of a galaxy cluster from the shapes of X-ray isophotes (Hayashi et al. 2007). However, we expect that the correlation between X and the equilibrium condition increases monotonically with the axis ratios. In other words, the closer the axis ratios are to unity, the more accurately X reflects the distance from equilibrium. Therefore, for the purpose of estimating the impact of asphericity on X in a certain radial bin, it is sufficient to consider a typical axis ratio that occurs in the line of sight integral corresponding to that bin.
The Aquarius simulations (Springel et al. 2008) of the evolution of pure dark matter halos were used to study subhalo shapes in Vera-Ciro et al. (2014). The authors then used the semianalytical method of Starkenburg et al. (2013) to determine which subhalos were likely to host local group satellite galaxies. The authors found that these had minor to major axis ratios ranging from 0.4 at about 1 kpc to 0.8 at the radius which the authors called r 95 , equal to about 2/3 of the Virial radius. In our case, at the smallest radii, X is dominated by binning effects and in particular the imprecision in our discrete derivative. At the radius r 95 the authors also estimated the intermediate length principal axis to major axis ratio and found about 0.9, suggesting that axisymmetry is a reasonable approximation. Tracers Figure 1. The panels represent a spherically symmetric (top) and two prolate axisymmetric configurations, with axis ratios of 0.85 (middle) and 0.8 (bottom). The two curves correspond to running GALIC with 0 (black) and 500 (red) steps/particle. The right panels count particles with specific energies beneath −1000 (km/sec) 2 . The left panels plot X for each 200 pc radial bin. The error bars reflect only statistical fluctuations. A value X ∼ 1 would correspond to a generic nonequilibrium configuration, the small values of X imply that the tracers already begin near equilibrium, and stay near equilibrium.

Time scales and dimensional analysis
As these systems are not in equilibrium, the velocity moments evolve with time. In this subsection we will use a dimensional analysis of the Jeans Eqs. (11) and (12) to relate the quantity X to the time scale t over which the velocity moments evolve.
Let us make the crude approximation that each term in Eq. (24) is equal to X/2. Now each term on the right hand sides of Eqs. (11) and (12) is of order v 4 /r where v is a characteristic stellar velocity, say 5 km/s for a UFD or 10 km/s for a classical dSph, and r is the radial coordinate, which is of order 100-300 pc for known dSphs. Thus the denominator of each term in Eq. (24) consists of four terms, each of which is about v 8 /r 2 . For simplicity we add these four terms in quadrature in what follows, as if they were of comparable magnitudes. On the other hand, the total right hand side is equal to the left hand side, which is of order v 3 /t where t is the characteristic timescale over which the nonequilibrium system evolves. For example, if the system were in equilibrium then t would be infinite.
One can then estimate the order of each term in X as where p is the period of a typical star, which we have roughly estimated to be 4r/v. We can now interpret X as giving the timescale over which the nonequilibrium system significantly evolves which for X ∼ 10 −3 would be of order one Gyr. This means that over a period of Gyr, if is sufficiently small then one would expect a fractional change in the odd velocity moments of order . As a result, if X is far from 0 but the first and third velocity moments are still quite small, this suggests that the system had been in equilibrium until it was disrupted roughly within the past Gyr. It is tempting to attempt to compare the disruption time estimated here with the relaxation time of a system relaxing in GALIC. However we know of no direct relation between the relaxation steps of the GALIC code and the real time evolution of a relaxing system, indeed it is unlikely that any such relation exists as the relaxation rate in GALIC is optimized to run quickly, rather than matched to the dynamical relaxation rate.

Application to All Particles
In the examples thus far we have considered small deviations from equilibrium. In these cases we have observed that the value of X is independent of the number of steps, and so appears to be independent of the distance from equilibrium. In Subsec. 4.6 we will see that this is not true for larger deviations from equilibrium. In this subsection, to show that X can indeed depend on the distance from equilibrium even for small perturbations from equilibrium, we will provide an example in which X does evolve with the number of steps. This example is the same as that presented in the previous section, except that instead of restricting our attention to tracers, we consider the full 10 5 particles. As these higher energy particles have a much greater spatial extent, we will use 4 kpc radial bins instead of 200 pc bins.
Our results are shown in Fig. 3. Here one can clearly observe that X decreases as the number of steps increases, and so as the system relaxes towards equilibrium. In the spherically symmetric case, X reduces from about 2 × 10 −2 to about 10 −3 at intermediate radii, where the errors are smallest. This improvement is less evident in the more prolate cases, suggesting that X saturates at some nonzero value, apparently between 2 × 10 −3 and 5 × 10 −3 at intermediate radii. Such a saturation to some extent can be attributed to the systematic error induced by falsely assuming spherical symmetry, but it may also result from the fact that the GALIC code relaxes prolate configurations more slowly than spherically symmetric configurations, although the later possibility seems unlikely as we have observed little GALIC evolution between steps 100 and 500. Particles Figure 3. As in Fig. 1 except that all 10 5 particles are considered, not only the tracers. Also the bin size has been increased to 4 kpc. The number of particles per bin is independent of the number of the steps as GALIC uses a made-to-measure algorithm to generate the distributions, and so the spatial coordinates of the particles are fixed.

Tracers Drawn From Known Distribution Functions
In this subsection we do not use the N -body code. Instead our tracers sample phase space distribution functions which we have fixed by hand. This will allow several tests which would not be feasible with the full code.
We apply the same analysis as above to tracer populations drawn from four well-defined phase space distribution functions. The first is a spherical, isotropic Hernquist profile with the same parameters as that used in the GALIC simulations. This is our only example of a true equilibrium system. The second and third are the same as the first but with one of component of the velocity at radial coordinate r multiplied by 1 + r/(2a) and 1 + 2a/r respectively, with a equal to 3.45 kpc. In the last, each velocity component is multiplied by 1 + 2sin (20r/a). This is our only example of a system which is far from equilibrium. In each case 10 5 particles were used with the same tracer condition, E < −1000 (km/sec) 2 , as in the simulated samples.
The results are shown in Fig. 4, where again the bin size is 200 pc. The first profile is an equilibrium configuration. The second configuration is not in equilibrium while the third is farther from equilibrium. The last is very far from equilibrium and so as expected X is of order unity. Thus X in the first case is lower than in the GALIC configurations, which are close to equilibrium but still not in equilibrium. One does not expect X to be zero in any case due to the finite sample size and the corresponding finite binning. As expected, the values of X are larger in the nonequilibrium second configuration and even larger in the third, which is yet further from equilibrium.
In practice the establishment of equilibrium requires a robust determination of the uncertainties entering X. The first curve in Fig. 4 shows the contribution to these errors arising from statistical fluctuations together with finite binning, but does not include asphericity.  Figure 4. This plot was generated by directly sampling 10 5 particles from predefined spherically symmetric distribution functions, not using an N -body simulation. In each case, the tracers are particles with specific energy less than −1000 (km/sec) 2 . The four curves correspond to an isotropic Hernquist distribution (black), the same distribution with one component of the velocity multiplied by 1 + r/(2a) (blue) and 1 + 2r/a (red) and also with all velocities multiplied by 1 + 2sin (20r/a) (green). The first curve samples an exact, spherical equilibrium distribution. Therefore the nonzero value of X is entirely due to statistical fluctuations and finite binning. The latter curves are further from equilibrium, and correspondingly X is larger. 10 -1 X Figure 5. As in Fig. 4 but only the Hernquist distribution is sampled. The black curve uses tracers from a 10 5 particle sample with 200 pc bins, the red and blue curves draw tracers from a 10 9 particle sample with 200 pc and 20 pc bins respectively. The red curve provides an estimate of the error due to finite binning alone.
As these configurations are simply drawn from known phase space distributions, and not generated by N -body simulations, it is not computationally expensive to create configurations with as many as 10 9 particles. In Fig. 5 we compare the 10 5 particle Hernquist configuration with two 10 9 particle configurations, with bin sizes of 200 pc and 20 pc. In each case, the same tracer condition is applied which selects 18% of the particles.
The 10 9 particle calculations, as expected, have much lower values of X. This of course is expected due to the small statistical fluctuations. The 200 pc bin case in fact has negligible statistical fluctuations and so it reflects the errors that are expected from binning alone, not only in the present case but throughout the paper where 200 pc bins are used. In fact, due to the very high number of tracers per bin, this 200 pc bin curve is the only curve in the paper whose statistical fluctuations are so small that they are not apparent on the plot. As expected, the binning error becomes greatest near the origin and also near the tracer edge at r = 4.5 kpc, where the bin to bin fractional changes are large and so the finite bin differences are a poor approximation for the derivatives.
On the other hand, the 20 pc bin case is clearly dominated by statistical fluctuations. While all X plots in the paper in 200 pc bin cases turn upwards beneath 1 kpc, this 20 pc bin plot turns upwards only beneath 100 pc. This confirms the fact that the high value of X beneath 1 kpc in the other plots is due to a combination of finite binning and shot noise. Here, with both of these contributions well below 10 −3 , one sees that the intrinsic value of X indeed is below 10 −3 . This is consistent with our central claim that in an equilibrium, spherical configuration, in the limit of a large number of tracers and perfect measurements, X tends to zero.

Projection
So far we have assumed that the full 3-dimensional positions and velocities are available to the observer. In practice only the projection along the line of sight is available. While the Thirty Meter Telescope will be able to measure the full 3-dimensional velocities for as many as 10 5 stars in the brightest dSphs, the positions will nonetheless be projected onto the celestial sphere.
We have also considered the projections of the GALIC configurations described above, with the same specific energy condition at −1000(km/sec) 2 . GALIC generates configurations with a rotational symmetry about the z-axis. For concreteness, we have considered observations with the line of sight along the x-axis. We then used the results of Sec. 3 to deproject these projected observables, yielding 3-dimensional observables which we used to calculate X as above. In the deprojection one needs to calculate derivatives of the projected moments. These derivatives, like all derivatives of moments in this note, are calculated by averaging the discrete difference from neighboring radial bins with the discrete difference from next-to-neighbor radial bins. An improved scheme for evaluating the derivatives would reduce the resulting fluctuations and as a result reduce X. To improve the precision of the derivatives is quite easy, one need only average information from nth neighbor bins where n > 2. However this improved precision will lead to reduced accuracy if the functions in question are not sufficiently linear on the distance scale equal to n times the bin size. Thus a higher precision comes at a cost of less sensitivity to the fine structure of the mass profile.
In Figs. 6 and 7 we present the equilibrium indicator X obtained using only the line of sight projected moments, restricting and not restricting our attention to the tracer subpopulation respectively. The main result of our analysis is that in general X, within the reported errors, is similar to or somewhat smaller using the projected moments. The cause of this somewhat surprising result is that the inverse Abel transforms, which are used to deproject the data, enforce spherical symmetry, even when the system is not spherically symmetric. This enforcement has little effect on the second moments. However the GALIC code begins by setting the second moments to their equilibrium values, it is the higher moments which evolve towards equilibrium in GALIC. And the deprojection has an averaging effect on the fourth moments which in these cases drive them closer to an equilibrium configuration, even in the case of a spherically symmetric configuration. This is because, when determining the moments from Eq. (22), one imposes all of the relations between the fourth moments which would be implied by spherical symmetry, such as Eq. (8). In our numerical calculations we have found that enforcing these relations greatly reduces the shot noise in X. Intuitively this is quite reasonable: statistical fluctuations generally do not respect the spherical symmetry, and so by spherically averaging these fluctuations are reduced.
How is it possible that with less information we obtained a better value of X? The reason is that X has simply become less sensitive to the details of the configuration, for example obviously X is insensitive to a distortion along the line of sight direction. Equilibrium implies X = 0, but the converse is not true and in the projected cases the equilibrium test provided by X is somewhat less powerful.

APPLICATIONS TO FUTURE OBSERVATIONS
To our knowledge, today there exists no astrophysical dataset to which these formulae can be usefully applied. The problem is that the nearest collisionless, dispersion-supported systems, the Milky Way's satellite galaxies, are too far away for the proper motions of individual stars to have been determined with a precision of order the velocity dispersion. However this will soon change.
In the next two or three years, the Gaia Space Telescope will determine the proper motions of more than 100 members of, for example, the Sculptor dwarf galaxy (Evslin 2015). The precision, better than 20 km/sec will be comparable to, although somewhat greater than, the stellar dispersion, and so will already provide interesting hints of a core or cusp, especially when combining the results for multiple satellite galaxies. However, in the present note we considered simulations with more than 10 4 tracers and even so we saw that the determination of equilibrium is challenging. Therefore with 10 2 tracers Gaia will only be sensitive to extremely large departures from equilibrium. However Gaia will not provide only proper motions, more importantly it has provided and is continuing to provide precise astrometric positions for hundreds of members of the Sculptor dwarf and other dSphs. These positions may be combined with a single epoch of precise astrometry using extremely large telescopes 10 years from now. The long baseline between these observations will yield 5 km/sec precision proper motions for the Gaia sample, with 10 km/sec for hundreds of other Sculptor members which were also observed by Gaia. At least in the spherically-symmetric case, this is already sufficient for a robust distinction of a core or cusp in a classical dSph if equilibrium can be shown to sufficient precision.
With a second epoch of observations, 5-7 years later, the Thirty Meter Telescope (TMT) will measure, with 3-5 km/sec precision, proper motions of 10 4 or more stars in the largest classical dSphs. As site selection for TMT is again open, with candidate sites in the northern and southern hemispheres, it is not yet clear which targets will be chosen. However the results of this note suggest that one needs of order 10 4 stars. To obtain the proper motion of a dSph member at 80 kpc with an uncertainty beneath the stellar dispersion, given less than a minute of observations per field per epoch, one requires an H magnitude of at most 22 (Skidmore et al. 2015). Therefore only the classical dSphs are suitable targets for such a study. Furthermore the Sculptor dwarf is quite large in the sky, and so would require several weeks of observing time per epoch. While Fornax's greater distance implies greater uncertainties per star with a fixed viewing time, given the greater stellar population it remains a promising target if a southern hemisphere site is chosen. The E-ELT is also capable of such measurements, but the systematics-limited astrometric precision will be 50 µas, as compared with 10-20 µas for TMT's first light IRIS instrument. This means that E-ELT will in general measure proper motions with errors comparable to the intrinsic velocity dispersions. This is critical, because if the measurement errors are smaller than the intrinsic velocity dispersions, then they will provide subdominant contributions to the total error budget in the determination of quantities such as the halo shape and also the equilibrium parameter X. Therefore we expect the above analysis of the uncertainties expected in X to be applicable to these observations. How might the results of this paper be applied to these datasets in practice? Our main result is that the matrices α ij and β ij , which may be summarised into X j , characterize the departure from equilibrium. More precisely, for a precisely spherical system with infinitely many observed tracers, X j = 0 in each radial bin j. In practice, X j will be nonzero even for an equilibrium system as a result of shot noise, asphericity and measurement uncertainties. In this paper we have studied the effects of shot noise and asphericity on X, while we believe that measurement uncertainties are unlikely to dominate over shot noise so long as these uncertainties are below the stellar dispersion, as is expected for the TMT observations described above. Our studies of these sources of error used the GALIC code, which produces dark matter halos which are close to Hernquist profiles and a rather simple tracer condition. To determine the uncertainties which are really expected in X, one needs to use a more realistic mass model for the dark matter and tracers, which of course needs to be determined by the data itself.
One way to implement this requirement is as follows. First, one may assume equilibrium and analyse the data using a Jeans analysis or orbit based methods to derive a 3-dimensional mass distribution for the dark matter and the stars in a dSph. Next, one can repeat an analysis like that done in this note to calculate the range of X j expected in an equilibrium configuration with such a distribution. The observed X j can be compared with this range to determine the confidence with which the equilibrium hypothesis may be excluded. This strategy may be applied to all of the observed stars, or to individual subpopulations.

CONCLUSIONS
We have provided a proof-in-principle that, for a spherically symmetric system, one may test equilibrium using the fourth order Jeans equations. We have introduced a test statistic X which tests the equilibrium hypothesis, such that for a spherically symmetric, equilibrium system X tends to zero as the number of particles observed tends to infinity. This property of X requires no assumptions on the form of the radial-dependence of the gravitational potential, nor on the velocity anisotropy.
We have calculated X, bin by bin, in spherically symmetric and axisymmetric configurations. We have considered both near-equilibrium configurations generated by GALIC and also samples drawn from exact equilibrium configurations, and their deformations. We have analyzed these directly from the full 6-dimensional phase space and also using only the projected observables together with inverse Abel transforms. In each case, we considered a sample of 10 5 particles of which 1.8 × 10 4 are tracers, with a radial and velocity distribution similar to stars in a very large dSph 3 . Due to our discrete derivatives and other binning effects, our uncertainties were large in the inner 800 pc. However, with this sample we were able to show that the indicator X works quite well, with a precision and stability below 10 −2 in the next 1 or 2 kpc. At larger distances, the small number of stars begins to lead to large statistical fluctuations.
Although less information is available, we found that the same configurations lead to similar or smaller values of X when only the projected data is used. However in this case X provides a less powerful test of the equilibrium hypothesis, as there are more nonequilibrium perturbations to which it is insensitive.
To actually conclude the significance with which equilibrium can be excluded, one may for example follow the procedure described at the end of Sec. 5. This can be done rather straightforwardly. However one must recall that this procedure is suboptimal, in the sense that it does not deliver the smallest possible uncertainties. This is because, as was described in Sec. 3, our formulas for the deprojection did not use the full available data in the observed moments. They were chosen because they are analytic formulas which are easily computed. However a more precise deprojection could be obtained by selecting the 3-dimensional distribution whose projection best matches the observed projected moments, as characterized for example by a χ 2 fit. More specifically, one may parametrize the various sphericalcoordinate functions and use (14) and (19) to determine the corresponding projected moments. Then a χ 2 , statistic can be used to find the best fit moments and their associated uncertainties. These can be inserted into Eqs. (11) and (12) to test the equilibrium of the system, with the uncertainties in the fixed time spherical coordinate moments propagated in the usual way through these two equations.
A quick and dirty approach to testing the equilibrium hypothesis is as follows. First, note that to determine the spherical second moments it suffices to use two of the three equations in (14) and three of the six in (19). Therefore there are 60 distinct ways to determine the moments, which using (11) and (12) yield 60 distinct estimates for the time derivatives. One can then estimate the uncertainty in X j to simply be the square root of the variance of the 60 values of X j given by these 60 permutations of our procedure. This source of error should be added in quadrature to other sources arising from shot noise and finite binning. Then one may test the null hypothesis of equilibrium by simply calculating the number of standard deviations between X and 0. This yields a very rough estimate of the significance with which the equilibrium hypothesis may be rejected.
To exclude equilibrium one needs a systematic investigation of the bias caused by the aspherical halo shape, for example via an analysis similar to that done in Laporte et al. (2013) or Pontzen et al. (2015). The main shortcoming with our approach in its current form is the assumption of spherical symmetry. The goal of our research program is to eventually generalize these results to aspherical systems. While we have seen that our spherical results may be applied to an aspherical system without an enormous increase in the uncertainties, we believe that a more precise incorporation of the asphericity will strengthen our test of equilibrium. To do this one will need to consider not only the radial dependence of the projected moments, but also the azimuthal dependence of these moments, as in Evslin (2016).