Table Of Content

Astronomy&Astrophysicsmanuscriptno.SormaniBertinRevised (cid:13)c ESO2013 January28,2013 Gravothermal Catastrophe: the dynamical stability of a fluid model M.C.Sormani1 andG.Bertin2 1 ScuolaNormaleSuperiore,PiazzadeiCavalieri7,56126Pisa,Italy(cid:63) e-mail:[email protected] 2 Universita`degliStudidiMilano,DipartimentodiFisica,viaCeloria16,I-20133Milano,Italy e-mail:[email protected] ReceivedXXX;acceptedYYY ABSTRACT 3 Are-investigationofthegravothermalcatastropheispresented.Bymeansofalinearperturbationanalysis,westudythedynamical 1 stabilityofasphericalself-gravitatingisothermalfluidoffinitevolumeandfindthattheconditionsfortheonsetofthegravothermal 0 catastrophe,underdifferentexternalconditions,coincidewiththoseobtainedfromthermodynamicalarguments.Thissuggeststhat 2 thegravothermalcatastrophemayreducetoJeansinstability,rediscoveredinaninhomogeneousframework.Wefindnormalmodes n andfrequenciesforthefluidsystemandshowthatinstabilitydevelopsonthedynamicaltimescale.Wethendiscussseveralrelated a issues.Inparticular:(1)Forperturbationsatconstanttotalenergyandconstantvolume,weintroduceasimpleheuristicterminthe J energybudgettomimictheroleofbinaries.(2)Weoutlinetheanalysisofthetwo-componentcaseandshowhowlinearperturbation analysiscanbecarriedoutalsointhismorecomplexcontextinarelativelystraightforwardway.(3)Wecomparethebehaviorofthe 5 fluidmodelwiththatofthecollisionlesssphere.Inthecollisionlesscasetheinstabilityseemstodisappear,whichisatvariancewith 2 thelinearJeansstabilityanalysisinthehomogeneouscase;wearguethatakeyingredienttounderstandthedifference(aspherical stellarsystemisexpectedtoundergothegravothermalcatastropheonlyinthepresenceofsomecollisionality,whichsuggeststhatthe ] A instabilityisdissipativeandnotdynamical)liesintheroleofthedetailedangularmomentuminacollisionlesssystem. Finally, we briefly comment on the meaning of the Boltzmann entropy and its applicability to the study of the dynamics of self- G gravitatinginhomogeneousgaseoussystems. . h Keywords.hydrodynamics–instabilities–gravitation–galaxies:clusters:general p - o r 1. Introduction models, collisionless models, and fluid models) admit equilib- t s rium configurations that are spatially truncated self-gravitating a Corecollapseinan N-bodysystemisaproblemrelevanttothe isothermalspheres,butitisnotcleartowhatextenttheycanbe [ dynamics of globular clusters, because these are considered to consideredasrepresentativeoftheequilibriumstatesofthepure betheonlystellarsystemsthatpossessthenecessarydegreeof 1 N-body problem. Several studies have addressed the stability collisionalitytorelaxthermally.Thediscoveryofsomeclusters v problemofisothermalspheresusingthermodynamicalmethods, withcuspycoreswasinterpretedasasignthattheyindeedexpe- 8 startingwithAntonov(1962)andLynden-Bell&Wood(1968) 3 riencedthegravothermalcatastrophe.Inthecontextofglobular andcontinuingwithalonglistofpapers(forexample,Hachisu 0 clusters,thethreemainphenomenacausedbycollisionalityare: & Sugimoto 1978; Nakada 1978; Katz 1978; Inagaki 1980; 6 core collapse, evaporation, and mass segregation. The focus of Padmanabhan 1989; Chavanis 2002, 2003, see also Thirring . thispaperisoncorecollapse. 1 1970).Inthethermodynamicalapproachthestudyofisothermal Consider the gravitational N-body problem, where N (cid:29) 1. 0 spheresisbasedonaformofentropyknownastheBoltzmann 3 Thatis,considerasetofN classicalpointmasses,eachofmass entropy: 1 m,mutuallyinteractingthroughNewtoniangravity.Theparticles : mayormaynotbeconfinedwithinasphericalvolumeofradius (cid:90) v R.Weareinterestedinthefollowingquestions: Sb[f]≡−k f(r,v)ln f(r,v)d3rd3v, (1) i X – Can the system reach some sort of equilibrium? Can this r where f is the one-particle distribution function,1 r and v the a equilibriumbecalledthermal? position and velocity vector respectively, and k the Boltzmann – Can a model (for example, collisionless or fluid) of the N- constant. bodyproblemreachequilibrium?Howdoestheequilibrium Unfortunately, the thermodynamics of self-gravitating sys- ofamodelrelatetotheequilibriumofthepureN-bodyprob- temsstilldependsonanumberofunresolvedissues,partlybe- lemandofarealstellarsystem? cause of the long-range nature of the force and partly because – Aretheseequilibriastable? ofthedivergentbehavioroftheforceatshortdistance(e.g.,see Thistopichasbeenstudiedbymanyauthors(forrecentreviews Padmanabhan1990;Katz2003;Chavanis2006;Mukamel2008; seeforexampleHeggie&Hut2003;Binney&Tremaine2008). Campaetal.2009;Bouchetetal.2010).Acriticalanalysisofthe Different models of the N-body problem (in particular gaseous useoftheBoltzmannentropyhasbeenmadebyMiller(1973). WewillbrieflycommentonthispointinSect.2. (cid:63) Present address: Rudolf Peierls Centre for Theoretical Physics, 1 (cid:82) KebleRoad,OxfordOX13NP 1 Normalizedinsuchawaythat fd3rd3v=N 1 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel Therefore,itwouldbeinterestingtostudythegravothermal Thepaperisorganizedasfollows.InSect.2wecommenton catastrophe by limiting as much as possible the use of thermo- the meaning of the Boltzmann entropy. In Sect. 3 we write the dynamicalarguments.Inthispaperwereconsiderthegravother- basic hydrodynamic equations in the Eulerian and Lagrangian mal catastrophe and analyze the dynamical stability of self- representation.InSect.4weshowtheresultsofthelinearanaly- gravitatingfluidsgovernedbytheEulerequationbyfindingthe sisofthehydrodynamicequationsbystudyingthepropertiesof normalmodesandfrequenciesofsphericalsystemsundervari- spherically-symmetric perturbations. We find the relevant nor- ousexternalconditions,inparticular:constanttotalenergyEand malmodesandfrequenciesunderdifferentboundaryconditions. volumeV(thiscasecorrespondstothegravothermalcatastrophe In the analysis of the gravothermal catastrophe, we propose a ofLynden-Bell&Wood1968),constanttemperatureT andvol- modifiedexpressionfortheenergy,whichcanbeconsideredas umeV(isothermalcollapse),constanttemperatureT andbound- a simple way to incorporate the energy generation from bina- ary pressure P (isobaric collapse, Bonnor 1956; Ebert 1955). ries in the linear regime; we show that the catastrophe can in- Note that, in contrast to thermodynamical arguments, the lin- deedbehaltedifweconsidersuchamodifiedexpression.Inthe earmodalanalysisautomaticallygivesthetimescaleforthede- lastsubsectionandinAppendixE,webrieflyconsiderthetwo- velopment of the instability. The constant {T,V} case has been componentcase.InSect.5.1wediscusstherelevanttimescales. addressedbySemelinetal.(2001),whostudiedthestabilityof In Sect. 5.2 we discuss the difference between the collisionless thesystemnumerically,andbyChavanis(2002),whofoundan model and the fluid model of the N-body problem. In Sect. 6 analytical solution for the marginally stable perturbations us- wedrawourconclusionsandidentifysomeopenquestionsand ing methods developed by Padmanabhan (1989). The constant issues. {T,P} case has been addressed previously with synthetic argu- mentsbyBonnor(1956)andEbert(1955),byYabushita(1968), 2. Thethermodynamicalapproachandthe whostudiedthestabilityofthesystemnumerically,byLombardi & Bertin (2001), who extended the analysis by Bonnor (1956) Boltzmannentropy to the nonspherically symmetric case, and by Chavanis (2003), Antonov (1962) and Lynden-Bell & Wood (1968) based their who gave an analytical solution for the case of marginally sta- stability analysis of self-gravitating spatially-truncated isother- bleperturbationsusingthesamemethodashehadusedforthe mal spheres on the use of the Boltzmann entropy (1). Starting constant {T,V} case. The constant {E,V} case has been studied fromakineticdescription,theylookedforstationarystatesofthe usingamodelbasedontheSmoluchowski-Poissonsystem(dif- Boltzmannentropywithrespecttothedistributionfunction f,at ferentfromtheEuler-Poissonsystemconsideredinthispaper), fixed2 total mass M, total energy E, and volume V = 4πR3/3, byChavanisetal.(2002). whereRistheradiusofasphericalbox.Thesestationarystates Inthispaper,weextendtheseanalysesoftheconstant{T,V} arespatiallytruncatedisothermalspheres.Itisfoundthat,fora and {T,P} cases to the general calculation of eigenfrequencies givensingle-particlemassm,eachspatiallytruncatedisothermal andeigenfunctionsforconditionsoutsidethoseofmarginalsta- sphere is identified by two dimensional scales (e.g., the central bility, and study in detail the constant {E,V} case. We use a densityρ (0)andthetemperatureT)andonedimensionlesspa- 0 EulerianorLagrangianrepresentationofhydrodynamicsassug- rameter{forexampleΞ=R/λ,whereλ=[kT/m4πGρ (0)]1/2is 0 gestedbytheboundaryconditionstobeimposedonthesystem reminiscentoftheJeanslength(Jeans1902);anotherequivalent and we provide a unified treatment for all cases. Surprisingly, choice for the dimensionless parameter is the density contrast wefindthatthesystembecomesdynamicallyunstableinallthe ρ (0)/ρ (R),thatis,theratioofthecentraldensitytothedensity 0 0 cases considered exactly at the same points found by Lynden- atthetruncationradiusR}.Conversely,eachchoiceofthethree Bell & Wood (1968) by means of a thermodynamical analysis, quantities{T,ρ(0),Ξ}identifiesoneparticularisothermalsphere. thatisforvaluesofthedensitycontrast(theratioofthecentral Thisone-to-onecorrespondencemaybelostforotherchoicesof densitytotheboundarydensity)≈14,32.1,and709respectively thequantitiescharacterizingthesystem. fortheconstant{T,P},{T,V},and{E,V}case.Theseresultssug- In this thermodynamical approach the necessary condition gest that the gravothermal catastrophe may reduce to the Jeans forstabilityisthatthestationarystatecorrespondstoa(localat (1902) instability, rediscovered in the inhomogeneous context. least)maximumoftheBoltzmannentropyfunctional.Antonov Thenwebrieflyoutlinetheproblemofthelineardynamicalsta- (1962) found that if Ξ < 34.4 [i.e., if ρ (0)/ρ (R) < 709] the 0 0 bilityofasphericalfluidgeneralizedtothetwo-componentcase, equilibrium configurations are local maxima, whereas for Ξ > byperturbingthetwo-componentspatiallytruncatedisothermal 34.4 they are saddle points. Therefore, he concluded that self- sphereconfigurationspreviouslyconsideredbyTaffetal.(1975), gravitatingisothermalspheresareunstableifΞ>34.4. Lightman (1977), Yoshizawa et al. (1978), de Vega & Siebert TheworkofAntonov(1962)wasextendedbyLynden-Bell (2002), who studied the stability of these systems with a ther- &Wood(1968).Theseauthorsstudiedthethermodynamicalsta- modynamicalapproach,andSopiketal.(2005),whoaddressed bility of the system under various conditions by applying the alsothedynamicalstabilitywithamodeldifferentfromtheone relevant thermodynamical potentials (based on the Boltzmann used in this paper. We show that the onset of thermodynamical entropy) and gave a physical interpretation of the instability anddynamicalstabilityoccursatthesamepointinthesimplest in terms of negative specific heats, as a characteristic feature case (constant {T,V}) and confirm that the component made of of self-gravitating systems, thus creating the paradigm of the heavierparticlesistheprimarydriveroftheinstability.Thisre- gravothermal catastrophe. They found that the critical value of sultislikelytoberelatedtoasimilarfindingbyBreen&Heggie Ξ that corresponds to the onset of instability depends on the (2012a,b)fortwo-componentgravothermaloscillations.Finally (cid:82) we focus on the following puzzling phenomenon: in the colli- 2 Intermsofthedistributionfunction f,wehaveM =m fd3xd3v sionlessmodeltheinstabilitydisappears.Thephenomenonisat and variancewith whathappens inhomogeneous systems(e.g., see (cid:90) mv2 1 (cid:90) f(r,v)f(r(cid:48),v(cid:48)) Bertin2000).Inthelastpartofthepaperwearguethatthecause E= f(r,v)d3rd3v− Gm2 d3rd3vd3r(cid:48)d3v(cid:48), ofthedifferenceliesintheroleoftheangularmomentumofthe 2 2 |r−r(cid:48)| individualparticles. wheremisthesingle-particlemass. 2 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel adoptedconditions.Forexample,Ξ=34.4holdsforasystemat temwhent→∞(tistime)hasaneffectiveequationofstate3not constanttotalenergyE andvolumeV,Ξ = 8.99forasystemat affected by the attractive nature of the gravitational force: thus constanttemperatureT andvolumeV,andΞ=6.45forasystem the Boltzmann entropy may not be applicable to find the true atconstanttemperatureT andexternalpressureP.Inthispaper, thermalequilibriumofaself-gravitatingN-bodysystem,sinceit wewillshowthatthethreepointscanbeidentifiedbymeansof assumestheequationofstateofanidealgas.Moreover,because adynamicalstabilityanalysisofthefluidmodel. gravitational forces have long range, each particle feels the in- fluenceofallotherdistantparticles.Thus,itmaybethat,strictly Now we briefly comment on the use of the Boltzmann ex- speaking, in the case of self-gravitating systems an equation of pression for the entropy (1). If we refer to systems for which statecannotbedefinedintermsoflocalquantities. the traditional thermodynamic limit (N → ∞ and V → ∞ at AdifferentwaytojustifytheuseoftheBoltzmannentropy N/V fixed)iswelldefinedandleadstohomogeneousequilibria, is to show that it can be obtained from the Gibbs microcanoni- itiswellknown(Jaynes1965)thattheBoltzmannentropycor- calentropybymeansoftheso-calledmeanfieldapproximation respondstotheentropydefinedinphenomenologicalthermody- (Katz 2003; Padmanabhan 1990). However, the range of appli- namicsonlywhentheinterparticleforcesdonotaffectthether- cability of this approximation is still not clear, also because a modynamic properties; that is, the Boltzmann entropy neglects small-scale cut-off is necessary to avoid divergences in the mi- the interparticle potential energy and the effect of the interpar- crocanonicalentropy(Padmanabhan1990).TheBoltzmannen- ticle forces on the pressure. If the equation of state is different tropyisalsothe H quantityoftheBoltzmann H-theorem;from from that of an ideal gas, the Boltzmann entropy is in error by thispointofview,acriticalanalysisoftheBoltzmannentropyin anon-negligibleamount.Thecorrectexpressionfortheentropy thecontextofthegravothermalcatastrophehasbeenperformed (corresponding to phenomenological thermodynamics) of sys- byMiller(1973). tems with a well-defined thermodynamic limit is the one given The discussion in this section supports the hypothesis that byGibbs[forthedefinitionofGibbsentropy,seeJaynes(1965)]. self-gravitating isothermal spheres are not true thermal equi- Thegeneralizationoftheaboveresulttosystemswithinho- librium states of the pure N-body problem (i.e., the states that mogeneous equilibria, such as self-gravitating systems, is that would be found in a spherical box containing N gravitating the Boltzmann entropy is valid (i.e., it corresponds to entropy particles eventually, after an infinite amount of time), but are asdefinedinphenomenologicalthermodynamics)ifandonlyif only metastable states of which the significance is still not (1) the equation of state of the ideal gas applies locally to the completely understood (for example, see Padmanabhan 1990, equilibriumstates;(2)hydrostaticequilibriumholds.Indeed,the Chavanis2006andreferencestherein). mostelementarywaytounderstandtheBoltzmannentropyinthe case of self-gravitating systems is to think of a fluid in hydro- staticequilibrium,withtheequationofstateofanidealgas p= 3. Basicequationsofthedynamicalapproach ρkT/m, subject to reversible transformations between equilibriumstates(i.e.,betweendifferenttruncatedisothermalspheres Inthissectionwederivethelinearizedhydrodynamicequations that govern the evolution of a fluid system for small deviations inthesphericallysymmetriccase).AsshownbyLynden-Bell& fromthetruncatedisothermalsphereequilibriumconfigurations. Wood (1968) in their Appendix I, the classical thermodynamic entropy defined from the relation dS = dQ/T and calculated We assume spherical symmetry, that is, we consider only ra- c dial perturbations. The configurations are known to be stable from such transformations coincides with S . For the general b against nonradial perturbations (Semelin et al. 2001; Chavanis case of a system of particles interacting through an arbitrary 2002;Binney&Tremaine2008). two-bodypotential(notnecessarilygravitational),itispossible Consider a self-gravitating fluid governed by the Navier- toshowthatthestationarystatesoftheBoltzmannentropyhave Stokes and continuity equations, together with the equation of thesamedensitydistributionasthatofafluidwiththeequation stateofanidealgas(forexample,seeLandau&Lifshitz1987): ofstateofanidealgasinhydrostaticequilibrium. The above considerations suggest that the validity of the ∂ρ +∇·(ρu)=0, Boltzmannentropyisstrictlyrelatedtothevalidityoftheequa- ∂t tionofstateofanidealgas.Inakineticdescription,theequation ∂u ∇p η ζ+ η ofstateofanidealgasisobtainedbyconsideringthelocalone- +(u·∇)u=− −∇Φ+ ∇2u+ 3∇(∇·u) , (2) particle distribution function to be a Maxwellian (i. e., of the ∂t ρ ρ ρ form f(r,v) = Aexp[−mv2/2kT(r)],whereAisanormalization kT (cid:82) p=ρ , constant)andbydefiningthepressureasp≡ mf(v2/3)d3rd3v. m Thisdefinitionofpressure,whichisobtainedbyconsideringpar- ticlesthatwouldreversetheirmomentumwhenhittinganimag- whereρisthedensity,uisthefluidvelocity, pisthelocalpres- inary wall, ignores the effects of interparticle forces. If, for ex- sure,Φisthegravitationalpotential,ηandζareviscositycoeffi- ample, particles repelling one another were confined in a box, cients,mistheone-particlemass,T isthetemperature,tistime, theywouldexertsomepressureonthewallsoftheboxevenifat andkistheBoltzmannconstant.Wewillkeeptrackofviscosity rest:thiscontributioniscompletelyneglectedbytheequationof termuntiltheendofSubsection3.1;then(fromSubsection3.2 stateoftheidealgas(theneglectedpressureissimilartothatof on), for simplicity, we will set η = ζ = 0. In this paper, we do rigid spheres when packed too closely; thus the Boltzmann en- not perform an analysisof the effects of viscosity. Inparticular tropycannotgiveacorrectresultforagasofalmostrigidspheres we do not discuss whether it has a stabilizing or destabilizing whentheirdensityistoohigh). effect:hopefully,theroleofviscositywillbestudiedindetailin asubsequentpaper. In the case of particles interacting through gravity, the ne- glected contribution is attractive and should therefore decrease 3 Byeffectiveequationofstatewemeananequationofstatethat,by thepressurecomparedtothatofanidealgas.Therefore,itisun- imposing hydrostatic equilibrium, would reproduce the density distri- likelythatthetruethermalequilibriumstateofan N-bodysys- butionoftheequilibriumstate. 3 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel Under the assumption of spherical symmetry, the gravita- Inthefollowing,weshalldropthesymbol˜tokeepthenotation tionalpotentialobeysthefollowingequations: simpler. After some manipulations (see Appendix B.1) and using the unperturbed density profile (4) we obtain the following GM(r) ∇Φ= rˆ, linearizedequation: r2 (cid:90) r (3) M(r)= ρ(s)4πs2ds, iω kT 0 Lf = Lf + 1ψ(cid:48)e−ψ (8) 4πGλ m whicThhearehyedqruoivstaalteinctetoquthileibProiaissoofnseuqcuhatfliouni.d system are spa- + iω m (cid:40)η 1 [ξ2(feψ)(cid:48)](cid:48)+(ζ+ η)(cid:34) 1 (cid:16)ξ2feψ(cid:17)(cid:48)(cid:35)(cid:48)(cid:41), ρ (0)kT ξ2 3 ξ2 tially truncated isothermal spheres (Chandrasekhar 1967) and 0 0 are briefly described in Appendix A. As mentioned in Sect. 2, where f(ξ)=ρ (ξ)u (ξ)istheunknownfunction, eachself-gravitatingtruncatedisothermalsphereisidentifiedby 0 1 two dimensional scales (for example the central density ρ (0) 0 ω2 and the temperature T) and one dimensionless parameter Ξ = L= (9) 4πGρ (0) R/λ, which is the value of the dimensionless radius ξ ≡ r/λ 0 at the truncation radius (the scale λ was defined at the begin- representsthedimensionless(squared)eigenfrequencyandLis ningofSect.2).Conversely,eachchoiceofthethreequantities thefollowingdifferentialoperator: {ρ (0),T,Ξ}identifiesaparticularisothermalsphere. 0 We shall denote unperturbed quantities by subscript 0 and d2 (cid:32)2 (cid:33) d (cid:32) 2 2ψ(cid:48) (cid:33) the perturbations by subscript 1. Thus the unperturbed density L≡− − +ψ(cid:48) + − −e−ψ . (10) dξ2 ξ dξ ξ2 ξ profile (the density profile of a truncated isothermal sphere) is, asafunctionofthedimensionlessradiusξ(seeAppendixA): Equation (8) is the equation to be solved to find the nor- (cid:40)ρ (ξ)=ρ (0)e−ψ(ξ) ifξ ≤Ξ mal modes and frequencies of the system, under the appropri- 0 0 (4) ate boundary conditions and constraints. The boundary condi- 0 ifξ >Ξ, tions for Eq. (8) are defined in the following way. The absence of sinks and sources of mass, together with the assumption of where ρ (0) is the unperturbed central density and ψ(ξ) is the 0 sphericalsymmetry,implies f(0)=0.SinceintheEulerianrep- regular solution to the Emden equation (the symbol (cid:48) denotes resentationweshallconsideronlysystemsinasphericalboxof derivativewithrespecttotheargumentξ): constant volume, the radial velocity at the edge must be zero, d whichimplies f(Ξ)=0.Thustheboundaryconditionsare: (ξ2ψ(cid:48))=ξ2e−ψ, ψ(0)=ψ(cid:48)(0)=0. (5) dξ f(0)= f(Ξ)=0. (11) When we perturb the hydrodynamic equations it is convenient Tosolvetheequations,westillneedtospecifythefunction to work in the Eulerian or Lagrangian representation of hydro- T (t).Thisspecificationdiscriminatesbetweentheconstanten- dynamics,dependingonwhichboundaryconditionsweimpose. 1 ergyandtheconstanttemperaturecase.Oncethisspecificationis IntheEulerianrepresentationtheindependentvariablesarethe made,Eq.(8)isaneigenvalueproblem:forfixedω,theequation timetandthepositionvectorr.IntheLagrangianrepresentation admitsasolutiononlyfordiscretevaluesof L.Thesesolutions theindependentvariablesarethetimetandtheoriginalposition aretheradialnormalmodesofthesystem. r that the fluid element under consideration had at the initial 0 The operator L has the following properties, which can be timet . 0 proveddirectlyfromtheEmdenequation(5): 3.1. Eulerianrepresentation L(ψ(cid:48))=−e−ψψ(cid:48) . (12) L(ξe−ψ)=−e−ψψ(cid:48) In this subsection we derive the linearized perturbation equations in the Eulerian representation. So we write each quantity Thesepropertiesallowustoobtainanalyticalsolutionsinsome as the sum of an unperturbed part (characterizing the truncated cases. They are equivalent to those found by Padmanabhan isothermalsphere)andaperturbation: (1989) in his review of the thermodynamical analysis of ρ(r,t)=ρ (r)+ρ (r,t) Antonov (1962) and later also used by Chavanis (2002, 2003) u(r,t)=u0(r,t) 1 (6) toobtainanalyticalsolutionsofthepresenthydrodynamicprob- T(t)=T 1+T (t), lemfortheconstant{T,V}(Sect.4.1.1)andconstant{T,P}(Sect. 0 1 4.2.1)cases. whereuistheradialcomponentofthefluidvelocity(recallthat Notethatviscositydisappearswhenω = 0,whichisthesit- weconsideronlyradialperturbations).WesubstituteEq. (6)in uationofmarginalstability.Thusviscositydoesnotmodifythe Eq.(2)andexpandtofirstorderintheperturbedquantities.We pointsoftheonsetofinstability(Semelinetal.2001;Chavanis allow the temperature to vary in time while remaining uniform 2002). inspace:thisconstraintwillbeusedinordertoimposethecon- ditionofconstanttotalenergy(seeSubsection4.1.2).Tofindthe 3.2. Lagrangianrepresentation normalmodesofthesystem,welookforsolutionsofthefollow- ingform: Inthissectionwepresentthelinearizedperturbationequationsin ρ (r,t)=ρ˜ (r)e−iωt theLagrangianrepresentationintheinviscidcase.Inthisrepre- 1 1 u1(r,t)=u˜1(r)e−iωt (7) sentation,theindependentvariableisthepositionr0 ofthefluid T1(t)=T˜1e−iωt . elementunderconsiderationattheinitialtimet0.ThusinEq.(2) 4 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel we have to perform the following change of independent vari- 4.1.1. Constant{T,V}case(isothermalcollapse) ables: (cid:40) HereweconsiderafluidatconstanttemperatureT containedina r→r (r,t) 0 (13) sphereoffixedradiusR.Theconditionofconstanttemperature t→t. is satisfied by imposing T = 0 in Eq. (8). The condition of 1 constant volume has been discussed in Sect. 3.1 and leads to IntheLagrangianrepresentation,eachquantityisafunctionof theboundarycondition f(Ξ) = 0.ThusEq.(8)fortheconstant thenewindependentvariablesr andt. 0 {T,V}casebecomes: ThecalculationsaresummarizedinAppendixB.2.1.Forlin- Lf = Lf (18) earperturbations,theresultingcontinuityandEulerequationsin theLagrangianrepresentationare(neglectingthetermquadratic withtheboundaryconditions(11). inthevelocity): Equation(18)isaneigenvalueequationthat,atgivenΞ,admits solutions only for discrete values of L. If the lowest value of LatfixedΞispositive,thenallmodesarestable(becauseω   ∂∂t − rr22ρρ u∂∂r ρ+ r12ρρ ∂∂r (cid:16)r2ρu(cid:17)=0 iastaalgwivayensrveaalul)eaonfdΞthisenseygsatetimvei,sthsteanbulen.sItfabthleemloowdeesstavraelupereosefnLt  0 0 0 0 0 0 (14) andthesystemisunstable. ∂∂t − rr22ρρ u∂∂r u=−rr22∂∂rρ ρ1 kmT − GMr2(r0). Bymeansofthestandardtransformation 0 0 0 0 0 0 (cid:34) (cid:90) ξ(cid:32)1 ψ(cid:48)(ξ(cid:48))(cid:33) (cid:35) f(ξ)= f˜(ξ)exp − + dξ(cid:48) , (19) AsfortheEuleriancaseweseparateeachquantityinanun- ξ¯ ξ(cid:48) 2 perturbedpartandaperturbation: whereξ¯isanarbitrarypointinthedomainof f,Eq.(18)canbe ρ(r ,t)=ρ (r )+ρ (r ,t) recastintheformofaSchro¨dingerequation: 0 0 0 1 0 u(r ,t)=u (r ,t) r(r0,t)=r1+0r (r ,t) (15) − f˜(cid:48)(cid:48)(ξ)+ f˜(ξ)U(ξ)= Lf˜(ξ), (20) 0 0 1 0 T(t)=T0+T1(t). whereU(ξ)istheeffectivepotentialgivenby: We substitute Eq. (15) in Eqs. (14) and expand to first order in U(ξ)≡ 2 + 1ψ(cid:48)(ξ)2− 2ψ(cid:48)(ξ)− 1e−ψ(ξ). (21) quantities with subscript 1. Then to find the normal modes we ξ2 4 ξ 2 take: ρ (r ,t)=ρ˜ (r )e−iωt Theboundaryconditionsbecome: 1 0 1 0 u1(r0,t)=u˜1(r0)e−iωt (16) f˜(0)= f˜(Ξ)=0. (22) r (r ,t)=r˜ (r )e−iωt . 1 0 1 0 The effective potential is shown in Fig. 1. From the boundary In the following we shall drop the symbol ˜ for simplicity of conditions(22)weseethatchoosingaspecificΞrequiresusto notation.Afterintroducingthedimensionlessradiusξ0 ≡ r0/λ, consider a potential that is infinite for ξ ≥ Ξ. The existence of using the density profile (4) of the unperturbed state and after negative eigenvalues implies that the system is unstable. From somemanipulationsweobtainthefollowingequationforρ1(ξ0) theformofthepotentialitisclearthatforsmallvaluesofΞneg- (seeAppendixB.2.2foranoutlineofthecalculations): ative eigenvalues (i.e., unstable modes) do not exist. Negative − ρ0ρ(10) + eξ−02ψddξ0 (cid:32)ddξξξdd40003(e+−ρ0ψρ1((0ξξ)002))ψ+L(cid:48)(ξ0TT)01(cid:33)=0, (17) teisnhnaiagefimmenpeniicrtvpeeaasollneiuannuepttmsspprbwaroepoahrpbceeolhreaefm(rnsnoeeetenghwlaiKystuiavonftoecsztcrae1ubsi9rguls7eeffi8nam)vtc.aioΞeldune≥etslsy8aap.lp9app9rege.aeaWrr,vihnapelrtunheeceΞsisthoe→elfyrΞm∞a.ot,Fdtahyone-r LetusnowconsiderthesolutionstoEq.(18)ingreaterdetail. whereL=ω2/4πGρ (0)asfortheEuleriancase.Boundarycon- Figure2showstheminimumvalueofLatgivendimensionless 0 radius Ξ as a function of Ξ. We see that L is negative, that is, ditions are discussed in Subsection 4.2.1. Solving Eq. (17) al- the system is unstable, for Ξ > 8.99, which is the same point lowsustofindthenormalmodesintheLagrangianrepresenta- found by Lynden-Bell & Wood (1968) in the thermodynamical tion.WehaveallowedthetemperaturetovaryintimeinEq.(17), approach.Highermodes,thatis,highervaluesof L atgivenΞ, but in the following we shall analyze only the isothermal case withT =0. would be represented by lines above the plotted curve. These 1 lineswouldintersecttheΞaxisatsomepoints,whicharetheze- rosoftheanalyticalsolutionG describedbelow[seeEq.(23)]. TV 4. Modalstabilityofaself-gravitatinginviscidfluid In Appendix D.1 a few density and velocity profiles of numerically calculated eigenfunctions are shown. Most of them are sphereunderdifferentboundaryconditions formodesofminimum LatgivenΞ.Densityprofilesofhigher Hereweanalyzetheequationsobtainedintheprevioussection modesexhibitoscillationsnotpresentinthelowestmode. byimposingdifferentboundaryconditions. For the case of marginal stability (L = 0), with the help of properties(12),therelevanteigenfunctioncanbeexpressedan- alytically(Chavanis2002).Thefunction 4.1. Eulerianrepresentation G (ξ)≡ψ(cid:48)(ξ)−ξe−ψ(ξ) (23) TV InthissubsectionweanalyzeEq.(8)byimposingtwokindsof boundary conditions: constant temperature T and volume V or is indeed a solution to Eq. (18) with L = 0, which satisfies constanttotalenergyEandvolumeV. G (0)=0.ThevaluesofΞforwhichG (Ξ)=0arethosefor TV TV 5 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel 1.0 It has two terms, which represent the thermal and gravitational contributions. The condition of constant energy is imposed in 0.8 thefollowingway.Whenthefluidisperturbed,itsgravitational 0.6 energychangesasaconsequenceoftheredistributionofmatter. We suppose that the temperature varies in time, while remain- 0.4 ing uniform in space, so as to keep the total energy (thermal U 0.2 plusgravitational)constant(foradifferentmodel,basedonthe Smoluchowski-Poissonsystemofequations,thesamenonstan- 0.0 dard assumption has been made by Chavanis et al. 2002). The (cid:45)0.2 thermalenergyexpressionattimetisthusgivenby3NkT(t)/2. Indoingso,weareassuminginfinitethermalconductivity(see (cid:45)0.4 also Sect. 5.1 for the relation of this fact to the relevant time 0 10 20 30 40 50 scales). Ξ ToreduceEq.(8)totheconstant{E,V}caseweneedtofind Fig.1.TheeffectivepotentialU [seeEqs.(20)and(21)]forthe the expression for the temperature as a function of the density constant {T,V} case. To find the frequencies of normal modes, distribution at fixed total energy, in the linear regime of small the Schro¨dinger equation (20) has to be solved with boundary perturbations. Starting from Eq. (24) for the total energy and conditions(22),whichmeansthattheeffectivepotentialisU(ξ) recallingthatΦ(r)=−G(cid:82)d3r(cid:48)ρ(r(cid:48))/|r−r(cid:48)|,wehave: for ξ ≤ Ξ and taken to be infinite for ξ > Ξ. States with negative energies, which exist for Ξ ≥ 8.99, correspond to unstable 3 G (cid:90) ρ(r)ρ(r(cid:48)) E = NkT − d3rd3r(cid:48). (25) modes. 2 2 |r−r(cid:48)| By substituting T = T +T , ρ = ρ +ρ in Eq. (25), keeping 0 1 0 1 0.20 onlyfirst-orderquantitiesandimposingthattheenergyremains constant,weobtainthefollowingexpressionforT : 1 0.15 (cid:82) ρ (r)Φ (r)d3r T =− 1 0 . (26) 0.10 1 3Nk 2 HereΦ isthegravitationalpotentialoftheunperturbeddensity 0.05 0 distribution ρ . After some manipulations (see Appendix B.3) 0 weobtain: 0.00 (cid:82)Ξ -0-.000-5-- -0.042 T1 =−i1ω 0{d3[ξρ2(f0(ξ)λ)]Ξ/d2ξψ}(cid:48)ψ(Ξ(ξ))dξT0. (27) 2 0 0 20 40 60 80 By substituting Eq. (27) in Eq. (8), recalling the definition λ = Fig.2. The minimum value of the eigenvalue L at given Ξ (the (cid:2)kT0/4πGρ0(0)m(cid:3)1/2,andneglectingviscosity,weobtain: dΞi,mfeonrstihoenlceossnsratadniuts{Tch,Var}acctaesreiz.inIfgtthheemsyisntiemmu)masLafisunncetgioantivoef Lf =Lf −ψ(cid:48)e−ψ 1 (cid:90) Ξψ(ξ) d (cid:104)ξ2f(ξ)(cid:105) dξ. (28) then the system is unstable. The system becomes unstable at 3Ξ2ψ(cid:48)(Ξ) dξ Ξ = 8.99, the same value obtained from the thermodynamical 2 0 approach.Fromthefigurewecanreadthetypicaltimescaleof AsdiscussedinSect. 3.1theboundaryconditionsaregivenby theinstability.AsΞ→∞,theasymptoticvalueL =−0.042is Eq.(11).ByintegratingthelasttermofEq.(28)bypartsunder ∞ thesameasfortheothercases(seeFigs.3and4). theboundaryconditions(11),weobtain: (cid:90) Ξ 1 Lf =Lf +ψ(cid:48)e−ψ ξ2ψ(cid:48)(ξ)f(ξ)dξ. (29) which the boundary conditions (11) are satisfied; thus they are 3Ξ2ψ(cid:48)(Ξ) thevaluesofΞatwhicheachnewunstablemodeappears.The 2 0 first zero ofG occurs where the first unstable mode appears, Note that Eq. (29) [or (28)] contains an integral global con- TV atΞ=8.99.Fromtheasymptoticbehaviorofψitispossibileto straint. Equation (29) is to be solved to find the normal modes obtainanasymptoticapproximationofthezeros:theyapproxi- andcorrespondstoEq.(18).Thenumericalprocedurefollowed √ mately follow a geometric progression of ratio e2π/ 7 (see also tosolveEq.(29)isrelativelystraightforwardandthusisnotre- Semelinetal.2001;Chavanis2002). portedhere. Figure3showstheminimumvalueof LatgivenΞ.Wesee that the system is unstable for Ξ > 34.36, which corresponds 4.1.2. Constant{E,V}case(gravothermalcatastrophe) to a density contrast ρ (0)/ρ (Ξ) > 709 (Antonov 1962). As Ξ 0 0 increases,newunstablemodesappear,similarlytothebehavior Here we consider a self-gravitating fluid sphere at constant to- observedintheconstant{T,V}case.AsΞ→∞,theasymptotic tal energy E and volume V. In this case the instability has valueoftheminimumLisfoundnumericallytobethesameas beennamedgravothermalcatastrophebyLynden-Bell&Wood intheconstant{T,V}case,L (cid:39)−0.042.Suchasymptoticvalue (1968).ThetotalenergyEofthefluidisdefinedas: ∞ isthesamealsointheconstant{T,P}case(seeSect.4.2.1). 3 1(cid:90) In the case of marginal stability (L = 0), Eq. (29) is equiv- E = NkT + ρ(r)Φ(r)d3r. (24) alent to that found and solved analytically by Padmanabhan 2 2 6 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel 0.5 instability: by numerically solving the algebraic equation (36), this minimum value is found to be Ξ = 34.36 (as shown by 0.4 Antonov 1962). As for the constant {T,V} case, the number of unstablemodesforeachΞcoincideswiththeresultsofthether- 0.3 modynamical analysis [see Katz (1978)]. Once a value of Ξ is (cid:76)0 obtained,thevalueofa/bandthenananalyticalsolutionG is min(cid:72)Ρ0 0.2 determined. EV 2ΩΠG Density profiles for modes of minimum L at given Ξ are 4 0.1 showninAppendixD.2.Sincetheequationsinvolvedareequiv- alent, the density profile for the marginally stable perturbation 0.0 (cid:88)=34.363 (L = 0) is the same as that found by Padmanabhan (1989). He pointed out that it has a “core-halo” structure, that is, an oscil- -0.042 lationinρ :thedensityperturbationispositiveintheinnerpart (cid:45)0.1 1 0 20 40 60 80 100 (nucleus),negativeinthemiddle(emptyingarea),andthenposi- (cid:88) tiveagain(halo).Thecore-halostructurehasbeenphysicallyin- Fig.3. Theminimumvalueoftheeigenvalue LatgivenΞ,the terpretedintheframeworkofthegravothermalcatastrophegiven dimensionlessradiuscharacterizingthesystem,asafunctionof by Lynden-Bell & Wood (1968), using the concept of negative Ξ,fortheconstant{E,V}case.WhentheminimumLisnegative specificheat.However,thisinterpretationisnotapplicabletothe the system is unstable. The system becomes unstable at Ξ = presentcontext. 34.36,whichcorrespondstoadensitycontrast(cid:39) 709,thesame WefoundthatformodesofminimumLatgivenΞthecore- value obtained from the thermodynamical approach. From the halo structure is present if L < 0.021, disappears between L = figure we can read the typical time scale of the instability. As 0.021andL = 0.022,andisabsentforL > 0.022,asillustrated Ξ→∞,theasymptoticvalueofL,L∞ =−0.042,isthesameas inthefigurespresentedinAppendixD.2.Highermodesalways fortheothercases(seeFigs.2and4). exhibitoneoremoreoscillations,asintheconstant{T,V}case. Wehavethusshownthatforthepresentcasethebehaviorof theeigenvaluesLasafunctionofΞissimilartothatofthecon- (1989) in the thermodynamical approach. Here, for complete- stant{T,V}case,withthedifferencethattheinstabilitythreshold ness, we record how the solution is derived, by adapting the is higher in the constant {E,V} case. The interpretation of this method of Padmanabhan (1989) to our choice of variables and factinthecontextofthefluidmodelisasfollows.Whenthefluid unknowns. is compressed, the gravitational potential energy decreases and TosolveanalyticallyEq.(29)for L = 0werewriteitinthe the temperature increases in order to maintain the total energy followingform: Lf =−Cψ(cid:48)e−ψ, (30) constant. Therefore, the tendency toward collapse is weakened and instability can take place only at higher values of Ξ (with where respecttothecaseinwhichthetemperatureremainsconstant). (cid:90) Ξ 1 Infact,ifinsteadofexpression(24)forthetotalenergywe C = ξ2ψ(cid:48)(ξ)f(ξ)dξ (31) 3Ξ2ψ(cid:48)(Ξ) consider a modified expression in such a way that the temper- 2 0 ature increase is greater, the collapse can be halted completely. isaconstant(withrespecttoξ)thatdependsgloballyon f.Using Considerthefollowingheuristicexpressionforthetotalenergy: theproperties(12),welookforasolutionoftheform (cid:90) 3 1 E = NkT + ρ(r)Φ(r)d3r−υσ2R3ρ(0), (37) GEV(ξ)=aψ(cid:48)(ξ)+bξe−ψ(ξ), (32) 2 2 0 where a and b are real numbers. Since Eq. (30) is linear in f, whichcontainsanadditionaltermproportionaltothedimension- onlytheratioa/bisrelevant.Susbtituting(32)in(30)weobtain less parameter υ and to the central density ρ(0). The quantity thefirstconditiononaandb: σ20 = kT0/mistheunperturbedthermalspeed.Werepeatedthe analysis of this subsection with such a modified expression for a+b=C. (33) theenergyandfoundthenewvalueofΞforthethresholdofthe instability. In this analysis the expression for T (27) obtained 1 UsingtheexpressionofC(31)anddividingbyb,weobtain: previously,tobesubstitutedinEq.(8),isreplacedwiththeex- a +1= 1 (cid:90) Ξξ2ψ(cid:48)(ξ)(cid:18)aψ(cid:48)(ξ)+ξe−ψ(ξ)(cid:19) dξ. (34) pTre=ssTio0n+Tfo1r,ρT1=tρh0a+t ρis1oanbdtabinyeidmfproosmingEqth.a(t3t7h)ebvyarisautbiosntitouftitnhge b 32Ξ2ψ(cid:48)(Ξ) 0 b totalenergyE iszero.Wefoundthatwhenυ>0(i.e.,whenthe fluid is compressed the new term contributes to make the tem- Asecondconditionona/bfollowsfromtheboundarycondition perature increase) the(linear) instability is postponed to higher f(Ξ)=0: valuesofΞforsmallυandiscompletelyhaltedforυ>5×10−4. a ψ(cid:48)(Ξ)+Ξe−ψ(Ξ) =0. (35) Wealsoverifiedthatifυ<0(i.e.,whenthefluidiscompressed b the new term contributes in the opposite direction), instability Substitutinga/bfrom(35)in(34)weobtain: occursatlowervaluesofΞ. In the case of the pure N-body problem, it is believed that 3Ξ2(cid:16)ψ(cid:48)(Ξ)−Ξe−ψ(Ξ)(cid:17)=(cid:90) Ξξ2ψ(cid:48)(ξ)(cid:32)ξe−ψ(ξ)−ψ(cid:48)(ξ)Ξe−ψ(Ξ)(cid:33) dξ collapse can be halted by energy “generation” through bina- 2 ψ(cid:48)(Ξ) ries, giving rise to a phenomenon called gravothermal oscil- 0 (36) lations (for a review, see Heggie & Hut 2003). Gravothermal The values of Ξ satisfying Eq. (36) are those for which oscillations were discovered by Sugimoto & Bettwieser (1983) Eq. (29) admits a solution for L = 0. In particular, the lowest using a gaseous model. These authors introduced in the model value of Ξ which satisfies Eq. (36) determines the threshold of of Lynden-Bell & Eggleton (1980) a phenomenological energy 7 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel generation term (to represent the role of binaries), proportional 0.5 toapowerofthedensityandfoundthatcollapsecanbehalted and reversed. Similarly, our υ term shows in a simple manner 0.4 thattheinstabilitycanbehaltedinthelinearregimewithasuit- ableterminthetotalenergybudgetthatmimicseffectspresumed 0.3 tobeassociatedwiththepresenceofbinaries. (cid:76)0 Gravothermal oscillations have been confirmed by N-body 2Ωmin(cid:72)ΡG0 0.2 simulations (Makino 1996). We will comment further on this Π4 topicinSect.5.2. 0.1 (cid:88)=6.45 0.0 4.1.3. Constant{T,V}:two-componentcase -0.042 (cid:45)0.1 Webrieflystudiedtheproblemofdynamicalstabilitybymeans 0 5 10 15 20 25 30 ofalinearmodalanalysisofatwo-componentidealfluid.Each (cid:88) component is assumed to interact with the other only through Fig.4. Theminimumvalueoftheeigenvalue LatgivenΞ,the thecommongravitationalpotential.Thehydrostaticequilibrium dimensionlessradiuscharacterizingthesystem,asafunctionof configurationsarespatiallytruncatedtwo-componentisothermal Ξ, for the constant {E,V} case. When the minimum L is neg- spheres, as considered by Taff et al. (1975), Lightman (1977), ative the system is unstable. The system becomes unstable at andYoshizawaetal.(1978);seealsodeVega&Siebert(2002) Ξ = 6.45,thesamevalueasfoundinthethermodynamicalap- andSopiketal.(2005). proach.Fromthefigurewecanreadthetypicaltimescaleofthe Wecallthesingle-particlemassesm andm ,withm >m . instability.AsΞ→∞,theasymptoticvalueL =−0.042isthe A B B A ∞ Two-component isothermal spheres are characterized by two sameasfortheothercases(seeFigs.2and3). additional dimensionless parameters (with respect to the one- componentcase):theratioofthesingle-particlemassesm /m B A andtheratioM /M ofthetotalmassesassociatedwiththetwo B A Boundary conditions are as follows. The condition u (0) = 0 components.ThereaderisreferredtoAppendixEforadescrip- 1 translates into the condition ρ(cid:48)(0) = 0 [using Euler equation tionoftheequationsused. 1 (14), under the assumption that ∂u/∂r does not diverge and We considered the constant {T,V} case, which is the sim- 0 ∂u/∂t =−iωu].Theconditionofconstantpressurerequiresthat plest from the mathematical point of view, and found that the afixedLagrangianfluidshell(whichfollowsthefluidduringthe equivalenceofthethermodynamicalanddynamicalapproaches motion and thus does not have a fixed position in space) feels stillholds.Itispossibletoshowanalyticallythattheonsetofdy- constant pressure. Constant pressure requires constant density namicalinstabilitytakesplaceatexactlythesamepointsasthose because of the ideal gas equation of state. Therefore, the rele- found with the thermodynamical approach. The analysis in the vantboundaryconditionsare: thermodynamicalapproachwasperformedbygeneralizingina straightforwardmannertheanalysisofChavanis(2002). ρ(cid:48)(0)=0 Aninterestingresultofthisanalysisisthattheinstabilityap- ρ1(Ξ)=0. (39) pears to be driven by the heavier component, in the following 1 sense.Considertheratiosρ /ρ andρ /ρ ofthedensityper- turbationofeachcomponen1tAto0thetota1lBloc0alunperturbedden- Frommasscontinuityandbyimposingu(cid:48)(0) (cid:44) 0(whichistrue 1 sity.Theratioreferringtotheheaviercomponentcanbehigher in the Eulerian representation and thus we expect to be true in even if the total mass of the heavier component is very small. thepresentcase),wealsohavetheconditionρ1(0)(cid:44)0. For example, for m /m = 3, we found that the ratios ρ /ρ , The numerical procedure to solve Eq. (38) is relatively B A 1A 0 ρ /ρ were comparable for M /M (cid:39) 0.066 (see Fig. E.1). straightforward and is thus not reported here. Figure 4 repre- 1B 0 B A Forhighervaluesof M /M theheaviercomponentdominates. sentstheminimumvalueof LatgivenΞ.Thesystembecomes B A Independent indications that the heavier component dominates unstableforΞ > 6.45,whichisthesameconditionasfoundby thecollapsehavebeenfoundbySopiketal.(2005)foradiffer- Bonnor(1956),Ebert(1955),andLynden-Bell&Wood(1968). entmodelbasedontheSmoluchowski-Poissonsystemofequa- As Ξ increases, other unstable modes appear, similarly to the tions. Breen & Heggie (2012a,b) found that the heavier com- casesdescribedpreviously.AsΞ → ∞,theasymptoticvalueof ponent dominates gravothermal oscillations in two-component the minimum L is found numerically to be the same as for the clusters.Thisislikelytoberelatedtothepresentanalysis. constant {E,V} and {T,V} cases, L∞ = −0.042, but, in contrast totheconstant{E,V}and{T,V}cases,itisreachedfrombelow. ThereisaminimumatL≈−0.043forΞ≈15. 4.2. Lagrangianrepresentation In Appendix D.3 density profiles of normal modes are shown. Most of them are for modes of minimum L at given 4.2.1. Constant{T,P}case(isobariccollapse) Ξ. Higher modes present oscillations. It would be interest- Inthissectionweanalyzeperturbationsatconstanttemperature ing to show whether the analysis of this subsection might be T and constant boundary pressure P. Therefore, T = 0 and generalized to the nonspherically-symmetric case by following 1 Eq.(17)becomes: Lombardi&Bertin(2001). Asforpreviouscases,wecanobtainanalyticalsolutionsfor  d ( ρ1 )  themarginallystableperturbations.ForL=0,Eq.(38)reads: − ρ0ρ(10) + e−ξψ02(ξ0)ddξ0 ξ403dd+ξξd00eξ−ρ020ψψ((0ξL(cid:48))0()ξ0)=0. (38) L(cid:32)ρ0ρ(10)(cid:33)=0, (40) 8 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel whereLisdefinedas: 5.2. Thedynamicsofacollisionlessself-gravitating  d  isothermalsphere L≡−1+ eξ−02ψddξ0  dξd0dξ4ξ0e30−ψ=0. (41) IigdsroeatavhliietzaretmidnagmliossdpoehthlesereormsftaahlreesppehuqereureiNlsib-cbraionudmbyepcsrtoounbdfilieegmdur.aaIstnisoptnaatsritoifncouarrlaysre,svtsaeetrleafs-l ofthecollisionlessBoltzmannequation: FromEmdenEq. (5),theoperatorLhasthefollowingproper- ties: ∂f(r,v,t) ∂f(r,v,t) ∂Φ(r,t) ∂f(r,v,t) L(ψ(cid:48)2)=−e−ψ +v· − · =0. (47) 2 (42) ∂t ∂r ∂r ∂v L(e−ψ)=−e−ψ . 4 Therefore,itisnaturaltoaskwhethersuchcollisionlessisother- Hence, an analytical solution of Eq. (40) is the following (for mal spheres are or can be unstable. It can be shown (Binney a derivation in the Eulerian representation, see also Chavanis & Tremaine 2008) that the unbounded collisionless isothermal 2003): G (ξ )=2e−ψ−ψ(cid:48)2. (43) sphere is linearly stable, in contrast to the results of the fluid TP 0 counterpart (in the fluid model we recover the unbounded case This solution satisfies the boundary conditions ρ1(0) = 2 and bytakingΞ → ∞).Itisgenerallybelievedthatthecollisionless ρ(cid:48)1(0) = 0, as required by Eq. (39). The zeros ofGTP(ξ0) allow isothermalsphereisalsostableinthenonlinearregime,although ustoidentifythemarginallystablenormalmodesthatsatisfythe toourknowledgearigorousproofofthisstatementisstilllack- correctboundaryconditions(39).ThefirstzeroofGTP(ξ0)isat ing. Moreover, if only spherically symmetric perturbations are Ξ = 6.45. Other zeros correspond to the values of Ξ at which considered,itispossibletoshow,throughanargumentbasedon new unstable modes appear and are found to be the same as in theconservationofthedetailedangularmomentum4 (Appendix thestandardthermodynamicalapproach. C), that a collisionless sphere (bounded or unbounded) cannot collapse,becauseeachparticlehasaminimumradiusitcanat- tain.Incontrast,thefluidsystemanalyzedinSects.3and4isun- 5. Thedifferentbehaviorofacollisionless stableonlywithrespecttosphericallysymmetricperturbations, self-gravitatingsphere with the subsequent nonlinear evolution presumably leading to 5.1. Timescales acollapse.Theaboveconsiderationssupportthehypothesisthat self-gravitatingcollisionlessisothermalspheres(boundedorun- Inthissubsectionwebrieflydiscussthetypicaltimescalesthat bounded)arestablewithrespecttoallkindsofperturbationsand characterize gaseous and fluid models and compare them to cannot collapse (Kandrup & Sygnet 1985; Kandrup 1990; Batt thoseofglobularclusters.Toalargeextent,ourdiscussionfol- etal.1995). lowsthatofInagaki(1980). The different behavior, between the collisionless and the Letτ bethelocalrelaxationtime,τ theglobalrelax- LTE GTE fluidcase,emergedinthispaperisactuallyquitesurprising,be- ationtime,andτ thedynamicaltimescale.Foragaseoussys- d causeforthehomogeneouscasethecollisionlessandfluidmod- tem,weidentifyτ withthetime(oftheorderoftheinverse LTE elsbehaveinthesamewaywithrespecttoJeansinstability(e.g., mean collision frequency) needed to reach a local Maxwellian see Bertin 2000), in the sense that both the fluid and the colli- distribution, τGTE with the time in which thermal conduction sionless systems are linearly unstable under the same criterion balances the temperatures of different parts of the system, and forinstability. τ withthesoundtraveltime.Gaseousmodelsgenerallyassume d If the self-gravitating collisionless isothermal sphere were thefollowingordering: unstable, its instability would develop on the dynamical time τ (cid:28)τ (cid:28)τ . (44) scale; in turn, it is commonly believed that the gravothermal LTE d GTE catastrophemayoccuronlyifthesystemisatleastweaklycol- For a globular cluster, because the mean free path is very large lisional and thus it is thought to develop on the collision time andstarscancrossthesystemmanytimesbeforeactuallycollid- scale. ing,theorderingoftimescalesisdifferent.Starsdonotcollide In the spherically symmetric case a difference between the significantly with neighboring stars, following the mechanisms twomodelsisthefollowing:whileinthefluidmodeleachfluid thatusuallycharacterizethermalconduction;instead,theytend element is sustained against gravity by pressure of the inner toreleasetheirenergythroughtheentirecluster(alsobecauseof parts,inthecollisionlessmodelstarsaresustainedbytheirindi- thelongrangenatureoftheforce).Thusforglobularclusterswe vidualangularmomentumrelativetothecenter,thatis,bytheir have: velocity dispersion.5 In moving from a kinetic description to a τd (cid:28)τLTE (cid:39)τGTE. (45) fluid description, all the information about microscopic veloci- In the fluid model analyzed in this paper we assumed infinite ties and detailed angular momentum is lost: in particular, each thermalconductivity,becausethetemperaturewasalwaystaken fluid element has zero angular momentum (see also Appendix to be and to remain uniform. Moreover, we implicitly assumed C).Inthisrespect,therealsituationofaglobularclusterresem- that the distribution function is locally Maxwellian. Therefore, bles more the collisionless case: strictly speaking, stars are not ourfluidmodelfollowstheordering: sustainedbypressure,butratherbytheirvelocitydispersion,because the mean free paths are long and stars cross the cluster τLTE (cid:28)τd manytimesbeforefeelingtheeffectsofcollisions. (46) τ (cid:28)τ . GTE d 4 Bydetailedangularmomentumwemeantheangularmomentumof Note that assumptions (44) and (46) are not applicable to the theindividualparticles. situation of a globular cluster (45). This difference and the re- 5 Evenifthetotalangularmomentumvanishes,thesumofthemag- latedlimitationsmustbekeptinmindifwewishtoapplythese nitudesoftheangularmomentaoftheindividualparticlesisdifferent modelstounderstandtheevolutionofglobularclusters. fromzero. 9 M.C.SormaniandG.Bertin:GravothermalCatastrophe:thedynamicalstabilityofafluidmodel If we come back to the original pure N-body problem, it is stablemodes.Indeed,asnotedintheIntroduction,someresults thereforenaturaltoask:whatisinthiscasetheroleofdetailed alongtheselineshavebeenobtainedpreviouslybyotherauthors. angularmomentum?Wemayarguethataquantityrelatedtothe Themainnewresultsobtainedinthispaperarethefollowing: detailed angular momentum should characterize the process of core collapse. If the collapse can happen in an N-body system, – UsingthefluidmodelbasedontheEulerequation,wehave it may be accompanied by a slow decrease of the sum of the extended previous studies of the constant {T,V} and {T,P} magnitudes of the angular momenta of the individual particles casestotheconstant{E,V}case(gravothermalcatastrophe), slowlysinkingtowardthecenter. proving that the onset of Jeans instability occurs exactly at For a particle with angular momentum J that moves in a the same point identified in the thermodynamical approach gravitational potential generated by a mass M a radius related also in this case. For this constant {E,V} case, we have in- toangularmomentum(whichistheradiusofthecircularorbitif troduced a heuristic term to incorporate effects akin to the Misapointmass)canbedefinedas: stabilizingroleofbinaries. – For all the three cases described above, we have calculated J2 numericallyeigenfrequenciesandeigenfunctionsoftherele- d ≡ . (48) vantmodesalsooutsidetheconditionsofmarginalstability. Gm2M The time scale for the instability that we have found is the Therefore, for a stellar system of N stars and total mass M we dynamical time scale. The excitation of higher modes has candefinearadiusrelatedtodetailedangularmomentuminthe beenillustratedinasimpleway,byreferringtoaneffective followingway: potentialthatgovernsthestructureofthelinearmodalanal- A ysis. R ≡ , (49) am NGm2M – We have found that for all the cases treated in our investi- whereAisthesumofthesquaresoftheangularmomentaofthe gation,asthedimensionlessradiusoftheisothermalsphere individualstars,thatis,thefollowingquantity: Ξbecomeslargerandlarger,thevalueofthedimensionless growth rate of the most unstable mode tends to a univer- (cid:88)N sal asymptotic constant value, independent of the adopted A≡ J2. (50) boundaryconditions. i i=1 – We have briefly shown that the correspondence between the stability in the dynamical and in the thermodynamical Thus we may argue that the typical time scale for core col- approach also holds for the two-component case and have lapseshouldcorrelatewith: foundindicationsthattheheaviercomponentisthemoreim- A portantdriveroftheinstability. tcc ≡ dA/dt . (51) – As a general discussion, we have commented on the meaning and applicability of the Boltzmann entropy for self- The main mechanism through which A varies with time is ex- gravitating systems and argued that the main difference be- pected to be that of two-body collisions. Thus t should be of tween the dynamical behavior of a fluid and a collisionless cc theorderofthetwo-bodyrelaxationtime. sphere,inrelationtotheirapplicationasmodelsofthepure By monitoring the quantity A(t) in N-body simulations, it N-body problem or a real weakly collisional stellar system would be interesting to test the role of detailed angular mo- suchasglobularclusters,istobeascribedtotheroleofthe mentuminthemechanismofgravothermaloscillations(Makino detailedangularmomentumbehaviorinthecollisionlessand 1996) [A(t) may reach an equilibrium value, around which the weaklycollisionalcases. systemgravothermallyoscillates;thiswouldbeconsistentwith The role of the viscosity and its consequences outside the con- the fact that the typical time scale of gravothermal oscillations dition of marginal stability have not been examined; hopefully, is the two-body relaxation time], its relevance to the studies of thisissuewillbeaddressedinafuturepaper. core-collapse in the gaseous model (Lynden-Bell & Eggleton Another interesting question is how the instability depends 1980; Sugimoto & Bettwieser 1983), and its connection with on the particle-particle interaction, for non-Newtonian cases thephenomenonofthegyro-gravothermalcatastrophe(Hachisu (Padmanabhan1989).Potentialsthatexhibitasofteningatsmall 1979). radii,suchas1/(r2+r2)1/2,wherer isaconstant(Chavanis& 0 0 Ispolatov2002,seealsoCasetti&Nardini2012),orthatdecline 6. Discussionandconclusions withadifferentpowerlawatlargeradii,suchas1/rα,withα(cid:44)1 (Ispolatov&Cohen2001),haveindeedbeenconsidered.Based Inthispaperwehavestudiedsystematicallythedynamicalsta- on the present article, we may argue that the results obtained bilityofaself-gravitatingisothermalfluidspherebymeansofa fromthethermodynamicalapproachwouldbereinterpretedand linearmodalanalysiswithrespecttosphericallysymmetricper- clarified as the analogue of the Jeans instability, by studying a turbations. In this sense, we have studied the Jeans instability fluid model for the potential considered, with the equation of intheinhomogeneouscontextofasphereoffinitesize.Within stateofaperfectgas. a unified framework, by imposing the boundary conditions of In general, this paper strengthens the view that the applica- constant{T,V}(isothermalcollapse),constant{E,V}(gravother- bility of different idealized models to describe the process of mal catastrophe), and constant {T,P} (isobaric collapse), we corecollapseinsystemsmadeofafinitenumberofstarsismore have proved that the onset of dynamical instability occurs ex- subtlethancommonlyreportedandthat,ingeneral,thestudyof actly at the points identified in the thermodynamical approach Jeansinstabilityofinhomogeneousstellarsystemsstillleavesa (see Antonov 1962; Lynden-Bell & Wood 1968; Bonnor 1956; numberofquestionsopen. Ebert 1955) and, by adapting derivations from other authors (Padmanabhan 1990; Chavanis 2002, 2003), we have provided Acknowledgements. We would like to thank Marco Lombardi, Francesco an analytic expression for the eigenfunctions of the marginally PegoraroandStevenN.Shoreformanyinterestingcommentsanddiscussions. 10