\documentclass[a4paper,11pt]{amsart}

\pdfoutput=1

\usepackage[text={400pt,660pt},centering]{geometry}

%\usepackage[utf8]{inputenc}
%\usepackage[T1]{fontenc}
%\usepackage{lmodern}
%\usepackage[french,english]{babel}
\usepackage{amsthm, amssymb, amsmath, amsfonts, mathrsfs}
\usepackage{mathtools}
\usepackage{scalerel} 
\usepackage{stackrel}
\usepackage[colorlinks=true, pdfstartview=FitV, linkcolor=blue, citecolor=blue, urlcolor=blue,pagebackref=false]{hyperref}
\usepackage{esint} % contains strokedint
\usepackage{MnSymbol} % contains notin inverted (which is \nowns)
\usepackage{booktabs} % nice tables
%\usepackage{enumitem} %hyperlinks to item in lists

%\usepackage{amscd, graphicx, epsfig, psfrag, enumerate, dsfont}

%\usepackage{color}

\usepackage{microtype}
\usepackage{cleveref}

\parskip= 2pt

\setcounter{tocdepth}{2}


%\usepackage[notcite,notref,color]{showkeys}
%\definecolor{labelkey}{gray}{.8}
%\definecolor{refkey}{gray}{.8}

%\definecolor{darkgreen}{rgb}{0,0.5,0}
%\definecolor{darkblue}{rgb}{0,0,0.7}
%\definecolor{darkred}{rgb}{0.9,0.1,0.1}

%\newcommand{\jccomment}[1]{\marginpar{\raggedright\scriptsize{\textcolor{darkgreen}{#1}}}}

%\newcommand{\gucomment}[1]{\marginpar{\raggedright\scriptsize{\textcolor{darkred}{#1}}}}

%\newcommand{\mncomment}[1]{\marginpar{\raggedright\scriptsize{\textcolor{magenta}{#1}}}}


%\newcommand{\jc}[1]{{\color{darkgreen}{#1}}}

%\newcommand{\gu}[1]{{\color{darkred}{#1}}}

%\newcommand{\mn}[1]{{\color{magenta}{#1}}}


\newtheorem{proposition}{Proposition}
\newtheorem{theorem}[proposition]{Theorem}
\newtheorem{lemma}[proposition]{Lemma}
\newtheorem{corollary}[proposition]{Corollary}
\newtheorem{claim}[proposition]{Claim}
\theoremstyle{remark}
\newtheorem{remark}[proposition]{Remark}
\newtheorem{example}[proposition]{Example}
\theoremstyle{definition}
\newtheorem{definition}[proposition]{Definition}
\newtheorem{conjecture}[proposition]{Conjecture}
\newtheorem{question}[proposition]{Question}
\newtheorem{assumption}[proposition]{Assumption}
\newtheorem{hypothesis}[proposition]{Hypothesis}
\newcommand{\blue}{\color{blue}}
\numberwithin{equation}{section}
\numberwithin{proposition}{section}
\numberwithin{table}{section}



\renewcommand{\le}{\leqslant}
\renewcommand{\ge}{\geqslant}
\renewcommand{\leq}{\leqslant}
\renewcommand{\geq}{\geqslant}
\newcommand{\ls}{\lesssim}
\newcommand{\les}{\ls}
\renewcommand{\subset}{\subseteq}
\renewcommand{\supset}{\supseteq}
\newcommand{\mcl}{\mathcal}
\newcommand{\msf}{\mathsf}
\newcommand{\msc}{\mathscr}
\newcommand{\A}{{\mathsf{A}}}
\newcommand{\G}{\mathcal{G}}
\renewcommand{\S}{\mathsf{S}}
\newcommand{\M}{\mathbf{M}}
\newcommand{\E}{\mathbb{E}}
\newcommand{\Er}{\mathbb{E}_{\rho}}
\renewcommand{\Pr}{\mathbb{P}_{\rho}}
\renewcommand{\a}{\mathbf{a}}
\newcommand{\aone}{\mathbf{a}^{\{1\}}}
\newcommand{\ab}{{\overbracket[1pt][-1pt]{\a}}}
\newcommand{\ahom}{{\overbracket[1pt][-1pt]{\a}}}
\newcommand{\ub}{{\overbracket[1pt][-1pt]{u}}}
\newcommand{\Ub}{{\overbracket[1pt][-1pt]{U}}}
\newcommand{\varb}{{\overbracket[1pt][-1pt]{\varphi}}}
\newcommand{\ut}{\tilde{u}}
\newcommand{\na}{\nabla}
\newcommand{\Ll}{\left}
\newcommand{\Rr}{\right}
\newcommand{\lhs}{left-hand side}
\newcommand{\rhs}{right-hand side}
\newcommand{\1}{\mathbf{1}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\Z}{\mathcal{Z}}
\newcommand{\Zd}{{\mathbb{Z}^d}}
\renewcommand{\P}{\mathbb{P}}
\newcommand{\ov}{\overline}
\renewcommand{\bar}{\overline}
\newcommand{\un}{\underline}
\newcommand{\td}{\widetilde}
\renewcommand{\tilde}{\widetilde}
\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}
\newcommand{\eps}{\varepsilon}
\newcommand{\ep}{\varepsilon}
\renewcommand{\epsilon}{\varepsilon}
\renewcommand{\d}{{\mathrm{d}}}
\newcommand{\var}{\mathbb{V}\!\mathrm{ar}}
\newcommand{\cov}{\mathbb{C}\mathrm{ov}}
\newcommand{\mfk}{\mathfrak}
\newcommand{\dr}{\partial}
\newcommand{\X}{\mathbb{X}}
\newcommand{\PP}{\mathbf{P}}
\newcommand{\EE}{\mathbf{E}}
\renewcommand{\O}{{\mathcal{O}}}
\newcommand{\Pb}{\bar{P}}
\newcommand{\Ee}{\mathcal{E}}
\newcommand{\D}{\mathcal{D}}
\newcommand{\la}{\left\langle}
\newcommand{\ra}{\right\rangle}
\newcommand{\rah}{\right\rangle_{\cH^{-1}, \cH^1_0}}


\newcommand{\cu}{{\scaleobj{1.2}{\square}}}
%\newcommand{\cud}{{\scaleobj{1.2}{\boxdot}}}

%\newcommand{\la}{\left\langle}
%\newcommand{\ra}{\right\rangle}
\newcommand{\vb}{\ \big \vert \ }

\renewcommand{\fint}{\strokedint}
\newcommand{\Rd}{{\mathbb{R}^d}}
\newcommand{\n}{\mathbf{n}}
\newcommand{\mmd}{\mathcal{M}_\delta}
\newcommand{\bmmd}{\bar{\mathcal{M}}_\delta}
%\newcommand{\mmb}{\mathcal{M}_\bullet}
%\newcommand{\ovN}{{\overrightarrow N}}
%\newcommand{\EUN}{\mathbb{E}_{U,N}}
%\newcommand{\EUNmu}{\mathbb{E}_{U,N,\mu_0}}
%\newcommand{\mt}{\mathscr{M}}
\newcommand{\fil}{\mathscr{F}}


%\DeclareMathOperator{\ins}{\msf{int}}
\DeclareMathOperator{\dist}{dist}
\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\diam}{diam}

%%% 
\newcommand{\cD}{\msc{D}}   % Smooth function in configuration space
\newcommand{\cC}{\msc{C}}   % Smooth function in configuration space
\newcommand{\cfC}{\mcl{F}\msc{C}} % Smooth, and finitely generated function in configuration space
\newcommand{\cL}{\msc{L}}   % L2 integrable for Poisson measure
\newcommand{\cH}{\msc{H}}   % Sobolev space in configuration space
\newcommand{\cB}{\dot{\mathbb{B}}}   % Besove space in configuration space
\newcommand{\cBs}{\cB^1_{2,\infty}}
\newcommand{\acL}{\underline{\msc{L}}}
\newcommand{\acH}{\underline{\msc{H}}}
\newcommand{\aL}{\underline{L}}
\newcommand{\aH}{\underline{H}}
\newcommand{\poi}{\text{Poi}}
\newcommand{\bin}{\text{Bin}}
\newcommand{\Ind}[1]{\mathbf{1}_{\left\{#1\right\}}}
\newcommand{\id}{\mathsf{Id}}
\newcommand{\norm}[1]{\left\Vert{#1}\right\Vert}
\newcommand{\bracket}[1]{\left\langle{#1}\right\rangle}
\newcommand{\dbint}[1]{[\![#1]\!]}
\newcommand{\xtag}{x_{\bullet}}
\renewcommand{\n}{\overrightarrow{\mathbf{n}}}
\newcommand{\op}[1]{#1^\mathrm{o}}
\newcommand{\inte}{\operatorname{int}}
\newcommand{\bmu}{\boldsymbol{\mu}}
\newcommand{\EEbmu}{\EE_{\bmu_0}}


%%%
%\newcommand{\red}{\color{red}}

\newcommand{\mres}{\mathbin{\vrule height 1.4ex depth 0pt width
		0.13ex\vrule height 0.13ex depth 0pt width 1ex}}

\newcommand{\mressmall}{\mathbin{\vrule height 1ex depth 0pt width
		0.1ex\vrule height 0.1ex depth 0pt width 0.7ex}}

% The averaged sum notation

%\newcommand{\avsum}{\mathop{\mathpalette\avsuminner\relax}\displaylimits}

%\makeatletter
%\newcommand\avsuminner[2]{%
	%  {\sbox0{$\m@th#1\sum$}%
		%   \vphantom{\usebox0}%
		%   \ooalign{%
			%     \hidewidth
			%     \smash{\rule[.5ex]{1.5ex}{.8pt} \relax}%
			%     \hidewidth\cr
			%     $\m@th#1\sum$\cr
			%   }%
		%  }%
	%}
%\makeatother




\title[Quantitative fluctuations for interacting particle systems]{Quantitative equilibrium fluctuations for interacting particle systems}
\author[C.\ Gu, 
J.-C.\ Mourrat, M.\ Nitzschner]{Chenlin Gu, Jean-Christophe Mourrat, Maximilian Nitzschner}


\address[Chenlin Gu]{Yau Mathematical Sciences Center, Tsinghua University,   China}

\address[Jean-Christophe Mourrat]{Ecole Normale Sup\'erieure de Lyon and CNRS, Lyon, France}

\address[Maximilian Nitzschner]{Department of Mathematics, The Hong Kong University of Science and Technology, Hong Kong}




\begin{document}
	%
	
	\begin{abstract}
		
		We consider a class of interacting particle systems in continuous space of non-gradient type, which are reversible with respect to Poisson point processes with constant density. For these models, a rate of convergence was recently obtained in \cite{bulk} for certain finite-volume approximations of the bulk diffusion matrix. Here, we show how to leverage this to obtain quantitative versions of a number of results capturing the large-scale fluctuations of these systems, such as the convergence of two-point correlation functions and the Green--Kubo formula. 
		
		
				\bigskip
		
				\noindent \textsc{MSC 2010:} 82C22, 35B27, 60K35.
		
				\medskip
		
				\noindent \textsc{Keywords:} interacting particle system, quantitative homogenization, equilibrium fluctuation, Green--Kubo formula, density field correlation, current-current correlation.
		
	\end{abstract}
	%
	\maketitle
	
	%\tableofcontents
	
	%
	%
	%
	%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	%
	%
	%
	\section{Introduction}
	We consider a system of particles $(X_i(t))_{i \in \N, t \ge 0}$ evolving on $\R^d$ with local interactions. We take particles to be indistinguishable, so we encode their positions by keeping track of the measure 
	\begin{equation}
		\label{e.def.mut}
		\bmu_t := \sum_{i \in \N} \delta_{X_i(t)},
	\end{equation}
	an element of the configuration space $\mcl M_\delta(\Rd)$ of locally finite measures that are sums of Dirac masses on $\R^d$. We denote by $\P_\rho$ the law of the Poisson point process with constant intensity $\rho > 0$, which is a probability measure on $\mcl M_\delta(\Rd)$, and we write $\E_\rho$ for the associated expectation. The model we consider is reversible with respect to $\P_\rho$ for every $\rho$, and is of non-gradient type. Precisely, its definition is based on the choice of a measurable function $\a_\circ$ from $\mcl M_\delta(\Rd)$ to the space $\R^{d \times d}_{\mathrm{sym}}$ of $d$-by-$d$ symmetric matrices that satisfies the following properties.
	\begin{itemize}
		\item \emph{Uniform ellipticity}: there exists $\Lambda < +\infty$ such that for every $\mu \in \mcl M_\delta(\Rd)$ and $\xi \in \Rd$,
		\begin{equation}
			\label{e.ass.unif.ell}
			|\xi|^2 \le \xi \cdot \a_\circ(\mu) \xi \le \Lambda |\xi|^2.
		\end{equation}
		
		
		\item \emph{Finite range of dependence}: for every $\mu \in \mcl M_\delta(\Rd)$, we have 
		\begin{equation}
			\label{e.finite.range}
			\a_\circ(\mu) = \a_\circ(\mu \mres B_1) .
		\end{equation}
	\end{itemize}
	In \eqref{e.finite.range} and throughout the paper, for every $r > 0$, we write $B_r$ to denote the open Euclidean ball of radius $r \ge 0$ centered at the origin, and for every measure $\mu$ on $\Rd$ and Borel set $A \subset \Rd$, we write $\mu \mres A$ to denote the restriction of the measure $\mu$ to $A$, that is, for every Borel set $B \subset \Rd$,
	\begin{equation*}  %\label{e.}
		(\mu \mres A)(B) = \mu(A \cap B). 
	\end{equation*}
	For every $x \in \Rd$, we write $\tau_x$ to denote the translation operator on measures, that is, for every Borel set $B \subset \Rd$,
	\begin{equation*}  %\label{e.}
		(\tau_{-x} \mu)(B) := \mu(x + B),
	\end{equation*}
	and we set, for every $\mu \in \mcl M_\de(\Rd)$ and $x \in \Rd$,
	\begin{equation}
		\label{e.def.amu}
		\a(\mu,x) := \a_\circ(\tau_{-x} \mu). 
	\end{equation}
	Roughly speaking, our model is built in such a way that each particle $X_i(t)$, $i \in \N$, follows a diffusive motion whose diffusion matrix at time $t$ is given by 
	\begin{equation*}  %\label{e.}
		\a(\bmu_t, X_i(t)) = \a_\circ\big(\tau_{-X_i(t)} \bmu_t\big) = \a_\circ \bigg( \sum_{j \in \N} \delta_{X_j(t) - X_i(t)}\bigg) .
	\end{equation*}
	As will be explained below, one convenient way to give a precise definition of this stochastic process and its associated semigroup $(P_t)_{t \ge 0}$ is to relate these to the Dirichlet form 
	\begin{equation}  
		\label{e.intro.dirichlet.form}
		f \mapsto \E_\rho \Ll[ \int_\Rd \frac 1 2 \nabla f \cdot \a \nabla f \, \d \mu \Rr] ;
	\end{equation}
	we refer to~\eqref{e.def.deriv} for a definition of the gradient $\nabla$, and to Section~\ref{subsec.Dirichletform} for a detailed construction. We denote by $\PP_{\bmu_0}$ the law of the cloud of particles $(\bmu_t)_{t \ge 0}$ started from $\bmu_0$, and by $\EE_{\bmu_0}$ the associated expectation. A classical object of study is the rescaled density fluctuation field
	\begin{equation}  
		\label{e.def.YtL}
		Y_t^N := N^{-\frac d 2}(\bmu_{N^2 t} - \rho m)(N\, \cdot),
	\end{equation}
	where $N \ge 1$ and $m$ denotes the Lebesgue measure on $\Rd$. We postulate that the initial configuration $\bmu_0$ is sampled according to $\P_\rho$; since we denote by~$\mu$ the canonical random variable under $\E_\rho$, we implicitly set $\mu = \bmu_0$ in expressions such that $\Er \EEbmu$. One can think of $Y_t^N$ as a (random) distribution over $\Rd$, and testing $Y_t^N$ against a smooth function $f \in C^\infty_c(\Rd;\R)$ of compact support results in
	\begin{equation}  
		\label{e.explicit.def.YtL}
		Y_t^N(f) = N^{-\frac d 2} \Ll( \int_\Rd f(N^{-1} x) \, \d \bmu_{N^2 t}(x) - \rho \int_\Rd f(N^{-1} x) \, \d x  \Rr) . 
	\end{equation}
	Using the invariance of the measure $\P_\rho$ under the dynamics, one can easily check that for each fixed $t \ge 0$, the distribution $Y_t^N$ converges in law under $\Pr \PP_{\bmu_0}$ to a white noise over~$\Rd$ with variance $\rho$ as $N$ tends to infinity. The space-time correlations are less direct and are governed by the bulk diffusion matrix $\ab \in \R^{d \times d}_{\mathrm{sym}}$, whose precise definition is given in~\eqref{eq:AbarDef}. For a discrete class of models analogous to our present setting, it is shown in \cite{cha96, fun96, lu94} that, as a space-time distribution, $(Y^N_t)_{t \ge 0}$ converges in law to~$(Y^\infty_t)_{t \ge 0}$ solution to
	\begin{equation}
		\label{e.eq.Yinfty}
		\dr_t Y^\infty_t = \frac 1 2 \nabla \cdot \ab \nabla Y^\infty_t + \nabla \cdot \Ll( \sqrt{\rho \ab} \xi_t \Rr) ,
	\end{equation}
	where $(\xi_t(x))_{t \ge 0, x \in \Rd}$ is a $d$-dimensional space-time white noise over $[0,+\infty) \times \Rd$. In particular, we expect the convergence of two-point density correlation functions, in the sense that for every $f, g \in C^\infty_c(\Rd;\R)$ and $t > s > 0$,
	\begin{equation}
		\label{e.intro.conv.two.point}
		\lim_{N \to + \infty} \E_\rho\EE_{\bmu_0} [Y^N_t(f) Y^N_s(g)]= \E[Y^\infty_t(f) Y^\infty_s(g)].
	\end{equation}
	The latter quantity can be rewritten explicitly as
	\begin{equation}
		\label{e.def.two.point.Yinfty}
		\E[Y^\infty_t(f) Y^\infty_s(g)] = \rho \int_{\Rd \times \Rd} f(x)\Psi_{t-s}(x-y)g(y)\, \d x  \d y,    
	\end{equation}
	where $(\Psi_t(x))_{t \ge 0, x \in \Rd}$ is the heat kernel for the operator $\partial_t - \frac 1 2 \nabla \cdot \ab \nabla$, that is,
	\begin{equation}
		\label{e.def.heat.kernel}
		\Psi_t(x) := \frac{1}{\sqrt{(2\pi t)^{{d}} \det \ab}}\exp\Ll(-\frac{x \cdot \ab^{-1} x}{2 t}\Rr).
	\end{equation}
	Our first main result gives a rate of convergence for \eqref{e.intro.conv.two.point}. We prefer to state it using the unrescaled centered field defined by
	\begin{equation}
		\label{e.def.Y}
		Y_t := Y^{1}_t = \bmu_t - \rho m,
	\end{equation}
	and for every $f \in L^1(\Rd) \cap L^2(\Rd)$, we write
	\begin{equation}
		\label{e.def.Y.f}
		Y_t(f) = \int_\Rd f \, \d \bmu_t - \rho \int_\Rd f \, \d m.
	\end{equation}
	\begin{theorem}[Quantitative asymptotics for two-point functions]
		\label{t.two.point}
		There exists an exponent $\beta(d,\Lambda,\rho) > 0$ and a constant $C(d,\Lambda,\rho) < + \infty$ such that for every $f, g \in L^1(\Rd) \cap L^2(\Rd)$ and $t > s > 0$, we have
		\begin{multline}
			\label{e.two.point}
			\Ll|\E_\rho \EEbmu \Ll[ Y_t(f) Y_s(g) \Rr] -\rho \int_{\Rd \times \Rd} f(x)\Psi_{t-s}(y-x)g(y)\, \d x  \d y \Rr|
			\\
			\le C (t-s)^{-\beta} \|f\|_{L^2} \, \|g\|_{L^2}. 
		\end{multline}	
	\end{theorem}
	By a simple scaling, Theorem~\ref{t.two.point} yields that for every $f, g \in L^1(\Rd) \cap L^2(\Rd)$, we have 
	\begin{equation}\label{e.two.pointScale}  %\label{e.}
		\Ll|  \E_\rho\EEbmu [Y^N_t(f) Y^N_s(g)]- \E [Y^\infty_t(f) Y^\infty_s(g)] \Rr| \le C (N^2(t-s))^{-\beta} \|f\|_{L^2} \, \|g\|_{L^2}. 
	\end{equation}
	
	
	As another illustration of our results, we show a quantitative version of the Green--Kubo formula. Standard versions of the Green--Kubo formula involve infinite-time and infinite-volume limits. Here we localize the time and space scales, with an estimate of the error. For convenience, we localize in space using boxes of the form
	\begin{equation}\label{eq.defcube}  %\label{e.}
		\cu_m := \Ll( -\frac{3^{m}}{2}, \frac{3^{m}}{2} \Rr) ^d, \quad m \in \N.
	\end{equation}
	The Green--Kubo formula involves an integral of the correlations of the current of the particle density, as for instance in \cite[Proposition~II.2.1]{spohn2012large}. Informally, for every $\mu \in \mmd(\Rd)$, the $\EE_\mu$-averaged instantaneous current at a point $x \in \supp \mu$ in the direction of $p \in \R^d$ should be given by
	\begin{equation*}  %\label{e.}
		\frac{1}{2}\nabla_x \cdot(\a(\mu,x)p) = \frac{1}{2}\nabla_x \cdot (\a(\mu' + \de_x, x) p), \qquad \text{ with } \mu' := \mu - \de_x.
	\end{equation*}
	This is a continuous analogue of the quantity 
	$j_l(\eta)$ in \cite[(II.2.23)]{spohn2012large}, with the correspondence of $(l,\eta,0)$ in~\cite{spohn2012large} to $(p,\mu,x)$ in our notation. Since we would like to state a version of the Green--Kubo formula in finite volume, it would be tempting to define the total current over the box $\cu_m$ as
	\begin{equation*}  %\label{e.}
		\int_{\cu_m} \frac{1}{2}\nabla_x \cdot (\a(\mu,x) p) \, \d \mu(x).
	\end{equation*}
	However, even if we leave aside the possible difficulties associated with a lack of smoothness of the mapping $x \mapsto \a(\mu' + \de_x, x)$, this definition fails to capture an important boundary contribution. Instead, we define the total current over the box $\cu_m$ by duality as follows. We let $\cL^2(\cu_m)$ be the space of square-integrable functions which are measurable with respect to the configuration in $\cu_m$; see also \eqref{eq.defL2}. We denote by $\cH^1_0(\cu_m)$ the space of functions with square-integrable gradients that can be approximated by smooth functions that only depend on $\mu \mres K$ for some compact set $K \subset \cu_m$; see below~\eqref{eq.defNormalH1} for a precise definition. We denote by $\cH^{-1}(\cu_m)$ the space dual to $\cH^{1}_0(\cu_m)$, with the canonical embedding
	\begin{equation*}  %\label{e.}
		\cH^{1}_0(\cu_m) \subset \cL^2(\cu_m) \subset \cH^{-1}(\cu_m),
	\end{equation*}
	and by $\la \cdot , \cdot \ra_{\cH^{-1}(\cu_m), \cH^1_0(\cu_m)}$ the duality pairing between $\cH^{-1}(\cu_m)$ and $\cH^1_0(\cu_m)$. 
	For every $p \in \Rd$ and $m \in \N$, we define $F_{p,m} \in \cH^{-1}(\cu_m)$, the integrated current over~$\cu_m$, by duality, so that for every $g \in \cH^1_0(\cu_m)$, we have
	\begin{equation}\label{eq.defFpm}  %\label{e.}
		\la F_{p,m} , g \ra_{\cH^{-1}(\cu_m), \cH^1_0(\cu_m)} := \Er \Ll[\int_{\cu_m} - \frac{1}{2} p \cdot \a \nabla g \, \d \mu\Rr].
	\end{equation}
	As is the case for diffusions on bounded subsets of $\Rd$, we need to choose a boundary condition in order to define the particle dynamics in finite volume, for instance of Dirichlet or Neumann type. In the context of particle systems, the ``Dirichlet''-type boundary condition is understood in the sense provided by the space $\cH^1_0(\cu_m)$. Heuristically, with this boundary condition, particles that hit the boundary instantaneously disappear, and some particles also appear spontaneously at the boundary to maintain detailed balance. Since this in effect places the finite-volume dynamics ``within a bath of particles at equilibrium'', we feel that this is a more appropriate dynamics for the estimation of the total flux than a Neumann-type boundary condition (which would correspond to dynamics with a conserved number of particles that reflect at the boundary). We denote by $(\bmu_t^{(m)})_{t \ge 0}$ the particle dynamics restricted to the box $\cu_m$ with ``Dirichlet'' boundary condition, and by $(P_t^{(m)})_{t \ge 0}$ the associated semigroup, so that $P_t^{(m)}(f)(\bmu_0) = \EEbmu [f(\bmu_t^{(m)})]$; these objects are defined more precisely in
Section~\ref{sec:5}. In analogy with \cite[(II.2.27)]{spohn2012large} (which we think should display a term $x_\alpha x_\beta$ in place of $\delta_{\alpha \beta} x_\alpha^2$), \cite[(1)]{sasada2018green}, or \cite[(2.2), (2.4), (8.7), (8.9)]{mourrat2019efficient}, our qualitative Green--Kubo formula takes the form 

	\begin{multline}
		\label{e.gk.quali}
		\frac{1}{2} p \cdot \ab p = \frac{1}{2}\Er[p\cdot \a_\circ(\mu + \delta_0) p] \\
		- \lim_{\lambda \to 0} \lim_{m \to \infty} \frac{1}{\rho |\cu_m|} \int_0^{+\infty}  e^{-\lambda t}  \la F_{p,m}, P_t^{(m)}(F_{p,m}) \rah \, \d t.
	\end{multline}
	%The factor of $\frac 1 2$ in front of the integral sign comes from the fact that we defined the semigroup with respect to the operator $\frac 1 2 \nabla \cdot \a \nabla$ and not $\nabla \cdot \a \nabla$. 
	Here we write $\la \cdot , \cdot \ra_{\cH^{-1}, \cH^1_0}$ as a shorthand notation for $\la \cdot , \cdot \ra_{\cH^{-1}(\cu_m), \cH^1_0(\cu_m)}$. Our second main result provides quantitative estimates that allow us to localize this formula in space and time. In particular, we can take $\lambda = 0$ there, in which case the integral of the current correlation becomes the classical case without any localization in time.
	\begin{theorem}[Quantitative Green--Kubo formula]
		\label{t.quanti.gk}
		There exist $\alpha(d,\Lambda,\rho) > 0$ and $C(d, \Lambda, \rho) < + \infty$ such that for every $m \in \N$, $p \in B_1$ and $\lambda \in [0, +\infty)$, we have 
		\begin{multline}\label{eq.quanti.gk}
			\Ll|\frac{1}{2} p \cdot \ab p - \frac{1}{2}\Er[p\cdot \a_\circ(\mu + \delta_0) p] + \frac{1}{\rho |\cu_m|} \int_0^{+\infty}  e^{-\lambda t}  \la F_{p,m}, P_t^{(m)}(F_{p,m}) \rah \, \d t \Rr|
			\\
			\leq C \Ll\{
			\begin{array}{ll}
				1 \qquad & \text{if } \lambda \in [1, +\infty),\\
				\lambda^{\frac{\alpha}{2(1+\alpha)}} \qquad & \text{if } \lambda \in (3^{-2(1+\alpha)m}, 1),\\
				3^{-\alpha m}&  \text{if } \lambda \in[0, 3^{-2(1+\alpha)m}].
			\end{array} \Rr. 
		\end{multline}
	\end{theorem}
	
	
	
	\subsection*{Related works} We now give a brief review of related works. We already mentioned \cite{cha96, fun96, lu94} which identified the scaling limit of the rescaled fluctuation field $Y^N$ for discrete models of non-gradient type similar to our continuous model. Similar results were also obtained in \cite{broxrost,cha94,chayau92, mpsw86, spo86} for other models. The hydrodynamic limit of some non-gradient models was obtained in \cite{klo94, quastel, varadhanII} using the entropy method introduced in \cite{gpv}, and in \cite{fuy} using the relative entropy method introduced in~\cite{yau}. Reference books on the topic include \cite{kipnis1998scaling, komorowski2012fluctuations, spohn2012large}. 
	
	We are not aware of any result addressing the out-of-equilibrium fluctuations of non-gradient models. While we focus here on the derivation of quantitative estimates for the fluctuations at equilibrium, we hope that our techniques will also be useful for at least some out-of-equilibrium situations. In relation to this, the fact  that the bulk diffusion matrix $\ab$ depends smoothly on the density $\rho$ has been shown in \cite{giunti2021smoothness}. For gradient models and small perturbations thereof, out-of-equilibrium fluctuation results have been derived in \cite{chayau92, mfl86, fpv88, jarmen18, prespo83, ravishankar1992fluctuations}.
	Recent progress on higher-order approximations and large deviations include \cite{cornalba2023dean, dirr2020conservative, fehrman2023nonequilibrium}. 
	
	The results of the present paper are based on the quantitative estimates obtained in \cite{bulk}. The approach taken up there is inspired by prior work on the homogenization of elliptic equations, as reviewed in \cite{armstrong2022elliptic, AKMbook, informal}. 
	
	Closely related to the problem investigated here is the question of obtaining quantitative estimates on the relaxation to equilibrium. One way to measure this is to estimate the rate of convergence to zero of quantities of the form
	\begin{equation*}  %\label{e.}
	\Er \Ll[ (P_t f)^2 \Rr] = \Er \Ll[ f \, (P_{2t} f) \Rr],
	\end{equation*}
	for suitable centered functions $f : \mmd(\Rd)\to \R$, and where we used the reversibility and the semigroup property of $(P_t)_{t \ge 0}$ to derive the identity above. For specific choices of the observable $f$, the statement of Theorem~\ref{t.two.point} gives an upper bound on the next-order correction to the leading-order behavior of $\Er \Ll[ (P_t f)^2 \Rr]$. For different classes of functions $f$, the leading-order behavior was investigated in \cite{gu2020decay} for the model we also consider here. Other models were considered in \cite{berzeg, ccr, deu94, jlqy, lanyau, lig91}; see also \cite{vardecay} that relates this sort of problem with the quantitative homogenization of elliptic equations. Heat-kernel estimates for a tagged particle evolving in a simple exclusion process were obtained in \cite{giunti2019heat}. 
	
	
	\subsection*{Organization of the paper} The remainder of this article is organized as follows. In Section~\ref{sec:2}, we introduce further notation and recall several results from~\cite{bulk}, as well as the semigroups on the configuration space. In Section~\ref{sec:3}, we develop a two-scale expansion for elliptic equations on the configuration space, which we find interesting in its own right and illustrates our proof approach in the simplest possible setting. Section~\ref{sec:4} establishes Theorem~\ref{t.two.point} using a two-scale expansion for parabolic equations. Finally, in Section~\ref{sec:5} we prove the quantitative error estimates in the Green--Kubo formula stated in Theorem~\ref{t.quanti.gk}.
	
	\section{Preliminaries}
	\label{sec:2}
	\subsection{Euclidean space}
	
	We denote by $\R^d$, $d\geq 1$, the standard Euclidean space, and as stated in \eqref{eq.defcube}, we denote by~$\cu_m$ the open hypercube of side-length~$3^m$ centered at the origin. For every $n \leq m \in \N$, we set
	\begin{equation}\label{e.Z.def}
			\Z_{m,n}  := 3^n \Zd \cap \cu_m, \qquad \Z_n   := 3^n  \Zd.
	\end{equation}
	Then we can partition a large cube $\cu_m$ into cubes of a smaller scale as $\bigcup_{z \in \Z_{m,n}} (z + \cu_n)$, up to a set of null Lebesgue measure, and similarly, partition $\Rd$ into $\bigcup_{z \in \Z_{n}} (z + \cu_n)$, up to a set of null Lebesgue measure.
	
	For any open set $U \subset \Rd$, we write $L^p(U)$ with $p \geq 1$ for the classical space of functions $f$ for which $|f|^p$ has a finite Lebesgue integral over $U$, and $H^k(U)$ with $k \geq 1$ for the classical Sobolev space of order $k$ on $U$. We let $H^1_0(U)$ stand for the closure in $H^1(U)$ of the set of smooth functions with compact support in $U$. For any non-empty and bounded open set $U$, and $f \in L^1(U)$, we also introduce shorthand notation for the Lebesgue integral of $f$, normalized by the Lebesgue measure of $U$, as
	\begin{equation*}  %\label{e.}
		\fint_U  f := \frac{1}{|U|} \int_U f.
	\end{equation*}
	
	
	
	\subsection{Configuration space}
	
	The configuration space $\mmd(\Rd)$ and the Poisson point process $\Pr$ are defined at the beginning of the introduction. For a Borel set $U\subset \R^d$,  we write $\mcl F_U$ for the $\sigma$-algebra generated by the mappings $\mu \mapsto \mu(V)$, for all Borel sets $V \subset U$, completed with all the $\P_{\rho}$-null sets. We also use the shorthand notation $\mcl F$ for $\mcl F_{\Rd}$. The $\sigma$-algebra $\mcl F$ is also the Borel $\sigma$-algebra associated with the topology of vague convergence on $\mmd(\Rd)$, as explained for instance in \cite[Exercise~5.2 and solution]{HJbook}.
	
	
	We introduce several function spaces on $\mmd(\Rd)$ that will be used below. On the configuration space $\mmd(\Rd)$, we define 
    \begin{align}\label{eq.defL2}
        \cL^2(U) := L^2(\mmd(\Rd), \mcl F_U, \Pr),
    \end{align}
    the complete space of the $\mcl F_U$-measurable functions with finite second moment under $\Pr$. We write $\cL^2 := \cL^2(\Rd)$ as a shorthand notation. It is a Hilbert space and we denote by $\bracket{\cdot , \cdot}_{\cL^2}$ its associated inner product.
	
	
	For every sufficiently smooth function $f : \mmd(\Rd) \to \R$, measure $\mu \in \mmd(\Rd)$, and $x \in \supp \mu$, the gradient $\nabla f(\mu,x)$ is defined by requiring, for every ${k \in \{1,\ldots,d\}}$, that
	\begin{equation}  
		\label{e.def.deriv}
		\mathrm e_k \cdot \nabla f(\mu, x) = \lim_{h \to 0} \frac{f(\mu - \delta_x + \delta_{x + h \mathrm e_k}) - f(\mu)}{h},
	\end{equation}
	where $(\mathrm e_1,\ldots, \mathrm e_d)$ denotes the canonical basis of $\Rd$. 
	%\jccomment{Do we ever use this piece in green or can I remove it?}
	%\jc{For definiteness, we extend the function $\nabla f(\mu,\cdot)$ (defined  naturally only on the support of $\mu$, which we denote by $\supp \mu$), as 
	%\begin{equation}  
	%	\label{e.nabla.convention}
	%	\text{for every } x \notin \supp \mu, \quad \nabla f(\mu,x) := 0. 
	%\end{equation}}
	We define the spaces of smooth functions $\cC^\infty(U)$ and $\cC^{\infty}_c(U)$ for an open set $U \subset \Rd$ as follows. For every open set $V \subset \Rd$ and any $\mcl F$-measurable function $f$, suppose the configuration outside of $V$ is fixed and the number of particles in $V$ is given by $\mu(V) = n \in \N$. This naturally gives rise to a measurable canonical projection $(x_1,...,x_n) \mapsto f_n(x_1,...,x_n, \mu \mres V^c) := f(\sum_{i = 1}^n \delta_{x_i} + \mu \mres V^c)$, defined on $V^n$. We define $\cC^\infty(U)$ as the space of $\mcl F$-measurable functions $f : \mmd(\Rd) \to \R$ such that for every bounded open set $V \subset U$, $\mu \in \mmd(\Rd)$, and $n \in \N$, the canonical projection $f_n(\cdot, \mu \mres V^c)$ is infinitely differentiable on $V^n$. The space $\cC^{\infty}_c(U)$ is the subspace of $\cC^{\infty}(U)$ of functions $f$ for which there exists a compact set $K \subseteq U$ such that $f$ is $\mcl F_K$-measurable.
	
	
	We now define $\cH^1(U, \Pr)$, an infinite dimensional analogue of the classical  Sobolev space $H^1$. For every $f \in \cC^\infty(U)$, we introduce the norm
	\begin{align}\label{eq.defH1}
		\norm{f}_{\cH^1(U, \Pr)} := \Ll(\Er[f^2(\mu)] + \Er\Ll[\int_U \vert \nabla f(\mu, x) \vert^2 \, \d \mu(x) \Rr]\Rr)^{\frac{1}{2}},
	\end{align}
	and we define the space $\cH^1(U, \Pr)$ as the completion of the set of functions $f \in \cC^\infty(U)$ such that $\norm{f}_{\cH^1(U, \Pr)}$ is finite, with respect to this norm. Since the density $\rho$ in this paper is kept fixed, we write $\cH^1(U)$ as a shorthand notation for $\cH^1(U, \Pr)$. We also denote by $\acH^1(U)$ a normalized version of this norm which is convenient in order to absorb the diameter, entering as a factor in Poincar\'e's inequality (see~\eqref{e.Poincare} below), 
	\begin{align}\label{eq.defNormalH1}
		\norm{f}_{\acH^1(U)} := \Ll(\vert U\vert^{-\frac{2}{d}} \Er[f^2(\mu)] + \Er\Ll[\int_U \vert \nabla f(\mu, x) \vert^2 \, \d \mu(x) \Rr]\Rr)^{\frac{1}{2}}.
	\end{align}
	We stress that {functions in $\cH^1(U)$ need not be $\mcl F_U$-measurable}. Indeed, the function $f$ can depend on $\mu \mres U^c$ in a relatively arbitrary (measurable) way, as long as $f \in \cL^2$. If $V \subset U$ is another open set, then $\cH^1(U) \subset \cH^1(V)$. 
	%Therefore, for a general function ${f \in \cH^1(U) \backslash \D(U)}$, its gradient $\nabla f$ should always be understood as the limit of $\nabla f^{(i)}$ such that $f^{(i)} \in \D(U), {f^{(i)} \xrightarrow{\cH^1(U)} f}$.
	
	
	
	We define the space $\cH^1_0(U)$ as the closure in $\cH^1(U)$ of the space of functions $f \in \cC^\infty_c(U)$   such that $\norm{f}_{\cH^1(U)}$ is finite. Notice in particular that, in contrast with functions in $\cH^1(U)$, a function in $\cH^1_0(U)$ does not depend on $\mu \mres U^c$. If $V \subset U$ is another open set, then $\cH^1_0(V) \subset \cH^1_0(U)$.% (the inclusion goes in the opposite direction to that for $\cH^1$ spaces).
	
	%% It is commented as this part is not so useful in this paper.
	\iffalse{
	We summarize the previous remarks in Table \ref{t.differences}.
	\begin{table}[ht]
		\begin{center}
			\begin{tabular}{ c|cccc } 
				\toprule
				\shortstack{Function \\ space} & \shortstack{$H^1$-regularity for \\  particles in $U$} & \shortstack{Compatibility\\ condition when \\ particles cross $\partial U$} & \shortstack{$\mcl F_U$- \\measurable} & \shortstack{For every\\ open set  $V \subset U$} \\
				\midrule
				$\cH^1(U)$ & Yes & No&  No& $\cH^1(U) \subset \cH^1(V)$\\ 
				$\cH^1_0(U)$ & Yes &  Yes &  Yes& $\cH^1_0(V) \subset \cH^1_0(U)$\\ 
				\bottomrule
			\end{tabular}
			\medskip
			\caption{Differences between $\cH^1(U)$ and $\cH^1_0(U)$.}
			\label{t.differences}
		\end{center}
	\end{table}
	}\fi
	
	We also recall a version of Poincar\'e's inequality on the configuration space, which was proved in~\cite[Proposition 3.3]{bulk}: there exists a constant $C(d) < + \infty$ such that for every bounded open set $U \subseteq \R^d$ and $f \in \cH_0^1(U)$, one has
	\begin{equation}
		\label{e.Poincare}
		\Er\left[(f - \Er[f])^2 \right] \leq C \diam(U)^2 \Er\left[\int_U |\nabla f|^2 \d \mu \right],
	\end{equation}
	where $\diam(U)$ denotes the Euclidean diameter of $U$.
	
	We denote by $\cH^{-1}(U)$ the dual space of $\cH^1_0(U)$, which is the space of all continuous linear forms on $\cH^1_0(U)$. We denote by $\bracket{\cdot,\cdot}_{\cH^{-1}(U), \cH^1_0(U)}$ the duality pairing, and define the norm $\norm{\cdot}_{\cH^{-1}(U)}$ as
	\begin{align}\label{eq.defHdual}
		\norm{g}_{\cH^{-1}(U)} := \sup_{f \in \cH^1_0(U), \norm{f}_{\cH^1(U)} \le 1} \bracket{g,f}_{\cH^{-1}(U), \cH^1_0(U)}.
	\end{align}

	%%%% The embedding is not used finally.
	\iffalse{
	A function $g \in \cL^2$ defines a natural linear functional in $\cH^{-1}(U)$ by embedding in the sense 
	\begin{align}\label{eq.embedding}
		\forall f \in \cH^1(U), \qquad \bracket{g, f}_{\cH^{-1}(U), \cH^1(U)} := \bracket{g, f}_{\cL^2}.
	\end{align}
	}\fi
	
	
	
	\subsection{Homogenization of the bulk diffusion matrix}
	
	In this subsection, we recall some results from \cite{bulk}, including the definition of the bulk diffusion matrix $\ab$ and its quantitative convergence from finite-volume approximations $\ab(\cu_m)$ and $\ab_\ast(\cu_m)$.
	
	For a bounded domain $U \subset \Rd$, vectors $p,q \in \Rd$, and $\rho > 0$, we consider the following optimization problems
	\begin{equation}\label{eq:defNu}
		\begin{split}
			\nu(U,p,\rho) &:= \inf_{v \in \ell_{p,U} + \cH^1_0(U)}\Er \Ll[ \frac{1}{\rho \vert U \vert} \int_{U} \frac{1}{2} \nabla v \cdot \a \nabla v \, \d \mu \Rr], \\
			\nu^*(U,q,\rho) &:=  \sup_{u \in \cH^1(U)} \Er \Ll[\frac{1}{\rho \vert U \vert}\int_{U} \Ll( -\frac 1 2 \nabla u \cdot \a \nabla u + q \cdot \nabla u \Rr) \, \d \mu \Rr],
		\end{split}
	\end{equation}
	with $\ell_{p,U} = \int_U p \, \cdot \, x \, \d \mu(x)$. By~\cite[Proposition 4.1]{bulk}, there exist symmetric $(d\times d)$-matrices $\ab(U,\rho), \ab_\ast(U,\rho)$ satisfying ${\id \leq \ab(U,\rho) \leq \Lambda \id}$ and ${\id \leq \ab_\ast(U,\rho) \leq \Lambda \id}$ such that for every $p,q \in \Rd$ and $\rho > 0$,
	\begin{equation}\label{eq:MatrixEq}
			\nu(U,p,\rho) = \frac{1}{2} p \cdot \ab(U,\rho) p, \quad  \quad \nu^\ast(U,q,\rho) = \frac{1}{2} q \cdot \ab_\ast^{-1}(U,\rho) q.
	\end{equation}
	The sequence $(\ab(\cu_m,\rho))_{m \in \N}$ is decreasing in $m$ (in the sense that $\ab(\cu_{m},\rho) - \ab(\cu_{m+1},\rho)$ is positive semi-definite for every $m\in \N$), and we denote its limit by
	\begin{equation}
		\label{eq:AbarDef}
		\ab(\rho) := \lim_{m \rightarrow \infty} \ab(\cu_m,\rho).
	\end{equation}
	We also refer to \cite[Appendix~B]{bulk} for equivalent definitions of $\ab(\rho)$. 
	The main result of~\cite{bulk} is a quantification of the speed of convergence in~\eqref{eq:AbarDef}. More precisely, by the proof of~\cite[Theorem 2.1]{bulk}, there exist constants $\alpha = \alpha(d,\Lambda,\rho) > 0$ and $C = C(d,\Lambda,\rho) < + \infty$ such that
	\begin{equation}
		\label{eq:QuantConvergence}
		|\ab(\cu_m,\rho) - \ab(\rho)|+ |\ab_\ast(\cu_m,\rho) - \ab(\rho)| \leq C 3^{-\alpha m}. 
	\end{equation}
	Concerning the optimization problems in \eqref{eq:defNu}, we write
	\begin{equation}
		\begin{split}
			v(\cdot,U,p) \in \ell_{p,U} + \cH^1_0(U) &\text{ for the minimizer of $\nu(U,p,\rho)$, and} \\
			u(\cdot,U,q) \in \cH^1(U) &\text{ for the maximizer of $\nu^\ast(U, q,\rho)$.}
		\end{split}
	\end{equation}
As observed in \cite[Proposition~4.1]{bulk}, these optimizers are unique provided that we also impose that $\Er[v(\mu,U,p)] = 0$ and that $\Er[u(\mu,U,q)| \mu(U), \mu \mres U^c] = 0$. 	
	By \cite[(72)]{bulk}, the first-order optimality conditions read as
	\begin{equation}
		\label{eq:VariationalFormula72}
		\begin{split}
			\forall v' \in \cH^1_0(U), \qquad \E_\rho\left[\int_U \nabla v(\cdot,U,p) \cdot \a \nabla v'  \d \mu\right] &= 0,\\
			\forall u' \in \cH^1(U), \qquad \E_\rho\left[\int_U \nabla u(\cdot,U,q) \cdot \a \nabla u'  \d \mu\right] &= \E_\rho\left[\int_U  q\cdot \nabla u' \d \mu \right].
		\end{split}
	\end{equation}
	By removing the affine part from $v(\cdot,U,p)$ and $u(\cdot,U,q)$, we define the \textit{approximate correctors}
	\begin{align}\label{eq.corrector}
		\phi_{p, U} := v(\cdot,U,p) - \ell_{p,U}, \qquad \phi^*_{p, U} := u(\cdot,U,\ab_\ast(U,\rho)p) - \ell_{p,U}.
	\end{align}
	Here the definition of $\phi^*_{p, U}$ is motivated as follows. As seen in \cite[(65)]{bulk}, we can combine \eqref{eq:VariationalFormula72} and \eqref{eq:MatrixEq} to see that the averaged slope of $ u(\cdot,U,q)$ is $\ab_\ast^{-1}(U,\rho) q$:
	\begin{align}\label{eq.slopeDuel}
		\forall q' \in \Rd, \qquad \E_\rho\left[ \frac{1}{\rho \vert U \vert}\int_U  q'\cdot \nabla u(\cdot,U,q) \d \mu \right] = q' \cdot \ab_\ast^{-1}(U,\rho) q.
	\end{align}
	Therefore, the function $u(\cdot,U,\ab_\ast(U,\rho)p)$ has an average slope of $p$. %We call these functions \emph{approximate} correctors because they are defined in finite volume.
	
	The following lemma subsumes important properties of correctors that will be instrumental in the following sections. Its proof essentially follows from~\cite[Sections~4 and~5]{bulk}.
	\begin{lemma}\label{lem.JBound}
		(1) For every $n \in \N$, $z,z' \in \mathcal{Z}_{n}$ with $\dist(z,z') > 3^n$, and $i \in \{1,\cdots,d\}$, the quantities $\phi_{e_i,z + \cu_n}$ and  $\phi_{e_i,z' + \cu_n}$ are independent random variables.
		
		(2) There exist constants $\alpha(d,\Lambda,\rho) > 0$ and $C(d,\Lambda,\rho) < + \infty$  such that for every $n\in \N$ and $i \in \{1,\cdots,d\}$, we have
		\begin{equation}\label{eq.JBound}
			 \Er\Ll[\frac{1}{\rho|\cu_n|}\int_{\cu_n} \vert \nabla \phi_{e_i, \cu_n} - \nabla \phi^*_{e_i, \cu_n} \vert^2 \, \d \mu\Rr] \leq C3^{-  2\alpha n}.
		\end{equation}
	\end{lemma}
	\begin{proof}
		The claim (1) follows directly from observing that $\phi_{e_i,z + \cu_n}$ and  $\phi_{e_i,z' + \cu_n}$ are in $\cH^1_0(z + \cu_n)$ and $\cH^1_0(z' + \cu_n)$ respectively.% and using the finite-range assumption~\eqref{e.finite.range}.
		
		Concerning the claim (2), for every $n \in \N$ and $p \in \R^d$, we introduce the quantity
		\begin{equation}\label{e.def.J}
			\begin{split}
				J(\cu_n,p,\ab_\ast(\cu_n,\rho) p) &:= \nu(\cu_n,p,\rho) + \nu^*(\cu_n,\ab_\ast(\cu_n,\rho)p,\rho ) - p \cdot \ab_\ast(\cu_n,\rho)p \\
				& = \frac{1}{2} p \cdot \ab(\cu_n,\rho)p - \frac{1}{2} p \cdot \ab_\ast(\cu_n,\rho) p.
			\end{split}
		\end{equation}
		By~\cite[Proposition 4.2, (86) and (87)]{bulk}, we see that upon writing
		\begin{equation}\label{eq.wDifference}
			\begin{split}
				w(\mu,\cu_n,e_i) &:= u(\mu,\cu_n,\ab_\ast(\cu_n,\rho)e_i) - v(\mu,\cu_n,e_i)\\
				&= \phi_{e_i, \cu_n} - \phi^*_{e_i, \cu_n},
			\end{split}
		\end{equation}
		the quantity $J$ admits the quadratic representation
		\begin{equation}\label{eq.JQuadra}
			J(\cu_n,e_i,\ab_\ast(\cu_n,\rho) e_i) = \Er\Ll[\frac{1}{\rho \vert \cu_n \vert} \int_{\cu_n} \frac{1}{2} \nabla  w(\mu,\cu_n,e_i)  \cdot \a  \nabla  w(\mu,\cu_n,e_i)  \, \d \mu \Rr].
		\end{equation}
		Here we remark that our $w$ is different from that in \cite[(86)]{bulk} up to a term that may only depend on $\mu(\cu_n)$ and $\mu \mres \cu_n^c$, but this does not change the value of the quadratic form on the right side of \eqref{eq.JQuadra}. Combining \eqref{eq.wDifference}, \eqref{eq.JQuadra}, and \eqref{e.ass.unif.ell},  we obtain that
		\begin{equation}\label{e.UpperBound.Corrector.Diff}
			 \Er\Ll[\frac{1}{\rho|\cu_n|}\int_{\cu_n} \vert \nabla \phi_{e_i, \cu_n} - \nabla \phi^*_{e_i, \cu_n} \vert^2 \, \d \mu\Rr]  \leq C J(\cu_n,e_i,\ab_\ast(\cu_n,\rho) e_i).
		\end{equation}
		On the other hand, by~\cite[Section 5.4]{bulk}, one knows that
		\begin{equation}
			\label{e.UpperBound.J}
			J(\cu_n,e_i,\ab_\ast(\cu_n,\rho) e_i)  \leq C3^{-2\alpha n}.
		\end{equation}
		Upon combining~\eqref{e.UpperBound.Corrector.Diff} and~\eqref{e.UpperBound.J}, the claim \eqref{eq.JBound} follows.
	\end{proof}
	
	\subsection{Semigroup and diffusion on configuration space}\label{subsec.Dirichletform}
	In this part, we summarize some properties of the semigroup and the associated diffusion process on the configuration space $\mmd(\Rd)$. One possible route to the definition of these objects would be to start by considering initial configurations with a finite number of particles. If $N$ particles are involved, the semigroup and stochastic process we aim to define are those associated with a divergence-form operator on $(\R^d)^N$, and standard methods apply: one can for instance construct the transition probabilities as in \cite[Appendix~E]{AKMbook} and then construct the stochastic process using the Kolmogorov extension and continuity theorems. Once this is done, one can then seek to define a stochastic process for sufficiently ``spread out'' initial point clouds (in particular for $\P_\rho$-almost every configuration) through a limit procedure. This is the method used in \cite{lang1977unendlich} for a slightly different model. 
	
	Since the stochastic process will not play a major role in this paper, we prefer to explain a more direct approach using the theory of Dirichlet forms. This approach has been explored in detail in \cite{albeverio1996differential, albeverio1996canonical, albeverio1998analysis, albeverio1998analysis2, kondratiev2003heat, ma2000construction, osada1996dirichlet, rockner1998stochastic, yoshida1996construction}, and the necessary conditions for the definition of the semigroup and stochastic process are easy to check for our model. 
	
	We endow the space $\mmd(\Rd)$ with the topology of vague convergence, and define the symmetric form on $\cL^2$ ($= L^2(\mmd(\Rd), \mcl F, \Pr)$) by setting
	\begin{align}\label{eq.Dirichlet}
		\Ll\{\begin{aligned}
			& \D(\Ee^{\a}) := \cH^1(\Rd),\\
			& \displaystyle{\forall u,v \in \D(\Ee^{\a}), \quad 			\Ee^{\a}(u,v) := \Er\Ll[\int_{\Rd} \frac{1}{2}\nabla u \cdot \a \nabla v \, \d \mu \Rr].}
		\end{aligned}\Rr.
	\end{align}
	The symmetric form $(\Ee^\a, \D(\Ee^\a))$ is closed: this means that if a sequence of elements of~$\D(\Ee^\a)$ is a Cauchy sequence with respect to the norm induced by the symmetric form $\Ee^{\a}_1(u,v) = \Ee^{\a}(u,v) + \bracket{u, v}_{\cL^2}$, then this sequence converges to an element of~$\D(\Ee^\a)$ with respect to this norm (see also \cite[Section~1.1]{fukushima2010dirichlet}). We can also verify that $\Ee^\a$ is Markovian, since for the unit contraction function $\phi_{\epsilon} \in C^\infty(\R)$ as in \cite[(1.1.5), Exercise~1.2.1]{fukushima2010dirichlet}, we have that $\vert \phi'_{\epsilon} \vert \leq 1$ and 
	\begin{align*}
		\forall u \in \D(\Ee^{\a}), \quad \Ee^\a(\phi_{\epsilon}(u), \phi_{\epsilon}(u)) = \Er\Ll[\int_{\Rd} \frac{1}{2}\vert \phi'_{\epsilon}(u)\vert^2 \nabla u \cdot \a \nabla u \, \d \mu \Rr] \leq \Ee^\a(u, u).
	\end{align*} 
 Therefore, $(\Ee^\a, \D(\Ee^\a))$ is a Dirichlet form. By the Riesz representation theorem (see \cite[Theorem~1.3.1]{fukushima2010dirichlet}), the Dirichlet form $(\Ee^\a, \D(\Ee^\a))$ defines a unique self-adjoint positive semi-definite operator $-A$ on $\cL^2$ such that 
	\begin{align}\label{eq.defA}
		\Ll\{\begin{array}{ll}
			\D(\Ee^{\a}) & = \D(\sqrt{-A}), \\
			\Ee^{\a}(u,v) & = \bracket{\sqrt{-A}u, \sqrt{-A}v}_{\cL^2},
		\end{array}\Rr.
	\end{align}
	where $\D(\sqrt{-A})$ denotes the domain of the operator $\sqrt{-A}$. Using the spectral decomposition, the self-adjoint operator $-A$ then defines a strongly continuous semigroup $P_t := \exp(tA)$ with $t \geq 0$ on $\cH^1(\Rd)$ (see \cite[Lemma~1.3.2]{fukushima2010dirichlet}). We will also exploit the close relationship between the semigroup and the resolvent of $-A$, which is the family of operators $(\lambda - A)^{-1}$ for $\lambda$ ranging in $(0,+\infty)$, and we refer to \cite[Section~1.3]{fukushima2010dirichlet} for more on this. We may use the informal notation $\frac{1}{2} \nabla \cdot \a \nabla$ to denote the operator $A$. 
	
	Since $\mmd(\Rd)$ is not locally compact with respect to the vague topology, we construct the associated Markov process using the notion of quasi-regular Dirichlet forms; see \cite[Chapter~IV.3]{ma2012introduction}, \cite{ma2000construction}, and here we follow the steps in \cite[Section~4]{rockner1998stochastic}. We denote by $\Ee$ the Dirichlet form with identity diffusion matrix
	\begin{align*}
		\Ll\{\begin{aligned}
			& \D(\Ee) := \cH^1(\Rd),\\
			& \displaystyle{\forall u,v \in \D(\Ee), \quad 			\Ee(u,v) := \Er\Ll[\int_{\Rd} \frac{1}{2}\nabla u \cdot  \nabla v \, \d \mu \Rr].}
		\end{aligned}\Rr.
	\end{align*}
	Our Dirichlet form $\Ee^{\a}$ is comparable to $\Ee$ because of the uniform ellipticity in \eqref{e.ass.unif.ell}, 
	\begin{align}\label{eq.Compare}
		\forall u \in \D(\Ee) = \D(\Ee^\a), \qquad \Ee(u,u) \leq \Ee^\a(u,u) \leq \Lambda \Ee(u,u).
	\end{align} 
	The quasi-regularity of $\Ee$ is proved in \cite[Theorem~4.12]{rockner1998stochastic} using the result in \cite{ma2000construction}. With the comparison in \eqref{eq.Compare}, the notions of $\Ee$-nest, $\Ee$-exceptional, and $\Ee$-quasi-continuous in \cite[Definitions~4.10 and~4.11]{rockner1998stochastic} are equivalent to those of $\Ee^\a$-nest, $\Ee^\a$-exceptional, and $\Ee^\a$-quasi-continuous, respectively. From this, we obtain  the quasi-regularity of $\Ee^\a$. The same argument also proves the locality of $\Ee^\a$ viewing that of $\Ee$ proved in \cite[Corollary~4.13]{rockner1998stochastic}. As a consequence of \cite[Theorems~IV.3.5, V.1.5 and~V.1.11]{ma2012introduction}, there exists an  $\mmd(\Rd)$-valued Markov process $(\bmu_t)_{t \geq 0}$ that is canonically associated with the Dirichlet form $\Ee^\a$, and we denote by $\PP_{\bmu_0}$ and $\EE_{\bmu_0}$ the probability and expectation associated with this process starting from the configuration $\bmu_0$. For every $u \in \cL^2$ and $\Ee^\a$-quasi every $\bmu_0 \in \mmd(\Rd)$, we have that 
	\begin{align*}
(P_t u)(\bmu_0) = \EE_{\bmu_0}[u(\bmu_t)].
	\end{align*}
	Moreover, $\Pr$ is a reversible measure for  $(\bmu_t)_{t \geq 0}$.
	

	In this paper, we mostly focus on the properties of the semigroup. 
	\begin{proposition}[Elementary properties of the semigroup]\label{prop.ElementSemigroup}
		For every $u \in \cL^2$, the function $u_t := P_t u \in \cH^1(\Rd)$, $t > 0$, satisfy the following properties.
		\begin{enumerate}
			\item (Energy solution) For every $T > 0$, the function $u_\cdot (\cdot) : [0,T] \times \mmd(\Rd) \to \R$ is the unique solution in $L^2([0,T], \cH^1(\Rd)) \cap C([0,T], \cL^2)$ with $\partial_t u_t \in L^2([0,T], \cH^{-1}(\Rd))$ such that for every $t \in [0,T]$,
			\begin{align}\label{eq.semigroupA}
				\forall v \in \cH^1(\Rd), \qquad \Er[u_t v] - \Er[u v] + \int_{0}^t \Ee^\a(u_s, v) \, \d s = 0.
			\end{align}	
			The function $u_\cdot(\cdot)$ also satisfies the energy evolution identity
			\begin{align}\label{eq.L2Energy}
				\forall \tau \in [0, t), \qquad  \frac{1}{2}\norm{u_\tau}^2_{\cL^2} = \frac{1}{2}\norm{u_t}^2_{\cL^2} + \int_{\tau}^t  \Ee^\a (u_s, u_s) \, \d s.
			\end{align}
			\item (Decay of $\cL^2$ and energy) The mappings $t \mapsto \Er[(u_t)^2]$ and $t \mapsto \Ee^{\a}(u_t, u_t)$ are decreasing. Moreover, we have
			\begin{align}\label{eq.L2Decay}
				\forall \tau \in [0, t), \qquad \Ee^\a(u_t, u_t)  \leq \frac{1}{2(t - \tau)}\Er[(u_\tau)^2].    
			\end{align}
			Finally, when $u \in \D(\Ee^\a)$, we have $\Ee^{\a}(u_t, u_t) \xrightarrow{t \to 0} \Ee^{\a}(u, u)$.
		\end{enumerate}
	\end{proposition}
	\begin{proof}
		Following the classical energy solution theory of the evolution equation, there exists a unique solution $\tilde u_t \in L^2([0,T], \cH^1(\Rd)) \cap C([0,T], \cL^2)$ and $\partial_t \tilde u_t \in L^2([0,T], \cH^{-1}(\Rd))$; see \cite[Sections~7.1.2.b and 7.1.2.c, Theorems~2, 3, 4]{evans2022partial}. Moreover, using the definition \eqref{eq.defA}, $u_t = P_t u = e^{tA} u$ satisfies that for every $v \in \cH^1(\Rd)$ 
		\begin{align*}
			\bracket{u_t, v}_{\cL^2} - \bracket{u, v}_{\cL^2} &= \bracket{e^{tA} u, v}_{\cL^2} - \bracket{u, v}_{\cL^2} \\
			&= \int_{0}^t \bracket{A e^{sA}u, v}_{\cL^2} \, \d s \\
			&= -\int_{0}^t \bracket{\sqrt{-A} e^{sA}u, \sqrt{-A} v}_{\cL^2} \, \d s \\
			&= -\int_{0}^t \Ee^\a(u_s, v) \, \d s.
		\end{align*}
		Then the uniqueness of the solution in \eqref{eq.semigroupA} implies that $u_t = \tilde u_t$.  Equation \eqref{eq.L2Energy} follows from \eqref{eq.semigroupA} by testing with $v = u_t$, and it also implies that  $t \mapsto \Er[(u_t)^2]$ is decreasing.
		

		To prove that $t \mapsto \Ee^{\a}(u_t, u_t)$ is also decreasing, we use the spectral measure decomposition for the Dirichlet form
		\begin{align}\label{eq.Spectral}
			\Ee^{\a}(u_t, u_t) = \bracket{e^{tA} u, - A e^{tA} u}_{\cL^2} = \int_{\lambda = 0}^\infty \lambda e^{-2 \lambda t} \,  \d \bracket{E_\lambda u, u}_{\cL^2}. 
		\end{align}
		Here $E_{\lambda}$ is the spectral projection operator with respect to $-A$; see \cite[Lemma~1.3.2]{fukushima2010dirichlet}. We see that the quantity above is decreasing with respect to $t$, as announced. We put this property into \eqref{eq.L2Energy}, which gives 
		\begin{align*}
			(t-\tau)\Ee^\a (u_t, u_t) \leq \int_{\tau}^t  \Ee^\a (u_s, u_s) \, \d s \leq  \frac{1}{2}\norm{u_\tau}^2_{\cL^2}. 
		\end{align*}
		This shows \eqref{eq.L2Decay}. When $u \in \D(\Ee^\a)$, we have 
		\begin{align*}
			\Ee^{\a}(u, u) = \bracket{\sqrt{- A} u,  \sqrt{- A} u}_{\cL^2} = \int_{\lambda = 0}^\infty \lambda \,  \d \bracket{E_\lambda u, u}_{\cL^2} < + \infty. 
		\end{align*}
Combining this with an application of the dominated convergence theorem in  \eqref{eq.Spectral}, we obtain that $\lim_{t \to 0}\Ee^{\a}(u_t, u_t)= \Ee^{\a}(u, u)$.
	\end{proof}
	
	
	
	\section{Two-scale expansion for elliptic equation}
	\label{sec:3}
	In this section, we take the Dirichlet problem of elliptic equations as a simple example to explain the two-scale expansion technique on the configuration space. Recall the homogenized matrix $\ab$ defined in \eqref{eq:AbarDef}, where we keep the dependence on $\rho$ implicit. Consider $U, \Ub \in \cH^1_0(\cu_m)$ centered solving
	\begin{align}\label{eq.defDirichlet}
		-\nabla \cdot (\a \nabla U) = {F}, \qquad -\nabla \cdot (\ab \nabla \Ub) =  F,
	\end{align}
	in the sense that for all $V \in \cH^1_0(\cu_m)$,
	\begin{align}
		\Er\Ll[\int_{\cu_m} \nabla V \cdot \a \nabla U \, \d \mu \Rr] &= \Er[V F],  \label{eq.Dirichlet1}\\
		\Er\Ll[\int_{\cu_m} \nabla V \cdot \ab \nabla \Ub \, \d \mu \Rr] &= \Er[V F].     \label{eq.Dirichlet2}       
	\end{align}
	We study the homogenization problem in the case when $F$ is a centered linear statistic, meaning that
	\begin{align}\label{eq.defF}
		F(\mu) := \int_{\cu_m} f \, \d \mu,  \qquad \text{ with } \qquad \int_{\cu_m} f \, \d m = 0.
	\end{align}
	Here $f \in L^2(\cu_m)$ is a function on Euclidean space $\Rd$ rather than the configuration space. This particular case is convenient because its homogenized solution has an explicit expression related to the homogenized solution in Euclidean space.
	
	\begin{lemma}
		Under the condition \eqref{eq.defF}, the unique solution of \eqref{eq.Dirichlet2} such that $\Er[\Ub] = 0$ is
		\begin{align}\label{eq.defUb}
			\Ub(\mu) = \int_{\cu_m} \ub \, \d \mu - \rho \int_{\cu_m} \ub \, \d m, 
		\end{align}
		with $\ub \in H^1_0(\cu_m)$ the unique solution of 
		\begin{equation}\label{eq.defub}
			\Ll\{\begin{array}{ll}
				- \nabla \cdot (\ab \nabla \ub) = f, & \text{ in } \cu_m, \\
				\ub = 0, & \text{ on } \partial \cu_m. 
			\end{array}\Rr. 
		\end{equation}
	\end{lemma}
	\begin{proof} It is clear that $\Ub \in \cH^1_0(\cu_m)$ as $\ub \in H^1_0(\cu_m)$. If we test \eqref{eq.Dirichlet2} with some $V \in \cH^1_0(\cu_m)$, then its {\lhs} is 
		\begin{multline}\label{eq.IPPLHS}
			\Er\Ll[\int_{\cu_m} \nabla V \cdot \ab \nabla \Ub \, \d \mu \Rr] \\
			= e^{-\rho \vert \cu_m \vert} \sum_{k=1}^{\infty} \frac{(\rho \vert \cu_m \vert)^k}{k!} \fint_{(\cu_m)^k} \sum_{i=1}^k \nabla_{x_i} \tilde{V}_k(x_1, \cdots, x_k) \cdot \ab \nabla \ub(x_i) \, \d x_1 \cdots \d x_k,
		\end{multline}
		where $\tilde{V}_k(x_1, \cdots, x_k) := V(\sum_{i=1}^k \delta_{x_i})$ is the projection of the function conditioned on the number of particles. We treat the integration with respect to $x_i$ for each term 
		\begin{equation}\label{eq.IPP}
			\begin{split}
				&\int_{\cu_m}  \nabla_{x_i} \tilde{V}_k(x_1, \cdots, x_k) \cdot \ab \nabla \ub(x_i) \, \d x_i \\
				&= \int_{\cu_m}   \nabla_{x_i} \Ll(\tilde{V}_k(x_1, \cdots, x_k) - \tilde{V}_{k-1}(x_1, \cdots, x_{i-1},  x_{i+1}, \cdots, x_k) \Rr) \cdot \ab \nabla \ub(x_i)  \, \d x_i \\
				&=  \int_{\cu_m}   \Ll(\tilde{V}_k(x_1, \cdots, x_k) - \tilde{V}_{k-1}(x_1, \cdots, x_{i-1},  x_{i+1}, \cdots, x_k) \Rr) f(x_i)  \, \d x_i \\
				&=  \int_{\cu_m}   \tilde{V}_k(x_1, \cdots, x_k)  f(x_i)  \, \d x_i. \\
			\end{split}
		\end{equation}	
		The equality between the expressions in the second and third lines is valid by the definition of $\bar u$ in \eqref{eq.defub} and since the mapping
		\begin{align*}
			x_i \mapsto \Ll(\tilde{V}_k(x_1, \cdots, x_k) - \tilde{V}_{k-1}(x_1, \cdots, x_{i-1},  x_{i+1}, \cdots, x_k) \Rr)
		\end{align*}
		is in $H^1_0(\cu_m)$ (see for instance \cite[(3.8)]{giunti2021smoothness}). From the third line to the fourth line, we make use of the property $\int_{\cu_m} f \, \d m = 0$ in \eqref{eq.defF} and $\tilde{V}_{k-1}(x_1, \cdots, x_{i-1},  x_{i+1}, \cdots, x_k)$ constant for $x_i \in \cu_m$. We put this back into \eqref{eq.IPPLHS} and obtain that 
		\begin{align*}
			&\Er\Ll[\int_{\cu_m} \nabla V \cdot \ab \nabla \Ub \, \d \mu \Rr] \\
			&= e^{-\rho \vert \cu_m \vert} \sum_{k=1}^{\infty} \frac{(\rho \vert \cu_m \vert)^k}{k!} \fint_{(\cu_m)^k} \tilde{V}_k(x_1, \cdots, x_k) \sum_{i=1}^k f(x_i) \, \d x_1 \cdots \d x_k \\
			&= \Er[V F]. 
		\end{align*}
		Here we remark that the case $\mu(\cu_m) = 0$ does not contribute in the definition of $F$ in~\eqref{eq.defF}. This finishes the proof.
	\end{proof}
	
	\bigskip
	
	
	Now, with the explicit expression \eqref{eq.defUb}, we propose the following two-scale expansion 
	\begin{align}\label{eq.defW}
		W := \Ub + \sum_{i=1}^d \sum_{z \in \Z_{m,n}} (\partial_i \ub)_{z+\cu_n} \phi_{e_i, z+\cu_n},
	\end{align}
	to approximate both $U$ and $\Ub$, i.e. $W \simeq U \simeq \Ub$ in $\cL^2$, and also with $W \simeq U$ in $\cH^1(\cu_m)$. The integer $n \le m$ will be chosen shortly depending on $m$. We recall that $\Z_{m,n} = 3^n \Zd \cap \cu_m$ as in~\eqref{e.Z.def} and $\phi_{e_i, z+\cu_n}$ is the $\cH^1_0(z+\cu_n)$ corrector defined in \eqref{eq.corrector}, and we define $(\partial_i \ub)_{z+\cu_n}$ as 
	\begin{align*}
		(\partial_i \ub)_{z+\cu_n} := \fint_{z+\cu_n} \partial_i \ub,
	\end{align*}
	which is a deterministic constant. Recall that $\ell_{e_i, z+\cu_n} + \phi_{e_i, z+\cu_n}$ is $\a$-harmonic on the domain $z+\cu_n$ by~\eqref{eq:VariationalFormula72} and~\eqref{eq.corrector}. We use $(\partial_i \ub)_{z+\cu_n}$ to represent the slope in the mesoscopic domain $z+\cu_n$, and use $ (\partial_i \ub)_{z+\cu_n} \phi_{e_i, z+\cu_n}$  to correct the solution in $z+\cu_n$. The following result can be seen as a generalization of two-scale expansion on the configuration space. Notice that the $\cH^1$ estimate and Poincar\'e inequality on \eqref{eq.defDirichlet} imply that 
	\begin{align}\label{eq.ClassicalH1}
		\norm{U}_{\acH^1(\cu_m)} \leq C 3^m \norm{F}_{\cL^2}, \qquad \norm{\Ub}_{\acH^1(\cu_m)} \leq C 3^m \norm{F}_{\cL^2},
	\end{align}
	so the error term in the following result is smaller than the bound above on large scales.
	
	\begin{proposition}[Two-scale expansion of elliptic equation]\label{prop.TwoScaleElliptic}
    The following holds for the exponent $\alpha(d,\lambda,\rho) > 0$  given by Lemma~\ref{lem.JBound} and for some constant $C(d,\lambda,\rho) < +\infty$. 
    Under the condition \eqref{eq.defF}, with the two-scale expansion \eqref{eq.defW} associated with \eqref{eq.Dirichlet1}, \eqref{eq.Dirichlet2} and with the choice of $n := \lfloor \frac {m} {1+\alpha} \rfloor$, we have
		\begin{align}\label{eq.TwoScaleElliptic}
			\norm{W - U}_{\acH^1(\cu_m)} + 3^{-m} \Ll(\norm{W-\Ub}_{\cL^2} + \norm{U-\Ub}_{\cL^2}\Rr)	\leq C 3^{\frac{m}{1+\alpha}} \norm{F}_{\cL^2}.
		\end{align}
	\end{proposition}
	\begin{proof} 			
		\textit{Step~1: estimate of $\nabla (W - U)$.} The main step of the proposition is to study $\nabla (W - U)$, because the two-scale expansion provides $\cH^1$ approximation for $U$. Using \eqref{eq.defDirichlet}, we take any $V \in \cH^1_0(\cu_m)$ and study 
		\begin{equation}\label{eq.ellipticDecom}
			\begin{split}
				\Er\Ll[\int_{\cu_m} \nabla V \cdot \a \nabla (W-U) \, \d \mu\Rr] &=\Er\Ll[\int_{\cu_m} \nabla V \cdot (\a \nabla W - \ab \nabla \Ub) \, \d \mu\Rr] \\
				&= \mathbf{I} + \mathbf{II},
			\end{split}
		\end{equation}
		which we decompose as the sum of the following terms 
		\begin{align*}
			\mathbf{I} &:= \sum_{z \in \Z_{m,n}}  \Er\Ll[\int_{z+\cu_n} \nabla V \cdot (\a-\ab)(\nabla \ub - (\nabla \ub)_{z+\cu_n}) \, \d \mu\Rr], \\
			\mathbf{II} &:= \sum_{i = 1}^d \sum_{z \in \Z_{m,n}}  (\partial_i \ub)_{z+\cu_n} \Er\Ll[\int_{z+\cu_n} \nabla V \cdot \Ll(\a(e_i + \nabla \phi_{e_i, z+\cu_n}) - \ab e_i\Rr) \, \d \mu\Rr].
		\end{align*}
		This decomposition is even simpler than the classical two-scale expansion in $\Rd$, because here the factor $(\partial_i \ub)_{z+\cu_n}$ is totally deterministic and no derivatives act on it. On the other hand, the corrector $\phi_{e_i, z+\cu_n}$ already contains all the influence from the interaction of particles in $z+\cu_n$. Here the term $\mathbf{I}$ can be seen as the error to fix the slope in the mesoscopic scale $3^n$, while the term $\mathbf{II}$ is the error of the flux replacement. We estimate them separately.
		
		\smallskip
		\textit{Step~2.1: term $\mathbf{I}$ as the error to fix the slope.} For the term $\mathbf{I}$, we apply the Cauchy--Schwarz inequality to get that
		\begin{align*}
			\vert \mathbf{I} \vert &\leq \Lambda\norm{V}_{\acH^1(\cu_m)}  \Ll(\sum_{z \in \Z_{m,n}}\Er\Ll[\int_{z+\cu_n} (\nabla \ub - (\nabla \ub)_{z+\cu_n})^2 \, \d \mu\Rr]\Rr)^{\frac{1}{2}}.
		\end{align*}
		Using Poincar\'e's inequality, we have
		\begin{equation}\label{eq.TermI}
			\begin{split}
				\sum_{z \in \Z_{m,n}}\Er\Ll[\int_{z+\cu_n} (\nabla \ub - (\nabla \ub)_{z+\cu_n})^2 \, \d \mu\Rr] &\leq  C 3^{2n} \rho \sum_{z \in \Z_{m,n}} \int_{z+\cu_n} \vert\nabla \nabla \ub\vert^2 \, \d m  \\
				&= C 3^{2n} \rho \int_{\cu_m} \vert\nabla \nabla \ub\vert^2 \, \d m.
			\end{split}
		\end{equation}
		
		
		
		\smallskip
		\textit{Step~2.2: term $\mathbf{II}$ as the error of flux replacement.} For the term $\mathbf{II}$, we would like to apply the flux replacement as the slope is already fixed by $(\partial_i \ub)_{z+\cu_n}$. However, as $V$ is not a local function belonging to $\cH^1_0(z+\cu_n)$, some manipulation is required to test against it. We replace $\phi_{e_i, z+\cu_n}$ by the dual corrector $\phi^*_{e_i, z+\cu_n}$ associated to the dual quantity $\nu^*(z+\cu, \ab_*(z+\cu_n) e_i)$, since $\phi^*_{e_i, z+\cu_n}$ admits more test functions from \eqref{eq:VariationalFormula72}.
		\begin{align*}
			\mathbf{II} &= \mathbf{II.1} + \mathbf{II.2},\\
			\mathbf{II.1} &:= \sum_{i = 1}^d \sum_{z \in \Z_{m,n}}  (\partial_i \ub)_{z+\cu_n} \Er\Ll[\int_{z+\cu_n} \nabla V \cdot \a\Ll( \nabla \phi_{e_i, z+\cu_n} - \nabla \phi^*_{e_i, z+\cu_n}\Rr) \, \d \mu\Rr],\\
			\mathbf{II.2} &:= \sum_{i = 1}^d \sum_{z \in \Z_{m,n}}  (\partial_i \ub)_{z+\cu_n} \Er\Ll[\int_{z+\cu_n} \nabla V \cdot \Ll(\a(e_i + \nabla \phi^*_{e_i, z+\cu_n}) - \ab e_i\Rr) \, \d \mu\Rr].
		\end{align*}
		For the term $\mathbf{II.1}$, we use the Cauchy--Schwarz inequality, which gives us
		\begin{equation}\label{eq.TermII1}
			\begin{split}
				\vert \mathbf{II.1}\vert &\leq C \norm{V}_{\acH^1(\cu_m)}\Ll(\sum_{i = 1}^d \sum_{z \in \Z_{m,n}}  \vert (\partial_i \ub)_{z+\cu_n} \vert^2  \Er\Ll[\int_{z+\cu_n} \vert \nabla \phi_{e_i, z+\cu_n} - \nabla \phi^*_{e_i, z+\cu_n} \vert^2 \, \d \mu\Rr]\Rr)^{\frac{1}{2}}\\
				&\leq C 3^{-\alpha n} \norm{V}_{\acH^1(\cu_m)}\Ll(\sum_{i = 1}^d \sum_{z \in \Z_{m,n}}  \rho  \vert \cu_n \vert\vert (\partial_i \ub)_{z+\cu_n} \vert^2 \Rr)^{\frac{1}{2}}\\
				&\leq C 3^{-\alpha n} \norm{V}_{\acH^1(\cu_m)} \norm{\Ub}_{\acH^1(\cu_m)}
			\end{split}
		\end{equation}
		Here the error from $\nabla (\phi_{e_i, z+\cu_n} - \phi^*_{e_i, z+\cu_n})$ %is that of $J$ quantity studied in \cite[Section 5]{bulk}
		could be bounded by employing Lemma~\ref{lem.JBound}.
		
		
		
		For the term $\mathbf{II.2}$, we use the variational formula~\eqref{eq:VariationalFormula72} for the maximizer of $\nu^*(z+\cu, \ab_*(z+\cu_n) e_i)$ to obtain
		\begin{align*}
			\Er\Ll[\int_{z+\cu_n} \nabla V \cdot \a (e_i + \nabla \phi^*_{e_i, z+\cu_n}) \, \d \mu\Rr] = \Er\Ll[\int_{z+\cu_n} \nabla V \cdot \ab_*(z+\cu_n) e_i  \, \d \mu\Rr].
		\end{align*}
		Then we obtain that 
		\begin{equation}\label{eq.TermII2}
			\begin{split}
				\vert \mathbf{II.2}\vert &= \Ll\vert\sum_{i = 1}^d \sum_{z \in \Z_{m,n}}  (\partial_i \ub)_{z+\cu_n} \Er\Ll[\int_{z+\cu_n} \nabla V \cdot (\ab_*(z+\cu_n) - \ab) e_i \, \d \mu\Rr]\Rr\vert\\
				&\leq C 3^{-\alpha n} \norm{V}_{\acH^1(\cu_m)} \norm{\Ub}_{\acH^1(\cu_m)},
			\end{split}
		\end{equation}
		because \cite[Theorem~2.1]{bulk} implies that $\vert \ab_*(z+\cu_n) - \ab \vert \leq C3^{-\alpha n}$, as recalled in \eqref{eq:QuantConvergence}.
		
		\smallskip
		\textit{Step~2.3: choice of parameters.}
		Combining the estimates \eqref{eq.TermI}, \eqref{eq.TermII1}, and \eqref{eq.TermII2}, we obtain 
		\begin{multline}\label{eq.TwoScaleMixed}
			\Ll\vert \Er\Ll[\int_{\cu_m} \nabla V \cdot \a \nabla (W-U) \, \d \mu\Rr]\Rr\vert  \\
			\leq C \norm{V}_{\acH^1(\cu_m)} \Ll( \rho  \int_{\cu_m} 3^{2n} \vert\nabla \nabla \ub\vert^2  + 3^{-2\alpha n} \vert \nabla \ub\vert^2 \, \d m \Rr)^\frac{1}{2}.
		\end{multline}
		In order to avoid confusion, here we write down all the integration with respect to $\ub$ explicitly. Then we apply $H^1$ and $H^2$ estimates for $\ub$ defined in \eqref{eq.defub}
		\begin{align*}
			\rho  \int_{\cu_m} 3^{2n} \vert\nabla \nabla \ub\vert^2  + 3^{-2\alpha n} \vert \nabla \ub\vert^2 \, \d m &\leq  C (3^{2n} + 3^{2m-2\alpha n}) \rho  \int_{\cu_m} f^2 \, \d m\\
			&= C (3^{2n} + 3^{2m-2\alpha n}) \norm{F}^2_{\cL^2}.
		\end{align*} 
		The validity of the $H^2$ estimate for $\ub$ is asserted for general convex domains in \cite[Theorem~3.1.2.1]{grisvard1985elliptic} (notice that the constant $C(\Omega)$ appearing there only depends on the diameter of $\Omega$); in the particular case of cubes we can also more directly adapt the proof of \cite[Lemma~B.19]{AKMbook} to the case of a Dirichlet boundary condition (the extension of the function $f$ appearing in this proof can be taken as $f(x_1,\ldots, x_d) := \mathrm{sgn}(x_1) \cdots \mathrm{sgn}(x_d) f(|x_1|,\ldots, |x_d|)$ in this case). 
		With the choice of $n = \lfloor \frac{m}{1+\alpha} \rfloor$, we obtain that 
		\begin{align*}
			\Ll\vert \Er\Ll[\int_{\cu_m} \nabla V \cdot \a \nabla (W-U) \, \d \mu\Rr] \Rr\vert \leq C 3^{\frac{m}{1+\alpha}} \norm{V}_{\acH^1}(\cu_m) \norm{F}_{\cL^2(\cu_m)}. 
		\end{align*}
		Viewing $W$ as the sum of $\Ub$ and a linear combination of some $\cH^1_0(\cu_m)$ functions,  it also belongs to $\cH^1_0(\cu_m)$. We choose $V = W-U \in \cH^1_0(\cu_m)$ and use the Poincar\'e inequality \eqref{e.Poincare} and the uniform ellipticity condition in \eqref{e.ass.unif.ell} to obtain that
		\begin{align}\label{eq.TwoScaleH-1}
			\|W-U\|_{\acH^1(\cu_m)} 
			%+ \Er\Ll[\int_{\cu_m} \vert  \nabla (W-U) \vert^2 \, \d \mu\Rr]^{\frac{1}{2}}
			 \leq C 3^{\frac{m}{1+\alpha}} \norm{F}_{\cL^2(\cu_m)}. 
		\end{align}
%		The  estimate of $\norm{W-U}_{\cL^2}$ is an application of Poincar\'e's inequality and \eqref{eq.TwoScaleH-1}.
		
		
		
		\smallskip
		\textit{Step 3: estimate of $\norm{W-\Ub}_{\cL^2}$.} The $L^2$  estimate of $\norm{W-\Ub}^2_{\cL^2}$  can be done directly
		\begin{align*}
			\norm{W-\Ub}^2_{\cL^2} &=  \Er\Ll[\Ll(\sum_{i=1}^d \sum_{z \in \Z_{m,n}} (\partial_i \ub)_{z+\cu_n} \phi_{e_i, z+\cu_n}\Rr)^2\Rr]\\
			&\leq  C \sum_{i=1}^d \sum_{z \in \Z_{m,n}} (\partial_i \ub)^2_{z+\cu_n} \Er\Ll[ \phi^2_{e_i, z+\cu_n}\Rr] \\
			&\leq  C 3^{2n} \sum_{i=1}^d \sum_{z \in \Z_{m,n}} (\partial_i \ub)^2_{z+\cu_n}  \Er\Ll[ \int_{z+\cu_n} \vert \nabla \phi_{e_i, z+\cu_n}\vert^2 \, \d \mu\Rr]\\
			&\leq C 3^{2n} \norm{\Ub}^2_{\acH^1(\cu_m)}.
		\end{align*}
		Here in the second line we use the fact that $\phi_{e_i, z+\cu_n}$ is local function,  thus for $\dist(z, z') > 3^n$, the associated correctors are independent; see (1) of Lemma~\ref{lem.JBound}. Then we use the basic $\cH^1$ estimate of $\Ub$ and Poincar\'e's inequality to conclude
		\begin{align}\label{eq.TwoScaleL2}
			3^{-m} \norm{W-\Ub}_{\cL^2} \leq C 3^{n} \norm{F}_{\cL^2}.
		\end{align}
		Recalling the choice of $n = \lfloor \frac{m}{1+\alpha} \rfloor$, we obtain the desired estimate.
		\end{proof}
	\begin{remark}
		One can also study the homogenization of elliptic equation \eqref{eq.defDirichlet} without the assumption \eqref{eq.defF}, and this will be an interesting question for future work. %Surprisingly, the main challenge for general case seems in the step to fix the slope rather than homogenization, because the particles far away can also have correlation even for the homogenized solution. The present version \eqref{eq.defW} can largely cutoff such dependence of particles thanks to the linear statistic.
	\end{remark}
	
	\section{Quantitative estimate of two-point function and semigroup}
	\label{sec:4}
	In this section, we use a parabolic two-scale expansion to prove Theorem~\ref{t.two.point}. Recall that $P_t$ is the parabolic semigroup associated with $(\Ee^\a, \D(\Ee^\a))$ defined in \eqref{eq.defA}, and let $\Pb_t$ be the one associated to $(\Ee^\ab, \D(\Ee^\ab))$ defined similarly. Theorem~\ref{t.two.point} is also related to the homogenization of the parabolic semigroup.
	\begin{proposition}\label{prop.HomoSemigroup}
		There exists an exponent $\beta(d,\Lambda,\rho) > 0$ and  $C(d,\Lambda,\rho) < + \infty$ such that for every $f \in L^1(\Rd) \cap L^2(\Rd)$, its associated linear statistic $F(\mu) := \int_{\Rd} f \, \d \mu$ satisfies, for every $t > 0$, 
		\begin{align}\label{eq.HomoSemigroup}
			\norm{(P_t  - \Pb_t)F}_{\cL^2} \leq C t^{- \frac{\beta}{2}} \norm{f}_{L^2}.
		\end{align}
%		Here $\beta$ is the same exponent as that in  Theorem~\ref{t.two.point}.
	\end{proposition}
	
	
	
	
	
	Our proof essentially follows the outline of \cite{zhikov2006estimates}. % by Zhikov and Pastukhova. 
	In the proof, there are two main ingredients: the two-scale expansion and the regularization by the parabolic semigroup. 
	%In our setting, since we only implement the two-scale expansion on linear statistics, the step of regularization should be done more carefully. 
	Since $(\Ee^\ab, \D(\Ee^\ab))$ defines the diffusion of independent Brownian motions, a very important observation is the explicit formula of $\Pb_t F(\mu)$ for the linear statistic $F(\mu) = \int_{\Rd} f \, \d \mu$, which reads as
	\begin{align}\label{eq.PbHomo}
		\Pb_t F (\mu) = \int_{\Rd}  \Psi_t \star f \, \d \mu,
	\end{align} 
	where we recall that $(\Psi_t(x))_{t \ge 0, x \in \Rd}$ is the heat kernel defined in \eqref{e.def.heat.kernel}, and $\star$ denotes the spatial convolution in $\Rd$. 
	We use the following classical heat kernel estimates on $\Rd$ (see e.g.\ \cite[(1.17)-(1.18)]{zhikov2006estimates}): for every $k \in \N$, there exists a constant $C_k(d, \Lambda, \rho) < +\infty$ such that for every $g \in L^2$ and with $g_t := \Psi_t \star g$,
	\begin{align}\label{eq.heatkernelClassic}
		\forall 0 < \tau < t < +\infty,  \qquad  \int_\tau^t \norm{\nabla^k g_s}^2_{L^2} \, \d s &\leq C_k \norm{\nabla^{k-1} g_\tau}^2_{L^2},\\
		\nonumber \norm{\nabla^k g_t}_{L^2} &\leq \frac{ C_k}{(t-\tau)^{k/2}}\norm{g_\tau}_{L^2}.
	\end{align}
	
	
	\subsection{Two-scale expansion of parabolic equation}
	We first give a proof of homogenization in the special case of a rather smooth initial condition, using a two-scale expansion. For any fixed $g \in L^1(\Rd) \cap L^2(\Rd)$, we denote its linear statistic by $G := \int_\Rd g \, \d \mu - \rho \int_{\Rd} g \, \d m$, and $G_t := P_t G$, $\bar{G}_t := \Pb_t G$. We propose its associated two-scale expansion
	\begin{align}\label{eq.TwoScalePara}
		\tilde G_t := \bar{G}_t + \sum_{i = 1}^d \sum_{z \in \Z_n} (\partial_i g_t)_{z+\cu_n} \phi_{e_i, z+\cu_n},
	\end{align} 
	where $\Z_n$ is defined in~\eqref{e.Z.def} and $g_t := \Psi_t \star g$. Our first estimate also assumes that $g$ belongs to $H^2(\Rd)$. 
	\begin{proposition}\label{prop.TwoScalePara}
		Let $\alpha(d, \Lambda, \rho) > 0$ be as given by Lemma~\ref{lem.JBound}. There exists a constant $C(d, \Lambda, \rho) < +\infty$ such that for every $g \in H^2(\Rd) \cap L^1(\Rd)$, $n \in \N_+$, and $t >0$, its linear statistic and associated two-scale expansion defined above satisfy 
		\begin{multline}\label{eq.TwoScaleParaBound}
			\norm{\tilde G_t - \bar G_t}_{\cL^2} + \norm{\tilde G_t - G_t}_{\cL^2} + \Ll(\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s\Rr)^{\frac{1}{2}}\\
			\leq C \Ll(\Ll(3^{-\alpha n} + t^{-\frac{1}{8}}\Rr)\norm{g}_{L^2} + 3^n \norm{\nabla g}_{L^2} + 3^n t^{\frac{5}{8}}\norm{\nabla^2 g}_{L^2}\Rr).
		\end{multline}
	\end{proposition}
	\begin{proof}		
		The part $\norm{\tilde G_t - \bar G_t}_{\cL^2}$ can be bounded directly
		\begin{equation}\label{eq.TwoScaleParaBound1}
			\begin{split}
				\norm{\tilde G_t - \bar G_t}^2_{\cL^2} &=  \Er\Ll[\Ll(\sum_{i=1}^d \sum_{z \in \Z_n} (\partial_i g_t)_{z+\cu_n} \phi_{e_i, z+\cu_n}\Rr)^2\Rr]\\
				&\leq  C \sum_{i=1}^d \sum_{z \in \Z_n} (\partial_i g_t)^2_{z+\cu_n} \Er\Ll[ \phi^2_{e_i, z+\cu_n}\Rr] \\
				&\leq  C 3^{2n} \sum_{i=1}^d \sum_{z \in \Z_n} (\partial_i g_t)^2_{z+\cu_n}  \Er\Ll[ \int_{z+\cu_n} \vert \nabla \phi_{e_i, z+\cu_n}\vert^2 \, \d \mu\Rr]\\
				&\leq C 3^{2n} \norm{\nabla g_t}^2_{L^2}\\
				&\leq C 3^{2n} \norm{\nabla g}^2_{L^2}.
			\end{split}
		\end{equation}		
		In the first three lines, we use the independence of the approximate correctors in (1) of Lemma~\ref{lem.JBound}. In the last line, we use the contraction from the convolution $\Psi_t$.  
		
		The main part of the proof is the estimation of $(\tilde G_t - G_t)$. Using the parabolic equations associated to $G_s$ and $\bar G_s$, we have
		\begin{align*}
			\Ll(\partial_s - \frac{1}{2} \nabla \cdot \a \nabla\Rr) G_s = \Ll(\partial_s - \frac{1}{2} \nabla \cdot \ab \nabla\Rr) \bar G_s = 0,
		\end{align*}
		and thus
		\begin{align*}
			\Ll(\partial_s - \frac{1}{2} \nabla \cdot \a \nabla\Rr)(\tilde G_s - G_s) = \Ll(\partial_s - \frac{1}{2} \nabla \cdot \a \nabla\Rr)\tilde G_s - \Ll(\partial_s - \frac{1}{2} \nabla \cdot \ab \nabla\Rr) \bar G_s .
		\end{align*}
		We test this with $(\tilde G_s - G_s)$ and integrate over $[0,t]$ 
		\begin{equation}\label{eq.paraDecom}
			\begin{split}
				&\norm{\tilde G_t - G_t}^2_{\cL^2} + \int_0^t \Er\Ll[\int_{\Rd} \vert \nabla(\tilde G_s - G_s) \vert^2 \, \d \mu \Rr] \, \d s \\
				&\leq \norm{\tilde G_0 - G_0}^2_{\cL^2} \\
				& \qquad + \int_{0}^t \Er\Ll[\int_{\cu_m} \nabla (\tilde G_s - G_s) \cdot (\a \nabla \tilde G_s - \ab \nabla \bar G_s) \, \d \mu \Rr] \, \d s \\
				& \qquad + \int_{0}^t   \Er\Ll[ (\tilde G_s - G_s) \sum_{i=1}^d \sum_{z \in \Z_n} (\partial_i \partial_s g_s)_{z+\cu_n}\phi_{e_i, z+\cu_n}\Rr]\, \d s. 
			\end{split}
		\end{equation}
		It suffices to estimate the three terms on the {\rhs} to obtain \eqref{eq.TwoScaleParaBound}. The first term can be estimated similarly to $\norm{\tilde G_t - \bar G_t}^2_{\cL^2}$ in \eqref{eq.TwoScaleParaBound1} as the evolution has not started and $G_0  = \bar G_0$, so
		\begin{align}\label{eq.paraDecom1}
			\norm{\tilde G_0 - G_0}^2_{\cL^2} \leq C 3^{2n} \norm{\nabla g}^2_{L^2}.
		\end{align}
		The third term in \eqref{eq.paraDecom} only depends on the elementary $\cL^2$ bound, as $\partial_i \partial_s g_s$ gains enough decay from \eqref{eq.heatkernelClassic}
		\begin{equation}\label{eq.paraDecom3}
			\begin{split}
				&\Ll\vert \int_{0}^t   \Er\Ll[ (\tilde G_s - G_s) \sum_{i=1}^d \sum_{z \in \Z_n} (\partial_i \partial_s g_s)_{z+\cu_n}\phi_{e_i, z+\cu_n}\Rr]\, \d s \Rr\vert \\
				&\leq  \int_{0}^t \frac{t^{-\frac{5}{4}}}{2} \norm{\tilde G_s - G_s}^2_{\cL^2} +  \frac{t^{\frac{5}{4}}}{2} \norm{\sum_{i=1}^d \sum_{z \in \Z_n}(\partial_i \partial_s g_s)_{z+\cu_n}\phi_{e_i, z+\cu_n}}^2_{\cL^2} \, \d s \\
				&\le C  \int_{0}^t t^{-\frac{5}{4}} \Ll(\norm{\bar G_s}^2_{\cL^2} + \norm{G_s}^2_{\cL^2} + \sum_{z \in \Z_n} \norm{\sum_{i = 1}^d (\partial_i g_t)_{z+\cu_n} \phi_{e_i, z+\cu_n}}^2_{\cL^2}\Rr)  \\
                &  \qquad \qquad +  t^{\frac{5}{4}} \sum_{z \in \Z_n}\norm{\sum_{i=1}^d (\partial_i \partial_s g_s)_{z+\cu_n}\phi_{e_i, z+\cu_n}}^2_{\cL^2} \, \d s \\
				&\leq C\int_{0}^t t^{-\frac{5}{4}}  (\norm{g}^2_{L^2} + 3^{2n}\|\nabla g\|_{L^2}^2)  + 3^{2n} t^{\frac{5}{4}}  \norm{\nabla \partial_s g_s}^2_{L^2}\, \d s\\
				&\leq C \Ll(t^{-\frac{1}{4}} \norm{g}^2_{L^2} + t^{-\frac{1}{4}} 3^{2n}\|\nabla g\|_{L^2}^2 + 3^{2n} t^{\frac{5}{4}}  \norm{\nabla^2 g}^2_{L^2}\Rr).
			\end{split}
		\end{equation}
		Here in the second line, we apply Young's inequality  and choose the weight $t^{\frac{5}{4}}$ to better balance each term (it suffices to take $t^{\theta}$ with an exponent $\theta \in (1,2)$). In the third line, we apply the triangle inequality to $\norm{\tilde G_s - G_s}_{\cL^2}$ with the expression \eqref{eq.TwoScalePara} at first, then we apply the independence of the approximate correctors in Lemma~\ref{lem.JBound} to the term of type $\norm{\sum_{i=1}^d \sum_{z \in \Z_n}(\partial_i \partial_s g_s)_{z+\cu_n}\phi_{e_i, z+\cu_n}}^2_{\cL^2}$. From the fourth line to the fifth line, we use \eqref{eq.heatkernelClassic} and the identity $\partial_s g_s = \frac{1}{2} \nabla \cdot (\ab \nabla g_s)$.  
		
		Using the explicit expression of $\bar G_s$ in \eqref{eq.PbHomo}, the second term in \eqref{eq.paraDecom} can be treated like the elliptic case \eqref{eq.ellipticDecom} by the following decomposition
		\begin{align}\label{eq.paraDecom2Decom}
			\int_{0}^t \Er\Ll[\int_{\Rd} \nabla (\tilde G_s - G_s) \cdot (\a \nabla \tilde G_s - \ab \nabla \bar G_s) \, \d \mu \Rr] \, \d s = \mathbf{I'} + \mathbf{II'},
		\end{align}
		where the two terms are respectively
		\begin{align*}
			\mathbf{I'} & := \int_{0}^t \sum_{z \in \Z_{n}}  \Er\Ll[\int_{z+\cu_n} \nabla (\tilde G_s - G_s) \cdot (\a-\ab)(\nabla g_s - (\nabla g_s)_{z+\cu_n}) \, \d \mu\Rr] \, \d s, \\
			\mathbf{II'} & := \int_{0}^t  \sum_{i = 1}^d \sum_{z \in \Z_{n}}  (\partial_i g_s)_{z+\cu_n} \Er\Ll[\int_{z+\cu_n} \nabla (\tilde G_s - G_s) \cdot \Ll(\a(e_i + \nabla \phi_{e_i, z+\cu_n}) - \ab e_i\Rr) \, \d \mu\Rr] \, \d s.
		\end{align*} 
		For the term $\mathbf{I'}$, the main issue is to use the Poincar\'e inequality to fix the slope, and the proof is similar to \eqref{eq.TermI}
		\begin{equation}\label{eq.paraDecom2Term1}
			\begin{split}
				\vert \mathbf{I'} \vert &\leq C  3^n\Ll(\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s\Rr)^{\frac{1}{2}} \Ll(\int_{0}^t \sum_{z \in \Z_{n}} \int_{z+\cu_n} \vert\nabla^2 g_s\vert^2 \, \d m  \, \d s\Rr)^{\frac{1}{2}}\\
				&=C 3^n\Ll(\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s\Rr)^{\frac{1}{2}} \Ll(\int_{0}^t \norm{\nabla^2 g_s}^2_{L^2}  \, \d s\Rr)^{\frac{1}{2}}\\
				&\leq C 3^n\Ll(\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s\Rr)^{\frac{1}{2}}  \norm{\nabla g}_{L^2}  \\
				&\leq \frac{1}{4}\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s + C 3^{2n}\norm{\nabla g}^2_{L^2}. 
			\end{split}
		\end{equation}
		In the third line, we apply the classical heat kernel decay \eqref{eq.heatkernelClassic}. In the fourth line, the Young's inequality is used, and the first term there can be compensated by part of the {\lhs} in \eqref{eq.paraDecom}. 
		
		For the term $\mathbf{II'}$, we apply the flux replacement argument like \eqref{eq.TermII1} and \eqref{eq.TermII2}
		\begin{equation}\label{eq.paraDecom2Term2}
			\begin{split}
				\vert \mathbf{II'}\vert &\leq C \int_0^t 3^{-\alpha n} \norm{\nabla g_s}_{L^2} \Er^{\frac{1}{2}}\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr]  \, \d s \\
				&\leq \int_0^t C 3^{- 2\alpha n} \norm{\nabla g_s}^2_{L^2}  \, \d s + \frac{1}{4}\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s\\
				&\leq C 3^{-2\alpha n} \norm{g}^2_{L^2}  + \frac{1}{4}\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s. 
			\end{split}
		\end{equation}
		We apply Young's inequality in the second line, and the classical heat kernel decay \eqref{eq.heatkernelClassic} in the third line. We put \eqref{eq.paraDecom2Term1}, \eqref{eq.paraDecom2Term2} back to \eqref{eq.paraDecom2Decom} and get 
		\begin{multline}\label{eq.paraDecom2}
			\Ll\vert \int_{0}^t \Er\Ll[\int_{\Rd} \nabla (\tilde G_s - G_s) \cdot (\a \nabla \tilde G_s - \ab \nabla \bar G_s) \, \d \mu \Rr]\Rr\vert \\
			\leq \frac{1}{2}\int_0^t \Er\Ll[\int_{\Rd} \vert \nabla (\tilde G_s - G_s)\vert^2 \, \d \mu \Rr] \, \d s + C 3^{-2\alpha n} \norm{g}^2_{L^2}  +  C 3^{2n}\norm{\nabla g}^2_{L^2}.
		\end{multline}
		Finally, we put \eqref{eq.paraDecom1}, \eqref{eq.paraDecom3}, \eqref{eq.paraDecom2} back to \eqref{eq.paraDecom} and obtain the desired result.
	\end{proof}
	
	
	\begin{remark}
		The result above is already enough for us to prove \eqref{e.two.pointScale} with the scaling and a good test function $g \in H^2 \cap L^1$. Let us define 
		\begin{align*}
			g^N(x) := N^{-\frac{d}{2}} g\Ll(\frac{x}{N}\Rr), \qquad G^N := \int_\Rd g \, \d \mu - \rho \int_{\Rd} g \d m.
		\end{align*}
		This scaling gives us 
		\begin{align*}
			\norm{G^N}^2_{\cL^2} = \rho \norm{g^N}^2_{L^2} = \rho \norm{g}^2_{L^2}.
		\end{align*}
		Then we calculate the density field correlation 
		\begin{align*}
			&\Ll\vert\Er \EEbmu\Ll[Y^N_t(g) Y^N_0(g)\Rr] - \rho \int_{\Rd \times \Rd} g(x)\Psi_{t}(x-y)g(y)\, \d x  \d y \Rr\vert\\
			&=  \Ll\vert \Er\Ll[G^N (P_{N^2t}-\Pb_{N^2t})G^N\Rr] \Rr\vert\\
			&\leq \norm{G^N}_{\cL^2} \norm{(P_{N^2t}-\Pb_{N^2t})G^N}_{\cL^2}\\
			&\leq C\Ll(3^{-\alpha n} + (N^2t)^{-\frac{1}{8}}\Rr)\norm{g}_{L^2} + 3^n N^{-1} \norm{\nabla g}_{L^2}+ 3^n (N^2t)^{\frac{5}{8}}N^{-2} \norm{\nabla^2 g}_{L^2}.
		\end{align*}
		With the scaling, we gain an extra factor $N^{-1}$ for the gradient. It suffices to choose $1 \ll 3^n \ll N^{\frac{3}{4}}$, for example, $3^n \simeq N^{\frac{1}{2}}$, then the error terms vanish when $N \to \infty$ with an explicit rate. 
	\end{remark}
	
	\subsection{Regularization effect}
	In this part, we will show that the diffusive regularization effect of the parabolic semigroup can help weaken the assumption $g \in H^2$ in Proposition~\ref{prop.TwoScalePara}.  
	
	\begin{lemma}\label{lem.regularization}
		For any $\tau, t > 0$, we have the following estimates
		\begin{align}\label{eq.regularization}
			\norm{\Pb_\tau P_t - P_t}_{\cL^2 \to \cL^2}  \leq \sqrt{\frac{\Lambda\tau}{t}} \quad \text{ and } \quad \norm{\Pb_\tau \Pb_t - \Pb_t}_{\cL^2 \to \cL^2} \leq \sqrt{\frac{\Lambda\tau}{t}} .
		\end{align}
	\end{lemma}
	\begin{proof}
		Using the parabolic equation \eqref{eq.semigroupA} associated to $\Pb_t$, for any $\tilde F \in \cL^2$
		\begin{align*}
			\norm{\Pb_\tau \tilde F}^2_{\cL^2} - \norm{\tilde F}^2_{\cL^2} + \int_0^\tau \Er\Ll[\int_{\Rd}  \nabla \Pb_s \tilde F \cdot \ab \nabla \Pb_s \tilde F \, \d \mu \Rr] \, \d s = 0,
		\end{align*}
		which implies that 
		\begin{align*}
			\norm{\Pb_\tau \tilde F}^2_{\cL^2} &\geq \norm{\tilde F}^2_{\cL^2} - \tau \Er\Ll[\int_{\Rd}   \nabla  \tilde F \cdot \ab  \nabla  \tilde F \, \d \mu \Rr] \\
			&\geq  \norm{\tilde F}^2_{\cL^2} - \Lambda \tau \Er\Ll[\int_{\Rd}  \vert \nabla  \tilde F \vert^2 \, \d \mu \Rr].
		\end{align*}
		Here we use the decay of the Dirichlet energy of the parabolic semigroup in Proposition~\ref{prop.ElementSemigroup}. Then we insert $\tilde{F} = P_t F$, and obtain 
		\begin{align}\label{eq.regularization1}
			\norm{\Pb_\tau P_t F}^2_{\cL^2} &\geq \norm{P_t F}^2_{\cL^2} - \Lambda \tau \Er\Ll[\int_{\Rd}  \vert \nabla  P_t F  \vert^2 \, \d \mu \Rr]\\
			\nonumber &\geq \norm{P_t F}^2_{\cL^2} - \frac{\Lambda \tau}{t} \norm{F}^2_{\cL^2}.
		\end{align}
		In the last line, we also use the semigroup property associated to $P_t$ in \eqref{eq.L2Decay}. This estimate results in that of $\norm{\Pb_\tau P_t - P_t}_{\cL^2 \to \cL^2}$
		\begin{align*}
			\norm{(\Pb_\tau P_t - P_t) F}^2_{\cL^2} &= \norm{\Pb_\tau P_t F}^2_{\cL^2} + \norm{P_t F}^2_{\cL^2} - 2 \norm{\Pb_{\frac{\tau}{2}} P_t F}^2_{\cL^2}\\
			&\leq 2 \norm{P_t F}^2_{\cL^2} - 2\Ll(\norm{P_t F}^2_{\cL^2} - \frac{\Lambda \tau}{2t}\norm{F}^2_{\cL^2} \Rr)\\
			&\leq \frac{\Lambda \tau}{t}\norm{F}^2_{\cL^2}.
		\end{align*}
		Here in the second line, we apply $\norm{\Pb_\tau P_t F}^2_{\cL^2} \leq \norm{P_t F}^2_{\cL^2}$ and the estimate \eqref{eq.regularization1} to $\norm{\Pb_{\frac{\tau}{2}} P_t F}^2_{\cL^2}$. A similar argument also works for $\norm{\Pb_\tau \Pb_t - \Pb_t}_{\cL^2 \to \cL^2}$.
	\end{proof}
	
	\begin{proof}[Proof of Theorem~\ref{t.two.point}]
		For two functions $f, g$ as in the statement, we define 
		\begin{align*}
			F(\mu) := Y_0(f) = \int_{\Rd} f \, \d \mu - \rho \int_{\Rd} f \, \d m,\\
			G(\mu) := Y_0(g) = \int_{\Rd} g \, \d \mu - \rho \int_{\Rd} g \, \d m.
		\end{align*}
		As $\Pr$ is a stationary measure for the process $(\bmu_t)_{t \geq 0}$, it suffices to treat the case $s = 0$, and \eqref{e.two.point} is equivalent to estimating $\Er[G (P_t -\Pb_t) F]$. We apply the regularization at first. Recall that we denote by $\bracket{\cdot, \cdot}_{\cL^2}$ the inner product of $\cL^2$ space
		\begin{align*}
			&\Ll\vert \bracket{G, (P_t -\Pb_t) F}_{\cL^2} \Rr\vert \\
			&\leq \Ll\vert \bracket{G, \Pb_\tau (P_t -\Pb_t) F}_{\cL^2} \Rr\vert + \norm{G}_{\cL^2} \norm{\Pb_\tau (P_t -\Pb_t) F - (P_t -\Pb_t) F}_{\cL^2}\\
			&\leq \Ll\vert \bracket{G, \Pb_\tau (P_t -\Pb_t) F}_{\cL^2} \Rr\vert + 2\sqrt{\frac{\Lambda\tau}{t}} \norm{G}_{\cL^2}  \norm{F}_{\cL^2}.
		\end{align*}
		In the second line, we apply Lemma~\ref{lem.regularization} to the second term. When $\tau \ll t$, the error paid is very small and it suffices to consider the regularized version. Then we use the self-adjoint property of the operator 
		\begin{align*}
			\bracket{G, \Pb_\tau (P_t -\Pb_t) F}_{\cL^2} = \bracket{ (P_t -\Pb_t) (\Pb_\tau G),   F}_{\cL^2}.
		\end{align*}
		The term $\Pb_\tau G$ gains some more regularity, and it still has an explicit expression 
		\begin{align*}
			\Pb_\tau G (\mu) = \int_{\Rd} g_\tau \, \d \mu - \rho \int_{\Rd} g_\tau \, \d m,
		\end{align*}
		with $g_\tau = \Psi_\tau \star g \in H^2(\Rd) \cap L^1(\Rd)$. Thus, we apply Proposition~\ref{prop.TwoScalePara} to it
		\begin{align*}
			\Ll\vert \bracket{G, \Pb_\tau (P_t -\Pb_t) F}_{\cL^2} \Rr\vert &\leq \norm{ (P_t -\Pb_t) (\Pb_\tau G)}_{\cL^2} \norm{F}_{\cL^2} \\
			&\leq C\Ll(\Ll(3^{-\alpha n} + t^{-\frac{1}{8}}\Rr)\norm{g_\tau}_{L^2} + 3^n \norm{\nabla g_\tau}_{L^2} + 3^n t^{\frac{5}{8}}\norm{\nabla^2 g_\tau}_{L^2}\Rr) \norm{f}_{L^2}\\
			&\leq C\Ll(3^{-\alpha n} + t^{-\frac{1}{8}} +  3^n \tau^{-\frac{1}{2}} + 3^n t^{\frac{5}{8}} \tau^{-1}\Rr)\norm{g}_{L^2} \norm{f}_{L^2}.
		\end{align*}
		From the second line to the third line, we apply the decay of classical heat kernel \eqref{eq.heatkernelClassic}. Combining all the estimates above, we obtain that 
		\begin{multline}\label{eq.ErrorSource}
			\Ll\vert \bracket{G, (P_t -\Pb_t) F}_{\cL^2} \Rr\vert \\
			\leq  C\Ll(3^{-\alpha n} + t^{-\frac{1}{8}} +  3^n \tau^{-\frac{1}{2}} + 3^n t^{\frac{5}{8}} \tau^{-1} + (\tau/t)^{\frac{1}{2}}\Rr)\norm{g}_{L^2} \norm{f}_{L^2}.
		\end{multline}
		With a choice of mesoscopic scales $1 \ll 3^n \ll t^{\frac{5}{8}} \ll \tau \ll t$, for example $\tau = t^{\frac{3}{4}}, 3^n \simeq t^{\frac{1}{16}}$ we obtain the desired result \eqref{e.two.point} with a parameter $\beta = \frac{\min\{\alpha, 1\}}{16}$.
	\end{proof}
	\begin{remark}
		The errors in \eqref{eq.ErrorSource} have clear interpretation as the regularization, the homogenization and the price to fix the local slope to homogenization.
	\end{remark}
	\begin{proof}[Proof of Proposition~\ref{prop.HomoSemigroup}] Since both $P_t$ and $\Pb_t$ are conservative semigroups, we denote by $\hat{F}:= F - \rho \int_{\Rd} f \, \d m$ and have $(P_t - \Pb_t )F = (P_t - \Pb_t )\hat{F}$. Then using the self-adjoint property of the semigroup, we have
		\begin{align*}
			\norm{(P_t - \Pb_t )F}^2_{\cL^2} = \bracket{\hat{F}, P_{2t} \hat{F}}_{\cL^2} + \bracket{\hat{F}, \Pb_{2t} \hat{F}}_{\cL^2} - 2\bracket{\Pb_t \hat{F}, P_t \hat{F}}_{\cL^2}.
		\end{align*}
		Theorem~\ref{t.two.point} applies to the first and third terms, then we obtain the estimate of the semigroup. 
	\end{proof}
	
	\section{Convergence rate of Green--Kubo formulas}
	\label{sec:5}
	In this section, we start by giving more details on how to make sense of the integral in the formula \eqref{eq.quanti.gk}, and then prove Theorem~\ref{t.quanti.gk}.
%	\subsection{Well-definedness of Green--Kubo formula}\label{subsec.WellDefined}
	We follow the discussion in Section~\ref{subsec.Dirichletform} to construct the semigroup $P_t^{(m)}$. We recall that $\cL^2(\cu_m)=L^2(\mmd(\Rd), \mcl F_{\cu_m}, \Pr)$ is the space of square-integrable $\mcl F_{\cu_m}$-measurable functions, and define the Dirichlet form $\Ee^{\a}_{\cu_m}$ by
	\begin{align}\label{eq.DirichletCube}
		\Ll\{\begin{aligned}
			& \D(\Ee^{\a}_{\cu_m}) := \cH^1_0(\cu_m),\\
			& \displaystyle{\forall u,v \in \D(\Ee^{\a}_{\cu_m}), \quad \Ee^{\a}_{\cu_m}(u,v) := \Er\Ll[\int_{\cu_m} \frac{1}{2}\nabla u \cdot \a \nabla v \, \d \mu \Rr].}
		\end{aligned}\Rr.
	\end{align}
	The Dirichlet form $(\Ee^{\a}_{\cu_m}, \D(\Ee^{\a}_{\cu_m}))$ is also associated to a unique self-adjoint positive semi-definite operator $-A^{(m)}$
	\begin{align}\label{eq.defAm}
		\Ll\{\begin{array}{ll}
			\D(\Ee^{\a}_{\cu_m}) & = \D(\sqrt{-A^{(m)}}), \\
			\Ee^{\a}_{\cu_m}(u,v) & = \bracket{\sqrt{-A^{(m)}}u, \sqrt{-A^{(m)}}v}_{\cL^2}.
		\end{array}\Rr.
	\end{align}
	We define the semigroup $P_t^{(m)} := e^{t A^{(m)}}$. It satisfies similar properties as in Proposition~\ref{prop.ElementSemigroup}, and in particular, for every $u \in \cL^2(\cu_m)$ and $t > 0$, we have $P_t^{(m)} u \in \cH^1_0(\cu_m)$.

	By the definition \eqref{eq.defFpm}, the functional $F_{p,m}$ is an element in $\cH^{-1}(\cu_m)$. We write $\bracket{\cdot, \cdot}_{\cH^{-1}, \cH^1_0}$ as shorthand for $\bracket{\cdot, \cdot}_{\cH^{-1}(\cu_m), \cH^1_0(\cu_m)}$, since $m$ is kept fixed throughout this section. For each $t > 0$, we define $P_t^{(m)} (F_{p,m}) \in \cL^2(\cu_m)$ by duality as the unique function in $\cL^2(\cu_m)$ such that 
	\begin{align}\label{eq.defPFpm}
		\forall f \in \cL^2(\cu_m), \qquad \bracket{P_t^{(m)} (F_{p,m}), f}_{\cL^2} = \bracket{F_{p,m}, P_t^{(m)} f}_{\cH^{-1}, \cH^1_0}.
	\end{align}
	This definition of $P_t^{(m)}(F_{p,m})$ relies on the Riesz representation theorem and on the observation that for each $t > 0$, the mapping $f \mapsto  \bracket{F_{p,m}, P_t^{(m)} f}_{\cH^{-1}, \cH^1_0}$ is a continuous linear form on $\cL^2(\cu_m)$, since by \eqref{eq.L2Decay},
	\begin{equation}\label{eq.L2FpmBound}
		\begin{split}
			\Ll\vert \bracket{F_{p,m}, P_t^{(m)} f}_{\cH^{-1}, \cH^1_0} \Rr\vert &\leq \norm{F_{p,m}}_{\cH^{-1}} \norm{P_t^{(m)} f}_{\cH^1} \\
			& \leq  t^{-\frac{1}{2}}\norm{F_{p,m}}_{\cH^{-1}} \norm{f}_{\cL^2}.
		\end{split}
	\end{equation}

	
	The next proposition collects additional properties satisfied by $F_{p,m}$ and $P_t^{(m)} (F_{p,m})$.
	\begin{proposition}\label{prop.Fpm}
		The following properties are valid.
		\begin{enumerate}
			\item With $\phi_{p,\cu_m}$ defined in \eqref{eq.corrector}, we have for every $ v$ in $\cH^1_0(\cu_m)$ that
			\begin{align}\label{eq.FpmCorrector}
				\bracket{F_{p,m}, v}_{\cH^{-1}, \cH^1_0} = \Ee^{\a}_{\cu_m}(\phi_{p,\cu_m}, v),
			\end{align}
			and 
			\begin{align}\label{eq.FpmContinuity}
				\lim_{t\to 0}\bracket{P_t^{(m)} (F_{p,m}), v}_{\cL^2} = \bracket{F_{p,m}, v}_{\cH^{-1}, \cH^1_0}. 
			\end{align}
			\item For every $t,s >0$, we have $P_{t+s}^{(m)} (F_{p,m}) = P_{t}^{(m)}( P_{s}^{(m)} (F_{p,m}))$.
			\item There exists $C(d) < +\infty$ such that for every $t >0$, we have that $P_t^{(m)} (F_{p,m})$ belongs to $\cH^1_0(\cu_m)$, is centered, and satisfies
			\begin{align}
				\norm{P_t^{(m)} (F_{p,m})}_{\cL^2} &\leq 4 \exp(-3^{-2m}t/C)  t^{-\frac{1}{2}}\norm{F}_{\cH^{-1}}, \label{eq.FpmL2Decay}\\
				\norm{P_t^{(m)} (F_{p,m})}_{\acH^1} &\leq 4 \exp(-3^{-2m}t/C)  t^{-1} \norm{F}_{\cH^{-1}}. \label{eq.FpmH1Decay}
			\end{align}
		\end{enumerate}
	\end{proposition}
	\begin{proof} 
            (1). Equation \eqref{eq.FpmCorrector} is a consequence of \eqref{eq.corrector} and \eqref{eq:VariationalFormula72}, which we combine with the definition of \eqref{eq.defFpm} and \eqref{eq.L2FpmBound} to obtain 
			\begin{align*}
				\bracket{P_t^{(m)} (F_{p,m}), v}_{\cL^2} = \bracket{F_{p,m}, P_t^{(m)} v}_{\cH^{-1}, \cH^1_0} = \Ee^{\a}_{\cu_m}(\phi_{p,\cu_m}, P_t^{(m)} v).
			\end{align*}
			The continuity of the Dirichlet energy in (2) of Proposition~\ref{prop.ElementSemigroup} implies that
			\begin{multline*}
				\lim_{t \to 0} \bracket{P_t^{(m)} (F_{p,m}), v}_{\cL^2} = \lim_{t \to 0} \Ee^{\a}_{\cu_m}(\phi_{p,\cu_m}, P_t^{(m)} v) \\
				=  \Ee^{\a}_{\cu_m}(\phi_{p,\cu_m}, v) = \bracket{F_{p,m}, v}_{\cH^{-1}, \cH^1_0}.
			\end{multline*} 
   
			(2). We test $P_{t}^{(m)}( P_{s}^{(m)} (F_{p,m}))$ with a function $f \in \cL^2(\cu_m)$
			\begin{align*}
				\bracket{P_{t}^{(m)}( P_{s}^{(m)} (F_{p,m})), f}_{\cL^2} &= \bracket{ P_{s}^{(m)} (F_{p,m}), P_{t}^{(m)} f}_{\cL^2}\\
				&= \bracket{F_{p,m}, P_{s}^{(m)} (P_{t}^{(m)} f)}_{\cL^2}\\
				&= \bracket{F_{p,m}, P_{t+s}^{(m)} f}_{\cL^2}.
			\end{align*}
			In the first line, we use the reversibility of $P_{t}^{(m)}$ with respect to $\Pr$, and then we apply the definition \eqref{eq.defFpm} from the first line to the second line. The result in the third line follows the semigroup property of $P_{t+s}^{(m)}$. Since it satisfies the characterization of $P_{t+s}^{(m)} (F_{p,m})$, the two quantities coincide. 
   
			(3). Using the property in (2), we have $P_t^{(m)} (F_{p,m}) = P_{t/2}^{(m)} \Ll(P_{t/2}^{(m)} (F_{p,m})\Rr)$, so the first semigroup $P_{t/2}^{(m)}$ maps the function from $\cL^2(\cu_m)$ to $\cH^1(\cu_m)$. $P_t^{(m)} (F_{p,m})$ is centered because we test \eqref{eq.defFpm} with constant function $f \equiv 1$ and it will give $0$ on the {\rhs}. Once we know the function is centered, we can apply the spectral gap inequality \eqref{e.Poincare} that
			\begin{align*}
				\frac{\d}{\d t} \norm{P_t^{(m)} (F_{p,m})}_{\cL^2}^2 &= - 2 \Ee^{\a}_{\cu_m}(P_t^{(m)} (F_{p,m}),P_t^{(m)} (F_{p,m}))\\ &\leq -\frac{3^{-2m}}{C} \norm{P_t^{(m)} (F_{p,m})}_{\cL^2}^2.
			\end{align*}
			This results in that $\norm{P_t^{(m)} (F_{p,m})}_{\cL^2}^2 \leq \exp(-3^{-2m}t/(2C)) \norm{P_{t/2}^{(m)} (F_{p,m})}_{\cL^2}^2$. Then we apply the estimate \eqref{eq.L2FpmBound} for $\norm{P_{t/2}^{(m)} (F_{p,m})}_{\cL^2}^2$ and obtain \eqref{eq.FpmL2Decay}. One more factor $t^{-1/2}$ can be gained when applying \eqref{eq.L2Decay} to the $\acH^1(\cu_m)$ estimate, which yields \eqref{eq.FpmH1Decay}.
	\end{proof}

	Using the properties above, we can justify that the integral in  \eqref{eq.quanti.gk} is well-defined. We write $\ell_{p,m}:= \int_{\cu_m} p \cdot x \, \d \mu(x)$ as the affine function in $\cu_m$.
	\begin{lemma}\label{lem.CurrentCorrelation}
		For every $m \in \N$ and $\lambda \geq 0$, the integral of the dynamic current-current correlation  
		\begin{multline}\label{eq.defCurrentCorrelation}
			\int_0^{+\infty}  e^{-\lambda t}  \la F_{p,m}, P_t^{(m)}(F_{p,m}) \rah \, \d t \\
			:= \lim_{\epsilon \to 0} \lim_{T \to \infty} \int_\epsilon^{T}  e^{-\lambda t}  \la F_{p,m}, P_t^{(m)}(F_{p,m}) \rah \, \d t,
		\end{multline}
		is well-defined. It is equal to 
        \begin{align}\label{eq.ValueCurrentCorrelation}
            \Er\Ll[\int_{\cu_m} \Ll(\frac{1}{2} p \cdot \a p -  \frac{1}{2} p \cdot \a \nabla U_\lambda\Rr) \, \d \mu \Rr],
        \end{align}
        where $U_\lambda$ is the unique centered solution in $\ell_{p,m} + \cH^1_0(\cu_m)$ such that   
		\begin{align}\label{eq.defUlambda}
			\forall v \in \cH^1_0(\cu_m), \qquad \lambda \bracket{U_{\lambda}, v}_{\cL^2} + \Ee^{\a}_{\cu_m}(U_{\lambda}, v) = \lambda \bracket{\ell_{p,m}, v}_{\cL^2}.
		\end{align} 
	\end{lemma}
	\begin{proof}
		Thanks to the exponential decay \eqref{eq.FpmL2Decay} and \eqref{eq.FpmH1Decay}, the integral 
		\begin{align*}
			V_{\lambda, \epsilon} := \lim_{T \to \infty} \int_\epsilon^{T} e^{-\lambda t} P_t^{(m)}(F_{p,m}) \, \d t,
		\end{align*}
		converges in $\cH^1_0(\cu_m)$, and the limit satisfies
		\begin{align}\label{eq.defVlambda}
			\forall v \in \cH^1_0(\cu_m), \qquad \lambda \bracket{V_{\lambda, \epsilon}, v}_{\cL^2} + \Ee^{\a}_{\cu_m}(V_{\lambda, \epsilon}, v) = \bracket{e^{-\lambda \epsilon}P_\epsilon^{(m)}(F_{p,m}), v}_{\cL^2}.
		\end{align}
		To see this, we test $V_{\lambda, \epsilon}$ with the {\lhs} of this equation, and apply \eqref{eq.defAm} 
		\begin{align*}
			\lambda \bracket{V_{\lambda, \epsilon}, v}_{\cL^2} + \Ee^{\a}_{\cu_m}(V_{\lambda, \epsilon}, v) &= \bracket{\int_\epsilon^\infty  (\lambda - A^{(m)}) e^{-t (\lambda - A^{(m)})}(F_{p,m}) \, \d t, v}_{\cL^2} \\
			&= \bracket{e^{-\lambda \epsilon}P_\epsilon^{(m)}(F_{p,m}), v}_{\cL^2}.
		\end{align*}		
		We test the equation \eqref{eq.defVlambda} with $V_{\lambda, \epsilon}$ and apply Young's inequality (for the case $\lambda=0$, we also need \eqref{e.Poincare}), then obtain that 
		\begin{align}\label{eq.VUniformBound}
			\sup_{\epsilon > 0}\norm{V_{\lambda, \epsilon}}_{\acH^1(\cu_m)} \leq C \norm{F_{p,m}}_{\cH^{-1}(\cu_m)}.
		\end{align}
		Therefore, the sequence admits at least one weak limit in the space $\cH^1_0(\cu_m)$. Moreover, when we pass $\epsilon \to 0$ in \eqref{eq.defVlambda} and apply \eqref{eq.FpmContinuity}, the weak limit $V_{\lambda}$ is unique which is characterized by 
		\begin{align}\label{eq.defV0}
			\forall v \in \cH^1_0(\cu_m), \qquad \lambda \bracket{V_{\lambda}, v}_{\cL^2} + \Ee^{\a}_{\cu_m}(V_{\lambda}, v) = \bracket{F_{p,m}, v}_{\cH^{-1}, \cH^1_0}.
		\end{align}
		Now we test $V_{\lambda, \epsilon}$ with $F_{p,m}$
		\begin{equation}\label{eq.DynamicLimit}
			\begin{split}
				\lim_{\epsilon \to 0} \lim_{T \to \infty} \int_\epsilon^{T}  e^{-\lambda t}  \la F_{p,m}, P_t^{(m)}(F_{p,m}) \rah \, \d t 
				&=\lim_{\epsilon \to 0} \bracket{F_{p,m}, V_{\lambda, \epsilon}}_{\cH^{-1}, \cH^1_0}\\
				&=\lim_{\epsilon \to 0} \Ee^{\a}_{\cu_m}(\phi_{p, \cu_m}, V_{\lambda, \epsilon})\\
				&= \Ee^{\a}_{\cu_m}(\phi_{p, \cu_m}, V_{\lambda})\\
				&= \bracket{F_{p,m},  V_{\lambda}}_{\cH^{-1}, \cH^1_0}.
			\end{split}
		\end{equation}
		Here we apply \eqref{eq.FpmCorrector} in the second line, and then pass to the limit using $V_{\lambda, \epsilon} \stackrel{\epsilon \to 0}{\rightharpoonup} V_{\lambda}$ in $\cH^1_0(\cu_m)$ from the second line to the third line. This justifies that \eqref{eq.defCurrentCorrelation} is well-defined. Finally, we define $U_\lambda := \ell_{p,m} + V_\lambda$, then \eqref{eq.defV0} becomes \eqref{eq.defUlambda} and we obtain the limit of the dynamic current correlation \eqref{eq.ValueCurrentCorrelation} using $U_\lambda$.
	\end{proof}

	
	We are now ready to prove the quantitative Green--Kubo formula in Theorem~\ref{t.quanti.gk}.

	\begin{proof}[Proof of Theorem~\ref{t.quanti.gk}]
		As discussed in Lemma~\ref{lem.CurrentCorrelation}, we use the stationarity  \eqref{e.def.amu} and insert \eqref{eq.ValueCurrentCorrelation} in the {\lhs} of \eqref{eq.quanti.gk} and it suffices to estimate 
		\begin{align}\label{eq.GreenKuboReduced} 
			\Ll\vert \frac{1}{2} p \cdot \ab p - \frac{1}{\rho |\cu_m|} \E_\rho \Ll[\int_{\cu_m} \frac{1}{2} p \cdot \a \nabla U_\lambda \, \d \mu \Rr]\Rr\vert,
		\end{align}
		with $U_\lambda$ defined in \eqref{eq.defUlambda}. We notice that the homogenized solution $\Ub \in \ell_{p,m} + \cH^1_0(\cu_m)$ solving  
		\begin{align}\label{eq.GKElliptic2}
			(\lambda - \nabla \cdot \ab \nabla )\Ub =  \lambda \ell_{p,m},
		\end{align}
		is exactly $\ell_{p,m}$ using a similar argument in \eqref{eq.IPP}. Then the main idea is to replace $U_\lambda$ in \eqref{eq.GreenKuboReduced} by its two-scale expansion to get the estimate.
		
		\smallskip
		\textit{Step~1: two-scale expansion.} We analyze the structure of $U_\lambda$ using the two-scale expansion technique \eqref{eq.defW} from Section~\ref{sec:3}, setting
		\begin{align}\label{eq.GKdefW}
			W := \ell_{p,m} + \sum_{i=1}^d \sum_{z \in \Z_{m,n}} p_i \phi_{e_i, z+\cu_n} = \ell_{p,m} + \sum_{z \in \Z_{m,n}} \phi_{p, z+\cu_n}. 
		\end{align}
		Here $p_i$ is the $i$-th coordinate of the vector $p$, and we use the linearity of the corrector as a function of $p$. Without loss of generality, we suppose $\vert p \vert = 1$. We can establish the following equation using \eqref{eq.defUlambda} and \eqref{eq.GKElliptic2} 
		\begin{align*}
			(\lambda - \nabla \cdot \a \nabla )(W - U_\lambda) &= 	(\lambda - \nabla \cdot \a \nabla )W - 	(\lambda - \nabla \cdot \ab \nabla )\Ub\\
			&= \lambda (W-\Ub) - \nabla \cdot (\a \nabla W - \ab \nabla \Ub).
		\end{align*}
		We test the equation above with $(W-U_\lambda)$ as their affine parts cancel and ${(W-U_\lambda)} \in \cH^1_0(\cu_m)$. This gives us 
		\begin{multline}\label{eq.GKEllipticDecom}
			\lambda \norm{W-U_\lambda}^2_{\cL^2} + \Er \Ll[ \int_{\cu_m} \nabla (W - U_\lambda) \cdot \a \nabla (W - U_\lambda) \, \d \mu\Rr]   \\
			\leq \lambda \vert \Er[(W-U_\lambda)(W-\Ub)]\vert + \Ll\vert \Er\Ll[\int_{\cu_m} \nabla (W-U_\lambda) \cdot  (\a \nabla W - \ab \nabla \Ub) \, \d \mu\Rr]\Rr\vert. 
		\end{multline}
		Its first term can be estimated by following Step~3 in the proof of Proposition~\ref{prop.TwoScaleElliptic}, 
		\begin{align}\label{eq.GKL2}
			\nonumber \lambda \vert \Er[(W-U_\lambda)(W-\Ub)]\vert  &\leq \frac{\lambda}{2} \norm{W-U_\lambda}^2_{\cL^2} + \frac{\lambda}{2} \norm{W-\Ub}^2_{\cL^2} \\
			&\leq \frac{\lambda}{2} \norm{W-U_\lambda}^2_{\cL^2} + \frac{C \lambda 3^{2n}}{2} \norm{\ell_{p,m}}^2_{\acH^1(\cu_m)}. 
		\end{align}
		For the second term on the {\rhs} of \eqref{eq.GKEllipticDecom}, we follow Step~2 in the proof of Proposition~\ref{prop.TwoScaleElliptic}, especially the estimate \eqref{eq.TwoScaleMixed}, to get
		\begin{equation}\label{eq.GKH1}
            \begin{split}
			&\Ll\vert \Er\Ll[\int_{\cu_m} \nabla (W-U_\lambda) \cdot  (\a \nabla W - \ab \nabla \Ub) \, \d \mu\Rr]\Rr\vert \\
            &\leq C 3^{-\alpha n}\norm{W-U_\lambda}_{\acH^1(\cu_m)} \norm{\ell_{p,m}}_{\acH^1(\cu_m)} \\
			& \leq \frac{C}{2}\norm{W-U_\lambda}^2_{\acH^1(\cu_m)}  + \frac{C 3^{-2\alpha n}}{2} \norm{\ell_{p,m}}^2_{\acH^1(\cu_m)}.
            \end{split}
		\end{equation}
		Notice that we do not have the error corresponding to the term $ 3^{2n} \vert\nabla \nabla \ub\vert^2$ in \eqref{eq.TwoScaleMixed}, because the second derivative of the affine function $\ell_{p,m}$ vanishes.
		
		We put the estimate \eqref{eq.GKH1} and \eqref{eq.GKL2} back to \eqref{eq.GKEllipticDecom}, which results in
		\begin{multline}\label{eq.GKTwoScale}
			\lambda \norm{W-U_\lambda}^2_{\cL^2} + \Er \Ll[ \int_{\cu_m} \nabla (W - U_\lambda) \cdot \a \nabla (W - U_\lambda) \, \d \mu\Rr] \\
			\leq C(\lambda 3^{2n} + 3^{-2\alpha n})\norm{\ell_{p,m}}^2_{\acH^1(\cu_m)}.
	\end{multline}
		
		
		\smallskip
		\textit{Step~2: replacement in the Green--Kubo formula.} Now, we replace $U_\lambda$ in \eqref{eq.GreenKuboReduced} by $W$ defined in \eqref{eq.GKdefW},
		\begin{align*}
			&\Ll\vert p \cdot \ab p - \frac{1}{\rho |\cu_m|} \E_\rho \Ll[\int_{\cu_m} p \cdot \a \nabla U_{\lambda}\, \d \mu \Rr]\Rr\vert \\
			& \leq \Ll\vert p \cdot \ab p - \frac{1}{\rho |\cu_m|} \E_\rho \Ll[\int_{\cu_m} p \cdot \a \nabla W \, \d \mu \Rr]\Rr\vert + \Ll\vert \frac{1}{\rho |\cu_m|} \E_\rho \Ll[\int_{\cu_m} p \cdot \a \nabla (W - U_\lambda) \, \d \mu \Rr]\Rr\vert \\
			& = \Ll\vert p \cdot \ab p - \frac{1}{\rho |\cu_n|} \E_\rho \Ll[\int_{\cu_n} p \cdot \a (p + \nabla \phi_{p, \cu_n})  \, \d \mu \Rr]\Rr\vert + \Ll\vert \frac{1}{\rho |\cu_m|} \E_\rho \Ll[\int_{\cu_m} p \cdot \a \nabla (W - U_\lambda) \, \d \mu \Rr]\Rr\vert. 
		\end{align*}
		Compared to $U_\lambda$, its two-scale expansion has a better structure and the first term above becomes the homogenization of the diffusion matrix. We apply \eqref{eq:QuantConvergence} to the first term, and \eqref{eq.GKTwoScale} to the second term, then obtain
		\begin{align}\label{eq.GreenKuboBoundMixed}
			\Ll\vert p \cdot \ab p - \frac{1}{\rho |\cu_m|} \E_\rho \Ll[\int_{\cu_m} p \cdot \a \nabla U_{\lambda}\, \d \mu \Rr]\Rr\vert \leq C(\lambda^{\frac{1}{2}} 3^{n} + 3^{-\alpha n}).
		\end{align}
		
		\smallskip
		\textit{Step~3: choice of the mesoscopic scale.} 
		Finally, we need a reasonable choice of the mesoscopic scale $3^n$ to quantify the error term in \eqref{eq.GreenKuboBoundMixed}.
		\begin{itemize}
			\item Regime $\lambda \in [1, \infty)$. This case is trivial, as the regularization from the resolvent is too strong to see any homogenization, thus our bound is of constant order. To analyze it, instead of applying \eqref{eq.GreenKuboBoundMixed}, we turn back directly to the last line \eqref{eq.DynamicLimit}, and use the $\cH^1$ estimate of $V_\lambda$ from \eqref{eq.defV0}.
			\item Regime $\lambda \in (3^{-2(1+\alpha)m}, 1)$. In the regime $(3^{-2(1+\alpha)m}, 1)$, the effect of regularization is still strong, but it allows us to see the homogenization in a mesoscopic scale roughly around $\lambda^{-\frac{1}{2}}$. In \eqref{eq.GreenKuboBoundMixed}, we choose $3^{n} \simeq \lambda^{- \frac{1}{2(1+\alpha)}}$, and obtain the convergence rate $\lambda^{\frac{\alpha}{2(1+\alpha)}}$. 
			\item Regime $\lambda \in [0, 3^{-2(1+\alpha)m}]$. In this regime, as the regularization is much smaller than the spectral gap of the domain, the error term from the regularization $\lambda^{\frac{1}{2}} 3^{n}$ in \eqref{eq.GreenKuboBoundMixed} is negligible compared to that from the homogenization $3^{-\alpha m}$, and the latter dominates the convergence rate of the Green--Kubo formula.
		\end{itemize}
	\end{proof}
	
	
	%
	%
	%
	%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	%
	%
	%
	
	%	\section{Variations on the Green--Kubo formula}
	%	\label{s.variations}
	%	
	%The goal of this section is to discuss the relationship between the formula for $\ab$ in Theorem~\ref{t.quanti.gk} and possibly more classical forms of the Green--Kubo formula referring to current-current correlations. We will seek neither full precision nor generality, as our aim is only to give some background and motivation for the formula in Theorem~\ref{t.quanti.gk}. For convenience, we assume throughout that for every $\mu\in \mmd(\Rd)$, the mapping 
	%\begin{equation*}  %\label{e.}
	%x \mapsto \a(\mu + \de_x, x)
	%\end{equation*}
	%is differentiable, and that the derivative is bounded uniformly over $\mu \in \mmd(\Rd)$. For every $p \in \Rd$, $\mu \in \mmd(\Rd)$ and $x \in \supp \mu$, we write
	%\begin{equation*}  %\label{e.}
	%\nabla_x \cdot(\a(\mu,x)p),
	%\end{equation*}
	%to denote the divergence in the $x$ variable of the mapping which, using more precise notation, one might alternatively want to denote as
	%\begin{equation*}  %\label{e.}
	%\nabla_x \cdot (\a(\mu' + \de_x, x) p), \qquad \text{ with } \mu' := \mu - \de_x.
	%\end{equation*}
	%The quantity $\nabla_x \cdot (\a(\mu,x) p)$ is meant to capture the $\EE_\mu$-averaged instantaneous current at a point $x \in \supp \mu$ at time zero in the direction of $p \in \Rd$. It is the continuous analogue of the quantity 
	%$j_l(\eta)$ in \cite[(II.2.23)]{spohn2012large}, with the correspondence of $(l,\eta,0)$ in~\cite{spohn2012large} to $(p,\mu,x)$ in our notation. 		For each $p \in \Rd$ and $m \in \N$, we let
	%		\begin{equation*}  %\label{e.}
		%			j_{p,m}(\mu) := \int_{\cu_m} \nabla_x \cdot (\a(\mu,x) p) \, \d \mu(x).
		%		\end{equation*}
	%We recall that we denote by $(P_t^{(m)})_{t \ge 0}$ the semigroup constructed \jc{somewhere}. Heuristically, the associated stochastic process is such that any particle that touches the boundary of $\cu_m$ instantaneously disappear, while particles spontaneously ``appear from the outside'' in order to maintain detailed balance. As will be seen, in order to account for the influx of particles ``coming from the outside'', we need to add some boundary term in the current-current correlation formula. To lighten the exposition we drop the quotation marks when speaking of particles ``coming from the outside'', but we stress that the process is such that they simply appear at the boundary of $\cu_m$ and do not literally come from the outside of the box, since the stochastic process associated with $(P_t^{(m)})_{t \ge 0}$ takes values in $\mmd(\cu_m)$. 
	
	%As a first step, let us ignore the contribution of particles coming from the outside. In analogy with \cite[(II.2.27)]{spohn2012large} or \cite[(2.2), (2.4), (8.7), (8.9)]{mourrat2019efficient} we may hope to represent $p \cdot \ab p$ as
	%\begin{equation}  %\label{e.}
	%p \cdot \ab p \stackrel{?}{=} \E_\rho[p \cdot \a p] - \frac 12 \int_0^{+\infty} \lim_{m \to \infty} \frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_t^{(m)}(j_{p,m})  \Rr] \, \d t.
	%\end{equation}
	%The factor of $\frac 1 2$ in front of the integral sign comes from the fact that we defined the semigroup with respect to the operator $\frac 1 2 \nabla \cdot \a \nabla$ and not $\nabla \cdot \a \nabla$. We localize this in time and space by considering finite $m \in \N$ and by introducing an exponential damping with parameter $\lambda > 0$ to the integral, leading to
	%\begin{equation}
	%\label{e.finitary}
	%\E_\rho[p \cdot \a p] - \frac 12 \frac{1}{\rho |\cu_m|}\int_0^{+\infty} e^{-\lambda t}  \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_t^{(m)}(j_{p,m})  \Rr] \, \d t.
	%\end{equation}
	%We have $P_0^{(m)}(j_{p,m}) = j_{p,m}$. One convenient way to think of $P_t^{(m)}(j_{p,m})$ (which would not require  any assumption on the regularity of $\a$) is as the time derivative of $2P_t^{(m)}(\ell_{p,m})$, where we recall that $\ell_{p,m}$ is defined in \eqref{e.def.ellp}. Integrating by parts, we can write
	%\begin{multline*}  %\label{e.}
	%\frac 12\int_0^{+\infty} e^{-\lambda t}  \E_\rho \Ll[ j_{p,m} P_t^{(m)}(j_{p,m})  \Rr] \, \d t 
	%\\
	%= \E_\rho \Ll[ j_{p,m} \Ll(P_0^{(m)}(\ell_{p,m}) - \int_0^{+\infty} \lambda {e^{-\lambda t}} P_t^{(m)}(\ell_{p,m}) \, \d t\Rr) \Rr] .
	%\end{multline*}
	
	
	
	%\bigskip\bigskip\bigskip
	%		
	%		From the definition, we may rewrite the Green--Kubo formula in the form of 
	%		\begin{equation}  %\label{e.}
		%			p \cdot \ab p = \E_\rho[p \cdot \a p] -  \lim_{t \to \infty} \lim_{m \to \infty}\frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) \big(P_t^{(m)}(\ell_{p,m}) - P_0^{(m)}(\ell_{p,m}) \big) \Rr] ,
		%		\end{equation}
	%		and maybe this can be rewritten as
	%		\begin{equation*}  %\label{e.}
		%			p\cdot \ab p = -  \lim_{t \to \infty}  \lim_{m \to \infty}\frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_t^{(m)}(\ell_{p,m})  \Rr].
		%		\end{equation*}
	%		This looks very much like \cite[(II.2.29)]{spohn2012large}. And now a formal integration by parts gives
	%		\begin{equation*}  %\label{e.}
		%			p \cdot \ab p = \lim_{t \to \infty} \lim_{m \to \infty}  \frac{1}{\rho |\cu_m|} \E_\rho \int_{\cu_m} p \cdot \a \nabla \big(P_t^{(m)}(\ell_{p,m})\big) \, \d \mu,
		%		\end{equation*}
	%		and this is finally some term that we can analyze without sweating too much about singular things like $\nabla_x \cdot (\a p)$ stuff. So maybe the true quantitative Green--Kubo formula we show is that we can control the difference
	%		\begin{equation}
		%			%\label{e.}
		%			\Ll|p \cdot \ab p -  \frac{1}{|\cu_m|} \E_\rho \int_{\cu_m} p \cdot \a \nabla \big(P_t^{(m)}(\ell_{p,m})\big) \, \d \mu\Rr|
		%		\end{equation}
	%		quantitatively in terms of $m$ and $t$. 
	%		
	%		
	%		
	%		
	%		For time-regularity purposes it would make more sense to try to estimate something like
	%		\begin{equation}\label{e.current}
		%			\Ll|p \cdot \ab p - \E_\rho[p \cdot \a p] + \frac 12 \frac{1}{\rho |\cu_m|}\int_0^{+\infty} e^{-\lambda t} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_t^{(m)}(j_{p,m})  \Rr] \, \d t\Rr|,
		%		\end{equation}
	%		for small $\lambda > 0$ and large $m \in \N$. Integrating this gives
	%		
	%		\begin{equation}  
		%			\label{e.integrated.gk}
		%			p \cdot \ab p = \E_\rho[p \cdot \a p] +  \lim_{\lambda \to 0} \lim_{m \to \infty}\frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) \Ll(P_0^{(m)}(\ell_{p,m}) - \int_0^{+\infty} \lambda {e^{-\lambda t}} P_t^{(m)}(\ell_{p,m}) \, \d t\Rr) \Rr] ,
		%		\end{equation}
	%		
	%		
	%		Letting $U := \int_0^{+\infty} \lambda e^{-\lambda t} P_t^{(m)}(\ell_{p,m}) \, \d t$, we have
	%		\begin{equation*}
		%			(\lambda - \nabla \cdot \ab \nabla )U =  \lambda \ell_{p,m}.
		%		\end{equation*}
	%	
	%	
	

\bigskip	
	\noindent \textbf{Acknowledgements}
	The research of CG is supported in part by the National Key R\&D Program of China (No. 2023YFA1010400) and NSFC (No. 12301166), and part of this project was developed during his visit at TSIMF and RIMS. We would like to thank Zhen-Qing Chen, Xiaodan Li and Xiangchan Zhu for helpful discussions.
	
	
	
	\iffalse{
	\section{Letovers}\label{s.informal}
	
	\jc{This is a temporary section that collects the leftovers we have left on our way. }
	
	
	
	
	%	\begin{equation*}
		%		\E_\rho \EEbmu \Ll[ Y_t(f) Y_s(g) \Rr] = \int f(x) S(t-s,y-x) g(y) \, \d x \, \d y.
		%	\end{equation*}
	%	
	\gu{I guess that if we treat the $L^1$ norm on the right-hand side in \eqref{e.two.point}, we could get the estimate like heat kernel of Nash
		\begin{multline}
			\label{e.two.pointL1}
			\Ll|\E_\rho \EEbmu \Ll[ Y_t(f) Y_s(g) \Rr] -\rho \int_{\Rd \times \Rd} f(x)\Psi_{t-s}(y-x)g(y)\, \d x  \d y \Rr|
			\\
			\le C (t-s)^{-\frac{d}{2}-\beta} \|f\|_{L^1} \, \|g\|_{L^1}. 
		\end{multline}	
		The proof follows Zhikov is correct, and the $L^2$ estimate cannot be much better as one can choose the function with arbitrary large support. However, when we treat the fluctuation, usually we are in the case of a \textbf{fixed} test function, thus we hope to gain more factor. After scaling, this extra time factor compensates the spatial normalization, and we still have 
		\begin{equation}  %\label{e.}
			\Ll|  \E_\rho\EE_{\bmu_0} [Y^N_t(f) Y^N_s(g)]- \E[Y^\infty_t(f) Y^\infty_s(g)] \Rr| \le C (N^2(t-s))^{-\beta} \|f\|_{L^1} \, \|g\|_{L^1}. 
		\end{equation}
		If all these estimates are true, we can even recover the general $L^p$ decay.
	}
	
	
	%\end{equation}
	%Formally, the two-point function $(S_t(x))_{t \ge 0, x \in \Rd}$ is such that, for every $f, g \in C^\infty_c(\Rd;\R)$, we have
	%\begin{equation}
	%\label{e.intro.def.two.point}
	%\E_\rho \EEbmu \Ll[ X_t(f) X_s(g) \Rr] = \int_{\Rd \times \Rd} f(x) S_{t-s}(y-x) g(y) \, \d x \, \d y.
	%\end{equation}
	% The \emph{Green--Kubo formula} states that, for every $p,q \in \Rd$, we have
	%\begin{equation}  
	%\label{e.intro.green.kubo}
	%p \cdot \ab q = \lim_{t \to +\infty} \frac{1}{2\rho t} \int_\Rd (p\cdot x) (q \cdot x) S_t(x) \, \d x.
	%\end{equation}
	
	%Here we show the following quantitative version of the Green--Kubo formula. 
	%\begin{theorem}[Quantitative Green--Kubo]
	%%\label{t.green.kubo}
	%There exist an exponent $\beta(d,\Lambda,\rho) > 0$ and a constant $C(d,\Lambda,\rho) < + \infty$ such that for every $p, q \in \Rd$,
	%\begin{equation}
	%\label{e.green.kubo}
	%\Ll|p \cdot \ab q - \frac{1}{2\rho t} \int_\Rd (p\cdot x) (q \cdot x) S_t(x) \, \d x\Rr| \le C t^{-\beta} \, |p| \, |q|. 
	%\end{equation}
	
	\jc{
		
		The following on Green--Kubo formulas is a description of my current best proposal, but can probably also be phrased in other ways. 
		
		While it's not immediate to make sense of the current in continuous space, the $\EEbmu$-averaged current at a point $x \in \supp \mu$, at time zero, and projected onto the direction $p \in \Rd$, is, I think, simply
		\begin{equation*}  %\label{e.}
			\nabla_x \cdot (\a(\mu,x) p) := \sum_{i, j = 1}^d \partial_{x_i} \big(\a_{i,j}(\mu,x) p_j\big).
		\end{equation*}
		I think this is the analogue of $j_l(\eta)$ in \cite[(II.2.12)]{spohn2012large}, with the correspondence $(l,\eta,0)$ in~\cite{spohn2012large} to $(p,\mu,x)$ in our notation. 
		
		For each $p \in \Rd$ and $m \in \N$, we let
		\begin{equation*}  %\label{e.}
			j_{p,m}(\mu) := \int_{\cu_m} \nabla_x \cdot (\a(\mu,x) p) \, \d \mu(x).
		\end{equation*}
		This does not look very well-defined, but I think we can make sense of $(P_t(j_{p,m}))_{t\ge 0}$ or of something closely related. One way to do so is to say that, with the notation
		\begin{equation*}  %\label{e.}
			\ell_{p,m}(\mu) := \int_{\cu_m} p\cdot x \, \d \mu(x),
		\end{equation*}
		we first make sense of $(P_{t}^{(m)}(\ell_{p,m}))_{t \ge 0}$, where $(P_t^{(m)})_{t \ge 0}$ is the semigroup for the evolution restricted to the box $\cu_m$ (say with Neumann boundary condition). And then the time derivative of this is $\frac 1 2 P_t^{(m)}(j_{p,m})$, I think (we may just take this as a definition) (the $\frac 1 2$ is again because of the $\frac 1 2$ in front of $\nabla \cdot \a \nabla$). For me the Green--Kubo formula is the claim that 
		\begin{equation}  %\label{e.}
			p \cdot \ab p = \E_\rho[p \cdot \a p] - \frac 12 \int_0^{+\infty} \lim_{m \to \infty} \frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_t^{(m)}(j_{p,m})  \Rr] \, \d t.
		\end{equation}
		The factor $\frac 1 2$ just before the limit in $m$ is because we use the semigroup with a factor of $\frac 1 2$ in front of $\nabla \cdot \a \nabla$. 
		The formula above looks similar to \cite[(II.2.27)]{spohn2012large}, or \cite[(2.2), (2.4), (8.7), (8.9)]{mourrat2019efficient} for instance. So one way to quantify the Green--Kubo formula would be to estimate the smallness of
		\begin{equation}
			%\label{e.}
			\Ll|p \cdot \ab p - \E_\rho[p \cdot \a p] + \frac 12 \frac{1}{\rho |\cu_m|}\int_0^{t} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_s^{(m)}(j_{p,m})  \Rr] \, \d s\Rr|
		\end{equation}
		in terms of $t$ and $m$. (I would guess that what is inside the absolute values above is non-positive in the limit of large $m$.)
		
		
		From the definition, we may rewrite the Green--Kubo formula in the form of 
		\begin{equation}  %\label{e.}
			p \cdot \ab p = \E_\rho[p \cdot \a p] -  \lim_{t \to \infty} \lim_{m \to \infty}\frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) \big(P_t^{(m)}(\ell_{p,m}) - P_0^{(m)}(\ell_{p,m}) \big) \Rr] ,
		\end{equation}
		and maybe this can be rewritten as
		\begin{equation*}  %\label{e.}
			p\cdot \ab p = -  \lim_{t \to \infty}  \lim_{m \to \infty}\frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_t^{(m)}(\ell_{p,m})  \Rr].
		\end{equation*}
		This looks very much like \cite[(II.2.29)]{spohn2012large}. And now a formal integration by parts gives
		\begin{equation*}  %\label{e.}
			p \cdot \ab p = \lim_{t \to \infty} \lim_{m \to \infty}  \frac{1}{\rho |\cu_m|} \E_\rho \int_{\cu_m} p \cdot \a \nabla \big(P_t^{(m)}(\ell_{p,m})\big) \, \d \mu,
		\end{equation*}
		and this is finally some term that we can analyze without sweating too much about singular things like $\nabla_x \cdot (\a p)$ stuff. So maybe the true quantitative Green--Kubo formula we show is that we can control the difference
		\begin{equation}
			%\label{e.}
			\Ll|p \cdot \ab p -  \frac{1}{|\cu_m|} \E_\rho \int_{\cu_m} p \cdot \a \nabla \big(P_t^{(m)}(\ell_{p,m})\big) \, \d \mu\Rr|
		\end{equation}
		quantitatively in terms of $m$ and $t$. 
		
		
		
		
		For time-regularity purposes it would make more sense to try to estimate something like
		\begin{equation}\label{e.current}
			\Ll|p \cdot \ab p - \E_\rho[p \cdot \a p] + \frac 12 \frac{1}{\rho |\cu_m|}\int_0^{+\infty} e^{-\lambda t} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) P_t^{(m)}(j_{p,m})  \Rr] \, \d t\Rr|,
		\end{equation}
		for small $\lambda > 0$ and large $m \in \N$. Integrating this gives
		
		\begin{equation}  
			\label{e.integrated.gk}
			p \cdot \ab p = \E_\rho[p \cdot \a p] +  \lim_{\lambda \to 0} \lim_{m \to \infty}\frac{1}{\rho |\cu_m|} \E_\rho \Ll[ P_0^{(m)}(j_{p,m}) \Ll(P_0^{(m)}(\ell_{p,m}) - \int_0^{+\infty} \lambda {e^{-\lambda t}} P_t^{(m)}(\ell_{p,m}) \, \d t\Rr) \Rr] ,
		\end{equation}
		
		
		Letting $U := \int_0^{+\infty} \lambda e^{-\lambda t} P_t^{(m)}(\ell_{p,m}) \, \d t$, we have
		\begin{equation*}
			(\lambda - \nabla \cdot \ab \nabla )U =  \lambda \ell_{p,m}.
		\end{equation*}
	}
	
	
	\gu{An alternative definition of Green--Kubo formula. We can try the definition in \cite{spohn2012large} (II.2.16) that 
		\begin{align*}
			(\ab)_{i,j} = \lim_{t \to \infty}\frac{1}{\rho t} \int_{\Rd} S(t,x) x_i x_j \d x.
		\end{align*}
		One can treat it as the limit of some $\frac{1}{\rho t}\Er[Y_0(f) Y_t(g)]$. However, in the original expression, the first function is too singular (Dirac mass), while the second function is not integrable. Thus, we propose other candidate.
		
		A possible idea is to use the equilibrium fluctuation with scaling 
		\begin{align*}
			\Ll|  \frac{1}{\rho t}\E_\rho\EE_{\bmu_0} [Y^N_0(f) Y^N_t(g)]- \frac{1}{t}\int_{\Rd \times \Rd} f(x)\Psi_{t}(y-x)g(y)\, \d x  \d y  \Rr| \\
			\le C t^{-1}(N^2 t)^{-\beta} \|f\|_{L^2} \, \|g\|_{L^2}.
		\end{align*}
		We can still choose different $g^{(\ell)} = \Ind{\ell \cu} x_i x_j$	
		\begin{align*}
			\Ll|  \frac{1}{\rho t_{\ell}}\E_\rho\EE_{\bmu_0} [Y^N_0(f) Y^N_t(g^{(\ell)})]- \frac{1}{t_{\ell}}\int_{\cu}\int_{\ell \cu} \Psi_{t_{\ell}} (y-x)y_i y_j\, \d x  \d y  \Rr| \\
			\le C t_{\ell}^{-1}(N^2 t_{\ell})^{-\beta} \|f\|_{L^2} \, \|g^{(\ell)}\|_{L^2}.
		\end{align*}
		Notice that $\|g^{(\ell)}\|_{L^2} \simeq \ell^{\frac{d+4}{2}}$, we need this quantity be overcome by $N^{-2\beta}$, thus we take $\ell = N^{\frac{4\beta}{d+4}}$. Moreover, we also need $\sqrt{t_\ell} \ll \ell$, for example we can choose $t_\ell = \ell$.
		
		
	}
	
	\jc{
		Using a previous result from Chenlin, it should be possible to argue that we can take $\ell = t^{\frac 1 2 + \ep}$ (in the sense that we can neglect what falls outside of the cube). We can control the error term provided that $t$ is a small power of $N$. I'm not sure how to interpret this. My guess is that we are seeing a limitation of our $L^2$ estimates. Personally I have a preference for the other Green--Kubo thing like in \eqref{e.integrated.gk}.
		
	}
	
	
	
	
	
	
	\subsection{Parabolic homogenization implies central fluctuation}
	Let $\a$ be the random conductance as defined in the previous paper. We denote by $(\mmd(\Rd), (\fil_t)_{t \geq 0}, \P_{\mu})$ the probability space for the diffusion process $(\mu_t)_{t \geq 0}$ of particles $\mu_t = \sum_{i=1}^\infty \delta_{X_i(t)}$ defined by the operator $- \frac{1}{2}\nabla \cdot \a \nabla$, starting from $\mu$, and $(\mmd(\Rd), \mcl{F}_{\Rd}, \Pr)$ the Poisson point process. We then use $(\mmd, \mcl{F}, \P)$ for the probability space by sampling at first $\mu$ under $\Pr$, then sampling the diffusion $(\mu_t)_{t \geq 0}$ under $(\mmd(\Rd), (\fil_t)_{t \geq 0}, \P_{\mu})$.
	
	We are interested in the \emph{fluctuation} theorem in interacting particles. Let $m$ be the standard Lebesgue measure, and we define a random field ${Y^N_t := N^{-d/2}(\mu_{N^2 t} - \rho m)(\cdot/N)}$, i.e.
	\begin{align}\label{eq.defFluctuation}
		Y^N_t(f) = N^{-\frac{d}{2}}\Ll(\int_\Rd f\Ll(\cdot/N\Rr) \, \d \mu_{N^2 t} - \rho \int_\Rd f\Ll(\cdot/N\Rr) \, \d m\Rr).
	\end{align}
	We let $(\xi_t(x))_{t \ge 0, x \in \Rd}$ be a $d$-dimensional space-time white noise over $\R_+ \times \R^d$. The following has been shown for discrete systems and should plausibly be standard to extend to our setting.
	
	\begin{theorem}
		When $N \to \infty$, $Y^N_t$ converges weakly to the limit $Y^\infty_t$ which is a generalized Ornstein-Uhlenbeck process
		\begin{align}
			\partial_t Y^\infty_t = \frac{1}{2} \nabla \cdot \ab \nabla Y^\infty_t + \nabla \cdot\big(\sqrt{\rho \ab} \xi_t\big).
		\end{align}
	\end{theorem}
	
	
	
	%	\begin{remark}
		%		\textcolor{red}{To be confirmed. I believe that here we do not have to define the Laplace term as $\frac{\ab}{\rho}$, because in our definition of $\ab$ we have already divided by $\rho$. It is the case in the finite volume approximation, and in the stationary corrector it is hidden in the Mecke's identity. Thus our definition by stationary corrector has a minor difference up to $\rho$ from the original one.} 
		%	\end{remark}
	
	Let use see how the homogenization $L^2$ estimate can prove this convergence of covariance. By polarization, without loss of generality we can suppose $g=f$ and we denote by $F_N(\mu_{N^2t}) := Y^N_t(f)$. Then we have 
	\begin{align*}
		\E[Y^N_t(f) Y^N_s(f)] &= \E[F_N(\mu_{N^2t}) F_N(\mu_{N^2s})] \\
		&= \E\Ll[\E_{\mu}[F_N(\mu_{N^2t}) F_N(\mu_{N^2s}) \vert \fil_{N^2 s}]\Rr] \\
		&= \E\Ll[F_N(\mu_{N^2s}) \E_{\mu}[F_N(\mu_{N^2t})  \vert \fil_{N^2 s}]\Rr] \\
		&= \E\Ll[F_N(\mu_{N^2s}) (P_{N^2(t-s)}F_N)(\mu_{N^2s}) \Rr].
	\end{align*}
	In the last line, we use the semigroup. Then observing the fact that Poisson measure is invariant measure for the dynamic, so we have 
	\begin{align*}
		\E[Y^N_t(f) Y^N_s(f)] &= \E\Ll[F_N(\mu) (P_{N^2(t-s)}F_N)(\mu) \Rr].
	\end{align*}
	We now replace the semigroup $P_t$ by $\bar{P}_t$ that of the diffusion following $-\frac{1}{2}\nabla \cdot \ab \nabla$, we have 
	\begin{multline}\label{eq.homo}
		\E[Y^N_t(f) Y^N_s(f)] \\
		=  \E\Ll[F_N(\mu) (\bar{P}_{N^2(t-s)}F_N)(\mu) \Rr] + \E\Ll[F_N(\mu) (P_{N^2(t-s)} - \bar{P}_{N^2(t-s)}F_N)(\mu) \Rr]. 
	\end{multline}
	In the first term on the {\rhs} of \eqref{eq.homo}, we can calculate $\bar{P}_{N^2(t-s)}F_N$ explicitly. In fact, the diffusion defined by  $-\frac{1}{2}\nabla \cdot \ab \nabla$ are independent Brownian motions of diffusion coefficient $\sigma^2$. Thus 
	\begin{align*}
		(\bar{P}_{N^2(t-s)}F_N)(\mu) = N^{-\frac{d}{2}}\Ll(\int_{\Rd} \Ll(\Psi_{N^2(t-s)} \ast f\Ll(\cdot/N\Rr)\Rr) \, \d \mu - \rho \int_\Rd f\Ll(\cdot/N\Rr) \, \d m\Rr).
	\end{align*}
	It is easy to verify that this quantity is still centered, so the first term on the {\rhs} of \eqref{eq.homo} becomes 
	\begin{align*}
		\E\Ll[F_N(\mu) (\bar{P}_{N^2(t-s)}F_N)(\mu) \Rr] &= N^{-d}\cov\left[\int_\Rd f\Ll(\cdot/N\Rr) \, \d \mu,\int_{\Rd} \Ll(\Psi_{N^2(t-s)} \ast f\Ll(\cdot/N\Rr)\Rr) \, \d \mu \right] \\
		&= \rho N^{-d} \int_{\Rd} f\Ll(\cdot/N\Rr) \Ll(\Psi_{N^2(t-s)} \ast f\Ll(\cdot/N\Rr)\Rr) \, \d m\\
		&= \rho N^{-d} \int_{\Rd \times \Rd} f\Ll(x/N\Rr) \Psi_{N^2(t-s)}(x-y)  f\Ll(y/N\Rr) \, \d x \d y.
	\end{align*}
	This is the desired result.
	For the second term on the {\rhs} of \eqref{eq.homo}, we can use the $L^2$ estimate from homogenization that 
	\begin{align*}
		&\Ll\vert \E\Ll[F_N(\mu) (P_{N^2(t-s)} - \bar{P}_{N^2(t-s)}F_N)(\mu) \Rr]\Rr\vert   \\
		& \leq \Ll(\E[(F_N)^2] \E[(P_{N^2(t-s)} - \bar{P}_{N^2(t-s)}F_N)^2]\Rr)^{\frac{1}{2}} \\
		& \leq C (N^2(t-s))^{-\alpha} \xrightarrow{N \to \infty} 0.
	\end{align*}
	
	\jc{
		
		We cannot directly apply the results of \cite{fukushima2010dirichlet} because the space $\mmd(\Rd)$ is not locally compact. We circumvent this by compactifying the space in the following way. We denote by $\bmmd(\Rd)$ the space of measures with countable support and taking values in $\N \cup \{+\infty\}$. We say that a sequence $(\mu_n)_{n \in \N}$ of elements of $\bmmd(\Rd)$ converges vaguely to $\mu \in \bmmd(\Rd)$ if for every $f \in C_c(\Rd;\R_+)$, we have
		\begin{equation*}  %\label{e.}
			\lim_{n \to +\infty} \Ll(1+\int f \, \d \mu_n\Rr)^{-1} = \Ll(1+\int f \, \d \mu\Rr)^{-1},
		\end{equation*}
		where we understand that $(+\infty)^{-1} = 0$.
		We let $\{f_n \ : \ n \in \N\}$ be the countable set of positive and negative parts of the functions in the set $\mcl F$ from \cite[Lemma~5.1]{HJbook}. For every $\mu, \nu \in \bmmd(\Rd)$, we define the distance between $\mu$ and $\nu$ as
		\begin{equation}
			\label{e.def.dist.bmmd}
			\msf d(\mu,\nu) := \sum_{n  = 0}^{+\infty} 2^{-n} \Ll| \Ll(1+\int f_n \, \d \mu\Rr)^{-1} - \Ll(1+\int f_n \, \d \nu\Rr)^{-1}  \Rr| .
		\end{equation}
		Essentially the topology for $\bmmd(\Rd)$ is the same as that of $[0,1]^\N$ for the product topology now. I believe that the space $\bmmd(\Rd)$ is compact. It is clear that if we give ourselves a sequence of measures $(\mu_n)_{n \in \N}$ in $\bmmd(\Rd)$, then we can extract a subsequence such that for each $k \in \N$, we have that $\int f_k \, \d \mu_n$ converges as $n$ tends to infinity along the subsequence. That the exists a measure $\mu \in \bmmd(\Rd)$ such that the limit of the thing above is $\int f_k \, \d \mu$ needs to be explained, but I think it is doable and resembles the argument for \cite[Lemma~5.3]{HJbook}, with the additional thing that we need to explain how to deal with points in the measure $\mu$ that have infinite mass (maybe we construct those in a first step). We still have our spaces $L^2(\P_\rho)$. Perhaps most convenient next is to say that our ``basic functions'' are of the form
		\begin{equation*}  %\label{e.}
			G \Ll( \Ll(1+ \int \phi_1 \, \d \mu\Rr)^{-1}, \ldots, \Ll(1+\int \phi_k \, \d \mu\Rr)^{-1} \Rr) ,
		\end{equation*}
		with $\phi_1, \ldots, \phi_k \in C^{\infty}(\Rd; \R_+)$ and $G \in C^{\infty}(\R^k;\R)$. Then we can define $\cH^1(\Rd)$ by taking the closure of the space of basic functions with respect to the usual norm. The space of basic functions will also be our core. I guess things work ok with this now!
		
		And once we have constructed the semigroup and the process we can observe that the measure $\P_\rho$ is invariant so really the stochastic process is actually taking values in $\mmd(\Rd)$. 
		
		
	}
	
	\jc{I have the impression that things work out ok if we use the usual basic functions and multiply them by some stuff to ensure compact support.
		The ``usual basic functions'' are those of the form
		\begin{equation*}  %\label{e.}
			f(\mu) = G \Ll( \int \phi_1 \, \d \mu, \ldots, \int \phi_k \, \d \mu \Rr) ,
		\end{equation*}
		with $G \in C^\infty(\R^k)$ bounded and with bounded derivatives, and $\phi_1, \ldots, \phi_k \in C^\infty_c(\Rd)$. 
		%We can take $k = 1$ in our discussion, and in this case
		%\begin{equation*}  %\label{e.}
		%\nabla f(\mu,x) = \nabla \phi(x) G' \Ll( \int \phi \, \d \mu \Rr) .
		%\end{equation*}
		The question is whether we can approximate $f$ in $\cH^1(\Rd)$ with a function of compact support. A function is of compact support if whenever you know that it is non-zero, you can deduce an upper bound on the number of points in any fixed compact subset of $\Rd$. We take $\phi \in C^\infty_c(\Rd;\R_+)$ with $\phi = 1$ on $\cu_0$, $\chi \in C^\infty_c(\R; \R_+)$ with $\chi = 1$ on $[0,1]$ and $0 \le \chi \le 1$ everywhere, and for every $M > 0$ $z,x \in \Rd$ and $r \in \R$, we define
		\begin{equation*}  %\label{e.}
			\phi_{M,z}(x) := \phi \Ll( \frac{x-z}{M} \Rr)  \quad \text{ and } \quad \chi_M(r) := \chi \Ll( \frac{r}{M} \Rr) .
		\end{equation*}
		How about we try to approximate $f$ with 
		\begin{equation*}  %\label{e.}
			\mu \mapsto f(\mu) \prod_{z \in M \mathbb{Z}^d} \chi_{2^{|z|} M^{2d}} \Ll( \int \phi_{M,z} \, \d \mu \Rr)  \quad ?
		\end{equation*}
		The gradient of this function has two terms, one with the gradient hitting $f$ and one with the gradient hitting the product. When the gradient hits $f$ we want to say that the product is essentially equal to $1$ for large $M$; this should be ok. When the gradient hits the product, we want to use that $f$ is bounded and that this product has a small $\cH^1$ semi-norm. Up to enlarging a bit the cubes, the norm of the gradient of this product is like
		\begin{equation*}  %\label{e.}
			\sum_{z \in M \mathbf{Z}^d} M^{-1-2d} 2^{-|z|} \1_{x \in z + Q_M} \1_{\mu(z + Q_M) \ge 2^{|z|} M^{2d}}.
		\end{equation*}
		$Q_M$ is the cube of side length $M$ (we could also use $\cu_m$). We square this (essentially we can pass the square inside the sum because of the indicator function with $x$), we integrate in $x$, take the expectation, and use some exponential tail estimate for Poisson random variables to bound this by 
		\begin{equation*}  %\label{e.}
			M^{-1-d} \sum_{z \in M \mathbf{Z}^d} 2^{-|z|} \exp \Ll( -c \rho^{-1} 2^{|z|} M^{d} \Rr) .
		\end{equation*}
		This tends to zero as $M$ tends to infinity, so all is well?
		
		
		
	}
	
	This guarantees that $(P_t)_{t \geq 0}$ is a Markovian semigroup; that is, for every $u \in \cL^2$ with $0 \le u \le 1$, we have $0 \le P_t u \le 1$ for every $t \ge 0$. In order to guarantee the existence of an associated stochastic process, by \cite[Theorem~7.2.1]{fukushima2010dirichlet}, it suffices to verify that the Dirichlet form $(\Ee^\a, \D(\Ee^\a))$ is regular. With $\cC_0(\Rd)$ denoting the space of continuous functions from $\mmd(\Rd)$ to $\R$ with compact support, the regularity of the Dirichlet form means that there exists a set $\cC \subset \D(\Ee^\a) \cap \cC_0(\Rd)$ that is dense in $\cC_0(\Rd)$ for the uniform norm, and dense in~$\D(\Ee^\a)$ for the norm induced by $\Ee^{\a}_1$. Denoting by $\msf d$ the distance on $\mmd(\Rd)$ defined in \cite[(5.7)]{HJbook}, we observe that a valid choice for $\cC$ is to take
	\begin{equation*}  %\label{e.}
		\cC := \Ll\{  \Rr\} 
	\end{equation*}
	
	\textcolor{red}{Some explanation for this.} 
	
	Although we will not make use of this, we mention that one can easily check that the Dirichlet form $(\Ee^\a, \D(\Ee^\a))$ is local (see \cite[Section~1.1]{fukushima2010dirichlet} for a definition), and thus \cite[Theorem~7.2.2]{fukushima2010dirichlet} guarantees that we can associate a diffusion process with this Dirichlet form; in other words, we can make sure that the stochastic process has continuous sample paths. 
	
	
	\subsection{Some discussion about the flux}

		\jc{I think we can leave it like that for the introduction. In Section~\ref{sec:5} we can explain something along the following lines.
		First, for every $t > 0$, we have that $P_t^{(m)}(F_{p,m})$ is well-defined and belongs to $\cH^1_0(\cu_m)$. So by definition,
		\begin{equation*}  %\label{e.}
			\la F_{p,m}, P_t^{(m)}(F_{p,m}) \rah = - \Er \int_{\cu_m}\a p \cdot \nabla (P_t^{(m)}(F_{p,m})) \, \d \mu. 
		\end{equation*}
		We set
		\begin{equation*}  %\label{e.}
			U_{\ep,\lambda} := \int_\eps^{+\infty}  e^{-\lambda t} P_t^{(m)}(F_{p,m}) \, \d t.
		\end{equation*}
		With informal notation, we have that
		\begin{equation}  
			\label{e.eq.for.Uepslambda}
			(\lambda - \frac 1 2 \nabla \cdot \a \nabla )U_{\ep,\lambda} =  e^{-\lambda \ep} P_\ep^{(m)}(F_{p,m}).
		\end{equation}
		When $\ep = 0$, this should mean that for every $\phi \in \cH^1_0(\cu_m)$,
		\begin{equation*}  %\label{e.}
			\lambda \Er[V U_{0,\lambda}] + \frac 1 2 \Er \int_{\cu_m} \nabla V \cdot \a (p + \nabla U_{0,\lambda}) \, \d \mu = 0.
		\end{equation*}
		In other words, $U_{0,\lambda}$ is the minimizer in $\cH^1_0(\cu_m)$ of the functional 
		\begin{equation*}  %\label{e.}
			U \mapsto \lambda \E[U^2] + \frac 1 2 \Er \int_{\cu_m} (p + \nabla U)\cdot \ab(p + \nabla U)\, \d \mu. 
		\end{equation*}
		This is essentially a corrector, up to the zero-order term perturbation. So in this case things work out. If we also know that $U_{\ep,\lambda}$ converges to $U_{0,\lambda}$ as $\ep$ tends to zero then we are all good. I think this should be a consequence of the fact that $P_\ep^{(m)}(F_{p,m})$ converges to $F_{p,m}$ in $\cH^{-1}(\cu_m)$ as $\ep$ tends to zero. 
		
		
		I suppose that there will be some functional analysis stuff to be discussed. For starters, let us denote by $-A_m$ the self-adjoint positive semi-definite operator corresponding to the Dirichlet form with domain $\cH^1_0(\cu_m)$ and bilinear form defined on $\cH^1_0(\cu_m)^2$ as
		\begin{equation*}  %\label{e.}
			(u,v) \mapsto \frac 1 2 \Er \int_{\cu_m} \nabla u \cdot \a \nabla v \, \d \mu. 
		\end{equation*}
		We call $-A_m^\circ$ the version with $\a$ replaced by $\id$.
		We have essentially by definition that 
		\begin{equation*}  %\label{e.}
			\cH^1_0(\cu_m) = \Ll\{f \in L^2(\Pr) \mid \sqrt{-A_m^\circ} f \in L^2(\Pr) \Rr\},
		\end{equation*}
		and $\|f\|_{\cH^1(\cu_m)} = \|\sqrt{\id -A_m^\circ} f\|_{\cL^2}$. From there on I forget about writing the $\circ$ exponent. 
		By the Riesz representation theorem, every element of the dual of $\cH^1_0(\cu_m)$ can be represented as the mapping
		\begin{equation*}  %\label{e.}
			f \mapsto \Er \Ll[ f \chi \Rr] +  \la \sqrt{-A_m} f, \sqrt{-A_m} \chi \ra_{\cL^2},
		\end{equation*}
		for some $\chi \in \cH^1_0(\cu_m)$. We consider the space
		\begin{equation*}  %\label{e.}
			D_m := \Ll\{ \phi \in L^2(\Pr) \mid \phi = (\id-A_m) \chi \quad \text{for some } \chi \in \cH^1_0(\cu_m) \Rr\} ,
		\end{equation*}
		which we equip with the norm $\phi \mapsto \|\phi\|_{\cH^{-1}(\cu_m)} := \|(\id-A_m)^{-1/2} \phi\|_{\cL^2}$. 
		One can check that for $\phi = - A_m \chi$, we have $\|\phi\|_{\cH^{-1}(\cu_m)} = \|\sqrt{\id-A_m}\chi\|_{\cL^2} = \|\chi\|_{\cH^1(\cu_m)}$. We can complete the space $D_m$ with respect to the norm $\|\cdot \|_{\cH^{-1}(\cu_m)}$, and call the completion $\cH^{-1}(\cu_m)$. The space $\cH^{-1}(\cu_m)$ really is the dual space to $\cH^1_0(\cu_m)$ (I guess we need to use a bit of spectral theory to assert this, but we know how to represent every element of the dual of $\cH^1_0(\cu_m)$ using the Riesz representation theorem as discussed above, and we can then truncate a bit the high frequencies of $\chi$ to justify the density thing), and when $\phi \in D_m = \cH^{-1}(\cu_m) \cap L^2(\Pr)$, we make sure that we have for every $f \in \cH^1_0(\cu_m)$ that
		\begin{equation*}  %\label{e.}
			\la \phi, f \ra_{\cH^{-1}, \cH^1_0} = \Er[\phi \, f].
		\end{equation*}
		
		
		For $G \in \cH^{-1}(\cu_m)$, when we say that we solve for $U \in \cH^1_0(\cu_m)$ solution to
		\begin{equation}  
			\label{e.ellipt.problem.with.G}
			\Ll(\lambda - \frac 1 2 \nabla \cdot \a \nabla \Rr)U = G,
		\end{equation}
		what we really mean is that we look for the minimizer in $\cH^1_0(\cu_m)$ of the functional
		\begin{equation*}  %\label{e.}
			U \mapsto  \lambda \Er[U^2] + \frac 1 2 \Er \int_{\cu_m} \nabla U \cdot \a \nabla U \, \d \mu - 2\la G, U \rah.
		\end{equation*}
		By writing the first variation, we see that
		\begin{align*}  %\label{e.}
			\lambda \Er[U^2] + \frac 1 2 \Er \int_{\cu_m} \nabla U \cdot \a \nabla U \, \d \mu 
			& = \la G, U \rah \\
			& \le  \|G\|_{\cH^{-1}(\cu_m)} \, \|U\|_{\cH^1(\cu_m)} \\
			& \le \frac{\lambda'}{2} \|G\|_{\cH^{-1}(\cu_m)}^2 + \frac 1 {2\lambda'} \E[U^2] + \frac 1 {2\lambda'} \Er \int_{\cu_m} |\nabla U|^2 \, \d \mu,
		\end{align*}
		where we write $\lambda' := \max(\lambda, \frac 1 2)$. Using also \eqref{e.ass.unif.ell}, we conclude that
		\begin{equation*}  %\label{e.}
			\lambda \Er[U^2] + \frac 1 2 \Er \int_{\cu_m} \nabla U \cdot \a \nabla U \, \d \mu  \le \lambda' \|G\|_{\cH^{-1}(\cu_m)}.
		\end{equation*}
		This shows that the mapping sending $G \in \cH^{-1}(\cu_m)$ to the solution $U \in \cH^1_0(\cu_m)$ of~\eqref{e.ellipt.problem.with.G} is continuous. Recalling the definition of $U_{\ep,\lambda}$ from \eqref{e.eq.for.Uepslambda}, if we show something for $U_{0,\lambda}$ that depends continuously on $U_{0,\lambda}$ in the sense of the $\cH^1(\cu_m)$ norm, then this something automatically transfers to $U_{\ep,\lambda}$ for small $\ep$. So I think we are good.
		
		
	}


	\gucomment{More discussions.}
	
			\jc{
		
		To be precise, the function $V := U - \ell_{p,m}$ is the unique minimizer in $\cH^1_0(\cu_m)$ of the functional
		\begin{equation*}  %\label{e.}
			V \mapsto  \frac \lambda 2 \Er [(\ell_{p,m} + V)^2] + \frac 1 2 \Er \int_{\cu_m} (p + \nabla V) \cdot \a (p + \nabla V) \, \d \mu - \lambda \E[\ell_{p,m} V],
		\end{equation*}
		or equivalently, of the functional
		\begin{equation*}  %\label{e.}
			V \mapsto \frac \lambda 2 \Er [V^2] + \frac 1 2 \Er \int_{\cu_m} (p + \nabla V) \cdot \a (p + \nabla V) \, \d \mu .
		\end{equation*}
		Equivalently, it is the unique $V \in \cH^1_0(\cu_m)$ such that for every $\phi \in \cH^1_0(\cu_m)$,
		\begin{equation*}  %\label{e.}
			\lambda \Er[\phi V] + \Er \int_{\cu_m} \nabla \phi \cdot \a(p+ \nabla V) \, \d \mu = 0.
		\end{equation*}
		Notice that the mapping
		\begin{equation*}  %\label{e.}
			F : \phi \mapsto -\Er \int_{\cu_m} \nabla \phi \cdot \a p  \, \d \mu
		\end{equation*}
		is an element of $\cH^{-1}(\cu_m)$, the dual space to $\cH^1_0(\cu_m)$. It is tempting to denote this guy by $\nabla \cdot (\a p)$, but I think this is missing some boundary term (which is relevant when discussing the Green--Kubo formula). Informally, we are solving for $V \in \cH^1_0(\cu_m)$ solution to
		\begin{equation*}  %\label{e.}
			(\lambda  - \nabla \cdot \a \nabla) V = F.
		\end{equation*}
		We said that we denote by $(P_t^{(m)})_{t \ge 0}$ the semigroup associated with the operator $-\nabla \cdot \a \nabla$ with domain $\cH^1_0(\cu_m)$; we can write $(R_\lambda^{(m)})_{\lambda > 0}$ for the corresponding resolvents. A priori these objects takes elements of $L^2(\Pr)$ as arguments. The right-hand side $F$ is not a priori an element of $L^2(\Pr)$, but if it was we would have liked to say that $V = R_\lambda^{(m)} F$. Maybe we can still make sense of this using some spectral decomposition stuff, I haven't thought about it. One point I want to stress though is that on the other hand, it is really not the case that $\ell_{p,m} + V = R_{\lambda}^{(m)}(\lambda \ell_{p,m})$. One way to see this clearly is that the image of $R_{\lambda}^{(m)}$ should be a subset of $\cH^1_0(\cu_m)$, but $\ell_{p,m} + V$ is not in there, the boundary conditions do not match. 
		
		Continuing with the informal discussion, we would like to represent $V$ as some exponentially weighted time integral of what we can informally denote as $t \mapsto P_t^{(m)}(F)$. I say informally because again a priori we are supposed to feed $L^2(\Pr)$ functions into $P_t^{(m)}$ but $F$ is not such a thing (even under the assumption that $\a$ is smooth I think it contains some boundary term). 
		This $u(t,\cdot) := P_t^{(m)}(F)$ might be the true guy that is supposed to be appearing in the Green--Kubo formula. I think it is a minimizer in the space
		\begin{equation*}  %\label{e.}
			\{u \in L^2([0,T]; \cH^1_0(\cu_m)) \ : \  \dr_t u \in L^2([0,T]; \cH^{-1}(\cu_m))\}
		\end{equation*}
		of the functional
		\begin{equation*}  %\label{e.}
			u \mapsto \frac 1 2 \E[u(T,\cdot)^2] + \frac 1 2  \E\int_0^T \int_{\cu_m} (p + \nabla u) \cdot \a (p + \nabla u) \, \d \mu \, \d t.
		\end{equation*}
		We may take some inspiration from \cite[Appendix~A]{armstrong2018parabolic} to show that the minimizer of this functional is indeed the object we want. Informally, we can think of
		\begin{equation*}  %\label{e.}
			\Ll\{
			\begin{aligned}
				& (\partial_t - \nabla \cdot \a \nabla )u = 0, \\
				& u(0,\cdot) = F.
			\end{aligned}
			\Rr.
		\end{equation*}
		I think the Green--Kubo formula becomes valid if we write it in the form of
		\begin{equation}  %\label{e.}
			p \cdot \ab p = \E_\rho[p \cdot \a p] - \frac 12 \int_0^{+\infty} \lim_{m \to \infty} \frac{1}{\rho |\cu_m|} \E_\rho \Ll[ F P_t^{(m)}(F)  \Rr] \, \d t.
		\end{equation}
		The notation is informal but, first of all $F$ depends on $m$ and $p$, and more importantly it's not really that we are computing $\E_\rho \Ll[ F P_t^{(m)}(F)  \Rr] $ since $F$ is not defined as a function but as an element of $\cH^{-1}(\cu_m)$, so really we mean something like 
		\begin{equation*}  %\label{e.}
			\E_\rho \Ll[ F P_t^{(m)}(F)  \Rr] \leadsto \langle F, P_t^{(m)}(F) \rangle_{\cH^{-1}(\cu_m), \cH^1_0(\cu_m)}.
		\end{equation*}
		This is not great notation but I mean the duality pairing between $\cH^{-1}(\cu_m)$ and $\cH^1_0(\cu_m)$. In short, by definition this is 
		\begin{equation*}  %\label{e.}
			-\E \int_{\cu_m} \nabla (P_t^{(m)}(F)) \cdot \a p \, \d \mu,
		\end{equation*}
		as we wanted from the beginning. The $F$ guy is taking care of the boundary term I was worried about on its own.
		
		(Incidentally, we can probably also write down a similar story with $\cH^1_0(\cu_m)$ replaced by $\cH^1(\cu_m)$ now.)
		
	}
	%		
	%\jc{I think the point is that we want to find $\bar V$, either in $\cH^1_0(\cu_m)$ or in $\cH^1(\cu_m)$ depending on the boundary condition we want, such that 
		%\begin{equation*}  %\label{e.}
		%(\lambda - \nabla \cdot \a \nabla )(\ell_{p,m} + \bar V) = \lambda \ell_{p,m}.
		%\end{equation*}
		%I'm not saying we should do both, we can stick to the version with $\cH^1_0$ if we want.
		%}		
	%

	}\fi


	\bibliographystyle{plain}
	\bibliography{KawasakiRef}
	
\end{document}

