Phase transitions at high energy vindicate negative microcanonical temperature

# Phase transitions at high energy vindicate negative microcanonical temperature

P. Buonsante QSTAR & CNR - Istituto Nazionale di Ottica, Largo Enrico Fermi 2, I-50125 Firenze, Italy.    R. Franzosi QSTAR & CNR - Istituto Nazionale di Ottica, Largo Enrico Fermi 2, I-50125 Firenze, Italy.    A. Smerzi QSTAR & CNR - Istituto Nazionale di Ottica, Largo Enrico Fermi 2, I-50125 Firenze, Italy.
July 25, 2019
###### Abstract

The notion of negative absolute temperature emerges naturally from Boltzmann’s definition of “surface” microcanonical entropy in isolated systems with a bounded energy density. Recently, the well-posedness of such construct has been challenged, on account that only the Gibbs “volume” entropy —and the strictly positive temperature thereof— would give rise to a consistent thermodynamics. Here we present analytical and numerical evidence that Boltzmann microcanonical entropy provides a consistent thermometry for both signs of the temperature. In particular, we show that Boltzmann (negative) temperature allows the description of phase transitions occurring at high energy densities, at variance with Gibbs temperature.

Our results apply to nonlinear lattice models standardly employed to describe the propagation of light in arrays of coupled waveguides and the dynamics of ultracold gases trapped in optical lattices. Optically induced photonic lattices, characterized by saturable nonlinearity, are particularly appealing because they offer the possibility of observing states —and phase transitions— at both signs of the temperature.

## I Introduction

Since the pioneering work of Purcell, Pound and Ramsey on nuclear spin systems Ramsey and Pound (1951); *Purcell_PR_81_279; Ramsey (1956), negative absolute temperature has been an established concept in Statistical Physics Landau and Lifschitz (1980); Kittel and Kroemer (1980). The ever growing control of ultracold atoms recently allowed the preparation of negative-temperature states for motional degrees of freedom of a bosonic gas loaded in an optical lattice Braun et al. (2013). Despite this remarkable result, the very notion of negative temperature has been challenged in a recent article Dunkel and Hilbert (2013), on account that the Boltzmann “surface” entropy it stems from would be inconsistent, and the only consistent picture would be based on the Gibbs “volume” entropy. This criticism spurred a lively and ongoing debate Sokolov (2014); Frenkel and Warren (2015); Dunkel and Hilbert (2014); Vilar and Rubi (2014); Treumann and Baumjohann (2014); Schneider et al. (2014); Hilbert et al. (2014); Dunkel and Hilbert (2014); Treumann and Baumjohann (2014); Swendsen and Wang (2014); Campisi (2015); Ferrari (2015); Cerino et al. (2015); Wang (2015); Hänggi et al. (2015); Poulter (2015).

Here we address the consistency of Boltzmann microcanonical temperature, focusing on a class of nonlinear lattice models standardly employed in the description of the propagation of light through arrays of waveguides, and the dynamics of ultracold bosons trapped in optical lattices. Our most interesting results apply to the case of optically induced photonic crystals Efremidis et al. (2002); Fleischer et al. (2003, 2003); Melvin et al. (2006). We find that, at variance with standard nonlinearity, the saturable nonlinearity characterizing these models supports both positive– and negative-temperature states in the same system. Even more interestingly, we show that the same physical system can undergo phase transitions for critical energies both in the lower and in the upper portion of the (bounded) energy spectrum. A finite-size scaling analysis shows that the Boltzmann picture provides a consistent description in both cases. While the former correspond to the standard situation, the critical temperature of the latter turns out to be finite and negative. Despite they typically manifest themselves in a clear ordering of the system, phase transitions at high energy densities are not consistently captured by the Gibbs picture. More in general, the Gibbs temperature does not appear to be a measurable quantity for the systems under concern.

As to the standard “cubic” nonlinearity typical of ultracold lattice bosons, we confirm that it does support negative-temperature states, as proposed in Refs. Mosk (2005); Rapp et al. (2010) and experimentally demonstrated in Refs. Braun et al. (2013, 2015). Owing to a pathological scaling of the energy density, states with positive or negative temperature turn out to be unstable for attractive or repulsive interactions, respectively. Nevertheless, phase transitions at negative temperatures should be observable in the former case, as heralded by the ordering phenomena observed experimentally Braun et al. (2013, 2015).

We support our conclusions with analytic arguments and extensive (microcanonical) numerical simulations. In particular, we provide independent and concordant tests of the fact that the considered dynamical states correspond to thermal equilibrium. We measure their (Boltzmann) temperature —which can be either positive or negative depending on the (conserved) energy density— either as a (time-averaged) function of the instantaneous configuration of the dynamical variables Franzosi (2011) or through a fit of the average distribution of the relevant modes of the system. Also, we show how, irrespective of the sign of the temperature, a large lattice acts as a thermostat for a small sublattice, thus confirming the equivalence between isolated and thermostated systems that is crucial for a consistent definition of temperature. From a different point of view, being able to support negative temperature states, the sublattice can be in principle used as a thermometer for the whole lattice.

The plan of this paper is the following. We briefly survey the concept of negative (Boltzmann) temperature in Sec. II, and introduce the nonlinear lattice models we focus on in Sec. III. After briefly addressing ensemble equivalence in Sec. IV, we discuss phase transitions —at both positive and negative temperatures— in Sec. V. The results of our numerical simulations are presented in Sec. VI. More detail about our results can be found in the appendix sections.

## Ii Negative Absolute Temperatures

In the microcanonical ensemble, the inverse temperature of the system is defined as

 β=1kB∂s∂h (1)

where is the entropy density corresponding to the energy density , and is Boltzmann’s constant. In general, two choices for are possible, corresponding to Boltzmann’s and Gibbs’ definitions. According to the former, , where is the number of degrees of freedom in the system, is the density of microstates at energy density and is a constant having the same dimension as . The Gibbs entropy is obtained by replacing with , i.e. the number of microstates having an energy density up to the chosen one. As it is well known, in the thermodynamic limit these two definitions are equivalent in “standard” systems lacking an upper bound to the energy density, for simple geometrical reasons (see e.g. Reichl (1998)). This comes about because is an increasing function of . However, some systems exist where the energy density has an upper bound and is a (non-negative) concave function featuring a maximum at some finite energy density . The logarithm of clearly has the same properties, which entails that is positive for , negative for , and vanishes at . A simple version of the lattice models introduced in Sec. III is employed to exemplify this behavior in Appendix A.

The key ingredient for the occurrence of negative Boltzmann temperatures is the existence of an upper bound to the available energy densities. Furthermore, the elements of the thermodynamical system must be in equilibrium, so that a temperature can be consistently defined. Finally, the system must be thermally isolated from any system that do not meet the previous requirements Ramsey (1956). In most cases the first condition fails because of the unboundedness of the kinetic energy term. The first experiments demonstrating negative Boltzmann temperatures Ramsey and Pound (1951); *Purcell_PR_81_279 involved spin systems, which are not affected by this problem owing to the lack of kinetic degrees of freedom. The kinetic energy density of a gas of particles can be effectively bounded from both above and below in the presence of a periodic potential inducing an energy gap sufficiently large that the physics of the system is dominated by states in the lowest energy band. Building on this observation, negative-temperature states for motional degrees of freedom in an ultracold bosonic gas have been demonstrated in a recent experiment Braun et al. (2013).

However, as observed in Ref Dunkel and Hilbert (2013), ensures that is a non-decreasing function of , and hence the Gibbs temperature is non negative even in systems with a bounded energy density. Therefore, plugging the Boltzmann or Gibbs entropy in Eq. (1) can produce very different temperatures, and the question arises as to which picture is the correct one.

The appearance of a negative sign in an absolute temperature might look disturbing enough to discard the Boltzmann framework at first glance. It should be noted, however, that a state at negative temperature is not colder than “the coldest possible state”, i.e. the ground state. Since its energy density exceeds , it is in fact hotter than the state attaining the maximum Boltzmann entropy, which has an infinite temperature. The fact that the concept of “hotter that ” sounds still somewhat disturbing can be merely ascribed to the traditional use of in the scale of temperatures. Using instead of restores the “correct order” of cold and hot in the whole range of Boltzmann temperatures111Here we do not introduce a minus sign in front of , in order to avoid confusion. Ramsey (1956).

In the following, we produce analytical and numerical evidence that the Boltzmann entropy does in fact provide a consistent thermodynamic picture. As we mention in the introduction, Ref. Dunkel and Hilbert (2013) advocates that only the Gibbs entropy results in a consistent thermostatics, and dismisses all previous claims about negative absolute temperatures. These arguments are further elaborated in Refs. Hilbert et al. (2014); Hänggi et al. (2015). While we refer the Reader to a different publication con () for a more systematic discussion of the points raised in Refs. Dunkel and Hilbert (2013); Hilbert et al. (2014); Hänggi et al. (2015), in the following we occasionally comment on some of them.

We start by observing that disturbing features also lurk behind the instinctively appealing non-negative temperatures characterizing the Gibbs formalism. Indeed, for lattice systems such as the ones we are going to address shortly, the Gibbs temperature corresponding to energy densities increases arbitrarily with the number of degrees of freedom in the system. Hence, in the thermodynamic limit, it is identically infinite on the whole finite interval of energy densities exceeding the one attaining the maximum Boltzmann entropy222This is the main criticisms leveled by Ref. Vilar and Rubi (2014) agains the Gibbs picture advocated by Ref. Dunkel and Hilbert (2013). We note that this weird feature of the Gibbs temperature is mentioned in Ref. Dunkel and Hilbert (2013) itself, albeit only in the Supplementary Information. . Correspondingly, the Gibbs heat capacity would be identically zero. We observe that, while this might be to some extent internally consistent, it presents at least two problems, as we illustrate in more detail in Sec. V and Appendix A. In the first place, it clearly makes the Gibbs temperature incapable of describing phenomena involving states with . Also, the Gibbs temperature cannot be measured as a microcanonical average by exploiting the usual formulation of the equipartition theorem advocated in Refs. Dunkel and Hilbert (2013); Hilbert et al. (2014).

We remark that the convexity properties of the density of states of the systems under investigation are a typical consequence of the large number of degrees of freedom and short-range interactions. The oscillating density of states discussed in Refs. Hilbert et al. (2014); Hänggi et al. (2015) can crop up in systems with a small number of degrees of freedom or long-range interactions. Standard thermodynamic relations can be —at least formally— used in such cases, although at the expense of features that appear crucial for a sensible definition of temperature Cerino et al. (2015). For instance, the equivalence of isolated and thermostated systems is lost. Also, when joined, two systems having the same temperature could equilibrate to a completely different temperature. The well-posedness and usefulness of the concept of temperature in such situations is at least arguable.

We finally note that some of the arguments against the existence of negative-temperature equilibrium states are based on the failure of one or more of the key conditions listed above Ramsey (1956). For instance, it has been argued that such states would not be stable if the system in which they occur is brought into contact with a system that is unable to sustain negative-temperature states, e.g. due to the lack of an upper bound to the energy density Dunkel and Hilbert (2013); Hilbert et al. (2014); Romero-Rochín (2013) . This, however, is basically a truism, since the composite system clearly does not meet the requirements Ramsey (1956) for sustaining negative-temperature equilibrium states.

The model we are going to address in the following has no doubt been devised as an extreme idealization of real physical systems, as remarked in Refs. Hilbert et al. (2014); Hänggi et al. (2015). On the other hand, over the last few years enormous steps towards the faithful experimental simulation of similarly ideal models have been made, the trailblazer being ultracold-atom physics Jaksch and Zoller (2005); Bloch et al. (2012). Lattice Hamiltonians originally conceived as rough yet extremely challenging toy models have been experimentally realized with an high degree of fidelity by loading a cold atomic gas into a “crystal of light”. As we mention, evidence of negative-temperature states has already been reported for these systems Braun et al. (2013). The observation of phase transitions occurring at negative critical temperature, such as the ones we discuss in Secs. V and VI involves a similar effort towards the faithful experimental simulation of the single-band lattice model described in the following Section. On a related note, we mention that a Mott insulator-superfluid quantum phase transition has been observed at negative temperature in the experimental realization of a Bose-Hubbard model with attractive interactions Braun et al. (2015). We once again remark that significant steps in the realization of synthetic nonlinear lattice models have been likewise made in the field of optical waveguides. These include an analysis of the effects of nonlinearity the Anderson localization Lahini et al. (2008) and the observation of a Berezinskii-Kosterlitz-Thouless phase transition at positive critical temperature Guohai and Fleischer (2012) in optical systems obeying the discrete nonlinear Schrödinger (DNLS) equation.

## Iii The model

We focus on nonlinear lattice models Kevrekidis (2009) of the form

 H=U∑ru(|zr|2)−J∑rr′z∗rArr′zr′ (2)

where denotes a site in a dimensional (D) lattice and is the relevant coordination matrix. The coordinates of the sites are integer numbers, , so that the total number of sites in the lattice is . Periodic boundary conditions are assumed, i.e. , where is the versor along the -th direction. We set the hopping amplitude to , so that the units for energy and time are and , respectively. As to the nonlinear term, we consider two cases

 u1(n)=−log(1+n),u2(n)=12n2. (3)

The former corresponds to the saturable nonlinearity typical of the equations describing the propagation of a light probe in an optically induced nonlinear photonic lattice Efremidis et al. (2002); Fleischer et al. (2003, 2003); Melvin et al. (2006), while the latter produces the cubic nonlinearity of standard DNLS equations. These are employed in the description of diverse phenomena Kevrekidis (2009); Lederer et al. (2008), including the dynamics of ultracold atoms loaded in optical lattices Cataliotti et al. (2001); Trombettoni and Smerzi (2001); Smerzi et al. (2002); Cataliotti et al. (2003); Trombettoni et al. (2005); Flach and Gorbach (2008); Pikovsky and Shepelyansky (2008); Kopidakis et al. (2008); Flach et al. (2009); Rumpf (2004); Small et al. (2011); Iubini et al. (2013); Franzosi et al. (2011) and the propagation of light in waveguide arrays Morandotti et al. (1999); *Morandotti_PRL_86_3296; Rumpf (2004); Small et al. (2011); Iubini et al. (2013); Franzosi et al. (2011); Lederer et al. (2008); Christodoulides et al. (2003); Lahini et al. (2008); *Lahini_PRL_103_013901.

The equations of motion generated by Hamiltonian (2) via the Poisson brackets have two first integrals, the energy and “particle” density,

 h=V−1H,a=V−1∑r|zr|2. (4)

The presence of a conserved quantity other than the energy density is important for the occurrence of negative temperatures, because it makes the configuration space of any finite system compact.

In the non interacting limit the Hamiltonian becomes linear, and the (thermo)dynamics is described exactly by the single-particle “plane-wave” eigenmodes , where is the quasimomentum along direction . The corresponding single-particle energies form a band bounded by .

It is easy to check that the “plane-wave” states are normal modes for the nonlinear equations as well, provided that the single-particle energy is replaced by the frequency . The corresponding energy density is . For repulsive interactions, —i.e. defocusing nonlinearity— the energy densities are bounded from below by . On finite lattices the energy densities also have an upper-bound, which however diverges in the thermodynamic limit for the standard nonlinearity, . Note indeed that the energy density of a state where only an individual site is occupied is . This means that negative-temperature equilibrium states are problematic for the standard defocusing nonlinearity Rumpf (2004), although metastable states at can persist for astronomically long times on 1D lattices Iubini et al. (2013); Franzosi et al. (2011).

In the case of saturable nonlinearity, , the upper bound of the energy density remains finite in the thermodynamic limit, and tends to . That is, the energy per particle is of the order of the maximum single-particle energy.

The situation for attractive interactions , —i.e. for self-focusing nonlinearity— is related to the previous case through the mapping

 (5)

This means that negative temperatures are well defined for the standard nonlinearity, , and positive ones are problematic. Note that switching the interaction strength to negative values is a crucial step for obtaining negative-temperature states in a bosonic gas loaded in an optical lattice Braun et al. (2013). The nonlinearity is typically self-focusing also in the case of waveguide arrays. Although the sign of can be reversed in photorefractive crystals Efremidis et al. (2002); Christodoulides et al. (2003), the defocusing case does not seem to lend itself to a tight-binding approach in current experimental realizations Efremidis et al. (2002). For these reasons, unless otherwise specified, in the following we fix our attention mainly on the self-focusing case, .

## Iv Ensembles

In general, the details of an experimental system determine the most appropriate choice for the statistical ensemble to be adopted in the description of its thermodynamic properties. The natural choice for an isolated system is the microcanonical ensemble. For ergodic systems, one expects that microcanonical thermodynamic quantities can be equivalently obtained as ensemble or temporal averages.

In the canonical and grand canonical ensembles, is a Lagrange multiplier fixing the total energy. It is possible to prove that in the thermodynamic limit this multiplier coincides with the microcanonical definition of temperature, Eq. (1) con (). Also, canonical —or, in the presence of additional conserved quantities, grand canonical— time averages can be obtained by considering a sufficiently macroscopic subset of an isolated, microcanonical system. In the absence of pathologies, one expects that the subsystem has the same thermodynamical properties as the whole system. More in general, one expects that all ensembles provide a consistent description of a given system.

In the non-interacting limit, , a detailed analysis of the statistical ensembles for the model in Eq. (2) can be carried out, at both the semiclassical and quantum level. In both cases, the relation between the energy density and the inverse temperature turns out to be the same for all the three ensembles con (). The easiest way of obtaining such relation is through the grand-canonical ensemble. The grand partition function for the model in Eq. (2) is

 Q=∫∏rdzre−βV[h({zr})−μa({zr})], (6)

where is the chemical potential, i.e. the Lagrange multiplier selecting the average density , and we omitted the dependance of the energy density on the parameters. In the non-interacting limit, the integral in Eq. (6) can be easily carried out. It is likewise easy to obtain the average occupation of the single particle modes333This is the “classical version” of the Bose-Einstein distribution , that comes about because the occupation of the single-particle modes of Eq. (2) is not restricted to integers.

 nq(β,μ)=⟨|~zq|2⟩=1β1εq−μ, (7)

where is the Fourier transform of the configuration of the system. For fixed , , and —where denotes the kinetic energy density per particle—, the chemical potential can be found by inverting the relations

 a=∑qnq,h=∑qεqnq. (8)

On a sufficiently large one-dimensional lattice this calculation can be carried out analytically, and gives

 β=−1a2κ4−κ2,μ=κ2+42κ. (9)

Thus, for an equilibrium thermodynamic state, if , and if .

As we mention, the first of Eqs. (9) accurately describes the relations between , and and that are found in the canonical and microcanonical ensembles con (). This means that, in the limit of a large number of sites , . Using the Laplace method it is possible to evaluate , and the Gibbs entropy thereof. As illustrated in Appendix A, it turns that the Gibbs inverse temperature is the same as in the Boltzmann picture in the lower energy interval, and vanishes identically on the whole upper interval. This causes the failure of the standard equipartition theorem, which underlies the measurability of the Gibbs temperature as a microcanonical average (see Sec. VI ).

The thermodynamics of model444Our choice for the Poisson brackets corresponds to the standard bosonic commutation rules when the C-number is interpreted as the expectation value of an on-site boson operator, e.g. in the Bose Hubbard model. Refs K.Ø. Rasmussen et al. (2000); Samuelsen et al. (2013) make a different choice for the same Poisson brackets. The two choices are connected by a simple mapping. We illustrate the results of Refs K.Ø. Rasmussen et al. (2000); Samuelsen et al. (2013) in the light of our choice. (2) on one-dimensional lattices has been addressed in Refs. K.Ø. Rasmussen et al. (2000) and Samuelsen et al. (2013) for defocusing standard and saturable nonlinearity, respectively. There, the grand-canonical partition function, Eq. (6) is calculated using a transfer-matrix approach, which allows the identification of the the region in the plane corresponding to positive temperatures. This is bounded from below by the ground-state energy, , and from above by a critical line that, in the case of standard nonlinearity, assumes the simple form (the superscript in the energy densities here refers to the grand canonical inverse temperature). In Ref. K.Ø. Rasmussen et al. (2000) the region is argued to correspond to negative temperatures, based on the change in the concavity of the probability distribution function of the amplitudes , as obtained from microcanonical dynamical simulations. Ref. Samuelsen et al. (2013) repeats basically the same analysis as in Ref. K.Ø. Rasmussen et al. (2000), but the only new insight it provides about the region is the observation that initial states picked in that region end up having one single very mobile localized excitation. This is contrasted with the larger number of pinned localized excitations characterizing the standard nonlinearity Rumpf (2004); Iubini et al. (2013); K.Ø. Rasmussen et al. (2000).

We find that, in view of the boundedness of the available energy densities, the transfer-matrix approach applies also for for the defocusing saturable nonlinearity considered in Ref. Samuelsen et al. (2013). Specifically it can be applied for , where in general depends on and . This is illustrated in Fig. 1, where we analyze the relation between between and the kinetic energy per particle in the interacting case, as provided by the transfer matrix approach. In particular, it is clear from panel a), that solutions exist at negative for the saturable nonlinearity. The leftmost symbol of each kind marks the largest negative we were able to analyze for the corresponding parameter choice. The failure of the transfer-matrix approach for larger negative ’s is related to a phase transition between an extended and a localized state, occurring at a finite negative . Indeed, as we discuss in the following, states with555We have verified that the energy thresholds K.Ø. Rasmussen et al. (2000); Samuelsen et al. (2013) also apply to two and three dimensional lattices ext (). exhibit a persistent and mobile localized excitation only for sufficiently large energies, i.e. for sufficiently small negative temperatures ext (). Panel b) in Fig. 1 illustrates similar results for the case of standard nonlinearities, where the transfer matrix approach clearly fails as soon as K.Ø. Rasmussen et al. (2000). It is interesting to note that, despite the non-negligible effective interaction, the data points closely follow the analytical relation derived in the non-interacting case, Eq. (9). Our microcanonical simulations confirm that this is the case also for two– and three–dimensional lattices ext ().

We once again remark that the above described situation is reversed for self-focusing nonlinearities.

## V Phase Transitions

In the non-interacting limit, the model in Eq. (2) is known to undergo Bose-Einstein condensation on a three-dimensional lattice. If , that is if the energy density is sufficiently small, a macroscopic fraction of the particle density occupies the ground state of the system, i.e. the kinetic-energy eigenstate corresponding to quasimomentum . On a finite-size lattice this second-order phase transition manifests itself as a crossover. For this transition to occur, it is crucial that the density of states in the vicinity of the ground state has a suitable behavior. Since this behavior is literally mirrored by the density of states in the vicinity of the highest-energy state666We remark that this symmetry in the system spectrum and in the relevant Boltzmann entropy does not imply the physical or thermodynamical equivalence of states having opposite energy density, as argued in Refs. Dunkel and Hilbert (2014); Hilbert et al. (2014). In fact, for one, the states are distinguished by the sign of the derivative of the Boltzmann entropy, and of the temperature thereof., it seems fair to expect that the system condenses into the highest-energy eigenstate for sufficiently large energy densities —i.e. at small negative temperatures—. Specifically, one expects that for a macroscopic density of particles occupies the state . In fact, this is what results from a simple grand-canonical calculation (see Fig. 5 in Appendix B). We mention that clear signatures of phase transitions at high-energies have been recently discussed for short-range ferromagnets Kastner and Pleimling (2009). Also, the emergence of order at small negative temperature in an isolated planar superfluid has been recently discussed in Ref. Simula et al. (2014), thus confirming an early prediction by L. Onsager Onsager (1949).

The condensation transition occurring at large is expected to survive the introduction of a defocusing standard nonlinear term, (see e.g. Ref. Ramakumar and Das (2005)). As we observed earlier, such nonlinearity has a dramatic effect on high-energy states. The upper bound to the energy density, and the negative temperatures thereof, are lost in the thermodynamic limit. In view of the mapping in Eq. (5), a condensation into the highest-energy state is expected at a negative critical temperature in the case of self-focusing standard nonlinearity.

The saturable nonlinearity, , produces less dramatic effects. It is not hard to check that for the ground state is localized. Specifically, the corresponding particle density features a single peak of finite width (corresponding to a breather), on top of a uniform background777This self-trapped state can be centered at any lattice site, and hence it is not unique. This simmetry breaking is a well known feature in nonlinear systems. . The highest-energy state is instead extended, and coincides with the uniform “plane-wave” state with . It is therefore tempting to envisage two phase transitions for this system: a condensation into the extended highest-energy state for , and a self-trapping transition —i.e. a “condensation” into the localized ground-state— for . In the defocusing case the localization properties of the extremal states are swapped, and hence one expects a condensation into an extended ground-state for , and a condensation into a localized highest-energy state for . According to the Mermin-Wagner theorem, condensation into an uniform state is not expected to occur for at a finite critical temperature, since it involves the breaking of the continuous symmetry in the phases of the dynamical variables . The localized ground-state instead breaks the discrete translational symmetry of the lattice, and does not exhibit long-range order. Therefore, the corresponding localization transition is not excluded for .

Also, the model in Eq. (2) is strictly related to the classical XY model, and it is therefore expected to undergo a Berezinskii-Kosterlitz-Thouless (BKT) transition at a finite temperature on a 2D lattice Trombettoni et al. (2005); Small et al. (2011). Thus a particularly intriguing scenario opens up for two-dimensional optically induced nonlinear photonic lattices, which realize model (2) for self-focusing saturable nonlinearity Fleischer et al. (2003); Efremidis et al. (2002). One could observe a localization transition at finite positive temperatures, and a BKT transition at finite negative temperatures. In fact, signatures of the latter transition have been reported for defocusing nonlinearity at positive temperatures Guohai and Fleischer (2012).

## Vi Thermalization and Thermometry

The numerical integration of the dynamical equations generated by Hamiltonian (2) reveals that, after a possibly long transient, the system reaches a stationary state in which the instantaneous value of observables characterizing the system performs small oscillations about an asymptotic value (see Appendix D). The observables we typically consider are for instance the kinetic and interaction energy per particle,

 κ =−1Va∑r,r′z∗rAr,r′zr′=1a∑q|~zq|2εq (10) I =a−1h−κ=UVa∑ru(|zr|2) (11)

where denotes the Fourier transform of . More importantly and interestingly, it is possible to give an estimate of the instantaneous microcanonical Boltzmann temperature as a function of the instantaneous configuration of the system Rugh (2001); Davis et al. (2002); Franzosi (2011) (see Appendix D for details). A time average excluding the initial transient,

 ⟨O⟩=1Δt∫t0+Δtt0dt′O({zr(t′)}) (12)

provides a “measure” of the generic observable , where we stress once again that this includes the Boltzmann temperature. Plotting vs reveals that the relation between these quantities remains remarkably close to the one applying in the non-interacting limit also for non-negligible nonlinearity, although some small deviations appear for the symmetry-broken phases ext ().

An additional evidence of thermalization is provided by the observation that the prediction in Eq. (7) is fulfilled by the “relevant modes” in the system, which we generically denote . That is, a plot of vs. the corresponding energies results into a straight line whose slope coincides with the time-averaged measure of the microcanonical Boltzmann temperature . In the absence of condensation the relevant modes are the single particle modes, i.e and . In the condensed phase the linear relation is fulfilled by the Bogoliubov quasiparticle modes (see Appendix D for details. Figures 7 and 8 contain examples of this behavior). Therefore, the interactions not only drive the system to equilibrium but, when this is reached, maintain it by acting as a “heat bath” for the relevant, effectively non-interacting modes of the system. One could argue that the above slope represents a sort of canonical or (even grand-canonical) measure of the temperature, since, while the total population of the modes may be strictly conserved (in the case of the kinetic modes), the relevant total energy is not. A grand-canonical description is also obtained by considering only the dynamical variables belonging to a sufficiently large sublattice of the whole lattice. Indeed, the first integrals of the motion, and , are not conserved when restricted to a portion of the whole sample. The fact that the relevant modes of the sublattice behave as those in the whole lattice (as apparent in Fig. 7 in Appendix D), is a further proof that the system is in equilibrium, and that the Boltzmann temperature has the properties expected of a temperature. If the whole lattice is much larger than the sublattice, then the former acts as a thermostat (and a “chemostat”) for the latter. We also checked that a grand-canonical Langevin approach generalizing the one introduced in Ref. Iubini et al. (2013), where is an external parameter, produces results in agreement with the above described observations for both positive and negative temperatures ext ().

A further fundamental test for a well defined temperature concerns the equilibration of two systems that are separately at equilibrium at different temperatures. One expects that, when these are brought into contact, energy —and, in our case, particles— flows in agreement with the intuitive notion of cold and hot, in such a way that eventually the inverse temperature of the composite system is intermediate between the two initial values. As discussed in Ref. con (), this is exactly the case for the Boltzmann temperature, irrespective of the sign it initially has in the separate systems. Note that some care must be taken when using a small system as a thermometer for a larger system. Indeed, both energy and particles will be in general exchanged to attain equilibrium con (); ext (). Of course, in a thermalization experiment where two systems with opposite signs of the temperature are brought into contact, the resulting composite system should support both positive– and negative–temperature states. For instance, if one of the two initially separate systems only supports positive temperatures and the other supports both, there is no chance that the equilibrium state of the composite system be negative, irrespective of the initial temperature of the latter subsystem Ramsey (1956). In this respect, in view of its ability to support both positive– and negative–temperature states, a nonlinear lattice system characterized by saturable nonlinearity seems to be an ideal setting for this kind of experiments, at least in principle.

The descriptive power of the Boltzmann microcanonical temperature becomes fully evident in the presence of phase transitions, as demonstrated in Figs. 2 and 3. Figure 2 refers to the 2D lattice system with self-focusing saturable nonlinearity modeling the propagation of light through an optically induced photonic lattice Efremidis et al. (2002); Fleischer et al. (2003, 2003). As we anticipated in Sec. V, two phase transitions occur in such a system. A suitably defined exponent , sensitive to the decay properties of the radial correlations, signals a BKT transition in the system Small et al. (2011) (see Appendix C). The data we obtain for different lattice sizes shown in Fig. 2 a), strongly suggest a transition at negative critical temperature, and cross at the value expected at the critical point Small et al. (2011); Polkovnikov et al. (2006) (see Appendix C for more detail). As we observed earlier, for sufficiently large the microcanonical state of the system develops a density peak. This suggests that the system is partially condensed into its ground-state, which is also characterized by a (taller) density peak. We quantify this condensation through the time-averaged projection of the microcanonical state of the system onto the ground-state888The occurrence of a density peak breaks the (discrete) symmetry of the lattice. Thus, we calculate the overlap of the instantaneous dynamical state and the ground-state after centering the relevant peaks at the same lattice site. . As illustrated by panel b) of Fig. 2, a plot of such quantity against the Boltzmann temperature of the corresponding microcanonical state clearly signals the occurrence of a phase transition.

Finally, Fig. 3 shows the average occupation of the highest-energy state in a three dimensional lattice system with standard self-focusing nonlinearity. The data points, obtained as temporal microcanonical averages according to Eq. (12), clearly signal a (condensation) phase transition, and nicely follow the (solid gray) curve obtained from a grand-canonical calculation based on the Bogoliubov approximation (see Appendix D). The same behavior is obtained for the saturable nonlinearity on a three dimensional lattice. In the self-focusing case one observes the condensation into the localized ground state at small positive temperatures, and the condensation into the extended highest-energy state state for small negative temperatures. Again, the relevant critical temperatures are finite ext ().

In the same physical situations the limitations of the Gibbs picture become evident. Although the microcanonical dynamics takes place in the phase-space “sheet” corresponding to a given value of the energy density (and, possibly, of other conserved quantities), the Gibbs entropy requires information about all energies below such value. It is therefore not immediately clear how the Gibbs (inverse) temperature could be obtained as a microcanonical ensemble average (and hence a time average). It is often argued Dunkel and Hilbert (2013); Hilbert et al. (2014) that this is made possible by the “equipartition theorem”

 β−1G=⟨ζj∂H∂ζj⟩, (13)

where denotes any element of the set of dynamical variables describing the microstate of the system, and the angle brackets denote the standard microcanonical average. However, Eq. (13) typically fails for system admitting negative Boltzmann temperatures. For instance, for a negative-temperature state of the self-focusing case of Eq. (2), the r.h.s. of Eq. (13) with is of the order of . As we have discussed, the corresponding tends instead to vanish as the system size increases, so that Eq. (13) cannot be possibly satisfied, since its l.h.s. diverges. This failure of the “equipartition theorem” stems from ignoring a surface term in the derivation of Eq. (13), which is legitimate only for systems that do not admit negative Boltzmann temperatures con (); ext (); Cerino et al. (2015) (where the Boltzmann and Gibbs pictures are equivalent, as we discussed before). Even if were actually measurable, its usefulness would be of very limited value in describing the phenomena involving the upper part of the spectrum, such as the phase transitions illustrated in Figs. 2 and 3. Indeed, while the critical energy density (and Boltzmann temperature) for such transitions becomes size-independent for sufficiently large systems, the corresponding Gibbs temperature indefinitely increases with the system size.

## Vii Discussion

We have addressed the statistical physics of nonlinear lattice models relevant in the description of the propagation of light in nonlinear media and the dynamics of ultracold atoms trapped in optical lattices. We have discussed how the Boltzmann picture provides a consistent description of equilibrium and equilibration processes, and that the absolute temperature characterizing the equilibrium states can have either sign. Negative temperatures correspond to high energy states, and come about due to the presence of an upper bound to the available energy density, and from the decreasing character of the entropy (density) in the vicinity of such bound. These features are already apparent in the noninteracting limit of the considered lattice models con (), and survive the introduction of interactions, provided that these do not give rise to pathological scaling in the thermodynamic limit. Interactions act as a heat bath for the relevant, effectively non interacting modes of the system, driving the system towards equilibrium. A large system can act as a thermostat for a smaller system, bringing it to a negative-temperature (i.e. higher-energy) state, provided that such a state can be supported by the composite system. Such a process might not fit the definition of “conventional heating” Dunkel and Hilbert (2013), but it is consistently described in terms of Boltzmann (inverse) temperature. A likewise consistent description is instead problematic in the Gibbs picture, where the whole upper interval of available energy densities corresponds to zero heat capacity and infinite temperature (which, in addition, is not measurable as a microcanonical average through the standard equipartition theorem). For the same reason, the description of the phase transitions taking place in the system at high energy densities is similarly problematic. As illustrated in Figs. 2 and 3, Boltzmann temperatures provide a consistent description of such transitions.

Optical systems represent an ideal testbed for our conclusions. Ordering phenomena related to the self-trapping transition discussed above have been observed in one– Fleischer et al. (2003) and two-dimensional Fleischer et al. (2003) lattices. Signatures of a BKT transition have been observed in a two-dimensional optically induced photonic lattice Guohai and Fleischer (2012), although they have been analyzed in terms of an “equipartition” effective temperature Small et al. (2011). Optically induced nonlinear photonic lattices Efremidis et al. (2002); Fleischer et al. (2003, 2003); Guohai and Fleischer (2012) are particularly intriguing, in that the relevant saturable nonlinearity preserves both the upper and lower bound characterizing the corresponding (single-band) linear lattice model, allowing the exploration of both positive- and negative-temperature states in the same system. Even more interestigly, the realizability of two-dimensional lattices Fleischer et al. (2003); Guohai and Fleischer (2012) opens up the possibility of observing phase transitions with critical temperature of both signs in the same (synthetic) physical system. States in the upper portion of the kinetic energy band can be excited by suitably tilting the input beam, as described in Refs. Christodoulides et al. (2003); Lahini et al. (2008). Phase transitions on lower dimensional systems could be engineered through the introduction of suitable “on-site” or topological defects Stamper-Kurn et al. (1998); Burioni et al. (2000).

The measure of the instantaneous microcanonical temperature requires the knowledge of the instantaneous configuration of the field, . A more feasible measure is obtained through Eq. (7). Indeed, as discussed above and demonstrated in Figs. 7 and 8 in Appendix D, a linear fit of the inverse average mode occupation versus the corresponding energy provides an estimate of the temperature that is in remarkable agreement with the time-average of the instantaneous microcanonical value.

###### Acknowledgements.
P.B. gratefully acknowledges helpful suggestions by R. Burioni, A. Vezzani and S. Wimberger.

## Appendix A Boltzmann Entropy

Figure 4 illustrates the concepts discussed in Sec. II, in the case of a uniform 1D lattice model comprising sites and containing non-interacting bosons, (, , ). The circles in the main panel correspond to the analytically calculated con () volume of the phase space relevant to the chosen energy and particle density

 ω(h,a)=∫∏qdzqδ(Lh−H)δ(La−N) (14)

where is the Hamiltonian in Eq. (2) and

 N=∑q|zq|2 (15)

Although the number of sites is not very large, the microcanonical result is very well described by the approximation derived from the grand canonical result in Eq. (9), where . The data in the inset refer to the Bose-Hubbard model obtained by changing the C-numbers and in Hamiltonian (2) into the lattice boson operators and , respectively. In view of the lack of interaction, a generic eigenstate of the system is a Fock state listing the number of bosons occupying each of the single-particle states. The behavior of can be therefore estimated by listing all the Fock states compatible with the chosen boson population and binning the corresponding energy densities. We divided the energy density interval into 183 bins but, for graphical reasons, plotted only a few of the corresponding numbers of microstates. The result is remarkably smooth owing to the very large number of Fock states (approximately ).

As discussed in Sec. II, the Boltzmann entropies are concave functions, and feature a maximum at . Therefore, according to Eq. (1) the corresponding Boltzmann inverse temperatures are positive for and negative for . Note that, owing to the small particle density, the classical and quantum results are quantitatively different, although qualitatively similar.

As we repeatedly mention, in the thermodynamic limit the Gibbs inverse temperature equals Boltzmann’s for and is identically zero for energy densities exceeding . This can be appreciated by plugging analytic function obtained from the grand canonical picture into

Using the Laplace method, for not too close to 0 we get

 Ω(h,a)≈⎧⎨⎩Ca22|h|[4−(a−1h)2]L+1,h∈[−2a,0)√πC22L+1a,h∈(0,2a], (17)

which, plugged into Eq. (1), gives

 βG(h,a)≈⎧⎨⎩L+1L−2a−2h4−(a−1h)2≈βB(h,a),h∈[−2a,0)0,h∈(0,2a].

where the subscripts in the inverse microcanonical temperatures refer to the Gibbs and Boltzmann pictures. Entirely similar results can be obtained numerically for higher dimensions or for the quantum case considered in the inset of Fig. 4.

## Appendix B Bose-Einstein condensation

As we mention in Sec. IV, the calculation of the grand partition function, Eq. (6), for the non-interacting version of the lattice model in Eq. (2) can be easily carried out analytically. The result is

 Q=∏qπβ(εq−μ) (18)

which gives rise to the average occupation distribution in Eq. (7). Despite the function in Eq. (7) is not the standard Bose-Einstein distribution, but rather its classical limit999See note 3., it still gives rise to condensation. The first of Eqs. (8) and Eq. (7) can be used to find the chemical potential corresponding to a given choice of the particle density and inverse temperature . Plugging the result into Eq. (7) gives the average occupation of each single-particle state. Note that the chosen can have either sign and, in view of the fact that , it must be for and for . This discontinuity in does not necessarily mean that the system undergoes a phase transition at , as argued in Ref. K.Ø. Rasmussen et al. (2000). Indeed, , so that . Thus the grand partition function is not singular for . The same is true in the presence of a non pathological interaction term, such as , as it is clear e.g. from the results of the transfer-matrix approach in Fig. 1 ext ().

The above sketched calculation can be easily carried out numerically. Figure 5 shows the behavior of the relative occupation of the extremal (kinetic) modes of a non-interacting discrete model on a 3D lattice, as the inverse temperature ranges from negative to positive temperatures. Expectedly, there exists a critical value above which the ground-state of the system is macroscopically occupied. As we discussed in Sec. V, owing to the symmetry of the energy spectrum, the highest-energy state is macroscopically occupied for .

In the presence of nonlinear interactions, only one of the extremal state maintains its extended character, while the other becomes localized. Since it corresponds to the breaking of a continuous (phase) symmetry, the condensation into the extended state occurs at finite only for , in agreement with the Mermin-Wagner theorem. As demonstrated in panel b) of Fig. 2, for saturable101010The analysis of the energy and temperature region in the vicinity of the extremal localized state is problematic in the case of the standard nonlinear term , because of the pathological scaling of the energy densities. nonlinearities the self-trapping transition occurs at finite also for .

The Bogoliubov approach allows the analysis of the condensation transition at finite interactions. In the following we sketch the simplest and most standard case, i.e. the condensation into the uniform ground state for defocusing nonlinearity and arbitrary interaction term. A more detailed and general calculation, including the stability of the excited plane-wave solutions, can be found elsewhere ext (). As usual, in view of the decoupling of quasimomenta expected in a uniform system, we consider a perturbation of the ground-state of the form

 zr(t)= +b∗qV∗qe−i(q⋅r−ϖqt)]}, (19)

and plug it into the equation of motion, retaining only the linear terms in and . The branch of the resulting excitation spectrum fulfilling the expected normalization relation, , corresponds to

 (20)

and

 (21)

When the Bogoliubov approximation applies, the interaction strength is incorporated into the spectrum in Eq. (20), and the system is governed by the effectively free Hamiltonian . One therefore expects that using in place of in Eq. (7), the average occupation of the Bogoliubov modes is obtained, i.e. . A procedure similar to the one sketched above for the noninteracting case allows to study the population of the macroscopically occupied (ground) state. This is how the solid curve in Fig. 3 has been obtained. Note that when the Bogoliubov spectrum coincides with the single-particle spectrum.

Fourier transforming the perturbed state in Eq. (19) we get

 ~zq =bqUqe−i(ε0+ϖq)t+b∗−qV∗−qe−i(ε0−ϖ−q)t (22)

and

 ⟨|~zq|2⟩=⟨|bq|2⟩(|Uq|2+|Vq|2) (23)

where we assumed that the time average of the terms containing the phase factors cancel out and that , on account that .

## Appendix C BKT transition

On 2D systems the BKT transition is signalled by a change in the decay properties of the radial correlations. Denoting , in the defocusing case one expects a power law decay, , for and an exponential decay, , for , with the decay exponent tending to as the critical point is approached Small et al. (2011); Polkovnikov et al. (2006). The quantity is a temperature-dependent correlation length. The quantity

 AΩ(β)=∫|r−r′|<√Ωdrdr′|Crr′|2∼Ω1+σβ (24)

where

 σβ={1−αββ>βBKT0β<βBKT

can be used as an indicator for the transition Small et al. (2011).

In view of the square modulus in the integrand of Eq. 24 one expect an entirely similar behavior in the self-focusing case,

 σβ={1−αββ<β′BKT0β>β′BKT

where .

This indeed is what we obtain in Fig. 2 a) for self-focusing saturable interactions. Note in particular that the exponent at the crossing point of the curves corresponding to different sizes is very close to the expected value .

We analyzed also the defocusing standard nonlinearity considered in Ref. Small et al. (2011), qualitatively confirming the findings therein discussed ext (). We recall that the temperature is expected to diverge as the energy density approaches the upper bound of the positive-temperature region. We observe that, conversely, the effective temperature defined in Ref. Small et al. (2011), tends to a finite value. For the considered standard nonlinearity we get .

## Appendix D Thermalization and Thermometry

On sufficiently ergodic systems, the microcanonical (Boltzmann) inverse temperature can be obtained as the time average of a suitable function of the dynamical variables. When the only first integral of the motion is the total energy, such function is related to the curvatures of the “energy sheet” involved in the dynamics Rugh (1997). Generalizing this approach to equations having one further first integral Franzosi (2011), the instantaneous microcanonical inverse temperature is obtained as

 β(t)=∥n∧h∥∇⋅¯v[∇⋅(¯v∥n∧h∥)−¯n⋅(¯n⋅∇)¯v∥n∧h∥] (25)

where is the gradient in the -dimensional Euclidean space of the real and imaginary parts of the complex dynamical variables and the boldface variables are vectors in the same space. Specifically

 v=¯h−(¯h⋅¯n)¯n,h=∇H,n=∇N (26)

where and are defined in Eqs. (2) and (15), respectively, and an overbar denotes a versor, i.e. and ( is the standard Euclidean norm). Clearly, all the quantities in the r.h.s. of Eq. (25) depend on the instantaneous value of the field, so that . This approach can be further generalized to equations having more than one additional first integral Franzosi (2012). See Refs. Rugh (2001); Davis et al. (2002) for similar approaches.

As we mention in Sec. VI, we find that the dynamics dictated by the lattice Hamiltonian in Eq. (2) brings the system to an asymptotic equilibrium state, characterized by well-defined values of observables such as the interaction or kinetic energy per particle, or the above-described instantaneous microcanonical temperature. Specifically, we observe that, after a transient whose duration depends on the initial state and the Hamiltonian parameters, the instantaneous value of said observables oscillates about an asymptotic value. This allows the definition of time averages as in Eq. (12). Fig. 6 shows some instances of the described equilibration process.

Typically, we initialize the dynamics on a suitably perturbed “plane-wave” state, , where and are random numbers uniformly chosen in , and control the magnitude of the random perturbations and is a normalization constant enforcing the desired particle density. It should be noted that unperturbed “plane wave” states can be either linearly stable or unstable depending on Smerzi et al. (2002); ext (). Plane-waves having an energy close to that of the ground-state are typically stable, and hence the relevant dynamics can be non ergodic. We checked that a suitable amount of noise destroys stability, so that equilibrium can be reached also for small energies. This may require long equilibration times. Conversely, for unstable plane-wave modes, a vanishingly small noise is sufficient to trigger a modulational instability that drives the system away from the initial state very quickly. We remark that these modulationally unstable states do not necessarily end up having a negative microcanonical temperature. Actually, for suitably large densities (or interaction strengths), the whole band of “plane-wave” modes can give rise to positive-temperature asymptotic states K.Ø. Rasmussen et al. (2000); ext ().

As we mention in Sec. VI the asymptotic average occupation distribution of the relevant modes of the dynamics provides a further proof that the system has reached equilibrium. This is demonstrated in Figs. 7 and 8. The histograms in the lower insets show that the instantaneous microcanonical temperature performs small oscillations around an asymptotic value, while the density plot in the upper insets show the quasimomentum average distribution. The scatter plots in the main figure show the average occupation of the lattice modes. In Fig. 7 the temperature is comparatively large, and neither extremal state is macroscopically occupied. Therefore, the prediction of Eq. (7) is fulfilled by single-particle states. As we discuss in Ref. VI, it is as if the interaction term simply acts as a heat bath maintaining the temperature non-interacting system. Note that the prediction of Eq. (7) applies for both positive and negative microcanonical temperatures.

Fig. 8 illustrates a case where Eq. (7) seems to fail. It refers to a two-dimensional lattice with strong defocusing nonlinearity, . The green scatter plot, obtained as in Fig. 7, shows that the occupation of the single particle modes is “larger than it should” at small energies, and bends towards the expected slope —i.e. the slope of the straight black line— only at high energies. This is because the system is in the condensed phase. Specifically, the average density is , while the condensate fraction is . As apparent from the orange scatter plot, the prediction of Eq. (7) is recovered when the Bogoliubov modes are considered. The same results are obtained for smaller interaction strengths, although the difference between the two scatter plots is not as dramatic as in the case presented in Fig. 8.

We observe that the scatter plot of versus deviates from a straight line also when the system condense into a localized state (not shown). The deviation is significant for energies close to that of the extremal localized state, while the scatter plot matches the expected linear behavior at the opposite end of the spectrum ext (). We expect that the linear behavior can be recovered on the whole energy spectrum if Bogoliubov modes are used instead of single-particle modes. However, owing to the localized character of the extremal state, the Bogoliubov approach is significantly more involved in this case. The localized character of the dynamical state might pose some problem with respect to the thermalization of the whole lattice and its subsystems. It is indeed clear that a large number of sublattices can be found which do not feature a localized density peak. In fact, most of the sublattices would have an average particle density much smaller than that of the whole system. We verified that the average distribution for the mode occupation agrees with Eq. (7) when calculated in a sublattice not containing the density peak ext ().

Finally, we checked that generalizing the grand-canonical Langevin approach introduced in Ref. Iubini et al. (2013), where is an external parameter, produces results in agreement with the prediction of Eq. (7) for both positive and negative temperatures ext (). This is one further evidence of the equivalence of the micronanonical and grand canonical ensemble for both signs of the temperature.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters