Finite-size effects in canonical and grand-canonical quantum Monte Carlo simulations for fermions

TL;DR

This paper presents a finite-temperature quantum Monte Carlo method for fermionic models in the canonical ensemble, analyzing finite-size effects and comparing them to grand-canonical results, with applications to nuclear and cold atom systems.

Contribution

The authors develop a canonical ensemble quantum Monte Carlo method and derive exact finite-size correction formulas relating canonical and grand-canonical observables.

Findings

Finite-size corrections in the canonical ensemble scale as inverse volume.

In the grand-canonical ensemble, corrections decay exponentially with system size.

Theoretical predictions are validated through simulations of the Hubbard and Ising models.

Abstract

We introduce a quantum Monte Carlo method at finite temperature for interacting fermionic models in the canonical ensemble, where the conservation of the particle number is enforced. Although general thermodynamic arguments ensure the equivalence of the canonical and the grand-canonical ensembles in the thermodynamic limit, their approach to the infinite-volume limit is distinctively different. Observables computed in the canonical ensemble generically display a finite-size correction proportional to the inverse volume, whereas in the grand-canonical ensemble the approach is exponential in the ratio of the linear size over the correlation length. We verify these predictions by quantum Monte Carlo simulations of the Hubbard model in one and two dimensions in the grand-canonical and the canonical ensemble. We prove an exact formula for the finite-size part of the free energy density, energy density and other observables in the canonical ensemble and relate this correction to a susceptibility computed in the corresponding grand-canonical ensemble. This result is confirmed by an exact computation of the one-dimensional classical Ising model in the canonical ensemble, which for classical models corresponds to the so-called fixed-magnetization ensemble. Our method is useful for simulating finite systems which are not coupled to a particle bath, such as in nuclear or cold atom physics.

Equations

$$C_V=\frac{d⟨\hat{H}⟩}{dT}=k_B{\beta }^2\left(⟨{\hat{H}}^2⟩-⟨\hat{H}{⟩}^2\right),$$

$${\Xi }_c=\frac{d⟨\hat{N}⟩}{d\mu }=\beta \left(⟨{\hat{N}}^2⟩-⟨\hat{N}{⟩}^2\right),$$

$$\underset{N\to \infty }{\lim }\frac{\sqrt{\left(⟨{\hat{H}}^2⟩-⟨\hat{H}{⟩}^2\right)}}{⟨\hat{H}⟩}=\underset{N\to \infty }{\lim }\frac{\sqrt{\left(⟨{\hat{N}}^2⟩-⟨\hat{N}{⟩}^2\right)}}{⟨\hat{N}⟩}=0$$

$$\lambda {\left(\hat{N}-N\right)}^2,$$

$$Z_{\mathrm{gc}}(h)=\sum _{\{S_k=\pm 1\}}\exp \left\{\beta J\sum _{<ij>}{S}_i{S}_j+h\sum _i{S}_i\right\}.$$

$$\begin{array}{rl}{Z}_{\mathrm{can}}(h,m)=\sum _{\{S_k=\pm 1\}}& \exp \left\{\beta J\sum _{<ij>}{S}_i{S}_j+h\sum _i{S}_i\right\}\\ & \cdot \delta \left(m,\frac{1}{V}\sum _i{S}_i\right),\end{array}$$

$$e^{-V{U}_{\mathrm{eff}}(m,V)}=\int [\mathcal{D}\phi ]e^{-S[\phi ]}\delta \left(m-\frac{1}{V}\int d^{d}x\phi (x)\right),$$

$$\begin{array}{rl}& \hat{H}=\hat{T}+\hat{V},\\ & \hat{T}\equiv \sum _{x,y}{\hat{c}}_x^{†}{T}_{x,y}{\hat{c}}_y^{},\\ & \hat{V}\equiv \sum _k{U}_k{\left({\hat{V}}^{(k)}+{\alpha }_k\right)}^2,{\hat{V}}^{(k)}\equiv \sum _{x,y}{\hat{c}}_x^{†}{V}_{x,y}^{(k)}{\hat{c}}_y^{},\end{array}$$

$$e^{\Delta \tau \lambda {\hat{A}}^2}=\frac{1}{4}\sum _{l=\pm 1,\pm 2}\gamma (l)e^{\sqrt{\Delta \tau \lambda }\eta (l)\hat{A}}+\mathcal{O}(\Delta {\tau }^4),$$

$$\begin{array}{rl}{e}^{-\beta \hat{H}}=& \sum _{\left\{l_{k,\tau }\right\}}{e}^{S_0\left\{l_{k,\tau }\right\}}\prod _{\tau =1}^{L_{\tau }}{e}^{-\Delta \tau \hat{T}}\prod _k{e}^{\sqrt{-\Delta \tau U_k}\eta (l_{k,\tau }){\hat{V}}^{(k)}}\\ & +i\Delta \tau \hat{R}+O(\Delta {\tau }^2).\end{array}$$

$$Z_{\mathrm{gc}}=\mathrm{Tr}\left\{e^{-\beta \hat{H}}\right\}=\sum _{\left\{l_{k,\tau }\right\}}{e}^{S_0\left\{l_{k,\tau }\right\}}\det (1+U(l_{k,\tau }))$$

$$U(l_{k,\tau })=\prod _{\tau =1}^{L_{\tau }}{e}^{-\Delta \tau T}\prod _k{e}^{\sqrt{-\Delta \tau U_k}\eta (l_{k,\tau })V^{(k)}}.$$

$$\begin{array}{rl}\det & (1+U)\\ & =1+\sum _{N=1}^{N_s}\sum _{x_N>x_{N-1}>\cdots >x_1}\det \left[\begin{array}{ccc}{U}_{x_1,x_1}& \dots & U_{x_1,x_N}\\ ⋮& \ddots & ⋮\\ U_{x_N,x_1}& \cdots & U_{x_N,x_N}\end{array}\right]\\ & =1+\sum _x{U}_{x,x}+\sum _{x_2>x_1}\det \left[\begin{array}{cc}{U}_{x_1,x_1}& U_{x_1,x_2}\\ U_{x_2,x_1}& U_{x_2,x_2}\end{array}\right]+\cdots \end{array}$$

$$Z_{\mathrm{can}}(N)=\frac{d^N}{d{z}^N}\sum _{\left\{l_{k,\tau }\right\}}{e^{S_0\left\{l_{k,\tau }\right\}}\det (1+zU(l_{k,\tau }))|}_{z=0}.$$

$$\begin{array}{rl}& Z_{\mathrm{can}}(N)\\ & =\frac{1}{N_s}\sum _{m=1}^{N_s}\sum _{\left\{l_{k,\tau }\right\}}{e}^{S_0\left\{l_{k,\tau }\right\}}{e}^{-i{\varphi }_{m}N}\det (1+e^{i{\varphi }_m}U(l_{k,\tau })),\end{array}$$

$$\begin{array}{rl}{Z}_{\mathrm{can}}(N)& =\mathrm{Tr}\left[{\delta }_{\hat{N},N}{e}^{-\beta \hat{H}}\right]\\ & =\frac{1}{N_s}\sum _{m=1}^{N_s}{e}^{-i{\varphi }_{m}N}\mathrm{Tr}\left[e^{i{\varphi }_m\hat{N}}{e}^{-\beta \hat{H}}\right].\end{array}$$

$$\begin{array}{rl}& \hat{H}(\lambda )=\hat{H}+{\hat{H}}_{\lambda },\\ & {\hat{H}}_{\lambda }\equiv \lambda {\left(\hat{N}-N_0\right)}^2,\end{array}$$

$$Z_{\mathrm{can}}(N_0)=\underset{\lambda \to \infty }{\lim }\mathrm{Tr}\left[e^{-\beta \hat{H}(\lambda )}\right].$$

$$e^{-\beta \hat{H}}=\prod _{\tau =1}^{L_{\tau }}\left[e^{-\Delta \tau \hat{T}}{e}^{-\Delta \tau \hat{V}}\underset{n_{\lambda }\text{-times }}{\underbrace{e^{-\frac{\Delta \tau }{n_{\lambda }}{\hat{H}}_{\lambda }}\cdots e^{-\frac{\Delta \tau }{n_{\lambda }}{\hat{H}}_{\lambda }}}}\right].$$

$${\chi }_c\equiv \frac{\beta }{V}\left(⟨{\hat{N}}^2⟩-⟨\hat{N}{⟩}^2\right).$$

$$\begin{array}{rl}{F}_{\mathrm{can}}& (n_0,V)-F_{\mathrm{gc}}(V)\\ & =\frac{1}{2V}\ln \left(2\pi V\right)+\frac{1}{2V}\ln \left(\frac{{\chi }_c}{\beta }\right)+O\left(\frac{1}{V^2}\right),\end{array}$$

$$e^{-V{F}_{\mathrm{can}}(n,V)}=\int [\mathcal{D}\Psi ]e^{-S[\Psi ]}\delta \left(n,\frac{1}{V}\hat{N}(\Psi )\right),$$

$$e^{-V{F}_{\mathrm{gc}}(V)}=\sum _{n=n_{\min }}^{n_{\max }}{e}^{-V{F}_{\mathrm{can}}(n,V)},$$

$$\begin{split}&e^{-VF_{\rm gc}(V)}=V\Bigg{[}\int_{n_{\rm min}}^{n_{\rm max}}dn\ e^{-VF_{\rm can}(n,V)}\\ &+\frac{e^{-VF_{\rm can}(n_{\rm min},V)}+e^{-VF_{\rm can}(n_{\rm max},V)}}{2V}+O\left(\frac{e^{-cV}}{V}\right)\Bigg{]},\end{split}$$

$$\begin{array}{rl}& e^{-V{F}_{\mathrm{gc}}(V)}/V=e^{-V{F}_{\mathrm{can}}(n_0,V)}\cdot \\ & {\left[\frac{2\pi }{V({\partial }^2{F}_{\mathrm{can}}/\partial n^2)(n_0,V)}\right]}^{1/2}\left[1+O\left(\frac{1}{V}\right)\right]\\ & +\frac{e^{-V{F}_{\mathrm{can}}(n_{\min },V)}+e^{-V{F}_{\mathrm{can}}(n_{\max },V)}}{2V}+O\left(\frac{e^{-cV}}{V}\right),\end{array}$$

$$\begin{array}{rl}{F}_{\mathrm{gc}}(V)& =F_{\mathrm{can}}(n_0,V)-\frac{1}{V}\ln V-\frac{1}{2V}\ln \left(\frac{2\pi }{V}\right)\\ & +\frac{1}{2V}\ln \left[\frac{{\partial }^2{F}_{\mathrm{can}}}{\partial n^2}(n_0,V)\right]+O\left(\frac{1}{V^2}\right),\end{array}$$

$${\chi }_c=\frac{\sum _{n=n_{\min }}^{n_{\max }}\frac{\beta }{V}(nV-n_{0}V{)}^2{e}^{-V{F}_{\mathrm{can}}(n,V)}}{\sum _{n=n_{\min }}^{n_{\max }}{e}^{-V{F}_{\mathrm{can}}(n,V)}}.$$

$${\chi }_c\underset{V\to \infty }{=}\frac{\beta }{\left({\partial }^2{F}_{\mathrm{can}}/\partial n^2\right)(n_0,V)}.$$

$$\begin{array}{rl}{F}_{\mathrm{can}}& (n_0,V)-F(V=\infty )\\ & =\frac{1}{2V}\ln \left(2\pi V\right)+\frac{1}{2V}\ln \left(\frac{{\chi }_c}{\beta }\right)+O\left(\frac{1}{V^2}\right).\end{array}$$

$$E_{\mathrm{can}}(V)-E(V=\infty )=\frac{\partial \left({\chi }_c/\beta \right)/\partial \beta }{2V\left({\chi }_c/\beta \right)}.$$

Full text

I Introduction

One of the central tenet of statistical mechanics is the notion of statistical ensembles. In thermal equilibrium, a system can be described by different statistical ensembles: the microcanonical, the canonical, and the grand-canonical ensemble. In the thermodynamic limit, and in the presence of short-ranged interactions, bulk properties do not generically depend on the choice of the ensemble. Such a property is known as ensemble equivalence.111Nevertheless, systems with long-ranged interactions exhibit violation of the ensemble equivalence Barré et al. (2001); Leyvraz and Ruffo (2002). In particular, a textbook argument for the equivalence between the canonical and the grand-canonical ensembles consists in the following observation. In the grand-canonical ensemble the particle number as well as the energy are sharp in the thermodynamic limit, i.e., their relative fluctuation vanishes in the thermodynamic limit. This stems from the fact that the specific heat,

[TABLE]

and the charge susceptibility

[TABLE]

are extensive quantities that measure energy and particle-number fluctuations. In Eqs. (1) and (2) $\hat{H}$ is the Hamiltonian of the system, $\hat{N}$ the particle-number operator, $\beta=1/k_{B}T$ the inverse temperature in units of the Boltzmann constant $k_{B}$, and $\mu$ the chemical potential. Thus,

[TABLE]

and the selection of the ensemble is merely a matter of convenience. Nonetheless, in many cases the choice of ensemble is dictated by the physical properties of the system under study. In fact, while the canonical ensemble requires the presence of a heat bath which fixes the temperature, the grand-canonical ensemble additionally needs a particle reservoir which allows to fix the chemical potential. Systems which lack such a particle bath, like those found in nuclear physics or in cold atoms, require a description in terms of the canonical ensemble. Moreover, in the case of mesoscopic systems with a finite particle number, a reliable comparison with experimental data needs a theoretical computation based on the canonical ensemble. In this context, we mention that, unlike the finite-temperature auxiliary field quantum Monte Carlo (QMC) method considered here, the so-called projective auxiliary field QMC, which targets the ground state of fermionic models, is a method which is intrinsically formulated within the canonical ensemble Assaad and Evertz (2008).

The aim of this paper is twofold. On one hand we will introduce a QMC method for fermionic models in the canonical ensemble, consisting in a simple formulation of the auxiliary field QMC method which enforces the conservation of the particle number. Our approach differs from that adopted in Ref. Koonin et al. (1997) and supplements the Hamiltonian that we simulate in the grand-canonical ensemble by the long-ranged interaction term

[TABLE]

such that in the infinite-$\lambda$ limit charge fluctuations are suppressed and the canonical ensemble is recovered. This type of interaction is easily incorporated in the auxiliary field QMC, especially in the formulation provided in Ref. Bercx et al. (2017). The advantage of such an approach is that $\lambda$ can be dynamically chosen. For instance, at low temperatures the charge susceptibility can vanish due to finite size or correlation-induced charge gaps. In this case $\lambda$ can be set to a very small number, or even to zero since both canonical and grand-canonical ensembles yield identical results. At high temperatures, where the grand-canonical ensemble exhibits significantly large charge fluctuations, bigger values of $\lambda$ are required to impose the constraint.

The second motivation of the paper is to look into finite-size corrections both in the canonical and grand-canonical ensembles, which we study in quantum and classical lattice models.

Concerning classical models on a lattice, it should be noted that in the literature the canonical ensemble is often defined by the usual partition function sum, where one considers all the configurations without any constraint. In the case of the standard Ising model, this corresponds to the usual partition function:

[TABLE]

However, through the mapping to the lattice gas, the magnetization of the model corresponds to the particle number, which in the ensemble of Eq. (5) is allowed to fluctuate. In order to provide a more meaningful comparison to quantum models, we refer to the lattice gas language and define the grand-canonical ensemble as the one where the magnetization is not fixed; in Eq. (5) we have anticipated this definition, such that the subscript gc refers the grand-canonical ensemble. Conversely, we define the canonical ensemble as the ensemble where the magnetization is fixed, so that the corresponding partition function of the Ising model is

[TABLE]

where the constraint is enforced by employing the Kronecker delta function $\delta(m,n)$. In the literature, the ensemble of Eq. (6) is often referred to as the fixed-magnetization ensemble. In three dimensions, the Ising model at fixed magnetization has been investigated by means of Monte Carlo (MC) simulations in Ref. Blöte et al. (2000), using the geometric cluster algorithm Heringa and Blöte (1998a, b).

In this work, we study the approach to the thermodynamic limit in the presence of a finite mass gap or, in the language of statistical physics, with a finite exponential correlation length. Generically for short-ranged Hamiltonians, on a finite volume with periodic boundary conditions and in the grand-canonical ensemble, the various observables are expected to show a finite-size correction which is proportional to $\exp(-L/\xi)$, where $L$ is the linear size of the system and $\xi$ is the exponential correlation length (or inverse mass gap). This expectation has been confirmed by explicit field theory calculations, both in the continuum Lüscher (1986); Neuberger (1989) and on a lattice Münster (1985); early numerical studies confirmed these prediction Montvay and Weisz (1987). An exponential approach to the thermodynamic limit is also verified, e.g., in the well-known solution of the one-dimensional Ising model, as well as in generic one-dimensional $O(N)-$invariant spin models Cucchieri et al. (1997). Nevertheless, it should be noted that, in the grand-canonical ensemble, some specific observables can exhibit a leading finite-size correction proportional to a power law of the system size. This is the case of the most common definition of the second-moment correlation length on a lattice, where finite-size corrections $\propto 1/L^{2}$ are due to the discretization of momenta on a finite lattice; see, e.g., the corresponding discussion in Ref. Caracciolo and Pelissetto (1998) and Appendix A of Ref. Parisen Toldin et al. (2015). We also remark that, in the presence of nontranslationally invariant boundary conditions, finite-size corrections polynomial in the inverse lattice size $1/L$ arise naturally, being related to subleading terms in the free energy; for instance, open boundary conditions result in the presence of a surface free energy which is depressed by a factor $1/L$ with respect to the bulk one and gives rise to finite-size corrections $\propto 1/L$ for bulk observables.

Conversely, in the canonical ensemble the prediction of exponentially decaying finite-size corrections fails, since the constraint introduces a long-ranged interaction, such that fluctuations in spatially separated regions (as measured by the correlation length) are not independent. Such a long-ranged (weak) correlation modifies also the high-temperature expansion of a model Iyer et al. (2015) and results in a slower approach to the thermodynamic limit of various observables, so that the leading finite-size correction is proportional to the inverse volume $V$. Several important properties of the free energy in the canonical ensemble have been, in fact, discussed in the literature, although under a different perspective and notation. In quantum field theory, the so-called constrained effective potential $U_{\rm eff}$, introduced in the context of scalar field theories in Ref. Fukuda and Kyriakopoulos (1975), is defined as

[TABLE]

where $S[\varphi]$ is the action of the theory and the right-hand side of Eq. (7) is a constrained path-integral over the field configurations where the volume-average value of $\varphi$ is fixed to $m$. In the language of statistical physics, the right-hand side of Eq. (7) is precisely a constrained partition function sum at fixed magnetization, i.e., the partition function in the canonical ensemble. Hence, $U_{\rm eff}(m,V)$ is the free energy per volume and $k_{B}T$ in the canonical ensemble. A detailed analysis of the constrained effective potential has shown that it admits an infinite-volume limit $U_{\rm eff}(m,V\rightarrow\infty)$ which coincides with the usual effective potential $\Gamma(m)$ of the theory O’Raifeartaigh et al. (1986). Moreover, as argued in Ref. Neuberger (1989), $U_{\rm eff}(m,V)$ exhibits finite-size corrections which are polynomial in $1/V$. This is because, as a consequence of the definition in Eq. (7), the grand-canonical average of any function of the magnetization $m$ is equivalent to an average over an effective probability measure $\propto\exp\{-VU_{\rm eff}(m,V)\}$, which for $V\rightarrow\infty$ can be evaluated by a saddle-point expansion, resulting in a series in $1/V$. On the other hand, the grand-canonical average converges exponentially to the limit $V\rightarrow\infty$. This is possible only if $U_{\rm eff}(m,V)$ displays finite-size corrections polynomial in $1/V$, which exactly cancel the expansion in $1/V$ originating from the saddle-point evaluation Neuberger (1989). A renormalized loop expansion for a $\phi^{4}$ theory on the lattice has confirmed the existence of finite-size corrections $\propto 1/V$ Palma (1992).

In this context, a recent study verified the existence of finite-size corrections $\propto 1/V$ in the canonical ensemble, and, conversely, of exponentially decaying finite-size corrections in the grand-canonical ensemble Iyer et al. (2015). In this paper we provide an exact formula for the leading finite-size corrections in the canonical ensemble of the free energy density, energy density and other observables. While our analysis is restricted to the case of a finite correlation length, we mention that the introduction of a constraint to a nonordering parameter results in the so-called Fisher renormalization, leading to a modification of the singularities associated with a critical point, such that the critical exponents differ from those observed in the corresponding unconstrained system Essam and Garelick (1967); Fisher (1968). The choice of ensemble is also relevant to the so-called critical Casimir force Fisher and de Gennes (1978), whose behavior in the canonical ensemble has been recently investigated within mean-field theory and MC simulations Gross et al. (2016).

This paper is organized as follows. In Sec. II we illustrate the QMC method that we use to generate numerical data for fermionic models in the canonical ensemble. In Sec. III we provide an exact determination of the leading finite-size corrections in the canonical ensemble. In Sec. IV we study the finite-size corrections of the Hubbard model in one and two dimensions. In Sec. V we summarize our results. In the Appendix we provide an exact solution of the one-dimensional classical Ising model in the canonical ensemble to the leading order in $1/V$, which confirms the general result of Sec. III.

II Canonical auxiliary field methods

II.1 General formulation

In this section we review various methods to achieve canonical auxiliary field QMC simulations at finite temperature. We will consider a Hamiltonian of the form

[TABLE]

that can be readily implemented in the ALF package Bercx et al. (2017). Here $x$ is a super-index encoding orbital and spin degrees of freedom, $\hat{c}^{\dagger}_{x}$ are fermion creation operators, $V^{(k)}$ and $T$ are Hermitian matrices, and $U_{k}$, $\alpha_{k}$ real numbers. To simplify the notation, in the following we assume that the chemical potential term $\mu\hat{N}$, with $\hat{N}\equiv\sum_{x}\hat{c}^{\dagger}_{x}\hat{c}^{\phantom{\dagger}}_{x}$, has been adsorbed into the Hamiltonian $\hat{H}$. Using the Trotter decomposition with $L_{\tau}\Delta\tau=\beta$, and a discrete version of the Hubbard-Stratonovich (HS) transformation,

[TABLE]

with $\gamma(\pm 1)=1+\sqrt{6}/3$, $\eta(\pm 1)=\pm\sqrt{\smash[b]{2(3-\sqrt{6})}}$, and $\gamma(\pm 2)=1-\sqrt{6}/3$, $\eta(\pm 2)=\pm\sqrt{\smash[b]{2(3+\sqrt{6})}}$, one can approximate the imaginary time propagator $e^{-\beta\hat{H}}$ as

[TABLE]

Here $S_{0}\left\{l_{k,\tau}\right\}=\sum_{l_{k,\tau}}\ln\left(\gamma(l_{k,\tau})\right)+\sqrt{-\Delta\tau U_{k}}\eta(l_{k,\tau})\alpha_{k}$ . It is easy to show that the contribution of the anti-Hermitian operator $i\hat{R}$ to the expectation value of an Hermitian oberservable is purely imaginary, so that the discretization error $\propto\Delta\tau$ can be filtered out, leading to a Trotter error $\propto\Delta\tau^{2}$. The systematic error involved in this discrete HS transformation is of a higher order than the one encountered in the Trotter decomposition so that it can be regarded as good as exact. At this point, one can integrate out the fermions so as to obtain the grand-canonical partition function:

[TABLE]

with

[TABLE]

Using the Leibniz formula for determinants, one can show that:

[TABLE]

Here $N_{s}$ corresponds to the number of single-particle states, and one can readily see that each term of the sum corresponds to the canonical trace of $N$ single-particle states. Thereby, the canonical partition function $Z_{\rm can}(N)$ is given by:

[TABLE]

A numerical implementation of the above equation reads:

[TABLE]

where $\phi_{m}=2\pi m/N_{s}$. An equivalent way to show the above result is to note that the total particle number $\hat{N}$ commutes with the Hamiltonian such that:

[TABLE]

By applying a Trotter decomposition and HS transformation to the right-hand side of Eq. (16), one can reproduce Eq. (15). Implementations of canonical simulations using the above results have been proposed in Refs. Ormand et al. (1994); Gilbreth and Alhassid (2015). In these approaches, the discrete Fourier transformation is computed exactly at each MC step. For the method to be successful, the chemical potential has to be chosen such that the average particle number is peaked around the desired value.

II.2 Constraint of the particle-number fluctuations

Here we follow a slightly different approach and modify the Hamiltonian as

[TABLE]

such that

[TABLE]

As discussed above, in Eqs. (17) and (18) the Hamiltonian $\hat{H}$ implicitly depends on the chemical potential $\mu$, which needs to be tuned such that $\langle\hat{N}\rangle=N_{0}$. In practice, this is done by computing $\langle\hat{N}\rangle$ as a function of $\mu$, for a suitable interval in $\mu$, by means of auxiliary field QMC and then fixing $\mu$ such that the equation $\langle\hat{N}\rangle=N_{0}$ is satisfied within the desired statistical accuracy; at half-filling one has exactly $\mu=0$.

Since $\hat{H}$ conserves the particle number, one can foresee rapid convergence in $\lambda$ because particle-number sectors with $\hat{N}=N\neq N_{0}$ have a statistical weight suppressed by a factor $e^{-\lambda\beta\left(N-N_{0}\right)^{2}}$. The latter also shows that the relevant parameter for the convergence is $\beta\lambda$ rather than $\lambda$ itself. The additional term is a perfect square term which is easily implemented within the ALF code Bercx et al. (2017). Since $\left(\hat{N}-N_{0}\right)^{2}$ effectively corresponds to a long-ranged interaction, one may face the issue that the acceptance rate of a single HS flip becomes excessively small on large lattices. To circumvent this problem we have used the following decomposition:

[TABLE]

Thereby, we need $n_{\lambda}$ fields per time slice to impose the constraint. For each field, the coupling constant is effectively suppressed by a factor $n_{\lambda}$, thus allowing to control the acceptance of the QMC algorithm.

In order to test the efficiency of our QMC method in the canonical ensemble, we computed the uniform intensive charge susceptibility $\chi_{c}$, defined as

[TABLE]

Note that compared with the extensive definition in Eq. (2), here the susceptibility is divided by the system volume $V$. In Fig. 1 we show $\chi_{c}$ for the one-dimensional (1D) Hubbard model as a function of $\beta\lambda$ and $n_{\lambda}$. As shown in Fig. 1(a), $\chi_{c}$ decays gradually from a finite value to zero on increasing $\beta\lambda$. The threshold in $\lambda$ for which $\chi_{c}$ converges to zero corresponds to the canonical ensemble. A comparison of the results for lattice sizes $L=4$, $8$, and $16$ suggests that the charge fluctuations are easier to suppress for larger system sizes. The dependence of $\chi_{c}$ on $n_{\lambda}$ defined in Eq. (19) is shown in Fig. 1(b), which illustrates the increased Trotter error for larger values of $\beta\lambda$.

Figure 2 shows the decay of charge susceptibility $\chi_{c}$ as a function of $\lambda$ in the two-dimensional (2D) Hubbard model, for $U=4.0$, $L=4$ and several inverse temperatures $\beta=0.5$, $2.0$, and $5.0$. Inspection of Fig. 2 reveals that in the grand-canonical ensemble the $\beta=2.0$ case exhibits charge fluctuations larger than the $\beta=0.5$ case, thereby requiring a larger value of $\beta\lambda$ to realize the canonical ensemble.

III Finite-size corrections in the canonical ensemble: exact results

In this section, by exploiting the relation between the canonical and the grand-canonical free energy, we determine the leading finite-size correction of the free energy in the canonical ensemble and relate it to a susceptibility. To be concrete, we consider a quantum model on a lattice, where in the canonical ensemble the number of particles is fixed, and we prove that on a finite volume $V$

[TABLE]

where $F_{\rm can}(n_{0},V)$ and $F_{\rm gc}(V)$ are the free energies per volume $V$ and in units of $k_{B}T$ in the canonical and grand-canonical ensembles, respectively, and $\chi_{c}$ is the charge susceptibility (in the grand-canonical ensemble), defined in Eq. (20); the filling fraction $n_{0}$ in $F_{\rm can}(n_{0},V)$ is fixed to the corresponding expectation value in the grand-canonical ensemble. Equation (21) provides the leading additional contribution to the free energy density due to the particle-number constraint. As discussed towards the end of this section, Eq. (21) allows also to determine the leading finite-size correction of observables in the canonical ensemble if, as expected, finite-size corrections in the grand-canonical ensemble decay faster than $1/V$.

In order to prove Eq. (21), we observe that the free energy density $F_{\rm can}(n,V)$ can be related to a path-integral formulation of the canonical partition function as

[TABLE]

where $\Psi$ indicates collectively the fields entering in the path integral, $S[\Psi]$ is the action of the model, $\hat{N}(\Psi)$ is the expression of the total number operator $\hat{N}$ in terms of the fields $\Psi$, and $n$ is the intensive filling fraction, which is fixed in the canonical ensemble. In Eq. (22), $S[\Psi]$, as well as $F_{\rm can}(n,V)$, additionally depend on the temperature and coupling constants, inessential for the present discussion. On a lattice, $\hat{N}$ is the sum of single-site and single-species number operators $\hat{N}(x)$, $\hat{N}=\sum_{x}\hat{N}(x)$, therefore $n$ can only take discrete values, separated by an interval of $1/V$. By summing over the allowed values of $n$, we obtain the grand-canonical free energy density $F_{\rm gc}(V)$

[TABLE]

where, as before, we have ignored the dependence of $F_{\rm gc}(V)$ on the various coupling constants, and $n_{\rm min}$, $n_{\rm max}$ indicate the minimum and maximum number of particles per volume that the model can host; usually $n_{\rm min}=0$, while $n_{\rm max}$ depend on the number and type of particle species. For $V\rightarrow\infty$ the sum in Eq. (23) can be approximated by the Euler-Maclaurin formula as

[TABLE]

where the next-to-leading term in the Euler-Maclaurin formula is $\propto(1/V^{2})\partial(e^{-VF_{\rm can}})/\partial n$ computed at the end points, hence it is of order $e^{-cV}/V$. In the limit $V\rightarrow\infty$, the integral on the right-hand side of Eq. (24) is dominated by the minimum $n_{0}$ of $F_{\rm can}(n,V)$. If $n_{0}$ is an interior point222The case of multiple saddle points, or a saddle point at an end point requires a separate analysis. of the integration interval $[n_{\rm min},n_{\rm max}]$, by using the saddle-point method we obtain

[TABLE]

where the factor $1+O(1/V)$ represents the next-to-leading term in the saddle-point expansion. The second term on the right-hand side of Eq. (24) is depressed by a factor $\propto\exp\{-V[F_{\rm can}(n_{\rm min},V)-F_{\rm can}(n_{0},V)]\}/V^{1/2}+\exp\{-V[F_{\rm can}(n_{\rm max},V)-F_{\rm can}(n_{0},V)]\}/V^{1/2}$ with respect to the first term, therefore, since $n_{0}$ is the minimum of $F_{\rm can}(n,V)$, it is subleading with respect to the first factor. Moreover, the convergence of the integral in Eq. (24) requires the last term on the right-hand side of Eq. (25) to be subleading with respect to the first factor. Thus, by factorizing the first term on the right-hand side of Eq. (25) and taking the logarithm, the last two terms give a contribution of order $\ln(1+\exp\{-cV\}/V^{1/2})\sim\exp\{-cV\}/V^{1/2}$, which is negligible with respect to the correction of order $1/V$ originating from the next-to-leading term of the saddle-point expansion. On taking the logarithm on both sides of Eq. (25) we find

[TABLE]

where subleading exponential corrections have been neglected. The second and third terms $\propto\ln V$ on the right-hand side of Eq. (26) represent an entropic contribution which is due to the larger configurational space of the grand-canonical ensemble as compared to the canonical one. In particular, the first constant originates from the discretization of the allowed values of $n$ [see the discussion after Eq. (22)] and is absent in continuous models. The saddle-point position $n_{0}$ appearing in the previous equations corresponds precisely to the grand-canonical expectation value of $\langle\hat{N}/V\rangle_{\rm gc}$. This is because, using Eq. (22) and Eq. (23), one can write $\langle\hat{N}/V\rangle_{\rm gc}e^{-VF_{\rm gc}(V)}=\sum_{n=n_{\rm min}}^{n_{\rm max}}ne^{-VF_{\rm can}(n,V)}$. Along the same line of reasoning as above, one finds that, as expected also from thermodynamic considerations, $\lim_{V\rightarrow\infty}\langle\hat{N}/V\rangle_{\rm gc}=n_{0}$. Thus, the quantity $F_{\rm can}(n_{0},V)$ on the right-hand side of Eq. (26) is precisely the free energy density with a particle number fixed to its expectation value in the grand-canonical ensemble, i.e., the thermodynamic quantity which is meaningful to compare with the grand-canonical free energy density. The fluctuation of the particle number, which determines the charge susceptibility $\chi_{c}$ defined in Eq. (20), can be related to the finite-size correction on the right-hand side of Eq. (26). By using Eq. (22), Eq. (23), and the definition of Eq. (20), one obtains

[TABLE]

The right-hand side of Eq. (27) can be evaluated for $V\rightarrow\infty$ using a saddle-point expansion as above, resulting in

[TABLE]

Finally, inserting Eq. (28) into Eq. (26), we obtain Eq. (21).

If finite-size corrections of $F_{\rm gc}(V)$ decay faster than $1/V$ (indeed, as discussed in Sec. I, we expect exponentially decaying finite-size corrections), we can replace $F_{\rm gc}(V)$ on the left-hand side of Eq. (21) with its ensemble-independent thermodynamic limit $F(V=\infty)=F_{\rm gc}(V=\infty)=F_{\rm can}(n_{0},V=\infty)$, such that the leading finite-size corrections in $F_{\rm can}(n_{0},V)$ are

[TABLE]

From Eq. (29) we can, e.g., determine the leading finite-size correction of the energy density in the canonical ensemble by taking the derivative with respect to $\beta$:

[TABLE]

It is useful to remark that the charge susceptibility $\chi_{c}$ appearing in Eqs. (21) and (27)-(30) is computed in the grand-canonical ensemble. Since $\chi_{c}$ has a finite thermodynamic limit and exponentially decaying finite-size corrections, it does not give rise to a further algebraic volume dependence. Equation  (30) can be generalized to other local observables and correlations thereof. To this end, one can supplement the action of the model with an external source term

[TABLE]

where the sum extends to the lattice sites and $O(x,\tau)$ is a local observable, to be expressed in terms of the fields $\Psi$ entering in the path integral of Eq. (22). Such an addition corresponds to the insertion of external lines in the Feynman diagram expansion, and hence one expects, in line with the analysis of Ref. Neuberger (1989), that in the presence of a finite mass gap correlations including $O(x,\tau)$ are characterized by exponentially decaying finite-size corrections. Under the substitution of Eq. (31), the charge susceptibility $\chi_{c}$ entering in Eq. (21) and Eq. (29) acquires a dependence on the external field $h$. Differentations of the free energy density with respect to $h$ provide the analogous of Eq. (30) for the finite-size corrections of the volume-average and susceptibility of $O$:

[TABLE]

[TABLE]

where we emphasize that the expectation values of $O(x,\tau)$ are computed in absence of the external field $h$. In particular, in a spinful model Eq. (33) implies a leading finite-size correction $\propto 1/V$ of the spin susceptibility in the canonical ensemble. We remark that the derivatives of $\chi_{c}(h)$ appearing in Eqs. (32) and (33) can be in principle directly computed by sampling a suitable observable, thus avoiding a numerical differentiation.

The results of Eqs. (21), (29), (30), (32), and (33) can be easily generalized to other correlations by considering a considering a space- and imaginary time-dependent source $h(x,\tau)$ in Eq. (31), or to other types of constrained models, along the same line of reasoning.

IV Fermionic simulations in the canonical ensemble

We performed QMC simulation of the $SU(2)$ Hubbard model in both the grand-canonical and canonical ensemble. The Hamiltonian of the Hubbard model is defined as:

[TABLE]

where $\hat{n}_{\bm{i},\sigma}\equiv\hat{c}^{\dagger}_{\bm{i},\sigma}\hat{c}_{\bm{i},\sigma}$. The canonical ensemble is realized by adding the constraint given in Eq. (17). For such a modified Hamiltonian, the total number of particles converges quickly to $N_{0}$ on increasing $\beta\lambda$.

Here we simulated both ensembles on a 1D lattice, as well as on the 2D square lattice at finite temperature, both of which are known to be disordered. We mainly considered the models at half filling ($N_{0}=N_{s}/2$, with $N_{s}=2L^{d}$) with zero chemical potential $\mu=0$ and carried out some test calculations for the two-dimensional doped Hubbard model. In all simulations we fixed $t=1$ and $U=4.0$. Our basic MC observables are as follows:

Energy density333Up to an inessential, filling-dependent, additive constant, $E$ corresponds to the energy part on the right-hand side of Eq (34).:

[TABLE] 2. 2.

Uniform spin susceptibility:

[TABLE]

IV.1 1D model

The QMC simulations of the one-dimensional Hubbard model are performed in both the grand-canonical and canonical ensembles at inverse temperature $\beta=0.5$, system sizes $L=4$, $8$, $12$, $16$, $20$, $24$, $28$, $32$, $40$, $48$, $64$, $72$, $80$, and at half-filling. A comparison of the size effect for the energy density $E(L)$ and for the uniform spin susceptibility $\chi_{s}(L)$ in the two ensembles is shown in Fig. 3 and Fig. 4, respectively. We observe that in the grand-canonical ensemble both $E$ and $\chi_{s}$ converge quickly to the thermodynamic limit for small system sizes. This indicates a small correlation length $\xi$ at this temperature.

On the other hand, except for the smallest system sizes, in the canonical ensemble both observables exhibit a linear-like behavior as a function of $1/L$. A fit of energy density in the canonical ensemble $E_{\rm can}(L)$ to $E_{\rm can}(L)=E(L\rightarrow\infty)+aL^{-1}$ exhibits a good $\chi^{2}/\text{DOF}$ (DOF denotes the number of degrees of freedom), when the data for the small sizes are discarded; the extrapolated value $E(L\rightarrow\infty)$ matches the grand-canonical result. Similar considerations hold for a fit of the spin susceptibility in the canonical ensemble $\chi_{s,{\rm can}}(L)$ to $\chi_{s,{\rm can}}(L)=\chi_{s}(L\rightarrow\infty)+aL^{-1}$.

Moreover, a fit of $E_{\rm can}(L)$ to $E(L\rightarrow\infty)+aL^{-d}$, leaving $d$ as a free parameter, gives $d=1.05(2)$ when the smallest lattice size taken into account for the fit is $L_{\rm min}=16$. An equivalent fit for $\chi_{s}(L)$ gives $d=1.04(2)$, when $L_{\rm min}=12$. This confirms that finite-size corrections of observables in the canonical ensemble are $\propto 1/L$.

IV.2 2D model

We simulated the Hubbard model on the square lattice for both ensembles at $\beta=0.5$ and $\beta=2.0$, lattice sizes $L=4$, $6$, $8$, $10$, $12$, $14$, $16$, and at half-filling.

Figure 5 and Fig. 6 show the size behavior of $E$ and $\chi_{s}$ for the two ensembles at $\beta=0.5$. The observed tiny size dependence of the observables in the grand-canonical ensemble suggests that the correlation length $\xi$ is smaller than the minimum lattice size $L=4$. On the other hand, in the canonical ensemble the energy density $E$ and the spin susceptibility $\chi_{s}$ show a linear-like behavior as function of $1/L^{2}$.

For a more quantitative check of the finite-size correction in the canonical ensemble, we fitted $E_{\rm can}(L)$ to $E_{\rm can}(L)=E(L\rightarrow\infty)+aL^{-1}+bL^{-2}+cL^{-3}$ and $\chi_{s,{\rm can}}(L)$ to an equivalent Ansatz, leaving $a,b$ and $c$ as free parameters. Fit results for both observables show a good $\chi^{2}/\text{DOF}$ when $L_{\rm min}=6$, and the coefficient $a$ vanishes within error bars, whereas $b$ acquires a finite value. On the other hand, a fit of $E_{\rm can}(L)$ to $E_{\rm can}(L)=E(L\rightarrow\infty)+bL^{-d}$, leaving $b$ and $d$ as free parameters, and of $\chi_{s}(L)$ to an equivalent Ansatz, gives $d=2.05(3)$ and $d=1.9(1)$ for $E$ and $\chi_{s}$, respectively, when $L_{\rm min}=6$. In line with the discussions of Sec. III, these fit results confirm that the leading finite-size correction in the canonical ensemble is $\propto 1/L^{2}$.

We also simulated the 2D Hubbard model at a lower temperature $\beta=2.0$. A corresponding comparison of the finite-size energy density for the grand-canonical and canonical ensemble is shown in Fig. 7. Generically, finite-size corrections in the canonical ensemble are expected to be temperature dependent. On the other hand, the exponential correction characterized by the correlation length in the grand-canonical ensemble may start to be relevant at a lower temperature, because of an increased correlation length.

The data shown in Fig. 7 exhibit a visible decay of the energy density in the grand-canonical ensemble $E_{\rm gc}$, on increasing the system size. As a guide to the eye, we fitted $E_{\rm gc}$ to $E_{\rm gc}(L)=E(L\rightarrow\infty)+b\cdot e^{-L/c}$. The finite-size values of $E$ in the canonical ensemble show a nonmonotonic behavior between $L=4$ and $6$, which might be due to a combination of various sources of finite-size corrections, such as the one $\propto 1/V$ originating from the particle-number constraint, the one related to the correlation length, and the residual correction term due to the regular part of the free energy. Nevertheless, a finite-size dependence $\propto 1/L^{2}$ can be clearly observed in Fig. 7 for $L>6$, with a smaller slope compared to the $\beta=0.5$ case (compare with Fig. 5).

We note that the auxiliary field QMC for the grand-canonical ensemble has a mild sign-problem under doping the system away from half filling, provided that the temperature is high enough. Here we also tested the efficiency of the canonical ensemble QMC method under doping. To this end, for every lattice size we tuned the chemical potential $\mu$ such that the expectation value of the number of particles in the grand-canonical ensemble matches the desired number $N_{0}$ of particles in the canonical ensemble. Subsequently, the canonical ensemble is realized by introducing a Lagrange multiplier, as discussed in Sec. II.2. In order to test the sign-problem, we also calculated the average sign during the MC simulation:

[TABLE]

where $S(C)$ is the action for the MC configuration $C$, so that the corresponding statistical weight is $\propto e^{-S(C)}$. The sign is not necessarily a real number (when the MC is sign-problem free, $S(C)$ is real and $\langle\text{sign}\rangle=1$). On the other hand, the expectation value of observables can be computed via a reweighting scheme only when $\langle\text{sign}\rangle$ is not too small. Figure 8 shows the average sign during MC of a doped 2D Hubbard model at $\frac{3}{8}$ filling, which in the grand-canonical ensemble system does not exhibit a significant sign-problem at an inverse temperature of $\beta=2.0$. Quite remarkably, when the number of particles is fixed in the canonical ensemble, the code still exhibits an average sign higher than $0.97$, for system sizes up to $L=20$. This confirms the feasibility of our QMC method for the canonical ensemble, even under doping. Figure 9 shows the finite-size behavior of energy density of the two ensembles at $\frac{3}{8}$ filling. Similar to the results at half-filling, and in line with the analysis of Sec. III, the energy density in the grand-canonical ensemble exhibits very small finite-size corrections when $L\geq 10$, whereas in the canonical ensemble we observe a finite-size correction approximately linear in $1/L^{2}$ for $L\geq 12$.

V Summary

In this paper we have introduced a method to simulate fermionic models in the canonical ensemble. It consists in an auxiliary field QMC simulation, where the Hamiltonian is supplemented by an additional Lagrange multiplier, which constraints the particle number. The method can implemented using the ALF package for fermionic simulations Bercx et al. (2017). In general, we find that canonical simulations are more computationally demanding than the corresponding ones in the grand-canonical ensemble. Although in the presence of short-ranged interactions the grand-canonical and the canonical ensemble are equivalent in the thermodynamic limit, their approach to the infinite-volume limit is distinctively different. In the canonical ensemble the observables are generically found to display a finite-size correction which is proportional to the inverse volume. In Sec. III we prove an exact formula for the leading finite-size correction of the free energy density, the energy density, and other observables. Such a correction is controlled by the charge susceptibility and is found to be proportional to the inverse volume. This result is further substantiated by an exact calculation for the one-dimensional Ising model reported in the Appendix. Our numerical simulations of the Hubbard model reported in Sec. IV confirm the presence of finite-size corrections proportional to the inverse volume in the canonical ensemble. In line with previous theoretical results, in the presence of a finite correlation length and for periodic boundary conditions, observables computed in the grand-canonical ensemble display a faster approach to the thermodynamic limit, such that the leading finite-size correction is exponential in the ratio of the linear size over the correlation length.

Note added: After completing this paper we became aware of related research presented in Ref. Gross et al. (2017), which investigates the effect of a constraint within statistical field theory.

Acknowledgements.

FFA acknowledges useful conversations with A. Sandvik. FPT acknowledges useful communications with M. Gross. This work was supported by the German Research Foundation (DFG) through Grant No. SFB 1170 ToCoTronics and Grant No. FOR 1807. We acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA Krause and Thörnig (2016) at the Jülich Supercomputing Centre.

Appendix: One-dimensional Ising model in the canonical ensemble: Exact results

In this Appendix we compute the leading finite-size correction of the free energy of the one-dimensional Ising model in the canonical ensemble. Although we mainly consider the antiferromagnetic Ising model, the results are also valid for the ferromagnetic model. Employing periodic boundary conditions, the Hamiltonian is

[TABLE]

where $L$ is the number of sites and $S_{L+1}\equiv S_{1}$. The partition function $Z_{\rm can}$ in the canonical ensemble with fixed magnetization $M=0$ is given by

[TABLE]

where we have defined $K\equiv\beta J$.

The constraint $M=0$ can be expressed by using an integral representation of the Kronecker delta function $\delta$ appearing in Eq. (A2), such that

[TABLE]

Inspecting Eq. (A3), we observe that $Z_{\rm can}$ is obtained as the integral over $\mu$ of the partition function for a one-dimensional Ising model in an external imaginary field $i\mu$. The trace over the configuration space can be computed using standard transfer-matrix techniques, yielding

[TABLE]

where the eigenvalues of the transfer matrix are

[TABLE]

and $\lambda_{\pm}(\mu)$ depends implicitly also on $K$. By noting that $\lambda_{\pm}(\mu+\pi)=-\lambda_{\mp}(\mu)$, Eq. (A4) can be cast in the form

[TABLE]

Equation (A6) shows that for $L$ odd the partition function vanishes exactly. This can be readily understood by the impossibility of imposing the constraint $\sum_{i}S_{i}=0$ with an odd number of spin variables $S_{i}$, which take values $\pm 1$. In the following we shall assume that $L$ is even, such that the two terms in Eq. (A6) are identical and we have

[TABLE]

For large $L$ the integral of Eq. (A7) is dominated by the saddle points of $\lambda_{\pm}(\mu)$ which are solutions of

[TABLE]

For both eigenvalues, Eq. (A8) has solutions for $\mu=0$ and $\mu=\pi$, which lie at the border of the integration domain in Eq. (A7). We observe that for $K>0$ (antiferromagnetic model), the eigenvalues are real. For $K<0$ (ferromagnetic model), $\lambda_{\pm}(\mu)$ given in Eq. (A5) are real for $|\mu|<\varepsilon_{0}\equiv\arcsin(e^{2K})$ and $|\mu-\pi|<\varepsilon_{0}$, i.e., in an interval around the saddle points. For this reason, without losing generality it is convenient to shift the domain of integration in Eq. (A7)

[TABLE]

where for $K>0$ one can take, e.g., $\varepsilon_{0}=\pi/2$, such that the single saddle point $\mu=0$ in the integration domain is an interior point. In order to determine the finite-size correction to the free energy, we need to compute the corrections around the saddle point. To this end, it is important to observe that for $K>0$ (antiferromagnetic model) both eigenvalues $\lambda_{\pm}(\mu)$ have a maximum around $\mu=0$; however, since $\lambda_{-}(\mu)<0$, with $L$ even the term $\lambda_{-}(\mu)^{L}$ in Eq. (A9) has a minimum at $\mu=0$, whereas $\lambda_{+}(\mu)>0$ and $\lambda_{+}(\mu)^{L}$ has a maximum. For $K<0$ (ferromagnetic model) $\lambda_{\pm}(\mu)$ are real and positive in an interval around $\mu=0$; however, while $\lambda_{+}(\mu)$ has a maximum at $\mu=0$, the other eigenvalue $\lambda_{-}(\mu)$ has instead a minimum around $\mu=0$. Therefore, for both cases $K>0$ and $K<0$ it is not possible to separate Eq. (A9) into a sum of two integrals to be evaluated for $L\rightarrow\infty$, but it is necessary to consider the behavior around $\mu=0$ of the sum of the two eigenvalues. To do so, we write the integrand of Eq. (A9) as

[TABLE]

A second-order Taylor expansion of $g(\mu,L)$ around $\mu=0$ gives

[TABLE]

Equation (A11) shows that, indeed, $g(\mu,L)$ exhibits a maximum around $\mu=0$. Moreover, the coefficient in front of $\mu^{2}$ remains finite in the limit $L\rightarrow\infty$. Inserting the expansion of Eq. (A11) in Eq. (A10), and using the resulting expression for $\lambda_{+}(\mu)^{L}+\lambda_{-}(\mu)^{L}$ in Eq. (A9) we obtain, after a Gaussian integration,

[TABLE]

The free energy per volume $L$, and in units of $k_{B}T$, $F_{\rm can}=-\ln Z_{\rm can}/L$ is

[TABLE]

Using Eq. (A5) and the known relation between the correlation length $\xi$ and the transfer-matrix eigenvalues

[TABLE]

Eq. (A13) can be written as

[TABLE]

where in the last equality we have expanded for $\xi/L\ll 1$, using the fact that $\xi$ is always finite. Equation (A15) agrees with the general result of Eq. (29), where subleading exponential finite-size corrections have been neglected. To confirm this, we observe that, under the mapping to the lattice gas model, the equivalent charge susceptibility $\chi_{c}$ (i.e., fluctuation of the particle number per volume and multiplied by $\beta$) is given by $\chi_{c}=\chi/4$, where $\chi$ is the usual spin susceptibility which, for the one-dimensional Ising model, is $\chi=\beta\exp\{-2K\}$; on substituting $\chi_{c}\rightarrow\chi/4=\beta\exp\{-2K\}/4$ in Eq. (29) we recover Eq. (A15). Alternatively, one can repeat the calculations of Sec. III, fixing in Eq. (22) the magnetization per volume $m$ instead of the filling fraction $n$. Then in the result of Eq. (29), $\chi_{c}$ is replaced by the fluctuations of the magnetization, i.e., the usual spin susceptibility $\chi$. Moreover, different than for $n$, the allowed values for $m$ are separated by an interval of $2/V$. This results in a factor $1/2$ in front of the right-hand side of Eq. (24), which in turns gives rise to an additional contribution $-(\ln 2)/V$ to the right-hand side of Eq. (29). Taking into account this additional term, and substituting $\chi_{c}\rightarrow\chi$ in Eq. (29) we recover Eq. (A15).

A comparison of Eq. (A15) with the corresponding result for the grand-canonical ensemble

[TABLE]

shows that, besides an irrelevant $L-$dependent constant, the free energy density in the canonical ensemble is affected by a finite $1/L$ correction to its thermodynamic limit $-\ln\left[2\cosh(K)\right]$, which is absent in the grand-canonical ensemble. We also notice that the constraint $M=0$ alters the coefficient in front of the subleading exponential correction $\exp\{-L/\xi\}$. From Eq. (A15) we can compute the energy density as

[TABLE]

which exhibits a leading finite-size correction $\propto 1/L$. Due to the fact that $\chi/\beta$ is exactly exponential in the one-dimensional Ising model, such a finite-size term is temperature independent [see Eq. (29)].

As emphasized in the derivation of the results, Eq. (A15) and Eq. (A17) are also valid for $J<0$, $K<0$, i.e., for a ferromagnetic model. We remark that it is not possible to take the limit $T\rightarrow 0$ in Eq. (A15) and Eq. (A17) because the calculation assumes a finite correlation length $\xi$. Indeed, for $T\rightarrow 0$ the coefficient of $\mu^{2}$ in Eq. (A11) either vanishes (for the antiferromagnetic model) or diverges in $L$ (for the ferromagnetic model), rendering the saddle-point expansion singular. In the ground state of the antiferromagnetic model $E_{\rm can}=-J$, with no size dependence.