Entropy Production in Collisionless Systems. III. Results from Simulations
aa r X i v : . [ a s t r o - ph . GA ] M a r Draft version August 4, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
ENTROPY PRODUCTION IN COLLISIONLESS SYSTEMS. III. RESULTS FROM SIMULATIONS
Eric I. Barnes
Department of Physics, University of Wisconsin — La Crosse, La Crosse, WI 54601
Colin P. Egerer
Department of Physics, University of Wisconsin — La Crosse, La Crosse, WI 54601
Draft version August 4, 2018
ABSTRACTThe equilibria formed by the self-gravitating, collisionless collapse of simple initial conditions havebeen investigated for decades. We present the results of our attempts to describe the equilibria formedin N -body simulations using thermodynamically-motivated models. Previous work has suggested thatit is possible to define distribution functions for such systems that describe maximum entropy states.These distribution functions are used to create radial density and velocity distributions for comparisonto those from simulations. A wide variety of N -body code conditions are used to reduce the chancethat results are biased by numerical issues. We find that a subset of initial conditions studied lead toequilibria that can be accurately described by these models, and that direct calculation of the entropyshows maximum values being achieved. Subject headings: galaxies:structure — galaxies:kinematics and dynamics INTRODUCTIONOver the past several decades, numerous investigationsof collisionless, self-gravitating systems have been under-taken. From early focus on the formation and evolutionof elliptical galaxies (e.g., van Albada 1982), to cosmo-logical simulations of dark matter structure formation(e.g., Navarro, Frenk, & White 1996; Moore et al. 1999;Springel et al. 2005), these works have taken advantageof ever-increasing levels of computing power. Simula-tions with higher and higher resolutions (numbers of par-ticles per unit volume) are constantly being performed.Our goal is not to attempt to replicate these state-of-the-art simulations, but rather to use more modest sim-ulations to investigate some basic questions. The sys-tems that we simulate are not direct analogues of puta-tive dark matter halos nor elliptical galaxies, but theydo share the fundamental physical conditions of self-gravitation and collisionless evolution. As many previ-ous works have shown evidence of “universal” behaviors[such as radial density (e.g., Navarro, Frenk, & White1997; Navarro et al. 2004) and power-law “phase space” ρ ( r ) /σ ( r ) (Taylor & Navarro 2001) profiles], an inter-esting possibility is that a basic physical mechanism mayunderlie the formation of these self-gravitating equilibria.Earlier work (Barnes & Williams 2011, 2012) presentsdescriptions of distribution functions of collisionless,self-gravitating systems that represent maximum en-tropy states. These results follow the seminal work ofLynden-Bell (1967), where it is argued that a fourth sta-tistical family is appropriate for describing the phase-space evolution of these kinds of systems. In a nut-shell, the familiar Maxwell-Boltzmann statistics describesystems in which particles are distinguishable and donot obey a phase-space exclusion principle – multi-ple particles can occupy a very small region of phase- [email protected]@uwlax.edu space. On the other hand, the Lynden-Bell statisti-cal family is appropriate for systems of distinguishableparticles that follow an exclusion principle – a classi-cal version of Fermi-Dirac statistics. By relaxing therequirement of large phase-space occupation numbers,Barnes & Williams (2012) show that finite-mass, max-imum entropy states exist for both Lynden-Bell andMaxwell-Boltzmann statistical families. The major goalof this work is to test whether or not any simulatedsystem will relax to such a state. It is certain thatthese models will fail for sufficiently strong collapses,as the assumed velocity isotropy underlying the mod-els has to disappear as the radial orbit instability be-gins to become important ( e.g., Merritt & Aguilar 1985;Barnes, Lanzel, & Williams 2009). As such, we also aimto identify ranges of conditions that can result in maxi-mum entropy states. A final goal is to monitor the be-havior of entropy in simulations to see if a maximumvalue is reached.We have created suites of simulations to test the use-fulness of the Lynden-Bell and Maxwell-Boltzmann fam-ilies of models. The publicly-available GADGET code(Springel 2005) has been employed to evolve simula-tions of collisonless systems comprised of N = 10 and N = 10 particles. These simulations approximate colli-sionless conditions by utilizing softened interactions. Asa point of comparison, we have also analyzed collisionalsystems with N = 2 ≈ particles using a version ofNBODY-6, enhanced with a Graphics Processing Unit(GPU) (Nitadori & Aarseth 2012). This code uses directNewtonian particle-particle interactions, but the largenumber of particles guarantees that two-body relaxationprocesses occur over timescales thousands of times longer(Binney & Tremaine 1987, Ch. 4) than any gravitationalpotential (or “violent”) relaxation processes, which occurover a few initial crossing times T .For any given simulation, we analyze spherically sym-metric density and velocity distributions by counting par- Barnes & Egererticles in spherical bins centered on the system center-of-mass. With estimates of the various uncertainties, theseradial profiles are then used as the data to be matched bythe Lynden-Bell and Maxwell-Boltzmann models. Chi-squared minimizations are used to indicate the appropri-ateness of each model. We will not insist that modelswith low chi-squared values are the only descriptions ofthese simulations, merely that such a model is consistentwith the data. For simulations that are well-describedby the thermodynamic models, we also investigate thebehavior of entropy during its evolution.We begin by describing simulation initial conditions,evolution code details, and analysis techniques in sec-tion 2. Two methods for describing entropy behavior inthese simulations are presented in Section 3. Section 4contains the findings inferred from minimizations for se-lected simulations, while section 5 outlines the entropybehaviors seen in the simulations. We summarize in Sec-tion 6. SIMULATION DETAILS2.1.
Initial Conditions
For all simulations discussed here, particles are as-sumed to be identical and system mass, system radiusand Newton’s gravitational constant are set to unity( m sys = R = G = 1). Each system is composed of N particles. Following van Albada (1982), particle posi-tions are chosen according to two different schemes.In “single” simulations, initial particle locations arerandomly chosen within the system according to a spec-ified density distribution – cuspy ( ρ ∝ /r ) or Gaussian( ρ ∝ exp − r ). A simple rejection scheme is used togenerate the distributions. The system center-of-massmust coincide with the center of the spherical boundaryto better than 1 / √ N to be acceptable.For “clumpy” simulations, centers-of-mass locationsare chosen for small clumps of particles according to theabove density distributions, and then particles are uni-formly distributed within a clump. The numbers of par-ticles in each clump are chosen from a Salpeter distribu-tion – the probability of generating a clump with N clump particles is proportional to N − α clump , where α = 2 / α )have not been scrutinized. The Salpeter form has beenadopted based on its simplicity. We demand that the sumof the volumes of all the clumps equals that of a spherewith radius R = 1. Individual clump radii are chosenproportional to the fraction of the total mass they con-tain, r clump ∝ ( N clump /N ) / . Clumps can overlap oneanother, leading these initial conditions to have regionsof higher-than-average density as well as nearly emptyregions where clumps fail to overlap.Particle velocities are given random orientations, guar-anteeing initial velocity isotropy. Speeds are chosen byadopting an initial virial ratio Q = 2 K / | W | that linksa system’s initial kinetic energy K to its initial potentialenergy W . Once particle positions have been selected,the virial ratio is used to define a scale speed. For singlesimulations, this is the speed given to every particle inthe system. Clumpy simulations distribute speeds in amore complicated manner. We define hot-clumpy sys- tems to be ones in which the clump centers-of-mass havezero initial velocities – particles in clumps are given ran-dom velocity directions with equal speeds. Cold-clumpysystems are composed of clumps in which all of the par-ticles move with the clump center-of-mass velocity. Thecenters-of-mass velocities are randomly oriented but havethe same magnitudes. As an intermediate case, warm-clumpy simulations split the kinetic energy equally be-tween individual particle motion and clump center-of-mass motion. Independent of the specifics of the setup,the systems are not initially in mechanical equilibrium,even though they may be in virial equilibrium (if Q = 1is adopted). 2.2. Evolution Code Details
As mentioned in the introduction, this work utilizestwo very different codes for evolving initial conditions.Our aim is to be able to identify any numerical effectsdue to particle number, softening parameters, and/orcode specifics. The GPU-enhanced NBODY-6 code hasan architecture designed for investigating globular clus-ter dynamics. GADGET has been designed to performcosmological simulations of structure formation. Find-ing agreement between our predictions and the results ofboth types of simulations will strengthen the assertionthat our analytical picture is relevant to collisionless sys-tems and is not biased by numerical issues.GADGET is a versatile tree-code that incorporatessoftened forces. In general, we have adopted the standardparameters for GADGET. However, we do not evolve inan expanding universe, and we adopt different softeninglengths, depending on the situation. For N = 10 simu-lations, we adopt a softening length ǫ = 10 − , about 100times smaller than the “optimal” softening length valuedescribed in Power et al. (2003). Test runs with soften-ing lengths between our adopted value and the optimalvalue have resulted in only minor differences in densityprofile shape. Smaller values lead to integration timesthat we judged to be unacceptably long, so our valueis as small as possible while keeping the wall-clock timemanageable. For N = 10 simulations, we adopt the op-timal softening length ǫ = 4 × − . Again, testing withsmaller softening lengths indicates that density profilesare largely unaffected, but integration times significantlylengthen.We have performed three types of GADGET simula-tions. For N = 10 , we have varied Q between 0.1and 1.0 for both types of initial density distributions andfor the three different velocity assignments when start-ing with clumpy initial conditions. We take these asthe standard set of simulations that provide zeroth-ordertests of our models. As the initial conditions are basedon random distributions of particles, we have also per-formed ensemble simulations of a subset of these initialconditions. Five independent realizations of initial con-ditions with 0 . ≤ Q ≤ . N = 10 particles. Both initial density distributions with0 . ≤ Q ≤ . N = 10 simulations.The GPU-NBODY-6 simulations evolve single andclumpy systems with N & , 0 . ≤ Q ≤ .
0, andboth initial density profiles. This code does not utilizesoftened forces, so it is a collisional code. However, thelarge number of particles provides a reasonable basis forthe assumption that two-body effects ( e.g., ejection) haverelatively minor impact on the simulations. For exam-ple, no particles gain escape speeds during any of oursimulations due to two-body collisions. Unfortunately,simulations with larger particle numbers could not becompleted with the current hardware available to the au-thors. 2.3.
Analysis
Independent of the evolution code, each simulation ex-tends for at least 20 initial crossing times. We have ob-served that this is generally a sufficient period for a sys-tem to reach a mechanical and virial equilibrium state.For clumpy initial conditions, individual clumps have“dissolved” by 10 initial crossing times, at the latest. Byconstruction, clumps in the hot-clumpy simulations be-gin to disperse almost immediately. Any simulation witha longer evolution will be noted in what follows. We in-fer mechanical equilibrium by observing that the radiusof the inner-most 90% of the particles stops varying andthat the average velocities of the particles are zero overthe radial range of the inner-most 90% of particles. Sim-ulations with Q & . ρ and rms speed v rms distribu-tions formed, we next compare to several model distri-butions. Our focus is on the comparison to the Lynden-Bell and Maxwell-Boltzmann models, but we also includetwo common analytical models, Plummer and de Vau-couleurs. The Plummer models we use have a density dis-tribution given by (Plummer 1911; Binney & Tremaine1987), ρ = ρ (cid:18) (cid:16) rr (cid:17) (cid:19) / , (1)where ρ is a scaling density and r is a scaling radius.This density profile is relatively constant in the core ofthe system and declines rapidly ( ∝ r − ) near the outeredge. The de Vaucouleurs profile (de Vaucouleurs 1948)has been used to fit the light profiles of elliptical galax-ies, and is a simple example of a broken power-law distri-bution (a cousin to commonly-discussed, cosmologically-motivated profiles such as, Navarro, Frenk, & White1996; Moore et al. 1999). It is well-established thatthe type of simulations discussed here result in struc-tures with outer density behavior that matches the de Vaucouleurs profile when Q . . ρ = ρ (cid:18) rr (cid:19) − δ (cid:18) rr (cid:19) δ − , (2)where δ = , and ρ and r are again a scaling densityand radius, respectively. A de Vaucouleurs density profilehas a central cusp, in contrast to the central density corebehavior of the Plummer model.As our models and simulated systems all have finitemasses, we demand that their density and v rms valuesmatch at the half-mass radius. With the connection be-tween model and simulated data values fixed, we use thereduced chi-squared statistic as the figure of merit forour fits, χ r = 1 N data N data X i =1 ( M i − D i ) ∆ i , (3)where N data = 100 for our 1% mass shells, D i is a massshell density or v rms value from a simulation, M i is amodel value corresponding to the same radial location,and ∆ i is an uncertainty estimate for the simulation value(as described in the first paragraph of this section). Oneshould expect, if the uncertainty estimates are appropri-ate, that a good model fit to the data produces χ r ≈ χ r values are determinedupon matching to simulation values at the half-mass ra-dius. For Lynden-Bell and Maxwell-Boltzmann models,a single parameter ( ν LB or ν MB ) determines the shapeof the density profile, and hence the match to the sim-ulation. Density profiles are determined by iterativelysolving the Poisson equation in straightforward fashion(Binney & Tremaine 1987, Sec. 4.4.2). The best-fit valueof ν is determined using an amoeba χ r minimization(Press et al. 1994). We have also performed MarkovChain Monte Carlo minimizations to corroborate theamoeba results.To determine model v rms distributions, we solve theJeans equation for a given density profile. This ap-proach demands a choice be made regarding the radialbehavior of the velocity anisotropy β ( r ). We utilize two v rms profiles: one assuming velocity isotropy β ( r ) = 0and one adopting β ( r ) from the simulated system. TheLynden-Bell and Maxwell-Boltzmann models are derivedassuming velocity isotropy, so their density profiles areconsistent only with the β ( r ) = 0 v rms profiles. In theabsence of a distribution function that incorporates themild tangential velocity anisotropy present in our sys-tems, we assume that the density derived from such afunction should be well-approximated by the density re-sulting from the isotropic version of the distribution func-tion. While not an exact description of the situation inour work, this assumption is consistent with the behaviorof the aniostropic Plummer model described in Merritt(1985). We choose the β ( r ) to be used in the Jeans equa-tion to be a smoothed version of the velocity anisotropypresent in the simulation. Specifically, a fourth-orderpolynomial fit to a simulation anisotropy profile is cre-ated. The parameters describing the Lynden-Bell and Barnes & EgererMaxwell-Boltzmann models are not allowed to vary dur-ing comparison to the v rms distribution, so the velocity χ r calculation for all models is straightforward once themodel and simulated v rms values are matched at the half-mass radius. ENTROPY BEHAVIOR3.1.
Microscopic Picture
The basis of the Barnes & Williams (2011, 2012)work is the phase-space counting approach outlined inLynden-Bell (1967). Here, we present a brief summaryof the chief ideas necessary for defining entropy fromthis microscopic viewpoint. Phase-space is imagined tobe sub-divided into two lattices; an array of nearly in-finitesimal micro-cells (each with volume ̟ ) that are ar-ranged into collections of macro-cells. Macro-cells con-tain ν micro-cells, and this value serves as the controlparameter for the models. Micro-cell occupation definesa fine-grained distribution function, while macro-cell oc-cupation defines a coarse-grained distribution functionthat can be realized through simulations. The occupa-tion of micro-cells determines the statistical propertiesof the system. Maxwell-Boltzmann statistics arise whenmultiple distinguishable particles can occupy a micro-cell. Lynden-Bell statistics describe situations in whicha classical exclusion principle disallows multiple particlesoccupying a single micro-cell. For collisionless systemslike the ones we investigate here, the fine-grained distri-bution function is a constant of motion. As such, oneexpects the Lynden-Bell statistical family to be most ap-propriate because the fine-grained distribution functionvalues cannot be increased through multiple occupancy.With these ideas, one can count the number of en-ergy states available to a system, and hence, define theentropy. The Lynden-Bell (1967) work relies on the Stir-ling approximation to simplify entropy calculations, butBarnes & Williams (2012) relax this assumption and al-low macro-cell occupation numbers to be small enoughthat the Stirling approximation fails. The end result isthat the Lynden-Bell entropy can be given as (Equation13 in Barnes & Williams 2012), S LB = S LB , − k B M X i =1 [( n i + 1 /
2) ln ( n i + 1)+( ν − n i + 1 /
2) ln ( ν − n i + 1) + λ ,n i + λ , ( ν − n i ) (cid:3) , (4)where S LB , = k B [ N ln N − N + M ( ν ln ν − ln 2 π )], N isthe number of phase-space elements/particles, and M isthe total number of macro-cells. The λ function arisesfrom the approximation,ln x ! = ( x + 12 ) ln ( x + 1) − x + ln 2 π λ ,x (5)where λ ,x = − ( x + 2 x + )( x + x + ) . (6)A similar expression can be found for theMaxwell-Boltzmann entropy (see Equation A2 inBarnes & Williams 2012). Adopting these modificationsleads to the possibility of finite-mass and energy systems that belong to the Maxwell-Boltzmann statistical family,in contrast to the findings of Lynden-Bell (1967). Like-wise, the discussion in Binney & Tremaine (1987, § f ln f term in the entropycalculation is now modified. Our Lynden-Bell expressiondoes not change the overall character of the associateddistribution function presented in Lynden-Bell (1967);it remains finite-mass and energy and closely resemblesa Fermi-Dirac distribution.We have calculated S LB using the results of GADGET N = 10 simulations. At every timestep, positions andvelocities are used to assign each particle to a macro-cell,giving n i . A fixed value of ν LB = 10 has been chosenfor these calculations, as that is always greater than themaximum n i value in these simulations. Tests varying ν LB show that the “zero point” of S LB is affected muchmore strongly than its time-dependent behavior . Resultsof these calculations are discussed in Section 5.3.2. Macroscopic Picture
As a complement to the microscopic approach, we fol-low the discussion of thermal non-equilibrium situationsgiven by de Groot & Mazur (1984). The behavior of self-gravitating systems composed of a large number of mas-sive particles can be described using equations that rep-resent macroscopic conservation laws. Expressions of theconservation of energy can be manipulated and combinedwith the first law of thermodynamics to provide insightinto the behavior of entropy. In particular, it is useful towrite the entropy production rate of a system as, ∂S ( c ) ∂t = Z V σ d x, (7)where σ is the entropy production per unit volume perunit time and the integral is taken over the system vol-ume. The specific entropy production rate in this pictureis, σ = − T q · ∇ T K − T K Π ↔ : ∇ v , (8)where T K is the kinetic temperature, q is the heat con-duction flux, Π ↔ is the anisotropic pressure tensor, and v is the mean velocity. As usual, the kinetic tempera-ture is a measurement of the random kinetic energy ina small region; T K = ( m/ k B ) h v p i , where k B is Boltz-mann’s constant and v p is the magnitude of the peculiarvelocity. The heat conduction flux q = ρ h v p v p i repre-sents the peculiar kinetic energy that is transported bypeculiar velocities in the system. Interested readers mayfind details of the calculation of σ in de Groot & Mazur(1984, § N = 10 particles broken into a spherical grid with 10 radial, 8ntropy Production 5polar, and 8 azimuthal bins, each will contain on the or-der of 100 particles. The need for gradients in quantitiesplaces some constraint on how coarse the spherical gridcan be made. We use simulations with N = 10 parti-cles and the aforementioned grid resolution to begin toinvestigate this macroscopic picture of entropy behaviorduring an evolution. Results of these calculations arediscussed in Section 5. DENSITY AND VELOCITY PROFILE FITTINGWe take the GADGET simulations with N = 10 asthe standards for the majority of our results. These simu-lations provide a wide-range of evolutions from which wedraw broad inferences. The other simulations (ensemble N = 10 GADGET, N = 10 GADGET, GPU-NBODY-6) are focused on testing relationships and behaviors sug-gested by the standard set.4.1.
Individual N = 10 GADGET Simulations
Lynden-Bell and Maxwell-Boltzmann models can pro-vide excellent fits to the density and velocity profiles ofthe final equilibrium states of single and hot-clumpy sys-tems when Q & .
6. However, independent of the typeof particle distribution (single or clumpy) and initial den-sity profile, systems with Q . . Q and/or cold-clumpy initial conditions, the fi-nal density distributions are best-described by a de Vau-couleurs profile, due to its cuspy nature. However, thespecifics of the central density cusp do not always agreewith the δ = 1 / δ vary in a few cases, finding that1 / . δ . /
4. The thermodynamic models investi-gated here simply cannot reproduce the central densitycusp. Interestingly, Plummer models provide good de-scriptions of the density profiles of warm-clumpy systemswhen Q & . Q = 1 .
0, while Figure 2 showsresults from evolutions with Q = 0 .
7. Analogous plotsfor simulations with Q = 0 . Q = 0 . Q ≥ .
9, with themodels reversing positions for Q ≤ .
8. Overall, LBmodels appear to do a better job of describing the den-sity distributions of the simulated equilibria.Corresponding plots for the v rms profiles are given inFigures 3 and 4. In these figures, the profiles createdassuming isotropic and anisotropic velocity distributions show quite different behaviors. Finding v rms ( r ) from theJeans equation under the assumption that β ( r ) = 0 re-sults in profiles that have flat cores. However, single sys-tems generically have v rms profiles with non-zero slopesnear their centers. Hot-clumpy systems are qualitativelysimilar to the isotropic model predictions, but are notterribly well described by the models. The set of thinlines below the data points illustrate the v rms ,r , v rms ,θ ,and v rms ,φ behavior of the simulation. The differencesbetween the line shapes indicates mild tangential veloc-ity anisotropy. Allowing this β ( r ) profile as input to theJeans equation results in v rms profiles that dramaticallyincrease the quality of the fits. We note that the LBmodels tend to provide better fits to the simulated v rms profiles than those produced by the MB models.4.2. Other Simulations
Fits to the standard simulations suggest that LB mod-els are generally better than MB models. We now beginto test the robustness of this observation using “aver-age” systems formed by ensembles of simulations withdifferent realizations of the same initial conditions using N = 10 particles. As with the standard simulations,GADGET has been used to evolve the initial conditions.Figures 5 and 6 are analogous to Figures 1 and 2. Themost significant changes one notices is that the simula-tion profiles are smoother and have smaller error bars.For nearly every simulation, LB model fits to the den-sity distributions are superior to those provided by theMB model (for Q = 1, the two models can producecomparable fits). The same is true for velocity profiles.Figures 7 and 8 show the improvement that LB modelsprovide over MB models. As with the standard simu-lations, the inclusion of velocity anisotropy significantlyimproves agreement between the model curves and thesimulation results.The equilibrium structures in simulations with N =10 particles (returning to only one realization per ini-tial condition) have also been analyzed in a similar man-ner. Figures 9 and 10 show how the LB and MB densityprofiles compare to the simulation results. The singlesystem profiles (particularly for the cuspy initial profile)display small-scale variations that are larger than theuncertainty estimates we have made. The systems arein virial equilibrium, and the near-zero average velocityvalues for all components indicate mechanical equilib-rium as well. Additional evolutions of these initial condi-tions with different (smaller) softening lengths have pro-duced very similar outcomes. Extending the evolutionsto longer times does reduce the variations somewhat, andwe have used the results of our longest evolutions (out to30 T ) to create the relevant figures. Somewhat surpris-ingly, we do not see similar structures in the profiles ofthe hot-clumpy simulations, which are very smooth. Asthese new features suggest that our naive uncertainty es-timates are questionable, we do not place much stock inthe actual values of χ r that we have found. However, weargue that since the LB and MB models are both beingcompared to the same data with the same uncertainties,their relative χ r values distinguish between which is themore appropriate model.Unlike the averaged simulations, the appropriateness ofthe LB model is not obvious here. In most cases, the LBand MB models provide comparable fits to the simulation Barnes & Egererdensity profiles; only for the hot-clumpy simulations withGaussian initial density profiles is the LB model clearlypreferred. Likewise, the differences between the LB andMB velocity profiles are more subdued in Figures 11 and12 compared to previous versions. It also appears thatthere are more significant differences between the dataand the LB model curves in these figures.Overall, the GADGET-based simulations support theidea that the LB model does a good job describing thedensity and velocity profiles of the equilibria of colli-sionless systems that have undergone mild gravitationalpotential relaxation. Our final test for this idea is touse the different evolution code GPU-NBODY-6. UnlikeGADGET, GPU-NBODY-6 does not incorporate grav-itational softening and instead treats particles as truepoint masses. While two-body encounters do occur, thelarge particle number ( N > ) minimizes the global im-pact of particle-particle interactions, leaving a basicallycollisionless evolution. The results of these simulationsare similar to those from the standard simulations. Ingeneral, LB and MB models provide comparable descrip-tions of the density profiles of systems with Q ≥ . Q ≤ . v rms profiles produced by thesesimulations. ENTROPY PRODUCTIONAs mentioned in section 3.1, the microscopic calcula-tion of the entropy has been carried out for the N = 10 GADGET simulations. We have created two differentmacro-cell grids; one with 10 macro-cell divisions perphase-space dimension ( M = 10 ), and one with 15 di-visions ( M ≈ ). We show a representative pair of S LB ( t ) curves in Figure 17a. The different curves corre-spond to the grid choices indicated. These are derivedfrom a single, cuspy Q = 0 . Q , Gaussian densityprofiles, and clumpy particle distributions. Most impor-tantly, it is clear that the entropy rises from an initialvalue and rapidly approaches a maximum, steady-statevalue. The most rapid increase (a roughly 5% change) inentropy occurs during the first couple of initial crossingtimes. It seems natural that the most rapid growth in S LB occurs during the period when the strongest gravi-tational potential relaxation (biggest variations in poten-tial) occur. For comparison, the variation in the virialratio Q = 2 K/ | W | as a function of time in the samesimulation is shown in Figure 17b.We also note that the finer macro-cell grid producessmaller variations in the value of S LB , a trend that alsooccurs as Q →
1. With the increasing smoothness of thecurves with higher Q , one can also discern that thereis slower growth of S LB occurring over tens of cross-ing times. We speculate that this slower growth is at-tributable to the phase mixing that continues after theinitial potential variations decrease.While the microscopic entropy picture dovetails nicelywith the overall scenario of entropy production in colli-sionless systems, the macroscopic picture results are notso clear. Figure 18 illustrates the behavior of Equation 7for the single, cuspy Q = 0 . N = 10 GADGET simulations. Given the relianceon derivative information required by this approach, thenoise present in the values is not surprising. It is possiblethat our particle number and grid resolutions are simplynot high enough to capture the true behavior of the en-tropy creation term. Preliminary testing with higher gridresolution did not produce appreciably different results.However, there could also be a fault in the assumptionof local thermodynamic equilibrium that underlies thederivation of Equation 8. CONCLUSIONSWe have evolved sets of N -body initial conditions todetermine if the thermodynamically-motivated Lynden-Bell or Maxwell-Boltzmann distribution functions candescribe the equilibrium distributions of particle loca-tions and velocities. Different evolution codes and nu-merical parameters have been adopted to reduce the like-lihood of spurious findings. These simulation results havealso been used to investigate the entropy behavior ofthese systems.Initially “cold” systems evolve in expected fashion,forming equilibria that show cuspy central density pro-files and outer density profiles that match the de Vau-couleurs form, not the thermodynamic models. How-ever, a subset of our simulations can be well-described bydistribution functions that maximize entropy. Initially“hot” systems form equilibria with density profiles thatare very similar to those produced by the Lynden-Belldistribution function. Unfortunately, the lack of velocityisotropy in these simulations seems to preclude the abilityof the LB distribution functions to predict v rms profiles.However, accounting for the mild tangential anisotropyproduces extremely good descriptions of simulation ve-locity results. Previous simulations involving mild ve-locity anisotropy indicate very similar global behaviorto fully isotropic systems ( e.g., Merritt & Aguilar 1985).We argue that the isotropic behavior is of central impor-tance, with velocity anisotropy determining higher-ordercorrections to density and v rms profiles ( e.g., Merritt1985). To be clear, introducing velocity anisotropy out-side the distribution function, as we have done, meansthat such models are not self-consistent. The anisotropic v rms profiles are not predictions of the maximum entropyargument, which assumes only constant system mass andenergy. Anisotropic models would require entropy max-imization including at least one other constraint, alongthe lines of (Trenti & Bertin 2005). Given these caveats,our results suggest that there are conditions under whicha collisionless self-gravitating system will evolve to amaximum entropy state. This conclusion is strengthenedby the fact that direct calculation of the entropy (usinga phase-space occupation approach) also shows a maxi-mum value being attained. An alternative constructionof the entropy behavior is inconclusive, presumably dueto resolution effects.There remain some unresolved issues with the idea ofthese equilibria representing maximum entropy states.One is that maximized entropy is normally taken to im-ply thermodynamic equilibrium, but these simulated sys-ntropy Production 7tems clearly have kinetic temperature gradients. Is itpossible that in self-gravitating systems, these two condi-tions are not equivalent? Another set of questions revolvearound the role of the ν parameter. This value representsthe number of micro-cells that occupy any macro-cell,but it does not have an ab initio value. In the directcalculation, it must be larger than the largest macro-cell occupation number n i in order for the entropy to bewell-defined. In fitting density profiles, we have placed nosuch restriction on its value. For the simulations focusedon in this work, we have found 5 . ν . Q values couple to lower ν values. On the otherhand, the slight positive concavity seen in the outer den- sity profiles of Q = 0 . ν & REFERENCESBarnes, E. I., Lanzel, P. A., Williams, L. L. R. 2009, ApJ, 704,372Barnes, E. I., Williams, L. L. R. 2011, ApJ, 728, 136Barnes, E. I., Williams, L. L. R. 2012, ApJ, 748, 144Binney, J., Tremaine, S. 1987, Galactic Dynamics, (Princeton,NJ:Princeton)de Groot, S. R., Mazur, P. 1984, Non-EquilibriumThermodynamics, (Mineola, NY:Dover)de Vaucouleurs, G. 1948, Ann. d’Astrophys., 11, 247Lynden-Bell, D. 1967, MNRAS, 136, 101Merritt, D. 1985, AJ, 90, 1027Merritt, D., Aguilar, L. A. 1985, MNRAS, 217, 787Moore, B., Quinn, T., Governato, F., Stadel, J., Lake, G., 1999,MNRAS, 310, 1147Navarro, J. F., Frenk, C. S., White, S. D. M. 1996, ApJ, 462, 563Navarro, J. F., Frenk, C. S., White, S. D. M. 1997, ApJ, 490, 493Navarro, J. F., Hayashi, E., Power, C., Jenkins, A. R., Frenk, C.S., White, S. D. M., Springel, V., Stadel, J., Quinn, T. R. 2004,MNRAS, 349, 1039 Nitadori, K., Aarseth, S. 2012, MNRAS, 424, 545Plummer, H. C. 1911, MNRAS, 72, 460Power, C., Navarro, J. F., Jenkins, A., Frenk, C. S., White, S. D.M., Springel, V., Stadel, J., Quinn, T. 2003, MNRAS, 338, 14Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P.1994, Numerical Recipes, Cambridge Univ. Press: New York,NYSpringel, V. 2005, MNRAS, 364, 1105Springel, V., White, S. D. M., Jenkins, A., Frenk, C. S., Yoshida,N., Gao, L., Navarro, J., Thacker, R., Croton, D., Helly, J.,Peacock, J. A., Cole, S., Thomas, P., Couchman, H., Evrard,A., Colberg, J., Pearce, F. 2005, Nature, 435, 629Sylos Labini, F. 2012, MNRAS, 423, 1610Taylor, J. E., Navarro, J. F. 2001, ApJ, 563, 483Trenti, M., Bertin, G. 2005, A&A, 429, 161van Albada, T. S. 1982, MNRAS, 201, 939
Barnes & Egerer
Fig. 1.—
Logarithmic density profiles for individual N = 10 GADGET simulations with Q = 1 .
0. The radial coordinate is scaledby the half-mass radius (indicated by the vertical dashed line) in all profiles. In each panel, the initial conditions are specified and theerrorbars indicate the data values from the simulations. The curves show the behaviors of the various models, specified in the legend. Asthey provide poor descriptions of the simulations, the Plummer and de Vaucouleurs models are only included in the single cuspy panel forreference. Lynden-Bell (LB) model fits produce density χ r values smaller or comparable to those with Maxwell-Boltzmann (MB) models(except for the hot-clumpy Gaussian simulation): single cuspy – χ = 0 . χ = 0 . χ = 0 . χ = 2 . χ = 1 . χ = 0 . χ = 5 . χ = 0 . ntropy Production 9 Fig. 2.—
Logarithmic density profiles for individual N = 10 GADGET simulations with Q = 0 .
7. The panels and linestyles are thesame as in Figure 1. In these simulations, the density profiles share a slight upward concavity in their outer profile. In general, the LBmodel fits are able to match this data behavior better than the MB models: single cuspy – χ = 1 . χ = 9 . χ = 0 . χ = 4 . χ = 0 . χ = 4 . χ = 3 . χ = 3 . Fig. 3.—
Logarithmic v rms profiles for individual N = 10 GADGET simulations with Q = 1 .
0. Again, the radial coordinate is scaled bythe half-mass radius (vertical dashed line). As before, the initial conditions are specified in each panel, and the errorbars indicate the datavalues from the simulations. The thin lines that appear below the errorbars show the behaviors of v rms ,r (the lowest line) and the nearlyidentical v rms ,θ and v rms ,φ . From these components, one can see the mild tangential velocity anisotropy that exists in these simulatedequilibria. Two sets of model curves are shown superimposed with the data. As indicated in the legend in the upper-left panel, there areisotropic and anisotropic model results. Isotropic model profiles reach nearly constant values near the centers of systems, while anisotropicmodel profiles provide better representations of the data for smaller r . In general, anisotropic LB models produce the smallest velocity χ r values: single cuspy – χ , iso = 1 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . ntropy Production 11 Fig. 4.—
Logarithmic v rms profiles for individual N = 10 GADGET simulations with Q = 0 .
7. The panels here are analogous tothose in Figure 3. The ineffectiveness of isotropic models remains evident, and the superiority of anisotropic LB models is even clearerthan in the Q = 1 . χ , iso = 2 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 1 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 3 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 1 . Fig. 5.—
Logarithmic density profiles for averaged N = 10 GADGET simulations with Q = 1 .
0. The panels are analogous to thosein Figure 1. The errorbars marking the data points are smaller and the data points show less point-to-point variation than in Figure 1.As before, the LB models tend to provide better descriptions of the data: single cuspy – χ = 4 . χ = 4 . χ = 0 . χ = 3 . χ = 0 . χ = 0 . χ = 1 . χ = 0 . ntropy Production 13 Fig. 6.—
Logarithmic density profiles for averaged N = 10 GADGET simulations with Q = 0 .
7. The panels are analogous to thosein Figure 2, and the data profiles present smoother versions of the outer-profile concavity features seen there. Again, the LB modelstend to describe the data behavior more completely than the MB models: single cuspy – χ = 7 . χ = 44 . χ = 3 . χ = 28 . χ = 0 . χ = 7 . χ = 2 . χ = 9 . Fig. 7.—
Logarithmic v rms profiles for averaged N = 10 GADGET simulations with Q = 1 .
0. The panels are analogous to those inFigure 3. These smoother versions again indicate that models incorporating the velocity anisotropy present in the simulations are bettersuited to describing the data. Of the models considered here, the anisotropic LB models provide the best representation of the data: singlecuspy – χ , iso = 5 . χ , aniso = 0 . χ , iso = 5 . χ , aniso = 0 . χ , iso = 6 . χ , aniso = 0 . χ , iso = 9 . χ , aniso = 3 . χ , iso = 2 . χ , aniso = 0 . χ , iso = 3 . χ , aniso = 0 . χ , iso = 3 . χ , aniso = 0 . χ , iso = 3 . χ , aniso = 0 . ntropy Production 15 Fig. 8.—
Logarithmic v rms profiles for averaged N = 10 GADGET simulations with Q = 0 .
7. These panels are analogous tothose in Figure 4. The smaller errorbars in these simulations result in some poorer fits between the anisotropic LB model and thedata ( e.g., in the clumpy Gaussian simulation, lower-right panel): single cuspy – χ , iso = 8 . χ , aniso = 0 . χ , iso = 3 . χ , aniso = 1 . χ , iso = 4 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 2 . χ , iso = 2 . χ , aniso = 0 . χ , iso = 6 . χ , aniso = 10 . χ , iso = 2 . χ , aniso = 0 . χ , iso = 2 . χ , aniso = 3 . Fig. 9.—
Logarithmic density profiles for individual N = 10 GADGET simulations with Q = 1 .
0. These panels are analogous to thosein Figure 1. The point-to-point variations seen in the single simulations can be relatively large, leading us to question the appropriatenessof our estimated data point uncertainties. However, the clumpy simulation results are quite smooth, and at least on a relative basis, theLB models describe the data better than the MB models: single cuspy – χ = 15 . χ = 16 . χ = 6 . χ = 7 . χ = 0 . χ = 0 . χ = 0 . χ = 1 . ntropy Production 17 Fig. 10.—
Logarithmic density profiles for individual N = 10 GADGET simulations with Q = 0 .
7. These panels are analogous tothose in Figure 2. Again, the data points for clumpy simulations show a smoothness not present in the single simulations. Except for thesingle cuspy simulation, the LB models clearly provide a superior description of the data: single cuspy – χ = 81 . χ = 105 . χ = 3 . χ = 32 . χ = 1 . χ = 28 . χ = 2 . χ = 22 . Fig. 11.—
Logarithmic v rms profiles for individual N = 10 GADGET simulations with Q = 1 .
0. These panels are analogous to thosein Figure 3. The models with velocity isotropy continue to be poor descriptors of the data, while the anisotropic LB models provide thebest fits to the simulation results: single cuspy – χ , iso = 16 . χ , aniso = 2 . χ , iso = 15 . χ , aniso = 3 . χ , iso = 12 . χ , aniso = 3 . χ , iso = 12 . χ , aniso = 3 . χ , iso = 7 . χ , aniso = 0 . χ , iso = 6 . χ , aniso = 0 . χ , iso = 3 . χ , aniso = 0 . χ , iso = 4 . χ , aniso = 1 . ntropy Production 19 Fig. 12.—
Logarithmic v rms profiles for individual N = 10 GADGET simulations with Q = 0 .
7. These panels are analogous to thosein Figure 4. Again, the anisotropic LB models provide the best descriptions of the data (even though it appears to under-predict thesingle cuspy simulation results for larger r ): single cuspy – χ , iso = 12 . χ , aniso = 1 . χ , iso = 9 . χ , aniso = 2 . χ , iso = 17 . χ , aniso = 0 . χ , iso = 7 . χ , aniso = 2 . χ , iso = 13 . χ , aniso = 0 . χ , iso = 5 . χ , aniso = 2 . χ , iso = 15 . χ , aniso = 0 . χ , iso = 7 . χ , aniso = 2 . Fig. 13.—
Logarithmic density profiles for individual N = 2 GPU-NBODY-6 simulations with Q = 1 .
0. These panels are analogous tothose in Figure 1. Evolving systems with this non-softened code produces equilibria with density profiles that are nearly indistinguishablefrom those discussed earlier. As with the GADGET results, LB models provide better representations of the data: single cuspy – χ =1 . χ = 2 . χ = 0 . χ = 1 . χ = 1 . χ = 0 . χ = 0 . χ = 0 . ntropy Production 21 Fig. 14.—
Logarithmic density profiles for individual N = 2 GPU-NBODY-6 simulations with Q = 0 .
7. These panels are analogousto those in Figure 2. The fact that the LB models can reproduce the slight concave features in the outer regions of the profiles again leadto them being preferred to the MB models: single cuspy – χ = 3 . χ = 17 . χ = 1 . χ = 13 . χ = 3 . χ = 13 . χ = 0 . χ = 8 . Fig. 15.—
Logarithmic v rms profiles for individual N = 2 GPU-NBODY-6 simulations with Q = 1 .
0. These panels are analogousto those in Figure 3. The trend for anisotropic LB models to provide the best description of the v rms data continues. Note that theclumpy cuspy simulation results in a nearly isotropic velocity distribution, and that the inner part of the profile is decently well-describedby the isotropic LB prediction: single cuspy – χ , iso = 1 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . ntropy Production 23 Fig. 16.—
Logarithmic v rms profiles for individual N = 2 GPU-NBODY-6 simulations with Q = 0 .
7. These panels are analogousto those in Figure 4. The anisotropic LB models continue to represent the data behavior better than the other models considered: singlecuspy – χ , iso = 1 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 1 . χ , aniso = 1 . χ , iso = 0 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 1 . χ , iso = 1 . χ , aniso = 0 . χ , iso = 0 . χ , aniso = 1 . Fig. 17.— (a) Lynden-Bell entropy versus time in the individual single, cuspy N = 10 GADGET simulation with Q = 0 .
7. The entropyis calculated from the phase-space macro-cell occupation values (Equation 4), with a value of ν LB = 10 . Two different macro-cell volumes(indicated in the legend) have been used to create comparison curves. Adopting a smaller macro-cell volume increases the zero point ofthe curve and reduces the small-scale variations in S LB but leaves the overall behavior unchanged. (b) The virial ratio Q as a function oftime for the same simulation. The initial increase in S LB appears to coincide with the initial growth in Q . Fig. 18.—
The macroscopic entropy production rate versus time in the individual single, cuspy N = 10 GADGET simulation with Q = 0 .
7. The production rate is calculated by determining gradients in the kinetic temperature and mean velocity field (Equation 8).The rather coarse grid used to determine the gradients may play an important role in the essentially null result shown. If the macroscopicentropy production could be properly calculated, one would expect a rather large positive spike in the interval 0 ≤ t/T ..