Non-Gibbs states on a Bose-Hubbard lattice

# Non-Gibbs states on a Bose-Hubbard lattice

Alexander Yu. Cherny Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, 141980, Dubna, Moscow region, Russia Center for Theoretical Physics of Complex Systems, Institute for Basic Science (IBS), Daejeon 34051, Republic of Korea    Thomas Engl Center for Theoretical Physics of Complex Systems, Institute for Basic Science (IBS), Daejeon 34051, Republic of Korea New Zealand Institute for Advanced Study, and Centre for Theoretical Chemistry and Physics, Massey University, Private Bag 102904 North Shore, Auckland 0745, New Zealand    Sergej Flach Center for Theoretical Physics of Complex Systems, Institute for Basic Science (IBS), Daejeon 34051, Republic of Korea New Zealand Institute for Advanced Study, and Centre for Theoretical Chemistry and Physics, Massey University, Private Bag 102904 North Shore, Auckland 0745, New Zealand
###### Abstract

We study the equilibrium properties of the repulsive quantum Bose-Hubbard model at high temperatures in arbitrary dimensions, with and without disorder. In its microcanonical setting the model conserves energy and particle number. The microcanonical dynamics is characterized by a pair of two densities: energy density and particle number density . The macrocanonical Gibbs distribution also depends on two parameters: the inverse nonnegative temperature and the chemical potential . We prove the existence of non-Gibbs states, that is, pairs which cannot be mapped onto . The separation line in the density control parameter space between Gibbs and non-Gibbs states corresponds to infinite temperature . The non-Gibbs phase cannot be cured into a Gibbs one within the standard Gibbs formalism using negative temperatures.

\RS@ifundefined

subsecref \newrefsubsecname = \RSsectxt \RS@ifundefinedthmref \newrefthmname = theorem  \RS@ifundefinedlemref \newreflemname = lemma

## I Introduction

Equipartition, ergodicity, and thermalization are essential properties, which an isolated many-body system has to have in order to qualify for the applicability of fundamental laws of statistical mechanics. Violation of the former, on the other side, could preserve coherence and might be of interest for, e.g., efficient information processing on classical and quantum levels. Surprising indications of transitions from ergodic to nonergodic dynamics have been reported in models of Josephson junction chains Pino et al. (2016) and Bose-Einstein condensates of ultracold atoms on optical lattices Mithun et al. (2018) upon “heating” the systems, i.e., upon increasing the average unbounded energy density.

For a macroscopic system with the energy being the only relevant conserved quantity — as in the case of the Josephson junction network Pino et al. (2016) — the canonical distribution function allows us to map any average energy density into a positive inverse temperature of the canonical distribution 111Note that this will hold independently of whether the microcanonical dynamics is ergodic (thermalizing) or nonergodic (nonthermalizing)..

In the presence of a second conserved quantity (e.g., the particle number) the situation changes. The additional constraint unfolds the existence of a non-Gibbs statistics of the one-dimensional Gross-Pitaevskii (GP) lattice also known as the discrete nonlinear Schrödinger equation with Hamiltonian

 H=∑l−J(ψ∗lψl+1+c.c.)+U2|ψl|4 (1)

and the Hamiltonian equations of motion , as obtained in Ref. Rasmussen et al. (2000). The GP dynamics preserves the energy and the total norm . The microcanonical dynamics at equilibrium — if existent — is defined by the two densities and , where is the number of lattice sites not (). In the macroscopic limit it follows that the Gibbs grand-canonical formalism becomes only applicable to the microcanonical dynamics for energy densities , despite the fact that the microcanonical dynamics can address states with . These states are therefore called non-Gibbs states Rasmussen et al. (2000). Subsequent studies addressed the question whether the microcanonical GP dynamics in the non-Gibbs phase is nonergodic Rumpf and Newell (2001); Rumpf (2004, 2007, 2008, 2009); Buonsante et al. (2017). Recent data Mithun et al. (2018) show that the dynamics stays ergodic; however, relaxation times quickly grow deep in the non-Gibbs phase, turning the system into a dynamical glass, which becomes quickly nonergodic for any practical purpose. The microscopic mechanism is related to the excitation of long-lived discrete breathers Flach and Willis (1998); Flach and Gorbach (2008). These discrete breathers become precise single-action excitations of the system in the infinite density limit, where the GP system turns into an integrable set of uncoupled anharmonic oscillators. The anomalously growing lifetime of these objects is the reason for the growth of relaxation times Mithun et al. (2018).

In this work we consider the corresponding quantum model, the Bose-Hubbard model on lattices with arbitrary dimension (note that the particular choice of lattice symmetry will not be of importance) Fisher et al. (1989). For simplicity we only show here its direct one-dimensional nearest-neighbour counterpart of the GP lattice (1):

 (2)

Here, and are the bosonic creation and annihilation operators on the th lattice site with the commutation relations , and is the site occupation number. The GP Hamiltonian (1) is the classical limit of the Bose-Hubbard (BH) Hamiltonian (2) for large occupation numbers , since quantum operators can be replaced by -numbers: , .

We will show that the non-Gibbs phase exists in the full quantum BH lattice model as well for energy densities in the - plane, where and are the average occupation number (filling factor) and energy per site, respectively, and is the on-site repulsion energy of the BH lattice Hamiltonian Fisher et al. (1989). We will also generalize to disorder potentials. Our general proof of the existence of nonthermal non-Gibbs states is relevant for a variety of experimental setups of Bose-Einstein condensates of ultracold atomic gases in optical lattice potentials, which are often used to study fundamental properties of matter Bloch et al. (2008).

This paper is organized as follows: in Sec. II, we obtain the line of infinite temperatures and, by means of the method of cumulant expansions, the thermodynamic susceptibilities for the Bose-Hubbard model. In the next section, we generalize these results to the Bose-Hubbard model with disorder (Sec. III.1) and with arbitrary local interactions (Sec. III.4). We consider the density of states of the system near the infinite temperature line and discuss the important role of boundedness of the single-particle spectrum. In the Conclusions, we review the obtained results and discuss possible applications and prospects.

## Ii Main results

### ii.1 The method

The limiting line , found in Ref. Rasmussen et al. (2000), corresponds to infinite temperatures of the canonical Gibbs ensemble. However, the non-Gibbs region above the line in the - parameter space is experimentally (numerically) accessible, since it is possible to initialize the lattice at any energy and total number of particles Rasmussen et al. (2000); Flach and Gorbach (2008). For instance, putting all bosons onto one site gives the energy per site , which is infinitely large in the thermodynamic limit. Dissipative boundaries were also proposed in order to drive the system into a non-Gibbs region Livi et al. (2006); Iubini et al. (2013). Here we address the question whether the limiting line also exists in the quantum Bose-Hubbard model and which equation it obeys.

In order to answer this question, we calculate the average energy in the canonical ensemble, that is, the internal energy as a function of the temperature, total number of particles, and the number of lattice sites, . The last equality is valid in the thermodynamic limit, because the internal energy and total number of particles are extensive properties of the system. The Gibbs canonical ensemble is, conceptually, the most convenient one to study this problem, because we know the full range of values of thermodynamic variables and . For a given density and for , the line goes through all possible values of in the - plane. The capacity

 (∂E/∂T)N,M=M(∂ε/∂T)n (3)

is positive, which implies that the energy reaches its maximum for , and, therefore, is the upper bound for the Gibbs region in the - plane, provided the limit exists.

The grand canonical potential is defined as with being the entropy, and is the chemical potential. Below we use two independent variables and instead of and to simplify calculations. For the grand canonical Gibbs ensemble, the partition function is given by

 Z(β,α,M)=Trexp(−β^H−α^N). (4)

Then the internal energy and total number of particles per site take the form

 ε(β,α) =−1M∂lnZ(β,α,M)∂β, (5) n(β,α) =−1M∂lnZ(β,α,M)∂α. (6)

The grand canonical ensemble is equivalent to the canonical one only for nonzero susceptibility , since otherwise the chemical potential is not a well-defined function of . The susceptibility vanishes at low temperatures for the Bose Hubbard model when the on-site repulsion is sufficiently large Fisher et al. (1989). As a consequence the system becomes a Mott insulator with zero compressibility. However, at high temperatures, the susceptibility is not zero, as shown below.

Excluding the variable from Eqs. (5) and (6) yields as a function of temperature and filling factor and hence the particular case at , which we are interested in.

### ii.2 Perturbation series at high temperatures

For the GP lattice, the value is shown to be finite and positive in the limit  Rasmussen et al. (2000). We assume that this property is also valid for the BH model and then show the consistency of the assumption. Thus, the fugacity (see, e.g., Ref. Huang (1987)) will tend to a constant value less than one in the same limit.

In order to obtain a perturbation series at high temperatures for the grand partition function, one can use the commutativity and and expand in a power series in the inverse temperature

 Z(β,α, M)=Tr[(1−β^H+β2^H2/2+⋯)exp(−α^N)] (7)

thus arriving at the expansion with the statistical moments . Here we denote and the brackets stand for the average over the density matrix in the zero-order approximation

 ^ρ0=exp(−α∑i^ni)Z0(α,M). (8)

The usage of the moment expansion is not convenient, since should be proportional to in the thermodynamic limit, while the moment is proportional to . In order to solve this problem, we reexpand in powers of thus arriving at the cumulant expansion with respect to Rushbrooke et al. (1974)

 lnZ(β,α,M)=λ0−βλ11!+β2λ22!−β3λ33!+⋯ (9)

where each cumulant depends on and , and . Cumulants can be obtained from the moments

 λ1=μ1, λ2=μ2−μ21, λ3=μ3−3μ1μ2+2μ31, (10)

and so on. Since and all its derivatives with respect to are proportional to in the thermodynamic limit, the terms of order and higher cancel each other, and each cumulant turns out to be proportional to .

When calculating the energy per site at infinite temperatures, it is sufficient to restrict ourselves to the first cumulant. For obtaining the capacity (3), we need the second order expansion

 lnZ(β,α,M)= (11)

### ii.3 The high-temperature expansion for the Bose-Hubbard model

We calculate the partition function in the zero-order approximation in the Fock basis for the operators and , where the density matrix is diagonal,

 Z0(α,M)=⎡⎣∞∑ni=0exp(−αni)⎤⎦M=1(1−e−α)M. (12)

Since the density matrix (8) is the exponential of a quadratic form of the bosonic operators and , the Wick-Bloch-De Dominicis theorem Bloch and Dominicis (1958) can be applied for calculating the average of a product of an arbitrary number of bosonic operators. According to the theorem, this average is given by the complete sum of all possible products of elementary binary averages (contractions)

 ⟨^a†i^aj⟩0 =δijn0. (13)

Here is the mean occupation number per site (lattice filling factor) at ,

 n0=⟨^a†i^ai⟩0=(1−e−α)∞∑k=0ke−αk=1eα−1, (14)

and the “nondiagonal” values are zero, since the average vanishes for .

The mean value of the energy (2) at is obtained as follows. The hopping terms in vanish, since they include only nondiagonal binary terms. We are left only with the interaction energy terms, for which in accordance with the Wick–Bloch–De Dominicis theorem, . Thus

 ⟨^H⟩0=MUn20. (15)

We substitute Eqs. (12) and (15) into Eq. (11) and use Eqs. (5) and (6) to derive

 ε =Un20+O(β), (16) n =n0[1−2βUn0(n0+1)]+O(β2). (17)

Here the symbol denotes terms of order and higher.

This system of equations allows us to obtain the upper border for the Gibbs region in the - plane. For infinite temperatures, we have and

 ε=Un2, (18)

which is the main result of this work. It coincides with the Gibbs–non-Gibbs border for the discrete nonlinear Schrödinger equation Rasmussen et al. (2000).

In the limit of infinite temperatures, the average energy per particle (18) is finite at fixed particle density. In the same limit, the density matrix (8) is independent of the energy, which makes all energies equiprobable at a fixed particle density. On the other hand, the upper bound of the spectrum of the Bose-Hubbard Hamiltonian grows as , which seemingly leads to a divergence of the average energy in the thermodynamic limit. The resolution of this apparent paradox is that the density of states decreases rapidly at high energies.

Equations (14) and (17) yield the parameter up to the first-order term proportional to :

 μ=−1βln(1+1/n)+2Un+O(β) (19)

 (∂μ∂n)β=1βn(n+1)+2U+O(β). (20)

The quantity (20) is directly related to the fluctuations of the total number of particles in the grand canonical ensemble, which remain finite in the limit of infinite temperatures,

 ⟨^N2⟩−⟨^N⟩2⟨^N⟩= 1βn(∂n∂μ)β=(n+1) ×[1−2βn(n+1)U+O(β2)]. (21)

In order to obtain the capacity (3), we have to consider the second-order terms in Eq. (11). In the same manner as for Eq. (15), we derive

 ⟨^H2⟩0−⟨^H⟩20=Mn0(n0+1)[2J2d+U2n0(5n0+1)], (22)

where is the dimension of the lattice. We rewrite the capacity per site in the more convenient form

 (∂ε∂T)n=−β2(∂ε∂β)n=−β2(∂n∂α)β∂(ε,n)∂(β,α), (23)

where we use the standard thermodynamic notation for the Jacobian determinant. By means of Eqs. (11)–(15) and (22), we arrive at the capacity per site in the lowest order in

 (∂ε∂T)n=β2n(n+1)[2J2d+U2n(n+1)]+O(β3), (24)

which is positive for any nonzero and turns zero at infinite temperature.

The positiveness of the expressions (21) and (24) implies that the system is stable at high temperatures.

The classical limit of the BH model, that is, the classical GP lattice can be evaluated in a similar manner for arbitrary lattice dimension by replacing sums by corresponding integrals. However, it is simpler to find the thermodynamic relations by considering the classical limit in the obtained quantum system equations. For instance, Eq. (14) becomes in this limit , and therefore, at infinite temperatures in the classical case (see the discussion in Sec. III.4.2 below). Similarly we can replace by in the equations for the susceptibility (20) and capacity (24), in order to find their classical limits. We note that the thermodynamic relations of the quantum and classical cases differ significantly for small particle densities.

Let us emphasize that the classical limit is not realized by a quantum system at infinite temperatures. Even on the infinite temperature line, the relation between the density of particles and fugacity () varies: in the quantum model , while in the classical one, . This is because the energy remains finite in the limit of infinite temperatures at a given particle density. Note that the classical case can be obtained from the quantum one in the limit , but an inverse scheme does not exist.

## Iii Generalizations

### iii.1 Adding disorder

Let us consider the Bose-Hubbard model with disorder by adding to the Hamiltonian (2) a disorder potential with random on-site energies , obeying some probability density distribution (see, e.g., Ref. Flach (2016)). Their average value is assumed to be zero, while the variance is finite.

The disorder potential does not change the upper bound for the Gibbs region (18). Its contribution to the energy per site is equal to zero in the limit of infinite temperatures: . However, the presence of disorder influences the second-order term in the high temperature expansion (11) into

 ⟨^H2⟩0−⟨^H⟩20= Mn0(n0+1) ×[σϵ+2J2d+U2n0(5n0+1)], (25)

and, hence, the capacity per site

 (∂ε∂T)n=β2f(n)+O(β3), (26)

where we use the notation

 f(n)=n(n+1)[σϵ+2J2d+U2n(n+1)]. (27)

The border of the Gibbs region and the capacity at high temperatures are independent of the correlations in disorder for . Therefore uncorrelated disorder with finite variance will result in the same border between the Gibbs and the non-Gibbs regime given by (LABEL:enborder). This result holds independent of the dimensionality of the lattice, and generalizes the previously obtained consideration of a one-dimensional disordered GP lattice Flach (2016).

### iii.2 The density of states in the vicinity of the Gibbs - non-Gibbs separation line

Let us calculate the density of states in the vicinity of the border separating the Gibbs and non-Gibbs states. Instead of following the route by calculating the entropy near by means of the grand canonical potential: , we use Eq. (26) for the capacity, since entropy and capacity per site are related by the equation . Then the expression (26) leads to and after integration over to

 s=s0(n)−β22f(n)+O(β3) (28)

with being the entropy per site at infinite temperature

 s0(n)=(n+1)ln(n+1)−nlnn. (29)

The last equation can be obtained directly from the density matrix (8) using the well-known relation for the entropy, .

In the same manner, we derive from Eqs. (18) and (26)

 ε=Un2−βf(n)+O(β2). (30)

Excluding the inverse temperature from Eqs. (28) and (30) yields

 s=(n+1)ln(n+1)−nlnn−(ε−Un2)22f(n), (31)

where is given by Eq. (27). By using the relation , we arrive at the density of states in the microcanonical ensemble near the Gibbs - non-Gibbs border

 W(E,N,M)=(M+N)M+NNNMMexp⎡⎣−M(ε−Un2)22f(n)⎤⎦, (32)

where and . This is the main contribution to the asymptotics of the density of states in the thermodynamic limit , . The density of states reaches its maximum precisely at the border .

The density of states for the discrete GP Hamiltonian (1)

 Wcl(E,N,M)=(NeM)Mexp⎡⎣−M(ε−Un2)22fcl(a)⎤⎦, (33)

can be inferred from Eq. (32) in the classical limit , where . In the particular case (the ideal Bose gas), the parabolic dependence of the entropy on the energy in the microcanonical ensemble was obtained up to a coefficient in Ref. Buonsante et al. (2017) by another method.

Note that for both quantum and classical cases, the entropy is monotonously increasing with at fixed , while it is nonmonotonous for increasing at fixed . This result is analogous to the one obtained in Ref. Rumpf (2004) for the one-dimensional classical GP lattice.

### iii.3 Bounded single particle spectra and non-Gibbs states

The Gibbs - non-Gibbs separation line (18) exists due to the single-particle dispersion being bounded, e.g., as for the one-band Bose-Hubbard model, (LABEL:BH). Note that generalizations to a finite number of bands generated by more complex lattice structures are straightforward. At variance, an infinite number of bands or the single-particle dispersion, or simply the case of free particle , leads to unbounded spectra. As a consequence non-Gibbs states are renormalized to infinite values in the - plane and vanish from any consideration. Indeed, in the particular case , we have in general an ideal Bose gas with the dispersion with being the band number. Then the total number of particles is given by the standard expression

 N=∑q,λ1/[exp(βJλ(q)+α)−1], (34)

where the quasimomentum runs over uniformly distributed points in the Brillouin zone. In the case of a one-band structure, only one band, say, contributes to the sum in (LABEL:multband). Then the limit yields which allows for a solution of as a function of . However, for an infinite number of bands, diverges. Therefore the assumption about the finiteness of when is not applicable anymore.

### iii.4 Further extensions

Non-Gibbs states can be obtained as well for the following generalized Bose-Hubbard models

 ^H=∑i≠jJ(i−j)^a†i^aj+∑iu(^ni). (35)

The sum in the hopping term includes all sites, and the single particle dispersion is assumed to be bounded. The standard Bose-Hubbard Hamiltonian (2) is a particular case of Eq. (35) with and .

If the condition

 limx→+∞u(x)x→+∞ (36)

is satisfied, then a non-Gibbs phase with anomalous scaling properties of the energy will emerge. This phase cannot be described by a grand partition function with negative temperature due to the divergence of the former. If on the other side Eq.(36) is not satisfied, i.e., with , then the line becomes the border between Gibbs states with positive and negative temperatures (see also Refs. Samuelsen et al. (2013); Buonsante et al. (2017)).

#### iii.4.1 The Gibbs–non-Gibbs transition line for the generalized Bose-Hubbard model

By complete analogy with Sec. II.3, one can show that at infinite temperatures the hopping terms do not contribute to the average energy. Its value is determined by the interaction terms and calculated with the density matrix (8)

 ε=(1−e−α)+∞∑m=0u(m)e−αm. (37)

The density of particles is still given by Eq. (14). Substituting the exponential into Eq. (37), we arrive at the curve of the energy density in the - plane

 (38)

above which non-Gibbs states for the Hamiltonian (35) emerge.

As an example, we calculate the transition line for the local interaction energy

 u(x)=Uδ[Γ(x+1)Γ(x+1−δ)−1Γ(1−δ)] (39)

with being Euler’s gamma function. This form of the potential generalizes the power-law interaction on each site with a positive integer exponent , for which . The last term in the right-hand side brackets in Eq. (39) is needed to enforce the condition , which the potential should satisfy for any positive . The long-range asymptotics of that potential is given by . The line of infinite temperatures is calculated by means of Eq. (38), which yields

 ε(n)=UδΓ(1−δ)[2F1(1,1;1−δ;nn+1)n+1−1], (40)

where is the hypergeometric function Abramowitz and Stegun (1974). As discussed above, this line separates the Gibbs and non-Gibbs phases when . For positive integer exponent , Eq. (40) is simplified

 ε(n)=UΓ(δ)nδ. (41)

This equation also gives the main term of the long-range asymptotics of Eq. (40) for arbitrary positive . For the standard Bose-Hubbard model with two-body interactions, we have and arrive at Eq. (18) (see blue line in Fig. 1). For the impact of the potential is reduced to a renormalization of the chemical potential of a noninteracting ideal Bose gas (see red line in Fig. 1). As a result, this case separates the appearance of non-Gibbs phases for from cases with , where the infinite temperature line separates Gibbs states with positive and negative temperatures (see black line in Fig. 1).

It is straightforward to extend the hopping network between sites. For instance, for binary interactions a one-band Hamiltonian takes the form

 ^H=∑i≠jJ(i−j)^a†i^aj+12∑i,jU(i−j)^a†2i^a2j. (42)

Using the same method as in Sec. II, we arrive at non-Gibbs states above the line with , which is assumed to be finite and positive.

Similar to Sec. III.1, we can show that the disorder does not modify the border line between Gibbs and non-Gibbs states for the both models, Eqs. (35) and (42).

#### iii.4.2 The Gibbs–non-Gibbs transition line for the discrete nonlinear Schrödinger equation

We also consider the generalizations of the GP lattice (1), whose Hamiltonian is given by Eq. (35) with the replacement , . One can calculate the Gibbs–non-Gibbs transition line in the same way as in Sec. II. However, the easiest way is to obtain it in the limit of large occupation numbers, when the quantum model approaches the classical one. Then the density of particles in Eq. (14) is and the sum in Eq. (37) tends to the integral . We finally obtain

 ε(n)=1n∫∞0u(x)e−x/ndx. (43)

We observe that the line of infinite temperatures is given by the Laplace transform of the interaction potential through . Similar to the quantum case, non-Gibbs states for the Hamiltonian (35) emerge for energy densities if the condition (36) is satisfied.

The classical analog of the potential (39) in the Bose-Hubbard model is

 u(x)=Uxδ/δ. (44)

Using Eq. (43), we find the Gibbs–non-Gibbs separation line given by Eq. (41). Thus, for an integer exponent , the infinite temperature line for the quantum and classical models is identical, but for noninteger it is not, see Fig. 1b.

One can see that for the interaction with being an entire function of , the line matches Eq. (43). Here the symbols stand for the normal ordering of the Bose operators and . In order to prove this equation, it is sufficient to expand an arbitrary function in a Taylor series and show the equation validity for each term of the series. Thus, we consider the particular case for any nonnegative integer , for which (here we omit the lattice site index for simplicity). By applying the methods of Sec. II and using the Wick-Bloch-De Dominicis theorem, we arrive at . The last equality is valid in zeroth order in . This coincides with the relation given by Eq. (43).

Samuelsen et al. Samuelsen et al. (2013) considered a one-dimensional chain with the saturable nonlinear potential (with negative ), and obtained the line of infinite temperatures using a transfer integral approach:

 ε=νexp(1/n)E1(1/n), (45)

where is the exponential integral Abramowitz and Stegun (1974). We note that our method is generating and therefore confirming that result using the above relation (43). We also note that our approach allows one to generalize the result (45) to arbitrary lattice dimensions.

### iii.5 The absence of the Gibbs–non-Gibbs transition line in a Josephson junction array model

Let us discuss why the Gibbs–non-Gibbs transition line is absent in a Josephson junction array

 ^H=∑iEC2^q2i−EJcos(^φi−^φi+1). (46)

This is a simple model describing the Josephson junction network of weakly coupled superconducting islands with the Josephson and charging energies (see, e.g., Ref. Pino et al. (2016)). The operators of phase and charge of the superconducting islands obey the commutation relations . They can be considered as the -component of momentum operators and on the ring of unit radius.

The Hamiltonian (46) is derived from the Bose-Hubbard model (2) in the limit of infinite densities Bruder et al. (2005). In the polar representation, the bosonic operators follow as

 ^ai=e−i^φi√n+δ^ni,^a†i=√n+δ^niei^φi (47)

with being the operator of deviation from the average occupation number per site. Since the density is assumed to be large, the operator takes now all integer values. The operators and obeys the same commutation relations as and .

Substituting the operators (47) into Eq. (2) yields

 ^H= ∑iU2δ^n2i−2Jncos(^φi−^φi+1)+U2(2n−1)δ^ni +MU2(n2−n). (48)

Here we neglect the operator in the hopping terms, because . The operator still commutes with the Hamiltonian (48). The substitution , and yields the Josephson junction array Hamiltonian (46) up to additive renormalizations of the energy and chemical potential, which do not influence the physical properties of the model. Note that we can relax the constraint and consider an arbitrary value of the charge , which is a conserved quantity.

Thus the model has two integrals of motion but lacks a line of infinite temperatures in the energy-charge plane. If we suppose that the infinite temperature line exists and try to construct the density matrix by analogy with Sec. II.2, we will obtain a divergence in the partition function , because takes all integer values in the Fock basis. The derivation of the model (48) from the Bose-Hubbard model reveals the reason why: the operators of local charge are supposed to be unbounded from below, while they should actually obey the inequality .

## Iv Conclusions

In this work we studied the equilibrium properties of the repulsive quantum Bose-Hubbard model in arbitrary lattice dimensions, with and without disorder. For the classical limit of a Gross-Pitaevskii lattice in one dimension, the existence of a non-Gibbs phase was proven for two-body interactions Rasmussen et al. (2000), and for a case with saturable nonlinearity Samuelsen et al. (2013). We extend these results to the full many-body quantum domain including its classical limit, and to arbitrary lattice dimensions, and to a generalized set of lattice couplings and interaction functions. We proved the existence of non-Gibbs states in the particle number and energy density control parameter space. The separation line in the density control parameter space between Gibbs and non-Gibbs states corresponds to the infinite temperature line where . We substantially extend these results both for the quantum and classical cases [see Eqs. (38) and (43) respectively] for a much wider class of interactions obeying asymptotic functional dependence on the particle density (36). This dependence tells us that the upper bound of the Hamiltonian spectrum at given particle density grows faster that the system size in the thermodynamic limit. For this reason, the non-Gibbs phase cannot be cured into a Gibbs one within the standard Gibbs formalism using negative temperatures (see Sec. III.4).

The existence of a non-Gibbs phase needs an infinite temperature line at finite densities in the density control parameter space in the first place. A prerequisite of such a line is the existence of a density control parameter space, which is at least two-dimensional. Then there needs to be at least one more conserved quantity in addition to the energy (since we consider Hamiltonian systems). The particle number in the Bose-Hubbard system is precisely the second integral of motion we need. Once the infinite temperature line is obtained, a second condition to be satisfied in order for a non-Gibbs phase to exist is the nonsaturability of the interaction potential (36). This property will apparently lead to a modification of finite-size scalings of energy and entropy densities. How precisely, remains to be addressed in future work.

From a practical perspective, the Bose-Hubbard lattice is a projection of an infinite band system with unbounded spectrum onto one (or a few) band(s) in the tight-binding approximation. The obtained results can be applied to a realistic system with a periodic potential with some caveats. In the case of one band, the limit of high temperatures actually means that but , where is the gap between the lowest bands and is the Boltzmann constant. On the other hand, the chemical potential should be restricted as well: . Then we arrive at the reliable temperature range

 α/Δ≪β≪α/J, α/U, (49)

where the parameter is directly related to the lattice filling factor: .

A natural question is about the impact of non-Gibbs states on the dynamics of the system. For the classical GP lattice case, the non-Gibbs states correlate with an entering of the system into a dynamical glass phase, in which relaxations of local fluctuations of the particle density are slowing down in a dramatic way Mithun et al. (2018). This happens despite the fact that the system shows nicely ergodic (chaotic) spots. It turns out that these spots are however strongly diluted, and the final time to reach ergodicity is much larger than a typical inverse Lyapunov exponent Kantz et al. (1994). In Ref. Pino et al. (2016) a similar scenario for a Josephson junction network is expected to lead to the celebrated many-body localization in the quantum case, which is a kind of Anderson localization in the Fock space. We therefore expect that the Bose-Hubbard model could also show a transition into a many-body localized regime, particularly, in the non-Gibbs phase.

## V Acknowledgement

This work was supported by the Institute for Basic Science, Project Code IBS-R024-D1. T.E. acknowledges financial support by the Alexander-von-Humboldt foundation through the Feodor-Lynen Research Fellowship Program No. NZL-1007394-FLF-P. A.Yu.Ch. thanks the hospitality of the IBS Center for Theoretical Physics of Complex Systems, where most of the work was conducted, and acknowledges support from the JINR–IFIN-HH projects.

## References

• Pino et al. (2016) Manuel Pino, Lev B. Ioffe,  and Boris L. Altshuler, “Nonergodic metallic and insulating phases of Josephson junction chains,” PNAS 113, 536–541 (2016).
• Mithun et al. (2018) Thudiyangal Mithun, Yagmur Kati, Carlo Danieli,  and Sergej Flach, “Weakly nonergodic dynamics in the gross-pitaevskii lattice,” Phys. Rev. Lett. 120, 184101 (2018).
• (3) Note that this will hold independently of whether the microcanonical dynamics is ergodic (thermalizing) or nonergodic (nonthermalizing).
• Rasmussen et al. (2000) K. Ø. Rasmussen, T. Cretegny, P. G. Kevrekidis,  and Niels Grønbech-Jensen, “Statistical mechanics of a discrete nonlinear system,” Phys. Rev. Lett. 84, 3740–3743 (2000).
• (5) In the literature about the discrete nonlinear Schrödinger equation, the notations , , , and are usually used.
• Rumpf and Newell (2001) Benno Rumpf and Alan C. Newell, “Coherent structures and entropy in constrained, modulationally unstable, nonintegrable systems,” Phys. Rev. Lett. 87, 054102 (2001).
• Rumpf (2004) Benno Rumpf, “Simple statistical explanation for the localization of energy in nonlinear lattices with two conserved quantities,” Phys. Rev. E 69, 016618 (2004).
• Rumpf (2007) B Rumpf, “Growth and erosion of a discrete breather interacting with Rayleigh-Jeans distributed phonons,” Europhys. Lett. 78, 26001 (2007).
• Rumpf (2008) Benno Rumpf, “Transition behavior of the discrete nonlinear Schrödinger equation,” Phys. Rev. E 77, 036606 (2008).
• Rumpf (2009) Benno Rumpf, “Stable and metastable states and the formation and destruction of breathers in the discrete nonlinear Schrödinger equation,” Physica D 238, 2067 – 2077 (2009).
• Buonsante et al. (2017) P. Buonsante, R. Franzosi,  and A. Smerzi, “Phase transitions at high energy vindicate negative microcanonical temperature,” Phys. Rev. E 95, 052135 (2017).
• Flach and Willis (1998) S. Flach and C.R. Willis, “Discrete breathers,” Phys. Rep. 295, 181 – 264 (1998).
• Flach and Gorbach (2008) Sergej Flach and Andrey V. Gorbach, “Discrete breathers — advances in theory and applications,” Phys. Rep. 467, 1 – 116 (2008).
• Fisher et al. (1989) Matthew P. A. Fisher, Peter B. Weichman, G. Grinstein,  and Daniel S. Fisher, “Boson localization and the superfluid-insulator transition,” Phys. Rev. B 40, 546–570 (1989).
• Bloch et al. (2008) Immanuel Bloch, Jean Dalibard,  and Wilhelm Zwerger, “Many-body physics with ultracold gases,” Rev. Mod. Phys. 80, 885–964 (2008).
• Livi et al. (2006) Roberto Livi, Roberto Franzosi,  and Gian-Luca Oppo, “Self-localization of bose-einstein condensates in optical lattices via boundary dissipation,” Phys. Rev. Lett. 97, 060401 (2006).
• Iubini et al. (2013) S Iubini, R Franzosi, R Livi, G-L Oppo,  and A Politi, “Discrete breathers and negative-temperature states,” New Journal of Physics 15, 023032 (2013).
• Huang (1987) K. Huang, Statistical mechanics (Wiley, 1987).
• Rushbrooke et al. (1974) G. S. Rushbrooke, G. A. Baker, Jr.,  and P. J. Wood, Phase Transitions and Critical Phenomena, edited by C.M. Domb and M.S. Green, V. 3 (Academic Press, New York, 1974) Chap. 5.
• Bloch and Dominicis (1958) Claude Bloch and Cyrano De Dominicis, “Un développement du potentiel de gibbs d’un système quantique composé d’un grand nombre de particules,” Nuclear Physics 7, 459 – 479 (1958).
• Flach (2016) Sergej Flach, “Spreading, nonergodicity, and selftrapping: A puzzle of interacting disordered lattice waves,” in Nonlinear Dynamics: Materials, Theory and Experiments: Selected Lectures, 3rd Dynamics Days South America, Valparaiso 3-7 November 2014, edited by Mustapha Tlidi and Marcel. G. Clerc (Springer International Publishing, Cham, 2016) pp. 45–57.
• Samuelsen et al. (2013) Mogens R. Samuelsen, Avinash Khare, Avadh Saxena,  and Kim Ø. Rasmussen, “Statistical mechanics of a discrete Schrödinger equation with saturable nonlinearity,” Phys. Rev. E 87, 044901 (2013).
• Abramowitz and Stegun (1974) Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables, (Dover Publications, Incorporated, 1974).
• Bruder et al. (2005) C. Bruder, R. Fazio,  and G. Schön, “The Bose-Hubbard model: from Josephson junction arrays to optical lattices,” Ann. Phys. (Leipzig) 14, 566–577 (2005).
• Kantz et al. (1994) H. Kantz, R. Livi,  and S. Ruffo, “Equipartition thresholds in chains of anharmonic oscillators,” Journal of Statistical Physics 76, 627–643 (1994).
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

371339

How to quickly get a good answer:
• Keep your question short and to the point
• Check for grammar or spelling errors.
• Phrase it like a question
Test
Test description