Nanothermodynamics: Stable Thermal Equilibrium and Nanoscale Fluctuations

1. Introduction

Thermodynamics and statistical mechanics provide two complementary approaches for understanding thermal processes in nature. Thermodynamics contains the fundamental physical laws governing what can or cannot occur, while statistical mechanics provides microscopic insight and quantitative tools for interpreting thermal behavior. Most physicists prefer the detailed pictures and predictive power of statistical mechanics. Indeed, Feynman has called Boltzmann’s factor a “fundamental physical law” that is the “the summit of statistical mechanics” [1]. But then he lists various assumptions needed to justify this “law.” Basically, for Boltzmann’s factor to fully describe any system, the system must have essentially instantaneous coupling to an effectively infinite and homogeneous heat bath. However, these assumptions are often unmet, especially for small systems that have finite-size thermal effects, including independently fluctuating small subsystems inside bulk samples. Therefore, a more basic approach is to first ensure that the system obeys the laws of thermodynamics before applying the tools of statistical mechanics. In practice, conservation of energy in thermal physics is usually found from the fundamental equation of thermodynamics, which combines the first and second laws.

Three people have made pioneering contributions to the fundamental equation of thermodynamics: Rudolf Clausius in 1850, Willard Gibbs in 1876, and Terrell Hill in 1962 (Figure 1) [2,3,4]. The first two names and their contributions are well known to essentially all scientists, whereas Hill’s contributions are relatively unknown. To put this into perspective, recall that it is legendary how Gibbs’ work was underappreciated for 15-20 years, until it was translated into German and used with great success by early physical chemists. It has now been 60-65 years since Hill introduced a mathematically similar (and comparably important!) contribution to thermodynamics, and yet it is still relatively unknown. One explanation comes from Joel Keiser: “It may be, as with much of Terrell’s work, that it was simply ahead of its time and that in future years much will be made of it,” but this statement was made about 40 years ago [5]. Here I intend to clarify the importance of Hill’s work by focusing on some key applications where statistical mechanics alone cannot explain the thermal behavior, so that one must go beyond Boltzmann’s factor to find the stable thermal equilibrium and resulting nanoscale fluctuations.

The outline of this paper is as follows. 2) Background reviews key concepts that aid in appreciating the power and utility of nanothermodynamics. 3) Measurements of Thermodynamic Heterogeneity presents various experimental results that establish the need to use nanothermodynamics for most types of materials, including liquids, glasses, spin glasses, polymers, and crystals. 4) Stable Thermal Equilibrium of Standard Models outlines the stable solution to Ising’s original model for finite chains of binary spins, and for the semiclassical ideal gas that gives a novel solution to Gibbs’ paradox. 5) Excess Energy Fluctuations in MD Simulations describes how and why a standard fluctuation relation must be modified to connect classical mechanics to statistical mechanics. 6) Nanoscale Fluctuations and the Second Law of Thermodynamics utilizes a simple model to test theoretical predictions for fluctuations of small systems, with the result that ideas based on thermodynamics yield better agreement with simulations than those based solely on statistical mechanics. 7) Summary contains the main conclusions from this review.

2. Background

Clausius first codified the two basic laws of thermodynamics: 1) that energy is always conserved and 2) that some of this energy is degraded by entropy. Specifically, in 1850, he recognized that the single internal energy of a system can be changed by two distinct processes involving either heat or work. Then he recognized that entropy limits the ability of internal energy to do external work. Combining these two laws yields Figure 1a, an early version of the fundamental equation of thermodynamics, written in modern form. The left side of the equation contains the change in total internal energy, $d{E}_t$, while the right side of Figure 1a contains two pairs of conjugate variables from two distinct processes that can change $E_t$. The heat term comes from temperature times the change in total entropy ($Td{S}_t$), while the work term comes from pressure times the change in total volume ($pd{V}_t$) with the minus sign from work done on the environment.

Gibbs’ greatest contributions to thermodynamics began in 1876 when he introduced a new pair of conjugate variables, the chemical potential times the change in total number of particles, $\mu d{N}_t$. The resulting fundamental equation for a simple system, such as a single-component ideal gas, is shown in Figure 1b. Gibbs also considered multiple types of particles, allowing the study of thermal equilibrium between different substances and between different phases of the same substance. This is the usual form of the fundamental equation found in textbooks. In general, the terms in Figure 1b apply only to the total quantities in effectively infinite systems, with an intensive variable ($T$, $p$, or $\mu$) times a change in an extensive variable ($d{S}_t$, $d{V}_t$, or $d{N}_t$). Because the equation in Figure 1b cannot systematically accommodate finite-size effects, non-extensive contributions to conservation of energy are missing from textbook treatments of thermodynamics and are similarly absent in standard statistical mechanics.

In 1962, Hill introduced the theory of small system thermodynamics [2]. As with Gibbs’ contribution, the foundation of Hill’s work involves a new pair of conjugate variables, the subdivision potential, $\mathcal{E}$, and the change in the number of subdivisions, $d\eta$. The resulting fundamental equation of thermodynamics is given in Figure 1c. In general, the number of small subsystems that are in a large system is given by $\eta +1$. Thus, $\eta$ varies inversely proportional to the average number of particles in each subsystem, $\eta \approx N_t/\overline{n}$, in contrast to the total number of particles in the system, $N_t$, that is purely extensive. Therefore, while Gibbs’ term dominates changes in large systems, ${dN}_t$, Hill’s term dominates changes in small systems, $\left|d\eta \right|~|d\overline{n}/{\overline{n}}^2|$ with $\overline{n}\ll N_t$. Indeed, Hill’s subdivision potential can be understood by comparison to Gibbs’ chemical potential. Gibbs’ $\mu$ is the change in (free) energy to take one particle from a bath of particles into the system, whereas Hill’s $\mathcal{E}$ is related to the change in (free) energy to take a cluster of interacting particles from a bath of clusters into the system. In general, a cluster of $n$ interacting particles does not have the same energy as $n$ isolated particles. Specifically surface terms, length-scale effects, and fluctuations all contribute to the energy of a finite cluster of interacting particles, so that Hill’s fundamental equation (Figure 1c) is required to systematically include these contributions for conservation of energy across all size scales. I will focus on small systems, especially small subsystems inside larger systems, where these contributions are difficult, if not impossible to include systematically in Boltzmann’s factor. In such examples, Hill’s small system thermodynamics is essential for understanding the stable thermal equilibrium and nanoscale fluctuations.

Hill developed his theory by considering large ensembles of small systems. Indeed, the foundation of his theory is that a large ensemble of small systems must behave like a large system, obeying the laws of standard thermodynamics. Starting in 1999, I adapted Hill’s theory by imagining the mathematically similar but conceptually inverse picture of a large system that subdivides into an ensemble of small subsystems [6]. In 2000, I first published the term “nanothermodynamics” for this inverse picture [7]. In 2001 Hill published an explanation that he used the term “nanothermodynamics” as “a shortened and more fashionable version of ‘small system thermodynamics’” [8]. Hill never fully accepted the inverse picture that I was working on, presumably because he did not recall his own mechanism favoring an ensemble of independent subsystems inside a large system. However, I stubbornly believed in the ubiquity of his theory, and continued to develop the inverse picture based on an abundance of experimental evidence showing thermodynamic heterogeneity inside bulk samples. One such experiment involves time-resolved specific heat (TRSH), especially measurements made by Meissner and collaborators [9]. Another experiment, that is effectively the inverse of TRSH, is non-resonant spectral hole burning (NSHB). NSHB was first proposed in 1994 and developed starting in 1996 [10,11,12]. These experiments provide some of the earliest and most direct evidence for independently fluctuating subsystems inside bulk samples. Now this thermodynamic heterogeneity is well established in essentially all types of materials, including liquids, glasses, spin glasses, polymers, and crystals. Some examples are given in Section 3, below.

Hill’s small system thermodynamics provides the theoretical foundation for the inverse process of nanothermodynamics, where a large system subdivides into a stable distribution of independent subsystems. This stable distribution is defined by setting Hill’s subdivision potential to zero, $\mathcal{E}=0$, as described in section 10-3 of Hill’s 1964 book [3]. Indeed, quoting from this section: “Although $\mathcal{E}$ is negligible for a macroscopic system (. . .), it is not equal to zero in the strict sense that we are using $\mathcal{E}=0$ above as an equilibrium condition. The macroscopic state is therefore not to be confused with the equilibrium state.” In other words, macroscopic systems that are constrained to be homogeneous so that they can be treated using standard thermodynamics are in an unstable equilibrium, having $\mathcal{E}\approx 0$ but not $\mathcal{E}=0$. Similarly, small systems with $\mathcal{E}\ne 0$ must also be constrained to prevent them from reaching their stable thermal equilibrium. With few exceptions, since 1964 no other researchers have used this $\mathcal{E}=0$ stability condition for the treatment of small systems or subsystems, not even Hill himself. Therefore, to emphasize the importance of the stable thermal equilibrium, and to reduce confusion [13] with the work of others, I will continue to restrict use of the term “nanothermodynamics” to the inverse problem that requires $\mathcal{E}=0$ for the stable equilibrium of independent subsystems inside larger systems, with “small system thermodynamics” used for the theory of small systems that are constrained to have $\mathcal{E}\ne 0$. It is difficult for me (as it was for Hill) to imagine how the inverse process of subdividing large systems could have $\mathcal{E}\ne 0$ for their internal subsystems.

Since the inception of nanothermodynamics for understanding how large systems subdivide into stable distributions of independent subsystems [6,7], a key goal has been to find applications that justify this extension of Hill’s theory to bulk materials. At first the applications were postulated, based primarily on experimental evidence for thermodynamic heterogeneity. Some examples include a mesoscopic mean-field theory for supercooled liquids and the glass transition [6,14], a mean-field cluster model for non-classical critical scaling [7], and maximum entropy as a mechanism for $1/f$-like noise [15]. However, over the past decade, several theoretical results and computer simulations have clarified the connections. Therefore, an emphasis in this review will be the need to go beyond standard thermodynamics and statistical mechanics to include finite-size thermal effects in experiments, theory and simulations.

It is useful to preview some key results from the rest of this paper. Section 3 focuses on experiments that show thermodynamic heterogeneity, providing direct evidence for the necessity of nanothermodynamics. Section 4 describes how nanothermodynamics is needed to obtain the stable solutions of two standard models. First is the stable equilibrium of Ising’s original (1925) model for finite chains of interacting binary spins. Second is a novel solution to Gibbs’ paradox, removing the residual non-extensive contributions to entropy that cause inconsistencies between statistical mechanics and thermodynamics. It is interesting to note that this solution to Gibbs’ paradox is contained in the equations Hill published in his 1964 book [3], but at that time he seems not to have recognized the physical significance. Section 5 and Section 6 present results from computer simulations used to study nanoscale thermal fluctuations. Specifically, Section 5 describes molecular dynamics (MD) simulations that address Boltzmann’s quest to find a connection between classical mechanics and statistical mechanics. The answer is that MD simulations give good agreement with statistical theories for thermal averages, but not for fluctuations. Instead, the MD simulations are quantitatively characterized by a theory that combines statistical mechanics with adiabatic modulations that are effectively decoupled from the heat bath, consistent with measurements of thermodynamic heterogeneity. Section 6 utilizes a one-dimensional (1-D) Creutz-like model [16,17] for tests of the 2nd law of thermodynamics. An appealing feature of this model is that it has an explicit heat bath that allows calculation of the exact total entropy ($S_t$) of the system and its bath. One result is that an intrinsically irreversible microscopic step is needed to maximize $S_t$, providing a clear counterexample to the claim that the 2nd law is an emergent property of reversible dynamics. Another result that is a consequence of strict adherence to the 2nd law is that fluctuations in $S_t$ deviate from standard theories of thermal fluctuations based on statistical mechanics, including detailed fluctuation theorems [18,19,20,21]. Instead, the fluctuations are consistent with Einstein’s reversal of Boltzmann’s relation and the 2nd law. Speculations are made about the relevance of such models to real systems, and what might cause the intrinsically irreversible dynamics needed for thermal-equilibrium behavior. Finally, Section 7 gives a summary and some conclusions.

3. Measurements of Thermodynamic Heterogeneity

3.1. Relevance of Heterogeneity to Nanothermodynamics

If nanothermodynamics is to apply to large systems, there must be thermodynamic heterogeneity from an ensemble of independent subsystems inside the system. Similarly, any system that shows thermodynamic heterogeneity cannot be described by standard thermodynamics that applies only to systems having a heat bath that is homogeneous and effectively infinite. In other words, nanothermodynamics is irrelevant for any system that is large and homogeneous, but it is essential for systems of any size that exhibit thermodynamic heterogeneity. The purpose of this section is to outline some of the experimental evidence for thermodynamic heterogeneity that is found in virtually all types of materials, including liquids, glasses, spin glasses, polymers, and crystals.

Dynamic heterogeneity, where separate responses come from independent subsystems, has been a theoretical concept at least since 1874 when Boltzmann introduced his superposition principle [22]. Thermodynamic heterogeneity is closely related. In fact, nanothermodynamics applies to any system having internal fluctuations that are statistically independent. Such statistical independence is needed so that neighboring subsystems have uncorrelated dynamics, yielding entropies that are additive for an ensemble of independent subsystems inside a bulk sample. Because local entropies are difficult to measure directly, distributions of effective local temperatures ($T_i$) usually provide the clearest evidence for thermodynamic heterogeneity. Therefore, we describe measurements that imply distinct values of $T_i$ for independent degrees of freedom inside individual samples.

3.2. Time-Resolved Specific Heat Measurements

Some of the earliest evidence for thermodynamic heterogeneity comes from time-resolved specific heat (TRSH) measurements at low temperatures, $T\ll 100$ K. The heterogeneity is deduced from time-dependent changes in effective temperature ($∆T_i$). Anomalies in $∆T_i$ reveal deviations from Debye’s theory for the specific heat, $C=C_D{T}^3$, where Debye’s theory yields values of $C_D$ from the elastic constants for harmonic modes. Specific heats that far exceed $C_D{T}^3$ were first found in amorphous materials [23,24] and attributed to tunneling defects [25,26]. However, ultra-pure single crystals also exhibit deviations from $C_D{T}^3$ [9]. Although the excess specific heat in crystals has been attributed to residual imperfections, the discovery of analogous excess energy fluctuations in MD simulations of crystals with no defects suggests that the excess specific heat may be intrinsic [27,28]. Evidence for excess specific heat is found in all materials that have been measured by TRSH at low enough temperatures. Symbols in Figure 2 show results from TRSH measurements on five different high-purity single crystals, given in the legend. The sketch in the insert helps to interpret the experiment. Each measurement is made by applying a heat pulse to one side of the sample and measuring the temperature as a function of time on the other side. Symbols in the insert show measurements from a quartz crystal at a base temperature of $T~0.2$ K. Note how $\mathsf{\Delta }{T}_i$ starts at zero, rises sharply until it reaches a peak at $t~0.01$ ms, then recovers towards zero in multiple steps. The peak at ~0.01 ms is due to the arrival of ballistic phonons from the heater to the thermometer. After the peak, as the added heat becomes dispersed throughout the heat bath, $\mathsf{\Delta }{T}_i$ stabilizes at a short-time plateau that is consistent with $C_D{T}^3$ for this crystal (black line). However, starting at $t~1$ ms, $\mathsf{\Delta }{T}_i$ relaxes again, reaching a long-time plateau (red dashed line) at $t>10$ ms that gives the equilibrium specific heat for the sample. Finally, at still longer times ($t>0.5$ s, not shown) $\mathsf{\Delta }{T}_i$ returns to zero as the excess heat flows out of the fine wires that connect the sample to the cryostat.

Note that specific heat is found from the inverse of $∆T_i$. Thus, the insert in Figure 2 shows that the equilibrium specific heat (red dashed line) is significantly higher than that predicted by Debye’s theory (black line). Furthermore, during the nearly exponential relaxation that occurs between 1 ms$<t<$ 10 ms the excess heat supplied to the sample is not flowing across the sample (which happens at $t<0.1$ ms), and it is not flowing out of the sample (which happens at $t>$ 100 ms), so it must be flowing into distinct degrees of freedom inside the sample. Specifically, the first plateau (that agrees with Debye’s theory) comes from the homogeneous heat bath of phonons, while the second plateau is due to some other degrees of freedom in the sample. From the zeroth law of thermodynamics, net heat flows only if there is a temperature difference. Thus, this experiment establishes that there are multiple temperatures inside a sample, which persist for times $t\gg 1$ ms. Furthermore, while these slow degrees of freedom are relaxing (or fluctuating), their exchange of energy with the heat bath is slow ($t\gg 1$ ms), thereby failing the fundamental requirement of fast thermalization needed for Boltzmann’s factor. This is direct evidence of thermodynamic heterogeneity. For such slow degrees of freedom, temperature is not a good thermodynamic variable, so that a full understanding requires going beyond Boltzmann’s factor to describe the fluctuations and relaxation. When measured by TRSH at low enough temperatures, similar slow relaxation is found in every sample studied, including dielectric crystals and amorphous materials.

Solid symbols in Figure 2 show the temperature dependence of the long-time (equilibrium) values of normalized specific heat from TRSH measurements of five different high-purity single crystals, given in the legend. The solid black line at $C_p/C_D{T}^3=1$ shows the $T$-dependence expected from Debye’s theory. The dashed lines come from power-law fits to the low-$T$ data for SiO₂ (red circles) and Si (blue triangles), showing how the equilibrium specific heat diverges from Debye behavior with an inverse power of $T$. This excess specific heat, that has been attributed to residual defects that are unavoidable in even ultra-pure crystals, could instead be due to excitations that are localized by anharmonic interactions. The anharmonicity decouples the interactions from the heat bath of phonon modes, yielding excitations that conserve local energy instead of having a single well-defined $T$. Qualitatively similar deviations at low temperatures are found in energy fluctuations of MD simulations, as described in Section 5, below.

TRSH measurements establish that there are at least two distinct temperatures in most (if not all) samples during their slow response that yields the excess specific heat at low temperatures. At higher temperatures, the excess specific heat becomes much smaller than standard Debye behavior. Specifically, extrapolating the dashed line for SiO₂ in Figure 2 to 1 K, the excess contribution is about 30% of the Debye value, and by room temperature the excess will be in the ppm range. However, at all temperatures, the excess specific heat gives 100% of the slow response that governs thermal equilibrium. Non-resonant spectral hole burning (NSHB) provides additional information about the slow degrees of freedom that cause the thermodynamic heterogeneity that justifies the need for nanothermodynamics.

3.3. Nonresonant Spectral Hole Burning in a Spin Glass

The sketch in the upper right of Figure 3 depicts the so called “box model” that is used to characterize the behavior of NSHB [10,11,12]. Thus, NSHB is essentially the inverse of TRSH (Figure 2). Specifically, in NSHB, a large amplitude low-frequency oscillation in an external field adds excess energy to slow degrees of freedom that are selected by the frequency ($f$) of the field. This oscillation can be in the electric field for dielectrics [10], in the magnetic field for magnetic materials [11], or in the stress for mechanical response [29]. The resulting excess energy changes the time scale of the response, which can be interpreted as a change in local temperature, $∆T_i$, due to “heating” of selected degrees of freedom in the sample. Usually, the slow degrees of freedom couple weakly to the heat bath (heat resistors in the sketch), yielding the distribution of slow relaxation times (${\tau }_i$) for the $∆T_i$ to return to zero.

The main part of Figure 3 shows NSHB measurements on a 5% Au:Fe spin glass [11,12]. This sample has a sharp peak in the field-cooled magnetization (not shown) that defines a spin-glass temperature of $T_g=21.6$ K. NSHB measurements are made utilizing a single oscillation “pump” in the magnetic field having a large-amplitude ($H_0$) at low $f$, followed by a small step in the field ($h=8$ Oe) to measure the linear response of the sample after being nonlinearly modified by the pump. The uppermost set of data in Figure 3 are from NSBH measurements at $T=$ 18.8 K. Symbols show the linear response with no pump (black squares), and after a pump oscillation of $H_0=$ 96 Oe at $f=$ 0.03 Hz (red circles), 1.0 Hz (green up triangles), and 30 Hz (blue down triangles). Even in the raw data the influence of the pump frequency is clear: $f=0.03$ Hz shifts the slow degrees of freedom at long times, $f=30$ Hz shifts the “fast” degrees of freedom at short times, while $f=1.0$ Hz shifts in the middle. Indeed, there is virtually no overlap in the degrees of freedom shifted by the lowest- and highest-frequency pumps, clearly and directly establishing the heterogeneous nature of the non-exponential response. In other words, there is negligible correlation between the slow and fast degrees of freedom. Solid curves show predictions of the box model with the distribution of relaxation times defined by the linear response (black), and hole burning behavior (red, green, blue) found from the known pump frequency and amplitude. The model gives a specific heat per spin of $\mathsf{\Delta }C/k=5x{10}^{-7}$ as the only free parameter.

Such a tiny specific heat per spin implies that relatively small magnetic fields can greatly increase the spin temperature above the bath temperature, consistent with the behavior shown in Figure 3 for NSHB at $T=$ 19.6 and at 20.6 K. First, at $T=$ 19.6 K, note that pump oscillations of $H_0\le 96$ Oe still produce spectral holes due to simply shifting the response times with no net loss in amplitude over the measurement time window. However, for $H_0=115$ Oe, a spectral “step” emerges where significant amplitude is lost from the time window. This spectral step indicates that many degrees of freedom are “heated” by at least 2 K to above $T_g$, where they respond on paramagnetic time scales ($t\ll 10 \mu s$), while the rest of the sample remains a spin glass. This interpretation is confirmed by measurements at $T=$ 20.6 K, where $H_0=77$ Oe is sufficient to heat about 30% of the sample into the paramagnetic phase. Measurements, modeling, and interpretations from Figure 3 indicate that spin temperatures can differ greatly from the bath temperature during the response. Such significant changes in spin temperature are well known in magnetism [30], including spin temperatures in magnetic resonance that can be negative. Thus, local temperatures that differ from the bath, with multiple values that depend on position in the sample are useful for interpreting the behavior, providing evidence for the thermodynamic heterogeneity that is needed for nanothermodynamics.

3.4. Nonresonant Spectral Hole Burning in Glass-Forming Liquids

Figure 4b shows dielectric NSHB measurements on a glass-forming liquid (propylene carbonate, PC) at $T=187.5$ K, slightly above the glass temperature, $T_g\approx 186$ K [10,12]. Again, there is negligible overlap in the degrees of freedom modified by the highest and lowest frequencies, directly establishing the heterogeneous nature of the slow response. For PC, the specific heat per molecule from the slow response is large enough to be measured directly, $\mathsf{\Delta }C/k\approx 7.85$ [31]. For theoretical analysis, knowing $∆C$ has the advantage that the box model (Figure 3) can be used with no adjustable parameters. Results utilizing a simplified box model (without heat resistors, ${\tau }_i\to \infty$) are shown by the smooth curves that mimic the measurements. Better agreement with the data, especially after the pump oscillation at lowest frequency, can be obtained (also with no free parameters) by using heat resistors with a distribution of recovery times (${\tau }_i$) that match the distribution of relaxation times in the equilibrium response [32]. Each ${\tau }_i$ governs how the associated effective local temperature returns to the bath temperature, $∆T_i\to 0$. This model quantifies what must happen for systems that obey Joule heating, Ohm’s law, and Fourier’s law. The main experimental challenge from the large $∆C$ is that very high fields are needed to yield measurable spectral holes, so that it is difficult to increase any effective local temperature to more than ~0.04% above the bath temperature. Still, the axes in Figure 4b show how the shift in response times after the pump oscillation (left axis) can be converted to a change in effective local temperatures (right axis). Again, it is useful to characterize the slow response using multiple internal temperatures, yielding the thermodynamic heterogeneity needed for nanothermodynamics.

Although NSHB and TRSH establish the ubiquity of thermodynamic heterogeneity in a wide variety of substances, other techniques are needed to determine the length scale of the heterogeneity. Perhaps the least invasive and most direct technique involves multi-dimensional nuclear magnetic resonance (NMR). Figure 4a is a picture that emerges from these measurements [33]. Analysis of the measurements yields typical length scales of 1-3 nm [34] (~10-100 molecules) for the thermal and dynamic heterogeneity of slow degrees of freedom in glass-forming liquids. Such nanoscale lengths and distributions of independent fluctuations inside bulk samples emphasize the need for nanothermodynamics.

4. Stable Thermal Equilibrium of Two Standard Models

4.1. Relevance to Nanothermodynamics

Perhaps the most important application of nanothermodynamics is to treat finite-size thermal effects in finding the stable equilibrium of systems, especially when a large system subdivides into subsystems. Subdividing into internal subsystems is usually favored due to an increase in entropy that lowers the net free energy. For sufficiently simple models, Hill’s subdivision potential can be calculated directly, allowing Hill’s fundamental equation to be used to determine the stable thermal equilibrium. Two such simple systems are Ising’s original model for magnetic spins, and the semiclassical ideal gas.

4.2. Stable Thermal Equilibrium of Ising’s Model for Finite Chains of Interacting Spins

Detailed histories of Ising and his model are given in [35,36]. Briefly, in 1925 Ernst Ising published his solution to the Lenz-Ising model for finite chains of interacting binary degrees of freedom (“spins”) [37]. His solution involved summing over all possible chain lengths, then forcing the chain to be infinitely long and homogeneous. Since 1925, most (but not all [38,39,40]) discussions of the Ising model have also assumed that the spin system is infinite and homogeneous, but this constraint prohibits the system from finding its true equilibrium if finite-size thermal effects are included. Here, we focus on the recent discovery of the stable thermal equilibrium of Ising’s original model [13,40]. Although Ising could not have found this equilibrium 40 years before Hill’s work, it is now a useful exercise that emphasizes the necessity of nanothermodynamics.

Most modern treatments of the Ising model start by assuming that the system is infinite and homogenous, as represented by the upper sketch in Figure 5. This limit is needed for standard thermodynamics. However, in Ising’s original work he treated finite chains, shown by the right-side sketch in Figure 5, which require nanothermodynamics. This right-side sketch shows several important features. Note that spin interactions are low energy between aligned spins (•), high energy between anti-aligned spins (X), or no energy when there is a break in the interaction (o). Breaks in the interaction subdivide the initially infinite and homogeneous system into finite chains, as originally treated by Ising. Also shown with the spins are two types of heat baths. The upper bath is idealized, assumed to be effectively infinite and homogeneous with essentially instantaneous coupling to every spin, thereby directly imposing a uniform temperature on the spins. The lower bath is explicit, comprised of a set of Einstein oscillators (Creutz “demons”), which accommodates heterogeneity in the spin system and facilitates testing the validity of Boltzmann statistics. For example, local temperatures can be extracted from the occupation probabilities in the equally spaced energy levels. Furthermore, this explicit heat bath provides information about the spins as they thermalize and allows the spins and bath to form a closed system for testing the 2nd law of thermodynamics (see Section 6, below).

The Lenz-Ising Hamiltonian for a 1-D chain of $n$ links ($n+1$ spins) in zero field is

$$U(n,n_x)=-J{\sum }_{i=0}^n{\sigma }_i{\sigma }_{i+1}.$$

(1)

Here, $J$ is the interaction energy between nearest-neighbor spins, with ${\sigma }_i=+1$ if the $i^{th}$ spin is up and ${\sigma }_i=-1$ if it is down. (Note that volume, $V$, is absent from Equation (1) because the energy is assumed to be independent of density.) Referring to the sketches of spins in Figure 5, this energy can be expressed in terms of the number of high-energy interactions ($n_x$),

$$U\left(n,n_x\right)=-J\left(n-2n_x\right).$$

(2)

Assume that the spin chain is ideally coupled to a heat bath having temperature $T$. Using the binomial coefficient for the multiplicity of these energies, the partition function for an Ising chain having a fixed number of spins ($n+1$) comes from summing over all energies,

$$Q(n,T)={\sum }_{n_x=0}^n{\displaystyle \frac{2n!}{\left(n-n_x\right)!n_x!}}{e}^{-U(n,n_x)/kT}=2^{n+1}{\left[\cosh \left(J/kT\right)\right]}^n.$$

(3)

Equation (3) is the partition function in the canonical ensemble. It is consistent with the zero-field limit of the size-dependent partition function that Ising gave in [37] (middle of pg. 256).

Now assume that the spin chain is ideally coupled to a bath of spins having chemical potential $\mu$ with $e^{\mu /kT}<1/\left[2\cosh (J/kT)\right]$ to ensure convergence. Because there are $n+1$ spins in the chain the chemical potential term is raised to the power of $n+1$. Summing the resulting series gives

$$\mathsf{{\rm Y}}\left(\mu ,T\right)={\sum }_{n=0}^{\infty }Q(n,T)e^{\mu (n+1)/kT}=2e^{\mu /kT}/\left[1-2e^{\mu /kT}\cosh (J/kT)\right].$$

(4)

Equation (4) is the partition function in the generalized ensemble. Although Equation (4) is formally equivalent to an equation in [37] (equation 7), Ising used the dummy variable $x\leftrightarrow e^{\mu /kT}$ as a tool to study the mathematics of 1-D chains. Ising’s mathematical treatment made sense in 1925 because Equation (4) (with no extensive environmental variables to fix the size of the system) is ill-defined in standard statistical mechanics. Small-system thermodynamics is needed to justify the physics of Equation (4). Specifically, Hill’s subdivision potential is

$$\mathcal{E}=-kT\ln \left(\mathsf{{\rm Y}}\right),$$

(5)

and the average number of spins in a chain is

$$\overline{n}\left(\mu ,T\right)+1={\displaystyle \frac{\partial \ln \left(\mathsf{{\rm Y}}\right)}{\partial (\mu /kT)}}={\left[1-2e^{\mu /kT}\cosh \left(J/kT\right)\right]}^{-1}.$$

(6)

Because Eqs. (4-6) are in the generalized ensemble, standard thermodynamics and statistical mechanics give no guidance on how to proceed. A key question is: what value of $\mu /kT$ is needed for Equation (6) to yield the equilibrium average chain length as a function of $T$? It is useful to contrast Equation (6) with the usual procedures for obtaining $\mu /kT$ and $\overline{n}(\mu ,V,T)$ in the grand-canonical ensemble. If there is a fixed density of particles, then $\mu /kT$ is chosen to give this $\overline{n}/V$. Familiar examples include the densities of charge carriers in conductors and bosons in superfluids. Alternatively, if the density of particles is not fixed, then $\mu /kT=0$ is used to minimize the free energy when the system itself controls the density of particles. Familiar examples include photons in black-body radiation and phonons in solids.

Using Equation (6) to write the chemical potential in terms of the average number of links in a chain, $e^{\mu /kT}=1/\left[2\left(1+1/\overline{n}\right)\cosh (J/kT)\right]$, the subdivision potential becomes

$$\mathcal{E}=-kT\ln \left(\mathsf{{\rm Y}}\right)=-kT\ln [\overline{n}/\cosh \left(J/kT\right)].$$

(7)

In standard statistical mechanics, $\mu =0$ is used if there are no external constraints on the density of particles. In nanothermodynamics, $\mathcal{E}=0$ [$\overline{n}=\cosh \left(J/kT\right)$] is used if there are no external constraints on subsystem size. Thus, the average number of spins per finite chain in the stable thermal equilibrium of Ising’s original model is given by

$$\overline{n}+1=1+\cosh (J/kT).$$

(8)

The free energy reduction as a function of average subsystem size, black curve in Figure 5, gives insight into understanding the thermal stability. This $∆F_t$ is obtained from Hill’s fundamental equation, Figure 1c, by integrating the subdivision potential over the number of subdivisions, $\mathsf{\Delta }{F}_t={\int }_0^{\overline{\eta }}\mathcal{E}d\eta$. For Figure 5, a change of variables is used to integrate over $n=N_t/(\eta +1)$. Starting with an effectively infinite system, $N_t\to \infty$, the change in free energy per particle to subdivide the system into subsystems of average size $\overline{n}$ is

$$∆F_t/\left(N_{t}kT\right)={\int }_0^{\overline{\eta }}\left(\mathcal{E}/kT\right)d\left(\eta /N_t\right)={\int }_{\infty }^{\overline{n}}\ln \left[n/\cosh \left(J/kT\right)\right]dn/n^2 \mathrm{br}-\mathrm{to}-\mathrm{break}=-\left\{\ln \left[\overline{n}/\cosh (J/kT)\right]+1\right\}/\overline{n}.$$

(9)

The black curve in Figure 5 shows Equation (9), using $J/kT=1$ for specific values. The minimum free energy occurs at $\partial \left(∆F_t\right)/\partial \overline{n}=0$, yielding $\mathcal{E}=0$ and $\overline{n}+1=1+\cosh (J/kT)=2.543\dots$. The resulting reduction in total free energy due to subdividing is significant at all nonzero temperatures. Thus, the large homogeneous solution is never in thermal equilibrium at any temperature, $T>0$. Indeed, the stable thermal equilibrium of Ising’s original model cannot be found from standard thermodynamics or statistical mechanics alone, instead insight from nanothermodynamics is required to obtain the minimum in free energy.

4.3. Equivalence Between Ising’s Original Model and a Model for Random Interactions

Equation (8) gives the average number of spins for the Ising chain in stable thermal equilibrium. It is obtained from the 1-D Ising Hamiltonian [Equation (1)] using the nanocanonical ensemble that is unique to nanothermodynamics. The same answer is obtained in the canonical ensemble from standard statistical mechanics by using a different Hamiltonian,

$$U\left(N,N_x,N_0\right)=-{\sum }_{i=1}^N{J}_i{\sigma }_i{\sigma }_{i+1}{J}_i=(J,0).$$

(10)

Here, $N$ is the total number of links between $N+1$ spins, with $N_x$ the number of high-energy links and $N_0$ the number of links having no interaction energy, $J=0$. The right-side sketch in Figure 5 shows $N_0=3$. Equation (10) differs from Equation (1) by having two values for the exchange interaction, $J$ or 0. If $J_i$ had a frozen distribution, Equation (10) would be analogous to a 1-D Edwards-Anderson spin glass model with a bimodal distribution [41]. Instead, we allow the interactions to change freely, $J\leftrightarrow 0$, yielding an equilibrium distribution of $J_i$. Using the number of each type of link, the energy of Equation (10) can be written as

$$U\left(N,N_x,N_0\right)=-J\left(N-2N_x-N_0\right).$$

(11)

Utilizing the trinomial coefficient for the multiplicity of each energy, the partition function for an Ising chain having a fixed number of spins ($N+1$) comes from summing over all energies,

$$Q(N,T)={\sum }_{N_0=0}^N{\sum }_{N_x=0}^{N-N_0}{\displaystyle \frac{2^{N_0+1}N!}{\left(N-N_x-N_0\right)!N_x!N_0!}}{e}^{-U(N,N_x,N_0)/kT}=2^{N+1}{\left[1+\cosh \left(J/kT\right)\right]}^N.$$

(12)

This is the partition function for a 1-D random-exchange Ising model in the canonical ensemble. It gives an average energy of

$$\overline{U}={\displaystyle \frac{\partial \ln Q}{\partial (-1/kT)}}=-NJ{\displaystyle \frac{\sinh \left(J/kT\right)}{1+\cosh \left(J/kT\right)}}.$$

(13)

This $\overline{U}$ differs from the usual expression for the average energy of the Ising model because many links between spins do not contribute to the energy. However, using the average number of breaks in the interaction, $\overline{N_0}=N/\left[1+\cosh \left(J/kT\right)\right]$ [13], the average energy per interacting link is identical to the standard Ising model,

$$\overline{U}/(N-\overline{N_0})=-J\tanh \left(J/kT\right).$$

(14)

Furthermore, using the probability of breaks in the interaction, $p=\overline{N_0}/N$, with percolation theory in the thermodynamic limit, the average number of interacting links between breaks in the interaction,

$$\overline{n}={\sum }_{n=0}^{\infty }n{\left(1-p\right)}^n{p}^2/{\sum }_{n=0}^{\infty }{\left(1-p\right)}^n{p}^2=(1-p)/p=\cosh \left(J/kT\right),$$

(15)

is equivalent to Equation (8). This emphasizes another key difference between standard statistical mechanics and nanothermodynamics. If finite-size thermal effects are included, the behavior depends on the choice of ensemble. Moreover, different Hamiltonians solved in their appropriate ensembles can yield identical results, but only if small subsystems are in the stable thermal equilibrium from the nanocanonical ensemble of nanothermodynamics.

4.4. Thermal Equilibrium of the Semi-Classical Ideal Gas

Nanothermodynamics provides a novel solution to Gibbs’ paradox that has several advantages over previous solutions [13,40]. First, the nanocanonical ensemble yields entropy that is exactly extensive for systems of all sizes, unlike previous solutions that are nonextensive for small systems. Specifically, entropy is precisely doubled when two systems having the same average number of subsystems are combined. Indeed, because previous solutions of Gibbs’ paradox neglect finite-size effects, nanothermodynamics is the only way to obtain extensive entropy for small systems. Second, as depicted in the sketch at the bottom of Figure 5, particles are indistinguishable only if they are within the same subsystems, allowing distant particles in large systems to be distinguishable by their locations when in distinct subsystems. Third, nanothermodynamics provides an explanation for the empirical accuracy of the Sackur-Tetrode equation. Specifically, this nanothermodynamic solution may be necessary to explain the absence of isotope effects in the measured entropies of real gases [42], thereby solving various questions about the entropy of mixing [43]. Fourth, as with Ising’s model for finite chains of spins, only the nanocanonical ensemble yields the minimum free energy of the semiclassical ideal gas. The remainder of this section outlines the derivation and advantages of the nanothermodynamic solution.

As with Ising’s original model, Hill’s subdivision potential can be calculated exactly for the semi-classical ideal gas in the generalized ensemble. In fact, Hill gives an answer in his 1964 book [3] (Equation 10-88) $\mathcal{E}=-kT\ln (1+\overline{N})$, where $\overline{N}$ is the average number of particles in a subsystem. However, this result assumes that the sum over particles in each subsystem starts at $N=0$, whereas Hill’s stability condition requires that each subsystem has at least one particle. Indeed, the sentence associated with the first footnote in [3] (section 10-3) states that “an empty ‘system’…is not counted as a system of the ensemble. This is the natural, realistic, and necessary point of view…in a solvent or gas.” Therefore, starting the summation at $N=1$ in Hill’s Equation 10-85 gives

$$\mathcal{E}=-kT\ln (\overline{N}-1).$$

(16)

This term lowers the total free energy when a large system subdivides into subsystems. Most previous studies of small-system thermodynamics assume $\mathcal{E}\ne 0$, which applies only if there are external constraints on the average subsystem size [44]. Instead, $\overline{N}=2$ is needed to yield $\mathcal{E}=0$ for the stable equilibrium of nanothermodynamics. Any value, $\overline{N}\ne 2$, fails to minimize the free energy, as indicated by the red curve in Figure 5. In fact, it is difficult to imagine how external constraints could control the distribution of subsystems inside bulk samples to yield $\mathcal{E}\ne 0$. An analogy from standard statistical mechanics would be if there were external constraints on a black body that controlled the distribution of photons by imposing $\mu /kT\ne 0$. Therefore, we focus on the stable equilibrium of nanothermodynamics.

The lower sketch in Figure 5 depicts how a semiclassical ideal gas can be subdivided into small subsystems. Ideal gases have no interactions, so that subsystems must be identified by another mechanism. We relate the dashed lines in the sketch with boundaries between subsystems of distinguishable particles. Specifically, all particles within a subsystem are indistinguishable from each other, but they are distinguishable from particles in other subsystems. Because distinguishable particles have higher entropy than indistinguishable particles, the total entropy of a semi-classical ideal gas is significantly increased by subdividing a large homogeneous system into a heterogeneous ensemble of small subsystems, which reduces the free energy. A similar concept also applies to Ising’s model. In real systems, the interaction energy $J_i=J$ between neighboring spins usually involves quantum mechanical exchange that requires particles to be statistically indistinguishable. In this case, $J_i=0$ would come from breaks in the interaction due to neighboring spins that are distinguishable, yielding the increase in net entropy that lowers the free energy and gives the stable thermal equilibrium.

As with Ising’s model, it is instructive to calculate the change in free energy per particle due to subdividing the semiclassical ideal gas. Starting from the thermodynamic limit ($N_t\to \infty$) the answer is

$$∆F_t/\left(N_{t}kT\right)={\int }_0^{\overline{\eta }}\left(\mathcal{E}/kT\right)(d\eta /N_t) =-{\int }_{\infty }^{\overline{N}}\left(\mathcal{E}/kT\right)dN/N^2\mathrm{br}\mathrm{to}\mathrm{break}=\left[\left(\overline{N}-1\right)\ln \left(\overline{N}-1\right)-\overline{N}\ln \left(\overline{N}\right)\right]/\overline{N} .$$

(17)

The red line in Figure 5 shows Equation (17) as a function of $\overline{N}$. The minimum at $\overline{N}=2$ occurs when most particles are distinguishable by belonging to separate subsystems, which significantly increases the entropy and lowers the free energy. Using $\overline{N}=2$, Equation (17) gives $∆F_t/\left(N_{t}kT\right)=-\ln \left(2\right)$, yielding the net reduction in free energy per particle when a semiclassical ideal gas subdivides into its stable thermal equilibrium of subsystems.

5. Excess Energy Fluctuations in MD Simulations

5.1. Equilibrium Averages and Fluctuations in Statistical Mechanics

Boltzmann often sought to connect classical mechanics to statistical mechanics, with mixed success. Now, molecular dynamics (MD) simulations provide a powerful tool for investigating Boltzmann’s goal, whether classical forces and Newton’s laws may yield thermal behavior in complex systems. In fact, the first MD simulation on a digital computer tested whether anharmonic interactions between idealized atoms would yield average energies that obey the equipartition theorem characteristic of thermal equilibrium. Specifically, in 1953, Fermi, Pasta, Ulam and Tsingou (FPUT) used what was then the world’s most powerful computer to simulate 1-D systems of up to 64 particles [45]. Although they knew that harmonic modes would not mix, they thought that adding anharmonic interactions would allow the system to evolve towards thermal equilibrium. Instead, they found that there was “very little, if any, tendency toward equipartition of energy among the degrees of freedom.”

The failure of the FPUT model to reach thermal equilibrium has led to various interpretations. One view is that the FPUT model is too small and simple to yield the stochastic dynamics needed for thermalization, at least when simulated at low energies [46,47,48]. Indeed, MD simulations of large systems often show good agreement with the equipartition theorem expected for thermal equilibrium. However, when Fermi’s idea of using MD simulations to study thermal averages is extended to study thermal dynamics, significant deviations are found from the standard fluctuation relation that is based on the equipartition theorem [27]. Thus, fluctuations are a more-stringent test of thermal behavior than equilibrium averages. One way to investigate fluctuations in small systems is to study large simulations that are subdivided into smaller subsystems or “blocks” [27,44,49,50]. Because there are interactions between the blocks, they cannot be treated as an ensemble of independent subsystems needed for Hill’s theory. Nevertheless, some ideas from small-system thermodynamics may still apply, such as thermal behavior that depends on the choice of ensemble and type of dynamics.

A basic test of whether MD simulations yield standard statistical mechanics comes from the connection between energy fluctuations and the equipartition theorem via specific heat ($C=\partial E/\partial T$). The canonical ensemble has a formal identity that makes this connection, which we call the $C\leftrightarrow {\left(∆E\right)}^2$ relation. It is useful to review the derivation of this relation. The terminology is from Ref. [28]. Let $Q$ denote either the kinetic ($K$) or potential ($U$) energy in a block inside a larger system. Average values of $Q$ can be obtained in two ways. From theory, the thermal average (denoted by angled brackets with subscript $th$, ${〈Q〉}_{th}$) is obtained in standard statistical mechanics by Boltzmann averaging. From simulations, the average is obtained by averaging over all time-steps then taking the arithmetic mean over all blocks in the system (denoted by angled brackets with subscripts $t,bl$, ${〈Q〉}_{t,bl}$). If the ergodic hypothesis holds, and if simulations and theory yield equivalent results, these two averages will equal each other,

$${〈Q〉}_{t,bl}={〈Q〉}_{th}.$$

(18)

The $C\leftrightarrow {\left(∆E\right)}^2$ relation is derived in canonical theory by taking the derivative of the thermal average with respect to $kT$, yielding a dimensionless heat capacity, then using the product rule of calculus to obtain the variance for the fluctuations,

$${\displaystyle \frac{\partial {〈Q〉}_{th}}{\partial \left(kT\right)}}={\displaystyle \frac{{\left({∆}_{th}Q\right)}^2}{{\left(kT\right)}^2}} ({\left({∆}_{th}Q\right)}^2\equiv {〈Q^2〉}_{th}-{〈Q〉}_{th}^2).$$

(19)

Equation (19) applies to any system with thermal-average energies that are determined solely by Boltzmann’s factor, making this $C\leftrightarrow {\left(∆E\right)}^2$ relation a very basic test for thermal behavior. We ask the question: if Equation (18) connects time-averages to thermal-averages, does Equation (19) hold for time-averages? Specifically, does the $C\leftrightarrow {\left(∆E\right)}^2$ relation apply to time-averages,

$${\displaystyle \frac{\partial {〈Q〉}_{t,bl}}{\partial \left(kT\right)}}={\displaystyle \frac{{\left({∆}_{t,bl}Q\right)}^2}{{\left(kT\right)}^2}} ? ({\left({∆}_{t,bl}Q\right)}^2\equiv {〈{\left(Q\left(t\right)-{〈Q〉}_t\right)}^2〉}_{t,bl}.$$

(20)

MD simulations of several models establish that Equation (18) is accurately obeyed for average energies, so that the ergodic hypothesis applies to time- and thermal-averages. However, at sufficiently low temperatures all these models show clear deviations from Equation (20).

5.2. Thermal Fluctuations and MD Simulations

MD simulations of the Lennard-Jones (L-J) model were made on systems of 442,368 atoms using parameters for crystalline argon [27,51]. The systems were subdivided into blocks containing 1-2048 atoms to study local thermal behavior. Each block gives excellent agreement with the equipartition theorem for average thermal behavior. One example is shown in Figure 6. From $T<10$ K, the black line (with error bars) near unity gives ${\delta }_T{〈U〉}_{t.bl}/{\delta }_T{〈K〉}_{t,bl}=0.9998\pm 0.0013$ for the ratio of partial specific heats from potential and kinetic energies. Here, ${\delta }_T$ denotes the temperature difference, e.g., ${\delta }_{T}U\equiv \left[U\left(T+\delta T\right)-U\left(T\right)\right]$. This ratio of partial specific heats is well within statistical uncertainty of unity, indicating that both $U$ and $K$ obey the equipartition of energy.

Symbols in the main part of Figure 6 show the temperature dependence of energy fluctuations for block sizes given in the legend. The blocks are a virtual construct, with an equal number of interacting neighbors for every atom. The $C\leftrightarrow {\left(∆E\right)}^2$ relation predicts that energy fluctuations should be proportional to specific heat, whereas the ratio ${〈{\left({∆}_{t}U\right)}^2〉}_{t,bl}/{〈{\left({∆}_{t}K\right)}^2〉}_{t,bl}$ tends to diverge as $T\to 0$, reaching at least 10 times the ratio of unity. These results establish that MD fluctuations deviate from canonical-ensemble behavior. The inset of Figure 6 shows why: time-dependent correlations within blocks (black) and between blocks (red) are strongly anticorrelated. Thus, most of the energy in each fluctuation comes from the surrounding shell of neighboring blocks, not from the large heat bath needed for the canonical ensemble. Specifically, the black line and symbols show auto correlations as a function of time, ${\mathsf{\Gamma }}_{\mathrm{0,0}}\left(\tau \right)/{\mathsf{\Gamma }}_{\mathrm{0,0}}\left(0\right)$, where ${\mathsf{\Gamma }}_{\mathrm{0,0}}\left(\tau \right)\equiv {〈U_0\left(t\right)U_0\left(t+\tau \right)〉}_{t,bl}$. Due to the normalization ${\mathsf{\Gamma }}_{\mathrm{0,0}}\left(\tau \right)/{\mathsf{\Gamma }}_{\mathrm{0,0}}\left(0\right)$ starts at 100%, then exhibits damped oscillations as a function of $\tau$. The red line and symbols show pair correlations between each central block and its six nearest-neighbor blocks, ${\mathsf{\Gamma }}_{\mathrm{0,1}}\left(\tau \right)\equiv {〈U_0\left(t\right)U_1\left(t+\tau \right)〉}_{t,bl}$, which are strongly anticorrelated with the autocorrelations. Indeed, at $\tau =0$, more than 50% of the energy of the fluctuation comes from the six nearest neighbor blocks. Thus, when the energy from all $3^3-1=26$ blocks in the local shell of blocks that touch the central block are included, essentially 100% of the initial energy in each fluctuation comes from this local shell, not the large heat bath needed for Boltzmann’s factor. A theory for these localized fluctuations based on fast energy modulations yields the dashed lines in the main part of Figure 6 that mimic the simulations with no adjustable parameters [28]. We now outline the theory.

The key assumption in the theory is that fast modulations of interactions between next-nearest-neighbor (nnn) atoms are effectively adiabatic, so that they are time-averaged before thermal averaging. The theory starts with the net potential energy of an atom in a crystal as a function of the displacement ($x$) from its equilibrium position ($x=0$) and as a function of time ($t$) from ground-state modulations due to the motion of nnn atoms,

$$u\left(x,t\right)=a{x}^2+bc\left(t\right).$$

(21)

Referring to Figure 7, the quadratic term ($a$) comes from interactions between nearest-neighbor (nn) atoms, and the linear term $b$ from interactions between nnn atoms. Now, assuming harmonic motion of the nnn atoms, the modulation term can be written as

$$c\left(t\right)=-\mathcal{A}\left[\sin (\omega t+\alpha )+\sin \omega t\right].$$

(22)

For classical MD simulations that have equipartition of energy, the modulation amplitude, $\sqrt{{〈{\left({\mathrm{A}}_{\pm }\right)}^2〉}_{\mathrm{t}\mathrm{h}}}\equiv \mathcal{A}$, is given by the thermal average amplitude

$$\mathcal{A}=\sqrt{k_{B}T/\left(2a\right)}.$$

(23)

The main mechanism causing deviations from the $C\leftrightarrow {\left(∆E\right)}^2$ relation can be traced to the fact that the time average of Equation (22) is zero, whereas the time average of ${\left[c\left(t\right)\right]}^2$ is not zero. This modifies the expression for potential energy fluctuations [28] so that Equation (20) becomes

$${\displaystyle \frac{{\left({∆}_{t,bl}U\right)}^2}{{\left(k_{B}T\right)}^2}}={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{b^2}{akT}} .$$

(24)

The interaction parameters for the L-J potential in its ground state fcc lattice are $a=144{ϵ}_{LJ}/\left(2^{1/3}{\sigma }^2\right)$ and $b=-63{ϵ}_{LJ}/\left(2^{14/3}\sigma \right)$ [28]. For argon atoms, the energy scale is ${ϵ}_{LJ}/k_B=118.2$ K and the length scale is $\sigma =0.3405$ nm [51]. Utilizing a proposed connection between anomalous diffusion and the block-size dependence for energy fluctuations, the theoretical expression for potential energy fluctuations of a block of $n$ atoms is

$${\displaystyle \frac{{\left({∆}_{t,bl}U\right)}^2}{{\left({∆}_{t,bl}K\right)}^2}}-1={\displaystyle \frac{{21}^2}{2^{13}}}{\displaystyle \frac{{ϵ}_{LJ}}{k_{B}T}}{n}^{1/4}\approx 0.053833{\displaystyle \frac{118.2}{T\left[K\right]}}{n}^{1/4}.$$

(25)

The dashed curves in Figure 6 show that Equation (25) quantitatively characterizes the excess fluctuations.

Note the qualitative similarity between simulations having excess fluctuations shown in Figure 6 and measurements having excess specific heat shown in Figure 2. This connection is quantified in Ref. [28]. Thus Equation (24) replaces the usual $C\leftrightarrow {\left(∆E\right)}^2$ relation for systems having fast dynamics that must be time-averaged before thermal-averaging. Similar separation of time scales and isolation from the heat bath are characteristic of measurements of thermodynamic heterogeneity, as described in Section 3, sketched in Figure 2 and Figure 3, and deduced from NMR measurements (Figure 4a). From these measurements, localized degrees of freedom often exhibit effective temperatures that differ significantly from the homogeneous heat bath, as found in liquids, glasses, spin glasses, polymers and crystals. Thus, one feature of nanothermodynamics that is crucial to the theory for fluctuations in MD simulations is the ubiquity of thermodynamic heterogeneity. Another feature is that appropriate ensembles must be used to ensure accurate behavior, especially for fluctuations.

6. Nanoscale Fluctuations and the Second Law of Thermodynamics

6.1. A Creutz-Like Model for Investigating Nanoscale Fluctuations

Nanothermodynamics impacts several aspects of nanoscale fluctuations. Two consequences of including finite-size thermal effects are that both the choice of ensemble [52] and the choice of dynamics [17] influence the fluctuations. Thus, if all sources of energy from finite-size thermal effects are included, fluctuations exhibit features that are not contained in standard theories for thermal fluctuations that are based on microscopic reversibility and ensemble equivalence [53].

A statement of the second law of thermodynamics is that the entropy of a closed system will most likely increase monotonically until it reaches its maximum value, which defines the thermal equilibrium. Thus, investigations of the 2nd law are ideally done using closed systems, having an explicit heat bath that thermalizes the system with no external influences [54]. Perhaps the simplest such model, adapted [17] from the work of Creutz [16], involves 1-D chains of Ising spins coupled to a bath of Einstein oscillators (that Creutz called “demons”). This simplified 1-D Creutz model was first used to study how the 2nd law applies not only to large and complex systems, but also to small and simple systems [17]. Sketches of the model are shown on the right-hand side of Figure 5 (above) and repeated in Figure 8b (below). An advantage of this model is that the exact entropy of the system and its bath can be calculated at every step. Moreover, analytic expressions for the stable equilibrium of the 1-D Ising model are known, as given in Section 4.2 (above). Here we will show that, at least for this model, the 2nd law requires an intrinsically irreversible microscopic step, providing a clear counterexample to the usual assumption that large systems have thermal behavior that emerges from intrinsically reversible microscopic dynamics. Another result shown by this model is that fluctuations in total entropy are best described by the 2nd law of thermodynamics, not by standard statistical theories based on Boltzmann’s factor.

6.2. Reversible and Irreversible Simulations of the 1-D Creutz Model

The simulations shown in Figure 8 come from a system of $N=\mathrm{250,000}$ spins with an equal number of demons. The closed system contained a fixed amount of energy ($1/2$ unit per spin), shared between the spins and demons, yielding an average temperature of $kT/J=1.30$. Again, a crucial feature of the 1-D Creutz model is that the exact entropy of the spin system and its bath can be calculated at every step. The solid lines in Figure 8 show the total entropy (from spins and demons, $S_t$) per spin as a function of time from simulations with two types of dynamics: reversible (black) and irreversible (red). Both simulations evolved from the same initial spin configuration that was highly ordered (low entropy). Details about the different dynamics are given in [17]. Briefly, reversible dynamics utilizes three randomly generated (but fixed) arrays: one array sets the sequence of spins to be updated while the other two arrays govern which two demons couple to each spin. In contrast, irreversible dynamics utilizes a new random number for each of these steps.

Figure 8a shows that at the start of the simulations, $S_t$ tends to increase with time, as expected from the 2nd law of thermodynamics for closed systems that are initially out of equilibrium. However, the nature of this increase is strongly dependent on the dynamics. Irreversible dynamics (red) yields a rapid and mostly monotonic increase, while reversible dynamics (black) shows slower relaxation with large oscillations. Figure 8c shows that at the end of the simulations the differences are even more pronounced. Of course, $S_t$ from irreversible dynamics (red) stays at or near its maximum ($S_m$), as expected from the 2nd law. However, for reversible dynamics (black) the simulation exhibits exact reversibility, with every step down in entropy that was taken up in entropy at the start, a clear demonstration of Loschmidt’s paradox for dynamics that is reversed at the midpoint of its evolution (not shown).

Most standard resolutions of Loschmidt’s paradox rely on the fact that for large and complex systems it is difficult to exactly reverse the dynamics at the midpoint of an evolution. In other words, a large isolated system has a vanishingly small fraction of initial conditions that cause $S_t$ to systematically decrease. However, as shown in Figure 8, it is relatively easy to find special sets of initial conditions for reversible systems [55] of any size by simulating from an initially low-entropy state before reversing the dynamics. More importantly, there are fundamental reasons why the standard argument about initial conditions cannot explain the 2nd law, at least for the 1-D Creutz model. First, from Figure 8b, $S_t$ from irreversible dynamics remains significantly higher than $S_t$ from reversible dynamics, regardless of the initial state. Second, as shown in [17] (Figure 2a), after an irreversible simulation $S_t$ always decreases when the dynamics is switched to reversible. Finally, the difference in entropy density between reversible and irreversible dynamics increases with decreasing system size, as shown in [17] (Figure 3). Thus, the usual assumption that only special initial conditions cause entropy to decrease in reversible systems cannot solve Loschmidt’s paradox, at least not for the 1-D Creutz model.

6.3. Comparison to Standard Theories for Fluctuations

Most theories of equilibrium fluctuations predict Gaussian distributions for thermal quantities. This, of course, comes from the central limit theorem for large numbers of uncorrelated, randomly distributed quantities with a finite variance. Standard derivations [56] start by expanding the entropy change as a function of the quantity. Specifically, for fluctuations in energy, $∆E$, assuming analyticity, from an expansion of the change in entropy, $p(∆E)\propto \exp \left(∆S/k\right)\approx \exp \left[\left(∆E/k\right)(\partial S/\partial E)+{\displaystyle \frac{1}{2}}(∆{E/k)}^2({\partial }^{2}S/\partial E^2)\right]$. Thermodynamic identities yield $\partial S/\partial E=1/T$ and ${\partial }^{2}S/\partial E^2=-1/\left(C{T}^2\right)$. Linear changes in energy in systems of any size are accommodated by the heat bath (subscript $b$), $∆E_b=-∆E$, so that the first term in the expansion becomes Boltzmann’s factor. The second term yields a Gaussian distribution for equilibrium energy fluctuations. This Gaussian distribution has many applications, including characterizing the slow equilibrium dynamics of supercooled liquids [14].

Fluctuation theorems have extended the range for theoretical treatments of thermal fluctuations to non-equilibrium conditions. Various fluctuation theorems are described in [53,57,58]. Applications of the theorems to experiments include characterizing the free energy change while stretching RNA [59] and the fluctuations of colloidal particles [60,61]. Of particular relevance for analyzing the 1-D Creutz model are detailed fluctuation theorems that apply to the total entropy of a closed system [18,19,20,21]. In general, these theorems yield the probability that the total entropy increases (positive entropy production) divided by the probability that it decreases (negative entropy production), as given by equation:

$$p(∆S_t)/p(-∆S_t)=\exp \left(\alpha ∆S_t\right).$$

(26)

Here, $\alpha =1$ for a driven system measured over a time interval [18] and along various stochastic trajectories [19,21]. In general, $\alpha$ may depend on a diffusion constant, $D$, and on the nonlinearity of the force [20], yielding $\alpha \ne 1$. Note that Gaussian distributions invariably obey Equation (26). Specifically, using $p(∆S_t)\propto \exp \left[-{\left(∆S_t-B\right)}^2/{\omega }^2\right]$ gives,

$$p(∆S_t)/p(-∆S_t)=\exp \left(4B∆S_t/{\omega }^2\right),$$

(27)

in agreement with Equation (26) if $\alpha =4B/{\omega }^2$. A remarkable feature of fluctuation theorems (that will not be addressed here) is that they may also apply to non-equilibrium conditions.

To compare theoretical treatments of fluctuations to simulations of the 1-D Creutz model it is useful to use histograms that convert entropy as a function of time from Figure 8 to distributions of multiplicity as a function of entropy, Figure 9. Two distinct types of distributions are found. Reversible dynamics (filled symbols) exhibit distributions that are approximately Gaussian, consistent with the central limit theorem when $S_t$ can fluctuate freely up and down. However, with decreasing system size the Gaussians become increasingly skewed towards low entropy, opposite to standard fluctuation theorems that favor high-entropy states. From Figure 9 It is evident that this asymmetry increases with decreasing system size because larger fluctuations allow $S_t$ to approach the maximum entropy ($S_m$) shown by the irreversible dynamics (open symbols). Indeed, this $S_m$ sets an effective upper limit for fluctuations in $S_t$. More importantly, equilibrium (irreversible) simulations have only asymmetric fluctuations, yielding one-sided distributions that are characterized by

$$p\left(S_t\right)=A\exp (S_t/k) (S_t\le S_m).$$

(28)

Here, the normalization constant for a continuous distribution is $A=1/\left[\exp \left(S_m/k\right)-1\right]$. The solid lines in Figure 9 show good agreement between Equation (28) and the open symbols from irreversible dynamics. Equation (28) is consistent with Einstein’s reversal of Boltzmann’s relation that he used to describe equilibrium fluctuations [62]. Equation (28) may also be consistent with Equation (26) if $\alpha =2$, except at the limit of $S_t=S_m$ where entropy can only decrease.

The asymmetry in the distribution of $S_t$ during irreversible dynamics gives insight into the physics of equilibrium thermal fluctuations. Standard fluctuation theories predict that entropy increases are more likely than entropy decreases. Equation (28) adheres to this prediction for most values of $S_t$, except at $S_m$ where $S_t$ can only decrease. This maximum in entropy is evident from the one-sided fluctuations shown by irreversible dynamics (red line) in Figure 8b, and by the one-sided distributions in Figure 9. Although the one-sided distributions differ from usual predictions of standard statistical mechanics, this asymmetry is needed to yield maximum entropy at the highest multiplicity, another indication that only irreversible dynamics is consistent with the 2nd law of thermodynamics.

6.4. Effective Temperatures from the 1-D Creutz Model

Because even von Neumann recognized that “nobody really knows what entropy is” [63], it is useful to go beyond the concept of maximum entropy and the 2nd law to further justify the need for irreversible dynamics. An appealing feature of the 1-D Creutz model is that it utilizes an explicit heat bath of demons to yield local temperature as an output of the simulation, not as an input, as indicated by the sketch in Figure 8c. One application is to use the distribution of energies in each oscillator for location-specific temperatures that are ideal for characterizing the thermodynamic heterogeneity measured by TRSH and NSHB described in Section 3. Another application is to use the relative occupation of adjacent energy levels to obtain an energy-specific temperature that is ideal for assessing the accuracy of the Boltzmann distribution. Indeed, because the demons that form the thermal bath have equal energy spacing, $J$, Boltzmann’s factor gives the relative probability of energy level $i$ as, $p_i\propto \exp (-iJ/kT)$. Thus, an energy-specific effective temperature ($T_i$) between adjacent levels is given by

$$k{T}_i/J=1/\ln \left(p_i/p_{i+1}\right).$$

(29)

Figure 10 shows the energy-level dependence of this effective temperature from simulations of the 1-D Creutz model with three system sizes and both types of dynamics. Irreversible dynamics gives a single well-defined temperature for all levels. Reversible dynamics gives an effective temperature that oscillates between even and odd levels. Again, reversible dynamics differs significantly from expected thermal behavior. At the very least, reversible dynamics yields a non-equilibrium steady state condition between the system and its bath, even when there are no external forces. Indeed, as shown in [17] (Figure 2), it is this failure of the heat bath to follow Boltzmann statistics that is the main cause of reduced $S_t$ for reversible dynamics, as shown here in Figure 8 and Figure 9.

6.5. Conclusions from Simulations of the 1-D Creutz Model

Three questions arise regarding the behavior shown in Figure 8, Figure 9 and Figure 10: what is the role of nanothermodynamics, can the simplistic 1-D Creutz model be connected to real-world systems, and what might be the source of intrinsically irreversible dynamics in nature? For the first question, all results presented here yielding maximum entropy require that the Ising model be able to access the stable thermal equilibrium of the nanocanonical ensemble. In fact, if the chain of Ising spins is constrained to be homogenous (depicted by the top sketch in Figure 5) then there is negligible difference between reversible and irreversible dynamics. Furthermore, total entropy is never maximized so that the distributions are not one-sided, and there are no odd energy levels with different effective temperatures for Figure 10. Thus, the stable thermal equilibrium of nanothermodynamics is essential for studying the influence of the 2nd law on fluctuations.

An overall answer to the second question is that Creutz demons play a crucial role in enhancing physical reality, which cannot be matched by any idealized heat bath such as the Metropolis algorithm or Langevin equation. On a practical level, random numbers are an inherent feature of idealized heat baths so that their simulations cannot be made reversible. Furthermore, as shown in [17] (Figure 2), entropy of the spin system decreases for irreversible dynamics: the increase in $S_t$ needed for the 2nd law (shown here in Figure 8 and Figure 9) requires the explicit heat bath. Similarly, an explicit heat bath and maximum entropy are also needed for a consistent local temperature that is not externally imposed on the system, as shown in Figure 10. On a fundamental level, the demons ensure that local energy is always conserved, unlike idealized heat baths where changes in local energy are assumed to be instantaneously absorbed by an infinite bath. This fundamental aspect also connects the 1-D Creutz model to MD simulations (Section 5) because both types of simulations are done in the most-basic microcanonical ensemble. Indeed, conservation of local energy is evident in the insert of Figure 6 showing that potential energy is traded back-and-forth between each central block and its surrounding shell of neighboring blocks. It is these local oscillations that insulate fast degrees of freedom from the heat bath, preventing them from having a well-defined temperature. Similarly, the oscillations in $S_t$ from the 1-D Creutz model with reversible dynamics (clearly visible in Figure 8a and Figure 8c) can be attributed to energy that is repeatedly traded back-and-forth between each spin and its heat bath of demons. This interpretation is consistent with the observation that the oscillation amplitude decreases with increasing number of demons that connect to each spin, as shown in [17] (Figure 1). Whereas, for irreversible dynamics, the random connection between spins and multiple demons efficiently disperses the initially localized energies. In contrast, reversible dynamics causes all simulations to exhibit persistent oscillations that inhibit the systems from reaching states of maximum entropy. Thus, the evidence indicates that an explicit heat bath more-accurately represents reality than idealized heat baths, and that intrinsically irreversible dynamics is a common requirement of any closed system that is to obey the 2nd law of thermodynamics.

So far there is no definitive answer to the third question. Although irreversible simulations of the 1-D Creutz model utilize a pseudo-random number generator, what is the source of randomness in real systems? In reference [17], it is speculated that randomness in nature could come from the quantum measurement process. Specifically, if the process connecting the system to its bath is irreversible, then the precise bath particle to be used is unknown until the connection is made. This connection is particularly impactful if the bath contains indistinguishable particles that become distinguishable when the measurement is made. Indeed, it is the nature of this connection between system and bath that is the only difference between reversible and irreversible simulations shown in Figure 8, Figure 9 and Figure 10. Presumably, it is the fact that random connections have higher entropy than reversible connections that allows the system to reach its maximum entropy. Future work will be necessary to understand how the connection between a system and bath controls the behavior of the combined system.

7. Summary

In this review it is emphasized that large systems usually favor subdividing into stable distributions of internal subsystems, which is the basis of nanothermodynamics. In Section 2, the history, foundations, and significance of nanothermodynamics are discussed. Figure 1 summarizes the progression of the fundamental equation of thermodynamics, culminating in Hill’s subdivision potential, $\mathcal{E}$, that is needed for the systematic treatment of finite-size thermal effects.

In Section 3, experimental evidence is given for thermodynamic heterogeneity inside bulk materials due to a distribution of effective local temperatures, $T_i$. The evidence for distinct values of $T_i$ during the response of most types of materials – including liquids, glasses, spin glasses, polymers, and crystals – justifies the need for nanothermodynamics. Figure 2 shows results from measurements of $∆T_i$ as a function of time on high-purity single crystals, yielding equilibrium specific heats that significantly exceed Debye’s theory. The essentially inverse process of nonresonant spectral hole burning establishes that there is negligible correlation between the slow and fast degrees of freedom in a spin glass (Figure 3) and a glass-forming liquid (Figure 4). Furthermore, Figure 3 shows that $T_i$ may exceed the bath temperature by 2 K, or more, heating a significant fraction of the sample into the paramagnet phase, while the rest of the sample remains a spin glass. Figure 4a is a picture that emerges from measurements of nuclear magnetic resonance indicating that the thermodynamic heterogeneity in glass-forming liquids occurs on length scales of nanometers.

Section 4 outlines the need for nanothermodynamics to solve two standard models in statistical physics. First, the stable thermal equilibrium of Ising’s original model for finite chains of spins is derived two ways, using two different ensembles starting from two different Hamiltonians. It is emphasized that standard thermodynamics and statistical mechanics give no guidance on how to find this equilibrium. Second, the stable equilibrium of the semi-classical ideal gas is described, which is needed to ensure that entropy is extensive for systems of any size. This novel solution to Gibbs’ paradox may also explain the lack of isotope effects in measured entropies of dilute gases. Figure 5 shows how minimizing the total free energy requires subdividing large systems into the stable distribution of small subsystems, which is the basis of nanothermodynamics.

In Section 5, MD simulations are used to investigate the excess energy fluctuations found at low temperatures (Figure 6) that are qualitatively similar to the excess specific heat measured in crystals (Figure 2). A theory is described based on fast (adiabatic) modulation of ground state energies due to localized fluctuations, as sketched in Figure 7. This theory is consistent with the anticorrelated local energies in the simulations (inset of Figure 6), and with measurements of thermodynamic heterogeneity. The theory gives quantitative agreement with the excess fluctuations, and a correction to the fluctuation relation for degrees of freedom that are too fast to fully couple to the heat bath, Equation (24).

In Section 6, simulations of a 1-D Creutz model are used to investigate nanoscale thermal fluctuations and the 2nd law of thermodynamics. This model is comprised of chains of Ising-like spins coupled to an explicit heat bath of Einstein oscillators. The model has several appealing features: the total energy is always exactly conserved, the spins can reach their stable thermal equilibrium, the heat bath yields the distribution of local temperatures, the dynamics can be made fully reversible, and the total entropy ($S_t$) can be calculated exactly at every step. As shown in Figure 8, maximizing $S_t$ requires an intrinsically irreversible microscopic step, providing a counterexample to the usual claim that thermal equilibrium can emerge from reversible dynamics. Figure 10 shows that this irreversible step is also needed for a well-defined temperature. Figure 9 shows that the distribution of $S_t$ from fluctuations is strongly dependent on the dynamics. Equilibrium fluctuations are best described by Einstein’s reversal of Boltzmann’s relation and the 2nd-law, not by any recent fluctuation theorem. Section 6 finishes with a discussion about the relevance of the 1-D Creutz model to real systems, and a possible source of intrinsic irreversibility in nature.

In conclusion, many measurements establish the necessity of nanothermodynamics for the thermal and dynamic behavior of most types of materials. Nanothermodynamics is crucial in finding the stable thermal equilibrium of standard models from statistical physics. Deviations of MD simulations from a fluctuation relation can be quantitatively characterized by including fast degrees of freedom from localized modes that are effectively decoupled from the heat bath. Nanoscale fluctuations in a simple model for a closed system are best described by Einstein’s reversal of Boltzmann’s relation and the 2nd law of thermodynamics, not by recent fluctuation theorems. It is anticipated that the stable thermal equilibrium introduced by Hill 62 years ago, and its adaptation to nanothermodynamics, will continue to find applications that give new insights into experiments, theory, and simulations.