A Stein’s Method Approach to the Linear Noise Approximation for Stationary Distributions of Chemical Reaction Networks

Systems composed of a set of chemical species interacting via a finite set of reactions are most often modeled as a Stochastic Chemical Reaction Network (SCRN), a continuous time Markov chain (CTMC) model for how the molecular counts of each species change over time. Such models have found widespread use in systems and synthetic biology as a method to capture the counts of relevant biomolecules within a single cell, where due to the small numbers of molecules present, the effect of stochasticity in the system’s evolution can be profound [67, 61, 52, 46, 6, 20, 65]. The time evolution of the probability distribution of molecular counts can be computed via the Kolmogorov forward equations, often called the Chemical Master Equation [27], or by a Monte Carlo approach using the Doob-Gillespie Algorithm [16, 26]. However, due to the very large, and often countably infinite, number of states in the Markov chain, simpler approximations have been proposed, the most prominent of which exploit the fact that for many physical systems the volume in which the reactions occur and the molecular counts of all species are large. These approximations are formalized by considering a family of Markov chains $\bm{X}_{\Omega}(t)$, where $\bm{X}_{\Omega}(t)$ is the vector of molecular counts at time $t$ and $\Omega$ is the volume. The Reaction Rate Equation (RRE) model is an Ordinary Differential Equation (ODE), whose solution $\bm{v}(t)$ is a deterministic approximation for $\bm{X}_{\Omega}(t)/\Omega$, rigorously justified by Kurtz on finite time intervals in the sense that under mild technical conditions, $\bm{X}_{\Omega}(t)/\Omega\rightarrow_{p}\bm{v}(t)$ uniformly on any finite time interval [39]. Due to $\bm{v}(t)$ being the solution to an ODE, the computational and analytical simplicity of the RRE model has lead to its widespread use in both modeling in systems biology [3] and design in synthetic biology [14, 24, 19]. In essence, the RRE model approximates the concentrations by a constant at each point in time, neglecting any stochastic fluctuation away from this value. In order to account for these fluctuations, Van Kampen expanded the Forward Kolmogorov Equations in powers of $1/\sqrt{\Omega}$ [35, 62]. The first order approximation he derived in this manner, called the Linear Noise Approximation (LNA), adds a correction term to the RREs, specifically a Gaussian approximation to fluctuations of $\bm{X}_{\Omega}(t)/\Omega$ about $\bm{v}(t)$. The LNA was rigorously justified by Kurtz in the sense of giving technical conditions for the central limit theorem result $\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}(t)-\Omega\bm{v}(t)\right)\rightarrow_{p}\mathcal{N}(0,P(t))$ to hold uniformly on finite time intervals with an appropriate covariance matrix $P(t)$. The LNA is similar to the RREs in that the approximating distribution can be easily computed from the solution to an ODE, and due to this simplicity it has been used widely to quantify noise in natural biological systems, [50, 51, 18, 60], as well as for design [56, 55], and inference [38, 22, 23] in synthetic biology.

However, a key weakness of the theoretical justification for the RREs and the LNA is that the results given in [39] do not say anything about the relationship between the stationary distribution of the Markov chain model and the steady state behavior of the RREs or LNA, even when the stationary distribution exists uniquely for all $\Omega$ and the RREs/LNA has a well defined limiting behavior. In a later work, Kurtz gave a sufficient condition for the uniform in time convergence of $\bm{X}_{\Omega}(t)/\Omega$ to $\bm{v}(t)$ and of $\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}(t)-\Omega\bm{v}(t)\right)$ to $\mathcal{N}(0,P(t))$ [40]. However, that result requires that the state space of $\bm{X}_{\Omega}(t)/\Omega$ be bounded uniformly in $\Omega$, which rules out the analysis of many biological systems of interest where the molecular counts of certain species can grow unbounded, and thus for all $\Omega$, $\bm{X}_{\Omega}(t)$ evolves on a countably infinite subset of the integer lattice. This technical difficulty in using the RREs or LNA to approximate the stationary distribution of an SCRN has long been understood [62], and in fact it is often the case that the Markov chain and RREs/LNA models have different long term behavior [43]. In this work we give sufficient conditions under which there exists a constant $C$ such that for all sufficiently large $\Omega$,

$$\mathcal{W}_{1}\left(\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}^{\infty}-\Omega\bm{v}^{*}\right),\bm{Y}^{\infty}\right)\leq C\frac{\ln{\Omega}}{\sqrt{\Omega}},$$

(1)

where $\bm{X}_{\Omega}^{\infty}$ is distributed according to the stationary distribution of $\bm{X}_{\Omega}(t)$, $\bm{v}^{*}$ is an appropriate equilbrium point of the RRE, $\bm{Y}^{\infty}\sim\mathcal{N}(0,P^{\infty})$ with $P^{\infty}$ a covariance matrix computed via the LNA, and $\mathcal{W}_{1}$ is the 1-Wasserstein distance between the laws of the two random variables. This result immediately implies both that the Markov chain’s stationary distribution converges to an equilibrium point of the RRE in the sense $\bm{X}_{\Omega}^{\infty}/\Omega\rightarrow_{d}\bm{v}^{*}$, where $\bm{v}^{*}$ is an appropriate equilbrium of the RRE, and that $\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}(t)-\Omega\bm{v}(t)\right)\rightarrow_{d}\mathcal{N}(0,P^{\infty})$ [25]. We give sufficient conditions in terms of second moment of $\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}^{\infty}-\Omega\bm{v}^{*}\right)$ for (1) to hold, and additionally give sufficient conditions on only the RRE model for (1) to hold, which include global exponential stability of $\bm{v}^{*}$. Interestingly, the convergence rate given by (1) is identical to the convergence rate shown by Kurtz on finite time intervals [41].

The primary technique we use to prove our results is Stein’s Method [13]. Stein’s Method is a powerful technique for proving central limit theorems introduced by Stein in 1972 [59, 58], and used in a wide variety of situations since then, such as when the limiting distribution is Poisson [7] or Binomial [17], as well as for convergence to a diffusion process [8]. The most closely related work to our study comes from queueing theory, where Braverman et al. using Stein’s Method, and Gurvich using a closely related technique, studied diffusion approximations of certain queueing systems in the heavy traffic regime [11, 10, 33]. Our application of Stein’s method follows that of [11] at a high level, where the key difference is that our family of CTMCs and limiting process is in $n$ dimensions instead of one. In order to handle this, we use results on Stein Factors from [30].

Although rigorous work on approximating the stationary distribution of SCRNs in the large volume limit has been limited, a variety of related ideas have been presented in the literature. In particular, when the microscopic rates are all affine, it is known that the LNA matches the first and second moments of $\bm{X}_{\Omega}(t)$ for all time [62, 43]. This idea has been extended to a restricted class of SCRNs with nonlinear propensities in [31]. Additionally, Sontag and Singh gave a special class of SCRNs where the moments can be exactly computed even when they may not match the moments of the LNA exactly [57]. However, in these cases it is still unclear if the stationary distribution of $\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}(t)-\Omega v^{*}\right)$ converges to the Gaussian given by the LNA in the large volume limit. In fact, the results we present allow one to rigorously go from convergence of $\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}^{\infty}-\Omega v^{*}\right)$ to $\mathcal{N}(0,P^{\infty})$ with respect to the first and second moments, to convergence in 1-Wasserstein distance. We can thus strengthen the results of e.g. [31] with limited additional work. We note that methods for establishing convergence of the stationary distributions of a family of Markov Processes to the stationary distribution of a limiting process are given in [21]. However, not only have these methods not been successfully applied to SCRNs, but also our Stein’s Method approach yields non-asymptotic bounds on the approximation error in 1-Wasserstein distance.

In certain special cases, such as detailed or complex balance, the stationary distribution of $\bm{X}_{\Omega}(t)$ can be computed in closed form, and thus in principle one could directly analyze the approximation error when using the LNA [66, 4, 48]. However, for most SCRNs of interest in biology, such results are not applicable. For example, a simple two step transcription-translation model for protein expression is not detailed or complex balanced [14], and thus the results in [66, 4] cannot be used to obtain a closed form expression for the stationary distribution, whereas the results of [48] are only exact for systems where $\bm{X}_{\Omega}(t)$ has finitely many states.

We note that the LNA is only one diffusion approximation for SCRNs, and in fact the Chemical Langevin Equation (CLE) is often a more accurate approximation to $\bm{X}_{\Omega}(t)/\Omega$ for moderately large volumes [28]. The convergence of $\bm{X}_{\Omega}(t)/\Omega$ to the CLE on finite time intervals was established by Kurtz [41], but using the CLE to approximate the stationary distribution of $\bm{X}_{\Omega}(t)/\Omega$ has remained mostly unstudied. In fact, the CLE often exhibits finite time breakdown and thus does not have a well defined stationary distribution, a problem that has been tackled by allowing complex values in the Stochastic Differential equation [54], as well as imposing the physical constraint that concentrations are nonnegative [42].

The rest of the paper is as follows: In Section 2 we give notation and mathematical background, and then introduce the three relevant models. In Section 3 we state our main results, the proofs of which appear in Section 4. In Section 5 we illustrate our results with several examples. In Section 6 we provide concluding remarks and future research directions.

We denote (column) vectors by bold symbols $\bm{x},\bm{Z}$, with the elements indexed as $\bm{x}=[x_{1},\dots,x_{n}]^{T}$. Let $\mathbb{Z}_{\geq 0}^{n}$ [$\mathbb{Z}_{>0}^{n}$] denote the nonnegative [positive] integer lattice in $n$ dimensions. Let $\mathbb{Q}$ denote the set of rational numbers. Let $\mathbb{R}^{n}$, $\mathbb{R}_{\geq 0}^{n}$, and $\mathbb{R}^{n}_{>0}$ denote n-dimensional Euclidian space, and the nonnegative and positive orthants of the same respectively. For $\bm{x}\in\mathbb{R}^{n}$ we denote the standard 2-norm by $\lVert\bm{x}\rVert$, and for a matrix $A\in\mathbb{R}^{n\times m}$ we denote by $\left\lVert A\right\rVert$ the induced 2-norm $\left\lVert A\right\rVert=\sup_{\bm{x}\neq 0}\frac{\left\lVert A\bm{x}\right\rVert}{\lVert\bm{x}\rVert}$, and denote by $\mathrm{cond}(A)$ the condition number of $A$. When all eigenvalues of $A\in\mathbb{R}^{n\times n}$ have strictly negative real parts we say that $A$ is Hurwitz. Given a function $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$, we denote by $\frac{\partial g}{\partial\bm{z}}$, $\nabla g$, and $\nabla^{2}g$ the row vector of partial derivatives, the gradient, and the Hessian of $g$ respectively. We denote the third (directional) derivative of $g$ by $\nabla^{3}g[\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}]$, along the directions $\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}$, and define its induced norm as $\sup_{\lVert\bm{x}_{1}\rVert=\lVert\bm{x}_{2}\rVert=\lVert\bm{x}_{3}\rVert=1}\lVert\nabla^{3}g[\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}]\rVert$. For two functions $f$ and $g$ from $\mathbb{R}^{n}$ to $\mathbb{R}$, we define the convolution of $f$ and $g$ as $(f\star g)(\bm{x})=\int_{\mathbb{R}^{n}}f(\bm{x}-\bm{s})g(\bm{s})\mu(d\bm{s})$, with $\mu$ the Lebesgue measure on $\mathbb{R}^{n}$. We denote the set of k-continuously differentiable functions on a domain $\mathcal{D}$ by $C^{k}(\mathcal{D})$, where we omit the domain if it is clear from context.

We consider Stochastic Chemical Reaction Networks parameterized by $\Omega\in\mathcal{O}\subseteq\mathbb{R}_{>0}$, composed of $n$ species $\mathrm{X}_{1},\mathrm{X}_{2},\dots,\mathrm{X}_{n}$, interacting according to $r$ reactions, each characterized by a nonzero reaction vector $\bm{\xi}_{i}\in\mathbb{Z}^{n}$ and a microscopic propensity function $\lambda_{i}(\bm{x})$, which depends on $\Omega$. The reaction vector describes the change in the counts on the species when reaction $i$ fires, and the microscopic propensity function describes the rate at which reaction $i$ fires. Thus, for each $\Omega\in\mathcal{O}$, the Stochastic Chenical Reaction Network describes a CTMC $\{\bm{X}_{\Omega}(t),t\geq 0\}$, where $\bm{X}_{\Omega}(t)=[X_{1}(t),X_{2}(t),\dots,X_{n}(t)]^{T}\in\mathcal{X}_{\Omega}\subseteq\mathbb{Z}_{\geq 0}^{n}$ represents the molecular counts of the species. The state space $\mathcal{X}_{\Omega}$ can be chosen as any subset of $\mathbb{Z}_{\geq 0}^{n}$ such that 1) $\lambda_{i}(\bm{x})\geq 0$ for all $\bm{x}\in\mathcal{X}_{\Omega}$ and 2) for all $\bm{x}\in\mathcal{X}_{\Omega}$ and all $i\in\{1,\dots,r\}$ such that $\bm{x}+\bm{\xi}_{i}\notin\mathcal{X}_{\Omega}$, $\lambda_{i}(\bm{x})=0$. The infinitesimal transition rates from $\bm{x}$ to $\bm{x}^{\prime}$ are given by $q_{\bm{x},\bm{x}^{\prime}}=\sum_{i:\bm{x}+\bm{\xi}_{i}=\bm{x}^{\prime}}\lambda_{i}(\bm{x})$. We refer the interested reader to [5] for a more thorough background on SCRNs. Throughout this work we will need to impose the following assumption on our SCRNs:

Under Assumption 2.1, we can define $\mathcal{O}=\left\{\Omega\in\mathbb{R}_{>0}\middle|\forall j,\;\Omega x_{j}^{max}\in\mathbb{Z}_{>0}\right\}$. Note that our assumption that $x_{j}^{max}\in\mathbb{Q}$ ensures that $\sup\mathcal{O}=+\infty$, i.e. our family of Markov chains includes those with arbitrarily large volume. Under Assumption 2.1 we will always choose to define the state space of our CTMC as $\mathcal{X}_{\Omega}=\mathbb{Z}_{\geq 0}^{n}\cap\Omega\mathcal{X}$. In order to apply our Stein’s Method technique, we will need the following assumption:

One form of $\lambda_{i}(\bm{x})$ which satisfies Assumptions 2.1 and 2.2 is mass action kinetics [62], which means that the propensities are of the form

$$\lambda_{i}(\bm{x})=k_{i}\Omega^{1-\sum_{j}\xi_{ri,j}}\prod_{j=1}^{n}x_{j}^{\xi_{ri,j}},$$

(3)

where $\bm{\xi}_{i}=\bm{\xi}_{pi}-\bm{\xi}_{ri}$ with $\bm{\xi}_{pi},\bm{\xi}_{ri}\in\mathbb{Z}_{\geq 0}^{n}$ is the reaction vector of the chemical reaction

$$\schemestart\subscheme{\sum_{j}\xi_{ri,j}\mathrm{X}_{j}}\arrow(x1--x1){->[k_{i}]}[0]\subscheme{\sum_{j}\xi_{pi,j}\mathrm{X}_{j}}\schemestop.$$

(4)

The requirement that $\sum\bm{\xi}_{pi}\leq 2$ must be imposed for $\lambda_{i}(\bm{x})$ to satisfy Assumption 2.2. See [27] for a discussion of the physical assumptions necessary to reach such propensities. We note that (3) requires that no reactant appears twice on the left hand side of the reaction. Having a reactant $\mathrm{X}_{j}$ appear twice on the left hand side of a reaction creates a propensity proportional to $x_{j}^{2}$ instead of the physically accurate $x_{j}(x_{j}-1)$ [27]. Here $k_{i}$ is the reaction rate constant of reaction $i$ and $\xi_{ri,j}$ is the $j$^th element of $\bm{\xi}_{ri}$. Another form of $\lambda_{i}(\bm{x})$ that satisfies Assumptions 2.1 and 2.2 is

$$\lambda_{i}(\bm{x})=\Omega\frac{k^{i}_{1}(x_{j_{i}}/\Omega)^{p}}{k^{i}_{2}+(x_{j_{i}}/\Omega)^{q}},$$

(5)

where $j_{i}\in\{1,\dots,n\}$ and $p,q\in\mathbb{Z}_{>0}$ such that $p+q\leq 3$ or $p=q$. Such “Hill-type propensities” can arise through timescale separation in systems of chemical reactions with mass action kinetics [36, 47, 34]. For simplicity, we will impose the following assumption:

Though Assumption 2.3 appears to rule out SCRNs that have conservation laws, i.e. a linear combination of species counts that remains unchanged by all reactions, we will show through an example in Section 5 that Assumption 2.3 does not fundamentally prevent the analysis of SCRNs with conservation laws, as long as a change of coordinates exist where the system can be represented in a lower dimensional lattice $\mathbb{Z}_{\geq 0}^{n^{\prime}}$ as an SCRN with stoichiometric subspace of dimension $n^{\prime}$. For $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ the generator of $\bm{X}_{\Omega}(t)$ is

$$G_{\bm{X}_{\Omega}}f(\bm{x})=\sum_{i=1}^{r}\lambda_{i}(\bm{x})\left(f(\bm{x}+\bm{\xi}_{i})-f(\bm{x})\right).$$

(6)

We will also consider two shifted and scaled version of $\bm{X}_{\Omega}(t)$ defined by

$$\tilde{\bm{X}}_{\Omega}(t)=\frac{1}{\sqrt{\Omega}}\left(\bm{X}_{\Omega}(t)-\Omega\bm{v}^{*}\right),$$

(7)

and

$$\tilde{\tilde{\bm{X}}}_{\Omega}(t)=\frac{1}{\Omega}\bm{X}_{\Omega}(t)-\bm{v}^{*},$$

(8)

where $\bm{v}^{*}$ is a specified point in $\mathcal{X}$. For $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ the generator of $\tilde{\bm{X}}_{\Omega}(t)$ is given by

$$G_{\tilde{\bm{X}}_{\Omega}}f(\bm{x})=\sum_{i=1}^{r}\lambda_{i}(\sqrt{\Omega}\bm{x}+\Omega\bm{v}^{*})\left(f(\bm{x}+\frac{1}{\sqrt{\Omega}}\bm{\xi}_{i})-f(\bm{x})\right),$$

(9)

and the generator of $\tilde{\tilde{\bm{X}}}_{\Omega}(t)$ is given by

$$G_{\tilde{\tilde{\bm{X}}}_{\Omega}}f(\bm{x})=\sum_{i=1}^{r}\Omega\bar{\lambda}_{i}(\bm{x}+\bm{v}^{*})\left(f(\bm{x}+\frac{1}{\Omega}\bm{\xi}_{i})-f(\bm{x})\right).$$

(10)

When $\bm{X}_{\Omega}(t)$ has a unique stationary distribution $\pi$, we denote by $\bm{X}_{\Omega}^{\infty}$ a random variable having law $\pi$. We likewise define $\tilde{\bm{X}}_{\Omega}^{\infty}$ and $\tilde{\tilde{\bm{X}}}_{\Omega}^{\infty}$.

The Reaction Rate Equations (RREs) are an Ordinary Differential Equation Model which is used to approximate $\bm{X}_{\Omega}(t)/\Omega$. Defining $\bar{\bm{\lambda}}(\bm{v})=\begin{bmatrix}\bar{\lambda}_{1}(\bm{v}),&\bar{\lambda}_{2}(\bm{v}),&\dots,&\bar{\lambda}_{r}(\bm{v})\end{bmatrix}^{T}$ and $S=\begin{bmatrix}\bm{\xi}_{1}&\bm{\xi}_{2}&\dots&\bm{\xi}_{r}\end{bmatrix}$, the reaction rate equations are

$\displaystyle\frac{d}{dt}\bm{v}(t)$

$\displaystyle=F(\bm{v}(t)),$

$\displaystyle\bm{v}(0)$

$\displaystyle=\bm{v}_{0}$

(11)

where $F(\bm{v})=S\bar{\bm{\lambda}}(\bm{v})$ and $\bm{v}\in\mathcal{X}\subseteq\mathbb{R}^{n}_{\geq 0}$. An equilibrium point of the RREs is any point $\bm{v}^{*}\in\mathcal{X}$ satisfying $0=F(\bm{v}^{*})$. We note that under Assumption 2.1, $\mathcal{X}$ will be positive invariant with respect to 11.

The Linear Noise Approximation (LNA) is a diffusion approximation obtained by expanding the Chemical Master Equation using the ansatz $\bm{X}_{\Omega}(t)\approx\Omega\bm{v}(t)+\sqrt{\Omega}\bm{Y}^{\prime}(t)$ [62]. The terms $\bm{v}(t)$ and $\bm{Y}^{\prime}(t)$ are given by

	$\displaystyle\frac{d}{dt}\bm{v}(t)$	$\displaystyle=F(\bm{v}(t)),$	$\displaystyle\bm{v}(0)$	$\displaystyle=\bm{v}_{0},$		(12a)
	$\displaystyle d\bm{Y}^{\prime}(t)$	$\displaystyle=\frac{\partial F(\bm{v})}{\partial\bm{v}}\bm{Y}^{\prime}(t)dt+S\operatorname{diag}\sqrt{\bar{\bm{\lambda}}(\bm{v})}d\bm{W}(t),$	$\displaystyle\bm{Y}^{\prime}(0)$	$\displaystyle=\bm{Y}^{\prime}_{0},$		(12b)

where evidently (12a) is the RRE (11). Here $\bm{W}(t)$ is a unit covariance Wiener Process. Let $\bm{v}^{*}$ be an equilibrium point of (12a). We note that $\bar{\lambda}_{i}(\bm{v})\geq 0$ must hold for (12b) to make sense, which is guaranteed by e.g. Assumption 2.1. We denote by $\bm{Y}(t)$ the solution to (12b) with $\bm{v}(t)=\bm{v}^{*}$ and $\bm{Y}^{\prime}(0)=\bm{Y}_{0}$. The generator of $\bm{Y}(t)$ is

$$G_{\bm{Y}}f(\bm{x})=\frac{\partial f(\bm{x})}{\partial\bm{x}}\cdot\frac{\partial F(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+\frac{1}{2}\sum_{i=1,j=1}^{n}D_{ij}\frac{\partial^{2}f(\bm{x})}{\partial x_{i}\partial x_{j}},$$

(13)

where $D=S\operatorname{diag}\left(\bar{\bm{\lambda}}(\bm{v}^{*})\right)S^{T}$. We remark that $\bm{Y}(t)$ is an Ornstein–Uhlenbeck process, and thus $\bm{Y}(t)$ is Gaussian as long as $\bm{Y}_{0}$ is. The stationary distribution of $\bm{Y}(t)$ is a zero mean Gaussian with covariance matrix $P^{\infty}\in\mathbb{R}^{n\times n}$ given as the solution to the Lyapunov Equation

$$\frac{\partial F(\bm{v}^{*})}{\partial\bm{v}}P^{\infty}+P^{\infty}\frac{\partial F(\bm{v}^{*})}{\partial\bm{v}}^{T}=-S\operatorname{diag}\bar{\bm{\lambda}}(\bm{v}^{*})S^{T},$$

(14)

which has a unique positive definite solution as long as $\frac{\partial F(\bm{v}^{*})}{\partial\bm{v}}$ is Hurwitz. In this case, we denote by $\bm{Y}^{\infty}$ a random variable distributed according to $\mathcal{N}(0,P^{\infty})$.

In this work we measure the distance between probability distributions using the 1-Wasserstein distance [25]. By a slight abuse of notation, we define the 1-Wasserstein distance between two random variables as the 1-Wasserstein distance between the laws of the random variables, formally:

We remark that this definition of 1-Wasserstein distance uses the dual formulation [63], and is a slight abuse of notation since the 1-Wasserstein distance is usually defined between two probability distributions. Our definition is equivalent to taking the 1-Wasserstein distance between the laws of $\bm{Z}_{1}$ and $\bm{Z}_{2}$.

We require the following Lemma from [29], which gives sufficient conditions for the the Basic Adjoint Relationship (BAR) to hold.

It will be convenient for us to consider the 1-Wasserstein distance computed by taking a supremum over not the 1-Lipschitz functions, but the 1-Lipschitz functions that are additionally in $C^{3}$ with bounded derivatives. To this end, let

$$\mathcal{H}=\{h\in C^{3}(\mathbb{R}^{n})|h\in\mathrm{Lip}(1),\;\forall k\in\{1,2,3\},\;\sup_{\bm{x}\in\mathbb{R}^{n}}\lVert\nabla^{k}h(\bm{x})\rVert<\infty\}.$$

(36)

We have the following Lemma, which is proved in Section 4.3.

We now present the derivative bounds that we will need. The proof of Lemma 4.3 uses results from [30], and is provided in 4.2.

Our goal is to analyze

$$\mathcal{W}_{1}\left(\tilde{\bm{X}}_{\Omega}^{\infty},\bm{Y}^{\infty}\right)=\sup_{h\in\mathrm{Lip}(1)}\left\lvert\mathbb{E}\left[h\left(\tilde{\bm{X}}_{\Omega}^{\infty}\right)\right]-\mathbb{E}\left[h\left(\bm{Y}^{\infty}\right)\right]\right\rvert.$$

(48)

We begin by observing that per Lemma 4.2 we can obtain $\mathcal{W}_{1}\left(\tilde{\bm{X}}_{\Omega}^{\infty},\bm{Y}^{\infty}\right)$ by taking the supremum over $h\in\mathcal{H}$ instead of $h\in\mathrm{Lip}(1)$, and thus we can analyze

$$\mathcal{W}_{1}\left(\tilde{\bm{X}}_{\Omega}^{\infty},\bm{Y}^{\infty}\right)=\sup_{h\in\mathcal{H}}\left\lvert\mathbb{E}\left[h\left(\tilde{\bm{X}}_{\Omega}^{\infty}\right)\right]-\mathbb{E}\left[h\left(\bm{Y}^{\infty}\right)\right]\right\rvert.$$

(49)

To do this we use Stein’s Method following [11]. The Stein Equation is

$$G_{\bm{Y}}f_{h}(\bm{x})=h(\bm{x})-\mathbb{E}\left[h(\bm{Y}^{\infty})\right],$$

(50)

which in our case is the Poisson equation

$$\frac{\partial f_{h}(\bm{x})}{\partial\bm{x}}\cdot\frac{\partial F(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+\frac{1}{2}\sum_{i=1,j=1}^{n}D_{ij}\frac{\partial^{2}f_{h}(\bm{x})}{\partial x_{i}\partial x_{j}}=h(\bm{x})-\mathbb{E}\left[h(\bm{Y}^{\infty})\right],$$

(51)

which is guaranteed to have a solution by Lemma 4.3. Throughout the rest of this proof, we will denote by $f_{h}$ a solution to (51) with $h\in\mathcal{H}$ that satisfies the bounds given by Lemma 4.3. We take the expectation of the Stein Equation with respect to $\tilde{\bm{X}}_{\Omega}^{\infty}$ to find

$$\mathbb{E}\left[G_{\bm{Y}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})\right]=\mathbb{E}\left[h(\tilde{\bm{X}}_{\Omega}^{\infty})\right]-\mathbb{E}\left[h(\bm{Y}^{\infty})\right].$$

(52)

Before proceeding, we show that $\mathbb{E}\left[G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})\right]=0$ for all $h\in\mathcal{H}$. We have that $G_{\tilde{\bm{X}}_{\Omega}}(\bm{x},\bm{x})=-\sum_{i=1}^{r}\lambda_{i}(\sqrt{\Omega}\bm{x}+\Omega\bm{v}^{*})$, and so by Condition 3.1, there exists $K>0$ such that

$$\sum_{\bm{x}}\mathbb{P}[\tilde{\bm{X}}_{\Omega}^{\infty}=\bm{x}]\lvert G_{\tilde{\bm{X}}_{\Omega}}(\bm{x},\bm{x})f_{h}(\bm{x})\rvert\leq\sum_{\bm{x}}\mathbb{P}[\tilde{\bm{X}}_{\Omega}^{\infty}=\bm{x}](K+\lVert\sqrt{\Omega}\bm{x}+\Omega\bm{v}^{*}\rVert^{\kappa})\lvert f_{h}(\bm{x})\rvert,$$

(53)

where the sums are taken over $\bm{x}$ such that $\sqrt{\Omega}\bm{x}+\Omega\bm{v}^{*}\in\mathcal{X}_{\Omega}$. By Lemma 4.3, $f_{h}$ is Lipschitz, and so $\lvert f_{h}(\bm{x})\rvert\leq\lvert f_{h}(0)\rvert+\lVert\bm{x}\rVert$. Therefore,

	$\displaystyle\sum_{\bm{x}}\mathbb{P}[\tilde{\bm{X}}_{\Omega}^{\infty}=\bm{x}]\lvert G_{\tilde{\bm{X}}_{\Omega}}(\bm{x},\bm{x})f_{h}(\bm{x})\rvert$	$\displaystyle\leq\sum_{\bm{x}}\mathbb{P}[\tilde{\bm{X}}_{\Omega}^{\infty}=\bm{x}](K+\lVert\sqrt{\Omega}\bm{x}+\Omega\bm{v}^{*}\rVert^{\kappa})\left(\lvert f_{h}(0)\rvert+\lVert\bm{x}\rVert\right),$		(54)
	$\displaystyle\begin{split}&\leq\sum_{\bm{x}}\mathbb{P}[\tilde{\bm{X}}_{\Omega}^{\infty}=\bm{x}]\left((K\lvert f_{h}(0)\rvert+K\lVert\bm{x}\rVert\vphantom{\max\left\{1,\left(\lVert\sqrt{\Omega}\bm{x}\rVert+\lVert\Omega\bm{v}^{*}\rVert\right)^{\kappa}\right\}}\right.\\ &\left.\quad\quad\quad+\lvert f_{h}(0)\rvert\lVert\sqrt{\Omega}\bm{x}+\Omega\bm{v}^{*}\rVert^{\kappa}\right.\\ &\quad\quad\quad\left.+\max\left\{1,\left(\lVert\sqrt{\Omega}\bm{x}\rVert+\lVert\Omega\bm{v}^{*}\rVert\right)^{\kappa}\right\}\lVert\bm{x}\rVert\right),\end{split}$			(55)

which is finite due to Condition 3.1. Therefore, from Lemma 4.1, $\mathbb{E}\left[G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})\right]=0$. Thus, for all $h\in\mathcal{H}$, (52) implies that

	$\displaystyle\left\lvert\mathbb{E}\left[h\left(\tilde{\bm{X}}_{\Omega}^{\infty}\right)\right]-\mathbb{E}\left[h\left(\bm{Y}^{\infty}\right)\right]\right\rvert$	$\displaystyle=\left\lvert\mathbb{E}\left[G_{\bm{Y}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})\right]-\mathbb{E}\left[G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})\right]\right\rvert,$		(56)
		$\displaystyle=\left\lvert\mathbb{E}\left[G_{\bm{Y}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})-G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})\right]\right\rvert,$		(57)
		$\displaystyle\leq\mathbb{E}\left\lvert G_{\bm{Y}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})-G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\tilde{\bm{X}}_{\Omega}^{\infty})\right\rvert,$		(58)

where $f_{h}$ is a solution to the Stein equation. The existence of $f_{h}$ is guaranteed by Lemma 4.3. For conciseness, let us define $\delta=\frac{1}{\sqrt{\Omega}}$. Let $\bm{l}=\bm{x}/\delta+\bm{v}^{*}/\delta^{2}$. For a fixed $\bm{l}$ we take the first order Taylor expansion of $G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\bm{x})$ in $\delta$ using the Lagrange form of the remainder.

	$\displaystyle G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\bm{x})$	$\displaystyle=\sum_{i=1}^{r}\lambda_{i}(\bm{l})\left(f(\bm{x}+\delta\bm{\xi}_{i})-f(\bm{x})\right),$		(59)
		$\displaystyle=\delta\sum_{i=1}^{r}\lambda_{i}(\bm{l})\frac{\partial f_{h}(\bm{x})}{\partial\bm{x}}\bm{\xi}_{i}+\delta^{2}\frac{1}{2}\sum_{i=1}^{r}\lambda_{i}(l)\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x}+\epsilon\bm{\xi}_{i})\bm{\xi}_{i},$		(60)
	$\displaystyle\begin{split}&=\delta\sum_{i=1}^{r}\lambda_{i}(\bm{l})\frac{\partial f_{h}(\bm{x})}{\partial\bm{x}}\bm{\xi}_{i}+\delta^{2}\frac{1}{2}\sum_{i=1}^{r}\lambda_{i}(l)\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\\ &\quad\quad+\delta^{2}\frac{1}{2}\sum_{i=1}^{r}\lambda_{i}(\bm{l})\left(\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x}+\epsilon\bm{\xi}_{i})\bm{\xi}_{i}-\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right),\end{split}$			(61)

where $0\leq\epsilon\leq\delta$, and $\epsilon$ can depend on $\bm{l}$. We have that $\lambda_{i}(\bm{l})=\frac{1}{\delta^{2}}\bar{\lambda}_{i}(\delta\bm{x}+\bm{v}^{*})$, and so from Taylor’s Theorem we have

$$\frac{1}{\delta^{2}}\bar{\lambda}_{i}(\delta\bm{x}+\bm{v}^{*})=\frac{1}{\delta^{2}}\bar{\lambda}_{i}(\bm{v}^{*})+\frac{1}{\delta}\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+R_{i}(\bm{x}),$$

(62)

where $\lvert R_{i}(\bm{x})\rvert\leq\frac{1}{2}\lVert\bm{x}\rVert^{2}\sup_{\bm{z}\in\mathcal{X}}\lVert\nabla^{2}{\bar{\lambda}_{i}}(\bm{z})\rVert$, in which the supremum over $\mathcal{X}$ is finite due to our assumption that $\bar{\lambda}_{i}(\bm{x})$ has bounded second derivative. We have that (LABEL:eq:taylor_expanded) becomes

$\displaystyle\begin{split}G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\bm{x})&=\sum_{i=1}^{r}\left(\frac{1}{\delta^{2}}\bar{\lambda}_{i}(\bm{v}^{*})+\frac{1}{\delta}\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+R_{i}(\bm{x})\right)\delta\frac{\partial f_{h}(\bm{x})}{\partial\bm{x}}\bm{\xi}_{i}\\ &+\sum_{i=1}^{r}\left(\frac{1}{\delta^{2}}\bar{\lambda}_{i}(\bm{v}^{*})+\frac{1}{\delta}\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+R_{i}(\bm{x})\right)\delta^{2}\frac{1}{2}\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\\ &+\sum_{i=1}^{r}\left(\frac{1}{\delta^{2}}\bar{\lambda}_{i}(\bm{v}^{*})+\frac{1}{\delta}\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}\right)\delta^{2}\frac{1}{2}\left(\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x}+\epsilon\bm{\xi}_{i})\bm{\xi}_{i}-\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right)\\ &+\sum_{i=1}^{r}R_{i}(\bm{x})\delta^{2}\frac{1}{2}\left(\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x}+\epsilon\bm{\xi}_{i})\bm{\xi}_{i}-\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right).\end{split}$

(63)

Collecting terms and using the fact that $F(\bm{v}^{*})=0$, we obtain

$\displaystyle\begin{split}&G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\bm{x})=\frac{\partial f_{h}(\bm{x})}{\partial\bm{x}}\sum_{i=1}^{r}\bm{\xi}_{i}\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+\frac{1}{2}\sum_{i=1}^{r}\bar{\lambda}_{i}(\bm{v}^{*})\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\\ &\quad\quad+\delta\frac{1}{2}\sum_{i=1}^{r}\left[2R_{i}(\bm{x})\frac{\partial f_{h}(\bm{x})}{\partial\bm{x}}\bm{\xi}_{i}+\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right]\\ &\quad\quad+\delta^{2}\frac{1}{2}\sum_{i=1}^{r}R_{i}(\bm{x})\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\\ &\quad\quad+\frac{1}{2}\sum_{i=1}^{r}\left(\bar{\lambda}_{i}(\bm{v}^{*})+\delta\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+\delta^{2}R_{i}(\bm{x})\right)\left(\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x}+\epsilon\bm{\xi}_{i})\bm{\xi}_{i}-\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right).\end{split}$

(64)

Observe that the first two terms are identical to $G_{Y}f_{h}(\bm{x})$, since

$$\frac{1}{2}\sum_{i=1}^{r}\bar{\lambda}_{i}(\bm{v}^{*})\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}=\frac{1}{2}\sum_{p=1,q=1}^{n}D_{pq}\frac{\partial f(\bm{x})}{\partial x_{p}\partial x_{q}},$$

(65)

with $D=S\operatorname{diag}(\bar{\bm{\lambda}}(\bm{v}^{*}))S^{T}$. Therefore, for all $\bm{x}$ such that $\sqrt{\Omega}\bm{x}+\Omega\bm{v}^{*}\in\mathcal{X}_{\Omega}$,

$\displaystyle\begin{split}&\left\lvert G_{\bm{Y}}f_{h}(\bm{x})-G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\bm{x})\right\rvert\leq\left\lvert\frac{1}{2}\delta\sum_{i=1}^{r}\left[2R_{i}(\bm{x})\frac{\partial f_{h}(\bm{x})}{\partial\bm{x}}\bm{\xi}_{i}+\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right]\right\rvert\\ &\quad\quad\quad+\left\lvert\frac{1}{2}\delta^{2}\sum_{i=1}^{r}R_{i}(\bm{x})\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right\rvert\\ &\quad\quad+\left\lvert\frac{1}{2}\sum_{i=1}^{r}\left(\bar{\lambda}_{i}(\bm{v}^{*})+\delta\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\bm{x}+\delta^{2}R_{i}(\bm{x})\right)\left(\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x}+\epsilon\bm{\xi}_{i})\bm{\xi}_{i}-\bm{\xi}_{i}^{T}\nabla^{2}{f_{h}}(\bm{x})\bm{\xi}_{i}\right)\right\rvert.\end{split}$

(66)

Letting $Q_{i}=\sup_{\bm{v}\in\mathcal{X}}\lVert\nabla^{2}{\bar{\lambda}_{i}}(\bm{v})\rVert$ and invoking Lemma 4.3 by the fact that Assumptions 2.1 and 2.3 ensure that $S\operatorname{diag}\sqrt{\bar{\bm{\lambda}}(\bm{v}^{*})}$ has a right inverse since $\bar{\bm{\lambda}}(\bm{v}^{*})>0$ and $S$ has full column rank, we have

$\displaystyle\begin{split}\left\lvert G_{\bm{Y}}f_{h}(\bm{x})-G_{\tilde{\bm{X}}_{\Omega}}f_{h}(\bm{x})\right\rvert&\leq\frac{1}{2}\delta\sum_{i=1}^{r}\left(C_{1}Q_{i}\left\lVert\bm{x}\right\rVert^{2}\lVert\bm{\xi}_{i}\rVert+C_{2}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert\left\lVert\bm{x}\right\rVert\lVert\bm{\xi}_{i}\rVert^{2}\right)\\ &+\frac{1}{4}\delta^{2}C_{2}\left\lVert\bm{x}\right\rVert^{2}\sum_{i=1}^{r}Q_{i}\left\lVert\bm{\xi}_{i}\right\rVert^{2}\\ &+\frac{1}{2}\delta^{2-\zeta}C_{3}(\zeta)\lVert\bm{x}\rVert\sum_{i=1}^{r}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}\\ &+\frac{1}{4}\delta^{3-\zeta}C_{3}(\zeta)\left\lVert\bm{x}\right\rVert^{2}\sum_{i=1}^{r}Q_{i}\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}\\ &+\frac{1}{2}\delta^{1-\zeta}C_{3}(\zeta)\sum_{i=1}^{r}\left\lvert\bar{\lambda}_{i}(\bm{v}^{*})\right\rvert\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}.\end{split}$

(67)

This bound holds for any $h\in\mathcal{H}$, with $f_{h}$ chosen as the solution to the Stein equation (51) specified by Lemma 4.3, and thus from (56) we have

$$\begin{split}\left\lvert\mathbb{E}\left[h\left(\tilde{\bm{X}}_{\Omega}^{\infty}\right)\right]-\mathbb{E}\left[h\left(\bm{Y}^{\infty}\right)\right]\right\rvert&\\ &\leq\mathbb{E}\left[\frac{1}{2}\delta\sum_{i=1}^{r}\left(C_{1}Q_{i}\left\lVert\tilde{\bm{X}}_{\Omega}^{\infty}\right\rVert^{2}\lVert\bm{\xi}_{i}\rVert+C_{2}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert\left\lVert\tilde{\bm{X}}_{\Omega}^{\infty}\right\rVert\lVert\bm{\xi}_{i}\rVert^{2}\right)\right.\\ +\frac{1}{4}\delta^{2}C_{2}\left\lVert\tilde{\bm{X}}_{\Omega}^{\infty}\right\rVert^{2}\sum_{i=1}^{r}Q_{i}\left\lVert\bm{\xi}_{i}\right\rVert^{2}\\ +\frac{1}{2}\delta^{2-\zeta}C_{3}(\zeta)\lVert\tilde{\bm{X}}_{\Omega}^{\infty}\rVert\sum_{i=1}^{r}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}\\ +\frac{1}{4}\delta^{3-\zeta}C_{3}(\zeta)\left\lVert\tilde{\bm{X}}_{\Omega}^{\infty}\right\rVert^{2}\sum_{i=1}^{r}Q_{i}\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}\\ \left.+\frac{1}{2}\delta^{1-\zeta}C_{3}(\zeta)\sum_{i=1}^{r}\left\lvert\bar{\lambda}_{i}(\bm{v}^{*})\right\rvert\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}\right].\end{split}$$

(68)

Using the moment bounds and the linearity of the expectation we have

$$\begin{split}\left\lvert\mathbb{E}\left[h\left(\tilde{\bm{X}}_{\Omega}^{\infty}\right)\right]-\mathbb{E}\left[h\left(\bm{Y}^{\infty}\right)\right]\right\rvert\leq\frac{1}{2}\delta\sum_{i=1}^{r}\left(C_{1}Q_{i}M_{2}\lVert\bm{\xi}_{i}\rVert+C_{2}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert M_{1}\lVert\bm{\xi}_{i}\rVert^{2}\right)\\ +\frac{1}{4}\delta^{2}C_{2}M_{2}\sum_{i=1}^{r}Q_{i}\left\lVert\bm{\xi}_{i}\right\rVert^{2}\\ +\frac{1}{2}\delta^{2-\zeta}C_{3}(\zeta)M_{1}\sum_{i=1}^{r}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}\\ +\frac{1}{4}\delta^{3-\zeta}C_{3}(\zeta)M_{2}\sum_{i=1}^{r}Q_{i}\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}\\ +\frac{1}{2}\delta^{1-\zeta}C_{3}(\zeta)\sum_{i=1}^{r}\left\lvert\bar{\lambda}_{i}(\bm{v}^{*})\right\rvert\left\lVert\bm{\xi}_{i}\right\rVert^{3-\zeta}.\end{split}$$

(69)

To prove the error bound proportional to $\frac{\ln\Omega}{\sqrt{\Omega}}$, we take $\zeta=-(\ln\delta)^{-1}$ for $\Omega>e^{2}$, in which case $C_{3}(\zeta)=-2\ln(\delta)C_{3}^{\prime}+e^{-\phi}C_{3}^{\prime}C_{1}$. Let $\iota>0$. For all $\Omega\geq\max\{\Omega_{0},e^{2}+\iota\}$ we have

$$\begin{split}\left\lvert\mathbb{E}\left[h\left(\tilde{\bm{X}}_{\Omega}^{\infty}\right)\right]-\mathbb{E}\left[h\left(\bm{Y}^{\infty}\right)\right]\right\rvert\leq\frac{1}{2}\delta\sum_{i=1}^{r}\left(C_{1}Q_{i}M_{2}\lVert\bm{\xi}_{i}\rVert+C_{2}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert M_{1}\lVert\bm{\xi}_{i}\rVert^{2}\right)\\ +\frac{1}{4}\delta^{2}C_{2}M_{2}\sum_{i=1}^{r}Q_{i}\left\lVert\bm{\xi}_{i}\right\rVert^{2}\\ +\frac{1}{2}\delta^{2}C_{3}^{\prime}e\left(2\ln(\delta^{-1})+e^{-\phi}C_{1}\right)M_{1}\sum_{i=1}^{r}\left\lVert\frac{\partial\bar{\lambda}_{i}(\bm{v}^{*})}{\partial\bm{v}}\right\rVert\max\{\left\lVert\bm{\xi}_{i}\right\rVert^{2},\left\lVert\bm{\xi}_{i}\right\rVert^{3}\}\\ +\frac{1}{4}\delta^{3}C_{3}^{\prime}e\left(2\ln(\delta^{-1})+e^{-\phi}C_{1}\right)M_{2}\sum_{i=1}^{r}Q_{i}\max\{\left\lVert\bm{\xi}_{i}\right\rVert^{2},\left\lVert\bm{\xi}_{i}\right\rVert^{3}\}\\ +\frac{1}{2}\delta C_{3}^{\prime}e\left(2\ln(\delta^{-1})+e^{-\phi}C_{1}\right)\sum_{i=1}^{r}\left\lvert\bar{\lambda}_{i}(\bm{v}^{*})\right\rvert\max\{\left\lVert\bm{\xi}_{i}\right\rVert^{2},\left\lVert\bm{\xi}_{i}\right\rVert^{3}\}.\end{split}$$

(70)

The $\delta\ln(\delta^{-1})=\frac{1}{2}\delta\ln(\delta^{-2})=\frac{1}{2}\frac{\ln\Omega}{\sqrt{\Omega}}$ term decays the slowest, which yields the desired result.