# Using Malliavin calculus to solve a chemical diffusion master equation

Alberto Lanconelli Dipartimento di Scienze Statistiche Paolo Fortunati,
Università di Bologna, Bologna, Italy. e-mail<EMAIL_ADDRESS>
###### Abstract

We propose a novel method to solve a chemical diffusion master equation of
birth and death type. This is an infinite system of Fokker-Planck equations
where the different components are coupled by reaction dynamics similar in
form to a chemical master equation. This system was proposed in [3] for
modelling the probabilistic evolution of chemical reaction kinetics associated
with spatial diffusion of individual particles. Using some basic tools and
ideas from infinite dimensional Gaussian analysis we are able to reformulate
the aforementioned infinite system of Fokker-Planck equations as a single
evolution equation solved by a generalized stochastic process and written in
terms of Malliavin derivatives and differential second quantization operators.
Via this alternative representation we link certain finite dimensional
projections of the solution of the original problem to the solution of a
single partial differential equations of Ornstein-Uhlenbeck type containing as
many variables as the dimension of the aforementioned projection space.

Key words and phrases: Particle-based reaction-diffusion models, Fock space,
Malliavin calculus.

AMS 2000 classification: 60H07; 60H30; 92E20.

##  1 Introduction and statement of the main result

Suppose we want to model the probabilistic evolution of a system that is
initially constituted by a single particle of a chemical species $A$, which is
located somewhere in the interval $[0,1]$ according to a given probability
density function $\zeta:[0,1]\to\mathbb{R}$, and that undergoes:

  * •

_degradation_ and _creation_ chemical reactions

$\displaystyle\mbox{(I)}\quad
A\xrightarrow{\lambda_{d}(x)}\varnothing\quad\quad\mbox{(II)}\quad\varnothing\xrightarrow{\lambda_{c}(x)}A,$
(1.1)

where $\lambda_{d}(x)$ denotes the propensity for reaction (I) to occur for a
particle located at position $x\in[0,1]$ (i.e., the probability per unit of
time for this particle to disappear) while $\lambda_{c}(x)$ is the propensity
for a new particle to be created at position $x\in[0,1]$ by reaction (II);

  * •

drift-less isotropic _diffusion_ in space.

While reactions (1.1) alone can be analysed via the standard chemical master
equation [5],[8],[11] and the sole diffusive motion of the particles through a
Fokker-Planck equation [1],[4], the combination of these two phenomena makes
the mathematical description quite challenging. This is due to the hybrid
nature of the considered reaction-diffusion process, i.e. discrete in the
evolution of the number of particles (and hence of the spatial dimension of
the problem) and continuous in the random movement of those particles. In the
recent paper [3], the authors proposed a set of equations for the functions

$\displaystyle\rho_{n}(t,x_{1},...,x_{n}):=p_{n}(t,x_{1},...,x_{n})\mathbb{P}(N(t)=n),\quad
n\geq 0,$

that aims to model the reaction-diffusion process described above. Here
$p_{n}(t,x_{1},...,x_{n})$ represents the joint probability density function,
conditioned to the event $\\{N(t)=n\\}$, i.e. the number of particles at time
$t$ is equal to $n$, for the positions of these particles at time $t$. The
model is named _chemical diffusion master equation_ of birth and death type
and it takes the form:

$\displaystyle\begin{split}\partial_{t}\rho_{0}(t)=&\int_{0}^{1}\lambda_{d}(y)\rho_{1}(t,y)dy-\int_{0}^{1}\lambda_{c}(y)dy\cdot\rho_{0}(t),\quad
t>0;\\\
\partial_{t}\rho_{n}(t,x_{1},...,x_{n})=&\sum_{i=1}^{n}\partial^{2}_{x_{i}}\rho_{n}(t,x_{1},...,x_{n})\\\
&+(n+1)\int_{0}^{1}\lambda_{d}(y)\rho_{n+1}(t,x_{1},...,x_{n},y)dy\\\
&-\sum_{i=1}^{n}\lambda_{d}(x_{i})\rho_{n}(t,x_{1},...,x_{n})\\\
&+\frac{1}{n}\sum_{i=1}^{n}\lambda_{c}(x_{i})\rho_{n-1}(t,x_{1},...,x_{i-1},x_{i+1},...,x_{n})\\\
&-\int_{0}^{1}\lambda_{c}(y)dy\cdot\rho_{n}(t,x_{1},...,x_{n}),\quad n\geq
1,t>0,(x_{1},...,x_{n})\in]0,1[^{n},\end{split}$ (1.2)

with initial and Neumann boundary conditions

$\displaystyle\begin{split}\rho_{0}(0)&=0;\\\
\rho_{1}(0,x_{1})&=\zeta(x_{1}),\quad x_{1}\in[0,1];\\\
\rho_{n}(0,x_{1},...,x_{n})&=0,\quad n>1,(x_{1},...,x_{n})\in[0,1]^{n};\\\
\partial_{\nu}\rho_{n}(t,x_{1},...,x_{n})&=0,\quad n\geq 1,t\geq
0,(x_{1},...,x_{n})\in\partial[0,1]^{n}.\end{split}$ (1.3)

This is an infinite system of Fokker-Planck equations where the components
have an increasing number of degrees of freedom (to account for all the
possible numbers of particles in the system) and are coupled through the
reaction mechanism (1.1). The term

$\displaystyle\sum_{i=1}^{n}\partial^{2}_{x_{i}}\rho_{n}(t,x_{1},...,x_{n})$

in (1.2) refers to the drift-less isotropic spatial diffusion; the terms

$\displaystyle(n+1)\int_{0}^{1}\lambda_{d}(y)\rho_{n+1}(t,x_{1},...,x_{n},y)dy-\sum_{i=1}^{n}\lambda_{d}(x_{i})\rho_{n}(t,x_{1},...,x_{n})$

formalize gain and loss, respectively, due to reaction (I), while

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\lambda_{c}(x_{i})\rho_{n-1}(t,x_{1},...,x_{i-1},x_{i+1},...,x_{n})-\int_{0}^{1}\lambda_{c}(y)dy\cdot\rho_{n}(t,x_{1},...,x_{n})$

relate to reaction (II). The functions
$\lambda_{d},\lambda_{c},\zeta:[0,1]\to\mathbb{R}$ in (1.2)-(1.3) are assumed
to be non negative, bounded and measurable; in addition, $\zeta$ is smooth
with $\zeta^{\prime}(0)=\zeta^{\prime}(1)=0$ and $\int_{[0,1]}\zeta(x)dx=1$.
The symbol $\partial_{\nu}$ in (1.3) stands for the directional derivative
along the outer normal vector at the boundary of $[0,1]^{n}$. The particles
are assumed to be indistinguishable thus entailing the symmetry of
$p_{n}(t,x_{1},...,x_{n})$, and hence of $\rho_{n}(t,x_{1},...,x_{n})$, in all
their space variables. We observe that by construction,

$\displaystyle\int_{[0,1]^{n}}p_{n}(t,x_{1},...,x_{n})dx_{1}\cdot\cdot\cdot
dx_{n}=1,$

and hence

$\displaystyle\int_{[0,1]^{n}}\rho_{n}(t,x_{1},...,x_{n})dx_{1}\cdot\cdot\cdot
dx_{n}=\mathbb{P}(N(t)=n).$ (1.4)

Therefore, since $\sum_{n\geq 0}\mathbb{P}(N(t)=n)=1$, the solution to 1.2
should fulfil

$\displaystyle\sum_{n\geq
0}\int_{[0,1]^{n}}\rho_{n}(t,x_{1},...,x_{n})dx_{1}\cdot\cdot\cdot dx_{n}=1.$
(1.5)

But this is indeed the case; in fact, if we assume the functions $\lambda_{d}$
and $\lambda_{c}$ to be constant and we integrate out the spatial degree of
freedom in (1.2), we see that the diffusive part vanishes by virtue of Gauss
theorem combined with the Neumann boundary condition in (1.3); moreover,
recalling (1.4) we get

$\displaystyle\begin{split}\partial_{t}\mathbb{P}(N(t)=n)=&(n+1)\lambda_{d}\mathbb{P}(N(t)=n+1)-n\lambda_{d}\mathbb{P}(N(t)=n)\\\
&+\lambda_{c}\mathbb{P}(N(t)=n-1)-\lambda_{c}\mathbb{P}(N(t)=n).\end{split}$
(1.6)

Equation (1.6) is exactly the chemical master equation as derived from the
classical law of mass action for reactions (1.1) in spatially well-mixed
systems [5]; this fact proves condition (1.5) and the consistency of model
(1.2)-(1.3) with classical equations of chemical kinetics.

Aim of the present paper is to propose a novel method to analytically solve
system (1.2)-(1.3). Our idea consists in transforming the deterministic
problem (1.2)-(1.3) into an equivalent stochastic one by using iterated Itô
integrals with respect to a one dimensional Brownian motion and the Wiener-Itô
chaos expansion theorem. Then, exploiting the Gaussian framework induced by
the aforementioned Itô integration, we obtain a single partial differential
equation of Ornstein-Uhlenbeck type in infinitely many variables which is
directly connected to the problem (1.2)-(1.3) via the so-called Stroock-Taylor
formula: this provides an analytical representation for the solution of the
original system.
The main steps of our approach are summarized as follows:

  * •

we assume the existence of a classical solution $\\{\rho_{n}\\}_{n\geq 0}$ for
(1.2)-(1.3); the continuity of $\rho_{n}(t,\cdot)$ together with its symmetry
in the spatial variables imply the membership of $\rho_{n}(t,\cdot)$ to
$L_{s}^{2}([0,1]^{n})$, the space of symmetric square integrable functions;

  * •

we Itô-integrate all the space variables of $\rho_{n}(t,\cdot)$ with respect
to a one dimensional Brownian motion $\\{B_{x}\\}_{x\in[0,1]}$; this produces
a sequence of multiple Itô integrals $\\{I_{n}(\rho_{n}(t,\cdot))\\}_{n\geq
0}$ which solves a new set of equations, equivalent to (1.2) but expressed in
terms of differential second quantization operators, Malliavin derivatives and
their adjoints (see equation (4.3) below);

  * •

we define $\Phi(t):=\sum_{n\geq 0}I_{n}(\rho_{n}(t;\cdot))$; this is a
generalized stochastic process, solution to a single infinite dimensional
differential equation, see (4.6) below, and whose kernels (from its Wiener-Itô
chaos expansion) are by construction the elements of the sequence
$\\{\rho_{n}\\}_{n\geq 1}$. In this Gaussian framework one can rewrite the
adjoint Malliavin derivative, which appears in the equation for
$\\{\Phi(t)\\}_{t\geq 0}$, as a difference between multiplication operator and
Malliavin derivative;

  * •

using this transformation, the equation for $\\{\Phi(t)\\}_{t\geq 0}$ is
reduced to a single parabolic equation of Ornstein-Uhlenbeck-type with
infinitely many variables; the kernels of the solution to suitable finite
dimensional projections of such equation provide, through the Stroock -Taylor
formula, an analytic representation for the corresponding projections of the
sequence $\\{\rho_{n}\\}_{n\geq 0}$ (see formula (1.10) below).

$\\{\rho_{n}\\}_{n\geq 0}$ solves
(1.2)-(1.3)$\\{I_{n}(\rho_{n}(t,\cdot))\\}_{n\geq 0}$ solves
(4.3)$\Phi(t):=\sum_{n\geq 0}I_{n}(\rho_{n}(t,\cdot))$ solves
(4.6)$\rho_{n}(t,x_{1},...,x_{n})=\frac{1}{n!}\mathbb{E}[D_{x_{1},...,x_{n}}u(t)]$$\Phi(t)=u(t,I_{1}(\xi_{1}),I_{1}(\xi_{2},...))$
with $u$ from (4.17) Figure 1: Representation of $\\{\rho_{n}\\}_{n\geq 0}$ as
Wiener-Itô kernels of the stochastic process $\\{\Phi(t)\\}_{t\geq 0}$

Problem (1.2)-(1.3) corresponds to one particular instance of the class of
problems formalized in the recent paper [3]. The aim of that work was to
develop a general framework for stochastic particle-based reaction-diffusion
processes in the form of an evolution equation for the probability density of
the open system. We refer readers to the nice account, presented in [3], of
the vaste literature concerning the mathematical modelling of chemical and
biochemical phenomena which combine diffusion and chemical reactions.
Compared to the model from [3] which system (1.2)-(1.3) refers to, we made
some simplifying assumptions. In [3] the authors assume that the single
particle can move in the open bounded region
$\mathbb{X}\subset\mathbb{R}^{3}$, so that $\rho_{n}$ is defined on
$[0,+\infty[\times\mathbb{X}^{\otimes n}$; in addition, they deal with a non
necessarily isotropic diffusion with drift term. Here, for simplicity we focus
on a one dimensional case, i.e. $\mathbb{X}=]0,1[$, and assume drift-less
isotropic diffusion but our technique readily extends to the general case.
Contrary to the $L^{1}$-framework adopted in [3] and utilized to set up the
Fock-space formalism for describing and analysing the problem, we restrict to
the $L^{2}$-space (which is legitimate by assuming the existence of a
classical solution for (1.2)-(1.3)) and fully exploit the Wiener-Itô-Segal
isomorphism between the symmetric Fock space based on a Hilbert space and the
family of square integral Brownian functionals.

To state our main result we introduce a few notations and a couple of standing
assumptions.

######  Assumption 1.1.

There exists an orthonormal basis $\\{\xi_{k}\\}_{k\geq 1}$ of $L^{2}([0,1])$
that diagonalizes the operator

$\displaystyle\mathcal{A}:=-\partial^{2}_{x}+\lambda_{d}(x),\quad x\in[0,1],$
(1.7)

with homogenous Neumann boundary conditions. This means that for all $j,k\geq
1$ we have

$\displaystyle\int_{0}^{1}\xi_{k}(y)\xi_{j}(y)dy=\delta_{kj},\quad\xi_{k}^{\prime}(0)=\xi_{k}^{\prime}(1)=0,$

and there exists a sequence of non negative real numbers
$\\{\alpha_{k}\\}_{k\geq 1}$ such that

$\displaystyle\mathcal{A}\xi_{k}=\alpha_{k}\xi_{k},\quad\mbox{ for all $k\geq
1$}.$

######  Remark 1.2.

We note that according to the classical Sturm-Liouville boundary value problem
the continuity of $\lambda_{d}$ is sufficient for Assumption 1.1 to hold true.

We now denote by $\Pi_{N}:L^{2}([0,1])\to L^{2}([0,1])$ the orthogonal
projection onto the finite dimensional space spanned by
$\\{\xi_{1},...,\xi_{N}\\}$, i.e.

$\displaystyle\Pi_{N}f(x):=\sum_{k=1}^{N}\langle
f,\xi_{k}\rangle_{L^{2}([0,1])}\xi_{k}(x),\quad x\in[0,1];$

we also set

$\displaystyle d_{k}:=\langle\lambda_{d},\xi_{k}\rangle_{L^{2}([0,1])},\quad
c_{k}:=\langle\lambda_{c},\xi_{k}\rangle_{L^{2}([0,1])},\quad\gamma:=\int_{0}^{1}\lambda_{c}(y)dy,\quad\mbox{
and }\quad\zeta_{k}:=\langle\zeta,\xi_{k}\rangle_{L^{2}([0,1])},$ (1.8)

where the functions $\lambda_{d}$, $\lambda_{c}$ and $\zeta$ are those from
(1.2)-(1.3).

######  Assumption 1.3.

There exists $N_{0}\geq 1$ such that $\Pi_{N_{0}}\lambda_{d}=\lambda_{d}$;
this is equivalent to say $\Pi_{N}\lambda_{d}=\lambda_{d}$ for all $N\geq
N_{0}$.

######  Remark 1.4.

Assumption (1.3) is readily fulfilled in the case of a constant function
$\lambda_{d}(x)=\lambda_{d}1(x),x\in[0,1]$; in fact, in this case

$\displaystyle(\mathcal{A}f)(x)=-f^{\prime\prime}(x)+\lambda_{d}f(x),\quad\quad\xi_{k}(x)=\cos((k-1)\pi
x),\quad k\geq 1,$

and

$\displaystyle\alpha_{k}=(k-1)^{2}\pi^{2}+\lambda_{d},\quad k\geq 1.$

This gives

$\displaystyle\xi_{1}(x)=1(x)\quad\mbox{ and hence
}\quad(\Pi_{1}\lambda_{d})(x)=\lambda_{d}(x),$

i.e. $N_{0}=1$.

We are now ready to state the main result of the present paper; its proof is
postponed to Section 4 but a direct verification of its validity is presented
in Section 2 for $n=0,1,2$. In the sequel we set $\Pi_{N}^{\otimes n}$ to be
the orthogonal projection from $L^{2}([0,1]^{n})$ to the linear space
generated by the functions
$\\{\xi_{i_{1}}\otimes\cdot\cdot\cdot\otimes\xi_{i_{n}},1\leq
i_{1},...,i_{n}\leq N\\}$.

######  Theorem 1.5.

Let Assumptions 1.1-1.3 be in force and denote by $\\{\rho_{n}\\}_{n\geq 0}$ a
classical solution of equation (1.2)-(1.3). Then, for any $N\geq N_{0}$ and
$t\geq 0$ we have the representation

$\displaystyle\rho_{0}(t)=\mathbb{E}[u(t,Z)],$ (1.9)

and for any $n\geq 1$ and $(x_{1},...,x_{n})\in[0,1]^{n}$,

$\displaystyle\Pi_{N}^{\otimes
n}\rho_{n}(t,x_{1},...,x_{n})=\frac{1}{n!}\sum_{j_{1},...j_{n}=1}^{N}\mathbb{E}\left[\left(\partial_{z_{j_{1}}}\cdot\cdot\cdot\partial_{z_{j_{n}}}u\right)(t,Z)\right]\xi_{j_{1}}(x_{1})\cdot\cdot\cdot\xi_{j_{n}}(x_{n}).$
(1.10)

Here,

$\displaystyle\mathbb{E}\left[\left(\partial_{z_{j_{1}}}\cdot\cdot\cdot\partial_{z_{j_{n}}}u\right)(t,Z)\right]=\int_{\mathbb{R}^{N}}\left(\partial_{z_{j_{1}}}\cdot\cdot\cdot\partial_{z_{j_{n}}}u\right)(t,z)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz,$
(1.11)

while $u:[0,+\infty[\times\mathbb{R}^{N}\to\mathbb{R}$ is a classical solution
of the partial differential equation

$\displaystyle\begin{split}\partial_{t}u(t,z)=&\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}u(t,z)+\sum_{k=1}^{N}\left(d_{k}-c_{k}-\alpha_{k}z_{k}\right)\partial_{z_{k}}u(t,z)+\left(\sum_{k=1}^{N}c_{k}z_{k}-\gamma\right)u(t,z)\\\
u(0,z)=&\sum_{k=1}^{N}\zeta_{k}z_{k},\quad t\geq
0,z\in\mathbb{R}^{N}.\end{split}$ (1.12)

The paper is organized as follows: in the next section we describe how to
verify the validity of formula (1.10) through a direct computation: this
should help the reader in understanding the mechanism that relates (1.2)-(1.3)
to (1.12); Section 3 describes the Gaussian setting needed to formalize our
approach: here we recall few basic ideas and tools from Malliavin calculus and
infinite dimensional Gaussian analysis; lastly, in Section 4 we prove formula
(1.10), passing through several intermediate steps that illustrate the main
ideas of our technique.

##  2 Verification of formula (1.9)-(1.10) for $n=0,1,2$

The aim of this section is to show via a direct verification the validity of
formula (1.9)-(1.10). This will be done only for $n=0,1,2$ and serves as an
illustration of the connection between (1.2)-(1.3) and (1.12). First of all,
using the notation (1.7) we rewrite (1.2) as

$\displaystyle\begin{split}\partial_{t}\rho_{0}(t)=&\int_{0}^{1}\lambda_{d}(y)\rho_{1}(t,y)dy-\gamma\rho_{0}(t);\\\
\partial_{t}\rho_{n}(t,x_{1},...,x_{n})=&-\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})\\\
&+(n+1)\int_{0}^{1}\lambda_{d}(y)\rho_{n+1}(t,x_{1},...,x_{n},y)dy\\\
&+\frac{1}{n}\sum_{i=1}^{n}\lambda_{c}(x_{i})\rho_{n-1}(t,x_{1},...,x_{i-1},x_{i+1},...,x_{n})\\\
&-\gamma\rho_{n}(t,x_{1},...,x_{n}).\end{split}$ (2.1)

Then, we find the equation solved by $\\{\Pi_{N}^{\otimes n}\rho_{n}\\}_{n\geq
0}$.

######  Proposition 2.1.

Let Assumptions 1.1-1.3 be in force and denote by $\\{\rho_{n}\\}_{n\geq 0}$ a
classical solution of equation (2.1)-(1.3). Then, for any $N\geq N_{0}$ the
sequence $\\{\Pi_{N}^{\otimes n}\rho_{n}\\}_{n\geq 0}$ solves

$\displaystyle\begin{split}\partial_{t}\rho_{0}(t)=&\int_{0}^{1}\lambda_{d}(y)\Pi_{N}\rho_{1}(t,y)dy-\gamma\rho_{0}(t);\\\
\partial_{t}\Pi_{N}^{\otimes
n}\rho_{n}(t,x_{1},...,x_{n})=&-\sum_{j_{1},...j_{n}=1}^{N}\left(\sum_{i=1}^{n}\alpha_{j_{i}}\right)\langle\rho_{n}(t,\cdot),\xi_{j_{1}}\otimes\cdot\cdot\cdot\otimes\xi_{j_{n}}\rangle_{L^{2}([0,1]^{n})}\xi_{j_{1}}(x_{1})\cdot\cdot\cdot\xi_{j_{n}}(x_{n})\\\
&+(n+1)\int_{0}^{1}\lambda_{d}(y)\Pi_{N}^{\otimes(n+1)}\rho_{n+1}(t,x_{1},...,x_{n},y)dy\\\
&+\frac{1}{n}\sum_{i=1}^{n}\Pi_{N}\lambda_{c}(x_{i})\Pi_{N}^{\otimes(n-1)}\rho_{n-1}(t,x_{1},...,x_{i-1},x_{i+1},...,x_{n})\\\
&-\gamma\Pi_{N}^{\otimes n}\rho_{n}(t,x_{1},...,x_{n}).\end{split}$ (2.2)

Here, we set $\Pi_{N}^{\otimes 0}\rho_{0}:=\rho_{0}$.

###### Proof.

Since $\\{\rho_{n}\\}_{n\geq 0}$ solves (2.1), we can write for $n\geq 1$ that

$\displaystyle\begin{split}\partial_{t}\Pi_{N}^{\otimes
n}\rho_{n}(t,x_{1},...,x_{n})=&\Pi_{N}^{\otimes
n}\partial_{t}\rho_{n}(t,x_{1},...,x_{n})\\\ =&-\Pi_{N}^{\otimes
n}\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})\\\
&+(n+1)\int_{0}^{1}\lambda_{d}(y)\Pi_{N}^{\otimes
n}\rho_{n+1}(t,x_{1},...,x_{n},y)dy\\\
&+\frac{1}{n}\sum_{i=1}^{n}\Pi_{N}^{\otimes
n}\left[\lambda_{c}(x_{i})\rho_{n-1}(t,x_{1},...,x_{i-1},x_{i+1},...,x_{n})\right]\\\
&-\gamma\Pi_{N}^{\otimes n}\rho_{n}(t,x_{1},...,x_{n})\\\ =&-\Pi_{N}^{\otimes
n}\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})\\\
&+(n+1)\int_{0}^{1}\lambda_{d}(y)\Pi_{N}^{\otimes
n}\rho_{n+1}(t,x_{1},...,x_{n},y)dy\\\
&+\frac{1}{n}\sum_{i=1}^{n}\Pi_{N}\lambda_{c}(x_{i})\Pi_{N}^{\otimes(n-1)}\rho_{n-1}(t,x_{1},...,x_{i-1},x_{i+1},...,x_{n})\\\
&-\gamma\Pi_{N}^{\otimes n}\rho_{n}(t,x_{1},...,x_{n}).\end{split}$ (2.3)

We now observe that

$\displaystyle\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})=$
$\displaystyle\sum_{i=1}^{n}\mathcal{A}_{i}\sum_{j_{1},...,j_{n}\geq
1}\langle\rho_{n}(t,\cdot),\xi_{j_{1}}\otimes\cdot\cdot\cdot\otimes\xi_{j_{n}}\rangle_{L^{2}([0,1]^{n})}\xi_{j_{1}}(x_{1})\cdot\cdot\cdot\xi_{j_{n}}(x_{n})$
$\displaystyle=$ $\displaystyle\sum_{j_{1},...,j_{n}\geq
1}\left(\sum_{i=1}^{n}\alpha_{j_{i}}\right)\langle\rho_{n}(t,\cdot),\xi_{j_{1}}\otimes\cdot\cdot\cdot\otimes\xi_{j_{n}}\rangle_{L^{2}([0,1]^{n})}\xi_{j_{1}}(x_{1})\cdot\cdot\cdot\xi_{j_{n}}(x_{n}),$

and hence

$\displaystyle\Pi_{N}^{\otimes
n}\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})$
$\displaystyle\quad=\sum_{j_{1},...,j_{n}=1}^{N}\left(\sum_{i=1}^{n}\alpha_{j_{i}}\right)\langle\rho_{n}(t,\cdot),\xi_{j_{1}}\otimes\cdot\cdot\cdot\otimes\xi_{j_{n}}\rangle_{L^{2}([0,1]^{n})}\xi_{j_{1}}(x_{1})\cdot\cdot\cdot\xi_{j_{n}}(x_{n}).$

Moreover, according to Assumption 1.3 we have

$\displaystyle\int_{0}^{1}\lambda_{d}(y)\Pi_{N}^{\otimes
n}\rho_{n+1}(t,x_{1},...,x_{n},y)dy=$
$\displaystyle\int_{0}^{1}\Pi_{N}\lambda_{d}(y)\Pi_{N}^{\otimes
n}\rho_{n+1}(t,x_{1},...,x_{n},y)dy$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\lambda_{d}(y)\Pi_{N}^{\otimes(n+1)}\rho_{n+1}(t,x_{1},...,x_{n},y)dy.$

If we employ these two last facts in (2.3), we arrive at (2.2) for $n\geq 1$.
Similarly, the equation for $n=0$ can be derived by virtue of Assumption 1.3:

$\displaystyle\partial_{t}\rho_{0}(t)=$
$\displaystyle\int_{0}^{1}\lambda_{d}(y)\rho_{1}(t,y)dy-\gamma\rho_{0}(t)$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\Pi_{N}\lambda_{d}(y)\rho_{1}(t,y)dy-\gamma\rho_{0}(t)$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\lambda_{d}(y)\Pi_{N}\rho_{1}(t,y)dy-\gamma\rho_{0}(t).$

∎

######  Remark 2.2.

The initial and boundary conditions for the sequence $\\{\Pi_{N}^{\otimes
n}\rho_{n}\\}_{n\geq 0}$ are easily deduced from (1.3) and the corresponding
boundary conditions for $\\{\xi_{k}\\}_{k\geq 1}$; more precisely,

$\displaystyle\begin{split}\rho_{0}(0)&=0;\\\
\Pi_{N}\rho_{1}(0,x_{1})&=\Pi_{N}\zeta(x_{1}),\mbox{ for all
$x_{1}\in[0,1]$};\\\ \Pi_{N}^{\otimes
n}\rho_{n}(0,x_{1},...,x_{n})&=0,\quad\mbox{for all $n>1$ and all
$(x_{1},...,x_{n})\in[0,1]^{n}$};\\\ \emph{grad }(\Pi_{N}^{\otimes
n}\rho_{n})(t,x_{1},...,x_{n})&=0,\quad\mbox{for all $n\geq 1$, $t\geq 0$ and
$(x_{1},...,x_{n})\in\partial[0,1]^{n}$}.\end{split}$ (2.4)

######  Remark 2.3.

It is useful to recall that the solution to the Cauchy problem

$\displaystyle\begin{split}\partial_{t}u(t,z)&=\sum_{i=1}^{N}\alpha_{i}\partial^{2}_{z_{i}}u(t,z)+\sum_{i=1}^{N}\left(d_{i}-c_{i}-\alpha_{i}z_{i}\right)\partial_{z_{i}}u(t,z)+\left(\sum_{i=1}^{N}c_{i}z_{i}-\gamma\right)u(t,z);\\\
u(0,z)&=\sum_{i=1}^{N}\zeta_{i}z_{i},\end{split}$ (2.5)

which is the key ingredient of formulas (1.9)-(1.10), admits the following
Feynman-Kac representation (see for instance [7]):

$\displaystyle
u(t,z)=\mathbb{E}\left[\left(\sum_{i=1}^{N}\zeta_{i}Z^{z_{i}}_{i}(t)\right)\exp\left\\{\int_{0}^{t}\left(\sum_{i=1}^{N}c_{i}Z^{z_{i}}_{i}(s)-\gamma\right)ds\right\\}\right],\quad
t\geq 0,z=(z_{1},...,z_{N})\in\mathbb{R}^{N}.$ (2.6)

Here, for $i\in\\{1,...,N\\}$, the stochastic process
$\\{Z_{i}^{z_{i}}(t)\\}_{t\geq 0}$ is the unique strong solution of the mean-
reverting Ornstein-Uhlenbeck stochastic differential equation

$\displaystyle
dZ_{i}^{z_{i}}(t)=\left(d_{i}-c_{i}-\alpha_{i}Z^{z_{i}}_{i}(t)\right)dt+\sqrt{2\alpha_{i}}dW_{i}(t),\quad
Z_{i}^{z_{i}}(0)=z_{i},$ (2.7)

with $\\{W_{1}(t)\\}_{t\geq 0}$,…, $\\{W_{N}(t)\\}_{t\geq 0}$ being
independent one dimensional Brownian motions. It is well known that the
solution to (2.7) can be explicitly written for $\alpha_{i}>0$ as

$\displaystyle Z^{z_{i}}_{i}(t)$
$\displaystyle=z_{i}e^{-\alpha_{i}t}+\frac{d_{i}-c_{i}}{\alpha_{i}}\left(1-e^{-\alpha_{i}t}\right)+\int_{0}^{t}e^{-\alpha_{i}(t-s)}\sqrt{2\alpha_{i}}dW_{i}(s),$

and simply

$\displaystyle Z^{z_{i}}_{i}(t)$ $\displaystyle=z_{i}+(d_{i}-c_{1})t,$

when $\alpha_{i}=0$. This shows that the function $z_{i}\mapsto
Z^{z_{i}}_{i}(t)$ is almost surely affine and hence that, according to
equation (2.6), for any $\tau>0$ there exist positive constant $m_{1}$ and
$m_{2}$ such that

$\displaystyle|u(t,z)|\leq m_{1}e^{m_{2}|z|},\quad\mbox{for all $t\in[0,\tau]$
and $z\in\mathbb{R}^{N}$}.$ (2.8)

This bound entails the finiteness of the expectation in (1.9); since the same
reasoning applies to the partial spatial derivatives of $u$ (they also satisfy
an equation of the form (2.5)), we conclude that the expectations in (1.10)
are well defined and finite as well.

We are now going to verify that the right hand sides of (1.9)-(1.10) solve
equation (2.2) for $n=0,1,2$; to ease the reference to formula (1.9)-(1.10),
we write down explicitly the first four terms:

$\displaystyle\begin{split}\rho_{0}(t)&=\mathbb{E}[u(t,Z)];\\\
\Pi_{N}\rho_{1}(t,x_{1})&=\sum_{j=1}^{N}\mathbb{E}\left[(\partial_{z_{j}}u)(t,Z)\right]\xi_{j}(x_{1});\\\
\Pi_{N}^{\otimes
2}\rho_{2}(t,x_{1},x_{2})&=\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\mathbb{E}\left[(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2});\\\
\Pi_{N}^{\otimes
3}\rho_{3}(t,x_{1},x_{2},x_{3})&=\frac{1}{6}\sum_{j_{1},j_{2},j_{3}=1}^{N}\mathbb{E}\left[(\partial_{z_{j_{3}}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})\xi_{j_{3}}(x_{3}).\end{split}$
(2.9)

###  2.1 Equation for $\rho_{0}$

We have:

$\displaystyle\partial_{t}\mathbb{E}[u(t,Z)]=$
$\displaystyle\mathbb{E}[(\partial_{t}u)(t,Z)]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}u(t,Z)+\sum_{k=1}^{N}\left(d_{k}-c_{k}-\alpha_{k}Z_{k}\right)\partial_{z_{k}}u(t,Z)+\left(\sum_{k=1}^{N}c_{k}Z_{k}-\gamma\right)u(t,Z)\right]$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{N}\alpha_{k}\mathbb{E}[\partial^{2}_{z_{k}}u(t,Z)]+\sum_{k=1}^{N}\mathbb{E}\left[\left(d_{k}-c_{k}-\alpha_{k}Z_{k}\right)\partial_{z_{k}}u(t,Z)\right]$
$\displaystyle+\mathbb{E}\left[\left(\sum_{k=1}^{N}c_{k}Z_{k}-\gamma\right)u(t,Z)\right].$

Observe that an integration by parts gives (the boundary term vanishes thanks
to the bound (2.8) applied to $\partial_{z_{k}}u(t,z)$)

$\displaystyle\mathbb{E}[\partial^{2}_{z_{k}}u(t,Z)]=\mathbb{E}[Z_{k}\partial_{z_{k}}u(t,Z)];$

therefore,

$\displaystyle\partial_{t}\mathbb{E}[u(t,Z)]=$
$\displaystyle\sum_{k=1}^{N}\alpha_{k}\mathbb{E}[Z_{k}\partial_{z_{k}}u(t,Z)]+\sum_{k=1}^{N}\mathbb{E}\left[\left(d_{k}-c_{k}-\alpha_{k}Z_{k}\right)\partial_{z_{k}}u(t,Z)\right]$
$\displaystyle+\mathbb{E}\left[\left(\sum_{k=1}^{N}c_{k}Z_{k}-\gamma\right)u(t,Z)\right]$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{N}\mathbb{E}\left[\left(d_{k}-c_{k}\right)\partial_{z_{k}}u(t,Z)\right]+\mathbb{E}\left[\left(\sum_{k=1}^{N}c_{k}Z_{k}-\gamma\right)u(t,Z)\right]$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{N}d_{k}\mathbb{E}\left[\partial_{z_{k}}u(t,Z)\right]-\sum_{k=1}^{N}c_{k}\mathbb{E}\left[\partial_{z_{k}}u(t,Z)\right]+\mathbb{E}\left[\left(\sum_{k=1}^{N}c_{k}Z_{k}-\gamma\right)u(t,Z)\right].$

An additional integration by parts yields

$\displaystyle\sum_{k=1}^{N}c_{k}\mathbb{E}\left[\partial_{z_{k}}u(t,Z)\right]=\mathbb{E}\left[\left(\sum_{k=1}^{N}c_{k}Z_{k}\right)u(t,Z)\right],$

and hence

$\displaystyle\partial_{t}\mathbb{E}[u(t,Z)]=$
$\displaystyle\sum_{k=1}^{N}d_{k}\mathbb{E}\left[\partial_{z_{k}}u(t,Z)\right]-\gamma\mathbb{E}\left[u(t,Z)\right]$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\lambda_{d}(y)\left(\sum_{j=1}^{N}\mathbb{E}\left[(\partial_{z_{j}}u)(t,Z)\right]\xi_{j}(y)\right)dy-\gamma\mathbb{E}[u(t,Z)].$

This corresponds to equation (2.2) for $n=0$, with
$\rho_{0}(t)=\mathbb{E}[u(t,Z)]$ and

$\displaystyle\Pi_{N}\rho_{1}(t,x_{1})=\sum_{j=1}^{N}\mathbb{E}\left[(\partial_{z_{j}}u)(t,Z)\right]\xi_{j}(x_{1}).$

###  2.2 Equation for $\Pi_{N}\rho_{1}$

We have

$\displaystyle\partial_{t}\sum_{j=1}^{N}\mathbb{E}\left[(\partial_{z_{j}}u)(t,Z)\right]\xi_{j}(x_{1})=\sum_{j=1}^{N}\mathbb{E}\left[(\partial_{z_{j}}\partial_{t}u)(t,Z)\right]\xi_{j}(x_{1}).$
(2.10)

Now,

$\displaystyle\mathbb{E}\left[(\partial_{z_{j}}\partial_{t}u)(t,Z)\right]=$
$\displaystyle\int_{\mathbb{R}^{N}}(\partial_{z_{j}}\partial_{t}u)(t,z)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{N}}\partial_{z_{j}}\left(\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}u(t,z)+\sum_{k=1}^{N}\left(d_{k}-c_{k}-\alpha_{k}z_{k}\right)\partial_{z_{k}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle+\int_{\mathbb{R}^{N}}\partial_{z_{j}}\left(\left(\sum_{k=1}^{N}c_{k}z_{k}-\gamma\right)u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{N}}\left(\sum_{k=1}^{N}\alpha_{k}\partial_{z_{j}}\partial^{2}_{z_{k}}u(t,z)-\alpha_{j}\partial_{z_{j}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle+\int_{\mathbb{R}^{N}}\left(\sum_{k=1}^{N}\left(d_{k}-c_{k}-\alpha_{k}z_{k}\right)\partial_{z_{j}}\partial_{z_{k}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle+\int_{\mathbb{R}^{N}}\left(c_{j}u(t,z)+\left(\sum_{k=1}^{N}c_{k}z_{k}-\gamma\right)\partial_{z_{j}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz.$

An integration by parts with respect to $\partial_{z_{k}}$ in the first term
of the fourth line above will produce the term

$\displaystyle\int_{\mathbb{R}^{N}}\left(\sum_{k=1}^{N}\alpha_{k}z_{k}\partial_{z_{j}}\partial_{z_{k}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz,$

which is identical to one of the terms from the fifth line but opposite in
sign. We can therefore write

$\displaystyle\begin{split}\mathbb{E}\left[(\partial_{z_{j}}\partial_{t}u)(t,Z)\right]=&\int_{\mathbb{R}^{N}}\left(-\alpha_{j}\partial_{z_{j}}u(t,z)+\sum_{k=1}^{N}\left(d_{k}-c_{k}\right)\partial_{z_{j}}\partial_{z_{k}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz\\\
&+\int_{\mathbb{R}^{N}}\left(c_{j}u(t,z)+\left(\sum_{k=1}^{N}c_{k}z_{k}-\gamma\right)\partial_{z_{j}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz.\end{split}$
(2.11)

Similarly, an integration by parts with respect to $\partial_{z_{k}}$ in the
term

$\displaystyle\int_{\mathbb{R}^{N}}\left(-\sum_{k=1}^{N}c_{k}\partial_{z_{j}}\partial_{z_{k}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$

from the first line of (2.11) will give

$\displaystyle-\int_{\mathbb{R}^{N}}\left(\sum_{k=1}^{N}c_{k}z_{k}\right)\partial_{z_{j}}u(t,z)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz.$

hence cancelling the corresponding term from the second line in (2.11).
Summing up,

$\displaystyle\mathbb{E}\left[(\partial_{z_{j}}\partial_{t}u)(t,Z)\right]=$
$\displaystyle\int_{\mathbb{R}^{N}}\left(-\alpha_{j}\partial_{z_{j}}u(t,z)+\sum_{k=1}^{N}d_{k}\partial_{z_{j}}\partial_{z_{k}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle+\int_{\mathbb{R}^{N}}\left(c_{j}u(t,z)-\gamma\partial_{z_{j}}u(t,z)\right)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle=$
$\displaystyle-\alpha_{j}\mathbb{E}[(\partial_{z_{j}}u)(t,Z)]+\sum_{k=1}^{N}d_{k}\mathbb{E}\left[(\partial_{z_{j}}\partial_{z_{k}}u)(t,Z)\right]$
$\displaystyle+c_{j}\mathbb{E}[u(t,Z)]-\gamma\mathbb{E}[(\partial_{z_{j}}u)(t,Z)].$

We can now plug the last expression in (2.10) to get

$\displaystyle\partial_{t}\sum_{j=1}^{N}\mathbb{E}\left[(\partial_{z_{j}}u)(t,Z)\right]\xi_{j}(x_{1})=$
$\displaystyle\sum_{j=1}^{N}\mathbb{E}\left[(\partial_{z_{j}}\partial_{t}u)(t,Z)\right]\xi_{j}(x_{1})$
$\displaystyle=$
$\displaystyle-\sum_{j=1}^{N}\alpha_{j}\mathbb{E}[(\partial_{z_{j}}u)(t,Z)]\xi_{j}(x_{1})+\sum_{j=1}^{N}\sum_{k=1}^{N}d_{k}\mathbb{E}\left[(\partial_{z_{j}}\partial_{z_{k}}u)(t,Z)\right]\xi_{j}(x_{1})$
$\displaystyle+\sum_{j=1}^{N}c_{j}\mathbb{E}[u(t,Z)]\xi_{j}(x_{1})-\gamma\sum_{j=1}^{N}\mathbb{E}[(\partial_{z_{j}}u)(t,Z)]\xi_{j}(x_{1})$
$\displaystyle=$
$\displaystyle-\sum_{j=1}^{N}\alpha_{j}\mathbb{E}[(\partial_{z_{j}}u)(t,Z)]\xi_{j}(x_{1})+2\int_{0}^{1}\lambda_{d}(y)\Pi_{N}^{\otimes
2}\rho_{2}(t,x_{1},y)dy$
$\displaystyle+\Pi_{N}\lambda_{c}(x_{1})\rho_{0}(t)-\gamma\Pi_{N}\rho_{1}(t,x_{1}).$

This corresponds to equation (2.2) for $n=1$ with the prescriptions (2.9).

###  2.3 Equation for $\Pi_{N}^{\otimes 2}\rho_{2}$

We start as before with

$\displaystyle\begin{split}&\partial_{t}\left(\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\mathbb{E}\left[(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})\right)\\\
&=\quad\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\mathbb{E}\left[(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{t}u)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2}).\end{split}$
(2.12)

Now,

$\displaystyle\mathbb{E}\left[(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{t}u)(t,Z)\right]$
$\displaystyle\quad=\int_{\mathbb{R}^{N}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{t}u(t,z)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad=\int_{\mathbb{R}^{N}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\left(\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}u(t,z)+\sum_{k=1}^{N}\left(d_{k}-c_{k}-\alpha_{k}z_{k}\right)\partial_{z_{k}}u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad\quad+\int_{\mathbb{R}^{N}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\left(\left(\sum_{k=1}^{N}c_{k}z_{k}-\gamma\right)u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad=\sum_{k=1}^{N}\alpha_{k}\int_{\mathbb{R}^{N}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial^{2}_{z_{k}}u(t,z)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad\quad+\int_{\mathbb{R}^{N}}\left(-\alpha_{j_{1}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u(t,z)-\alpha_{j_{2}}\partial_{z_{j_{1}}}\partial_{z_{j_{2}}}u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad\quad+\int_{\mathbb{R}^{N}}\left(\sum_{k=1}^{N}\left(d_{k}-c_{k}-\alpha_{k}z_{k}\right)\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{z_{k}}u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad\quad+\int_{\mathbb{R}^{N}}\left(c_{j_{1}}\partial_{z_{j_{2}}}u(t,z)+c_{j_{2}}\partial_{z_{j_{1}}}u(t,z)+\left(\sum_{k=1}^{N}c_{k}z_{k}-\gamma\right)\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz.$

Integration by parts with respect to $\partial_{z_{k}}$ in

$\displaystyle\int_{\mathbb{R}^{N}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial^{2}_{z_{k}}u(t,z)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$

and

$\displaystyle\int_{\mathbb{R}^{N}}\left(-\sum_{k=1}^{N}c_{k}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{z_{k}}u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$

will simplify the expression to

$\displaystyle\mathbb{E}\left[(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{t}u)(t,Z)\right]$
$\displaystyle\quad=\int_{\mathbb{R}^{N}}\left(-\alpha_{j_{1}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u(t,z)-\alpha_{j_{2}}\partial_{z_{j_{1}}}\partial_{z_{j_{2}}}u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad\quad+\sum_{k=1}^{N}d_{k}\int_{\mathbb{R}^{N}}\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{z_{k}}u(t,z)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad\quad+\int_{\mathbb{R}^{N}}\left(c_{j_{1}}\partial_{z_{j_{2}}}u(t,z)+c_{j_{2}}\partial_{z_{j_{1}}}u(t,z)-\gamma\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u(t,z)\right)(2\pi)^{-\frac{N}{2}}e^{-\frac{|z|^{2}}{2}}dz$
$\displaystyle\quad=-\left(\alpha_{j_{1}}+\alpha_{j_{2}}\right)\mathbb{E}\left[\left(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u\right)(t,Z)\right]+\sum_{k=1}^{N}d_{k}\mathbb{E}\left[\left(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{z_{k}}u\right)(t,Z)\right]$
$\displaystyle\quad\quad+c_{j_{1}}\mathbb{E}\left[(\partial_{z_{j_{2}}}u)(t,Z)\right]+c_{j_{2}}\mathbb{E}\left[(\partial_{z_{j_{1}}}u)(t,Z)\right]-\gamma\mathbb{E}\left[\left(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u\right)(t,Z)\right]$

Therefore,

$\displaystyle\partial_{t}\Pi_{N}^{\otimes 2}\rho_{2}(t,x_{1},x_{2})=$
$\displaystyle\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\mathbb{E}\left[(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{t}u)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\left(\alpha_{j_{1}}+\alpha_{j_{2}}\right)\mathbb{E}\left[\left(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u\right)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})$
$\displaystyle+\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\sum_{k=1}^{N}d_{k}\mathbb{E}\left[\left(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}\partial_{z_{k}}u\right)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})$
$\displaystyle+\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\left(c_{j_{1}}\mathbb{E}\left[(\partial_{z_{j_{2}}}u)(t,Z)\right]+c_{j_{2}}\mathbb{E}\left[(\partial_{z_{j_{1}}}u)(t,Z)\right]\right)\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})$
$\displaystyle-\frac{\gamma}{2}\sum_{j_{1},j_{2}=1}^{N}\mathbb{E}\left[\left(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u\right)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N}\left(\alpha_{j_{1}}+\alpha_{j_{2}}\right)\mathbb{E}\left[\left(\partial_{z_{j_{2}}}\partial_{z_{j_{1}}}u\right)(t,Z)\right]\xi_{j_{1}}(x_{1})\xi_{j_{2}}(x_{2})$
$\displaystyle+3\int_{0}^{1}\lambda_{d}(y)\Pi_{N}^{\otimes
3}\rho_{3}(t,x_{1},x_{2},y)dy+(\Pi_{N}\lambda_{c}\hat{\otimes}\Pi_{N}\rho_{1}(t,\cdot))(x_{1},x_{2})$
$\displaystyle-\gamma\Pi_{N}^{\otimes 2}\rho_{2}(t,x_{1},x_{2}).$

This corresponds through identities (2.9) to equation (2.2) for $n=2$.

##  3 Preliminary material

In this section we introduce the framework utilized for proving our main
theorem. For more details on these topics, we refer the reader to one of the
books [2], [6] and [10].

###  3.1 Wiener-Itô chaos expansion

Let $(\Omega,\mathcal{B},\mathbb{P})$ be the classical Wiener space over the
interval $[0,1]$, i.e. $\Omega$ is the space of continuous functions defined
on the interval $[0,1]$ and null at zero, $\mathcal{B}$ is the Borel
$\sigma$-algebra of $\Omega$ induced by the supremum norm and $\mathbb{P}$ the
Wiener measure on $(\Omega,\mathcal{B})$. We denote by

$\displaystyle\begin{split}B_{x}&:\Omega\to\mathbb{R}\\\ &\quad\omega\mapsto
B_{x}(\omega):=\omega(x),\quad x\in[0,1],\end{split}$

the coordinate process which by construction is a one dimensional Brownian
motion under $\mathbb{P}$. According to the Wiener-Itô chaos expansion
theorem, any random variable $\Phi$ in $\mathbb{L}^{2}(\Omega)$ can be
uniquely represented as

$\displaystyle\Phi=\sum_{n\geq 0}I_{n}(h_{n}),$ (3.1)

where

  * •

$h_{0}:=\mathbb{E}[\Phi]$;

  * •

for $n\geq 1$, $h_{n}\in L_{s}^{2}([0,1]^{n})$, the space of square integrable
symmetric functions;

  * •

$I_{0}(h_{0}):=h_{0}=\mathbb{E}[\Phi]$;

  * •

for $n\geq 1$, $I_{n}(h_{n})$ stands for the $n$-th order multiple Itô
integral defined as

$\displaystyle
I_{n}(h_{n}):=n!\int_{0}^{1}\int_{0}^{x_{1}}\cdot\cdot\cdot\int_{0}^{x_{n-1}}h_{n}(x_{1},...,x_{n})dB_{x_{n}}\cdot\cdot\cdot
dB_{x_{2}}dB_{x_{1}}.$

The series in (3.1) provides an orthogonal decomposition of $\Phi$ that
converges in $\mathbb{L}^{2}(\Omega)$; in fact, multiple Itô integrals possess
the following general properties:

  * •

for all $n\geq 1$, $\mathbb{E}[I_{n}(h_{n})]=0$;

  * •

if $n\neq m$, then $\mathbb{E}[I_{n}(h_{n})I_{m}(h_{m})]=0$;

  * •

for $n\geq 1$,
$\mathbb{E}[I_{n}(h_{n})^{2}]=n!|h_{n}|^{2}_{L^{2}([0,1]^{n})}$.

From the last two identities we get

$\displaystyle\mathbb{E}[\Phi^{2}]=\sum_{n\geq
0}n!|h_{n}|^{2}_{L^{2}([0,1]^{n})},$

and, for $\Psi\in\mathbb{L}^{2}(\Omega)$ with

$\displaystyle\Psi=\sum_{n\geq 0}I_{n}(g_{n}),$ (3.2)

that

$\displaystyle\mathbb{E}[\Phi\Psi]=\sum_{n\geq 0}n!\langle
h_{n},g_{n}\rangle_{L^{2}([0,1]^{n})}.$ (3.3)

Two notable subsets of $\mathbb{L}^{2}(\Omega)$ are

$\displaystyle\mathtt{F}:=\left\\{\sum_{n=0}^{M}I_{n}(h_{n}),\mbox{ for some
$M\in\mathbb{N}\cup\\{0\\}$, $h_{0}\in\mathbb{R}$ and $h_{n}\in
L_{s}^{2}([0,1]^{n})$, $n=1,...,M$}\right\\},$

which collects the random variables with a finite order chaos expansion, and

$\displaystyle\mathtt{E}:=\left\\{\mathcal{E}(f):=\sum_{n\geq
0}I_{n}\left(\frac{f^{\otimes n}}{n!}\right),\mbox{ for some $f\in
L^{2}([0,1])$}\right\\},$

which is the family of the so-called _stochastic exponentials_. It is well
known that

$\displaystyle\mathcal{E}(f)=\exp\left\\{I_{1}(f)-\frac{1}{2}|f|^{2}_{L^{2}([0,1])}\right\\}$

and that $\mathtt{F}$ and the linear span of $\mathtt{E}$ are both dense in
$\mathbb{L}^{2}(\Omega)$. In particular,

$\displaystyle\mathbb{E}[\Phi Z]=\mathbb{E}[\Psi Z],\quad\mbox{for all
$Z\in\mathtt{F}$},$

or

$\displaystyle\mathbb{E}[\Phi\mathcal{E}(f)]=\mathbb{E}[\Psi\mathcal{E}(f)],\quad\mbox{for
all $f\in L^{2}([0,1])$ (or some dense subset of $L^{2}([0,1])$)},$

implies $\Phi=\Psi$, $\mathbb{P}$-a.s.. Note also that, according to (3.3), we
can write

$\displaystyle\mathbb{E}[\Phi\mathcal{E}(f)]=\sum_{n\geq 0}\langle
h_{n},f^{\otimes n}\rangle_{L^{2}([0,1]^{n})},$

whenever $\Phi=\sum_{n\geq 0}I_{n}(h_{n})$. We recall in addition that, by
virtue of the _Hu-Meyer formula_

$\displaystyle I_{n}(h_{n})\cdot I_{m}(h_{m})=\sum_{r=0}^{n\wedge
n}r!{n\choose r}{m\choose r}I_{n+m-2r}(h_{n}\hat{\otimes}_{r}h_{m}),$ (3.4)

the linear space $\mathtt{F}$ is an algebra with respect to the point-wise
multiplication. Here, $h_{n}\otimes_{r}h_{m}$ stands for the _$r$ -th order
contraction_ of $h_{n}$ and $h_{m}$, i.e.

$\displaystyle\begin{split}&(h_{n}\otimes_{r}h_{m})(x_{1},....,x_{n+m-2r})\\\
&:=\int_{[0,1]^{r}}h_{n}(x_{1},...,x_{n-r},y_{1},...,y_{r})h_{m}(y_{1},...,y_{r},x_{n-r+1},...,x_{n+m-2r})dy_{1}\cdot\cdot\cdot
dy_{r},\end{split}$ (3.5)

while $h_{n}\hat{\otimes}_{r}h_{m}$ denotes the _symmetrization_ of
$h_{n}\otimes_{r}h_{m}$, i.e.

$\displaystyle\begin{split}&(h_{n}\hat{\otimes}_{r}h_{m})(x_{1},....,x_{n+m-2r})\\\
&:=\frac{1}{(n+m-2r)!}\sum_{\sigma\in
S_{n+m-2r}}(h_{n}\otimes_{r}h_{m})(x_{\sigma(1)},...,x_{\sigma(n+m-2r)}),\end{split}$
(3.6)

with $S_{n+m-2r}$ being the group of permutations on $\\{1,...,n+m-2r\\}$.

###  3.2 Malliavin derivative

The _Malliavin derivative_ of $\Phi=\sum_{n=0}^{M}I_{n}(h_{n})\in\mathtt{F}$,
denoted $\\{D_{x}\Phi\\}_{x\in[0,1]}$, is the element of
$L^{2}([0,1];\mathtt{F})$ defined by

$\displaystyle D_{x}\Phi:=\sum_{n=0}^{M-1}(n+1)I_{n}(h_{n+1}(\cdot,x)),\quad
x\in[0,1].$

For $l\in L^{2}([0,1])$ and $\Phi=\sum_{n=0}^{M}I_{n}(h_{n})\in\mathtt{F}$, we
also write

$\displaystyle\begin{split}D_{l}\Phi:=\langle
D\Phi,l\rangle_{L^{2}([0,1])}&=\sum_{n=0}^{M-1}(n+1)I_{n}\left(\int_{0}^{1}h_{n+1}(\cdot,y)l(y)dy\right)\\\
&=\sum_{n=0}^{M-1}(n+1)I_{n}\left(h_{n+1}\otimes_{1}l\right)\end{split}$ (3.7)

for the _directional Malliavin derivative_ of $\Phi$ along $l$ (in the last
member above we utilized the notation (3.5)). We remark that $D_{l}\Phi$ is
also a member of $\mathtt{F}$.
If we now take $l\in L^{2}([0,1])$,
$\Phi=\sum_{n=0}^{M}I_{n}(h_{n})\in\mathtt{F}$ and
$\Psi=\sum_{n=0}^{K}I_{n}(g_{n})\in\mathtt{F}$, we can write

$\displaystyle\mathbb{E}[D_{l}\Phi\cdot\Psi]$
$\displaystyle=\sum_{n=0}^{(M-1)\wedge K}n!(n+1)\langle
h_{n+1}\otimes_{1}l,g_{n}\rangle_{L^{2}([0,1]^{n})}$
$\displaystyle=\sum_{n=0}^{(M-1)\wedge K}(n+1)!\langle h_{n+1},l\otimes
g_{n}\rangle_{L^{2}([0,1]^{n+1})}$ $\displaystyle=\sum_{n=0}^{(M-1)\wedge
K}(n+1)!\langle h_{n+1},l\hat{\otimes}g_{n}\rangle_{L^{2}([0,1]^{n+1})}$ (3.8)
$\displaystyle=\sum_{n=1}^{M\wedge(K+1)}n!\langle
h_{n},l\hat{\otimes}g_{n-1}\rangle_{L^{2}([0,1]^{n})}$
$\displaystyle=\mathbb{E}[\Phi\cdot D_{l}^{\star}\Psi]$ ,

where

$\displaystyle
D^{\star}_{l}\Psi:=\sum_{n=1}^{K+1}I_{n}(l\hat{\otimes}g_{n-1})$

and

$\displaystyle(l\hat{\otimes}g_{n-1})(x_{1},...,x_{n}):=\frac{1}{n}\sum_{i=1}^{n}f(x_{i})g_{n-1}(x_{1},...,x_{i-1},x_{i+1},...,x_{n})$

(compare with definitions (3.5) and (3.6) for $r=0$). We remark that in (3.2)
we utilized the symmetry of $h_{n}$ and the fact that the symmetrization
operator is idempotent and self-adjoint in $L^{2}([0,1]^{n})$.
It is clear that $D^{\star}_{l}\Psi$ also belongs to $\mathtt{F}$; moreover,
$\mathbb{E}[D^{\star}_{l}\Psi]=0$, for all $l\in L^{2}([0,1])$ and
$\Psi\in\mathtt{F}$.
If in the Hu-Meyer formula (3.4) we take $m=1$ and set $g_{1}=l\in
L^{2}([0,1])$, we get

$\displaystyle I_{n}(h_{n})\cdot I_{1}(l)$
$\displaystyle=I_{n+1}(h_{n}\hat{\otimes}l)+I_{n-1}(h_{n}\otimes_{1}l)$
$\displaystyle=D_{l}^{\star}I_{n}(h_{n})+D_{l}I_{n}(h_{n}).$

Summing over $n$ and using the linearity of the operators $D^{\star}_{l}$ and
$D_{l}$, we deduce the identity

$\displaystyle D_{l}^{\star}\Psi+D_{l}\Psi=\Psi\cdot I_{1}(l),$ (3.9)

which is valid for $\Psi\in\mathtt{F}$ and $l\in L^{2}([0,1])$. One can also
introduce the adjoint of $D_{x}$, denoted $\delta$: if
$\Phi(x)=\sum_{n=0}^{M}I_{n}(h_{n}(\cdot;x))$ is a stochastic process in
$\mathtt{F}$, then

$\displaystyle\delta(\Phi(\cdot)):=\sum_{n=0}^{M}I_{n+1}(\tilde{h}_{n})\in\mathtt{F},$

where $\tilde{h}_{n}$ stands for the symmetrization of $h_{n}$ with respect to
the $n+1$ variables $x_{1},...,x_{n},x$.
It is useful to mention that the definition of Malliavin derivative can also
be extended to the members of the family $\mathtt{E}$: for any $f,l\in
L^{2}([0,1])$, we have

$\displaystyle D_{x}\mathcal{E}(f)=f(x)\mathcal{E}(f),x\in[0,1]\quad\mbox{ and
}\quad D_{l}\mathcal{E}(f)=\langle f,l\rangle_{L^{2}([0,1])}\mathcal{E}(f).$
(3.10)

###  3.3 Second quantization operators

Let $A:L^{2}([0,1])\to L^{2}([0,1])$ be a bounded linear operator; for
$\Phi=\sum_{n=0}^{M}I_{n}(h_{n})\in\mathtt{F}$ we define the _second
quantization operator_ of $A$ as

$\displaystyle\Gamma(A)\Phi:=\sum_{n=0}^{M}I_{n}\left(A^{\otimes
n}h_{n}\right),$

and the _differential second quantization operator_ of $A$ as

$\displaystyle
d\Gamma(A)\Phi:=\sum_{n=1}^{M}I_{n}\left(\sum_{i=1}^{n}A_{i}h_{n}\right),$

where $A_{i}$ stands for the operator $A$ acting on the $i$-th variable of
$h_{n}$. The boundedness of $A$ implies that both $\Gamma(A)\Phi$ and
$d\Gamma(A)\Psi$ also belong to $\mathtt{F}$; note in addition that for $A$
being the identity, we recover from $d\Gamma(A)$ the well known _number
operator_ :

$\displaystyle\mathcal{N}\Phi=\sum_{n=1}^{M}nI_{n}\left(h_{n}\right).$

Via a simple verification on can see that for all $\Phi$ and $\Psi$ in
$\mathtt{F}$ the following identities hold true:

$\displaystyle\begin{split}&\mathbb{E}[\Gamma(A)\Phi]=\mathbb{E}[\Phi];\quad\mathbb{E}[d\Gamma(A)\Phi]=0;\\\
&\mathbb{E}[\Gamma(A)\Phi\cdot\Psi]=\mathbb{E}[\Phi\cdot\Gamma(A^{\star})\Psi];\quad\mathbb{E}[d\Gamma(A)\Phi\cdot\Psi]=\mathbb{E}[\Phi\cdot
d\Gamma(A^{\star})\Psi].\end{split}$

Here, $A^{\star}$ denotes the adjoint of $A$ in $L^{2}([0,1])$. As for the
Malliavin derivative, the actions of second quantization and differential
second quantization operators can be extended to the class $\mathtt{E}$ of
stochastic exponentials:

$\displaystyle\Gamma(A)\mathcal{E}(f)=\mathcal{E}(Af)\quad\mbox{ and }\quad
d\Gamma(A)\mathcal{E}(f)=D^{\star}_{Af}\mathcal{E}(f).$

The differential second quantization operator can also be represented as a
composition of the Malliavin derivative and its adjoint; more precisely,

$\displaystyle d\Gamma(A)\Phi=\delta\left(AD_{\cdot}\Phi\right).$ (3.11)

###  3.4 Space of generalized random variables

We remark that the convergence of the series in (3.3) is implied via the
Cauchy-Schwartz inequality by the conditions

$\displaystyle\sum_{n\geq
0}n!|h_{n}|^{2}_{L^{2}([0,1]^{n})}<+\infty\quad\mbox{ and }\quad\sum_{n\geq
0}n!|g_{n}|^{2}_{L^{2}([0,1]^{n})}<+\infty.$ (3.12)

However, if $\Psi$ has a finite order expansion, i.e.

$\displaystyle\Psi=\sum_{n=0}^{M}I_{n}(g_{n}),\mbox{ for some
$M\in\mathbb{N}\cup\\{0\\}$},$ (3.13)

then one can drop the first condition in (3.12) and get a still well-defined
pairing between the _generalized_ random variable $\Phi$, represented by the
formal series $\sum_{n\geq 0}I_{n}(h_{n})$ (which in general will not
convergence in $\mathbb{L}^{2}(\Omega)$), and the _regular_ or _test_ random
variable $\Psi$ with finite order expansion (3.13). Let

$\displaystyle\mathtt{F}^{\star}:=\left\\{\sum_{n\geq 0}I_{n}(h_{n}),\mbox{
for some $h_{0}\in\mathbb{R}$ and $h_{n}\in L_{s}^{2}([0,1]^{n})$, $n\geq
1$}\right\\}$

be a family of _generalized_ random variables. The action of $T=\sum_{n\geq
0}I_{n}(h_{n})\in\mathtt{F}^{\star}$ on
$\varphi=\sum_{n=0}^{M}I_{n}(g_{n})\in\mathtt{F}$ is defined as

$\displaystyle\langle\langle T,\varphi\rangle\rangle:=\sum_{n=0}^{M}n!\langle
h_{n},g_{n}\rangle_{L^{2}([0,1]^{n})}.$

By construction, we have the inclusions

$\displaystyle\mathtt{F}\subset\mathbb{L}^{2}(\Omega)\subset\mathtt{F}^{\star}$

with

$\displaystyle\langle\langle T,\varphi\rangle\rangle=\mathbb{E}[T\varphi],$

whenever $T\in\mathbb{L}^{2}(\Omega)$. We will say that $T=U$ in
$\mathtt{F}^{\star}$ if

$\displaystyle\langle\langle T,\varphi\rangle\rangle=\langle\langle
U,\varphi\rangle\rangle,\quad\mbox{for all }\varphi\in\mathtt{F}.$

The _generalized expectation_ of $T=\sum_{n\geq
0}I_{n}(h_{n})\in\mathtt{F}^{\star}$ is $\mathbb{E}[T]:=\langle\langle
T,1\rangle\rangle=h_{0}$. It is also important to observe that, according to
the Hu-Meyer formula (3.4), the vector space $\mathtt{F}$ is closed with
respect to the point-wise product between random variables. Therefore, if
$T\in\mathtt{F}^{\star}$ and $\psi\in\mathtt{F}$, the product $T\cdot\psi$ is
well defined and corresponds to the element of $\mathtt{F}^{\star}$ given by
the prescription

$\displaystyle\langle\langle T\cdot\psi,\varphi\rangle\rangle:=\langle\langle
T,\psi\cdot\varphi\rangle\rangle,\quad\varphi\in\mathtt{F}.$

The definitions of Malliavin derivative, its adjoint and (differential) second
quantization operators can be lifted from $\mathtt{F}$ to $\mathtt{F}^{\star}$
by duality:

$\displaystyle\langle\langle D_{l}T,\varphi\rangle\rangle:=\langle\langle
T,D^{\star}_{l}\varphi\rangle\rangle,\quad\langle\langle
D^{\star}_{l}T,\varphi\rangle\rangle:=\langle\langle
T,D_{l}\varphi\rangle\rangle$
$\displaystyle\langle\langle\Gamma(A)T,\varphi\rangle\rangle:=\langle\langle
T,\Gamma(A^{\star})\varphi\rangle\rangle,\quad\langle\langle
d\Gamma(A)T,\varphi\rangle\rangle:=\langle\langle
T,d\Gamma(A^{\star})\varphi\rangle\rangle.$

Lastly, we recall a generalized version of the so-called Stroock-Taylor
formula: if $T=\sum_{n\geq 0}I_{n}(h_{n})\in\mathtt{F}^{\star}$, then

$\displaystyle
h_{n}(x_{1},...,x_{n})=\frac{1}{n!}\mathbb{E}[D_{x_{1},...,x_{n}}T],\quad(x_{1},...,x_{n})\in[0,1]^{n}.$
(3.14)

Here, $\mathbb{E}[D_{x_{1},...,x_{n}}T]$ stands for the the generalized
expectation of the $n$-th order Malliavin derivative of $T$.

######  Remark 3.1.

The space $\mathtt{F}^{\star}$ has been already utilized for solving some
stochastic partial differential equations which admits only generalized
solutions. See for instance the paper [9] where the authors investigate an
unbiased version of the stochastic Navier-Stokes equations.

##  4 Proof of Theorem 1.5

In this section we present the rigorous derivation of formula (1.9)-(1.10).
Firstly, we devote our attention to the original system, i.e. the one without
projection operators $\Pi_{N}^{\otimes n}$, that we report here for easiness
of reference:

$\displaystyle\begin{split}\partial_{t}\rho_{0}(t)=&\int_{0}^{1}\lambda_{d}(y)\rho_{1}(t,y)dy-\gamma\rho_{0}(t);\\\
\partial_{t}\rho_{n}(t,x_{1},...,x_{n})=&-\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})\\\
&+(n+1)\int_{0}^{1}\lambda_{d}(y)\rho_{n+1}(t,x_{1},...,x_{n},y)dy\\\
&+\frac{1}{n}\sum_{i=1}^{n}\lambda_{c}(x_{i})\rho_{n-1}(t,x_{1},...,x_{i-1},x_{i+1},...,x_{n})\\\
&-\gamma\rho_{n}(t,x_{1},...,x_{n}),\end{split}$ (4.1)

with initial and boundary conditions

$\displaystyle\begin{split}\rho_{0}(0)&=0;\\\
\rho_{1}(0,x_{1})&=\zeta(x_{1}),\quad x_{1}\in[0,1];\\\
\rho_{n}(0,x_{1},...,x_{n})&=0,\quad n>1,(x_{1},...,x_{n})\in[0,1]^{n};\\\
\partial_{\nu}\rho_{n}(t,x_{1},...,x_{n})&=0,\quad n\geq 1,t\geq
0,(x_{1},...,x_{n})\in\partial[0,1]^{n}.\end{split}$ (4.2)

The first fundamental step of our analysis consists in integrating all the
spatial variables of $\rho_{n}(t,x_{1},...,x_{n})$ with respect to the one
dimensional Brownian motion $\\{B_{x}\\}_{x\in[0,1]}$; this procedure will
produce a sequence of time-dependent multiple Itô integrals which satisfies a
stochastic counterpart of equation (4.1)-(4.2).

######  Proposition 4.1.

Let $\\{\rho_{n}\\}_{n\geq 0}$ be a classical solution to (4.1)-(4.2). Then,
the sequence of random variables $\\{I_{n}(\rho_{n}(t,\cdot))\\}_{n\geq 0}$
satisfies the equations

$\displaystyle\begin{split}\partial_{t}I_{n}(\rho_{n}(t,\cdot))=&d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))+D_{\lambda_{d}}I_{n+1}\left(\rho_{n+1}(t,\cdot)\right)\\\
&+D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma
I_{n}(\rho_{n}(t,\cdot)),\quad t>0,n\geq 0;\\\
I_{1}(\rho_{1}(0,\cdot))=&I_{1}(\zeta);\\\ I_{n}(\rho_{n}(0,\cdot))=&0,\mbox{
for all $n\neq 1$},\end{split}$ (4.3)

with probability one. Here, we agree on setting $I_{-1}(\cdot)\equiv 0$.

###### Proof.

Using the first order contraction, see (3.5), and symmetrized tensor product,
see (3.6), we can reformulate (4.1) as

$\displaystyle\begin{split}\partial_{t}\rho_{0}(t)=&\lambda_{d}\otimes_{1}\rho_{1}(t,\cdot)-\gamma\rho_{0}(t);\\\
\partial_{t}\rho_{n}(t,x_{1},...,x_{n})=&-\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})+(n+1)(\lambda_{d}\otimes_{1}\rho_{n+1}(t,\cdot))(x_{1},...,x_{n})\\\
&+(\lambda_{c}\hat{\otimes}\rho(t,\cdot))(x_{1},...,x_{n})-\gamma\rho_{n}(t,x_{1},...,x_{n}).\end{split}$
(4.4)

The continuity and symmetry of the functions $\rho_{n}(t,x_{1},...,x_{n})$,
$\partial_{t}\rho_{n}(t,x_{1},...,x_{n})$ and
$\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,x_{1},...,x_{n})$ entail their
membership to $L^{2}_{s}([0,1]^{n})$; this allows us to perform, for any
$n\geq 1$, an $n$-th order multiple Itô integral on both sides of the equation
(4.4) to get

$\displaystyle\begin{split}I_{n}(\partial_{t}\rho_{n}(t,\cdot))=&-I_{n}\left(\sum_{i=1}^{n}\mathcal{A}_{i}\rho_{n}(t,\cdot)\right)+(n+1)I_{n}\left(\lambda_{d}\otimes_{1}\rho_{n+1}(t,\cdot)\right)\\\
&+I_{n}\left(\lambda_{c}\hat{\otimes}\rho_{n-1}(t,\cdot)\right)-\gamma
I_{n}(\rho_{n}(t,\cdot)),\end{split}$

or equivalently,

$\displaystyle\begin{split}\partial_{t}I_{n}(\rho_{n}(t,\cdot))=&d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))+D_{\lambda_{d}}I_{n+1}\left(\rho_{n+1}(t,\cdot)\right)\\\
&+D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma
I_{n}(\rho_{n}(t,\cdot)).\end{split}$

The last identity holds for all $n\geq 1$, $t>0$, $\mathbb{P}$-almost surely.
The initial conditions in (4.3) are readily checked. ∎

Our next step is to construct a generalized stochastic process
$\\{\Phi(t)\\}_{t\geq 0}$ out of the sequence
$\\{I_{n}(\rho_{n}(t,\cdot))\\}_{n\geq 0}$ in the spirit of the Wiener-Itô
chaos expansion.

######  Proposition 4.2.

Let $\\{\rho_{n}\\}_{n\geq 0}$ be a classical solution to (4.1)-(4.2). Then,
the stochastic process

$\displaystyle\Phi(t):=\sum_{n\geq 0}I_{n}(\rho_{n}(t,\cdot)),\quad t\geq 0,$
(4.5)

belongs to $\mathtt{F}^{\star}$ and solves the differential equation

$\displaystyle\begin{split}\partial_{t}\Phi(t)=&d\Gamma(-\mathcal{A})\Phi(t)+D_{\lambda_{d}}\Phi(t)+D^{\star}_{\lambda_{c}}\Phi(t)-\gamma\Phi(t),\quad
t>0,\\\ \Phi(0)=&I_{1}(\zeta).\end{split}$ (4.6)

###### Proof.

The continuity of $\rho_{n}(t,\cdot)$ together with its partial derivatives up
to the second order ensure that the dual pairing
$\langle\langle\Phi(t),\varphi\rangle\rangle$ between $\Phi(t)$ from (4.5) and
any element $\varphi=\sum_{n=0}^{M}I_{n}(h_{n})$ of $\mathtt{F}$ is well
defined thus entailing the membership of $\Phi(t)$ to $\mathtt{F}^{\star}$,
for all $t\geq 0$. To prove that $\\{\Phi(t)\\}_{t\geq 0}$ solves (4.6), we
sum equation (4.3) over $n\geq 0$ with the convention that $I_{-1}(\cdot):=0$
(recall also that $D_{\lambda_{d}}I_{0}(\cdot)=0$): this will lead to

$\displaystyle\partial_{t}\Phi(t)=$ $\displaystyle\partial_{t}\sum_{n\geq
0}I_{n}(\rho_{n}(t,\cdot))=\sum_{n\geq 0}\partial_{t}I_{n}(\rho_{n}(t,\cdot))$
$\displaystyle=$ $\displaystyle\sum_{n\geq
0}d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))+\sum_{n\geq
0}D_{\lambda_{d}}I_{n+1}\left(\rho_{n+1}(t,\cdot)\right)$
$\displaystyle+\sum_{n\geq
0}D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma\sum_{n\geq
0}I_{n}(\rho_{n}(t,\cdot))$ $\displaystyle=$ $\displaystyle
d\Gamma(-\mathcal{A})\sum_{n\geq
0}I_{n}(\rho_{n}(t,\cdot))+D_{\lambda_{d}}\sum_{n\geq
0}I_{n+1}\left(\rho_{n+1}(t,\cdot)\right)$
$\displaystyle+D^{\star}_{\lambda_{c}}\sum_{n\geq
0}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma\sum_{n\geq
0}I_{n}(\rho_{n}(t,\cdot))$ $\displaystyle=$ $\displaystyle
d\Gamma(-\mathcal{A})\sum_{n\geq
0}I_{n}(\rho_{n}(t,\cdot))+D_{\lambda_{d}}\sum_{n\geq
1}I_{n}\left(\rho_{n}(t,\cdot)\right)$
$\displaystyle+D^{\star}_{\lambda_{c}}\sum_{n\geq
1}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma\sum_{n\geq
0}I_{n}(\rho_{n}(t,\cdot))$ $\displaystyle=$ $\displaystyle
d\Gamma(-\mathcal{A})\Phi(t)+D_{\lambda_{d}}\Phi(t)+D^{\star}_{\lambda_{c}}\Phi(t)-\gamma\Phi(t).$

Observe that the smoothness of the functions $\\{\rho_{n}\\}_{n\geq 0}$ allows
for the interchange between the series and operators $\partial_{t}$,
$d\Gamma(-\mathcal{A})$, $D_{\lambda_{d}}$ and $D^{\star}_{\lambda_{c}}$
(since we are working in the space $\mathtt{F}^{\star}$). The initial
condition in (4.6) is trivially verified and the proof is complete. ∎

######  Remark 4.3.

The stochastic process $\\{\Phi(t)\\}_{t\geq 0}$ and differential equation
(4.6) provide a concise reformulation of the sequence $\\{\rho_{n}\\}_{n\geq
0}$ and system of partial differential equations (4.1)-(4.2). In principle,
one may start solving equation (4.6) and then identify, via the generalized
Stroock-Taylor formula (3.14), the kernels of this solution as the sequence
$\\{\rho_{n}\\}_{n\geq 0}$ fulfilling (4.1)-(4.2).

We now proceed with the investigation of the projected sequence
$\\{\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot)\\}_{n\geq 0}$.

######  Proposition 4.4.

Let $\\{\rho_{n}\\}_{n\geq 0}$ be a classical solution to (4.1)-(4.2). Then,
for any $N\geq N_{0}$ the sequence of random variables
$\\{I_{n}(\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot))\\}_{n\geq 0}$ satisfies for
any $t\geq 0$ the equations

$\displaystyle\begin{split}\partial_{t}I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))=&d\Gamma(-\mathcal{A})I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))+D_{\lambda_{d}}I_{n+1}\left(\Pi_{N}^{\otimes(n+1)}\rho_{n+1}(t,\cdot)\right)\\\
&+D^{\star}_{\Pi_{N}\lambda_{c}}I_{n-1}\left(\Pi_{N}^{\otimes(n-1)}\rho_{n-1}(t,\cdot)\right)-\gamma
I_{n}(\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot));\\\
I_{1}(\Pi_{N}\rho_{1}(0,\cdot))=&I_{1}(\Pi_{N}\zeta);\\\
I_{n}(\Pi_{N}^{\otimes n}\rho_{n}(0,\cdot))=&0,\mbox{ for all $n\neq
1$},\end{split}$ (4.7)

with probability one.

###### Proof.

We know from Proposition 4.1 that the sequence of random variables
$\\{I_{n}(\rho_{n}(t,\cdot))\\}_{n\geq 0}$ satisfies equation (4.3); moreover,
according to the definition of second quantization operator we can write

$\displaystyle I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))=\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot)).$

Therefore,

$\displaystyle\begin{split}\partial_{t}I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))=&\partial_{t}\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))\\\
=&\Gamma(\Pi_{N})\partial_{t}I_{n}(\rho_{n}(t,\cdot))\\\
=&\Gamma(\Pi_{N})\left[d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))+D_{\lambda_{d}}I_{n+1}\left(\rho_{n+1}(t,\cdot)\right)\right]\\\
&+\Gamma(\Pi_{N})\left[D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma
I_{n}(\rho_{n}(t,\cdot))\right]\\\
=&\Gamma(\Pi_{N})d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))+\Gamma(\Pi_{N})D_{\lambda_{d}}I_{n+1}\left(\rho_{n+1}(t,\cdot)\right)\\\
&+\Gamma(\Pi_{N})D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot)).\end{split}$
(4.8)

We now have to investigate the commutation relations of $\Gamma(\Pi_{N})$ with
$d\Gamma(-\mathcal{A})$, $D_{\lambda_{d}}$ and $D^{\star}_{\lambda_{c}}$. Let
$h\in C^{2}([0,1])$; then,

$\displaystyle\begin{split}\mathbb{E}[\Gamma(\Pi_{N})d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]&=\mathbb{E}[d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(\Pi_{N}h)]\\\
&=\mathbb{E}[I_{n}(\rho_{n}(t,\cdot))\cdot
d\Gamma(-\mathcal{A})\mathcal{E}(\Pi_{N}h)]\\\
&=\mathbb{E}[I_{n}(\rho_{n}(t,\cdot))\cdot
D^{\star}_{-\mathcal{A}\Pi_{N}h}\mathcal{E}(\Pi_{N}h)]\\\
&=\mathbb{E}[I_{n}(\rho_{n}(t,\cdot))\cdot
D^{\star}_{-\Pi_{N}\mathcal{A}h}\mathcal{E}(\Pi_{N}h)].\end{split}$ (4.9)

On the other hand,

$\displaystyle\begin{split}\mathbb{E}[d\Gamma(-\mathcal{A})\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]&=\mathbb{E}[\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))\cdot
d\Gamma(-\mathcal{A})\mathcal{E}(h)]\\\
&=\mathbb{E}[\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))\cdot
D^{\star}_{-\mathcal{A}h}\mathcal{E}(h)]\\\
&=\mathbb{E}[I_{n}(\rho_{n}(t,\cdot))\cdot\Gamma(\Pi_{N})D^{\star}_{-\mathcal{A}h}\mathcal{E}(h)]\\\
&=\mathbb{E}[I_{n}(\rho_{n}(t,\cdot))\cdot
D^{\star}_{-\Pi_{N}\mathcal{A}h}\mathcal{E}(\Pi_{N}h)].\end{split}$ (4.10)

Comparing the first and last members of (4.9) and (4.10) we deduce that

$\displaystyle\mathbb{E}[\Gamma(\Pi_{N})d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]=\mathbb{E}[d\Gamma(-\mathcal{A})\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)],$

for all $h\in C^{2}([0,1])$ and hence that

$\displaystyle\Gamma(\Pi_{N})d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))=d\Gamma(-\mathcal{A})\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot)),\quad\mathbb{P}-\mbox{a.s}.$
(4.11)

We now consider the commutation between $\Gamma(\Pi_{N})$ and
$D_{\lambda_{d}}$:

$\displaystyle\mathbb{E}[\Gamma(\Pi_{N})D_{\lambda_{d}}I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]$
$\displaystyle=\mathbb{E}[D_{\lambda_{d}}I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(\Pi_{N}h)]$
$\displaystyle=\mathbb{E}[nI_{n-1}(\lambda_{d}\otimes_{1}\rho_{n}(t,\cdot))\cdot\mathcal{E}(\Pi_{N}h)]$
$\displaystyle=n\langle\lambda_{d}\otimes_{1}\rho_{n}(t,\cdot),(\Pi_{N}h)^{\otimes{n-1}}\rangle_{L^{2}([0,1]^{n-1})}$
$\displaystyle=n\langle\lambda_{d}\otimes_{1}\Pi_{N}^{\otimes(n-1)}\rho_{n}(t,\cdot),h^{\otimes{n-1}}\rangle_{L^{2}([0,1]^{n-1})},$

while

$\displaystyle\mathbb{E}[D_{\lambda_{d}}\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]$
$\displaystyle=\mathbb{E}[D_{\lambda_{d}}I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]$
$\displaystyle=\mathbb{E}[nI_{n-1}(\lambda_{d}\otimes_{1}\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]$
$\displaystyle=\mathbb{E}[nI_{n-1}(\Pi_{N}\lambda_{d}\otimes_{1}\Pi_{N}^{\otimes(n-1)}\rho_{n}(t,\cdot))\cdot\mathcal{E}(h)]$
$\displaystyle=n\langle\Pi_{N}\lambda_{d}\otimes_{1}\Pi_{N}^{\otimes(n-1)}\rho_{n}(t,\cdot),h^{\otimes{n-1}}\rangle_{L^{2}([0,1]^{n-1})}.$

Therefore,

$\displaystyle\Gamma(\Pi_{N})D_{\lambda_{d}}I_{n}(\rho_{n}(t,\cdot))=D_{\lambda_{d}}\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot)),$
(4.12)

if $\lambda_{d}=\Pi_{N}\lambda_{d}$; but this is the case since we are
assuming $N\geq N_{0}$ (recall Assumption 1.3). Lastly, we investigate the
commutation between $\Gamma(\Pi_{N})$ and $D^{\star}_{\lambda_{c}}$:

$\displaystyle\mathbb{E}[\Gamma(\Pi_{N})D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)\cdot\mathcal{E}(h)]=$
$\displaystyle\mathbb{E}[D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)\cdot\mathcal{E}(\Pi_{N}h)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)\cdot
D_{\lambda_{c}}\mathcal{E}(\Pi_{N}h)]$ $\displaystyle=$
$\displaystyle\mathbb{E}[I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)\cdot\langle\lambda_{c},\Pi_{N}h\rangle\mathcal{E}(\Pi_{N}h)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)\cdot\langle\Pi_{N}\lambda_{c},h\rangle\mathcal{E}(\Pi_{N}h)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\Gamma(\Pi_{N})I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)\cdot\langle\Pi_{N}\lambda_{c},h\rangle\mathcal{E}(h)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[D^{\star}_{\Pi_{N}\lambda_{c}}\Gamma(\Pi_{N})I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)\cdot\mathcal{E}(h)];$

comparing the first and last members we can conclude that

$\displaystyle\Gamma(\Pi_{N})D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)=D^{\star}_{\Pi_{N}\lambda_{c}}\Gamma(\Pi_{N})I_{n-1}\left(\rho_{n-1}(t,\cdot)\right),$
(4.13)

almost surely. Now, using (4.11)-(4.12)-(4.13) in (4.8) we obtain

$\displaystyle\partial_{t}I_{n}(\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot))=$
$\displaystyle\Gamma(\Pi_{N})d\Gamma(-\mathcal{A})I_{n}(\rho_{n}(t,\cdot))+\Gamma(\Pi_{N})D_{\lambda_{d}}I_{n+1}\left(\rho_{n+1}(t,\cdot)\right)$
$\displaystyle+\Gamma(\Pi_{N})D^{\star}_{\lambda_{c}}I_{n-1}\left(\rho_{n-1}(t,\cdot)\right)-\gamma\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))$
$\displaystyle=$ $\displaystyle d\Gamma(-\mathcal{A})I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))+D_{\lambda_{d}}I_{n+1}\left(\Pi_{N}^{\otimes(n+1)}\rho_{n+1}(t,\cdot)\right)$
$\displaystyle+D^{\star}_{\Pi_{N}\lambda_{c}}I_{n-1}\left(\Pi_{N}^{\otimes(n-1)}\rho_{n-1}(t,\cdot)\right)-\gamma
I_{n}(\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot)),$

which corresponds to (4.7). ∎

The last step of our construction consists in rewriting equations (4.7) from a
standard partial differential equation’s perspective.

######  Proposition 4.5.

For any $n\geq 1$,

$\displaystyle I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))=\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right),$
(4.14)

where $\varphi_{n,N}:[0,+\infty[\times\mathbb{R}^{N}\to\mathbb{R}$ is a
polynomial of degree $n$ in the variables $I_{1}(\xi_{1})$,…,
$I_{1}(\xi_{N})$. With this representation, equations (4.7) read

$\displaystyle\begin{split}(\partial_{t}\varphi_{n,N})\left(t,z\right)=&\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}\varphi_{n,N}\left(t,z\right)-\sum_{k=1}^{N}\alpha_{k}z_{k}\partial_{z_{k}}\varphi_{n,N}\left(t,z\right)+\sum_{k=1}^{N}d_{k}\partial_{z_{k}}\varphi_{n+1,N}\left(t,z\right)\\\
&+\varphi_{n-1,N}\left(t,z\right)\cdot\left(\sum_{k=1}^{N}c_{k}z_{k}\right)-\sum_{k=1}^{N}c_{k}\partial_{z_{k}}\varphi_{n-1,N}\left(t,z\right)-\gamma\varphi_{n,N}(t,z),\end{split}$
(4.15)

for $n\geq 1$, $t\geq 0$ and $z=(z_{1},...,z_{N})\in\mathbb{R}^{N}$. Here we
set $\varphi_{-1,N}\equiv 0$.

###### Proof.

We start observing (see for instance [6]) that

$\displaystyle I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))=\Gamma(\Pi_{N})I_{n}(\rho_{n}(t,\cdot))=\mathbb{E}[I_{n}(\rho_{n}(t,\cdot))|\sigma(I_{1}(\xi_{1}),...,I_{1}(\xi_{N}))],$

where the right hand side above stands for the conditional expectation of the
random variable $I_{n}(\rho_{n}(t,\cdot))$ with respect to the sigma-algebra
generated by the random variables $\\{I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\\}$.
It is also well known from the theory of multiple Itô integrals that
$I_{n}(\rho_{n}(t,\cdot))$ can be written as an infinite linear combination of
polynomials of degree $n$ in the variables
$I_{1}(\xi_{1}),I_{1}(\xi_{2}),...$; the action of the conditional expectation
above reduce that linear combination to a finite number of terms in the
variables $\\{I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\\}$. This proves identity
(4.14). Let us now derive equation (4.15); according to formula (3.11) we have

$\displaystyle d\Gamma(-\mathcal{A})I_{n}(\Pi_{N}^{\otimes
n}\rho_{n}(t,\cdot))=$ $\displaystyle
d\Gamma(-\mathcal{A})\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle=$
$\displaystyle-\delta\left(\mathcal{A}D_{\cdot}\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)\right)$
$\displaystyle=$
$\displaystyle-\delta\left(\sum_{k=1}^{N}(\partial_{z_{k}}\varphi_{n,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)\mathcal{A}\xi_{k}(\cdot)\right)$
$\displaystyle=$
$\displaystyle-\sum_{k=1}^{N}\alpha_{k}\delta\left((\partial_{z_{k}}\varphi_{n,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)\xi_{k}(\cdot)\right)$
$\displaystyle=$
$\displaystyle-\sum_{k=1}^{N}\alpha_{k}\delta(\xi_{k})\cdot(\partial_{z_{k}}\varphi_{n,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle+\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right);$

moreover, using the chain rule for Malliavin derivatives we get

$\displaystyle
D_{\lambda_{d}}I_{n+1}\left(\Pi_{N}^{\otimes(n+1)}\rho_{n+1}(t,\cdot)\right)=$
$\displaystyle
D_{\lambda_{d}}\varphi_{n+1,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{N}(\partial_{z_{k}}\varphi_{n+1,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)\langle\lambda_{d},\xi_{k}\rangle_{L^{2}([0,1])}$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{N}d_{k}(\partial_{z_{k}}\varphi_{n+1,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right),$

and

$\displaystyle
D^{\star}_{\Pi_{N}\lambda_{c}}I_{n-1}\left(\Pi_{N}^{\otimes(n-1)}\rho_{n-1}(t,\cdot)\right)=$
$\displaystyle
D^{\star}_{\Pi_{N}\lambda_{c}}\varphi_{n-1,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle=$
$\displaystyle\varphi_{n-1,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)\cdot
I_{1}(\Pi_{N}\lambda_{c})$ $\displaystyle-
D_{\Pi_{N}\lambda_{c}}\varphi_{n-1,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle=$
$\displaystyle\varphi_{n-1,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)\cdot\left(\sum_{k=1}^{N}c_{k}I_{1}(\xi_{k})\right)$
$\displaystyle-\sum_{k=1}^{N}(\partial_{z_{k}}\varphi_{n-1,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)c_{k}.$

Note that in the second equality above we made use of identity (3.9).
Combining all the previous identities we can rewrite equations (4.7) as

$\displaystyle(\partial_{t}\varphi_{n,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)=$
$\displaystyle\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle-\sum_{k=1}^{N}\alpha_{k}\delta(\xi_{k})\cdot(\partial_{z_{k}}\varphi_{n,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle+\sum_{k=1}^{N}d_{k}(\partial_{z_{k}}\varphi_{n+1,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)$
$\displaystyle+\varphi_{n-1,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)\cdot\left(\sum_{k=1}^{N}c_{k}I_{1}(\xi_{k})\right)$
$\displaystyle-\sum_{k=1}^{N}(\partial_{z_{k}}\varphi_{n-1,N})\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)c_{k}$
$\displaystyle-\gamma\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right).$

This corresponds to (4.15) upon replacing $I_{1}(\xi_{k})$ with $z_{k}$, for
$k=1,...,N$. ∎

We are now ready to prove Theorem 1.5. Let

$\displaystyle u(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})):=\sum_{n\geq
0}\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right).$ (4.16)

According to the last proposition, summing over $n\geq 0$ equations (4.15) we
see that the function $u$ in (4.16) solves

$\displaystyle\begin{split}\partial_{t}u(t,z)&=\sum_{k=1}^{N}\alpha_{k}\partial^{2}_{z_{k}}u(t,z)+\sum_{k=1}^{N}\left(d_{k}-c_{k}-\alpha_{k}z_{k}\right)\partial_{z_{k}}u(t,z)+\left(\sum_{k=1}^{N}c_{k}z_{k}-\gamma\right)u(t,z)\\\
u(0,z)&=\sum_{k=1}^{N}\zeta_{k}z_{k}.\end{split}$

On the other hand, by construction

$\displaystyle\sum_{n\geq
0}\varphi_{n,N}\left(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N})\right)=\sum_{n\geq
0}I_{n}\left(\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot)\right);$

this gives

$\displaystyle u(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N}))=\sum_{n\geq
0}I_{n}\left(\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot)\right).$

Thus, the kernels in the Wiener-Itô chaos expansion of
$u(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N}))$, which in general can be represented
via the Stroock-Taylor formula (3.14), coincide with the sequence
$\\{\Pi_{N}^{\otimes n}\rho_{n}(t,\cdot)\\}_{n\geq 0}$. Therefore,

$\displaystyle\Pi_{N}^{\otimes n}\rho_{n}(t,x_{1},...,x_{n})=$
$\displaystyle\frac{1}{n!}\mathbb{E}\left[D_{x_{1},...,x_{n}}u(t,I_{1}(\xi_{1}),...,I_{1}(\xi_{N}))\right]$
$\displaystyle=$
$\displaystyle\frac{1}{n!}\sum_{j_{1},...j_{n}=1}^{N}\mathbb{E}\left[\left(\partial_{z_{j_{1}}}\cdot\cdot\cdot\partial_{z_{j_{n}}}u\right)(t,Z)\right]\xi_{j_{1}}(x_{1})\cdot\cdot\cdot\xi_{j_{n}}(x_{n}),$

where

$\displaystyle\mathbb{E}\left[\left(\partial_{z_{j_{1}}}\cdot\cdot\cdot\partial_{z_{j_{n}}}u\right)(t,Z)\right]=\int_{\mathbb{R}^{N}}\left(\partial_{z_{j_{1}}}\cdot\cdot\cdot\partial_{z_{j_{n}}}u\right)(t,z)(2\pi)^{-N/2}e^{-\frac{|z|^{2}}{2}}dz.$

This completes the proof of our main result (recall that the orthonormality of
the functions $\xi_{1},...,\xi_{N}$ implies that
$I_{1}(\xi_{1}),...,I_{1}(\xi_{N})$ are independent standard Gaussian random
variables).

######  Remark 4.6.

Following the previous construction, we could formally associate equation
(4.6) with

$\displaystyle\begin{split}\partial_{t}u(t,z)&=\sum_{k\geq
1}\alpha_{k}\partial^{2}_{z_{k}}u(t,z)+\sum_{k\geq
1}\left(d_{k}-c_{k}-\alpha_{k}z_{k}\right)\partial_{z_{k}}u(t,z)+\left(\sum_{k\geq
1}c_{k}z_{k}-\gamma\right)u(t,z)\\\ u(0,z)&=\zeta(z),\end{split}$ (4.17)

where now $u(t,z)=u(t,z_{1},z_{2},....)$ is a function of infinitely many
variables, and obtain the identity

$\displaystyle\Phi(t)=u\left(t,I_{1}(\xi_{1}),I_{1}(\xi_{2}),...\right).$

From this point of view, the use of the projection operator $\Pi_{N}$ is
needed for reducing (4.17) to a standard partial differential equation with a
finite number of spatial variables.

## References

  * [1] E. Allen, _Modelling with Itô Stochastic Differential Equations_ , Springer-Verlag, London, 2007.
  * [2] V. I. Bogachev, _Gaussian Measures_ , American Mathematical Society, Providence, 1998.
  * [3] M. J. del Razo, D. Frömberg, A. V. Straube, C. Schütte, F. Höfling, and S. Winkelmann. A probabilistic framework for particle-based reaction-diffusion dynamics using classical Fock space representations, arXiv:2109.13616, 2021.
  * [4] C. Gardiner, _Stochastic Methods: A Handbook for the Natural and Social Sciences - IV edition_ , Springer Series in Synergetics, Springer, Berlin Heidelberg, 2009.
  * [5] D. T. Gillespie, A rigorous derivation of the chemical master equation, _Physica A: Statistical Mechanics and its Applications_ , 188 (1992) pp. 404-425.
  * [6] S. Janson, _Gaussian Hilbert spaces_ , Cambridge Tracts in Mathematics, 129. Cambridge University Press, Cambridge, 1997.
  * [7] I. Karatzas and S. E. Shreve, _Brownian motion and stochastic calculus_ , Springer-Verlag, New York, 1991.
  * [8] P. Lecca, I. Laurenzi, I. and F. Jordan, _Deterministic Versus Stochastic Modelling in Biochemistry and Systems Biology_ , Oxford: Woodhead Publishing, 2013.
  * [9] R. Mikulevicius and B. L. Rozovskii, On unbiased stochastic Navier-Stokes equations, _Probability Theory and Related Fields_ 154 (2012) pp. 787-834.
  * [10] D. Nualart, _Malliavin calculus and Related Topics - II Edition_ , Springer, New York, 2006\.
  * [11] N. G. Van Kampen, _Stochastic Processes in Physics and Chemistry_ , Elsevier, Amsterdam, 1992.