Efficient Numerical Evaluation of Thermodynamic Quantities on Infinite (Semi-)classical Chains

Abstract

This work presents an efficient numerical method to evaluate the free energy density and associated thermodynamic quantities of (quasi) one-dimensional classical systems, by combining the transfer operator approach with a numerical discretization of integral kernels using quadrature rules. For analytic kernels, the technique exhibits exponential convergence in the number of quadrature points. As demonstration, we apply the method to a classical particle chain, to the semiclassical nonlinear Schrödinger (NLS) equation and to a classical system on a cylindrical lattice. A comparison with molecular dynamics simulations performed for the NLS model shows very good agreement.

1 Introduction

The partition function Z is a cornerstone of statistical mechanics [1, 2], in particular since thermodynamics quantities can be obtained as derivatives of Z. Nevertheless, computing Z for large systems is a challenging task if a closed-form solution cannot be found. In this work we are concerned with an efficient numerical method for evaluating Z (corresponding to the Boltzmann-Gibbs ensemble) and the associated free energy density for (semi-)classical chains, in the thermodynamic limit of infinite system size. In more detail, we consider a Hamiltonian H on a (quasi) one-dimensional lattice of size L, containing solely nearest neighbor terms:

$$\begin{aligned} H(z^1, \dots , z^L) = \sum _{\ell =1}^L h(z^\ell , z^{\ell +1}), \end{aligned}$$

(1)

with $z^\ell \in \Omega \subset {\mathbb {R}}^n$ for all $\ell$ and some fixed $n \in {\mathbb {N}}$. We assume periodic boundary conditions and thus set $z^{L+1} = z^1$. A canonical example is a classical particle chain, with $z^\ell = (p_\ell , q_\ell ) \in {\mathbb {R}}^{2d}$ the momentum and position of the $\ell$-th particle ($d \in \{1, 2, 3\}$), and

$$\begin{aligned} H_{\text {pc}}(p_1, \dots , p_L, q_1, \dots , q_L) = \sum _{\ell =1}^L \left( \tfrac{1}{2} p_\ell ^2 + V(q_\ell , q_{\ell +1}) \right) \end{aligned}$$

(2)

consisting of kinetic energy terms $\tfrac{1}{2} p_\ell ^2$ and the potential energy described by V.

The partition function of the physical system is defined as

$$\begin{aligned} Z_L(\beta ) = \int _{\Omega } \cdots \int _{\Omega } {{\,\mathrm{e}\,}}^{-\beta H(z^1, \dots , z^L)} \mathrm {d}z^1 \cdots \mathrm {d}z^L, \end{aligned}$$

(3)

with $\beta = 1/(k_{\text {B}} T) \in {\mathbb {R}}^{+}$ the “inverse temperature”. Our goal is to take the thermodynamic limit $L \rightarrow \infty$ and compute the canonical free energy density

$$\begin{aligned} F(\beta ) = -\beta ^{-1} \lim _{L \rightarrow \infty } \frac{1}{L} \log Z_L(\beta ). \end{aligned}$$

(4)

2 Method

To evaluate the partition function numerically, we combine the well-known transfer-matrix method [3,4,5] with a numerical discretization of integral kernels based on numerical quadrature methods. The latter idea traces back to Nyström [6] and was recently applied by the second author to the computation of Fredholm determinants [7]. Our approach is conceptually similar to [4], but generalizes the method there by considering general quadrature formulas thus making it computationally significantly more efficient.

2.1 Assumptions

To keep technical preliminaries as simple as possible, we make the following assumptions, which cover the applications in the next section: Let $\Omega$ be a (not necessarily finite) measurable subset of ${\mathbb {R}}^n$ and $\nu$ be a finite positive measure on $\Omega$ with a strictly positive, continuous density $\omega$ (weight function). Now we rewrite the partition function in the form

$$\begin{aligned} Z_L(\beta ) = \int _{\Omega } \cdots \int _{\Omega } \prod _{\ell =1}^L k_{\beta }(z^\ell , z^{\ell +1}) \, \mathrm {d}\nu (z^1) \cdots \mathrm {d}\nu (z^L) \end{aligned}$$

(5)

with $z^{L+1} = z^1$ as before. Here the integral kernel

$$\begin{aligned} k_{\beta }: \Omega \times \Omega \rightarrow {\mathbb {R}}, \quad k_{\beta }(z, z') = \frac{{{\,\mathrm{e}\,}}^{-\beta h(z, z')}}{\sqrt{\omega (z) \omega (z')}} \end{aligned}$$

(6)

is assumed to be continuous and bounded; it is strictly positive by construction. We further assume that $h(z, z') = h(z', z)$, i.e., $k_\beta$ is a symmetric kernel. These assumptions imply, in particular,

$$\begin{aligned} \int _{\Omega } \int _{\Omega } |k_{\beta }(z, z') |^2 \, \mathrm {d}\nu (z) \, \mathrm {d}\nu (z') < \infty , \end{aligned}$$

(7)

so that the symmetric kernel $k_\beta$ induces an integral operator [8, §16.1]

$$\begin{aligned} ({\mathcal {T}}_{\beta } u)(z) = \int _{\Omega } k_{\beta }(z, z') u(z') \, \mathrm {d}\nu (z'), \quad u \in L^2(\Omega , \nu ). \end{aligned}$$

(8)

on the Hilbert space $L^2(\Omega ,\nu )$ that is self-adjoint and Hilbert–Schmidt (and hence, compact). Since the product of two Hilbert-Schmidt operators is of trace class and the trace class operators form an ideal within the bounded ones, Eq. (5) for $L \ge 2$ may be represented as¹

$$\begin{aligned} Z_L(\beta ) = {{\,\mathrm{Tr}\,}}\!\left( {\mathcal {T}}_{\beta }^L\right) . \end{aligned}$$

(9)

2.2 Dominant Eigenvalue and Free Energy

Following the mathematical theory of compact operators on Hilbert spaces [8, §21.2], the non-zero elements of the spectrum of ${\mathcal {T}}_{\beta }$ (which is real since the operator is selfadjoint) are at most countably many eigenvalues of finite multiplicity that accumulate only at 0 (which belongs to the spectrum).

Generalizing the Perron–Frobenius theory in matrix analysis [9, §8.2], Jentzsch’s Theorem [10, §8.7, Satz 3] states that the Hilbert-Schmidt operator ${\mathcal {T}}_\beta$ with strictly positive kernel $k_\beta$ has a simple, dominant, strictly positive eigenvalue. That is, all the non-zero eigenvalues can be ordered as (with each eigenvalue listed as often as given by its multiplicity)

$$\begin{aligned} \lambda _1({\mathcal {T}}_\beta )> |{\lambda _2({\mathcal {T}}_\beta )}| \ge |{\lambda _3({\mathcal {T}}_\beta )}| \ge \cdots > 0. \end{aligned}$$

Hence, evaluating Eq. (9) together with Lidskii’s theorem [8, §30.3], that is,

$$\begin{aligned} {{\,\mathrm{Tr}\,}}\!\left( {\mathcal {T}}_{\beta }^L\right) = \sum _{j}\lambda _j({\mathcal {T}}_{\beta })^L \qquad (L\ge 2), \end{aligned}$$

leads to

$$\begin{aligned} \lim _{L \rightarrow \infty } \frac{1}{L} \log Z_L(\beta ) = \lim _{L \rightarrow \infty } \frac{1}{L} \log \sum _j \lambda _j({\mathcal {T}}_{\beta })^L = \log \lambda _1({\mathcal {T}}_{\beta }),\qquad F(\beta ) = -\beta ^{-1} \log \lambda _1({\mathcal {T}}_{\beta }). \end{aligned}$$

(10)

Thus the calculation of the free energy amounts to computing the dominant eigenvalue of ${\mathcal {T}}_{\beta }$.

2.3 Nyström Method for Computing the Dominant Eigenvalue

For simplicity, we restrict the discussion of the numerical method to $n=1$, that is, $\Omega \subset {\mathbb {R}}$. For the multi-dimensional case $n>1$, see the discussion of the example in Sect. 3.3.

We make use of a numerical quadrature rule of the form

$$\begin{aligned} \int _\Omega f(z)\,d\nu (z) \approx \sum _{i=1}^m w_i f(z_i) \end{aligned}$$

(11)

with positive weights $w_i > 0$ and nodes $z_i$. The approximation is of order $p\ge 1$ if the quadrature rule is exact for polynomials of degree at most $p-1$. It is possible to construct such a unique quadrature rule of maximum order $p=2m$, called Gauss quadrature rule; see [11] for the classical construction of Gaussian weights and nodes from the tridiagonal Jacobi matrix of the orthogonal polynomials associated with the measure $\nu$. Error estimates depend on the regularity of the integrand f, e.g., [12, §4.8]: the error for $m\rightarrow \infty$ is of the form $O(m^{-k})$ if $f \in C^{k-1,1}$ and of the form $O(e^{-c m})$ with some $c>0$ if f extends analytically to an ellipse or (semi)-strip in the complex plane containing $\Omega$. The former estimate is called algebraic convergence, the latter exponential convergence.

Inserting the quadrature rule into Eq. (8) results in

$$\begin{aligned} ({\mathcal {T}}_{\beta } u)(z_i) \approx \sum _{j=1}^m k_{\beta }(z_i, z_j) w_j u(z_j) \end{aligned}$$

(12)

for all $i = 1, \dots , m$. Setting $u_i = \sqrt{w_i} u(z_i)$, we have thus discretized the integral operator by the symmetric $m \times m$ matrix (Nyström method)

$$\begin{aligned} T_{\beta } = \big (k_{\beta }(z_i, z_j)\,\sqrt{w_i\,w_j}\big )_{i,j=1}^m. \end{aligned}$$

(13)

By the Perron–Frobenius theory [9, Cor. 8.2.6], this (element-wise) strictly positive matrix has a simple, dominant, strictly positive eigenvalue $\lambda _1(T_\beta )$. The following theorem shows that $\lambda _1(T_{\beta }) \approx \lambda _1({\mathcal {T}}_{\beta })$, and thus for the free energy density in Eq. (4)

$$\begin{aligned} F(\beta ) \approx -\beta ^{-1} \log \lambda _1(T_{\beta }), \end{aligned}$$

(14)

with a speed of convergence, algebraic or exponential, depending on the smoothness of $k_\beta$.

Theorem 1

Under the assumptions in Sect. 2.1 the error $\epsilon _m = \lambda _1(T_{\beta }) - \lambda _1({\mathcal {T}}_{\beta })$ of the Nyström method for computing the dominant eigenvalue behaves as follows as $m\rightarrow \infty$:

• if $k_\beta \in C^{k-1,1}$ then there is algebraic convergence $\epsilon _m=O(m^{-k})$;

• if $k_\beta$ extends analytically to $E\times E$, where $\Omega \subset E\subset {\mathbb {C}}$ is an ellipse or (semi-)strip, there is exponential convergence $\epsilon _m=O(e^{-c m})$ with some constant $c>0$.

The classical, though long and quite involved proof goes by the theory of collectively compact operators and can be found in many books on the numerical treatment of integral equations, e.g., [13, Thm. 4.8.20]. Note, that though only the case of algebraic convergence is covered in this reference, the proof extends literally to the case of exponential convergence by adjusting the consistency assumptions accordingly.

An alternative, conceptually much simpler novel proof based on the theory of the Fredholm determinant can be found in the Appendix of this paper.

Remark

While the Nyström method for the dominant eigenvalue essentially inherits the convergence properties of an underlying cubature formula also in the multidimensional case $n>1$ (see the proof in the Appendix for guidance), we refrained from formulating a general theorem since the characterization of convergence properties of general cubature formulae is much more involved in the first place. If, however, a Gaussian quadrature is applied coordinate-wise, Theorem 1 extends in a straightforward fashion (in the case of exponential convergence $E\subset {\mathbb {C}}^n$ has then to be chosen as an poly-ellipse or poly-(semi-)strip), cf. the example in Sect. 3.3.

3 Applications

We demonstrate the range of applicability of the method via a diverse selection of model systems.

3.1 Classical Particle Chain

Consider a classical particle chain governed by a Hamiltonian of the form (2) above. For evaluating the partition function, the integration over the momentum variables can be performed in closed form, such that

$$\begin{aligned} Z_L(\beta ) = \left( \frac{2\pi }{\beta }\right) ^{L/2} {\tilde{Z}}_L(\beta ) \end{aligned}$$

(15)

with

$$\begin{aligned} {\tilde{Z}}_L(\beta ) = \int _{-\infty }^{\infty } \cdots \int _{-\infty }^{\infty } \prod _{\ell =1}^L {{\,\mathrm{e}\,}}^{-\beta V(q_\ell , q_{\ell +1})} \mathrm {d}q_1 \cdots \mathrm {d}q_L. \end{aligned}$$

(16)

As specific example for the following, we choose

$$\begin{aligned} V(q, q') = V_{\text {loc}}(q) + \tfrac{1}{2} \gamma (q - q')^2 \end{aligned}$$

(17)

with anharmonic on-site potential

$$\begin{aligned} V_{\text {loc}}(q) = \tfrac{1}{2} \eta q^2 + \tfrac{1}{6} \mu q^3 + \tfrac{1}{24} \lambda q^4 \end{aligned}$$

(18)

and coefficients $\lambda , \mu , \eta , \gamma \in {\mathbb {R}}$, $\eta > 0$, $|{\lambda }| \ge |{\mu }|$. The $\sim q^4$ term ensures that V grows asymptotically to infinity as $|{q}| \rightarrow \infty$, and the above $V(q, q')$ is equivalent to its symmetrized version $\frac{1}{2} (V_{\text {loc}}(q) + V_{\text {loc}}(q')) + \frac{1}{2} \gamma (q - q')^2$.

We now assign the term $\frac{1}{2} \eta q^2$ from the on-site potential to the weight function:

$$\begin{aligned} \omega : {\mathbb {R}}\rightarrow {\mathbb {R}}^+, \quad \omega (q) = \frac{{{\,\mathrm{e}\,}}^{-\frac{1}{2} \beta \eta q^2}}{\sqrt{2\pi /(\beta \eta )}}. \end{aligned}$$

(19)

This leads to a (rescaled) Gauss-Hermite quadrature rule [14, §3.5(v)], and we denote the weights by $w_i$ and nodes by $z_i$, $i = 1, \dots , m$, as above.² The particular choice of the weight function is somewhat arbitrary—we could have also included the cubic and quartic terms from the on-site potential into the weight function; at the expense, though, of a less straightforward computation of the weights and nodes of the quadrature rule. The general reasoning is to capture most of the local weight while retaining a well-behaved kernel for the genuine inter-particle potential [see Eq. (21] below).

Combining the weight function (19) with Eq. (16) leads to

$$\begin{aligned} {\tilde{Z}}_L(\beta ) = \left( \frac{2\pi }{\beta \eta }\right) ^{L/2} \int _{-\infty }^{\infty } \cdots \int _{-\infty }^{\infty } \prod _{\ell =1}^L k_{\beta }(q_\ell , q_{\ell +1}) \, \mathrm {d}\nu (q_1) \cdots \mathrm {d}\nu (q_L) \end{aligned}$$

(20)

with the symmetrized kernel

$$\begin{aligned} k_{\beta }(q, q') = {{\,\mathrm{e}\,}}^{-\frac{1}{12} \beta \mu (q^3 + {q'}^3) - \frac{1}{48} \beta \lambda (q^4 + {q'}^4) - \frac{1}{2} \beta \gamma (q - q')^2}. \end{aligned}$$

(21)

We then assemble the symmetric matrix in Eq. (13). Finally, after taking all prefactors in to account, the numerical approximation of the free energy density reads

$$\begin{aligned} - \beta F(\beta ) = \lim _{L \rightarrow \infty } \frac{1}{L} \log Z_L(\beta ) \approx \log \frac{2\pi }{\beta } - \frac{1}{2} \log \eta + \log \lambda _1(T_{\beta }). \end{aligned}$$

(22)

The special case $\gamma = 0$ serves as reference, since the partition function factorizes in this case, i.e., $Z_L(\beta )\vert _{\gamma = 0} = (Z_1(\beta )\vert _{\gamma = 0})^L$ with

$$\begin{aligned} \log Z_1(\beta )\vert _{\gamma = 0} = \frac{1}{2} \log \frac{2\pi }{\beta } + \log \int _{-\infty }^{\infty } {{\,\mathrm{e}\,}}^{-\beta V_{\text {loc}}(q)} \mathrm {d}q. \end{aligned}$$

(23)

For the following numerical examples, we set $\lambda = \mu = \frac{1}{5}$ and $\eta = 1$. Figure 1 visualizes the kernel in Eq. (21). The factor ${{\,\mathrm{e}\,}}^{-\frac{1}{2} \beta \gamma (q - q')^2}$ localizes the kernel around the line $q = q'$, which poses a challenge for accurately “sampling” it using a limited number of points $z_i$ in Eq. (13).

Fig. 1

Integration kernel (21) of the particle chain, for the factorized case $\gamma = 0$ and generic $\gamma = 1$. The remaining parameters are $\lambda = \mu = \frac{1}{5}$ and $\beta = 5$

Fig. 2

a Free energy of the classical particle chain based on a Hamiltonian of the form (2) with potential (17), for parameters $\lambda = \mu = \frac{1}{5}$ and $\eta = 1$. b Corresponding convergence plot of the free energy computation in dependence of the number of quadrature points m, for $\beta = 5$. The reference values for $\gamma = 0$ have been obtained via Eq. (23), and for $\gamma > 0$ by using a large $m = 30$

Figure 2a shows the free energy as function of $\beta$, for several values of $\gamma$, and Fig. 2b the relative error depending on the number of quadrature points m. One observes exponential convergence. The less favorable shape of the kernel with increasing $\gamma$, as mentioned above (see also Fig. 1), translates to a slower convergence rate.

Based on the free energy one can obtain averages and higher-order cumulants following the well-known procedure based on derivatives of F. For example, the average squared particle distance and energy per site are

$$\begin{aligned} \left\langle \tfrac{1}{2} (q_\ell - q_{\ell +1})^2 \right\rangle = \partial _{\gamma } F, \qquad \langle e_\ell \rangle = \partial _{\beta } (\beta F), \end{aligned}$$

(24)

independent of $\ell$ by translation invariance. In practice, a higher-order finite difference scheme on a fine grid is well suited to calculate the derivatives. Figure 3 shows these averages, for the same parameters as before. One notices that the average energy hardly depends on $\gamma$.

Fig. 3

Average quantities in Eq. (24) based on derivatives of the free energy and plotted on a logarithmic scale in dependence of $\beta$, using $m = 30$ quadrature points. The derivatives are numerically approximated via a finite difference scheme of order 6

As a remark, for the case of vanishing on-site potential, $V_{\text {loc}} \equiv 0$, the model conserves momentum, and the statistical mechanics description changes accordingly [16]. Numerically computing the free energy is less challenging in this case since the partition function factorizes after introducing the “stretch” $r_{\ell } = q_{\ell +1} - q_\ell$.

3.2 Discrete Nonlinear Schrödinger Equation

The method described in Sect. 2 has been employed in the work [17] on the discrete nonlinear Schrödinger equation. Here we present and elaborate on the numerical aspects in more detail, and compare thermodynamic expectation values with results of molecular dynamics (MD) simulations.

To be self-contained, we first restate the physical setup: the central object is a complex-valued wave field $\psi _\ell$ ($\ell = 1, \dots , L$) governed by the semiclassical Hamiltonian

$$\begin{aligned} H(\psi _1, \dots , \psi _L) = \sum _{\ell =1}^L \big ( \tfrac{1}{2} |{\psi _{\ell +1} - \psi _\ell }|^2 + \tfrac{1}{2}\,g\,|{\psi _\ell }|^4\big ), \end{aligned}$$

(25)

with parameter $g > 0$ (so-called defocusing case). The corresponding partition function reads

$$\begin{aligned} Z_L(\mu , \beta ) = \int {{\,\mathrm{e}\,}}^{-\beta (H - \mu N)} \, \mathrm {d}\psi _1 \cdots \mathrm {d}\psi _L, \end{aligned}$$

(26)

where we have introduced the chemical potential $\mu$ as additional parameter, which is dual to the total particle number $N = \sum _{\ell =1}^L |{\psi _\ell }|^2$.

A symplectic change of variables to polar coordinates leads to the representation

$$\begin{aligned} \psi _\ell = \sqrt{\rho _\ell }\, {{\,\mathrm{e}\,}}^{\mathrm {i} \varphi _\ell } \end{aligned}$$

(27)

with $\rho _\ell \in \Omega = [0, \infty )$, and $\varphi _\ell \in S^1$ (unit circle). The Hamiltonian in these variables reads

$$\begin{aligned} H = \sum _{\ell =1}^L \big ( {-} \sqrt{\rho _{\ell +1}\,\rho _\ell }\, \cos (\varphi _{\ell +1} - \varphi _\ell ) + \rho _\ell + \tfrac{1}{2}\,g\,\rho _\ell ^2 \big ). \end{aligned}$$

(28)

It depends only on phase differences, which implies the invariance under the global shift $\varphi _\ell \mapsto \varphi _\ell + \phi$. For evaluating the partition function, the $\varphi _\ell$ integrals can be calculated in closed form [18]. This leads to

$$\begin{aligned} Z_L(\mu , \beta ) = {{\,\mathrm{e}\,}}^{\beta \frac{1}{2} \mu ^2 L/g} \int _{\Omega } \cdots \int _{\Omega } \prod _{\ell =1}^L k_{\beta }(\rho _\ell , \rho _{\ell +1}) \, {{\,\mathrm{e}\,}}^{-\beta \frac{1}{2} g \left( \rho _\ell - \frac{\mu }{g}\right) ^2} \, \mathrm {d}\rho _1 \cdots \mathrm {d}\rho _L \end{aligned}$$

(29)

with the kernel

$$\begin{aligned} k_{\beta }(\rho , \rho ') = 2\pi I_0\big (\beta \sqrt{\rho \,\rho '}\big )\, {{\,\mathrm{e}\,}}^{-\beta \frac{1}{2} (\rho + \rho ')}. \end{aligned}$$

(30)

$I_0$ is the modified Bessel function of the first kind. For this example, we use the second factor of the integrand in Eq. (29) as weight function:

$$\begin{aligned} \omega : \Omega \rightarrow {\mathbb {R}}^+, \quad \omega (z) = c {{\,\mathrm{e}\,}}^{-a (z - b)^2/2} \end{aligned}$$

(31)

with $a = \beta g$, $b = \mu /g$ and the normalization constant

$$\begin{aligned} c = \frac{2 \sqrt{\frac{a}{2 \pi }}}{1 + {{\,\mathrm{erf}\,}}(b \sqrt{\frac{a}{2}})}. \end{aligned}$$

(32)

After constructing a Gauss quadrature rule on $\Omega = [0, \infty )$ as described in Sect. 2.3, we form the symmetric matrix in Eq. (13), here denoted $T_{\mu , \beta }$ since it implicitly also depends on $\mu$ via the quadrature points and weights. Then, taking the prefactor in Eq. (29) and the normalization constant (32) into account, one arrives at

$$\begin{aligned} - \beta F(\mu , \beta ) = \lim _{L \rightarrow \infty } \frac{1}{L} \log Z_L(\mu , \beta ) \approx \beta \,\tfrac{1}{2} \tfrac{\mu ^2}{g} + \log \lambda _1(T_{\mu , \beta }) - \log c(\mu , \beta ). \end{aligned}$$

(33)

Fig. 4 shows the free energy as function of $\beta$, for various values of $\mu$.

Fig. 4

Free energy of the discrete nonlinear Schrödinger equation governed by the Hamiltonian (25) on an infinite chain, for parameter $g = 1$

Numerically, we again observe exponential convergence with respect to the number of quadrature points, see Fig. 5. At $\beta = 15$ and $\mu = 1$ for example, $m = 16$ points suffice for double precision accuracy. The reference data stems from a calculation with $m = 20$.

Fig. 5

Relative error of the free energy calculation for the discrete nonlinear Schrödinger equation depending on the number of quadrature points m, exemplified for $\beta = 1$ and $\beta = 15$

Fig. 6

Average density and energy for the discrete nonlinear Schrödinger equation. Solid curves are computed via Eq. (34), and the dots show the results of MD simulations of a system with $L = 4096$ sites, using overdamped Langevin dynamics for equilibration. The error bars (of the size of line thickness) visualize the standard deviation based on 128 independent realizations

One can obtain thermodynamic averages based on derivatives of $F(\mu , \beta )$. For example, the average density and energy per lattice site are

$$\begin{aligned} \langle \rho _\ell \rangle = - \partial _{\mu } F(\mu , \beta ), \qquad \langle e_\ell \rangle = \partial _{\beta } (\beta \,F(\mu , \beta )) + \mu \,\langle \rho _\ell \rangle . \end{aligned}$$

(34)

Fig. 6 shows these quantities as function of $\beta$, and compares them with molecular dynamics simulations of the microscopic model (system size $L = 4096$). As described in [17], we equilibrate the system based on overdamped Langevin dynamics [19], using 1024 Langevin steps. One observes very good agreement, as expected.

3.3 Classical Oscillators on a Cylindrical Lattice

The numerical method is in principle also applicable to two-dimensional lattices, by using periodic boundary conditions in one direction and reducing the setting to a quasi-one dimensional problem. Specifically, we consider the lattice $\Gamma = {\mathbb {Z}}/(L_x) \otimes {\mathbb {Z}}/(L_y)$ for $L_x, L_y \in {\mathbb {N}}$, i.e., starting with periodic boundary conditions both in x- and y-direction, but eventually sending $L_x \rightarrow \infty$ while keeping $L_y$ finite. Thus we arrive at a cylindrical lattice, as visualized in Fig. 7.

Fig. 7

Cylindrical lattice topology, using periodic boundary conditions in y-direction. Each small dot is a lattice site

We identify a lattice site by the index $\ell = (\ell _x, \ell _y) \in \Gamma$, and consider for simplicity scalar spatial variables $q_\ell \in {\mathbb {R}}$; these could be displacements from the reference positions in one fixed direction, for example. $p_\ell$ denotes the momentum of the $\ell$-th particle.

As demonstration, let the system be governed by the Hamiltonian

$$\begin{aligned} H = \sum _{\ell \in \Gamma } \left( \tfrac{1}{2} p_\ell ^2 + V_{\text {loc}}(q_\ell )\right) + \sum _{\langle \ell , \ell ' \rangle } V_{\ell ' - \ell }(q_\ell , q_{\ell '}), \end{aligned}$$

(35)

consisting of site-local kinetic and potential energy terms (first sum) as well as nearest neighbor interactions (second sum). Specifically, we consider a local quadratic potential $V_{\text {loc}}(q) = \frac{1}{2} \eta q^2$, $\eta > 0$, and an interaction potential $V_{\Delta \ell }(q, q') = \frac{1}{2} a_{\Delta \ell } (q - q')^2$ with two coefficients $a_{(\pm 1, 0)} = a_x$ and $a_{(0, \pm 1)} = a_y$.

To cast the Hamiltonian (35) into the form of Eq. (1), we subsume the particles contained in one lattice “ring” into the vectors

$$\begin{aligned} \mathbf {p}_{\ell _x}&= \left( p_{(\ell _x, 1)}, \dots , p_{(\ell _x, L_y)}\right) \in {\mathbb {R}}^{L_y}, \end{aligned}$$

(36)

$$\begin{aligned} \mathbf {q}_{\ell _x}&= \left( q_{(\ell _x, 1)}, \dots , q_{(\ell _x, L_y)}\right) \in {\mathbb {R}}^{L_y} \end{aligned}$$

(37)

for $\ell _x = 1, \dots , L_x$.

Similar to the particle chain in Sect. 3.1, the momentum integration for evaluating the partition function can be performed explicitly. This results in

$$\begin{aligned} Z_{(L_x,L_y)}(\beta ) = \left( \frac{2\pi }{\beta }\right) ^{L_x L_y /2} {\tilde{Z}}_{(L_x,L_y)}(\beta ) \end{aligned}$$

(38)

with

$$\begin{aligned} {\tilde{Z}}_{(L_x,L_y)}(\beta ) = \int _{{\mathbb {R}}^{L_y}} \cdots \int _{{\mathbb {R}}^{L_y}} \prod _{\ell _x=1}^{L_x} {{\,\mathrm{e}\,}}^{-\beta \left( {\tilde{V}}_{\text {loc}}(\mathbf {q}_{\ell _x}) + {\tilde{V}}_{\text {int}}(\mathbf {q}_{\ell _x}, \mathbf {q}_{\ell _x+1})\right) } \mathrm {d}^{L_y}q_1 \cdots \mathrm {d}^{L_y}q_{L_x}, \end{aligned}$$

(39)

where

$$\begin{aligned} {\tilde{V}}_{\text {loc}}(\mathbf {q}) = \sum _{\ell _y=1}^{L_y} V_{\text {loc}}(q_{\ell _y}) \end{aligned}$$

(40)

and

$$\begin{aligned} {\tilde{V}}_{\text {int}}(\mathbf {q}, \mathbf {q}') = \sum _{\ell _y=1}^{L_y} \left( \tfrac{1}{2} a_x (q_{\ell _y} - q_{\ell _y}')^2 + \tfrac{1}{4} a_y (q_{\ell _y} - q_{\ell _y+1})^2 + \tfrac{1}{4} a_y (q_{\ell _y}' - q_{\ell _y+1}')^2\right) . \end{aligned}$$

(41)

Note that ${\tilde{V}}_{\text {int}}$ is symmetric in its arguments, i.e., ${\tilde{V}}_{\text {int}}(\mathbf {q}, \mathbf {q}') = {\tilde{V}}_{\text {int}}(\mathbf {q}', \mathbf {q})$ for any $\mathbf {q}, \mathbf {q}' \in {\mathbb {R}}^{L_y}$.

It suggests itself to use the factor ${{\,\mathrm{e}\,}}^{-\beta {\tilde{V}}_{\text {loc}}(\mathbf {q})}$ in Eq. (39) as integration measure, resulting in a tensor product of normal distributions:

$$\begin{aligned} \omega : {\mathbb {R}}^{L_y} \rightarrow {\mathbb {R}}^+, \quad \omega (\mathbf {q}) = \frac{{{\,\mathrm{e}\,}}^{-\frac{1}{2} \beta \eta \Vert {\mathbf {q}}\Vert ^2}}{\big (\frac{2\pi }{\beta \eta }\big )^{L_y/2}} = \prod _{\ell _y=1}^{L_y} \frac{{{\,\mathrm{e}\,}}^{-\frac{1}{2} \beta \eta q_{\ell _y}^2}}{\sqrt{2\pi /(\beta \eta )}}. \end{aligned}$$

(42)

We use a rescaled Gauss-Hermite quadrature rule along each coordinate direction, as in Sect. 3.1. Theorem 1 extends straightforwardly to this choice. An alternative, which is less affected by the inherent curse of dimensionality, is a cubature rule dedicated to multidimensional integration [20,21,22,23], or sparse grid methods. The convergence properties of such cubature rules are more involved, but they would essentially be inherited by the Nyström method for the dominant eigenvalue. We leave an exploration of these ideas for future work.

Using $\mathrm {d}\nu (q) = \omega (q)\mathrm {d}q$ as measure, Eq. (39) becomes

$$\begin{aligned} {\tilde{Z}}_{(L_x,L_y)}(\beta ) = \left( \frac{2\pi }{\beta \eta }\right) ^{L_x L_y/2} \int _{{\mathbb {R}}^{L_y}} \cdots \int _{{\mathbb {R}}^{L_y}} \prod _{\ell _x=1}^{L_x} k_{\beta }(\mathbf {q}_{\ell _x}, \mathbf {q}_{\ell _x+1}) \, \mathrm {d}\nu (q_1) \cdots \mathrm {d}\nu (q_{L_x}) \end{aligned}$$

(43)

with the kernel

$$\begin{aligned} k_{\beta }(\mathbf {q}, \mathbf {q}') = {{\,\mathrm{e}\,}}^{-\beta {\tilde{V}}_{\text {int}}(\mathbf {q}, \mathbf {q}')}. \end{aligned}$$

(44)

Following the factorized quadrature rule, the symmetric matrix in Eq. (13) takes the form

$$\begin{aligned} T_{\beta } = \big (k_{\beta }(\mathbf {q}_{\mathbf {i}}, \mathbf {q}_{\mathbf {j}})\,\sqrt{w_{\mathbf {i}}\,w_{\mathbf {j}}}\big )_{\mathbf {i},\mathbf {j}} \end{aligned}$$

(45)

with multi-indices $\mathbf {i}, \mathbf {j} \in \{1, \dots , m_0\}^{L_y}$ and the definitions $\mathbf {q}_{\mathbf {i}} = (q_{i_1}, \dots , q_{i_{L_y}})$, $\omega _{\mathbf {i}} = \omega _{i_1} \cdots \omega _{i_{L_y}}$ and $w_i$, $q_i$, $i = 1, \dots , m_0$ the weights and points of the one-dimensional rescaled Gauss-Hermite quadrature rule. Thus the overall number of weights and points is $m = m_0^{L_y}$.

The numerical approximation of the free energy per lattice site is then

$$\begin{aligned} - \beta F(\beta ) = \lim _{L_x \rightarrow \infty } \frac{1}{L_x L_y} \log Z_{(L_x,L_y)}(\beta ) \approx \log \frac{2\pi }{\beta } - \frac{1}{2} \log \eta + \frac{1}{L_y} \log \lambda _1(T_{\beta }), \end{aligned}$$

(46)

with $L_y$ kept fixed. For the following examples, we set $L_y = 3$, such that the overall number of quadrature points $m = m_0^{L_y}$ remains manageable up to $m_0 = 8$. Figure 8a visualizes the free energy as function of $\beta$, for several combinations of $a = (a_x, a_y)$. One notices that the curve for $a = (\frac{1}{2}, \frac{1}{5})$ is visually indistinguishable from the case with interchanged parameters $a_x \leftrightarrow a_y$, pointing to the conclusion that the influence of a finite $L_y$ compared to the “infinite” $L_x$ on the free energy is quite small.

Fig. 8

a Free energy of the cylindrical lattice model governed by the Hamiltonian (35), for $L_y = 3$ and $\eta = 1$. b Corresponding convergence plot of the free energy computation in dependence of the number of quadrature points $m_0$ along one dimension, for $\beta = 5$. The reference values for $a_x = 0$ have been obtained via Eq. (47), and for $a_x \ne 0$ by using $m_0 = 8$

Figure 8b shows the corresponding relative error depending on the number of quadrature points $m_0$ along one dimension. In the special case $a_x = 0$ (without coupling in x-direction), the partition function factorizes, such that, analogous to Sect. 3.1, ${\tilde{Z}}_{(L_x,L_y)}(\beta )\vert _{a_x = 0} = ({\tilde{Z}}_{(1,L_y)}(\beta )\vert _{a_x = 0})^{L_x}$ with

$$\begin{aligned} {\tilde{Z}}_{(1,L_y)}(\beta )\vert _{a_x = 0} = \int _{{\mathbb {R}}^{L_y}} {{\,\mathrm{e}\,}}^{-\beta \left( \frac{1}{2} \eta \Vert {\mathbf {q}}\Vert ^2 + \sum _{\ell _y=1}^{L_y} \frac{1}{2} a_y (q_{\ell _y} - q_{\ell _y+1})^2 \right) } \mathrm {d}^{L_y} q. \end{aligned}$$

(47)

We evaluate this integral numerically and use it as reference for computing the relative error in Fig. 8b for $a_x = 0$. The relative error is still rather large for $a = (\frac{1}{2}, \frac{1}{5})$ and $a = (\frac{1}{5}, \frac{1}{2})$; as before, this observation can be explained by the difficulty of accurately sampling the kernel (44) via (45) using a small number of quadrature points along each coordinate. To mitigate this issue for the present example, one could associate the $a_y$ terms in the Hamiltonian to the integration measure $\omega$ instead of the kernel, at the expense of a more complicated quadrature rule.

In summary, this application example demonstrates that our method can in principle handle two-dimensional lattice topologies as well, although the large number of required quadrature points (when interpreting the problem as quasi one-dimensional) limits the size of the periodic dimension $L_y$ in practice.

4 Conclusion and outlook

The convergence plots of the example applications illustrate the validity of Theorem 1, which states that approximating the dominant eigenvalue of the discretized kernel inherits the favorable exponential convergence properties of the underlying quadrature rule in the case of kernels extending analytically into the complex domain.

To optimize the numerical performance of the method for the cylindrical topologies, one could further exploit the factorized structure of $k_{\beta }$ in (44) along coordinate directions, or use sparse grid methods for the quadrature as mentioned above.

Concerning transfer operator techniques in general, it could be fruitful to adopt ideas from quantum mechanics (see e.g. [24]), or set oriented numerical methods [25].

5 Appendix

We give here an alternative, conceptually much simpler proof³ of Theorem 1 based on the theory of the Fredholm determinant of the kernel $k_\beta$, namely

$$\begin{aligned} D(\mu ) = 1 + \sum _{n=1}^\infty \frac{(-\mu )^n}{n!} \int _{\Omega }\cdots \int _{\Omega }\det _{i,j=1}^n k_\beta (z_i,z_j) \,d\nu (z_1)\cdots d\nu (z_n). \end{aligned}$$

Given the assumptions in Sect. 2.1, D is an entire function whose roots are exactly the reciprocal non-zero eigenvalues of ${\mathcal {T}}_\beta$, see [10, §6.2, Satz 3]. The Weierstrass product [10, §6.4, Satz 1]

$$\begin{aligned} D(\mu ) = e^{\alpha \mu + \beta \mu ^2} \prod _n \left( 1-\mu \lambda _n({\mathcal {T}}_\beta )\right) e^{\mu \lambda _n({\mathcal {T}}_\beta )} \end{aligned}$$

shows that the multiplicities of the roots of D and the multiplicities of the non-zero eigenvalues of ${\mathcal {T}}_\beta$ agree. In particular, $1/\lambda _1({\mathcal {T}}_\beta )$ is a simple root and therefore $D'(\lambda _1({\mathcal {T}}_\beta )^{-1})\ne 0$. Also, [7, Thm. 6.2] (whose proof can literally be extended to the current assumptions) gives

where, uniformly for bounded $\mu$, the error is given by $\epsilon _m(\mu ) = O(m^{-k})$ or $\epsilon _m(\mu ) = O(e^{-cm})$ according to whether $k_\beta \in C^{k-1,1}$ or $k_\beta$ extends analytically into the complex plane. By Perron–Frobenius $\lambda _1(T_\beta )$ is the simple, dominant, strictly positive eigenvalue of the entry-wise positive matrix $T_\beta$, which by the argument principle of complex analysis must satisfy

$$\begin{aligned} 1/\lambda _1(T_\beta ) \rightarrow 1/\lambda _1({\mathcal {T}}_\beta ) \qquad (m\rightarrow \infty ). \end{aligned}$$

Hence, by inserting $\mu =1/\lambda _1(T_\beta )$ into (#), followed by a Taylor expansion, we get

$$\begin{aligned} \lambda _1(T_\beta ) = \lambda _1({\mathcal {T}}_\beta ) + \frac{\lambda _1({\mathcal {T}}_\beta )^2}{D'(\lambda _1({\mathcal {T}}_\beta )^{-1})} \epsilon _m(\lambda _1({\mathcal {T}}_\beta )^{-1})) + O(\epsilon _m(\lambda _1({\mathcal {T}}_\beta )^{-1}))^2), \end{aligned}$$

which completes the proof.