2.1. Itô diffusion process

This section presents a class of SDEs which is related to the representation of our extended RM equation on a domain D in $\mathbb{R}^n$ . Let $T>0$ be an arbitrary fixed time horizon and D a domain in $\mathbb{R}^n$ , i.e., an open connected subset of $\mathbb{R}^n$ . Consider continuous functions $b\,:\,D\to\mathbb{R}^{n\times1}$ and $\sigma\,:\,D\to\mathbb{R}^{n\times m}$ , where $m\geq1$ . We then introduce the following second-order differential operator acting on $C^2(D)$ :

(1) \begin{align}{\mathcal{A}}f(x) &\,:\!=\, b(x)^{\top}\nabla_xf(x) + \frac{1}{2}{\textrm{tr}}[\sigma\sigma^{\top}(x)\nabla_x^2f(x)],\qquad f\in C^2(D), x\in D,\end{align}

where $\nabla_x=(\partial_{x_1},\ldots,\partial_{x_n})^{\top}$ denotes the gradient operator, $\nabla_x^2$ is the corresponding Hessian matrix, and tr denotes the trace operator. The operator ${\mathcal{A}}$ acting on $C^2(D)$ is in fact the infinitesimal generator of the (Markovian) solution of the SDE given, for $t\in[0,T]$ , by

(2) \begin{align}dX_s^{t,x} &= b(X_s^{t,x})ds + \sigma(X_s^{t,x})dW_s, \qquad s\in[t,T],\nonumber\\X_t^{t,x}&=x\in D,\end{align}

where $W=(W_t)_{t\in[0,T]}$ is an m-dimensional Brownian motion on the filtered probability space $(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ with the filtration $\mathbb{F}=(\mathcal{F}_t)_{t\in[0,T]}$ satisfying the usual conditions.

We impose the following assumptions on the coefficients of the SDE (2) to guarantee the existence a unique strong solution of the equation:

1. (A _b,σ) Either (i) $b\,:\,D\to\mathbb{R}^{n\times1}$ and $\sigma\,:\,D\to\mathbb{R}^{n\times m}$ are locally Lipschitz continuous; or (ii) for $m=n=1$ , $b\,:\,D\to\mathbb{R}$ is locally Lipschitz continuous and $\sigma\,:\,D\to\mathbb{R}$ is Hölder continuous with exponent $\gamma\in[\frac{1}{2},1)$ .

2. (A _X ) For all $(t,x)\in[0,T]\times D$ , the solution $X^{t,x}=(X_s^{t,x})_{s\in[t,T]}$ of (2) neither explodes nor leaves D before T, i.e.,

   \begin{equation*}\mathbb{P}\left(\sup_{s\in[0,T]}|X_s^{t,x}|<+\infty\right)=1\end{equation*}
   and
   \begin{equation*}\mathbb{P}(X_s^{t,x}\in D, \forall s\in[t,T])=1.\end{equation*}

By Theorem 5.2.5 of Karatzas and Shreve [Reference Karatzas and Shreve12], the assumption (A $_{b,\sigma}\!$ ) (i) implies that the SDE (2) admits a unique strong solution $X^{t,x}$ up to a possibly finite random explosion time. Then, the assumption ${({\textbf{A}}_{X})}$ yields that this explosion time should be greater than T, $\mathbb{P}$ -almost surely, and hence $X^{t,x}$ is well-defined on [t, T]. For the case with $m=n=1$ (i.e., the SDE (2) is a one-dimensional equation), by Proposition 5.2.13 in [Reference Karatzas and Shreve12], the assumption (A $_{b,\sigma}\!$ ) (ii) gives that the SDE (2) admits a unique strong solution $X^{t,x}$ up to a possibly finite random explosion time. Therefore, the assumption ${({\textbf{A}}_{X})}$ yields that this explosion time is greater than T, almost surely, and hence $X^{t,x}$ is well-defined on [t, T].

2.2. Extended replicator–mutator equations on $[0,\textbf{\textit{T}}]\times \textbf{\textit{D}}$

We introduce a class of extended RM equations as a non-local Cauchy problem on $[0,T]\times D$ . More precisely, let ${\mathcal{A}}^*$ be the adjoint operator of ${\mathcal{A}}$ defined by (1). Then the extended RM equation considered in this paper is given by

(3) \begin{align}\left\{ \begin{array}{ll} \displaystyle \partial_tu(t,x)={\mathcal{A}}^*u(t,x)+\left(g(x)-\int_{D}g(y)u(t,y)dy\right)u(t,x),\\[16pt] \displaystyle u(0,x)=u_0(x),\quad (t,x)\in[0,T]\times D, \end{array}\right.\end{align}

where $u_0\,:\,D\to\mathbb{R}_+$ is a probability density function, i.e., $\int_D u_0(x)dx=1$ , and $g\,:\,D\to\mathbb{R}$ is referred to as the fitness function (see, e.g. Kimura [Reference Kimura13] and Tsimring et al. [Reference Tsimring, Levine and Kessler25]). The condition satisfied by the fitness function g will be imposed later.

We next give the definition of a weak solution of the above extended RM equation (3). To this end, let ${\mathcal{D}}\,:\!=\,C_b^{2}(D)$ be the space of real-valued functions with bounded and continuous partial derivatives up to order 2. Then a function $u\,:\,[0,T]\times D\to\mathbb{R}$ is called a weak solution of the RM equation (3) if it satisfies the following variational form:

(4) \begin{align}\langle u(t),f\rangle &=\langle u_0,f\rangle +\int_0^t\langle u(s),({\mathcal{A}}+g)f\rangle ds\quad-\int_0^t\langle u(s),g\rangle \langle u(s),f\rangle ds,\qquad f\in{\mathcal{D}},\end{align}

where the integral $\langle u(t),f\rangle\,:\!=\,\int_D u(t,x)f(x)dx$ . One of the objectives of this paper is to establish a closed-form probabilistic representation of the weak solution to the RM equation (3) in the above distributional sense. We provide sufficient conditions on the fitness function g and the initial function $u_0$ under which the solution of (4) admits an explicit form which can be expressed in terms of the transition density function of some Itô diffusion processes. The main strategy for achieving this aim is to provide a probabilistic representation of the solution of (3) by using the theory of FPK equations and a martingale extraction approach.

The second objective of this paper is to establish the propagation of chaos of the RM equation (3) for a relatively large class of fitness functions g. It is in fact equivalent to the convergence of the empirical measure of a (homogeneous) particle system to the solution of the FPK equation related to (4) with respect to a Wasserstein-like distance adapted to our probabilistic framework. We also establish an explicit rate of convergence in the propagation of chaos.

3.1. Probabilistic solution of FPK equation

Denote by ${\mathcal{P}}(D)$ the set of probability measures on ${\mathcal{B}}_{D}$ (i.e., the $\sigma$ -algebra generated by the open subsets of D). For $p\geq1$ , let ${\mathcal{P}}_p(D)\subset{\mathcal{P}}(D)$ be the set of probability measures on ${\mathcal{B}}_{D}$ with finite pth-order moment. We introduce the FPK equation associated with (4), which is given by

(5) \begin{align}\langle\mu_t,f\rangle = \langle\rho_0,f\rangle+\int_0^t \langle\mu_s,({\mathcal{A}}+g)f\rangle ds - \int_0^t\langle\mu_s,f\rangle\langle\mu_s,g\rangle ds, \qquad f\in{\mathcal{D}},\end{align}

where the initial data $\rho_0\in{\mathcal{P}}(D)$ and the integral $\langle\mu_t,f\rangle\,:\!=\,\int_D fd\mu_t$ for $t\in[0,T]$ .

The following lemma provides a probabilistic representation of the ${\mathcal{P}}(D)$ -valued solution $\mu=(\mu_t)_{t\in[0,T]}$ for the FPK equation.

Lemma 1. Let the assumptions ( $\textbf{\textrm{A}}_{b,\sigma}\!$ ) and ${({\textbf{A}}_{X})}$ hold. Assume that the fitness function $g\,:\,D\to\mathbb{R}$ satisfies

1. (A _g ) $(t,x)\to\mathbb{E}[g(X_t^x)\exp(\int_0^tg(X_s^{x})ds)]$ is in $L^1([0,T]\times D;\,dt\otimes\rho_0(dx))$ , where $X_t^{x}\,:\!=\,X_{t}^{0,x}$ for $t\in[0,T]$ .

Let us make the definition

(6) \begin{align}\mu_t= \frac{\int_{D}\mathbb{E}\left[\delta_{X_t^y}\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)}{\int_{D}\mathbb{E}\left[\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)},\quad {\textrm{on}}\ {\mathcal{B}}_{D},\end{align}

where $\delta$ denotes the Dirac delta measure. Then $\mu=(\mu_t)_{t\in[0,T]}$ as defined by (6) is a ${\mathcal{P}}(D)$ -valued solution of the FPK equation (5).

Proof. It is obvious that $\mu_t$ as defined by (6) is a probability measure on ${\mathcal{B}}_D$ for any $t\in[0,T]$ , where $\mu_0=\rho_0$ . We define

\begin{align*}h_t\,:\!=\,\int_D\mathbb{E}\left[\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy),\qquad t\in[0,T],\end{align*}

which is well-defined due to the assumption ( $\textbf{\textrm{A}}_{g}$ ). By virtue of the representation (6), we have that

\begin{align*}h_s\langle\mu_s,g\rangle=\int_{D}\mathbb{E}\left[g(X_s^y)\exp\left(\int_0^sg(X_r^y)dr\right)\right]\rho_0(dy),\qquad s\in[0,T].\end{align*}

Integrating on both sides of the above equality with respect to s from 0 to $t\in(0,T]$ , it follows from Fubini’s theorem and ( $\textbf{\textrm{A}}_{g}$ ) that

(7) \begin{align}\int_0^t h_s\langle\mu_s,g\rangle ds&=\int_{D}\mathbb{E}\left[\int_0^tg(X_s^y)\exp\left(\int_0^sg(X_r^y)dr\right)ds\right]\rho_0(dy)\nonumber\\[5pt] &=\int_D\mathbb{E}\left[\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)-1\\[5pt] &=h_t-1.\nonumber\end{align}

This yields the ODE given by $h^{\prime}_t=h_t\langle\mu_t,g\rangle$ with $h_0=1$ . It then holds that

(8) \begin{align}\exp\left(\int_0^t\langle\mu_s,g\rangle ds\right)=\int_{D}\mathbb{E}\left[\exp\left(\int_0^t g(X_s^{y})ds\right)\right]\rho_0(dy),\quad t\in[0,T].\end{align}

It follows from (6) and (8) that

(9) \begin{align}\mu_t &= \int_{D}\mathbb{E}\left[\delta_{X_t^{y}}\exp\left(\int_0^t(g(X_s^{y})-\langle\mu_s,g\rangle) ds\right)\right]\rho_0(dy),\quad {\textrm{on}}\ {\mathcal{B}}_{D}.\end{align}

Then, for any test function $f\in{\mathcal{D}}$ , by applying integration by parts, it follows that

\begin{align*}&f\left(X_t^y\right)\exp\left(\int_0^t(g(X_s^{y})-\langle\mu_s,g\rangle) ds\right)=f(y) \nonumber\\[3pt] &\qquad\qquad+ \int_0^t\exp\left(\int_0^s(g(X_r^{y})-\langle\mu_r,g\rangle) dr\right){\mathcal{A}}f(X_s^y)ds\nonumber\\[3pt] &\qquad\qquad+\int_0^t\exp\left(\int_0^s(g(X_r^{y})-\langle\mu_r,g\rangle) dr\right)\nabla_xf(X_s^y)^{\top}\sigma(X_s^y)dW_s\nonumber\\[3pt] &\qquad\qquad+\int_0^t\exp\left(\int_0^s(g(X_r^{y})-\langle\mu_r,g\rangle) dr\right)f(X_s^y)(g(X_s^{y})-\langle\mu_s,g\rangle)ds.\end{align*}

This yields that

\begin{align*}&\mathbb{E}\left[f\left(X_t^y\right)\exp\left(\int_0^t(g(X_s^{y})-\langle\mu_s,g\rangle) ds\right)\right]= f(y)\nonumber\\[3pt] &\qquad+\int_0^t \mathbb{E}\left[\exp\left(\int_0^s(g(X_r^{y})-\langle\mu_r,g\rangle) dr\right){\mathcal{A}}f(X_s^y)\right]ds\nonumber\\[3pt] &\qquad+\int_0^t \mathbb{E}\left[\exp\left(\int_0^s(g(X_r^{y})-\langle\mu_r,g\rangle) dr\right)f(X_s^y)(g(X_s^{y})-\langle\mu_s,g\rangle)\right]ds.\end{align*}

Integrate on both sides of the above equality with respect to $\rho_0(dy)$ . Then, by virtue of (9), we arrive at

\begin{align*}\langle\mu_t,f\rangle&=\langle\rho_0,f\rangle+\int_0^t\langle\mu_s,{\mathcal{A}}f\rangle ds+\int_0^t \langle\mu_s,gf-\langle\mu_s,g\rangle f\rangle ds.\end{align*}

This yields the FPK equation (5).

3.2. Closed-form representation of $\mu$ given by (6)

We study a closed-form representation of the ${\mathcal{P}}(D)$ -valued function $\mu$ given by (6) in Lemma 1 by applying a martingale approach. We find that the obtained explicit form of $\mu$ can be in fact represented in terms of the transition density function of a class of Itô diffusion processes. Moreover, we provide examples in which the expressions for $\mu$ and its density admit closed-form representations.

By (6), the representation of $\mu$ can be reduced to identifying the following function: for $f\in C_b(D)$ , i.e., the space of bounded continuous functions on D, define

(10) \begin{align}\Phi^{f}(t,x)\,:\!=\, \mathbb{E}\left[f(X_t^{x})\exp\left(\int_0^t g(X_s^{x})ds\right)\right],\quad (t,x)\in[0,T]\times D,\end{align}

where the process $X^{x}=(X_t^{0,x})_{t\in[0,T]}$ is the unique strong solution of the SDE (2). We next apply a martingale approach to give a closed-form representation of $\mu$ when the coefficients of the equation $(b,\sigma,g)$ satisfy an additional constraint.

Lemma 2. Let the assumptions $(\textbf{\textrm{A}}_{b,\sigma}\!)$ , ${({\textbf{A}}_{X})}$ , and $(\textbf{\textrm{A}}_{g})$ hold. In addition, assume the following:

1. (i) The fitness function $g\in C^2(D)$ and satisfies that, for all $x\in D$ ,

   (11) \begin{align} {\mathcal{A}}g(x)=C_1\in\mathbb{R}, \text{and} \nabla_xg(x)^{\top}\sigma(x)=C_2\in\mathbb{R}^{1\times m}, \end{align}
   where the operator ${\mathcal{A}}$ acting on $C^2(D)$ is given by (1).

2. (ii) For $y\in D$ , let $\overline{X}^y=\big(\overline{X}_t^y\big)_{t\in[0,T]}$ be the strong solution of the following SDE:

   (12) \begin{align} \overline{X}_t^y = y+ \int_0^t\overline{b}\left(s,\overline{X}_s^y\right)ds + \int_0^t\sigma\left(\overline{X}_s^y\right)dW_s,\end{align}
   where the time-dependent function $\overline{b}(t,x)\,:\!=\,b(x)-t\sigma(x)C_2^{\top}$ for $(t,x)\in[0,T]\times D$ .

Then, for any $T>0$ , the solution $\mu$ given by (6) of the FPK equation (5) admits the following form:

(13) \begin{align}\mu_t(dx) =\frac{e^{tg(x)}\int_{D}\mathbb{P}\left(\overline{X}_t^y\in dx\right)\rho_0(dy)}{\int_{D}e^{tg(z)}\int_{D}\mathbb{P}\left(\overline{X}_t^y\in dz\right)\rho_0(dy)},\quad t\in[0,T].\end{align}

Moreover, if $\rho_0(dx)=u_0(x)dx$ and the density function $\overline{p}(t,y;\,x)$ of $\overline{X}_t^y$ for $t\in[0,T]$ exists, i.e.

\begin{equation*}\mathbb{P}\left(\overline{X}_t^y\in dx\right)=\overline{p}(t,y;\,x)dx,\end{equation*}

then, for all $t\in[0,T]$ , the probability measure $\mu_t(dx)=u(t,x)dx$ . Here u(t, x) satisfies the extended RM equation (4), and it has the following representation:

(14) \begin{align}u(t,x)=\frac{e^{tg(x)}\int_{D}\overline{p}(t,y;\,x)u_0(y)dy}{\int_{D}e^{tg(z)}\left(\int_{D}\overline{p}(t,y;\,z)u_0(y)dy\right)dz},\quad (t,x)\in[0,T]\times D.\end{align}

Proof. Integration by parts yields that, for $(t,y)\in[0,T]\times D$ ,

\begin{align*}\int_0^t g(X_s^y)ds &= tg(X_t^y)-\int_0^t s{\mathcal{A}}g(X_s^y)ds-\int_0^t s \nabla_xg(X_s^y)^{\top}\sigma(X_s^y)dW_s\nonumber\\&=tg\left(X_t^y\right)+\int_0^t\left\{\frac{1}{2}s^2\nabla_xg\left(X_s^y\right)^{\top}\left(\sigma\sigma^{\top}\right)\left(X_s^y\right)\nabla_xg\left(X_s^y\right)-s{\mathcal{A}}g\left(X_s^y\right)\right\}ds\nonumber\\& -\int_0^t s \nabla_xg(X_s^y)^{\top}\sigma(X_s^y)dW_s-\frac{1}{2}\int_0^ts^2\nabla_xg(X_s^y)^{\top}\left(\sigma\sigma^{\top}\right)(X_s^y)\nabla_xg(X_s^y)ds.\end{align*}

Note that, by the condition (11) in (i), we have that ${\mathcal{A}}g(x)\equiv C_1\in\mathbb{R}$ and $\nabla_xg(x)^{\top}\sigma(x)=C_2\in\mathbb{R}^{1\times m}$ . Then, by verifying Novikov’s condition, we can define a $(\mathbb{P},\mathbb{F})$ -martingale as follows:

\begin{align*}N_t\,:\!=\,\exp\left(-\int_0^t s C_2dW_s-\frac{C_2C_2^{\top}}{6}t^3\right),\quad t\in[0,T].\end{align*}

Therefore, it follows from (10) that

(15) \begin{align}\Phi^f(t,y)=\mathbb{E}\left[N_tf(X_t^{y})\exp\left(tg(X_t^y)+\frac{C_2C_2^{\top}}{6}t^3-\frac{C_1}{2}t^2\right)\right], \qquad (t,y)\in[0,T]\times D.\end{align}

Let us define

\begin{equation*}\frac{d\mathbb{Q}}{d\mathbb{P}}|_{\mathcal{F}_t}=N_t,\end{equation*}

for $t\in[0,T]$ . Then

\begin{equation*}B_t\,:\!=\,W_t+\frac{C_2^{\top}t^2}{2}\end{equation*}

for $t\in[0,T]$ is an m-dimensional Brownian motion under the probability measure $\mathbb{Q}$ . Moreover, under the new probability measure $\mathbb{Q}$ , the dynamics of the process $X^y=(X_t^y)_{t\in[0,T]}$ is given by

\begin{align*}X_t^y = y+ \int_0^t\{b(X_s^y)-s\sigma(X_s^y)C_2^{\top}\}ds + \int_0^t\sigma(X_s^y)dB_s,\quad y\in D.\end{align*}

By applying (15), we arrive at

\begin{align*}\Phi^f(t,y)&=\exp\left(\frac{C_2C_2^{\top}}{6}t^3-\frac{C_1}{2}t^2\right)\mathbb{E}^{\mathbb{Q}}\left[f(X_t^{y})e^{tg(X_t^y)}\right]\nonumber\\&=\exp\left(\frac{C_2C_2^{\top}}{6}t^3-\frac{C_1}{2}t^2\right)\int_{D}f(x)e^{tg(x)}\mathbb{Q}(X_t^y\in dx),\end{align*}

where $\mathbb{E}^{\mathbb{Q}}$ denotes the expectation operator under $\mathbb{Q}$ . Consequently, the solution representation (6) yields that

(16) \begin{align}\int_{D} f(x)\mu_t(dx) = \frac{\int_{D}f(x)e^{tg(x)}\int_{D}\mathbb{Q}(X_t^y\in dx)\rho_0(dy)}{\int_{D}e^{tg(x)}\int_{D}\mathbb{Q}(X_t^y\in dx)\rho_0(dy)},\quad t\in[0,T].\end{align}

Note that $\mathbb{Q}(X_t^y\in dx)=\mathbb{P}(\overline{X}_t^y\in dx)$ . Therefore, the desired results (13) and (14) follow from the equality (16).

Lemma 2 works well for the case where the drift and volatility functions $(b,\sigma)$ are constant matrices and the fitness function g is a linear map. This is considered in the following example.

Example 1. Consider $m=n$ , $b(x)\equiv b\in\mathbb{R}^{n}$ , and $\sigma(x)\equiv\sigma\in\mathbb{R}^{n\times n}$ . Let $\sigma\in\mathbb{R}^{n\times n}$ be invertible, with its inverse matrix given by $\sigma^{-1}$ . In this example, the SDE (2) is reduced to a drift-Brownian motion described as

\begin{align*}X_t^{x}=x+bt+\sigma W_t,\quad t\in[0,T].\end{align*}

Accordingly, the domain is $D=\mathbb{R}^n$ . For any $C_2\in\mathbb{R}^{1\times n}$ , we consider the following fitness function:

(17) \begin{align}g(x) = C_2\sigma^{-1}x,\quad x\in\mathbb{R}^{n\times1}.\end{align}

Then a direct calculation shows that if the initial data $\rho_0$ of the FPK equation (5) satisfies

\begin{equation*}\int_0^T\int_D(|x|+1)e^{tC_2\sigma^{-1}x}\rho_0(dx)dt<\infty,\end{equation*}

then the assumption $(\textbf{\textrm{A}}_{g})$ holds. Thus, we have from (17) that

\begin{equation*}{\mathcal{A}}g\equiv C_1\,:\!=\,C_2\sigma^{-1}b\in\mathbb{R}.\end{equation*}

Using the condition (ii) of Lemma 2, the process $\overline{X}^y=\big(\overline{X}_t^y\big)_{t\in[0,T]}$ is given by

\begin{align*}\overline{X}_t^y = y+ bt -\frac{\sigma C_2^{\top}}{2}t^2 +\sigma W_t,\quad (t,y)\in[0,T]\times\mathbb{R}^n.\end{align*}

Therefore, under the original probability measure $\mathbb{P}$ , we have that

\begin{equation*}\overline{X}_t^y\sim N\left(y+ bt -\frac{\sigma C_2^{\top}}{2}t^2,\sigma\sigma^{\top}t\right),\quad t\in(0,T].\end{equation*}

This implies that the density function $\overline{p}(t,y;\,x)$ of $\overline{X}_t^y$ has the following closed-form representation: $\overline{p}(0,y;\,x)=\delta_{x-y}$ , and for $(t,x,y)\in(0,T]\times\mathbb{R}^{2n}$ ,

(18) \begin{align}\overline{p}(t,y;\,x)&=\frac{1}{(2\pi)^{\frac{n}{2}}}\frac{1}{\sqrt{{\textrm{det}}(\sigma\sigma^{\top})t}}\\[5pt] & \times\exp\left\{-\frac{1}{2t}\left(x-y-bt +\frac{\sigma C_2^{\top}}{2}t^2\right)^{\top}(\sigma\sigma^{\top})^{-1}\left(x-y-bt +\frac{\sigma C_2^{\top}}{2}t^2\right)\right\}.\nonumber\end{align}

It follows from (14) in Lemma 2 that the density function of $\mu_t$ is given by

(19) \begin{align}u(t,x)=\frac{e^{C_2\sigma^{-1}xt}\int_{\mathbb{R}^n}\overline{p}(t,y;\,x)u_0(y)dy}{\int_{\mathbb{R}^n}e^{C_2\sigma^{-1}zt}\left(\int_{\mathbb{R}^n}\overline{p}(t,y;\,z)u_0(y)dy\right)dz},\quad (t,x)\in[0,T]\times\mathbb{R}^n.\end{align}

By Remark 1, it follows that u(t, x) is a weak solution of the RM equation (3). Note that the generator ${\mathcal{A}}$ of $X^x=(X_t^{x})_{t\in[0,T]}$ in this example is given by

\begin{equation*}{\mathcal{A}}=b^{\top}\nabla_x+\frac{1}{2}{\textrm{tr}}\left[\sigma\sigma^{\top}\nabla_x^2\right].\end{equation*}

Then the adjoint operator of ${\mathcal{A}}$ is given by

\begin{equation*}{\mathcal{A}}^*=\frac{1}{2}{\textrm{tr}}\left[\sigma\sigma^{\top}\nabla_x^2\right]-b^{\top}\nabla_x.\end{equation*}

Let $\rho_0(dx)=u_0(x)dx$ ; i.e., the initial $\rho_0$ admits a density function $u_0(x)$ for $x\in D$ . In view of (3), we in fact establish the probabilistic representation of the weak solution to the following RM equation with the fitness space $D=\mathbb{R}^n$ :

(20) \begin{align}\left\{ \begin{array}{ll} \displaystyle \partial_tu(t,x)=\frac{1}{2}{\textrm{tr}}\left[\sigma\sigma^{\top}\nabla_x^2\right]u(t,x)-b^{\top}\nabla_xu(t,x)\\[4pt] \displaystyle \qquad\qquad\quad+\left(g(x)-\int_{\mathbb{R}^n}g(y)u(t,y)dy\right)u(t,x),\\[10pt] \displaystyle u(0,x)=u_0(x), \qquad x\in\mathbb{R}^n. \end{array}\right.\end{align}

A special case of Example 1 is the one-dimensional RM equation with $b\equiv0$ , $\sigma=\sqrt{2}$ and the fitness function $g(x)=x$ for $x\in\mathbb{R}$ . In this case, (20) becomes that

(21) \begin{align}\left\{ \begin{array}{ll} \displaystyle \partial_tu(t,x)=\Delta u(t,x)+\left(x-\int_{\mathbb{R}}yu(t,y)dy\right)u(t,x),\\[4pt] \displaystyle u(0,x)=u_0(x), x\in\mathbb{R}. \end{array}\right.\end{align}

The explicit solution to the RM equation (21) has been studied by Alfaro and Carles [Reference Alfaro and Carles1]. We next verify that the probabilistic solution (19) of the RM equation (21) coincides with the solution form given by (2.3) of Theorem 2.3 in [Reference Alfaro and Carles1], which uses a tricky transform of the solution based on the Avron–Herbst formula. In fact, in view of Lemma 2, from (17) and (18), it results in $C_2=\sqrt{2}$ , and hence

(22) \begin{align}\overline{p}(t,y;\,x)=\frac{1}{2\sqrt{\pi t}}\exp\left(-\frac{(x-y+t^2)^2}{4t}\right),\quad (t,x,y)\in[0,T]\times\mathbb{R}^2.\end{align}

This implies that the probabilistic solution (19) is given by

(23) \begin{align}u(t,x)=\frac{e^{tx}\int_{\mathbb{R}}\exp\left(-\frac{(x-y+t^2)^2}{4t}\right)u_0(y)dy}{\int_{\mathbb{R}}e^{tz}\left(\int_{\mathbb{R}}\exp\left(-\frac{(z-y+t^2)^2}{4t}\right)u_0(y)dy\right)dz},\quad (t,x,y)\in[0,T]\times\mathbb{R}^2.\end{align}

By (6) and the assumption $(\textbf{\textrm{A}}_{g})$ , it is not difficult to verify that

\begin{align*} \int_{\mathbb{R}}xu(t,x)dx=\int_{\mathbb R}x\mu_t(dx)= \frac{\int_{\mathbb{R}}\mathbb{E}\left[X_t^y\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)}{\int_{\mathbb{R}}\mathbb{E}\left[\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)}<\infty, \quad t\in[0,T].\end{align*}

A direct calculation gives that

\begin{align*}&\int_{\mathbb{R}}e^{tz}\left(\int_{\mathbb{R}}\exp\left(-\frac{(z-y+t^2)^2}{4t}\right)u_0(y)dy\right)dz\nonumber\\[4pt] &\qquad=\int_{\mathbb{R}^2}\exp\left(\frac{-(z-y-t^2)^2}{4t}+ty\right)u_0(y)dzdy=\sqrt{4\pi t}\int_{\mathbb{R}}e^{ty}u_0(y)dy.\end{align*}

As a consequence,

(24) \begin{align}u(t,x)&=\frac{e^{tx}\int_{\mathbb{R}}\exp\left(-\frac{(x-y+t^2)^2}{4t}\right)u_0(y)dy}{\int_{\mathbb{R}}e^{tz}\left(\int_{\mathbb{R}}\exp\left(-\frac{(z-y+t^2)^2}{4t}\right)u_0(y)dy\right)dz}\nonumber\\[8pt] &=\frac{\frac{e^{tx}}{\sqrt{4\pi t}}\int_{\mathbb{R}}\exp\left(-\frac{(x-y+t^2)^2}{4t}\right)u_0(y)dy}{\int_{\mathbb{R}}e^{ty}u_0(y)dy},\quad (t,x)\in[0,T]\times\mathbb{R},\end{align}

which is the solution representation given by (2.3) of Theorem 2.3 in [Reference Alfaro and Carles1].

However, the condition (11) on the coefficients $(b,\sigma,g)$ in Lemma 2 is very restrictive. In order to bypass this constraint on $(b,\sigma,g)$ , we next apply a martingale extraction technique to identify the function $\Phi^f$ defined by (10). More precisely, recall the operator ${\mathcal{A}}$ acting on $C^2(D)$ , which is defined by (1). Let $(\lambda,\phi)$ be a (nonzero) eigenpair (if it exists) of the following characteristic equation:

(25) \begin{align}({\mathcal{A}}+ g(x))\phi(x) = -\lambda \phi(x),\quad x\in D.\end{align}

Then, it follows from (25) and Itô’s formula that

(26) \begin{align}M_t^y \,:\!=\, \frac{1}{\phi(y)}\exp\left(\lambda t+\int_0^tg(X_s^y)ds\right)\phi(X_t^y),\qquad t\in[0,T],\end{align}

is a local $(\mathbb{P},\mathbb{F})$ -martingale with $M_0^y=1$ . Then, by (10), it holds that

(27) \begin{align}\Phi^{f}(t,y)= e^{-\lambda t}\phi(y)\mathbb{E}\left[M_t^y\frac{f(X_t^y)}{\phi(X_t^y)}\right],\quad (t,y)\in[0,T]\times D.\end{align}

A key observation is that the expectation in (27) does not depend on the path of the process $X^y$ . This yields the following result.

Theorem 1. Let the assumptions $(\textbf{\textrm{A}}_{b,\sigma}\!)$ , ${({\textbf{A}}_{X})}$ , and $(\textbf{\textrm{A}}_{g})$ hold. In addition, for all $y\in D$ , assume the following:

1. (i) There exists an eigenpair $(\lambda,\phi)$ of the characteristic equation (25) such that the process $M^y=(M_t^y)_{t\in[0,T]}$ defined by (26) is a $(\mathbb{P},\mathbb{F})$ -martingale.

2. (ii) Let $\overline{X}^y=\big(\overline{X}_t^y\big)_{t\in[0,T]}$ be the strong solution of the following SDE:

   (28) \begin{align} \overline{X}_t^y = y+ \int_0^t\left(b\left(\overline{X}_s^y\right)+\sigma\left(\overline{X}_s^y\right)\overline{\sigma}\left(\overline{X}_s^y\right)^{\top}\right)ds + \int_0^t\sigma\left(\overline{X}_s^y\right)dW_s,\quad t\in[0,T],\end{align}
   where the function
   \begin{equation*}\overline{\sigma}(x)\,:\!=\,\frac{\nabla_x\phi(x)^{\top}}{\phi(x)}\sigma(x)\end{equation*}
   for $x\in D$ .

Then, for any $T>0$ , the solution $\mu$ given by (6) of the FPK equation (5) admits the following form:

(29) \begin{align}\mu_t(dx) =\frac{\frac{1}{\phi(x)}\int_{D}\phi(y)\mathbb{P}\left(\overline{X}_t^y\in dx\right)\rho_0(dy)}{\int_{D}\frac{\int_{D}\phi(y)\mathbb{P}(\overline{X}_t^y\in dz)\rho_0(dy)}{\phi(z)}},\quad t\in(0,T].\end{align}

Moreover, if $\rho_0(dx)=u_0(x)dx$ and the density function $\overline{p}(t,y;\,x)$ of $\overline{X}_t^y$ for $t\in(0,T]$ exists, i.e., $\mathbb{P}\big(\overline{X}_t^y\in dx\big)=\overline{p}(t,y;\,x)dx$ , then, for all $t\in[0,T]$ , the probability measure $\mu_t(dx)=u(t,x)dx$ . Here, u(t, x) satisfies the extended RM equation (4), and it admits the following representation:

(30) \begin{align}u(t,x)=\frac{\frac{1}{\phi(x)}\int_{D}\phi(y)\overline{p}(t,y;\,x)u_0(y)dy}{\int_{D}\frac{\int_{D}\phi(y)\overline{p}(t,y;\,z)u_0(y)dy}{\phi(z)}dz},\quad (t,x)\in[0,T]\times D.\end{align}

Proof. Using (26) and (25), the dynamics of the martingale $M^y$ defined by (26) is given by

(31) \begin{align}dM_t^y &=M_t^y\frac{\nabla_x\phi\left(X_t^y\right)^{\top}\sigma\left(X_t^y\right)}{\phi\left(X_t^y\right)}dW_t=M_t^y\overline{\sigma}\left(X_t^y\right) dW_t,\quad M_0^y=1.\end{align}

For fixed $y\in D$ , define

\begin{equation*}\frac{d\mathbb{Q}}{d\mathbb{P}}|_{\mathcal{F}_t}=M_t^y\end{equation*}

for $t\in[0,T]$ . Then, from Girsanov’s theorem, it follows that $B_t\,:\!=\,W_t-\int_0^t\overline{\sigma}(X_s^y)^{\top}ds$ for $t\in[0,T]$ is an m-dimensional Brownian motion under the new probability measure $\mathbb{Q}$ . Thus, under $\mathbb{Q}$ , the dynamics of the process $X^y=(X_t^y)_{t\in[0,T]}$ is given by

\begin{align*}X_t^y = y+ \int_0^t\left(b\left(X_s^y\right)+\sigma\left(X_s^y\right)\overline{\sigma}\left(X_s^y\right)^{\top}\right)ds + \int_0^t\sigma\left(X_s^y\right)dB_s,\quad (t,y)\in[0,T]\times D.\end{align*}

By virtue of (27), we have that

\begin{align*}\Phi^{f}(t,y)&= e^{-\lambda t}\phi(y)\mathbb{E}^{\mathbb{Q}}\left[\frac{f(X_t^y)}{\phi(X_t^y)}\right]\nonumber\\[6pt] &=e^{-\lambda t}\phi(y)\int_{D}\frac{f(x)}{\phi(x)}\mathbb{Q}(X_t^y\in dx), (t,y)\in[0,T]\times D.\end{align*}

Apply Lemma 1, we obtain that

(32) \begin{align}\int_{D}f(x)\mu(dx) = \frac{\int_{D}\frac{f(x)}{\phi(x)}\int_{D}\phi(y)\mathbb{Q}(X_t^y\in dx)\rho_0(dy)}{\int_{D}\frac{\int_{D}\phi(y)\mathbb{Q}(X_t^y\in dx)\rho_0(dy)}{\phi(x)}},\quad \forall f\in{\mathcal{D}}.\end{align}

Note that $\mathbb{Q}(X_t^y\in dx)=\mathbb{P}(\overline{X}_t^y\in dx)$ . The desired result then follows from (32).

To apply Theorem 1, the key point is to identify an eigenpair $(\lambda,\phi)$ of the characteristic equation (25) so as to satisfy the condition (i) of Theorem 1. We next provide examples illustrating how to verify this condition (i) and derive the probability solution (29)–(30) explicitly. We first consider the one-dimensional case, i.e., $n=m=1$ . In this case, let $D=(\alpha,\beta)$ with $-\infty\leq\alpha<\beta\leq+\infty$ . To this end, for a function $V\,:\,D\to\mathbb{R}$ which is locally Hölder continuous, we define the following operator acting on $C^2(D)$ :

(33) \begin{align}{\mathcal{H}}_V \,:\!=\, {\mathcal{A}} + V.\end{align}

We call a positive function $\phi\in C^2(D)$ a positive harmonic function of ${\mathcal{H}}_V$ if ${\mathcal{H}}_V\phi=0$ . A function $\phi\in C^2(D)$ is said to be an invariant function of the semigroup generated by ${\mathcal{H}}_V$ if

\begin{equation*}\mathbb{E}\left[e^{\int_0^tV(X_s^x)ds}\phi\left(X_t^x\right)\right]=\phi(x)\end{equation*}

for $(t,x)\in[0,T]\times D$ . The following result provides a necessary and sufficient condition under which a positive harmonic function of ${\mathcal{H}}_V$ on D is also an invariant function of the semigroup generated by this operator.

Lemma 3. (Theorem 5.1.8 in Pinsky [Reference Pinsky17]) Let the assumptions $(\textbf{\textrm{A}}_{b,\sigma}\!)$ and ${({\textbf{A}}_{X})}$ hold, and let $x_0\in D$ . Then a positive harmonic function of ${\mathcal{H}}_V$ is also an invariant function of the semigroup generated by ${\mathcal{H}}_V$ if and only if the following conditions hold:

(34) \begin{align}&\int_{\alpha}^{x_0}\left[\frac{1}{\phi^2(x)}\exp\left(-\int_{x_0}^x\frac{2b(z)}{\sigma^2(z)}dz\right)\int_{x}^{x_0}\frac{\phi^2(y)}{\sigma^2(y)}\exp\left(\int_{x_0}^y\frac{2b(z)}{\sigma^2(z)}dz\right)dy\right]dx=+\infty,\nonumber\\[6pt] &\int_{x_0}^{\beta}\left[\frac{1}{\phi^2(x)}\exp\left(-\int_{x_0}^x\frac{2b(z)}{\sigma^2(z)}dz\right)\int_{x_0}^{x}\frac{\phi^2(y)}{\sigma^2(y)}\exp\left(\int_{x_0}^y\frac{2b(z)}{\sigma^2(z)}dz\right)dy\right]dx=+\infty.\end{align}

Lemma 3 can be used to verify the martingale property of $M^y=(M_t^y)_{t\in[0,T]}$ defined by (26) in some one-dimensional cases. We next provide examples in the one-dimensional case by applying Lemma 3.

Example 2. (The Vasicek model) Consider $m=n=1$ , $b(x)=\kappa(\theta-x)$ , and $\sigma(x)=\sigma$ for $x\in D\,:\!=\,\mathbb{R}$ , where $\theta,\kappa\in\mathbb{R}$ with $\kappa\neq0$ and $\sigma>0$ . Then the assumptions ( $\textbf{\textrm{A}}_{b,\sigma}\!$ ) and ${({\textbf{A}}_{X})}$ are satisfied. In this case, the process $X^{x}=(X_t^{x})_{t\in[0,T]}$ given by (2) becomes the Vasicek model, i.e.,

(35) \begin{align}dX_t^{x} = \kappa(\theta-X_t^{x})dt + \sigma dW_t,\quad X_0^{x}=x\in\mathbb{R}.\end{align}

The fitness function is $g(x)=-x$ for $x\in\mathbb{R}$ . We next compute $\mathbb{E}\left[\exp\left(-\int_0^tX_s^{x}ds\right)\right]$ for $(t,x)\in[0,T]\times\mathbb{R}$ . In fact, we have that

\begin{align*}\int_0^tX_s^{x}ds&=\frac{1-e^{-\kappa t}}{\kappa}(x-\theta)+\theta t + \frac{\sigma}{\kappa}W_t - \frac{\sigma}{\kappa}\int_0^te^{-\kappa(t-s)}dW_s.\end{align*}

This yields that $\int_0^tX_s^{x}ds\sim N\big({\mu}(t,x),\sigma^2(t)\big)$ for $t\in(0,T]$ , where

\begin{align*}\mu(t,x) &\,:\!=\,\frac{1-e^{-\kappa t}}{\kappa}(x-\theta)+\theta t,\nonumber\\\sigma^2(t)&\,:\!=\,\frac{\sigma^2t}{\kappa^2}+\frac{\sigma^2}{2\kappa^3}(1-e^{-2\kappa t})-\frac{2\sigma^2}{\kappa^3}(1-e^{-\kappa t}).\end{align*}

Therefore

(36) \begin{align}\mathbb{E}\left[\exp\left(-\int_0^tX_s^{x}ds\right)\right]&=e^{-\mu(t,x)+\frac{1}{2}\sigma^2(t)},\quad t\in[0,T].\end{align}

It follows from (36) that there exists a constant $C_T>0$ which is independent of x such that

\begin{align*}&\mathbb{E}\left[\left|X_t^{x}\right|\exp\left(-\int_0^tX_s^{x}ds\right)\right]\leq xe^{-\kappa t}\mathbb{E}\left[\exp\left(-\int_0^tX_s^{x}ds\right)\right]\\[6pt] & +e^{-\mu(t,x)}\mathbb{E}\left[\left|\theta(1-e^{-\kappa t})+\sigma\int_0^te^{-\kappa(t-s)}dW_s\right|\exp\left(-\frac{\sigma}{\kappa}W_t + \frac{\sigma}{\kappa}\int_0^te^{-\kappa(t-s)}dW_s\right)\right]\\[6pt] & \leq C_T(x+1)e^{-\frac{1-e^{-\kappa t}}{\kappa}x}.\end{align*}

This implies that if the initial data $\rho_0$ of the FPK equation (5) satisfies

(37) \begin{align} \int_0^T\int_{\mathbb{R}}(x+1)e^{-\frac{1-e^{-\kappa t}}{\kappa}x}\rho_0(dx)dt<\infty,\end{align}

then the assumption $(\textbf{\textrm{A}}_{g})$ is satisfied.

Define $V(x)\,:\!=\,g(x)+\lambda=\lambda-x$ for $x\in\mathbb{R}$ . Thus, the characteristic equation (25) becomes ${\mathcal{H}}_V\phi=0$ . For $\lambda\in\mathbb{R}$ , let us define

(38) \begin{align}\mu(\lambda)\,:\!=\,\left\{ \begin{array}{ll} \displaystyle \frac1\kappa\left(\lambda-\theta+\frac{\sigma^2}{2\kappa^2}\right), & \kappa>0; \\[9pt] \displaystyle \frac1{|\kappa|}\left(\lambda-\theta+\frac{\sigma^2}{2\kappa^2}+\kappa\right), & \kappa<0. \end{array}\right.\end{align}

We take $\lambda$ so that $\mu(\lambda)\leq0$ . Then, by Lemma 3 together with Proposition E.1 and Theorem E.2 in the Supplementary Material of Qin and Linetsky [Reference Qin and Linetsky18], the associated semigroup generated by ${\mathcal{H}}_V$ (depending on $\lambda$ ) has invariant functions given as follows:

1. (i) For $\mu(\lambda)<0$ , the invariant functions are of the form

   (39) \begin{align} \phi_\lambda(x)=C_1e^{\frac{z^2(x)\epsilon}4}E_\mu(z(x)-\alpha)+C_2e^{\frac{z^2(x)\epsilon}4}E_\mu(\alpha-z(x)),\quad C_1,C_2>0,\end{align}
   where
   \begin{equation*}\alpha\,:\!=\,\sigma\sqrt{\frac2{|\kappa|^3}}, \qquad z(x)\,:\!=\,\frac{\sqrt{2|\kappa|}}\sigma(\theta-x),\end{equation*}
   and $\epsilon\,:\!=\,{\textrm{sign}}(\kappa)$ .

2. (ii) For $\mu(\lambda)=0$ , the invariant functions are of the form

   (40) \begin{align} \phi_\lambda(x)=C_1e^{\frac{z(x)^2\epsilon}4}E_\mu(z(x)-\alpha),\quad C_1>0.\end{align}

Here, $E_\mu(\cdot)$ denotes the Weber parabolic cylinder function. Note that the following recurrence relations on the Weber parabolic cylinder function hold:

\begin{align*}E_{\mu+1}(z)-zE_{\mu}(z)+\mu E_{\mu-1}(z)=0,\qquad E^{\prime}_\mu(z)+\frac12zE_\mu(z)-\mu E_{\mu-1}(z)=0.\end{align*}

Then we arrive at the following relations:

(41) \begin{align}\begin{cases}\displaystyle \frac d{dz}\left(e^{\frac {z^2}4}E_{\mu}(z-\alpha)\right)=\frac\alpha2e^{\frac{z^2}4}E_\mu(z-\alpha)+\mu e^{\frac{z^2}4}E_{\mu-1}(z-\alpha);\\[12pt] \displaystyle \frac d{dz}\left(e^{\frac {z^2}4}E_{\mu}(\alpha-z)\right)=\frac\alpha2e^{\frac{z^2}4}E_\mu(\alpha-z)-\mu e^{\frac{z^2}4}E_{\mu-1}(\alpha-z);\\[12pt] \displaystyle \frac d{dz}\left(e^{-\frac {z^2}4}E_{\mu}(z-\alpha)\right)=-\frac\alpha2e^{-\frac {z^2}4}E_{\mu}(z-\alpha)-e^{\frac {z^2}4}E_{\mu+1}(z-\alpha);\\[12pt] \displaystyle \frac d{dz}\left(e^{-\frac {z^2}4}E_{\mu}(\alpha-z)\right)=-\frac\alpha2e^{-\frac {z^2}4}E_{\mu}(\alpha-z)+e^{\frac {z^2}4}E_{\mu+1}(\alpha-z).\end{cases}\end{align}

We next derive the dynamics of $\overline{X}^y=\big(\overline{X}^y_t\big)_{t\in[0,T]}$ given in Theorem 1(ii). Then it is enough to compute the function $b(\cdot)+\sigma\overline{\sigma}(\cdot)$ with

\begin{equation*}\overline{\sigma}(\cdot)\,:\!=\,\frac{\phi^{\prime}_{\lambda}(\cdot)}{\phi_{\lambda}(\cdot)}\sigma.\end{equation*}

We compute this function in terms of the sign of $\kappa$ by defining

\begin{equation*}b^{+}(\cdot)\,:\!=\,b(\cdot)+\sigma^2\frac{\phi^{\prime}_\lambda(\cdot)}{\phi_\lambda(\cdot)}\end{equation*}

when $\kappa>0$ and

\begin{equation*}b^{-}(\cdot)\,:\!=\,b(\cdot)+\sigma^2\frac{\phi^{\prime}_\lambda(\cdot)}{\phi_\lambda(\cdot)}\end{equation*}

when $\kappa<0$ . First, consider $\kappa>0$ ; then we have $\epsilon=1$ in (39). This yields that

(42) \begin{align} \phi^{\prime}_\lambda(x)&=C_1\frac{dz(x)}{dx}\cdot\frac{d}{dz(x)}\left(e^{\frac{z^2(x)}4}E_\mu(z(x)-\alpha)\right)\nonumber\\ &\quad+C_2\frac{dz(x)}{dx}\cdot\frac{d}{dz(x)}\left(e^{\frac{z^2(x)}4}E_\mu(\alpha-z(x))\right).\end{align}

It follows from the first and the second equations in (41) that

\begin{align*} &\frac{dz(x)}{dx}\cdot\frac{d}{dz(x)}\left(e^{\frac{z^2(x)}4}E_\mu(z(x)-\alpha)\right)\\[4pt] &\qquad=-\frac{\sqrt{2|\kappa|}}\sigma e^{\frac{z^2(x)}4}\left[\frac\alpha2E_\mu(z(x)-\alpha)+\mu E_{\mu-1}(z(x)-\alpha)\right],\\[4pt] &\frac{dz(x)}{dx}\cdot\frac{d}{dz(x)}\left(e^{\frac{z^2(x)}4}E_\mu(\alpha-z(x))\right)\\[4pt] &\qquad=-\frac{\sqrt{2|\kappa|}}\sigma e^{\frac{z^2(x)}4}\left[\frac\alpha2E_\mu(\alpha-z(x))-\mu E_{\mu-1}(\alpha-z(x))\right].\end{align*}

Plugging the above equalities into (42) and rearranging terms, it follows that

\begin{align*} \phi^{\prime}_\lambda(x)&=-\frac{\sqrt{2|\kappa|}}\sigma\cdot\frac\alpha2\left[C_1 e^{\frac{z(x)^2}4}E_\mu(z(x)-\alpha)+C_2e^{\frac{z(x)^2}4}E_\mu(\alpha-z(x))\right]\notag\\[4pt] &\quad-\frac{\sqrt{2|\kappa|}}\sigma e^{\frac{z(x)^2}4}\mu\left[C_1 E_{\mu-1}(z(x)-\alpha)-C_2E_{\mu-1}(\alpha-z(x))\right]\notag\\[4pt] &=-\frac1{|\kappa|}\phi_\lambda(x)-\frac{\sqrt{2|\kappa|}}\sigma e^{\frac{z(x)^2}4}\mu\left[C_1 E_{\mu-1}(z(x)-\alpha)-C_2E_{\mu-1}(\alpha-z(x))\right].\notag\end{align*}

As a consequence, for $\kappa>0$ ,

(43) \begin{align} \bar{b}^+(x) &=\kappa(\theta-x)-\frac{\sigma^2}{|\kappa|}-\sigma\mu\sqrt{2|\kappa|}\frac{C_1 E_{\mu-1}(z(x)-\alpha)-C_2E_{\mu-1}(\alpha-z(x))}{C_1 E_{\mu-1}(z(x)-\alpha)+C_2E_{\mu-1}(\alpha-z(x))}\notag\\[4pt] &=\kappa\theta-\frac{\sigma^2}{\kappa}-\kappa x+\sigma\mu\sqrt{2\kappa}\frac{C_2E_{\mu-1}(\alpha-z(x))-C_1 E_{\mu-1}(z(x)-\alpha)}{C_1 E_{\mu-1}(z(x)-\alpha)+C_2E_{\mu-1}(\alpha-z(x))}.\end{align}

Similarly, for the case $\kappa<0$ , we have $\epsilon=-1$ in (39), and the invariant function is given by

\begin{align*} \phi_\lambda(x)=C_1e^{-\frac{z^2(x)}4}E_\mu(z(x)-\alpha)+C_2e^{-\frac{z^2(x)}4}E_\mu(\alpha-z(x)),\quad C_1,C_2>0.\end{align*}

Then, in a similar fashion to the derivation of (43), we can apply the third and the fourth equations in (41) to obtain

(44) \begin{align} \bar{b}^{-}(x)=\kappa\theta+\frac{\sigma^2}\kappa-\kappa x+\sqrt{2|\kappa|}\sigma \frac{C_1E_{\mu+1}\left(z(x)-\alpha\right)-C_2E_{\mu+1}\left(\alpha-z(x)\right)} {C_2E_\mu\left(\alpha-z(x)\right)+C_1E_\mu\left(\alpha-z(x)\right)}.\end{align}

Therefore, we have the dynamics of $\overline{X}^y=\big(\overline{X}^y_t\big)_{t\in[0,T]}$ in terms of the sign of $\kappa\neq0$ :

\begin{align*}\begin{cases}\displaystyle d\overline{X}_t^y=b^{+}(\overline{X}_t^y)+\sigma dW_t, & \kappa>0,\\[5pt] \displaystyle d\overline{X}_t^y=b^{-}(\overline{X}_t^y)+\sigma dW_t, & \kappa<0.\end{cases}\end{align*}

We next derive the density function $\overline{p}(t,y;\,x)$ of $\overline{X}_t^{y}$ for $t\in(0,T]$ . For this purpose, define $W^y_t\,:\!=\,\sigma^{-1}y+W_t$ ,

\begin{equation*}N_{t}^{y,\pm}\,:\!=\,\exp\left(\int_0^t\sigma^{-1}b^{\pm}(\sigma W^y_s)dW^y_s-\frac{1}{2\sigma^2}\int_0^tb^{\pm}(\sigma W^y_s)^2ds \right),\end{equation*}

and $\overline{Y}_t^y\,:\!=\,\sigma^{-1}\overline{X}_t^y$ for $t\in[0,T]$ . From Exercise 5.5.38 in Karatzas and Shreve [Reference Karatzas and Shreve12], it follows that

\begin{align*}\mathbb{P}\left(\overline{Y}_t^{y}\leq z\right)&=\mathbb{E}\left[N_t^{y,{\textrm{sgn}}(\kappa)}{\mathbf{1}}_{W^y_t\leq z}\right]\\&=\int_{-\infty}^{z}\mathbb{E}\left[N_{t}^{y,{\textrm{sgn}}(\kappa)}\Big|W^y_t=r\right]p_0\left(t,\frac{y}{\sigma};\,r\right)dr,\quad (t,z)\in(0,T]\times\mathbb{R},\end{align*}

where

\begin{equation*}p_0(t,y;\,x)\,:\!=\,\frac{1}{\sqrt{2\pi t}}\exp\Big(-\frac{(x-y)^2}{2t}\Big)\end{equation*}

for $(t,x,y)\in(0,T]\times\mathbb{R}^2$ . Therefore

\begin{align*}\mathbb{P}\left(\overline{X}_t^{y}\leq x\right)&=\mathbb{P}\left(\overline{Y}_t^{y}\leq \frac{x}{\sigma}\right)=\int_{-\infty}^{\frac{x}{\sigma}}\Phi^{{\textrm{sgn}}(\kappa)}(t,r,y)p_0\left(t,\frac{y}{\sigma};\,r\right)dr, (t,x,y)\in(0,T]\times\mathbb{R}^2.\end{align*}

Here $\Phi^{\pm}(t,r,y)\,:\!=\,\mathbb{E}\big[N_{t}^{y,\pm}|W^y_t=r\big]$ for all $(t,r,y)\in(0,T]\times\mathbb{R}^2$ . The density function $\overline{p}(t,y;\,x)$ of $\overline{X}_t^{y}$ for $t\in(0,T]$ is then given by

(45) \begin{align}\overline{p}(t,y;\,x)=\frac{1}{\sigma}\Phi^{{\textrm{sgn}}(\kappa)}\left(t,\frac{x}{\sigma},y\right)p_0\left(t,\frac{y}{\sigma};\,\frac{x}{\sigma}\right),\quad (t,x,y)\in(0,T]\times\mathbb{R}^2.\end{align}

Then, by the solution representation (30) in Theorem 1, and (37), for any probability density function $u_0$ satisfying

\begin{equation*}\int_0^T\int_{\mathbb{R}}(x+1)e^{-\frac{1-e^{-\kappa t}}{\kappa}x}u_0(x)dxdt<\infty,\end{equation*}

we have the solution of the extended RM equation (4) given by

\begin{align*}u(t,x)=\frac{\frac{1}{\phi_{\lambda}(x)}\int_{-\infty}^{\infty}\phi_{\lambda}(y)\overline{p}(t,y;\,x)u_0(y)dy}{\int_{-\infty}^{\infty}\frac{\int_{-\infty}^{\infty}\phi(y)\overline{p}(t,y;\,z)u_0(y)dy}{\phi_{\lambda}(z)}dz},\quad (t,x)\in[0,T]\times\mathbb{R}.\end{align*}

Example 3. (The CIR model) Consider $m=n=1$ , $b(x)=a+bx$ , and $\sigma(x)=\sigma\sqrt{x}$ for $x\in D\,:\!=\,(0,\infty)$ , where $a,\sigma>0$ and $b\in\mathbb{R}$ . Let the Feller condition hold, i.e., $2a\geq\sigma^2$ . Then the assumptions $(\textbf{\textrm{A}}_{b,\sigma}\!)$ and ${({\textbf{A}}_{X})}$ are satisfied. In this case, the process $X^{x}=(X_t^{x})_{t\in[0,T]}$ given by (2) becomes the CIR model, i.e.,

(46) \begin{align}dX_t^{x} = (a+bX_t^{x})dt + \sigma\sqrt{X_t^{x}}dW_t,\quad X_0^{x}=x\in D.\end{align}

The fitness function is $g(x)=-x$ for $x\in D$ . Then

\begin{align*} \mathbb{E}\left[\left|g(X_t^{x})\right|\exp\left(-\int_0^tX_s^{x}ds\right)\right]\leq\mathbb{E}\left[X_t^{x}\right]=\left\{ \begin{array}{ll} \displaystyle xe^{bt}+\frac ab(e^{bt}-1), & b\neq0,\\ \displaystyle x+at, & b=0, \end{array} \right.\end{align*}

and hence the assumption $(\textbf{\textrm{A}}_{g})$ holds if $\rho_0\in{\mathcal{P}}_1(D)$ . Define $V(x)\,:\!=\,g(x)+\lambda=\lambda-x$ for $x\in D$ . Thus, the characteristic equation (25) becomes ${\mathcal{H}}_V\phi=0$ . By Lemma 3 together with Proposition E.1 and Theorem E.2 in the Supplementary Material of Qin and Linetsky [Reference Qin and Linetsky18], for

\begin{equation*}\lambda\leq\lambda_0\,:\!=\,\frac{a(\sqrt{b^2+2\sigma^2}+b)}{\sigma^2},\end{equation*}

the semigroup generated by ${\mathcal H}_V$ (depending on $\lambda$ ) has positive invariant functions $\phi_\lambda$ with the form given by

(47) \begin{align} \phi_\lambda(x)=C_1e^{\frac{\kappa-\gamma}{\sigma^2}x}K\left(\alpha,\beta,\frac{2\gamma}{\sigma^2}x\right),\quad x>0, C_1>0.\end{align}

Here $\kappa=-b$ , $\gamma=\sqrt{\kappa^2+2\sigma^2}$ , $\alpha=\frac{\lambda-\lambda_0}\gamma$ , $\beta=\frac{2a}{\sigma^2}$ , and $K(\cdot,\cdot,\cdot)$ is the Kummer confluent hypergeometric function. By plugging (47) into (28) and (30), we obtain the dynamic of $\overline{X}^y_t$ given in the condition (ii) of Theorem 1 as follows:

(48) \begin{align} d\overline{X}_t^y=\left\{a-\gamma \overline{X}^y_t+\frac{2\alpha\gamma}\beta\frac{K\left(\alpha+1,\beta+1,\frac{2\gamma}{\sigma^2}\overline{X}^y_t\right)} {K\left(\alpha,\beta,\frac{2\gamma}{\sigma^2}\overline{X}^y_t\right)}\overline{X}^y_t\right\}dt+\sigma\sqrt{\overline{X}^y_t}dW_t.\end{align}

For $(t,x,y)\in[0,T]\times(0,\infty)^2$ , let $\overline{p}(t,y;\,x)$ be the transition density function of $\overline{X}_t^y$ (note that $\overline{p}(t,y;\,x)$ does not admit a closed form in general). Then, by (30) in Theorem 1, for any probability density function $u_0(x)$ on $x\in(0,\infty)$ , we have the solution of the extended RM equation (4) given by

\begin{align*}u(t,x)=\frac{\frac{1}{\phi_{\lambda}(x)}\int_{0}^{\infty}\phi_{\lambda}(y)\overline{p}(t,y;\,x)u_0(y)dy}{\int_{0}^{\infty}\frac{\int_{0}^{\infty}\phi_{\lambda}(y)\overline{p}(t,y;\,z)u_0(y)dy}{\phi_{\lambda}(z)}dz},\quad (t,x)\in[0,T]\times(0,\infty).\end{align*}

Example 4. Consider $m=n=1$ , $b(x)\equiv0$ , and $\sigma(x)\equiv\sqrt{2}\sigma$ for $x\in D=\mathbb{R}$ , where $\sigma>0$ . Then the assumptions $(\textbf{\textrm{A}}_{b,\sigma}\!$ ) and ${({\textbf{A}}_{X})}$ are satisfied. In this case, the process $X^{x}=(X_t^{x})_{t\in[0,T]}$ given by (2) satisfies

(49) \begin{align}X_t^{x} = x + \sqrt{2}\sigma W_t,\quad X_0^x=x\in\mathbb{R}.\end{align}

Here, we consider the polynomial confining fitness function of degree 2q which is introduced in Alfaro and Veruete [Reference Alfaro and Veruete3]; i.e., for a given $q\in\mathbb{N}$ ,

(50) \begin{align}g(x) = -x^{2q} + \sum_{l=0}^{2q-1}\alpha_lx^l,\quad x\in\mathbb{R},\end{align}

where $(\alpha_l)_{l=0,\ldots,2q-1}$ is a sequence of finite real numbers. Then the resulting equation (3) becomes the integro-differential Cauchy problem studied in [Reference Alfaro and Veruete3]. Note that the operator ${\mathcal{A}}$ acting on $C^2(\mathbb{R})$ is given by

\begin{equation*}{\mathcal{A}}f(x)=\sigma^2f^{\prime\prime}(x),\quad x\in\mathbb{R}.\end{equation*}

The eigenpair $(\lambda,\phi)$ of the characteristic equation (25) is related to the spectral basis of Schrödinger operator ${\mathcal{H}}\,:\!=\,-{\mathcal{A}}-g$ with confining potential $-g$ . By virtue of (50), it can be seen that $g\,:\,\mathbb{R}\to\mathbb{R}$ is continuous and $\lim_{|x|\to\infty}g(x)=-\infty$ . Then, by Proposition 2.1 in [Reference Alfaro and Veruete3], ${\mathcal{H}}$ has discrete spectrum, i.e., there is an orthonormal basis $(\phi_k)_{k\geq0}\subset L^2(\mathbb{R})$ consisting of eigenfunctions of ${\mathcal{H}}$ ,

\begin{align*}{\mathcal{H}}\phi_k=\lambda_k\phi_k,\quad \|\phi_k\|_{L^2(\mathbb{R})}=1,\end{align*}

with corresponding eigenvalues

\begin{equation*}\lambda_0<\lambda_1\leq\lambda_2\leq\cdots\leq\lambda_k\to\infty\end{equation*}

of finite multiplicity. Moreover, Proposition 2.6 in [Reference Alfaro and Veruete3] gives that

\begin{equation*}\|\phi_k\|_{L^{\infty}(\mathbb{R})}\leq Ck^{\frac{q}{2(q+1)}},\quad k\geq1.\end{equation*}

Up to subtracting a constant from g, we can assume without loss of generality that $g\leq0$ . Note that $X_t^x\in N\big(x,2\sigma^2t\big)$ . This yields the assumption $(\textbf{\textrm{A}}_{g}$ ) if $\rho_0\in{\mathcal{P}}_{2q}(D)$ . Consider $(\lambda_0,\phi_0)$ , i.e., the smallest eigenvalue and the corresponding eigenfunction. According to Sturm–Liouville theory, $\phi_0$ has no zeros on the real line, and hence we may assume that $\phi_0>0$ . Recall the local martingale $M^y$ defined by (26) with the eigenpair $(\lambda_0,\phi_0)$ . Then it holds that $M_0^y=1$ and

\begin{align*} dM_t^y =\frac{\sqrt2\sigma}{\phi_0(y)}\exp\left(\lambda_0 t+\int_0^tg(X_s^y)ds\right)\phi^{\prime}_0(X_t^y)dW_t,\quad t\in(0,T].\end{align*}

Note that the normalized eigenfunction $\phi_k$ for $k\geq0$ satisfies

(51) \begin{align} \sigma^2\int_{\mathbb{R}}|\phi^{\prime}_k(x)|^2dx-\int_{\mathbb{R}}g(x)|\phi_k(x)|^2dx=\lambda_k;\end{align}

this yields that $\phi^{\prime}_k(\pm\infty)=0$ . Therefore $\phi^{\prime}_k$ for $k\geq0$ is bounded on $\mathbb{R}$ . This implies that the local martingale $M^{y}$ defined by (26) with the eigenpair $(\lambda_0,\phi_0)$ is a $(\mathbb{P},\mathbb{F})$ -martingale. In other words, the condition (i) of Theorem 1 holds.

Consider the following SDE:

(52) \begin{align} \overline{Y}_t^{y} = \frac{y}{\sqrt{2}\sigma}+ \int_0^t\overline{\sigma}\left(\sqrt{2}\sigma\overline{Y}_s^{y}\right)ds + W_t,\quad t\in[0,T],\end{align}

where

\begin{equation*}\overline{\sigma}(x)\,:\!=\,\sqrt{2}\frac{\phi^{\prime}_0(x)}{\phi_0(x)}\sigma\end{equation*}

for $x\in\mathbb{R}$ . Then the process $\overline{X}^{y}=\big(\overline{X}^{y}_t\big)_{t\in[0,T]}$ in Theorem 1(ii) is given by

\begin{equation*}\overline{X}_t^{y}=\sqrt{2}\sigma\overline{Y}_t^{y},\quad t\in[0,T].\end{equation*}

By applying Exercise 5.5.38 in Karatzas and Shreve [Reference Karatzas and Shreve12], we arrive at

\begin{align*}\mathbb{P}\left(\overline{Y}_t^{y}\leq z\right)&=\mathbb{E}\left[N_t^{y}{\mathbf{1}}_{W^y_t\leq z}\right]\\[6pt] &=\int_{-\infty}^{z}\mathbb{E}\left[N_t^{y}\big|W^y_t=r\right]p_0\left(t,\frac{y}{\sqrt{2}\sigma};\,r\right)dr,\quad (t,z)\in(0,T]\times\mathbb{R},\end{align*}

where

\begin{equation*}W^y_t=\frac{y}{\sqrt{2}\sigma}+W_t\end{equation*}

and

\begin{equation*}N_t^{y}\,:\!=\,\exp\bigg(\int_0^t\overline{\sigma}\left(\sqrt{2}\sigma W^y_s\right)dW^y_s-\frac{1}{2}\int_0^t\overline{\sigma}^2\Big(\sqrt{2}\sigma W^y_s\Big)ds\bigg)\end{equation*}

for $(t,y)\in[0,T]\times{\mathbb R}$ . Here

\begin{equation*}p_0(t,y;\,x)\,:\!=\,\frac{1}{\sqrt{2\pi t}}\exp\bigg(-\frac{(x-y)^2}{2t}\bigg)\end{equation*}

for $(t,x,y)\in(0,T]\times\mathbb{R}^2$ . Therefore, for all $(t,x,y)\in(0,T]\times\mathbb{R}^2$ ,

\begin{align*}\mathbb{P}\left(\overline{X}_t^{y}\leq x\right)&=\mathbb{P}\left(\overline{Y}_t^{y}\leq \frac{x}{\sqrt{2}\sigma}\right)=\int_{-\infty}^{\frac{x}{\sqrt{2}\sigma}}\Phi(t,r,y)p_0\left(t,\frac{y}{\sqrt{2}\sigma};\,r\right)dr.\end{align*}

Here $\Phi(t,r,y)\,:\!=\,\mathbb{E}[N_t^{y}|W^y_t=r]$ for all $(t,r,y)\in(0,T]\times\mathbb{R}^2$ . Then the density function $\overline{p}(t,y;\,x)$ of $\overline{X}_t^{y}$ for $t\in(0,T]$ exists, and it is given by

(53) \begin{align}\overline{p}(t,y;\,x)=\frac{1}{\sqrt{2}\sigma}\Phi\left(t,\frac{x}{\sqrt{2}\sigma},y\right)p_0\left(t,\frac{y}{\sqrt{2}\sigma};\,\frac{x}{\sqrt{2}\sigma}\right),\quad (t,x,y)\in(0,T]\times\mathbb{R}^2.\end{align}

Thus, the function

(54) \begin{align}u(t,x)=\frac{\frac{1}{\phi_0(x)}\int_{\mathbb{R}}\phi_0(y)\overline{p}(t,y;\,x)u_0(y)dy}{\int_{\mathbb{R}}\frac{\int_{\mathbb{R}}\phi_0(y)\overline{p}(t,y;\,z)u_0(y)dy}{\phi_0(z)}dz},\quad (t,x)\in[0,T]\times\mathbb{R},\end{align}

is a (weak) solution of the extended RM equation (4).

An example with multidimensional fitness space and quadratic fitness function is provided below.

Example 5. Consider $m=n\geq1$ , $b(x)=b+ Bx$ , and $\sigma(x)\equiv\sigma$ for $x\in D\,:\!=\,\mathbb{R}^n$ , where $b\in\mathbb{R}^{n}$ , $B\in\mathbb{R}^{n\times n}$ , and $\sigma\in\mathbb{R}^{n\times n}$ is an invertible matrix, so that $a\,:\!=\,\sigma\sigma^{\top}$ is strictly positive definite. The fitness function is given by $g(x)=-r(x)$ for $x\in\mathbb{R}^n$ , and

(55) \begin{align}r(x)=\alpha+ \delta^{\top}x+ x^{\top}Gx,\quad x\in\mathbb{R}^n,\end{align}

where $\alpha\in\mathbb{R}$ , $\delta\in\mathbb{R}^{n}$ , and $G\in\mathbb{R}^{n\times n}$ is a symmetric positive semi-definite matrix. These quantities are taken to be such that r(x) is nonnegative for all $x\in\mathbb{R}^n$ . Then the assumption ( $\textbf{\textrm{A}}_{g}$ ) holds if $\rho_0\in{\mathcal{P}}_2(D)$ . Let $V=g(x)+\lambda=\lambda-r(x)$ for $x\in\mathbb{R}^n$ ; then the characteristic equation (25) becomes ${\mathcal{H}}_V\phi=0$ . Let $H\in\mathbb{R}^{n\times n}$ satisfy the following continuous-time algebraic Riccati equation:

(56) \begin{align}2HaH-B^{\top}H-HB-G=0.\end{align}

Let $v\in\mathbb{R}^n$ satisfy the following linear equation:

(57) \begin{align}2Hav -B^{\top}v-2Hb-\delta=0.\end{align}

By Section H in the Supplementary Material of Qin and Linetsky [Reference Qin and Linetsky18], the eigenpair $(\lambda,\phi)$ of the characteristic equation (25) which satisfies the condition (i) of Theorem 1 is given by

(58) \begin{align}(\lambda,\phi(x))=\left(\alpha-\frac{1}{2}v^{\top}av+{\textrm{tr}}(aV)+v^{\top}b,\exp\left(-v^{\top}x-x^{\top}Hx\right)\right),\quad x\in\mathbb{R}^n.\end{align}

Define $\overline{\sigma}(x)=-(v+2Hx)^{\top}\sigma$ for $x\in\mathbb{R}^n$ . Then the process $\overline{X}^y=\big(\overline{X}_t^y\big)_{t\in[0,T]}$ given in the condition (ii) of Theorem 1 satisfies

(59) \begin{align} \overline{X}_t^y = y+ \int_0^t\left(b+B\overline{X}_s^y+\sigma\overline{\sigma}\left(\overline{X}_s^y\right)^{\top}\right)ds + \sigma W_t,\quad (t,y)\in[0,T]\times\mathbb{R}^n.\end{align}

Note that the process $\overline{X}^y$ is an n-dimensional affine process; hence the transition density function $\overline{p}(t,y;\,x)$ of $\overline{X}^y$ exists, and the characteristic function of $\overline{p}(t,y;\,\cdot)$ has exponential-affine dependence on $y\in\mathbb{R}^n$ (see Theorem 2.7 of Duffie et al. [Reference Duffie, FilipoviĆ and Schachermayer7]). Then (59) reduces to the following linear SDE:

(60) \begin{align} \overline{X}_t^y = y+ \int_0^t\big\{b-\sigma\sigma^\top v+(B-2\sigma\sigma^\top H)\overline{X}_s^y\big\}ds + \sigma W_t,\quad (t,y)\in[0,T]\times\mathbb R^n.\end{align}

By Proposition IX-2.3 on p. 378 of Revuz and Yor [Reference Revuz and Yor19], it follows that the linear SDE (60) has the following closed-form solution:

\begin{align*} \overline{X}_t^y=e^{\Gamma t}y+\int_0^te^{\Gamma(t-s)}\beta ds+\int_0^te^{\Gamma(t-s)}\sigma dW_s,\quad (t,y)\in[0,T]\times\mathbb{R}^n,\end{align*}

where $\beta\,:\!=\,b-\sigma\sigma^\top v$ and $\Gamma\,:\!=\,B-2\sigma\sigma^\top H$ . This yields that

\begin{align*}\overline{X}^y_t\sim N\left(\mu_t,\Sigma_t\right),\quad t\in(0,T],\end{align*}

where

\begin{equation*}\mu_t\,:\!=\,e^{\Gamma t}y+\int_0^te^{\Gamma(t-s)}\beta ds\end{equation*}

and

\begin{equation*}\Sigma_t\,:\!=\,\int_0^te^{\Gamma(t-s)}\sigma\sigma^\top\left(e^{\Gamma(t-s)}\right)^\top ds.\end{equation*}

Therefore, the transition density function of $\overline{X}_t^y$ is given by

(61) \begin{align}\overline{p}(t,y;\,x)=\frac1{\sqrt{(2\pi)^n{\textrm{det}}(\Sigma_t)}}\exp\left(-\frac12(x-\mu_t)^\top\Sigma_t^{-1}(x-\mu_t)\right),\end{align}

for $ (t,x,y)\in(0,T]\times\mathbb{R}^{2n}$ . Moreover, by (30) in Theorem 1, for any probability density function $u_0(x)$ on $x\in\mathbb{R}^n$ , we have the solution of the extended RM equation (4) given by

\begin{align*}u(t,x)=\frac{\frac{1}{\phi(x)}\int_{\mathbb{R}^n}\phi(y)\overline{p}(t,y;\,x)u_0(y)dy}{\int_{\mathbb{R}^n}\frac{\int_{\mathbb{R}^n}\phi(y)\overline{p}(t,y;\,z)u_0(y)dy}{\phi(z)}dz},\quad (t,x)\in[0,T]\times\mathbb{R}^n.\end{align*}

4.1. Particle system

Building upon the probabilistic solution form of the extended RM equation discussed in Section 3, we construct a suitable coupling between a particle system and the realization of the solution to the FPK equation (5).

Before introducing the state dynamics of the particle system proposed in this section, we first impose the following assumption on the fitness function used in this section:

1. (A _gl ) (i) The fitness function $g\,:\,D\to\mathbb{R}$ is continuous and bounded from above. (ii) There exists a polynomial $Q_g\,:\,\mathbb{R}_+\to\mathbb{R}_+$ such that

   (62) \begin{align} |g(x)-g(y)|\leq Q_g(|x|+|y|)|x-y|,\quad \forall x,y\in D. \end{align}

Note that the condition (i) on the boundedness from above of g in the assumption ${({\textbf{A}}_{gl})}$ with ( $\textbf{\textrm{A}}_{X_0}$ ) below implies the assumption ( $\textbf{\textrm{A}}_{g}$ ) imposed in Lemma 3. It can be seen that the fitness functions in the examples presented in Section 3 satisfy the condition (62).

We next describe a particle system. More precisely, applying the solution representation (6) of the FPK equation given in Lemma 1, we introduce the following particle system: let the number of particles be given by $N\geq1$ . For $i=1,\ldots,N$ , the dynamics of the state process of the ith particle is given by

(63) \begin{align}dX_t^i = b\left(X_t^i\right)dt + \sigma\left(X_t^i\right)dW_t^i,\quad X_0^i\in D, t\in[0,T],\end{align}

where $W^i=\big(W_t^i\big)_{t\in[0,T]}$ , $i=1,\ldots,N$ , and $W=(W_t)_{t\in[0,T]}$ are independent (m-dimensional) Brownian motions under the filtered probability space $(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ . Here, the initial values $(X_0^i)_{i\geq1}$ are assumed to satisfy the following

1. (A _{X 0} ) For $q\geq2$ , the sequence of random variables $\left(X_0^i\right)_{i\geq1}$ is independent and identically distributed (i.i.d.) according to the probability distribution

   \begin{equation*}\rho_0\in{\mathcal{P}}_{({\textrm{deg}}(Q_g)+1)q}(D),\end{equation*}
   where we recall from Section 3.1 that
   \begin{equation*}{\mathcal{P}}_{({\textrm{deg}}(Q_g)+1)q}(D)\end{equation*}
   is the space of probability measures on ${\mathcal{B}}_D$ with finite $({\textrm{deg}}(Q_g)+1)q$ th-order moment, and ${\textrm{deg}}(Q_g)$ denotes the degree of the polynomial $Q_g$ which is given in ${({\textbf{A}}_{gl})}$ .

We also strengthen the assumption on the coefficients $(b,\sigma)$ of the SDE (2) imposed in ( $\textbf{\textrm{A}}_{b,\sigma}\!$ ), used in the previous section, as follows: ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ The functions $b\,:\,D\to\mathbb{R}^{n\times1}$ and $\sigma\,:\,D\to\mathbb{R}^{n\times m}$ are Lipschitz continuous with linear growth.

The current assumption ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ implies the assumption ${({\textbf{A}}_{X})}$ imposed in Section 2.1, and moreover, under the assumption ( $\textbf{\textrm{A}}_{X_0}$ ), for $q\geq2$ , there exists a constant $C=C_{g,q,b,\sigma,T}>0$ such that

(64) \begin{align}\mathbb{E}\left[\sup_{t\in[0,T]}\left|X_t^i\right|^{({\textrm{deg}}(Q_g)+1)q}\right]\leq C\left\{1+\mathbb{E}\left[|X_0^i|^{({\textrm{deg}}(Q_g)+1)q}\right]\right\}<\infty;\end{align}

see also Chapter 5 in Karatzas and Shreve [Reference Karatzas and Shreve12].

Motivated by the solution representation (6) of the FPK equation (5) given in Lemma 1, we introduce the following sequence of ${\mathcal{P}}(D)$ -valued processes:

(65) \begin{align}\mu_t^N = \frac{\frac{1}{N}\sum_{i=1}^N\delta_{X_t^i}\exp\left(\int_0^tg(X_s^i)ds\right)}{\frac{1}{N}\sum_{i=1}^N\exp\left(\int_0^tg(X_s^i)ds\right)},\qquad {\textrm{on}}\ {\mathcal{B}}_{D}.\end{align}

Similarly to (9), for $N\geq1$ , the ${\mathcal{P}}(D)$ -valued process $\mu^N=(\mu_t^N)_{t\in[0,T]}$ also satisfies

(66) \begin{align}\mu_t^N= \frac{1}{N}\sum_{i=1}^N\delta_{X_t^i}\exp\left(\int_0^t(g(X_s^i)-\langle\mu_s^N,g\rangle)ds\right),\quad {\textrm{on}}\ {\mathcal{B}}_{D}.\end{align}

By the representation (66), $(\mu^N)_{N\geq1}$ is in fact a sequence of empirical measure-valued processes with an exponential weight.

4.2. Metric between $\mu^N$ and $\mu$

This section introduces an appropriate metric between $\mu^N$ and $\mu$ for establishing the propagation of chaos. We first give the following spaces. Let ${\mathcal{P}}^s(D)$ be the set of sub-probability measures on ${\mathcal{B}}_{D}$ ; i.e., for any $\mu\in{\mathcal{P}}^s(D)$ , $\mu$ is a finite measure on ${\mathcal{B}}_{D}$ such that $\mu(D)\leq1$ . We next introduce the (Alexandroff) one-point compactification. Add to D one point which is outside of D, called ‘ $\star$ ’, and define $D_{\star}\,:\!=\,D\cup\{\star\}$ . Let D be topologized by a topology ${\mathcal{T}}$ ; we can then define a topology ${\mathcal{T}}^{\star}$ for $D_{\star}$ as follows: (i) each open subset of D is also in ${\mathcal{T}}^{\star}$ , i.e., ${\mathcal{T}}\subset{\mathcal{T}}^{\star}$ ; (ii) for each compact set $C\subset D$ , define an element $U_C\in{\mathcal{T}}^{\star}$ by $U_C\,:\!=\,(D\setminus C)\cup\{\star\}$ . Let us define a bijection $\iota\,:\,{\mathcal{P}}^s(D)\to{\mathcal{P}}(D_{\star})$ by

(67) \begin{align}(\iota\mu)(A) \,:\!=\, \mu(A\cap D)+(1-\mu(D))\delta_{\star}(A),\quad \forall A\in{\mathcal{B}}(D_{\star}), \mu\in{\mathcal{P}}^s(D).\end{align}

Then the integral of $\mu\in{\mathcal{P}}^s(D)$ with respect to a measurable function $f\,:\,D_{\star}\to\mathbb{R}$ is defined as follows (if it is well-defined):

(68) \begin{align}\int_{D_{\star}}f(x)(\iota\mu)(dx)=\int_{D}f(x)\mu(dx)+f(\star)(1-\mu(D))=\langle\mu,f\rangle+f(\star)(1-\mu(D)).\end{align}

Consider $\mu=(\mu_{t})_{t\in[0,T]}$ as an arbitrary ${\mathcal{P}}(D)$ -valued solution of the FPK equation (5). We then define

(69) \begin{align}\mu_t^g\,:\!=\,\exp\left(\int_0^t\langle\mu_s,g\rangle ds\right)\mu_t,\quad t\in[0,T].\end{align}

This implies that

(70) \begin{align}\langle\mu_t^g,f\rangle=\langle\rho_0,f\rangle+\int_0^t\langle\mu_s^g,({\mathcal{A}}+g)f\rangle ds,\quad\forall f\in{\mathcal{D}}.\end{align}

In view of (1.2) and (1.4) in Manita et al. [Reference Manita, Romanov and Shaposhnikov16], we have the following equivalent representation of (70): for all $\phi\in C_b^{1,2}([0,T]\times D)$ ,

(71) \begin{align} \left\langle\mu_t^g,\phi(t,\cdot)\right\rangle=\left\langle\mu_0^g,\phi(0,\cdot)\right\rangle+\int_0^t\left\langle\mu_s^g,(\partial_t+{\mathcal{A}}+g)\phi(s,\cdot)\right\rangle ds.\end{align}

For the sequence of empirical measure-valued processes $(\mu^N)_{N\geq1}$ defined by (65), similarly to (69), we can also define

(72) \begin{align} \mu_t^{g,N}\,:\!=\,\frac{1}{N}\sum_{i=1}^N\delta_{X_t^i}\exp\left(\int_0^tg(X_s^i)ds\right),\quad t\in[0,T].\end{align}

Hereafter, without loss of generality, we assume that $g(x)\leq0$ for all $x\in D$ by the assumption ${({\textbf{A}}_{gl})}$ . Then, for any $t\in(0,T]$ , it follows from (69) and (72) that $\mu_t^g$ and $\mu_t^{g,N}$ are sub-probability measures on ${\mathcal{B}}_{D}$ , i.e. $\mu_t^g,\mu_t^{g,N}\in{\mathcal{P}}^s(D)$ (note that $\mu_0^{g,N}=\mu^N$ and $\mu_0^g=\mu_0$ , and hence $\mu_0^{g,N},\mu_0^g\in{\mathcal{P}}(D)$ ).

For the parameter $q\geq2$ given in the assumption ( $\textbf{\textrm{A}}_{X_0}$ ), we next establish a metric $d_{q,T}$ between $\mu=(\mu_t)_{t\in[0,T]}$ and $\mu^N=(\mu_t^N)_{t\in[0,T]}$ as follows (see also Lucon and Stannat [Reference LuÇon and Stannat14]):

(73) \begin{align}d_{q,T}(\mu,\mu^N)\,:\!=\,\sup_{t\in[0,T]}d_{BL}\left(\iota\mu_t^{g,N},\iota\mu_t^{g}\right),\end{align}

where $d_{BL}$ is a distance for ${\mathcal{P}}(D)$ -valued random variables, which is defined as

(74) \begin{align}d_{BL}\left(\iota\mu_t^{g,N},\iota\mu_t^{g}\right)=\sup_{\psi\in{\mathcal{R}}_1}\mathbb{E}\left[\left|\int_{D_{\star}}\psi(x)\left(\iota\mu_t^{g,N}-\iota\mu_t^{g}\right)(dx)\right|^q\right]^{\frac1q}.\end{align}

Here, ${\mathcal{R}}_1$ is the set of (bounded) Lipschitz continuous functions $\psi\,:\,D_{\star}\to\mathbb{R}$ satisfying $\|\psi\|_{\infty}+\|\psi\|_{\textrm{Lip}}\leq1$ (where

\begin{equation*}\|\psi\|_{\infty}\,:\!=\,\sup_{x\in D_{\star}}|\psi(x)|\end{equation*}

and $\|\psi\|_{\textrm{Lip}}$ denotes the Lipschitz coefficient of $\psi$ ). Let $|\cdot|$ be the Euclidean norm. Then we can define a metric $d_{\star}$ on $D_{\star}$ as in Mandelkern [Reference Mandelkern15]: fix $x_0\in D$ and define

\begin{equation*}l(x)\,:\!=\,\frac{1}{1+|x-x_0|}\end{equation*}

for $x\in D$ . For $x_1,x_2\in D$ , define $d_{\star}(x_1,x_2)\,:\!=\,|x_1-x_2|\wedge(l(x_1)+l(x_2))$ , $d_{\star}(x,\star)\,:\!=\,l(x)$ for $x\in D$ , and $d_{\star}(\star,\star)=0$ . Then the Lipschitz coefficient of $\psi$ (as a seminorm) is given by

\begin{equation*}\|\psi\|_{\textrm{Lip}}=\sup_{x_1\neq x_2,x_1,x_2\in D_{\star}}\frac{\left|\psi(x_1)-\psi(x_2)\right|}{d_{\star}(x_1,x_2)}.\end{equation*}

This implies that, for any $\psi\in{\mathcal{R}}_1$ ,

(75) \begin{align}\left|\psi(x_1)-\psi(x_2)\right|\leq |x_1-x_2|\wedge(l(x_1)+l(x_2))\leq |x_1-x_2|,\quad x_1,x_2\in D.\end{align}

In other words, $\psi\in{\mathcal{R}}_1$ is also a (bounded) Lipschitz continuous function on D with the Lipschitz coefficient being less than one.

4.3. Propagator of FPK equation

This section establishes the propagator corresponding to the transformed FPK equation (70), which is defined as follows: for $(s,x)\in[t,T]\times D$ ,

(76) \begin{align}P_{t,s}^gf(x)\,:\!=\,\mathbb{E}\left[f(X_{s}^{t,x})\exp\left(\int_t^sg(X_r^{t,x})dr\right)\right],\quad f\in C_b(D),\end{align}

where the process $X^{t,x}=(X_s^{t,x})_{s\in[t,T]}$ is the unique strong solution of the SDE (2), i.e.,

\begin{align*}X_s^{t,x} = x + \int_t^s b(X_r^{t,x})dr + \int_t^s\sigma(X_r^{t,x})dW_r,\quad s\in[t,T].\end{align*}

The following lemma can be proved by verifying the conditions of Theorem 1 and Lemma 2 in Health and Schweizer [Reference Health and Schweizer10].

Lemma 4. Let the assumptions $(\textbf{\textrm{A}}_{b,\sigma}\!)$ , ${({\textbf{A}}_{X})}$ , and ${({\textbf{A}}_{gl})}$ (i) hold. Also suppose the following:

1. (A _D,σ) There exists a sequence $(D_k)_{k\in\mathbb{N}}$ of bounded domains with $\overline{D}_k\subset D$ such that $\bigcup_{k=1}^{\infty}D_k=D$ , each $D_k$ has a $C^2$ -boundary, and for each $k\geq1$ , $\sigma\sigma^{\top}(x)$ is uniform elliptic on $\mathbb{R}^n$ for $(t,x)\in[0,T)\times D_k$ .

Then the propagator $P_{t,s}^gf$ defined by (76) satisfies the following: for $(s,x)\in[t,T]\times D$ ,

(77) \begin{align}\partial_tP_{t,s}^gf(x) + ({\mathcal{A}}+g)P_{t,s}^gf(x)=0,\quad P_{s,s}^gf(x)=f(x).\end{align}

The operator ${\mathcal{A}}$ is defined by (1). Moreover, $P_{\cdot,s}^gf\in C^{1,2}((0,s]\times D)\cap C([0,s]\times D)$ , and there exists a unique classical solution of the Cauchy problem (77).

We next estimate the gradient of the propagator $P_{t,s}^g$ defined by (76) acting on test functions. Before doing so, in view of (68), we have that, for all $\psi\in{\mathcal{R}}_1$ ,

(78) \begin{align}\int_{D_{\star}}\psi(x)\left(\iota\mu_t^{g,N}-\iota\mu_t^{g}\right)(dx)=\int_{D}\psi_{\star}(x)\left(\mu_t^{g,N}-\mu_t^{g}\right)(dx)=\big\langle\mu_t^{g,N}-\mu_t^{g},\psi_{\star}\big\rangle,\end{align}

where $\psi_{\star}(x)\,:\!=\,\psi(x)-\psi(\star)$ for $x\in D$ . It follows from (75) that $\psi_{\star}$ is a (bounded) Lipschitz continuous function on D with the Lipschitz coefficient being less than one. We then give the estimates of the gradient $\nabla_xP_{t,s}\psi_{\star}$ and $\partial_tP_{t,s}\psi_{\star}$ as follows.

Lemma 5. Let the assumptions ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ , ${({\textbf{A}}_{D,\sigma})}$ , and ${({\textbf{A}}_{gl})}$ hold. Then there exist polynomials $\hat{Q}_g^{(1)}=\hat{Q}^{(1)}_{g,b,\sigma,T}\,:\,\mathbb{R}_+\to\mathbb{R}_+$ and $\hat{Q}_g^{(2)}=\hat{Q}^{(2)}_{g,b,\sigma,T}\,:\,\mathbb{R}_+\to\mathbb{R}_+$ (depending only on $g,b,\sigma,T$ ), with ${\textrm{deg}}(\hat{Q}^{(1)}_g)={\textrm{deg}}(Q_g)$ and ${\textrm{deg}}(\hat{Q}^{(2)}_g)={\textrm{deg}}(Q_g)+1$ , such that for all $\psi\in{\mathcal{R}}_1$ ,

(79) \begin{align}\left|\nabla_xP_{t,s}\psi_{\star}(x)\right|\leq \hat{Q}^{(1)}_g(|x|),\quad \left|\partial_tP_{t,s}\psi_{\star}(x)\right|\leq\hat{Q}^{(2)}_g(|x|),\quad (s,x)\in[t,T]\times D.\end{align}

Proof. Note that $\|\psi_{\star}\|_{\infty}\leq 2$ for all $\psi\in{\mathcal{R}}_1$ . Then, in view of (76) and the fact that g is bounded from above (without loss of generality, we have assumed that $g\leq0$ on D), we have that, for all $x_1,x_2\in D$ ,

(80) \begin{align}&\left|P_{t,s}^g\psi_{\star}(x_1)-P_{t,s}^g\psi_{\star}(x_2)\right|^2\nonumber\\[3pt] &\quad\leq2\mathbb{E}\left[\left|\psi_{\star}(X_s^{t,x_1})\exp\left(\int_t^sg(X_r^{t,x_1})dr\right)-\psi_{\star}(X_s^{t,x_2})\exp\left(\int_t^sg(X_r^{t,x_1})dr\right)\right|^2\right]\notag\\[3pt] &\quad\quad+2\mathbb{E}\left[\left|\psi_{\star}(X_s^{t,x_2})\exp\left(\int_t^sg(X_r^{t,x_1})dr\right)-\psi_{\star}(X_s^{t,x_2})\exp\left(\int_t^sg(X_r^{t,x_2})dr\right)\right|^2\right]\notag\\[3pt] &\quad\leq 2\mathbb{E}\left[\left|\psi_{\star}(X_s^{t,x_1})-\psi_{\star}(X_s^{t,x_2})\right|^2\exp\left(2\int_t^sg(X_r^{t,x_1})dr\right)\right]\nonumber\\[3pt] &\qquad\quad+2\mathbb{E}\left[\left|\psi_{\star}(X_s^{t,x_2})\right|^2\left|\exp\left(\int_t^s g(X_r^{t,x_1})dr\right)-\exp\left(\int_t^s g(X_r^{t,x_2})dr\right)\right|^2\right]\nonumber\\[3pt] &\qquad\quad\leq2\mathbb{E}\left[\left|X_s^{t,x_1}-X_s^{t,x_2}\right|^2\right]+8\mathbb{E}\left[\left|\int_t^s g(X_r^{t,x_1})dr-\int_t^s g(X_r^{t,x_2})dr\right|^2\right]\notag\\[3pt] &\qquad\quad\leq2\mathbb{E}\left[\left|X_s^{t,x_1}-X_s^{t,x_2}\right|^2\right]+8T\int_t^s\mathbb{E}\left[\left|g(X_r^{t,x_1})-g(X_r^{t,x_2})\right|^2\right]dr.\end{align}

We first note that it follows from the assumption ${({\textbf{A}}_{gl})}$ (ii) that, for the polynomial $Q_g\,:\,\mathbb{R}_+\to\mathbb{R}_+$ and for all $r\in[t,s]$ ,

(81) \begin{align}\mathbb{E}\left[\left|g(X_r^{t,x_1})-g(X_r^{t,x_2})\right|^2\right]&\leq\mathbb{E}\left[Q^2_g\left(\left|X_r^{t,x_1}\right|+\left|X_r^{t,x_2}\right|\right)\left|X_r^{t,x_1}-X_r^{t,x_2}\right|^2\right]\\[5pt] &\leq\left\{\mathbb{E}\left[Q^4_g\left(\left|X_r^{t,x_1}\right|+\left|X_r^{t,x_2}\right|\right)\right]\right\}^{\frac12} \left\{\mathbb{E}\left[\left|X_r^{t,x_1}-X_r^{t,x_2}\right|^4\right]\right\}^{\frac12}.\notag\end{align}

On the other hand, by the assumption ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ , it follows from Burkholder–Davis–Gundy (BDG) inequality that, for $k\geq1$ ,

(82) \begin{align}\mathbb{E}\left[\sup_{r\in[t,s]}\left|X_r^{t,x_1}-X_r^{t,x_2}\right|^{2k}\right]\leq |x_1-x_2|^{2k}+C_k\int_t^s \mathbb{E}\left[\left|X_r^{t,x_1}-X_r^{t,x_2}\right|^{2k}\right]dr,\end{align}

where $C_k=C_{b,\sigma,T,k}$ is a positive constant which depends on $b,\sigma,T,k$ . Then Gronwall’s lemma yields that

(83) \begin{align}\mathbb{E}\left[\sup_{r\in[t,s]}\left|X_r^{t,x_1}-X_r^{t,x_2}\right|^{2k}\right]\leq|x_1-x_2|^{2k}e^{TC_k}.\end{align}

It then follows from (80), (81), and (83) that for all $(s,x_1,x_2)\in[t,T]\times D^2$ ,

\begin{align*} &\left|P_{t,s}^g\psi_{\star}(x_1)-P_{t,s}^g\psi_{\star}(x_2)\right|^2\\ &\quad\leq2\mathbb{E}\left[\left|X_s^{t,x_1}-X_s^{t,x_2}\right|^2\right]+8T\int_t^s\mathbb{E}\left[\left|g(X_r^{t,x_1})-g(X_r^{t,x_2})\right|^2\right]dr\notag\\ &\quad\leq2|x_1-x_2|^2e^{TC_{1}}+8T\int_t^s\left\{\mathbb{E}\left[Q^4_g(\left|X_r^{t,x_1}\right|+\left|X_r^{t,x_2}\right|)\right]\right\}^{\frac12} \left\{\mathbb{E}\left[\left|X_r^{t,x_1}-X_r^{t,x_2}\right|^4\right]\right\}^{\frac12}dr\notag\\ &\quad\leq2e^{TC_{1}}|x_1-x_2|^2+8Te^{\frac{TC_{2}} {2}}|x_1-x_2|^2\int_t^s\left\{\mathbb{E}\left[Q^4_g\left(\left|X_r^{t,x_1}\right|+\left|X_r^{t,x_2}\right|\right)\right]\right\}^{\frac12}dr.\end{align*}

It is not difficult to show from (64) that there exists a polynomial ${\tilde Q}_{g}\,:\!=\,{\tilde Q}_{g,b,\sigma,T}\,:\,\mathbb{R}_+\to\mathbb{R}_+$ (which depends only on $g,b,\sigma,T$ ) such that ${\textrm{deg}}({\tilde Q}_{g})=2{\textrm{deg}}(Q_{g})$ , and for all $(s,x_1,x_2)\in[t,T]\times D^2$ ,

\begin{align*} 2e^{TC_{1}}+8Te^{\frac{TC_{2}}{2}}\int_t^s\left\{\mathbb{E}\left[Q^4_g\left(\left|X_r^{t,x_1}\right|+\left|X_r^{t,x_2}\right|\right)\right]\right\}^{\frac12}dr \leq{\tilde Q}_{g}(|x_1|+|x_2|).\end{align*}

Therefore, we arrive at, for all $(s,x_1,x_2)\in[t,T]\times D^2$ ,

(84) \begin{align} \left|P_{t,s}^g\psi_{\star}(x_1)-P_{t,s}^g\psi_{\star}(x_2)\right|^2\leq{\tilde Q}_{g}(|x_1|+|x_2|)|x_1-x_2|^2.\end{align}

Thanks to Lemma 4, we have that $P_{\cdot,s}^g\psi_{\star}\in C^{1,2}((0,s]\times D)\cap C([0,s]\times D)$ by the fact that $\psi_{\star}\in C_b(D)$ . By letting $x_2\rightarrow x_1$ in (84), it follows that

(85) \begin{align} |\nabla_xP_{t,s}\psi_{\star}(x)|\leq{\tilde Q}^{\frac12}_{g}(2|x|),\quad (s,x)\in[t,T]\times D.\end{align}

The estimate (79) then follows from (85) and the fact that

\begin{equation*}{\tilde Q}^{\frac12}_{g}(2|x|)\leq\hat{Q}_{g}^{(1)}(|x|)\end{equation*}

for some polynomial $\hat{Q}_g^{(1)}\,:\!=\,\hat{Q}^{(1)}_{g,b,\sigma,T}\,:\,\mathbb{R}_+\to\mathbb{R}_+$ with

\begin{equation*}{\textrm{deg}}\left({\hat Q}_{g}^{(1)}\right)=\frac12 {\textrm{deg}}({\tilde Q}_{g})={\textrm{deg}}(Q_{g}).\end{equation*}

This completes of the proof of the first estimate in (79).

We next prove the second estimate in (79). Let $t\in[0,T]$ , $\delta\leq T-t$ , and $(s,x)\in$ $[t+\delta,T]\times D$ . It follows from the strong Markov property of $X^{t,x}=(X_s^{t,x})_{s\in[t,T]}$ that

\begin{align*}P_{t+\delta,s}^g\psi_{\star}(x)&=\mathbb{E}\left[\psi_{\star}(X_{s}^{t+\delta,x})\exp\left(\int_{t+\delta}^{s}g(X_r^{t+\delta,x})dr\right)\right]\notag\\&=\mathbb{E}\left[\psi_{\star}(X_{s-\delta}^{t,x})\exp\left(\int_{t}^{s-\delta}g(X_{r}^{t,x})dr\right)\right]=P_{t,s-\delta}^g\psi_{\star}(x).\end{align*}

Therefore, using the fact that $g\leq0$ on D, we obtain

(86) \begin{align}&\left|P_{t+\delta,s}^g\psi_{\star}(x)-P_{t,s}^g\psi_{\star}(x)\right|^2=\left|P_{t,s-\delta}^g\psi_{\star}(x)-P_{t,s}^g\psi_{\star}(x)\right|^2\notag\\&\qquad\quad\leq 2\mathbb{E}\left[\left|\psi_{\star}(X_s^{t,x})-\psi_{\star}(X_{s-\delta}^{t,x})\right|^2\exp\left(2\int_t^sg(X_r^{t,x})dr\right)\right]\nonumber\\&\qquad\qquad+2\mathbb{E}\left[\left|\psi_{\star}(X_{s-\delta}^{t,x})\right|^2\left|\exp\left(\int_t^s g(X_r^{t,x})dr\right)-\exp\left(\int_t^{s-\delta}g(X_r^{t,x})dr\right)\right|^2\right]\nonumber\\&\qquad\quad\leq2\mathbb{E}\left[\left|X_s^{t,x}-X_{s-\delta}^{t,x}\right|^2\right]+8\mathbb{E}\left[\left|\int_{s-\delta}^{s}g(X_r^{t,x})dr\right|^2\right]\notag\\&\qquad\quad\leq2\delta^2\mathbb{E}\left[\sup_{r\in[t,s]}\left|X_r^{t,x}\right|^2\right]+8\delta^2\mathbb{E}\left[\sup_{r\in[t,s]}\left|g(X_r^{t,x})\right|^2\right].\end{align}

Thanks to Lemma 4, the similar argument used in the estimate of (84) and (85), together with (86), (64), and the assumption ${({\textbf{A}}_{gl})}$ (ii), yield the validity of the second estimate in (79). This completes the proof of the lemma.

We next discuss the estimate (79) in Lemma 5 for the one-dimensional case (i.e., $m=n=1$ ), but when the volatility function $\sigma\,:\,D\to\mathbb{R}$ of the SDE (2) is only Hölder continuous. For this purpose, we impose the following assumption:

1. (A _1d) For $m=n=1$ , $b\,:\,D\to\mathbb{R}$ is Lipschitz continuous, and $\sigma\,:\,D\to\mathbb{R}$ is Hölder continuous with exponent $\gamma\in[\frac{1}{2},1)$ .

Proof. Note that under the assumption ( $\textbf{\textrm{A}}_{1\textrm{d}}$ ), the estimate (82) in the proof of Lemma 5 does not hold. However, the assumption (A $_{1\textrm{d}}$ ) implies that the SDE (2) satisfies the Yamada–Watanabe condition given in Proposition 5.2.13 of Karatzas and Shreve [Reference Karatzas and Shreve12]. This proposition therefore yields that, for all $s\in[t,T]$ , there exists a constant $C=C_{b,\sigma,T}>0$ such that

\begin{align*}\mathbb{E}\left[\left|X_s^{t,x_1}-X_s^{t,x_2}\right|\right]\leq |x_1-x_2|+C\int_t^T \mathbb{E}\left[\left|X_s^{t,x_1}-X_s^{t,x_2}\right|\right]ds.\end{align*}

Then Gronwall’s lemma yields that, for all $s\in[t,T]$ ,

(87) \begin{align}\mathbb{E}\left[\left|X_s^{t,x_1}-X_s^{t,x_2}\right|\right]\leq|x_1-x_2|e^{(T-t)C}.\end{align}

Note that under the assumption ${({\textbf{A}}_{gl})}$ with ${\textrm{deg}}(Q_g)=0$ , the fitness function $g\,:\,D\to\mathbb{R}$ is Lipschitz continuous. Similarly to (80), we obtain from (87) that

\begin{align*}\left|P_{t,s}^g\psi_{\star}(x_1)-P_{t,s}^g\psi_{\star}(x_2)\right|&\leq\mathbb{E}\left[\left|X_s^{t,x_1}-X_s^{t,x_2}\right|\right]+2\int_t^T\mathbb{E}\left[\left|g(X_r^{t,x_1})-g(X_r^{t,x_2})\right|\right]dr\nonumber\\&\leq \left\{1+2T\|g\|_{\textrm{Lip}}\right\}e^{TC}|x_1-x_2|.\end{align*}

This proves the first estimate in (79) with $\hat{Q}_g^{(1)}\equiv K\,:\!=\,\{1+2T\|g\|_{\textrm{Lip}}\}e^{TC}$ . We can apply the similar argument used in (86) to prove the validity of the second estimate in (79).

We next show the uniqueness of solutions to the FPK equation (5) using the propagator $P_{t,s}^g$ introduced in (76).

Lemma 7. Let the assumptions $(\textbf{\textrm{A}}_{X_0})$ , ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ , ${({\textbf{A}}_{D,\sigma}\!)}$ , and ${({\textbf{A}}_{gl})}$ (or the assumptions $(\textbf{\textrm{A}}_{X_0})$ , $(\textbf{\textrm{A}}_{1\textrm{d}})$ , ${({\textbf{A}}_{D,\sigma}\!)}$ , and ${({\textbf{A}}_{gl})}$ with ${\textrm{deg}}(Q_g)=0$ ) hold. Let $\mu=(\mu_t)_{t\in[0,T]}$ be a solution of the FPK equation (5) satisfying $\mu_0=\rho_0$ and the integrability condition

(88) \begin{align}\sup_{t\in[0,T]}\int_D |x|^{{\textrm{deg}}(Q_g)+1}\mu_t(dx)<+\infty.\end{align}

Then $\mu$ is unique.

Proof. Let $\hat{\mu}=(\hat{\mu}_t)_{t\in[0,T]}$ be another solution of the FPK equation (5) with $\hat{\mu}_0=\rho_0$ . Similarly to (69), we define $\hat{\mu}_t^g\,:\!=\,\exp\big(\int_0^t\langle\hat{\mu}_r,g\rangle dr\big)\hat{\mu}_t$ for $t\in[0,T]$ . We first show that

(89) \begin{align} \langle\mu_t^g,\psi\rangle=\langle\hat\mu_t^g,\psi\rangle, t\in[0,T], \psi\in{\mathcal{R}}_1.\end{align}

For this purpose, let $s\in[0,T]$ be an arbitrary (fixed) time and define $h(t,x)\,:\!=\,P_{t,s}^g\psi(x)$ with $(t,x)\in[0,s]\times D$ and $\psi\in{\mathcal{R}}_1$ . Note that h is not in $C^{1,2}_b([0,s]\times D)$ . We therefore introduce the cutoff function $\xi\in C^\infty_0(\mathbb R^n)$ which satisfies $\xi(x)\in[0,1]$ for all $x\in\mathbb{R}^n$ , $\xi(x)=1$ for $|x|\leq1$ , and $\xi(x)=0$ for $|x|>2$ . Moreover, for $N\geq1$ , define $\xi_N(x)\,:\!=\,\xi(x/N)$ for $x\in\mathbb{R}^n$ . By Lemma 4, it follows that $h_N(t,x)\,:\!=\,\xi_N(x)h(t,x)$ is in $C^{1,2}_b([0,s]\times D)$ . Therefore, using (71), we have, for all $t\in[0,s]$ ,

\begin{align*} \big\langle\mu_t^g,h_N(t,\cdot)\big\rangle=\big\langle\mu_0^g,h_N(0,\cdot)\big\rangle+\int_0^t\big\langle\mu_r^g,(\partial_t+{\mathcal{A}}+g)h_N(r,\cdot)\big\rangle dr.\end{align*}

In terms of the definition of $h_N$ , one gets from (77) in Lemma 4 that

\begin{align*} (\partial_t+{\mathcal{A}}+g)h_N(t,x)&=\xi_N(x)(\partial_t+{\mathcal{A}}+g)h(t,x)+h(t,x){\mathcal{A}}\xi_N(x)\notag\\ &\quad+\nabla_x\xi_N(x)^\top\sigma\sigma^{\top}\nabla_xh(t,x)\notag\\ &=h(t,x){\mathcal{A}}\xi_N(x)+\nabla_x\xi_N(x)^\top\sigma\sigma^{\top}\nabla_xh(t,x).\end{align*}

Note that both $\nabla^2_x\xi_N$ and $\nabla_x\xi_N$ are bounded (the boundedness is independent of N) and supported on $\{x;\, N\leq|x|\leq2N\}\cap D$ . Moreover, for $N\leq|x|\leq2N$ , there exists a constant C (which is independent of N) such that

\begin{align*} \left|\nabla_x\xi_N(x)\right|=\frac{1}{N}\left|\nabla_x\xi(x/N)\right|\leq\frac C{|x|}.\end{align*}

Let $\alpha\,:\!=\,{\textrm{deg}}(Q_g)+1$ . Then it follows from Lemma 4 that there exists a constant C (which is independent of N) such that

(90) \begin{align}&\left|\int_D\left\{h(t,x){\mathcal{A}}\xi_N(x)+\nabla_x\xi_N(x)^\top\sigma\sigma^{\top}\nabla_xh(t,x)\right\}\mu^g_t(dx)\right|\nonumber\\ &\qquad\leq C\int_{D_N}\left\{|x|^\alpha+1\right\}\mu^g_t(dx),\end{align}

where $D_N\,:\!=\,\{x\,:\,N\leq|x|\leq2N\}\cap D$ . Then, by the assumption (88), i.e.

\begin{equation*}\sup_{t\in[0,T]}\int_D|x|^{\alpha}\mu_t(dx)<\infty,\end{equation*}

and the assumption ${({\textbf{A}}_{gl})}$ (i), there exists a constant C (which is independent of N) such that

\begin{align*} \sup_{t\in[0,T]}\int_{D_N}\left\{|x|^\alpha+1\right\}\mu^g_t(dx)\leq C\sup_{t\in[0,T]}\int_{D}\left\{|x|^\alpha+1\right\}\mu_t(dx)<+\infty,\end{align*}

which implies that the right-hand side of (90) tends to 0 as $N\to\infty$ . Hence, it holds that

\begin{align*}\lim_{N\to\infty}\int_D\left\{h(t,x){\mathcal{A}}\xi_N(x)+\nabla_x\xi_N(x)^\top\sigma\sigma^{\top}\nabla_xh(t,x)\right\}\mu^g_t(dx)=0,\quad t\in[0,T].\end{align*}

This yields that, for all $t\in[0,s]$ ,

(91) \begin{align} \langle\mu_t^g,h(t,\cdot)\rangle=\lim_{N\rightarrow+\infty}\big\langle\mu_t^g,h_N(t,\cdot)\big\rangle=\lim_{N\rightarrow+\infty}\big\langle\mu_0^g,h_N(0,\cdot)\big\rangle =\big\langle\mu_0^g,h(0,\cdot)\big\rangle.\end{align}

Note that $\langle\mu_s^g,h(s,\cdot)\rangle=\langle\mu_s^g,\psi\rangle$ , and hence (91) gives that

(92) \begin{align}\langle\mu_s^g,\psi\rangle=\big\langle\mu_0^g,h(0,\cdot)\big\rangle=\big\langle\rho_0,h(0,\cdot)\big\rangle.\end{align}

The same reasoning yields that

(93) \begin{align} \langle\hat\mu_s^g,\psi\rangle=\big\langle\hat\mu_0^g,h(0,\cdot)\big\rangle=\big\langle\rho_0,h(0,\cdot)\big\rangle.\end{align}

The equality (89) then follows from (92) and (93) with the arbitrariness of $s\in[0,T]$ . If we choose $\psi\equiv1$ , it follows from (89) that

\begin{align*} \exp\left(\int_0^t\langle\mu_s,g\rangle ds\right)=\langle\mu_t^g,1\rangle=\langle\hat\mu_t^g,1\rangle=\exp\left(\int_0^t\langle\hat\mu_s,g\rangle ds\right),\quad t\in[0,T].\end{align*}

Therefore

(94) \begin{align} \langle\mu_t,\psi\rangle=\langle\hat\mu_t,\psi\rangle, t\in[0,T], \psi\in{\mathcal{R}}_1.\end{align}

This gives the uniqueness of solutions of the FPK equation (5).

Proof. It follows from Lemma 1 that $\mu=(\mu_t)_{t\in[0,T]}$ given by (6) is a solution of the FPK equation (5). To prove the corollary, by Lemma 7, it suffices to show that this $\mu$ satisfies the integrability condition (88). In fact, without loss of generality, we have assumed that $g\leq0$ on D. Then, using (64), for $\alpha={\textrm{deg}}(Q_g)+1$ , we have

\begin{align*} \sup_{t\in[0,T]}\int_D|x|^{\alpha}\mu_t(dx)&= \sup_{t\in[0,T]}\frac{\int_{D}\mathbb{E}\left[|X_t^y|^{\alpha}\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)} {\int_{D}\mathbb{E}\left[\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)}\notag\\ &\leq \frac{\int_{D}\mathbb{E}\left[\sup_{t\in[0,T]}|X_t^y|^{\alpha}\right]\rho_0(dy)} {\int_{D}\mathbb{E}\left[\exp\left(\int_0^Tg(X_s^y)ds\right)\right]\rho_0(dy)}\notag\\ &\leq \frac{C+C\int_{D}|y|^{\alpha}\rho_0(dy)} {\int_{D}\mathbb{E}\left[\exp\left(\int_0^Tg(X_s^y)ds\right)\right]\rho_0(dy)}, \end{align*}

where the constant $C>0$ depends only on $g,b,\sigma,T$ . Moreover, the assumption (A $_{X_0}$ ) yields that $\int_S|y|^{\alpha}\rho_0(dy)<\infty$ , and obviously $\int_{D}\mathbb{E}\big[\exp\big(\int_0^Tg\big(X_s^y\big)ds\big)\big]\rho_0(dy)$ is finite since $g\leq0$ on D. This proves that $\mu$ satisfies the integrability condition (88).

Building upon the previous results, the propagator defined by (76) can be used to establish the following relation satisfied by $\mu_T^{g,N}-\mu_T^g$ for any fixed $T>0$ .

Lemma 8. Let the conditions of Lemma 7 hold. Then, for any fixed $T>0$ , it holds that

(95) \begin{align}\big\langle\mu_T^{g,N}-\mu_T^g,f\big\rangle&= \big\langle\mu_0^{g,N}-\mu_0^g,P_{0,T}^gf\big\rangle\nonumber\\&\quad+\frac{1}{N}\sum_{i=1}^N\int_0^T\exp\left(\int_0^sg\left(X_r^i\right)dr\right)\nabla_xP_{s,T}^gf\left(X_s^i\right)^{\top}\sigma\left(X_s^i\right)dW_s^i,\end{align}

where $P_{t,T}^gf$ for $t\in[0,T]$ is the propagator defined by (76).

Proof. Recall the state process of the particle system $X^i=(X_t^i)_{t\in[0,T]}$ defined by (63) for $i\geq1$ . Lemma 4 allows us to apply Itô’s formula to $P_{t,T}^gf(X_t^i)$ with $t\in[0,T]$ , and we have that

\begin{align*}P_{t,T}^gf(X_t^i)&=P_{0,T}^gf(X_0^i)+\int_0^t(\partial_s+{\mathcal{A}})P_{s,T}^gf(X_s^i)ds+\int_0^t\nabla_xP_{s,T}^gf(X_s^i)^{\top}\sigma(X_s^i)dW_s^i.\end{align*}

Thus, by (77) in Lemma 4, we have

\begin{align*}P_{t,T}^gf(X_t^i)\exp\left(\int_0^tg(X_s^i)ds\right)&=P_{0,T}^gf(X_0^i)\nonumber\\&\quad+\int_0^t\exp\left(\int_0^sg(X_r^i)dr\right)\nabla_xP_{s,T}^gf(X_s^i)^{\top}\sigma(X_s^i)dW_s^i.\end{align*}

Recall (72). Then, using (66), it follows that

\begin{align*}\big\langle\mu_t^{g,N},P_{t,T}^gf\big\rangle&= \exp\left(\int_0^t\langle\mu_s^N,g\rangle ds\right)\big\langle\mu_t^N,P_{t,T}^gf\rangle=\frac{1}{N}\sum_{i=1}^NP_{t,T}^gf(X_t^i)\exp\left(\int_0^t g(X_s^{i})ds\right)\nonumber\\ &=\big\langle\mu_0^{g,N},P_{0,T}^gf\big\rangle+\frac{1}{N}\sum_{i=1}^N\int_0^t\exp\left(\int_0^sg(X_r^i)dr\right)\nabla_xP_{s,T}^gf(X_s^i)^{\top}\sigma(X_s^i)dW_s^i.\end{align*}

We first claim that

(96) \begin{align}\partial_t\langle\mu_t^g, P_{t,T}^gf\rangle=0,\quad {\forall} t\in[0,T].\end{align}

In fact, by Corollary 1, we can use the probabilistic representation (6) to validate (96). Thanks to Lemma 4, it follows from Itô’s formula that, for all $(t,y)\in[0,T]\times D$ ,

\begin{align*}&\exp\left(\int_0^tg(X_s^y)ds\right)P_{t,T}^gf(X_t^y)\\&\quad =P_{0,T}^gf(y)+\int_0^t\exp\left(\int_0^sg(X_r^y)dr\right)\left(\partial_t+\mathcal A+g\right)P_{s,T}^gf(X_s^y)ds\\&\qquad+\int_0^t\exp\left(\int_0^sg(X_r^y)dr\right)\nabla_xP_{s,T}^gf(X_s^y)^{\top}\sigma(X_s^y)dW_s\\&\quad=P_{0,T}^gf(y)+\int_0^t\exp\left(\int_0^sg(X_r^y)dr\right)\nabla_xP_{s,T}^gf(X_s^y)^{\top}\sigma(X_s^y)dW_s.\end{align*}

It follows from Lemma 4 again that

\begin{align*} \mathbb E\left[\exp\left(\int_0^tg(X_s^y)ds\right)P_{t,T}^gf(X_t^y)\right]=P_{0,T}^gf(y),\quad (t,y)\in[0,T]\times D.\end{align*}

Then, by virtue of (6), we have

\begin{align*}\displaystyle\langle\mu_t,P_{t,T}^gf\rangle= \frac{\int_{D}\mathbb{E}\left[\exp\left(\int_0^tg(X_s^y)ds\right)P_{t,T}^gf(X_t^y)\right]\rho_0(dy)}{\int_{D}\mathbb{E}\left[\exp\left(\int_0^tg(X_s^y)ds\right)\right]\rho_0(dy)},\quad t\in[0,T].\end{align*}

This yields from (8) and (69) that, for all $t\in[0,T]$ ,

\begin{align*} \langle\mu_t^g, P_{t,T}^gf\rangle=\int_{D}\mathbb E\left[\exp\left(\int_0^tg(X_s^y)ds\right)P_{t,T}^gf(X_t^y)\right]\rho_0(dy) =\int_D P_{0,T}^gf(y)\rho_0(dy).\end{align*}

Since the right-hand side of the above equation is independent of t, (96) holds. Consequently, $\langle \mu_t^g,P_{t,T}^gf\rangle=\langle \mu_0^g,P_{0,T}^gf\rangle$ for all $t\in[0,T]$ . This gives

\begin{align*}\big\langle\mu_t^{g,N}-\mu_t^g,P_{t,T}^gf\big\rangle&= \big\langle\mu_0^{g,N}-\mu_0^g,P_{0,T}^gf\big\rangle\nonumber\\&\quad+\frac{1}{N}\sum_{i=1}^N\int_0^t\exp\left(\int_0^sg(X_r^i)dr\right)\nabla_xP_{s,T}^gf(X_s^i)^{\top}\sigma(X_s^i)dW_s^i.\end{align*}

In the equality, taking $t=T$ and using $P_{T,T}^gf=f$ given in Lemma 4, we arrive at (95). This completes the proof of the lemma.

Building upon the auxiliary results established by the previous lemmas, we next prove the main result of this section on the propagation of chaos of the FPK equation (5) with respect to the metric defined by (73).

Theorem 2. Let the assumptions $(\textbf{\textrm{A}}_{X_0})$ , ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ , ${({\textbf{A}}_{D,\sigma}\!)}$ , and ${({\textbf{A}}_{gl})}$ (or the assumptions $(\textbf{\textrm{A}}_{X_0})$ , $(\textbf{\textrm{A}}_{1\textrm{d}})$ , ${({\textbf{A}}_{D,\sigma}\!)}$ , and ${({\textbf{A}}_{gl})}$ with ${\textrm{deg}}(Q_g)=0$ ) hold. Let $\mu=(\mu_{t})_{t\in[0,T]}$ be the ${\mathcal{P}}(D)$ -valued solution of the FPK equation (5), and let $\mu^N=(\mu_t^N)_{t\in[0,T]}$ be the ${\mathcal{P}}(D)$ -valued process given by (65) in Section 4.1. Then, for any $T>0$ and $N\geq1$ , there exists a constant $C>0$ which is independent of N such that, for any $p\geq1$ ,

(97) \begin{align}d_{q,T}\left(\mu,\mu^N\right)\leq C\left(\alpha(p,q,n,N)+\frac{1}{N^{q-1}}\right),\quad q\geq2,\end{align}

where the metric $d_{q,T}(\cdot,\cdot)$ is defined by (73), and the first rate of convergence $\alpha(p,q,n,N)$ in (97) is given by

(98) \begin{align}\alpha(p,q,n,N)\,:\!=\,\left\{ \begin{array}{ll} \displaystyle N^{-\frac{1}{2}}+N^{-\frac{p-q}{p}}, & q>\frac{n}{2}, p\neq2q;\\[5pt] \displaystyle N^{-\frac{1}{2}}\ln(1+N)+N^{-\frac{p-q}{p}}, & q=\frac{n}{2}, p\neq 2q;\\[5pt] \displaystyle N^{-\frac{q}{n}}+N^{-\frac{p-q}{p}}, & q<\frac{n}{2}, p\neq\frac{n}{n-q}. \end{array} \right.\end{align}

Proof. Using the equality (95) in Lemma 8, we get, for all $\psi\in{\mathcal{R}}_1$ ,

(99) \begin{align}\big\langle\mu_T^{g,N}-\mu_T^g,\psi_{\star}\big\rangle&= \big\langle\mu_0^{g,N}-\mu_0^g,P_{0,T}^g\psi_{\star}\big\rangle\nonumber\\&\quad+\frac{1}{N}\sum_{i=1}^N\int_0^T\exp\left(\int_0^sg(X_r^i)dr\right)\nabla_xP_{s,T}^g\psi_{\star}(X_s^i)^{\top}\sigma(X_s^i)dW_s^i.\end{align}

We next consider the estimate of the first term of the right-hand side of the equality (99). Note that

(100) \begin{align}\mu_0^{g,N}=\mu_0^N=\frac{1}{N}\sum_{i=1}^N\delta_{X_0^i},\quad \mu_0^g=\rho_0.\end{align}

Let $\phi_0\,:\!=\,P_{0,T}\psi_{\star}$ with $\psi\in{\mathcal{R}}_1$ . It follows from the assumption (A $_{X_0}$ ) that $(\phi_0(X_0^i))_{i\geq1}$ are i.i.d. random variables. Let $\tilde{\rho}_0\,:\!=\,{\mathcal{L}}(\phi_0(X_0^1))$ , i.e., the law of the random variable $\phi_0(X_0^1)$ , and define

\begin{equation*}\tilde\mu^N\,:\!=\,\frac{1}{N}\sum_{i=1}^N\delta_{\phi_0(X_0^i)}.\end{equation*}

Note that, under the assumption (A $_{X_0}$ ), by the Glivenko–Cantelli theorem, the empirical measure $\tilde{\mu}^N$ tends weakly to $\tilde{\rho}_0$ as $N\to\infty$ . On the other hand, by (76) and the fact that $\|\psi_{\star}\|_{\infty}\leq2$ , we have from the assumption ${({\textbf{A}}_{gl})}$ (i) (without loss of generality, we have assumed that $g\leq0$ ) that $\|\phi_0\|_{\infty}\leq2$ . Therefore $\int_{D}|x|^p\tilde\rho_0(dx)=\int_{D}|\phi_0(x)|^p\rho_0(dx)\leq 2^p$ for any $p\geq1$ . Then, by Theorem 1 of Fournier and Guillin [Reference Fournier and Guillin9], there is a constant C depending only on n, p, q such that

(101) \begin{align}\mathbb{E}\left[{\mathcal{W}}_q(\tilde{\mu}_0^N,\tilde{\rho}_0)^q\right]&\leq C\left(\int_{D}|x|^p\tilde{\rho}_0(dx)\right)^{\frac{q}{p}}\alpha(p,q,n,N)\leq 2^qC\alpha(p,q,n,N),\end{align}

where ${\mathcal{W}}_q$ denotes the Wasserstein metric with order $q\geq2$ on the qth-order Wasserstein space ${\mathcal{P}}_q(D)$ (see Villani [Reference Villani26]), and the rate of convergence is given by (98). Moreover, for $I(x)=x$ for $x\in D$ , we obtain from (101) and the Kantorovich–Rubinstein dual formula (see Villani [Reference Villani26]) that, for $q\geq2$ ,

(102) \begin{align}\mathbb{E}\left[\left|\big\langle\mu_0^{g,N}-\mu_0^g,P_{0,T}^g\psi_{\star}\big\rangle\right|^q\right]&=\mathbb{E}\left[\left|\big\langle\mu_0^{g,N}-\mu_0^g,\phi_0\big\rangle\right|^q\right]=\mathbb{E}\left[\left|\big\langle\tilde\mu^N-\tilde\rho_0,I\big\rangle\right|^q\right]\nonumber\\&\leq\mathbb{E}\left[{\mathcal{W}}_1(\tilde\mu^N,\tilde\rho_0)^q\right]\leq\mathbb{E}\left[{\mathcal{W}}_q(\tilde\mu^{N},\tilde\rho_0)^q\right].\end{align}

Then, in view of (102) and (101), there is a constant C depending only on n, p, q such that

(103) \begin{align}\mathbb{E}\left[\left|\big\langle\mu_0^{g,N}-\mu_0^g,P_{0,T}^g\psi_{\star}\big\rangle\right|^q\right]\leq C\alpha(p,q,n,N).\end{align}

On the other hand, we can apply the estimate (79) in Lemma 5 (or in Lemma 6). Then, for the estimate of the second term of the right-hand side of the equality (99), under the assumption (A $_{X_0}$ ), the BDG inequality yields that, for $q\geq2$ ,

(104) \begin{align}&\mathbb{E}\left[\sup_{s\in[0,t]}\left|\frac{1}{N}\sum_{i=1}^N\int_0^s\exp\left(\int_0^rg(X_u^i)du\right)\nabla_xP_{r,t}^g\psi_{\star}(X_r^i)^{\top}\sigma(X_r^i)dW_r^i\right|^q\right]\nonumber\\ &\qquad \leq \frac{C}{N^q}\sum_{i=1}^N\mathbb{E}\left[\left(\int_0^t\exp\left(2\int_0^sg(X_r^i)dr\right)\left|\nabla_xP_{s,t}^g\psi_{\star}(X_s^i)\right|^2\left|\sigma(X_s^i)\right|^2ds\right)^{q/2}\right]\nonumber\end{align}

\begin{align} &\qquad\leq\frac{C}{N^q}\sum_{i=1}^N\mathbb{E}\left[\left(\int_0^t\exp\left(2\int_0^sg(X_r^i)dr\right)\hat{Q}^{(1)}_{g}(|X_s^i|)^2\left|\sigma(X_s^i)\right|^2ds\right)^{q/2}\right]\nonumber\\ &\qquad\leq\frac{C}{N^q}\sum_{i=1}^N\mathbb{E}\left[\left(\int_0^t\hat{Q}^{(1)}_{g}(|X_s^i|)^2\left|\sigma(X_s^i)\right|^2ds\right)^{q/2}\right]\notag\\ &\qquad=\frac{C}{N^{q-1}}\mathbb{E}\left[\left(\int_0^t\hat{Q}^{(1)}_{g}(|X_s^1|)^2\left|\sigma(X_s^1)\right|^2ds\right)^{q/2}\right],\end{align}

for some constant $C=C_{g,b,\sigma,q,T}>0$ which may be different from line to line. Here we have used $g\leq0$ in the assumption ${({\textbf{A}}_{gl})}$ , and $\hat{Q}_g^{(1)}$ is the polynomial with ${\textrm{deg}}(\hat{Q}_g^{(1)})={\textrm{deg}}(Q_g)$ given in Lemma 5. In addition, by the linear growth condition in the assumption ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ or the Hölder continuity given in the assumption (A $_{1\textrm{d}}$ ), there exists a constant $C=C_{g,b,\sigma,T}>0$ such that, $\mathbb{P}$ -almost surely,

(105) \begin{align} \hat{Q}^{(1)}_{g}(|X_s^1|)^2\left|\sigma(X_s^1)\right|^2\leq {\textrm{P}}_1\left(\sup_{t\in[0,T]}|X^1_t|^2\right),\quad s\in[0,T],\end{align}

where ${\textrm{P}}_1\,:\,\mathbb{R}_+\to\mathbb{R}_+$ is a polynomial with ${\textrm{deg}}({\textrm{P}}_1)={\textrm{deg}}(Q_g)+1$ . By (64) under the assumption (A $_{X_0}$ ), it follows that

(106) \begin{align} \mathbb{E}\left[\sup_{t\in[0,T]}\left|X^1_t\right|^{{\textrm{deg}}({\textrm{P}}_1)q}\right]\leq C\left\{\mathbb{E}\left[\left|X^1_0\right|^{{\textrm{deg}}({\textrm{P}}_1)q}\right]+1\right\}<\infty,\end{align}

for some constant $C=C_{g,b,\sigma,T}>0$ . Then, by (104), (105), and (106), there exists a constant $C=C_{g,b,\sigma,T,q}>0$ , which may be different from line to line, such that

(107) \begin{align}&\mathbb{E}\left[\sup_{s\in[0,t]}\left|\frac{1}{N}\sum_{i=1}^N\int_0^s\exp\left(\int_0^rg(X_u^i)du\right)\nabla_xP_{r,t}^g\psi_{\star}(X_r^i)^{\top}\sigma(X_r^i)dW_r^i\right|^q\right]\nonumber\\&\qquad\qquad\leq\frac{C}{N^{q-1}}\mathbb{E}\left[\left(\int_0^t\hat{Q}^2_{g}(|X_s^1|)\left|\sigma(X_s^1)\right|^2ds\right)^{q/2}\right]\nonumber\\&\qquad\qquad\leq\frac{C}{N^{q-1}}\left\{{\mathbb E}\left[\left|X^1_0\right|^{{\textrm{deg}}({\textrm{P}}_1)q}\right]+1\right\}.\end{align}

We deduce from (99), (103), and (107) that

(108) \begin{align} &\mathbb{E}\left[\left|\big\langle\mu_t^{g,N}-\mu_t^g,\psi_{\star}\big\rangle\right|^q\right]\leq2^q\mathbb{E}\left[\left|\big\langle\mu_0^{g,N}-\mu_0^g,P_{0,t}^g\psi_{\star}\big\rangle\right|^q\right]\notag\\ &\qquad\quad+2^q\mathbb{E}\left[\sup_{s\in[0,t]}\left|\frac{1}{N}\sum_{i=1}^N\int_0^s\exp\left(\int_0^rg(X_u^i)du\right)\nabla_xP_{r,t}^g\psi_{\star}(X_r^i)^{\top}\sigma(X_r^i)dW_r^i\right|^q\right]\notag\\&\qquad\leq C\alpha(p,q,n,N)+\frac{C}{N^{q-1}}\left\{{\mathbb E}\left[\left|X^1_0\right|^{{\textrm{deg}}({\textrm{P}}_1)q}\right]+1\right\}.\end{align}

Note that the right-hand side of (108) is independent of $\psi_{\star}$ and t. We can obtain (97) by taking the supremum with respect to $\psi\in{\mathcal{R}}_1$ and $t\in[0,T]$ . This completes the proof of the theorem.

Recall that Example 3 relates to the one-dimensional RM equation with $b(x)=a+bx$ , $\sigma(x)=\sigma\sqrt{x}$ , and $g(x)=-x$ for $x\in D=(0,\infty)$ ; Example 4 relates to the one-dimensional RM equation with $b(x)\equiv0$ , $\sigma(x)=\sigma\sqrt{2}$ , and

\begin{equation*}g(x)=-x^{2q} + \sum_{l=0}^{2q-1}\alpha_lx^l\end{equation*}

given by (50) for $x\in D=(-\infty,\infty)$ ; Example 5 relates to the n-dimensional RM equation with $b(x)=b+ Bx$ , $\sigma(x)\equiv\sigma$ , and $g(x)=-r(x)$ given by (55) for $x\in D=\mathbb{R}^n$ . It is not difficult to verify the validity of the assumptions (A $_{1\textrm{d}}$ ), ${({\textbf{A}}_{D,\sigma})}$ , and ${({\textbf{A}}_{gl})}$ with ${\textrm{deg}}(Q_g)=0$ for Example 3, and the validity of the assumptions ${({\textbf{A}}^{\prime}_{b,\sigma}\!)}$ , ${({\textbf{A}}_{D,\sigma})}$ , and ${({\textbf{A}}_{gl})}$ for Examples 4 and 5. Therefore, the propagation of chaos with the rate of convergence (97) in Theorem 2 under the assumption (A $_{X_0}$ ) on the initial data of the particle system holds for these examples.