- Email: [email protected]
Monte Carlo simulations including energy from an entropic force
Physica A 391 (2012) 5384–5391
Contents lists available at SciVerse ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Monte Carlo simulations including energy from an entropic force Ralph V. Chamberlin ∗ Department of Physics, Arizona State University, Tempe, AZ 85287-1504, USA
article
info
Article history: Received 8 August 2011 Available online 27 June 2012 Keywords: Entropic force Heterogeneous dynamics Ising model Monte Carlo simulations Small-system thermodynamics
abstract Several experimental techniques have shown that the primary response of many materials comes from a heterogeneous distribution of independently relaxing nanoscale regions; but most Monte Carlo simulations have homogeneous correlations. Resolving this discrepancy may require including the energy needed to change the configurational entropy, which is often used in theoretical treatments of thermal fluctuations, but not in computer simulations. Here the local configurational entropy is shown to give a nonlinear correction to the Metropolis algorithm that restores conservation of energy, maintains maximum entropy, and yields heterogeneous correlations. The nonlinear correction also improves agreement between Monte Carlo simulations of the Ising model and measurements of specific heat and structural correlations from the Jahn–Teller distortion in LaMnO3 . © 2012 Elsevier B.V. All rights reserved.
1. Introduction In most computer simulations of classical models, the energy fluctuations in each local subsystem, and resulting local entropy (SL ), fail to increase linearly with increasing subsystem size, NL [1,2]. Although superficially similar to the usual Gibbs’ paradox for non-extensive entropy in large systems of distinguishable particles [3,4], the non-extensive energy fluctuations in most computer simulations arise from homogeneous correlations. However, several experimental techniques have shown that the primary response of many materials comes from a heterogeneous distribution of independently relaxing nanometer-sized ‘‘regions’’, with negligible correlations between the regions [5–8]. These experiments often show that the heterogeneity is dynamical: occurring transiently in the thermal-equilibrium fluctuations of single-component substances having homogeneous average density and uniform interactions. For most computer simulations of independent regions in these materials, either the dynamics needs a correction to make the entropy additive, or the foundation of statistical mechanics needs to be modified to accommodate non-extensive entropy [9,10]. Key insight comes from a quantitative analysis of non-resonant spectral hole burning measurements showing that the heat capacity of these regions is extensive [11,12]. Thus, simulations exhibiting non-extensive energy fluctuations and uniform correlations cannot explain the measured response. The cellular method [13] is a theoretical tool for treating fluctuations by subdividing a bulk sample into small subsystems, so that each subsystem couples to a self-consistent ensemble of similar subsystems. The usual procedure is to include uniform coupling between subsystems to maintain homogeneous correlations, unlike the measurements showing heterogeneous dynamics. Here I explore a type of dynamics that includes a local configurational entropy, giving a nonlinear correction to the Boltzmann factor that decreases the coupling between subsystems, allowing otherwise homogeneous systems to exhibit response from an ensemble of effectively independent regions. The energy from this entropy is emergent [14], hence it does not appear in the interaction Hamiltonian, yet it is needed maintain conservation of energy from nanometer-sized fluctuations. Because the correction is dynamical, the system will eventually become
∗
Tel.: +1 480 965 3922; fax: +1 480 965 7954. E-mail address: [email protected].
0378-4371/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2012.06.016
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
5385
homogeneous, but at longer times and lower energies than the uncorrected system due to the heterogeneous response at intermediate times. Entropic forces are often used in the theory of Gaussian fluctuations [15,16], but most computer simulations focus on interparticle interactions and statistical probabilities, not statistical forces. Even when entropic forces are simulated [17,18], finite-size effects from the local entropy are usually neglected. Recently, a configurational entropy from the local order parameter ML has been shown to give a Gaussian correction to the Metropolis algorithm that lowers the energy, makes the entropy extensive, and improves agreement between computer simulations and measurements of ferromagnetic materials and critical fluids [19,20]. Here I add a similar Gaussian correction to Monte Carlo simulations of the antiferromagnetic (AF) Ising model. I show that away from the critical temperature Tc , the correction reduces the correlations between regions, yielding heterogeneous dynamics in the simulation. Moreover, the correction improves agreement between the measured specific heat and structural correlations found directly from the Jahn–Teller distortion in LaMnO3 . 2. Background 2.1. Thermal properties of small regions Let NL be the number of particles in a subsystem having interaction energy UL . In most computer simulations, local energy fluctuations are non-extensive, i.e. the thermal average ⟨UL2 ⟩ − ⟨UL ⟩2 increases nonlinearly with increasing NL [1,2]. Thus, the resulting heat capacity CV = (⟨UL2 ⟩ − ⟨UL ⟩2 )/(kB T 2 ) [21] and local entropy difference between two temperatures
T
SL (T2 ) − SL (T1 ) = T 2 (CV /T ) dT , are also non-extensive. Consequently, the entropy of a combined subsystem a + b does 1 not equal the entropy of its constituents: Sa+b ̸= Sa + Sb . For homogeneous systems, this non-additive entropy is expected from uniform correlations [9]. Specifically, using the Boltzmann–Gibbs expression for entropy S = −kB Σi pi ln(pi ) with kB Boltzmann’s constant and pi the probability of each state, correlations prohibit separating the probability into a product: pa+b ̸= pa ∗ pb . However, the primary response of many materials is now known to come from independently relaxing nanometer-sized regions that are effectively uncorrelated [5–8]. For this response, to a good approximation pa+b = pa ∗ pb , so that Sa+b = Sa + Sb , unlike most simulations. Fundamental transition rates in thermostatistics are given by w ∝ e1St /kB , where 1St is the change in total entropy. First consider two independent contributions to entropy, one from the change of a single particle (1S1 ) and one from an infinite bath (1S ∗ ), so that 1St = 1S1 + 1S ∗ . The Boltzmann factor is obtained by assuming that the single particle changes from a known initial state to a known final state, so that 1S1 = 0. Then the total entropy change is assumed to come only from the change in energy of the infinite bath 1U ∗ , so that using the fundamental equation of thermodynamics (combined first and second laws) gives 1St = 1U ∗ /T . The minus sign in the Boltzmann factor comes from replacing 1U ∗ by the change in energy of the single particle (1U1 ) using conservation of energy 1U ∗ = −1U1 , giving w ∝ e−1U1 /kB T . Strictly speaking, the standard Boltzmann factor assumes instantaneous thermal contact to an infinite heat bath. Here I consider some consequences of treating an independently relaxing region as a small local bath, and including configurational entropy as another source of heat. A more complete expression for the change in total entropy is 1St = 1S1 + 1S ∗ + 1SL [20]. For a single particle changing between known states, 1S1 = 0 as before. For an effectively infinite heat bath, 1S ∗ = −1U1 /T as before. The new term 1SL is the local offset in entropy from the region containing the single particle. Note that adding 1SL to 1S ∗ is justified because the region relaxes independently, even if the rest of the sample forms a heat bath for the region. For large regions, 1SL may include entropy from particles that do not interact directly with the single particle; whereas for small regions, 1U1 often includes interactions that extend into neighboring regions. In any case, the change in energy rarely defines the total change in entropy; many different configurations in a region may yield the same 1U1 . In other words, it is unlikely that a single scalar quantity will define the locations and alignments of all particles in a region. Thus, if the statistical probabilities are to be based on the total entropy: a local configurational entropy must be added to the entropy of the infinite bath, giving a nonlinear correction to the Boltzmann factor. Another way to think of this correction is to include conservation of energy from work done on entropic forces. Entropic forces generally favor maximum entropy, such as the unaligned state of independent dipoles, or the uniform density of ideal gases. Even for ideal systems with no interactions, nonzero work is needed to displace the system from its maximum entropy. One example is the entropic restoring force that opposes stretching a polymer [16]. Another example is Fick’s law, which gives a restoring force proportional to the gradient in the particle density, so that the restoring force increases linearly with how far the particles have moved from a uniform distribution. The resulting emergent energy depends quadratically on the average offset of particles from equilibrium, even if there is no interaction between particles. Quadratic energies from entropic forces often form the basis of theoretical treatments of thermal fluctuations [15,16], but they are usually neglected in computer simulations. For example, the Metropolis algorithm weights the probability of configurations using only the interaction energy, whereas including the local configurational entropy can shift critical temperatures by nearly 50% ([20] and shown below). In systems of independently relaxing regions, a nonlinear correction to the Boltzmann factor is necessary if energy is to be conserved and total entropy is to be maximized during normal thermal fluctuations.
5386
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
2.2. Nanothermodynamics Although small-system thermodynamics [22] was originally developed to treat isolated small systems, similar nonextensive contributions to energy occur in ensembles of independently fluctuating regions. In this ‘‘nanothermodynamics’’, energy conservation is restored by adding a new pair of conjugate variables to the fundamental equation of thermodynamics. Thus, for magnetic systems:
1U = T 1S + H 1M + µ1N + E1η.
(1)
Here η ∼ N /NL is the total number of regions, E = ∂ U /∂η is the subdivision potential, and the other symbols have the usual meaning (T —temperature, S—entropy, H—magnetic field, and M—magnetic moment). This E can be understood by comparison to the chemical potential (µ): µ is the change in energy to add a single particle to the system, whereas E is the extra energy from adding a region of interacting particles to the system; and in general a region of NL interacting particles does not have the same energy as NL isolated particles. In other words, E contains all of the non-extensive and emergent energy terms that occur in finite-sized systems of interacting particles. Classical thermodynamics is based on the assumption that both energy and entropy are extensive, but finite-size effects alter the behavior of small systems. Two standard sources of non-extensive energy for internal regions are: interface terms (breaking correlations between regions costs energy proportional to the area of internal surfaces) and fluctuation effects. In general, the average energy of an independently fluctuating region is lowered inversely proportional to NL , which can be understood from mesoscopic mean-field theory [23–25] as follows. In thermal equilibrium at high temperatures, the meanfield energy of an infinite system is at its maximum because the entropy is also maximized. Subdividing the system into an ensemble of independent regions lowers the mean-field energy due to fluctuations. Mean-square fluctuations of intensive variables vary inversely proportional to the number of particles, so that the mean-field energy per particle is reduced inversely proportional to NL . Thermal equilibrium is found by balancing this energy reduction from increased fluctuations with the entropy reduction from finite-sized regions, thus minimizing the free energy of the generalized ensemble (µ, H , T ). Note that treating independent regions inside a bulk system using a canonical or grand-canonical ensemble artificially fixes the volume of each region. Also note that although this mean-field energy reduction is in the interaction, it comes from an entropic force that allows regions to fluctuate independently, so that the energy reduction emerges only for finite-sized regions. A similar 1/NL reduction in the average interaction energy is seen in computer simulations of the Ising model using the nonlinear correction to the Boltzmann factor that allows the regions to fluctuate independently [19]. It is often more difficult to justify non-extensive entropy. Furthermore, because the nonlinear correction modifies the dynamics, not the Hamiltonian, the instantaneous interaction energy remains additive at every step. Boltzmann’s definition of entropy involves a sum over all microstates, not a weighted average, so that the entropy of uncorrelated regions should be additive at all times. 2.3. Beyond the Boltzmann factor Now return to the fundamental expression for transition rates w ∝ e1St /kB . Here 1St = −1U1 /T + 1SL , from Eq. (1) for a fixed number of particles (1N = 0) in zero field (H = 0). Recall that 1U1 is the energy borrowed from the infinite bath to balance the change in interaction energy during single-particle dynamics. Eq. (1) gives the connection between the emergent energy and the change in entropy for the small regions 1SL = −E1η/T . Although this 1SL may contain several factors, the main contribution in many systems comes from the quadratic term 1SL = 0.5(∂ 2 SL /∂ ML2 )(1ML )2 , where 1ML is the offset from equilibrium in the local order parameter. I emphasize that because thermal equilibrium is required for Eq. (1), this 1ML is the offset from equilibrium due to normal thermal fluctuations, not the difference between two out-of-equilibrium states. Note that in a large enough external field, H 1M may dominate conservation of energy in Eq. (1), but if H ≈ 0 the entropic force that yields 1SL is the main source of magnetic work. Also note that this magnetic work exists even if H = 0, because ∂ 2 SL /∂ ML2 < 0 is necessary for stable fluctuations. The usual Boltzmann factor neglects all such nonlinear terms. If nonlinear magnetic work is included, the fundamental rate for changing the energy in a local region of a bulk sample comes from the usual Boltzmann factor e−1U1 /kB T , times a nonlinear correction from the entropy of the region wL = e1SL /kB . For realistic systems, nonlinear interactions and kinetic energy may also influence this correction. For the Ising model, using the binomial coefficient for the multiplicity of local alignments, the full correction to all orders is wL = (0.5NL )!(0.5NL )!/{[0.5(NL + 1ML )]![0.5(NL − 1ML )]!}. This nonlinear correction can also be understood by maximizing the entropy of a system and its heat bath. Using Lagrange multipliers, α for normalized probabilities and 1/kB T for heat flow from the thermal bath, the Boltzmann–Gibbs entropy of a uniform system is maximized by maximizing −kB Σi pi [ln(pi ) + α + Ui /kB T ], yielding pi ∝ e−Ui /kB T for the probability of a state of energy Ui . Subdividing the system into an ensemble of independently fluctuating regions lowers its entropy, but its energy is also lowered inversely proportionally to NL . This extra energy added to the heat bath usually increases the combined entropy of the system and its bath, as shown by the reduction in free energy [19]. The Boltzmann factor and canonical ensemble neglect this contribution from local configurational entropy, yielding thermal equilibrium only if the system can be thought of as isolated states that couple individually to a large heat bath. The total entropy from local configurations of interacting particles plus heat bath can be increased by subdividing the system into a generalized ensemble of independently fluctuating regions, with transition rates that include a nonlinear correction factor w ∝ e−1U1 /kB T e1SL /kB .
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
5387
Fig. 1. Schematic representations of binary particles that may be up (white) or down (black). Region (a) has an equal number of up and down particles giving 1SL = 0. Region (b) has fluctuated into a highly aligned state with 1SL /kB = ln{(8!8!)/[14!2!]} = ln(1/107.25). Both regions have a particle (marked by 1) with identical interaction energies; but (b) has a nonlinear correction factor wL = e1SL /kB , causing its inversion rate to be 107.25 times slower than that predicted by the Boltzmann factor alone. (c) 2 × 2 regions containing, from left to right: four-up, three-up, and two-up particle alignments. The statistical probabilities of each alignment are: p4 = 1/16, p3 = 4/16, and p2 = 6/16. At high temperatures, the classical Metropolis inversion rates between each pair of alignments are w4+ = 4/4 and w3− = 1/4, with w3+ = 3/4 and w2− = 2/4. Including the nonlinear correction factors, which depend only on the offset of entropy from its maximum, the inversion rates become ⟨w4+ ⟩ = 1/6 and ⟨w3− ⟩ = 1/6, with ⟨w3+ ⟩ = 1/2 and ⟨w2− ⟩ = 1/2. Thus, detailed balance gives equal probabilities ⟨p4 ⟩ = ⟨p3 ⟩ = ⟨p2 ⟩, and the alignment entropy −kB ⟨pL ⟩ ln⟨pL ⟩ does not rise and fall during equilibrium fluctuations.
The correction factor is not used in Onsager’s solution of the two-dimensional Ising model, which assumes homogenous correlations in the canonical ensemble. If energy reduction from independent nanometer-sized fluctuations is included, the Onsager solution fails to give the thermal equilibrium of the Ising model. 2.4. Examples The sketches in Fig. 1(a) and (b) show a specific example that emphasizes the importance of the nonlinear correction. Note that both (a) and (b) have an up site (marked by 1) with all four nearest neighbors down. For nearest-neighbor interactions, inverting these sites yields the same energy change 1U1 , so that the Boltzmann factor alone predicts identical inversion rates for particle 1 in both (a) and (b). Whereas the more fundamental expression w ∝ e−1U1 /kB T wL has wL = 1/107.25 for 1ML = 12, so that the inversion rate of particle 1 in Fig. 1(b) is two orders of magnitude slower. This slow inversion rate can be attributed to the highly aligned region having less entropy available to invert spins, similar to how a low-temperature heat bath has less entropy available to increase energy. Alternatively, if total entropy is to remain maximized during equilibrium fluctuations, the lower-entropy aligned region must be balanced by a higher-entropy heat bath, due to heat transferred to the bath from work done by entropic forces in reaching the low-entropy alignment. This increased entropy of the heat bath is similar to, but in addition to any increased entropy of the bath from work done by interaction forces in reaching lowenergy states. In any case, because the inversion rate governs lifetimes, the configuration in Fig. 1(b) is 107.25 times more likely than that predicted by the Boltzmann factor alone; the Boltzmann factor neglects configurational entropy from the alignment of the local region. Fig. 1(c) shows three 2 × 2 regions with various configurations of up and down particles. The binomial coefficient gives the statistical probability of the alignments, ranging from p4 = 1/16 for the all-up state, to p2 = 6/16 for no net alignment. If SL = −kB pL ln(pL ) is the entropy of each alignment, then this entropy will rise and fall during equilibrium fluctuations. Although seeming to violate the second law of thermodynamics, there are several possible explanations. (1) The net alignment of each region is not measureable, so that the entropy must always include the sum over all alignments. (2) The second law may apply only to macroscopic systems [26,27], allowing violations even for isolated small regions. (3) Particle inversion rates are altered by nonlinear corrections to Boltzmann’s factor, equalizing the probabilities, so that local configurational entropies also obey the second law. Because nonlinear corrections are expected for most small systems, this third mechanism may exist regardless of other explanations. Details are as follows. The master equation for the time dependence of the all-up state is: dp4 /dt = p3 w3+ − p4 w4− . At high temperatures where Boltzmann’s factor is e−1U1 /kB T ≈ 1, the standard Metropolis inversion rates are w4− = (4/4) and w3+ = (1/4), found from the fact that every simulation step converts 4 → 3, while only one step in four converts 3 → 4. Now include the nonlinear correction factor (wL ), so that ⟨w4− ⟩ = (4/4) ∗ 2!2!/(4!0!) = 1/6 and ⟨w3+ ⟩ = (1/4) ∗ 2!2!/(3!1!) = 1/6. Using detailed balance, the ratio of the thermal-average probabilities of each alignment is ⟨p4 ⟩/⟨p3 ⟩ = ⟨w3+ ⟩/⟨w4− ⟩ = 1. Thus, the probability of the all-up state now equals the probability of the three-up configurations, not because there are additional all-up states, but because there is additional heat transferred to the thermal bath to reach the low-entropy configuration, analogous to how low-energy interactions are favored by the Boltzmann factor at low temperatures. Similarly, for the other two configurations in Fig. 1(c): ⟨p3 ⟩/⟨p2 ⟩ = ⟨w2+ ⟩/⟨w3− ⟩ = [(2/4) ∗ 2!2!/(2!2!)]/[(3/4) ∗ 2!2!/(3!1!)] = 1. Each
5388
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
net alignment is now equally likely, allowing the alignment entropy-kB ⟨pL ⟩ ln⟨pL ⟩ to remain constant during equilibrium fluctuations, so that the second law of thermodynamics may also apply to an isolated small system. One interpretation is that the nonlinear correction is a simplistic way to allow semi-classical models to simulate indistinguishable particles, again reminiscent of the standard explanation of the usual Gibbs’ paradox. Specifically, if all the particles in a region are described by a single wave function, then each region would contain a superposition of identical particles, so that the statistics would depend only on the net alignment of the region, not on the location of each particle. Thus, the region would define a fundamental limit of experimental resolution, so that individual particles could not be identified inside the region. Another explanation is if particles change positions rapidly within each region, so that slower dynamics occurs between regions having time-averaged behavior. This separation of time scales could come from conservation laws, such as conservation of angular momentum during spin exchange [28], or conservation of particle density during exchange of local structures. Then the loss of correlation between regions would come from time averaging the interactions across the interfaces between independent regions. In any case, regions that contain an ensemble of timeaveraged states may provide the coarse-graining necessary for entropy increase as equilibrium is approached [29]. 3. Monte Carlo simulations I investigate nonlinear corrections to the Boltzmann factor using Monte Carlo simulations of the three-dimensional Ising model, on a simple-cubic lattice with periodic boundary conditions, containing a total of N = (24)3 or (48)3 particles (‘‘spins’’). The spin at site i may be up (σi = +1) or down (σi = −1), which simulates the two states of a uniaxial magnet. Other systems that map to the Ising model include lattice gases and binary alloys. For the Jahn–Teller distortion that I focus on here, the crystal unit cell may elongate in the x- or y-direction, with alternating orientations at low temperatures N corresponding to the AF Ising model. The interaction energy is U = (1/2)J i=1 σi Hi , where the exchange constant is J > 0 for AF interactions. The internal field at each site is Hi = Σ⟨ij⟩ σj , with the sum over all six nearest neighbors in the simplecubic lattice. The specific heat is ⟨cV ⟩ = (⟨U 2 ⟩ − ⟨U ⟩2 )/(NkB T 2 ). The maximum in ⟨cV ⟩ gives Tc . Thermal averages are found by averaging over the final 30% of 105 –106 Monte Carlo sweeps. Finite-size effects are studied by subdividing the system into cubic regions containing NL spins, ranging in size from N NL = (3)3 to (16)3 . The interaction energy of each region is UL = (1/2)J i=L 1 σi Hi . The sum is over all spins in the region, with Hi from all six nearest-neighbor spins, even if the neighbor is across an interface in another region. Thus, this subdivision does not change the interaction Hamiltonian; the nonlinear correction simply slows down the dynamics of regions that have fluctuated into low-entropy states. For AF behavior, the local order parameter comes from the sublattice alignment, NL x+y+z ML = σi , where x, y, and z are the integers that define each lattice site in Cartesian coordinates. The offset i=1 (−1) in the local order parameter is 1ML = ML − ML , where ML is the instantaneous average over all regions in the sample. Although the binomial coefficient gives the exact values for wL , for computational efficiency I use a quadratic (Gaussian) approximation that is accurate to third order. Specifically, the particle inverts if:
−1U1
> rand[0, 1) AND kB T (1ML )2 exp −g (1 − δ1U1 0 ) > rand[0, 1)
exp
(2)
2NL
where each rand[0, 1) is a pseudo-random number between 0 and 1. Here g is a constraint parameter that controls the correction factor, e.g. g = 0 gives the standard Metropolis algorithm, while g = 1 gives the expected Gaussian correction. The Kronecker delta (δ1U1 0 ) cancels the correction when inverting the spin does not change the interaction energy. Although this δ1U1 0 alters the equilibrium properties somewhat, it provides a much faster pathway towards equilibrium for the Ising model. Classical models with continuous variables should not need this factor to facilitate the dynamics. Furthermore, I argue that if 1U1 = 0, so that the spin has no net interaction with its neighbors, then E = ∂ U /∂η = 0 should give no correction. In other words, non-interacting particles are adequately described by the Boltzmann factor alone. 4. Results and discussion The main part of Fig. 2 is a comparison between measurements and simulations of the specific heat as a function of normalized temperature. The solid line is from data for LaMnO3 , after removing the smooth (non-critical) background [30]. The dashed line comes from simulations of the AF Ising model using the Metropolis algorithm, g = 0 in Eq. (2). The square symbols are from simulations of the same model with the Gaussian correction factor, g = 1. For these simulations, the region size is fixed at 3 × 3 × 3 (NL = 27), with no other adjustable parameters. The inset of Fig. 2 shows that this small region size significantly reduces the interaction energy at all temperatures, with a phase transition at kB Tc /J = 6.59 ± 0.03, 46% above the usual critical point of kB Tc /J = 4.51. Also note the abruptness of the transition, consistent with the sharp maximum in ⟨cV ⟩ shown in the main part of Fig. 2. The Ornstein–Zernike picture yields a pair-correlation function that decreases uniformly with increasing distance between particles [31], consistent with most MC simulations, but inconsistent with experiments showing dynamical
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
5389
Fig. 2. (Color online). Main plot shows specific heat as a function of normalized temperature. Solid line is from measurements of LaMnO3 [30] with the non-critical background removed. Dashed line is from simulations of the AF Ising model using the standard Metropolis algorithm (g = 0 in Eq. (2)). Solid squares (with error bars) are from the same model, but with a Gaussian correction factor (g = 1) from the local entropy of cubic (3 × 3 × 3) regions of size NL = 27. Inset shows the average energy per particle as a function of temperature. Dashed line is from the standard Metropolis algorithm (g = 0). Symbols include the Gaussian correction factor (g = 1) with the local entropy of cubic regions of size NL = 163 - J, 123 -, 83 -H, 63 -N, 43 -•, and 33 -. Because of negligible correction for large regions, some sets of symbols lie underneath the dashed line.
Fig. 3. (Color online). Semi-log plot of AF spin–spin correlation function versus distance from the central site. Symbols are from simulations of the AF Ising model at eight temperatures (given in the legend), which include a Gaussian correction (g = 1 in Eq. (2)) from 3 × 3 × 3 regions. At high temperatures, the spin–spin correlation function decreases exponentially over three lattice sites, then drops sharply to the next three sites. Solid lines are from linear fits to the data over the nearest-neighbor regions between r = 3 and 5 lattice parameters, giving the correlation length within regions λ = −1/[slope ∗ ln(10)]. Note the abrupt change in λ around kB T /J = 6.6, marking Tc .
heterogeneity. For the AF Ising model that is subdivided into cubic regions, I obtain a spin–spin correlation function from the alignment of each spin in the plane at x = 0, multiplied by the alignment of each spin at distance r in the x-direction, with an alternating sign for the sublattice: σ (0)(−1)r σ (r ) = Σy,z [σ (0, y, z )(−1)r σ (r , y, z )]/Σy,z 1. Fig. 3 shows this spin–spin correlation function versus distance, from simulations of the AF Ising model with the Gaussian correction (g = 1) from 3×3×3 cubic regions. The influence of the regions is most obvious at high temperatures, where the average ⟨σ (0)(−1)r σ (r )⟩ decreases smoothly over three lattice sites, then the correlation drops abruptly to the next three sites. Close to Tc there is no such drop in ⟨σ (0)(−1)r σ (r )⟩ between regions, allowing long-range correlations near the transition. Solid lines are from linear fits to the data over the nearest-neighbor regions between r = 3 and 5 unit cells, yielding an exponential decrease in the average correlation ⟨σ (0)(−1)r σ (r )⟩ ∝ e−r /λ , where λ is the correlation length within each region. The main part of Fig. 4 shows 1/λ as a function of normalized temperature. Circles are from neutron-scattering measurements of LaMnO3 [32]. Squares are from fits to simulations of the AF Ising model over the nearest-neighbor regions, as shown in Fig. 3. Conversion to physical units is made using the measured value of 0.40 nm for the average lattice parameter of LaMnO3 , with no adjustable parameters. Note that because of the enhanced correlation within each region, and loss of correlation between regions, this 1/λ is significantly smaller than the inverse correlation length from the homogeneous model (g = 0) at the same reduced temperature, and therefore closer to the data. The upper part of the inset shows two very similar pair distribution functions (pdfs) obtained from neutron scattering measurements at two temperatures, 1.2 Tc
5390
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
Fig. 4. (Color online). Main plot shows inverse correlation length as a function of normalized temperature. Circles are from neutron-scattering measurements of LaMnO3 [32]. Squares (with error bars) are from simulations of the AF Ising model including a Gaussian correction (g = 1 in Eq. (2)) from 3 × 3 × 3 regions, as found from the solid lines in Fig. 3. The upper part of the inset shows two nearly overlapping pair-distribution functions obtained from neutron scattering measurements at two temperatures, 1.2 Tc and 0.9 Tc . (Note that a common functional dependence has been removed from both pdfs.) The lower part of the inset shows the difference between these two pdfs. A sharp increase in this difference at distances beyond ∼0.9 nm is consistent with the sharp drop-off in correlation across the interface between regions shown in Fig. 3.
and 0.9 Tc . The lower part of the inset shows the difference between these two pdfs. A sharp rise in this difference at distances beyond ∼0.9 nm differs from the smooth loss in correlation of the Ornstein–Zernike picture, but is consistent with the sharp drop in correlation across the interface between regions shown in the Monte Carlo simulations of Fig. 3. 5. Conclusions I describe a nonlinear correction to the Metropolis algorithm that allows Monte Carlo simulations to show nonextensive energy instead of non-extensive entropy. Although neither scenario occurs in the usual thermodynamic limit, for independently relaxing regions in heterogeneous materials I present several arguments favoring the former: nonextensive energy is expected from interface terms and fluctuation effects, total entropy is maximized by including the region and its thermal bath, total energy is conserved by adding workdone by entropic forces, and the second law of thermodynamics is obeyed during normal thermal fluctuations. Then I show that the nonlinear correction improves agreement between simulations of the AF Ising model and measurements of heat capacity and structural correlations in LaMnO3 . More sophisticated models may further improve quantitative agreement with data; but if total energy is to be conserved, simulations of these models should also include the emergent energy from entropic forces. Moreover, the qualitative feature of thermodynamic heterogeneity is difficult to obtain from the usual Boltzmann factor that assumes homogenous thermal contact to an infinite heat bath. Acknowledgments I thank S.J.L. Billinge for providing me with the neutron scattering data, and for his helpful comments. I thank G.H. Wolf for several detailed discussions. I also thank R.A. Anthenien, N. Bernhoeft, K. Ghosh, R. Richert, T.D. Sewell, R.H. Swendsen, and D.L. Thompson for additional comments. I thank the Army Research Office for funding this research, and the ASU Advanced Computing Center for technical support. References [1] K. Binder, Finite size scaling analysis of Ising model block distribution functions, Z. Phys. B 43 (1981) 119–140. [2] M. Rovere, D.W. Hermann, K. Binder, Block density distribution function analysis of two-dimensional Lennard-Jones fluids, Europhys. Lett. 6 (1988) 585–590. [3] D. Hestenes, Entropy and indistinguishability, Amer. J. Phys. 38 (1970) 840–845. [4] R.H. Swendsen, Statistical mechanics of colloids and Boltzmann’s definition of entropy, Amer. J. Phys. 74 (2006) 187–190. [5] R.V. Chamberlin, Experiments and theory of the nonexponential relaxation in liquids, glasses, polymers and crystals, Phase Transit. 65 (1998) 169–209. [6] R. Böhmer, R.V. Chamberlin, G. Diezemann, B. Geil, A. Heuer, G. Hinze, S.C. Kuebler, R. Richert, B. Schiener, H. Sillescu, H.W. Spiess, U. Tracht, M. Wilhelm, Nature of the non-exponential primary relaxation in structural glass-formers probed by dynamically selective experiments, J. Non-Cryst. Solids 235–237 (1998) 1–9. [7] M.D. Ediger, Spatially heterogeneous dynamics in supercooled liquids, Annu. Rev. Phys. Chem. 51 (2000) 99–128. [8] R. Richert, Heterogeneous dynamics in liquids: fluctuations in space and time, J. Phys.: Condens. Matter 14 (2002) R703–R738. [9] R. Balian, From Microphysics to Macrophysics, Volume 1, Springer-Verlag, Berlin, 1992. Translated by D. ter Haar and J.F. Gregg (Chapter 3). [10] C. Tsallis, What should a statistical mechanics satisfy to reflect nature? Physica D 193 (2004) 3–34. [11] B. Schiener, R. Böhmer, A. Loidl, R.V. Chamberlin, Nonresonant spectral hole burning in the slow dielectric response of supercooled liquids, Science 274 (1996) 752–754.
R.V. Chamberlin / Physica A 391 (2012) 5384–5391 [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]
5391
K. Duvvuri, R. Richert, Dielectric hole burning in the high frequency wing of supercooled glycerol, J. Chem. Phys. 118 (2003) 1356–1363. M.J. Klein, L. Tisza, Theory of critical fluctuations, Phys. Rev. 76 (1949) 1861–1868. E. Verlinde, On the origin of gravity and the laws of Newton, J. High Energy Phys. 4 (2011) 029. L.D. Landau, E.M. Lifshitz, Statistical Physics, third ed., in: Course of Theoretical Physics, vol. 5, Pergamon Press, Oxford, 1980 (Chapter XII). H.B. Callen, Thermodynamics and An Introduction to Thermostatistics, second ed., Wiley, New York, 1985 (Section 15.4 and Chapter 19). R. Dickman, P. Attard, V. Simonian, Entropic forces in binary liquids, J. Chem. Phys. 107 (1997) 205–213. I.M. Sokolov, Statistical mechanics of entropic forces: disassembling a toy, European J. Phys. 31 (2010) 1353–1367. R.V. Chamberlin, G.H. Wolf, Fluctuation theory constraint for extensive entropy in Monte Carlo simulations, Eur. Phys. J. B 67 (2009) 495–499. R.V. Chamberlin, J.V. Vermaas, G.H. Wolf, Beyond the Boltzmann factor for corrections to scaling in ferromagnetic materials and critical fluids, Eur. Phys. J. B 71 (2009) 1–6. T.L. Hill, R.V. Chamberlin, Fluctuations in energy in completely open small systems, Nano Lett. 2 (2002) 609–613. T.L. Hill, Thermodynamics of Small Systems, Dover, New York, 1994. R.V. Chamberlin, Mesoscopic mean-field theory for supercooled liquids and the glass transition, Phys. Rev. Lett. 82 (1999) 2520–2523. R.V. Chamberlin, Critical behavior from Landau theory in nanothermodynamic equilibrium, Phys. Lett. A 315 (2003) 313–318. M.R.H. Javaheri, R.V. Chamberlin, A free-energy landscape picture and Landau theory for the dynamics of disordered materials, J. Chem. Phys. 125 (2006) 154503. D.J. Evans, E.G.D. Cohen, G.P. Morriss, Probability of second law violations in shearing steady states, Phys. Rev. Lett. 72 (1993) 2401–2404. G.M. Wang, E.M. Sevick, E. Mittag, D.J. Searles, D.J. Evans, Experimental demonstration of violations of the second law of thermodynamics for small systems and short time scales, Phys. Rev. Lett. 89 (2002) 050601. R.V. Chamberlin, K.J. Stangel, Monte Carlo simulation of supercooled liquids using a self-consistent local temperature, Phys. Lett. A 350 (2006) 400–404. A. Wehrl, General properties of entropy, Rev. Modern Phys. 50 (1978) 221–260. H. Satoh, M. Takagi, K. Kinukawa, N. Kamegashira, Heat capacity of LaMnO3 , Therm. Acta 299 (1997) 123–126. C. Domb, Critical phenomena: a brief historical survey, Contemp. Phys. 26 (1985) 49–72. X. Qiu, T. Proffen, J.F. Mitchell, S.J.L. Billinge, Orbital correlations in the pseudocubic O and rhombohedral R phases of LaMnO3 , Phys. Rev. Lett. 94 (2005) 177203.
Contents lists available at SciVerse ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Monte Carlo simulations including energy from an entropic force Ralph V. Chamberlin ∗ Department of Physics, Arizona State University, Tempe, AZ 85287-1504, USA
article
info
Article history: Received 8 August 2011 Available online 27 June 2012 Keywords: Entropic force Heterogeneous dynamics Ising model Monte Carlo simulations Small-system thermodynamics
abstract Several experimental techniques have shown that the primary response of many materials comes from a heterogeneous distribution of independently relaxing nanoscale regions; but most Monte Carlo simulations have homogeneous correlations. Resolving this discrepancy may require including the energy needed to change the configurational entropy, which is often used in theoretical treatments of thermal fluctuations, but not in computer simulations. Here the local configurational entropy is shown to give a nonlinear correction to the Metropolis algorithm that restores conservation of energy, maintains maximum entropy, and yields heterogeneous correlations. The nonlinear correction also improves agreement between Monte Carlo simulations of the Ising model and measurements of specific heat and structural correlations from the Jahn–Teller distortion in LaMnO3 . © 2012 Elsevier B.V. All rights reserved.
1. Introduction In most computer simulations of classical models, the energy fluctuations in each local subsystem, and resulting local entropy (SL ), fail to increase linearly with increasing subsystem size, NL [1,2]. Although superficially similar to the usual Gibbs’ paradox for non-extensive entropy in large systems of distinguishable particles [3,4], the non-extensive energy fluctuations in most computer simulations arise from homogeneous correlations. However, several experimental techniques have shown that the primary response of many materials comes from a heterogeneous distribution of independently relaxing nanometer-sized ‘‘regions’’, with negligible correlations between the regions [5–8]. These experiments often show that the heterogeneity is dynamical: occurring transiently in the thermal-equilibrium fluctuations of single-component substances having homogeneous average density and uniform interactions. For most computer simulations of independent regions in these materials, either the dynamics needs a correction to make the entropy additive, or the foundation of statistical mechanics needs to be modified to accommodate non-extensive entropy [9,10]. Key insight comes from a quantitative analysis of non-resonant spectral hole burning measurements showing that the heat capacity of these regions is extensive [11,12]. Thus, simulations exhibiting non-extensive energy fluctuations and uniform correlations cannot explain the measured response. The cellular method [13] is a theoretical tool for treating fluctuations by subdividing a bulk sample into small subsystems, so that each subsystem couples to a self-consistent ensemble of similar subsystems. The usual procedure is to include uniform coupling between subsystems to maintain homogeneous correlations, unlike the measurements showing heterogeneous dynamics. Here I explore a type of dynamics that includes a local configurational entropy, giving a nonlinear correction to the Boltzmann factor that decreases the coupling between subsystems, allowing otherwise homogeneous systems to exhibit response from an ensemble of effectively independent regions. The energy from this entropy is emergent [14], hence it does not appear in the interaction Hamiltonian, yet it is needed maintain conservation of energy from nanometer-sized fluctuations. Because the correction is dynamical, the system will eventually become
∗
Tel.: +1 480 965 3922; fax: +1 480 965 7954. E-mail address: [email protected].
0378-4371/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2012.06.016
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
5385
homogeneous, but at longer times and lower energies than the uncorrected system due to the heterogeneous response at intermediate times. Entropic forces are often used in the theory of Gaussian fluctuations [15,16], but most computer simulations focus on interparticle interactions and statistical probabilities, not statistical forces. Even when entropic forces are simulated [17,18], finite-size effects from the local entropy are usually neglected. Recently, a configurational entropy from the local order parameter ML has been shown to give a Gaussian correction to the Metropolis algorithm that lowers the energy, makes the entropy extensive, and improves agreement between computer simulations and measurements of ferromagnetic materials and critical fluids [19,20]. Here I add a similar Gaussian correction to Monte Carlo simulations of the antiferromagnetic (AF) Ising model. I show that away from the critical temperature Tc , the correction reduces the correlations between regions, yielding heterogeneous dynamics in the simulation. Moreover, the correction improves agreement between the measured specific heat and structural correlations found directly from the Jahn–Teller distortion in LaMnO3 . 2. Background 2.1. Thermal properties of small regions Let NL be the number of particles in a subsystem having interaction energy UL . In most computer simulations, local energy fluctuations are non-extensive, i.e. the thermal average ⟨UL2 ⟩ − ⟨UL ⟩2 increases nonlinearly with increasing NL [1,2]. Thus, the resulting heat capacity CV = (⟨UL2 ⟩ − ⟨UL ⟩2 )/(kB T 2 ) [21] and local entropy difference between two temperatures
T
SL (T2 ) − SL (T1 ) = T 2 (CV /T ) dT , are also non-extensive. Consequently, the entropy of a combined subsystem a + b does 1 not equal the entropy of its constituents: Sa+b ̸= Sa + Sb . For homogeneous systems, this non-additive entropy is expected from uniform correlations [9]. Specifically, using the Boltzmann–Gibbs expression for entropy S = −kB Σi pi ln(pi ) with kB Boltzmann’s constant and pi the probability of each state, correlations prohibit separating the probability into a product: pa+b ̸= pa ∗ pb . However, the primary response of many materials is now known to come from independently relaxing nanometer-sized regions that are effectively uncorrelated [5–8]. For this response, to a good approximation pa+b = pa ∗ pb , so that Sa+b = Sa + Sb , unlike most simulations. Fundamental transition rates in thermostatistics are given by w ∝ e1St /kB , where 1St is the change in total entropy. First consider two independent contributions to entropy, one from the change of a single particle (1S1 ) and one from an infinite bath (1S ∗ ), so that 1St = 1S1 + 1S ∗ . The Boltzmann factor is obtained by assuming that the single particle changes from a known initial state to a known final state, so that 1S1 = 0. Then the total entropy change is assumed to come only from the change in energy of the infinite bath 1U ∗ , so that using the fundamental equation of thermodynamics (combined first and second laws) gives 1St = 1U ∗ /T . The minus sign in the Boltzmann factor comes from replacing 1U ∗ by the change in energy of the single particle (1U1 ) using conservation of energy 1U ∗ = −1U1 , giving w ∝ e−1U1 /kB T . Strictly speaking, the standard Boltzmann factor assumes instantaneous thermal contact to an infinite heat bath. Here I consider some consequences of treating an independently relaxing region as a small local bath, and including configurational entropy as another source of heat. A more complete expression for the change in total entropy is 1St = 1S1 + 1S ∗ + 1SL [20]. For a single particle changing between known states, 1S1 = 0 as before. For an effectively infinite heat bath, 1S ∗ = −1U1 /T as before. The new term 1SL is the local offset in entropy from the region containing the single particle. Note that adding 1SL to 1S ∗ is justified because the region relaxes independently, even if the rest of the sample forms a heat bath for the region. For large regions, 1SL may include entropy from particles that do not interact directly with the single particle; whereas for small regions, 1U1 often includes interactions that extend into neighboring regions. In any case, the change in energy rarely defines the total change in entropy; many different configurations in a region may yield the same 1U1 . In other words, it is unlikely that a single scalar quantity will define the locations and alignments of all particles in a region. Thus, if the statistical probabilities are to be based on the total entropy: a local configurational entropy must be added to the entropy of the infinite bath, giving a nonlinear correction to the Boltzmann factor. Another way to think of this correction is to include conservation of energy from work done on entropic forces. Entropic forces generally favor maximum entropy, such as the unaligned state of independent dipoles, or the uniform density of ideal gases. Even for ideal systems with no interactions, nonzero work is needed to displace the system from its maximum entropy. One example is the entropic restoring force that opposes stretching a polymer [16]. Another example is Fick’s law, which gives a restoring force proportional to the gradient in the particle density, so that the restoring force increases linearly with how far the particles have moved from a uniform distribution. The resulting emergent energy depends quadratically on the average offset of particles from equilibrium, even if there is no interaction between particles. Quadratic energies from entropic forces often form the basis of theoretical treatments of thermal fluctuations [15,16], but they are usually neglected in computer simulations. For example, the Metropolis algorithm weights the probability of configurations using only the interaction energy, whereas including the local configurational entropy can shift critical temperatures by nearly 50% ([20] and shown below). In systems of independently relaxing regions, a nonlinear correction to the Boltzmann factor is necessary if energy is to be conserved and total entropy is to be maximized during normal thermal fluctuations.
5386
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
2.2. Nanothermodynamics Although small-system thermodynamics [22] was originally developed to treat isolated small systems, similar nonextensive contributions to energy occur in ensembles of independently fluctuating regions. In this ‘‘nanothermodynamics’’, energy conservation is restored by adding a new pair of conjugate variables to the fundamental equation of thermodynamics. Thus, for magnetic systems:
1U = T 1S + H 1M + µ1N + E1η.
(1)
Here η ∼ N /NL is the total number of regions, E = ∂ U /∂η is the subdivision potential, and the other symbols have the usual meaning (T —temperature, S—entropy, H—magnetic field, and M—magnetic moment). This E can be understood by comparison to the chemical potential (µ): µ is the change in energy to add a single particle to the system, whereas E is the extra energy from adding a region of interacting particles to the system; and in general a region of NL interacting particles does not have the same energy as NL isolated particles. In other words, E contains all of the non-extensive and emergent energy terms that occur in finite-sized systems of interacting particles. Classical thermodynamics is based on the assumption that both energy and entropy are extensive, but finite-size effects alter the behavior of small systems. Two standard sources of non-extensive energy for internal regions are: interface terms (breaking correlations between regions costs energy proportional to the area of internal surfaces) and fluctuation effects. In general, the average energy of an independently fluctuating region is lowered inversely proportional to NL , which can be understood from mesoscopic mean-field theory [23–25] as follows. In thermal equilibrium at high temperatures, the meanfield energy of an infinite system is at its maximum because the entropy is also maximized. Subdividing the system into an ensemble of independent regions lowers the mean-field energy due to fluctuations. Mean-square fluctuations of intensive variables vary inversely proportional to the number of particles, so that the mean-field energy per particle is reduced inversely proportional to NL . Thermal equilibrium is found by balancing this energy reduction from increased fluctuations with the entropy reduction from finite-sized regions, thus minimizing the free energy of the generalized ensemble (µ, H , T ). Note that treating independent regions inside a bulk system using a canonical or grand-canonical ensemble artificially fixes the volume of each region. Also note that although this mean-field energy reduction is in the interaction, it comes from an entropic force that allows regions to fluctuate independently, so that the energy reduction emerges only for finite-sized regions. A similar 1/NL reduction in the average interaction energy is seen in computer simulations of the Ising model using the nonlinear correction to the Boltzmann factor that allows the regions to fluctuate independently [19]. It is often more difficult to justify non-extensive entropy. Furthermore, because the nonlinear correction modifies the dynamics, not the Hamiltonian, the instantaneous interaction energy remains additive at every step. Boltzmann’s definition of entropy involves a sum over all microstates, not a weighted average, so that the entropy of uncorrelated regions should be additive at all times. 2.3. Beyond the Boltzmann factor Now return to the fundamental expression for transition rates w ∝ e1St /kB . Here 1St = −1U1 /T + 1SL , from Eq. (1) for a fixed number of particles (1N = 0) in zero field (H = 0). Recall that 1U1 is the energy borrowed from the infinite bath to balance the change in interaction energy during single-particle dynamics. Eq. (1) gives the connection between the emergent energy and the change in entropy for the small regions 1SL = −E1η/T . Although this 1SL may contain several factors, the main contribution in many systems comes from the quadratic term 1SL = 0.5(∂ 2 SL /∂ ML2 )(1ML )2 , where 1ML is the offset from equilibrium in the local order parameter. I emphasize that because thermal equilibrium is required for Eq. (1), this 1ML is the offset from equilibrium due to normal thermal fluctuations, not the difference between two out-of-equilibrium states. Note that in a large enough external field, H 1M may dominate conservation of energy in Eq. (1), but if H ≈ 0 the entropic force that yields 1SL is the main source of magnetic work. Also note that this magnetic work exists even if H = 0, because ∂ 2 SL /∂ ML2 < 0 is necessary for stable fluctuations. The usual Boltzmann factor neglects all such nonlinear terms. If nonlinear magnetic work is included, the fundamental rate for changing the energy in a local region of a bulk sample comes from the usual Boltzmann factor e−1U1 /kB T , times a nonlinear correction from the entropy of the region wL = e1SL /kB . For realistic systems, nonlinear interactions and kinetic energy may also influence this correction. For the Ising model, using the binomial coefficient for the multiplicity of local alignments, the full correction to all orders is wL = (0.5NL )!(0.5NL )!/{[0.5(NL + 1ML )]![0.5(NL − 1ML )]!}. This nonlinear correction can also be understood by maximizing the entropy of a system and its heat bath. Using Lagrange multipliers, α for normalized probabilities and 1/kB T for heat flow from the thermal bath, the Boltzmann–Gibbs entropy of a uniform system is maximized by maximizing −kB Σi pi [ln(pi ) + α + Ui /kB T ], yielding pi ∝ e−Ui /kB T for the probability of a state of energy Ui . Subdividing the system into an ensemble of independently fluctuating regions lowers its entropy, but its energy is also lowered inversely proportionally to NL . This extra energy added to the heat bath usually increases the combined entropy of the system and its bath, as shown by the reduction in free energy [19]. The Boltzmann factor and canonical ensemble neglect this contribution from local configurational entropy, yielding thermal equilibrium only if the system can be thought of as isolated states that couple individually to a large heat bath. The total entropy from local configurations of interacting particles plus heat bath can be increased by subdividing the system into a generalized ensemble of independently fluctuating regions, with transition rates that include a nonlinear correction factor w ∝ e−1U1 /kB T e1SL /kB .
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
5387
Fig. 1. Schematic representations of binary particles that may be up (white) or down (black). Region (a) has an equal number of up and down particles giving 1SL = 0. Region (b) has fluctuated into a highly aligned state with 1SL /kB = ln{(8!8!)/[14!2!]} = ln(1/107.25). Both regions have a particle (marked by 1) with identical interaction energies; but (b) has a nonlinear correction factor wL = e1SL /kB , causing its inversion rate to be 107.25 times slower than that predicted by the Boltzmann factor alone. (c) 2 × 2 regions containing, from left to right: four-up, three-up, and two-up particle alignments. The statistical probabilities of each alignment are: p4 = 1/16, p3 = 4/16, and p2 = 6/16. At high temperatures, the classical Metropolis inversion rates between each pair of alignments are w4+ = 4/4 and w3− = 1/4, with w3+ = 3/4 and w2− = 2/4. Including the nonlinear correction factors, which depend only on the offset of entropy from its maximum, the inversion rates become ⟨w4+ ⟩ = 1/6 and ⟨w3− ⟩ = 1/6, with ⟨w3+ ⟩ = 1/2 and ⟨w2− ⟩ = 1/2. Thus, detailed balance gives equal probabilities ⟨p4 ⟩ = ⟨p3 ⟩ = ⟨p2 ⟩, and the alignment entropy −kB ⟨pL ⟩ ln⟨pL ⟩ does not rise and fall during equilibrium fluctuations.
The correction factor is not used in Onsager’s solution of the two-dimensional Ising model, which assumes homogenous correlations in the canonical ensemble. If energy reduction from independent nanometer-sized fluctuations is included, the Onsager solution fails to give the thermal equilibrium of the Ising model. 2.4. Examples The sketches in Fig. 1(a) and (b) show a specific example that emphasizes the importance of the nonlinear correction. Note that both (a) and (b) have an up site (marked by 1) with all four nearest neighbors down. For nearest-neighbor interactions, inverting these sites yields the same energy change 1U1 , so that the Boltzmann factor alone predicts identical inversion rates for particle 1 in both (a) and (b). Whereas the more fundamental expression w ∝ e−1U1 /kB T wL has wL = 1/107.25 for 1ML = 12, so that the inversion rate of particle 1 in Fig. 1(b) is two orders of magnitude slower. This slow inversion rate can be attributed to the highly aligned region having less entropy available to invert spins, similar to how a low-temperature heat bath has less entropy available to increase energy. Alternatively, if total entropy is to remain maximized during equilibrium fluctuations, the lower-entropy aligned region must be balanced by a higher-entropy heat bath, due to heat transferred to the bath from work done by entropic forces in reaching the low-entropy alignment. This increased entropy of the heat bath is similar to, but in addition to any increased entropy of the bath from work done by interaction forces in reaching lowenergy states. In any case, because the inversion rate governs lifetimes, the configuration in Fig. 1(b) is 107.25 times more likely than that predicted by the Boltzmann factor alone; the Boltzmann factor neglects configurational entropy from the alignment of the local region. Fig. 1(c) shows three 2 × 2 regions with various configurations of up and down particles. The binomial coefficient gives the statistical probability of the alignments, ranging from p4 = 1/16 for the all-up state, to p2 = 6/16 for no net alignment. If SL = −kB pL ln(pL ) is the entropy of each alignment, then this entropy will rise and fall during equilibrium fluctuations. Although seeming to violate the second law of thermodynamics, there are several possible explanations. (1) The net alignment of each region is not measureable, so that the entropy must always include the sum over all alignments. (2) The second law may apply only to macroscopic systems [26,27], allowing violations even for isolated small regions. (3) Particle inversion rates are altered by nonlinear corrections to Boltzmann’s factor, equalizing the probabilities, so that local configurational entropies also obey the second law. Because nonlinear corrections are expected for most small systems, this third mechanism may exist regardless of other explanations. Details are as follows. The master equation for the time dependence of the all-up state is: dp4 /dt = p3 w3+ − p4 w4− . At high temperatures where Boltzmann’s factor is e−1U1 /kB T ≈ 1, the standard Metropolis inversion rates are w4− = (4/4) and w3+ = (1/4), found from the fact that every simulation step converts 4 → 3, while only one step in four converts 3 → 4. Now include the nonlinear correction factor (wL ), so that ⟨w4− ⟩ = (4/4) ∗ 2!2!/(4!0!) = 1/6 and ⟨w3+ ⟩ = (1/4) ∗ 2!2!/(3!1!) = 1/6. Using detailed balance, the ratio of the thermal-average probabilities of each alignment is ⟨p4 ⟩/⟨p3 ⟩ = ⟨w3+ ⟩/⟨w4− ⟩ = 1. Thus, the probability of the all-up state now equals the probability of the three-up configurations, not because there are additional all-up states, but because there is additional heat transferred to the thermal bath to reach the low-entropy configuration, analogous to how low-energy interactions are favored by the Boltzmann factor at low temperatures. Similarly, for the other two configurations in Fig. 1(c): ⟨p3 ⟩/⟨p2 ⟩ = ⟨w2+ ⟩/⟨w3− ⟩ = [(2/4) ∗ 2!2!/(2!2!)]/[(3/4) ∗ 2!2!/(3!1!)] = 1. Each
5388
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
net alignment is now equally likely, allowing the alignment entropy-kB ⟨pL ⟩ ln⟨pL ⟩ to remain constant during equilibrium fluctuations, so that the second law of thermodynamics may also apply to an isolated small system. One interpretation is that the nonlinear correction is a simplistic way to allow semi-classical models to simulate indistinguishable particles, again reminiscent of the standard explanation of the usual Gibbs’ paradox. Specifically, if all the particles in a region are described by a single wave function, then each region would contain a superposition of identical particles, so that the statistics would depend only on the net alignment of the region, not on the location of each particle. Thus, the region would define a fundamental limit of experimental resolution, so that individual particles could not be identified inside the region. Another explanation is if particles change positions rapidly within each region, so that slower dynamics occurs between regions having time-averaged behavior. This separation of time scales could come from conservation laws, such as conservation of angular momentum during spin exchange [28], or conservation of particle density during exchange of local structures. Then the loss of correlation between regions would come from time averaging the interactions across the interfaces between independent regions. In any case, regions that contain an ensemble of timeaveraged states may provide the coarse-graining necessary for entropy increase as equilibrium is approached [29]. 3. Monte Carlo simulations I investigate nonlinear corrections to the Boltzmann factor using Monte Carlo simulations of the three-dimensional Ising model, on a simple-cubic lattice with periodic boundary conditions, containing a total of N = (24)3 or (48)3 particles (‘‘spins’’). The spin at site i may be up (σi = +1) or down (σi = −1), which simulates the two states of a uniaxial magnet. Other systems that map to the Ising model include lattice gases and binary alloys. For the Jahn–Teller distortion that I focus on here, the crystal unit cell may elongate in the x- or y-direction, with alternating orientations at low temperatures N corresponding to the AF Ising model. The interaction energy is U = (1/2)J i=1 σi Hi , where the exchange constant is J > 0 for AF interactions. The internal field at each site is Hi = Σ⟨ij⟩ σj , with the sum over all six nearest neighbors in the simplecubic lattice. The specific heat is ⟨cV ⟩ = (⟨U 2 ⟩ − ⟨U ⟩2 )/(NkB T 2 ). The maximum in ⟨cV ⟩ gives Tc . Thermal averages are found by averaging over the final 30% of 105 –106 Monte Carlo sweeps. Finite-size effects are studied by subdividing the system into cubic regions containing NL spins, ranging in size from N NL = (3)3 to (16)3 . The interaction energy of each region is UL = (1/2)J i=L 1 σi Hi . The sum is over all spins in the region, with Hi from all six nearest-neighbor spins, even if the neighbor is across an interface in another region. Thus, this subdivision does not change the interaction Hamiltonian; the nonlinear correction simply slows down the dynamics of regions that have fluctuated into low-entropy states. For AF behavior, the local order parameter comes from the sublattice alignment, NL x+y+z ML = σi , where x, y, and z are the integers that define each lattice site in Cartesian coordinates. The offset i=1 (−1) in the local order parameter is 1ML = ML − ML , where ML is the instantaneous average over all regions in the sample. Although the binomial coefficient gives the exact values for wL , for computational efficiency I use a quadratic (Gaussian) approximation that is accurate to third order. Specifically, the particle inverts if:
−1U1
> rand[0, 1) AND kB T (1ML )2 exp −g (1 − δ1U1 0 ) > rand[0, 1)
exp
(2)
2NL
where each rand[0, 1) is a pseudo-random number between 0 and 1. Here g is a constraint parameter that controls the correction factor, e.g. g = 0 gives the standard Metropolis algorithm, while g = 1 gives the expected Gaussian correction. The Kronecker delta (δ1U1 0 ) cancels the correction when inverting the spin does not change the interaction energy. Although this δ1U1 0 alters the equilibrium properties somewhat, it provides a much faster pathway towards equilibrium for the Ising model. Classical models with continuous variables should not need this factor to facilitate the dynamics. Furthermore, I argue that if 1U1 = 0, so that the spin has no net interaction with its neighbors, then E = ∂ U /∂η = 0 should give no correction. In other words, non-interacting particles are adequately described by the Boltzmann factor alone. 4. Results and discussion The main part of Fig. 2 is a comparison between measurements and simulations of the specific heat as a function of normalized temperature. The solid line is from data for LaMnO3 , after removing the smooth (non-critical) background [30]. The dashed line comes from simulations of the AF Ising model using the Metropolis algorithm, g = 0 in Eq. (2). The square symbols are from simulations of the same model with the Gaussian correction factor, g = 1. For these simulations, the region size is fixed at 3 × 3 × 3 (NL = 27), with no other adjustable parameters. The inset of Fig. 2 shows that this small region size significantly reduces the interaction energy at all temperatures, with a phase transition at kB Tc /J = 6.59 ± 0.03, 46% above the usual critical point of kB Tc /J = 4.51. Also note the abruptness of the transition, consistent with the sharp maximum in ⟨cV ⟩ shown in the main part of Fig. 2. The Ornstein–Zernike picture yields a pair-correlation function that decreases uniformly with increasing distance between particles [31], consistent with most MC simulations, but inconsistent with experiments showing dynamical
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
5389
Fig. 2. (Color online). Main plot shows specific heat as a function of normalized temperature. Solid line is from measurements of LaMnO3 [30] with the non-critical background removed. Dashed line is from simulations of the AF Ising model using the standard Metropolis algorithm (g = 0 in Eq. (2)). Solid squares (with error bars) are from the same model, but with a Gaussian correction factor (g = 1) from the local entropy of cubic (3 × 3 × 3) regions of size NL = 27. Inset shows the average energy per particle as a function of temperature. Dashed line is from the standard Metropolis algorithm (g = 0). Symbols include the Gaussian correction factor (g = 1) with the local entropy of cubic regions of size NL = 163 - J, 123 -, 83 -H, 63 -N, 43 -•, and 33 -. Because of negligible correction for large regions, some sets of symbols lie underneath the dashed line.
Fig. 3. (Color online). Semi-log plot of AF spin–spin correlation function versus distance from the central site. Symbols are from simulations of the AF Ising model at eight temperatures (given in the legend), which include a Gaussian correction (g = 1 in Eq. (2)) from 3 × 3 × 3 regions. At high temperatures, the spin–spin correlation function decreases exponentially over three lattice sites, then drops sharply to the next three sites. Solid lines are from linear fits to the data over the nearest-neighbor regions between r = 3 and 5 lattice parameters, giving the correlation length within regions λ = −1/[slope ∗ ln(10)]. Note the abrupt change in λ around kB T /J = 6.6, marking Tc .
heterogeneity. For the AF Ising model that is subdivided into cubic regions, I obtain a spin–spin correlation function from the alignment of each spin in the plane at x = 0, multiplied by the alignment of each spin at distance r in the x-direction, with an alternating sign for the sublattice: σ (0)(−1)r σ (r ) = Σy,z [σ (0, y, z )(−1)r σ (r , y, z )]/Σy,z 1. Fig. 3 shows this spin–spin correlation function versus distance, from simulations of the AF Ising model with the Gaussian correction (g = 1) from 3×3×3 cubic regions. The influence of the regions is most obvious at high temperatures, where the average ⟨σ (0)(−1)r σ (r )⟩ decreases smoothly over three lattice sites, then the correlation drops abruptly to the next three sites. Close to Tc there is no such drop in ⟨σ (0)(−1)r σ (r )⟩ between regions, allowing long-range correlations near the transition. Solid lines are from linear fits to the data over the nearest-neighbor regions between r = 3 and 5 unit cells, yielding an exponential decrease in the average correlation ⟨σ (0)(−1)r σ (r )⟩ ∝ e−r /λ , where λ is the correlation length within each region. The main part of Fig. 4 shows 1/λ as a function of normalized temperature. Circles are from neutron-scattering measurements of LaMnO3 [32]. Squares are from fits to simulations of the AF Ising model over the nearest-neighbor regions, as shown in Fig. 3. Conversion to physical units is made using the measured value of 0.40 nm for the average lattice parameter of LaMnO3 , with no adjustable parameters. Note that because of the enhanced correlation within each region, and loss of correlation between regions, this 1/λ is significantly smaller than the inverse correlation length from the homogeneous model (g = 0) at the same reduced temperature, and therefore closer to the data. The upper part of the inset shows two very similar pair distribution functions (pdfs) obtained from neutron scattering measurements at two temperatures, 1.2 Tc
5390
R.V. Chamberlin / Physica A 391 (2012) 5384–5391
Fig. 4. (Color online). Main plot shows inverse correlation length as a function of normalized temperature. Circles are from neutron-scattering measurements of LaMnO3 [32]. Squares (with error bars) are from simulations of the AF Ising model including a Gaussian correction (g = 1 in Eq. (2)) from 3 × 3 × 3 regions, as found from the solid lines in Fig. 3. The upper part of the inset shows two nearly overlapping pair-distribution functions obtained from neutron scattering measurements at two temperatures, 1.2 Tc and 0.9 Tc . (Note that a common functional dependence has been removed from both pdfs.) The lower part of the inset shows the difference between these two pdfs. A sharp increase in this difference at distances beyond ∼0.9 nm is consistent with the sharp drop-off in correlation across the interface between regions shown in Fig. 3.
and 0.9 Tc . The lower part of the inset shows the difference between these two pdfs. A sharp rise in this difference at distances beyond ∼0.9 nm differs from the smooth loss in correlation of the Ornstein–Zernike picture, but is consistent with the sharp drop in correlation across the interface between regions shown in the Monte Carlo simulations of Fig. 3. 5. Conclusions I describe a nonlinear correction to the Metropolis algorithm that allows Monte Carlo simulations to show nonextensive energy instead of non-extensive entropy. Although neither scenario occurs in the usual thermodynamic limit, for independently relaxing regions in heterogeneous materials I present several arguments favoring the former: nonextensive energy is expected from interface terms and fluctuation effects, total entropy is maximized by including the region and its thermal bath, total energy is conserved by adding workdone by entropic forces, and the second law of thermodynamics is obeyed during normal thermal fluctuations. Then I show that the nonlinear correction improves agreement between simulations of the AF Ising model and measurements of heat capacity and structural correlations in LaMnO3 . More sophisticated models may further improve quantitative agreement with data; but if total energy is to be conserved, simulations of these models should also include the emergent energy from entropic forces. Moreover, the qualitative feature of thermodynamic heterogeneity is difficult to obtain from the usual Boltzmann factor that assumes homogenous thermal contact to an infinite heat bath. Acknowledgments I thank S.J.L. Billinge for providing me with the neutron scattering data, and for his helpful comments. I thank G.H. Wolf for several detailed discussions. I also thank R.A. Anthenien, N. Bernhoeft, K. Ghosh, R. Richert, T.D. Sewell, R.H. Swendsen, and D.L. Thompson for additional comments. I thank the Army Research Office for funding this research, and the ASU Advanced Computing Center for technical support. References [1] K. Binder, Finite size scaling analysis of Ising model block distribution functions, Z. Phys. B 43 (1981) 119–140. [2] M. Rovere, D.W. Hermann, K. Binder, Block density distribution function analysis of two-dimensional Lennard-Jones fluids, Europhys. Lett. 6 (1988) 585–590. [3] D. Hestenes, Entropy and indistinguishability, Amer. J. Phys. 38 (1970) 840–845. [4] R.H. Swendsen, Statistical mechanics of colloids and Boltzmann’s definition of entropy, Amer. J. Phys. 74 (2006) 187–190. [5] R.V. Chamberlin, Experiments and theory of the nonexponential relaxation in liquids, glasses, polymers and crystals, Phase Transit. 65 (1998) 169–209. [6] R. Böhmer, R.V. Chamberlin, G. Diezemann, B. Geil, A. Heuer, G. Hinze, S.C. Kuebler, R. Richert, B. Schiener, H. Sillescu, H.W. Spiess, U. Tracht, M. Wilhelm, Nature of the non-exponential primary relaxation in structural glass-formers probed by dynamically selective experiments, J. Non-Cryst. Solids 235–237 (1998) 1–9. [7] M.D. Ediger, Spatially heterogeneous dynamics in supercooled liquids, Annu. Rev. Phys. Chem. 51 (2000) 99–128. [8] R. Richert, Heterogeneous dynamics in liquids: fluctuations in space and time, J. Phys.: Condens. Matter 14 (2002) R703–R738. [9] R. Balian, From Microphysics to Macrophysics, Volume 1, Springer-Verlag, Berlin, 1992. Translated by D. ter Haar and J.F. Gregg (Chapter 3). [10] C. Tsallis, What should a statistical mechanics satisfy to reflect nature? Physica D 193 (2004) 3–34. [11] B. Schiener, R. Böhmer, A. Loidl, R.V. Chamberlin, Nonresonant spectral hole burning in the slow dielectric response of supercooled liquids, Science 274 (1996) 752–754.
R.V. Chamberlin / Physica A 391 (2012) 5384–5391 [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]
5391
K. Duvvuri, R. Richert, Dielectric hole burning in the high frequency wing of supercooled glycerol, J. Chem. Phys. 118 (2003) 1356–1363. M.J. Klein, L. Tisza, Theory of critical fluctuations, Phys. Rev. 76 (1949) 1861–1868. E. Verlinde, On the origin of gravity and the laws of Newton, J. High Energy Phys. 4 (2011) 029. L.D. Landau, E.M. Lifshitz, Statistical Physics, third ed., in: Course of Theoretical Physics, vol. 5, Pergamon Press, Oxford, 1980 (Chapter XII). H.B. Callen, Thermodynamics and An Introduction to Thermostatistics, second ed., Wiley, New York, 1985 (Section 15.4 and Chapter 19). R. Dickman, P. Attard, V. Simonian, Entropic forces in binary liquids, J. Chem. Phys. 107 (1997) 205–213. I.M. Sokolov, Statistical mechanics of entropic forces: disassembling a toy, European J. Phys. 31 (2010) 1353–1367. R.V. Chamberlin, G.H. Wolf, Fluctuation theory constraint for extensive entropy in Monte Carlo simulations, Eur. Phys. J. B 67 (2009) 495–499. R.V. Chamberlin, J.V. Vermaas, G.H. Wolf, Beyond the Boltzmann factor for corrections to scaling in ferromagnetic materials and critical fluids, Eur. Phys. J. B 71 (2009) 1–6. T.L. Hill, R.V. Chamberlin, Fluctuations in energy in completely open small systems, Nano Lett. 2 (2002) 609–613. T.L. Hill, Thermodynamics of Small Systems, Dover, New York, 1994. R.V. Chamberlin, Mesoscopic mean-field theory for supercooled liquids and the glass transition, Phys. Rev. Lett. 82 (1999) 2520–2523. R.V. Chamberlin, Critical behavior from Landau theory in nanothermodynamic equilibrium, Phys. Lett. A 315 (2003) 313–318. M.R.H. Javaheri, R.V. Chamberlin, A free-energy landscape picture and Landau theory for the dynamics of disordered materials, J. Chem. Phys. 125 (2006) 154503. D.J. Evans, E.G.D. Cohen, G.P. Morriss, Probability of second law violations in shearing steady states, Phys. Rev. Lett. 72 (1993) 2401–2404. G.M. Wang, E.M. Sevick, E. Mittag, D.J. Searles, D.J. Evans, Experimental demonstration of violations of the second law of thermodynamics for small systems and short time scales, Phys. Rev. Lett. 89 (2002) 050601. R.V. Chamberlin, K.J. Stangel, Monte Carlo simulation of supercooled liquids using a self-consistent local temperature, Phys. Lett. A 350 (2006) 400–404. A. Wehrl, General properties of entropy, Rev. Modern Phys. 50 (1978) 221–260. H. Satoh, M. Takagi, K. Kinukawa, N. Kamegashira, Heat capacity of LaMnO3 , Therm. Acta 299 (1997) 123–126. C. Domb, Critical phenomena: a brief historical survey, Contemp. Phys. 26 (1985) 49–72. X. Qiu, T. Proffen, J.F. Mitchell, S.J.L. Billinge, Orbital correlations in the pseudocubic O and rhombohedral R phases of LaMnO3 , Phys. Rev. Lett. 94 (2005) 177203.