Thermodynamics of the self-gravitating ring model Takayuki Tatekawa1,2, Freddy Bouchet1, Thierry Dauxois1, Stefano Ruffo1,3 1. Laboratoire de Physique, UMR-CNRS 5672, ENS Lyon, 46 All´ee d’Italie, 69364 Lyon c´edex 07, France 2. Department of Physics, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo, 169-8555, Japan 3. Dipartimento di Energetica, “S. Stecco” and CSDC, 5 Universit`a di Firenze, and INFN, via S. Marta, 3, 50139 Firenze, Italy 0 (Dated: 2nd February2008) 0 2 We present the phase diagram, in both the microcanonical and the canonical ensemble, of the n Self-Gravitating-Ring (SGR) model, which describes the motion of equal point masses constrained a on a ring and subject to 3D gravitational attraction. If the interaction is regularized at short J distances by the introduction of a softening parameter, a global entropy maximum always exists, 5 and thermodynamics is well defined in themean-field limit. However, ensembles are not equivalent 2 andaphaseofnegative specific heatinthemicrocanonical ensembleappearsinawideintermediate energyregion,ifthesofteningparameterissmallenough. Thephasetransitionchangesfromsecond ] to first order at a tricritical point, whose location is not the same in the two ensembles. All these h features make of the SGR model the best prototype of a self-gravitating system in one dimension. c Inordertoobtainthestablestationarymassdistribution,weapplyanewiterativemethod,inspired e by a previous one used in 2D turbulence, which ensures entropy increase and, hence, convergence m towards an equilibrium state. - t a t I. INTRODUCTION s . t a There are many objects in our universe whose behavior can be understood considering only the gravitational m interaction. Examples are globular clusters, galaxies, clusters of galaxies, molecular clouds [1]. Different theoretical - approaches have been proposed to explain the peculiar statistical properties of self-gravitating systems. The main d n difficulty is thatthese systemscannotapproachstatisticalequilibriumbecause ofthe short-distancedivergenceofthe o potential and of the evaporation at the boundaries. Even if one puts the system in a box with adiabatic walls, thus c eliminating evaporation, still gravity causes the well-known phenomenon of gravothermal catastrophe [2, 3, 4]. The [ introduction of a small-scale softening of the interaction potential avoids such a catastrophe, so that self-gravitating 1 systems can approach the final (thermal) equilibrium state. However, such a state may have a negative specific heat. v Moreover,afirst-orderphasetransitionfromthe highenergygasphase to the lowenergyclusteredphase appears[3]. 3 Direct studies of the full three-dimensional N-body gravitational dynamics are particularly heavy [5] and even 8 special purpose machines have been built to this aim [6]. Therefore, lower dimensional models have been introduced 5 to describe gravitational systems with additional symmetries. For instance, the gravitational sheet model describing 1 the motion of infinite planar mass distributions perpendicularly to their surface has been considered [7]. Although 0 this model shows interesting behaviors [8, 9], the specific heat is always positive and no phase transition is present. 5 0 Recently, another one-dimensionalmodel has been introduced [10] where particle motion is confined on a ring, but / the interaction is the true Newtonian 3D one. At short distances, the potential is regularized, so that the particles t a do not interact. This model has been called the Self-Gravitating Ring model (SGR) and will be the subject of the m studydiscussedinthispaper. Ithasbeenshowninnumericalsimulations[10],thatthis modelmaintainsthe peculiar - features ofthe 3DNewtonianpotential, showinganegative specific heatphase anda phase transitionif the softening d parameter is small enough. Moreover, for large softening, this model reduces to the Hamiltonian Mean-Field model n (HMF) [11], which has been recently extensively studied as a prototype system with long-range interactions. This o latter model, however, although it displays a second order phase transition, does not have a negative specific heat c : phase at equilibrium. v In this paper, we derive the equilibrium thermodynamics of the SGR model both in the canonical and in the i X microcanonical ensemble. For all non vanishing softening parameter values, this model has a thermal equilibrium r state. If the softening parameter is small enough, the model shows ensemble inequivalence [12, 13] with a negative a specific heat phase in the microcanonical ensemble and a first order phase transition. Therefore, the SGR model displays severalfeatures of the true 3D Newtonian interaction,and canserve as a better prototype ofself-gravitating systems in one dimension than all previously introduced models. The paper is organized as follows. In Sec. II, we briefly introduce the SGR model and we discuss the essential features of previous numerical simulations [10]. In Sec. III, we show the general scheme for deriving all stationary density distributions which maximize Boltzmann-Gibbs entropy at fixed total energy and mass. Section IV presents a new iterative method which ensures entropy increase and leads in a unique way towards the stable equilibrium single particle distribution function. The method is inspiredby a similar one used to compute entropymaxima in 2D turbulence [14]. In Sec. V, we describe in full detail how to implement the iterative algorithmin a numerical scheme. 2 Figure 1: Self-Gravitating Ring model with a fixed unitary radius. Particles are constrained to move on a ring and therefore theirlocation isspecifiedbytheangles measuredwith respect toafixeddirection. Each pairofparticlesatθi andθj interacts throughtheinverse-squarethree-dimensionalgravitationalforce. Thedistanceismeasuredbythechord,asshowninthefigure. In Sec. VI, we calculate the thermodynamic quantities of the SGR model using the iterative method. We also show how, reducing the softening parameter, one enters into a region of ensemble inequivalence, where a tricritical point exists which is not the same in the twoensembles [15]. Finally, in Sec. VII, we discuss the dynamicalevolutionof the SGR model, emphasizing the properties of relaxation to equilibrium. II. SELF-GRAVITATING RING MODEL In this section, we briefly present the Self-Gravitating Ring (SGR) model [10]. In this model, particle motion is constrained on a ring and particles interact via a true 3D Newtonian potential (Fig. 1). The Hamiltonian of the SGR model is N 1 1 H = p2+ V (θ θ ), (1) 2 i 2N ε i− j i=1 i,j X X 1 1 V (θ θ ) = , (2) ε i j − −√2 1 cos(θ θ )+ε i j − − where ε is the softening parameter, which is introduced,pas usual, in order to avoidthe divergence of the potential at short distances. Taking the large ε limit, the potential becomes 1 1 cos(θ θ ) V = − i− j 1 +O(ε−2), (3) ε √2ε 2ε − (cid:20) (cid:21) which is the one of the Hamiltonian Mean-Field (HMF) model [11]. It is well known that the HMF model [11] has a second order phase transition, separating a low energy phase, where the particles form a single cluster, from a high energy gas phase where kinetic energy dominates and the particles are homogeneously distributed on the circle. One usually draws the so-called caloric curve, where temperature, given by twice the averagedkinetic energy per particle T β−1 = 2 K /N, is plotted against the total energy per particle U H/N. In a situation close to that of the ≡ h i ≡ HMF model, e.g. for ε = 10, the caloric curve determined from microcanonical numerical simulations is reported in Fig. 2(a). In the homogeneous phase U > U (ε), the caloric curve is almost linear, while in the clustered phase c U < U (ε), it is bent downward. Nonetheless, temperature always grows with energy and one does not observe any c negative specific heat energy range. However, as it happens for 3D Newtonian gravity simulations [5], when one reducesthesofteningparameter,anegativespecificheatphasedevelops. Forinstance,inFig.2(b),weshowtwocases at small ε where three phases can be identified [10]: a low-energy clustered phase for U <U (ε), where U is defined as the energy at which ∂T/∂U =0. top top • an intermediate-energy phase, U (ε)<U <U (ε), with negative specific heat. top c • 3 Figure2: CaloriccurvesoftheSelfGravitatingRing(SGR)modelobtainedfromnumericalsimulationsofHamiltonian1. Panel (a)referstothesofteningparametervalueε=10,forwhichasecondorderphasetransition appearsatUc. Nobackbendingof thecaloric curve,indicating anegativespecificheat,is present. Simulations wereperformed for N =100. Panel (b)showsthe caloriccurvesfortwodifferentvaluesofthesofteningparameter,ε1 =1.010−6 andε2 =2.510−7,andN =100. Thetransition isherefirst orderin themicrocanonical ensemble(see Sec.VIfora discussion). Thetwotransition energies Uc(ε1)and Uc(ε2) are prettyclose, suggesting a slow variation of thecritical energy with thesoftening parameter ε. On thecontrary,Utop(ε1) is significantly smaller than Utop(ε2), indicating that this characteristic energy value diminishes with ε. A negative specific heat phase appears for Utop<U <Uc, and expandsas thesoftening parameter is reduced. a high-energy gaseous phase for U (ε)<U. c • The clustered phase is created by the presence of softening ε, without which the particles would fall into the zero distance singularity. In the gas phase, the particles are hardly affected by the potential and behave as almost free particles. The intermediate phase is expected to show the characters of gravity, persisting and even widening in the ε 0 limit. → In the following, several of these features will be given a theoretical explanation and we will detail the analysis of the nature of the phase transition (first or second order) when ε is varied. 4 III. STATIONARY DENSITY DISTRIBUTION In the mean-field limit (N with fixed length [16]), one can introduce the single particle distribution function →∞ f(p,θ)suchthatf(p,θ)dpdθ isthefractionofparticlesinthe domain[θ,θ+dθ][p,p+dp]. Intermsoff,thepotential energy can be written as 1 ′ ′ E [f] = dθ dφ dp dp f(p,θ)V (θ φ)f(p,φ) (4) P ε 2 − Z 1 = dφ dθ ρ(θ)ρ(φ)V (θ φ) (5) ε 2 − Z where ρ(θ)= dp f(p,θ) (6) Z is the mass density. The kinetic energy is 1 E [f]= dθ dp p2f(p,θ) (7) K 2 Z and the total energy E[f]=E [f]+E [f]. (8) K P The equilibrium distribution in the microcanonical ensemble is determined by maximizing entropy S[f]= dθ dp flogf (9) − Z under the constraints of fixed total energy, momentum and mass. In the following, we fix the total energy E[f]=U, the total mass M[f]= ρ dθ =1 (10) Z and the total momentum p[f]= pf(p,θ)dθ dp=0. (11) Z A necessary condition to get an entropy maximum is to require that the free energy F[f] S[f] βE[f] α f dp dθ γp[f], (12) ≡ − − − Z where α, β and γ are Lagrange multipliers, is stationary δF[f] p2 = logf 1 β +W(θ) α γp=0, (13) δf − − − 2 − − (cid:18) (cid:19) where W(θ) is defined as +π W(θ) ρ(φ)V (θ φ)dφ. (14) ε ≡Z−π − Since p[f]=0, the Lagrange multiplier γ vanishes. From Eq. (13), the normalized stationary distribution function can be written as p2 f(p,θ)=Aexp β +W(θ) , (15) − 2 (cid:20) (cid:18) (cid:19)(cid:21) 5 where A=exp( 1 α) is the normalization constant and the mass density is given by − − ρ(θ)=A e−βW(θ), (16) where A=A 2π/β. When Eq. (14) and (16) are combined, we obtain the consistency equations e p +π e W(θ) = A e−βW(φ)V (θ φ)dφ, (17) ε Z−π − and the equilibrium density equation e +π ρ(θ) = Aexp βA ρ(φ)V (θ φ)dφ , (18) ε (cid:20)− Z−π − (cid:21) which are solved numerically in the followineg. Once tehe stationary mass distributions ρ and the function W are obtained for each value of ε, the full single particle distribution function f(θ,p) is derived from Eq. (15). The potential energy and the kinetic energy are determined by Eq. (5) and Eq. (7) respectively, allowing to draw the caloric curve by plotting T β−1 =2E against the total energy U =E +E . K K P ≡ IV. AN ITERATIVE METHOD TO SOLVE THE EQUILIBRIUM DENSITY EQUATION The inverse temperature β can be expressed in terms of the energy +π +π −1 β = 2U ρ(θ)ρ(φ)V (θ φ)dθdφ . (19) ε (cid:26) −Z−π Z−π − (cid:27) Onceaninitialdensitydistributionρ (θ)ischosen,onecandetermineaninitialinversetemperatureβ usingEq.(19), 0 0 andthensolvetheconsistencyequation(18)iteratively(asdoneforinstanceinRef.[17]). However,wewillfollowhere a different iterative method, which ensures entropy increase and, hence, convergence of the algorithm. The method is inspired by a similar one used by Turkington and Whittaker [14] to compute entropy maxima for two dimensional turbulence. The functional to maximize S[f] is strictly concave and we must fix both a linear constraint M[f] = 1 and a nonlinearoneE[f]=U. Itis this latternonlinearconstraintwhichmakesthe variationalproblemmoredifficult than usual. Thetricktosolvethisnonlinearproblemconsistsinconsideringalinearizationoftheenergyconstraintaround the distribution function resulting from the previous step in the iterative process. One begins with the normalized distribution f obtained at the kth step of the algorithm. From that, one can k compute the mass density ρ and the average potential W . k k ρ (θ) = dpf (p,θ) (20) k k Z +π W (θ) = dφρ (φ)V (θ φ). (21) k k ε Z−π − The distribution at the next step f will be then determined by solving the following variational problem k+1 δE max S[f] M[f]=1,E[f ]+ (f f )dpdθ U , (22) k k | δf − ≤ (cid:26) Z (cid:12)fk (cid:27) (cid:12) where the functional derivative of the energy is (cid:12) (cid:12) δE p2 = +W (θ). (23) k δf 2 (cid:12)fk (cid:12) This variational problem has a unique solution f (cid:12) , since it corresponds to the maximization of a strictly concave k+(cid:12)1 functional with linear constraints. Thisiterativeprocessensuresconvergenceoftheentropy. Letusproveit. ByusingageneralizationoftheLagrange multiplier rule for our inequality constrained variational problem [18, 19] δS δE =α +β (24) k+1 k+1 δf δf (cid:12)fk+1 (cid:12)fk (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 6 with the additional requirement δE β E fk + (f f )dpdθ U =0, (25) k+1 k+1 k " Z δf (cid:12)fk − − # (cid:2) (cid:3) (cid:12) (cid:12) where β 0 is the multiplier associated with th(cid:12)e energy constraint and α , the one associated with mass k+1 k+1 ≥ conservation. When solving Eq. (25), we have either β = 0, which removes the energy constraint, or β > 0, k+1 k+1 and an equality for the linearized energy constraint. In order to prove convergence of the entropy, let us first prove that the energy functional E[f] is concave. Since the kinetic part is linear in f, the second variation of E[f] is δ2E = dφdθδρ(θ)δρ(φ)V (θ φ) (26) ε − Z = V δρ 2, (27) ε,k k | | k X where the second equality is obtained using the Fourier series expansion for both the mass density variation δρ and the potential V ε 1 +π δρ = dϕ exp(ikϕ)δρ(ϕ) (28) k 2π Z−π 1 +π V = dϕ exp(ikϕ)V (ϕ). (29) ε,k ε 2π Z−π SinceV isevenintheargumentφ,V isarealnumber. Moreover,sinceV 0andV (ϕ)isstrictlyincreasingfor ε ε,k ε ε ≤ 0 ϕ π, it is easy to prove that for any k, V is strictly negative. Hence, from formula (27) the second variation ε,k ≤ ≤ of the energy functional is negative and this functional is strictly concave. On the other hand, the entropy is strictly concave. We have δS 1 (δf)2 S[f +δf] S[f]+ dθdp δf dθdp , (30) ≤ δf − 2 f Z (cid:12)f Z (cid:12) (cid:12) where in the derivation we have used ln(1+x) x x2/2 for(cid:12)x> 1. Applying this latter inequality with f =f k+1 ≥ − − and δf =f f , and using both condition (24) and (25), we obtain k k+1 − 1 (f f )2 k+1 k S[f ] S[f ] β (U E[f ])+ dθdp − , (31) k+1 k k+1 k − ≥ − 2 f Z k where the term proportionalto α vanishes because of mass conservation. k+1 We will now use the concavity of the energy functional E[f]. For k >1, δE E[fk] E[fk−1]+ (fk fk−1)dpdθ. (32) ≤ δf − Z (cid:12)fk−1 (cid:12) (cid:12) As β 0 and E[f ] U, directly from the variational(cid:12)problem (22), Eq. (31) implies that k+1 k ≥ ≤ 1 (f f )2 k+1 k S[f ] S[f ] − dθdp 0. (33) k+1 k − ≥ 2 f ≥ Z k Hence,the entropyhasto increaseforalliteratesafterthesecond. Sincethe entropyisboundedfromabove,ithasto converge. UsingEqs.(31)and(33),onederivesthattheenergyE[f ]convergestoU frombelow. Moreover,assuming k thatf convergestowardf,onecanprovetheconvergenceofthemultiplierstolimitvaluesαandβ 0,whichimplies k ≥ that f verifies Eq. (24) for equilibrium states. Although mathematically one cannot prove the convergence of f, in allpracticalcaseswe willanalyze,it appearsto be verified. For a more thoroughdiscussionof the convergencein the similar case of the Euler equation, see Sec. IV in Ref. [14]. 7 V. IMPLEMENTATION OF THE ALGORITHM We describe in this section the practical implementation of an algorithm which allows the calculation of the stable distribution, using the method described in the previous section. From (24), we obtain p2 f =A exp β +W (θ) , (34) k+1 k+1 k+1 k − 2 (cid:18) (cid:18) (cid:19)(cid:19) where A =exp( α 1) and β are unknown at this stage. Using (20), we get k+1 k+1 k+1 − − ρ (θ)=A e−βk+1Wk(θ). (35) k+1 k+1 where A =A 2π/β . This equation allows to compute W (θ) from Eq. (21) and k+1 k+1 k+1 e k+1 p 1 1 +π e E E[f ] = + ρ (θ)W (θ)dθ. (36) k+1 k+1 k+1 k+1 ≡ 2βk+1 2Z−π Then the multipliers α and β must be computed from Eqs (10) and (25) and, from these, one gets A . In k+1 k+1 k+1 order to compute numerically these Lagrange parameters, let us define the Lagrangian[19] e δE L [f](β,α)= S[f]+β E + (f f )dpdθ U +α(M[f] 1). (37) k k k − " Z δf (cid:12)fk − − # − (cid:12) (cid:12) From this, one further defines (cid:12) L⋆(β,α)=inf L [f](β,α) . (38) k f { k } One can prove on a general ground [19] that L⋆ is concave and that α and β are the unique maxima of L⋆. k k+1 k+1 k Using condition (34) for the extrema of L [f](β,α), we can compute L⋆. We obtain, using for practical reasons the k k variable A instead of α, 1 1 +π e L⋆(β,A) = logA+ logβ β U +E A dθe−βWk(θ), (39) k 2 − (cid:18) k− 2βk(cid:19)− Z−π Necessary conditions for theeconcave fuenction L⋆ to be maximal are e k ∂L⋆ 1 +π k = dθe−βWk(θ) =0, (40) ∂A A −Z−π ∂L⋆ 1 1 +π k = + U E +A dθW (θ)e−βWk(θ) =0. (41) e e k k ∂β 2β 2βk − − Z−π Substituting Eq. (40) into Eq. (41), one gets the condition e +π dθ W (θ)e−βk+1Wk(θ) k 1 1 +U +E Z−π =0, (42) − 2β − 2β k− +π k+1 k dθe−βk+1Wk(θ) Z−π which,sinceL⋆ isconcave,hasauniquesolution. Numerically,thesolutionβ isfoundbyusingaNewtonalgorithm k k+1 for Eq.(42). Then, fromEq.(40), we get A . Finally, we cancalculate the new density distribution fromEq.(35). k+1 e VI. DISCUSSION OF THE RESULTS Using the iterative method described in the previous section, we are able to derive the stable mass density ρ(θ) solution of Eq. (18) and, from that, all thermodynamic functions in the microcanonicalensemble. In the first part of this section, we will show the numerical solution obtained for ρ(θ), and its dependence on energy for a small value of the softening parameter ε. In the second part, we will discuss the phase diagram of the SGR model, both in the microcanonical and in the canonical ensemble, when ε is varied. 8 A. Mass density, entropy and caloric curves For energies above a certain critical value U (ε), the stable mass density solution is uniform. In this case, one can c compute the entropy from Eq. (9) 1 S = (3log(2π)+1 logβ), (43) 2 − and the inverse temperature from Eq. (19) −1 β = 2U 2E , (44) p − where (cid:0) (cid:1) 1 1 +π +π E = dθdφV (θ φ) (45) p 2(2π)2 Z−π Z−π ε − 1 1 2 = , (46) −π√2√2+ε K ε+2 (cid:18) (cid:19) where is the complete elliptic integral of the first kind (x) π/2dθ/ 1 xsin2θ. K K ≡ 0 − RemarkthatEq.(44)impliesthatthehomogenousstatecannotbecontinuedbelowU =E ,becausethislatter R p hom p energy corresponds to zero temperature. For U <U (ε), the stable mass distribution must be determined numerically. We have checked in this case, that a c direct iterative method of solution of the consistency Eqs (17) and (18) does not always converge. On the contrary, the novel algorithm presented in Sec. V ensures convergence as shown in Fig. 3 for the entropy. 3.6 3.5 3.4 3.3 3.2 S 3.1 3 2.9 2.8 2.7 0 10 20 30 40 50 60 70 80 90 100 k Figure 3: Convergence of the entropyusing the algorithm of Sec. V for ε=10−5,U = 1. − In Fig. 4 we show both entropy and temperature T = β−1 as a function of energy U. The most striking feature is the presence of a negative specific heat region for U U U . For U U U , the entropy does not top c low high ≤ ≤ ≤ ≤ coincide with its convex envelope. Hence, microcanonical and canonical ensemble do not give the same predictions. Indeed, the mainpeculiarityofthe microcanonicalensemble isthatmacroscopicstates withinthis intervalarestable, while they would be either metastable or unstable in the canonical ensemble. The mass density is uniform above U , while, below this value, it is localized. The appropriate order parameter to characterize this localization is the c “magnetization” +π B = dθeiθρ(θ), (47) Z−π which vanishes if the mass distribution is uniform while it reaches the value B =1 when the mass is concentrated in only one point. Intermediate degrees of localizationgive intermediate values of B. The “magnetization”is plotted in Fig.5asafunctionofU. ItisadecreasingfunctionofU,uptoU ,whereithasajumptothelimitingvalue0. Hence, c we havehere a firstordermicrocanonicalphase transition. The firstordernature ofthe phase transitionis confirmed 9 Figure 4: Temperature (panel (a)) and entropy (panel (b)) versus energy U for the softening parameter value ε=10−5. Four valuesoftheenergy,indicatedbytheshort-dashedverticallines,canbeidentifiedfromthispicture: U 93andU 6 low high ≃− ≃ bound from below and above the region of inequivalence of ensembles. Uc 0 is the transition energy in the microcanonical ≃ ensemble. Utop 66 limits from below the negative specific heat region, where temperature decreases as energy increases. ≃ − Tcan 15, represented with a dashed line in panel (a), is the canonical transition temperature and corresponds to the inverse ≃ slope of the entropy, both at U and U , as represented by the straight dashed line in panel (b). The full lines represent low high theanalyticalsolutionsofthetemperatureandofentropyintheuniformcase(seeformulas (43)and(44)). Theyareextended slightly below Uc, in the metastable phase, in order to identify them. The insets in panels (a) and (b) show a zoom of the temperature and of the entropy around Uc, revealing a temperature jump at Uc and different slopes of the entropy above and below Uc, which emphasizes thefirst order natureof thephase transition. zoomingtheentropyaroundU (seetheinsetinpanel(b)ofFig.4). Thisrevealsthatthisfirstorderphasetransition c isofthe convex-concavetype(seeRef.[13]). The canonicalensembleisobtainedbytakingtheconvexenvelopeofthe microcanonicalentropy. The transitionis firstorderinthecanonicalensembleandthe transitiontemperatureT is can given by the inverse slope of the entropy at U and U . No canonical macrostate is present in the energy range low high [U ,U ]. low high Atypicallocalizedmassdensity distributionisshowninFig.6. Itcorrespondstoanenergywherethe specific heat 10 Figure 5: “Magnetization” B versus energy U for ε=10−5, which emphasizes the microcanonical first order phase transition at Uc 0 by showing a jump in theorder parameter. ≃ is negative. 100 10 1 0.1 0.01 0 Figure 6: A typical mass density distribution ρ(θ) for ε=10−5 and U = 20.0, in thenegative specific heat region. − The firstorderphase transitionisassociatedwith the existenceofmetastablestates. Using acontinuationmethod, we have been able to compute them. Their entropy is represented in Fig. 7 around the transition energy U for the c specific case ε = 10−5. The inhomogeneous metastable state turns out to exist for U U U with U 0.16, c in in while the homogeneous metastable state exists for U U U , with U =E ε=≤10−≤5 1.19. ≃ hom c hom p ≤ ≤ ≃− (cid:0) (cid:1) B. Behaviour as the softening parameter ε is varied Let us first examine a situation where the softening parameter is much larger than previously, ε = 10−2. In the microcanonical ensemble, Fig. 8 shows that a concavity change still occurs at U 0.8, and that a phase top ≃ − transition exists at U = U 0.3. However, the temperature being now a continuous function of the energy but c ≃ − with discontinuous derivative at U , the phase transition is of second order, and is associated with the symmetry c breaking of the order parameter. The caloric curve shows that this second order phase transition is of the convex- concave type. As it is necessary for this type of microcanonical second order phase transition [13], we observe a