The monomer-dimer problem and moment Lyapunov exponents of homogeneous Gaussian random fields

We consider an"elastic"version of the statistical mechanical monomer-dimer problem on the n-dimensional integer lattice. Our setting includes the classical"rigid"formulation as a special case and extends it by allowing each dimer to consist of particles at arbitrarily distant sites of the lattice, with the energy of interaction between the particles in a dimer depending on their relative position. We reduce the free energy of the elastic dimer-monomer (EDM) system per lattice site in the thermodynamic limit to the moment Lyapunov exponent (MLE) of a homogeneous Gaussian random field (GRF) whose mean value and covariance function are the Boltzmann factors associated with the monomer energy and dimer potential. In particular, the classical monomer-dimer problem becomes related to the MLE of a moving average GRF. We outline an approach to recursive computation of the partition function for"Manhattan"EDM systems where the dimer potential is a weighted l1-distance and the auxiliary GRF is a Markov random field of Pickard type which behaves in space like autoregressive processes do in time. For one-dimensional Manhattan EDM systems, we compute the MLE of the resulting Gaussian Markov chain as the largest eigenvalue of a compact transfer operator on a Hilbert space which is related to the annihilation and creation operators of the quantum harmonic oscillator and also recast it as the eigenvalue problem for a pantograph functional-differential equation.


1.
Introduction. It is not uncommon that deterministic problems may reveal connections with stochastic counterparts, and it is particularly so if the original formulation inherently involves a pseudo-randomness production mechanism or rich combinatorics. When this happens, probability theoretic methods can provide an additional insight into the deterministic setting. Such are, for example, the problems concerned with chaotic dynamical systems under increasingly refined spatial discretization which models the finite computer arithmetic. Here, the effect of the round-off errors, magnified in the long run by positive Lyapunov exponents, leads to a subtle interplay between the ergodic behaviour of the underlying continuous system and the computer discretization. This can distort the invariant measure to the extent that, in place of its original absolute continuity, a persistent attracting centre emerges, thus destroying resemblance between the long-term behaviour of the discretised system and its continuous predecessor. As reported in a series of papers by A.V.Pokrovskii, V.S.Kozyakin and their collaborators (see, for example, [39] and references therein), the attracting centre can be modelled through a finite Markov chain with an absorbing state. Moreover, a strikingly accurate numerical reproduction of statistics of discretised systems has been achieved with the subsidiary Markov chains obtained by composing random maps with a common fixed point on the finite state space, sampled independently according to a specially designed probability measure. 1 The idea of random composite maps on a countable set can be extended to the continuum, for example, by considering the compositions ϕ 1 • . . . • ϕ k of jointly Gaussian real-valued random processes ϕ k (x) of a real parameter x. If these processes have differentiable sample paths and share a nonrandom fixed point x * , then its stability under the maps ϕ k , as a random dynamical system, is determined (in the linear approximation) by the asymptotic behaviour of the products of jointly Gaussian random variables ξ k := ϕ k (x * ) as N → +∞. Assuming that the processes ϕ k are all identically distributed, the standard ergodicity argument under a suitable mixing condition (on the decay of correlations between the processes ϕ k with distant values of k) leads to the Lyapunov exponent Here, the expectation E(·) involves only the common probability law of the processes ϕ k as if the latter were independent copies of a "mother" process, in which case (2) would follow immediately from the strong law of large numbers [43]. Similar products of appropriately rescaled random processes have attracted an active research interest due to their applicability to multifractal modelling of turbulence cascades and other complex phenomena; see, for example, [1]. The problem of investigating the asymptotic behaviour of the products ψ N from (1) becomes considerably harder if it is concerned with the moment Lyapunov exponent (MLE) which we define here as provided this limit exists (otherwise, the upper limit or other limit points can be considered instead). In contrast to the conventional Lyapunov exponent (2), the MLE (3) is easy to compute only in the case of independent random variables ξ k in (1). If they are dependent, their joint probability law nontrivially enters the MLE even if correlations in the Gaussian random sequence ξ 1 , ξ 2 , ξ 3 , . . . decay fast enough to ensure mixing. Although the problem of computing the MLE as a modification to the usual Lyapunov exponent (2) may seem artificial in the context of random dynamics stability mentioned above, it has a remarkable connection with a statistical mechanical problem of computing the equilibrium free energy of a system of monomers and dimers on a lattice, and this connection is the main theme of the present paper.
In the classical monomer-dimer problem [3,11,12,13,20,25,29], which is essentially solved in the planar case and remains unsolved in higher dimensions, the monomers are represented by singletons, and the dimers are two-element subsets of the lattice formed by nearest neighbours, with the energy of a dimer depending on its orientation. In particular, in a monomer-free "isotropic" setting on the plane, the partition function counts the number of "domino" tilings of a given finite subset of the planar lattice [3, p. 125]. It turns out that, for any dimension n, the number of such partitions of a finite subset Λ ⊂ Z n of the n-dimensional integer lattice coincides with the product moment of a specially constructed homogeneous Gaussian random field (GRF) ξ := (ξ x ) x∈Z n . Moreover, this GRF is of moving average type which is easy to simulate. This, in principle, provides an alternative probabilistic way (different from the probabilistic algorithm proposed in [27]) for approximate computation of the partition function in the monomerdimer problem by modelling the auxiliary GRF and statistically estimating its product moments instead of the direct Monte-Carlo simulation of the underlying particle system which is complicated by the tilability issue. Due to this connection, the free energy of the monomer-dimer system per lattice site in the thermodynamic limit is proportional to the MLE of the auxiliary GRF, defined, similarly to (3), by where #(·) is the number of elements in a finite set, and Λ → ∞ is understood in the sense of van Hove [41]. This connection between the MLEs of homogeneous GRFs and the statistical mechanical partition functions remains valid for a wider class of monomerdimer systems.
In the present paper, we consider an "elastic" version of the monomer-dimer problem, which extends the original "rigid" formulation by allowing each dimer to consist of two particles at arbitrarily distant sites of the lattice, with the energy of interaction between the particles depending on their relative position. Thus, the geometric constraint that the dimers are formed only by nearest neighbours is removed, thereby introducing the longrange interactions; see Fig. 1. In this more general setting, the partition function of the FIGURE 1. Examples of configurations of the elastic dimer-monomer system in dimensions one, two and three (left to right). Monomers are depicted as "•"s. Dimer particles are represented by "•"s, whilst the bonds between them are shown as arcs in the one-dimensional case and as line segments in the higher dimensions.
elastic dimer-monomer (EDM) system and its thermodynamic limit can be encoded in the product moments and the MLE of an auxiliary homogeneous GRF whose mean value and the covariance function are the Boltzmann factors associated with the monomer energy and the dimer potential, respectively. The structure of the covariance function singles out a class of EDM systems where the dimer potential is a weighted 1 -norm (also known as "Manhattan norm") of the dimer vector. A remarkable feature of the Manhattan dimer potential is that the corresponding auxiliary GRF is Markov and, moreover, belongs to the class of Pickard random fields (PRFs) [37,38,45]; see also, [6,7]. PRFs were developed originally for the planar case and subsequently extended to three dimensions [22]. Unlike general Markov random fields, they are amenable to a unilateral single-pass simulation, for example, using the level-sets algorithm of [8]. The simplicity of the autoregressive-like modelling allows the Gaussian PRF (GPRF) to be utilized for the statistical estimation of the partition function of the Manhattan EDM system, similar to the role of moving average GRFs for the classical monomer-dimer problem. However, in addition to this utility for the alternative Monte-Carlo simulation, we take advantage of the specific Markov structure of the auxiliary GPRF to establish a recurrence relation for its conditional product moments which involves local transitions from one level set to another. Although the exposition is focused mainly on the two-dimensional case for better visualizability, a recurrence equation for the conditional product moments of the GPRF can also be established in three dimensions. Note that similar recurrence relations are present in the statistical mechanical transfer matrix method for systems with short-range interactions [2,3]. Since the interactions in the Manhattan EDM system are of long-range nature, the possibility of recursive representation of its partition function seems striking.
For the one-dimensional Manhattan EDM system, where the dimer potential is proportional to the standard Euclidean distance, the auxiliary GRF becomes a homogeneous Gaussian Markov sequence satisfying an autoregressive equation. Due to the Markov property, the MLE is reduced to the largest eigenvalue of a compact transfer operator on an invariant cone in a Hilbert space which governs a recurrence equation for the conditional product moments. We obtain a matrix representation of the transfer operator in the basis of Hermite polynomials and provide numerical results on the dependence of the MLE on two parameters of the univariate Manhattan EDM system. Although this eigenanalysis lies within the general framework of the transfer matrix method, we employ probability theoretic techniques and constructs such as measure change and conditional Gaussian distributions. This allows the transfer operator to be linked with the annihilation and creation operators of the quantum harmonic oscillator [42, p. 91], and the eigenvalue problem to be recast into a pantograph functional-differential equation [5,10,26,31,44] which involves scaling of the independent variable.
Note that the approach to the monomer-dimer problem through the product moments and MLEs of auxiliary GRFs, proposed in the present paper, is different from, and in fact, more general than, the classical methods based on Pfaffians of matrices (see [3] and references therein) and the links of the problem with the determinants of random matrices [30]. Our machinery differs also from the Gaussian integrals in Grassmann variables and their connections with the Pfaffians [18].
The paper is organised as follows. Section 2 describes the class of EDM systems being considered and formulates the problem of computing the equilibrium free energy in the thermodynamic limit. In Section 3, this problem is related to the product moments and the MLE of a homogeneous GRF. Section 4 shows that the auxiliary GRF for the classical rigid dimer problem is of moving average type. Section 5 relates the family of Manhattan EDM systems with GPRFs and focuses on the two-dimensional case. Section 7 obtains a recurrence relation for their conditional product moments associated with aggregated level sets of Section 6. Section 8 proceeds to one-dimensional Manhattan EDM systems. Section 9 derives the transfer operator which governs the recurrence equation for the conditional product moments of the auxiliary Gaussian Markov chain. Section 10 establishes an inequality for the norm of this operator and a related upper bound for the MLE. Section 11 considers the matrix representation of the transfer operator in the basis of Hermite polynomials and discusses its links with the ladder operators of the quantum harmonic oscillator. Section 12 relates the MLE for the univariate Manhattan EDM system with the eigenvalue problem for a pantograph functional-differential equation. Section 13 provides concluding remarks.
2. Elastic dimer-monomer problem. Consider a system of identical particles residing at sites of a nonempty finite subset Λ of the n-dimensional integer lattice Z n , with each site of Λ being occupied by exactly one particle. For some disjoint two-element subsets {x, y} ⊂ Λ, the particles at the lattice sites x and y are (chemically) bonded and form a dimer. The particles which do not participate in dimers are interpreted as monomers.
The monomers and dimers are endowed with particular values of energy. The individual energies of monomers are all equal and their common value is denoted by V . The energy W (z) associated with a dimer with endpoints x = y is that of the interaction between the particles in it and depends on their relative position z := x − y. Here, W : Z n \ {0} → R {+∞} is a symmetric function which we will refer to as the dimer potential.
The classical monomer-dimer formulation [12,13,25] allows dimers to consist of nearest neighbours only, so that all dimers have a fixed length and the energy of a dimer can only depend on its orientation. The "elastic" dimer-monomer (EDM) setting, outlined above, is more general. Indeed, it still leaves room for geometric constraints which can be imposed by putting W (z) := +∞ for prohibited values of the dimer vector z. On the other hand, for admissible z, where W (z) is finite, the dimer potential W plays the role of a scoring function which favors dimers with smaller energy, though, in general, it does not forbid arbitrarily distant particles to form a bond. As in the classical rigid setting, the sites occupied by monomers can also be interpreted as vacancies in which case the EDM system provides a lattice model for a gas of dimers. We will not, however, employ this alternative interpretation in what follows.
A configuration ω of the EDM system is completely specified by a subset D of evenly many #D = 2r sites of Λ, and a partition {{x 1 , y 1 }, . . . , {x r , y r }} of D into two-element subsets which represent dimers. The class of pair partitions of D is denoted by P D , so that #P D = (2r − 1)!!. With the complementary set Λ \ D being occupied by monomers, the total energy of the particle system in the configuration ω is where N := #Λ. The class Ω Λ := {ω} of all possible configurations of the EDM system associated with the set Λ consists of partitions of Λ into two-element subsets and singletons. Here, · is the floor function, ζ is a standard normal (that is, zero mean unit variance Gaussian) random variable with the probability density function (PDF) and is the kth modified Hermite polynomial. In (7), use is made of the even order moments E(ζ 2r ) = (2r − 1)!! and the moment-generating function Ee λζ = e λ 2 /2 for the standard normal distribution in combination with the generating function of the Hermite polynomials [33, Theorem 2.5.6 on p. 7]: By the variational postulate of Statistical Mechanics describing the canonical ensemble [21], the particle system, in contact with a heat bath at absolute temperature T > 0, acquires an equilibrium probability distribution over the configurational space Ω Λ that delivers the minimum value to the free energy functional over all probability mass functions (PMFs) P : Ω Λ → [0, 1] satisfying ω∈ΩΛ P (ω) = 1. By the finiteness of the set Ω Λ and the strict convexity of the functional being minimized in (11) over the simplex of PMFs, the equilibrium value F Λ of the free energy is achieved at the unique Gibbs- where the units of the temperature T are chosen so as to "absorb" the Boltzmann constant, and is the canonical partition function [21,34] which encodes the equilibrium properties of the system. For example, the average total energy of the particle system is recovered from (13) as ω∈ΩΛ U (ω)P Λ (ω) = −∂ β ln Z Λ . We will be concerned with the statistical mechanical problem of computing the partition function Z Λ and the thermodynamic limit of the equilibrium free energy (11) per site of an unboundedly increasing fragment Λ of the lattice Z n : lim Although the convergence Λ → ∞ is usually understood in the sense of van Hove [41], we will often assume, for simplicity, that Λ is a discrete hypercube {0, . . . , N − 1} n (or another simply shaped subset of the lattice), where N is an integer parameter which tends to +∞.
3. Auxiliary Gaussian random field. Suppose the temperature T is fixed and the dimer potential W is extended to the origin so as to satisfy the condition W (0) = ∞ and the inequality e −βW (0) To make this possible, we assume that the sum on the right-hand side of (15) is finite. Then the function C : Z n → R + , defined as the Boltzmann factor associated with the dimer potential W by C(z) := e −βW (z) , (16) with the shorthand notation (12), is positive definite [15]. That is, for any positive integer m and any points z 1 , . . . , z m in Z n , the matrix (C(z j − z k )) 1 j,k m is positive semi-definite. The fact that (15) implies the positive definiteness of C in (16) is established by a diagonal dominance argument [19, p. 349]. Hence, by Bochner's theorem, C is the covariance function of a real-valued homogeneous random field on Z n which can be chosen to be Gaussian (without additional assumptions on C). Therefore, the fulfillment of (15) ensures the existence of a homogeneous Gaussian random field (GRF) ξ := (ξ x ) x∈Z n satisfying for all x, y ∈ Z n . Here, cov(·, ·) is the covariance of random variables, V is the monomer energy, and C is given by (16). In what follows, such a random field ξ will be referred to as the auxiliary GRF. Its significance for the EDM system, described in Section 2, is clarified by the theorem below.
Theorem 3.1. Suppose the dimer potential W satisfies (15). Then for any finite set Λ ⊂ Z n , the partition function (13) coincides with the product moment (4) of the auxiliary GRF ξ over Λ: Proof. The product moment of the auxiliary GRF ξ over the set Λ can be computed as where N := #Λ, and η := (η x ) x∈Z n is a GRF obtained by centering ξ as η in conformance with (17). By the Isserlis-Wick theorem [23], which relates mixed moments of evenly many zero mean Gaussian random variables with their cross-covariances, 2 it follows that Here, the sum extends over the class P D of partitions {{x 1 , y 1 }, . . . , {x r , y r }} of the set D into two-element subsets, with r := #D/2. Now, upon substitution of (16) into (21), (19) takes the form In view of (6) and the definition of the configurational space Ω Λ of the EDM system from Section 2, the right-hand side of (22) coincides with the partition function Z Λ in (13), thus proving (18).
In view of Theorem 3.1, the rightmost limit in (14) coincides with the moment Lyapunov exponent (MLE) of the auxiliary GRF ξ: which is identical to (5) since the product moments of the GRF being considered are all nonnegative. This limit is completely specified by the mean µ, given by (20), and the spectral density function S : [−π, π] n → R + defined as the Fourier transform of the covariance function (16) by Here, u T v is the inner product in R n , vectors are organised as columns unless indicated otherwise, and (·) T denotes the transpose. Computing the MLE m(ξ), with M Λ (ξ) in (23) being replaced with |M Λ (ξ)|, is of considerable interest in its own right for a general homogeneous GRF ξ. However, the salient feature of this problem in the context of the EDM system is that, in view of (16), both the mean value and the covariance function of the auxiliary GRF in (17) are nonnegative. The nonnegativeness of C implies that the spectral density S in (24) is not only symmetric and nonnegative, but is also positive definite. Hence, S is the covariance function of yet another homogeneous GRF, this time, on the ndimensional torus [−π, π) n , whose probability distributions are invariant under the group of entrywise additions modulo 2π.
Remark 1. The fact that the choice of W (0), which affects C(0) in (16), has no influence on the product moments of the auxiliary GRF ξ is explained by their invariance under adding a ξ-independent "white noise" GRF ζ := (ζ x ) x∈Z n whose values ζ x are independent Gaussian random variables with zero mean and variance γ 2 . The resulting GRF ξ + ζ := (ξ x + ζ x ) x∈Z n has spectral density S + γ 2 and its product moment over a finite set Λ ⊂ Z n is computed as Here, ξ Λ := (ξ x ) x∈Λ is the restriction of ξ to the set Λ, and, in extending (4), use is made of the conditional product moment M Λ (η | S) := E x∈Λ η x | S of a random field η := (η x ) x∈Z n over the set Λ with respect to a σ-subalgebra S (or a conditioning random element which generates S). In (25), we have also used the tower property of iterated conditional expectations [43], and the assumption that ζ consists of independent zero mean random variables, independent of ξ. Therefore, under the transformation µ → ρµ and S → ρ 2 S + γ 2 , which corresponds to ξ → ρξ + ζ with a constant factor ρ > 0, the MLE m(ξ) in (23), as a functional of µ and the spectral density S from (24), is transformed as m → m + ln ρ.

4.
The classical monomer-dimer problem and moving average Gaussian random fields. Consider the structure of the auxiliary GRF for the classical monomer-dimer problem [12,13,25] as a particular case of EDM systems. Suppose the dimer potential W on the set Z n \ {0} is given by where α k is the energy ascribed to dimers which are parallel to the kth coordinate axis, and is the kth basis vector in R n . Note that (26) forbids dimers {x, y} with |x − y| > 1 and describes the dependence of the energy of a dimer of Euclidean length 1 on its orientation. In particular, if V := +∞ (so that monomers are energetically forbidden) and α 1 = . . . = α n = 0, then the partition function Z Λ in (13) in this monomer-free isotropic case counts the number of "domino" tilings [3, p. 125] of the set Λ, that is, partitions of Λ into pairs of nearest neighbours, assuming that #Λ is even. In general, the dimer orientations with smaller values of α k are energetically favoured. The resulting series on the right-hand side of (15) consists of finitely many nonzero terms, and the auxiliary GRF ξ for this classical "rigid" dimer system can be explicitly constructed as follows. Using (27), for every k = 1, . . . , n, we denote by the set of midpoints of the edges between neighbouring sites in Z n parallel to the kth coordinate axis. For any point x ∈ Z n , the nearest sites of E k are x ± e k /2, both at the Euclidean distance 1/2 from x. Accordingly, those midpoints from the combined set which are nearest to a given x ∈ Z n , are x ± e k /2, with k = 1, . . . , n. We will now associate independent standard normal random variables η y with the sites y ∈ E, and ascribe, to every x := (x k ) 1 k n ∈ Z n , a weighted sum over the 2n sites of the set (29), nearest to x, where ρ 1 , . . . , ρ n are nonnegative parameters given by The resulting GRF ξ := (ξ x ) x∈Z n is homogeneous and has mean µ and covariance function which, in view of (31), is related with the dimer potential (26) by (16), with W being extended to the origin as W (0) = −T ln(2 n k=1 ρ k ). Therefore, (30) indeed produces an auxiliary GRF ξ for the classical dimer system. It is a moving average random field which is easy to simulate. The simulation boils down to generating the standard normal random variables at sites of the set E given by (28)- (29). This, in principle, provides an alternative probabilistic way for approximate computation of the partition function Z Λ for the monomer-dimer problem over a finite set Λ ⊂ Z n by modelling the auxiliary GRF ξ and statistically estimating its product moment M Λ (ξ) instead of the Monte-Carlo simulation of the underlying particle system which is more complicated.
In view of (32), the covariance matrix cov(ξ Λ ) of the restriction of ξ to a finite set Λ ⊂ Z n , with ξ Λ being considered as a Gaussian random vector of dimension #Λ, has a (2n+1)-diagonal structure. A similar sparsity structure is known for the inverse covariance matrices (that is, precision matrices) of finite fragments of Markov GRFs [40]. The Markov GRFs describe the equilibrium states of systems of linear harmonic oscillators with nearest neighbour interaction and are, in general, hard to simulate in contrast to the moving average random fields above; see also [4,16]. An important subclass of Markov GRFs which is amenable to efficient unilateral simulation in a single pass through the simulation domain is provided by Pickard random fields (PRFs) [6,7,22,37,38,45] which behave in space like autoregressive processes do in time. 5. Manhattan dimers and Gaussian Pickard random fields. Suppose the dimer potential is given by a weighted 1 -norm (also called "Manhattan norm") of the dimer vector z := (z k ) 1 k n ∈ Z n : Here, α 1 , . . . , α n are positive parameters which quantify the projections of the force of attraction between the particles in a dimer onto the coordinate axes. Unlike ideal springs, the hypothetical Manhattan dimer does not develop a linearly increasing force when being stretched, albeit the attraction between the particles in it persists. The covariance function (16), associated with the Manhattan dimer potential (33), is where 0 < ρ 1 , . . . , ρ n < 1 are computed according to (31). The resulting auxiliary GRF ξ is a Gaussian PRF (GPRF) in dimensions two [45] and three [22] and lends itself to the unilateral simulation mentioned above. Such simulation of the GPRF ξ is made possible by the following enhancement of the spatial Markov property for any x := (x k ) 1 k n ∈ Z n : Here, [[[η | ζ]]] is a shorthand notation 3 for the conditional probability distribution of η given ζ, and The sets Λ − x and Λ + x are the "past" and "future" of the lattice with respect to the site x in the spatial sense [7]. Accordingly, the set Λ 0 x , which consists of 2 n − 1 points from Λ − x , plays the role of the nearest past for x; see x and the spatial future Λ + x for a given site x ∈ Z 2 , defined by (36), are represented by "•"s and "•"s, respectively. The sites of the nearest past Λ 0 x are marked by " "s.
which we will study in Sections 8-12, the relation (35) between the conditional probability laws reduces to the standard Markov property, thus allowing ξ (which becomes a Gaussian random sequence) to be generated by an autoregressive equation.
A remarkable feature of the Markov structure of the GPRF in the multivariate case n > 1 is that its restrictions ξ L to the level sets L 0 , L 1 , L 2 , . . ., described below, form a Markov chain of order n. This property is employed by the level-sets algorithm of [8] for the unilateral simulation of such random fields on the nonnegative orthant Z n + of the lattice, where Z + is the set of nonnegative integers. For example, in the two-dimensional case n = 2 elucidated by Fig. 2, the level sets are given by and the restrictions ξ L of the GPRF, which are organised as Gaussian random vectors (with ξ L having dimension + 1), form a Markov chain of order two. More precisely, for any > 0, the conditional probability distribution . Furthermore, the values of ξ at sites of the next level set L +1 are conditionally independent given ξ L −1 and ξ L , with their conditional joint probability distribution being the product of + 2 conditional Gaussian distributions. Here, ξ jk at site (j, k) is conditioned on the values of ξ in the nearest past Λ 0 jk which is restricted to the nonnegative quadrant Z 2 + of the lattice as in conformance with (36) for the two-dimensional case. In terms of the centered random field σ := (σ x ) x∈Z 2 , defined by with µ associated with the monomer energy by (20) as before, the conditional Gaussian distributions on the right-hand side of (38) are given by where ρ 1 , ρ 2 are related to the attraction force parameters α 1 , α 2 of the Manhattan dimer potential (33) by (31), and the conditional variances γ 2 1 , γ 2 2 of the last two distributions are computed as This follows from the well-known results on conditional distributions for jointly Gaussian random variables [43]. The coefficients a, b, c and the conditional variance d 2 of the Gaussian distribution for j, k > 0 in (41) are computed as with var(·) the variance of a random variable. Here, the conditioning restriction σ Λ 0 jk of the centered random field (40) to the set (39) is organized as the three-dimensional row-vector jk := σ j−1,k−1 , σ j,k−1 , σ j−1,k whose covariance (3 × 3)-matrix cov( jk ) and the row-vector cov(σ jk , jk ) of cross-covariances with σ jk are appropriate submatrices of the covariance matrix obtained by evaluating (34) for the square cell {j − 1, j} × {k − 1, k} of the lattice; see also [38,45].
6. Aggregation of level sets. We will now employ the technique of reducing the order of a Markov chain by augmenting its state space. For any > 0, consider the restriction of the centered GPRF σ, defined by (40), to the union of two adjacent level sets from (37) which consists of 2 +1 lattice sites. Also, let Θ 0 := L 0 be the singleton formed by the origin of Z 2 . The resulting sequence of random vectors σ Θ , labeled by 0, is a Markov chain of order one. In view of the reflection symmetry of the covariance function (34) with respect to the coordinate axes, the restrictions σ −Θ of σ to the sets −Θ , which are centrally symmetric to the aggregated level sets (45) about the origin, also form a Markov chain. Furthermore, the Markov property remains valid upon removal of the "endpoints" (0, − ) and (− , 0) from the set −Θ . More precisely, the restrictions σ Γ of the GPRF σ to the sets Γ 0 , Γ 1 , Γ 2 , . . . defined by also form a Markov chain. We label the values of σ at sites of the set Γ N by integers 0 to 2N as shown in Fig. 3, which allows σ Γ N to be identified with a (2N + 1)-dimensional Gaussian random vector Due to this numbering, Σ N is obtained by sampling the GPRF σ along a monotone path in the planar lattice. Hence, by the results of [38,45], for every given N 0, the entries  (34) for the zigzag shaped set Γ N , it follows that Σ N is an alternating autoregressive sequence [36] governed by the recurrence equation where the initial element Σ N,0 and the innovations w N,0 , w N,1 , w N,2 , . . . are independent standard normal random variables, and use is made of (42). For any N > 0, the conditional Gaussian distribution of Σ N −1 given Σ N is degenerate since these random vectors have common entries: The other entries of Σ N −1 (namely, Σ N −1,2k for 0 k < N ) are conditionally independent given Σ N . More precisely, where each of the conditional distributions on the right-hand side is described by the first of the Gaussian distributions from (41): involves at most three entries of the random vector Σ N for every 0 j 2N − 2. This property is inherited from the localness of the conditional probability measures describing the transitions of PRFs between the level sets.

7.
Recurrence relation for conditional product moments. Assuming N 0 and using the sets Γ 0 , . . . , Γ N from (46), we will now consider the two-dimensional Manhattan EDM system on the set which consists of #∆ N = N (N + 5)/2 + 1 lattice sites; see Fig. 3. The product moment of the auxiliary GPRF ξ over ∆ N is representable as Here, the function R N : R 2N +1 → R is defined as the conditional product moment of ξ over the set ∆ N with respect to the random vector Σ N from (47): In particular, since the set ∆ 0 is a singleton which consists of the origin of Z 2 , the definition (54) yields R 0 (s) = E(ξ 00 | σ 00 = s) = µ + s.
Theorem 7.1. For any N > 0, the conditional product moment function R N of the auxiliary GPRF ξ for the two-dimensional Manhattan EDM system satisfies the recurrence equation with the initial condition (55). Here, the conditional expectation is taken over the condi- Proof. By partitioning the set ∆ N from (52) into two disjoint subsets ∆ N −1 and ∆ N \ ∆ N −1 = Γ N \ ∆ N −1 , with Γ N given by (46) (see also Fig. 3), it follows that the function (54) can be factorized as Here, is a polynomial function of the appropriate entries of the random vector Σ N in (47). Furthermore, where, in addition to the tower property of iterated conditional expectations, we have used which follows from the measurability of ξ ∆ N −1 with respect to Σ 0 , . . . , Σ N −1 and the Markov property of the sequence Σ 0 , . . . , Σ N . The recurrence equation (56) is now obtained by substitution of (58) and (59) into (57).
By induction, (55) and (56) imply that the function R N is a polynomial (in 2N + 1 scalar variables) of degree #∆ N . Therefore, calculating the partition function of the twodimensional Manhattan EDM system for the set ∆ N through the product moment of the auxiliary GPRF ξ according to (53) consists in evaluating the expectation of the recursively computed polynomial R N on the segment Σ N of the alternating autoregressive sequence governed by (48). For example, by recalling the conditional distribution (41)- (44), it follows that is a quartic polynomial. Further iterations of the recurrence relation (56) involve the conditional expectations of higher order polynomials applied to conditionally independent Gaussian random variables; cf. (50), (51). Although this recurrence is not easy to implement, its significance is that the above discussed localness of the conditional probability , which specifies the linear operator R N −1 → R N in (56), appears to be a striking reduction of the long-range nature of interactions in the underlying Manhattan EDM system. 4 Remark 3. The "localization" of the original problem above has been achieved through its connection with the product moments of the auxiliary GPRF and the possibility of their recursive computation (in terms of the conditional product moments) due to the specific Markov structure of the random field. A similar line of reasoning, employing the aggregated level sets and associated conditional product moments, is applicable in three dimensions for the trivariate GPRFs from [22]. This suggests that such an indirect approach to the monomer-dimer problem and its extensions may contain a hidden resource which is worth exploring. 8. One-dimensional Manhattan dimer-monomer system. Consider the Manhattan EDM system in the one-dimensional case, where the dimer potential (33) reduces to W (z) := α|z| for z ∈ Z, with α > 0 the attraction force parameter. In this case, the auxiliary GRF, which becomes a sequence ξ := (ξ k ) k∈Z , is a homogeneous Gaussian Markov chain with mean µ from (20) and covariance function cov(ξ j , ξ k ) = ρ |j−k| , where, in accordance with (31) and (34), ρ := e −αβ , (60) and the shorthand notation (12) is used. In view of (23), the MLE of the sequence ξ takes the form where is the product moment over the discrete interval 0 : N − 1, with a : b denoting the set Z [a, b] of integers from a to b. For convenience of calculations, we will use a centered sequence σ := (σ k ) k∈Z , defined by σ k := ξ k − µ and satisfying an autoregressive equation for k 0. Here, the initial condition σ 0 and the innovations w 0 , w 1 , w 2 , . . . are independent standard normal random variables, and The autoregressive sequence σ can be obtained by sampling an Ornstein-Uhlenbeck diffusion process Y := (Y t ) t∈R+ with time step := αβ as σ k := Y k . The process Y is governed by an Itô stochastic differential equation dY t = −Y t dt + √ 2dW t , where W := (W t ) t∈R+ is a standard Wiener process, independent of the standard normal initial value Y 0 . 9. Conditional product moments and transfer operator. For the one-dimensional Manhattan EDM system of the previous section, the product moment (62) of the auxiliary Gaussian Markov chain ξ is representable as Here, for any N 0, the function Q N : R → R is defined by the conditional product moment with Q 0 ≡ 1. Note that Q N is a univariate counterpart of the function R N associated with the two-dimensional Manhattan EDM system by (54). Thus, by employing a line of reasoning similar to Theorem 7.1, and using the Markov property and homogeneity of the sequence ξ, it follows that the function Q N +1 is expressed in terms of Q N from (66) through a linear transfer operator G as Here, the last expectation is taken over the transition probability law of the sequence σ which, in view of (63) and (64), is described by that is, the conditional Gaussian distribution with mean E(σ 1 | σ 0 ) = ρσ 0 , variance var(σ 1 | σ 0 ) = γ 2 and PDF By induction, (67) and (68) imply that the function Q N is a polynomial of degree N with the leading coefficient The latter also employs the property that if X ∼ N (m, s 2 ) with a fixed variance s 2 , and L is a polynomial, then EL(X) is a polynomial in the mean value m with the same leading term as L. In particular, the first three of the polynomials Q N are computed as which shows that the recursive calculation of their coefficients, which are needed for (65), in the standard basis of monomials 1, x, x 2 , . . . quickly becomes complicated. In Section 11, we will perform these calculations in a more suitable basis of Hermite polynomials after establishing the boundedness of the transfer operator G which governs the recurrence equation (67).

Remark 4.
A similar recursive computation of conditional product moments can also be developed for a wider class of homogeneous Gaussian random sequences with rational spectral densities which are realizable as a hidden Markov chain with a higher dimensional state space.
10. An upper bound for the moment Lyapunov exponent. Let H denote the Hilbert space of functions h : R → R which are square integrable over the standard normal PDF ν from (8) and are endowed with the norm h := h, h and inner product f, g := , where X is a standard normal random variable. The Hermite polynomials (9) form an orthogonal basis of the space H, with H k = √ k!; see, for example, [24, pp. 19-20]. The following theorem establishes an upper bound for the · -induced norm of the transfer operator G defined by (67): Theorem 10.1. The norm (70) of the transfer operator G in (67), associated with the onedimensional Manhattan EDM system, satisfies Proof. Let h ∈ H. Then, by employing a change of measure technique, it follows that for any real x, Here, with p 1|0 the transition PDF of the sequence σ described by (69); ν is the standard normal PDF from (8), and p 01 (x, y) := e (2ρxy−x 2 −y 2 )/(2γ 2 ) (2πγ) (74) is the joint PDF of σ 0 and σ 1 . Application of the Cauchy-Bunyakovsky-Schwarz inequality to (72) yields (75) An upper bound for the transfer operator norm (70) can now be obtained by integrating both sides of (75) with respect to the standard normal distribution and dividing the result by h 2 : Here, in view of (73), (74) and (64), s(x, y) := (q(x, y)) 2 ν(x)ν(y) = (p 01 (x, y)) 2 ν(x)ν(y) The function γ 2 s : R 2 → R + is recognizable as a bivariate Gaussian PDF with zero mean and covariance matrix whose determinant is equal to one and the diagonal entries describe the common marginal variance (1 + ρ 2 )/γ 2 . It is the latter quantity that enters the right-hand side of (76) through (77) as thus proving (71). Now, by using the Cauchy-Bunyakovski-Schwarz inequality and the submultiplicativity of the operator norm ||| · |||, it follows that the product moment (65) satisfies for any N 0, where 1 denotes the identically unit function. Hence, the MLE (61) admits an upper bound m(ξ) ln |||G|||, which, in combination with (71) and (64), yields Remark 5. The last term in (79) is recognizable as Shannon's mutual information I(σ 0 , σ 1 ) between the random variables σ 0 and σ 1 (see [9,17]), which, in the Markov case being considered, coincides with I(σ 1 , σ 0 ), where σ 0 := (σ k ) k 0 is the past history of the sequence σ up until time 0.
Remark 6. The proof of Theorem 10.1 shows that the inequality (76) remains valid for any homogeneous Markov random sequence ξ, and hence, its MLE, in this more general case, admits an upper bound where ν and p 01 are understood as the invariant PDF of ξ and the joint PDF of ξ 0 and ξ 1 , respectively.
11. Matrix representation in the basis of Hermite polynomials. As an orthonormal basis of the Hilbert space H from Section 10, we will use the functions h 0 , h 1 , h 2 , . . . obtained by scaling the Hermite polynomials (9) as Their significance for computing the product moment (62) is that, in view of the leftmost equality from (78), where q N,k := Q N , h k (82) are the Fourier coefficients of the polynomial Q N with respect to the orthonormal basis (80), so that We will show that the coefficients q N,0 , . . . , q N,N satisfy a recurrence equation in N 0 which can be "encoded" in an equivalent matrix representation of the transfer operator G defined by (67). To this end, we will employ a lemma below which, although it follows from the properties of the Ornstein-Uhlenbeck operator [32,Proposition 1.5.4(v) on p. 233], is provided with a proof for completeness of exposition.
Lemma 11.1. Let X be a Gaussian random variable with mean m and variance s 2 < 1. Then, for any nonnegative integer k, Proof. By combining the generating function (10) Both A and A † are defined (and are mutually adjoint) on the linear manifold of those functions h ∈ H whose Fourier coefficients in the basis (80) are square summable with an appropriate weight: k 1 k h, h k 2 < +∞. The action of the operators A and A † in (85) on the basis functions is identical to that of the annihilation and creation operators on the eigenfunctions of the Schrödinger equation for the quantum harmonic oscillator; see, for example, [35, pp. 52-53] and [42, p. 91]. In this sense, A + A † (which is a symmetric operator) corresponds to the quantum-mechanical position operator [24, pp. 210-211].
Although A and A † are unbounded, their compositions AR and A † R with the operator R, as well as R itself, are bounded operators on the whole space H, since 0 < ρ < 1 in view of (60).
Theorem 11.2. The transfer operator G in (67) can be expressed in terms of the operators R, A and A † from (85), and its matrix representation in the basis (80) is given by Proof. In view of (68), application of Lemma 11.1 with m := ρσ 0 and s := γ from (64) yields the identity by which the random sequence (ρ −jk H k (σ j )) j 0 is a martingale with respect to the natural filtration of the homogeneous Gaussian Markov chain σ for any given k 0. Hence, the action of the transfer operator G on the kth Hermite polynomial is described by where we have also used the recurrence relation H k+1 (x) = xH k (x) − kH k−1 (x). Upon division of both parts of (87) by √ k!, the result is representable in terms of the basis functions (80) and the operators (85) as Since the last relation holds for any k 0, the transfer operator G is indeed expressed in terms of R, A and A † as in (86). The first of the equalities (88) yields the kth column of the matrix representation of G in the basis (80) on the right-hand side of (86).
Since 0 < ρ < 1, the matrix representation (86) implies that G is a compact operator. With √ R denoting the square root of the operator R from (85) defined by √ R(h k ) = ρ k/2 h k , the spectrum of G coincides with that of a self-adjoint operator √ R(µI + A + A † ) √ R and is all real. Moreover, since µ > 0, the convex cone H + of functions h ∈ H with all nonnegative Fourier coefficients h, h k with respect to the orthonormal basis (80) is invariant under G. Hence, by the Ruelle-Perron-Frobenius theorem, the operator G has a unique (up to a positive scalar factor) eigenfunction in H + , which corresponds to the largest eigenvalue λ of G. By applying (86) to the Fourier coefficients (82) of the polynomial Q N , it follows that for any N 0, they satisfy the recurrence equation Here, q N,k = 0 for k > N and k < 0 in view of (83), and the recursion is initialized by q 0,0 = 1 and q 0,k = 0 for any k = 0. Since the cone H + is invariant under G, a straightforward induction yields the nonnegativeness q N,k 0. Furthermore, the ratio of the free terms (81) of the polynomials Q N and Q N +1 converges to the largest eigenvalue of the transfer operator G as lim which can be used for finding this eigenvalue numerically as it is done in the power method. The dependence of the MLE m(ξ) = ln λ on µ and ρ, obtained by numerical computation of the eigenvalue λ, is presented in Fig. 4.
where the second sum consists of finitely many nonzero terms (with k N ), and X is a standard normal random variable. Here, use has also been made of the generating function of Hermite polynomials (10). The definition of T N implies that it is a polynomial of degree N .
13. Concluding remarks. We have established a novel connection between the product moments and moment Lyapunov exponents of homogeneous Gaussian random fields on a multidimensional integer lattice with the partition function of an elastic dimer-monomer system which extends the classical monomer-dimer problem of equilibrium statistical mechanics by allowing for long-range interactions. A salient feature of the auxiliary GRF, which encodes the thermodynamic properties of the EDM system, is that its mean value and the covariance function are the Boltzmann factors, associated with the energetics of the underlying particle system, and, therefore, are both nonnegative. Any insight into how the MLE of such a GRF depends on its spectral density function in an arbitrary dimension would shed light on the solution of the EDM problem, including its particular case, the classical monomer-dimer problem, which remains unsolved in dimensions three and higher.
To this end, we have outlined a research paradigm which builds on the observation that a specific spatial Markov structure of the GRF, present, for example, in Gaussian Pickard random fields, facilitates recursive representation of the product moments of such GRFs. We have discussed the possibility to recursively compute the product moments with a local transition from one level set to another, for a Gaussian PRF, which is the auxiliary GRF for the Manhattan EDM system. This paradoxical "localization" of the effect of longrange interactions suggests that the product moment approach, proposed in the present study, is a promising resource towards the solution of the monomer-dimer problem and its generalizations.
In addition to its links with the statistical mechanical lattice models of interacting particle systems, the problem of computing the product moments and MLEs for homogeneous GRFs appears to be of interest in its own right from the probability theoretic point of view. In this regard, even the deceptively simple case of univariate Gaussian Markov chains reveals rich algebraic and probabilistic structures. In particular, we have linked the product moments and MLEs in this case to the ladder operators of the quantum harmonic oscillator and the pantograph functional-differential equation by using a change of measure technique.