text:
\section{Introduction}

The motion of flagellated bacteria can be modeled as a sequence of run phases, during which a bacterium moves in a straight line at constant speed. Between two run phases, the bacterium changes direction in a tumble phase, which takes much less time than the run phase and acts as a reorientation. To bias movement towards regions with high concentration of chemoattractant, the bacterium adjusts its turning rate by increasing, resp.~decreasing, the probability of tumbling when moving in an unfavorable, resp.~favorable, direction \cite{Alt:1980p8992,Stock:1999p8984}. Since many species are unable to sense chemoattractant gradients reliably due to their small size, this adjustment is often done via an intracellular mechanism that allows the bacterium to retain information on the history of the chemoattractant concentrations along its path \cite{Bren:2000p7499}. The resulting model, which will be called the ``internal state''  or ``fine-scale'' model in this text, can be formulated as a velocity-jump process, combined with an ordinary differential equation (ODE) that describes the evolution of an internal state that incorporates this memory effect \cite{Erban:2005p4247}.

The probability density distribution of the velocity-jump process evolves according to a kinetic equation, in which the internal variables appear as additional dimensions.  A direct deterministic simulation of this equation is therefore prohibitively expensive, and one needs to resort to a stochastic particle method. Unfortunately, a direct fine-scale simulation using stochastic particles presents a large statistical variance, even in the diffusive asymptotic regime, where the behavior of the bacterial density is known explicitly to satisfy a Keller-Segel advection-diffusion equation. Consequently, it is difficult to assess accurately how the solutions of the fine-scale model differ from their advection-diffusion limit in intermediate regimes. We refer to \cite{Horst1,Horst2,KelSeg70} for numerous historical references on the Keller-Segel equation, to \cite{ErbOthm04,Erban:2005p4247} for formal derivations (based on moment closures) of convergence of the the velocity-jump process with internal state to the Keller-Segel equation, and to the companion paper \cite{limits} for a proof of this convergence based on probabilistic arguments.

In this paper, we propose and analyze a numerical method to simulate  individual-based models for chemotaxis of bacteria with internal dynamics with reduced variance. The variance reduction is based on a coupling technique (control variate): the main idea is to simultaneously simulate, using the same random numbers, a simpler, ``coarse'' process where the internal dynamics is replaced by a direct ``gradient sensing'' mechanism (see \cite{Alt:1980p7984,Othmer:1988p7986,Patlak:1953p7738} for references on such gradient sensing models). The probability density of the latter satisfies a kinetic equation without the additional dimensions of the internal state, and converges to a similar advection-diffusion limit, see e.g.~\cite{Chalub:2004p7641,Hillen:2000p4751,Othmer:2002p4752,limits}.

More precisely, we consider the coupling of two velocity-jump processes with exactly the same advection-diffusion limit (for both processes, a proof of the latter convergence has been obtained in the companion paper~\cite{limits}). The first one is the process with internal dynamics described by the system of equations~\eqref{eq:process_noscale}, with jump rates satisfying \eqref{eq:lin_rate} and internal dynamics satisfying the two main assumptions \eqref{eq:tau_eps_def} and \eqref{eq:tau_scales} as well as other technical assumptions contained in Section~\ref{sec:asympt}. The second one is the process with direct gradient sensing~\eqref{eq:cprocess}, with jump rates satisfying~\eqref{eq:control_rate} with~\eqref{eq:Adef}. The same random numbers are used for the two processes in the definition of both the jump times and the velocity directions.  The main contributions of the present paper are then twofold~:
\begin{itemize}
\item From a numerical point of view, we couple two systems of many particles consisting of different realizations of the fine-scale (internal state) and coarse (gradient sensing) processes, simulated with similar discretization schemes (see Section~\ref{sec:coupl_princ}) and the same random numbers. Additionally, for the coarse model, also the continuum description is simulated on a spatial grid. Evolution of the fine-scale model is then evaluated on the grid at hand by adding the difference between the two coupled particle descriptions to the evolution of the coarse continuum description. This can be seen as using the coarse process as a control variate. As coupling is lost over time, a regular reinitialization of the two coupled particle systems is included in the algorithm. The asymptotic variance reduction method can be depicted using Table~\ref{tab:table} (see also Figure~\ref{fig:pic} in Section~\ref{sec:method}).
\begin{table}\label{tab:table}
\begin{center}
  \begin{tabular}{|l|l|}
\hline
Model & Numerics \\
\hline
Internal state model--Individual description & Monte-Carlo method \\
Gradient sensing model--Individual description & Monte-Carlo method (coupled)\\
Gradient sensing model--Density description & Grid method \\
\hline
\end{tabular}
\end{center}
\caption{Models and numerical methods used in the variance reduction algorithm.}
\end{table}
	\item From an analysis point of view, we rigorously prove, using probabilistic arguments, $L^p$-bounds on the position difference at diffusive times of both processes in terms of appropriate powers of the scaling parameter (Theorem~\ref{thm:coupling}). This analysis shows that the numerical variance reduction obtained by the method is \emph{asymptotic} (the remaining statistical variance vanishes with the small scaling parameter), and provides upper bounds on the rate of the latter convergence. We argue in Section~\ref{sec:sharp} that this bound should be sharp in relevant regimes.
\end{itemize}

The idea of asymptotic variance reduction is a general and very recent idea in scientific computing that appears when using hybrid Monte-Carlo/PDE (partial differential equation) methods. While the use of control variates is fairly common in Monte Carlo simulation, the only explicit attempt to develop asymptotic variance reduction in hybrid methods, up to our knowledge, can be found in \cite{DimPar08} and related papers in the context of the Boltzmann equation, and is based on an importance sampling approach (see also~\cite{10.1063/1.1899210} and related papers).

This paper is organized as follows. In Section \ref{sec:model}, we briefly review the mathematical models that we will consider, and summarize the results on their advection-diffusion limit that were obtained in the companion paper \cite{limits}. In Section \ref{sec:method}, we proceed to describe the coupled numerical method in detail.
Section \ref{sec:variance} is dedicated to a detailed analysis of the variance reduction of the coupled method. The main result is the following: up to a time shift perturbation, the variance of the difference between the two coupled process scales according to the small parameter of the diffusive asymptotics.
We proceed with numerical illustrations in Sections \ref{sec:illustr} and \ref{sec:appl}, and conclude in Section \ref{sec:concl} with an outlook to future research.

\section{Particle-based models for bacterial chemotaxis \label{sec:model}}

In this paper, we directly present the models of interest in a nondimensional form, incorporating the appropriate space and time scales.  For a  discussion on the chosen scaling, we refer to the companion paper \cite{limits}.

\subsection{Model with internal state\label{sec:model-internal}}
We consider bacteria that are sensitive to the concentration of $m$ chemoattractants $\pare{\rho_i(x)}_{i=1}^m$, with $\rho_i(x) \geq 0$ for $x \in \mathbb{R}^d$. While we do not consider time dependence of chemoattractant via production or consumption  by the bacteria, a generalization to this situation is straightforward, at least for the definition of the models and the numerical method.  Bacteria move with a constant speed $v$ (run), and change direction at random instances in time (tumble), in an attempt to move towards regions with high chemoattractant concentrations.  As in \cite{ErbOthm04}, we describe this behavior by a velocity-jump process driven by some internal state $y \in \mathbb{Y}\subset\mathbb{R}^{n}$ of each individual bacterium. The internal state models the memory of the bacterium and is subject to an evolution mechanism attracted by a function  $\psi: \mathbb{R}^m \to \mathbb{R}^n$ of the chemoattractants concentrations,
\[
S(x):= \psi(\rho_1(x),\ldots,\rho_m(x)) \in \mathbb{R}^n,
\]
where $x$ is the present position of the bacterium. $S$ is assumed to be smooth with bounded derivatives up to order~$2$.
 We refer to Section~2.2 in~\cite{limits} for some typical choices; in the concrete example in this paper, we choose $m=n=1$, and $S(x)=\rho_1(x)$. By convention, the gradient of $S(x):\mathbb{R}^d \to \mathbb{R}^n$ is a matrix with dimension
\[
\nabla S(x) \in \mathbb{R}^{n\times d}.
\]
We then denote the evolution of each individual bacterium position by~
$
  t \mapsto X_t,
$
with normalized velocity
\[
  \frac{\d X_t}{\d t} = \epsilon V_t, \qquad V_t\in \mathbb{V}=\mathbb{S}^{d-1},
\]
with $\mathbb{S}^{d-1}$ the unit sphere in $\mathbb{R}^d$. Hence, $V_t$ represents the direction and the parameter $\epsilon$ represents the size of the velocity.  The evolution of the internal state is denoted by
$t \mapsto Y_t$. The internal state adapts to the local chemoattractant concentration through an ordinary differential equation (ODE),
\begin{equation}\label{eq:internal}
  \frac{\d Y_t}{\d t}= F_\epsilon(Y_t,S(X_t)),
\end{equation}
which is required to have a unique fixed point $y^*=S(x^*)$ for every fixed value $x^* \in \mathbb{R}^d$.
We also introduce the deviations from equilibrium $z=S(x)-y$.  The evolution of these deviations is denoted as
\[
t \mapsto Z_t = S(X_t) -Y_t.
\]
The velocity of each bacterium is switched at random jump times
 $(T_n)_{n \geq 1}$
that are generated via a Poisson process with a time dependent rate given by $\lambda(Z_t)$, where $z \mapsto \lambda(z)$ is a smooth function satisfying
\begin{equation}
  \label{eq:ratebound}
  0 < \lambda_{\rm min} \leq \lambda \pare{z} \leq \lambda_{\rm max},
\end{equation}
as well as,
\begin{equation}\label{eq:lin_rate}
\lambda(z) = \lambda_0 - b^T  z + c_\lambda\mathcal{O} \pare{ \abs{z}^k } ,
\end{equation}
with $b \in \mathbb{R} ^ n $, $k \geq 2$, and $c_\lambda >0$ is used to keep track of the non-linearity in the analysis. The new velocity at time $T_n$ is generated at random according to a centered probability distribution
$\mathcal M (dv)$
with $\int v \, \mathcal M (dv) =0 $, typically
\[
\mathcal M (dv) = \sigma_{\mathbb{S}^{d-1}}(dv),
\]
where $\sigma_{\mathbb{S}^{d-1}}$ is the uniform distribution on the unit sphere.

The resulting fine-scale stochastic evolution of a bacterium is then described by a left continuous with right limits (lcrl) process
 \[
 t \mapsto \pare{X_t, V_t, Y_t} ,
 \]
which satisfies the following differential velocity-jump equation~:
\begin{equation} \label{eq:process_noscale}
  \begin{cases}
   \displaystyle \dfrac{\d X_t}{\d t} = \epsilon  V_t,  \\[8pt]
    \displaystyle \dfrac{\d Y_t}{\d t} =  F_\epsilon(Y_t,S(X_t)), \\[8pt]
    \displaystyle \int_{T_n}^{T_{n+1}} \lambda(Z_t) \d t = \theta_{n+1}, \qquad \text {with }
    Z_t := S(X_t)-Y_t ,\\[8pt]
    \displaystyle V_{t}  =  \mathcal{V}_{n} \quad \text{ for $t \in [T_{n},T_{n+1}]$ } ,
  \end{cases}
\end{equation}
with initial condition $X_0, \mathcal{V}_0 \in \mathbb{R}^d$, $Y_0 \in \mathbb{R}^n$ and $T_0 = 0$.
In \eqref{eq:process_noscale}, $\pare{\theta _n }_{n \geq 1}$ denote i.i.d. random variables with normalized exponential distribution, and $\pare{ \mathcal{V}_n }_{n \geq 1}$ denote i.i.d. random variables with distribution $ \mathcal M (dv)$.

\begin{example}
For concreteness, we provide a specific example, adapted from \cite{ErbOthm04}, which will also be used later on to illustrate our results numerically. We consider $m=1$, i.e., there is only one chemoattractant $S(x)$ and the internal dynamics~\eqref{eq:internal} reduces to a scalar equation
\begin{equation}\label{e:scalar-y}
	\frac{\d y(t)}{\d t} = \frac{S(x)-y}{\tau}=\frac{z}{\tau}.
\end{equation}
For the turning rate $z\mapsto \lambda(z)$, we choose the following nonlinear strictly decreasing smooth function
\begin{equation}\label{eq:rate}
\lambda(z)=2\lambda_0 \pare{\frac{1}{2}-\frac{1}{\pi}\arctan \pare{\frac{\pi}{2\lambda_0}z}}.
\end{equation}
\end{example}

\subsection{Asymptotic regimes and linearizations\label{sec:asympt}}

In the adimensional models~\eqref{eq:process_noscale}-\eqref{eq:cprocess}, a small parameter $\epsilon$ is introduced, which can be interpreted as the ratio of the typical time of random change of velocity direction to  the typical time associated with chemoattractant variations as seen by the bacteria.  In the present section, we present the main scaling assumptions on the internal dynamics that are necessary to consider the asymptotic regime~$\epsilon \ll 1$, as well as the corresponding notation.  Most of these assumptions were introduced in the companion paper~\cite{limits}, to which we refer for more details and motivation.

\begin{rem}[Notation]
Throughout the text, the Landau symbol $\mathcal{O}$ denotes a \emph{deterministic} and globally Lipschitz function satisfying $\mathcal{O}(0)=0$. Its precise value \emph{may vary} from line to line, may depend on all parameters of the model, but its Lipschitz constant is uniform in $\epsilon$. In the same way, we will denote generically by $C > 0$ a deterministic constant that may depend on all the parameters of the model except $\epsilon$.
\end{rem}

We assume that the ODE~\eqref{eq:internal} driving the internal state is well approximated by a near equilibrium evolution equation in following sense~:
\begin{ass}\label{a:1}
 We have
\begin{equation}\label{eq:tau_eps_def}
 {F_{\eps}}(y,s) =  - {\tau_{\eps}}^{-1} (y-s) + \epsilon^{1-\delta}c_F\mathcal{O}(\abs{s-y}^2),
\end{equation}
where $\delta >0$ and ${\tau_{\eps}}  \in\mathbb{R}^{n\times n}$ is an invertible constant matrix.
\end{ass}
The constant $c_F$ is used to keep track of the dependance on the non-linearity of $F_\epsilon$ in the analysis. For technical reasons that are specific to the analysis of the time shift in the coupling (Lemma~\ref{lem:DTestimbis}), we need the following assumption on the scale of the non-linearity (this assumption was not necessary in the companion paper~\cite{limits}):
\begin{ass}\label{a:1bis}
 Assume that $\abs{y-s} = \mathcal{O}(\epsilon^\delta)$, then there is a $\gamma > 0$ such that:
\begin{equation}\label{eq:taueps_hypbis}
\abs { (y-s) - {\tau_{\eps}}{F_{\eps}}(y,s) } =  c_F\mathcal{O}( \epsilon^{1 + \gamma} ).
\end{equation}
\end{ass}

We now formulate the necessary assumptions on the scale of the matrix ${\tau_{\eps}}$~:
\begin{ass}\label{a:2}
There is a constant $C >0$ such that for any $t \geq 0$, one has
\begin{equation}\label{eq:tau_scales}
  \norm{{\rm e}^{- t {\tau_{\eps}} ^{-1} }} \leq C {\rm e}^{- t \epsilon^{1-\delta}/C}.
\end{equation}
\end{ass}
\begin{ass}\label{a:3}
There is a constant $C >0$ such that
\[
 \sup_{t \geq 0} \norm{t{\tau_{\eps}} ^{-1}{\rm e}^{- t {\tau_{\eps}} ^{-1} } } \leq C  \epsilon^{-1}.
\]
\end{ass}
Assumption~\ref{a:2} is used to ensure exponential convergence of the linear ODE with time scale at least of order $\epsilon^{1-\delta}$. Assumption~\ref{a:3} is necessary for technical reasons in the proof of Lemma~\ref{lem:DTestim} (which is given in~\cite{limits}). Note that if ${\tau_{\eps}}^{-1}$ is symmetric and strictly positive with a lower bound of order $\mathcal{O}(\epsilon^{1-\delta})$ on the spectrum, then Assumption~\ref{a:2} and Assumption~\ref{a:3} are satisfied. This implies that the slowest time scale of the ODE~\eqref{eq:internal} is at worse of order $\epsilon^{1-\delta}$.  Indeed, integration over time in Assumption~\ref{a:2} implies $\norm{{\tau_{\eps}}} \leq C \epsilon^{\delta-1}$.
Moreover, we  assume that the solution of the ODE driving the internal state in \eqref{eq:process_noscale} satisfies the following long time behavior~:
\begin{ass}
   \label{ass:ODE}
Consider a path $t \mapsto S_t$ such that $ \sup_{t \in [0,+\infty] } \abs{ \dfrac{ d S_t}{ dt} }= \mathcal{O}( \epsilon )$, and assume $\abs{Y_0 - S_0} = \mathcal{O}(\epsilon^\delta)$. Then the solution of
\[
\frac{ d Y_t}{ dt} = {F_{\eps}}(Y_t,s_t),
\]
satisfies
\[
\sup_{t \in [0,+\infty] } \abs{ Y_t - S_t } = \mathcal{O} \pare{ \epsilon^\delta  }.
\]
\end{ass}
Assumption~\ref{ass:ODE} is motivated by the case of linear ODEs; indeed, in that case, Assumption~\ref{ass:ODE} is a consequence of Assumption~\ref{a:2} (see Section~3.3 in~\cite{limits}).
Recalling that $Z_t=S_t-Y_t$, we will assume in the remainder of the paper that the solution of~\eqref{eq:process_noscale} satisfies,
\[
\abs{Z_0} = \mathcal{O}(\epsilon^\delta),
\]
as well as
\[
\abs{{\tau_{\eps}}^{-1} Z_0} = \mathcal{O}( 1 ).
\]
Then, under Assumption~\ref{a:1},  Assumption~\ref{a:2},  Assumption~\ref{a:3}, and Assumption~\ref{ass:ODE}, we obtain the bounds (see Section~3.3 in~\cite{limits})
\begin{equation}
    \label{eq:techbound}
    \sup_{t \in [0,+\infty] } \abs{ {\tau_{\eps}}^{-1}{\rm e}^{- {\tau_{\eps}}^{-1} t } (S_t - Y_t) } = \mathcal{O} \pare{ 1 } ,
  \end{equation}
\[
\sup_{t \in [0,+\infty]} \pare{ \abs{Z_t} } =  \mathcal{O}(\epsilon^{\delta}),
\]
as well as
\[
\sup_{t \in [0,+\infty]} \pare{ \abs{ {\tau_{\eps}}^{-1} {\rm e}^{- {\tau_{\eps}}^{-1} t} Z_t} } =  \mathcal{O}(1) .
\]

Finally, we will make use of the following notation. Denoting
$c_S = \norm{\nabla^{2}  S}_{\infty},$
with $\nabla^2$ the Hessian,
the error terms due to non-linearities of the problem will be handled through
\begin{equation}
  \label{eq:nl_err}
  \begin{cases}
      {\rm nl} (\epsilon) := c_{F} \epsilon^{\delta} + c_S \epsilon + c_{\lambda} \epsilon^{ k \delta-1}  = o(1),\\
      {\rm nl_{2}} (\epsilon) := c_{F}( \epsilon^{\delta} + \epsilon^\gamma)+ c_S( \epsilon + \epsilon^\delta) + c_{\lambda} \epsilon^{ k \delta-1}  = o(1) .
  \end{cases}
\end{equation}
 Note that $  {\rm nl} (\epsilon) \leq  {\rm nl_{2}} (\epsilon)$ and $  {\rm nl} (\epsilon) \sim  {\rm nl_{2}} (\epsilon)$ if $\delta \leq 1$ and $\delta \leq \gamma$. Recall that $k \geq 2$ is defined in~\eqref{eq:lin_rate}.

Note that all the assumptions of the present section hold in the following case: (i) the ODE is linear ($F_\epsilon$ is linear), (ii) the associated matrix $\tau_{\epsilon}$ is symmetric positive and satisfies
\[
\tau_{\epsilon}^{-1} \geq C \epsilon^{1-\delta},
\]
for some $\delta > 0$, (iii) the initial internal state is close to equilibrium in the sense that ${\tau_{\eps}}^{-1} Z_0=\mathcal{O}(1)$. The assumptions of the present section can be seen as technical generalizations to non-linear ODEs with non-symmetric linear part.

\subsection{Model with direct gradient sensing (control process)}

We now turn to a simplified, ``coarse'' model, in which the internal process (\ref{eq:internal}), and the corresponding state variables, are eliminated.  Instead, the turning rate  depends directly on the chemoattractant gradient. This process will be called the \emph{control process} since it will be used in Section~\ref{sec:method} as a control variate to perform variance reduced simulations of (\ref{eq:process_noscale}).

The control process is a Markov process in position-velocity variables
\[
t \mapsto (X^c_t,V^c_t) ,
\]
which evolves according to the following differential velocity-jump equations~:
\begin{equation} \label{eq:cprocess}
\system{
& \frac{\d X_t^c}{\d t} = \epsilon V_t^c , &\\
& \int_{T_n^c}^{T_{n+1}^c} \lambda^c_\epsilon(X_t^c , V_t^c) \d t = \theta_{n+1} , & \\
& V_{t}^c  =  \mathcal{V}_n  \quad \text{ for $t \in [T_{n}^c,T^c_{n+1}]$ } , &
}
\end{equation}
with initial condition $X_0, \mathcal{V}_0 \in \mathbb{R}^d$.  In (\ref{eq:cprocess}), $\pare{\theta _n }_{n \geq 1}$ denote i.i.d. random variables with normalized exponential distribution, and $\pare{ \mathcal{V}_n }_{n \geq 1}$ denote i.i.d. random variables with distribution $ \mathcal M (dv)$.
The turning rate of the control process is assumed to satisfy
\begin{equation}
  \label{eq:rateboundc}
  0 < \lambda_{\rm min} \leq \lambda^c_{\epsilon}(x,v) \leq \lambda_{\rm max},
\end{equation}
as well as
\begin{equation}\label{eq:control_rate}
\lambda^c_\epsilon(x,v) := \lambda_0 - \epsilon\; A^T_\epsilon(x)  v  + \mathcal{O}(\epsilon^2).
\end{equation}

Typically, $\lambda^c_{\epsilon}(x,v)$ is a function of $\nabla S(x)$, so that the model \eqref{eq:control_rate} may describe a large bacterium that is able to directly sense chemoattractant gradients. When $m=n=1$, and the turning rate (\ref{eq:control_rate}) is proportional to $ \nabla S(x)  v \in \mathbb{R}$, it can be interpreted as follows: the rate at which a bacterium will change its velocity direction depends on the alignment of the velocity with the gradient of the chemoattractant concentration $\nabla S(x)$, resulting in a transport towards areas with higher chemoattractant concentrations.
In the companion paper \cite{limits}, it is shown that, when the vector field $A_\epsilon(x) \in \mathbb{R}^{d}$ is given by the formula
\begin{equation}
  \label{eq:Adef}
  A_\epsilon(x) = b^T \frac{{\tau_{\eps}}}{\Id + \lambda_0 {\tau_{\eps}}} \nabla S(x),
\end{equation}
asymptotic consistency with the internal state model is obtained in the limit when $\epsilon\to 0$ (see Section~\ref{sec:model-hydro}).

\subsection{Kinetic formulation and advection-diffusion limits\label{sec:model-hydro}}

We now turn to the kinetic description of the probability distribution of the processes introduced above.
The probability distribution density of the fine-scale process with internal state at time $t$ with respect to the measure $\d x \, \mathcal M ( \d v ) \, \d y$ is denoted as
$p(x,v,y,t)$, suppressing the dependence on $\epsilon$ for notational convenience, and evolves according to the Kolomogorov forward evolution equation (or master equation). In the present context, the latter is the following kinetic equation
\begin{equation} \label{e:kinetic}
\partial_t p + \epsilon v \cdot \nabla_x p + \div_y  \pare{ F_\epsilon(x,y) p } = \lambda \pare{S(x)-y} \pare{   R(p)    - p   },
\end{equation}
where
\[
R(p) := \int_{v \in \mathbb{S}^{d-1}} p(\cdot,v,\cdot) \, \mathcal M(dv)
\]
is the operator integrating velocities with respect to $\mathcal M$, and $F(x,y)$ and $\lambda(z)$ are defined as in Section \ref{sec:model-internal}.
Similarly, the distribution density of the control process is denoted as
$p^c(x,v,t)$, and evolves according to
the kinetic/master equation:
\begin{equation} \label{e:kinetic_control}
\partial_t p^c + \epsilon v \cdot \nabla_x p^c = \pare{   R(\lambda^c_\epsilon p)    - \lambda^c_\epsilon p   },
\end{equation}
with $\lambda^c_\epsilon$ defined as in \eqref{eq:control_rate}.
The reader is referred to \cite{EthKur86} for the derivation of master equations associated to Markov jump processes.

In the companion paper \cite{limits}, we show, using probabilistic arguments, that, in the limit $\epsilon \to 0$, both the equation for the control process~\eqref{e:kinetic_control} and the equation for the process with internal state~\eqref{e:kinetic} converge to an advection-diffusion limit on diffusive time scales.  Convergence is to be understood \emph{pathwise}, in the sense of convergence of probability distribution on paths endowed with the uniform topology.

Denote
diffusive times by
\[
{\bar{t}}:=t  \epsilon^2,
\]
and the processes considered on diffusive time scales as
\[
X^{\epsilon}_{\bar{t}}:=X_{{\bar{t}}/\epsilon^2}, \qquad X^{c,\epsilon}_{\bar{t}}:=X^c_{{\bar{t}}/\epsilon^2}.
\]
For the control process, we then have the following theorem ~:
\begin{theorem}
For $\epsilon \to 0$, the process ${\bar{t}} \mapsto X^{c,\epsilon}_{\bar{t}}$, solution of~\eqref{eq:cprocess}, converges towards an advection-diffusion process, satisfying the stochastic differential equation (SDE)
\begin{equation}\label{eq:cprocess_hydro}
\d X_{\bar{t}}^{c,0} =\pare{\frac{D A_0(X_{\bar{t}}^{c,0})}{\lambda_0} \d {\bar{t}} + \pare{\frac{2D}{\lambda_0}}^{1/2} \d W_{\bar{t}}},
\end{equation}
where ${\bar{t}} \mapsto  W_{\bar{t}}$ is a standard Brownian motion, the parameters $\lambda_0$ and $A_0:=\lim_{\epsilon\to 0} A_\epsilon$ originate from the turning rate~\eqref{eq:control_rate}, and the diffusion matrix is given by the covariance of the Maxwellian distribution~:
\begin{equation}
  \label{eq:D}
  D =  \int_{\mathbb{S}^{d-1}} v \otimes v \, \mathcal M(dv) \in \mathbb{R}^{d \times d}.
\end{equation}
\end{theorem}
In particular, this result implies that, at the level of the Kolomogorov/master evolution equation, the evolution of the position bacterial density,
\begin{equation}\label{eq:cdensity}
	n^{c,\epsilon}(x,{\bar{t}}):= n^c(x,{\bar{t}}/\epsilon^2):=\int_{\mathbb{V}}p^c(x,v,{\bar{t}}/\epsilon^2) \,  \mathcal M(\d v)
\end{equation}
converges to the advection-diffusion equation \begin{equation}\label{eq:cdens_hydro}
\partial_{\bar{t}} n^{c,0} = \frac{1}{\lambda_0}\div_x \pare{ D \nabla_x n^{c,0} - D A_0(x) n^{c,0}}
\end{equation}
on diffusive time scales as $\epsilon \to 0$.

In the same way, a standard probabilistic diffusion approximation argument can be used to derive the pathwise diffusive limit of the process with internal state~\eqref{eq:process_noscale}~:
\begin{theorem}
For $\epsilon \to 0$, the process ${\bar{t}} \mapsto X^{\epsilon}_{\bar{t}}$, solution of~\eqref{eq:process_noscale}, converges towards an advection-diffusion process, satisfying the stochastic differential equation (SDE)~\eqref{eq:cprocess_hydro}, where $A_0$ originates from
\begin{equation}\label{eq:control_field}
A_0(x) = b^T \lim_{\epsilon \to 0}  \frac{{\tau_{\eps}}}{\lambda_0 {\tau_{\eps}} + \Id } \nabla S(x),
\end{equation}
in which, $b$, ${\tau_{\eps}}$, and $\lambda_0$ were introduced in \eqref{eq:tau_eps_def}-\eqref{eq:lin_rate} as parameters of the process with internal state, and $\Id \in \mathbb{R}^{n\times n}$ is the identity matrix. Again, the diffusion matrix $D$ is given by the covariance of the Maxwellian distribution~\eqref{eq:D}.
\end{theorem}
Also here, convergence needs to be understood in terms of probability distribution on paths endowed with the uniform topology.

Introducing the bacterial density of the process with internal state as
\begin{equation}
	n(x,t)= \int_{\mathbb{Y}}\int_{\mathbb{V}}p(x,v,y,t) \mathcal M(\d v) \d y,
\end{equation}
this implies that the evolution of $n$ converges to \eqref{eq:cdens_hydro} on diffusive time scales in the limit of $\epsilon \to 0$.

\section{Numerical method\label{sec:method}}

To simulate the process with internal state, solving the kinetic equation \eqref{e:kinetic} over diffusive time scales can be cumbersome, due to the additional dimensions associated with the internal state. The alternative is to use to stochastic particles. However, a particle-based simulation of equation \eqref{e:kinetic}  is subject to a large statistical variance of the order $\mathcal{O}(N^{-1/2})$, where $N$ is the number of simulated particles. The asymptotic analysis shows that the position bacterial density approaches an advection-diffusion limit \eqref{eq:cdens_hydro} when $\epsilon\to 0$. Consequently, to accurately assess the deviations of the process with internal state \eqref{eq:process_noscale} as compared to its advection-diffusion limit (for small but non-vanishing (intermediate) values of $\epsilon$), the required number of particles needs to increase substantially with decreasing $\epsilon$, which may become prohibitive from a computational point of view.

In this section, we therefore propose a hybrid method, based on the principle of control variates, that couples the process with internal dynamics with the control process, which is simulated simultaneously using a grid-based method. We first describe the variance reduction technique (Section \ref{sec:coupl_princ}).  The analysis in Section~\ref{sec:variance} will reveal that the variance reduction is \emph{asymptotic}, in the sense that the variance vanishes in the diffusion limit.  To ensure this asymptotic variance reduction during actual simulations, one needs to ensure that the time discretization preserves the diffusion limits of the time-continuous process.  An appropriate time discretization is discussed in Section \ref{sec:discr}.

\subsection{Coupling and asymptotic variance reduction\label{sec:coupl_princ}}

The proposed variance reduction technique is based on the introduction of a control variate that exploits a coupling between the process with internal state and the control process.
While the idea of control variates for Monte Carlo simulation is already well known, see e.g.~\cite{KloPla92} and references therein, the coupling that is proposed here is particular (we call it \emph{asymptotic}), since the difference of the coupled processes, and hence the variance, vanishes with an estimable rate in the diffusion limit $\epsilon\to 0$, as will be shown in Section \ref{sec:variance}.

\subsubsection{The control variates}
Let us first assume that we are able to compute the exact solution of the kinetic equation for the control process, \eqref{e:kinetic_control}, with infinite precision in space and time.

The algorithm of asymptotic variance reduction is based on a coupling between an ensemble of realizations evolving according to the process with internal state \eqref{eq:process_noscale}, denoted as
\[ \left\{X^i_t,V^i_t,Y^i_t\right\}_{i=1}^N,
\]
and an ensemble of realization of the control process \eqref{eq:cprocess}, denoted as
\[ \left\{X^{i,c}_t,V^{i,c}_{t}\right\}_{i=1}^N.
\]
We denote the empirical measure of the particles with internal state in position-velocity space as
\[
\mu_{\bar{t}}^{N}(x,v) = \frac{1}{N} \sum_{i=1}^{N} \delta_{X^i_{{\bar{t}}/\epsilon^2},V^i_{{\bar{t}}/\epsilon^2}},
\]
and, correspondingly, the empirical measure of the control particles as
\[
\mu_{\bar{t}}^{c,N} = \frac{1}{N} \sum_{i=1}^{N} \delta_{X^{c,i}_{{\bar{t}}/\epsilon^2},V^{c,i}_{{\bar{t}}/\epsilon^2}}.
\]

A coupling between the two ensembles is obtained by ensuring that both simulations use \emph{the same random numbers} $(\theta_n)_{n\ge 1}$ and $(\mathcal{V}_n)_{n\ge 0}$, which results in a strong correlation between $(X^i_{t},V^{i}_t)$ and $(X_t^{i,c},V^{c,i}_t)$ for each realization.
Simultaneously, the kinetic equation for the control process \eqref{e:kinetic_control} is also solved using a deterministic method (which, for now, is assumed to be exact). We formally denote the corresponding semi-group evolution as
\[
{\rm e}^{{\bar{t}} L^c}, \qquad \text{ with } L^c(p^c)= - \epsilon v \cdot \nabla_x p^c + \pare{   R(\lambda^c_\epsilon p^c)    - \lambda^c_\epsilon p^c   }.
\]
Besides the two particle measures $\mu^N_{\bar{t}}$ and $\mu^{c,N}_{\bar{t}}$, we denote by $\overline{\mu}^N_{\bar{t}}$ the variance reduced measure, which will be defined by the algorithm below. Since, with increasing diffusive time, the variance of the algorithm increases due to a loss of coupling between the particles with internal state and the control particles, the variance reduced algorithm will also make use of a reinitialization time step $\overline{\delta t}_{ri}$, which is defined on the diffusive time scale. The corresponding time instances are denoted as $\bar{t}_n=n\overline{\delta t}_{ri}$ on the diffusive time scale, or equivalently, on the original time scale as $t_n=n\overline{\delta t}_{ri}/\epsilon^2$.

Starting from an initial probability measure $\mu_0$ at time $t=0$, we sample $\mu_0$ to obtain the ensemble $\left\{X^i_t,V^i_t,Y^i_t\right\}_{i=1}^N$, corresponding to $\mu^N_0$, and then set $\mu^{c,N}_0:=\mu^N_0$, i.e., $X^{i,c}_{0}=X^{i}_{0}$ and $V^{i,c}_{0}=V^{i}_{0}$ for all $i=1,\ldots,N$.  Furthermore, we set the variance reduced estimator as $\overline{\mu}_0^N:=\mu_0=\E(\mu_0^N)$. We then use the following algorithm to advance from $\bar{t}_n$ to $\bar{t}_{n+1}$, (see also Figure~\eqref{fig:pic})~:
\begin{figure}

  \includegraphics[scale=0.5]{figuren/fig-pic.pdf}
	\caption{\label{fig:pic} A schematic description of Algorithm \ref{algo}. The dashed line represent the evolution of $N$ bacteria with internal state. The dotted line represent the coupled evolution of $N$ bacteria with gradient sensing, subject to regular reinitializations. The dashed-dotted line is computed according to a deterministic method simulating the density of the model with gradient sensing, and subject to regular reinitializations. The solid line is the variance reduced simulation of the internal state dynamics, and is computed by adding the difference between the particle computation with internal state, and the particle simulation with gradient sensing to the deterministic gradient sensing simulation.  At each reinitialization step, the two simulations (deterministic and particles) of the gradient sensing dynamics are reinitialized to the values of their internal state simulation counterpart (as represented by the arrows).}
\end{figure}
\begin{alg}\label{algo}
At time $t_n$, we have that the particle measure $\mu_{\bar{t}_n}^{c,N}= \mu_{\bar{t}_n}^{N}$, and the variance reduced measure is given by  $\overline{\mu}^N_{\bar{t}_n}$.
To advance from time $ \bar{t}_n$ to $\bar{t}_{n+1}$, we perform the following steps~:
  \begin{itemize}
  \item Evolve the particles $\left\{X^i_t,V^i_t,Y^i_t\right\}_{i=1}^N$ from $t_n$ to $t_{n+1}$, according to \eqref{eq:process_noscale},
   \item  Evolve the particles $\left\{X^{i,c}_t,V^{i,c}_{t}\right\}_{i=1}^N$ according to \eqref{eq:cprocess}, using the same random numbers as for the process with internal state,
 \item Compute the variance reduced evolution
   \begin{equation}
     \label{eq:var_reduced_estimation}
     \overline{\mu}^N_{\bar{t}_{n+1}} = \overline{\mu}^N_{\bar{t}_n} {\rm e}^{\overline{\delta t}_{ri} / \epsilon^2 L^c}  + \pare{\mu_{\bar{t}_{n+1}}^{N} - \mu_{\bar{t}^{-}_{n+1}}^{c,N} }.
   \end{equation}
\item Reinitialize the control particles by setting
\[
X^{i,c}_{t_{n+1}}=X^{i}_{t_{n+1}}, \qquad V^{i,c}_{t_{n+1}}=V^{i}_{t_{n+1}}, \qquad i=1,\ldots,N,
\]
i.e., we set the state of the control particles to be identical to the state of the particles with internal state.
\end{itemize}
\end{alg}
In formula \eqref{eq:var_reduced_estimation}, we use the symbol $\bar{t}^{-}_{n+1}$ to emphasize that the involved particle positions and velocities are those obtained \emph{before} the reinitialization.
An easy computation shows that the algorithm is unbiased in the sense that for any $n \geq 0$,
\[
\E \pare{ \overline{\mu}^N_{\bar{t}_n}} = \E \pare{\mu^N_{\bar{t}_n}},
\]
since the particles with internal dynamics are unaffected by the reinitialization, and, additionally,
\[
\E \mu^{c,N}_{\bar{t}_{n+1}}= \E \overline{\mu}^N_{{\bar{t}}_n}{\rm e}^{\overline{\delta t}_{ri}/ \epsilon^2 L^c}.
\]
Moreover, the variance is controlled by the coupling between the two processes. Indeed, using the independence of the random numbers between two steps of Algorithm~\ref{algo}, and introducing $\varphi$ as a position and velocity dependent test function), we get (according to the main result of Theorem~\ref{thm:coupling})
\begin{align}
  \var( \overline{\mu}^N_{\bar{t}_n}(\varphi)  ) &= \sum_{k=1}^n \E \pare{\mu_{}^{N}(\varphi) - \mu_{\bar{t}_k^{-}}^{c,N}(\varphi) }^2 \nonumber \\
 & \leq  \sum_{k=1}^n \frac{\norm{\nabla \varphi}_\infty}{N} \E \pare{  \abs{X_{\bar{t}_k/\epsilon^2} - X_{\bar{t}_k/\epsilon^2}^c  }^2 } , \nonumber
\end{align}
and thus
 \begin{align}\label{eq:MC_bound}
  \var( \overline{\mu}^N_{\bar{t}_n}(\varphi)  )
& \leq C n\frac{ \epsilon + \epsilon^\delta +  {\rm nl_{2}} (\epsilon) }{N},
\end{align}
where in the last line, $C$ is independant of $n$, $\epsilon$, and $N$.

In some generic situations (see Section~\ref{sec:sharp}), we can argue that the statistical error in Algorithm~\ref{algo} coming from the coupling is ``sharp'' with respect to the order in $\epsilon$. This means that the difference between the probability distribution of the model with internal state and the probability distribution of the model with gradient sensing is of the same order.  This would imply that, with the asymptotic variance reduction technique, one is able to reliably assess the true deviation of the process with internal variables from the control process using a number of particles $N$ that is independent of $\epsilon$.

\subsubsection{Time-space discretization of the kinetic equation}

For both processes, the probability density distribution of the position and velocity can then be computed via binning in a histogram, or via standard kernel density estimation \cite{WScott:1992p8096,WSilverman:1986p8151} using a kernel~$K_h$, in which the chosen bandwidth $h$ can be based on the grid used for the deterministic simulation and on the data, for instance using the Silverman heuristic \cite{WSilverman:1986p8151}. This yields
\begin{equation}\label{eq:kde}
	\hat{p}_N(x,v,t) = \frac{1}{N h} \sum_{i=1}^{N} K_h(x - X^i_{t},v-V^i_{t}),
\end{equation}
and similarly
\begin{equation}\label{eq:kdec}
	\hat{p}^c_N(x,v,t) = \frac{1}{N h} \sum_{i=1}^{N} K_h(x - X^{c,i}_{t},v-V^{c,i}_{t}).
\end{equation}

The important point is that the solution $p^c(x,v,t)$ for the control process may be accurately approximated with a deterministic grid-based method, which ensures that
\[
\abs{p^c(x,v,{\bar{t}} / \epsilon^2)-p^c_{\rm grid}(x,v,{\bar{t}} / \epsilon^2)} = \mathcal{O}(\delta x^l) +\mathcal{O}(\delta t^k),
\]
for some integers $k,l\ge 1$.
\begin{remark}
In this text, we will perform the simulations in one space dimension using a third-order upwind-biased scheme, and perform time integration using the standard fourth order Runge--Kutta method, i.e., $l=3$ and $k=4$.
\end{remark}
Then, the unbiased nature of the variance-reduced estimator is conserved up to $\mathcal{O}(\delta x^l)+\mathcal{O}(\delta t^k)+\mathcal{O}(h)$ discretization errors.

\subsection{Time discretization of velocity-jump processes \label{sec:discr}}

When simulating equation~\eqref{eq:process_noscale}, a time discretization error originates from the fact that the equation for the evolution of the internal state, and hence the evolution of the fluctuations $Z_t$, is discretized in time.  This results in an approximation of the jump times $(T_n)_{n\ge 1}$, and hence of $X_t$. To retain the diffusion limits of the time-continuous process, special care is needed. We now briefly recall the time discretization procedure that was proposed and analyzed in the companion paper~\cite{limits}. For ease of exposition, we consider the scalar equation \eqref{e:scalar-y} for the internal state; generalization to nonlinear systems of equations is briefly discussed in~\cite{limits}.

Consider first the linear turning rate \eqref{eq:lin_rate}. In that case, we define a numerical solution  $(\disc{X}_t,\disc{V}_t,\disc{Y}_t)$ as follows.
Between jumps, we discretize the simulation in steps of size $\delta t$ and denote by $(\disc{X}_{n,k},\disc{Z}_{n,k})$ the solution at $t_{n,k}=\disc{T}_n+k\delta t$.  The numerical solution for $t\in[t_{n,k},t_{n,k+1}]$ is given by
\begin{equation} \label{e:discr-nonlin}
  \begin{cases}
    \disc{X}_t = \disc{X}_{n,k}+\epsilon \mathcal{V}_{n}\;(t-t_{n,k})  \\
    \disc{Z}_t = \exp(-(t-t_{n,k}) {\tau_{\eps}}^{-1}) \disc{Z}_{n,k}  + \epsilon {\tau_{\eps}} \pare{\Id-\exp(-(t-t_{n,k}){\tau_{\eps}}^{-1})} \nabla S(\disc{X}_{n,k}) \, \mathcal{V}_n.
  \end{cases}
\end{equation}
We denote by $K\geq 0$ the integer such that the simulated jump time $\disc{T}_{n+1}\in [t_{n,K},t_{n,K+1}]$. To find $\disc{T}_{n+1}$, we first approximate the integral $ \int_{\disc{T}_n}^{\disc{T}_{n+1}}\lambda(\disc{Z}_t)\d t$ using
\begin{align}\label{eq:time-disc-lin}
	\int_{\disc{T}_n}^{\disc{T}_{n+1}}\lambda(\disc{Z}_t)\d t &=\sum_{k=0}^{K-1}\int_{t_{n,k}}^{t_{n,k+1}}\lambda(\disc{Z}_t)\d t
	+ \int_{t_{n,K}}^{\disc{T}_{n+1}}\lambda(\disc{Z}_t)\d t ,
\end{align}
and then compute~:
\begin{multline}
	\int_{t_{n,k}}^{t_{n,k+1}}\lambda(\disc{Z}_t)\d t=
	\lambda_0 \delta t - b^T \pare{\Id - {\rm e}^{- {\delta t}{{\tau_{\eps}}^{-1}}} }{{\tau_{\eps}}} \disc{Z}_{T_n}  \\
 - \epsilon b^T \pare{\delta t {\tau_{\eps}} -(\Id-{\rm e}^{-\delta t {\tau_{\eps}}^{-1} } ){\tau_{\eps}}^2   }  \nabla S(\disc{X}_{n,k}) \mathcal{V}_n .
\end{multline}
The jump time $\disc{T}_{n+1}$ can then be computed as the solution of
\begin{equation}\label{e:lin-Newton}
\int_{t_{n,K}}^{\disc{T}_{n+1}}\lambda(\disc{Z}_t)\d t=\theta_{n+1}- \sum_{k=0}^{K-1}\int_{t_{n,k}}^{t_{n,k+1}}\lambda(\disc{Z}_t)\d t,
\end{equation}
using a Newton procedure. It is shown in \cite{limits} that the results on the diffusion limit (as outlined in Section \ref{sec:model-hydro}) are not affected by the discretization.

We now consider a general nonlinear turning rate. We again discretize in time to obtain the time-discrete solution \eqref{e:discr-nonlin}.  The jump time $\disc{T}_{n+1}$ is now computed by linearizing \eqref{eq:rate} in each time step,
\begin{equation}\label{eq:nonlinrate_disc}
\disc{\lambda}\pare{\disc{Z}_t,\disc{Z}_{n,k}}=\lambda(\disc{Z}_{n,k})+ \frac{\d \lambda(\disc{Z}_{n,k})}{\d z}\pare{\disc{Z}_t-\disc{Z}_{n,k}}.
\end{equation}
One then replaces equation \eqref{eq:time-disc-lin} by
\begin{align}
	\int_{\disc{T}_n}^{\disc{T}_{n+1}}\disc{\lambda}(\disc{Z}_t)\d t &=\sum_{k=0}^{K-1}\int_{t_{n,k}}^{t_{n,k+1}}\disc{\lambda}(\disc{Z}_t,\disc{Z}_{n,k})\d t
	+ \int_{t_{n,K}}^{\disc{T}_{n+1}}\disc{\lambda}(\disc{Z}_t,\disc{Z}_{n,K})\d t,
\end{align}
and proceeds in the same way as for the linear case.  Also in this case, the diffusive limit is recovered in an exact fashion for the time-discretized process.

\begin{remark}[Discretization of the control process]
For the control process, there is no internal state. Hence, the only time dependence of turning rate $\lambda^c$ is due to the spatial variation of $\nabla S(x)$, which can be treated by discretizing the integral of $\lambda$ in the same way as above.
\end{remark}

\section{Asymptotic variance reduction of the coupling\label{sec:variance}}

In this section, we show how the difference between the two coupled processes on diffusive time scales $\bar{t}=t/\epsilon^2$ behaves in the limit of $\epsilon\to 0$.  We first recall some notation from the companion paper~\cite{limits} (Section~\ref{sec:not}), after which we state and prove the main theorem (Section~\ref{sec:proof}).

\subsection{Notations and asymptotic estimates\label{sec:not}}

The main theorem in Section~\ref{sec:proof} relies on following auxiliary definitions and lemmas that were given in the companion paper~\cite{limits}:
\begin{defn}
  We denote by $m:\mathbb{R} \to \mathbb{R}^{n \times n}$ the function
  \begin{equation}
    \label{eq:m(t)}
    m(t) := t {\tau_{\eps}} - \pare{\Id - {\rm e}^{-t {\tau_{\eps}}^{-1}}} {\tau_{\eps}}^2,
  \end{equation}
whose derivative is given by
\[
m'(t) = {\tau_{\eps}}\pare{\Id - {\rm e}^{-t {\tau_{\eps}}^{-1}} }.
\]
\end{defn}
This function satisfies the following lemmas~:
\begin{lem}
  Let $\theta$ be an exponential random variable of mean $1$. Then:
  \begin{equation}
    \label{eq:av_m}
    \E \pare{ m\left( \frac{\theta}{\lambda_0} \right)  } = \frac{1}{\lambda_0} \frac{{\tau_{\eps}}}{\Id + \lambda_0 {\tau_{\eps}}}=A_\epsilon(x),
  \end{equation}
see equation \eqref{eq:Adef}.
\end{lem}
\begin{lem}
For all $t \in \mathbb{R}$, we have
$\norm{ m'(t) } \leq t$, as well as
$\norm{ m(t) } \leq t^2/2$.
\end{lem}

The difference between jump times is denoted as
\begin{equation*}
\Delta T_{n+1}^{c} :=  T_{n+1}^{c}- T_{n}^{c}, \qquad \Delta T_{n+1} :=  T_{n+1}- T_{n}.
\end{equation*}

The proofs of the asymptotic variance reduction make use of the following asymptotic expansions of the jump time differences of both processes~:
\begin{lem}\label{lem:dt-control}
The difference between two jump times of the control process satisfies
\begin{equation}\label{eq:deltaTc}
\Delta T_{n+1}^{c} = \frac{\theta_{n+1}}{\lambda_0} + \epsilon \frac{\theta_{n+1}}{\lambda_0^2} A_\epsilon^T(X_{T^c_n}^c) \mathcal{V}_n + \theta_{n+1} \mathcal{O} (\epsilon^2  ).
\end{equation}
\end{lem}
and
\begin{lem}\label{lem:DTestim}
The jump time variations of the process with internal state can be written in the following form ($ {\rm nl} (\epsilon)$ being defined by~\eqref{eq:nl_err})~:
\begin{equation}\label{eq:deltaT}
  \Delta T_{n+1} = \Delta T_{n+1}^0 + \epsilon \Delta T_{n+1}^1 + (\theta_{n+1}^6+\theta_{n+1})\mathcal{O} \pare{\epsilon^2+ \epsilon\, {\rm nl} (\epsilon)} ,
\end{equation}
where
\begin{align}\label{eq:estimDT0}
\Delta T_{n+1}^0 &= \frac{\theta_{n+1}}{\lambda_0} + \frac{b^T}{\lambda_0} m'(\Delta T_{n+1}^0 ) Z_{T_n} \nonumber \\
&= \frac{\theta_{n+1}}{\lambda_0} + \theta_{n+1} \mathcal{O} ( \epsilon^\delta)
\end{align}
and, correspondingly,
\begin{align}\label{eq:estimDT1}
  \Delta T_{n+1}^1 &= \frac{1}{\lambda_0 - b^T {\rm e}^{- {\Delta T_{n+1}^0}{\tau^{-1}}} Z_{T_n} } b^T m(\Delta T_{n+1}^0) \nabla  S(X_{T_n})  \mathcal{V}_{n}\nonumber  \\
& = \frac{1}{\lambda_0}   b^T m\pare{\frac{\theta_{n+1}}{\lambda_0}} \nabla S(X_{T_n}) \mathcal{V}_{n} +  \theta_{n+1} \mathcal{O}(\epsilon^\delta) .
\end{align}
\end{lem}
It is also useful to recall that according to \eqref{eq:ratebound}-\eqref{eq:rateboundc}, the following hold~:
\[
\abs{\Delta T_{n+1}} \leq C \theta_{n+1}, \qquad \abs{\Delta T_{n+1} ^c} \leq C \theta_{n+1}.
\]
Finally, we will also need a different expansion of the jump times of the process with internal state~:
\begin{lem}\label{lem:DTestimbis}
When using $ {\rm nl} _2(\epsilon)$ as defined by~\eqref{eq:nl_err}, the jump time variations of the process with internal state can be written in the following form~:
\begin{equation}\label{eq:deltaTbis}
  \begin{split}
&     \Delta T_{n+1} = \frac{\theta_{n+1}}{\lambda_0} - \frac{b^T {\tau_{\eps}}}{\lambda_0}(Z_{T_{n+1}}-Z_{T_{n}})+ \epsilon \Delta T _{0 } b^T {\tau_{\eps}} \nabla S (X_{T_n}) \mathcal{V}_n \\
& + \pare{\theta_{n+1}^3+\theta_{n+1}}\mathcal{O} \pare{ \epsilon\, {\rm nl_{2}} (\epsilon) + \epsilon^{1+\delta} }.
  \end{split}
 \end{equation}
\end{lem}
\begin{proof}
  The proof is based on Duhamel integration of the ODE~\eqref{eq:process_noscale} on $[T_n,T_{n+1}]$, as in the proof of Lemma 4.7 of the companion paper~\cite{limits}. Following equation~(4.24) in~\cite{limits}, we get for $t\in [T_n,T_{n+1}]$~:
\begin{align}\label{eq:key_Duhamel}
  Z_t = {\rm e}^{- \pare{t-T_n}{{\tau_{\eps}}^{-1}}} Z_{T_n} + \epsilon \pare{\Id - {\rm e}^{- \pare{t-T_n}{{\tau_{\eps}}^{-1}}} }{{\tau_{\eps}}} \nabla S(X_{T_n})  \mathcal{V}_{n} + (\theta_{n+1}^2+\theta_{n+1})\mathcal{O} \pare{c_S \epsilon^2+ c_F \epsilon^{1+\delta} }.
\end{align}
Multiplying both sides by $\tau_\epsilon$, recalling that Assumption~\ref{a:2} implies $\|\tau_\epsilon\|=\mathcal{O}(\epsilon^{\delta-1})$, and using Assumption~\ref{a:1bis}, we get for $t\in [T_n,T_{n+1}]$~:
\begin{align}\label{e:help}
 {\tau_{\eps}} Z_t = {\tau_{\eps}} {\rm e}^{- \pare{t-T_n}{{\tau_{\eps}}^{-1}}} Z_{T_n} + \epsilon \pare{\Id - {\rm e}^{- \pare{t-T_n}{{\tau_{\eps}}^{-1}}} }{{\tau_{\eps}}^2} \nabla S(X_{T_n})  \mathcal{V}_{n} + (\theta_{n+1}^2+\theta_{n+1})\mathcal{O} \pare{c_S \epsilon^{1+\delta}+ c_F \epsilon^{1+\gamma} }.
\end{align}
Moreover, integrating~\eqref{eq:key_Duhamel} on $t\in [T_n,T_{n+1}]$ yields
\begin{align}\label{eq:key_Duhamelint}
  \begin{split}
   \int_{T_n}^{T_{n+1}} Z_t \d t &= \pare{\Id - {\rm e}^{- \pare{T_{n+1}-T_n}{{\tau_{\eps}}^{-1}}} }{{\tau_{\eps}}} \pare{Z_{T_n} -\epsilon{\tau_{\eps}} \nabla S(X_{T_n}) \mathcal{V}_{n} } \\
& \quad +\epsilon (T_{n+1}-T_n){\tau_{\eps}} \nabla S(X_{T_n})  \mathcal{V}_{n} + (\theta_{n+1}^3+\theta_{n+1}^2)\mathcal{O} \pare{c_S \epsilon^2+ c_F \epsilon^{1+\delta} }.
\end{split}
\end{align}
Evaluating equation~\eqref{e:help} at $t=T_{n+1}$, add adding equation~\eqref{eq:key_Duhamelint}, we get
\begin{align*}
  \int_{T_n}^{T_{n+1}} Z_t \d t + {\tau_{\eps}} Z_{T_{n+1}}  = & {\tau_{\eps}} Z_{T_n} +\epsilon (T_{n+1}-T_n){\tau_{\eps}} \nabla S(X_{T_n})  \mathcal{V}_{n}  \\
& \quad  + \pare{\theta_{n+1}^3+\theta_{n+1}}\mathcal{O} \pare{c_S(\epsilon^{1+\delta}+\epsilon^2 ) + c_F(\epsilon^{1+\delta}+\epsilon^{1+\gamma}) } .
\end{align*}
The estimates in Lemma~$4.8$ from~\cite{limits} (i.e., $\abs{ \Delta T_{n+1}-\Delta T_{n+1}^0} = (\theta_{n+1}^4+\theta_{n+1})\mathcal{O}(\epsilon) $) yields
\begin{align*}
  \int_{T_n}^{T_{n+1}} Z_t \d t + {\tau_{\eps}} Z_{T_{n+1}}  = & {\tau_{\eps}} Z_{T_n} +\epsilon \Delta T_{n+1}^0 {\tau_{\eps}} \nabla S(X_{T_n})  \mathcal{V}_{n}  \\
& \quad  + \pare{\theta_{n+1}^3+\theta_{n+1}}\mathcal{O} \pare{c_S(\epsilon^{1+\delta}+\epsilon^2 ) + (c_F+1)\epsilon^{1+\delta}+c_F\epsilon^{1+\gamma}) } .
\end{align*}

Finally, plugging the last equation in the estimate~\eqref{eq:lin_rate} yields the result.
\end{proof}

\subsection{Analysis of variance of the coupling\label{sec:proof}}
The present section is devoted to the proof of the following theorem:
\begin{thm}\label{thm:coupling}
Assume the assumptions of Section~\ref{sec:asympt} hold, and that $k \delta > 1$ where $k\geq 2$ is defined in~\eqref{eq:lin_rate}. Then the difference between the process with internal state and the coupling process satisfies~:
\begin{equation}
  \label{eq:coupl_estim}
   \E \pare{\pare { X_{{\bar{t}}/\epsilon^2} - X^c_{{\bar{t}}/\epsilon^2}  } ^ p}^{1/p} = \mathcal{O}( \epsilon + \epsilon^\delta +  {\rm nl_{2}} (\epsilon)  ), \qquad p\ge 1,
\end{equation}
where $ {\rm nl_{2}} (\epsilon)$ is defined in~\eqref{eq:nl_err}.
\end{thm}

The proof of Theorem~\ref{thm:coupling} relies on a number of steps that will be detailed by a series of lemmas. We will make use of the following random sequences, defined as the position of both processes after $n$ jumps,
\begin{equation}
  \label{eq:Xi}
  \Xi_n := X_{T_n}= X_0 + \epsilon \sum_{m=0}^{n-1}\Delta T _{m+1} \mathcal{V}_m, \qquad n\geq 0,
\end{equation}
as well as
\begin{equation}
  \label{eq:Xic}
  \Xi^c_n := X_{T^c_n}^c= X_0 + \epsilon \sum_{m=0}^{n-1}\Delta T _{m+1}^c \mathcal{V}_m, \qquad n\geq 0.
\end{equation}
We  will also make use of the following random integers~:
\begin{defn}\label{def:N_Nc}
  The random integers $N \geq 1$ and $N^c \geq 1$ are uniquely defined by:
\begin{equation}
  T_N \leq {\bar{t}} / \epsilon^2 < T_{N+1}, \qquad T^c_{N^c} \leq {\bar{t}} / \epsilon^2 < T^c_{N^c+1} .
\end{equation}
\end{defn}

The first lemma bounds the difference at time ${\bar{t}} / \epsilon^2$ between the two processes $X_{{\bar{t}}/\epsilon^2}$ and $X^c_{{\bar{t}}/\epsilon^2}$ by expressing it in terms of  differences of positions and jump times of both processes after the same random number of jumps, $\Xi_n$ and $\Xi_n^c$~:
\begin{lem}
  The difference between the rescaled process with internal state $X^\epsilon_{{\bar{t}}}:= X_{{\bar{t}}/\epsilon^2}$ and the rescaled coupling process $X^{c,\epsilon}_{{\bar{t}}}:= X^c_{{\bar{t}}/\epsilon^2}$ satisfies:
  \begin{align}
\abs{X^\epsilon_{{\bar{t}}} - X^{\epsilon,c}_{{\bar{t}}}} & \leq  \abs{\Xi_{N}-\Xi_{N}^c} + \epsilon \abs{T_{N}-T_{N}^c}+\pare{\theta_{N+1}+\theta_{N_c+1}} \mathcal{O}(\epsilon) \label{eq:coupl_11} \\
&  \leq \abs{\Xi_{N^c}-\Xi_{N^c}^c} + \epsilon \abs{T_{N^c}-T_{N^c}^c}+\pare{\theta_{N+1}+\theta_{N_c+1}} \mathcal{O}(\epsilon) \label{eq:coupl_12}
\end{align}
\end{lem}
Only one of the two estimates \eqref{eq:coupl_11}-\eqref{eq:coupl_12} is necessary in the remainder of the proof, but we detail both to highlight the symmetry.
\begin{proof}
  By definition, we have
  \begin{equation*}
    \begin{cases}
      X^\epsilon_{{\bar{t}}}= \Xi_N+ \epsilon({\bar{t}}/\epsilon^2-T_N)\mathcal{V}_{N}, \\
      X^{c,\epsilon}_{{\bar{t}}}= \Xi_{N^c}^c+ \epsilon({\bar{t}}/\epsilon^2-T_{N^c}^c)\mathcal{V}_{N},
    \end{cases}
  \end{equation*}
so that by Definition~\ref{def:N_Nc} of $N$ and $N^c$, and by realizing that $\abs{{\bar{t}}/\epsilon^2-T_{N^c}^c} \leq C \theta_{N_c+1}$ and $\abs{{\bar{t}}/\epsilon^2-T_{N}} \leq C \theta_{N+1}$,  we get
\begin{equation}\label{eq:diffX}
    \begin{cases}
     \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}}\leq \abs{\Xi_N-\Xi^c_N}+\underbrace{\abs{\Xi^c_{N^c}-\Xi^c_{N}}}_{(a)}+\pare{\theta_{N+1}+\theta_{N_c+1}} \mathcal{O}(\epsilon)  , \\
     \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}} \leq \abs{\Xi_{N^c}-\Xi^c_{N^c}}+\underbrace{\abs{\Xi_{N^c}-\Xi_{N}}}_{(b)}+ \pare{\theta_{N+1}+\theta_{N_c+1}} \mathcal{O}(\epsilon).
    \end{cases}
  \end{equation}
To analyze $(a)$, we consider three cases:
\begin{itemize}
\item[$N=N_c$.] Then $(a)=0$.
\item[$N > N_c$.] Then
  \begin{align*}
    \frac{(a)}{\epsilon} & = \abs{\sum_{n = N_c}^{N-1} \Delta T _{n+1}^c \mathcal{V}_n } \\
& \leq T _{N}^c-T_{N_c}^c=\Delta T^c_{N^c+1} + T^c_ N - T^c_{N^c+1}  \\
&\leq T^c_ N - T_{N} + C \theta_{N^c+1},
  \end{align*}
where in the last line we have used that, $\Delta T^c_{N^c+1}\le C\theta_{N^c+1}$, and, by definition, $ T_{N} \leq {\bar{t}}/ \epsilon^2 < T^c_{N^c+1}$.
\item[$N_c > N$.] Then
  \begin{align*}
    \frac{(a)}{\epsilon} & = \abs{\sum_{n = N}^{N_c-1} \Delta T _{n+1}^c \mathcal{V}_n }  \\
& \leq T _{N_c}^c-T_{N}^c \\
&\leq T_{N+1} - T_{N}^c  = \Delta T_{N+1} + T_N-T_N^c \leq T_ N - T^c_{N} + C \theta_{N+1} ,
  \end{align*}
where in the last line we have used that by definition
 $ T^c_{N^c} \leq {\bar{t}}/ \epsilon^2 \leq T_{N+1}$.
\end{itemize}
This proves~\eqref{eq:coupl_11}. By symmetry, we get:
\[
\frac{(b)}{\epsilon} \leq \abs{T_{N^c} - T^c_{N^c}} + C\pare{\theta_{N+1}+\theta_{N^c+1}},
\]
which yields~\eqref{eq:coupl_12}.
\end{proof}

In the next lemma, we estimate the supremum of the difference of the position of both processes after $n$ jumps for $n\in [0,n_\epsilon]$.
\begin{lem}
  Let $n_\epsilon \geq 1$ be a deterministic integer verifying $n_\epsilon = \mathcal{O}(\epsilon^{-2})$. Then
  \begin{equation}
    \label{eq:estim_jumps}
    \E \pare{ \sup_{0 \leq n \leq n_\epsilon}\abs{\Xi_n-\Xi_n^c}^p } ^{1/p} = \mathcal{O}(\epsilon + \epsilon^\delta +  {\rm nl} (\epsilon) ).
  \end{equation}
\end{lem}
\begin{proof}
  First, using the estimates of jump times~\eqref{eq:deltaTc}-\eqref{eq:deltaT}, we can decompose the differences of positions of both processes as follows~:
  \begin{align*}
    \Xi_{n}-\Xi_{n}^c  &= \sum_{m=0}^{n-1}\epsilon (\Delta T _ {m+1}-\Delta T _ {m+1} ^c) \mathcal{V}_{m} \\
& = \sum_{m=0}^{n-1} \alpha_{m+1} + d_{m+1} + r_{m+1},
  \end{align*}
where, by definition,
\begin{equation*}
  \begin{cases}
    \alpha_{m+1} := \epsilon^2 \frac{\mathcal{V}_m}{\lambda_0^2} \pare{A_\epsilon(\Xi_{m})^T-A_\epsilon(\Xi_{m}^c)^T } \mathcal{V}_m, \\
       d_{m+1} :=\epsilon^2\frac{\mathcal{V}_m}{\lambda_0^2} \pare{1-\theta_{m+1}}A_\epsilon(\Xi_{m}^c)^T  \mathcal{V}_m \\
\hspace{1.5cm} + \epsilon^2 \frac{\mathcal{V}_{m}}{\lambda_0^2} \pare{ \lambda_0  b^T m\pare{\frac{\theta_{m+1}}{\lambda_0}} \nabla S(\Xi_m) - A_\epsilon(\Xi_{m})^T }  \mathcal{V}_{m} \\
\hspace{1.5cm} + \epsilon b^T m'(\Delta T_{m+1}^0 ) Z_{T_m} \mathcal{V}_m , \\
r_{m+1}:= \pare{\theta_{m+1}^6+\theta_{n+1}} \mathcal{O}(\epsilon^{3} + \epsilon^2  {\rm nl} (\epsilon) ).
     \end{cases}
  \end{equation*}

Since $A_\epsilon$ is Lipschitz, we have that $\alpha_{m+1} \leq C \epsilon^2 \abs{\Xi_{m}-\Xi_{m}^c}$, and we can write
\begin{equation*}
  \abs{\Xi_{n}-\Xi_{n}^c} \leq  C\epsilon^2 \sum_{m=0}^{n-1} \abs{\Xi_{m}-\Xi_{m}^c} + \sup_{0 \leq l \leq n_{\epsilon}}\abs{\sum_{m=1}^{l} d_{m}}+ \sum_{m=1}^{n_\epsilon}\abs{r_{m}},
\end{equation*}
for $n \leq n_{\epsilon}$, which yields
\begin{align}
  \abs{\Xi_{n}-\Xi_{n}^c} &\leq \pare{\sup_{0 \leq l \leq n_{\epsilon}}\abs{\sum_{m=1}^{l} d_{m}}+ \sum_{m=1}^{n_\epsilon}\abs{r_{m}}} \pare{1+C \epsilon^2+ \dots + (C \epsilon^2)^{n_\epsilon}}, \label{eq:key_xi}.
\end{align}
Note that
\[
\pare{1+C \epsilon^2+ \dots + (C \epsilon^2)^{n_\epsilon}} \leq  C .
\]
Now, we can remark that the discrete time process
$
\pare{M_l}_{l \geq 0}:=\pare{ \sum_{m=1}^{l} d_{m} }_{l \geq 0}
$ is a martingale with respect to the filtration $\mathcal F _{l} = \sigma\pare{(\theta_{m},\mathcal{V}_{m-1}) \vert 1 \leq m \leq l }$ with $M_0=0$, so that we can apply the Burkholder-Davies-Gundy upper bound (which is a simple consequence of Doob maximal inequality here, see~\cite{Karatzasbook}), for some $p \geq 1$:
\begin{equation}
  \label{eq:BDG}
  \E \pare{\sup_{0 \leq l \leq n_\epsilon } \abs{M_l}^p} \leq C_p n_{\epsilon}^{p/2 - 1} \sum_{l=1}^{n_{\epsilon}} \E \pare{\abs{M_l-M_{l-1}}^p}.
\end{equation}
In the present case, since $\abs{Z_{T_n}}=\mathcal{O}(\epsilon^\delta)$:
\[
\E \pare{\abs{d_l}^p} \leq C_p \epsilon^{p(1+\delta)}.
\]
Moreover,
\begin{align*}
  \E \pare{\abs{\sum_{m=1}^{n_\epsilon}\abs{r_{m}}}^p} &\leq n_{\epsilon}^{p - 1} \sum_{m=1}^{n_\epsilon}\E \abs{r_{m}}^p \\
& \leq C \pare{ \epsilon +  {\rm nl} (\epsilon) }^p,
\end{align*}
and the result follows from~\eqref{eq:key_xi}.
\end{proof}
In the same way, we can estimate the difference of the $n_\epsilon$-th jump time of both processes~:
\begin{lem}
  Let $n_\epsilon \geq 1$ a deterministic integer. Then
  \begin{equation}
    \label{eq:estim_jumps_T}
    \epsilon \, \E \pare{ \sup_{0 \leq n \leq n_\epsilon}\abs{T_n-T_n^c}^p }^{1/p} \leq C_p \pare{ \epsilon^\delta +  {\rm nl_{2}} (\epsilon)  }
  \end{equation}
\end{lem}
\begin{proof}
   First, using the estimates of jump times~\eqref{eq:deltaTc}-\eqref{eq:deltaTbis}, we can decompose the differences of jump times of both processes as follows~:
  \begin{align*}
    T_{n}-T_{n}^c  &= \sum_{m=0}^{n-1} \Delta T _ {m+1}-\Delta T _ {m+1} ^c\\
& = \sum_{m=0}^{n-1} \delta_{m+1} + d_{m+1} + r_{m+1},
  \end{align*}
where, by definition,
\begin{equation*}
  \begin{cases}
    \delta_{m+1} := - \frac{b^T {\tau_{\eps}}}{\lambda_0}(Z_{T_{m+1}}-Z_{T_{m}}), \\
       d_{m+1} :=\epsilon \Delta T _{0 } b^T {\tau_{\eps}} \nabla S (X_{T_n}) \mathcal{V}_n \\
\hspace{1.5cm} - \epsilon \frac{\theta_{n+1}}{\lambda_0^2} A_\epsilon^T(X_{T^c_n}^c) \mathcal{V}_n  \\
r_{m+1}:= \pare{\theta_{m+1}^3+ \theta_{n+1} } \epsilon \mathcal{O}(\epsilon^\delta +  {\rm nl_{2}} (\epsilon)).
     \end{cases}
  \end{equation*}
We can write for $n \leq n_{\epsilon}$:
\begin{align*}
  \abs{T_{n}-T_{n}^c} \leq \abs{\frac{b^T {\tau_{\eps}}}{\lambda_0}(Z_{T_{n}}-Z_{T_{0}}) } + \sup_{0 \leq l \leq n_{\epsilon}}\abs{\sum_{m=1}^{l} d_{m}} + \sum_{m=1}^{n_\epsilon}\abs{r_{m}}.
\end{align*}
Now, we can remark that the discrete time process
$
\pare{M_l}_{l \geq 0}:=\pare{ \sum_{m=1}^{l} d_{m} }_{l \geq 0}
$ is a martingale with respect to the filtration $\mathcal F _{l} = \sigma\pare{(\theta_{m},\mathcal{V}_{m-1}) \vert 1 \leq m \leq l }$ with $M_0=0$. Using the Burkholder-Davies-Gundy inequality~\eqref{eq:BDG}, together with:
\[
\E \pare{\abs{m_{m}}^p} \leq C \epsilon^{p\delta},
\]
yields
\[
\E \pare{ \sup_{0 \leq l \leq n_{\epsilon}}\abs{\sum_{m=1}^{l} d_{m}}^p  }^{1/p} \leq C \epsilon^{\delta-1}.
\]
Using:
\[
\abs{\frac{b^T {\tau_{\eps}}}{\lambda_0}(Z_{T_{n}}-Z_{T_{0}}) } = \mathcal{O}(\epsilon^{2 \delta-1}) < \mathcal{O}(\epsilon^{\delta-1}),
\]
and plugging in the rest term, we finally get~\eqref{eq:estim_jumps_T}.
\end{proof}

We can now conclude with the proof of Theorem~\ref{thm:coupling}~:

\begin{proof}[Proof of Theorem~\ref{thm:coupling}]
We consider a given integer $n_{\epsilon} \geq 1$. We decompose the difference between both processes as
\begin{equation}\label{eq:proof1}
  \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}} = \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}} {\bf 1 }_{N\leq n_{\epsilon}}+\abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}} {\bf 1 }_{N > n_{\epsilon}},
\end{equation}
and use~\eqref{eq:coupl_11} to write
\begin{equation}
  \label{eq:proof2}
  \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}} {\bf 1 }_{N\leq n_{\epsilon}} \leq \sup_{0 \leq n \leq n_\epsilon}\abs{\Xi_n-\Xi_n^c} +\epsilon \sup_{0 \leq n \leq n_\epsilon}\abs{T_n-T_n^c} + \pare{\theta_{N+1}+\theta_{N_c+1}} \mathcal{O}(\epsilon).
\end{equation}
Using now the estimates~\eqref{eq:estim_jumps} and \eqref{eq:estim_jumps_T} in~\eqref{eq:proof2}, we get
\begin{align}\label{eq:estim_jumps_tot}
  \E \pare{ \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}}^p {\bf 1 }_{N \leq n_{\epsilon}} }^{1/p}& \leq C  ( \epsilon + \epsilon^\delta +  {\rm nl_{2}} (\epsilon)  )\\
& \qquad + C\epsilon \E (\theta_{N}^p+\theta_{N^c}^p)^{1/p} \nonumber \\
& = \mathcal{O}( \epsilon + \epsilon^\delta +   {\rm nl_{2}} (\epsilon)     ),
\end{align}
where in the last line we have used a technical lemma (Lemma~\ref{lem:tech1}).

It remains to control the probability of the event $\{ N > n_\epsilon\}$. To this purpose, we consider $t \mapsto N_t \in \mathbb N$ the Poisson process uniquely defined by:
\[
\sum_{n=1}^{N_t} \theta_n \leq t < \sum_{n=1}^{N_t +1} \theta_n.
\]
By definition of $N$ and $N_c$, there is a constant $C$ such that $N_{C{\bar{t}} / \epsilon^2} \geq N$ and $N_{C{\bar{t}} / \epsilon^2} \geq N^c$, so that we can write:
\begin{equation}
  \label{eq:proof3}
  \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}}  \leq \epsilon C \pare{\sum_{n=1}^{N_{C{\bar{t}} / \epsilon^2}+1}\theta_{n}}{\bf 1 }_{N_{C{\bar{t}} / \epsilon^2} > n_{\epsilon}}.
\end{equation}
Let us choose $n_{\epsilon}  \geq  4 {\rm e} C{\bar{t}} / \epsilon^2 $, a choice that  satisfies $n_{\epsilon}=\mathcal{O}(\epsilon^{-2})$. A second technical lemma (Lemma~\ref{lem:tech2}) implies that
\begin{equation}
  \label{eq:proof4}
  \E \pare{ \abs{ X^\epsilon_{{\bar{t}}}- X^{c,\epsilon}_{{\bar{t}}}}^p {\bf 1 }_{N > n_{\epsilon}} } \leq \epsilon C_p 2^{-1/(C\epsilon^2)},
\end{equation}
so that the contribution of the latter to the final estimate~\eqref{eq:coupl_estim} is negligible, since it is exponentially small with respect to $\epsilon^{-2}$. The proof is complete.
\end{proof}

\subsection{Discussion about the sharpness of the coupling}\label{sec:sharp}
One may ask about the ``sharpness'' of the precise estimates of the coupling, given mainly by Theorem~\ref{thm:coupling}; (see also~\eqref{eq:estim_jumps} (position shift), and \eqref{eq:estim_jumps_T} (time shift)). Let us discuss the case where $c_F=0$, $ \delta \in [1/k , 1/(k-1)]$. In this case the dominant error term in the jump times expansion of the internal state model in Lemma~$4.8$ in~\cite{limits} is due to the non-linearity of the turning rate: $c_\lambda \epsilon^{k \delta -1}$. Since this term is due to the internal state mechanism, it can be conjectured that the difference between the density of the internal state model and gradient sensing model will be at least of this order. On the other hand, the coupling estimate in Theorem~\ref{thm:coupling} is controlled by the same term: $c_\lambda \epsilon^{k \delta -1}$ , suggesting the sharpness of the latter.

\section{Numerical illustration\label{sec:illustr}}

In this section, we  demonstrate the validity of the analysis above. For the numerical experiments, we restrict ourselves to one space dimension.
For the process with internal state, we use internal dynamics that are given by the scalar cartoon model \eqref{e:scalar-y}. In all experiments, $\tau$ is chosen independently of $\epsilon$ so that $\delta=1$. The corresponding turning rate is given by \eqref{eq:rate} with $\lambda_0=1$. For the control process, we choose the linear turning rate \eqref{eq:control_rate}, with parameters \eqref{eq:control_field} such that the two processes have the same diffusion limit; in particular, $b=1$. The physical domain is $x\in [0,20]$, and we use reflecting boundary conditions, {\it i.e}. the bacterial velocity is reversed when $x=0$ or $x=20$. We fix a scalar bimodal chemoattractant concentration field
\begin{equation}\label{eq:chemoattractant}
	S(x)=\alpha\left(\exp\left(-\beta\left(x-\xi\right)^2\right)+\exp\left(-\beta\left(x-\eta\right)^2\right)\right),
\end{equation}
in which $\alpha=5$, $\beta=1$, $\xi=7.5$ and $\eta=12.5$.

\paragraph{Single bacterium as a function of time. } In a first experiment, we simulate a single bacterium evolving according to the fine-scale model (\ref{eq:process_noscale}), in which we set the parameters to $\epsilon=0.2$ and $\tau=1$.  We compare this evolution to that of a bacterium that satisfies the corresponding control process \eqref{eq:cprocess}. Both simulations are performed using the same random numbers starting from the initial condition $X_0=X_0^c=8$ and $V_0=V_0^c=+1$; the bacterium with internal state starts from $Y_0=S(X_0)$. The time step $\delta t=0.1$. The results are shown in figure \ref{fig:one-particle}.
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.8]{figuren/fig-one-particle}
	\caption{\label{fig:one-particle}Evolution of a single bacterium evolving according to the fine-scale process (\ref{eq:process_noscale}) (solid line) and the corresponding control process (\ref{eq:cprocess}) (dashed) on short (left) and long (right) time-scales.}
\end{center}
\end{figure}
We see a very good coupling initially, which degrades over time. Note the time shift in the short time picture. In the long time picture, coupling is completely lost at some point.

\paragraph{Expectation and variance as a function of $\epsilon$ and ${\bar{t}}$.} Next, we repeat the experiment using $N=10000$ particles and compute the empirical mean and variance of the coupling, i.e.,\
\[
E\left(\left|X_{{\bar{t}}}-X_{{\bar{t}}}^c\right|\right)=\frac{1}{N}\sum_i^N \left(\left|X^i_{{\bar{t}}}-X^{i,c}_{{\bar{t}}}\right|\right), \;\;\; \text{resp., }  \;\;\; E\left((X_{{\bar{t}}}-X_{{\bar{t}}}^c)^2\right)=\frac{1}{N}\sum_i^N \left(X^i_{{\bar{t}}}-X^{i,c}_{{\bar{t}}}\right)^2.
\]
As fine-scale parameters, we choose $\tau=1$, $\lambda_0=b=1$ and several values of $\epsilon$. The chemoattractant concentration is again given as \eqref{eq:chemoattractant}, now with $\alpha=\beta=1$, $\xi=7.5$ and $\eta=12.5$.  The time interval for the computation is $t\in [0,30/\epsilon^2]$.
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.68]{figuren/fig-eps-long}
	\caption{\label{fig:eps-long} Empirical mean (left) and variance (right) of the difference between the fine-scale process (\ref{eq:process_noscale}) and the corresponding control process (\ref{eq:cprocess}) as a function of $\epsilon$ for different values of the reporting time $T$.  The theoretical slope is indicated with a dashdotted line. The sample size is $N=10000$.}
\end{center}
\end{figure}
Figure \ref{fig:eps-long} shows the dependence in $\epsilon$ of the coupling, by plotting the empirical mean and variance defined above as a function of $\epsilon$ for different values of the reporting time $T$. The results shown in figure \ref{fig:eps-long} are in clear accordance with the theoretical slope predicted by the asymptotic analysis.
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.68]{figuren/fig-T-long}
	\caption{\label{fig:T-long}Evolution of the empirical mean (left) and variance (right) of the difference between the fine-scale process (\ref{eq:process_noscale}) and the corresponding control process (\ref{eq:cprocess}) as a function of ${\bar{t}}=t/\epsilon^2$ for different values of $\epsilon$. The sample size is $N=10000$. }
\end{center}
\end{figure}
Figure \ref{fig:T-long} shows the mean and variance of the coupling difference as a function of time.  The time dependence of the variance has not been analyzed mathematically in Section~\ref{sec:variance}; the specific behaviour viewed in figure \ref{fig:T-long} and is probably due to: (i) sufficiently short diffusive times; (ii) the specific double-well form of the chemoattractant potential.

\paragraph{Expectation and variance as a function of $\tau$.} Finally, in a last experiment, we illustrate the dependence on $\tau$. To this end, we again simulate $N=10000$ particles choosing $\lambda_0=1$, $\epsilon=0.1$ and $X_0=7.5$, for different values of $\tau$. The results in figure \ref{fig:tau-long} show that the variance quickly increases with $\tau$ until it reaches a plateau for $\tau>1$.
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.68]{figuren/fig-tau-long}
	\caption{\label{fig:tau-long} Empirical mean (left) and variance (right) of the difference between the fine-scale process (\ref{eq:process_noscale}) and the corresponding control process (\ref{eq:cprocess}) as a function of $\tau$ for different values of the reporting time $T$. The sample size is $N=10000$.}
\end{center}
\end{figure}

\section{Simulation on diffusive time scales\label{sec:appl}}

In this section, we consider a simulation of the density of an ensemble of particles, with and without variance reduction and/or reinitialization, as described in Section~\ref{sec:coupl_princ}. We again restrict ourselves to one space dimension, with domain $x\in [0,20]$ and periodic boundary conditions. In this case, the kinetic equation corresponding to the control process reduces to the system
\begin{equation} \label{eq:cprocess_hydro1d}
\system{
	 & \partial_t  p^c_++\epsilon\partial_x p^c_+=-\frac{\lambda^c(x,+1)}{2}p_+^c+ \frac{\lambda^c(x,-1)}{2}p_-^c &\\
	 & \partial_t  p^c_--\epsilon\partial_x p^c_-=\frac{\lambda^c(x,+1)}{2}p_+^c -\frac{\lambda^c(x,-1)}{2}p_-^c &
	 }.
 \end{equation}
of two PDEs, which is straightforward to simulate using finite differences.

We fix the chemoattractant concentration field as \eqref{eq:chemoattractant}, with parameters $\alpha=2$, $\beta=1$, $\xi=7.5$ and $\eta=12.5$. For the internal dynamics, the same model (\eqref{e:scalar-y}-\eqref{eq:rate}, \eqref{eq:control_rate}-\eqref{eq:control_field}) is used. The parameters are $\epsilon=0.5$, $\lambda_0=1$, $\tau=1$, $\delta t=0.1$.

All simulations are performed with $N=5000$ particles. The initial positions are uniformly distributed in the interval $x\in[13,15]$; the initial velocities are chosen uniformly, i.e., each particle has an equal probability of having an initial velocity of $\pm\epsilon$. The initial condition for the internal variable is chosen to be in local equilibrium, i.e., $Y_0^i=S(X^i_0)$. The initial positions and velocities of the control particles are chosen to be identical.

We discretize the continuum description \eqref{eq:cprocess_hydro1d} on a mesh with $\Delta x = 0.1$ using a third-order upwind-biased scheme, and perform time integration using the standard fourth order Runge--Kutta method with time step $\delta t_{pde}= 10^{-1}$. The initial condition is given as
\begin{equation}
	p^+(x,0)=p^-(x,0)=\begin{cases}
	0.25 , & {x\in [13,15],} \\
	0, & \text{otherwise.}
\end{cases}
\end{equation}

\paragraph{Simulation without variance reduction.} First, we simulate both stochastic processes up to time $\bar{t}=50$ ($t=50/\epsilon^2$) and estimate the density of each of these processes $\hat{n}_N(x,{\bar{t}})$, resp.~$\hat{n}^c_N(x,{\bar{t}})$, without variance reduction.  The density is obtained via binning in a histogram, in which the grid points of the deterministic simulation are the centers of the bins. Figure \ref{fig:no-variance-reduction} (left)
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.8]{figuren/fig-no-variance-reduction}
	\caption{\label{fig:no-variance-reduction} Bacterial density as a function of space at $t=50/\epsilon^2$ without variance reduction. Left: one realization. Right: mean over $100$ realizations and $95\%$ confidence interval. The solid line is the estimated density from a particle simulation using the process with internal state; the dashed line is estimated from a particle simulation using the control process. Both used $N=5000$ particles. The dotted line is the solution of the deterministic PDE (\ref{eq:cdensity}).   }
\end{center}
\end{figure}
shows the results for a single realization.  We see that, given the fluctuations on the obtained density, it is impossible to conclude on differences between the two models.  This observation is confirmed by computing the average density of both processes over $100$ realizations.  The mean densities are shown in figure \ref{fig:no-variance-reduction} (right), which also reveals that the mean density of the control process is within the $95\%$ confidence interval of the process with internal state.  Both figures also show the density that is computed using the continuum description, which coincides with the mean of the density of the control particles.

\paragraph{Simulation with variance reduction.} Next, we compare the variance reduced estimation \eqref{eq:var_reduced_estimation} with the density of the control PDE. We reinitialize the control particles after each coarse-scale step, i.e., each $k$ steps of the particle scheme, where $k \delta t = \delta t_{pde}$, (here $k=1$). The results are shown in figure \ref{fig:var-red-reinit}.
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.8]{figuren/fig-variance-reduction}
	\caption{\label{fig:var-red-reinit} Bacterial density as a function of space at $t=50/\epsilon^2$ with variance reduction and reinitialization. Left: variance reduced density estimation of one realization with $N=5000$ particles (solid) and deterministic solution for the control process (\ref{eq:cdensity}) (dashed). Right: mean over $100$ realization and $95\%$ confidence interval (solid) and the deterministic solution for the control process (\ref{eq:cdensity}) (dashed).}
\end{center}
\end{figure}
We see that, using this reinitialization, the difference between the behaviour of the two processes is visually clear from one realization (left figure).  Also, the resulting variance is such that the density of the control PDE is no longer within the $95\%$ confidence interval of the variance reduced density estimation (right figure). We see that there is a significant difference between both models: the density corresponding to the control process is more peaked, indicating that bacteria that follow the control process are more sensitive to sudden changes in chemoattractant gradient. This difference can be interpreted from the fact that the bacteria with internal state do not adjust themselves instantaneously to their environment, but instead with a time constant $\tau$.

\paragraph{Simulation with variance reduction but without reinitialization. } We also compare the variance reduced estimation \eqref{eq:var_reduced_estimation} with the density of the control PDE without performing any reinitialization of the control process to restore the coupling.
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.8]{figuren/fig-var-red-noreinit}
	\caption{\label{fig:var-red-noreinit} Bacterial density as a function of space at $t=50/\epsilon^2$ with variance reduction. Left: variance reduced density estimation of one realization with $N=5000$ particles (solid) and deterministic solution for the control process (\ref{eq:cdensity}) (dashed). Right: mean over $200$ realization and $95\%$ confidence interval (solid) and the deterministic solution for the control process (\ref{eq:cdensity}) (dashed).}
\end{center}
\end{figure}
Figure \ref{fig:var-red-noreinit} (right) shows the mean estimation of the density over $200$ realizations, as well as the density of the control process. As evidenced by the error bars, simulating a single realization of the process, even with variance reduction, is not able to reliably reveal this difference. This phenomenon is due to the degeneracy of the coupling on long diffusive times. This is also illustrated in figure \ref{fig:var-red-noreinit} (left), which compares the variance reduced density estimation of a single realization with the density of the control PDE. Note that, as predicted by the analysis, the variance is much larger in regions where $\nabla S(x)$ is large.

Finally, we repeat the experiment with a larger chemoattractant gradient. We again choose the chemoattractant field as \eqref{eq:chemoattractant}, with the same parameters as above, except that we now take $\alpha=5$. We also use the same discretization parameters as above. The variance reduced estimation \eqref{eq:var_reduced_estimation}, as well as the density of the control PDE, are shown in figure \ref{fig:var-red-reinit-alpha5}.
\begin{figure}
	\begin{center}
	\includegraphics[scale=0.8]{figuren/fig-variance-reduction-alpha5}
	\caption{\label{fig:var-red-reinit-alpha5} Bacterial density as a function of space at $t=50/\epsilon^2$ with variance reduction and reinitialization. Left: variance reduced density estimation of one realization with $N=5000$ particles (solid) and deterministic solution for the control process (\ref{eq:cdensity}) (dashed). Chemoattractant given by \eqref{eq:chemoattractant} with $\alpha=5$. Right: mean over $100$ realization and $95\%$ confidence interval (solid) and the deterministic solution for the control process (\ref{eq:cdensity}) (dashed).}
\end{center}
\end{figure}
We see that, although the difference between the process with internal state and the control process becomes much larger, the variance is still very well controlled by the coupling.

\section{Conclusions and discussion\label{sec:concl}}

In this paper, we studied the simulation of stochastic individual-based models for chemotaxis of bacteria with internal state. We have used a coupling with a simpler, direct gradient sensing model with the same diffusion limit to obtain an ``asymptotic'' variance reduction.  Throughout this work, we have assumed that the computation of the bacterial density using the kinetic description of the direct gradient sensing model can be performed accurately with a grid-based method. Note that in higher spatial dimensions, this may become cumbersome, and the grid-based computation should be carried out at the level of the advection-diffusion description (or any desired moment system) only, using a second level of coupling between the velocity-jump process associated with the kinetic description, and the stochastic differential equation associated with the advection-diffusion limit. This is left for future work.

The current algorithm allows to explore the differences between the fine-scale process with internal state and the simpler, coarse process \emph{using a fixed number of particles, independently of the small parameter $\epsilon$}.  However, the computational cost to simulated an individual particle over diffusive time-scales becomes very large for $\epsilon\to 0$.  In future work, we will therefore also study the development of truly ``asymptotic preserving'' schemes in the diffusion asymptotics, in the sense that the computational cost of the simulation of processes is independent of the small parameter $\epsilon$.
This will require to deal with two kinds of difficulties: (i)
 We will need to use an asymptotic preserving method to solve the density evolution of the control model on the grid; and (ii) We will need to extrapolate forward in time the state of the fine-scale simulation.

The first difficulty implies that, instead of solving the full kinetic equation associated with the control process, we only solve its diffusion limit. This may imply the use of a second level of variance reduction, coupling the control velocity jump process and its limiting drift-diffusion process.

The second difficulty, extrapolation in time, is related to the equation-free \cite{KevrGearHymKevrRunTheo03,Kevrekidis:2009p7484} and HMM \cite{EEng03,E:2007p3747} types of methodologies; ideas of this type can be traced back to Erhenfest \cite{Ehrenfest:1990p9192}. One approach is to use a \emph{coarse projective integration} method \cite{GearKevrTheo02,KevrGearHymKevrRunTheo03}.  One then extrapolates the bacterial density over a projective time step on diffusive time scales, after which the projected density needs to be \emph{lifted} to an ensemble of individual bacteria.  For such methods, besides the effect of extrapolation on the variance of the obtained results, also the effects of reconstructing the velocities and internal variables have to be systematically studied.

\section*{Acknowledgements}

The authors thank Radek Erban, Thierry Goudon, Yannis Kevrekidis and Tony Leli\`evre for interesting discussions that eventually led to this work.
This work was performed during a research stay of GS at SIMPAF (INRIA - Lille). GS warmly thanks the whole SIMPAF team for its hospitality. GS is a Postdoctoral Fellow of the
Research Foundation -- Flanders.
This work was partially supported by the Research Foundation -- Flanders
through Research Project G.0130.03 and by the Interuniversity
Attraction Poles Programme of the Belgian Science Policy Office through grant
IUAP/V/22 (GS). The scientific responsibility rests with its authors.

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv


text:
James Madison University College of Business
A business school is a university-level institution that confers degrees in business administration or management. Such a school can also be known as school of management, school of business, school of business administration, or, colloquially, b-school or biz school. A business school teaches topics such as accounting, administration, strategy, economics, entrepreneurship, finance, human resource management, information systems, logistics, marketing, organizational psychology, organizational behavior, public relations, research methods and real estate among others.
There are several forms of business schools, including school of business, business administration, and management.
Most of the university business schools are faculties, colleges, or departments within the university, and teach predominantly business courses.
In North America, a business school is often understood to be a university program that offers a graduate Master of Business Administration degrees and/or undergraduate bachelor's degrees.
This page contains text from Wikipedia, the Free Encyclopedia -	https://wn.com/Business_school
The College of Business is the business school of James Madison University in Harrisonburg, Virginia. It is a fully accredited business school that offers undergraduate degrees in accounting, computer information systems and business analytics, economics, finance and business law, international business, management, marketing, and quantitative finance. Additionally, the College of Business offers two distinctive master of business administration programs – an Innovation MBA program and an Information Security MBA program. The college also offers a master of science in accounting It has received recognitions and awards from multiple organizations, including being listed in BusinessWeek's Top 5% Undergraduate Programs and the 40th best business school in the United States. Its facilities are located in Zane Showker Hall on the southwestern part of JMU's campus.
The College of Business at JMU was first accredited in the 80s by The Association to Advance Collegiate Schools of Business. The graduate and undergraduate programs in accounting and business were reaccredited in 2012.
This page contains text from Wikipedia, the Free Encyclopedia -	https://wn.com/James_Madison_University_College_of_Business
Speaking on Business: LDS Business College
source: Building Utah
Running a Small Business in College
source: Thriftin' Ain't Easy
From College to Business
source: Loyola School of Music Industry Forum Series
Speaking on Business: Nightingale College
Speaking on Business: Westminster College
Bridging the College to Business Gap
source: Tandem Radio
Business Bros – Episode 75 – To College or Not To College?
source: Business Bros
Jason Aldrich with Robinson College of Business
source: GSU ENI
Ep. 24: The Business of College Athletics
source: Connect with KB
Simha Magal with Robinson College of Business
Balancing College, Business & Job???
source: Lens Caps
17 Business ideas for college students
source: Self Improvement Options
RaiseTheW - Samantha McGrath, Lee Honors College and Haworth College of Business Student
source: Raise The W
Business Beat: Navigating Student Debt as a College Graduate
source: Business Beat
Is College Worth It? - A Business School Grad Perspective
source: John Marty
Red Business - Episode 55 - University College Cork
source: Red Business | Cork's RedFM

redpajama_set_name: RedPajamaCommonCrawl


Chris Claremont is the man who made the X-Men what they are. He took the comic book series in the seventies and, with Dave Cockrum, John Byrne, Paul Smith, John Romita Jr, Marc Silvestri, Jim Lee, and Alan Davis among others, made those characters as we know them today. And in the process co-created Rogue, Psylocke, Shadowcat, Phoenix, The Brood, Lockheed, Shi'ar, Shi'ar Imperial Guard, Mystique, Destiny, Selene, Reverend William Stryker, Lady Mastermind, Emma Frost, Tessa, Siryn, Jubilee, Rachel Summers, Madelyne Pryor, Moira MacTaggert, Lilandra, Shadow King, Cannonball, Warpath, Mirage, Wolfsbane, Karma, Cypher, Sabretooth, Empath, Sebastian Shaw, Donald Pierce, Avalanche, Pyro, Legion, Nimrod, Gateway, Strong Guy, Proteus, Mister Sinister, Marauders, Purifiers, Captain Britain, Sunspot, Forge and Gambit amongst many. Plenty of who have made it or are making it into film.
Well, now Chris Claremont is looking for a film agent. He wrote or co-wrote the comics that much of the X-Men movies are based on, especially Days of Future Past. But he's written many other comics and novels – and now seems to like the idea of branching out into film.
If you are a film agent and would like to talk to Chris, if only to get the lay of the land, he's waiting for your message right now on Facebook.
Mondo Music Release of the Week: Big Trouble in Little China!


Sugar girl: all about sugaring bikini (+REVIEWS) | Fashion ... - Click here!
Waxing intimate areas at home wax, ... Owners of fine hairs in the bikini area, ... How to do waxing bikini area. Woman can do both conventional and deep depilation.


This NEWEST and most stylish home in beautiful Brookhaven is a true masterpiece! Located on a cul-de-sac street, this turnkey Manning model offers 5 bedrooms, 4.5 baths, a guest suite on main, a bright open floor plan and 3 levels of usable space. This is a GREAT house for entertaining and everyday living. The gorgeous and well-equipped chef's kitchen opens to the light-filled family room with coffered ceilings. Let the parties spill out onto the deck and below in your own private backyard oasis featuring a salt water, gunite pool. Retire to your luxe master bedroom retreat with spa-like bath and enormous custom walk-in closet. 2nd & 3rd floor bonus rooms! This gem offers tons of storage, beautiful built-ins throughout, 5 inch ebony plank flooring downstairs, plantation shutters, updated lighting and more. If you're looking for an amenity filled neighborhood, you've hit the jackpot. Drool-worthy house + resort style amenities + top rated schools + low Union County taxes = WOW!


2014 v. 2015- who wore Golden Globes fashion best? Find out here.
Pastels for Spring? Groundbreaking. No, our favourite leading ladies steered well clear of subtler hues and instead opted for punchy head to toe brights to get the cameras flashing- and we loved it. Major props to Girls' Allison Williams in Armani Prive, Jenna Dewan Tatum in Emilio Pucci and Amy Poehler in Stella McCartney.
The only way to wear trousers on the red carpet? Like a boss.
Just when we thought we couldn't love Emma Stone any more, she goes and pulls this Lanvin lovely out of the bag, perfectly accessorised with killer points and just rolled out of bed hair. Can trousers ever beat a killer dress on the red carpet? With the backing of Miss Stone, along with Lorde in Narciso Rodriguez, we'd say the answer's a resounding yes.
Need a show stopping look? No probs. When the chips are down and you can't decide which look to go for, always go big. Big silhouette, big thigh high split, big embellishment… big headline. Just ask J-Lo and Jemima Kirke in all their shimmering caped awesomeness.


text:
Tech/Modernization
Email: admin@navalnews.net
Naval News
Technology/Modernization
Russia to build five more Project 23130 logistic support tankers
Russian Navy's "Gremyashchiy" Class corvette completes its sea trials
India successfully test-fires K-4 Ballistic Missile
Royal Malaysian Navy's first littoral mission ship (LMS) arrives Malaysia
Theodore Roosevelt Carrier Strike Group departs for Indo-Pacific deployment
Philippines Navy sends a naval force to Mid-East to evacuate OFWs
U.S. Navy seizes suspected Iranian missile parts set for Yemen
CTF150 (FS NIVOSE) Achieves Major Drug Bust in the Arabian Sea
US Coast Guard Seized Submarine While Smuggling Tons of Drugs
Egypt conducts "QADER 2020" exercise with the participation of combined units
China-Pakistan Sea Guardians 2020 joint exercise concludes fleet drills
Russia and Japan navies to hold drills in the Gulf of…
China, Pakistan kick off "Sea Guardians-2020" naval exercise in Karachi
Iran, Russia and China completed trilateral exercise
Abraham Lincoln Carrier Strike Group Arrives in San Diego after 295-day…
BrahMos equipped Sukhoi-30 MKI to keep eye on Indian Ocean
US Navy announces aircraft carrier CVN81 will be named USS DORIS…
"NAVIRIS" JV consists of Fincantieri and Naval Group is fully operational
HMAS Toowoomba leaves home for sixth Middle East deployment
Home Miscellaneous Wreck of the German sunken cruiser found off the Falkland Islands
Wreck of the German sunken cruiser found off the Falkland Islands
Dorian Archus
We have been witnessing succesful researches of wreckages of the ships sank during battles at previous century. These are very precious efforts because those wreckages help us to discover historical details. Last week researchers shared new photos and videos of a sunken warship off coast of the Falkland Islands, which could be considered as the latest link of the chain of researches.
The wreck of SMS Scharnhorst, a sunken armored battlecruiser of the German Navy, has been located off the Falkland Islands, the Falklands Maritime Heritage Trust said.
The wreck, lying some 1,600 metres (5,280ft) down, was discovered by marine archaeologist Mensun Bound, the first success in a search by the Falklands Maritime Heritage Trust to locate the warship.
It set out to locate all four German ships – heavy cruisers Scharnhorst and Gneisenau, plus light cruisers Nürnberg and Leipzig – initially to mark the centenary of the sinkings in 2014; the battle is a seminal moment in Falklands history, second only to the 1982 conflict, with December 8 commemorated each year as 'Battle Day'.
The search five years ago drew a blank, but just three days into a renewed attempt in 2019, four deep-sea robot submarines and high-tech side-scan sonar on their mother ship, the search vessel Seabed Constructor, located the Scharnhorst as an area of seabed the size of Somerset was methodically scoured.
"The moment of discovery was extraordinary," said Mr Bound. "We are often chasing shadows on the seabed, but when the Scharnhorst first appeared in the data flow, there was no doubt that this was one of the German fleet. You could even see the impact crater. We sent down a remotely-operated vehicle to explore and almost straight away we were into a debris field that said 'battle'. Suddenly she just came out of the gloom with great guns poking in every direction."
The Battle of the Falkland Islands followed the Battle of Coronel, fought off the coast of Chile in November 1914, where Graf von Spee's fleet overpowered the Royal Navy and in which 1,600 British sailors perished. A month later a British squadron under the command of Vice-Admiral Doveton Sturdee defeated Graf von Spee's squadron, comprising the Scharnhorst, Gneisenau, Nürnberg and Leipzig. The action was particularly important because as a consequence of the battle, the German East Asia Squadron, Germany's only permanent overseas naval formation, effectively ceased to exist, bringing an end to commerce raiding by German warships.
The Scharnhorst, built in Hamburg in 1905, was the first to be sunk, having sustained substantial damage inflicted by HMS Invincible and HMS Inflexible. Tragically, 2,200 German sailors died, including Graf von Spee himself and his two sons – Heinrich aboard the Gneisenau, and Otto aboard the Nürnberg.
"It is less than a month since Remembrance Day, when we commemorated the millions who died in the First World War and subsequent conflicts. One episode in that conflict was the Battle of the Falklands in 1914," Donald Lamont, Chairman of The Falklands Maritime Heritage Trust, noted.
The search we organised had as its aim the locating of all ships of the German squadron, so that we may learn more about the Battle and commemorate all who perished in it. The site of the wrecks can now be protected," Lamont added.
The Falkland Maritime Heritage Trust is now seeking to have the site formerly protected in law. The wreck was not touched or in any way disturbed during the operation.
Article bottom
Falklands Maritime Heritage Trust
German Cruiser
SMS Scharnhorst
World War-I
Calender of Archives (Click to See Daily News)
Which is the most powerful frigate in the world?
Admiral Grigorovich class (Russia)
Type 054A class (China)
Shivalik class (India)
Carlo Bergamini class (Italy)
Aquitaine class (France)
Alvaro de Bazan class (Spain)
Iver Huitfeldt class (Denmark)
Sachsen class (Germany)
Admiral Gorshkov class (Russia)
Follow us on Instagram @naval_news
"Naval News" is an online defense media company, focused on Warships, Submarines, Naval Air, Marines, Coast Guard and Naval Missiles.
Contact us: admin@navalnews.net
© Copyright 2016 - 2019 The Coral Media/Naval News
All new Russian ships to have capability to carry Zircon Missiles
How South Korea will respond the North Korean missile attack
... https://navalnews.net/wp-content/uploads/2020/01/naval.ico https://navalnews.net/wp-content/uploads/2020/01/88NN.png ...