Factorized Duality, Stationary Product Measures and Generating Functions

We find all self-duality functions of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} D(\xi , \eta )= \prod _{x} d(\xi _x, \eta _x) \end{aligned}$$\end{document}D(ξ,η)=∏xd(ξx,ηx)for a class of interacting particle systems. We call these duality functions of simple factorized form. The functions we recover are self-duality functions for interacting particle systems such as zero-range processes, symmetric inclusion and exclusion processes, as well as duality and self-duality functions for their continuous counterparts. The approach is based on, firstly, a general relation between factorized duality functions and stationary product measures and, secondly, an intertwining relation provided by generating functions. For the interacting particle systems, these self-duality and duality functions turn out to be generalizations of those previously obtained in Giardinà et al. (J Stat Phys 135:25–55, 2009) and, more recently, in Franceschini and Giardinà (Preprint, arXiv:1701.09115, 2016) . Thus, we discover that only these two families of dualities cover all possible cases. Moreover, the same method discloses all simple factorized self-duality functions for interacting diffusion systems such as the Brownian energy process, where both the process and its dual are in continuous variables.


Introduction
Duality and self-duality are very useful and powerful tools that allow to analyze properties of a complicated system in terms of a simpler one. In case of self-duality for particle systems, the dual system is the same and the simplification arises because in the dual one considers only a finite number of particles (e.g. [3]).
Several methods are available to construct dual processes and duality relations. In the context of population dynamics, the starting point to find dualities is to consider the coalescent, the simplest example here being the duality between Kingman's coalescent block-counting process and the Wright-Fisher diffusion (for an overview of this kind of dualities in more general contexts, see [2]). Another method is provided by the pathwise dualities based on graphical constructions and time reversal, see e.g. [15,19].
In the context of conservative particle systems such as the exclusion process and its generalizations, zero range processes, etc. (cf. [3,15]), the algebraic method developed in [9] offers a general framework to construct self-duality functions starting from a reversible product measure by using symmetries of the generator, i.e. operators commuting with the generator. Additionally, if these symmetries are in product form, i.e. of the form x S x , where the action of S x depends only on the variables associated to site x, then the self-duality functions produced by these product-like symmetries and a reversible product measure also factorize over the sites, i.e. are a product over the sites of functions that depend only on the variables associated to that site. In this paper, we call such duality functions "simple factorized self-duality functions".
A complete picture of how to obtain all simple factorized self-duality functions for such particle systems is missing (except in the simplest case of symmetric simple exclusion, see e.g. [18] and references therein). Natural associated questions are: which of these conservative particle systems allow self-duality and is it possible then to obtain all simple factorized selfduality functions for these systems? One of the useful applications of disposing of all simple factorized self-duality functions is that, depending on the target, one can choose appropriate ones: e.g. in the hydrodynamic limit and the study of the structure of the stationary measures, the "classical" duality functions are the appropriate ones (see e.g. [3]), whereas in the study of (stationary and non-stationary) fluctuation fields and associated Boltzmann-Gibbs principles [12,Chap. 11], as well as in the study of speed of relaxation to equilibrium in L 2 or in the study of perturbation theory around models with duality (cf. [3]), "orthogonal" duality functions turn out to be very useful.
In this paper, we develop an approach to answer the above questions and systematically determine all simple factorized self-duality functions for a class of conservative interacting particle systems. As a consequence, by considering many-particle limits, we also obtain simple factorized dualities and self-dualities for a class of conservative diffusion processes, such as the Brownian energy process.
In this route, starting from examples, we first investigate a general connection between stationary product measures and simple factorized duality functions. This shows, in particular, that for infinite systems with simple factorized self-duality functions, the only stationary measures which are ergodic (w.r.t. either space-translation or time) are in fact product measures. Then we use this connection between stationary product measures and simple factorized duality functions to recover all possible candidate simple factorized duality functions from the stationary product measures. More precisely, we show that, given the first duality function, i.e. the duality function with a single dual particle, all other simple factorized duality functions are determined. This provides a simple machinery to obtain all simple factorized self-duality functions in processes such as Symmetric Exclusion Process (SEP), Symmetric Inclusion Process (SIP) and Independent Random Walkers (IRW). In particular, we recover via this method all orthogonal polynomial duality functions obtained in [8].
Moreover, we prove that in the context of conservative particle systems where the rates for particle hopping depend only on the number of particles of the departure and arrival sites, the processes SEP, SIP and IRW are the only systems which have self-duality with simple factorized self-duality functions and that the first duality function is necessarily an affine function of the number of particles. Next, in order to prove that the only candidate simple factorized self-duality functions derived via the method described above are actual self-duality functions, we develop a method based on generating functions. This method, via an intertwining relation, allows to go from discrete systems (particle hopping dynamics) to continuous systems (such as diffusion processes or deterministic dynamics) and back, and also allows to pass from self-duality to duality and back. The proof of a self-duality relation in a discrete system then reduces to the same property in a continuous system, which is much easier to check directly. The generating function method also provides new examples of self-duality for processes in the continuum such as the Brownian Energy Process (BEP), which intertwines with the SIP via the generating function. In fact, we show equivalence between self-duality of SIP, duality between SIP and BEP and self-duality of BEP. Finally, this method based on generating functions generalizes the concept of obtaining dualities from symmetries to intertwinings, being a symmetry an intertwining of the generator with itself.
Notice that we restrict in this paper to conservative particle systems (and associated diffusion processes) in discrete space, although there are many other possible contexts such as spin-systems, reaction diffusion systems, Fleming-Viot processes, etc. as well as continuumspace limits of such systems where similar duality questions could be addressed. Moreover, most of the emphasis will be on self-duality, as dualites in this context in fact follow from self-dualities.
Organization of the Rest of the Paper In Sect. 2 we introduce the basic definitions of duality and systems considered. Additionally, in Theorem 2.1 we prove which particle systems out of those considered admit simple factorized self-duality. In Sect. 3, we investigate a general relation between simple factorized duality functions and stationary product measures. We treat separately the finite and infinite contexts in which this relation arises; in the latter case, we exploit this connection to draw some conclusions on the product structure of ergodic measures. Section 4 is devoted to the derivation of all possible simple factorized self-duality and duality functions. Here Theorem 2.1 and the relation in the previous section are the two key ingredients. In Sect. 5, after an introductory example and a brief introduction on the general connection between duality and intertwining relations, we establish an intertwining between the discrete and the continuum processes. This intertwining relation is then used to produce all the self-duality functions for the Brownian Energy Process.

Setting
We start defining what we mean by duality for Markov processes. Then, we introduce a general class of Markov interacting particle systems with associated interacting diffusion systems arising as many-particle limits.
Main Result of the Section Theorem 2.1 states that the only conservative particle systems described by the infinitesimal generator (7) below which admit a "non-trivial" factorized selfduality are necessarily IRW, SIP and SEP-type of processes. Moreover, in the same statement, we find the general form of the first single-site self-duality function for such systems.

Duality
Given two (Polish) state spaces andˆ and two Markov processes {ξ(t), t ≥ 0} and {η(t), t ≥ 0} evolving on them, we say that they are dual with duality function D :ˆ × → R (where D is a measurable function) if, for all t > 0, ξ ∈ˆ and η ∈ , we have the so-called duality relationÊ If the laws of the two processes coincide, we speak about self-duality. More generally, we say that two semigroups {S(t), t ≥ 0} and {Ŝ(t), t ≥ 0} are dual with duality function D if, for all t ≥ 0, where "left" (resp. "right") refers to action on the left (resp. right) variable. In the case that these semigroups are Markov semigroups, (1) is exactly the same as (2). Even more generally, we say that two operators L andL are dual to each other with duality function D if In the context of Markov processes, the operators L andL which we have in mind here are the generators of the Markov processes. Moreover, we refer to [11] for more technical aspects of duality relations, e.g. when generator duality implies semigroup duality or which are the exact restrictions on the state spaces needed.

Lattice and Factorization over Sites
The underlying geometry of all systems that we will look at consists of a set of sites V either finite or V = Z d . Moreover we are given a family of transition rates p : V × V → R + , satisfying the following conditions: for all x, y ∈ V , (1) Vanishing diagonal p(x, x) = 0, (2) Irreducibility there exist x 1 = x, x 2 , . . . , x m = y such that m−1 l=1 p(x l , x l+1 ) > 0. In case of infinite V , we further require the following: Uniform bound on total jump rate sup x∈V y∈V p(x, y) < ∞.
Notice that when p is finite-range and translation invariant, then the uniform bound on total jump rate follows automatically. Notice also that the finite-range assumption is not necessary and can be relaxed; however, we assume it here for simplicity in order to avoid existence problems of the processes in the context when V = Z d .
To each site x ∈ V we associate a variable η x ∈ E = N, {0, . . . , N } or R + , with the interpretation of either the number of particles or the amount of energy associated to the site x. Configurations are denoted by η ∈ = E V .
As discussed in the introduction, we look at duality functions which factorize over sites, i.e. of the form We then call the functions d(ξ x , η x ) the singe-site duality functions and further assume d(0, ·) ≡ 1.
The above condition (5) is related to the fact that we want to have duality functions which make sense for infinite systems when the dual configuration has a finite total mass. A typical example is when η, ξ ∈ N Z d , where η is an infinite configuration while ξ is a finite configuration, so that in the product (4) there are only a finite number of factors different from d(0, η x ). In this sense, the choice d(0, ·) ≡ 1 is the only sensible one for infinite systems. When V is finite and E = N or {0, . . . , N }, this condition is not necessary and e.g. if a (positive) reversible product measure μ = ⊗ x∈V ν exists, then the so-called cheap self-duality function

Interacting Particle Systems with Simple Factorized Self-duality
The class of interacting particle systems we consider is described by the (formal) infinitesimal generator acting on local functions f : → R as follows: where L x,y , the single-edge generator, is defined as and η x,y denotes the configuration arising from η by removing one particle at x and putting it at y, i.e. η x,y Note the conservative nature of the system and the form of the particle jump rates in (7) for all λ > 0 for which the normalizing constant Z λ < ∞ and with ϕ(

Remark 2.1
The existence of the processes with formal generator (6) poses no problem if V is finite. If V is infinite, further growth conditions on the functions u(n) and v(n) are required in order to ensure non-explosion. For the processes that we will be considering in the next sections, which have the self-duality property, existence can be proved via self-duality as in [3, Sect. 2.2.4].

Basic Examples and Characterization of Self-dual Particle Systems
We recall here the basic examples of self-dual interacting particle systems and corresponding simple factorized self-duality functions known in literature (cf. e.g. [9]). Next, in Theorem 2.1 below we prove that these are the only particle systems (within our defined class) that are self-dual with simple factorized self-duality functions.
(I) Independent random walkers (IRW) In the following theorem we show that the only processes with generator of the type (6) which have non-trivial simple factorized self-duality functions are of one of the types described in the examples above, i.e. IRW, SIP or SEP. Here by "non-trivial" we mean that the first single-site self-duality function d(1, n) is not a constant (as a function of n). (7) is self-dual with simple factorized self-duality function D(ξ, η) = d(ξ x , η x ) in the form (4) with d(0, ·) ≡ 1 as in (5). If d (1, n) is not constant as a function of n, then

Theorem 2.1 Assume that the process with generator
and the first single-site self-duality function is of the form for some a ∈ R and b = 0.
Proof Using the self-duality relation for ξ x = 1 and no particles elsewhere, together with u(0) = 0, we obtain the identity Setting η x = η y = n ≥ 1, this yields anytime u(n)v(n) = 0 from which we derive d(1, n) = a + bn. Because d(1, n) is not constant as a function of n, we must have b = 0. Inserting d(1, n) = a + bn in (11) we obtain from which, by setting η x = n and η y = 0 we obtain the first in (9), while via η x = n and η y = 1 we get the second condition.

Remark 2.2
More generally, if we replace (5) with d(0, n) = 0 in the above statement, we analogously obtain (9) and for some constants a, b, c ∈ R, b, c = 0.

Interacting Diffusion Systems as Scaling Limits
Conservative interacting diffusion processes arise as scaling limits of the particle systems in Sect. 2.4 (cf. [9]). More in details, by "scaling limit" we refer to the limit process of the par- In case of IRW, one obtains a deterministic (hence degenerate diffusion) process {z(t), t ≥ 0} whose evolution is described by a first-order differential operator. In case of SIP(α), the scaling limit is a proper Markov process of interacting diffusions known as Brownian Energy Process (BEP(α)) (cf. [9]). For the SEP(γ ), this limit cannot be taken in the sense of Markov processes, but we can extend the SEP generator to a larger class of functions defined on a larger configuration space and take the many-particle limit. The limiting second-order differential operator is then not a Markov generator, but still a second order differential operator. We will explain this more in detail below.
The limiting differential operators in the case of IRW and SIP(α) can be described as acting on smooth functions f : E V → R as follows: with single-edge generators L x,y given, respectively, by and For the SEP(γ ) we proceed as follows. For each N ∈ N, consider the operator L N working on functions f : where This operator is not a Markov generator anymore, because the factors η x (γ −η y ) can become negative. With this operator, we consider the limit where η N = ( N z x ) x∈V and f : E V → R is a smooth function. This then gives the differential operator L which is the analogue of (16) in the context of SEP(γ ). This differential operator L , with single-edge operators does not generate a Markov process but it is still useful because, as we will see in Sect. 5 below, via generating functions, it is intertwined with the operator (19) for the choice N = 1.

Remark 2.3
The existence and ergodic properties of diffusion processes with generator of type (16) in the context of infinite volume V = Z d has been addressed in [10, Sect. 3, Proposition 3.1, Lemma 3.2, Example 3.4) (this is the proof of the existence of the so-called BMP process in infinite volume, the BEP process is then obtained from a transformation of the BMP process [9, Theorem 6.1]. In the finite-volume case, the process is a multi-type Wright-Fisher diffusion with mutation, of which the existence is well-known. Naturally, as we can see for the case of SIP(α) and BEP(α), when going to the scaling limit, some properties concerning stationary measures and duality pass to the limit. Indeed, BEP(α) admits a one-parameter family of stationary product measures {⊗ ν λ , λ > 0}, where and is dual to SIP(α) with simple factorized duality function After noting that property (5) holds also in this situation, we show that the first single-site duality functions d (1, x) between SIP(α) and BEP(α) are affine functions of z ∈ R + , as we found earlier for single-site self-dualities in Theorem 2.1.

Proposition 2.1 Assume that SIP(α) and BEP(α)'s single-edge generators are dual with simple factorized duality function D(ξ, z)
for some a, b ∈ R.
Proof The duality relation for ξ x = 1 and no particles elsewhere, by using (5), reads If we set z x = z y = z, then z 2 d 2 dz 2 d(1, z) = 0 leads to (25) as unique solution.

Relation Between Simple Factorized Duality Functions and Stationary Product Measures
In the examples of duality that we have encountered in the previous section, we have a universal relation between the stationary product measures and the simple factorized duality (4)-(5) and stationary product measures {μ λ = ⊗ ν λ , λ > 0}, then there is a relation between these measures and these functions, namely there exists a function θ (λ) such that This function θ (λ) is then simply the expectation of the first single-site duality function, i.e.
In the examples of Sect. 2.4, we have θ (λ) = λ for IRW, θ (λ) = λ 1−λ for SIP(α) and θ (λ) = λ 1+λ for SEP(γ ). In this section we first investigate under which general conditions this relation holds, and further use it in Sects. 3.2.1-3.2.2 as a criterion of characterization of all extremal measures. We refer to Sects. 2.1-2.2 for the general setting in which these results hold.
Later on, we will see that this relation (27) is actually a characterizing property of the simple factorized duality functions, meaning that all duality functions are determined once the first single-site duality function is fixed.

Finite Case
We start with the simplest situation in which V is a finite set.
First, we assume that the total number of particles/the total energy of the dual process is the only conserved quantity. More precisely, we assume the following property, which we refer to as harmonic triviality of the dual system: then H (ξ ) is only a function of |ξ | := x∈V ξ x .
Moreover, let us assume the simple factorized form of the duality function D(ξ, η) as in (4)- (5). Of the single-site functions d we may require the following additional property: MD The function d is measure determining, i.e. for two probability measures ν * , ν on E such that for all x ∈ V and ξ x ∈Ê it follows that ν * = ν .
Then we have the following. (1) with simple factorized duality function (4) satisfying condition (5). Moreover, assume that [HT] holds and that μ is a probability measure on . We distinguish two cases: (1) Interacting particle system case. IfÊ is a subset of N, then we assume that the self-duality functions D(ξ, ·) are μ-integrable for all ξ ∈ˆ . (2) Interacting diffusion case. IfÊ = R + , then we assume the following integrability condition: for each ε > 0, there exists a μ-integrable function f ε such that Then (a) μ is a stationary product measure for the process where δ x denotes the configuration with a single particle at x ∈ V and no particles elsewhere.
By μ-integrability in the interacting particle system case (resp. (31) in the interacting diffusion case), self-duality and invariance of μ we have Therefore by [HT] we conclude that H (ξ ) = ψ(|ξ |). By using d(0, ·) ≡ 1 and the factorization of the duality functions, we have that ψ(0) = 1. For the particle case, we obtain that In particular, we obtain that the l.h.s. does not depend on x. Next, for n ≥ 2, put then we have for x = y ∈ V , using the simple factorized duality function and the product form of the measure, from which it follows that ψ(n) = ψ(1) n . Via an analogous reasoning that uses the factorization of D(ξ, η) and the product form of μ, for the diffusion case we obtain, for all ε, ρ ≥ 0, and hence, by measurability of ψ(ε), we get ψ(ε) = ψ(1) ε .
To prove the other implication, put We then have by assumption and so it follows that μ is stationary by self-duality, μ-integrability, the conservation of the number of particles and the measure-determining property. Indeed, From the factorized form of D(ξ, η), (32) implies that for all x ∈ V and ξ x ∈Ê and also therefore μ is a product measure by the fact that d is measure determining.

Infinite Case
If V = Z d , then one needs essentially two extra conditions to state an analogous result in which a general relation between duality functions and corresponding stationary measures can be derived. In this section we will assume that the dual process is a discrete particle system, i.e.Ê is a subset of N, in which the number of particles is conserved. In this case we need an additional property ensuring that for the dynamics of a finite number of particles there are no bounded harmonic functions other than those depending on the total number of particles. Therefore, we introduce the condition of existence of a successful coupling for the discrete dual process with a finite number of particles. This is defined below. Definition 3. 1 We say that the discrete dual process {ξ(t), t ≥ 0} has the successful coupling property when the following holds: if we start with n particles then there exists a labeling such that for the corresponding labeled process {X 1 (t), . . . , X n (t), t ≥ 0} there exists a successful coupling. This means that for every two initial positions x = (x 1 , . . . , x n ) and y = (y 1 , . . . , y n ), there exists a coupling with path space measure P x,y such that the coupling time is finite P x,y almost surely.
Notice that the successful coupling property is the most common way to prove the following equivalent property (cf. [15]), which is the analogue of [HT], referred here to as bounded harmonic triviality of the dual process: [BHT] If H is a bounded harmonic function, then H (ξ ) = ψ(|ξ |) for some bounded ψ : E → R.

Remark 3.1
The condition of the existence of a successful coupling (and the consequent bounded harmonic triviality) is quite natural in the context of interacting particle systems, where we have that a finite number of walkers behave as independent walkers, except when they are close and interact. Therefore, the successful coupling needed is a variation of the Ornstein coupling of independent walkers, see e.g. [3,5,14].
Furthermore, we need a form of uniform μ-integrability of the duality functions which we introduce below and call uniform domination property of D w.r.t. μ (note the analogy with condition (31)): [UD] Given μ a probability measure on , the duality functions {D(ξ, ·), |ξ | = n} are uniformly μ-integrable, i.e. for all n ∈ N there exists a function f n such that f n is μ-integrable and such that for all η ∈ sup ξ ∈ˆ , |ξ |=n Under these conditions, the following result holds, whose proof resembles that of Theorem 3.1. (1) that {η(t), t ≥ 0} is dual to the discrete process {ξ(t), t ≥ 0} with simple factorized duality function as in (4)- (5). Moreover, assume [BHT] in place of [HT] for the dual process and that μ is a probability measure on such that [UD] holds. Then the same conclusions as in Theorem 3.1 follow, where (32) holds for all finite configurations ξ ∈ˆ .

Translation Invariant Case
In this section, we show that under the assumption of simple factorized duality functions, minimal ergodicity assumptions on a stationary probability measure μ on are needed to ensure (32) and, as a consequence, that μ is product measure.
Here we restrict to the case V = Z d because we will use spatial ergodicity.  η). By the Birkhoff ergodic theorem, we have that μ-a.s. as N → ∞. Using this, together with (44), we have Iterating this argument gives (32).

Remark 3.2
As follows clearly from the proof, the condition of factorization of the duality function can be replaced by the weaker condition of for all η, x 1 , . . . , x n . We note that this approximate factorization of the duality function leads to (32), though μ is not necessarily a product measure.

Non-translation Invariant Case
We continue here with V = Z d but drop the assumption of translation invariance. Indeed, equality (32) is also valid in contexts where one cannot rely on translation invariance. Examples include spatially inhomogeneous SIP(α) and SEP(γ ), where the parameters α = (α x ) x∈V and γ = (γ x ) x∈V in Sect. 2.3 may depend on the site accordingly (cf. e.g. [17]). Also in this inhomogeneous setting the self-duality functions factorize over sites and the stationary measures are in product form, with site-dependent single-site duality functions, resp. site-dependent marginals. We will show that the relation (32) between the self-duality functions and any ergodic stationary measure still holds, and as a consequence this ergodic stationary measures is in fact a product measure. The idea is that the averaging over space w.r.t. μ, used in the proof of Theorem 3.3 above, can be replaced by a time average. If we start with a single dual particle, the dual process is a continuous-time random walk on V , for which we denote by p(t; x, y) the transition probability to go from x to y in time t > 0. A basic assumption will then be lim t→∞ p(t; x, y) = 0 for all x, y ∈ V .

Theorem 3.4
In the setting of Theorem 3.2 with D(ξ, η) duality function and μ probability measure on , if μ is an ergodic stationary measure for the process {η(t), t ≥ 0} and (48) holds for the dual particle, then we have (32) for all finite configurations ξ ; as a consequence, μ is a product measure.

Proof
The idea is to replace the spatial average in the proof of Theorem 3.3 by a Cesaro average over time, which we can deal by combining assumption (48) with the assumed temporal ergodicity. Fix x 1 , . . . , x n ∈ V , y ∈ V . Define It is sufficient to obtain that H n+1 (x 1 , . . . , x n , y) = H 1 (y)H n (x 1 , . . . , x n ). We already know by the bounded harmonic triviality that H n only depends on n and not on the given locations x 1 , . . . , x n . Therefore, we have By assumption (48), this implies where in the last step we used the assumed temporal ergodicity of μ and Birkhoff ergodic theorem.

From Stationary Product Measures to Duality Functions
As we have just illustrated in Sects As a consequence, knowing the first single-site duality function d(1, ·) and the explicit expression of the marginal ν λ is enough to recover the l.h.s. in (54). However, rather than obtaining d(k, ·), at this stage the l.h.s. has still the form of an "integral transform"-type of expression for d(k, ·).
In the next two subsections, we show how to recover d(k, η x ) from (54) and the knowledge of θ (λ). This then leads to the characterization of all possible simple factorized (self-)duality functions. (27), together with the knowledge of the first singlesite self-duality function, determines all candidate simple factorized (self)-duality functions. This is shown in Section 4.1 for particle systems (self-duality for Exclusion, Inclusion and Independent Random Walkers) and in Sect. 4.2 for diffusion processes (duality between Inclusion and Brownian Energy processes).

Particle Systems and Orthogonal Polynomial Self-duality Functions
Going back to the interacting particle systems introduced in Sect. 2.3 with infinitesimal generator (6) and product stationary measures with marginals (8), the integral relation (54) rewrites, for each k ∈ N and λ > 0 for which Z λ < ∞, as where as the Taylor series expansion around λ = 0 of the function θ (λ) k Z λ , we can re-obtain the explicit formula of d(k, n)ϕ(n) as its nth order derivative evaluated at λ = 0, namely and hence, anytime ϕ(n) > 0, Together with the full characterization obtained in Theorem 2.1 of the first single-site selfduality functions d(1, ·) -and θ (λ) in turn -we obtain via this procedure a full characterization of all single-site self-duality functions. Beside recovering the "classical" dualities illustrated in Sect. 2.3, we find the single-site self-duality functions in terms of orthogonal polynomials { p k (n), k ∈ N} of a discrete variable (cf. e.g. [16]) recently discovered via a different approach in [8]. We add the observation that all these new single-site self-duality functions can be obtained from the classical ones via a Gram-Schmidt orthogonalization procedure w.r.t. the correct probability measures on N, namely the marginals of the associated stationary product measures (cf. [8]).
We divide the discussion in three cases, one suitable for processes of IRW-type, the other for SIP and SEP and the last one for the remaining particle systems.

Independent Random Walkers
We recall that the IRW-case corresponds to the choice of values in (9) satisfying the relation v(1) = v(0). If we compute θ (λ) for the general first single-site self-duality function d(1, n) = a + bn obtained in (10), we get and, in turn via relation (58), we recover all functions d(k, ·), for k > 1: In case a = 0, d(k, n) = 0 for n < k, while for n ≥ k, in the summation all terms but the one corresponding to r = k vanish, thus d(k, n) = n! (n−k)! b k . In case a = 0, In conclusion, for the choice a · b < 0, where {C k (n; μ), k ∈ N} are the Poisson-Charlier polynomials -orthogonal polynomials w.r.t. the Poisson distribution of parameter μ > 0 (cf. [16]).

Inclusion and Exclusion Processes
For SIP and SEP we are in the case v(1) = v(0), and hence we abbreviate where for SIP(α) we choose σ = 1 and β = α, while for SEP(γ ) we set σ = −1 and β = γ . If we compute θ (λ) for d(1, n) = a + bn in (10), we have By applying formula (58), we obtain all functions d(k, n) for k > 1 as follows: In case a = 0, clearly d(k, n) = 0 for k < n, while for n ≥ k only the term for r = k is nonzero in the summation: In case a = 0, by using the known relation (cf. [16], p. 51), If σ = 1, β = α > 0 and if a · b < 0, we recognize in (67) the Meixner polynomials as defined in [16], i.e.
where in our case μ = a a−bα and {M k (n; μ, α), k ∈ N} are the Meixner polynomialsorthogonal polynomials w.r.t. the discrete Gamma distribution of scale parameter μ and shape parameter α.
Furthermore, if σ = −1, β = γ with γ ∈ N and given the additional requirements a · b < 0 and − a b ≤ γ , from the expression in (67) we have a representation of the single-site duality functions in terms of the Kravchuk polynomials as defined in [16], i.e.
where p = − a bγ in our case and {K k (n; p, γ ), k ∈ N} the Kravchuk polynomials -orthogonal polynomials w.r.t. the Binomial distribution of parameters γ and p.
As a conclusion of this procedure, we note that all factorized self-duality functions for independent random walkers, inclusion and exclusion processes satisfying (5) are either in the "classical" form of Sect. 2.3 (case a = 0) or consist of products of rescaled versions of orthogonal polynomials (case a = 0). Other simple factorized self-duality functions for the systems IRW, SIP and SEP do not exist.

Remark 4.1
It is interesting to note that, apart from the leading factor a k , the remaining polynomials in the expressions of d(k, n) for a = 0 are "self-dual" in the sense of the orthogonal polynomials literature, i.e. p n (k) = p k (n) in our context (cf. Definition 3.1, [13]). Henceforth, if d(k, n) is interpreted as a countable matrix, the value a ∈ R is the only responsible for the asymmetry of d(k, n): upper-triangular for a = 0 while symmetric for a = 1.

Trivial Factorized Self-duality
To conclude, for the sake of completeness, we can implement the same machinery to cover all factorized self-dualities with property (5) for all discrete processes of type (6).
Indeed, from the proof of Theorem 2.1, if the process is neither of the types IRW, SIP and SEP, then the only possible choice is d (1, n) = a for some a ∈ R, i.e. it is not depending on n. From this we get θ (λ) = a, and d(k, n) = a k from formula (58). Hence, the self-duality functions must be of the form i.e. depending only on the total number of dual particles (and not on the configuration η). Hence, the duality relation in that case reduces to the trivial relation, for all t ≥ 0 and ξ ∈ˆ , which is just conservation of the number of particles in the dual process. No other self-duality relation with simple factorized self-duality functions can exist.

Interacting Diffusions and Orthogonal Polynomial Duality Functions
As shown in Theorem 3.1, relation (54) still holds whenever the discrete right-variables n ∈ N are replaced by continuous variables z ∈ R + and sums by integrals. With this observation in mind, we provide a second general method to characterize all simple factorized duality functions between the continuous process BEP(α) and its discrete dual SIP(α). More precisely, if d(k, z) is a single-site duality function with property (5) between BEP(α) and SIP(α), and ν λ is the stationary product measure marginal for BEP(α) as in (23), then, from the analogue of relation (54) for k = 1, namely we necessarily have by Theorem 3.1 that As a consequence, the function d(k, z) z α−1 (α) is the inverse Laplace transform of θ(λ) k λ −α . Given the first single-site duality function d (1, z) in (25), from (72) we obtain As a consequence, the r.h.s. in (73) becomes and there exist explicit expressions for the inverse Laplace transform of this function. We split the computation in two cases. In case a = 0, since the inverse Laplace transform of λ −α−k is z α+k−1 (α+k) , we immediately obtain i.e. the "classical" single-site duality function as in Sect. 2.5, up to set b = 1 α . In case a = 0, the inverse Laplace transform of (75) is more elaborated: As the above expression must equal d(k, z) z α−1 (α) , it follows that As a final consideration, we note that for the choice a · b < 0, where μ = − bα a here and {L k (z; α − 1, μ), k ∈ N} are the generalized Laguerre polynomials-orthogonal polynomials w.r.t. to the Gamma distribution of scale parameter 1 μ and shape parameter α as defined in [16].

Intertwining and Generating Functions
In this section, we introduce the generating function method (cf. [9]), which allows to go from a self-duality of a discrete process towards duality between a discrete and continuous process, and further towards a self-duality of a continuous process, and back. This then allows e.g. to simplify the proof of a discrete self-duality by lifting it to a continuous self-duality, which is usually easier to verify. The key of this method is an intertwining between discrete multiplication and derivation operators and their continuous analogues, via an appropriate generating function.
To reduce issues of well-definition of operators and their domains, in this section we restrict our discussion to the case of finite vertex set V .
We start with the introductory example of Independent Random Walkers in Sect. 5.1, showing its generator intertwines with a first order differential operator, and from that recover in an easy way the self-duality of independent random walkers with the "classical" self-duality functions.

Main Results of the Section
After an introduction to intertwining in Sect. 5.2 and its relation with duality in Theorem 5.1, we find in Propositions 5.1-5.2 (product) intertwiners between the discrete particle system generators and their diffusion counterparts. As a corollary of Propositions 5.1-5.2, in Sect. 5.4 we prove that all the candidate simple factorized selfduality functions produced in Sect. 4.1 (from the stationary product measures via (27)) are actual self-duality functions. Seemingly, we also produce several (self-)duality functions for the diffusion counterparts of the particle systems.

Introductory Example
To make the method clear, let us start with a simple example of independent random walkers on a single edge. The generator is with n 1 , n 2 ∈ N. Define now the (exponential) generating function Then it is easy to see that where Now assume that we have a self-duality function for the particle system, i.e.
then we have where In words, a self-duality function of L is "lifted" to a duality function between the independent random walk generator L and its continuous counterpart L by applying the generating function to the n-variables. Conversely, given a duality function between the independent random walk generator L and its continuous counterpart L , its Taylor coefficients provide a self-duality function of L.
We can then also take the generating function w.r.t. the k-variables in the function D to produce a self-duality function for L , i.e. defining we have For the classical self-duality function D(k 1 , k 2 ; n 1 , n 2 ) = n 1 ! (n 1 −k 1 )! n 2 ! (n 2 −k 2 )! 1{k 1 ≤ n 1 }1{k 2 ≤ n 2 }, we find that Beside the factor e z 1 +z 2 which depends only on the conserved quantity z 1 + z 2 , to check the self-duality relation for L w.r.t. the function e v 1 z 1 +v 2 z 2 is rather straightforward, the computation involving only derivatives of exponentials. By looking at the Taylor coefficients w.r.t. both v and z-variables of this self-duality relation, we obtain the self-duality relation for L w.r.t. D where we started from.
In conclusion, all these duality relations turn out to be equivalent, and the proof of selfduality for particle systems requiring rather intricate combinatorial arguments (cf. e.g. [3]) is superfluous once the more direct self-duality for diffusion systems is checked.

Intertwining and Duality
The notion of intertwining between stochastic processes was originally introduced by Yor in [20] in the context of Markov chains and later pursued in [4] and [6] as an abstract framework, in discrete-time and continuous-time respectively, for the problem of Markov functionals, i.e. finding sufficient and necessary conditions under which a random function of a Markov chain is again Markovian.
For later purposes, we adopt a rather general definition of intertwining, in which {η(t), t ≥ 0} and {ζ(t), t ≥ 0} are continuous-time stochastic processes on the Polish spaces and , respectively, whose expectations read E, E resp., and M( ) denotes the space of signed measures on . We say that {ζ(t), t ≥ 0} is intertwined on top of {η(t), t ≥ 0} if there exists a mapping Λ : → M( ) such that, for all t ≥ 0, ζ ∈ and f : → R, Working at the abstract level of semigroups, we say that {S (t), t ≥ 0} on a space of functions f : → R denoted by F ( ), is intertwined on top of {S(t), t ≥ 0}, a semigroup on a space of functions f : → R denoted by F ( ), with intertwiner Λ if Λ is a linear operator from F ( ) into F ( ) and if, for all t ≥ 0 and f : → R, Similarly, operators L with domain D(L ) and L with domain and Notice that with a slight abuse of notation we used the same symbol Λ for an abstract intertwining operator as for the intertwining mapping. In other words, in case the intertwining mapping as in (90) is given byΛ, then the corresponding operator is An intertwining mapping Λ has a probabilistic interpretation if it takes values in the subset of probability measures on . Indeed, in (90) the process {ζ(t), t ≥ 0} may be viewed as an added structure on top of {η(t), t ≥ 0} or, alternatively, the process {η(t), t ≥ 0} as a random functional of {ζ(t), t ≥ 0}, in which Λ provides this link.
where L † denotes the transpose of L, and Intertwiners as Λ in (95) may be also interpreted as natural generalizations of symmetries of generators, indeed (95) with L = L just means that commutes with L, which is the definition of a symmetry of L. As outlined in [9, Theorem 2.6], the knowledge of symmetries of a generator and dualities of this generator leads to the construction of new dualities. The following theorem presents the analogue procedure in presence of intertwiners: a duality and an intertwining lead to a new duality.
and a duality function D :ˆ × → R forL and L, namely D(ξ, ·) ∈ D(L) for all ξ ∈ˆ , D(·, η) ∈ D(L) for all η ∈ and L left D = L right D. (97) Then, if Λ right D(ξ, ·) ∈ D(L ) for all ξ ∈ˆ and Λ right D(·, ζ ) ∈ D(L) for all ζ ∈ , Λ right D is a duality function forL and L , i.e. ProofL Here in the first equality we used that left and right actions commute, in the second equality we used the assumed duality ofL and L, and in the third equality we used the assumed intertwining.

Intertwining Between Continuum and Discrete Processes
In this section we prove the existence of an intertwining relation between the interacting diffusion processes presented in Sect. 2.5 and the particle systems of Sect. 2.3. This intertwining relation provides a second connection, besides the scaling limit procedure (cf. Sect. 2.5), between continuum and discrete processes, which proves to be better suited for the goal of establishing duality relations among these processes. Indeed, the characterization of all possible simple factorized self-dualities for particle systems obtained in Sect. 4 and the intertwining relation below, via the application of Theorem 5.1, produces a characterization of all possible dualities, resp. self-dualities, between the discrete and the continuum processes, resp. of the continuum process. In the following proposition, we prove the intertwining relation for operators L σ,β and L σ,β defined, respectively, on functions f : where and, on real analytic functions f : where L σ,β x,y f (z) Note that L σ,β in (100) is a special instance of the generator L u,v in (6) with conditions (9), while L σ,β above, matches on a common sub-domain, for particular choices of the parameters σ and β, those in (16).
that intertwines L σ,β and L σ,β , in this order. The natural candidate for H is the "inverse operator" of G, whenever this is well-defined. In general, this "inverse intertwiner" lacks any probabilistic interpretation, but indeed establishes a second intertwining relation useful in the next section.

Proposition 5.2
Let H be the differential operator mapping real analytic functions g : R → R into functions H g : N → R as Then H is the inverse operator of G, namely, for all f : Moreover, the tensorized operator H ⊗ = ⊗ x∈V H (x) is an intertwiner for L σ,β and L σ,β , i.e. for all real analytic g : Before giving the proof, we need the following lemma.

Lemma 1 Let A be the operator acting on functions f
Then, the tensorized operator A ⊗ = ⊗ x∈V A (x) is a symmetry for the generator L σ,β , i.e. for all f : Proof Instead of going through tedious computations, we exploit the fact that the operator A ⊗ has the form where J (x) is an operator defined for functions f : N V → R which acts only on the xth variable as Since all these operators {J (x) , x ∈ V } commute over the sites we have We conclude the proof by noting that the operator is a symmetry for the generator L σ,β , cf. e.g. [1,9].
form of the self-duality functions. Indeed, while the single-site self-duality functions for the SIP(α) process, for instance, have the generic form of an hypergeometric function the single-site duality functions between discrete and continuum processes involve in their expressions hypergeometric functions and those for the self-duality of continuum processes are even simpler, namely as the number of arguments of the hypergeometric function drops. The tables below schematically report all single-site (self)-duality functions for the operators L σ,β in (100) and L σ,β in (102). Recall that the parameters a, b ∈ R in (9)   More in detail, on the left-most column we place the single-site self-duality functions d(k, n) for the particle systems of Sect. 2.3: while the top-left functions are those already appearing in [1,9], cf. also Sect. 2.3 and (66)-and, thus, for this reason denoted here as the "classical" ones -the second-to-the-top functions are those derived in Sect. 4.1 in (61)-(67) and being related to suitable families of orthogonal polynomials. While these two classes of single-site self-duality functions satisfy condition (5) (they are the only ones doing so), the bottom-left single-site self-duality functions correspond to the "cheap" self-duality (cf. end of Sect. 2.2), namely the detailed-balance condition w.r.t. the measures {⊗ x∈V ν λ , λ > 0} with marginals (8).
On the mid-column, we find the single-site duality functions between the difference operators L σ,β and the differential operators L σ,β , obtained from their left-neighbors by a direct application of the operator G in (104) on the n-variables. The new functions will depend hence on the two variables k ∈ N and z ∈ R + .
A second application w.r.t. the k-variables of the same operator G on the functions just obtained gives us back the right-most column, functions depending now on variables v, z ∈ R + . These functions represent the single-site self-duality functions for the differential operator L σ,β . As an immediate consequence of Proposition 5.2, we could also proceed from right to left by applying the inverse intertwiner H in (114).
Note that the single-site self-duality functions for L σ,β on the right-most columns, though they have been derived from different discrete analogues, i.e. classical, orthogonal and cheap single-site functions, within the same table they differ only for a factor which depends only on the conserved quantities |z| = x∈V z x and |v| = x∈V v x . Henceforth, when proving the self-duality relation, this extra-factor does not play any role and it is enough to check that the functions for constants c ∈ R, are single-site self-duality functions for the operators L 0,1 , L 1,α and L −1,γ , respectively. This final computation is the content of the next proposition.
Proof To prove that is a single-site self-duality function for the differential operator L 0,1 , we first observe that Hence, the self-duality relation for the single-edge generator L 1,0 x,y rewrites which indeed holds. For the second proof of self-duality for the single-site function we use the following shortcut: for x ∈ V , Additionally, we recall a formula for the z x -derivative of F x , namely and a recurrence identity Hence, the self-duality relation for L 1,α x,y w.r.t. the function F x (α)F y (α) rewrites by using (141) as cz x v y F x (α + 1)F y (α) + cv x z y F x (α)F y (α + 1) = cv x z y F x (α + 1)F y (α) + cv y z x F x (α)F y (α + 1) By substituting (142), the identity holds. The proof for the operator L −1,γ follows the same lines and we omit it.