Density of states of the XY model: an energy landscape approach
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec Density of states of the XY model: an energy landscape approach Cesare Nardini ∗ Laboratoire de Physique, ´Ecole Normale Sup´erieure de Lyon and CNRS, 46 all´ee d’Italie, F-69007 Lyon, France
Rachele Nerattini † and Lapo Casetti ‡ Dipartimento di Fisica e Astronomia and Centro per lo Studio delle Dinamiche Complesse (CSDC),Universit`a di Firenze, and Istituto Nazionale di Fisica Nucleare (INFN),Sezione di Firenze, via G. Sansone 1, I-50019 Sesto Fiorentino (FI), Italy (Dated: December 17, 2013)Among the stationary configurations of the Hamiltonian of a classical O( n ) lattice spin model,a class can be identified which is in one-to-one correspondence with all the the configurations ofan Ising model defined on the same lattice and with the same interactions. Starting from thisobservation it has been recently proposed that the microcanonical density of states of an O( n )model could be written in terms of the density of states of the corresponding Ising model. Later, ithas been shown that a relation of this kind holds exactly for two solvable models, the mean-field andthe one-dimensional XY model, respectively. We apply the same strategy to derive explicit, albeitapproximate, expressions for the density of states of the two-dimensional XY model with nearest-neighbor interactions on a square lattice. The caloric curve and the specific heat as a function of theenergy density are calculated and compared against simulation data, yielding a very good agreementover the entire energy density range. The concepts and methods involved in the approximationspresented here are valid in principle for any O( n ) model. I. INTRODUCTION
A quite recent point of view on the study of the thermodynamic properties of a classical N − body Hamiltoniansystem, commonly referred to as the “energy landscape” approach [1], implies the study of the properties of thegraph of the energy function H : Γ N → R , where Γ N is phase space of the system. Examples of applications includedisordered systems and glasses [2, 3], clusters [1], biomolecules and protein folding [4–6]. In this context a specialrˆole is played by the stationary points of the energy function. The latter are the points p c in phase space such that ∇ H ( p c ) = 0. Based on knowledge about the stationary points of the energy function, landscape methods allow toestimate both dynamical and equilibrium properties of a system. In the present paper the attention will be focusedonly on equilibrium properties.The natural setting to understand the connection between the stationary points of the Hamiltonian and equilibriumstatistical properties is the microcanonical one [7]. This is the statistical description in which the system is consideredas isolated; the control parameter is the energy density ε = E/N , and the relevant thermodynamic potential is theentropy density s N given by s N ( ε ) = 1 N log ω N ( ε ) , (1)where ω N ( ε ) is the density of states of the system. The density of states is given by ω N ( ε ) = Z Γ N d Γ N δ ( H − N ε ) = Z Γ N ∩ Σ ε d Σ ε |∇ H | , (2)where Σ ε is the hypersurface of constant energy N ε and d Σ is the Hausdorff measure; the rightmost integral stemsfrom a co-area formula [8]. When a stationary configuration p c is considered, the integrand in (2) diverges. In finite- N systems this divergence is compensated by the measure d Σ that shrinks in such a way that the integral in Eq. (2)remains finite. This notwithstanding, the density of states is non-analytic at the corresponding value of the energydensity ε c = H ( p c ) /N . Although these nonanalyticities typically disappear in the thermodynamic limit N → ∞ , insome special cases they may survive and give rise to equilibrium phase transitions [9–11]. Stationary points of the ∗ [email protected] † rachele.nerattini@fi.infn.it ‡ lapo.casetti@unifi.it energy correspond to topology changes in the accessible phase space, and this has suggested a relation between thesetopology changes and equilibrium phase transitions (see Refs. [12–14] and references therein).In [15] an approximate form for the density of states of a paradigmatic class of classical spin models, the O( n )models, was conjectured on the basis of energy landscape considerations. According to the conjecture, the density ofstates ω ( n ) of a classical O( n ) spin model on a lattice can be approximated in terms of the density of states ω (1) ofthe corresponding Ising model, i.e., an Ising model defined on the same lattice and with the same interactions. Thecrucial observation is that all the configurations of the corresponding Ising model are stationary configurations of aO( n ) model for any n . The density of states would then be given by ω ( n ) ( ε ) ≈ ω (1) ( ε ) g ( n ) ( ε ) , (3)where g ( n ) ( ε ) is an unknown function representing the volume of a neighborhood of the Ising configuration in the phasespace of the O( n ) model. Since the function g ( n ) comes from the evaluation of local integrals over a neighborhoodof the phase space, one expects that it is regular in all the energy density range. Equation (3) will be discussed indetail in Sec. II. The assumptions made in [15] to derive Eq. (3) are rather crude and uncontrolled. However, werethis relation exact, there would be an interesting consequence: the critical energy densities of the phase transitions ofall the O( n ) models on a given lattice would be the same and equal to that of the corresponding Ising model.Despite the fact that Eq. (3) is approximate[16] the critical energy densities are indeed, if not equal, very close toeach other, whenever a phase transitions is known to take place, at least for ferromagnetic models defined on regular d − dimensional hypercubic lattices. In particular, the critical energy densities are the same and equal to the Ising onefor all the O( n ) models with long-range interactions, as shown by the exact solution [17], and the same happens forall the O( n ) models on an one-dimensional lattice with nearest-neighbor interactions. For what concerns O( n ) modelswith nearest-neighbor interactions defined on (hyper)cubic lattices with d >
1, only numerical results are availablefor the critical energy densities. In d = 2 a Bereˇzinskij-Kosterlitz-Thouless (BKT) phase transition is present in the n = 2 case, occurring at a value of the energy density that differs of about 2% from the Ising case (see [15, 18] fora deeper discussion[19]). In d = 3, the difference between the critical energy densities of the O( n ) models and thoseof the Ising model is smaller than 3% for any n , including the n = ∞ case, as recently shown in [20]. In the mostfrequently considered cases, that is for n = 2 (the XY model), n = 3 (the Heisenberg model) and for the O(4) model,it becomes less than 1%.In two particular cases, where the prediction on the critical energy densities is exact, the derivation of Eq. (3) can befollowed rigorously, that is for the mean-field and for the 1 d ferromagnetic XY models. In these cases, an expressionvery similar to Eq. (3), and which reduces to the latter when ε → ε c , can be derived exactly in the thermodynamiclimit [18]. The technical aspects of the derivation strongly rely on the peculiarities of the two models and especiallyon the fact that they are exactly solvable in the microcanonical ensemble. This feature is crucial for the analysispresented in [18] and its generalization to O( n ) models with n > d > ω ( n ) of nearest-neighbor ferromagnetic O( n ) models in d > d XY models andwill be referred to as “first-principle” approximation in the following; the second one will be named “ansatz-based”approximation, its starting point being the ansatz on the form of the density of states given by Eq. (3) supposedto be valid in all the energy density range. These techniques can be applied to estimate the density of states of inprinciple any O( n ) model in d >
1. However we have explicitly performed the calculations only for n = 2, i.e., for the XY model in d = 2; the technical aspects of the generalization to other members of the O( n ) class in d > g ( n ) . A possible way in which this can bedone will be presented in Sec. III A. In Secs. IV and V the “first-principle” and the “ansatz-based” approximationswill be, respectively, discussed in details through their application to the XY model in d = 2. Some remarks on thetechnical limits of the “first-principle” approach will be included in Sec. IV A. The paper ends with some concludingremarks in Sec. VI. II. STATIONARY POINTS AND DENSITY OF STATES
In our analysis we are going to consider some special cases of classical O( n ) spin models defined on a regular squarelattice in d = 2 and with periodic boundary conditions. In these models, to each lattice site i an n -component classicalspin vector S i = ( S i , . . . , S ni ) of unit length is assigned. The energy of the model is given by the Hamiltonian H ( n ) = − J X h i,j i S i · S j = − J X h i,j i n X a =1 S ai S aj , (4)where the angular brackets denote a sum over all pairs of nearest-neighbor lattice sites. The exchange coupling J willbe assumed positive, resulting in ferromagnetic interactions, and without loss of generality we shall set J = 1 in thefollowing. The Hamiltonian (4) is globally invariant under the O ( n ) group. In the case n = 1, the symmetry group O (1) ≡ Z is a discrete one and the Hamiltonian (4) becomes the Ising Hamiltonian H (1) = − X h i,j i σ i σ j , (5)where σ i = ± ∀ i . In all the other cases n ≥
2, the O ( n ) group is a continuous one. For n = 2 we obtain the XY model that will be the subject of our study. In this model the spins are constrained on the unit circle S (1) and thecomponents of the i th spin can be parametrized by a single angular variable ϑ i ∈ [0 , π ) such that ( S i = cos ϑ i ,S i = sin ϑ i . (6)The Hamiltonian of the XY model can thus be conveniently written as H (2) = − N X i =1 X j ∈N ( i ) cos ( ϑ i − ϑ j ) , (7)where N ( i ) denotes the set of nearest neighbors of lattice site i . The energy density ε = H (2) /N lies in the energyrange [ − d, d ] where d is the lattice dimension.Let us now recall the derivation of Eq. (3) made in [15]. The stationary configurations of H ( n ) for n ≥ S = ( ¯ S , . . . , ¯ S N ) of the N vector equations ∇ H ( n ) = 0 , (8)such that the constraint P a =1 n ( S ai ) = 1 is satisfied ∀ i = 1 , . . . , N . To find explicitly all the stationary configurationsis an essentially impossible task but in some special cases, see e.g. [21], or for very small systems, see e.g. [22] andreferences therein. However, as shown in Ref. [15], a particular class of solutions can be found, given by all theconfigurations in which the spins are parallel or anti-parallel to a fixed direction (say the n th direction): S i = . . . = S n − = 0 ∀ i . Indeed, in this case, the constraint P a =1 n ( S ai ) = 1 can be fulfilled by choosing S ni = σ i ∀ i andthe stationary points equations (8) are satisfied by any of the 2 N possible choices of the σ ’s. The Hamiltonian (4)becomes the Ising Hamiltonian (5) on this class of stationary configurations. Therefore a one-to-one correspondencebetween a class of stationary configurations of the Hamiltonian (4) and all the configurations of the Ising model[23];the corresponding stationary values are the energy levels of the Ising Hamiltonian. We shall refer to the class ofstationary configurations ¯ S = (0 , . . . , , σ i ) ∀ i = 1 , . . . , N as “Ising stationary configurations”. Although this classdoes not include all the stationary points of H ( n ) , see e.g. [24], the 2 N Ising ones are a non-negligible fraction ofthe whole, especially for large N . The total number of stationary configurations is expected to be O ( e N ) [21, 25].When ferromagnetic O( n ) models defined on regular d − dimensional lattices are considered, as in our case, in thethermodynamic limit N → ∞ the energy density levels of the Ising Hamiltonian (5) become dense and cover thewhole energy density range of all the O( n ) models. This latter observation, together with the above mentionedproperties of the Ising points, suggests that Ising stationary configurations are the most important ones, so that thedensity of states ω ( n ) of an O( n ) model may be approximated in terms of these configurations.Let us then rewrite the density of states ω ( n ) as a sum of integrals over a partition of the phase space, ω ( n ) ( ε ) = X p Z U p δ ( H ( n ) − N ε ) d Γ , (9)where p runs over the 2 N Ising stationary configurations, U p is a neighborhood of the p -th Ising configuration suchthat { U p } N p =1 is a proper partition of the configuration space Γ N , that coincides with phase space for spin modelsdefined by the Hamiltonian H ( n ) in Eq. (4).Two approximations were introduced in [15] to derive Eq. (3) from Eq. (9): ( i ) it was assumed that the integralsin Eq. (9) depend only on ε , i.e., the neighborhoods U p can be chosen, or deformed, such as g ( n ) ( ε, p ) = Z U p δ ( H ( n ) − N ε ) d Γ = g ( n ) ( ε, q ) = Z U q δ ( H ( n ) − N ε ) d Γ = g ( n ) ( ε ) , (10)for any p, q such that H ( n ) ( p ) = H ( n ) ( q ) = N ε ; ( ii ) Since non-Ising stationary configurations have been neglected inthis analysis, only neighborhoods centered around stationary configurations at energy density ε have been retainedin the sum (9). These two assumptions immediately lead to Eq. (3), that we rewrite for convenience: ω ( n ) ( ε ) ≈ ω (1) ( ε ) g ( n ) ( ε ) . (11)Both assumptions are needed to derive Eq. (11), and are strictly related to each other. However, these two assumptionsmight well play a very different rˆole. In [18] it has been shown that in the two analytically tractable special cases,the mean-field and the one-dimensional XY models, assumption ( ii ) does not hold in general: it holds only when ε → ε ( n ) c . As a consequence, one should include also stationary configurations with energy ε ′ = ε in the sum. Byintroducing the continuous function G ( n ) ( ε, ε ′ ) = Z U p | H ( n )( p ) N = ε ′ δ (cid:16) H ( n ) − N ε (cid:17) d Γ (12)such that G ( n ) ( ε, ε ) = g ( n ) ( ε ) , (13)the density of states can be written as ω ( n ) ( ε ) = ω (1) (˜ ε ) G ( n ) ( ε, ˜ ε ) . (14)In the above expression ˜ ε is a suitable function of ε . If ˜ ε = ε , then using Eq. (13) one recovers Eq. (11). This preciselyhappens in the mean-field and 1- d XY models when ε → ε ( n ) c , that is at the phase transition.To compute G ( n ) ( ε, ˜ ε ) the microcanonical solutions of the models are needed, as shown in [18]. For this reason Eq.(14) becomes of little use when O( n ) models in d > ε = ε ( n ) c , aswe are going to discuss in detail in the following Sections. III. FIRST-PRINCIPLE AND “ANSATZ-BASED” APPROXIMATIONS: AN OVERVIEW
Let us now discuss how the results recalled in Sec. II can be generalized to O( n ) spin models in d >
1. Our resultscan be grouped into two different categories according to the approach involved in their derivation.The first attempt in the generalization is to set up an approximation procedure that allows to estimate each integralappearing in the sum in Eq. (9), that yields the density of states ω ( n ) ( ε ). This approach is the direct generalization ofthe techniques already applied to the mean-field and to the one-dimensional XY models to generic short-range O( n )models except that, for short range systems, only an approximate evaluation of each of these integrals is possible.This approximation protocol will be denoted as “first-principles” approximation, its starting point being simply thedensity of states ω ( n ) expressed in terms of a suitable partition of the phase space Γ of the system. The procedurewill be presented in Sec. IV and the two-dimensional XY model will be considered as a test model of our analysis.The second approach starts from the ansatz on the form of the density of states given by Eq. (11) with g ( n ) ( ε ) as inEq. (10). In this approach it is assumed that the integral in Eq. (10) does not depend on the specific point consideredbut only on its energy density ε and that only Ising points with energy density ε ′ = ε contribute to the density ofstates ω ( n ) ( ε ) in Eq. (11). We will denote this kind of analysis as “ansatz-based” approximation. The general aspectsof the method will be presented in Sec. V and its application to the XY model in two dimensions will be discussedin detail.Both in the “ansatz-based” and “first-principle” approximations, the main point is to estimate the quantity g ( n ) ( ε, p ) = Z U p d Γ δ (cid:16) H ( n ) − N ε (cid:17) = Z U p d Γ δ − N X i =1 X j ∈N ( i ) n X a =1 S ai S aj − N ε , (15) FIG. 1. Index convention for the nearest-neigbors spins ϑ (1) i and ϑ (2) j of ϑ i for the XY model on a square lattice with edge L = √ N and periodic boundary conditions. The spins in the lattice sites i + 1 and i + L are said to be second neighbors. where the expression of H ( n ) given in (4) has been explicitly introduced. Among the O( n ) class of models, Eq. (15)can be analytically evaluated only in the cases of the mean-field and one-dimensional XY models; its computationin the general case being as difficult as to find an exact solution of the models. However some computationalprocedures can be set up to carry on the calculations, albeit approximate. In [26] an approximation scheme has beenintroduced, named the Local-Mean-Field (LMF) approximation. This is a model-dependent procedure that allows toreduce the N − dimensional integral in Eq. (15) over the configuration space Γ to N one-dimensional integrals overuncoupled variables. The uncoupling procedure becomes possible once suitable model-dependent collective variablesare defined[27]. The LMF procedure will be presented in Sec. III A. We shall restrict the presentation to the case n = 2and to two-dimensional square lattices with periodic boundary conditions. The generalization to higher-dimensionallattices and different values of n should be possible along similar lines, but we did not work this out in detail althoughsomething similar (but in the canonical ensemble) has been done in [26] for the case n = 2 and d = 3.Other schemes than LMF for approximating the g ( n ) ’s may be used. However, their implementation may bequite complicated in the case of the “first-principle” approximation, while feasible in the “ansatz-based” one. Thereason for this will become clear in Sec. IV A. We will give an example of an alternative scheme, corresponding to aharmonic approximation of the Hamiltonian around each Ising stationary point, while discussing the “ansatz-based”approximation in Sec. V. A. The Local Mean-Field (LMF) approximation
In the case n = 2 and d = 2, i.e. for the XY model in two spatial dimensions, the parametrization (6) can bechosen such that the integral in Eq. (15) becomes Z U p dϑ . . . dϑ N δ − N X h i,j i cos ( ϑ i − ϑ j ) − N ε == Z U p dϑ . . . dϑ N δ − N X i =1 cos ϑ i X j =1 cos ϑ ( j ) i + sin ϑ i X j =1 sin ϑ ( j ) i − N ε == Z ∞−∞ dk π e − ikNε Z U p dϑ . . . dϑ N e − ik P Ni =1 h cos ϑ i P j =1 cos ϑ ( j ) i +sin ϑ i P j =1 sin ϑ ( j ) i i (16)where we have introduced the integral representation of the δ − function, δ ( x ) = π R ∞−∞ e − ikx dk ; ϑ ( j ) i denotes thenearest-neighbor spin j of the lattice site i considered according to the convention reported in Fig. 1. The two spinsidentified by ϑ (1) i and ϑ (2) i are said to be second-neighbor spins. A generalization of Eq. (16) to n > d > p denotes any Ising stationary point; the LMF approximation can then beintroduced as follows.Let us consider a particular Ising point, i.e., ¯ p = (cid:0) ¯ ϑ , . . . , ¯ ϑ N (cid:1) with ¯ ϑ i ∈ { , π } ∀ i = 1 , . . . , N . For each angularvariable ϑ i , in Eq. (16) we replace the N − ϑ j with j = i with their values ¯ ϑ j at point ¯ p . The variable ϑ i is left free to vary in all the range specified by U ¯ p ( ϑ i ) . In this way the angular variables in Eq. (16) become uncoupledand Z ∞−∞ dk π e − ikNε Z U p dϑ . . . dϑ N e − ik P Ni =1 h cos ϑ i P j =1 cos ϑ ( j ) i +sin ϑ i P j =1 sin ϑ ( j ) i i == Z ∞−∞ dk π e − ikNε N Y i =1 Z U ¯ p ( ϑi ) dϑ i e − ik cos ϑ i (cid:16) cos ¯ ϑ (1) i +cos ¯ ϑ (2) i (cid:17) (17)with H (¯ p ) = N ε ′ . Eq. (17) clarifies the choice of “Local Mean-Field” as the name for the approximation. Indeed, afterLMF approximation is applied, the angular variables ϑ i ’s in Eq. (17) become independent variables in U ¯ p ( ϑ i ) ; eachdegree of freedom interacts only with a sort of local mean-field generated by the spins located in the nearest-neighborlattice sites. Remarkably, with this approximation the contribution of the sines in the exponent of Eq. (17) vanishes,since sin ϑ i = 0 when ϑ i ∈ { , π } , and the expression in Eq. (17) simplifies.Since the analysis is restricted to Ising stationary configurations, in Eq. (17) only two cases are possible:(a) The second-neighbor spins are equal. In this case cos ¯ ϑ (1) i = cos ¯ ϑ (2) i = ±
1, and so Z U ¯ p ( ϑi ) dϑ i e − ik cos ϑ i (cid:16) cos ¯ ϑ (1) i +cos ¯ ϑ (2) i (cid:17) = Z U ¯ p ( ϑi ) dϑ i e ∓ ik cos ϑ i ; (18)(b) The second-neighbor spins are opposite. In this case cos ¯ ϑ (1) i = − cos ¯ ϑ (2) i , and so Z U ¯ p ( ϑi ) dϑ i e − ik cos ϑ i (cid:16) cos ¯ ϑ (1) i +cos ¯ ϑ (2) i (cid:17) = Z U ¯ p ( ϑi ) dϑ i . (19)To evaluate Eqs. (18) and (19), the neighborhoods U ¯ p ( ϑ i ) have to be defined. A priori there are no particular clues onthe best choice of U ¯ p ( ϑ i ) , the only request being that { U p ( ϑ ) , . . . , U p ( ϑ N ) } N p =1 has to be a partition of the phase spaceof the system. In the following, two different choices of U p ( ϑ i ) will be considered: U p ( ϑ i ) = ( [ − π, π ] , if ¯ ϑ i = 0,[0 , π ] , if ¯ ϑ i = π . (20) U p ( ϑ i ) = ((cid:2) − π , π (cid:3) , if ¯ ϑ i = 0, (cid:2) π , π (cid:3) , if ¯ ϑ i = π . (21)In principle the choice in Eq. (20) should be avoided, the neighborhoods U p ( ϑ i ) ’s being not disjoint. We will anyhowconsider it in the following, since it is the easiest choice that can be done a priori for these systems. For the “first-principle” approximation only the first choice for the U p ( ϑ i ) will be considered. The technical reasons for this choicewill be discussed in Sec. IV A. For the “ansatz-based” approximations, instead, both (20) and (21) will be consideredand the effect of neighborhood superposition will be explicitly discussed. IV. “FIRST-PRINCIPLE” APPROXIMATION
Let us consider the form of the density of states ω ( n ) ( ε ) given by Eq. (9). In the following we are going to apply ourprocedure to derive an approximate form of the density of states ω (2) ( ε ) for the XY model in d = 2. A generalizationof these techniques to O( n ) models with n > d > XY model in d = 2 is considered, Eq. (9) can be written as ω (2) N ( ε ) ≃ X p ∈ Γ π Z ∞−∞ dk e − iNkε Z U p dϑ . . . dϑ N e − ik P Ni =1 h cos ϑ i P j =1 cos ϑ ( j ) i +sin ϑ i P j =1 sin ϑ ( j ) i i , (22)where p is any Ising stationary configuration. The integral over the angular variables ϑ i can be evaluated by applyingthe LMF approximation presented in Sec. III A. In this case we choose a definition of the integration neighborhoodsas in Eq. (20). Similar results are supposed to hold also for a choice of U p ( ϑ i ) as in Eq. (21) but the calculationshave not been carried out in this case. The generalization of Eq. (22) to other O( n ) models would depend onthe parametrization chosen to describe the spin variables and on the number of nearest neighbors (that is, on thedimensionality of the lattice).From Eqs. (18) and (19) we have Z π − π dϑ e ∓ ik cos ϑ = 2 πJ (2 | k | ) (23)whenever a couple of equal second-neighbors spins is present in the system, and Z π dϑ π (24)whenever a couple of opposite second-neighbors spins is present in the system; J ( x ) is the zero-order Bessel functionof the first kind [28]. We will denote by N c the number of couples of equal second-neighbor spins[29] and by n c = N c /N its number density.For any given value of N c , a particular family of Ising stationary configurations is selected and the integral in Eq.(22) becomes the product of two different contributions: (2 πJ (2 | k | )) N c due to couples of equal second neighborsspins, and (2 π ) N − N c due to the opposite couples. In this way Eq. (22) becomes ω (2) N ( ε ) ≃ N X N c =0 ν ( N c ) (2 π ) N (1 − n c ) π Z ∞−∞ dk e N [ − i k ε + n c log ( πJ (2 √ k ) )] , (25)where ν ( N c ) is a degeneracy factor counting the number of Ising configurations with a given value N c of the collectivevariable. Due to periodic boundary conditions, determining ν ( N c ) in Eq. (25) reduces to a combinatorial problemanalogous to disposing N c distinct elements over N possible empty spaces; we then have ν ( N c ) = (cid:0) NN c (cid:1) = N ! N c !( N − N c )! .The evaluation of the degeneracy factor ν ( N c ) is crucial for this kind of analyses and we will come back on this pointin Sec. IV A.Since we are interested in the large- N behavior of the system, the integration over k in Eq. (25) can be computedusing the saddle-point method; the saddle-point equation is given by k = iτ with τ satisfying the self-consistencyequation I (2 τ ) I (2 τ ) = ε n c . (26)Eq. (25) can then be written as ω (2) N ( ε ) ≃ N X Nn c =0 N ! N c !( N − N c )! e N [ − τ ε + log 2 π + n c log( I (2 τ ))] ≃ N Z dn c e N [ − τ ε − n c log n c − (1 − n c ) log(1 − n c )+log 2 π + n c log( I (2 τ ))] , (27)where we have replaced the sum over N c with an integration over n c , we have neglected the term − N log 2 π inthe exponent since its contribution vanishes for N → ∞ , and we have introduced the Stirling approximation of thefactorial terms in the binomial coefficient to get N ! N c !( N − N c )! ≃ e N [ − n c log n c − (1 − n c ) log(1 − n c )] . (28)For each value of the energy density ε , we can approximate Eq. (27) as ω (2) N ≃ N e N max nc ∈ [0 , [ − τ ε − n c log n c − (1 − n c ) log(1 − n c )+log 2 π + n c log( I (2 τ ))] . (29)The entropy density s (2) ( ε ) in the thermodynamic limit is finally given by: s (2) ( ε ) ≃ max n c ∈ [0 , f ( n c , τ ) , (30) ò ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò ò ò òò ò ò - - - - Ε T FIG. 2. Temperature T as a function of the energy density ε as derived from Eq. (30) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than symbol sizes. The dashed blue line and the solid red line connecting the two sets of data aremeant as guide to the eyes. with f ( n c , τ ) = − τ ε − n c log n c − (1 − n c ) log(1 − n c ) + + log 2 π + n c log ( I (2 τ )) (31)and τ numerically determined from Eq. (26). The maximization procedure in Eq. (30) can be performed numerically.The temperature T and the specific heat c as a function of the energy density can be computed from Eq. (30) bynumerical differentiation and are shown as red circles in Figs. 2 and 3, respectively. As a comparison, data for T and c as functions of the energy density ε have been computed with a Monte Carlo simulation of the XY model with edge L = 32 in d = 2, performed with the optimized cluster algorithm spinmc provided by the ALPS project [30]. Thenumerical data are plotted in Figs. 2 and 3 as blue squares together with the results obtained from our approximation.Fig. 2 shows that Eq. (30) correctly reproduces the asymptotic behavior of the function T ( ε ) both in the low andin the high energy regime, at a semiquantitative level. In particular, for ε ≃ − ε ≥ − . ε p, ≃ − .
495 markedby the vertical red dot-dashed line, at a slightly lower energy density value than ε p ≃ − .
24 where the peak occursin the simulation data (vertical dashed blue line). The overall shape of the specific heat sketched by our results is inqualitative agreement with the numerical results for ε ∈ [ ε p, , − . ε ≥ − . ε < ε p, the agreement becomes worse. Both from theoretical and numerical results, we know that c ( ε ) → . ε →
0, see i.e. [31, 32]. Our results show an abrupt increase for ε ≃ − .
85. This is a shortcoming of ourapproximation that is still under investigation.
A. The degeneracy factor ν and a different approximation for ω ( n ) The analysis presented in Sec. IV shows that the “first-principle” procedure provides a practical method to derivean approximate form of the density of states ω (2) in two dimensions. The thermodynamic functions derived from Eq.(30) are in reasonably good agreement with the simulations, as shown in Figs. 2 and 3. This suggests that our ideaof considering only Ising stationary points in the derivation of ω (2) ( ε ) is trustworthy and provides a good strategy toapproximate the thermodynamic properties of continuous O( n ) models, in principle for any value of n and d .An important feature of the “first-principle” approximation is the natural emergence of collective variables, like N c , in terms of which the stationary configurations can be parametrized and the density of states re-written as in ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò òò òò ò ò ò ò - - - - Ε c FIG. 3. Specific heat c as a function of the energy density ε as derived from Eq. (30) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than symbol sizes. The dashed blue line and the solid red line connecting the two sets of data aremeant as guide to the eyes. Eq. (25). The number and the type of the collective variables depend on several aspects: the model considered,the dimensionality of the lattice, the definition of the integration neighborhoods U p ( ϑ i ) , the specific Ising point p considered, and the approximation strategy applied to evaluate the integral in Eq. (15). Indeed, as we are going toshow in Sec. V C, instead of applying the LMF approximation other strategies could have been adopted to computethe integral in Eq. (15), like e.g. a harmonic expansion of H ( n ) around the Ising stationary points. In this case otherquantities, as the determinant D ( p ) of the Hessian matrix of H ( n ) or the density of index ι ( p ) = I N ( p ) would haveemerged in the analysis[33].Let us denote by z ( n, d, p ) the vector of collective variables needed in the evaluation. In the “first-principle”approach the density of states ω ( n ) N ( ε ) can always be reduced to a form of the type ω ( n ) N ( ε ) = N X z ( n,d,p ) ν ( z ( n, d, p ) , N ) f ( z ( n, d, p ) , N ) , (32)as shown in Eq. (25). In the above expression ν ( z ( n, d, p ) , N ) is the degeneracy factor associated to the collectivevector of parameters z counting the number of Ising configurations for given values z ( n, d, p ) and N of the collectivevariables and of the number of degrees of freedom of the system, respectively.The degeneracy factor ν can not be evaluated analytically but in some specific cases like the one discussed in Sec.IV and those discussed in [18]. In the analysis presented in Sec. IV only integration neighborhoods as in Eq. (20)have been considered. Indeed, as will be shown in Sec. V B, a choice of the integration neighborhoods as in (21)would have produced the emergence of two different collective variables, N c and N with N denoting the number oftriplets of equal Ising spins. To compute ω (2) , it would then be necessary to estimate the degeneracy factor ν ( N c , N )counting the number of Ising configurations with given values of N, N c and N . Up to our knowledge, this quantityis not analytically known. To carry on the calculation one may then estimate ν ( z ( n, d, p )) numerically. This can bedone for instance by performing a Monte Carlo simulation in which the Ising configurations are grouped and countedaccording to their value of z . This method, although possibly correct, would require a strong computational effortand we preferred to leave it to future work.0On the other hand the ansatz on the form of the density of states ω ( n ) given by Eq. (11) is able to reproduce withunexpected accuracy both the emergence of the phase transitions in the O( n ) models and even the critical energydensity values at which the transitions are located [15]. Then, one may consider Eq. (11) as the new starting pointto approximate the thermodynamic functions of the O( n ) system in the whole energy density range [ − d, d ]. In thiskind of approach, called “ansatz-based” approach, the main point remains the estimation of g ( n ) ( ε ); this implies theemergence of collective variables z ( n, p, d ), as before. This notwithstanding, it is now reasonable to assume that,given a particular Ising point p with energy density ε ′ = H ( n ) ( p ) /N , the possible values of the collective variables z ( n, d, p ) would narrow around a typical value ˜ z ( n, d, p ) when N → ∞ (see e.g. Ref. [24]). In this limit the typicalvalue ˜ z ( n, d, p ) would not depend on p anymore but only on its energy density ε ′ . We then have ˜ z ( n, d, p ) ≃ ˜ z ( n, d, ε ′ ).Since the ansatz in Eq. (11) imposes that only stationary points with energy density ε ′ = ε have to be considered inthe evaluation of ω ( n ) ( ε ), we have that g ( n ) ( z ( d, p )) → g ( n ) ( z ( d, ε )) = g ( n ) ( ε ) in d dimensions, with z ( d, ε ) suitablefunctions that can be easily estimated by fits of numerical simulation data.All these concepts will be clarified in the next Sections where the “ansatz-based” approximation will be applied tothe XY model in d = 2. V. “ANSATZ-BASED” APPROXIMATION
Let us consider the density of states ω ( n ) as given by Eq. (11); then, our purpose is to estimate g ( n ) ( ε ). This canbe done for instance by applying the LMF approximation introduced in Sec. III A. Let us consider the XY model in d = 2 as test model of our procedure so that all the technical tools presented in Sec. III A can be immediately appliedto this case; the generalization to other O( n ) models should be straightforward. A. LMF approximation for g (2) ( ε ) and U p ( ϑ i ) given by Eq. (20). We start considering the LMF approximation with the first choice of the neighborhoods U p ( ϑ i ) given by Eq. (20).In this case the calculations proceed on the same lines as in the case of the “first-principles” approximation, the onlydifference being that only Ising points with energy density ε are considered and the collective variable z (2 , , ε ) = N c ( ε )does not depend on the specific Ising point p anymore but only on its energy density ε .The density of states can then be written as ω (2) N ( ε ) = ω (1) N ( ε ) (2 π ) N (1 − n c ( ε )) π Z ∞−∞ e N [ − ikε + n c ( ε ) log ( πJ ( √ k ))] (33)where ω (1) N ( ε ) plays the rˆole of ν ( N c , N ) in Eq. (25) and is analytically known thanks to the exact solution of the 2- d Ising model. On the other hand, n c = n c ( ε ) is an unknown function that has to be determined. This has been doneinterpolating the numerical data obtained by a Monte Carlo simulation of the two-dimensional XY model with edge L = 32. The result is shown as a solid blue line in Fig. 4, together with the simulation data.The integral in Eq. (33) can be computed with the saddle-point method and the saddle-point equation is given by k = iλ ; λ satisfies a self-consistency equation analogous to Eq. (26) given by I (2 λ ) I (2 λ ) = ε n c ( ε ) . (34)Eq. (33) can then be written as ω (2) N ( ε ) ≃ ω (1) N ( ε ) e N [ − λε +log 2 π + n c ( ε ) log( I (2 λ ))] = e N h s (1) N ( ε ) − λε +log 2 π + n c ( ε ) log( I (2 λ )) i (35)valid for N ≫ s (1) N ( ε ) represents is the entropy density of the two-dimensional Ising model. Dividing by N thelogarithm of the above expression, letting N → ∞ and neglecting the sub-leading terms in N , we finally get thefollowing expression for the entropy density of the XY model in d = 2 s (2) ( ε ) ≃ s (1) ( ε ) + log 2 π − λε + n c ( ε ) log [ I (2 λ )] , (36)where λ satisfies the self-consistency equation (34). In Fig. 5 we plot the temperature T as a function of the energydensity ε as obtained from numerical differentiation of Eq. (36) (red points). As in Fig. 2, the values obtained witha Monte Carlo simulation are shown for comparison (blue squares).1 FIG. 4. Data for n c obtained by a Monte Carlo simulation of the two-dimensional XY model with edge L = 32 (blue solidcircles). Errorbars of the simulation data are smaller than the symbol sizes. The solid blue line is the interpolating function n c ( ε ). ò ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò ò ò òò ò ò - - - - Ε T FIG. 5. Temperature T as a function of the energy density ε as derived from Eq. (36) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than the symbol sizes. The dashed blue line and the solid red line connecting the two sets of dataare meant as guide to the eyes. Fig. 5 shows that the asymptotic behavior of the function T ( ε ) is well reproduced at a semiquantitative level by ourapproximation in the harmonic regime (very low-energy), in the low-energy regime and in the high-energy limit; theagreement is extremely good for low energies. For ε & − . T is about 50%.In Fig. 6 the values of the specific heat c obtained with a numerical differentiation of Eq. (36) are plotted as afunction of the energy density ε . As in Fig. 5 the theoretical results are displayed as red circles and are plottedtogether with the values of c computed by Monte Carlo simulation (blue squares). The specific heat shows a peak for ε p, ≃ − .
258 marked by the vertical red dot-dashed line. This value of the energy density is very close to ε p ≃ − . c is about 20% for ε < ε p and smaller for ε > ε p . The trend of the theoretical results up to ε ≃ − . c ( − ∈ [0 . , . ε ≃ − ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò òò òò ò ò ò ò - - - - Ε c FIG. 6. Specific heat c as a function of the energy density ε as derived from Eq. (36) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than the symbol sizes. The dashed blue line and the solid red line connecting the two sets of dataare meant as guide to the eyes. observed as in Fig. 3. Again, this is due to a shortcoming of the approximation.The calculations presented in this Section have been repeated using the expression of U p ( ϑ i ) given in Eq. (21). Theresults are presented below. B. LMF approximation for g (2) ( ε ) and U p ( ϑ i ) given by Eq. (21). We now consider the LMF approximation with the neighborhoods U p ( ϑ i ) given by Eq. (21). Given an Ising config-uration p , from Eq. (18) we have two possibilities.( i ) if ¯ ϑ (1) i = ¯ ϑ (2) i = 1 (resp. −
1) and ¯ ϑ i = 1 (resp. − · ↑ · · ↓ ·· ↑ ↑ or respectively · ↓ ↓ , · · · · · · (37)then Z π − π e − ik cos ϑ i (cos ¯ ϑ (1) i +cos ¯ ϑ (2) i ) dϑ i = Z π π e − ik cos ϑ i (cos ¯ ϑ (1) i +cos ¯ ϑ (2) i ) dϑ i = π ( J (2 k ) − iH (2 k )) (38)regardless of whether ¯ ϑ i = 0 or π ; H ( x ) denotes the zero-order Struve function.( ii ) If ¯ ϑ (1) i = ¯ ϑ (2) i = 1 (resp. −
1) and ¯ ϑ i = − · ↑ · · ↓ ·· ↓ ↑ or respectively · ↑ ↓ , · · · · · · (39)3then Z π − π e − ik cos ϑ i (cos ¯ ϑ (1) i +cos ¯ ϑ (2) i ) dϑ i = Z π π e − ik cos ϑ i (cos ¯ ϑ (1) i +cos ¯ ϑ (2) i ) dϑ i = π ( J (2 k ) + iH (2 k )) (40)regardless of whether ¯ ϑ i = 0 or π .On the other hand, Eq. (19) is simply given by Z π − π dϑ i Z π π dϑ i π . (41)We will denote by n = N /N the density of triplets of equal spins forming a local configuration as in Eq. (37).Combining Eqs. (38) and (40) with Eq. (17) we get g (2) ( ε ) = π N π Z ∞−∞ dk e ikNε e N [ n ( ε ) log( J (2 k ) − iH (2 k ))] × e N [( n c ( ε ) − n ( ε )) log( J (2 k )+ iH (2 k ))] . (42)In this case, the vector of parameters z (2 , , ε ) is given by z (2 , , ε ) = n c ( ε ) ∪ n ( ε ).The integral in Eq. (42) can be computed with the saddle point method in the large- N limit. The saddle pointequation is given by k = − iζ ; making use of the properties of the Bessel and of the Struve functions and performingsome algebra, Eq. (42) can be written as g (2) ( ε ) ≃ e N [ − ζε +log π + n ( ε ) log( I (2 ζ ) − L (2 ζ ))] e N [( n c ( ε ) − n ( ε )) log( I (2 k )+ L (2 ζ ))] (43)with ζ satisfying the self-consistency equation − ε + n ( ε ) I (2 ζ ) − L (2 ζ ) (cid:18) I (2 ζ ) − (cid:18) π + L − (2 ζ ) + L (2 ζ ) (cid:19)(cid:19) ++ n c ( ε ) − n ( ε ) I (2 ζ ) + L (2 ζ ) (cid:18) I (2 ζ ) + (cid:18) π + L − (2 ζ ) + L (2 ζ ) (cid:19)(cid:19) = 0 . (44)As in the case of n c ( ε ), the function n ( ε ) can be obtained interpolating the numerical data arising from a MonteCarlo simulation of the system; the simulation data for n are not shown here.We can now insert Eq. (43) in Eq. (11). By taking the logarithm of the resulting expression and neglecting thesub-leading term in N , we finally arrive to the following expression for the entropy density s (2) ( ε ) ≃ s (1) ( ε ) − ζ ε + log π + n ( ε ) log [ I (2 ζ ) − L (2 ζ )] + ( n c ( ε ) − n ( ε )) log [ I (2 ζ ) + L (2 ζ )] (45)valid in the N → ∞ limit; ζ has to be determined numerically from Eq. (44).Fig. 7 shows T as function of the energy density ε obtained by numerical differentiation from Eq. (45) (red points).As in Fig. 5 the theoretical values are plotted together with the data obtained by a Monte Carlo simulation of thesystem (blue squares).For Fig. 7 the same comments can be done as for Fig. 5. The theoretical results are in good agreement with thenumerics in the entire energy density range. In comparison with the theoretical caloric curve in Fig. 5, the caloriccurve resulting from this approximation and shown in Fig. 7 is closer to the numerical results: the largest differencebetween theory and simulation is here around 20%. This fact is the effect of the different choice of the integrationneighborhoods. In particular, if the integration neighborhoods are superposed, as in Eq. (20), g (2) ( ε ) is overestimatedand the difference between the theoretical and the numerical results is larger than in the case of a choice of theintegration neighborhoods as in Eq. (21).This fact is even more evident for the specific heat. The theoretical results are plotted in red in Fig. 8. With achoice of the integration neighborhoods of Eq. (42) as in Eq. (21), the energy value of the peak of the specific heatderived with our approximation is ε p, ≃ − . ≃ ε p ≃ − .
24; moreover, the entire high energy regime for ε > ε p is in good quantitative agreement with the numerics. On the other hand, for ε < ε p the two sets of data separatethemselves and in the low energy regime the same considerations as for Figs. 3 and 6 can be done. C. Harmonic approximation for g (2) ( ε ) In this last Section we present an alternative way to estimate the function g (2) appearing in Eq. (11), using aharmonic expansion of the Hamiltonian around each Ising stationary point. This approach can be seen as the natural4 ò ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò ò ò òò ò ò - - - - Ε T FIG. 7. Temperature T as a function of the energy density ε as derived from Eq. (45) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than the symbol sizes. The dashed blue line and the solid red line connecting the two sets of dataare meant as guide to the eyes. ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò òò òò ò ò ò ò - - - - Ε c FIG. 8. Specific heat c as a function of the energy density ε as derived from Eq. (36) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than the symbol sizes. The dashed blue line and the solid red line connecting the two sets of dataare meant as guide to the eyes. g (2) ( ε, ¯ p ) in Eq. (15). The approximation consists in expanding theHamiltonian up to the harmonic order around the Ising configuration ¯ p = ( ¯ ϑ , ... ¯ ϑ N ): g (2) ( ε, ¯ p ) ≃ g (2) har ( ε, ¯ p ) = Z U ¯ p dϑ ...dϑ N δ (cid:18)
12 ( ϑ − ¯ p ) H (¯ p ) ( ϑ − ¯ p ) T − N ( ε − ε ¯ p ) (cid:19) , (46)where H (¯ p ) is the Hessian matrix of H (2) evaluated on the Ising stationary point ¯ p with energy H (2) (¯ p ) = N ε ¯ p . Wenow shift the coordinates to move the stationary point ¯ p to the origin of the coordinate system, a change of variablesto diagonalize the Hessian matrix and a rescaling of the integration variables according to the eigenvalues of theHessian. With this procedure, we have g (2) ( ε, ¯ p ) ≃ g (2) har ( ε, ¯ p ) = 1 p N D (¯ p ) Z U ¯ p dy ...dy N + dx ...dx N − δ N + (¯ p ) X i =1 y i − N − (¯ p ) X i =1 x i − N ( ε − ε ¯ p ) (47)where D (¯ p ) is the determinant of H (¯ p ), N − (¯ p ) is the index of ¯ p , N + (¯ p ) = N − N − (¯ p ). Observe that in the previousexpression, y i with 1 ≤ i ≤ N + (¯ p ) are the variables associated with the directions in phase space corresponding tothe positive eigenvalues and x i with 1 ≤ i ≤ N − (¯ p ) are the variables associated with the directions in the phase spacecorresponding to the negative eigenvalues.To go further, we have to specify in a suitable way the integration neighborhoods U ¯ p . As in the LMF case, thischoice determines the feasibility of the calculation. In the following, we consider U ¯ p as the union of two balls [36]centered on the stationary point. The first one is associated with the directions y i and has radius α , and the secondone to the directions x i and has radius β . With this choice, we have g (2) har ( ε, ¯ p ) = 1 p N D (¯ p ) Z ∞−∞ dy ...dy N + dx ...dx N − δ (cid:0) Y − X − N ( ε − ε ¯ p ) (cid:1) Θ (cid:0) − Y + N + α (cid:1) Θ (cid:0) − X + N − β (cid:1) . (48)where we have introduced the following variables: Y = N + ( p ) X i =1 y i and X = N − ( p ) X i =1 x i . (49)In principle, α and β should be considered as functions of the specific Ising configuration ¯ p . This would howeverintroduce in the theory a huge number of free parameters. Thus, we will assume that α and β are fixed real numbers.Moreover, g (2) har depends only very weakly on α : this can be easily understood by recalling that, in the canonicalensemble, large values in the positive directions are damped in a Gaussian way. For this reason, we can also assumethat α ≫ β .As explained in Section IV A, in the “ansatz-based” approach the continuous weight does not really depend on thespecific Ising stationary point ¯ p but only on its energy density, which corresponds to assume that (cid:26) D (¯ p ) N , N + (¯ p ) N , N − (¯ p ) N (cid:27) ≃ (cid:26) D ( ε ) N , N + ( ε ) N , N − (¯ ε ) N (cid:27) . (50)These functions can be computed numerically as in the case of n c ( ε ) and n ( ε ) discussed in Secs. IV and V, respec-tively (data not shown). We recall that in the “ansatz-based” approach, only continuous weights corresponding toIsing stationary points with energy ε ¯ p = ε contribute to ω (2) . We refer the reader to Appendix A for the detailedcomputation of g (2) . Under these conditions, the approximate density of states becomes ω (2) ( ε ) ≃ exp n N h s (1) ( ε )+1 − log 2 −LD ( ε )+log π − n + ( ε ) log n +( ε )2 − n − ( ε ) log n − ( ε )2 +log( n − ( ε ))+log( β ) io (51)where LD = lim N →∞ N log D , n − = N − /N , n + = 1 − n − . Two observations are mandatory at this level. First, wenote that Eq. (51) does not depend anymore on α ; second, even if the entropy does depend on β , thermodynamicfunctions such as temperature or specific heat do not depend on the free parameter β . Thus, thermodynamic functionsdo not depend on any free parameters.The behavior of the temperature and of the specific heat as a function of the energy density ε derived from Eq.(51) are reported in Fig. 9 and 10, respectively. For the latter figures, similar considerations as in the case of Fig. 5and 6 can be done.6 ò ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò ò ò òò ò ò - - - - Ε T FIG. 9. Temperature T as a function of the energy density ε as derived from Eq. (51) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than the symbol sizes. The dashed blue line and the solid red line connecting the two sets of dataare meant as guide to the eyes. ò ò òòòòòòòòòòòòòòòòòòòòòòòò òòòòòò òò òò ò ò ò ò - - - - Ε c FIG. 10. Specific heat c as a function of the energy density ε as derived from Eq. (51) for the XY model in d = 2. Ourresults (red circles) are plotted together with the data obtained by a Monte Carlo simulation (blue triangles). Errorbars of thesimulation data are smaller than the symbol sizes. The dashed blue line and the solid red line connecting the two sets of dataare meant as guide to the eyes. VI. CONCLUDING REMARKS
In this paper, we have presented two different semi-analytic ways to approximate the density of states, hencethermodynamic functions, of short-range O( n ) models on a regular lattice. We have performed explicit calculations inthe case of the XY model in d = 2, but the generalization of our results to higher dimensions and/or to O( n ) modelswith n > d and the mean-field XY models in [18]. Thesecond one, called “ansatz-based” approximation, assumes Eq. (3) as an ansatz on the form of the density of statesof O( n ) models with n ≥
2. In both cases, the delicate point is to approximate the weights g given by Eq. (15). Twodifferent approximations have been used, the local mean field (LMF) and a harmonic approximation around eachIsing configuration.The XY model in d = 2 is not exactly solvable and any approximation scheme has to be compared with the resultscoming from numerical simulations. Calculated thermodynamic functions always show a qualitative agreement withthe numerical data and in some cases also a quantitative agreement. Particularly interesting is the presence of thepeak in the specific heat reported in Figs. 3, 6 and 8. Despite the approximations involved in its derivation, thespecific heat correctly shows a peak and not a divergence as it happens, instead, for the Ising model in d = 2 at ε (1) c = √
2. For this particular value of the energy density our numerical procedure correctly produces a finite valueof c although the numerical convergence is more delicate as highlighted by the scattered data present in Figs. 6 and 8for ε ≃ ε (1) c . In case of Fig. 8 the agreement is also quantitative as far as the location of the peak of the specific heatis concerned.As stressed throughout the paper, the concepts presented here are valid in principle for any O( n ) model in anyspatial dimensions. Hence the generalization of the calculations carried out here for the XY model on a square latticeto other O( n ) models in d dimensions should be possible, with possibly some technical complications, and may give ahint towards the development of approximation techniques that may be valid for estimating the density of states oflarge classes of Hamiltonian systems. Appendix A: “Ansatz-based” approximation with harmonic expansion
In this Appendix, we compute exactly the integrals in Eq. (48) under the assumption ε ¯ p = ε . We observe that thecomputation in this appendix can be extended also to the case ε p = ε ; however, we do not report the details of thismore general case here, as we do not need it for the development of the “ansatz-based” approximation.Inserting the change of measure in Eq. (48) due to the change of variables in Eq. (49), and setting ε p = ε , we have g (2) har ( ε, ¯ p ) = C ( N + ) C ( N − ) p N D (¯ p ) Z ∞ dY dX Y N + − X N − − δ (cid:0) Y − X (cid:1) Θ (cid:0) − Y + N + α (cid:1) Θ (cid:0) − X + N − β (cid:1) , (A1)where C ( k ) is the surface of the unitary ( k − S k − : C ( k ) = kπ k Γ (cid:0) k + 1 (cid:1) . (A2)Using the delta functional to express X as a function of Y and with simple computations, we get g (2) har ( ε, ¯ p ) = C ( N + ) C ( N − ) p N D (¯ p ) Z ∞ dY Y N − Θ (cid:0) − Y + N + α (cid:1) Θ (cid:0) − Y + N − β (cid:1) == C ( N + ) C ( N − ) p N D (¯ p ) Z √ N n − β dY Y N − = (A3)= 1 N − C ( N + ) C ( N − ) p N D (¯ p ) ( N n − β ) N − . To pass from the first to the second line we assumed that α ≫ β , as discussed in Section V C, so that min { n + α, n − β } = n − β . The above expression for g (2) har can be simplified in the large N limit by retaining only those factors that areexponential in N . We thus obtain g (2) har ( ε, ¯ p ) ∼ N ≫ exp (cid:26) N (cid:20) − LD (¯ p ) + log π − n + (¯ p ) log n + (¯ p )2 − n − (¯ p ) log n − (¯ p )2 + log ( n − (¯ p )) + log ( β ) (cid:21)(cid:27) , (A4)8where we have used the classical asymptotic expansion of the Gamma function. Inserting the last expression in theansatz (11), making use of Eq. (50), we obtain the approximate density of states given in Eq. (51). [1] D. J. Wales, Energy Landscapes (Cambridge University Press, Cambridge, 2004).[2] P. G. Debenedetti and F. H. Stillinger, Nature, , 259 (2001).[3] F. Sciortino, Journal of Statistical Mechanics: Theory and Experiment, , P05015 (2005).[4] J. N. Onuchic, Z. Luthey-Schulten, and P. G. Wolynes, Annual Review of Physical Chemistry, , 545 (1997).[5] L. N. Mazzoni and L. Casetti, Phys. Rev. Lett., , 218104 (2006).[6] L. N. Mazzoni and L. Casetti, Phys. Rev. E, , 051917 (2008).[7] M. Kastner, Physica A: Statistical Mechanics and its Applications, , 447 (2006).[8] H. Federer, Geometric Measure Theory (Springer, New York, 1969).[9] L. Casetti and M. Kastner, Phys. Rev. Lett., , 100602 (2006).[10] M. Kastner, O. Schnetz, and S. Schreiber, Journal of Statistical Mechanics: Theory and Experiment, , P04025 (2008).[11] L. Casetti, M. Kastner, and R. Nerattini, Journal of Statistical Mechanics: Theory and Experiment, , P07036 (2009).[12] L. Casetti, M. Pettini, and E. G. D. Cohen, Physics Reports, , 237 (2000).[13] M. Pettini, Geometry and topology in Hamiltonian dynamics and statistical mechanics (Springer, New York, 2007).[14] M. Kastner, Rev. Mod. Phys., , 167 (2008).[15] L. Casetti, C. Nardini, and R. Nerattini, Phys. Rev. Lett., , 057208 (2011).[16] In this form, Eq. (3) cannot be exact neither for finite systems nor in the thermodynamic limit [18] because it wouldimply wrong values of the specific heat critical exponent α ( n ) , that would be given by α ( n ) = − α (1) for any n . However,it correctly reproduces the sign of the critical exponent that entails a cusp-like behavior of the specific heat at criticalityrather than a divergence, as happens in the Ising case.[17] A. Campa, A. Giansanti, and D. Moroni, Journal of Physics A: Mathematical and General, , 6897 (2003).[18] C. Nardini, R. Nerattini, and L. Casetti, Journal of Statistical Mechanics: Theory and Experiment, , P02007 (2012).[19] As to a possible relation between the BKT transition and the 2- d Ising transition, it is interesting to notice that they doshare some “weak” universality as defined by Suzuki [37], despite being very different transitions. The critical exponentratio β/ν and the exponent δ , for instance, are equal in the two universality classes [38].[20] R. Nerattini and L. Casetti, in preparation (2013).[21] D. Mehta and M. Kastner, Annals of Physics, , 1425 (2011).[22] D. Mehta, Phys. Rev. E, , 025702 (2011).[23] These considerations are valid also for O( n ) models defined on generic graphs and with a generic interaction matrix J ij inEq. (4).[24] R. Nerattini, M. Kastner, D. Mehta, and L. Casetti, Phys. Rev. E, , 032140 (2013).[25] R. Schilling, Physica D: Nonlinear Phenomena, , 157 (2006).[26] R. Nerattini, Paesaggio energetico dei modelli XY , Master’s thesis, Universit`a di Firenze (2010).[27] Their rˆole reminds the one played by N π and N d in [18] for the mean-field and the one-dimensional XY models, respectively.[28] M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical Functions (Dover, New York, 1965).[29] The number of couples of opposite spins will be simply given by N D = N − N c .[30] Web page: http://alps.comp-phys.org/ .[31] J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics, , 1181 (1973).[32] V. L. Bereˇzinskij, Sov. Phys. JETP, , 493 (1971).[33] More precisely, the quantity that would appear is D ( p ) N and will be referred to as the “reduced determinant” of theHessian matrix of H ( n ) in p ; the index I of a stationary point p is the number of negative eigenvalues of the Hessian matrixof the Hamiltonian evaluated at p .[34] This behavior is analogous to the one obtained in [26] with a similar procedure applied to the partition function of thesystem.[35] F. H. Stillinger and T. A. Weber, Science, , 983 (1984).[36] A different choice of the U ¯ p given by a single ball, may look more natural at first sight. This choice, however, lead to theinconsistent result that if N − = 0, the minimum energy density value available to the system in the set U ¯ p is lower than ε ¯ p .[37] M. Suzuki, Progress of Theoretical Physics, , 1992 (1974).[38] P. Archambault, S. T. Bramwell, and P. C. W. Holdsworth, Journal of Physics A: Mathematical and General,30