text:
# I Introduction

## I Introduction

The physics of non-equilibrium systems includes a broad class of phenomena, such as the physics of steady states, relaxation and dynamics far from equilibrium. Dynamical processes which range from those in the early universe, to ultra-relativistic heavy ion collisions and the formation of the quark-gluon plasma all involve non-equilibrium physics in an essential manner, perhaps requiring an understanding of the physics beyond linear response. Physics of steady states contains interesting non-equilibrium phenomena, such as transport, spatially varying observables, and is also important to the study of systems where the time scales for local equilibration are smaller than macroscopic time scales which might describe some global evolution of the system. Non-equilibrium steady states can be realized in many ways, such as the placement of a systems in a thermal gradient, or in an environment which provides shearing, pressure gradients and so forth. In this article, we will examine the statistical mechanics of the non-equilibrium steady states of a classical field theory in thermal gradients. This will allow us to understand the behavior of the theory under these non-equilibrium conditions and to consider problems related to the range or validity of local equilibrium, linear response, equilibrium thermodynamics and statistical mechanics. In particular, the thermal gradients we study can be quite strong and in such situations, it is natural to ask whether the system always relaxes to local equilibrium.

When systems are near equilibrium, one expects linear response to provide a description of the transport coefficients. However, there is no means to address its regime of validity within the theory itself. Furthermore, even in this regime, in low dimensions, $`d2`$, it has been argued using kinetic theory that linear response often does not hold; the Green-Kubo autocorrelation functions are expected to behave as $`t^{d/2}`$, leading to divergent transport coefficients. Such divergences of the Green-Kubo integral have been observed in certain low dimensional systems such as the FPU model  or the diatomic Toda lattice , and seem to be endemic to Hamiltonians which conserve total momentum. Other studies have also found thermal transport in the linear regime to diverge  or have focused on somewhat exotic models , while strong thermal gradients have been studied in cellular automata.

In this work, we attempt to address many of the basic questions regarding the non-equilibrium properties of the $`\varphi ^4`$ theory on the lattice in (1+1) dimensions, by studying the model both near and far from equilibrium. We choose the $`\varphi ^4`$ theory since it is a prototypical model which appears in a variety of contexts, including particle physics. From the outset, we should point out that our work has two limitations, namely that it is classical and that it is a lattice field theory. On the other hand, we make no further approximations and we analyze the model from first principles without any dynamical assumptions. This will allow us to answer interesting physical questions that cannot yet be addressed in the full quantum case. Furthermore, the approach we adopt can be generalized to other classical lattice field theories in a straightforward manner. Our main objective is to develop a comprehensive understanding of the underlying dynamics of the scalar field theory in thermal gradients and to lay the ground work for further analysis. As such, we shall provide the necessary details of our methods for further work and in such a manner so that can be easily generalized to other models.

Classical field theory is relevant to high temperature behavior of quantum field theories: For instance, recently, it was used to derive properties of the standard electroweak model at finite temperature . Properties of the quantum $`\varphi ^4`$ theory in equilibrium has also been studied previously . However, the relation between the physical quantities in the classical lattice field theory and the quantum theory is far from trivial and it is beyond the scope of this work. The dynamics of classical field theories is of interest on its own right, which in our case can be thought of as the dynamics of an anharmonic chain. While we discuss the results of the (1+1) dimensional theory in this work, it is not an essential limitation of our model or approach; we do find that the basic physics understanding developed in (1+1) dimensions carries over to results in (3+1) dimensions, although the latter will be discussed elsewhere. As such, our analysis in (1+1) dimensions is not specific to 1-d systems. It should be noted that classical $`\varphi ^4`$ theory has been considered in various contexts in the past: The Lyapunov spectra has been computed in the microcanonical ensemble in massless and massive models. Equilibration of the model was studied in .

In addition to elucidating the physics underlying some of the non–equilibrium phenomena of field theories, we believe that our results shed light on the nature of non-equilibrium statistical ensembles. The extension of the Gibbs ensemble to the non-equilibrium steady state remains unanswered. However, limited results to date suggest that the theory is far from trivial, including divergent Gibbs entropy and singular steady-state measures. While some approaches exist, such as maximum entropy or projection techniques to construct non-equilibrium operators, assumptions must be made about the non-equilibrium state in order to compute its properties. By using classical field theory as our starting point, we can use existing techniques to construct non-equilibrium steady states without assumptions on the dynamics of the model which are symptomatic to other approaches. A seemingly simple question concerns the thermal profile $`T(x)`$ which develops in a system when it is in a thermal gradient. In other approaches, often some form is assumed for the profile $`T(x)`$, which in principle should be dynamically obtained, as we shall do here. Furthermore, in our study, whether local equilibrium is achieved is not an assumption but is determined dynamically by the system. We will see that there are qualitative differences between the behavior of the system in the various regimes. These are characterized in Table 1.

We would like to emphasize once again the motivation for this analysis: While the lattice model we work with does have well defined thermal transport, a single component scalar continuum field theory does not support thermal conductivity in the usual sense. It would be interesting to generalize to more complex theories with more than one conserved current, such as the two component scalar field theory, gauge theory with matter and so on. On the other hand, rather than introduce additional degrees of freedom and observables such as charge or matter density, we prefer first to focus on some important questions in the lattice theory. By focusing on the simplest lattice theory, we are able to clearly present an approach to the non-equilibrium statistical mechanics of classical lattice models from first principles, which can then be applied to other models like the ones mentioned above. Our results are also of interest to the statistical mechanics of many body systems in non-equilibrium. So as is, we do not attempt to make claims concerning the continuum limit but rather concentrate on elucidating the physics of the lattice theory near and far from equilibrium.

In section 2, we describe the model we study and how we analyze the theory, particularly paying attention to the way the temperature boundary conditions are implemented. In section 3, we discuss thermal transport in our theory. In section 4, we analyze the equilibrium physics of the model and in section 5, the non–equilibrium physics. In particular, we study the thermal conductivity and its temperature dependence. We analyze the thermal profiles and establish that the linear response theory works even for visibly curved profiles, up to certain thermal gradients. We further examine the relations between the various physical quantities such as entropy, speed of sound, heat capacity and the thermal conductivity. We end with a discussion in section 6.

Table 1: Behavior of the $`\varphi ^4`$ theory under varying thermal gradients.

Here, $`\mathrm{}`$ is the mean free path as explained in section V.

| Regime |  | Properties |
| --- | --- | --- |
| Global Equilibrium (GE) |  | $`T=0`$, $`f(\pi ,\varphi )\mathrm{exp}[H/T]`$. |
| Local Equilibrium I (LE-I) |  | $`\mathrm{\Delta }T/T1`$; $`T=`$constant; |
|  |  | Fourier’s law holds globally; |
|  |  | Agreement with linear response theory. |
| Local Equilibrium II (LE-II) |  | $`\mathrm{}T/T1/10`$; $`T`$constant; |
|  |  | Fourier’s law holds locally; |
|  |  | Small deviations from linear response theory; |
|  |  | Existence of boundary temperature jumps. |
| Local Non-Equilibrium (LNE) |  | $`\mathrm{}T/T1/10`$; $`T`$constant; |
|  |  | Local equilibrium description inadequate; |
|  |  | Definition of temperature ambiguous; |
|  |  | No suitable definition for transport coefficients. |

## II The Model

We start with the $`\varphi ^4`$ Lagrangian (with the metric convention $`(,+)`$),

$$=\frac{1}{2}\left(\frac{\stackrel{~}{\varphi }(\stackrel{~}{x})}{\stackrel{~}{x}_\mu }\right)^2+\frac{1}{2}\stackrel{~}{m}^2\stackrel{~}{\varphi }(\stackrel{~}{x})^2+\frac{\stackrel{~}{g}^2}{4}\stackrel{~}{\varphi }(\stackrel{~}{x})^4.$$
(1)
We discretize and perform the rescaling

$$\varphi _x(t)=a\stackrel{~}{g}\stackrel{~}{\varphi }(\stackrel{~}{x},\stackrel{~}{t}),t=\stackrel{~}{t}/a,x=\stackrel{~}{x}/a,m^2=\stackrel{~}{m}^2a^2,$$
(2)
where $`a`$ is the lattice spacing. We then obtain the corresponding Hamiltonian where the lattice spacing is scaled out

$$H(\pi ,\varphi )=\frac{1}{2}\underset{i}{}\left[\pi _i^2+\left(\varphi _i\right)^2+m^2\varphi _i^2+\frac{1}{2}\varphi _i^4\right].$$
(3)
Here $`k=1,2,\mathrm{}L`$ runs over all sites in the lattice, the lattice derivative is $`\varphi _k\varphi _{k+1}\varphi _k`$. The resulting equations of motion are:

$$\dot{\varphi }_i=\pi _i,\dot{\pi }_i=(^2\varphi )_im^2\varphi _i\varphi _i^3.$$
(4)
Here, we defined the lattice Laplacian as $`(^2\varphi )_k\varphi _{k+1}2\varphi _k+\varphi _{k1}`$.

### A  Finite Temperature Equilibrium Dynamics: $`T=0`$

Starting from the microcanonical dynamics (4), we can develop a realization of the constant temperature dynamics using the global demons of Ref. . In this approach, auxiliary variables are added to the systems which dynamically emulate the presence of a heat bath. This type of dynamics, while not as rigorously understood, has been shown to converge much faster than the optimized hybrid Monte Carlo methods. Further, it does not suffer as much from critical slowing down.

When we are interested in studying the statistical properties of a system described by an action $`S(\phi )`$, where $`\phi =(\phi _1,\phi _2,\mathrm{}\phi _n)`$ are the degrees of freedom, we usually start with the definition of a statistical measure, such as

$$f𝒟\mu (\phi )\mathrm{exp}[S(\phi )/T]𝒟\mu (\phi ).$$
(5)
Here, $`𝒟\mu (\phi )`$ might include constraints in the dynamical space, as in the case of motion on curved manifolds, such as Lie groups. Because we know the measure, steady-state values of observables are readily determined. For an arbitrary observable $`𝒪`$, we have

$$𝒪=\frac{1}{Z}𝒟\mu (\phi )e^{S(\phi )/T}𝒪,Z=𝒟\mu (\phi )e^{S(\phi )/T}.$$
(6)
While the approach we discuss is suited to general measures, in this article we use the measure $`𝒟\mu (\phi )=𝒟\phi `$ over the phase space, where $`\phi `$ will typically represent canonically conjugate coordinates and momenta, $`\phi =(\varphi ,\pi )`$, and $`S(\phi )`$ is taken as a Hamiltonian, $`S(\varphi ,\pi )=H(\varphi ,\pi )`$. While the dynamics of the model may now be easily implemented using equations of motion, (4) — often referred to as the molecular dynamics method — we would like to add finite temperature constraints to the equations. For this to happen, we must no longer evolve on the constant energy surface, so that $`S(\phi )`$ should no longer be conserved. The method we discuss here is reminiscent of Parisi and Wu’s stochastic quantization, although the one adopted in our work is deterministic and time-reversal invariant. It is also a versatile approach in that it has been applied to systems with non-trivial measures, such as Lie algebraic Hamiltonians, equilibrium and non-equilibrium quantum systems, atomic clusters and molecules, magnetic materials and lattice models.

There are many formulations of this dynamics, initially motivated by the approaches of Nosé and Hoover. Consider the following equations of motion for a thermostatted site labeled by $`k`$:

$$\dot{\varphi }_k=\pi _k,\dot{\pi }_k=\frac{S}{\varphi _k}\frac{dG(w_k)}{dw_k}F(\pi _k)\frac{dG^{}(w_k^{})}{dw_k^{}}F^{}(\pi _k).$$
(7)
We have added two additional degrees of freedom, $`w_k,w_k^{}`$, which couple through forces indicated above. These extra degrees of freedom, so called ‘demons’, may be coupled either to the fields, $`\varphi _k`$, or to their ‘momenta’, $`\pi _k`$. Here, we choose to couple them only to $`\pi _k`$’s in order to have the ability to apply thermostats locally at any one site. The microcanonical limit is recovered when these extra degrees of freedom are decoupled. In this extended space, $`\phi =(\varphi _i,\pi _i,w_k,w_k^{})`$ , we define a new action which is the old one plus additional terms for the demons:

$$f(\varphi ,\pi ,w,w^{})=\mathrm{exp}\left(\left[S(\varphi ,\pi )+\underset{\genfrac{}{}{0pt}{}{k:\mathrm{thermostatted}}{\mathrm{sites}}}{}\left(G(w_k)+G^{}(w_k^{})\right)\right]/T\right).$$
(8)
In contrast to microcanonical dynamics, where the Hamiltonian is a constant of the motion, $`f`$ is not preserved by the constant temperature dynamics. While the choice of forces in the equations of motion as well as that of $`f`$ are seemingly arbitrary, steady state expectation values will be independent of these under reasonably general conditions. Consequently, they are chosen to optimize convergence of the physical variables.

To find the dynamics associated with the demons $`w_i,w_i^{}`$, we simply require that $`f`$ satisfy a continuity (Liouville) equation in the configuration space $`\phi =(\varphi _i,\pi _i,w_k,w_k^{})`$:

$$0=\frac{f}{t}+\underset{i}{}\frac{(\dot{\phi }_if)}{\phi _i}.$$
(9)
This is equivalent to requiring that the master equation, enforcing conservation of probability under evolution of the ensemble, be satisfied. By substituting the equations of motion into the continuity equation, and using the definition of $`f`$, we can derive solutions for $`\dot{w}_k,\dot{w}_k^{}`$:

$$\dot{w_k}=\pi _kF(\pi _k)T\frac{dF(\pi _k)}{d\pi _k},\dot{w^{}}_k=\pi _kF^{}(\pi _k)T\frac{dF^{}(\pi _k)}{d\pi _k}.$$
(10)
By construction, this dynamics preserves the measure Eq. (8), so that time averages of observables on a given trajectory will converge to the configuration space average over the canonical measure. There is clearly some freedom in defining the dynamics; namely, the functions $`G(w),G^{}(w^{})`$ and $`F(\pi ),F^{}(\pi )`$. The only restriction on $`G(w),G^{}(w^{})`$ is that the measure Eq. (8) leads to a finite integral; in general the auxiliary variables $`w`$ can have any desired measure. In practice, highly non-linear functions are impractical since they will require small integration time steps. For these reasons, it is convenient to take $`G(w),G^{}(w)`$ to be $`\mu w^2/2`$ or $`\mu ^{}w^4/4`$, where $`\mu ,\mu ^{}`$ are positive constants. The constant couplings $`\mu ,\mu ^{}`$ of the demons to the physical degree of freedom are in principle arbitrary as long as the phase space integration is finite. Choosing these couplings to be too weak will make the convergence slow or choosing them to be too strong will lead to small time steps in the evolution, so that the couplings are chosen to optimize the convergence of physical observables. A necessary condition for $`F(\pi ),F^{}(\pi )`$, on the other hand, is for it to be at least linear in its argument, the minimal requirement for the existence of the fluctuations in the phase space volume which allows for the exploration of the canonical measure. The precise relation to the fluctuations in a phase space volume $`𝒱`$, or equivalently, the instantaneous phase space compressibility, can be found using the divergence theorem

$$\frac{1}{𝒱}\frac{d𝒱}{dt}=_𝒱𝒟\phi \underset{k}{}\frac{dG(w_k)}{dw_k}\frac{dF(\pi _k)}{d\pi _k}+\frac{dG^{}(w_k^{})}{dw_k^{}}\frac{dF^{}(\pi _k)}{d\pi _k},$$
(11)
where the index $`k`$ runs over the thermostatted sites. In this paper, we do not explore the effect of different choices of $`G(w),G^{}(w),F(\pi ),F^{}(\pi )`$. Such studies have been done on various other systems .

Finally, we note that the linearized equations of motion are evolved by the stability matrix, $`\dot{\phi }_i/\phi _j`$. The eigenvalues of this matrix are the Lyapunov exponents for the system. Hence we have the relation that

$$\underset{i}{}\lambda _i=\underset{i}{}\frac{\dot{\phi }_i}{\phi _i}.$$
(12)
For the canonical ensemble, the Liouville equation gives $`\lambda _i=0`$, while in steady state non-equilibrium systems, $`\lambda _i<0`$.

One specific realization of the finite temperature equilibrium dynamics we use couples the thermostats only at the endpoints of the systems, $`k=1,L`$. With the choice of $`G(w)=w^4/4,G^{}(w^{})=w^2/2,F(\pi )=\pi /T,F^{}(\pi )=\pi ^3/T`$, we obtain

$$\begin{array}{ccc}\dot{\varphi }_k\hfill & =& \pi _kk=1,2,\mathrm{}L\hfill \\ \dot{\pi }_k\hfill & =& \{\begin{array}{cc}(^2\varphi )_km^2\varphi _k\varphi _k^3\hfill & k=2,3,\mathrm{},L1\hfill \\ (^2\varphi )_km^2\varphi _k\varphi _k^3w_k^3\pi _k/Tw_k^{}\pi _k^3/T\hfill & k=1,L\hfill \end{array}\hfill \\ \dot{w}_k\hfill & =& \pi _k^2/T1,\dot{w}_k^{}=\pi _k^4/T3\pi _k^2k=1,L.\hfill \end{array}$$
(13)
Just thermostatting the two boundary points is sufficient to thermalize the entire system. We then can examine the interior system far from the boundaries to study the finite temperature theory. Either free or fixed boundary conditions were used for the $`\varphi `$ field with no significant effect on the physics behavior of the theory. To compute observables, we use the fact that in this type of dynamics, the time averages converge to the ensemble average using the desired ensemble (8):

$$\overline{𝒪}=\underset{t\mathrm{}}{lim}\frac{1}{t}_0^t𝑑t𝒪(\varphi (t),\pi (t))=𝒪_{EQ}=\frac{𝒟\phi 𝒪f}{𝒟\phi f}.$$
(14)
### B Non-equilibrium Boundary Conditions: $`T0`$

One of the problems immediately encountered in the study of non-equilibrium systems is the nature of the steady state statistical distribution. While equilibrium statistical mechanics is well understood, once we apply thermal gradients to the system, much less is known about the system. To set up the non-equilibrium molecular dynamics, we use the demons to thermostat the end-points of our system in the same way we did in the equilibrium simulation. The only difference is that we now control the two endpoint temperatures separately. A consequence will be that the phase space distribution, $`f(\varphi ,\pi ,t)`$, will evolve to a non-smooth function that describes the non-equilibrium steady state.

We take the equations of motion at finite $`T`$ and now introduce two temperatures $`T_1^0`$ and $`T_2^0`$. The superscript is needed to distinguish the thermostatted temperatures from those measured just inside the system, which can suffer from boundary jumps which we will analyze in detail. One set of equations, similar to the equilibrium case in Eq. (13), we have used are

$$\begin{array}{cccc}\dot{\varphi }_k\hfill & =& \pi _kk=1,2,\mathrm{}L\hfill & \\ \dot{\pi }_k\hfill & =& \{\begin{array}{cc}(^2\varphi )_km^2\varphi _k\varphi _k^3\hfill & k=2,3,\mathrm{},L1\hfill \\ (^2\varphi )_1m^2\varphi _1\varphi _1^3w_1^3\pi _1/T_1^0w_1^{}\pi _1^3/T_1^0\hfill & \\ (^2\varphi )_Lm^2\varphi _L\varphi _L^3w_L^3\pi _L/T_2^0w_L^{}\pi _L^3/T_2^0\hfill & \end{array}\hfill & \\ \dot{w}_1\hfill & =& \pi _1^2/T_1^01,\dot{w}_1^{}=\pi _1^4/T_1^03\pi _1^2\hfill & \\ \dot{w}_L\hfill & =& \pi _L^2/T_2^01,\dot{w}_L^{}=\pi _L^4/T_2^03\pi _L^2.\hfill & \end{array}$$
(15)
One can see that the thermostats are applied to the endpoints $`k=1`$ and $`k=L`$. It should be noted that inside the boundaries, the dynamics of the system is that of the $`\varphi ^4`$ theory itself with no other degrees of freedom. We have also considered both variations of these thermostats, such as different forms of the interactions or increasing the number of sites at each end where we apply the demons. We will comment when these distinctions are relevant.

The equations of motion are solved on a spatial grid, using two methods: fifth and sixth order Runge-Kutta, and leap-frog algorithms. We used from $`10^6`$ to $`10^9`$ time steps of $`dt`$ from $`0.1`$ to $`0.001`$, with observables being sampled every $`\mathrm{\Delta }t=20100dt`$. The lattice size was varied from $`L=20`$ to 8000.

A consequence of the non-equilibrium steady state is that the measure becomes singular with respect to the Liouville measure . We do not rely explicitly on the singular nature of the non-equilibrium measure. Rather, we use it to interpret certain observables in the non-equilibrium state. To see this we start with the two main equations for $`f(\varphi ,\pi ,t)`$; the continuity equation and the expression for the total derivative of a phase space valued function. If we denote the vector $`\phi `$ to include all degrees of freedom, $`\phi =(\varphi ,\pi ,w)`$, we have

$$\frac{df}{dt}=\frac{f}{t}+\underset{i}{}\dot{\phi }_i\frac{f}{\phi _i}.$$
(16)
By combining this equation with the continuity equation (9), which holds both in equilibrium and in non-equilibrium, we derive

$$\frac{df}{dt}=f\underset{i}{}\frac{\dot{\phi }_i}{\phi _i}.$$
(17)
Solving for $`f`$, we obtain

$$f(t)=f(0)\mathrm{exp}\left(_0^t𝑑t\underset{i}{}\frac{\dot{\phi }_i}{\phi _i}\right)=f(0)\mathrm{exp}\left(t\underset{i}{}\frac{\dot{\phi }_i}{\phi _i}_{NE}\right)=f(0)\mathrm{exp}\left(t\underset{i}{}\lambda _i\right).$$
(18)
In these steps we have replaced the time average with the non-equilibrium ensemble average. The average of the divergence of the equations of motion is nothing more than the sum of Lyapunov exponents. In non-equilibrium steady states, whether from thermal gradients, or shearing, and so forth, the sum of the exponents is observed to become negative, signaling the presence of a fractal dimension. In the steady state,

$$f(t)t\mathrm{}\mathrm{}.$$
(19)
As the continuity equation is satisfied, the allowed phase space volume must be shrinking onto a set of measure zero, with respect to the original measure. A consequence of the divergence of the distribution function in the non-equilibrium steady state is that Gibbs entropy will also diverge,

$$S_G=\mathrm{log}f\mathrm{}.$$
(20)
Although we do not know how to properly define the fractal measure for our non-equilibrium steady state, we are able to compute non-equilibrium expectation values with respect to it using time averages:

$$𝒪_{NE}=\overline{𝒪}=\underset{t\mathrm{}}{lim}\frac{1}{t}_0^t𝑑t𝒪(\varphi (t),\pi (t)).$$
(21)
We will verify this in the linear response regime where near-equilibrium results can be compared to thermal equilibrium predictions obtained using the linear response theory.

Simulations of many-body systems in non-equilibrium steady states have found that the steady state measure is typically singular with respect to the original equilibrium measure. As one moves further from thermal equilibrium, the available phase space contracts onto an ergodic fractal: the accessible points are dense in the phase space, but fractal in nature. The resulting loss of dimension is related to the transport coefficient. One can see that in our dynamics as well since $`f`$ satisfies the continuity equation, so that total probability is conserved, yet Eq. (19) is satisfied. This means that the accessible phase space volume is contracted onto a set of measure zero. This type of ‘dimensional loss’ in the steady state can be more rigorously seen in low dimensional systems like the Lorentz gas.

From the point of view of dynamical systems theory, it is possible to understand the steady state measures under certain special conditions, namely that the dynamics is hyperbolic or Anosov. Unfortunately the conditions for hyperbolicity are not satisfied for our system since the number of positive Lyapunov exponents can vary along a trajectory. Nevertheless, it is useful to note that for hyperbolic systems described by some flow, $`x(t)=S_tx`$, for time evolution operator $`S_t`$ and initial point $`x`$, the Sinai-Ruelle-Bowen theorem provides the existence of the steady state measure, denoted $`\mu _{SRB}`$. It can be shown that for a continuous function $`f(x)`$, there exists a unique measure $`\mu _{SRB}`$ such that

$$\underset{t\mathrm{}}{lim}\frac{1}{t}_0^t𝑑t^{}f(x(t))=\frac{f(x)𝑑\mu _{SRB}}{𝑑\mu _{SRB}}.$$
(22)
This measure is known to be fractal for various non-equilibrium systems, and can be explicitly constructed for certain maps such as the modified Baker’s map. Here a basis of fractal Takagi functions has been implemented.

For the Lorentz gas, it was shown numerically under general conditions, and rigorously proven in the linear response regime, that the steady state measure is singular. Here one can explicitly see that the dimensional loss $`\mathrm{\Delta }D`$ is proportional to the transport coefficient, which in this case is the electrical conductivity. While there is still some contention concerning the existence of singular measures in the extension of the Gibbs ensemble to the far from equilibrium steady state, the compelling evidence suggests that whatever its nature is, it is far from trivial.

## III Transport

One of the relevant observables is the stress-energy tensor,

$$𝒯^{\mu \nu }=\frac{}{(_\mu \varphi )}_\nu \varphi +\eta _{\mu \nu }.$$
(23)
From the continuity equation $`0=_\mu 𝒯^{\mu \nu }`$, we have

$`0={\displaystyle \frac{}{x^\mu }}𝒯^{0\mu }={\displaystyle \frac{𝒯^{00}}{t}}+{\displaystyle \frac{𝒯^{0i}}{x^i}}.`$ (24)
We can identify the heat flux as $`𝒯^{0i}`$. On the lattice and in one spatial dimension,

$$𝒯^{01}(x)=\pi (x)\varphi (x)𝒯_k^{01}=\pi _k(\varphi _{k+1}\varphi _k).$$
(25)
Defining the lattice energy density consistently as

$$𝒯^{00}(x)=\frac{1}{2}\pi ^2+\frac{1}{2}(\varphi )^2+V(\varphi )𝒯_k^{00}=\frac{1}{2}\pi _k^2+\frac{1}{2}(\varphi _{k+1}\varphi _k)^2+V(\varphi _k),$$
(26)
we find that the discrete version of the continuity equation is satisfied:

$$\frac{}{t}𝒯_k^{00}+\left(𝒯_k^{01}𝒯_{k1}^{01}\right)=0.$$
(27)
This establishes that $`𝒯_k^{01}`$ is a constant in space and time in steady state. On the other hand, there is no way to satisfy the spatial component of the continuity equation, $`_t𝒯_k^{01}+𝒯^{11}=0`$. The reason for this can be understood as follows; while the translation invariance in the time direction is preserved, the translation invariance in the spatial direction has been lost due to the lattice discretization.

It is interesting to contrast this with the so called FPU $`\beta `$ model, which has divergent thermal conductivity in one spatial dimension. In this case the Hamiltonian is

$$H_{FPU}=\frac{1}{2}\underset{k}{}\left[p_k^2+(q_{k+1}q_k)^2+\frac{\beta }{2}(q_{k+1}q_k)^4\right],$$
(28)
so that the heat flux reads

$$𝒯^{01}=\pi _k(q_{k+1}q_k)\left(1+\beta (q_{k+1}q_k)^2\right).$$
(29)
If we think of this as originating from a continuum field theory, we would have the following (non-relativistic) Lagrangian and heat flux:

$$=\frac{1}{2}\left(\frac{\varphi }{t}\right)^2+\frac{1}{2}\left(\frac{\varphi }{x}\right)^2+\frac{\beta }{4}\left(\frac{\varphi }{x}\right)^4,𝒯^{01}=\frac{\varphi }{t}\frac{\varphi }{x}\left(1+\beta \left(\frac{\varphi }{x}\right)^2\right).$$
(30)
It is straightforward to check that the heat current satisfies the continuity equation (24). In this case, the thermal conductivity diverges as $`L^{2/5}`$ .

In the simulations we will define the temperature locally using the concept of an ideal gas thermometer, which measures the second moment of the momentum distribution. Specifically, at a site $`k(x=ka)`$ we define in both equilibrium and non-equilibrium

$$T(x)=T_k=\pi _k^2.$$
(31)
We will see that this definition is sensible in these regimes in the context of both statistical mechanics and thermodynamics.

### A Green-Kubo Approach

The standard approach to near equilibrium or local equilibrium transport is to apply linear response theory or Green-Kubo formulas. In this approach, external fields are added to a Hamiltonian which generate the desired transport processes and expressions for the transport coefficients can be derived in the weak field limit. For thermal conduction, the arguments are somewhat heuristic since one does not have well defined external fields that produce heat flow. Nevertheless, the procedure empirically results in sensible expressions. The Green-Kubo approach expresses transport coefficients in terms of equilibrium autocorrelation functions. For the thermal conductivity, one has

$`\kappa (T)`$ $`=`$ $`{\displaystyle \frac{1}{T^2}}{\displaystyle _0^{\mathrm{}}}𝑑t{\displaystyle 𝑑x𝒯^{01}(x,t)𝒯^{01}(x_0,0)_{EQ}}`$ (32)
$`=`$ $`{\displaystyle \frac{1}{NT^2}}{\displaystyle _0^{\mathrm{}}}𝑑t{\displaystyle \underset{k,k^{}=1}{\overset{N}{}}}𝒯^{01}(x_k,t)𝒯^{01}(x_k^{},0)_{EQ}.`$ (33)
$`N(<L)`$ is the number of sites in the region inside the boundaries used in the computation. We typically choose the region to be as large as possible while excluding the boundary effects. While this expression is expected to hold near equilibrium, the region of its validity cannot be determined within the linear response theory itself.

Autocorrelation functions, such as Eq. (32), have been argued to decay algebraically, rather than exponentially. Originally observed in early molecular dynamics simulations, the velocity-velocity autocorrelation function was found to decay as $`v(t)v(0)_{EQ}t^{d/2}`$ up to times on the order of a few ten times the mean free time. This behavior, if it persists, leads to divergence of the diffusion constants for $`d<3`$. Using kinetic theory and some general assumptions, it was later argued that this is a generic feature of dynamical systems for long times. These long-time power law tails cause divergences in transport coefficients in a large class of low dimensional systems including the FPU $`\beta `$ model (28)—(30). It is believed, however, that for theories without strict momentum conservation, such a divergence can be absent. This is the case for our lattice model due to our ‘on-site’ nature of the potential, as we shall see below.

### B Non-Equilibrium Approach

We also compute the thermal conductivity directly by constructing ensembles near equilibrium which lead to constant gradient thermal profiles. This is arguably more fool-proof than the Green-Kubo formula since no assumptions are necessary for the computation. Here we first confirm Fourier’s law and then use it to compute the conductivity using

$$\kappa (T)=\frac{𝒯^{01}_{NE}}{T(x)}.$$
(34)
While Fourier’s law is believed to be valid near equilibrium, what constitutes the ‘linear regime’ is not known. We will see that it eventually breaks down as we move far from equilibrium, while still in the steady state. It is not clear at this point how to make a sensible definition of the thermal conductivity. One approach would be to try to the characterize the ‘nonlinear’ response by incorporating the dependence of the conductivity on higher order terms in the derivatives, such as $`(T)^3,(^2T)(T)`$ and so on. We will try to clarify these issues below.

## IV  Equilibrium Ensemble

When the boundary temperatures are equal, we recover equilibrium physics for our system. In the remainder of the paper we specialize to the $`m^2=0`$ case.

### A Thermalization

A test of the thermalization of the system is readily performed by studying the distribution functions of various quantities. If we take a single trajectory, $`(\varphi (t),\pi (t))`$, and histogram $`\pi _k(t)`$ at $`t=0`$, $`\mathrm{\Delta }t`$, $`2\mathrm{\Delta }t`$,…, where $`\mathrm{\Delta }t`$ is some time step which ensures the points are reasonably well decorrelated, we will converge to the thermal distribution function

$$f(\pi _k)\mathrm{exp}[\pi _k^2/2T].$$
(35)
In Fig. 1 (left column), we show computed equilibrium distributions for the momenta and heat flux at the center of a lattice with $`L=163`$. Here we have taken the endpoint temperatures to be $`T_1^0=T_2^0=1`$. The measured temperature in the middle is $`T=0.995(8)`$, and one can see that the measured momentum distribution (histogram, top left) agrees with the predicted value (solid, top left). In the bottom left we show the measured thermal distribution of $`𝒯^{01}`$, for which we do not have a theoretical prediction. Since we are in equilibrium, there is no heat flow, so $`𝒯^{01}=0`$, which is evident from the symmetric nature of the thermal distribution $`f(𝒯^{01})`$. We have verified that when $`T_1^0=T_2^0`$, these boundary conditions reproduce the equilibrium canonical measure $`f_{eq}(\pi ,\varphi )\mathrm{exp}[H(\pi ,\varphi )/T_1]`$ at all points, including the thermostatted sites. The measured thermal profiles satisfy $`T(x)=T_1^0(=T_2^0)`$ for any $`x`$ within numerical error over the temperature range we investigated: $`T=0.01`$ to $`T=10`$. Hence the boundary conditions (15) do produce the desired physics for this theory.

In our presentation of equilibrium and non-equilibrium steady state expectation values, we verify that stationary results are obtained. This is done by examining the time evolution of observables. For instance, if we want to measure the heat flow through the system, we measure at all sites and verify that the time averages converge to the same value. A typical result is shown in Fig. 2, where the average $`𝒯^{01}`$ is shown as a function of time at three sites: $`L/4`$, $`L/2`$ and $`3L/4`$ for $`L=800`$. One can see that the values eventually converge to the ensemble average. In this case the endpoint temperatures are distinct, so that there is a net heat flux. Equilibrium figures are similar, with convergence to $`𝒯^{01}=0`$.

### B Auto-Correlation Functions

We compute auto-correlation functions for the components of the stress-energy tensor, which includes the Green-Kubo formula for the thermal conductivity. We define the normalized correlation functions

$$C(t)=\frac{1}{(𝒯^{\mu \nu }(0))^2}\left[𝒯^{\mu \nu }(t)𝒯^{\mu \nu }(0)𝒯^{\mu \nu }(0)^2\right].$$
(36)
In Fig. 3, we show the time dependence of $`C(t)`$ for $`𝒯^{00}`$ (dashes), $`𝒯^{01}`$ (solid) and $`𝒯^{11}`$ (dots) as a function of time, where the time is normalized by the mean free time $`\tau `$. (As we discuss below, $`\tau `$ will be of order of the thermal conductivity.) One can see that these functions decay to half their initial values on the order of the mean free time. For larger times $`t>10\tau `$, these functions oscillate about zero. This behavior is seen at all temperatures studied.

To compute the thermal conductivity, we are interested in the auto-correlation function for $`𝒯^{01}`$:

$$\kappa (T,t)=\frac{1}{NT^2}_0^t𝑑t\underset{k,k^{}}{}𝒯^{01}(x_k,t)𝒯^{01}(x_k^{},0)_{EQ},\kappa (T)=\underset{t\mathrm{}}{lim}\kappa (T,t).$$
(37)
By studying the convergence $`\kappa (T,t)`$ to the thermal conductivity, we can explore the problems associated with the long-time tails. In Fig. 4 we plot $`\kappa (T,t)`$ as a function of time, where $`t`$ is normalized by the mean free time $`\tau `$. According to the predictions based on kinetic theory, in $`d=1`$ we would expect $`\kappa (T,t)t^{1/2}`$, leading to an infinite conductivity. At the temperatures shown ($`T=1/50`$, 1/10, 1, 2), $`\kappa (T,t)`$ does apparently display behavior similar to $`t^{1/2}`$ (dashes) on time scales up to $`t10\tau `$, but on longer time scales, the results converge to well defined values. Consequently, divergences such as those found in the FPU model are not present here, and long-time tails are at best a transient aspect of this dynamics.

The behavior of $`\kappa (T)`$ found from the Green-Kubo approach is summarized in Fig. 5 (crosses). This analysis will be continued when we discuss the direct measurements below.

### C Speed of Sound

To better understand the kinetic theory aspects of this finite temperature theory, we would like to know the thermal behavior of the ‘speed of sound’. This can be estimated in a number of ways. A convenient approach is to use the $`𝒯^{01}`$ auto-correlation function. The sound speed, $`c_s`$, is defined here as how fast excitations travel through the system and is the relevant velocity for the transport theory in our particular model. We note, however, that strictly speaking, this is not ‘sound’ in the hydrodynamic sense. We define

$$G(x,t;x_0,t_0)𝒯^{01}(x_0,t_0)𝒯^{01}(x,t)_{EQ}.$$
(38)
Instead of the time dependence, we consider the spatial dependence of $`G(x,t;x_0,t_0)`$. Not only will it decorrelate in time, as we saw previously in the Green-Kubo integrals, but it will also decorrelate over space. A typical behavior of $`G(x,t;x_0,t_0)`$ is shown in Fig. 6. Here we choose $`x_0=0`$ to be the center of the lattice, and $`t_0=0`$. The temperature is $`T=1/10`$ and $`L=160`$. The autocorrelation function is shown for times $`t=0,30,60,90`$ where one can see the regions of high correlation separate. By measuring the rate at which the peaks separate, we can obtain an estimate for $`c_s`$ at that temperature. In Fig. 7, we plot the sound speed extracted in this manner against the temperature. There is little temperature dependence until $`T1/10`$, after which the speed begins to decrease with $`T`$, as one might naively expect. So although it shows some temperature dependence, $`c_s`$ is generally of order unity over this temperature range.

## V  Non-Equilibrium Steady States

By controlling the boundary temperatures, we can begin to develop an understanding of the physics of hot scalar field theory, as it moves increasingly further from thermal equilibrium. We discuss the regimes we have found as we move away from equilibrium.

### A Near Equilibrium: $`T_1T_2`$

In the near equilibrium limit, the physics is consistent with linear response and Fourier’s law. Here $`𝒯^{01}`$ will develop a non-zero expectation value. As the temperatures $`T_1`$ and $`T_2`$ begin to separate, a constant $`T`$ develops. This is illustrated in Fig. 8 for a lattice of $`L=8000`$ where a small gradient is developed around a temperature $`T=1/4`$ using boundary temperatures $`T_1^0=0.15`$ and $`T_2^0=0.35`$ (solid). A linear profile, expected from Fourier’s law, is superimposed (dashes). While such linear profiles in principle provide a direct measure of the conductivity by dividing the measured heat flux by the average gradient, we perform a more careful analysis by taking increasingly smaller gradients about the same average temperature, thereby explicitly verifying Fourier’s law. This is shown in Fig. 9 for the temperature $`T=1/4`$. All these points correspond to thermal profiles which are linear. As the boundary temperatures approach $`T=1/4`$, we see that the heat flux also decreases. The slope, $`\kappa (T)=𝒯^{01}/T`$, corresponds to the thermal conductivity at this temperature.

The measurements of $`\kappa (T)`$ in the range $`1/100<T<20`$ are summarized in Fig. 5 for both direct measurements ($`\mathrm{}`$) and Green-Kubo ($`\times `$). We find that both methods agree and that the general behavior is well described by a power law:

$$\kappa (T)=\frac{A}{T^\gamma },\gamma =1.35(2),A=2.83(4).$$
(39)
This type of behavior is reminiscent of lattice phonons at high temperature. The fact that the direct measurements agree with the Green-Kubo integrals of the auto-correlation functions is sufficient to dispel the notion of asymptotic long-time tails in this system.

Another important aspect to verify is that we have achieved a bulk limit in the values of $`\kappa (T)`$ reported. In Fig. 10, we show the dependence of the direct measurements of $`\kappa `$ on the size $`L`$ of the system for several temperatures. The dashed lines are the predictions from the power law fit. One can see that a bulk limit is realized for relatively small lattices. On closer inspection, we find that the bulk limit is reached for smaller lattices when the temperature is higher. This can be understood as follows; the mean free path of the system is of the order of the thermal conductivity, as we shall see when we analyze the kinetic theory aspects of this problem. The bulk limit is reached for lattices larger than the mean free path, roughly speaking.

### B Far from Equilibrium: $`T_1T_2:`$

We now consider the behavior of the system as we move it further from equilibrium. One of the first characteristics to emerge is the development of curvature in the temperature profile. In Fig. 11 we plot a succession of steady state temperature profiles as we change the boundary temperatures from $`(T_1^0,T_2^0)=(0.3,0.7)`$ (dots) to $`(0.2,0.8)`$ (dashes) and finally to $`(0.1,2)`$ (solid). As the system moves away from the constant gradient profiles, it starts to feel the temperature dependence of the thermal conductivity, which decreases with increasing temperatures. As a consequence, the hot end of the system cannot conduct heat as well as to lower temperature end. (Of course the converse of this argument would hold if the power law behavior has $`\gamma <0`$, in which case the curvature would have an opposite sign.)

Fig. 2, which was already mentioned when we discussed the convergence properties, displays the time evolution of the heat flux at three sites within a system thermostatted at $`(T_1^0,T_2^0)=(0.3,0.7)`$ (dots in Fig. 11). Regardless of how far the system is from equilibrium, the heat flux must be independent of position in steady state, since there are no sources or sinks for heat inside the boundaries. We see that the values at the three sites converges to the same value in the steady state, namely, a constant non-equilibrium heat flux independent of $`x`$.

Another property of the non-equilibrium system is the appearance of jumps in the temperature at the boundaries. This is illustrated in Fig. 12; in Fig. 12 (top), we show a typical non-equilibrium thermal profile with curvature (solid — the dashes will be discussed below). In the lower panels we examine the low and high temperature ends. $`T_1^0`$ and $`T_2^0`$ are the temperatures enforced by our boundary temperatures. One can see that there is a difference between these temperatures and what one obtains by smoothly extrapolating the temperature profile near the edge. These are physical phenomena associated with the dynamics at the interface, which can be readily understood quantitatively, as discussed in section V D below. We will define the points obtained through smooth extrapolation by $`(T_1,T_2)`$. If we now focus on the thermal profile away from the edges, we can understand the curvature in the temperature profile using Fourier’s law. At or near equilibrium, we have found the power law behavior $`\kappa (T)=AT^\gamma `$, as in Eq. (39). If we integrate Fourier’s law and re-express the result in terms of extrapolated temperatures $`(T_1,T_2)`$, we find

$$T(x)=\{\begin{array}{cc}T_1\left[1\left(1\left(\frac{T_2}{T_1}\right)^{1\gamma }\right)\frac{x}{L}\right]^{\frac{1}{1\gamma }},\hfill & \gamma 1\hfill \\ T_1\left(\frac{T_2}{T_1}\right)^{x/L},\hfill & \gamma =1.\hfill \end{array}$$
(40)
Some previous efforts to model the temperature profiles in systems exist , although these were generally model fits and furthermore, a full understanding requires a description of the boundary effects, which up to now have been lacking. We have found that this formula provides an excellent description, but ultimately breaks down as we see below. In Fig. 12, the dashed line is the description given by Eq. (40) with $`\gamma =1.35`$, which offers a fit to within a few percent. In should be noted that $`\gamma `$ used here was obtained from systems at or near equilibrium as in Eq. (39), independently of the systems far from equilibrium. One can see that at the high temperature end there is a tendency to overshoot the measured behavior as in Fig. 12 (bottom right). This analytic expression provides a sufficiently good description of the physics as we move away from equilibrium, until Fourier’s law fails to hold locally. This includes the first three regimes listed in Table 1. Eventually, Fourier’s law will break down locally as well, as we develop steady states which are locally non-equilibrium (LNE).

The analytic understanding of the thermal profile allows us to understand the behavior of the heat flux for systems not too far from equilibrium. Using Fourier’s law and Eq. (40), we have for

$$𝒯^{01}_{NE}=\frac{A}{L(1\gamma )}(T_1^{1\gamma }T_2^{1\gamma })(\gamma 1).$$
(41)
For $`\gamma =1`$, $`𝒯^{01}_{NE}=(A/L)\mathrm{log}(T_1/T_2)`$. When compared to measured heat flux using the extrapolated temperatures from the thermal profile, this formula is found to provide a very good description of the heat flux, typically to within a few percent. In the near equilibrium regime where the temperature profile is visibly linear, the boundary jumps vanish, and we expect it to behave as

$$𝒯^{01}_{NE}^0=\kappa (T_{av})\frac{T_2T_1}{L}$$
(42)
where the superscript denotes the linear response, or constant gradient limit. In general, these are related by

$$𝒯^{01}_{NE}=𝒯^{01}_{NE}^0\left\{1+\frac{\gamma (\gamma +1)}{24}\left(\frac{\mathrm{\Delta }T}{T_{av}}\right)^2+𝒪\left(\left[\frac{\mathrm{\Delta }T}{T_{av}}\right]^4\right)\right\},$$
(43)
where $`T_{av}=(T_2+T_1)/2`$ and $`\mathrm{\Delta }T=(T_2T_1)`$. Note that these temperatures are the extrapolated ones and not the boundary temperatures, $`(T_1^0,T_2^0)`$. In this way we see that the temperature profile will no longer be linear when $`\mathrm{\Delta }T/T_{av}LT(x)/T(x)1`$.

In this regime of curved thermal profiles, we see some differences between the thermal distributions measured in equilibrium. In Fig. 1 we display on the right side the measured steady-state non-equilibrium statistical distributions for momenta (histogram, upper right) compared to the equilibrium thermal distribution expected at that temperature based on the notion of local equilibrium (solid). We find that as we more further from equilibrium, the momentum distributions are typically sharper than the gaussian expected based on the ideal gas thermometer. In the lower right, we show the non-equilibrium distribution of heat flux which develops the asymmetry needed to give a non-vanishing expectation value $`𝒯^{01}<0`$ for $`T_1^0<T_2^0`$.

An illustration of a system very far from equilibrium is the non-equilibrium steady state shown in Fig. 13. Here the ratio of the endpoint temperatures is 100, and the thermal profile is shown on a log scale (solid). The fit of Eq. (40) is shown by the dashed line. While we still measure temperature locally through the second moment, the measured distributions are no longer gaussian. In Fig. 14 (top) we show $`f(\pi )`$ in the center of the system. In the bottom panel, the ratio of Eq. (40) to the measured behavior in Fig. 13 is shown. Deviations up to 40% are evident on the low temperature end. Here, there is no accepted formalism to describe the dynamics and we are in the LNE regime of Table 1. We will return to this issue when we discuss local equilibrium in more detail in section V F.

### C Entropy and Kinetic Theory

While the Gibbs entropy for the theory diverges as it contracts onto a set of measure zero, not all measures of entropy behave in this manner. Irreversible thermodynamics provides a description of systems on scales much larger than the mean free path. When one has heat flow, one defines the local rate of entropy production as

$`\sigma (x)`$ $`=`$ $`𝒯^{01}_{NE}{\displaystyle \frac{d}{dx}}{\displaystyle \frac{1}{T(x)}}`$ (44)
$`=`$ $`{\displaystyle \frac{A}{L^2(1\gamma )^2T_1^\gamma }}\left[1\left({\displaystyle \frac{T_2}{T_1}}\right)^{1\gamma }\right]^2\left[1\left(1\left({\displaystyle \frac{T_2}{T_1}}\right)^{1\gamma }\right){\displaystyle \frac{x}{L}}\right]^{\frac{(2\gamma )}{(\gamma 1)}}>0.`$ (45)
Integrating this formula, we can compute the net rate of entropy production to be

$`\dot{S}_{irr}`$ $`=`$ $`{\displaystyle _0^L}𝑑x\sigma (x)`$ (46)
$`=`$ $`𝒯^{01}_{NE}\left({\displaystyle \frac{1}{T_2}}{\displaystyle \frac{1}{T_1}}\right)`$ (47)
$`=`$ $`{\displaystyle \frac{A}{L(\gamma 1)T_1^\gamma }}\left(1{\displaystyle \frac{T_1}{T_2}}\right)\left[1\left({\displaystyle \frac{T_1}{T_2}}\right)^{\gamma 1}\right]>0.`$ (48)
This expression can be simply interpreted as saying that the global entropy production rate is due to the difference in entropy production rates at the boundaries coming from the demons. On the other hand, irreversible thermodynamics predicts a constant local rate of entropy production, $`\sigma (x)`$, which is at odds with the coarse-grained local Boltzmann entropy computed below. The latter calculation is microscopically based and does not rely on any hydrodynamic limit of the theory.

If we envision the expansion of hot systems similar to RHIC collisions, the local entropy is an important quantity. From the statistical mechanics point of view, when one discusses the notion of entropy in a local frame, it is not Gibbs entropy of the entire systems that is needed. Hence we would like to consider Boltzmann’s entropy. To understand the behavior of this entropy in increasingly non-equilibrium environments, we consider the $`n`$-body Boltzmann entropy, defined through the $`n`$-body distribution functions $`f_B^{(n)}`$. The best we can do, of course, is the coarse-grained limit of these quantities. In the ‘local frame’ at $`x=x^{}`$, we define them as integrals of the full phase space distribution $`f`$ over all quantities except the arguments of $`f_B^{(n)}`$:

$`f_B^{(1)}(\varphi (x^{}),\pi (x^{}))`$ $`=`$ $`f_B^{(1)}(\varphi _k,\pi _k)`$ (49)
$`f_B^{(2)}(\varphi (x^{}),\varphi (x^{\prime \prime }),\pi (x^{}),\pi (x^{\prime \prime }))`$ $`=`$ $`f_B^{(2)}(\varphi _k,\varphi _{k+1},\pi _k,\pi _{k+1})`$ (50)
$`\mathrm{}`$ . (51)
The distributions $`f_B^{(k)}`$ are obtained by histogramming the corresponding degrees of freedom ( $`(\varphi _k,\pi _k)`$ for $`f_B^{(1)}`$, $`(\varphi _k,\pi _k,\varphi _{k+1},\pi _{k+1})`$ for $`f_B^{(2)}`$,…), and is readily constructed in equilibrium and non-equilibrium conditions. In the equilibrium or non-equilibrium steady states, we then compute

$`S_B^{(1)}`$ $`=`$ $`{\displaystyle 𝑑\pi (x)𝑑\varphi (x)f_B^{(1)}\mathrm{log}f_B^{(1)}}`$ (52)
$`S_B^{(2)}`$ $`=`$ $`{\displaystyle 𝑑\pi (x)𝑑\varphi (x)𝑑\pi (x^{})𝑑\varphi (x^{})f_B^{(2)}\mathrm{log}f_B^{(2)}}`$ (53)
$`\mathrm{}`$ (54)
These satisfy the inequalities

$$S_G\mathrm{}\frac{L}{2}S_B^{(2)}LS_B^{(1)}.$$
(55)
To compute $`S_B`$ we must coarse-grain the phase space. Hence we are actually computing the coarse-grained 1- and 2-body densities, which we denote $`f_\mathrm{\Delta }^{(k)}`$. These are related by

$$S_B^{(1)}f_\mathrm{\Delta }\mathrm{log}f_\mathrm{\Delta }\mathrm{log}(\mathrm{\Delta }\pi _k\mathrm{\Delta }\varphi _k),f_\mathrm{\Delta }=1.$$
(56)
We have computed these entropies and find that $`S_B^{(1)}`$ does not shift noticeably from its equilibrium value regardless of how far the system is from equilibrium (see Fig. 16). Further, $`S_B^{(2)}`$ $`(2S_B^{(1)})`$ is only slightly less than its upper limit $`2S_B^{(1)}`$ and remains so even far from equilibrium. So unlike $`S_G`$ $`\left(LS_B^{(1)}\right)`$, $`S_B`$ is rather insensitive to the non-equilibrium nature of the system. Since this is coarse-grained, it is not too surprising that this is so. It would be more revealing to consider the behavior of $`f_B^{(k)}`$ as the size of the bins decreases, but this is currently too computationally intensive.

Thermodynamic quantities, such as the entropy, allow us to investigate the underlying dynamics of the theory and also probe the possible deviations of the physical observables from their equilibrium values under non-equilibrium conditions. We shall analyze these questions below. Let us first obtain the specific heat, $`C_V`$, from the equilibrium ensembles; $`C_V`$ may be obtained using the standard formula,

$$C_V=\frac{E^2_{EQ}E_{EQ}^2}{T^2},$$
(57)
where $`E`$ is the energy per site. We find that the specific heat has a weak temperature dependence which may be fitted by a simple temperature dependence,

$$C_V=C_0T^\alpha ,C_0=0.86(2),\alpha =0.025(6).$$
(58)
This entails that the energy of the lattice per site has the following behavior

$$E=\frac{C_0}{1\alpha }T^{1\alpha }\alpha 0E=C_0T.$$
(59)
Since the temperature dependence is weak, it is natural to compare $`C_V`$ against $`E/T`$ obtained locally both in equilibrium and in non-equilibrium. Such a comparison of $`C_V`$, $`E/T`$ for equilibrium, near equilibrium and far from equilibrium is shown in Fig. 15. We find that they all agree quite well. There seems to be an intriguing tendency for the non-equilibrium values of $`E/T`$ to be larger than their equilibrium counterparts. However, this cannot be delineated within the errors and more investigation is necessary to see if this is a real effect.

Let us now move onto the computation of the Gibbs entropy, $`S_G`$, in equilibrium. From the application of the first law of thermodynamics, we obtain

$$S_G=^T\frac{C_V}{T}𝑑T.$$
(60)
$$S_G(T)=S_{G,0}\frac{C_0}{\alpha }\left(T^\alpha 1\right)\alpha 0S_{G,0}+C_0\mathrm{ln}T.$$
(61)
In Fig. 16 we plot the computed Boltzmann entropy $`S_B^{(1)}`$ as a function of temperature for systems in equilibrium, as well as for systems near and far from equilibrium and find a similar form, $`S_B^{(1)}=S_{B,0}+S_{B,1}\mathrm{log}T`$. While we cannot measure $`S_G`$ away from equilibrium, $`S_{B,1}C_0`$ so that the behavior of $`S_B^{(1)}`$ in equilibrium and in non-equilibrium is similar to that of $`S_G`$ in equilibrium. We also find that $`S_B^{(1)}`$ is completely insensitive to the departure from equilibrium. Consequently, Boltzmann’s entropy and some of the thermodynamic concepts such as temperature, still retain their relationships in these local states far from equilibrium, even when the momentum distributions are no longer gaussian. Other physical quantities were also studied for effects of the system being in non-equilibrium, including $`𝒯^{11}`$,

$$𝒯^{11}=\frac{1}{2}\pi ^2+\frac{1}{2}\left(\varphi \right)^2\frac{\varphi ^4}{4}.$$
(62)
We find that the results for $`𝒯^{11}`$ in equilibrium, near equilibrium and far from equilibrium are consistent with each other in a manner similar to $`𝒯^{00}`$ in Fig. 15.

### D Boundary Temperature Jumps

Boundary temperature jumps in non-equilibrium steady states, such as systems experiencing thermal gradients are well known, although it seems that their behavior has never had a suitable explanation. Systems experiencing shearing due to moving walls also display jumps in the velocity profile, with the fluid inside the wall moving slightly slower than the wall velocity, for large velocities. Such effects are known to be sources of error in experiments that measure transport coefficients . We find that it is possible to achieve a quantitative understanding of these effects using simple kinetic arguments.

We have found in this system that $`c_s`$ and $`C_V`$ are of order unity, and have at most a weak temperature dependence (see Fig. 7 and Fig. 15). The thermal conductivity is related to the mean free path by $`\kappa C_Vc_s\mathrm{}`$ by a standard kinetic theory argument. Hence we expect $`\kappa \mathrm{}`$. Strictly speaking, $`\mathrm{}`$ is the mean free path in equilibrium, but we will refer to this loosely as the mean free path also away from equilibrium, to use as a natural length scale in the system. It corresponds to the mean free path when the thermal gradients are not too strong (up to LE-II regime in Table 1), but the standard kinetic theory arguments presumably break down when local equilibrium no longer holds. The boundary temperature jumps are due to the mean free path being non-zero and the deviation of the system from equilibrium. When $`\mathrm{}L`$,

$$T_i^0T_i=\eta \frac{T}{n}|_{boundary},$$
(63)
where $`n`$ denotes the normal to the boundary . This formula should apply when the jump $`|T_i^0T_i|/T_i`$ is relatively small. On dimensional grounds, the coefficient $`\eta `$ should be on order of the mean free path $`\mathrm{}`$, so that $`\eta \kappa `$. We have verified this relation by plotting $`\eta `$ measured directly from thermal profiles, as a function of the extrapolated temperature $`T`$ at the boundary. These are summarized in Fig. 18 (top). In the figure we display only the jumps on the low temperature end since the errors are much smaller. (The large gradient data on the high temperature end, while consistent with the low end data, is quite noisy.) We also show a fit to the data (dashes), which gives

$$\eta (T)=(6.1\pm 0.5)T^{1.5\pm 0.1}.$$
(64)
On the same figure we also show the behavior of $`\kappa (T)`$, which indicates that $`\eta \kappa `$ is consistent with the observed behavior.

An independent verification of the behavior of the jumps can be made by studying how the jumps depend on the heat flux. We let $`\eta =\alpha \kappa `$, where $`\alpha `$ is a constant to be determined. From Eq. (63), we can then associate the heat flux with the right side of Eq. (63):

$$T_iT_i^0\alpha 𝒯^{01}_{NE}\alpha (T_2^0T_1^0)\frac{\kappa (T_{av})}{L}+\mathrm{}.$$
(65)
By plotting the boundary jumps directly as a function of the heat flux in non-equilibrium steady states both near and far from equilibrium, we obtain the data in Fig. 18 (bottom). One can see that there is a simple relationship, where the slope gives $`\alpha =2.6(1)`$, consistent with the assumption that $`\eta =\alpha \kappa `$. The understanding of these jumps together with that of the temperature profile (40) provides a complete description of $`T(x)`$ in terms of the boundary temperatures, $`(T_1^0,T_2^0)`$. As was alluded to earlier, these simple relations (63),(65) break down in the LNE region of Table 1, where the thermal gradients are too strong.

A note on boundary conditions is in order: There is reasonable freedom in how the thermostats are implemented. We can vary the number of thermostatted sites or how the demons are couple to the physical degrees of freedom. We find that the changes in the way thermostats are implemented bring about changes in the boundary jumps when we are sufficiently far from equilibrium. This will also result in a different $`𝒯^{01}`$ in a manner consistent with Eq. (41). Different thermostats will correspond to different values of $`\alpha `$, which is not an intrinsic parameter of the $`\varphi ^4`$ theory but rather reflects different manners in which a heat bath might efficiently couple to the scalar field theory. In other words, as expected, different thermostats do not change our understanding of the physics at all.

In understanding the generality of Eq. (65), we note that the FPU $`\beta `$ model (28) displays temperature profiles that depend on $`L`$. In Refs. , it has been shown that the heat flux varies with the size as $`𝒯^{01}L^{0.55}`$. Consequently, one would expect that the boundary jumps seen there behave as $`\delta TL^{0.55}`$. A cursory analysis of those thermal profiles suggests that this is indeed the case.

### E  Non-Equilibrium Distributions

It is worth making a few remarks on the non-equilibrium distribution functions we have seen. Many usual approaches to the non-equilibrium statistical properties of systems make from the outset certain model dependent choices. In some thermo-field approaches, the temperature profile $`T(x)`$ is chosen to have a specific form, resulting in a statistical operator of the type $`\mathrm{exp}(\stackrel{~}{H}(\pi ,\varphi )/T(x))`$, where $`\stackrel{~}{H}(\pi ,\varphi )`$ might include contributions due to transport. Employing the approach of Zubarev, one assumes a particular behavior for the steady state distribution, which is then used to solve the dynamical equations.

The local non-equilibrium distribution, $`f_B^{(1)}(\varphi _k,\pi _k)`$ is shown in Fig. 17 for a system very far from equilibrium. While the distribution $`f`$ in the full configuration space is expected to be fractal, $`f_B^{(1)}(\varphi _k,\pi _k)`$ is only a projection on a lower dimensional space, and should be smooth. Such a non-equilibrium statistical operator will differ from those assumed in non-equilibrium approaches in their correlations. One aspect of the fractal nature of $`f`$ is that the dynamical space is of reduced dimension, so that additional correlations will exist that are not present in those standard non-equilibrium approaches.

Another limit in which one obtains unusual thermal profiles is a ‘ballistic’ limit in which the mean free path is on the order or larger than the size of the system. Because the excitations pass through the system rather readily, one has large boundary jumps on both ends. As a consequence, the thermal profile is almost flat, in spite of the values of the endpoint temperatures. In this case the system is still thermalized; the distributions of momenta at the endpoints and inside are gaussian. These types of results are included in our figures of entropy, specific heat and so forth, in the near-equilibrium data set, and they fall in line with non-ballistic results.

### F Local Equilibrium

While we do not yet have any precise criteria regarding when local equilibrium fails to be a good approximation, we can analyze various physical quantities and compare them to their values when local equilibrium holds. This should provide us with a measure of the dependence of the physics on local equilibrium. While field dependent observables, such as moments of $`\varphi (x)`$ will depend on the model under investigation, the momenta should be more robust. A natural place to start is with the analysis of cumulants. If local equilibrium is expected holds, we should have

$$\pi _k^2=T,\pi _k^4=3T^2,\pi _k^6=15T^3$$
(66)
and so forth. Equivalently, we can say that

$$\pi _k^4=\pi _k^43\pi _k^2^2,\pi _k^6=\pi _k^615\pi _k^4\pi _k^2+30\pi _k^2^3.$$
(67)
vanish in local equilibrium. Hence one measure of the breaking of local equilibrium when we are far from equilibrium is

$$\frac{\pi _k^4}{3\pi _k^2^2}=\frac{\pi _k^4}{3\pi _k^2^2}1.$$
(68)
Similar expressions for higher cumulants are also possible.

Non-equilibrium behavior begins to emerge in certain quantities as soon as $`T_1T_2`$. For instance, $`\pi ^4(x)/\pi ^2(x)^23`$, the equality holding in equilibrium, and the value growing as one departs from equilibrium. One can also see that the non-equilibrium measure is not locally Boltzmann since quantities such as $`\pi (x)\varphi (x^{})0`$ for $`xx^{}`$, as is expected for a system supporting some type of transport. However, the origins in this case we attribute not to additional terms in the statistical distribution $`f`$, but rather to additional correlations in non-equilibrium measure due to its reduced dimensionality, as expected from dynamical systems theory. Unfortunately, there is no way to explicitly estimate the dimensional loss without direct measurement of the entire Lyapunov spectrum, which is a rather numerically intensive task.

As one moves away from equilibrium, the linear response theory is expected to eventually breakdown. Such behavior is confirmed for our system; when the temperature gradient becomes too large, the formula for the heat flow (41) ceases to be valid. We display the relative difference of the measured current to the current obtained from the linear response theory (41) in Fig. 19 (bottom). It can be seen that the relative deviations can be of order one for large thermal gradients signaling the breakdown of the linear response theory. The deviations are plotted against $`\kappa (T)T/T`$ which seems to be the natural scale, since $`\kappa `$ is the mean free path, roughly speaking, which is the natural length scale in the problem as discussed in section V D. One obvious possibility would be to interpret this as a nonlinear response of the thermal conductivity, $`\kappa (T,𝒯^{01})`$, which may, for instance, be parametrized as

$$\kappa (T,𝒯^{01})=\kappa _0(T)+\kappa _2(T)\left(𝒯^{01}\right)^2+𝒪\left(\left(𝒯^{01}\right)^4\right).$$
(69)
Such approaches have been discussed in the literature .

However, such an interpretation presupposes, perhaps tacitly, that the concept of local equilibrium holds and that the standard notion of temperature applies, amongst other things. Therefore, it is imperative to first check that the local equilibrium is achieved in this ‘non-linear’ regime and that it is this point that we now wish to investigate with care. Such questions have been asked previously and in those situations, the local equilibrium was seen to be valid . First, we need to ask ourselves what constitutes local equilibrium? The concept of local equilibrium has been defined previously  and it is not our intention here to discern the possible differences in the various definitions of local equilibrium. In the least, local equilibrium assumes that we have a Maxwellian distribution for the momentum for the class of Hamiltonians we work with, leading to the usual concept of temperature, which will be our point of investigation.

In Fig. 19 (top), we plot the fourth cumulant of $`\pi `$, (68), against $`\kappa (T)T/T`$. The cumulant, which is defined locally, quite clearly deviates from the equilibrium value under strong thermal gradients. Even in these situations, the thermostatted boundary sites are in local equilibrium, as they should be, as shown in Fig. 20. It can be seen that the deviations from local equilibrium start to occur at roughly the same value of $`\kappa (T)T/T1/10`$ where linear response theory breaks down. In Table 1, we identify this as a steady state which is locally non-equilibrium (LNE). We have verified that higher cumulants display similar behavior. Therefore, at least in our model, the breakdown of local equilibrium needs to be considered when nonlinearity of the response is to be analyzed. Of course, this does not preclude the possibility of the nonlinearity of the response as discussed above, but that the nonlinearity needs to be disentangled from the deviations from local equilibrium with care.

We emphasize here that à priori, this needs not be the case. Namely, it is in principle possible that there is a region where the concept of local equilibrium is still valid and yet the linear response theory breaks down. In such a situation, ‘non–linear response’ theory should be quite appropriate for analyzing the situation. However, in the cases we studied, such regimes do not exist and the failure of the linear response theory occurs simultaneously with the breakdown of local equilibrium.

## VI Conclusions

We have constructed non-equilibrium steady states for classical $`\varphi ^4`$ lattice field theory in one dimension, under conditions near and far from equilibrium. We obtained the behavior of the thermal conductivity with respect to the temperature in the linear regime and found that the results were consistent with the linear response regime. The underlying dynamics of the theory was investigated and physical quantities such as the speed of sound, heat capacity, Boltzmann’s entropy and their temperature dependence, were obtained. The results could consistently understood using the kinetic theory approach. This understanding was further used to clarify the dynamics behind the temperature gaps that arise at the boundaries of the system. We also found that for temperature gradients that are not too large, the linear response law is adequate for understanding the behavior of the system, even though the temperature profile might be visibly non-linear. For even larger gradients, even though the system is in a steady state, we found that the linear response law eventually ceases to hold, but the local equilibrium is also violated.

We have classified the steady states in Table 1, which identify distinct dynamical regimes of the theory. It would be nice to develop more precise measures for these dynamical regimes, but it is clear that even a one component classical, lattice field theory does contain a means to understand the non-equilibrium physics of many-body systems. It would be interesting to extend these results to theory with phase transitions, multi-component theories which would also allow the analysis of Onsager reciprocity relations and so on. The additional degrees of freedom will provide additional measurable quantities, but we expect many of the qualitative features of this simple model to persist.

We acknowledge support through the grants from Keio University and DOE grant DE-FG02-91ER40608. We would like to thank Guy Moore and Larry Yaffe for enlightening discussions and the the Institute for Nuclear Physics at University of Washington for hospitality, where some of the work was conducted.


format: text

text:
# Logarithmic correction to the Bekenstein-Hawking entropy

## Abstract

The exact formula derived by us earlier for the entropy of a four dimensional non-rotating black hole within the quantum geometry formulation of the event horizon in terms of boundary states of a three dimensional Chern-Simons theory, is reexamined for large horizon areas. In addition to the semiclassical Bekenstein-Hawking contribution proportional to the area obtained earlier, we find a contribution proportional to the logarithm of the area together with subleading corrections that constitute a series in inverse powers of the area.

The derivation of the Bekenstein-Hawking (BH) area law for black hole entropy from the quantum geometry approach  (and also earlier from string theory  for some special cases), has lead to a resurgence of interest in the quantum aspects of black hole physics in recent times. However, the major activity has remained focussed on confirming the area law for large black holes, which, as is well-known, was obtained originally on the basis of arguments of a semiclassical nature. The question arises as to whether any essential feature of the bona fide quantum aspect of gravity, beyond the domain of the semiclassical approximation, has been captured in these assays. Indeed, as has been most eloquently demonstrated by Carlip , a derivation of the area law alone seems to be possible on the basis of some symmetry principle of the (semi)classical theory itself without requiring a detailed knowledge of the actual quantum states associated with a black hole. The result seems to hold for arbitrary number of spatial dimensions, so long as a particular set of isometries of the metric is respected. That quantum gravity has a description in terms of spin networks (or for that matter, in terms of string states in a fixed background) appears to be of little consequence in obtaining the area law, although these proposed underlying structures also lead to the same behaviour via alternative routes, in the semiclassical limit of arbitrarily large horizon area.

Although there is as yet no complete quantum theory of gravitation, one would in general expect key features uncovered so far to lead to modifications of the area law which could not have been anticipated through semiclassical reasoning. Thus, the question as to what is the dominant quantum correction due to these features of quantum gravity becomes one of paramount importance. Already in the string theory literature  examples of leading corrections to the area law, obtained by counting D-brane states describing special supersymmetric extremal black holes (interacting with massless vector supermultiplets) have appeared. This has received strong support recently from semiclassical calculations in $`N=2`$ supergravity  supplemented by ostensible stringy higher derivative corrections which are incorporated using Wald’s general formalism describing black hole entropy as Noether charge . However, the geometrical interpretation of these corrections remains unclear. Further, there are subtleties associated with direct application of Wald’s formalism which assumes a non-degenerate bifurcate Killing horizon, to the case of extremal black holes which have degenerate horizons. Moreover, the string results do not pertain to generic (i.e., non-extremal) black holes of Einstein’s general relativity, and are constrained by the unphysical requirement of unbroken spacetime supersymmetry.

In this paper, we consider the corrections to the semiclassical area law of generic four dimensional non-rotating black holes, due to key aspects of non-perturbative quantum gravity (or quantum geometry) formulated by Ashtekar and collaborators . In , appropriate boundary conditions are imposed on dynamical variables at the event horizon considered as an inner boundary. These boundary conditions require that the Einstein-Hilbert action be supplemented by boundary terms describing a three dimensional $`SU(2)`$ Chern-Simons theory living on a finite ‘patch’ of the horizon with a spherical boundary, punctured by links of the spin network bulk states describing the quantum spacetime geometry interpolating between asymptopia and the horizon. On this two dimensional boundary there exists an $`SU(2)`$ Wess Zumino model whose conformal blocks describe the Hilbert space of the Chern-Simons theory modelling the horizon. An exact formula for the number of these conformal blocks has been obtained by us  two years ago, for arbitrary level $`k`$ and number of punctures $`p`$. It has been shown that in the limit of large horizon area given by arbitrarily large $`k`$ and $`p`$, the logarithm of this number duly yields the area law. Here we go one step further, and calculate the dominant sub-leading contribution, as a function of the classical horizon area, or what is equivalent, as a function of the BH entropy itself.

On purely dimensional grounds, one would expect the entropy to have an expansion, for large classical horizon area, in inverse powers of area so that the BH term is the leading one,

$$S_{bh}=S_{BH}+\underset{n=0}{\overset{\mathrm{}}{}}C_nA_H^n$$
(1)
where, $`A_H`$ is the classical horizon area and $`C_n`$ are coefficients which are independent of the horizon area but dependent on the Planck length (Newton constant). Here the Barbero-Immirzi parameter  has been ‘fitted’ to the value which fixes the normalization of the BH term to the standard one. However, in principle, one could expect an additional term proportional to $`lnA_H`$ as the leading quantum correction to the semiclassical $`S_{BH}`$. Such a term is expected on general grounds pertaining to breakdown of naïve dimensional analysis due to quantum fluctuations, as is common in quantum field theories in flat spacetime and also in quantum theories of critical phenomena. We show, in what follows, that such a logarithmic correction to the semiclassical area law does indeed arise from the formula derived earlier  and derive its coefficient.

We first briefly recapitulate the derivation  of the general formula for the number of conformal blocks of the $`SU(2)_k`$ Wess Zumino model on a punctured 2-sphere appropriate to the black hole situation. This number can be computed in terms of the so-called fusion matrices $`N_{ij}^r`$

$$N_𝒫=\underset{\{r_i\}}{}N_{j_1j_2}^{r_1}N_{r_1j_3}^{r_2}N_{r_2j_4}^{r_3}\mathrm{}\mathrm{}N_{r_{p2}j_{p1}}^{j_p}$$
(2)
Diagrammatically, this can be represented as shown in fig. 1 below.

Here, each matrix element $`N_{ij}^r`$ is $`1or0`$, depending on whether the primary field $`[\varphi _r]`$ is allowed or not in the conformal field theory fusion algebra for the primary fields $`[\varphi _i]`$ and $`[\varphi _j]`$ ($`i,j,r=0,1/2,1,\mathrm{}.k/2`$):

$$[\varphi _i][\varphi _j]=\underset{r}{}N_{ij}^r[\varphi _r].$$
(3)
Eq. (2) gives the number of conformal blocks with spins $`j_1,j_2,\mathrm{},j_p`$ on $`p`$ external lines and spins $`r_1,r_2,\mathrm{},r_{p2}`$ on the internal lines.

We then use the Verlinde formula  to obtain

$$N_{ij}^r=\underset{s}{}\frac{S_{is}S_{js}S_s^r}{S_{0s}},$$
(4)
where, the unitary matrix $`S_{ij}`$ diagonalizes the fusion matrix. Upon using the unitarity of the $`S`$-matrix, the algebra (2) reduces to

$$N_𝒫=\underset{r=0}{\overset{k/2}{}}\frac{S_{j_1r}S_{j_2r}\mathrm{}S_{j_pr}}{(S_{0r})^{p2}}.$$
(5)
Now, the matrix elements of $`S_{ij}`$ are known for the case under consideration ($`SU(2)_k`$ Wess-Zumino model); they are given by

$$S_{ij}=\sqrt{\frac{2}{k+2}}sin\left(\frac{(2i+1)(2j+1)\pi }{k+2}\right),$$
(6)
where, $`i,j`$ are the spin labels, $`i,j=0,1/2,1,\mathrm{}.k/2`$. Using this $`S`$-matrix, the number of conformal blocks for the set of punctures $`𝒫`$ is given by

$$N_𝒫=\frac{2}{k+2}\underset{r=0}{\overset{k/2}{}}\frac{\underset{l=1}{\overset{p}{}}sin\left(\frac{(2j_l+1)(2r+1)\pi }{k+2}\right)}{\left[sin\left(\frac{(2r+1)\pi }{k+2}\right)\right]^{p2}}.$$
(7)
Eq. (7) thus gives the dimensionality of the $`SU(2)`$ Chern-Simons states corresponding to a three-fold bounded by a two-sphere punctured at $`p`$ points. The black hole microstates are counted by summing $`N_𝒫`$ over all sets of punctures $`𝒫`$, $`N_{bh}=_{\{𝒫\}}N_𝒫`$. Then, the entropy of the black hole is given by $`S_{bh}=\mathrm{log}N_{bh}`$.

We are however interested only in the leading correction to the semiclassical entropy which ensues in the limit of arbitrarily large $`A_H`$. To this end, recall that the eigenvalues of the area operator  are given by

$$A_H=8\pi \beta l_{Pl}^2\underset{l=1}{\overset{p}{}}[j_l(j_l+1)]^{\frac{1}{2}},$$
(8)
where, $`l_{Pl}`$ is the Planck length, $`j_l`$ is the spin on the $`l`$th puncture on the 2-sphere and $`\beta `$ is the Barbero-Immirzi parameter . Clearly, the large area limit corresponds to the limits $`k\mathrm{},p\mathrm{}`$. Now, from eq. (8), it follows that the number of punctures $`p`$ is largest for a given $`A_H`$ provided all spins $`j_l=\frac{1}{2}`$. Thus, for a fixed classical horizon area, we obtain the largest number of punctures $`p_0`$ as

$$p_0=\frac{A_H}{4l_{Pl}^2}\frac{\beta _0}{\beta },$$
(9)
where, $`\beta _0=1/\pi \sqrt{3}`$. In this approximation, the set of punctures $`𝒫_0`$ with all spins equal to one-half dominates over all other sets, so that the black hole entropy is simply given by

$$S_{bh}=lnN_{𝒫_0},$$
(10)
with $`N_{𝒫_0}`$ being given by eq. (7) with $`j_l=1/2`$.

Observe that $`N_{𝒫_0}`$ can now be written as

$$N_{𝒫_0}=\frac{2^{p_0+2}}{k+2}\left[F(k,p_o)F(k,p_0+2)\right]$$
(11)
where,

$$F(k,p)=\underset{\nu =1}{\overset{[\frac{1}{2}(k+1)]}{}}cos^p\left(\frac{\nu \pi }{k+2}\right).$$
(12)
The sum over $`\nu `$ in eq. (12) can be approximated by an integral in the limit $`k\mathrm{},p_0\mathrm{}`$, with appropriate care being taken to restrict the domain of integration; one obtains

$$F(k,p_0)\left(\frac{k+2}{\pi }\right)_0^{\pi /2}𝑑xcos^{p_0}x,$$
(13)
so that,

$$N_{𝒫_0}\frac{2^{p_0+2}}{\pi (p_0+2)}B(\frac{p_0+1}{2},\frac{1}{2}),$$
(14)
where, $`B(x,y)`$ is the standard $`B`$-function . Using well-known properties of this function, it is straightforward to show that

$`lnN_{𝒫_0}`$ $`=`$ $`p_0ln2{\displaystyle \frac{3}{2}}lnp_0ln(2\pi )`$ (15)
$``$ $`{\displaystyle \frac{5}{2}}p_0^1+O(p_0^2).`$ (16)
Substituting for $`p_0`$ as a function of $`A_H`$ from eq. (9) and setting the Barbero-Immirzi parameter $`\beta `$ to the ‘universal’ value $`\beta _0ln2`$ , one obtains our main result

$$S_{bh}=S_{BH}\frac{3}{2}ln\left(\frac{S_{BH}}{ln2}\right)+const.+\mathrm{},$$
(17)
where, $`S_{BH}=A_H/4l_{Pl}^2`$, and the ellipses denote corrections in inverse powers of $`A_H`$ or $`S_{BH}`$.

Admittedly, the above calculation is restricted to the leading correction to the semiclassical approximation. It has been done for a fixed large $`A_H`$ by taking the spins on all the punctures to be 1/2 so that we have the largest number of punctures. But it is not difficult to argue that the coefficient of the $`lnA_H`$ term is robust in that inclusion of spin values higher than 1/2 do not affect it, although the constant term and the coefficients of sub-leading corrections with powers of $`O(A_H^1)`$ might get affected. The same appears to be true for values of the level $`k`$ away from the asymptotic value which we have assigned it above: the coefficient of the $`lnA_H`$ is once again unaffected. Thus, the leading logarithmic correction with coefficient -3/2 that we have discerned for the black hole entropy is in this sense universal. Moreover, although we have set $`\beta =\beta _0ln2`$ in the above formulae, the coefficient of the $`lnA_H`$ term is independent of $`\beta `$, a feature not shared by the semiclassical area law.

It is therefore clear that the leading correction (and maybe also the subleading ones) to the BH entropy is negative. One way to understand this could be the information-theoretic approach of Bekenstein : black hole entropy represents lack of information about quantum states which arise in the various ways of gravitational collapse that lead to formation of black holes with the same mass, charge and angular momentum. Thus, the BH entropy is the ‘maximal’ entropy that a black hole can have; incorporation of leading quantum effects reduces the entropy. The logarithmic nature of the leading correction points to a possible existence of what might be called a ‘non-perturbative fixed point’. That this happens in the physical world of four dimensions is perhaps not without interest.

Recently, the zeroth and first law of black hole mechanics have been derived for situations with radiation present in the vicinity of the horizon, using the notion of the isolated horizon . Our conclusions above for the case of non-rotating black holes hold for such generalizations  as well. Note however that while, the foregoing analysis involves $`SU(2)_k`$ Chern Simons theory, for large $`k`$ this reduces to a specific $`U(1)`$ theory presumably related to the ‘gauge fixed’ classical theory discussed in . The charge spectrum of this $`U(1)`$ theory is discrete and bounded from above by $`k`$. The $`SU(2)`$ origin of the theory thus provides a natural ‘regularization’ for calculation of the number of conformal blocks.

Note Added: After the first version of this paper appeared in the Archives, it has been brought to our attention that corrections to the area law in the form of logarithm of horizon area have been obtained earlier  for extremal Reissner-Nordstrom and dilatonic black holes. These corrections are due to quantum scalar fields propagating in fixed classical backgrounds appropriate to these black holes. The coefficient of the $`lnA_H`$ term that appears in ref.  is different from ours. This is only expected, since in contrast to ref. , our corrections originate from non-perturbative quantum fluctuations of spacetime geometry (for generic non-rotating black holes), in the absence of matter fields. Thus, this correction is finite and independent of any arbitrary ‘renormalization scale’ associated with divergences due to quantum matter fluctuations in a fixed classical background.

We thank Prof. A. Ashtekar for many illuminating discussions and Prof. R. Mann for useful correspondence.