text:
\section{Introduction}

Diffusion processes have won an increasing interest in the past years in the statistical physics community, to model physical phenomena and to sample the underlying probability measure characterizing the state of the system~\cite{Balian07}. The average value of a thermodynamic function~$R$ (as the energy, the pressure, a length, a flux, ...) under this probability distribution is given by an integral over the very high-dimensional configurational space. An important motivation for this work is the averaging of mean properties for systems subject to an external driving. In this case the invariant probability measure is often not known explicitly. The goal can be to compute a transport coefficient, a free energy or more generally the response to a non-equilibrium forcing. From a practical point of view, the unknown probability measure is sampled by integrating a stochastic dynamics~\cite{Allen87,Frenkel02,Tuckerman10,Leimkuhler16}
\begin{equation}
\label{eq:process intro}
	\dd X_t = b(X_t) \, \dd t + \sigma(X_t) \, \dd W_t.
\end{equation}
Two prototypical dynamics in molecular simulation are the Langevin and overdamped Langevin dynamics~\cite{Allen87,Leimkuhler16_book}. At equilibrium, the Langevin dynamics evolves positions $q$ and momenta $p$ as
\begin{equation}
  \label{eq:Langevin_introduction}
\left\{
\begin{aligned}
  \dd q_t &= \frac {p_t} m \, \dd t, \\
  \dd p_t &= -\nabla V(q_t) \, \dd t - \frac \gamma m p_t \, \dd t + \sqrt{2 \gamma \beta^{-1}} \, \dd W_t,
\end{aligned}
\right.
\end{equation}
where $\gamma>0$ is the friction coefficient, $m>0$ is the mass of a particle and $\beta > 0$ is proportional to the inverse temperature. The potential energy function is denoted by $V$ and $W_t$ is a multi-dimensional standard Brownian motion. In the limit of large frictions $\gamma$, this equation becomes after proper rescaling the overdamped Langevin dynamics~\cite{Freidlin04}:
\begin{equation}
\label{eq:cv overdamped}
\dd q_t = -\nabla V(q_t) \, \dd t + \sqrt{2 \beta^{-1}} \, \dd W_t.
\end{equation}
Nonequilibrium versions of the above dynamics are obtained for instance by considering non-gradient forces rather than $-\nabla V$.

For any ergodic dynamics~\eqref{eq:process intro}, macroscopic properties are computed via averages over a trajectory as
\begin{equation*}
  \mathbb{E}[R] := \lim_{T \to \infty} \widehat R_T \quad \mbox{a.s.}, \qquad \widehat R_T = \frac 1 T \int_0^T R(X_t) \, \dd t.
\end{equation*}
The statistical error for these estimators is characterized by the asymptotic variance:
\begin{equation}
\label{eq:def variance}
	\sigma_R^2 = \lim_{T \to \infty} T \, \mathrm{Var}\left[\widehat R_T \right].
\end{equation}
In many cases of interest ergodic means converge very slowly, requiring the use of variance reduction techniques to speed up the computation. Two types of phenomena can lead to large statistical errors: first, the metastability arising from multimodal potentials, which can greatly increase the correlation of the trajectory in time and lead to large variances; second, a high signal-to-noise ratio, which is typical when averaging small linear responses as for the computation of transport coefficients~\cite{Evans13_book, Tuckerman10}.

When the system is at equilibrium, by which we mean that detailed balance holds, the invariant probability measure is often known and it is possible to use standard variance reduction techniques~\cite{Rubinstein81, Fishman96, Caflisch98, Liu01, Lapeyre03} such as importance sampling~\cite{Chen97, Liu01} or stratification~\cite{Geyer94, Meng96, Kong03, Tan04, Shirts08}. This allows to address both metastability issues and high noise-to-signal ratios.

For non-equilibrium systems, and more generally when the invariant probability measure is not known, reducing the variance is challenging since standard variance reduction methods cannot be used. Note that reducing the metastability would require modifying the dynamics while keeping the invariant probability measure unchanged, or at least knowing how it changes (see for instance~\cite[Section~3.4]{Lelievre16}). When the latter is not known, this task is hard to perform. On the other hand increasing the signal-to-noise ratio is feasible even for non-equilibrium dynamics. This is the goal of the present work, where we rely on control variates. Control variates laying on the concept of "zero-variance" principle have been already used in molecular simulation~\cite{Assaraf99}. This approach was also studied in Bayesian inference simulations~\cite{Henderson97, Mira03, Dellaportas12, Mira13, Oates17}, where configurations are sampled with Markov chains rather than diffusion processes. This type of techniques has however been restricted to cases where the invariant probability measure is known, except for specific settings such as~\cite{Goodman09}, where a coupling strategy is described.

In the present work, which relies on ideas announced in~\cite[Section 3.4.2]{Lelievre16}, control variates are constructed without any knowledge of the expression of the invariant probability measure. We build an unbiased modified observable $R+\xi=R+\mathcal{L} \Phi$, where $\mathcal{L}$ is the generator of the dynamics, which is of smaller variance (at least in some asymptotic regime):
\begin{equation*}
	\mathbb{E}[R + \xi] = \mathbb{E}[R] \quad \mbox{and} \quad \sigma_{R + \xi}^2 < \sigma_R^2.
\end{equation*}
The optimal choice for the control variate $\xi$ is $\xi = \mathcal{L} \Phi$ where $\Phi$ is the solution of the following Poisson equation:
\begin{equation*}
	-\mathcal{L} \Phi = R - \mathbb{E}[R].
\end{equation*}
The general strategy we consider consists in approximating this partial differential equation (PDE) by a simplified one, with an operator $\calL_0$, for which the solution ${\Phi_0}$ can be analytically computed or numerically approximated with a good precision. Theorems~\ref{th:var order 2 eta exact} and~\ref{th:var order 2 eta} provide an analysis of the asymptotic variance $\sigma_{R + \xi}^2$.

We present numerical results illustrating the general method in three practical cases. In particular we provide in each case a simplified process, associated to a simplified Poisson problem. In these applications we are interested in averaging the linear response of an observable with respect to a non-equilibrium perturbation. This is a challenging class of problems since the average quantity is small and thus the relative statistical error is large. We present the problems we consider by increasing complexity of the setup. We start with the computation of the mobility of a particle in a periodic two-dimensional potential. The control variate can be approached with a very high precision by a numerical method based on a spectral basis, allowing to illustrate Theorem~\ref{th:var order 2 eta}. We next estimate the conductivity of an atom chain~\cite{Bonetto00, Lepri16, Dhar08}. The number of state variables is much larger (up to several hundreds of degrees of freedom in our simulations) but the geometrical setting is one-dimensional. The control variate can be computed analytically when taking a harmonic model as a reference, and thus it does not require any additional numerical procedure. The third application is a dimer in a solvent, whose mean length is estimated under an external shearing force. In the latter case the difficulty comes from the fact that the system is high-dimensional and not as structured as the atom chain.

This article is organized as follows. We present in Section~\ref{s:strategy} the general strategy for building control variates and state a result making precise how control variates behave in a perturbative framework. We then turn to the case of a single particle in a one-dimensional periodic potential under a non-gradient forcing in Section~\ref{s:sinus potential}; the computation of the thermal flux passing through a chain in Section~\ref{s:quasi harmonic}; and the estimation of the mean length of a dimer in a solvent under an external shearing stress in Section~\ref{s:dimer}. Some technical results are gathered in the appendices.

\section{General strategy}
\label{s:strategy}

The definition of the asymptotic variance of time averages along a trajectory requires to introduce more precisely the generator of the process and some associated functional spaces, which is done in Section~\ref{ss:asyvar}. The concept of control variate is then explained Section~\ref{ss:ideal cv}, as well as the so-called "zero-variance principle". We give in Section~\ref{ss:perturbative} our perturbative construction of control variate in an abstract setting and state the main theorem quantifying the variance reduction in a limiting regime. Finally Section~\ref{ss:numerics poisson} provides a generalized version of this theorem in the case when an approximate solver is used.

\subsection{Asymptotic variance}
\label{ss:asyvar}

\newcommand{X_t^\eta}{X_t^\eta}

The state space $\mathcal{X}$ is typically the full space $\mathbb{R}^d$ or a bounded domain with periodic boundary conditions $\mathbb{T}^d$. For some dynamics such as Langevin dynamics, auxiliary variables with values in $\mathbb{R}^d$ are added, so that in this case $\mathcal{X} = \mathbb{R}^d \times \mathbb{R}^d$ or $\mathcal{X} = \mathbb{T}^d \times \mathbb{R}^d$. As suggested in the introduction, we decompose the generator of the process~\eqref{eq:process intro} as a sum $\mathcal{L} = \calL_0 + \widetilde{\mathcal{L}}$ of a reference generator $\calL_0$ and a perturbation $\widetilde{\mathcal{L}}$. In order to study the asymptotic regime corresponding to small perturbations, we use a parameter $\eta \in \mathbb{R}$ to interpolate smoothly between $\calL_0$ and $\mathcal{L}$, and define
\begin{equation*}
	\calL_\eta = \calL_0 + \eta \widetilde{\mathcal{L}}.
\end{equation*}
We suppose that these operators $\calL_\eta$ write:
\begin{equation*}
	\calL_\eta = b_\eta\cdot \nabla + \frac 1 2 \sigma_\eta \sigma_\eta^\top : \nabla^2 = \sum_{i=1}^d \left( b_\eta \right)_i \partial_{x_i} + \frac 1 2 \sum_{i,j=1}^d \left(\sigma_\eta \sigma_\eta^\top \right)_{i,j} \partial_{x_i x_j},
\end{equation*}
with $b_\eta$ and $\sigma_\eta$ are smooth. They are then the generators of the following stochastic processes on $\mathcal{X}$, indexed by $\eta$:
\begin{equation}
\label{eq:process}
	\dd X_t^\eta = b_\eta(X_t^\eta) \, \dd t + \sigma_\eta(X_t^\eta) \, \dd W_t,
\end{equation}
where $W_t$ is a $d$-dimensional standard Brownian motion. Let us assume that $b_\eta$ and $\sigma_\eta$ are such that the following holds.
\begin{assumption}
\label{as:pieta}
The dynamics~\eqref{eq:process} admits a unique invariant probability measure $\pi_\eta$ for any $\eta \in \mathbb{R}$. Moreover, trajectorial ergodicity holds: for any observable $R \in \mathrm{L}^1(\pi_\eta)$,
\begin{equation*}
  \mathbb{E}_\eta[R] := \int_\mathcal{X} R(x) \, \dd \pi_\eta(x) = \lim_{T \to \infty} \frac 1 T \int_0^T R(X_t^\eta) \, \dd t \quad \mathrm{a.s.}
\end{equation*}
\end{assumption}
Sufficient conditions for this to hold are discussed after Assumption~\ref{as:lyapunov in L2}. Let us now make the functional spaces precise. We denote by $(\mathcal{K}_n)_{n \in \mathbb{N}}$ a family of so-called Lyapunov functions with values in $[1, +\infty)$. The associated weighted $\mathrm{L}^\infty$ spaces are:
\begin{equation*}
\begin{aligned}
  \forall n \in \mathbb{N},\quad {\rmL_n^\infty} &= \left\{ \varphi \ \mbox{measurable} \ \left| \ \|\varphi\|_{\rmL_n^\infty} < \infty \right. \right\} , \qquad \|\varphi\|_{\mathrm{L}_n^\infty} = \left\| \frac \varphi {\mathcal{K}_n} \right\|_{\rmL^\infty}.
\end{aligned}
\end{equation*}
We make the following assumption on the Lyapunov functions.
\begin{assumption}
\label{as:lyapunov in L2}
For any $\eta \in \mathbb{R}$, the function $\mathcal{K}_n$ belongs to~${\rmL^2(\pieta)}$. In particular,
\begin{equation*}
	\forall n \in \mathbb{N},\, \forall \eta \in \mathbb{R}, \quad {\rmL_n^\infty} \subset {\rmL^2(\pieta)}.
\end{equation*}
We also assume that for any $n,n' \in \mathbb{N}$, there exists $m \in \mathbb{N}$ such that $\mathcal{K}_n \mathcal{K}_{n'} \in {\rmL_m^\infty}$.
\end{assumption}
The first part of Assumption~\ref{as:lyapunov in L2} is typically obtained by using a family of Lyapunov functions satisfying conditions of the form
\begin{equation}
\label{eq:lyapunov condition}
	\calL_\eta \mathcal{K}_n \leqslant -\alpha_{n,\eta} \mathcal{K}_n + b_{n,\eta},
\end{equation}
for some $\alpha_{n,\eta}>0$ and $b_{n,\eta} \in \mathbb{R}$. Indeed, after integration against $\pi_\eta$,
\begin{equation*}
  \begin{aligned}
    0 = \int_\mathcal{X} \calL_\eta \mathcal{K}_n \dd \pi_\eta \leqslant -\alpha_{n,\eta} \int_\mathcal{X} \mathcal{K}_n \dd \pi_\eta + b_{n,\eta},
\end{aligned}
\end{equation*}
so that
\begin{equation*}
  1 \leqslant \int_\mathcal{X} \mathcal{K}_n \dd \pi_\eta \leqslant \frac {b_{n,\eta}}{\alpha_{n,\eta}}.
\end{equation*}
We can then conclude with the second part of Assumption~\ref{as:lyapunov in L2} since, for any $n \geqslant 1$, there exist $C_n > 0$ and $m \geqslant 1$ such that $1 \leqslant \mathcal{K}_n^2 \leqslant C_n \mathcal{K}_m$. The condition~\eqref{eq:lyapunov condition} also implies the existence of an invariant probability measure for any $\eta \in \mathbb{R}$ when a minorization condition holds~\cite{Hairer11_chap}. A typical choice for the Lyapunov functions are the polynomials $\mathcal{K}_n(x) = 1+|x|^n$. This choice satisfies the second part of Assumption~\ref{as:lyapunov in L2} when the invariant probability measure has moments of any order. Unless otherwise mentioned, we always consider this choice in the sequel (which is standard for Langevin and overdamped Langevin dynamics, see~\cite{Talay02,Mattingly02,Kopec14,Kopec15}). As for Assumption~\ref{as:pieta}, trajectorial ergodicity holds when the generator $\calL_\eta$ is elliptic or hypoelliptic, and there exists an invariant probability measure with positive density with respect to the Lebesgue measure~\cite{Kliemann87}. The latter condition follows if the measure which appears in the minorization condition has a positive density with respect to the Lebesgue measure.

We denote for any function $\varphi \in \mathrm{L}^1(\pi_\eta)$ the projection on the space of mean zero functions by:
\begin{equation*}
  \Pi_\eta \varphi = \varphi - \mathbb{E}_\eta[\varphi].
\end{equation*}
For any operator $A \in \mathcal{B}(E)$ (bounded on the Banach space $E$), the operator norm is defined as
\begin{equation*}
	\| A \|_{\mathcal{B}(E)} = \sup_{\| \varphi \|_E = 1} \| A \varphi \|_E.
\end{equation*}
Let us now make the following assumption.
\begin{assumption}
\label{as:generator eta}
For any $n \in \mathbb{N}$, the $\mathrm{L}^2(\pi_\eta)$ norms of the Lyapunov functions are uniformly bounded on compact sets of $\eta$: for any $\eta_*>0$ there exists a constant $C_{n,\eta_*}$ such that
\begin{equation}
\label{eq:uniform lyapunov measure}
	\forall |\eta| \leqslant \eta_*, \quad \| \mathcal{K}_n \|_{\rmL^2(\pieta)} \leqslant C_{n,\eta_*}.
\end{equation}
Moreover $\calL_\eta$ is invertible on $\Pi_\eta {\rmL_n^\infty}$. Finally the inverse generator is bounded uniformly on compact sets of $\eta$:
\begin{equation}
\label{eq:uniform inverse generator}
	\forall |\eta| \leqslant \eta_*, \quad \left\| -{\calL_\eta^{-1}} \right\|_{\mathcal{B} \left(\Pi_\eta {\rmL_n^\infty} \right)} \leqslant C_{n,\eta_*}.
\end{equation}
\end{assumption}
The invertibility of $\calL_\eta$ on $\Pi_\eta {\rmL_n^\infty}$ is a standard result which follows typically from the Lyapunov conditions~\eqref{eq:lyapunov condition} and a so-called minorization condition~\cite{Hairer11_chap}. It has been proved for a large variety of problems~\cite{Eckmann99,Talay02,Mattingly02,Kopec14,Kopec15}. Conditions~\eqref{eq:uniform lyapunov measure} and~\eqref{eq:uniform inverse generator} are needed to prove Theorems~\ref{th:var order 2 eta exact} and~\ref{th:var order 2 eta} to come. Condition~\ref{eq:uniform lyapunov measure} can be obtained by showing uniform bounds on the coefficients which appear in the Lyapunov conditions, while condition~\eqref{eq:uniform inverse generator} additionally requires some uniformity on the minorization condition.

When Assumptions~1 to 3 hold, the asymptotic variance introduced in~\eqref{eq:def variance} is finite for any $\varphi \in {\rmL_n^\infty}$ and the following formula holds~\cite{Lelievre16}:
\begin{equation}
  \label{eq:variance eta}
  \sigma_{\varphi, \eta}^2 = 2 \left\langle \varphi, -{\calL_\eta^{-1}} \Pi_\eta \varphi \right\rangle_{\eta},
\end{equation}
where $\left\langle \cdot, \cdot \right\rangle_{\eta}$ denotes the canonical scalar product on~${\rmL^2(\pieta)}$. We refer to Appendix~\ref{ap:variance} for more details on the numerical estimation of the asymptotic variance.
\begin{remark}
We choose to work directly with weighted $\mathrm{L}^\infty$ spaces as this is the relevant setting for Theorems~\ref{th:var order 2 eta exact} and~\ref{th:var order 2 eta}. Note that the asymptotic variance of an observable $\varphi \in {\rmL^2(\pieta)}$ can also be defined using perturbative arguments relatively to an equilibrium reference dynamics~\cite{Dolbeault15, Iacobucci17}. Contrarily to the ${\rmL_n^\infty}$ framework one would however be restricted in this case to small non-equilibrium perturbations.
\end{remark}

\subsection{Ideal control variate}
\label{ss:ideal cv}

We recall in this section what a control variate is in our context and show how the construction of an optimal control variate can be reformulated as solving a Poisson problem. The functional framework is made precise in a second step using Assumption~\ref{as:generator eta}. We say that a function $\xi$ is a control variate of the observable $R$ for the process $X_t^\eta$ with generator $\calL_\eta$ if
\begin{equation*}
	\mathbb{E}_\eta[R + \xi] = \mathbb{E}_\eta[R] \quad \mbox{and} \quad \sigma_{R + \xi, \eta}^2 < \sigma_{R, \eta}^2.
\end{equation*}
The principle of our method, already explained in~\cite{Lelievre16}, is based on the equation which characterizes the invariance of the measure $\pi_\eta$: for any function $\Phi$,
\begin{equation*}
	\mathbb{E}_\eta[\calL_\eta \Phi] = 0.
\end{equation*}
This shows that control variates $\xi$ of the form $\xi = \calL_\eta \Phi$ automatically ensure that $\mathbb{E}_\eta[R + \xi] = \mathbb{E}_\eta[R]$, whatever the choice of $\Phi$. In order for $\xi$ to be a good control variate, the modified observable $R + \xi = R + \calL_\eta \Phi$ should however be of small asymptotic variance. The optimal choice, denoted by $\Phi_\eta$ and referred to as the "zero-variance principle"~\cite{Assaraf99, Mira03}, is to make the modified observable constant. This constant is then necessarily equal to $\mathbb{E}_\eta[R]$, and $\Phi_\eta$ is the solution of the Poisson problem:
\begin{equation}
\label{eq:cv poisson problem}
	-\calL_\eta \Phi_\eta = R - \mathbb{E}_\eta[R].
\end{equation}
Assuming that $R \in {\rmL_n^\infty}$ for some $n\in \mathbb{N}$, the problem~\eqref{eq:cv poisson problem} admits a unique solution $\Phi_\eta \in {\rmL_n^\infty}$ when Assumption~\ref{as:generator eta} holds.

In practice two problems arise when trying to solve~\eqref{eq:cv poisson problem}. First, the equation~\eqref{eq:cv poisson problem} is a very high-dimensional PDE for most purposes and the complexity of such problems scales exponentially with the dimension. Second, $\mathbb{E}_\eta[R]$ is not known since it is precisely the quantity we are trying to compute. We discuss in the next section how to approximate the solution of~\eqref{eq:cv poisson problem}, at least for small $\eta$.

\subsection{Perturbative control variate}
\label{ss:perturbative}

The key assumption in our approach is to assume that we can compute the solution $\Phi_0$ of the reference Poisson problem corresponding to $\eta=0$:
\begin{equation}
\label{eq:poisson problem simple}
-\calL_0 \Phi_0 = R - \mathbb{E}_0[R].
\end{equation}
In practice $\calL_0$ is the generator of a simplified dynamics (depending on the problem), and $\calL_\eta$ is the generator of the problem at hand (say for $\eta=1$). Let us emphasize that the dynamics associted with $\mathcal{L}_0$ need not be an equilibrium dynamics (see the example discussed in Section~\ref{s:quasi harmonic}). The small parameter $\eta$ is used to quantify the discrepancy between the optimal function $\Phi_\eta$ and its approximation ${\Phi_0}$ in a perturbative framework. We refer to Section~\ref{ss:numerics poisson} for a discussion on the numerical resolution of~\eqref{eq:poisson problem simple}, and to the end of this section for further comments on the decomposition.

We define the so-called core space $\mathcal{S}$ as the set of all $\mathcal{C}^\infty$ functions which grow at infinity at most like $K_n$ for some $n$, and whose derivatives also grow at most like $K_n$ for some $n$. Such a space was considered in~\cite{Talay02} for instance. More precisely,
\begin{equation}
\label{eq:core}
\begin{aligned}

	\mathcal{S} &= \big\{ \varphi \in \mathcal{C}^\infty(\mathcal{X}) \big|\ \forall k \in \mathbb{N},\, \exists n \in \mathbb{N}, \quad\partial_k \varphi \in {\rmL_n^\infty} \big\}.
\end{aligned}
\end{equation}
The space $\mathcal{S}$ is dense in ${\rmL^2(\pieta)}$ under Assumption~\ref{as:lyapunov in L2}, since $\mathcal{C}^\infty$ functions with compact support are included in $\mathcal{S}$. We need an additional assumption in our analysis to ensure that ${\Phi_0} \in \mathcal{S}$.
\begin{assumption}
\label{as:generator simple}
The space $\mathcal{S}$ is stable by the generator $\calL_0$ and $\calL_0$ is invertible on the space $\Pi_0 \mathcal{S}$ composed of functions with average $0$ with respect to the invariant probability measure~$\pi_0$. This means that, for any $\varphi \in \Pi_0 \mathcal{S}$, there exists a unique solution $\psi \in \Pi_0 \mathcal{S}$ to the Poisson equation
\begin{equation*}
  -\calL_0 \psi = \varphi.
\end{equation*}
\end{assumption}
Assumption~\ref{as:generator simple} can be proved to hold for Langevin and overdamped Langevin dynamics at equilibrium under certain assumptions on the potential $V$, see~\cite{Talay02, Kopec14, Kopec15}. The generator of the perturbation should also satisfy the following condition.
\begin{assumption}
\label{as:Lrep}
The generator $\widetilde{\calL}$ of the perturbation is such that $\mathcal{S}$ is stable by $\widetilde{\calL}$ and $\widetilde{\mathcal{L}}^* \mathbf{1} \in {\rmL^2(\pinot)}$.
\end{assumption}
Here and in the following we denote by $B^*$ the adjoint of a closed operator $B$ on the functional space ${\rmL^2(\pinot)}$. Assumption~\ref{as:Lrep} is easy to check, since $\widetilde{\mathcal{L}}$ is typically a differential operator with coefficients in $\mathcal{S}$.

Let us define the following modified observable involving ${\Phi_0} \in \Pi_0 \mathcal{S}$, defined in~\eqref{eq:poisson problem simple}:
\begin{equation*}
\begin{aligned}
  \phi_\eta &:= R + \calL_\eta \Phi_0. \\

\end{aligned}
\end{equation*}
The following theorem makes precise the main properties of this modified observable.
\begin{theorem}
\label{th:var order 2 eta exact}
Fix $R \in \mathcal{S}$. Under Assumptions~1 to 5, $\phi_\eta \in \mathcal{S}$ is well defined for any $\eta \in \mathbb{R}$, and $\mathbb{E}_\eta[\phi_\eta] = \mathbb{E}_\eta[R]$. Moreover, for any $\eta_*>0$, there exists $C_{R, \eta_*} > 0$ such that, for any $| \eta | \leqslant \eta_*$, the asymptotic variance satisfies
\begin{equation}
\label{eq:var order 2 eta exact}
	\sigma_{\phi_\eta, \eta}^2 = 2 \eta^2 \left\langle A R, -\calL_0^{-1} A R \right\rangle_0 + \eta^3 E_{R, \eta},
\end{equation}
with $A= - \widetilde{\mathcal{L}} \calL_0^{-1} \Pi_0$ and $| E_{R, \eta} | \leqslant C_{R, \eta_*}$.
\end{theorem}

The scalar products involved in the previous theorem are well defined since $\mathcal{S}$ is stable by $A$, and $\mathcal{S} \subset {\rmL^2(\pieta)}$ for any $\eta \in \mathbb{R}$. Equation~\eqref{eq:var order 2 eta exact} shows that the standard error $\sqrt{\frac{\sigma_{\phi_\eta, \eta}^2} T}$ committed on the empirical estimator $\widehat \varphi_T$ of $\mathbb{E}_\eta[R]$ after a time $T$ is of leading order $\eta/\sqrt{T}$.

The scaling $\eta^2$ of the asymptotic variance formally comes from the fact that the modified observable writes $\phi_\eta = \mathbb{E}_\eta[R] + \mathrm{O}(\eta)$. Indeed,
\begin{equation}
  \label{eq:Leta Linv}
  \calL_\eta \calL_0^{-1} \Pi_0 = \Pi_\eta (\calL_0 + \eta \widetilde{\mathcal{L}}) \calL_0^{-1} \Pi_0 = \Pi_\eta + \eta \Pi_\eta \widetilde{\mathcal{L}} \calL_0^{-1} \Pi_0 = \Pi_\eta (1-\eta A),
\end{equation}
so that the modified observable can be rewritten as:
\begin{equation}
\begin{aligned}
\label{eq:modified obs}
\phi_{\eta} &= R + \calL_\eta \Phi_0 = R - \calL_\eta \calL_0^{-1} \Pi_0 R = \mathbb{E}_\eta[R] + \eta \Pi_\eta A R.
\end{aligned}
\end{equation}
In particular,
\begin{equation}
\label{eq:projected modified obs}
\Pi_\eta \phi_\eta = \eta \Pi_\eta A R.
\end{equation}
The remainder of the proof consists in carefully estimating remainders in some truncated series expansion of $-\mathcal{L}_\eta^{-1}\Pi_\eta$; see Appendix~\ref{ap:proof Th2}.

\begin{remark}
  The formula~\eqref{eq:var order 2 eta exact} can in fact be replaced by an expansion in powers of $\eta$ with a truncation at an arbitrarily high order and a remainder controlled uniformly in $|\eta| \leqslant \eta_*$. This can be proved by an immediate generalization of the proof we provide in Appendix~\ref{ap:proof Th2}.
\end{remark}

In the following applicative sections we cannot always prove that Assumptions~\ref{as:generator eta} and~\ref{as:generator simple} hold true, but the scaling of the variance predicted by Theorem~\ref{th:var order 2 eta exact} is nevertheless numerically observed to hold. More importantly, let us emphasize that the modeling process is crucial to write the generator as the sum of a reference generator~$\mathcal{L}_0$ for which we can solve Poisson equations, and an additional term. Such a decomposition can always be performed, but the quality of the resulting control variate will strongly depend on the choice of the decomposition. Typically, the control variate method works well when the additional term is a perturbation of the reference generator (which corresponds to the regime where Theorem~\ref{th:var order 2 eta exact} can be applied). In practice, this can be checked a posteriori, once the simulation has been carried out, by comparing the asymptotic variance of the two estimators, one with and one without the control variate. Of course, one should keep the one with the smallest variance.

\subsection{Numerical resolution of the reference Poisson problem}
\label{ss:numerics poisson}

We discuss in this section a strategy to compute the solution to~\eqref{eq:poisson problem simple} when this equation cannot be analytically solved. We rely for this on a Galerkin strategy, and look for an approximation of the solution ${\Phi_0}$ to the Poisson problem~\eqref{eq:poisson problem simple} in a subspace $V_M \subset {\rmL^2(\pinot)}$ of finite dimension~$M$. For simplicity we suppose that $V_M \subset \Pi_0 {\rmL^2(\pinot)}$ (which corresponds to a conformal approximation). This implies in particular that $\Phi_M \in \Pi_0 {\rmL^2(\pinot)}$ has mean zero with respect to $\pi_0$. We also assume that $V_M \subset \mathrm{H}^2(\pi_0)$ to avoid regularity issues. The optimal choice for the approximation $\Phi_{0,M}$ is to minimize the variance of the modified observable $R+\calL_0 \Phi_{0,M}$ for the reference dynamics:
\begin{equation}
  \label{eq:minimization_variance_Galerkin}
  \min_{\varphi \in V_M} \sigma_{R+\calL_0 \varphi,0}^2 = \sigma_{R,0}^2 + \min_{\varphi \in V_M} \frac 1 2 \left\langle (\calL_0 + \calL_0^*) \varphi, 2 \calL_0^{-1} \Pi_0 R + \varphi \right\rangle_0.
\end{equation}
The latter equality follows from Lemma~\ref{lemma:variance LU} in Appendix~\ref{ap:proof Th2}, in the particular case when $\eta=0$. In the case when $\calL_0$ is not the generator of a stochastic differential equation the quantity $\sigma_{R+\calL_0 \varphi,0}^2$ cannot be interpreted as a variance but the analysis we provide here remains valid. The necessary optimality condition for a minimizer $\Phi_{0,M}$ of~\eqref{eq:minimization_variance_Galerkin} is given by the following Euler-Lagrange equation:
\begin{equation*}
	\forall \psi \in V_M, \quad \left\langle (\calL_0 + \calL_0^*) (\Phi_{0,M} + \calL_0^{-1} \Pi_0 R), \psi \right\rangle_0 = 0.
\end{equation*}
Introducing the orthogonal projector $\Pi_M$ on $V_M$ (with respect to the scalar product on ${\rmL^2(\pinot)}$), the latter equation can be rewritten as
\begin{equation}
\label{eq:euler lagrange numerics}
\Pi_M (\calL_0 + \calL_0^*) (\Phi_{0,M} + \calL_0^{-1} \Pi_0 R) = 0.
\end{equation}
In practice we distinguish two cases.
\begin{itemize}
\item[(i)] For reversible dynamics such as the overdamped Langevin dynamics~\eqref{eq:cv overdamped} $\calL_0 = \calL_0^*$, so the equation reduces to
\begin{equation}
\label{eq:minimizer reversible}
	-\Pi_M \calL_0 \Phi_{0,M} = \Pi_M R.
\end{equation}
\item[(ii)] For Langevin dynamics at equilibrium (see~\eqref{eq:Langevin_introduction}), we consider a tensorized basis involving Hermite elements as in Section~\ref{s:sinus potential} or in~\cite{Roussel17}. The symmetric part of the generator:
\begin{equation*}
	\frac 1 2 (\calL_0 + \calL_0^*) = -\gamma \beta^{-1} \nabla_p^* \nabla_p,
\end{equation*}
diagonalizes the Hermite polynomials so we have the commutation rule $\Pi_M (\calL_0 + \calL_0^*) = (\calL_0 + \calL_0^*) \Pi_M$. Moreover the kernel of $\calL_0 + \calL_0^*$ is composed of functions depending only on the position variables. The condition~\eqref{eq:euler lagrange numerics} then implies that there exists $g= g(q)$ in ${\rmL^2(\pinot)}$ such that
\begin{equation}
\label{eq:minimizer langevin}
	\Phi_{0,M} = - \Pi_M \calL_0^{-1} \Pi_0 R + g.
\end{equation}
\end{itemize}

The solution to~\eqref{eq:minimizer reversible} coincides with the result provided by the Galerkin method on the approximation space $V_M$. For~\eqref{eq:minimizer langevin} the optimal solution, for $g=0$, is given by the Galerkin method apart from a consistency error (see~\cite{Roussel17} for a detailed analysis). This justifies the use of Galerkin methods to determine a good approximation $\Phi_{0,M}$ of ${\Phi_0}$ in the general case.

For the Langevin equation the operator $\mathcal{L}$ is not coercive on ${\rmL_0^2(\pi)}$ so the associated rigidity matrix is not automatically invertible. The existence of a unique solution $\Phi_{0,M} \in V_M$ converging to ${\Phi_0} = -{\calL^{-1}} \Pi_0 R$ when $M \to +\infty$, as well as error estimates and a discussion on non-conformal approximations, can be found in~\cite{Roussel17}.

When a Galerkin method (or any other approximation method) is used, the error committed on ${\Phi_0}$ induces an error on the modified observable $\phi_\eta$. The modified asymptotic variance is then the sum of terms coming from~\eqref{eq:var order 2 eta exact} (depending on $\eta$) and terms coming from the approximation error (of order $\varepsilon$) committed on ${\Phi_0}$, as made precise in the following result.

\begin{theorem}
\label{th:var order 2 eta}
Fix $R \in \mathcal{S}$ and assume that $\Phi_0$ is approximated by $\Phi_{0,\varepsilon} = \Phi_0 + \varepsilon f$ with $f \in \mathcal{S}$ and $\varepsilon \geqslant 0$. Denote by $\phi_{\eta,\varepsilon} = R + \calL_\eta \Phi_{0,\varepsilon}$ the modified observable. Under Assumptions~1 to~5, for any $\eta_*,\varepsilon_* > 0$, there exists $E_{R,\eta_*,\varepsilon_*} > 0$ such that, for any $|\eta | \leqslant \eta_*$ and $| \varepsilon | \leqslant \varepsilon_*$,
\begin{equation}
\begin{aligned}
\label{eq:variance with cv}
\sigma_{\phi_{\eta, \varepsilon}, \eta}^2 &= 2 \varepsilon^2 \left\langle -\calL_0 f, f \right\rangle_0
- 2 \varepsilon \eta \left\langle (\calL_0 + \calL_0^*) f, \calL_0^{-1} \Pi_0 A R \right\rangle_0 \\
&\quad +2 \eta^2 \left\langle A R, -\calL_0^{-1} \Pi_0 A R \right\rangle_0
+ (\eta^3 + \varepsilon^3) C_{R, \eta, \varepsilon},
\end{aligned}
\end{equation}
with $| C_{R, \eta, \varepsilon} | \leqslant E_{R,\eta_*,\varepsilon_*}$.
\end{theorem}
The proof of this result can be read in Appendix~\ref{ap:proof Th2}. It shows that the variance is globally of order 2 with respect to both $\eta$ and $\varepsilon$. This suggests to take $\eta$ and $\varepsilon$ of the same order.
\begin{remark}
The dependence of $\widetilde C_{R,\eta_*,\varepsilon_*}$ with respect to $R$ can be made more explicit (see for instance the discussion in~\cite{Leimkuhler16}).
\end{remark}
The error committed on ${\Phi_0}$ can also arise from additional approximations on the right hand side of the Poisson problem, in situations when the observable $R_\eta = R_0 + \eta \widetilde R$ depends on $\eta$. A result similar to Theorem~\ref{th:var order 2 eta} can be obtained upon assuming that $\widetilde R \in \mathcal{S}$, where ${\Phi_0}$ is the solution to~\eqref{eq:poisson problem simple} with $R$ replaced by~$R_0$.

\begin{remark}
\label{rmk:variance irreversible}
Note that in the expression \eqref{eq:variance with cv} the error term $f$ only appears through $(\calL_0 + \calL_0^*) f$. This term may vanish even if $f$ is not identically zero. For example, for a Langevin process at equilibrium, $(\calL_0+\calL_0^*) f$ vanishes when $f$ is a function depending only on the positions.
\end{remark}

\section{One-dimensional Langevin dynamics}
\label{s:sinus potential}

We construct in this section a control variate for a one-dimensional system by solving a simplified Poisson equation using a spectral Galerkin method. The simplification consists in neglecting a non-equilibrium perturbation, the small parameter $\eta$ being the amplitude of this perturbation. We first present in Section~\ref{ss:sinus full dynamics} the model and define the quantity of interest, namely the mobility. We next construct in Section~\ref{ss:sinus simplified} the approximate control variate and conclude in Section~\ref{ss:sinus numerics} with some numerical results.

\subsection{Full dynamics}
\label{ss:sinus full dynamics}
We consider the following Langevin process on the state space $\mathcal{X} = 2 \pi \mathbb{T} \times \mathbb{R}$:
\begin{equation}
\label{eq:dynamics sin}
\left\{
\begin{aligned}
	\dd q_t &= \frac {p_t} m \, \dd t, \\
	\dd p_t &= (-v'(q_t)  + \eta ) \, \dd t - \frac \gamma m p_t \, \dd t + \sqrt{2 \gamma \beta^{-1}} \, \dd W_t,
\end{aligned}
\right.
\end{equation}
where $\gamma, m, \beta > 0$ and $v$ is a smooth $2 \pi$-periodic potential. The particle experiences a constant external driving of amplitude $\eta \in \mathbb{R}$. This force is not the gradient of a periodic function, so the system is out of equilibrium and the invariant measure is not known. We are interested in the average velocity $R(q,p)=\frac p m$ induced by the non-gradient force $\eta$, which can also be seen as a mass flux. The linear response of the average velocity with respect to the external driving is characterized by the mobility of the particle~\cite{Rodenhausen89}:
\begin{equation*}
	D = \lim_{\eta \to 0} \frac{\mathbb{E}_\eta[R]} \eta.
\end{equation*}
The generator of~\eqref{eq:dynamics sin} is the sum of the generator associated with the Langevin dynamics at equilibrium and of a non-equilibrium perturbation:
\begin{equation*}
\begin{aligned}
	\calL_\eta = \calL_0 + \eta \widetilde{\mathcal{L}},
\end{aligned}
\end{equation*}
where
\begin{equation*}
	\calL_0 = -v'(q) \partial_p + \frac p m \partial_q - \frac \gamma m p \partial_p + \gamma \beta^{-1} \partial_p^2, \qquad \widetilde{\mathcal{L}} = \partial_p.
\end{equation*}
In this setting the Lyapunov functions are defined for all $n\in \mathbb{N}$ as:
\begin{equation*}
	\forall (q,p) \in \mathcal{X}, \quad \mathcal{K}_n(q,p) = 1 + |p|^n,
\end{equation*}
and Assumptions~\ref{as:pieta}, \ref{as:lyapunov in L2}, \ref{as:generator eta} and \ref{as:generator simple} correspond to standard results for Langevin dynamics~\cite{Talay02, ReyBellet06_chap1, Kopec15, Leimkuhler16, Lelievre16}. Assumption~\ref{as:Lrep} trivially holds: the core space $\mathcal{S}$ is stable by $\widetilde{\mathcal{L}} = \partial_p$ by definition, while $\widetilde{\mathcal{L}}^* \mathbf{1} = (-\partial_p + \frac \beta m p) \mathbf{1} = \frac \beta m p \in {\rmL^2(\pinot)}$.

\subsection{Simplified dynamics and control variate}
\label{ss:sinus simplified}

Solving the Poisson problem $-\calL_\eta \Phi_\eta = R-\Pi_\eta[R]$ associated with the nonequilibrium dynamics is not practical because $\mathbb{E}_\eta[R]$ is not know. For this simple one-dimensional example it would still be technically doable. Since our purpose is however to illustrate both Theorems~\ref{th:var order 2 eta exact} and~\ref{th:var order 2 eta}, we do not follow this path. We therefore consider the control variate associated to a reference Poisson problem, namely
\begin{equation}
  \label{eq:simplified langevin 1D}
  -\calL_0 \Phi_0 = R.
\end{equation}
Note that the average drift vanishes at equilibrium: $\mathbb{E}_0[R] = 0$. Equation~\eqref{eq:simplified langevin 1D} cannot be solved analytically, but it is possible to approach its solution by a Galerkin method as explained in Section~\ref{ss:numerics poisson}. The modified observable is
\begin{equation*}
  \phi_{\eta,M} = R + \calL_\eta \Phi_{0,M},
\end{equation*}
with the notation of Theorem~\ref{th:var order 2 eta} (the error committed when estimating ${\Phi_0}$ is indexed by $M$ instead of $\varepsilon$). In practice we construct a basis $(e_m)_m$ of $V_M$ and write $\Phi_{0,M} = \sum_{m=1}^M a_m e_m$ where the coefficients $(a_m)_m \in \mathbb{R}^M$ are the solution of a linear system obtained from~\eqref{eq:euler lagrange numerics} (see~\cite{Roussel17}). The modified observable is then $\phi_{\eta,M} = R + \sum_{m=1}^M a_m \calL_\eta e_m$ where the functions $\calL_\eta e_m$ are explicit for appropriate choices of basis functions $(e_m)_m$; see Section~\ref{ss:sinus numerics}.

The mobility is estimated with an ergodic average of $R$ along a trajectory during a time~$T$, for a forcing $\eta$, by$\widehat D_{\eta, T} = \frac 1 \eta \widehat R_T$. This estimator has an expectation of order 1 and a large variance when $\eta$ is small and $T$ is large~\cite{Lelievre10, Lelievre16}:
\begin{equation*}
	\mathbb{E}[\widehat D_{\eta, T}] = D + \mathrm{O} \left( \eta + \frac 1 T\right), \qquad \mathrm{Var}[\widehat D_{\eta, T}] \sim \frac {\sigma_R^2} {\eta^2 T},
\end{equation*}
so that the relative statistical error scales as
\begin{equation*}
	\frac {\sqrt{\mathrm{Var}[\widehat D_{\eta, T}]}} {\mathbb{E}[\widehat D_{\eta, T}]} \sim \frac {\sigma_R} {D \eta \sqrt T}.
\end{equation*}
In order for this quantity to be small, the simulation time should be taken of order $T \sim \frac 1 {\eta^2}$, which is very large since $\eta$ is small.

When the Poisson problem is exactly solved, the modified observable is
\begin{equation*}
	\phi_{\eta} = R + \calL_\eta \Phi_0 = \eta \widetilde{\calL} \Phi_0,
\end{equation*}
which is proportional to $\eta$, so that the associated asymptotic variance scales as $\sigma_{\phi_\eta}^2 \sim \eta^2$. The relative statistical error for the mobility estimator $\widetilde D_{\eta, T} = \frac 1 \eta \widehat \phi_{\eta,T}$ is then bounded with respect to $\eta$:
\begin{equation*}
	\frac {\mathrm{Var}[\widetilde D_{\eta, T}]} {\mathbb{E}[\widetilde D_{\eta, T}]} \sim \frac 1 {D \sqrt T},
\end{equation*}
and the simulation time can be fixed independently of the value of $\eta$. Now, if an error of order $\varepsilon_M$ is committed on ${\Phi_0}$, the asymptotic variance of $\phi_{\eta, M}$ scales as $\eta^2 + \varepsilon_M^2$ so the relative statistical error on $\widetilde D_{\eta, T}$ is of order
\[
\frac {|\eta| + \varepsilon_M} {\sqrt{T} {|\eta|}} = \frac 1 {\sqrt{T}} \left( 1 + \frac {\varepsilon_M} {|\eta|} \right).
\]
This implies that the simulation time $T$ can be taken of order $\displaystyle{1+\left( \frac{\varepsilon_M}{\eta} \right)^2}$ instead of $\eta^{-2}$.

\subsection{Numerical results}
\label{ss:sinus numerics}

In order to simplify the numerical resolution of the Galerkin problem (see Section~\ref{ss:numerics poisson}) we consider the simple potential:
\begin{equation*}
  \forall q \in 2 \pi \mathbb{T}, \quad v(q) = 1 - \cos(q).
\end{equation*}
We construct $V_M$ using a tensorized basis made of weighted Fourier modes in position and Hermite modes in momenta. The particular weights of the Fourier modes are chosen so that the basis is orthogonal for the ${\rmL^2(\pinot)}$ scalar product. Obtaining error estimates on ${\Phi_0}-\Phi_{0,M}$ requires some work, see~\cite[Section~4]{Roussel17} for a detailed analysis and a precise expression of the basis. In the following we take either $M=15 \times 10$ basis elements ($15$ Fourier modes and $10$ Hermite modes), $M=7 \times 5$ or $5 \times 3$ basis elements. Estimating ${\Phi_0}$ allows to construct a control variate, and also to compute directly the mobility since the Green--Kubo formula~\cite{Kubo91} states that
\begin{equation}
  \label{eq:mobility algebra}
  D = \beta \left\langle R, -\calL_0^{-1} R \right\rangle_0 = \beta \left\langle R, \Phi_0 \right\rangle_0.
\end{equation}

In order to compute the mobility using Monte-Carlo simulations, we fix $\eta >0$ small and rely on the estimators $\widehat D_{\eta, T}$ or $\widetilde D_{\eta, T}$ defined in Section~\ref{ss:sinus simplified}. We are interested in the reduction of the asymptotic variance provided by our control variate, \textit{i.e.} comparing $\mathrm{Var}[\widetilde D_{\eta, T}]$ and $\mathrm{Var}[\widehat D_{\eta, T}]$. The Langevin dynamics is integrated over a time $T=2 \times 10^4$ with time steps $\Delta t = 0.02$ for an external forcing $\eta$ ranging from $0.01$ to $2.56$. The numerical integration is done with a Geometric Langevin Algorithm~\cite{BouRabee10}. This scheme ensures that the invariant probability measure is correct up to terms of order $\mathrm{O}({\Delta t}^2)$ at equilibrium. Moreover the transport coefficients estimated by linear response are also correct up to terms of order ${\Delta t}^2$ (see~\cite{Leimkuhler16}). The scheme writes:
\begin{equation}
  \label{eq:gla}
  \left\{ \begin{aligned}
    p^{k+{\nicefrac{1}{2}}} &= p^k + \left(- v'(q^k)+\eta\right) \frac {\Delta t} 2, \\
    q^{k+1} &= q^k + \frac {p^{k+{\nicefrac{1}{2}}}} m \Delta t, \\
    \widetilde p^{k+1} &= p^{k+{\nicefrac{1}{2}}} + \left(- v'(q^{k+1})+\eta\right) \frac {\Delta t} 2, \\
    p^{k+1} &= \alpha_{{\Delta t}} \widetilde p^{k+1} + \sqrt{m \beta^{-1} (1 - \alpha_{{\Delta t}}^2)} \ G^k,
	\end{aligned} \right.
\end{equation}
where the superscript $k$ is the iteration index, $\alpha_{{\Delta t}} = \exp(-\frac \gamma m \Delta t)$ and the $(G^k)_{k \in \mathbb{N}}$ are independent and identically distributed (i.i.d.) standard one-dimensional Gaussian random variables. The results which are reported are obtained for $m = \gamma = \beta = 1$.

\paragraph{Linear response.}

The results presented in Figure~\ref{fig:langevin mobility} (Left) show that the average velocity scales linearly with respect to the forcing for $\eta$ small, as predicted by linear response theory. The slope, which is the mobility $D$, matches the one computed using~\eqref{eq:mobility algebra} for $M$ large. An effective mobility is obtained by dividing the average velocity by the forcing, see Figure~\ref{fig:langevin mobility} (Right). We are interested in its limiting value for a small forcing. On the one hand the result is biased if $\eta$ is too large, but on the other hand the statistical error scales like $\frac 1 \eta$. For the standard observable we see that the optimal trade-off value of $\eta$ is around $0.2$ for the chosen simulation time $T$.  When using a control variate the variance is much smaller in the small forcing regime, so $\eta$ can be taken very small to reduce the bias while keeping the statistical error under control. We discuss next the estimation of the error bars plotted on Figure~\ref{fig:langevin mobility}.

\begin{figure}
\caption{Linear response for the standard MC simulation (black squares) compared to the version with control variate (blue, red) and to the asymptotic response $D \eta \approx 0.48 \eta$ (black line).}
\begin{center}
\includegraphics[scale=.58]{drift.pdf}
\includegraphics[scale=.58]{mobility.pdf}
\end{center}
\label{fig:langevin mobility}
\end{figure}

\paragraph{Correlation profiles.}

The asymptotic variance of time averages for an observable $\varphi \in \mathcal{S}$ writes, using the Green--Kubo formula~\cite{Kubo91},
\begin{equation*}
	\sigma_\varphi^2 = 2 \int_0^\infty C_\varphi(t) \, \dd t, \qquad C_\varphi(t) = \mathbb{E}[(\varphi(X_0)-\mathbb{E}[\varphi]) (\varphi(X_t)-\mathbb{E}[\varphi])],
\end{equation*}
where the autocorrelation function $C_\varphi$ involves an expectation over all initial condition $X_0 = (q_0,p_0) \in \mathcal{X}$ distributed according to the invariant probability measure~$\pi_0$ at equilibrium, and for all realizations of the dynamics with generator~$\mathcal{L}_0$. The integrability of $C_\varphi(t)$ can be guaranteed when the semi-group $\mathrm{e}^{t \mathcal{L}_0}$ decays sufficiently fast in ${\rmL_n^\infty}$~\cite{Lelievre16}. The function $C_\varphi$ is characterized by three major features explaining the value of the asymptotic variance; see Figure~\ref{fig:correlation cartoon} for an illustration.
\begin{enumerate}
\item[(i)] The first one is the amplitude of the signal $\| \varphi - \mathbb{E}[\varphi] \|_{\mathrm{L}^2(\pi_0)}^2 = C_\varphi(0)$ corresponding to the value of the autocorrelation at $t=0$.
\item[(ii)] The second one is the characteristic decay time $\tau$ of the autocorrelation, which can be related to the decay of its exponential envelope.
\item[(iii)] The last one is the presence of anticorrelations, which arise only for non-reversible dynamics such as Langevin dynamics.
\end{enumerate}
A proper estimation of the asymptotic variance requires to compute autocorrelation profiles on a sufficiently long time interval $[0,t_{\rm deco}]$ (here $t_{\rm deco}=6$). One can check a posteriori that this time is sufficient by looking at the convergence of the cumulated autocorrelation $t \mapsto \int_0^t C_\varphi$ toward its limit $\frac {\sigma_\varphi^2} 2 = \int_0^\infty C_\varphi$ (see Appendix~\ref{ap:variance} for more details on the variance estimators, and the computation of error bars for these quantities).

\begin{figure}
\caption{Illustrative autocorrelation profile. The dashed line is the exponential envelope of the correlation function.}
\begin{center}
\includegraphics[scale=.58]{correlation_cartoon}
\end{center}
\label{fig:correlation cartoon}
\end{figure}

Figure~\ref{fig:langevin autocorrelation} compares the autocorrelation profile of the velocity with the ones for the modified observables, for two different Galerkin basis sizes~$M$ and two different forcing amplitudes~$\eta$. For a small forcing $\eta=0.08$ the two modified observables have an amplitude and a decorrelation time which are both much smaller than for the standard velocity. Note that the two modified observables do not exhibit any anti-correlation, contrarily to the velocity observable. The cumulated plots show that the control variates drastically reduce the asymptotic variance in this case, especially for the one based on a more accurate Galerkin approximation. For a larger value $\eta=1.28$ the modified observables have a significantly larger amplitude (\textit{i.e.} $C_\varphi(0)$ is larger), especially in the case of a low accuracy~$M$. However the decorrelation times are small and there is anti-correlation, resulting in a reasonable variance reduction in both cases.

\begin{figure}
\caption{Left: Autocorrelation profile of the velocity compared to the one of the modified observables for two different accuracies, either for a small forcing $\eta=0.08$ (top) or a larger one $\eta=1.28$ (bottom). Right: Corresponding cumulated autocorrelations. The limit value is half the asymptotic variance of the observable.}
\begin{center}
\includegraphics[scale=.58]{correlations_eta0.08}
\includegraphics[scale=.58]{cumulated_eta0.08}\\
\includegraphics[scale=.58]{correlations_eta1.28}
\includegraphics[scale=.58]{cumulated_eta1.28}
\end{center}
\label{fig:langevin autocorrelation}
\end{figure}

\paragraph{Asymptotic variances.}

Let us now compute the asymptotic variance for a whole range of Galerkin accuracies~$M$ and forcing amplitudes $\eta$. The results presented in Figure~\ref{fig:langevin variance} confirm that for a very accurate Galerkin resolution the variance of the modified observable scales as $\eta^2$ with the prefactor $\alpha = \left\langle \Pi_0 A R, -\calL_0^{-1} \Pi_0 A R \right\rangle_0$ predicted theoretically in Theorem~\ref{th:var order 2 eta exact}. This prefactor has been computed independently by solving~\eqref{eq:simplified langevin 1D} using a Galerkin method and plugging this approximation in~\eqref{eq:mobility algebra}. When the Galerkin discretization is not sufficiently accurate, the variance reaches a plateau in the region of small forcings as predicted by Theorem~\ref{th:var order 2 eta}.

\begin{figure}
\caption{Asymptotic variance of the velocity (black squares) compared to its counterpart when using a control variate (blue, grey, red) and to the reduced variance (black line) predicted theoretically ($\alpha \approx 0.53$ computed with a Galerkin discretization).}
\begin{center}
\includegraphics[scale=.58]{variance.pdf}
\end{center}
\label{fig:langevin variance}
\end{figure}

\section{Thermal transport in atom chains}
\label{s:quasi harmonic}

Thermal transport in one-dimensional systems has been the topic of many investigations, both from theoretical and numerical points of view~\cite{Bonetto00, Lepri16, Dhar08}. Determining which microscopic ingredients influence the scaling of the conductivity with respect to the length of the chain is still an active line of research. Studying numerically this scaling requires to simulate chains of thousands of particles. In these systems the temperature gradient and the thermal flux are both very small, which induces large statistical errors when estimating the conductivity. Introducing variance reduction techniques not requiring the knowledge of the invariant probability measure could alleviate (at least partly) these difficulties.

\subsection{Full dynamics}

\subsubsection{Equations of motion}

We consider a chain composed of $N$ particles interacting through a nearest-neighbor potential~$v$ (see Figure~\ref{fig:thermal chain}). The evolution is dictated by a Hamiltonian dynamics and a thermalization mechanism at the boundaries, where the first and the last particles are submitted to Ornstein-Uhlenbeck processes, at temperatures $T_{\rmL}$ and $T_{\rmR}$ respectively. The unknowns are the momenta~$p=(p_n)_{1 \leqslant n \leqslant N}$ of the particles and the interparticle distances $r=(r_n)_{1 \leqslant n \leqslant N-1}$. In these variables, the dynamics reads:
\begin{equation}
\label{eq:dynamics chain}
\left\{ \begin{aligned}
	\dd r_n &= \frac 1 m (p_{n+1} - p_n) \, \dd t, \\
	\dd p_1 &= v'(r_1) \, \dd t - \frac \gamma m p_1 \dd{t} + \sqrt{2 \gamma T_{\rmL}} \dd W_t^L, \\
	\dd p_n &= (v'(r_n) - v'(r_{n-1})) \, \dd t, \\
	\dd p_N &= -v'(r_{N-1}) \, \dd t - \frac \gamma m p_N \dd{t} + \sqrt{2 \gamma T_{\rmR}} \dd W_t^R,
\end{aligned} \right.
\end{equation}
where $m>0$ is the mass of a particle, $\gamma>0$ is the friction coefficient, $T_{\rmL} \geqslant T_{\rmR}$ and $W_t^L$, $W_t^R$ are two independent standard one-dimensional Brownian motions. Notice that the ends of the chain are free. The Hamiltonian of the system is the sum of the potential and kinetic energies:
\begin{equation*}
	H(r,p) = V(r) + \sum_{n=1}^N \frac {p_n^2}{2 m}, \quad V(r) = \sum_{n=1}^{N-1} v(r_n).
\end{equation*}
The infinitesimal generator of the dynamics~\eqref{eq:dynamics chain} reads
\begin{equation*}
\begin{aligned}
	\mathcal{L} &= \frac 1 m \sum_{n=1}^{N-1} (p_{n+1} - p_n) \partial_{r_n} + \sum_{n=1}^N (v'(r_n) - v'(r_{n-1})) \partial_{p_n} \\
	&\quad - \frac \gamma m p_1 \partial_{p_1} + \gamma T_L \partial_{p_1}^2 - \frac \gamma m p_N \partial_{p_N} + \gamma T_R \partial_{p_N}^2,
\end{aligned}
\end{equation*}
using the convention $v'(r_0) = v'(r_N) = 0$.
\begin{figure}
\begin{center}
\tikzsetnextfilename{scales}
\begin{tikzpicture}[scale=.4]
\input{Img/chain_free.tex}
\end{tikzpicture}
\caption{Heat transport in a one-dimensional chain.}
\label{fig:thermal chain}
\end{center}
\end{figure}

\subsubsection{Properties of the dynamics}
\label{sss:chain properties}

\newcommand{\rmL_\theta^\infty}{\mathrm{L}_\theta^\infty}
\newcommand{\rmL_{\theta,0}^\infty}{\mathrm{L}_{\theta,0}^\infty}

Let us recall some properties of the dynamics~\eqref{eq:dynamics chain} which hold under the following assumption.
\begin{assumption}
\label{as:potential chain}
The interaction potential $v$ is $\mathcal{C}^\infty$ and there exist $k \geqslant 2$ and $a > 0$ such that
\begin{equation*}
	\forall r_1 \in \mathbb{R}, \quad \lim_{\tau \to +\infty} \tau^{-k} v(\tau r_1) = a |r_1|^k, \quad \lim_{\tau \to +\infty} \tau^{1-k} v'(\tau r_1) = k a |r_1|^{k-1} \mathrm{sign}(r_1).
\end{equation*}
Moreover the interaction potential is not degenerate: for any $r_1 \in \mathbb{R}$ there exists $m=m(r_1) \geqslant 2$ such that $\partial_m v(r_1) \neq 0$.
\end{assumption}
These conditions hold for the potentials we use in the numerical simulations reported in Section~\ref{ss:chain numerics}. When Assumption~\ref{as:potential chain} holds, the dynamics admits a unique invariant probability measure~$\pi$~(see \cite{Carmona07}). This invariant probability measure is explicit when the chain is at equilibrium ($T_{\rmL} = T_{\rmR} = \beta^{-1}$), in which case it has the tensorized form
\begin{equation}
\label{eq:measure chain}

	\pi_{\rm eq}(\dd r \, \dd p) = Z_\beta^{-1} \exp(- \beta \left( \frac {|p|^2} {2 m} + V(r) \right)) \, \dd r \, \dd p,
\end{equation}
where $Z_\beta^{-1}$ is a normalization constant. Let us emphasize that the reference system considered later on in Section~\ref{ss:chain simplified} is not at equilibrium. Following the framework considered in~\cite{Carmona07} (which is itself based on~\cite{ReyBellet06_chap2}), we consider in this section the Lyapunov functions $\mathcal{K}_\theta = \mathrm{e}^{\theta H}$. There exist $\theta_* > 0$ such that $\mathcal{K}_\theta \subset \mathrm{L}^2(\pi)$ for any $\theta \in [0,\theta_*)$. The functional spaces we use are also indexed by the continuous parameter $\theta \in [0,\theta_*)$:
\begin{equation*}
  \mathrm{L}_\theta^\infty = \left\{ \varphi \ {\rm measurable} \ \big| \ \| \varphi \, \mathrm{e}^{-\theta H} \|_{\mathrm{L}^\infty} < + \infty \right\},
\end{equation*}
and the space $\mathcal{S}$ is defined similarly to~\eqref{eq:core}. For $\theta \in [0,\theta_*)$, we also define the vector space $\rmL_{\theta,0}^\infty$ of functions of $\rmL_\theta^\infty$ with mean zero with respect to $\pi$. One can prove the exponential decay of the semi-group on the associated functional space $\rmL_{\theta,0}^\infty$ (see~\cite{Carmona07}): for any $\theta \in [0,\theta_*)$, there exist $C, \lambda > 0$ such that, for any $\varphi \in \rmL_{\theta,0}^\infty$,
\begin{equation*}
	\forall t \geqslant 0, \quad \| \mathrm{e}^{t \mathcal{L}} \varphi \|_{\mathrm{L}_\theta^\infty} \leqslant C \mathrm{e}^{-\lambda t} \| \varphi \|_{\mathrm{L}_\theta^\infty}.
\end{equation*}
This implies that $\mathcal{L}$ is invertible on $\rmL_{\theta,0}^\infty$, and that its inverse is bounded.

\paragraph{Validity of Assumptions 1 and 2.}
Assumption~\ref{as:pieta} holds true since there exist a unique invariant probability measure with positive density with respect to the Lebesgue measure, and the generator of the dynamics is hypoelliptic~\cite{Carmona07,Kliemann87,ReyBellet06_chap2}. The first part of Assumption~\ref{as:lyapunov in L2} is also satisfied for $\theta \in [0,\theta_*)$. Note that at equilibrium ($T_{\rmL}=T_{\rmR}=\beta^{-1}$) the invariant probability measure is explicit and $\theta^* = \beta/2$. The product of two Lyapunov functions $\mathcal{K}_\theta$ and $\mathcal{K}_{\theta'}$ is in a Lyapunov space only if $\theta + \theta' < \theta_*$, so the second part of Assumption~\ref{as:lyapunov in L2} is not satisfied.

\subsubsection{Heat flux and conductivity}

When studying heat transport in atom chains the typical quantity of interest is the thermal flux through the chain:
\begin{equation}
\label{eq:heat flux}
	\forall n \in [1, N-1], \quad j_n(r,p) = -\frac {p_n + p_{n+1}} 2  v'(r_n),
\end{equation}
see for example the review~\cite{Lepri03} on thermal transport in low-dimensional lattices for further background material. We also make use of the two boundary elementary fluxes:
\begin{equation}
\label{eq:heat flux boundaries}
	j_0(r,p) = \frac \gamma m \left(T_{\rmL} - \frac {p_1^2} m\right), \qquad j_N(r,p) = \frac \gamma m \left( \frac {p_N^2} m - T_{\rmR} \right).
\end{equation}
The definition of the elementary fluxes $j_n$ is motivated by the local energy balance, centered on particle~$n$:
\begin{equation}
\label{eq:energy balance}
	\forall n \in [1, N], \quad \mathcal{L} \varepsilon_n = j_{n-1} - j_n, \qquad \varepsilon_n(r,p) = \frac {v(r_{n-1})} 2 + \frac {p_n^2}{2 m} + \frac {v(r_n)} 2.
\end{equation}
The quantity $\mathcal{L} \varepsilon_n$ is of mean zero since it is in the image of the generator. Therefore the elementary fluxes $j_n$ all have the same stationary values:
\begin{equation}
\label{eq:equality fluxes}
	\mathbb{E}_\pi[j_0] = \mathbb{E}_\pi[j_1] = \cdots = \mathbb{E}_\pi[j_N].
\end{equation}
Any linear combinations of such fluxes, namely
\begin{equation}
\label{eq:linear combination flux}
	J_\lambda = \sum_{n=0}^N \lambda_n j_n, \qquad \sum_{n=0}^N \lambda_n = 1,
\end{equation}
has the same stationary value. The most common choice is the spatial mean
\begin{equation}
\label{eq:heat flux sum}
	\widetilde R = \frac 1 {N-1} \sum_{n=1}^{N-1} j_n.
\end{equation}
We call the latter observable "standard heat flux" in the sequel. Notice that it does not depend on the boundary fluxes $j_0$ and $j_N$. The linear response of $\widetilde R$ (or of any flux $J_\lambda$) with respect to the temperature gradient $\frac {T_{\rmL}-T_{\rmR}}{N-1}$ defines the effective conductivity:
\begin{equation}
\label{eq:def conductivity}
	\kappa = \frac{N-1}{T_{\rmL}-T_{\rmR}}\mathbb{E}_\pi[\widetilde R],
\end{equation}
which is here the transport coefficient of interest. There exist infinitely many observables having the same expectation as $\widetilde R$, see \eqref{eq:linear combination flux} for example. Let us first discuss the choice of observable for the heat flux, before trying to reduce its variance by constructing a control variate.

\paragraph{Asymptotic variance of the heat fluxes at equilibrium.}

The chain is supposed to be at equilibrium in all this paragraph ($T_{\rmL}=T_{\rmR}=\beta^{-1}$). The conclusions remain unchanged for nonequilibrium systems when $T_{\rmL}-T_{\rmR}$ is small since the results are only perturbed to first order with respect to this quantity. In the remainder of Section~\ref{s:quasi harmonic}, an index 'eq' refers to the equilibrium dynamics and to the equilibrium invariant probability measure~$\pi_{\rm eq}$. In this setting the asymptotic variance $\sigma_{\widetilde R, {\rm eq}}^2$ for the standard heat flux $\widetilde R$ is not smaller than the one associated with any elementary flux $(j_n)_{1\leqslant n \leqslant N-1}$. These two variances are in fact equal, as made precise in the following proposition (similar in spirit to Remark~\ref{rmk:variance irreversible}).

\begin{prop}
\label{prop:variance invariance}
Consider an observable $\varphi \in \mathcal{S}$ and a function $U \in \mathcal{S}$ which does not depend on $p_1$ and $p_N$. Then adding $\mathcal{L} U$ to the observable does not modify the variance:
\begin{equation}
  \label{eq:equality_variances}
  \sigma_{\varphi + \mathcal{L} U, {\rm eq}}^2 = \sigma_{\varphi, {\rm eq}}^2.
\end{equation}
\end{prop}
\begin{proof}
At equilibrium, the invariant probability measure $\pi_{\rm eq}$ is explicit. The symmetric part of the generator can then be computed and it corresponds to the fluctuation-dissipation part of the process:
\begin{equation}
  \label{eq:def_LFD}
  \frac 1 2 (\mathcal{L} + \mathcal{L}^*) = \calL_\mathrm{FD} := - \frac{\gamma}{\beta} \left(  \partial_{p_1}^* \partial_{p_1} + \partial_{p_N}^* \partial_{p_N} \right),
\end{equation}
where adjoints are considered on ${\rmL^2(\pieq)}$. When $U$ does not depend on $p_1$ nor $p_N$, it holds $\calL_\mathrm{FD} U = 0$. The claimed result then follows from Lemma~\ref{lemma:variance LU}.
\end{proof}

The equality~\eqref{eq:equality_variances} is perturbed by terms of order $T_{\rmL}-T_{\rmR}$ for out of equilibrium dynamics according to linear response theory. Upon taking for $U$ a linear combination of the energies $(\varepsilon_n)_{2 \leqslant n \leqslant N-1}$, we directly obtain, thanks to~\eqref{eq:energy balance}, that all the fluxes of the form~\eqref{eq:linear combination flux} which do not depend on the boundary fluxes ($\lambda_0=\lambda_N=0$) share the same asymptotic variance at equilibrium; in particular
\begin{equation}
  \label{eq:flux same asymptotic variance}
  \forall 1 \leqslant n \leqslant N-1, \quad \sigma_{j_n, {\rm eq}}^2 = \sigma_{\widetilde R, {\rm eq}}^2.
\end{equation}

\begin{remark}
\label{rmk:flux eq}
Linear response theory indicates that the previous asymptotic variances are related to the conductivity through the Green--Kubo formula~\cite{Kubo91}:
\begin{equation}
  \label{eq:variance and conductivity}
  \sigma_{\widetilde R, {\rm eq}}^2 = \frac {2 \kappa}{\beta^2 (N-1)}.
\end{equation}
We show in Appendix~\ref{ap:proof variance end} that, in this equilibrium situation, the variance of the two boundary fluxes $j_0$ and $j_N$ is also related to the conductivity as:
\begin{equation}
  \label{eq:variance boundary}
  \sigma_{j_0, {\rm eq}}^2 = \frac \gamma {m \beta^2} - \frac {2 \kappa}{\beta^2 (N-1)}, \qquad
  \sigma_{j_N, {\rm eq}}^2 = \frac \gamma {m \beta^2} - \frac {2 \kappa}{\beta^2 (N-1)}.
\end{equation}
Apart from special cases such as integrable systems (as the harmonic system considered in Appendix~\ref{ap:harmonic}), we generically observe numerically that $\frac {\kappa(N)} {N-1} \xrightarrow[N \to \infty]{} 0$. The boundary fluxes therefore have (asymptotically in~$N$) a larger variance than bulk fluxes. Since the variances are perturbed to first order in $T_{\rmL}-T_{\rmR}$ in nonequilibrium situations, the same conclusion holds for temperature differences which are not too large.
\end{remark}

In the next section we construct a modified observable by adding a control variate to the reference observable $R = \frac 1 2 (j_0 + j_N)$. This particular heat flux does not depend on the potential energy function~$v$, which simplifies the computation of the control variate (see Appendix~\ref{ap:harmonic}). In Appendix~\ref{ap:justification flux} we prove using Proposition~\ref{prop:variance invariance} that, at equilibrium, the asymptotic variance of the resulting modified observable does not depend on the choice of the reference observable~$R$.

\subsection{Simplified dynamics and control variate}
\label{ss:chain simplified}

We split the interaction potential into a harmonic part with parameters $\hat \omega > 0$ and $\hat r \in \mathbb{R}$, and an anharmonic part~$w$:
\begin{equation}
  \label{eq:potential decomposition}
  v(r_1) = v_0(r_1) + w(r_1), \qquad v_0(r_1) = \frac 1 2 m \hat \omega^2 (r_1-\hat r)^2.
\end{equation}
The potential energy is then decomposed as
\begin{equation*}
  V(r) = V_0(r) + W(r), \qquad V_0(r) = \sum_{n=1}^{N-1} v_0(r_n), \qquad W(r) = \sum_{n=1}^{N-1} w(r_n).
\end{equation*}
Following the general strategy outlined in Section~\ref{s:strategy} we decompose the generator as
\begin{equation*}
	\mathcal{L} = \calL_0 + \widetilde{\calL},
\end{equation*}
where $\calL_0$ is the generator of the harmonic chain corresponding to the harmonic interaction potential~$v_0$ and $\widetilde{\calL}$ is the generator of the anharmonic perturbation:
\begin{equation*}
	\widetilde{\calL} = \sum_{n=1}^N \left( w'(r_n) - w'(r_{n-1}) \right) \partial_{p_n},
\end{equation*}
with the same convention $w'(r_0) = w'(r_N) = 0$ as for the potential $v$. We simplify the Poisson problem for the optimal control variate
\begin{equation*}
	-\mathcal{L} \Phi = R - \mathbb{E}[R],
\end{equation*}
into the harmonic Poisson problem
\begin{equation}
\label{eq:simplified poisson chain}
	-\calL_0 {\Phi_0} = R - \mathbb{E}_0[R].
\end{equation}
Note that the observable $R$ does not depend on the potential, contrarily to other heat fluxes such as $\widetilde R$ in~\eqref{eq:heat flux sum}, so that the right hand side of the Poisson equation needs not be changed when looking for an approximate control variate. Equation~\eqref{eq:simplified poisson chain} can be solved analytically for ${\Phi_0}$. In Appendix~\ref{ap:harmonic} we show that
\begin{equation}
\begin{aligned}
\label{eq:cv chain}
	\mathbb{E}_0[R] &= \frac {\nu^2} {1 + \nu^2} \frac {\gamma (T_{\rmL}-T_{\rmR})} {2 m}, \\
	{\Phi_0}(r,p) &= \frac m {2 \gamma \left( 1 + \nu^2 \right)} \left[-\hat \omega^2 \sum_{n=1}^{N-1} (r_n-\hat r) (p_n + p_{n+1}) + \frac \gamma {2m^2} \left( p_N^2 - p_1^2 \right) \right] + C,
\end{aligned}
\end{equation}
where $\nu = \frac {m \hat \omega} \gamma$ and $C \in \mathbb{R}$ is such that $\mathbb{E}_0[{\Phi_0}] = 0$. The modified observable is therefore
\begin{equation}
\label{eq:modified obs chain}
\begin{aligned}
	&\quad (R+\mathcal{L} {\Phi_0})(r,p) = \mathbb{E}_0[R] + \widetilde{\calL} {\Phi_0}(r,p) \\
		&= \frac 1 {2 ( 1 + \nu^2 )} \left[ \nu \hat \omega (T_{\rmL}-T_{\rmR}) - \nu \hat \omega \sum_{n=1}^{N-1} (r_n-\hat r) \left( w'(r_{n+1}) - w'(r_{n-1}) \right)
		 - \left( \frac {p_1} m w'(r_1) + \frac {p_N} m w'(r_{N-1}) \right) \right] \\
		&= \frac 1 {2(1+\nu^2)} \left[ \nu \hat \omega (T_{\rmL}-T_{\rmR}) - \sum_{n=1}^{N-1} \left( \widetilde v_{n+1}(r,p) - \widetilde v_{n-1}(r,p) \right) w'(r_n) \right],
\end{aligned}
\end{equation}
where
\begin{equation*}
\widetilde v_n(r,p) = \left\{
\begin{array}{cl}
  -\frac {p_1} m \quad &\mbox {if} \quad n=0, \\
  -\nu \hat \omega (r_n-\hat r) \quad &\mbox {if} \quad 1 \leqslant n \leqslant N-1, \\
  \frac {p_N} m \quad &\mbox {if} \quad  n=N.
\end{array}
\right.
\end{equation*}
Notice that, by construction, this observable is constant when the chain is harmonic (\emph{i.e.} $w=0$).

\paragraph{Harmonic fitting.}

For a given pair potential $v=v(r)$, there is some freedom in the decomposition~\eqref{eq:potential decomposition}, namely the choice of the parameters $\hat \omega$ and $\hat r$. The optimal choice would be such that the variance of the modified observable~\eqref{eq:modified obs chain} is minimal, but this condition is not practical. A possible (and simpler) heuristic is to choose these coefficients in order to minimize the $\mathrm{L}^2(\pi_{\rm eq})$ norm of the anharmonic force $-\nabla W$ at equilibrium, namely when $T_{\rmL}=T_{\rmR}=\beta^{-1}$. In view of the tensorized form~\eqref{eq:measure chain} of the invariant probability measure at equilibrium,
\begin{equation*}
  \| \nabla W(r) \|_{\rm eq}^2 = (N-1) z_\beta^{-1} \int_\mathbb{R} \left( v'(r_1) - m \hat \omega^2 (r_1-\hat r) \right)^2\, \mathrm{e}^{- \beta v(r_1)} \dd r_1,
\end{equation*}
where $z_\beta = \int_\mathbb{R} \mathrm{e}^{- \beta v(r_1)} \dd r_1$. Therefore the minimization problem defining $\hat r$ and $\hat \Omega = m \hat \omega^2$ writes
\begin{equation}
\label{eq:minimization}
	\argmin_{\hat \omega, \hat r} \int_\mathbb{R} \left( v'(r_1) - \hat \Omega (r_1-\hat r) \right)^2\, \mathrm{e}^{- \beta v(r_1)} \dd r_1.
\end{equation}
There exists a minimizer $(\hat r, \hat \Omega)$ since the function to be minimized is continuous and coercive; uniqueness is proved in Appendix~\ref{ap:proof minimizer fpu}. Define the moments of the marginal measure for inter-particle distances as:
\begin{equation*}
	\mathcal{M}_k = \int_\mathbb{R} r_1^k \mathrm{e}^{- \beta v(r_1)} \dd r_1.
\end{equation*}
The Euler-Lagrange equation associated with~\eqref{eq:minimization} provides the expression of the minimizer (see Appendix~\ref{ap:proof minimizer fpu}):
\begin{equation}
  \label{eq:rhar omegahar values}
  \hat r = \frac{\mathcal{M}_1}{\mathcal{M}_0}, \qquad \hat \Omega = m \hat \omega^2 = \beta^{-1} \frac{\mathcal{M}_0^2}{\mathcal{M}_0 \mathcal{M}_2 - \mathcal{M}_1^2},
\end{equation}
where, by the Cauchy-Schwarz inequality, $\mathcal{M}_0 \mathcal{M}_2 - \mathcal{M}_1^2 > 0$ for any continuous potential $v$ which is not constant.

\paragraph{Validity of Assumptions 3 to 5.}

The standard way to verify Assumption~\ref{as:generator eta} would be to show that the  coefficients in the Lyapunov condition exhibited in~\cite{Carmona07} depend continuously on perturbations of the potential $v$. This is not straightforward, especially if the exponent~$k$ introduced in Assumption~\ref{as:potential chain} is discontinuous with the perturbation amplitude at $\eta=0$ (which corresponds to the harmonic chain). For example, for the FPU potential~\eqref{eq:fpu}, $k=2$ for the harmonic chain ($\eta=0$) and $k=4$ for $\eta >0$.  We show that Assumption~\ref{as:generator simple} holds under Assumption~\ref{as:potential chain} in Appendix~\ref{ap:ass 4}. Moreover it is clear that $\mathcal{S}$ is stable by $\widetilde{\mathcal{L}}$. A simple computation shows that $\pi_0$ is a Gaussian probability measure, which implies that $\widetilde{\mathcal{L}}^* \mathbf{1} \in {\rmL^2(\pinot)}$. Therefore Assumption~\ref{as:Lrep} holds as well.

\subsection{Numerical results}
\label{ss:chain numerics}

We consider a Fermi-Pasta-Ulam (FPU) potential
\begin{equation}
  \label{eq:fpu}
  v(r_1) = \frac a 2 r_1^2 + \frac b 3 r_1^3 + \frac c 4 r_1^4,
\end{equation}
where $c=\frac {b^2}{3 a}$ is such that $v''(r_1) = a + 2 b r_1 + 3 c r_1^2 = \frac 1 a (a + br_1)^2$ is positive except at a single point where it vanishes. This choice makes the potential both asymmetric and convex. Symmetric potentials indeed exhibit special behaviors~\cite{Spohn14}, whereas we want to be as general as possible. On the other hand, non-convex potentials are not typical in the literature on the computation of thermal conduction in one-dimensional chains. In the following we fix $a=1$ and vary $b$ only. The parameters $\hat r$ and $\hat \omega$ are given by~\eqref{eq:rhar omegahar values}, where the moments $\mathcal{M}_k$ are computed using one-dimensional numerical quadratures.

The system is simulated for a time $T=10^8$ with time steps $\Delta t = 10^{-2}$, and with $m=\gamma=1$. The atom chain is simulated either at equilibrium with $T_{\rmL}=T_{\rmR}=2$, or for $T_{\rmL}=3$ and $T_{\rmR}=1$. The dynamics is discretized using a Geometric Langevin Algorithm scheme as in~\eqref{eq:gla}. The estimator of the asymptotic variance is made precise in Appendix~\ref{ap:variance}. The decorrelation time is set to $t_{\rm deco} = 3 N$ for the standard flux $R$ and to $t_{\rm deco} = 32$ for the modified observable.

We plot in Figure~\ref{fig:autocorrelation chain} the autocorrelation profiles of the heat flux for a chain of size $N=128$ at equilibrium, for two different anharmonicities. Let us comment this picture in greater detail. The results first show that the chosen decorrelation time $t_{\rm deco}$ is sufficiently large. Similar plots were used to check this is also the case for all the range of anharmonicities $b$ and numbers of particles $N$ we consider. They can also be used to understand the eventual variance reduction granted by the control variate~\eqref{eq:cv chain}. For a small anharmonicity $b=0.08$ we see that the asymptotic variance, which is twice the limit of the cumulated autocorrelation (right plot) is greatly reduced for two reasons. First, the signal amplitude, which is related to the autocorrelation value at $t=0$ (left plot), is slightly smaller, and (right plot) the contribution of the times $0\leqslant t \leqslant 2$ is twice smaller for the modified flux. Second, and that is the actual reason for the variance reduction, there is anticorrelation for $2 \leqslant t \leqslant 5$. For a larger anharmonicity it appears that the amplitude of the modified observable is much larger, but this is compensated by a long-time anticorrelation. The resulting asymptotic variance of the modified flux is slightly smaller than the one of the standard flux. The plots are essentially the same out of equilibrium, when $T_{\rmL} = 3$ and $T_{\rmR}=1$ (numerical results not presented here).

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{nice_correlation_N128}
\includegraphics[scale=0.5]{cumulated_nice_correlation_N128}
\caption{Comparison of results obtained either with the elementary flux $j_n$~\eqref{eq:heat flux} or with the modified observable~\eqref{eq:modified obs chain}, for two different anharmonicities $b=0.08, 0.64$, and for a chain of $N=128$ particles at equilibrium. Left: Autocorrelation profile. Right: Cumulated autocorrelation $t \mapsto \int_0^t C$ at longer times.}
\label{fig:autocorrelation chain}
\end{center}
\end{figure}

The asymptotic variances, with associated error bars (see Appendix~\ref{ap:variance}), are plotted on Figure~\ref{fig:variance chain} for a whole range of anharmonicities $b$ and numbers of particles $N$. The two left plots are at equilibrium ($T_{\rmL}=T_{\rmR} = 2$) while the two right plots are out of equilibrium ($T_{\rmL}=3$ and $T_{\rmR} = 1$). We check that the variances are extremely similar in both cases, which is expected by linear response theory. We observe that the asymptotic variance of the modified flux scales as $b^2$ for $b \ll 1$, as expected from Theorem~\ref{th:var order 2 eta exact}, providing an excellent variance reduction in this case. Note that, in the limit $b \to 0$, the variance of the standard flux tends to $\frac {\gamma \nu^2}{m \beta^2(1+\nu^2)}=2$ (since $\nu=\frac{m \omega}{\gamma}=1$ here), which is the theoretical value predicted at equilibrium in view of the expression of the mean flux for a harmonic chain (see Equations~\eqref{eq:cv chain},~\eqref{eq:variance and conductivity}, and~\eqref{eq:def conductivity}). The modified flux can sometimes have a larger variance than the standard one, for example in the regime $b=\mathrm{O}(1)$ and $N$ large. Note that for the particular choice $\hat \omega=0$ the modified flux is $\frac 1 2 (j_1 + j_{N-1})$, which has the same asymptotic variance as the standard flux~$\widetilde R$ according to~\eqref{eq:flux same asymptotic variance}. Therefore, for any set of parameters, there exist an optimal choice of the coefficients $\hat \omega, \hat r$ providing a modified flux whose asymptotic variance is smaller or equal to its counterpart without control variate. In the present application these two coefficients are instead chosen according to the heuristic~\eqref{eq:rhar omegahar values}, leading to a degradation of the asymptotic variance in certain cases.

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{niceVarianceFctb_DT0}
\includegraphics[scale=0.5]{niceVarianceFctb_DT1}\\
\includegraphics[scale=0.5]{nice_variance_ratio_DT0}
\includegraphics[scale=0.5]{nice_variance_ratio_DT1}
\caption{Comparison of variances for the standard and modified fluxes. The reported variance ratio corresponds to the modified variance divided by the standard one. Left: Chain at equilibrium ($T_{\rmL}=T_{\rmR}=2$). Right: Chain out of equilibrium ($T_{\rmL}=3$ and $T_{\rmR}=1$).}
\label{fig:variance chain}
\end{center}
\end{figure}

\section{Solvated dimer under shear}
\label{s:dimer}

\newcommand{\nabla_{q_i}}{\nabla_{q_i}}
\newcommand{\Delta_{q_i}}{\Delta_{q_i}}

A solvated dimer is a pair of bonded particles in a bath constituted of many other particles. It serves as a prototypical model of a molecule in solvent (e.g. peptide in water). This model has been used in the context of free energy computations~\cite{Lelievre10}. We apply to this system an external shearing force as in~\cite{Joubaud12}, also coined sinusoidal transverse field in the physics literature (see~\cite[Section 9.1]{Todd17} and \cite{Evans13_book, Todd07}), so that the invariant measure of the system is not known. A typical question is the influence of the shear force on the average bond length of the dimer.

\subsection{Full dynamics}

We consider $N$ particles in a two-dimensional box of length $L$ with periodic boundary conditions, with positions $q=(q_1,\cdots,q_N) \in \mathcal{D} = (L \mathbb{T})^{2N}$. Two of these particles, with positions $q_1$ and $q_2$, form a dimer whereas the other $N-2$ particles, with positions $q_3,\cdots,q_N$, constitute the solvent. The potential energy of the system is composed of three parts:
\begin{equation*}
  \begin{aligned}
    V(q) &= v(|q_1-q_2|) + \sum_{i \in \{1,2\}} \sum_{j=3}^N v_{\rm sol} (|q_i-q_j|) + \sum_{3 \leqslant i<j \leqslant N} v_{\rm sol} (|q_i-q_j|) \\

    &=: V_{\rm dim}(q) + V_{\rm inter}(q) + V_{\rm sol}(q),
  \end{aligned}
\end{equation*}
where $V_{\rm dim}$ is the potential energy of the dimer, $V_{\rm inter}$ is the interaction energy between the dimer and the solvent and $V_{\rm sol}$ is the potential energy of the solvent. The two particles forming the dimer interact via a double-well potential: denoting by $r=|q_1-q_2|$ the bond length,
\begin{equation}
  \label{eq:double well}
  v(r) = h \left[ 1 - \left( \frac {r-r_0}{\Delta r} \right)^2 \right]^2,

\end{equation}
where $r_0$, $\Delta r > 0$ (see Figure~\ref{fig:potentials}, Left). The potential $v$ presents two minima: one associated with a compact state of length $r=r_0-\Delta r$ and one associated with a stretched state of length $r=r_0+\Delta r$. These minima are separated by a potential barrier of height $h$. The particles of the solvent interact both with the other particles of the solvent and the particles of the dimer through a purely repulsive potential. In the following we consider two types of potentials with compact support (see Figure~\ref{fig:potentials}, Right): a soft repulsion potential (used in~\cite{Groot97} for example)
\begin{equation}
\label{eq:soft}
\forall r > 0, \quad v_{\rm sol}(r) = \varepsilon \left( 1 - \frac r {r_{\rm cut}} \right) ^2  \mathds{1}_{r \leqslant r_{\rm cut}},
\end{equation}
where $\varepsilon, r_{\rm cut} > 0$; and a singular Coulomb-like potential:
\begin{equation}
  \label{eq:coulomb}

  \forall r > 0, \quad v_{\rm sol}(r) = \varepsilon \left( \frac {\frac 1 {\sqrt r} - \frac 1 {\sqrt {r_{\rm cut}}}} {\frac 1 {\sqrt \sigma} - \frac 1 {\sqrt {r_{\rm cut}}}} \right)^2 \mathds{1}_{\rm r < r_{\rm cut}} = \varepsilon \frac \sigma r \left( \frac{1 - \sqrt{\frac r {r_{\rm cut}}}} {1 - \sqrt{\frac \sigma {r_{\rm cut}}}} \right)^2 \mathds{1}_{r \leqslant r_{\rm cut}}.

\end{equation}
This potential behaves like $\frac 1 r$ for $r \to 0$, reaches the value $\varepsilon$ at $r=\sigma$ and vanishes at $r=r_{\rm cut}$,where its derivative also vanishes. Note that we recover the Coulomb potential $\varepsilon \frac \sigma r$ in the limit $r_{\rm cut} \to +\infty$.
\begin{figure}
\caption{Pairwise potentials for the solvated dimer model. Left: Potential for the dimer and associated free energy in vacuum (see~\eqref{eq:1D dimer poisson problem}). Right: Potentials for the solvent interaction.}
\begin{center}
\includegraphics[scale=.58]{dimer_potentials.pdf}
\includegraphics[scale=.58]{sol_potentials.pdf}
\end{center}
\label{fig:potentials}
\end{figure}

The system is driven out of equilibrium by a shearing force of amplitude $\nu$. More precisely, a particle located at a position $(q_{i,x}, q_{i,y})$ experiences the force~\cite{Joubaud12, Todd17}:
\begin{equation*}
	(0,f(q_{i,x}))= \left( 0,\nu \sin(2\pi\frac{q_{i,x}}L) \right).
\end{equation*}
This force is in the $y$ direction and depends only on $x$. It therefore induces a non-equilibrium forcing since it is not of gradient type. We are interested in computing the mean length of the dimer $R(q,p) = |q_1 - q_2|$ as a function of this external forcing. The corresponding average is denoted by $\mathbb{E}[|q_1 - q_2|]$.

For simplicity we study the overdamped dynamics associated with~$V$, but everything can be adapted to the Langevin case. Since the space $\mathcal{D}$ is compact and the noise in the dynamics is non degenerate, there exist a unique invariant probability measure $\pi$ by the Doeblin condition when the potentials under consideration are smooth. This invariant measure depends on~$\nu$, and is not explicit. Proving a similar result for singular potentials such as the Coulomb-like potential~\eqref{eq:coulomb} would require more work.

The generator can be decomposed as:
\begin{equation*}
\begin{aligned}
  \mathcal{L} &= -\nabla V(q) \cdot \nabla + \beta^{-1} \Delta + \nu \sum_{i=1}^N f(q_{i,x}) \partial_{q_{i,y}}
  = \calL_{\rm dim} + \calL_{\rm inter} + \calL_{\rm sol} + \nu \calL_\mathrm{pert},
\end{aligned}
\end{equation*}
where
\begin{alignat*}{3}
	\calL_{\rm dim} &= \sum_{i=1,2} \left( -\nabla_{q_i} V_{\rm dim}(q) \cdot \nabla_{q_i} + \beta^{-1} \Delta_{q_i} \right),& \qquad
	\calL_{\rm inter} &= - \nabla V_{\rm inter}(q) \cdot \nabla,\\
	\calL_{\rm sol} &= \sum_{i=3}^N \left( -\nabla_{q_i} V_{\rm sol}(q) \cdot \nabla_{q_i} + \beta^{-1} \Delta_{q_i} \right), &\qquad
	\calL_\mathrm{pert} &= \sum_{i=1}^N f(q_{i,x}) \partial_{q_{i,y}}.
\end{alignat*}
Note that $\calL_{\rm dim}$ is the generator of the dynamics of the dimer at equilibrium in vacuum and $\calL_{\rm sol}$ is the generator of the dynamics of the solvent at equilibrium and without dimer.

\subsection{Simplified dynamics and control variate}

We consider the following reference Poisson equation where the system is at equilibrium and the interaction between the dimer and the solvent has been switched off:
\begin{equation}
\label{eq:dimer ref poisson problem}
	-\calL_0 {\Phi_0} = R - \mathbb{E}_0[R],
\end{equation}
where
\begin{equation*}
	\calL_0 = \calL_{\rm dim} + \calL_{\rm sol}.
\end{equation*}
Let us show that this equation admits a solution ${\Phi_0}$ depending only on the length $|q_1-q_2|$ of the dimer. In order to highlight the dependence on the dimension of the underlying space, let us denote by $d=2$ this dimension. Assume that ${\Phi_0}$ is defined for any $q \in (L \mathbb{T})^{dN}$ by ${\Phi_0}(q) = \frac 1 2 \psi(|q_1-q_2|)$ for some smooth function $\psi$. The Laplacian of ${\Phi_0}$ can be rewritten using spherical coordinates as:
\begin{equation}
\label{eq:spherical laplacian}
	\Delta {\Phi_0}(q) = \psi''(|q_1-q_2|) + \frac {d-1}{|q_1-q_2|} \psi'(|q_1-q_2|),
\end{equation}
where $d=2$ is the dimension of the underlying physical space. We obtain by substituting ${\Phi_0}$ into~\eqref{eq:dimer ref poisson problem} that $\psi$ satisfies the following one-dimensional differential equation:
\begin{equation}
\begin{aligned}
\label{eq:1D dimer poisson problem}
	\forall r >0, \quad v_*'(r) \psi'(r) - \beta^{-1} \psi''(r) &= r - r^*,
\end{aligned}
\end{equation}
where $v_*(r) = v(r) - \frac {d-1} \beta \ln(r)$ and $r^* = \mathbb{E}_*[r]$ is the expectation of the length $r$ with respect to the probability measure $\pi_*(\dd r) = Z_*^{-1} \mathrm{e}^{-\beta v_*(r)} \, \dd r$. Note the additional term $- \frac {d-1} \beta \ln(r)$ in the expression of $v_*$ coming from~\eqref{eq:spherical laplacian}, which can be interpreted as an entropic contribution.

Let us first discuss the well-posedness of~\eqref{eq:1D dimer poisson problem}. The double-well potential~\eqref{eq:double well} considered here is such that $v_*$ is a bounded perturbation of a convex function. Therefore $\pi_*(\dd r)$ satisfies a log-Sobolev inequality and thus a Poincaré inequality by the Holley-Stroock theorem~\cite{Holley87} and the Bakry-Emery criterion~\cite{Bakry85}. This implies that the one-dimensional Poisson problem~\eqref{eq:1D dimer poisson problem} then admits a unique solution in
\[
\mathrm{H}^1(\pi_*) \cap \mathrm{L}_0^2(\pi_*) = \left\{ \varphi \in H^1(\pi_*), \ \int_0^{+\infty} \varphi \, \dd\pi_* = 0 \right\}
\]
by the Lax-Milgram theorem for the variational formulation:
\begin{equation*}
  \forall u \in \mathrm{H}^1(\pi_*) \cap \mathrm{L}_0^2(\pi_*), \quad \beta^{-1} \int_0^\infty \psi'(r) u'(r) \, \pi_*(\dd r) = \int_0^\infty (r-r_*) u(r) \, \pi_*(\dd r).
\end{equation*}
We discuss precisely in Appendix~\ref{ap:resolution ode} how we numerically solve~\eqref{eq:1D dimer poisson problem}. Knowing the solution $\psi$, the corresponding modified observable then writes
\begin{equation*}
\begin{aligned}
  (R + \mathcal{L} {\Phi_0})(q) &= |r_{12}| + \beta^{-1} \psi''(|r_{12}|) \\
  & + \left[ \frac 1 2 \big( \nabla_{q_1} V(q) - \nabla_{q_2} V(q) - \nu (f(q_{1,x}) - f(q_{2,x}))e_y \big) \cdot \frac {r_{12}}{|r_{12}|} +  \frac {d-1}{\beta |r_{12}|} \right] \psi'(|r_{12}|),
\end{aligned}
\end{equation*}
where $r_{12} = q_2 - q_1$ and $e_y=(0,1)$. Note that $\nabla_{q_1} V$ and $\nabla_{q_2} V$ are the forces that apply on particles $1$ and $2$ respectively, which depend also on the solvent variables.

\subsection{Numerical results}

We simulate a system of $N=64$ particles in $d=2$ dimensions, using periodic boundary conditions. We fix $L=8$ (so that the particle density is $1$), and $\beta = 1$. The parameters of the potentials are set to $r_{\rm cut}=2.5$, $\varepsilon=1$, $h=1$, $r_0=3$ and $\Delta r = 1$ (see Figure~\ref{fig:potentials}). For the finite difference method used to solve the Poisson equation~\eqref{eq:1D dimer poisson problem} we use a mesh size $\Delta r = 10^{-3}$ on an interval $[0,r_{\rm max}]$ with $r_{\rm max} = 10$ (see Appendix~\ref{ap:resolution ode}).

The influence of the shearing on the average dimer length is plotted on Figure~\ref{fig:dimer length}. We see that a shear force of amplitude $\nu=1$ increases the mean length by roughly $1\%$, and that the response of the mean length to the nonequilibrium forcing is of order $2$. The response is small thus difficult to estimate accurately, hence the need for control variates to alleviate this issue.

In the case of an unsolvated dimer, Figure~\ref{fig:dimer variance reduction} (Left) shows that the variance of the modified observable scales like $\nu^2$, as predicted by Theorem~\ref{th:var order 2 eta exact}. Note that in the limit $\nu \to 0$ the modified observable is the constant $\mathbb{E}_0[R] = r_*$, which is computed by a numerical quadrature, so that the variance converges to zero.

When the solvent interacts with the dimer the variance of the modified observable plateaus at a certain value when $\nu \to 0$ , as expected from Theorem~\ref{th:var order 2 eta}. For the soft potential~\eqref{eq:soft}, the variance scales like $\nu$ for a forcing amplitude of order $1$, which is expected from Theorem~\ref{th:var order 2 eta} (see Figure~\ref{fig:dimer variance reduction}, Right). The variance stabilizes at a value which is ten times smaller than the initial one. For the Coulomb-like potential the influence of the solvent on the dimer is stronger and the control variate does not perform as well, as seen on Figure~\ref{fig:dimer variance reduction} (Bottom). For a small shearing the variance is however reduced by a factor 4.

\begin{figure}
\caption{Left: Mean length of a dimer, either unsolvated (in vacuum) or in a solvent with soft or Coulomb-like potential. Right: Relative variation of this mean length induced by the shearing. The solid line represents the reference scaling $\nu^2$.}
\begin{center}
\includegraphics[scale=.58]{nice_length_wrt_nu.pdf}
\includegraphics[scale=.58]{relative_length_wrt_nu.pdf}
\end{center}
\label{fig:dimer length}
\end{figure}

\begin{figure}
\caption{Asymptotic variance of the length of the dimer, with or without control variate. Left: Unsolvated dimer. Right: Solvent with the soft potential~\eqref{eq:soft}. Bottom: Solvent with the Coulomb-like potential~\eqref{eq:coulomb}.}
\begin{center}
\includegraphics[scale=.58]{none_variance_wrt_nu.pdf}
\includegraphics[scale=.58]{softdpd_variance_wrt_nu.pdf}
\includegraphics[scale=.58]{coulomb_variance_wrt_nu.pdf}
\end{center}
\label{fig:dimer variance reduction}
\end{figure}

\paragraph{Generalization.}

The variance reduction strategy discussed here can be easily adapted to similar systems. For example Langevin dynamics can be treated by replacing~\eqref{eq:1D dimer poisson problem} by a two-dimensional PDE where the variables are the dimer length and the radial part of the momentum associated to this length. The Poisson equation should then be solved using a Galerkin approximation similar to what is done Section~\ref{s:sinus potential}. One could also consider a solvated molecule more complex than a dimer. In this case the PDE~\eqref{eq:1D dimer poisson problem} would be posed in several dimensions and thus become rapidly impossible to solve in practice. In general one has to reduce the system to a few relevant variables corresponding to a simplified Poisson equation in order to use the control variate approach developed here. This is connected to coarse-graining, \emph{i.e.} finding a few (nonlinear) functions of the degrees of freedom providing some macroscopic information on the system -- think of identifying an appropriate molecular backbone for proteins. An alternative route, which does not require an a priori physical knowledge of the system, would be to resort to greddy methods~\cite{Mokdad07, Temlyakov08, Cances11, Figueroa12}. If the system possesses in addition a specific symmetry or structure, one can make profit of dedicated tensor formats~\cite{Hackbusch12} as done for the Schrödinger equation in~\cite{Yserentant10}. This approach would be particularly adapted when studying an isotropic system composed of identical particles for example. Let us also mention recent advances on the resolution of Poisson equations based on deep convolutional neural networks~\cite{Shan18, Avrutskiy17}, which offer the promise of a better scalability with respect to the dimension of the system.

\section*{Acknowledgements}

The idea of using control variates came out of discussions with Antonietta Mira (USI) while Gabriel Stoltz was participating to the workshop "Free-energy calculations: a mathematical perspective" at Oaxaca. We thank Antoine Levitt (ENPC), Greg Pavliotis (Imperial College), Stefano Lepri (ISC) and Jonathan Goodman (NYU) for helpful discussions. This work is supported by the Agence Nationale de la Recherche under grant ANR-14-CE23-0012 (COSMOS); as well as the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) -- ERC Grant Agreement number 614492. We also benefited from the scientific environment of the Laboratoire International Associ\'e between the Centre National de la Recherche Scientifique and the University of Illinois at Urbana-Champaign. Part of this work was done during the authors' stay at the Institut Henri Poincar\'e - Centre Emile Borel during the trimester "Stochastic Dynamics Out of Equilibrium" (April-July 2017). The authors warmly thank this institution for its hospitality.

redpajama_set_name: RedPajamaArXiv


text:
Q: Can't prevent page reload after form submit/file upload i want to upload photo with ajax to Express, but every time i submit a form it doesn't generate html that i have in code, it just redirects me and shows picture data as an object. is it even possible to upload a picture and if it succeeds, generate a desired response as html?
here's the Express code:
const storage = multer.diskStorage({
   destination: function (req, file, cb) {
      cb(null, 'public/uploads/basic_uploads/')
   },
   filename: function (req, file, cb) {
      cb(null, Date.now() + '-' + file.originalname)
   }
});
const upload = multer({ storage: storage })

router.post('/ajax_upload', upload.single('ajaxFile'), (req, res) => {
   let data = req.file;
   res.status(200).json(data)
});

everything is ok with uploading a picture.
and then (Ajax):
$(document).ready(function (event){
      $("#uploadForm").on('submit',( function (event) {
         event.preventDefault();
         $.ajax({
            url: "/ajax_upload",
            type: "POST",
            data:  new FormData(this),
            contentType: false,
            cache: false,
            processData: false,
            success: successHandler,
            error: function(error){
               $("#targetLayer").html(`<span>` + error + `</span>`)
            }
         });
      }));

      function successHandler(data) {
         // let html = `<img src="/uploads/basic_uploads/` + data.filename + `" >`
         // $("#targetLayer").html(html)
         $("#targetLayer").html(`<p>AT LAAAAAST</p>`)

      }
});

and here is html:

                <form id="uploadForm" action="/ajax_upload" method="POST" enctype="multipart/form-data" onsubmit="submitFormAjax();return false;">

                        Ajax file Upload with Multer
                        <input type="file" name="ajaxFile" id="exampleInputFile">

                        <input class="btn btn-primary" type="submit" value="Basic file Upload">
                </form>

A: after successHandler function add return false statement and try again,
$(document).ready(function (event){
  $("#uploadForm").on('submit',( function (event) {
     event.preventDefault();
     $.ajax({
        url: "/ajax_upload",
        type: "POST",
        data:  new FormData(this),
        contentType: false,
        cache: false,
        processData: false,
        success: successHandler,
        error: function(error){
           $("#targetLayer").html(`<span>` + error + `</span>`)
        }
     });
  }));

  function successHandler(data) {
     // let html = `<img src="/uploads/basic_uploads/` + data.filename + `" >`
     // $("#targetLayer").html(html)
     $("#targetLayer").html(`<p>AT LAAAAAST</p>`)

  }

 return false;
});

redpajama_set_name: RedPajamaStackExchange


text:
1953 Kurtis-Kraft 500S Murphy Special
One of a handful of authentic 500 S sports cars known to exist. Documented race history including at Pebble Beach and Sebring 12 hours. A rare opportunity to acquire one of America's finest home-grown sports racing cars.
Among the approximately 25 Kurtis Kraft 500 S roadsters built, chassis number 015, also known as The Murphy Special, is one of the most fascinating and unique. This car, presented here with a Pebble-Beach class-winning restoration, has an impressive record of race wins and podium finishes that came primarily with its first owner, Mr. Bill Murphy. Based in Culver City, California, Bill Murphy was a life-long car guy who ran one of the most successful Buick dealerships in America. He was already a successful race car owner, running a series of Kurtis-Kraft Midgets for local hero Sam Hanks. However, this new Kurtis would be Bill's personal steed to use in the rapidly growing West Coast sports car racing scene.
In 1953, Bill Murphy ordered his 500S directly from Kurtis-Kraft in Los Angeles. Considering Murphy already owned several Kurtis cars, it is not surprising to learn the 500S he had built was no off-the-shelf effort. It is believed that Murphy bought a used Indy Car to obtain the trick quick-change rear-end and aircraft-inspired disc brakes made by the legendary Conze Brothers Machine Shop in California. Of course, Buick power was the only way to go for Murphy, who supplied a brand new 322 cubic-inch Nailhead V8, which he later claimed made upwards of 300 horsepower on fuel injection and was good for an estimated 180 miles per hour.

redpajama_set_name: RedPajamaCommonCrawl


text:
Home News Op-Ed: Newport Democratic City Committee expresses support for the leadership of Congressman...
Op-Ed: Newport Democratic City Committee expresses support for the leadership of Congressman David Cicilline following the attack on the Capitol
In response to the events which unfolded last week, the Newport Democratic City Committee expresses its support for the leadership of our Congressman, the Honorable David N. Cicilline, to hold President Trump accountable for his role in inciting the mob attack on our nation's Capitol and the legitimate Constitutional process underway there on January 6th.
The images were shocking to those who respect our Constitutional democracy and the rule of law. The desecration of the house where our duly elected representatives meet to deliberate and legislate was saddening, heartbreaking, and the cause of legitimate outrage. The loss of life was tragic.
The President, and each person involved with this seditious attack on our democracy, must be held responsible. The President must face consequences for his role in encouraging the extreme far-right agitators who attacked last Wednesday, as well as for his speech that day and his subsequent expression of support for these domestic terrorists and their actions.
Congressman Cicilline has taken the lead since Wednesday, first in calling for the invocation of the procedure defined in the 25th Amendment, and today, with no response from the Republican administration, in introducing articles of impeachment. Congressman Cicilline recognizes the importance of our democratic form of government and the attendant need to take decisive action against this lawless and dangerous President. We are proud of you, Congressman. You have our support.
Newport Democratic City Committee
Susan D. Taylor, Chair
J. Clement Cicilline, outgoing Chair
Previous articleRISCA to commission public art for Garrahy Judicial Complex
Next articlePortsmouth Free Public Library to present 'Slavery in Colonial Rhode Island' with Fred Zilian
This content has been contributed to What's Up Newp. The views and opinions included within are not necessarily those of What's Up Newp, our contributors, or our advertisers. We welcome letters to the editor on current local topics. Email them to Ryan@whatsupnewp.com.
Aquidneck Island Climate Caucus to host a virtual conversation on climate change with Senator Whitehouse
COVID-19 in Rhode Island: The latest data
Governor Raimondo to deliver State of the State Address on February 3
LIVE ON WUN: Newport City Council COVID-19 Status Update (Jan. 21 at 4:30 pm)
What's Up in Newport County today: Tuesday, January 19