Higher order fluctuation fields and orthogonal duality polynomials

Inspired by the works in [1] and [8] we introduce what we call $k$-th-order fluctuation fields and study their scaling limits. This construction is done in the context of particle systems with the property of orthogonal self-duality. This type of duality provides us with a setting in which we are able to interpret these fields as some type of discrete analogue of powers of the well-known density fluctuation field. We show that the weak limit of the $k$-th order field satisfies a recursive martingale problem that formally corresponds to the SPDE associated with the $k$th-power of a generalized Ornstein-Uhlenbeck process.


Introduction
In the context of interacting particle systems with a conserved quantity (such as the number of particles) in [5,9] one studies the time-dependent density fluctuation field X (n) (ϕ, η(n 2 t)) = 1 n d/2 x∈Z d ϕ(x/n)(η x (n 2 t) − ρ).
Here ϕ denotes a test-function, and η x the number of particles at site x ∈ Z d . The quantity X (n) (ϕ, η(n 2 t)) is then considered as a random time-dependent (Schwartz) distribution. In a variety of models with particle number conservation (such as zero range processes, simple exclusion process, etc.), this time-dependent field is proved to converge, at equilibrium, to a stationary infinite dimensional Ornstein-Uhlenbeck process. This scaling limit behavior of the density fluctuation field can be thought of as a generalized central limit theorem, and, as such, as a correction or refinement of the hydrodynamic limit, which can be thought of as a law of large numbers.
The usual strategy of proof (see e.g. Chapter 11 of [9]) is to start from the Dynkin martingale associated to the density field and prove convergence of the drift term via the Boltzmann-Gibbs principle (the drift term becomes in the scaling limit a function of the density field), and convergence of the noise term via characterization of its quadratic variation (which becomes deterministic in the scaling limit). This then eventually leads to the informally written SPDE dX t = D∆X t + σ(ρ)∇dW t where ρ is the parameter of the invariant measure associated to the density, ∆ denotes the Laplacian, and where σ(ρ)∇dW t is an informal notation for Gaussian white noise with variance σ 2 (ρ) (∇ϕ) 2 dx.
In interacting particle systems with (self-)duality, the drift term in the equation for the density field is already microscopically (i.e., without rescaling) a (linear) function of the density field. As a consequence, closing the equation and proving convergence to the limiting Ornstein-Uhlenbeck process, is, for self-dual systems, particularly simple and don't require the use of a Boltzmann-Gibbs principle. This simplification suggests that, in that context, we can obtain more detailed results about fluctuation fields of more general observables. Orthogonal polynomial duality is a useful tool in the study of fluctuation fields, and associated Boltzmann-Gibbs principles, as we have seen in [2].
The density fluctuation field can be viewed as the lowest (i.e., first) order of a sequence of fields associated to orthogonal polynomials. Indeed, in all the models with orthogonal polynomial self-duality, the function (η x − ρ) is the first order orthogonal polynomial up to a multiplicative constant. Orthogonal polynomials are indexed by finite particle configurations, i.e., the dual configurations. If we denote by D(x 1 , . . . x k ; η) the orthogonal polynomial associated to the dual configuration n i=1 δ x i , then a natural field generalizing the density fluctuation field is X (n,k) (Φ, η) = n −kd/2 x i ∈Z d D(x 1 , . . . , x k ; η) · Φ x 1 n , . . . , x k n .
In the context of exclusion process the case k = 2 (orthogonal polynomial of order 2) has been studied in [8], where this field, called the quadratic fluctuation field, is shown to converge, in the limit n → ∞, to the solution of a martingale problem. The quadratic variation of this 2nd-order field is proven to be a function of the 1st order field (the density field). From the result on the quadratic (k=2) field one can conjecture the existence of a more general structure where the kth-order orthogonal polynomials field satisfies, in the scaling limit, a martingale problem with quadratic variation depending on the k − 1-order field. In this paper we show exactly the emergence of a scenario of this type: within a general class of models with orthogonal polynomial self-duality we consider the fluctuation fields associated to orthogonal polynomials and prove that they converge, in the scaling limit, to the solution of a recursive system of martingale problems. We believe that this can also be a first step in the direction of defining non-linear fields, such as the square of the density field, via approximation of the identity, i.e., via a singular linear observable (cf. [8]) of the field constructed in our paper.
The rest of our paper is organized as follows. In Section 2 we define the basic models, and introduce orthogonal polynomial duality. In Section 3 we define the fluctuation fields, in Section 4 we introduce a coordinate version of the dual process, a technical tool that will prove to be useful later on. In Section 5 we state the main result, Theorem 5.1 below, and outline a strategy of its inductive proof. Finally, the rest of the sections are devoted to the proof of Theorem 5.1.

The models 2.1 The infinite configuration process
We consider an interacting particle system where an infinite number of particles randomly hop on the lattice Z d . Configurations are denoted by η, ξ, ζ and are elements of Ω ⊆ N Z d (where N denotes the natural numbers including zero). We denote by η x the number of particles at x in the configuration η ∈ Ω. We have in mind symmetric processes of the type independent random walkers, inclusion or exclusion. We fix two parameters (σ, α) ∈ {0, 1} × [0, ∞) ∪ {−1} × N and we define the generator working on local functions f : Ω → R as where η i,i+r denotes the configuration obtained from η by removing a particle from i and putting it at i + r. The state space Ω has to be defined and its form depends on the choice of the parameters α and σ. We assume that p(r) is a symmetric, finite range, irreducible Markov transition function on Z d : 1. Symmetry. The function p : R d → [0, ∞) is of the form: p(r 1 , . . . , r d ) = p(|r 1 |, . . . , |r d |) and such that p(r σ(1) , . . . , r σ(d) ) = p(r 1 , . . . , r d ) for all σ ∈ P(d), the set of permutations of {1, . . . , d}.
2. Finite-range. There exists a finite subset R ⊂ Z d of the form R = [−R, R] d ∩ Z d , for some R ∈ N, R > 1, such that p(r) = 0 for all r / ∈ R.
3. Irreducibility. For all x, y ∈ Z d there exists a sequence i 1 = x, . . . , i n = y such that We will also assume, without loss of generality, that p(0) = 0, and denote by χ the second moment: For the associated Markov processes on Ω, we use the notation {η(t) : t ≥ 0}, η x (t) denoting the number of particles at time t at location x ∈ Z d . These particle systems have a one-parameter family of homogeneous (w.r.t. translations) reversible and ergodic product measures ν ρ , ρ > 0, indexed by the density of particles, i.e., The nature of the underlying dynamics and the type of reversible measure we obtain is regulated by the parameter σ ∈ Z as follows.
Independent random walkers (IRW): This particle system corresponds to the choice σ = 0 and the intensity parameter α ∈ R regulates the rate at which the walkers move. The reversible measures ν ρ , ρ > 0 are products of Poisson distributions with parameter ρ, ν ρ = ⊗ i∈Z d Pois(ρ), i.e. the marginals are given by Symmetric exclusion process (SEP(α)): The choice σ = −1 results in exclusion interaction. For this process the parameter α takes values in the set of natural number, α ∈ N, as it determines the maximum number of particles allowed per site. This system is well known to have reversible measures ν ρ , ρ ∈ (0, α), that are products of Binomial distributions: ν ρ = ⊗ i∈Z d Binom α, ρ α whose marginals are given by Symmetric inclusion process (SIP(α)): The choice σ = 1 gives rise to an interaction of inclusion-type consisting of particles attracting each other. The SIP is known to have products of Negative-Binomial distributions as reversible measures, i.e. ν ρ , ρ > 0 with The definition of the state space Ω is different on each case, depending on whether there are restrictions or not on the total number of particles allowed per site. This is finite for the exclusion process, thus, for SEP(α), we have Ω = {0, 1, . . . , α} Z d . The situation is different in the cases of IRW and SIP, for which, in principle, there are no restrictions. Nevertheless, one has to avoid explosions of the number of particles in a given site. For this reason the characterization of Ω in these cases (i.e. for σ ≥ 0) is a more subtle problem whose treatment is beyond the scope of this paper. Here we will restrict to implicitly define Ω as the set of configurations in N Z d whose evolution η(t) is well-defined and belonging to Ω for all subsequent times t ≥ 0. For what we are concerned, a precise characterization of Ω is not necessary as we will mostly work in L 2 (ν ρ ) spaces, and the set Ω (implicitly defined above) has ν ρ measure 1 for all ρ.

The finite configuration processes
The process introduced in Section 2.1 can also be realized with a fixed finite number of particles. For a process with k ∈ N particles we denote by Ω k its state space, more precisely: We will then denote by {ξ(t) : t ≥ 0} the Ω k -valued Markov process, with infinitesimal generator given by working on functions f : Ω k → R.
We now define the measure Λ(·), which is a product measure not depending on k: , m ∈ N for σ = 1 SIP(α) Notice that by detailed balance we can verify that the measure Λ(·) is reversible. As a consequence of this reversibility we can infer that the k-particles generator L (k) is self-adjoint with respect to the inner product ·, · Λ , i.e. for all f, g ∈ L 2 (Ω k , Λ) we have where the inner product is given by

Orthogonal polynomial self-duality
The processes defined in Section 2.1 share a self-duality property that will be crucial in our analysis. Define the set of configurations with a finite number of particles, the self-duality functions that we consider in this paper are functions D ρ : Ω f × Ω → R parametrized by the density ρ > 0 satisfying the following properties.
3. Orthogonality: For proofs of self-duality with these functions we refer to [7] and [11]. REMARK 2.1. Notice that, as a consequence of the orthogonality property (13), we have that where p t (·, ·) is the transition probability function of the dual process {ξ(t) : t ≥ 0}. Moreover, if we use the reversibility of the measure ν ρ on the LHS of (15) we obtain which, by detailed balance, implies the reversibility of the measure µ ρ (ξ). This in turn implies that there exists a constant c(k, ρ) such that From now on we will often suppress this dependence on the parameter ρ, D(·, ·) = D ρ (·, ·), in order not to overload the notation. For each of the processes we are considering, the orthogonal duality polynomials are given as follows.
IRW: Charlier polynomials. The duality polynomials are given by where C (m, ·) is the Charlier polynomial of degree m that we characterize by means of its generating function: SEP(α): Krawtchouk polynomials. Strictly speaking these polynomials do not satisfy a selfduality relation. However, under a proper normalization we can find a duality function in terms of them. The duality polynomials are given by where K(m, ·) is the Krawtchouk polynomial of degree m whose generating function is SIP(α): Meixner polynomials. In this case the polynomials satisfying the self-duality relation are given by the following normalization of the Meixner polynomials where M (m, ·) is the Meixner polynomial of degree m with generating function We refer the reader to [10] and [4] for more details on these polynomials and their generating functions.

Fluctuation fields
The density fluctuation field Y is the stochastic object usually defined to study small deviations from the hydrodynamic limit. This field corresponds to a central limit type of rescaling of the density field, i.e. X (n) Fields of this type have been intensively studied in the literature. For different models, the sequence X (n) t is proven to converge to a limiting field X t that is identified as the distributionvalued random variable satisfying the following martingale problem: for any ϕ ∈ S(R d ) the process is a square integrable continuous martingale of quadratic variation given by the expression being also a martingale.
Formally speaking, (23) and (24) imply that the limiting field X t satisfies the Ornstein-Uhlenbeck equation We refer the reader to [5] for a precise statement on the convergence for the case of the exclusion process, corresponding, in our setting, to the case α = 1 and σ = −1.
The density field (22) can be written, in our context, in terms of our orthogonal polynomial dualities D ρ (ξ, η) by choosing ξ ∈ Ω 1 . Indeed, in all models considered we have that there exists a constant c σ,α,ρ such that Later on, in order not to overload notation we will suppress the dependence on ρ and α and just write c σ . From (26) we observe that the field (22) can be re-written (modulo a multiplicative constant) as X where the superindex (n, 1) suggests that, in some sense, this is the first-order density field.
Using (26) and (25) the formal limiting SPDE for X t is The observation that the field (22) can be expressed in terms of duality polynomials opens the possibility of defining higher order fields and study their scaling limits. For k ∈ N, k ≥ 1 we define the k-th order field as is a test function, Λ is as in (7), and Notice that the only difference between X (n,k) (ϕ (k) , η) and Y (n,k) (Φ, η) is that the latter works on test functions over configuration space, i.e., Φ ∈ S(Ω k ), while the former works on test functions ϕ (k) ∈ S(R kd ). Then, using the notation we can rewrite the k-th order field (30) as and define: Y The choice of multiplying the duality function by the measure Λ(·) in (33) is dictated simply by computational convenience that, even if obscure at the moment, will be made clearer in the course of the paper.

The coordinate process
Thinking of k ∈ N as the number of particles in our process, we want to introduce a family of permutation-invariant coordinate processes {X (k) (t) : t ≥ 0} compatible with the finite configuration processes {ξ(t) : t ≥ 0} on Ω k . Here the coordinate process is a Markov process on Z dk with X i (t) being the position of the i-th particle at time t ≥ 0. For a further explanation of the notion of compatibility we refer the reader to [3].
Denote by x ∈ Z kd the coordinate vector x : is defined by means of its infinitesimal generator: where x i,i+r denotes x after moving the particle in position x i to position x i + r ∈ Z d . Notice that for x ∈ Z kd the compatible configuration ξ(x) ∈ Ω k is given by

Product σ-finite reversible measures
It is possible to verify, by means of detailed balance, that the coordinate-process {X (k) (t) : t ≥ 0} admits a reversible σ-finite measure that is given by then we can rewrite Π in the product form: with π given as follows: Given the measures Π, we now consider the spaces of permutation-invariant functions: with P(k) denoting the set of all possible permutations of the set {1, 2, 3, . . . k}. We endowed the spaceL 2 (Z kd , Π) with the inner product given by: Notice that any function f ∈L 2 (Z kd , Π) can be interpreted also as a function on the configuration space. In this work we will extensively use this fact by changing between interpretations sometimes from one line to another in the same derivation.
As a consequence of reversibility of the measures Π, we can infer that the kparticles generator L (k) is self-adjoint with respect to the inner product ·, · Π , i.e.

The fluctuation fields in coordinate notation
It is possible to rewrite the fluctuation field (30) in the coordinate variables. Notice that in this context the test function Φ defined in (32) becomes a tensor function: i.e. it is the homogeneous k-tensor test function ϕ (k) ∈ S(R kd ) of the form then, after a change of variable in the sum we can rewrite the k-th field as follows Notice that we can also let the field X act on a general f ∈ S(R kd ) as expected, i.e., Because we deal with unlabeled particle systems it is natural to define the higher order fluctuation fields acting on symmetric test functions Φ i.e. on elements of the Schwartz space S(R kd ) that are permutation-invariant: Φ(x σ(1) , . . . , x σ(k) ) = Φ(x 1 , . . . , x k ) for all σ ∈ P(k), the set of permutations of {1, . . . , k}. By polarization and the fact that the orthogonal polynomials are functions of the unlabeled particle configurations, it is therefore sufficient to restrict our analysis to the set of symmetric tensor products of S(R d ). Indeed the linear combinations of those generate the symmetric test functions, see [12].

Heuristics: macroscopic dynamics
The goal of this section is to provide some intuitions on the type of limiting field that we should expect for fields of order greater than one. We will start by considering the cases k = 1, 2 and, inspired by the results obtained in [8], we will propose a heuristic interpretation of the two SPDEs obtained as scaling limits and their relation. Based on this interpretation we will conjecture a possible generalization to the kth-order case. In Section 5.2 we will give the rigorous result confirming the validity of the conjecture.
Here we will informally use the notation Y (k) t and X (k) t for the distributional limits of Y (n,k) and X (n,k) respectively.
Recall that from (29) we know that formally the distribution valued first order field X where for x ∈ R d , W t (x) is a space-time white noise. Additionally, from the martingale problem given in [8], we can deduce that the distribution-valued second-order field X where W t (x) is the white noise in (51) and ∆ (2) is the R 2d dimensional Laplacian, which is the sum of the Laplacian in the x variable plus the Laplacian in the y variable.
The key idea to extrapolate these relations to higher orders is to interpret the non-linearity on the RHS of (52) as a fields product, that we denote by ⋄, that satisfies the Leibniz rule of differentiation. This interpretation suggests that the second-order field X (2) t (x, y) is, in turn, a second power of the first-order field X (1) t (x). More precisely conjecturing since the product ⋄ follows the Leibniz rule we would have that which indeed agress with (52).
After the discussion above, it seems natural to expect that the kth-order field is a kth ⋄-power of the first order one. More precisely we conjecture that a relation of the type is satisfied. If this holds true, computations analogous to (53) would imply the formal SPDE dX (k) where ∆ (k) is the R kd -dimensional Laplacian, defined as the sum of the Laplacians at each coordinate and x −j is the (k −1)d-dimensional vector obtained from x by removing its coordinate x j .
The equation that we thus informally obtained will be later on justify in the main result of the next section.

Main theorem
Let us spend one paragraph to introduce the probability notions which are relevant for our main result. As we already mentioned, the kth-order fluctuation field can be considered as taking values in S ′ (R k ), the space of tempered distributions which is dual to S(R k ). Our original process η n 2 t has state space Ω (n) corresponding to the rescaled lattice 1 n Z. We then denote by P n , respectively E n , the probability measure, respectively expectation, induced by the measure ν ρ and the diffusively rescaled process η n 2 t on D([0, T ]; Ω (n) ). Finally, we denote by Q   Recursive martingale problem: for any symmetric ϕ (k) ∈ S(R kd ) the process is a continuous square integrable martingale of quadratic variation given by the solution of (51). REMARK 5.1. This recursive martingale problem is the rigorous counter part of the formal SPDE (54) that we heuristically obtained.

Strategy of the proof
We will show Theorem 5.1 by using induction on k. In the proof we will take advantage of the fact that the base case, k = 1, is already proved in the literature. On the other hand, the inductive step will be proven by means of an approach based on the natural Dynkin martingales: and where Γ is the so-called carré-du-champ operator given by: Notice that the Dynkin martingales can also be expressed in terms of the fields X (n,k) t .
Roughly our approach consists of the following steps: 1. we express the integrand term of equation (57) in terms of the kth-order fluctuation field Y (n,k) using duality (Section 6.1); 2. we close the equation (58) by expressing the integrand in the RHS in terms of the (k−1)thorder fluctuation field Y (n,k−1) (Section 6.2); 3. we show tightness for the sequence of probability measures Q

Inductive argument
The proof is done by induction over the order of the field k. The base case k = 1, corresponding to the density fluctuation field (22), is assumed to be true. Indeed, as mentioned in Section 3, a proof of Theorem 5.1 for exclusion dynamics and zero-range processes (of which independent random walkers are a particular case) is given in [5] and [9] respectively. By similar arguments the result can be extended to the case of inclusion process.
To implement the inductive argument we formalize the following inductive hypothesis that will be referred to several times in the course of the proof of Theorem 5.1. Martingale problem: for any symmetric ϕ (k 0 ) ∈ S(R k 0 d ) the process is a continuous square integrable martingale of quadratic variation 6 Proof of Theorem 5.1

Closing the equation for the drift term: k ≥ 2
In order to close the equation for the drift term thanks to Remark 4.2 we can just proceed as follows We proceed evaluating the action of the k-particles generator on Φ n . We then have (62) REMARK 6.1. Notice that the contribution coming from the second term in the RHS of (62) does not appear in the case k = 1.
First of all we prove that To prove this we use the Taylor expansion: and then

Now we have
It is now convenient to pass to the coordinate notation to treat sums of the type: for some ψ ∈ S(R d ). First of all we notice that summing over ξ ∈ Ω k is the same as summing over x ∈ Z kd : where the last identity follows using the expression of the field acting on more general (i.e., non-symmetric) test functions (50). Then, substituting in (69) we get where we used the fact that ϕ is uniformly bounded on Z. From this we can see that it is possible to close the equation for the second order fluctuation field, modulo an error term that we define as follows then we have that has to be estimated. Analogously to the previous computation we have where in the last step we used the symmetry of p(·). Then Hence we have It remains to show that the L 2 (P n ) norm of E (n,k) (ϕ, η(n 2 t)) vanishes in the limit as n goes to infinity. This is done in the following lemma: LEMMA 6.1. Let E (n,k) (ϕ, η) be given by (69), then, for every test function ϕ ∈Ŝ(R d ) there exists C > 0 such that, for all t ≥ 0 and n ∈ N, PROOF. Using the fact that ϕ is bounded and that p(·) has finite range we can conclude that there exists an M > 0 such that We recall here that the duality function is parametrized by the density parameter ρ, i.e. D(·, ·) = D ρ (·, ·) and that {D ρ (ξ, ·), ξ ∈ Ω} is a family of products of polynomials that are orthogonal with respect to the reversible measure ν ρ . From the stationarity of ν ρ we have The fact that we can exchange expectations and integral is a consequence of Proposition 6.1 in Section 6.2.2, which does not use any results of the current section.
Let us denote by V n (ϕ) the integrand in (76), then, using (15), we have where we used (15) in the second identity, (40) and (17) in the third identity (with c = c(k, ρ)) and (75) in the fourth line. From (73) we have Using (64) we have that the first term in the r.h.s. of (78) is bounded by a constant times n −1 .
For what concerns the second term, we have: Now, from the Taylor expansion (65) we know that there exists a sequence of functions where, using the fact that the range of p(·) is R = [−R, R] d , and the Taylor expansion (65) we have that there exists a smooth function ψ ∈ S(R d ) such that, for all x ∈ Z d , as a consequence we obtain the upper bound 1 n 2d where the inequality holds for a suitable c ′ > 0. In conclusion we have that there exists a constant C > 0 such that from which the statement follows.
As a consequence of Lemma 6.1 we can close the drift term, i.e.

Closing the equation for the carré-du-champ
In this section we will show that the integrand in the RHS of equation (58) can be expressed in terms of the (k −1)th-order fluctuation field Y (n,k−1) . To achieve this we consider the expression for the carré-du-champ given by (141) in the Appendix. For the case of our kth-order fluctuation field this becomes Notice that here we multiplied by a factor n d/2+1 the squared term in order to cancel the n 2 in front of the carré-du-champ and get a general factor n −d in front of the sum.
In the next section we find some recursion relations for duality polynomials. The main application of these relations consists in allowing us to rewrite any polynomial depending on η x,x+r in terms of polynomials depending on the unmodified η.

Recursion relation for duality polynomials
In this section we obtain a recurrence relation for the single-site orthogonal polynomials. Before giving the result it is convenient to summarize the expression for the self-duality generating function by defining the function then f σ can be written in the form with c σ given by (27). Then we define the functions g σ ,g σ : N → R given by for m ≥ 1 and g σ (0) =g σ (0) := 1 (86) that are exactly computable: for m ≥ 1, that can be rewritten as notice, in particular, thatg We have the following result. PROOF. From (85) we have that then, from the generating function definition (84), we deduce that hence, the recurrence relation (91) and an application of Leibniz product rule for differentiation in the RHS above give where in the second equality we used (92). This concludes the proof of (89). Equation (90) can be proved from the same reasoning, with the difference that we now have the inverse relation This change results, after the application of Leibniz rule, in the relation that concludes the proof.

Controlling the moments of the fields
The objective of this section is to take advantage of the ergodic properties of our process to introduce a result that will allow us to make multiple replacements, in the appropriate sense, inside the expression of the carré-du-champ given in (83). Let us start first with a uniform estimate for moments of the fields Y (n,l) (Φ, η).
PROOF. As claimed in the statement of the proposition, this result holds for any finite natural number m. Nevertheless for simplicity we will only show how to obtain the estimates for m ∈ {2, 4} (which indeed are the only two uses that we make of this result). Let us start with the simplest non-trivial case, m = 2, for which the result comes directly from orthogonality where in the second line we used (13) and K is given by Notice that the previous estimate was possible due the fact that orthogonality, in the form of expression (13), allowed us to reduce the summation in the RHS of (95) from a 2ld dimensional sum to an ld dimensional sum in (96).
For the case m = 4 we have For this case the sum in the RHS of (98) is 4ld-dimensional. Given the factor n −2ld in front of the RHS, in order to obtain a uniform estimate, we would like this summation to be 2ld dimensional instead. In order to see that this is indeed the case, we analyze the non-zero contribution coming from E νρ D(ξ (1) , η)D(ξ (2) , η)D(ξ (3) , η)D(ξ (4) , η) .
By the product nature of the measure ν ρ and the duality polynomials we have (99) Notice that for every x for which ξ (j) x = 0 for all j ∈ {1, 2, 3, 4}, the corresponding contribution in the RHS of (99) is equal to 1 and therefore negligible. This is precisely the reason why the summation in the RHS of (98) is at most 4ld-dimensional. We have indeed that the maximum number of x ∈ Z d contributing to the product in the RHS of (99) is at most 4l, i.e. one for each of the 4l particles that all the ξ (j) have in total. In reality we can see that there are less xs giving a non-zero contribution. In order to see is, consider an x ∈ Z d such that there exists a unique j ∈ {1, 2, 3, 4} for which ξ (j) x = 0. In this case, because of the zero mean of the single-site duality function we have E νρ d(ξ (1) x , η)d(ξ (2) x , η)d(ξ (3) x , η)d(ξ (4) x , η) = 0 (100) this means that whenever x ∈ Z d is such that there exists a j ∈ {1, 2, 3, 4} for which ξ (j) x = 0 there must be another j ′ ∈ {1, 2, 3, 4} for which ξ (j ′ ) x = 0. In other words we only have a possibility of 2l particles to distribute freely, and hence the summation in the RHS of (98) is at most 2ld-dimensional. PROPOSITION 6.2. Let f : R d → R be a test function, and {M n : Ω × R → R : n ∈ N} be a sequence of uniformly bounded cylindrical functions of the form where only a finite number of b j are different from zero. Let also {a n : n ∈ N} be a sequence of real numbers converging to 0, we then have for all l ∈ {1, 2, . . . , k − 1}, and m ∈ N.
PROOF. By Cauchy-Schwartz we have where in the last line we used Proposition 6.1, the boundedness of the single-site duality polynomials d(b j , η x ) and the smoothness of f in the representation (101). The result then follows from the convergence a n → 0.

The gradient of the fluctuation fields
Our goal for this section is to rewrite the square inside the RHS of (83) in terms of lower order fluctuation fields. We will see that this can be expressed, in agreement with (56) , only in terms of the field of order k − 1. Let us then denote by ∇ i,i+r d the d-dimensional gradient Notice that, by linearity of the k-th order field, we have with D(·, ·) as in (34). We define now, for i, j ∈ Z d , ℓ ≤ k, the auxiliary field then we have the following formula for the gradient of the fluctuation field.
Using the product nature of the polynomials D(·, η) and of Φ n (·) we get and now, calling b = κ and m = a we get where the first identity follows from the change of variable ℓ = κ + a. Then then This concludes the proof.
The advantage that Proposition 6.3 gives us is that we now have an expression in terms of the auxiliary field (105): Recall that we claimed that we are able to close the carré-du-champ in an expression depending only on the field of order k − 1. In order to achieve this it remains to: 1. replace the first term in the RHS of (108) by some expressions depending on the field of order k − 1; 2. show that the second term in the RHS of (108) vanishes as n → ∞.
Thanks to Proposition 6.2 we conclude the proof.
For what concerns the second step, let us denote by G (n,k) i,j (Φ, η) the second term in the RHS of (108), i.e. x,x+r (ϕ, η(n 2 s)) 2 , the statement follows from applying multiple times Propositions 6.4 and 6.2.
PROOF. We start with the proof (129) which is more involved. Thanks to stationarity, the expectation does not depend on time, and then we can ignore the supremum over time in (129). From (82) we already have an expression for the integrand of (129): n 2 L X (n,k) (ϕ (k) , η) = αk · χ 2 · X (n,k) (ϕ (k−1) ⊗ ∆ϕ, η) + O(n −1 ) recall that here again we are using the fact that the field X (n,k) can be also thought as acting on general ( not necessarily symmetric ) test functions. Because of stationarity it is enough to estimate then the desired bound is obtained by applying Proposition 6.1. In the same spirit we can use Proposition 6.1 to bound (128).