Squeezed ensemble for systems with first-order phase transitions

TL;DR

This paper introduces a new class of statistical ensembles, called squeezed ensembles, that accurately describe systems with first-order phase transitions, especially finite systems, improving computational efficiency and reducing finite-size effects.

Contribution

The paper proposes and derives a new class of ensembles that remain valid in the first-order transition region and are practical for finite systems, with formulas relating different ensembles.

Findings

Squeezed ensembles accurately capture equilibrium states in first-order transitions.

Formulas relate different ensembles, reducing finite-size effects.

Application to models confirms advantages of the new ensembles.

Abstract

All ensembles of statistical mechanics are equivalent in the sense that they give the equivalent thermodynamic functions in the thermodynamic limit. However, when investigating microscopic structures in the first-order phase transition region, one must choose an appropriate statistical ensemble. The appropriate choice is particularly important when one investigates finite systems, for which even the equivalence of ensembles does not hold. We propose a class of statistical ensembles, which always give the correct equilibrium state even in the first-order phase transition region. We derive various formulas for this class of ensembles, including the one by which temperature is obtained directly from energy without knowing entropy. Moreover, these ensembles are convenient for practical calculations because of good analytic properties. We also derive formulas which relate statistical-mechanical quantities of different ensembles, including the conventional ones, for finite systems. The formulas are useful for obtaining results with smaller finite-size effects, and for improving the computational efficiency. The advantages of the squeezed ensembles are confirmed by applying them to the Heisenberg model and the frustrated Ising model.

Figures

Equations

$${\sigma }_N(u)\equiv \frac{1}{N}\log g_N(u)$$

$$\underset{N\to \infty }{\lim }{\sigma }_N^{(n)}(u)=s^{(n)}(u).(n=0,1,2)$$

$${\hat{\rho }}_N^{\eta }\equiv \frac{e^{-N\eta (\hat{h})}}{{\Phi }_N^{\eta }},$$

$${\Phi }_N^{\eta }\equiv \mathrm{Tr}\left[e^{-N\eta (\hat{h})}\right].$$

$$\mathrm{Tr}\left[f(\hat{h})e^{-N\eta (\hat{h})}\right]$$

$${\beta }_N({\upsilon }_N^{\eta })\equiv {\sigma }_N^{\prime }({\upsilon }_N^{\eta })={\eta }^{\prime }({\upsilon }_N^{\eta }).$$

$${\xi }_N^{\eta }(u)={\xi }_N^{\eta }({\upsilon }_N^{\eta })-$$

$${{\xi }_N^{\eta }}^{\prime \prime }({\upsilon }_N^{\eta })={\sigma }_N^{\prime \prime }({\upsilon }_N^{\eta })-{\eta }^{\prime \prime }({\upsilon }_N^{\eta })<0.$$

$$u_N^{\eta }\equiv \mathrm{Tr}\left[\hat{h}{\hat{\rho }}_N^{\eta }\right]$$

$$\mathrm{Tr}\left[\hat{A}{\hat{\rho }}_N^{\mathrm{mic}}(u_N^{\eta })\right]=\mathrm{Tr}\left[\hat{A}{\hat{\rho }}_N^{\eta }\right]\left(1+O(N^{-1})\right).$$

$$\mathrm{Tr}\left[{\hat{a}}^2{\hat{\rho }}_N^{\eta }\right]-\mathrm{Tr}{\left[\hat{a}{\hat{\rho }}_N^{\eta }\right]}^2=o(N^0),$$

$$\frac{1}{N}\log {\Phi }_N^{\eta }$$

$$u_N^{\eta }$$

$${\sigma }_N(u_N^{\eta })=\frac{1}{N}\log {\Phi }_N^{\eta }+\eta \left(u_N^{\eta }\right)+O(N^{-1}\log N).$$

$${\sigma }_N(u_N^{\eta })=s_N^{\mathrm{vN}}({\hat{\rho }}_N^{\eta })+O(N^{-1}\log N),$$

$$s_N^{\mathrm{vN}}(\hat{\rho })\equiv -\frac{1}{N}\mathrm{Tr}\left[\hat{\rho }\log \hat{\rho }\right].$$

$${\beta }_N(u_N^{\eta })={\eta }^{\prime }\left(u_N^{\eta }\right)+O(N^{-1}).$$

$$c_N(u_N^{\eta })$$

$$=\frac{{\left({\beta }_N(u_N^{\eta })\right)}^2}{\frac{1}{N\mathrm{Tr}\left[{\left(\hat{h}-u_N^{\eta }\right)}^2{\hat{\rho }}_N^{\eta }\right]}-{\eta }^{\prime \prime }(u_N^{\eta })+O(N^{-1})}.$$

$$\frac{\partial {\upsilon }_N^{\eta }}{\partial \kappa }(\kappa )$$

$$\underset{\kappa }{\inf }{u}_N^{\eta }(\kappa )={ϵ}^{\min },\underset{\kappa }{\sup }{u}_N^{\eta }(\kappa )=\overline{ϵ},$$

$${\psi }_N^{\eta }(\kappa )\equiv -\frac{1}{N}\log {\Phi }_N^{\eta }(\kappa ),$$

$$\frac{\partial {\psi }_N^{\eta }}{\partial \kappa }(\kappa )=\mathrm{Tr}\left[\frac{\partial \eta }{\partial \kappa }\left(\kappa ,\hat{h}\right){\hat{\rho }}_N^{\eta }(\kappa )\right].$$

$${\psi }^{\eta }(\kappa )\equiv \underset{N\to \infty }{\lim }{\psi }_N^{\eta }(\kappa ).$$

$$s(u)=\underset{\kappa }{\inf }\left\{\eta (\kappa ,u)-{\psi }^{\eta }(\kappa )\right\}.$$

$${\psi }^{\eta }(\kappa )=\underset{u}{\inf }\left\{\eta (\kappa ,u)-s(u)\right\}.$$

$$s_N^{\eta }(\kappa )\equiv \eta (\kappa ,u_N^{\eta }(\kappa ))-{\psi }_N^{\eta }(\kappa ).$$

$$s(u_{\infty }^{\eta }(\kappa ))=s_N^{\eta }(\kappa )+\delta s_N^{\eta }(\kappa ).$$

$${\psi }_N^{\eta }$$

$$u_N^{\eta }$$

$${\psi }_{\tilde{\eta }}$$

$$u_N^{\tilde{\eta }}$$

Full text

I Introduction

Using statistical mechanics, one can obtain not only the thermodynamic functions but also the density operator, by which one can investigate microscopic structures of equilibrium statesGibbs (1902); Landau and Lifshiz (1980). To obtain thermodynamic functions in the thermodynamic limit, one can employ any statistical ensemble because all statistical ensembles give the equivalent thermodynamic functions, i.e., the functions are Legendre transformations of each otherRuelle (1969). This useful property, called the equivalence of ensembles, holds even for the system which undergoes the first-order phase transition, where the thermodynamic functions exhibit the strongest singularities. By contrast, to investigate microscopic structures in the first-order phase transition region, one must choose an appropriate statistical ensemble.

For example, consider the liquid-gas phase transition of water at a pressure of $1\ \mathrm{atm}$. If one uses the canonical ensemble specified by temperature $T$, the phase transition takes place at a single point $T=100\ {}^{\circ}\mathrm{C}$. However, at this phase transition point, the molar ratio of the liquid and the gas phases can take various valuesGibbs (1902). Consequently, the equilibrium state changes discontinuously in temperature, and it is impossible to obtain equilibrium states for each molar ratio of liquid and gas phases using the canonical ensemble. By contrast, if one employs the microcanonical ensemble specified by energy, the phase transition takes place in a finite region of energy, called the phase transition region or coexisting regionShimizu (2008); Kastner and Pleimling (2009); Alder and Wainwright (1962). At every point in this region, the ensemble gives the correct equilibrium state, in which all macroscopic variables including the liquid-gas molar ratio are uniquely determined. When gradually heating up water in experiments one obtains a sequence of such equilibrium states in the transition region.

As seen from this example, one must choose an appropriate statistical ensemble, such as the microcanonical ensemble, to obtain a correct density operator and thereby investigate microscopic structures in the first-order phase transition region.

The appropriate choice of the ensemble is particularly important when one investigates finite systems because even the equivalence of ensembles (that holds in the thermodynamic limit) does not hold for finite systems Stump and Hetherington (1987); Gross (2001). That is, even when one is only interested in thermodynamic functions, one must choose an appropriate statistical ensemble in order to derive correct properties of finite systems from the functions. For example, even with short-range interactions, finite systems which undergo first-order phase transitions exhibit thermodynamic anomalies such as a negative specific heat Bixon and Jortner (1989); Labastie and Whetten (1990); Gross (1990, 1997). In fact, the recent technical development enabled the experimental realization of first-order phase transitions in small systems, and evidence of the negative specific heat was observed D’Agostino et al. (2000); Schmidt et al. (2001). Nevertheless, the canonical ensemble always gives positive specific heat. Moreover, it gives double peaks of the energy distributionJanke (1998), i.e., an unphysical state which is a classical mixture of macroscopically distinct states Penrose and Lebowitz (1971); Binder and Kalos (1980). To correctly obtain the negative specific heat and a physical equilibrium state, one must use another ensemble such as the microcanonical ensemble Stump and Hetherington (1987); Gross (2001); Junghans et al. (2006, 2008); Chen et al. (2009); Penrose and Lebowitz (1971); Binder and Kalos (1980); Tröster et al. (2012).

The use of an appropriate ensemble is also important for numerical studies of the systems which undergo the first-order phase transition. However, conventional ensembles are not appropriate enough. If one employs the canonical ensemble its energy distribution has double peaks, which are separated by exponentially suppressed phase coexisting states, near the transition pointSchierz et al. (2016). This degrades greatly the efficiency of the Monte Carlo calculations using the importance sampling with local update algorithms or the replica exchange method Hukushima and Nemoto (1996) (also called parallel tempering) Martin-Mayor (2007); Kim et al. (2010). Furthermore, the singularities of thermodynamic quantities at the first-order transition point are smeared significantly in the canonical ensemble for finite systems due to the large fluctuation Stump and Hetherington (1987); Hüller (1992, 1994). This makes it difficult to identify the phase transition and to determine its order Hüller (1994). By contrast, if such finite systems are studied using the microcanonical ensemble, the phase transitions are directly detected Challa and Hetherington (1988a, b); Behringer et al. (2005); Behringer and Pleimling (2006). Unfortunately, however, the microcanonical ensemble has technical difficulties in practical calculations. It is difficult, especially for quantum systems, to construct a microcanonical ensemble. Moreover, one needs to differentiate the entropy in order to calculate the temperature, but it gives very noisy results in numerical calculationKanki et al. (2005).

Several attempts were made to overcome these problems. For example, the Gaussian ensemble Hetherington (1987); Challa and Hetherington (1988b, a); Johal et al. (2003) and the dynamical ensemble Gerling and Hüller (1993) were conceived as elaborate numerical methods for classical systems. Furthermore, the generalized canonical ensemble Costeniuc et al. (2005, 2006); Toral (2006) was introduced, which gives the entropy in the thermodynamic limit via the Legendre transformation even when the entropy is not concave. While its mathematical aspects were studied, the physical aspects were not discussed, such as the physical properties of the state described by that ensemble.

In this paper, we propose a class of statistical ensembles, which we call the squeezed ensembles. They always give the correct equilibrium state even in the first-order phase transition region. In particular, thermodynamic anomalies, such as negative specific heat, are correctly obtained, which appear generally in the transition region for finite systems with short-range interactions. We derive various formulas for this class of ensembles, including the one by which temperature is obtained directly from energy without knowing entropy.

Moreover, the squeezed ensembles are convenient for practical calculations because of good analytic properties. They can be numerically constructed more easily than the microcanonical ensemble, and the construction is even easier than that of the canonical ensemble in some cases. Furthermore, efficient numerical methods, such as the replica exchange method, are applicable in almost the same manner as in the canonical ensemble.

We also derive formulas which relate statistical-mechanical quantities of different ensembles, including the conventional ones, for finite systems. The formulas are useful for obtaining results with smaller finite size effects, and for improving the computational efficiency. The advantages of the squeezed ensembles and these formulas are confirmed by applying them to the Heisenberg model and the frustrated Ising model.

Various ensembles, including the Gaussian and the dynamical ensembles of the previous works, are included in the class of squeezed ensembles. One can choose an appropriate squeezed ensemble depending on the purpose, without losing the above advantages. By contrast, the conventional ensembles, such as the canonical and microcanonical ensembles, are understood as certain limiting cases of the squeezed ensembles so that some of the advantages are lost by the limiting procedure.

II Squeezed ensemble

II.1 Definition

We consider a quantum system which has $N$ degrees of freedom and the Hamiltonian $\hat{H}$. To take the thermodynamic limit, we use $\hat{h}\equiv\hat{H}/N$ and the energy density $u\equiv\mathrm{energy}/N$ 111Precisely speaking, $u$ is the mean energy. It agrees with the energy density only when a single uniform phase is realized in the equilibrium state. For simplicity, we use the term “energy density” throughout this paper even when several phases coexist.. We assume that all quantities are nondimensionalized with an appropriate scale. We denote the minimum and the maximum eigenvalues of $\hat{h}$ by $\epsilon^{\mathrm{min}}$ and $\epsilon^{\mathrm{max}}$, respectively.

We assume that the equilibrium state is specified by the energy density for each value of $N$. In other words, we assume that the state described by the microcanonical ensemble is not a classical mixture of macroscopically distinct states.

We also assume that the system is consistent with thermodynamics in the sense that

[TABLE]

converges to an $N$-independent concave function (entropy density, $s$) as $N\to\infty$, where $g_{N}(u)$ denotes the density of microstates222 More precisely, $\sigma_{N}$ is a twice continuously differentiable function which closely approximates $\frac{1}{N}\log g_{N}$ and satisfies

$\displaystyle\lim_{N\to\infty}\sigma_{N}^{(n)}(u)=s^{(n)}(u).\qquad(n=0,1,2)$

(2)

We assume the existence of such $\sigma_{N}$. The validity of analysis based on this assumption will be checked numerically in Section VII.. For the moment, we assume that $\sigma_{N}$ is also a concave function for finite $N$. Later on, it turns out that this assumption is unnecessary. Our formulation is valid even for systems whose concavity of $\sigma_{N}$ is broken.

We introduce the squeezed ensemble. Let $\eta$ be a convex function on $\left[\epsilon^{\mathrm{min}},\epsilon^{\mathrm{max}}\right]$. We define the squeezed ensemble associated with $\eta$ by

[TABLE]

where

[TABLE]

When $\eta(u)=\beta u$, $\hat{\rho}_{N}^{\eta}$ gives the canonical ensemble. When $\eta(u)=\left(\frac{u-\varepsilon}{\delta}\right)^{2n}$, $\hat{\rho}_{N}^{\eta}$ approaches the microcanonical ensembles as $n\to\infty$. We will show that other, appropriate forms of $\eta(u)$ give better ensembles.

II.2 Requirements on $\eta$

As discussed in Section I, the canonical ensemble gives an unphysical state at the first-order phase transition point. To make the squeezed ensembles free from such deficiency, we require

1. (A)

$\eta$ is a strongly convex function.

To calculate temperature easily (using Eq. (17) below) and to simplify the analysis, we also assume

1. (B)

$\eta$ is a twice continuously differentiable function.

These conditions ensure that the squeezed ensembles give the correct equilibrium state even in the first-order phase transition region, as follows.

We examine how the energy density distributes in $\hat{\rho}_{N}^{\eta}$. Let $f$ be an $N$-independent function. Then

[TABLE]

where $\xi_{N}^{\eta}(u)\equiv\sigma_{N}(u)-\eta(u)$. $\xi_{N}^{\eta}(u)$ takes the maximum at $u=\upsilon_{N}^{\eta}$ which satisfies

[TABLE]

Expanding $\xi_{N}^{\eta}$ around $\upsilon_{N}^{\eta}$ and noting ${\xi_{N}^{\eta}}^{\prime\prime}(\upsilon_{N}^{\eta})<0$, we get

[TABLE]

Hence, in the vicinity of $\upsilon_{N}^{\eta}$, $e^{N\xi_{N}^{\eta}}$ behaves as the Gaussian distribution, peaking at $\upsilon_{N}^{\eta}$, with the small variance $\frac{1}{N\left|{\xi_{N}^{\eta}}^{\prime\prime}(\upsilon_{N}^{\eta})\right|}$ (Fig. 1).

Unlike the canonical ensemble, $e^{N\xi_{N}^{\eta}}$ has a sharp peak even when $\sigma_{N}^{\prime\prime}(\upsilon_{N}^{\eta})=0$ in the first-order phase transition region thanks to the strong convexity of $\eta$, i.e.,

[TABLE]

Therefore, as proven in Appendix A, $\hat{\rho}_{N}^{\eta}$ represents the equilibrium state specified by the energy density

[TABLE]

in the following senses:

1. (i)

We can obtain the expectation value of any “mechanical variable” $\hat{A}$ in the microcanonical ensemble $\hat{\rho}_{N}^{\mathrm{mic}}$ with energy density $u_{N}^{\eta}$ as 333 We take the energy width of $\hat{\rho}_{N}^{\mathrm{mic}}$ as specified in Appendix A.

[TABLE]

Here, by mechanical variable, we mean a local observable (i.e., an observable on a continuous $O(N^{0})$ sites), such as the two-point correlation functions, or the sum of local operators, such as the total magnetization. Hence, both ensembles give the same result in the thermodynamic limit, even in the first-order phase transition region. 2. (ii)

Any macroscopic additive observable $\hat{A}$ has small variance as

[TABLE]

where $\hat{a}\equiv\hat{A}/N$.

Obviously, the state described by this ensemble depends on $\eta$. In order to obtain the states specified by a series of energies, it is convenient for practical calculations if $\eta$ depends on a parameter. We will discuss this parameter dependence in Section III.

II.3 Genuine thermodynamic variables

We can also obtain genuine thermodynamic variables such as the entropy and temperature. Using Eqs. (5)-(8) and applying Laplace’s method, we have

[TABLE]

Therefore, we obtain

[TABLE]

It is sometimes convenient to rephrased this relation as

[TABLE]

where $s_{N}^{\mathrm{vN}}$ is the von Neumann entropy density,

[TABLE]

Using Eqs. (6) and (13), we also obtain

[TABLE]

Here, since $\eta^{\prime}=\sigma_{N}^{\prime}$ only at $\upsilon_{N}^{\eta}$, replacing $\upsilon_{N}^{\eta}$ with $u_{N}^{\eta}$ yields the difference of $O(N^{-1})$.

In the thermodynamic limit, $\sigma_{N}$ and $\beta_{N}$ converge to the thermodynamic entropy density and inverse temperature, respectively, which are well-defined in the thermodynamic limit. Therefore, one can obtain the temperature of the equilibrium state specified by $u_{N}^{\eta}$ just by calculating $u_{N}^{\eta}$, via Eq. (17). By contrast, in order to calculate the temperature using the microcanonical ensemble, one needs to differentiate the entropy, and it gives very noisy results in numerical calculationKanki et al. (2005).

In a similar manner, we obtain

[TABLE]

In the thermodynamic limit, $c_{N}$ converges to the thermodynamic specific heat.

II.4 Interpretation

Although the above results have been derived naturally using Laplace’s method, the following scenario may be more intuitive from the viewpoint of the principle of equal weight.

Let us consider the equilibrium state of the target system which is in weak thermal contact with an external system. We note the ratio of the degrees of freedom of the external system to that of the target system. Using the principle of equal weight, one can obtain the canonical ensemble as the ensemble of the target system if the external system is much larger than the target system, and the microcanonical ensemble if the external system is much smaller than the target system. Then, we now consider the case where the external system has about the same degrees of freedom as the target system. In this case, the details of the external system affect the state of the target system and there are infinitely many ensembles as the state of the target system that converge the same equilibrium state as $N\to\infty$.

Suppose that the target system is in weak thermal contact with the external system which has the same degrees of freedom as the target system and that the restriction to the interval $[-\epsilon^{\mathrm{max}},-\epsilon^{\mathrm{min}}]$ of the $\sigma_{N}$ of the external system is equal to $-\eta(-u)$. Assume that the target system plus the external system together are isolated, with fixed total energy [math]. Then, applying the principle of equal weight to the total system, we get $\hat{\rho}_{N}^{\eta}$ as an ensemble of the target system. We can extract all statistical-mechanical quantities about the target system because we are familiar with $\sigma_{N}$ of the external system.

III Parameter of squeezed ensemble

Suppose that $\eta$ depends on a certain parameter. The equilibrium states described by the squeezed ensembles are specified by the parameter. We now examine the parameter dependencies of the physical quantities defined by the squeezed ensemble.

III.1 Energy density

Let $K$ be a real interval for the parameter $\kappa$, and $\eta$ be a function on $K\times[\epsilon^{\mathrm{min}},\epsilon^{\mathrm{max}}]$. Assume that $\eta(\kappa,\cdot)$ satisfies conditions (A)-(B) for all $\kappa$. Then, for each $\kappa$, $\eta(\kappa,\cdot)$ corresponds to a squeezed ensemble. Thus, a quantity related to the squeezed ensemble can be regarded as a function of $\kappa$. To examine the parameter dependence, we assume

1. (C)

$\eta$ is a twice continuously differentiable function of two variables.

For treating the systems which undergo the first-order phase transition, it is necessary that the energy density $u_{N}^{\eta}(\kappa)$ in the squeezed ensemble takes every possible value of the energy density of the system. The canonical ensemble does not satisfy this condition because $\kappa$ of the canonical ensemble corresponds to the inverse temperature $\beta$. In the thermodynamic limit, $u_{N}^{\mathrm{can}}(\beta)$ changes discontinuously in $\beta$ at the first-order phase transition point. Hence, as discussed in Section I, the canonical ensemble is unable to describe the equilibrium states in the first-order phase transition region, in which the energy density takes continuous values. [As we will discuss in Section IV, the correct state is not obtained even if the system is finite.] By contrast, in the case of the squeezed ensemble, the energy density changes continuously in $\kappa$ thanks to the strong convexity of $\eta(\kappa,\cdot)$ for each $\kappa$. In fact,

[TABLE]

is finite even in the first-order phase transition region. Furthermore, we take $\eta$ such that

1. (D)

[TABLE]

where $\bar{\epsilon}$ is the arithmetic mean of the eigenvalues of $\hat{h}$. Then, $u_{N}^{\eta}(\kappa)$ takes every possible value of $u$ in the physical region $\epsilon^{\mathrm{min}}<u<\bar{\epsilon}$ 444 We say the region $\epsilon^{\mathrm{min}}<u<\bar{\epsilon}$ is physical because temperature is positive in this region..

III.2 Thermodynamic function

As in the case of the conventional ensembles, we consider the logarithm of the “partition function”

[TABLE]

which is a function not of a physical quantity (such as $\beta$) but of our parameter $\kappa$. In numerical calculations (using, e.g., the Monte Carlo calculation), $\psi_{N}^{\eta}(\kappa)$ can be obtained easily by integrating

[TABLE]

Here, the right hand side is obtained simply by calculating the expectation value of $\displaystyle\frac{\partial\eta}{\partial\kappa}\left(\kappa,\hat{h}\right)$.

Let us define the thermodynamic function associated with $\eta$ as the thermodynamic limit of $\psi_{N}^{\eta}$:

[TABLE]

As proven in Appendix B, $\psi^{\eta}$ is equivalent to the thermodynamic entropy density in the following sense:

[TABLE]

Using this relation, one can obtain the thermodynamic entropy from $\psi^{\eta}$ without knowing $u_{N}^{\eta}$. We can also invert this relation as

[TABLE]

From a physical point of view, these relations are a generalization of the equivalence of the entropy density and the canonical free energy density. From a mathematical point of view, this is a generalization of the Legendre transformation, to which it reduces for $\eta(\beta,u)=\beta u$.

It is worth mentioning that Eqs. (23)-(24) are valid even when the concavity of the entropy is broken, as long as conditions (E)-(F) of Section IV are satisfied. Even in such a case, we can obtain entropy from $\psi^{\eta}$, as schematically explained in Fig. 2-2. By contrast, the Legendre transformation does not preserve the information in the non-convex function (Fig. 2).

IV Non-concave $\sigma_{N}$

So far we have assumed that $\sigma_{N}$ is concave. Since we consider systems with short-range interactions, this assumption is valid in the thermodynamic limit. However, for finite $N$, it was suggested that the concavity of $\sigma_{N}$ is broken in a system which undergoes a first-order phase transition with a phase separationBixon and Jortner (1989); Labastie and Whetten (1990); Gross (1990, 1997), as illustrated schematically in Fig. 3. In Appendix C, we prove this under reasonable conditions.

In order to obtain the correct equilibrium state even in such a case, we require

1. (E)

$\sigma_{N}-\eta$ has a single peak

and

1. (F)

$\sigma_{N}^{\prime\prime}-\eta^{\prime\prime}$ is negative at the peak.

In fact, under these conditions, the energy density distribution in $\hat{\rho}_{N}^{\eta}$ can still be well approximated by the Gaussian distribution, and all results in Sections II-III are valid.

Therefore, even for a system whose $\sigma_{N}$ is not concave, we can obtain the equilibrium states and the statistical-mechanical quantities using a squeezed ensemble associated with an appropriately chosen $\eta$. By contrast, using the canonical ensemble, the equilibrium states with the energy density in the region other than the black dotted lines in Fig. 3 cannot be obtained (see Appendix D).

V Physical natures of a system with non-concave $\sigma_{N}$

The concavity breaking of $\sigma_{N}$ causes thermodynamic anomalies. In this section, we discuss physical natures of equilibrium states with non-concave $\sigma_{N}$.

V.1 Realization of equilibrium states with non-concave $\sigma_{N}$

In general, an isolated system with sufficiently complex dynamics evolves into (relaxes to) an equilibrium state. During this “thermalization” process, the energy density keeps the same value as that of the initial state, which may be a non-equilibrium state such as a local equilibrium state. In this way, one can obtain equilibrium states of any possible value of the energy density by tuning the energy density of the initial state.

In particular, one can obtain an equilibrium state in the region where concavity of $\sigma_{N}$ is broken by tuning the energy density of the initial state in such a region. Such an experiment will be possible in various finite systems, such as cold atoms.

A more practical way is to heat the system slower enough (across the transition temperature) so that a quasistatic process is realized.

V.2 Measurability of statistical-mechanical quantities

As illustrated schematically in Fig. 3, $\beta_{N}$ has an S-shaped curve. It is a physical quantity that can be measured as follows.

Let place the target system in weak thermal contact with an external system. We assume the external system is sufficiently small so that its effect on the target system is negligible. That is, the external system works as a thermometer and the target system as a heat bath. After the thermal equilibrium of the total system is reached, the thermometer is in the canonical Gibbs state with $\beta=\beta_{N}(u)$, where $u$ is the energy density of the target system, because the canonical typicality Sugita (2006, 2007); Popescu et al. (2006); Goldstein et al. (2006) holds regardless of the concavity of $\sigma_{N}$ of the heat bath. Hence, one can read the value of $\beta_{N}(u)$ from the thermometer. Using this setup with various values of $u$, one can measure the function $\beta_{N}$ and $c_{N}$ of the target system. Then one will find not only the S-shaped $\beta_{N}$ but also the anomalous behaviors of $c_{N}$. The latter is illustrated schematically in Fig. 3, where $c_{N}$ takes negative values between the two singular points. See Refs. Stodolsky (1995); Chomaz and Gulminelli (1999); Schmidt et al. (1997) for related discussions.

Evidence of the concavity breaking of $\sigma_{N}$ due to the first-order phase transition has been observed both in experimentsD’Agostino et al. (2000); Schmidt et al. (2001) and in dynamical simulationsLitz et al. (1992); Nielsen et al. (1994); Reyes-Nava et al. (2003).

V.3 Thermodynamic stability

Even when $\sigma_{N}$ is not concave the equilibrium states of an isolated finite system is stable because of the energy conservations. In this subsection, we discuss whether the equilibrium states are stable when the system is in thermal contact with an infinitely large heat bath.

Suppose that the system is initially isolated and in the equilibrium state with energy density $u$. Then let the system, which we now call the target system, interact with a large heat bath of inverse temperature $\beta$ that is equal to $\beta_{N}(u)$. According to the canonical typicalitySugita (2006, 2007); Popescu et al. (2006); Goldstein et al. (2006), the target system is driven by the heat bath toward the canonical Gibbs state. However, as discussed in Appendix D, if $u$ lies in the red solid or blue dashed region of Fig. 3 the equilibrium state of the target system is not realized as a canonical Gibbs state. Therefore, for such $u$, the equilibrium state of the target system is no longer stable when attached to a heat bath.

More specifically, if $u$ lies in the blue dashed region of Fig. 3 the target system quickly evolves into the canonical Gibbs state of the same temperature by absorbing or emitting energy from or to the heat bath. That is, the equilibrium state for such $u$ is thermodynamically unstable. On the other hand, if $u$ lies in the red solid region the equilibrium state of the target system is thermodynamically metastable. That is, the state is maintained for a certain macroscopic time scale because $\sigma_{N}(u)$ is locally concave in such a region Antoni et al. (2004).

We can thus classify three types of regions in Fig. 3, for the case where the target system is in thermal contact with a heat bath, as the stable equilibrium states (black dotted line), the metastable states (red solid line), and the unstable states (blue dashed line).

VI Relation with other ensembles in finite systems

In finite systems, different ensembles are not completely equivalent. The inequivalence is most prominent in the first-order transition region. Let us investigate the effects of the inequivalence.

VI.1 Inequivalence of ensembles in finite systems

As an example, we consider the entropy density. What one can calculate using statistical mechanics is a sequence indexed by $N$ such as

[TABLE]

As $N\to\infty$, this $s_{N}^{\eta}(\kappa)$ converges to thermodynamic entropy density $s(u)$ at $u=u_{\infty}^{\eta}(\kappa)$. In other words, there exists $o(N^{0})$ quantity $\delta s_{N}^{\eta}(\kappa)$ such that

[TABLE]

On the other hand, for another pair of $(\tilde{\eta},\tilde{\kappa})$ which satisfies $u_{\infty}^{\tilde{\eta}}(\tilde{\kappa})=u_{\infty}^{\eta}(\kappa)$, $s_{N}^{\tilde{\eta}}(\tilde{\kappa})$ also converges to $s(u_{\infty}^{\eta}(\kappa))$, and there exists $o(N^{0})$ quantity $\delta s_{N}^{\tilde{\eta}}(\tilde{\kappa})$ such that $s(u_{\infty}^{\eta})=s_{N}^{\tilde{\eta}}(\tilde{\kappa})+\delta s_{N}^{\tilde{\eta}}(\tilde{\kappa})$. Generally, $s_{N}^{\eta}(\kappa)$ and $s_{N}^{\tilde{\eta}}(\tilde{\kappa})$ are different in finite systems. The difference is $o(N^{0})$ because $s_{N}^{\eta}(\kappa)-s_{N}^{\tilde{\eta}}(\tilde{\kappa})=\delta s_{N}^{\eta}(\kappa)-\delta s_{N}^{\tilde{\eta}}(\tilde{\kappa})$. Note that the rates at which $\delta s_{N}^{\eta}(\kappa)$ and $\delta s_{N}^{\tilde{\eta}}(\tilde{\kappa})$ converge to [math] are different in general.

The same is true for other thermodynamic quantities. Suppose that one calculates a thermodynamic quantity of $O(N^{0})$. What one can calculate using statistical mechanics is a function sequence indexed by $N$ that converges to the thermodynamic quantity as $N\to\infty$. Therefore, there is an arbitrariness in the choice of the function sequence up to $o(N^{0})$. In order to emphasize the distinction from a thermodynamic quantity, we call the function sequence that converges to the thermodynamic quantity a “statistical-mechanical quantity”.

VI.2 Formulas relating different ensembles

One might think that it is meaningless to calculate quantities with such an arbitrariness accurately. However, it is important in numerical calculation, where one has to deal with finite systems inevitably, for the following reasons. The rates of convergence of the statistical-mechanical quantities to the thermodynamic quantities with increasing $N$ are different among ensembles. If we can choose an ensemble which converges very quickly, it is advantageous for practical calculations. On the other hand, there may also be an ensemble that gives thermodynamic quantities that converge very slowly, even though it is equivalent to other ensembles in the thermodynamic limit. Using the following formulas, we can evaluate the difference between the statistical-mechanical quantity calculated using a squeezed ensemble and that calculated using a conventional ensemble (such as the canonical ensemble and microcanonical ensemble). If necessary, we can correct the difference from the conventional ensemble.

It is known that, for a system without phase transition, the statistical-mechanical quantities in the canonical ensemble, among those in various Gibbs ensembles, are closest to thermodynamic quantities of the infinite systemIyer et al. (2015). Therefore, it is useful to derive formulas by which a squeezed ensemble gives statistical-mechanical quantities in the canonical ensemble. Such formulas are obtained as follows.

First, let us derive formulas by which a squeezed ensemble gives statistical-mechanical quantities in another squeezed ensemble. Below in this section, we assume that $\eta$ is differentiable as many times as necessary. Using Eqs. (5)-(7) and (20), we have

[TABLE]

for any choice of $\eta$ that satisfies conditions (A)-(B). Suppose that two squeezed ensembles associated with $\eta$ and $\tilde{\eta}$ satisfy $\upsilon_{N}^{\eta}=\upsilon_{N}^{\tilde{\eta}}$. Then, we have

[TABLE]

These formulas relate statistical-mechanical quantities between different squeezed ensembles.

VI.3 Formulas relating the squeezed ensemble to the canonical ensemble

Next, using these formulas, let us derive formulas by which the squeezed ensemble gives the canonical entropy density $s^{\mathrm{can}}_{N}$ and the canonical energy density $u^{\mathrm{can}}_{N}$. Assume that $\sigma_{N}^{\prime\prime}(\upsilon_{N}^{\eta})<0$. This assumption is satisfied in many cases except when $\upsilon_{N}^{\eta}$ is in the phase transition region. Set $\tilde{\eta}=\beta u$, where

[TABLE]

Then, the squeezed ensemble associated with $\tilde{\eta}$ gives the canonical ensemble at inverse temperature $\beta$. Using the relation

[TABLE]

we get

[TABLE]

These are the desired formulas.

Though $\upsilon_{N}^{\eta}$ is unknown, using Eq. (28), we find

[TABLE]

One can calculate $\xi_{N}^{\eta}$’s derivatives at $\upsilon_{N}^{\eta}$ from the central moments of energy distribution. For example, as the second and third order derivatives, we have

[TABLE]

Using Eqs. (35)-(37), we can calculate $\upsilon_{N}^{\eta}$ and $\xi_{N}^{\eta}$’s derivatives at $\upsilon_{N}^{\eta}$ with accuracy up to $O(N^{-1})$. Substituting these into Eqs. (31)-(34), we can calculate the statistical-mechanical quantities in the canonical ensemble from the statistical-mechanical quantities in the squeezed ensemble with accuracy up to $O(N^{-1})$.

VII Application to Heisenberg model

Before studying a system which undergoes a first-order phase transition, we confirm the validity of our results by applying our formulation to the Heisenberg chain, defined by the Hamiltonian

[TABLE]

where $J=+1$ (ferromagnetic). The exact results at finite temperature have been derived for $N\to\infty$Nakamura and Takahashi (1994).

We can freely choose $\eta$ for practical convenience according to the physical situation. Here, we choose $K=(0,\infty)$ and $\eta(\kappa,u)=-2\kappa\log(l-u)$, which is particularly convenient for quantum systems, because $\hat{\rho}_{N}^{\eta}$ and $\Phi_{N}^{\eta}$ are obtained by simply multiplying the Hamiltonian repeatedly. In this model, $-0.75\leq\epsilon^{\mathrm{min}}<\bar{\epsilon}<\epsilon^{\mathrm{max}}\leq+0.75$. Hence, we take $l=1(>\epsilon^{\mathrm{max}})$, then $\eta(\kappa,\cdot)$ satisfies conditions (A)-(B) for all $\kappa>0$. Furthermore, $\eta$ satisfies conditions (C)-(D).

We calculate statistical-mechanical quantities in the squeezed ensemble using thermal pure quantum formulationSugiura and Shimizu (2012). Then, using Eqs. (12)-(17), we calculate $\psi_{N}^{\mathrm{can}}$, with an error of $O(N^{-1}\log N)$, from the statistical-mechanical quantities in the squeezed ensemble as

[TABLE]

which corresponds to the result obtained in the previous work Sugiura and Shimizu (2013). The result at $\beta\simeq 2.82$ is plotted by the red crosses in Fig. 4.

Furthermore, we correct the difference from the canonical free energy density using Eqs. (31)-(37), which reduces an error to $O(N^{-2})$. The result is plotted by the blue crosses in Fig. 4. It is confirmed that the difference from the canonical free energy density is proportional to $N^{-2}$.

VIII Application to Frustrated Ising model

To confirm the advantages of our formulation for studying the systems which undergo first-order phase transitions, we apply our formulation to the two-dimensional frustrated Ising model on the square lattice, defined by the Hamiltonian

[TABLE]

Here, $\hat{\sigma}^{z}_{i}=\pm 1$ is a classical variable, and $\langle i,j\rangle$ and \mathopen{\hbox{\set@color{\langle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\langle}}}i,j\mathclose{\hbox{\set@color{\rangle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\rangle}}} denote the nearest and the next-nearest neighbors, respectively. We take both couplings $J_{1},J_{2}$ positive, and measure energy density and temperature in units of $J_{1}$. It is known that the model undergoes a weak first-order phase transition for $0.5<J_{2}/J_{1}\lesssim 0.67$Jin et al. (2012). We take $J_{2}/J_{1}=0.6$. The numerical simulations are performed for systems of size $L\times L$ with periodic boundary conditions.

We choose $K=(-\infty,\bar{\epsilon})$ and $\eta(\kappa,u)=\frac{1}{2}\lambda(u-\kappa)^{2}$. This $\eta$ satisfies conditions (A)-(D) for positive $\lambda$ and defines the so-called Gaussian ensemble introduced by HetheringtonHetherington (1987). As proven in Appendix C, the concavity of $\sigma_{N}$ is broken for finite $N$. Since the second order derivative of $\eta(\kappa,\cdot)$ is given by $\lambda$, $\eta$ satisfies conditions (E)-(F) for sufficiently large $\lambda$ even in the case of the non-concave $\sigma_{N}$. In the classical systems, it is easy to calculate the energy of a given configuration and the energy change due to a change in the local configuration. Therefore, this choice of $\eta$ is convenient for classical systems which undergo first-order phase transitions.

We calculate the expectation values of mechanical variables in the squeezed ensemble using the Monte Carlo calculations. The acceptance probability can be easily computed using the Metropolis algorithmMetropolis et al. (1953). As argued in Section III, the equilibrium state changes continuously in $\kappa$. Therefore the replica exchange methodHukushima and Nemoto (1996) works wellKim et al. (2010) even in the first-order phase transition region, unlike the canonical ensemble. Other advantages of the Gaussian ensemble for studying the phase transitions were already studied in Ref. Challa and Hetherington (1988b, a). These advantages are enjoyed also by other choices of $\eta$ as long as conditions (E)-(F) are satisfied.

Fig. 5 shows the relation between $u$ and $\beta_{N}$ for $L=64$. We calculate $\beta_{N}$ using the canonical ensemble and the squeezed ensembles with various values of $\lambda$. When the canonical ensemble is used, $u$ is given as a function of $\beta$, which is a single-valued function, changing monotonically and continuously, even when $\sigma_{N}$ is not concave. Hence, the effect due to the first-order phase transition is greatly diminished. This should be contrasted with the results of the squeezed ensemble, which show that, for sufficiently large $\lambda$, $\beta_{N}$ becomes a multi-valued, S-shaped function in the transition region. Consequently, one can correctly obtain the negative specific heat due to the first-order phase transition by using the squeezed ensemble.

It is also seen from Fig. 5 that the functional form of $\beta_{N}(u)$ without correction (full circles) becomes insensitive to the magnitude of $\lambda$ for sufficiently large $\lambda$. This is because the error in $\beta_{N}$ by Eq. (17) scales as $O(N^{-1}\lambda^{-1})$. In fact, Eq. (28) gives

[TABLE]

Therefore, the larger $\lambda$ gives the smaller error, when using Eq. (17). On the other hand, the computational efficiency decreases with $\lambda$ because the acceptance probability for configurations with higher energy is $e^{-O(N\lambda)}$. To compromise these conflicting demands, we can take $\lambda$ small (to get a good efficiency) such as $\lambda=0.5$ and correct the result using the formulas derived in Section VI (to decrease the error). In fact, we can obtain the relation between $u$ and $\beta_{N}$ with an error of $O(N^{-2})$ from Eqs. (6) and (28) using the squeezed ensemble with small $\lambda$ (open circles in Fig. 5), and it agrees well with $\beta_{N}$ calculated using the squeezed ensemble with larger $\lambda$.

As discussed in Sec. IV, $c_{N}$ takes negative values in the transition region for finite $N$. To confirm this fact, we calculate the relation between $u$ and $c_{N}$ for $L=64$ from Eq. (18) using the squeezed ensemble with $\lambda=5$. The results are plotted in Fig. 6. It is confirmed that $c_{N}$ takes negative values between the two singular points.

Finally, we investigate how the result approaches that of the infinite system. Fig. 7 shows the relation between $\beta_{N}$ and $u$ for various values of $L$, where $\beta_{N}$ is calculated using the squeezed ensemble with $\lambda=5$. With increasing $L$, the S-shaped region shifts towards the transition temperature of the infinite system, and this makes it possible to precisely extrapolate the transition temperature Challa et al. (1986); Borgs and Kotecký (1992); Borgs et al. (1991); Challa and Hetherington (1988a). As $L$ increases, the amplitude of the S-shape decreases. This implies that the concavity of $\sigma_{N}$ is broken due to the finite system size.

IX Summary

To summarize, we have proposed the squeezed ensembles (Sections II-III). They give the correct equilibrium state even in the first-order phase transition region, as the microcanonical ensemble does. In particular, for finite systems, one can correctly obtain thermodynamic anomalies such as negative specific heat, which appear generally in the transition region even when interactions are short-ranged (Sections IV-V and Appendix C). Moreover, the squeezed ensembles have good analytic properties, which yield useful analytic formulas (Section II.3) including the one by which temperature is obtained directly from energy without knowing entropy (Eq. (17)). The squeezed ensembles are convenient for practical calculations because they can be numerically constructed easily, and efficient numerical methods, such as the replica exchange method, are applicable straightforwardly. We also derive formulas which relate statistical-mechanical quantities of different ensembles for finite systems (Section VI). The formulas are useful for obtaining results with smaller finite-size effects, and for improving the computational efficiency. We have confirmed the advantages of the squeezed ensembles and these formulas by applying them to the Heisenberg model (Section VII) and the frustrated Ising model (Section VIII).

Acknowledgements.

We thank Y. Chiba, K. Hukushima, H. Tasaki, H. Hakoshima, and S. Sugiura for discussions. This work was supported by The Japan Society for the Promotion of Science, KAKENHI No. 15H05700 and No. 19H01810.

Appendix A Proof of Eqs. (10)-(11)

Consider systems with the translationally invariance. We make a reasonable assumption that, by taking an appropriate constant $C$ depending only on $N$, the expectation value of $\hat{a}\equiv\hat{A}/C$ in the microcanonical ensemble converges to an $N$-independent continuous function of $u$ as $N\to\infty$. Then, by symmetry, the expectation value approaches a constant multiple of the expectation value of an appropriate additive observable as $N\to\infty$.

We denote by $\hat{\rho}_{N}^{\mathrm{mc}}(u-\delta(u),u]$ the microcanonical ensemble with the energy shell $(u-\delta(u),u]$. There is an arbitrariness in the choice of $\delta$. Here, we take $\delta$ as follows. Let $\delta_{0}$ be a positive constant of $O(N^{0})$ such that $\epsilon^{\mathrm{min}}+\delta_{0}<\upsilon_{N}^{\eta}$, and let $\delta$ be a positive valued function such that

[TABLE]

for all $u^{\prime}\in[\epsilon^{\mathrm{min}}+\delta_{0},\epsilon^{\mathrm{max}}]$. Then it holds

[TABLE]

Using Eq. (45) and applying Laplace’s method, we have

[TABLE]

Therefore, we have Eq. (10).

Since we assumed that the equilibrium state is specified by the energy density, we have

[TABLE]

for all macroscopic additive observable $\hat{A}$. Then using Eqs. (46) and applying Laplace’s method, we have

[TABLE]

for all macroscopic additive observable $\hat{A}$. Therefore, we have Eq. (11).

Appendix B Proof of Eqs. (23)-(24)

Eq. (24) can be immediately obtained from Eq. (12). Eq. (23) is proved as follows. Using Eq. (24), we have

[TABLE]

for all $\kappa$. Then we have

[TABLE]

On the other hand, using Eq. (24), we have

[TABLE]

for all $\kappa$. The continuity of $\upsilon_{\infty}^{\eta}$ in $\kappa$ (Eq. (19)) implies that there exists $\kappa$ such that $\upsilon_{\infty}^{\eta}(\kappa)=u$, then it holds

[TABLE]

Therefore we have Eq. (23).

Appendix C Generality of Concavity breaking of $\sigma_{N}$

We prove that the concavity breaking always occurs for finite $N$ under reasonable conditions. For proper handling of a system which undergoes a first-order phase transition, we consider a more general case where the equilibrium state is specified using a general set of extensive variables.

C.1 Concavity breaking in liquid-gas systems

As an example, we take a $d$-dimensional system which undergoes a first-order liquid-gas phase transition (whereas general systems will be discussed in Appendix C.2). Its equilibrium state is assumed to be specified by the energy $U$, volume $V$ and number of particles $N$. Hence, for a given value of $N$, $\sigma_{N}$ is a function of the energy per particle, $u\equiv U/N$, and the volume per particle, $v\equiv V/N$, i,e.,

[TABLE]

We investigate properties of $\sigma_{N}(u,v)$, for each fixed value of $N$, in the state space spanned by $u$ and $v$. Note that in this state space the liquid-gas coexisting region is a two-dimensional region, whereas it is a one-dimensional region (line) in the state space spanned by temperature $T$ and pressure $P$.

Suppose that the system is in an equilibrium state where liquid and gas phases coexist. We assume that the thickness and surface area of phase boundaries are $O(N^{0})$ and $O(N^{1-1/d})$, respectively. Furthermore, we assume that finite size effects in each phase are smaller than those due to the phase boundaries. This condition is natural and reasonable, as discussed in Appendix C.4. Under the above conditions, we show that the concavity of $\sigma_{N}$ is always broken for finite $N$ in such an equilibrium state.

Let us compare the equilibrium states with the same values of $(u,v)$ for two cases where 8(a) the effects of the phase boundaries are ignored and 8(b) not (Fig. 8).

First, we consider case 8(a). We denote a physical quantity in this case by a tilde over its symbol, such as $\tilde{\sigma}_{N}$. Since phase boundaries are ignored, one can define quantities (such as the number of particles) in each phase without ambiguity. Let $\tilde{N}^{p}$, $\tilde{u}^{p}$, and $\tilde{v}^{p}$ be the number of particles, energy per particle, and volume per particle, respectively, in phase $p$. Since phase boundaries are ignored, the fraction $\tilde{\nu}^{p}\equiv\tilde{N}^{p}/N$ satisfies

[TABLE]

and $u$, $v$ and $\tilde{\sigma}_{N}$ agree with the weighted arithmetic means as

[TABLE]

Next, we consider case 8(b) where effects of phase boundaries are not ignored. For $U,V$ and $N$ of the whole system including the phase boundaries, we take their values same as those in case 8(a). Then the fraction and the state of the bulk of each phase change slightly (if possible) as compared with the case where the phase boundaries can be neglected. Hence, compared with $\tilde{\sigma}_{N}$, the $\sigma_{N}$ changes according to two main causes: the number of possible configurations of the phase boundaries and the changes in the fraction and the state of the bulk of each phase.

Let us compare $\sigma_{N}(u,v)$ and $\tilde{\sigma}_{N}(u,v)$. For this purpose, we consider an equilibrium state with the number of particles $M$ and the same values of $(u,v)$, where $1\ll M/N$. Let us compare the microstates for two cases where 9(c) surface area of the phase boundaries is $O(M^{1-1/d})$ and 9(d) $M/N$ times that in the equilibrium state with the number of particles $N$ (Fig. 9). The density of microstates 9(c) agree with $e^{M\left(\tilde{\sigma}_{N}(u,v)+O(M^{-1/d})\right)}$. On the other hand, the density of microstates 9(d) is larger than $e^{M\sigma_{N}(u,v)}$. Therefore, If $\sigma_{N}(u,v)$ were larger than $\tilde{\sigma}_{N}(u,v)$, in the sufficiently large system, the density of microstates 9(c) would be exponentially smaller than the density of microstates 9(d). Hence, the macrostate with the phase boundaries whose surface area is $O(M^{1-1/d})$ would not be realized as an equilibrium state which is a typical macrostate with the largest number of microstates. Since such an equilibrium state contradicts with our assumption, we conclude

[TABLE]

From Eqs. (63)-(67), we find

[TABLE]

for certain values of $\tilde{u}^{p},\tilde{v}^{p}$ and $\tilde{\nu}^{\mathrm{Gas}}\in(0,1)$. This shows that the concavity of $\sigma_{N}(u,v)$ is broken.

This concavity breaking occurs in equilibrium states where liquid and gas phases coexist, i.e., in the phase-coexisting region in the state space spanned by $(u,v)$.

C.2 Concavity breaking in general systems

The above discussions can easily be extended to general systems with short-range interactions, as follows.

Consider a $d$-dimensional system whose equilibrium states is specified by extensive variables $X_{1},X_{2},\cdots,X_{M}$ and $N$. Hence, for a given value of $N$, $\sigma_{N}$ is a function of $(x_{1},x_{2},\cdots,x_{M})$, where $x_{i}\equiv X_{i}/N$. Suppose that the system which undergoes a first-order phase transition and that several phases coexist in the phase transition region. We assume that the thickness and surface area of phase boundaries are $O(N^{0})$ and $O(N^{1-1/d})$, respectively. Furthermore, we assume that finite size effects in each phase are smaller than those due to the phase boundaries. Then concavity of $\sigma_{N}$ is always broken for finite $N$.

C.3 Concavity restoration in the thermodynamic limit

The anomalous behavior of $\sigma_{N}$ is peculiar to the finite systems and the concavity of $\sigma_{N}$ recovers in the thermodynamic limitLynden-Bell (1995); Ispolatov and Cohen (2001). Since the surface area of the phase boundaries is $O(N^{1-1/d})$. the contribution of the phase boundaries to $x_{i}$ and $\sigma_{N}$ is $O(N^{-1/d})$. On the other hand, in the model with the long-range interaction obtained using the mean-field approximation, the concavity of $\sigma_{N}$ is broken even in the thermodynamic limit. Therefore, an S-shaped caloric curve of such the model (so-called van der Waals loop) is an unphysical artifact of the approximation and must be distinguished from S-shaped $\beta_{N}$ in the finite system with short-range interaction.

C.4 Validity of the condition

In Appendices C.1-C.2, we have assumed that finite size effects can be neglected except for those due to the phase boundaries. Here we discuss the validity of this condition.

Let us compare the asymptotic behavior of the finite size effects due to the phase boundaries and those due to other causes at large $N$. We consider the case where surface effects can be neglected.

As mentioned in Appendix C.2, the finite size effects due to the phase boundaries on $\sigma_{N}$ are $O(N^{-1/d})$. With short-range interactions, the spatial dimension of the system which undergoes a first-order phase transition at finite temperature satisfies $d>1$Araki (1969, 1975). On the other hand, the finite size effects excluding those due to the phase boundaries on $\sigma_{N}$ is expected to be $O(N^{-1})$. In fact, the difference between the transition points in finite and infinite systems is $O(N^{-1})$Challa et al. (1986).

Therefore, for sufficiently large $N$, finite size effects can be neglected except for those due to the phase boundaries.

Appendix D Inapplicability of the canonical ensemble

We examine whether the canonical ensemble is applicable to the system with non-concave $\sigma_{N}$ (Fig. 3). The canonical ensemble corresponds to the case where $\eta=\beta u$, and the equilibrium state is specified by $\beta$. Since neither condition (A) nor (E)-(F) are satisfied by this choice of $\eta$, the following difficulty arises for systems with non-concave $\sigma_{N}$.

If $\sigma_{N}$ were strictly convex, the energy density distribution in the canonical ensemble would have a single peak. The peak position $\upsilon_{N}^{\mathrm{can}}$ would be uniquely determined as the solution to $\beta_{N}(\upsilon_{N}^{\mathrm{can}})=\beta$ (Eq. (6)). However, this property is lost in the case where $\sigma_{N}$ is not concave.

For example, suppose that one wants to investigate microscopic structures of the equilibrium state with the energy density $u^{*}$ on the blue dashed line in Fig. 3. This is impossible if one uses the canonical ensemble. In fact, if one takes $\beta=\beta^{*}\equiv\beta_{N}(u^{*})$ in the canonical ensemble, there are three points of $u$ (purple squares in Fig. 3) that give the same value of $\beta_{N}$. Consequently, the energy distribution becomes bimodal, as shown in Fig. 10, and takes the local minimum at $u^{*}$. Hence the desired equilibrium state is not obtained, and one cannot investigate its microscopic structures. This should be contrasted with a squeezed ensemble, whose energy distribution has the global maximum at $u=u^{*}$, giving the desired equilibrium state correctly.

Similarly, for $u$ on the red solid line, the energy density distribution in the canonical ensemble at $\beta=\beta_{N}(u)$ has the local (but not global) maximum at $u$.

To sum up, the equilibrium states with the energy density in the region other than the black dotted line in Fig. 3 cannot be obtained using the canonical ensemble.