\section{Introduction}\label{intro}
We consider the kinetic Fokker-Planck operator $\mathscr{L}_{\scriptscriptstyle{\rm FP}}:=\nabla_v\cdot\left(\nabla_{v}+v\right)$ and the following spatially inhomogeneous nonlinear drift-diffusion model, 
\begin{equation}\label{nk}
\left\{ 
\begin{aligned}
\ &\left(\partial_t+v\cdot\nabla_x\right)f(t,x,v)=\rho^{\beta}_f(t,x)\,\mathscr{L}_{\scriptscriptstyle{\rm FP}} f(t,x,v),  \\
\ &\;f(0,x,v)=f_{\rm in} (x,v), \\
\end{aligned}
\right. 
\end{equation} 
for an unknown $f(t,x,v)\ge0$ with $t\in\mathbb{R}_+$, $x\in\mathbb{T}^d$ or $\mathbb{R}^d$, $v\in\mathbb{R}^d$, where $\mathbb{T}^d$ denotes the $d$-dimensional torus with unit volume, the constant $
{\beta}\in[0,1]$ and 
\begin{equation*}
\rho_f(t,x) := \int_{\mathbb{R}^{d}} f(t,x,v) \mathop{}\!\mathrm{d} v. 
\end{equation*}
Given a constant ${\epsilon}\in(0,1)$, the equation under the low field scaling $t\mapsto{\epsilon}^2 t$, $x\mapsto{\epsilon} x$ reads 
\begin{equation}\label{epfnk}
\left\{ 
\begin{aligned}
\ &\left({\epsilon}\partial_t+v\cdot\nabla_x\right)f_{\epsilon}(t,x,v)=\frac{1}{{\epsilon}}\rho^{\beta}_{f_{\epsilon}}(t,x)\,\mathscr{L}_{\scriptscriptstyle{\rm FP}} f_{\epsilon}(t,x,v),  \\
\ &\;f_{\epsilon}(0,x,v) = f_{{\epsilon},\rm in} (x,v). \\
\end{aligned}
\right. 
\end{equation} 
Our aim is to show the global well-posedness and the trend to equilibrium with smoothness a priori estimates for the equation~\eqref{nk}, and the quantitative asymptotic dynamics of the equation~\eqref{epfnk} as ${\epsilon}$ tends to zero. 

\subsection{Main results}
Let us recall that a classical solution of an evolution equation is a nonnegative function verifying the equation pointwise everywhere and 
matching the initial data continuously.
Unless otherwise specified, any solution we consider below is intended in the classical sense. 
For $k\in\mathbb{N}$, $\mathcal{C}^k({\Omega})$ is the set of functions having all derivatives of order less than or equal to $k$ continuous in the domain ${\Omega}$. 
For ${\alpha}\in(0,1)$, $\mathcal{C}^{\alpha}({\Omega})$ is the classical H\"older space on ${\Omega}$ with exponent ${\alpha}$. 
Besides, we write the measure $\mathop{}\!\mathrm{d} m:=\mathop{}\!\mathrm{d} x\mathop{}\!\mathrm{d}{\mu}$, where
\begin{equation*}
 {\mu}(v):=(2\pi)^{-\frac{d}{2}}e^{-\frac{|v|^2}{2}}  {\quad\rm and\quad}   
\mathop{}\!\mathrm{d}{\mu}:={\mu}\mathop{}\!\mathrm{d} v 
\end{equation*}
denote the Gaussian function and the Gaussian measure, respectively. 
A function takes the form of $C{\mu}(v)$ for some constant $C>0$ is called a Maxwellian.


\begin{theorem}\label{wellpose}
Let the space domain ${\Omega}_x=\mathbb{T}^d$ or $\mathbb{R}^d$, and the constants $0<\lambda<\Lambda$ be given. 
\begin{item}
\item [\bf {\ \,(i)}] 
If $f_{\rm in}\in{\mathcal C}^0({\Omega}_x\times\mathbb{R}^d)$ satisfies $0\lef_{\rm in} \le \Lambda {\mu}$ in ${\Omega}_x\times\mathbb{R}^d$, then there exists a solution $f$ to the Cauchy problem~\eqref{nk} such that $0\le f \le \Lambda {\mu}$ in $\mathbb{R}_+\times{\Omega}_x\times\mathbb{R}^d$. 
Moreover, for any ${\nu}\in(0,1)$, $k\in\mathbb{N}$, and any compact subset $K\subset(0,T]\times{\Omega}_x$, there is some constant $C_{T,{\nu},k,K}>0$ depending only on $d,{\beta},\lambda,\Lambda,T,{\nu},k,K$ and the initial data such that
\begin{equation*}\label{golbalregu}
\|{\mu}^{-{\nu}} f\|_{{\mathcal C}^{k}(K\times\mathbb{R}^d)} \le C_{T,{\nu},k,K}.  
\end{equation*}
Additionally, if $f_{\rm in}$ is H\"older continuous and $\rho_{f_{\rm in}}\ge\lambda$ in ${\Omega}_x$, then the solution that lies below any Maxwellian is unique. 
			
\medskip
\item [\bf {\;(ii)}]
For ${\Omega}_x=\mathbb{T}^d$, if the initial data satisfies $\lambda{\mu}\lef_{\rm in}\le \Lambda{\mu}$ in $\mathbb{T}^d\times\mathbb{R}^d$, 
then for any $k\in\mathbb{N}$, there exists some constant $c>0$ depending only on $d,{\beta},\lambda,\Lambda$ and some constant $C_k>0$ depending additionally on $k$ such that for any $t\ge1$,
\begin{equation*}\label{longtimethm}
\left\|\frac{f-{\mu}\int f_{\rm in}\mathop{}\!\mathrm{d} x\mathop{}\!\mathrm{d} v}{\sqrt{{\mu}}}\right\|_{{\mathcal C}^k(\mathbb{T}^{d}\times\mathbb{R}^d)} 
\le C_k e^{-ct}. 
\end{equation*}
For ${\Omega}_x=\mathbb{R}^d$, if the initial data satisfies $\lambda{\mu}\lef_{\rm in}\le \Lambda{\mu}$ in $\mathbb{R}^d\times\mathbb{R}^d$ and $f_{\rm in}-M_1{\mu}\in L^1(\mathbb{R}^d\times\mathbb{R}^d)$ for some constant $M_1>0$, then there is some constant $C'>0$ depending only on $d,{\beta},\lambda,\Lambda,M_1$ such that
\begin{equation*}
\left\|\frac{f-M_1{\mu}}{\sqrt{{\mu}}}\right\|_{L^2(\mathbb{R}^d\times\mathbb{R}^d)}
\le C'\left(1+\|f_{\rm in}-M_1{\mu}\|_{L^1(\mathbb{R}^d\times\mathbb{R}^d)}\right) t^{-\frac{d}{4}}. 
\end{equation*}
\end{item}
\end{theorem}


\begin{remark}\label{wellposer}
For general measurable initial data $f_{\rm in}$, if it satisfies $f_{\rm in} \le \Lambda \mu$ and an extra locally uniform lower bound assumption (see \eqref{lowerdd} below for a precise description), the existence of solutions still holds in some weak sense, as pointed out in Remark~\ref{remark-weak} below. 
\end{remark}


In order to describe the diffusion asymptotics of the equation~\eqref{epfnk}, we introduce the (Bregman) distance characterized by the relative phi-entropy functional $\mathcal{H}_{\beta}$. 
\begin{definition}\label{phi-entropy}
Let ${\beta}\in[0,1]$. For any measurable functions $h_1\ge0$ and $h_2>0$ defined in $\mathbb{T}^d\times\mathbb{R}^d$, the \emph{relative phi-entropy} of $h_1$ with respect to $h_2$ is defined by
\begin{align*}
\mathcal{H}_{\beta}(h_1|h_2):=&
\int_{\mathbb{T}^d\times\mathbb{R}^d}\left( \varphi_{\beta}(h_1)-\varphi_{\beta}(h_2)-\varphi'_{\beta}(h_2)(h_1-h_2) \right)\mathop{}\!\mathrm{d} m, 
\end{align*}
where $\varphi_{\beta}:\mathbb{R}_+\rightarrow\mathbb{R}$ is defined by $\varphi_{\beta}(z):=\frac{1}{1-{\beta}}\left(z^{2-{\beta}}-(2-{\beta})z+1-{\beta}\right)$ for ${\beta}\in[0,1)$ and $\varphi_1(z):=z\log z-z+1$. 
\end{definition}
	
\begin{theorem}\label{limit}
Let the constants ${\alpha}_0\in(0,1)$, $0<\lambda<\Lambda$ be given and consider a sequence of functions  $\left\{f_{{\epsilon},\rm in}\right\}_{{\epsilon}\in(0,1)}\subset 
{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)$ satisfying $0\le f_{{\epsilon},\rm in} \le \Lambda{\mu}$ in $\mathbb{T}^d\times\mathbb{R}^d$ and ${\rho}_{f_{{\epsilon},\rm in}}\ge\lambda$ in $\mathbb{T}^d$. 
Let $f_{\epsilon}$ be the solution to \eqref{epfnk} associated with the initial data $f_{{\epsilon}, \rm in}$. 
\begin{item}
\item [\bf {\ \,(i)}]
If there exists some constant ${\epsilon}'\in(0,1)$ and some function $\rho_{\rm in}\in {\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)$ valued in $[\lambda,\Lambda]$ such that
\begin{equation*}
	\mathcal{H}_{\beta}\left(\,{\mu}^{-1}f_{{\epsilon},\rm in}\,|\,\rho_{\rm in}\,\right)\le {\epsilon}', 
\end{equation*}
then there exist some constants $M,m>0$ depending only on $d,{\beta},\lambda,\Lambda,{\alpha}_0, \left\|{\rho}_{\rm in}\right\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}$ and $\left\|f_{{\epsilon},\rm in}\right\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)}$ such that for any $T>0$, 
\begin{equation*}\label{limit1}
\left\|{\mu}^{-1}f_{\epsilon}-\rho\right\|_{L^\infty\left([0,T];\,L^2\left(\mathbb{T}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m\right)\right)}
\le Me^{MT}({\epsilon}+{\epsilon}')^{m},  
\end{equation*}
where $\rho(t,x)$ for $(t,x)\in\mathbb{R}_+\times\mathbb{T}^{d}$ is the solution to the following fast diffusion equation, 
\begin{equation}\label{limiteq}
\left\{ 
\begin{aligned}
\ &\partial_t\rho(t,x)=\nabla_x\cdot\left(\rho^{-{\beta}}(t,x)\nabla_x\rho(t,x)\right), \\
\ &\;\rho(0,x) = \rho_{\rm in} (x).  \\
\end{aligned}
\right. 
\end{equation} 

	
\medskip
\item [\bf {\;(ii)}]
If we additionally assume that $f_{{\epsilon},\rm in}\ge\lambda{\mu}$ in $\mathbb{T}^d\times\mathbb{R}^d$, 
then there exist some constants $M',m'>0$ with the same dependence as $M,m$ such that 
\begin{equation*}\label{limit2}
\left\|{\mu}^{-1}f_{\epsilon}-\rho\right\|_{L^\infty\left(\mathbb{R}_+;\,L^2\left(\mathbb{T}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m\right)\right)}
\le M'({\epsilon}+{\epsilon}')^{m'}. 
\end{equation*}
\end{item}
\end{theorem}



\subsection{Strategy and background}\label{background}

\subsubsection{Cauchy problem of the nonlinear model}
The study of the well-posedness of the nonlinear model~\eqref{nk} was firstly addressed
in \cite{IM} mixing H\"older and Sobolev spaces in the torus, and in \cite{LWYang} under the regime of perturbation to the 
global equilibrium in the whole space. 
We develop it with rough initial data by means of the combination of hypoelliptic analog of the parabolic theory with a positivity-spreading result; 
in particular, the technique we employ allows us to drop the smallness and lower bound assumptions as asserted in Theorem~\ref{wellpose}. 
Besides, the global behavior of solutions to \eqref{nk} is derived under the assumption of upper and lower bounds on the initial data only. 

When the drift-diffusion coefficient $\rho^{\beta}_f$ in \eqref{nk} is proportional to the local mass of the solution, that is when ${\beta}=1$, the equations \eqref{nk}, \eqref{epfnk} have the same quadratic homogeneity as the Landau equation, but simpler global bounds and conservation laws. 
Due to the complex structure of the Landau equation, most of the existing results for its classical solutions are about the global theory under the near Maxwellian equilibrium regime \cite{Guo,KGH}, and about the local well-posedness associated with low-regularity and non-perturbative initial data \cite{HST2,HST1}. 
By contrast, the boundedness from above and from below by Maxwellians of the initial data will be preserved along 
time for the solutions to \eqref{nk}, \eqref{epfnk}, 
and the lack of conservation of momentum and energy of \eqref{epfnk} reduces its hydrodynamic limit to the fast diffusion 
flow~\eqref{limiteq} rather than the Navier-Stokes dynamics of scaling limit of the Landau equation, 
which makes its Cauchy problem and global behavior more tractable in a very general setting. 


The method we propose to address the nonlinear Cauchy problem with the only requirement of the initial data lying below a Maxwellian involves several ingredients. 
First of all, in Section~\ref{linearsection} we carry out a preliminary study on the linear counterpart of \eqref{nk}, that is the Cauchy problem associated to the Kolmogorov operator, 
\begin{equation}\label{k}
\mathscr{L}_1:=\partial_t+v\cdot\nabla_x-{\text{\rm tr}}\left(A(t,x,v)D_v^{2}\cdot\right)+B(t,x,v)\cdot\nabla_v,
\end{equation} 
where the coefficients including the entries of the positive definite $d\times d$ real symmetric matrix $A$ and the 
$d$-dimensional vector $B$ are H\"older continuous ($B$ is not necessarily bounded over $v\in\mathbb{R}^d$). 
Even if the well-posedness theory for the Cauchy problem associated to the linear operator~\eqref{k} was already 
well-developed in some sense in the existing literature (see \cite{M}, as well as the survey paper \cite{APsurvey} and the references therein), 
the H\"older spaces (see Definition~\ref{kineticholder}) considered in these works are instead different from the one studied in \cite{ISkinetic,IM} (see Definition~\ref{holdercontinuous}) which we will use. 
Indeed, in contrast to \cite{IM}, the (Schauder type) apriori estimates proved in the previous literature are weaker and not appropriate to bootstrap higher regularity for nonlinear problems (see Subsection~\ref{smoothsection}). 


Secondly, the treatment of the existence issue for \eqref{nk} in H\"older spaces is based on a fixed point argument, where the compactness is provided by hypoelliptic regularization results; see Subsection~\ref{wellsection}. 
A breakthrough on such a priori estimates for spatially inhomogeneous kinetic equations with a quasilinear diffusive structure in velocity was obtained in the works \cite{GIMV} and \cite{HS,IM}, where the authors proved the kinetic (hypoelliptic) counterparts of the De Giorgi-Nash-Moser theory and the Schauder theory for classical elliptic equations (see for instance \cite{GilbargTrudinger}), respectively. One may refer to \cite{Mouhot} for a summary. 
Armed with the Schauder estimate developed in \cite{IM} in kinetic H\"older spaces and the bootstrap procedure developed in \cite{IS} adapted to our case, we are then able to derive instantaneous ${\mathcal C}^\infty$ regularization for the solutions to the equation~\eqref{nk} in Subsection~\ref{smoothsection}, provided that the solution is bounded from above and bounded away from vacuum, which guarantees the ellipticity in the velocity variable for \eqref{nk}. 


Thirdly, in order to remove the lower bound assumption on the initial data, in Subsection~\ref{lowersection} we establish a self-generating lower bound result showing that the positivity of solutions spreads everywhere instantaneously. 
Its proof is based on repeated applications of the spreading of positivity forward in time (see Lemma~\ref{lowerlemma1}) and the spreading for all velocities (see Lemma~\ref{lowerlemma2}), as proposed in \cite{HST3}. 
On the one hand, the barrier function argument will be used in the same spirit as \cite{HST3} to show the former one. 
Indeed, a \emph{lower (resp.\! upper) barrier} for a certain equation is a  subsolution (resp.\! supersolution) of the equation which bounds its solution from below (resp.\! above) on the boundary; it then follows from the maximum principle that the barrier function performs as a lower (resp.\! upper) bound of the solution. 
On  the  other  hand, combining the local Harnack inequality obtained in \cite{GIMV} with the construction of a Harnack chain yields the latter one. We remark that the idea of the Harnack chain was firstly used in \cite{Moser} and an example of its application to Kolmogorov equations can be found in \cite{AEP}. 
Essentially, the spreading of positivity can be seen as a lower bound estimate of the fundamental solution, which is thus related to the result in \cite{HST2}, where they applied a probabilistic method. 

A subtle point of the lower bound result lies in the possibilities of the degeneracy of solutions as $t\rightarrow0^+$ or 
$t\rightarrow\infty$, which leads to two delicate issues. 
First, with the same difficulty as mentioned in \cite{HST1}, in order to prove the uniqueness of the Cauchy problem~\eqref{nk}
the nondegeneracy of diffusion up to the initial time is required so that the a priori estimates can be still applicable. 
We remark that, generally speaking, deriving uniqueness of solutions to nonlinear equations in rough spaces is always a classical difficulty, and the presence of a vacuum sometimes gives rise to non-uniqueness phenomenon even for the limiting equation \eqref{limiteq} (see for instance \cite{DK}). 
Under the additional assumptions of absence of vacuum on the local mass and H\"older continuity of initial data, we achieve the uniqueness by using the scaling argument and Grönwall's lemma, since the H\"older estimate around the initial time implies that the integrand in the inequality of Grönwall's type is improved to be integrable with respect to the time variable; see the proof of Theorem~\ref{unique} below for more details. 
Second, we are only able to show the convergence to equilibrium if the drift-diffusion coefficient $\rho^{\beta}_f$ decays slower than $t^{-1}$ as $t\rightarrow\infty$ in Proposition~\ref{longtime}. Therefore, an additional lower Maxwellian bound on the initial data is imposed in Theorem~\ref{wellpose}~(ii) and in Theorem~\ref{limit}~(ii) to ensure the solutions will be away from the vacuum uniformly along time. 
It would be expected that such additional lower bound assumption could be removed, especially when ${\beta}$ is small. 


\subsubsection{Long time behavior}
The drift-diffusion operator $\mathscr{L}_{\scriptscriptstyle{\rm FP}}$ acts only on the velocity variable and ceases to be dissipative on its unique steady state ${\mu}$, which also ensures that the null space of $\mathscr{L}_{\scriptscriptstyle{\rm FP}}$ is spanned by ${\mu}$ and the conservation law of mass is satisfied. Consequently, the convergence to equilibrium is to be expected.  
With the help of the global smoothness a priori estimates, we are able to pass from the exponential convergence to equilibrium in $L^2$-framework to the uniform convergence in ${\mathcal C}^\infty$ in Subsection~\ref{longtimesection}, when the spatial domain is compact, that is the periodic box $\mathbb{T}^d$. 
Therein, the $L^2$-convergence is obtained by the $L^2$-hypocoercivity under a macro-micro (fluid-kinetic) decomposition scheme, which suggests to construct some proper entropy (Lyapunov) functional that would provide an equivalent $L^2$-norm for solutions. 
The key ingredient is to control the macroscopic part by means of the microscopic part in view of the decomposition. 
This hypocoercive theory was studied in \cite{EGKM,DMS,He} via different approaches, while their ideas are essentially the same. 
In \cite{EGKM}, the authors intended to develop the nonlinear energy estimate in an $L^2$-$L^\infty$ framework. 
In \cite{DMS} and \cite{He}, the authors studied the $L^2$-hypocoercivity theory in an abstract setting and in the framework of pseudo-differential calculus, respectively.  
Besides, if the spatial domain is $\mathbb{R}^d$, meaning that it is not confining in a compact region, then the convergence rate slows down to an algebraic decay, whose hypocoercive theory was captured by \cite{BDM}. We remark that the $L^2$-framework allows us to avoid some difficulties from the nonlinearity of the operator $\rho^{\beta}_f\mathscr{L}_{\scriptscriptstyle{\rm FP}} f$, in contrast with $H^1$-entropic hypocoercivity methods (see for instance the memoir~\cite{Villanihypo}). 


\subsubsection{Diffusion asymptotics}
The diffusion approximation serves as a simplification of collisional kinetic equations when the mean free path is much smaller than the typical length of 
observation in a long time scale. This approximation for linear Fokker-Planck models can be traced back to \cite{Degond}, where the authors applied the Hilbert expansion method. 
One is also able to achieve the diffusion limit for \eqref{epfnk} in some weak sense by applying a similar strategy to the one given in \cite{ELM}. However, weak convergence is sometimes ineffective for application, as precise description of the convergence is not given. 
Still the nonlinearity of the term $\rho^{\beta}_{f_{\epsilon}}\mathscr{L}_{\scriptscriptstyle{\rm FP}} f_{\epsilon}$ in \eqref{epfnk} associated with non-perturbative initial data reveals some difficulties to derive a quantitative convergence. 

In order to overcome this difficulty, in Subsection~\ref{finitetimesection} we will rely on the phi-entropy of solutions relative to their limit to see the finite time asymptotics on the torus.
The relative entropy method, that heavily relies on the regularity of solutions to the target equation,  has become an effective tool in the study of hydrodynamic limits since \cite{YauHZ,BGL} (see also \cite{LSR}). 
The method applied to the diffusion asymptotics of the kinetic Fokker-Planck equation of the type with linear diffusion can be found in \cite{Markou}. 
The so-called phi-entropy (relative to the global equilibrium) was used to study the convergence of certain kinds of Fokker-Planck equations; see for instance \cite{AMTU,DL}. 
Finally, combining the barrier function method with a careful treatment of the regularity estimate of the target equation enables us to deal with the asymptotic dynamics for the cases associated with general H\"older continuous initial data.  


	
\subsection{Physical motivation}
The spatially inhomogeneous Fokker-Planck equation~\eqref{nk} arises from modeling the evolution of some system of a large number of interacting particles from the statistical mechanical point of view. 
These models appear for instance in the study of plasma physics and biological dynamics; see \cite{Villani,Cha2}. 
Its solution can be interpreted as the probability density of the particles lying at the position $x$ at time $t$ with velocity $v$. The scaled model \eqref{epfnk} for small ${\epsilon}$ describes the evolution of the particle density in the small mean free path and long time regime, where the nondimensional parameter ${\epsilon}\in(0,1)$ designates the ratio between the mean free path (microscopic scale) and the typical macroscopic length. 
The limiting equation~\eqref{limiteq} characterizes its macroscopic dynamics. 

From the perspective of a stochastic process $\{(X_t,V_t):t\ge0\}$ driven by a Brownian motion $\{\mathcal{B}_t:t\ge0\}$, 
\begin{equation*}
\left\{ 
\begin{aligned}
\ &\mathop{}\!\mathrm{d} X_t= V_t\mathop{}\!\mathrm{d} t, \\   
\ &\mathop{}\!\mathrm{d} V_t= {\rho}^{\beta}_f(t,X_t)V_t\mathop{}\!\mathrm{d} t+\sqrt{2{\rho}^{\beta}_f(t,X_t)}\mathop{}\!\mathrm{d}\mathcal{B}_t, \\
\end{aligned}
\right. 
\end{equation*} 
the dual equation describing the dynamics of $\{(X_t,V_t):t\ge0\}$ is given by \eqref{nk}; see the review paper \cite{Chand}. 
Indeed, the nonlinear term $\rho^{\beta}_f\mathscr{L}_{\scriptscriptstyle{\rm FP}} f$ models the collisional interaction of the particles, where the mobility of these particles is hampered by their aggregation. 
More precisely, the nonlinear dependence on the drift-diffusion coefficient $\rho^{\beta}_f$ translates the fact that the friction effect in the interaction is positively correlated to the local mass of particles occupying the position $x$ at time $t$. 
Moreover, the low field scaling $t\mapsto{\epsilon}^2 t$, $x\mapsto{\epsilon} x$ of the equation~\eqref{nk} formally implies the equation~\eqref{epfnk}. As ${\epsilon}$ tends to zero, its spatial diffusion phenomena are characterized by the equation~\eqref{limiteq}. 

Regarding its physical interpretation, we point out that the factor multiplying the time derivative in \eqref{epfnk} takes into account the long time scale. The inverse of the factor multiplying $\rho^{\beta}_{f_{\epsilon}}\mathscr{L}_{\scriptscriptstyle{\rm FP}} f_{\epsilon}$ stands for the scaled average distance traveled by particles between each collision and it is usually referred to as mean free path. 
In the small mean free path regime, it was noticed in \cite{Chand} that the spatial variation occurs significantly only under the long time scale that is consistent with the particle motion. 
In such an overdamped process, named as low field limit or diffusion limit, the kinetic model is reduced to a macroscopic model. 

	
Finally, we recall that the associated phi-entropy introduced in Definition~\ref{phi-entropy} is also known as Tsallis entropy in the physics community, which generalizes the Boltzmann–Gibbs entropy (the  phi-entropy with $\beta=1$) in nonextensive statistical mechanics~\cite{tsallis1}.  
It gives some hints for the formulation of the correlated diffusion, where the index $\beta$ measures the degree of nonextensivity and nonlocality of the system; see \cite{tsallisbook} and the references therein. 

	
\subsection{Organization of the paper}
The article is organized as follows. In Section~\ref{preliminaries}, we recall some basic notions related to kinetic H\"older spaces that are adapted to the Fokker-Planck equations. 
Section~\ref{linearsection} is devoted to the study of the linear Fokker-Planck equation with H\"older continuous coefficients. 
The well-posedness result Theorem~\ref{wellpose}~(i) is proved in Section~\ref{nonlinear}. 
The asymptotic behaviors, including Theorem~\ref{wellpose}~(ii) and Theorem~\ref{limit}, are proved in Section~\ref{diffusionlimit}. 
	
\medskip\noindent\textbf{Acknowledgement.}
The authors are grateful to Cyril Imbert for suggesting the question, and both François Golse and Cyril Imbert for the helpful discussions, and the anonymous referees for their careful reading and comments.
\includegraphics[height=0.275cm]{EU_flag.png} 
YZ's research has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 754362. 
	
	
	
\section{Preliminaries}\label{preliminaries}
This section is devoted to the basic notations, including the invariant structure and the kinetic H\"older space for the equations we are concerned with. 
Instead of the usual parabolic scaling and translations, the  invariant scaling and transformation associated with the Kolmogorov operator $\mathscr{L}_1$ (see \eqref{k}) is replaced by kinetic scaling and Galilean transformations, respectively. 
It then turns out that the appropriate H\"older space as well as its norm should be adapted to the new scaling and transformation. 
Let us now precisely state these notions. 


\subsection{The geometry associated to Kolmogorov operators}
Let $z:=(t,x,v)\in\mathbb{R}\times\mathbb{R}^d\times\mathbb{R}^d$. 
We define the \emph{kinetic scaling}, 
\begin{equation*}
S_r(t,x,v) := (r^2 t,\, r^{3} x,\, r v), {\quad\rm for\ } r>0,  
\end{equation*}
and the \emph{Galilean transformation}, 
\begin{equation*}
(t_0,x_0,v_0)\circ(t,x,v):=(t_0+t,\,x_0+x+tv_0,\,v_0+v), 
	{\quad\rm for\ } (t_0,x_0,v_0)\in\mathbb{R}\times\mathbb{R}^d\times\mathbb{R}^d.
\end{equation*}
With respect to the product $\circ$, we are able to define the inverse of $z$ as given by $z^{-1}:=(-t,-x+tv,-v)$. 
In view of this structure of scaling and transformation, it is natural to define the cylinder centered at the origin of
radius $r>0$ as 
$$Q_r:= (-r^2, 0] \times B_{r^3}(0) \times B_r(0).$$
More generally, the cylinder  centered at $z_0=(t_0,x_0,v_0)$ with radius $r$ is defined by 
\begin{align*}
Q_r(z_{0}):=&\left\{ z_{0} \circ S_r(z): \, z \in Q_1\right\}\\
=&\left\{ (t,x,v): \, t_{0} - r^{2} < t \le t_{0}, \, |x - x_{0} - (t - t_{0})v_{0}| < r^{3},\, |v - v_{0}| < r\right\}.
\end{align*}
Roughly speaking, for fixed $z_0\in\mathbb{R}^{1+2d}$, the Kolmogorov operator $\mathscr{L}_1$ is invariant under the kinetic scaling and left invariant under the Galilean transformation. 
It means that if $f$ is a solution to the equation $\mathscr{L}_1f=0$ in $Q_{r}(z_{0})$, then $f\left(z_0\circ S_r(\cdot)\right)$ solves an equation of the same ellipticity class in $Q_1$. 

Besides, the associated quasi-norm $\|\cdot\|$ is defined by 
\begin{align*}
\|z\|:=\max\Big\{ |t|^{\tfrac12}, |x|^{\tfrac13}, |v| \Big\}, 
\end{align*}   
as we notice that $\|S_r(z)\|=r\|z\|$ and $\|z_0\circ z\|\le3\left(\|z_0\|+\|z\|\right)$. 
For further information on the non-Euclidean geometry associated to Kolmogorov operators, one may refer to \cite{APsurvey, ISkinetic} and the references therein. 

\subsection{Kinetic H\"older spaces and differential operators}
The proper kinetic H\"older space and kinetic degree of basic differential operators should be adapted to the above definitions, so that they are homogeneous under these transformations. Their definitions were introduced in \cite{ISkinetic} (see also \cite{IM}), and recalled below. 

Given a monomial $m (t,x,v) =t^{k_0} x_1^{k_1} \ldots x_d^{k_d} v_{1}^{k_{d+1}} \ldots v_d^{k_{2d}}$, we define its \emph{kinetic degree} as 
\begin{equation*}
{\rm deg}_{\rm kin}(m)= 2 k_{0} + 3 \sum\nolimits_{j=1}^{d} k_{j} + \sum\nolimits_{j=d+1}^{2d} k_{j}.
\end{equation*}
Any polynomial $p\in\mathbb{R}[t,x,v]$ can be uniquely written as a linear combination of monomials and its kinetic degree ${\rm deg}_{\rm kin}(p)$ is defined by the maximal kinetic degree of the monomials appearing in $p$. 
This definition is justified by the fact that $p(S_r(z)) = r^{{\rm deg}_{\rm kin}(p)} p(z)$. 


\begin{definition}\label{holdercontinuous}
Let the constant ${\alpha}>0$ and the open subset ${\Omega}\subset\mathbb{R}\times\mathbb{R}^d\times\mathbb{R}^d$ be given. We say a function
$f : {\Omega} \rightarrow \mathbb{R}$ is ${\mathcal C}_l^{\alpha}$-continuous at a point $z_0\in{\Omega}$, 
if there exists a polynomial $p_{0}\in \mathbb{R}[t,x,v]$ with ${\rm deg}_{\rm kin} (p_{0})<{\alpha}$ and a constant $C>0$ such that, for any $z\in{\Omega}$ with $z_0\circ z\in{\Omega}$,  
\begin{equation}\label{holdercontinuousineq}
| f(z_0\circ z) - p_{0}(z) | \le C \| z\|^{\alpha}. 
\end{equation}
If this property holds for any $z_0,z$ on each compact subset of ${\Omega}$, then we say $f \in {\mathcal C}_l^{{\alpha}}({\Omega})$. 
If the constant $C$ in \eqref{holdercontinuousineq} is uniformly bounded for any $z_0,z\in{\Omega}$, 
we define the smallest value of $C$ as the semi-norm $[f]_{{\mathcal C}_l^{{\alpha}}({\Omega})}$, 
and the norm $\|f\|_{{\mathcal C}_l^{{\alpha}}({\Omega})}:=[f]_{{\mathcal C}_l^0({\Omega})}+[f]_{{\mathcal C}_l^{{\alpha}}({\Omega})}$, where we additionally define ${\mathcal C}_l^0({\Omega}):={\mathcal C}^0({\Omega})$, the space of continuous functions on ${\Omega}$, with the norm $\|f\|_{{\mathcal C}_l^0({\Omega})}:=[f]_{{\mathcal C}_l^0({\Omega})}:=\|f\|_{{\mathcal C}^0({\Omega})}=\|f\|_{L^\infty({\Omega})}$. 
\end{definition}
\begin{remark}
For ${\alpha}\in[0,1)$, this ${\mathcal C}_l^{\alpha}$-continuity is equivalent to the standard definition of 
${\mathcal C}^{\alpha}$-continuity with respect to the distance $\|\cdot\|$.
The subscript ``\textit{l}'' of ${\mathcal C}_l$ stems from the definition of H\"older continuity above, that it is given in terms of a  left-invariant distance with respect the group structure of $\circ$. 
\end{remark}


We also mention another kind of H\"older space suitable for the study of Kolmogorov operators that was first used in \cite{M}. 
\begin{definition}\label{kineticholder}
Let ${\alpha}\in[0,1)$ and ${\Omega}\subset\mathbb{R}\times\mathbb{R}^d\times\mathbb{R}^d$ be given. 
The space ${\mathcal C}_{\rm kin}^{2+{\alpha}}({\Omega})$ consists of functions $f\in{\mathcal C}_l^0({\Omega})$ such that $D^2_vf$, $(\partial_t+v\cdot\nabla_x)f\in{\mathcal C}_l^{\alpha}({\Omega})$, equipped with the norm 
 \begin{equation*}
 \|f\|_{{\mathcal C}_{\rm kin}^{2+{\alpha}}({\Omega})}:=\|f\|_{{\mathcal C}_l^0({\Omega})}+\|D^2_vf\|_{{\mathcal C}_l^{\alpha}({\Omega})}+\|(\partial_t+v\cdot\nabla_x)f\|_{{\mathcal C}_l^{\alpha}({\Omega})}. 
 \end{equation*}
\end{definition}

The consistence between these two definitions is given by \cite[Lemma 2.7]{ISkinetic} (see also \cite[Lemma 2.4]{IM}), a result that we state here as follows. 
\begin{lemma}\label{kinetic-degree}
Let ${\alpha}\in(0,1)$ and $f\in{\mathcal C}_l^{2+{\alpha}}(Q_2)$. Then, there exists some constant $C>0$ depending only on the dimension $d$ such that  
\begin{equation*}
\|\nabla_vf\|_{{\mathcal C}_l^{{\alpha}}(Q_1)} \le C \|f\|_{{\mathcal C}_l^{1+{\alpha}}(Q_2)},\quad
\|D^2_vf\|_{{\mathcal C}_l^{\alpha}(Q_1)}+\|(\partial_t+v\cdot\nabla_x)f\|_{{\mathcal C}_l^{\alpha}(Q_1)} \le C \|f\|_{{\mathcal C}_l^{2+{\alpha}}(Q_2)}. 
\end{equation*}
\end{lemma}

\begin{remark}\label{holderremark}
	For ${\alpha}>2$, one can easily check that the polynomial $p_0$ in \eqref{holdercontinuousineq} has the form
	\begin{equation*}
	p_0(t,x,v)= f(z_0)+(\partial_t+v_0\cdot\nabla_x)f(z_0)\,t+\nabla_vf(z_0)\cdot v + \frac{1}{2} D^2_vf(z_0)\,v\cdot v+\ldots. 
	\end{equation*}
	In particular, if ${\alpha}\in(2,3)$, the polynomial expansion is independent of the $x$-variable. 
\end{remark}

\begin{remark}
	A subtle difference between ${\mathcal C}_l^2$ and ${\mathcal C}_{\rm kin}^2$ comes from the fact that for $f\in{\mathcal C}_l^2$, $D^2_vf$, $(\partial_t+v\cdot\nabla_x)f$ are lying in $L^\infty$ rather than ${\mathcal C}^0$. 
\end{remark}

We will also employ the following notions of weighted H\"older norms in Section~\ref{linearsection}. 

\begin{definition}\label{weighted-norm}
Let $z=(t,x,v)\in{\Omega}:=(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$ with $T\in\mathbb{R}_+$. 
For $f\in{\mathcal C}_l^{{\alpha}}({\Omega})$ with ${\alpha}>0$ and ${\sigma}\in\mathbb{R}$, we define 
\begin{equation*}
[f]_{0}^{(\sigma)}:=
\sup\nolimits_{z\in{\Omega}}\iota^{\sigma}[f]_{{\mathcal C}_l^0(Q_\iota(z))},\ \ 
[f]_{{\alpha}}^{(\sigma)}:=
\sup\nolimits_{z\in{\Omega}}\iota^{{\alpha}+\sigma}[f]_{{\mathcal C}_l^{\alpha}(Q_\iota(z))},\ \ 
\|f\|_{\alpha}^{(\sigma)}:= [f]_0^{(\sigma)}+[f]_{\alpha}^{(\sigma)},  
\end{equation*} 
where $\iota:=\min\big\{1,t^\frac{1}{2}\big\}$ measures the distance between $z$ and the (parabolic) boundary of ${\Omega}$.
\end{definition}


\subsection{Other notations}
Throughout the article, $B_R$ denotes the Euclidean ball in $\mathbb{R}^d$ centered at the origin with radius $R>0$. 
We employ the Japanese bracket defined as $\langle\cdot\rangle:=\left(1+|\cdot|^2\right)^\frac{1}{2}$. 
By abuse of notation, $\langle\cdot\rangle$ will also denote the velocity mean in Section~\ref{diffusionlimit}. 

Moreover, we assume $0<\lambda<\Lambda$. 
We denote by $C$ a \emph{universal} constant, that is to say a constant depending only on ${\beta},d,\lambda,\Lambda,{\alpha},{\sigma},{\alpha}_0$ specified in context. 
Finally, we write $X\lesssim Y$ to denote that $X\le CY$ for some universal constant $C>0$, 
and $X\lesssim_{q} Y$ to denote that $X\le C_qY$ for $C_q>0$ depending only on universal constants and the quantity $q$. 


	
	
	
\section{Kolmogorov-Fokker-Planck equation}\label{linearsection}
This section is devoted to the study of the Cauchy problem associated to the operator~\eqref{k}, 
\begin{equation}\label{linearequation}
\left\{ 
\begin{aligned}
\ &\mathscr{L}_1 f :=\left(\partial_t+v\cdot\nabla_{x}\right)f-{\text{\rm tr}}(AD_v^{2}f)-B\cdot\nabla_{v}f=s \quad{\rm in\ } (0,T] \times\mathbb{R}^{d}\times\mathbb{R}^{d}, \\
&\; f(0,x,v)=f_{\rm in}(x,v) \quad {\rm in\ } \mathbb{R}^{d}\times\mathbb{R}^{d}, \\
\end{aligned}
\right. 
\end{equation} 
where the $d\times d$ symmetric matrix $A(t,x,v)$ and the $d$-dimensional vector $B(t,x,v)$ satisfy the following condition, 
\begin{equation}\label{schaudercondition}
\left\{ 
\begin{aligned}
\ & A\,\xi\cdot\xi\ge\lambda|\xi|^2{\quad\rm\ for\ any\ }\xi\in\mathbb{R}^d,\\
& \|A\|_{{\mathcal C}_l^{\alpha}} +\|B\|_{{\mathcal C}_l^{\alpha}}\le\Lambda, \\
\end{aligned}
\right. 
\end{equation} 
where ${\alpha}\in(0,1)$ and the norm $\|\cdot\|_{{\mathcal C}_l^{\alpha}({\Omega})}$ of matrix denotes the summation of the norm of each entry. The boundedness condition at infinity means that the solution shall be bounded, which is intended for the validity of maximum principle; see the proof of Lemma~\ref{max} below.

The aim of this section is to solve the Cauchy problem~\ref{linearequation} by virtue of the weighted H\"older norm (Definition~\ref{weighted-norm}) and by means of the standard continuity method combined with Schauder type estimates. 
One may refer to \cite[Subsection 6.5]{GilbargTrudinger} for the corresponding treatment in classical elliptic theory.  

Throughout this section we work with the domain ${\Omega}:=(0,T] \times\mathbb{R}^{d} \times\mathbb{R}^{d}$, with $T\in\mathbb{R}_+$. We shed light on the fact that all of the results below can be restricted to $(0,T] \times \mathbb{T}^{d} \times\mathbb{R}^{d}$ whenever required. 




\subsection{Schauder estimates}
In order to apply the continuity method, first of all one needs to prove a global a priori estimate for solutions to \eqref{linearequation} with respect to the weighted H\"older norm. 
In the kinetic setting, we have at our disposal the following interior Schauder estimates proved in \cite[Theorem 3.9]{IM}. 


\begin{proposition}[Interior Schauder estimate]\label{interiorschauder} 
Let the constant ${\alpha}\in(0,1)$ be given and the cylinder $Q_{2r}(z_0)\subset{\Omega}$ with $r\in(0,1]$. 
If $f$ verifies the equation~\eqref{k} with the condition~\eqref{schaudercondition} in $Q_{2r}(z_0)$ and $s\in{\mathcal C}_l^{\alpha}(Q_{2r}(z_0))$, then we have 
\begin{equation}\label{schauderIM}
r^{2+{\alpha}} [f]_{{\mathcal C}_l^{2+{\alpha}}(Q_r(z_0))}
\lesssim \|f\|_{L^\infty(Q_{2r}(z_0))}+ r^{2+{\alpha}} [s]_{{\mathcal C}_l^{{\alpha}}(Q_{2r}(z_0))}; 
\end{equation}
in particular, the right hand side controls $r^2\|\left(\partial_t+v\cdot\nabla_x\right)f\|_{L^\infty(Q_r(z_0))} +r^2\|D_v^2f\|_{L^\infty(Q_r(z_0))}$. 
\end{proposition}

First of all, we enhance this result to a global estimate for the Cauchy problem \eqref{linearequation} under a vanishing condition for the initial data. 


\begin{proposition}[Global Schauder estimate]\label{globalschauder} 
Let ${\Omega}=(0,T] \times\mathbb{R}^{d} \times\mathbb{R}^{d}$, and the constants ${\alpha}\in(0,1)$, ${\sigma}\in(0,2)$ be universal, $s\in{\mathcal C}_l^{\alpha}({\Omega})$ such that $\|s\|_{\alpha}^{(2-\sigma)}<\infty$ and $f$ be a bounded solution to the Cauchy problem~\eqref{linearequation} under the condition \eqref{schaudercondition} in ${\Omega}$. 
If the initial data $f_{\rm in}=0$, then we have 
\begin{equation*}
\|f\|_{2+\alpha}^{(-\sigma)}\lesssim\|s\|_{\alpha}^{(2-\sigma)}. 
\end{equation*}
\end{proposition}
	
\begin{proof}
In view of Proposition~\ref{interiorschauder}, it suffices to deal with the estimates around the initial time. Without loss of generality, we assume $T\le1$. 

Let $z_0=(t_0,x_0,v_0)\in{\Omega}$ and $2r=t_0^\frac{1}{2}$. 
Applying the interior Schauder estimate~\eqref{schauderIM} yields 
\begin{equation*}
r^{2+\alpha}[f]_{{\mathcal C}_l^{2+{\alpha}}(Q_{r}(z_0))}
\lesssim \|f\|_{L^\infty(Q_{2r}(z_0))}+ r^{2+\alpha}[s]_{{\mathcal C}_l^\alpha(Q_{2r}(z_0))}. 
\end{equation*}
It then follows from the arbitrariness of $z_0$ that, for any ${\sigma}<2$ such that $ [f]_{0}^{(-{\sigma})}<\infty$, 
\begin{equation}\label{schauderIM0}
[f]_{2+{\alpha}}^{(-{\sigma})} 
\lesssim [f]_{0}^{(-{\sigma})}  + [s]_{{\alpha}}^{(2-{\sigma})}. 
\end{equation}
With ${\sigma}\in(0,2)$, observing that 
\begin{align*}
\mathscr{L}_1\left(\frac{2}{{\sigma}}[s]_{0}^{(2-{\sigma})}t^{\frac{{\sigma}}{2}}\pm f\right)=[s]_{0}^{(2-{\sigma})}t^{\frac{{\sigma}}{2}-1}\pm s\ge0 
{\quad\rm in\ }{\Omega}\\
 \frac{2}{{\sigma}}[s]_{0}^{(2-{\sigma})}t^\frac{{\sigma}}{2}\pm f=0 {\quad\rm on\ }\{t=0\}&, 
\end{align*} 
we apply the maximum principle (Lemma~\ref{max}) to the function $\frac{2}{{\sigma}}[s]_{0}^{(2-{\sigma})}t^{\frac{{\sigma}}{2}}\pm f$ 
to deduce that 
$\pm t^{-\frac{{\sigma}}{2}} f\le \frac{2}{{\sigma}}[s]_{0}^{(2-{\sigma})},$ 
which means 
\begin{equation*}
[f]_{0}^{(-\sigma)}\lesssim [s]_{0}^{(2-\sigma)}.
\end{equation*} 
Combining this estimate with \eqref{schauderIM0}, we get the desired result. 
\end{proof}
	
	
	
\subsection{Cauchy problem for the linear equation}
The goal of this subsection is to prove the well-posedness of the Cauchy problem~\eqref{linearequation} with H\"older continuous coefficients. 
\begin{proposition}\label{variable}
Let ${\Omega}=(0,T] \times\mathbb{R}^{d} \times\mathbb{R}^{d}$, and the constants $\alpha\in(0,1)$, ${\sigma}\in(0,2)$ be universal. Assume that 
\begin{equation}\label{schaudercondition'}
\left\{ 
\begin{aligned}
\ & A\,\xi\cdot\xi\ge\lambda|\xi|^2{\quad\rm\ for\ any\ }\xi\in\mathbb{R}^d,\\
& \|A\|_{{\mathcal C}_l^{\alpha}({\Omega})} +\|\langle v\rangle^{-1}B\|_{{\mathcal C}_l^{\alpha}({\Omega})}\le\Lambda.\\
\end{aligned}
\right. 
\end{equation}

\noindent
Then, for any $s\in{\mathcal C}_l^{\alpha}({\Omega})$ such that $\|s\|_{\alpha}^{(2-\sigma)}<\infty$ and $f_{\rm in}\in {\mathcal C}^0(\mathbb{R}^d\times\mathbb{R}^d)$, 
there exists a unique bounded solution $f\in {\mathcal C}_l^{2+{\alpha}}({\Omega})$ to the Cauchy problem~\eqref{linearequation}. 
\end{proposition}

\begin{remark}
In contrast with \eqref{schaudercondition}, the condition~\eqref{schaudercondition'} is weaker, which allows the coefficients of $B$ to not necessarily be bounded globally. This fact will be applied to the Ornstein-Uhlenbeck operator $\mathscr{L}_{\scriptscriptstyle{\rm OU}}=\left(\nabla_v-v\right)\cdot\nabla_v$ in Subsection~\ref{wellsection}. 
\end{remark}


The simplest possible setting of \eqref{linearequation} under the condition~\eqref{schaudercondition'} is recovered by choosing $A=I$ and $B=0$, which turns out to be the classical Kolmogorov operator $\mathscr{L}_0:=\partial_t+v\cdot\nabla_x-\Delta_v$. 
Such operator was first studied in \cite{K1}, where its fundamental solution was calculated explicitly, 
\begin{equation}\label{fun-sol}
{\Gamma}(t,x,v) = 
\begin{cases}
\ \left( \frac{\sqrt{3}}{2 \pi t^{2}} \right)^{d} 
 e^{-\frac{3\left|x+\frac{t}{2}v\right|^{2}}{t^{3}}-\tfrac{|v|^{2}}{4t}}&{\rm for\ }t>0,\\
\ \ 0 &{\rm for\ }t\le 0.
\end{cases}
\end{equation}
One can easily see that ${\Gamma}$ is smooth outside of its pole (the origin). 
In fact, in this latter case the following result holds true.

\begin{lemma}\label{constant}
Let ${\Omega}=(0,T] \times\mathbb{R}^{d} \times\mathbb{R}^{d}$ and $\alpha\in(0,1)$. For any $s\in{\mathcal C}_l^{\alpha}({\Omega})$ such that $\|s\|_{\alpha}^{(2-\sigma)}<\infty$, the function 
\begin{equation}\label{representation}
f(t,x,v)=\int_{\mathbb{R}^d\times\mathbb{R}^d} {\Gamma}\left(({\tau},{\xi},{\eta})^{-1}\circ(t,x,v)\right)s({\tau},{\xi},{\eta}) \mathop{}\!\mathrm{d}{\tau} \mathop{}\!\mathrm{d}\xi\mathop{}\!\mathrm{d}{\eta}
\end{equation}
is the unique bounded solution in ${\mathcal C}_l^{2+{\alpha}}({\Omega})$ to \eqref{linearequation} with $\mathscr{L}_1$ replaced by $\mathscr{L}_0$ and $f_{\rm in}=0$. 
\end{lemma} 
\begin{remark}
When the spatial domain is $\mathbb{T}^d$, one can apply the Green function 
${G}(t,x,v):=\sum\nolimits_{n\in{\mathbb{Z}}^d}{\Gamma}(t,x+n,v)$, which is well-defined due to the decay of ${\Gamma}$. 
\end{remark}

We are now in a position to apply the standard continuity method to derive Proposition~\ref{variable}. 
\begin{proof}[Proof of Proposition~\ref{variable}]
We split the proof in three steps. 
In the first step, we establish the case for vanishing initial data under the stronger assumption~\eqref{schaudercondition}. 
We point out that the assumption on the coefficient $B$ can be weakened in the second step. 
Finally, we deal with general continuous initial data. 


\medskip\noindent{\textbf{Step 1.}}	
Assume that $f_{\rm in}=0$ and the condition~\eqref{schaudercondition} holds. 
Let the constant ${\sigma}\in(0,2)$ be fixed and consider the Banach space 
$\mathcal{Y}:=\left({\mathcal C}_l^{2+\alpha}({\Omega}),\|\cdot\|_{2+\alpha}^{(-\sigma)}\right)$. 
In particular, every function lying in $\mathcal{Y}$ vanishes at $t=0$. 

For ${\tau} \in [0,1]$, we define the operator
$\mathscr{L}_{\tau}:=(1-{\tau})\mathscr{L}_0+{\tau}\mathscr{L}_1$, which can be written in the form of 
\begin{equation*}
\mathscr{L}_{\tau}=\partial_t+v\cdot\nabla_{x}-{\text{\rm tr}}(A_{\tau} D_v^{2}\cdot)-{\tau} B\cdot\nabla_{v},
\end{equation*}
where its coefficients $A_{\tau}:=(1-{\tau})I+{\tau} A$ and ${\tau} B$ still satisfy the condition~\eqref{schaudercondition} (with $\lambda$, $\Lambda$ replaced by $\min\{1,\lambda\}$, $\max\{1,\Lambda\}$, respectively). 
For any $w\in \mathcal{Y}$, we have
\begin{equation}\label{contraction0}
\|\mathscr{L}_{\tau} w\|_{\alpha}^{(2-\sigma)}
\lesssim \left(1+\|A_{\tau}\|_{\alpha}^{(2)}\right)\|w\|_{2+\alpha}^{(-\sigma)}
+\|B\|_{\alpha}^{(2)}\|w\|_{1+\alpha}^{(-\sigma)}
\lesssim\|w\|_{2+\alpha}^{(-\sigma)}.
\end{equation}

Let the set $\mathcal{I}$ be the collection of ${\tau}\in[0,1]$ such that
the Cauchy problem~\eqref{zerocauchy} is solvable for any $s\in {\mathcal C}_l^\alpha({\Omega})$ with $\|s\|_{\alpha}^{(2-\sigma)}<\infty$: 
there is a unique bounded solution $f\in \mathcal{Y}$ verifying 
\begin{equation}\label{zerocauchy}
\left\{ 
\begin{aligned}
\ & \mathscr{L}_{\tau} f=s  {\quad\rm in\ } {\Omega}, \\
\ &\;f(0,x,v)=0 {\quad\rm in\ } \mathbb{R}^{d}\times\mathbb{R}^{d}. 
\end{aligned}
\right. 
\end{equation}  
By Lemma~\ref{constant}, we see that $0\in\mathcal{I}$; in particular, $\mathcal{I}$ is not empty. 

It now suffices to show that $1\in\mathcal{I}$.
Pick ${\tau}_0\in\mathcal{I}$, then the global Schauder estimate provided by Proposition~\ref{globalschauder} implies that, for any $s\in{\mathcal C}_l^{\alpha}({\Omega})$ with $\|s\|_{\alpha}^{(2-\sigma)}<\infty$, $f=\mathscr{L}^{-1}_{{\tau}_0} s$ satisfies 
\begin{equation}\label{contraction}
\|\mathscr{L}^{-1}_{{\tau}_0} s\|_{2+\alpha}^{(-\sigma)}
\lesssim\|s\|_{\alpha}^{(2-\sigma)}. 
\end{equation}
For any $w\in \mathcal{Y}$, since ${\tau}_0\in\mathcal{I}$ and \eqref{contraction0} holds, 
the following Cauchy problem is solvable for any $s\in {\mathcal C}_l^\alpha({\Omega})$ with $\|s\|_{\alpha}^{(2-\sigma)}<\infty$, 
\begin{equation*}
\left\{ 
\begin{aligned}
\ & \mathscr{L}_{{\tau}_0} f=s+({\tau}-{\tau}_0)\left(\mathscr{L}_0-\mathscr{L}_1\right)w {\quad\rm in\ } {\Omega},\\
\ &\; f(0,x,v)=0 {\quad\rm in\ } \mathbb{R}^{d}\times\mathbb{R}^{d}. 
\end{aligned}
\right. 
\end{equation*} 
Thus, we can define the mapping $F:\mathcal{Y}\rightarrow \mathcal{Y}$ by setting $F(w)=f$.
Armed with \eqref{contraction} and \eqref{contraction0}, there exists a universal constant $C>0$ such that, for any $u,w\in \mathcal{Y}$, 
\begin{equation*}
\|F(u)-F(w)\|_{2+\alpha}^{(-\sigma)} 
\le C|{\tau}-{\tau}_0| \left\|\left(\mathscr{L}_0-\mathscr{L}_1\right)(u-w)\right\|_{{\alpha}}^{(2-\sigma)}
\le C|{\tau}-{\tau}_0|\|u-w\|_{2+\alpha}^{(-\sigma)}. 
\end{equation*}
Hence, $F$ is a contraction mapping, provided that $|{\tau}-{\tau}_0|\le \delta:=\frac{1}{2C}$. 
Then, $F$ gives a unique fixed point $f\in \mathcal{Y}$, that is the unique bounded solution to the Cauchy problem \eqref{zerocauchy} in $\mathcal{Y}$. 
By dividing the interval $[0,1]$ into subintervals of length less than $\delta$, we conclude that $1\in\mathcal{I}$. 




\medskip\noindent{\textbf{Step 2.}}
If $f_{\rm in}=0$ and the condition~\eqref{schaudercondition'} holds, we approximate the coefficient $B$ by $B_n:=B\varrho_n$, 
where $\varrho_n(v):=\varrho_1\big(\frac{v}{n}\big)$ for $v\in\mathbb{R}^d$, $n\in\mathbb{N}_+$, and $\varrho_1\in{\mathcal C}_c^\infty(B_2)$ is valued in $[0,1]$ such that $\varrho_1\equiv1$ in $B_1$. 
Then, for each $n\in\mathbb{N}_+$, the result obtained in the previous step provides a bounded solution $f_n$ to the equation~\eqref{linearequation} with $B$ replaced by $B_n$. 
Indeed, applying the maximum principle (Lemma~\ref{max}) to the function $\pm f-e^t\sup_{\Omega} |s|$ implies that $\sup_{\Omega} |f_n|\le e^T\sup_{\Omega} |s|$. 
Thanks to the interior Schauder estimate (Proposition~\ref{interiorschauder}), for any compact subset $K\subset{\Omega}$, $\{f_n\}_{n\ge N}$ is precompact in ${\mathcal C}_{\rm kin}^2(K)$, provided that $N$ (depending on $K$) is large enough. 
Sending $n\rightarrow\infty$ in the equation satisfied by $f_n$ yields that the limit function $f\in{\mathcal C}_l^{2+{\alpha}}({\Omega})$ is a bounded solution to \eqref{linearequation}, which satisfies $\sup_{\Omega} |f|\le e^T\sup_{\Omega} |s|$. 




\medskip\noindent{\textbf{Step 3.}}		
For general $f_{\rm in} \in {\mathcal C}^0(\mathbb{R}^d\times\mathbb{R}^d)$, we approximate $f_{\rm in}$ uniformly as $\varepsilon\rightarrow0$ by a sequence of smooth functions $\{f_{\rm in}^\varepsilon\}$ on $\mathbb{R}^d\times\mathbb{R}^d$. 
Thus, the function $f-f_{\rm in}^\varepsilon$ is a solution to \eqref{linearequation}, with the source term equal to $s-\mathscr{L}_1f_{\rm in}^\varepsilon$, and associated with the vanishing initial data. 
The procedure presented in the previous steps ensures a unique bounded solution $f^\varepsilon$ to \eqref{linearequation} for each $f_{\rm in}^\varepsilon$. 

The uniform convergence of $\{f_{\rm in}^\varepsilon\}$ and the maximum principle (Lemma~\ref{max}) implies the uniform convergence of $\{f^\varepsilon\}$. 
We may denote its limit by $f\in{\mathcal C}^0(\overline{{\Omega}})$, which satisfies $f(0,x,v)=f_{\rm in}(x,v)$ on $\mathbb{R}^d\times\mathbb{R}^d$. 
The interior Schauder estimate again implies that $\{f^\varepsilon\}$ is precompact in ${\mathcal C}_{\rm kin}^2(K)$ for any compact subset $K\subset{\Omega}$; 
then sending $\varepsilon\rightarrow0$ gives the solution $f\in{\mathcal C}_l^{2+{\alpha}}({\Omega})$ to the equation~\eqref{linearequation}. 
Its uniqueness is again given by the maximum principle. 
This concludes the proof. 
\end{proof}
	
	
	
	
	
\section{Well-posedness of the nonlinear model}\label{nonlinear}
This section is devoted to the proof of Theorem \ref{wellpose} (i), including a self-generating lower bound given in Subsection~\ref{lowersection},  the existence and uniqueness given in Subsection~\ref{wellsection} and smoothness a priori estimate given in Subsection~\ref{smoothsection}.

First of all, the Cauchy problem \eqref{nk} can be recasted in terms of the unknown function $g(t,x,v):= {\mu}(v)^{-\frac{1}{2}}f(t,x,v)$ with $g_{\rm in}(x,v):={\mu}(v)^{-\frac{1}{2}}f_{\rm in}(x,v)$ as follows, 
\begin{align}\label{gnk}
&\left\{ 
\begin{aligned}
\ &\left(\partial_t+ v\cdot\nabla_x\right)g=\mathcal{R}[g]\, \mathcal{U}[g], \\
\ &\; g(0,x,v)= g_{\rm in}(x,v), \\
\end{aligned}
\right. 
\end{align} 
where $\mathcal{R}[g]$ and $\mathcal{U}[g]$ on the right hand side are defined by  
\begin{equation*}
\mathcal{R}[g]:= \left(\int_{\mathbb{R}^{d}} g{\mu}^{\frac{1}{2}} \mathop{}\!\mathrm{d} v\right)^{\beta} {\quad\rm and\quad}
\mathcal{U}[g]:={\mu}^{-\frac{1}{2}}\nabla_v\left({\mu}\nabla_{v}\left({\mu}^{-\frac{1}{2}}g\right)\right)
=\Delta_vg+\left(\frac{d}{2}-\frac{|v|^{2}}{4}\right)g. 
\end{equation*}
The main advantage of this formulation is that it allows us to get rid of the first order term in $v$ and the zero order term is bounded, since $g$ is bounded from above by a Maxwellian. 

For convenience, we are also concerned with the substitution $h(t,x,v):= {\mu}(v)^{-1}f(t,x,v)$ and the Ornstein-Uhlenbeck operator $\mathscr{L}_{\scriptscriptstyle{\rm OU}}:=\left(\nabla_v-v\right)\cdot\nabla_v$. The equation~\eqref{nk} is then equivalent to 
\begin{equation}\label{rnk}
\left(\partial_t+v\cdot\nabla_x\right) h(t,x,v) = \mathscr{R}_h(t,x)\,\mathscr{L}_{\scriptscriptstyle{\rm OU}} h(t,x,v),\quad 
\mathscr{R}_h(t,x):= \left(\int_{\mathbb{R}^{d}} h(t,x,v)\mathop{}\!\mathrm{d}{\mu}\right)^{\beta}. 
\end{equation} 
In contrast with \eqref{nk}, the zero order term disappears. 
Let us begin by exhibiting the global bounds of solutions to \eqref{rnk} in $(0,T)\times\mathbb{T}^d\times\mathbb{R}^d$, which is a variant of \cite[Lemma 4.1]{IM}. 



\begin{lemma}[Global bounds]\label{Gaussian}
Let $a(t,x)\in L^\infty((0,T)\times\mathbb{T}^d)$ be nonnegative. 
Assume that $h(t,x,v)\in L^\infty \big((0,T);H^1\big(\mathbb{T}^d\times\mathbb{R}^d,\mathop{}\!\mathrm{d} m\big)\big)$ verifies the equation $\left(\partial_t+v\cdot\nabla_x\right)h=a\mathscr{L}_{\scriptscriptstyle{\rm OU}} h$ in $(0,T)\times\mathbb{T}^d\times\mathbb{R}^d$ in the sense of distributions. 
If $h(0,\cdot,\cdot)\le\Lambda$ in $\mathbb{T}^d\times\mathbb{R}^d$, then $h\le\Lambda$ in $(0,T)\times\mathbb{T}^d\times\mathbb{R}^d$; 
if $h(0,\cdot,\cdot)\ge\lambda$ in $\mathbb{T}^d\times\mathbb{R}^d$, then $h\ge\lambda$ in $(0,T)\times\mathbb{T}^d\times\mathbb{R}^d$.
\end{lemma} 
\begin{proof}
Integrating the equation $\left(\partial_t+v\cdot\nabla_x\right)(h-\Lambda)=a\mathscr{L}_{\scriptscriptstyle{\rm OU}}(h-\Lambda)$ against $(h-\Lambda)_+$ yields
\begin{align*}
\frac{1}{2}\int_{\mathbb{T}^d\times\mathbb{R}^d}\left[(h(t,\cdot,\cdot)-\Lambda)_+^2-(h(0,\cdot,\cdot)-\Lambda)_+^2\right]\mathop{}\!\mathrm{d} m&\\
=-\int_{[0,t]\times\mathbb{T}^d\times\mathbb{R}^d}a\left|\nabla_v(h-\Lambda)_+\right|^2 \mathop{}\!\mathrm{d} t \mathop{}\!\mathrm{d} m\le0&, 
\end{align*}
for any $t\in[0,T]$. 
This means that the upper bound is preserved along time. 
Similarly, the lower bound can be obtained by integrating the equation $\left(\partial_t+v\cdot\nabla_x\right)(\lambda-h)=a\mathscr{L}_{\scriptscriptstyle{\rm OU}}(\lambda-h)$ against $(\lambda-h)_+$. 
\end{proof}

In particular, the above global bounds preserving result holds for solutions to the equations~\eqref{rnk} and \eqref{epnk} in $(0,T)\times\mathbb{T}^d\times\mathbb{R}^d$. 
We will also apply such result to the substitution $g={\mu}^\frac{1}{2}h$ appearing in the following Subsection~\ref{wellsection} below. 
Unless otherwise specified, throughout this section we set the domain  ${\Omega}:=(0,T]\times\mathbb{T}^d\times\mathbb{R}^d$ with $T\in\mathbb{R}_+$. 
Nevertheless, as it is specified in Remark~\ref{lowerwholeremark}, Corollary~\ref{existwhole} and Proposition~\ref{unique} below, 
the results of this section also hold if the spatial domain is $\mathbb{R}^d$. 


\subsection{Self-generating lower bound}\label{lowersection}
Throughout this subsection, we assume that the bounded solution $h$ of \eqref{nk} lies below the universal constant $\Lambda$, which is guaranteed by Lemma~\ref{Gaussian} if the initial data lies below $\Lambda$. 
The aim of this subsection is to show the following positivity-spreading result. We remark that this proposition only relies on the mixing structure of the classical parabolic-type maximum principle and the transport operator, but not on the structure of the mass conservation. 

\begin{proposition}[Lower bound]\label{lower}
Let $\delta>0$, $\underline{T}\in(0,T)$ and $h$ be a bounded solution to \eqref{rnk} in ${\Omega}$ satisfying
\begin{equation}\label{lowers00}
h(0,x,v)\ge\delta\mathbbm{1}_{\left\{|x-x_0|<r,\,|v-v_0|<r\right\}},
\end{equation}
for some $(x_0,v_0)\in\mathbb{T}^d\times\mathbb{R}^d$. 
Then, there exist two positive continuous functions $\eta_1(t),\eta_2(t)$ on $(0,T]$ depending only on universal constants, $T,\delta,r$ and $v_0$ such that for any $(t,x,v)\in{\Omega}$, 
\begin{equation}\label{lowers0}
h(t,x,v)\ge \eta_1(t) e^{-\eta_2(t)|v|^2}. 
\end{equation}
\end{proposition}
\begin{remark}
In particular, the functions $\eta_1(t),\eta_2(t)$ are positive and bounded on any compact subset of $(0,T]$, but $\eta_1$ might degenerate to zero and $\eta_2$ may go to infinity as $t$ tends to zero or infinity.
\end{remark}

\begin{remark}\label{lowerwholeremark}
If one is concerned with the problem in the whole space, that is ${\Omega}=(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$, we can proceed along the same lines of the proof in Appendix~\ref{lowerapp} to see that \eqref{lowers00} implies the following lower bound 
\begin{equation}\label{lowerwhole}
h(t,x,v)\ge \eta_1(t,x)^{-1} e^{-\eta_2(t,x)|v|^4}, 
\end{equation}
where the functions $\eta_1(t,x),\eta_2(t,x)$ on $(0,T]\times\mathbb{R}^d$ are positive, continuous and only depend on universal constants, $T,\delta,r$ and $v_0$. 
Compared with \eqref{lowers0}, $\eta_1(t,x),\eta_2(t,x)$ lose the uniformity in $x$ as $\mathbb{R}^d$ is not compact (see Step~3 of the proof of the proposition in Appendix~\ref{lowerapp}). 
Besides, the exponential tail with respect to $v$ cannot be improved to a Gaussian type, since there is no uniform-in-$x$ lower bound on the local mass $\int h\mathop{}\!\mathrm{d}{\mu}$ so that Step~4 in Appendix~\ref{lowerapp} fails. 
\end{remark}
	
	
Note that the proof of the proposition is composed mainly of two lemmas. 
On the one hand, Lemma~\ref{lowerlemma1} below extends the lower bounds forward a short time from a neighborhood of any given point in $\mathbb{T}^d\times\mathbb{R}^d$ and at any given time. 
On the other hand, Lemma~\ref{lowerlemma2} below is used to spread the lower bound for all velocities. 
The spreading of the lower bound in space is given by selecting proper velocity to transport the positivity which is guaranteed by Lemma \ref{lowerlemma1}. 
By applying these two lemmas repeatedly, as proposed in \cite{HST3}, we are able to control the solution from below for any finite time. 
We postpone the full proof of Proposition~\ref{lower}, 
obtained by combining the previously introduced lemmas, 
until Appendix~\ref{lowerapp}. 


	
\begin{lemma}[Lower bound forward in time]\label{lowerlemma1}
Let $\delta,\tau,r\in(0,1]$ and $h$ be a bounded solution to \eqref{rnk} in ${\Omega}$ with 
\begin{equation*}
h(0,x,v)\ge\delta\mathbbm{1}_{\left\{|x-x_0|<r,\, |v-v_0|<\frac{r}{\tau}\right\}},
\end{equation*}
for some $(x_0,v_0)\in\mathbb{T}^d\times\mathbb{R}^d$. Then, there exists some universal constant $c_0>0$ such that 
\begin{equation*}
h\ge\frac{\delta}{8}\mathbbm{1}_\mathcal{P},{\quad}
\mathcal{P}:=\left\{t\le \min\big\{T,\tau,c_0\left\langle\tau r^{-1}\right\rangle^{-2}\!\langle v_0\rangle^{-2}\big\},\, |x-x_0-tv|<\frac{r}{2},\, |v-v_0|<\frac{r}{2\tau}\right\}.
\end{equation*}
\end{lemma}
	
\begin{proof}
Let us consider the barrier function $$\underline{h}(t,x,v):=-C_0t+\frac{\delta}{2}\left(1-r^{-2}|x-x_0-tv|^2-\tau^2r^{-2}|v-v_0|^2\right),$$ 
where the constant $C_0>0$ is to be determined. 
The region $\mathcal{Q}:=\big\{t\le\min\{T,\tau\},\,|x-x_0-tv|^2+\tau^2|v-v_0|^2<r^2\big\}$ contains $\mathcal{P}$. 
A direct computation yields that
\begin{equation*}
|\mathscr{L}_{\scriptscriptstyle{\rm OU}}\underline{h}|\le |\Delta_v\underline{h}|+ |v\cdot\nabla_v\underline{h}|
\lesssim \delta\left\langle\tau r^{-1}\right\rangle^{2}\!\langle v_0\rangle^2 {\quad\rm in\ }\mathcal{Q}. 
\end{equation*}
By choosing $C_0:=\frac{1}{8c_0}\delta\left\langle\tau r^{-1}\right\rangle^{2}\!\langle v_0\rangle^2$ for some (small) universal constant $c_0>0$, we have
\begin{equation}\label{hunder}
\left(\partial_t+v\cdot\nabla_x-\mathscr{R}_h\,\mathscr{L}_{\scriptscriptstyle{\rm OU}}\right)\underline{h}\le -C_0+\Lambda^{\beta}|\mathscr{L}_{\scriptscriptstyle{\rm OU}}\underline{h}|<0 
{\quad\rm in\ }\mathcal{Q}.  
\end{equation}
Besides, $\underline{h}(t,x,v)\ge\frac{\delta}{8}$ in $\left\{t\le c_0\left\langle\tau r^{-1}\right\rangle^{-2}\!\langle v_0\rangle^{-2},\; |x-x_0-tv|^2+\tau^2|v-v_0|^2<\frac{r^2}{2}\right\}\supset\mathcal{P}$. 

Then, observing that $\underline{h}-h\le0$ on the parabolic boundary $\{t=0\}\cap\mathcal{Q}$,  $\big\{t\le\min\{T,\tau\},\,|x-x_0-tv|^2+\tau^2|v-v_0|^2=r^2\big\}$, and applying the classical maximum principle to $\underline{h}-h$ in $\mathcal{Q}$ yields the result. 
\end{proof}


The spreading of lower bound to all velocities relies on the construction of a Harnack chain through iterative application of the local Harnack inequality \cite[Theorem 1.6]{GIMV}. 
Although some coefficients of the equation~\eqref{rnk} are unbounded globally over $v\in\mathbb{R}^d$, we remark that their local boundedness is sufficient for us to achieve the result through a careful study on the rescaling during the construction of the Harnack chain. 

\begin{lemma}[Lower bound for all velocities]\label{lowerlemma2}
Let $\delta>0$, $T,R\in(0,1]$, $\underline{T}\in(0,T)$ and $h$ be a bounded solution to \eqref{rnk} in ${\Omega}$ such that, for any $t\in[0,T]$,  
\begin{equation}\label{lower2assumption}
h(t,x,v)\ge\delta\mathbbm{1}_{\left\{|x-x_0-tv_0|<R,\,|v-v_0|<R\right\}},
\end{equation}
for some $(x_0,v_0)\in\mathbb{T}^d\times\mathbb{R}^d$. Then, there exists some (large) constant $\underline{C}>0$ depending only on universal constants, $\underline{T},\delta,R$ and $v_0$ such that, for any $t\in[\underline{T},T]$, we have 
\begin{equation}\label{lower2bound}
h(t,x,v)\ge \underline{C}^{-1} e^{-\underline{C}|v|^4}\mathbbm{1}_{\left\{|x-x_0-tv_0|<\frac{R}{2}\right\}}. 
\end{equation}
\end{lemma}
	
\begin{proof}
For any $z:=(t,x,v)\in\big\{t\in[\underline{T},T],\,|x-x_0-tv_0|<\frac{R}{2},\,v\in\mathbb{R}^d\big\}$, we will construct a finite sequence of points to reach $z$ from the region $\big\{t\le{T},\,|x-x_0-tv_0|<R,\,|v-v_0|<R\big\}$ where the solution is positive by the assumption. 
In particular, $x$ does not exit this region. 
The nonlocality of the coefficient $\mathscr{R}_h$, with the assumption \eqref{lower2assumption}, implies the nondegeneracy of the diffusion in velocity so that the positivity of the solution $h$ propagates over $v\in\mathbb{R}^d$ in a localized space region. 

\medskip\noindent{\textbf{Step 1.}}	Iterate the Harnack inequality. \\
For $i\in\{1,2,...,N+1\}$ with $N\in\mathbb{N}$, we define $z_{N+1}:=z$ and $z_{i}:=(t_i,x_i,v_i)$ by the relation 
$$z_i=z_{i+1}\circ S_r\left(-{\tau}_1,0,-{\tau}_2\frac{v-v_0}{|v-v_0|}\right),$$  
where the parameters $N,r,\tau_1,\tau_2>0$ are to be determined next. 	
Consider the function for $\tilde{z}:=(\tilde{t},\tilde{x},\tilde{v})\in Q_1$, 
\begin{equation*}
h_i(\tilde{z}):=h\left(z_i\circ S_r(\tilde{z})\right)
=h(t_i+r^2\tilde{t},\,x_i+r^3\tilde{x}+r^2\tilde{t}v_i,\,v_i+r\tilde{v}).
\end{equation*} 
We observe that if the following is true for any $\tilde{z}\in Q_1$ 
\begin{align}
&\;t_{i+1}+r^2\tilde{t}\in[0,T],\quad Nr\tau_2\le |v-v_0|, \label{htxv21}\\
&\left|x_{i+1}+r^3\tilde{x}+r^2\tilde{t}v_{i+1}-x_0-\left(t_{i+1}+r^2\tilde{t}\right)v_0\right|<R,\label{htxv22}
\end{align}
then the function $h_{i+1}(\tilde{z})$, for $1\le i\le N$, verifies the equation 
\begin{equation*}
\left(\partial_{\tilde{t}}+\tilde{v}\cdot\nabla_{\tilde{x}}\right)h_i
=\mathscr{R}_h\left(z_i\circ S_r(\tilde{z})\right)\, \left(\Delta_{\tilde{v}}h_i-r(v_i+r\tilde{v})\cdot\nabla_{\tilde{v}}h_i\right) {\rm\quad in\ } Q_1, 
\end{equation*}
where the coefficients satisfy 
		\begin{equation*}
		\delta^{\beta} R^{d{\beta}}\lesssim\mathscr{R}_h\lesssim 1  {\quad\rm and\quad}
		|r(v_i+r\tilde{v})|\le r(1+|v_0|+|v-v_0|)\le 1, 
		\end{equation*}
		provided that $r\le(1+|v_0|+|v-v_0|)^{-1}$. 
Applying the Harnack inequality \cite[Theorem 1.6]{GIMV} to $h_i$ implies that there exist constants $c_0,\tau_1\in(0,1)$ depending only on universal constants, $\delta$ and $R$ such that, for any $\tau_2\in[0,1-{\tau}_1]$ and $1\le i\le N$, we have
\begin{equation}\label{hharnack}
h(z_{i+1})=h_{i+1}(0,0,0)\ge c_0h_{i+1}\Big(-\tau_1,\,0,\,-\tau_2\frac{v-v_0}{|v-v_0|}\Big)=c_0h(z_i).  
\end{equation}
Hence, it remains to determine the chain $\{z_i\}_{1\le i\le N+1}$ such that the conditions \eqref{htxv21}, \eqref{htxv22} hold and the point $z_1$ stays in the region $\big\{(t,x,v):t\le{T},\,|x-x_0-tv_0|<R,\,|v-v_0|<R\big\}$. 

\medskip\noindent{\textbf{Step 2.}}	Determine the Harnack chain (including $N,r,\tau_2$) from a proper starting time $t_1$.\\ 
For $M>0$, we set 
\begin{equation*}
t_1:=\max\left\{\frac{\underline{T}}{2},\,t-\frac{R}{8}(1+|v_0|+|v-v_0|)^{-1}\right\}
{\quad\rm and\quad}
r:=\frac{R}{M}(1+|v_0|+|v-v_0|)^{-2}.
\end{equation*}
Recalling that $T,R\in(0,1]$, by choosing 
$M\ge \frac{2}{\underline{T}} +\frac{{\tau}_1}{1-{\tau}_1}\left(8+\frac{2}{\underline{T}}\right)$, we have 
\begin{equation*}
r^2\le\frac{\underline{T}}{2} {\quad \rm and\quad}
\tau_2:=\frac{r\tau_1|v-v_0|}{t-t_1}\le 1-{\tau}_1. 
\end{equation*}
To determine the parameter $M>0$, 
we point out that there exists some constant $\overline{C}$ depending only on universal constants, $\underline{T},\delta,R$ and $v_0$ such that $M\le\overline{C}$ and 
\begin{equation*}
N:=\frac{t-t_1}{r^2\tau_1}\in\mathbb{N}^+. 
\end{equation*}
Thus, $Nr\tau_2=|v-v_0|$. This setting then guarantees the condition \eqref{htxv21}. 

It also follows from the iteration relation that $v_1=v_0$, and for $1\le i\le N+1$,  
\begin{equation}\label{iterationrelation}
t_i=t_1+(i-1)r^2{\tau}_1,\quad 
v_{i}=v_0+(i-1)r{\tau}_2\frac{v-v_0}{|v-v_0|},\quad 
x_{i}=x-r^2\tau_1\sum\nolimits^{N}_{j=i}v_{j+1}. 
\end{equation} 


\medskip\noindent{\textbf{Step 3.}} Determine the starting point $x_1$.\\	
For any $1\le i\le N$, we estimate the departure distance from the expression \eqref{iterationrelation}, 
\begin{equation*}
\left|x_{i+1}-x_1-(t_{i+1}-t_1)v_0\right| =\frac{i(i+1)}{2} r^3\tau_1\tau_2
\le N^2 r^3\tau_1\tau_2 =(t-t_1)|v-v_0|\le\frac{R}{8}. 
\end{equation*}
Therefore, for any $x\in B_\frac{R}{2}(x_0+tv_0)$, there exists some $x_1\in B_{\frac{5R}{8}}(x_0+t_1v_0)$ such that $x_{N+1}=x$. 
In this setting, for any $1\le i\le N$, we also have 
\begin{align*}
&\left|x_{i+1}+r^3\tilde{x}+r^2\tilde{t}v_{i+1}-x_0-\left(t_{i+1}+r^2\tilde{t}\right)v_0\right|\\
&\le \left|x_{i+1}-x_1-(t_{i+1}-t_1)v_0\right|+|x_1-x_0-t_1v_0| +r^2|r\tilde{x}+\tilde{t}v_{i+1}-\tilde{t}v_0|\\
&\le \frac{R}{8}+\frac{5R}{8}+r^2(1+|v-v_0|)<\frac{3R}{4}+\frac{R^2}{M^2}<R. 
\end{align*}
Thus, the condition~\eqref{htxv22} ensuring the inequality~\eqref{hharnack} is satisfied for $1\le i\le N$, which yields 
\begin{equation*}
h({t},{x},{v})\ge c_0^Nh(t_1,x_1,v_0) \ge {\delta} e^{-N\log \frac{1}{c_0}}. 
\end{equation*}
Recalling that $c_0\in(0,1)$ appears in \eqref{hharnack} and  $N\le\frac{T\overline{C}^2\left(1+|v_0|+|v-v_0|\right)^4}{{\tau}_1R^2}$, we obtain the desired result \eqref{lower2bound}. 
\end{proof}



	
	
\subsection{Existence and uniqueness}\label{wellsection}
Let us begin by summarizing some basic a priori estimates for solutions to the  equation~\eqref{gnk}. 
\begin{lemma}[H\"older estimates]\label{boundschauder} 
Let ${\Omega}_x=\mathbb{T}^d$ or $\mathbb{R}^d$, and $g$ be a solution to \eqref{gnk} in $[0,T]\times{\Omega}_x\times\mathbb{R}^d$ satisfying 
\begin{equation*}
\mathcal{R}[g]\ge\lambda  {\quad\rm in\ } [0,T]\times{\Omega}_x {\quad\rm and\quad}
0\le g_{\rm in}\le\Lambda{\mu}^{\frac{1}{2}} {\quad\rm in\ }{\Omega}_x\times\mathbb{R}^d. 
\end{equation*}
\begin{item}
\item [ {\bf\ \,(i)}] 
Let $\underline{T}\in(0,T)$ and ${\delta}\in\big(0,\frac{1}{2}\big)$. 
There exists some universal constant ${\alpha}\in(0,1)$ such that, 
for any $Q_{2r}(z_{0})\subset[\underline{T},T]\times{\Omega}_x\times\mathbb{R}^d$, we have 
\begin{equation}\label{schauder1}
\| g \|_{{\mathcal C}_l^{2+{\alpha}}(Q_{r}(z_{0}))} \lesssim_{\underline{T},{\delta}} {\mu}^{{\delta}} (v_{0}).
\end{equation}

\item [ {\bf\;(ii)}] 
If $g_{\rm in}\in {\mathcal C}^{{\alpha}_0}({\Omega}_x\times\mathbb{R}^d)$ with (universal) ${\alpha}_0\in(0,1)$, 
then for any ${\delta}\in\big(0,\frac{1}{2}\big)$, there exists some universal constant $\alpha\in(0,1)$ such that 
\begin{equation*}
\|g\|_{{\mathcal C}_l^{\alpha}\left([0,T]\times{\Omega}_x\times B_1(v_0)\right)}
\lesssim_{\delta} \big(1+[g_{\rm in}]_{{\mathcal C}^{{\alpha}_0}({\Omega}_x\times\mathbb{R}^d)}\big) {\mu}^{\delta}(v_0).
\end{equation*}
\end{item}
\end{lemma}

We remark that, armed with Lemma~\ref{Gaussian}, the assertions (i) and (ii) in the above lemma directly follow from \cite[Proposition 4.4]{IM} and \cite[Corollary 4.6]{YZ}, respectively. 
	
	
\begin{proposition}[Existence]\label{existence}
For any $g_{\rm in}\in{\mathcal C}^0(\mathbb{T}^d\times\mathbb{R}^d)$ such that $0\le g_{\rm in}\le\Lambda{\mu}^\frac{1}{2}$ in $\mathbb{T}^d\times\mathbb{R}^d$, there exists a (classical) solution $g$ to the Cauchy problem~\eqref{gnk} satisfying $0\le g\le\Lambda{\mu}^\frac{1}{2}$ in ${\Omega}$. 
\end{proposition}

\begin{remark}\label{remark-weak}
For any given nonnegative continuous function $g_{\rm in}$ that is not identically zero, there is some point $(x_0,v_0)\in\mathbb{T}^d\times\mathbb{R}^d$ and some constants $\delta,r>0$ such that 
\begin{equation}\label{lowerdd}
g_{\rm in}(x,v)\ge\delta\mathbbm{1}_{\{|x-x_0|<r,|v-v_0|<r\}} 
{\quad\rm in\ }\mathbb{T}^d\times\mathbb{R}^d. 
\end{equation} 
We will see that the upper bound $g_{\rm in}\le\Lambda{\mu}^\frac{1}{2}$ and the lower bound \eqref{lowerdd} assumptions on the initial data $g_{\rm in}$ (could be discontinuous) are sufficient to ensure the existence of solution $g\in{\mathcal C}_{\rm kin}^2({\Omega})$ in the weak sense that for any $\phi\in {\mathcal C}^\infty_c([0,T)\times\mathbb{T}^d\times\mathbb{R}^d)$, 
\begin{equation}\label{existence0}
	\int_{\mathbb{T}^d\times\mathbb{R}^d} g_{\rm in}\,\phi\big|_{t=0}
	=\int_{{\Omega}}\left\{-g\left(\partial_t+v\cdot\nabla_x\right)\phi
	+\mathcal{R}[g]\nabla_vg\cdot\nabla_v\phi-\mathcal{R}[g]\left(\frac{d}{2}-\frac{|v|^2}{4}\right)g\phi\right\}.
\end{equation}
As solutions become regular instantaneously, the difference between the weak solution and the classical one lies only in the continuity around the initial time. 
\end{remark}


\begin{proof}
Let us assume that $g_{\rm in}$ satisfies \eqref{lowerdd}, for some point $(x_0,v_0)\in\mathbb{T}^d\times\mathbb{R}^d$ and some constants $\delta,r>0$. 
By Proposition~\ref{lower}, for any solution $g$ to \eqref{gnk} and for any $\underline{T}\in(0,T)$, there is some $\lambda_*>0$ depending only on universal constants, $\underline{T},T,\delta,r$ and $v_0$ such that, 
\begin{equation}\label{exist-lower}
\mathcal{R}[g](t,x)\ge \lambda_*
{\quad\rm in\ }[\underline{T},T]\times\mathbb{T}^d. 
\end{equation}	

\medskip\noindent{\textbf{Step 1.}}	
We first approximate the initial data $g_{\rm in}$ by 
$g_{\rm in}^{\epsilon}:=g_{\rm in}*\varrho_\varepsilon+\varepsilon{\mu}^{\frac{1}{2}}$, where
$\varrho_\varepsilon(x,v):=\frac{1}{\varepsilon^{2d}}\varrho_1\left(\frac{x}{\varepsilon},\frac{v}{\varepsilon}\right)$ 
with $(x,v)\in\mathbb{T}^d\times\mathbb{R}^d$, $\varepsilon\in(0,1]$, and $\varrho_1\in {\mathcal C}^\infty_c\left(B_1\times B_1\right)$ is a nonnegative bump function such that $\int_{\mathbb{R}^{2d}}\varrho_1=1$. 
Then, we have $\varepsilon\mu^{\frac{1}{2}}\le g_{\rm in}^\varepsilon\le(1+\Lambda)\mu^{\frac{1}{2}}$ in $\mathbb{T}^d\times\mathbb{R}^d$.

Let us fix $\varepsilon\in(0,1]$. 
In order to establish the existence of solution to \eqref{gnk} associated with the initial data $g_{\rm in}^\varepsilon$, 
we are going to find a fixed point of the mapping $F:w\mapsto{g}$ defined by solving the Cauchy problem, 
\begin{align}\label{nonexun}
&\left\{ 
\begin{aligned}
\ &\left(\partial_t+v\cdot\nabla_x\right)g
= \mathcal{R}[w]\,\mathcal{U}[g]
{\ \ \rm in\ } {\Omega}, \\
\ &\; {g}(0,\cdot,\cdot)=g_{\rm in}^\varepsilon {\ \ \rm in\ }\mathbb{T}^d\times\mathbb{R}^d,\\
\end{aligned}
\right. 
\end{align} 
on the closed convex subset $\mathcal{K}$ of the Banach space ${\mathcal C}_l^{\gamma}(\overline{{\Omega}})$, 
\begin{equation*}
\mathcal{K}:=\Big\{
w\in{\mathcal C}_l^{{\gamma}}(\overline{{\Omega}})
:\,\|w\|_{{\mathcal C}_l^{{\gamma}}({\Omega})}\le\mathcal{N},\ \;
\varepsilon\mu^{\frac{1}{2}}\le w\le(1+\Lambda)\mu^{\frac{1}{2}}{\ \ \rm in\ }\overline{{\Omega}}
\Big\}, 
\end{equation*}	
where the constants ${{\gamma}}\in(0,1)$ and $\mathcal{N}>0$ are to be determined. 
We remark that the equation~\eqref{nonexun} is equivalent to 
\begin{equation*}
\left(\partial_t+v\cdot\nabla_x\right)\big({\mu}^{-\frac{1}{2}}g\big)
= \mathcal{R}[w]\,\mathscr{L}_{\scriptscriptstyle{\rm OU}}\big({\mu}^{-\frac{1}{2}}g\big). 
\end{equation*}
By Lemma~\ref{Gaussian} and the fact that $\mathcal{R}\left[w\right]\ge \varepsilon$, we have  $\varepsilon\mu^{\frac{1}{2}}\le{g}\le(1+\Lambda)\mu^{\frac{1}{2}}$ in $\overline{{\Omega}}$. 
In particular, the lower order term $\big|\mathcal{R}\left[w\right]\big(\frac{d}{2}-\frac{|v|^{2}}{4}\big)g\big|\lesssim1$ for any $w\in\mathcal{K}$. 
Thus, the global H\"older estimate \cite[Corollary 4.6]{YZ} implies that
there exist some constants ${\gamma}\in(0,1)$ and $\mathcal{N}>0$ depending only on universal constants and $\varepsilon$ such that $\|{g}\|_{{\mathcal C}_l^{2{\gamma}}({\Omega})}\le\mathcal{N}$, which also implies that the lower order term $\big|\mathcal{R}\left[w\right]\big(\frac{d}{2}-\frac{|v|^{2}}{4}\big)g\big|$ is bounded in ${\mathcal C}_l^{2{\gamma}}({\Omega})$. 
It then follows from Proposition~\ref{variable} with the interior Schauder estimate (Proposition~\ref{interiorschauder}) that the mapping $F:\mathcal{K}\rightarrow\mathcal{K}\cap{\mathcal C}_l^{2{\gamma}}(\overline{{\Omega}})\cap{\mathcal C}_l^{2+2{\gamma}}({\Omega})$ is well-defined. 
Besides, with the help of the Arzelà–Ascoli theorem, we know that $F(\mathcal{K})$ is precompact in ${\mathcal C}_l^{\gamma}(\overline{{\Omega}})$. 

As far as the continuity of $F$ is  concerned, we take a sequence $\{w_n\}$ converging to $w_\infty$ in ${\mathcal C}_l^{\gamma}(\overline{{\Omega}})$. 
Since $\{F(w_n)\}$ is precompact in ${\mathcal C}_l^{\gamma}(\overline{{\Omega}})$, there exists a converging subsequence whose limit is ${g}_\infty\in{\mathcal C}_l^{\gamma}(\overline{{\Omega}})$ which satisfies 
${g}_\infty(0,\cdot,\cdot)=g_{\rm in}^\varepsilon$ in $\mathbb{T}^d\times\mathbb{R}^d$. 
In view of the interior Schauder estimate (Proposition~\ref{interiorschauder}), $\{F(w_n)\}$ is precompact in ${\mathcal C}_{\rm kin}^2(K)$ for any compact subset $K\subset{\Omega}$ and ${g}_\infty\in{\mathcal C}_{\rm kin}^2({\Omega})\cap\mathcal{C}^0(\overline{{\Omega}})$. 
Sending $n\rightarrow\infty$ in \eqref{nonexun} satisfied by $(w,{g})=(w_n,F(w_n))$, we see that the equation \eqref{nonexun} also holds for the couple of limits $(w,{g})=(w_\infty,{g}_\infty)$. 
Then, applying the maximum principle (Lemma~\ref{max}) to 
\begin{align*}
&\left\{ 
\begin{aligned}
\ &\left(\partial_t+v\cdot\nabla_x\right)\left({\mu}^{-\frac{1}{2}}(g_\infty-F(w_\infty))\right)
=\mathcal{R}[w_\infty]\,\mathscr{L}_{\scriptscriptstyle{\rm OU}} \left({\mu}^{-\frac{1}{2}}\left(g_\infty-F(w_\infty)\right)\right)
{\quad\rm in\ }{\Omega}, \\
\ &\; (g_\infty-F(w_\infty))(0,\cdot,\cdot)=0 {\ \ \rm in\ }\mathbb{T}^d\times\mathbb{R}^d,\\
\end{aligned}
\right. 
\end{align*} 
we arrive at ${g}_\infty=F(w_\infty)$. 

Then for every $\varepsilon\in(0,1]$, we are allowed to apply the Schauder fixed point theorem (see for instance \cite[Corollary 11.2]{GilbargTrudinger}) to get $g^\varepsilon\in{\mathcal C}_{\rm kin}^2({\Omega})\cap\mathcal{C}^0(\overline{{\Omega}})$ such that $F(g^\varepsilon)=g^\varepsilon$, which is a (classical) solution to \eqref{gnk} associated with the initial data $g_{\rm in}^\varepsilon$. 
		
\medskip
\noindent{\textbf{Step 2.}}	Passage to the limit.\\
Recalling the lower bound~\eqref{exist-lower} on the coefficient and the higher order H\"older estimate given by Lemma~\ref{boundschauder}~(i), we point out that for any $\underline{T}\in(0,T)$, $\{g^\varepsilon\}$ is uniformly bounded in ${\mathcal C}_l^{2+{\alpha}_*}([\underline{T},T]\times\mathbb{T}^d\times\mathbb{R}^d)$, for some constant ${\alpha}_*\in(0,1)$ with the same dependence as $\lambda_*$. 
Hence, $g^\varepsilon$ converges uniformly to $g$ in ${\mathcal C}_{\rm kin}^2([\underline{T},T]\times\mathbb{T}^d\times\mathbb{R}^d)$, up to a subsequence. 
		
Write the equation satisfied by $g^\varepsilon$ in the weak formulation, that is, for any $\phi\in {\mathcal C}^\infty_c(\overline{{\Omega}})$, 
\begin{equation}\label{epexistence}
\begin{split}
&\int_{\mathbb{T}^d\times\mathbb{R}^d}\left[g^\varepsilon(T,x,v)\phi(T,x,v)-g^\varepsilon_{\rm in}(x,v)\phi(0,x,v)\right]\\
&=\int_{{\Omega}}\left\{g^\varepsilon\left(\partial_t+v\cdot\nabla_x\right)\phi-\mathcal{R}[g^\varepsilon]\nabla_vg^\varepsilon\cdot\nabla_v\phi+\mathcal{R}[g^\varepsilon]\left(\frac{d}{2}-\frac{|v|^2}{4}\right)g^\varepsilon\phi\right\}.
\end{split}
\end{equation}
Combining the energy estimate derived by choosing $\phi=g^\varepsilon$ above with the upper bound of $g^\varepsilon$ provided by Lemma~\ref{Gaussian}, we have  
\begin{equation*}
\int_{{\Omega}}\left|\mathcal{R}\left[g^\varepsilon\right]\nabla_vg^\varepsilon\right|^2
\lesssim \int_{{\Omega}}\mathcal{R}\left[g^\varepsilon\right]|\nabla_vg^\varepsilon|^2
\le \frac{1}{2}\int_{\mathbb{T}^d\times\mathbb{R}^d}|g_{\rm in}^\varepsilon|^2
 +\int_{{\Omega}}\mathcal{R}[g^\varepsilon]\left(\frac{d}{2}-\frac{|v|^2}{4}\right)|g^\varepsilon|^2
\lesssim 1.
\end{equation*}
Therefore, after passing to a subsequence, $\mathcal{R}[g^\varepsilon]\nabla_vg^\varepsilon$ converges weakly in $L^2({\Omega})$. On account of its local uniform convergence, we know that its weak limit is $\mathcal{R}[g]\nabla_vg$. 
Besides, since ${\mu}^{-\frac{1}{2}}g^\varepsilon$ is uniformly bounded, by their local uniform convergence, we can also derive that the sequences $g^\varepsilon$ and $\mathcal{R}[g^\varepsilon]\left(\frac{d}{2}-\frac{|v|^2}{4}\right) g^\varepsilon$ converge to  $g$ and $\mathcal{R}[g]\left(\frac{d}{2}-\frac{|v|^2}{4}\right) g$, respectively, weakly in $L^2({\Omega})$, up to a subsequence. 		
Then, for any $\phi\in {\mathcal C}^\infty_c([0,T)\times\mathbb{T}^d\times\mathbb{R}^d)$, sending $\varepsilon\rightarrow0$ in \eqref{epexistence} gives \eqref{existence0}. 

Furthermore, if the initial data $g_{\rm in}$ is continuous, then the barrier function method shows that the continuity around the initial time depends only on the upper bound of the solution and the continuity of $g_{\rm in}$; see the derivation of the estimate~\eqref{coninitialbarrier} of a general type in Subsection~\ref{finitetimesection} below. 
Indeed, by \eqref{coninitialbarrier} (with ${\epsilon}=1$, $R=|v_1|$, $h_{\epsilon}={\mu}^{-\frac{1}{2}}g$ and $h_{{\epsilon},\rm in}={\mu}^{-\frac{1}{2}}g_{\rm in}$), we see that 
for any fixed ${\delta}\in(0,1)$, $(x_1,v_1)\in\mathbb{T}^d\times\mathbb{R}^d$ and for any $(t,x,v)\in \big[0,\frac{{\delta}}{4(1+|v_1|)}\big]\times B_{\delta}(x_1,v_1)$, 
\begin{equation}\label{barrierinitial}
\begin{split}
|g(t,x,v)-g_{\rm in}(x_1,v_1)|
\lesssim&\; {\delta}^{-2}\langle v_1\rangle^2{\mu}^\frac{1}{2}(v_1)t + {\delta}^{-2}{\mu}^\frac{1}{2}(v_1)\left(|x-x_1-tv|^2+|v-v_1|^2\right)\\
 &+{\mu}^\frac{1}{2}(v_1)\sup\nolimits_{B_{\delta}(x_1,v_1)}|g_{\rm in}(x,v)-g_{\rm in}(x_1,v_1)|\\
\lesssim&\; {\delta}^{-2}\left(t+|x-x_1|^2+{\mu}^\frac{1}{2}(v_1)|v-v_1|^2\right)\\
  &+\sup\nolimits_{B_{\delta}(x_1,v_1)}|g_{\rm in}(x,v)-g_{\rm in}(x_1,v_1)|. 
\end{split}
\end{equation}
It implies the continuity of the solution $g$ around $t=0$. 
This finishes the proof. 
\end{proof}




One may extend the above existence result to the case where the spacial domain $\mathbb{T}^d$ is replaced by $\mathbb{R}^d$. 
\begin{corollary}\label{existwhole}
For any $g_{\rm in}\in{\mathcal C}^0(\mathbb{R}^d\times\mathbb{R}^d)$ such that $0\le g_{\rm in}\le\Lambda{\mu}^\frac{1}{2}$ in $\mathbb{R}^d\times\mathbb{R}^d$, there exists a solution $g$ to the Cauchy problem~\eqref{gnk} satisfying $0\le g\le\Lambda{\mu}^\frac{1}{2}$ in $(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$. 
If additionally $h_{\rm in}={\mu}^{-\frac{1}{2}}g_{\rm in}$ satisfies $h_{\rm in}\ge\lambda$ and $h_{\rm in}-M_1\in L^1(\mathbb{R}^d\times\mathbb{R}^d,\mathop{}\!\mathrm{d} m)$ for some constant $M_1>0$, then $h={\mu}^{-\frac{1}{2}}g$ satisfies $h\ge\lambda$  in $(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$ and 
\begin{equation}\label{l2l1}
\|h-M_1\|_{L_t^\infty([0,T]; L^1(\mathbb{R}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m))}\le \|h_{\rm in}-M_1\|_{L^1(\mathbb{R}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m)}.
\end{equation}
\end{corollary}

\begin{proof}
For $R>1$, we set $g_{\rm in}^R=g_{\rm in}\mathbbm{1}_{[-R+R^{-1},R-R^{-1}]^d}$ for $x\in[-R,R]^d$ with periodic extension to $\mathbb{R}^d$.
In the light of Proposition~\ref{existence}, we take a solution $g^R$ to the equation~\eqref{gnk} associated with the initial data $g_{\rm in}^R$ in $(0,T]\times[-R,R]^d\times\mathbb{R}^d$, where $[-R,R]^d$ is considered as a periodic box. 
After extracting a subsequence, we define the function $g:=\lim_{R\rightarrow\infty}g^R$ in $(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$ pointwisely; 
furthermore, since $0\le {\mu}^{-\frac{1}{2}}g^R\le\Lambda$ in $(0,T]\times[-R,R]^d\times\mathbb{R}^d$, 
we know that the limiting function satisfies $0\le {\mu}^{-\frac{1}{2}}g\le\Lambda$ in $(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$. 
Similarly, ${\mu}^{-\frac{1}{2}}g_{\rm in}\ge\lambda$ in $\mathbb{R}^d\times\mathbb{R}^d$ implies that ${\mu}^{-\frac{1}{2}}g\ge\lambda$ in $(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$. 

Since the initial data is continuous, unless it is identically zero, we may assume that $g_{\rm in}\ge\delta\mathbbm{1}_{\{|x-x_0|<r,|v-v_0|<r\}}$, for some point $(x_0,v_0)\in\mathbb{R}^d\times\mathbb{R}^d$ and some constants $\delta,r>0$. 
Consider $R>|x_0|+r$. 
Applying the the lower bound of the solution given by \eqref{lowerwhole} yields that, for any compact subset $K\subset(0,T]\times\mathbb{R}^d\times\mathbb{R}^d$, the coefficient $\mathcal{R}[g^R]\ge\lambda_*$, where the constant $\lambda_*>0$ only depends on universal constants,  $\delta,r,v_0$ and $K$.  
In view of the higher order H\"older estimate given by Lemma~\ref{boundschauder}~(i), we know that $g^R$ uniformly converges to $g$ in ${\mathcal C}_{\rm kin}^2(K)$, up to a subsequence. 
Besides, due to the estimate derived in \eqref{barrierinitial}, the limiting function $g$ is a solution to~\eqref{gnk} that matches the initial data $g_{\rm in}$ continuously. 


As for \eqref{l2l1}, we notice that the function $\big(h^R-M_1\big)_\pm$ with $h^R:={\mu}^{-\frac{1}{2}}g^R$ verifies 
\begin{equation*}
\left(\partial_t+v\cdot\nabla_x\right)\big(h^R-M_1\big)_\pm \le \mathscr{R}_{h^R}\,\mathscr{L}_{\scriptscriptstyle{\rm OU}} \big(h^R-M_1\big)_\pm  {\quad \rm in\ } (0,T]\times[-R,R]^d\times\mathbb{R}^d. 
\end{equation*} 
Integrating the equation against the function $\int_{[-R,R]^d}\big(h^R-M_1\big)_\pm\mathop{}\!\mathrm{d} x$ yields that 
\begin{align*}
\int_{[-R,R]^d\times\mathbb{R}^d}\big(h^R(t,\cdot,\cdot)-M_1\big)_\pm\mathop{}\!\mathrm{d} m -\int_{[-R,R]^d\times\mathbb{R}^d}\big(h^R(0,\cdot,\cdot)-M_1\big)_\pm\mathop{}\!\mathrm{d} m\le0. 
\end{align*}
Sending $R\rightarrow\infty$, we acquire
\begin{equation*}
\|(h-M_1)_\pm\|_{L_t^\infty([0,T]; L^1(\mathbb{R}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m))}
\le \|(h_{\rm in}-M_1)_\pm\|_{L^1(\mathbb{R}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m)}, 
\end{equation*}
which implies the estimate \eqref{l2l1} as asserted. 
The proof is complete. 
\end{proof}





The following proposition concerned with the uniqueness of the Cauchy problem~\eqref{gnk} is derived from a Grönwall-type argument. 
The standard scaling technique and the H\"older estimate up to the initial time given by Lemma~\ref{boundschauder}~(ii) can improve the integrability with respect to $t$ in the energy estimate so that Grönwall's inequality becomes admissible; see \eqref{gronwall} below for the precise expression. This kind of phenomena was also noticed in \cite{HST1} (see the remarks in \S1.4.2). 
Besides, the global energy estimate of the equation~\eqref{gnk} is not available when the spatial domain is unbounded, since there is no decay of the solution as $|x|\rightarrow\infty$. 
To work it out, we take advantage of the idea originated from the uniformly local space that was used in \cite{Kato,HST2}. We also remark that such a technique is not necessary when working with the periodic box $\mathbb{T}^d$. 


\begin{proposition}[Uniqueness]\label{unique}
Let the domain ${\Omega}_x=\mathbb{T}^d$ or $\mathbb{R}^d$, the constants ${\alpha}_0\in(0,1)$, and the functions $0\le g_1,g_2\lesssim{\mu}^\frac{1}{2}$ be two solutions to \eqref{gnk} in $(0,T]\times{\Omega}_x\times\mathbb{R}^d$ associated with the same initial data $g_{\rm in}\in {\mathcal C}^{{\alpha}_0}({\Omega}_x\times\mathbb{R}^d)$ such that 
\begin{equation*}
\int_{\mathbb{R}^{d}} g_{\rm in}{\mu}^{\frac{1}{2}}\mathop{}\!\mathrm{d} v\ge\lambda {\ \ \rm in\ }{\Omega}_x{\quad\rm and\quad}
0\le g_{\rm in}\lesssim{\mu}^\frac{1}{2}{\ \ \rm in\ }{\Omega}_x\times\mathbb{R}^d.
\end{equation*} 
Then,  $g_1=g_2$ in $[0,T]\times{\Omega}_x\times\mathbb{R}^d$. 
\end{proposition}

\begin{proof}
In view of the lower bound given by Lemma~\ref{lowerlemma1} and Proposition~\ref{lower}, we know that there is some constant $\lambda_*\in(0,1)$ depending only on universal constants, $T$ and the initial data such that
\begin{equation}\label{uniquelower}
\int_{\mathbb{R}^{d}} g_i{\mu}^{\frac{1}{2}}\mathop{}\!\mathrm{d} v\ge\lambda_* {\ \ \rm in\ } [0,T]\times{\Omega}_x, \quad i=1,2.
\end{equation} 
Therefore, we may assume $T=\Lambda^{-1}$ with $\Lambda>1$. 
Let us set the difference $\tilde{g}:=e^{-\frac{|v|^2t}{8}}(g_1-g_2)$. We have to show that $\tilde{g}$ is identically zero. 

In view of the equation~\eqref{gnk}, a direct computation yields that the function $\tilde{g}$ satisfies 
\begin{equation}\label{uniqueiden0}
\begin{split}
&\left(\partial_t+v\cdot\nabla_x\right)\tilde{g} +\frac{|v|^2}{8}\tilde{g}\\
&=e^{-\frac{|v|^2t}{8}}\left(\mathcal{R}[g_1]-\mathcal{R}[g_2]\right)\mathcal{U}[g_1]
+\mathcal{R}[g_2]\left(\mathcal{U}[\tilde{g}] +\frac{t}{2}\,v\cdot\nabla_v\,\tilde{g} +\left(\frac{dt}{4}+\frac{|v|^2t^2}{16}\right)\tilde{g}\right), 
\end{split}
\end{equation}
with the initial condition $\tilde{g}(0,x,v)=0$ in ${\Omega}_x\times\mathbb{R}^d$. 

Let $y\in\mathbb{R}^d$. We introduce a cut-off function $\phi_y(x):=\phi(x-y)$, where $\phi\in C^\infty_c(\mathbb{R}^d)$ is valued in $[0,1]$ such that $\phi|_{B_1}\equiv 1$, $\phi|_{B_2^c}\equiv 0$ and $|\nabla\phi|\lesssim1$ in $\mathbb{R}^d$. 
For any $t\in(0,T]$, integrating the equation~\eqref{uniqueiden0} against $\phi_y^2\tilde{g}$ in ${\Omega}_x\times\mathbb{R}^d$ and applying integration by parts yields that 
\begin{align*}
\frac{1}{2}\int_{{\Omega}_x\times\mathbb{R}^d}\phi_y^2\tilde{g}^2(t)
=\int_0^t\!\int_{{\Omega}_x\times\mathbb{R}^d} \bigg\{(v\cdot\nabla\phi_y)\,\phi_y\tilde{g}^2 -\frac{|v|^2}{8}\phi_y^2\tilde{g}^2
 + e^{-\frac{|v|^2t}{8}}\left(\mathcal{R}[g_1]-\mathcal{R}[g_2]\right)\mathcal{U}[g_1]\,\phi^2\tilde{g} &\nonumber\\
-\mathcal{R}[g_2]\left(\left|\nabla_{v}\left({\mu}^{-\frac{1}{2}}\,\tilde{g}\right)\right|^2\!\mu\phi_y^2
+\frac{dt}{4}\phi_y^2\tilde{g}^2
-\left(\frac{dt}{4}+\frac{|v|^2t^2}{16}\right)\phi_y^2\tilde{g}^2\right)\bigg\}&.
\end{align*}
Since $\mathcal{R}[g_2]\in[0,\Lambda]$, for any $t\in(0,T]$, we have 
\begin{equation*}\label{uniqueiden}
\frac{1}{2}\int_{{\Omega}_x\times\mathbb{R}^d}\phi_y^2\tilde{g}^2(t)
\le \int_0^t\!\int_{{\Omega}_x\times\mathbb{R}^d} \bigg\{|v||\nabla\phi_y|\,\phi_y\tilde{g}^2 -\frac{|v|^2}{16}\phi_y^2\tilde{g}^2
+ {\mu}^{-\frac{1}{4}}\left|\mathcal{R}[g_1]-\mathcal{R}[g_2]\right||\mathcal{U}[g_1]|\phi_y^2|\tilde{g}|\bigg\}. 
\end{equation*}
Due to the elementary inequality $|z^\beta-1|\le|z-1|$ $(z\in\mathbb{R}^+)$ and the lower bound estimate in \eqref{uniquelower}, as well as the boundedness assumption on $g_1,g_2$, we have 
\begin{align*}
\left|\mathcal{R}[g_1]-\mathcal{R}[g_2]\right|\le \mathcal{R}[g_1]^\frac{{\beta}-1}{{\beta}}|\mathcal{R}[\tilde{g}]|^\frac{1}{{\beta}}\le \frac{1}{\lambda_*}\int_{\mathbb{R}^d}|\tilde{g}(t,x,\cdot)|{\mu}^\frac{1}{2}
\lesssim_{\lambda_*}1  {\quad\rm in\ }[0,1]\times{\Omega}_x. 
\end{align*}
It then follows that for any $t\in(0,T]$, 
\begin{align}\label{uniquenon0}
\frac{1}{2}\int_{{\Omega}_x\times\mathbb{R}^d}\phi_y^2\tilde{g}^2(t)
\le& \int_0^t\!\int_{{\Omega}_x\times\mathbb{R}^d} \left(|v||\nabla\phi_y|\,\phi_y\tilde{g}^2 -\frac{|v|^2}{16}\phi_y^2\tilde{g}^2\right)\nonumber\nonumber\\
  &+\frac{1}{\lambda_*}\int_0^t\big\|\mu^{-\frac{3}{8}}\mathcal{U}[g_1]\big\|_{L^\infty_{x,v}}\int_{{\Omega}_x\times\mathbb{R}_v^d} \phi_y^2|\tilde{g}(t,x,v)|\mu^{\frac{1}{8}}
   \int_{\mathbb{R}_\xi^d}|\tilde{g}(t,x,\xi)|{\mu}^\frac{1}{2}\mathop{}\!\mathrm{d}\xi \nonumber\\  
\lesssim&_{\lambda_*} \int_0^t\!\int_{{\Omega}_x\times\mathbb{R}^d} |\nabla\phi_y|^2\tilde{g}^2 + 
\int_0^t\big\|\mu^{-\frac{3}{8}}\mathcal{U}[g_1]\big\|_{L^\infty_{x,v}}\int_{{\Omega}_x\times\mathbb{R}^d} \phi_y^2\tilde{g}^2,   
\end{align}
where we used the Cauchy–Schwarz inequality and Hölder's inequality in the last line. 
Recalling that $\phi_y(x)=\phi(x-y)\in C^\infty_c(\mathbb{R}^d)$ and $|\nabla\phi|\lesssim1$ in $\mathbb{R}^d$, we have 
\begin{align*}
\sup_{y\in\mathbb{R}^d}\int_{{\Omega}_x\times\mathbb{R}^d}|\nabla\phi_y|^2\tilde{g}^2
\lesssim \sup_{y\in\mathbb{R}^d}\int_{{\Omega}_x\times\mathbb{R}^d} \phi_y^2\tilde{g}^2. 
\end{align*}
By the definition of $\mathcal{U}[g_1]$ and the upper bound $g_1\lesssim {\mu}^\frac{1}{2}$, 
\begin{align*}
\big\|\mu^{-\frac{3}{8}}\mathcal{U}[g_1]\big\|_{L^\infty_{x,v}}
\lesssim 1+\big\|\mu^{-\frac{3}{8}}\Delta_{v}g_1\big\|_{L^\infty_{x,v}}. 
\end{align*}
Hence, for any $t\in(0,T]$, taking supremum over $y\in\mathbb{R}^d$ in \eqref{uniquenon0}, we obtain 
\begin{align}\label{uniquenon}
\sup_{y\in\mathbb{R}^d}\int_{{\Omega}_x\times\mathbb{R}^d}\phi_y^2\tilde{g}^2(t)
\lesssim_{\lambda_*} \int_0^t\left(1+\big\|\mu^{-\frac{3}{8}}\Delta_{v}g_1\big\|_{L^\infty_{x,v}}\right) \sup_{y\in\mathbb{R}^d}\int_{{\Omega}_x\times\mathbb{R}^d} \phi_y^2\tilde{g}^2. 
\end{align}

Now we have to consider the pointwise estimate on $D_v^2g_1$. Let $z_0=(t_0,x_0,v_0)\in(0,T]\times{\Omega}_x\times\mathbb{R}^d$ and $2r=t_0^\frac{1}{2}$. 
In view of \eqref{uniquelower}, Lemma~\ref{boundschauder}~(ii) implies that there exists some constant ${{\alpha}_*}\in(0,1)$ with the same dependence as $\lambda_*$ such that 
\begin{align*}
\|g_1\|_{{\mathcal C}_l^{{\alpha}_*}\left([0,T]\times{\Omega}_x\times B_1(v_0)\right)}
\lesssim_{\lambda_*} 1+[g_{\rm in}]_{{\mathcal C}^{{\alpha}_0}({\Omega}_x\times\mathbb{R}^d)}. 
\end{align*}
Then, applying the interior Schauder estimate (Proposition~\ref{interiorschauder}) and the upper bound $g_1\lesssim {\mu}^\frac{1}{2}$ yields that 
\begin{equation*}
\begin{split}
\big\|D_v^2g_1\big\|_{L^\infty(Q_{r}(z_0))}
&\lesssim_{\lambda_*} r^{-2}\| g_1-g_1(z_0)\|_{L^\infty(Q_{2r}(z_0))}
+r^{{\alpha}_*}\left[\mathcal{R}[g_1]\left(\frac{d}{2}-\frac{|v|^{2}}{4}\right)g_1\right]_{{\mathcal C}_l^{{\alpha}_*}\left(Q_{2r}(z_0)\right)}\\
&\lesssim_{\lambda_*} r^{-2+\frac{{{\alpha}_*}}{4}}{\mu}^\frac{3}{8}(v_0)[g_1]_{{\mathcal C}_l^{{\alpha}_*}(Q_{2r}(z_0))}^\frac{1}{4}
 +{\mu}^\frac{3}{8}(v_0)[g_1]_{{\mathcal C}_l^{{\alpha}_*}(Q_{2r}(z_0))}^\frac{1}{4}\\
& \lesssim_{\lambda_*}t_0^{-1+\frac{{\alpha}_*}{8}} {\mu}^\frac{3}{8}(v_0) \Big(1+[g_{\rm in}]_{{\mathcal C}^{{\alpha}_0}({\Omega}_x\times\mathbb{R}^d)}^\frac{1}{4}\Big). 
\end{split}
\end{equation*}
By the arbitrariness of $z_0$, we know that for any $s\in(0,T]$,  
\begin{equation*}
\big\|\mu^{-\frac{3}{8}}\Delta_{v}g_1(s)\big\|_{L^\infty({\Omega}_x\times\mathbb{R}^d)}
\lesssim_{\lambda_*} \Big(1+[g_{\rm in}]^\frac{1}{4}_{{\mathcal C}^{{\alpha}_0}({\Omega}_x\times\mathbb{R}^d)}\Big)\,s^{-1+\frac{{\alpha}_*}{8}}. 
\end{equation*} 
Dragging this estimate into \eqref{uniquenon} yields that for any $t\in(0,T]$,
\begin{align}\label{gronwall}
\sup_{y\in\mathbb{R}^d}\int_{{\Omega}_x\times\mathbb{R}^d}\phi_y^2\tilde{g}^2(t)
\le C_*\int_0^t\mathop{}\!\mathrm{d} s\left(1+s^{-1+\frac{{\alpha}_*}{8}}\right) \sup_{y\in\mathbb{R}^d}\int_{{\Omega}_x\times\mathbb{R}^d}\phi_y^2\tilde{g}^2(s),  
\end{align}
where the constant $C_*>0$ depends only on universal constants and the initial data. 
The desired result is then given by Grönwall's inequality. 
\end{proof}
	
	

	
	

\subsection{Global regularity}\label{smoothsection}
The instantaneous smoothness a priori estimate in Theorem~\ref{wellpose}~(i) is made up of the lower bound 
given by Proposition~\ref{lower} and the following proposition. 

\begin{proposition}\label{smoothness}
Let ${\Omega}_x=\mathbb{T}^d$ or $\mathbb{R}^d$,  
$\underline{T}\in(0,T)$, and the function$g$ be a solution to \eqref{gnk} in $(0,T)\times\mathbb{T}^d\times\mathbb{R}^d$ such that 
\begin{equation}\label{lowerunif}
\mathcal{R}[g]\ge\lambda {\ \ \rm in\ }\left[{\underline{T}}/{4},T\right]\times{\Omega}_x {\quad\rm and\quad}
	0\le g\le\Lambda{\mu}^\frac{1}{2}{\ \ \rm in\ }[0,T]\times{\Omega}_x\times\mathbb{R}^d. 
\end{equation}
Then, for any ${\nu}\in\big(0,\frac{1}{2}\big)$ and $k\in\mathbb{N}$, we have 
\begin{equation*}
	\|{\mu}^{-{\nu}}g\|_{{\mathcal C}^{k}\left([\underline{T}, T]\times{\Omega}_x\times\mathbb{R}^d\right)} \le C_{\underline{T},{\nu},k}, 
\end{equation*}
for some constant $C_{\underline{T},{\nu},k}>0$ depending only on universal constants, $\underline{T},{\nu}$ and $k$. 
\end{proposition}


Generally speaking, if $g$ be a solution to \eqref{gnk} in $(0,T]\times{\Omega}_x\times\mathbb{R}^d$ constructed by Proposition~\ref{existence} (${\Omega}_x=\mathbb{T}^d$) or Corollary~\ref{existwhole} (${\Omega}_x=\mathbb{R}^d$), then the uniform positivity  assumption \eqref{lowerunif} should be replaced by 
\begin{equation*} 
\mathcal{R}[g]\ge\lambda_{t,x} {\ \ \rm in \ } (0,T]\times{\Omega}_x
{\quad\rm and\quad} 0\le g\le\Lambda{\mu}^\frac{1}{2}	{\ \ \rm in \ } {\Omega}_x\times\mathbb{R}^d,
\end{equation*}
where $\lambda_{t,x}>0$ may degenerate to zero as $t\rightarrow0$, or $t+|x|\rightarrow\infty$; see Proposition~\ref{lower} and Remark~\ref{lowerwholeremark}. 
As an immediate consequence of the above proposition, for any ${\nu}\in\big(0,\frac{1}{2}\big)$, $k\in\mathbb{N}$ and for any compact subset $K\subset(0,T]\times{\Omega}_x$, there exists some constant $C_{{\nu},k,K}>0$ depending only on universal constants, ${\nu},k$ and $K$ such that 
\begin{equation*}
	\|{\mu}^{-{\nu}}g\|_{{\mathcal C}^{k}\left(K\times\mathbb{R}^d\right)} \le C_{{\nu},k,K}, 
\end{equation*}
which is exactly the assertion in Theorem~\ref{wellpose} (i). 


In order to show the higher regularity, we will apply the bootstrap procedure developed in \cite{IS} which was intended for the non-cutoff Boltzmann equation. 
The classical bootstrap iteration proceeds by differentiating the equation, using a priori estimates to the new equation to improve the regularity of solutions, and repeating such procedure. 
Nevertheless, since ${\mathcal C}_l^{2+{\alpha}}\not\subset{\mathcal C}_x^1$ for any ${\alpha}\in(0,1)$ by their definitions, the hypoelliptic structure of the equation~\eqref{gnk} lacks gain of enough regularity in $x$-variable which disables the $x$-differentiation at each iteration. 
Indeed, the Schauder type estimate provided by Lemma~\ref{boundschauder}~(i) only shows that the solution to \eqref{gnk} belongs to ${\mathcal C}^\frac{2+{\alpha}}{3}$ with respect to $x$-variable. 
In order to overcome it, we have to apply estimates to increments of the solution to recover a full derivative. 
From now on, for $y \in \mathbb{R}^{d}$ and $w\in\mathbb{R}\times\mathbb{R}^d\times\mathbb{R}^d$, we denote the spatial increment 
\begin{equation*}
{\delta}_yg(z):= g(w\circ (0,y,0))-g(w).
\end{equation*}

Let us proceed with the proof of the regularity estimate. 
\begin{proof}[Proof of Proposition~\ref{smoothness}]
We are going to show that for any multi-index $k:=(k_t,k_x, k_v)\in\mathbb{N}\times\mathbb{N}^d\times\mathbb{N}^d$ and ${\nu}\in\big(0,\frac{1}{2}\big)$, 
there exists some constant ${\alpha}_k\in(0,1)$ depending only on $|k|$ such that for any $Q_{r}(z_0)\subset[\underline{T}/2,T] \times {\Omega}_x\times\mathbb{R}^d$, 
\begin{equation}\label{inequality}
\|\partial_{t}^{k_t}\partial_{x}^{k_x}\partial_{v}^{k_v}g\|_{{\mathcal C}_l^{2+{\alpha}_k}(Q_r(z_{0}))} 
\lesssim_{\underline{T},{\nu},k} \mu^{\nu}(v_0). 
\end{equation}

For simplicity, we will omit the domain in estimates below, since the estimates can be always localized around the center $z_0$. 
		
\medskip\noindent{\textbf{Step 0.}} The case of $k=(0,0,0)$ in \eqref{inequality} is a direct consequence of Lemma \ref{boundschauder}~(i). 
		
		
\medskip\noindent{\textbf{Step 1.}} We will establish that \eqref{inequality} holds for any differential operators of the type $\partial^{k_x}_x$. It suffices to show that for any $n\in\mathbb{N}$, $k_x\in\mathbb{N}^d$ with $|k_x|=n$, ${\nu}\in\big(0,\frac{1}{2}\big)$ and $y \in B_{\frac{r^3}{4}}$,
\begin{equation}\label{xincre}
\|\delta_y\partial^{k_x}_xg\|_{{\mathcal C}_l^{2+{\alpha}_n}} \lesssim_{\underline{T},{\nu},n} |y|\mu^{\nu}(v_0).
\end{equation}
Indeed, sending $y \rightarrow 0$ in \eqref{xincre} will complete this step. 
		
Based on an induction on $|k_x|=n$, we suppose that \eqref{xincre} holds for any $|k_x|\le n-1$, which implies for any $k_x\in\mathbb{N}^d$ with $|k_x|\le n$, 
\begin{equation}\label{step1hypo}
\|\partial^{k_x}_xg\|_{{\mathcal C}_l^{2+{\alpha}_n}} \lesssim_{\underline{T},{\nu},n} \mu^{\nu}(v_0).
\end{equation}
We remark that the induction here begin with \eqref{step1hypo} for $|k_x|=0$, which holds due to the previous step. 
		
Let $q:=\delta_{y}\partial^{k_x}_x g$ with $|k_x|=n$. 
Lemma~\ref{incrementi-IS} and \eqref{step1hypo} gives 
\begin{equation}\label{step1lemma}
\|q\|_{{\mathcal C}_l^{{\alpha}_n}} 
\lesssim \|\partial^{k_x}_x g\|_{{\mathcal C}_l^{2+{\alpha}_n}}\|(0,y,0)\|^{2}
\lesssim_{\underline{T},{\nu},n} |y|^{\frac{2}{3}}\mu^{\nu}(v_0).
\end{equation}
Therefore, we have to enhance the exponent $\frac{2}{3}$ on the right hand side to $1$; as a sacrifice, the H\"older exponent on the left hand side will decrease. 
		
Set ${\tau}_yg(w):=g(w\circ (0,y,0))$ for $y\in\mathbb{R}^d$ and $w\in\mathbb{R}\times\mathbb{R}^d\times\mathbb{R}^d$. 
A direct computation shows that $q$ verifies the equation, 
\begin{equation}\label{step1eq}
(\partial_{t} + v \cdot \nabla_{x}) q = \mathcal{R}[g]\, \mathcal{U}[q] 
+ \sum\nolimits_{|i| \le n \atop i \le k_{x}} \delta_{y}\hat{D_{i}}\mathcal{R}[g]\, \mathcal{U}[{\tau}_yD_{i} g]
+\sum\nolimits_{|i| \le n-1 \atop i \le k_{x}}\hat{D_{i}}\mathcal{R}[g]\, \mathcal{U}[\delta_{y} D_{i} g],
\end{equation}
where the multi-indices such that $i \le k$ means each component of $i$ is lower or equal than the corresponding component of $k$
and $\hat{D_{i}}$ denotes the differential operator satisfying $\partial_{t}^{k_t} D_{x}^{k_{x}} = \hat{D_{i}} \circ D_{i}$. 
		
In view of \eqref{step1hypo} \eqref{step1lemma} and the induction hypothesis, each term in two summations on the right hand side of \eqref{step1eq} is bounded in ${\mathcal C}_l^{{\alpha}_n}$ by $C_n\|(0,y,0)\|^2\mu^{\nu'}(v_0)$ for any $\nu'\in(0,\nu)$. 
Then, by the interior Schauder estimate (Proposition~\ref{interiorschauder}), 
\begin{equation}\label{gain-x}
\|q\|_{{\mathcal C}_l^{2+{\alpha}_n}}\lesssim_{\underline{T},{\nu}',n}\|(0,y,0)\|^2\mu^{\nu'}(v_0).
\end{equation}
Combining Lemma~\ref{mio-lemma} with \eqref{step1lemma} and \eqref{gain-x}, we obtain \eqref{xincre}. 
		
		
\medskip\noindent{\textbf{Step 2.}} For \eqref{inequality} in the case of $k_v=0$, we proceed a bidimensional induction on $(m,n)=(k_t,|k_{x}|)$ such that for any ${\nu}\in\big(0,\frac{1}{2}\big)$, 
\begin{equation}\label{inductive-t}
\| \partial_{t}^{k_{t}} D^{k_{x}}_{x} g \|_{{\mathcal C}_l^{2+{\alpha}_{m,n}}} \lesssim_{\underline{T},{\nu},m,n}\mu^{\nu}(v_0).
\end{equation}
Based on the previous step ($m=0$), we have to show that \eqref{inductive-t} holds for $k_t=m\ge1$, $|k_{x}| =n$, under the induction hypothesis that \eqref{inductive-t} holds for any $k_t\le m-1$ and $|k_x|\le n+1$, 
		
With $k_t=m>0$, $|k_{x}| =n$, set $q:=\partial_{t}^{k_t} D_{x}^{k_{x}} g$. Then, there holds
\begin{equation}\label{eqstep2}
(\partial_{t} + v \cdot \nabla_{x}) q = \mathcal{R}[g]\, \mathcal{U}[q] + \sum\nolimits_{i \le (k_t, k_{x}, 0) \atop i \ne (k_t, k_{x}, 0)}\hat{D_{i}}\mathcal{R}[g]\, \mathcal{U}[D_{i} g],
\end{equation}
where we use the notation $\hat{D_{i}}$ such that $\partial_{t}^{k_t} D_{x}^{k_{x}} = \hat{D_{i}} \circ D_{i}$. 
		
By the induction hypothesis, each term in the remainder (the summation on the right hand side of \eqref{eqstep2}) with $i \ne (0,0,0)$ can be controlled in ${\mathcal C}_l^{{\alpha}_{m,n}}$. 
It now suffices to deal with the exceptional term $\partial_{t}^{k_t} D_{x}^{k_{x}} \mathcal{R}[g]$ so that the whole remainder can be controlled in ${\mathcal C}_l^{{\alpha}_{m,n}}$; and then \eqref{inductive-t} follows from the interior Schauder estimate (Proposition~\ref{interiorschauder}). 
To this end, using Lemma~\ref{kinetic-degree} and the induction hypothesis with the pair $(m-1,n)$ yields 
\begin{equation}\label{1}
\| (\partial_{t} + v \cdot \nabla_{x}) \partial_{t}^{m-1} D^{k_{x}}_{x} g \|_{{\mathcal C}_l^{{\alpha}_{m,n}}} 
\lesssim_{\underline{T},{\nu},m,n}\mu^{{\nu}}(v_0).
\end{equation}
Due to the induction hypothesis with the pair $(m-1,n+1)$, for any $\nu'\in(0,\nu)$, 
\begin{align}\label{2}
\mu^{-\nu'}(v_0)\| (v \cdot \nabla_{x}) \partial_{t}^{m-1} D^{k_{x}}_{x} g \|_{{\mathcal C}_l^{2+{\alpha}_{m,n}}} 
\lesssim_{{\nu},{\nu}'} \mu^{-\nu}(v_0)\| \partial_{t}^{m-1} \nabla_{x} D^{k_{x}}_{x} g \|_{{\mathcal C}_l^{2+{\alpha}_{m,n}}} 
\lesssim_{\underline{T},{\nu},m,n} 1.
\end{align}
Then, \eqref{1} and \eqref{2} produce the bound on $\mu^{-\nu'}(v_0)\|q\|_{{\mathcal C}_l^{{\alpha}_{m,n}}}$. 
		
		
\medskip\noindent{\textbf{Step 3.}} Similarly, to show \eqref{inequality} for any differential operator $\partial_{t}^{k_{t}} D^{k_{x}}_{x} D^{k_{v}}_{v}$, 
we proceed a bidimensional induction on $(m,n)=(k_t+|k_x|,k_{v})$ such that for any ${\nu}\in\big(0,\frac{1}{2}\big)$,
\begin{equation}\label{inductive-v}
\| \partial_{t}^{k_{t}} D^{k_{x}}_{x} D^{k_{v}}_{v} g \|_{{\mathcal C}_l^{2+{\alpha}_{m,n}}} 
\lesssim_{\underline{T},{\nu},m,n}\mu^{\nu}(v_0).
\end{equation}
The case $n=0$ is treated in the previous step. By Lemma~\ref{kinetic-degree} and the induction hypothesis \eqref{inductive-v} with $k_{t} + |k_{x}| = m$ and $|k_v|=n-1$, $n \ge 1$, we have
\begin{equation*}
\| \partial_v  \partial_t^{k_{t}} \partial_x^{k_{x}} \partial_v^{k_v} g \|_{{\mathcal C}_l^{{\alpha}_{m,n}}} 
\lesssim \| \partial_t^{k_{t}} \partial_x^{k_{x}} \partial_v^{k_v} g \|_{{\mathcal C}_l^{1+{\alpha}_{m,n}}}  
\lesssim_{\underline{T},{\nu},m,n} {\mu}^{{\nu}}(v_0).
\end{equation*}
Computing the equation satisfied by $\partial_v \partial_{t}^{k_{t}} \partial^{k_{x}}_{x} \partial^{n-1}_{v} g$ and proceeding the argument as in the previous step, we conclude the proof.
\end{proof}
	




\section{Diffusion asymptotics}\label{diffusionlimit}
This section is devoted to the study of the global in time quantitative diffusion asymptotics which consists of the (uniform-in-${\epsilon}$) convergence towards the equilibrium over long times and of the finite time asymptotics, including the results of Theorem~\ref{wellpose}~(ii) and Theorem~\ref{limit}. 

First of all, let us introduce the required notation. 
For any scalar or vector valued function $\Psi\in L^1(\mathbb{R}^d,\mathop{}\!\mathrm{d}\mu)$, we denote its velocity mean by
	\begin{equation*}
	\langle \Psi\rangle:=\int_{\mathbb{R}^d} \Psi(v) \mathop{}\!\mathrm{d}\mu. 
	\end{equation*}
For any couple of functions (scalar, vectors or $d\times d$-matrices) $\Psi_1,\Psi_2\in L^2\left(\mathbb{T}^d\times\mathbb{R}^d,\mathop{}\!\mathrm{d} m\right)$, we denote their $L^2$ inner product with respect to the measure $\mathop{}\!\mathrm{d} m$ by
\begin{equation*}
\left(\Psi_1, \Psi_2\right):=\int_{\mathbb{T}^d\times\mathbb{R}^d} \Psi_1(x,v)\Psi_2(x,v)\mathop{}\!\mathrm{d} m, 
\end{equation*}
where the multiplication between the couple in the integrand is replaced by scalar contraction product, if $\Psi_1,\Psi_2$ is a couple of vectors or matrices.  
	
Recalling our notation for the Ornstein-Uhlenbeck operator  $\mathscr{L}_{\scriptscriptstyle{\rm OU}}=\left(\nabla_v-v\right)\cdot\nabla_v$, we apply the substitutions $f_{\epsilon}={\mu} h_{\epsilon}$, $f_{{\epsilon},\rm in}=\mh_{{\epsilon},\rm in}$ to \eqref{epfnk} and obtain 
\begin{equation}\label{epnk}
\left\{ 
\begin{aligned}
\ &\left({\epsilon}\partial_t+v\cdot\nabla_x\right) h_{\epsilon}(t,x,v) = \frac{1}{{\epsilon}}\langle h_{\epsilon}\rangle^{\beta}(t,x)\,\mathscr{L}_{\scriptscriptstyle{\rm OU}} h_{\epsilon}(t,x,v), \\
\ &\;h_{\epsilon}(0,x,v) = h_{{\epsilon},\rm in} (x,v), \\
\end{aligned}
\right. 
\end{equation} 
In this setting, by applying integration by parts, for any $h_1,h_2\in  \mathcal{C}_c^\infty(\mathbb{T}^d\times\mathbb{R}^d)$ we get
\begin{equation*}
\left(h_1,\, \mathscr{L}_{\scriptscriptstyle{\rm OU}} h_2\right)=-\left(\nabla_vh_1,\, \nabla_vh_2\right). 
\end{equation*}
We will use this identity repeatedly in the computation below. 	
Then, the operator $\mathscr{L}_{\scriptscriptstyle{\rm OU}}$ is self-adjoint with respect to the inner product $(\cdot,\cdot)$ and the bracket $\langle \cdot\rangle$ is a projection on the null space of $\mathscr{L}_{\scriptscriptstyle{\rm OU}}$. 
Moreover, as the total mass is conserved, we define
	\begin{equation}\label{FP2cond}
	M_0:=\int_{\mathbb{T}^d\times\mathbb{R}^d} h_{\epsilon}\mathop{}\!\mathrm{d} m=\int_{\mathbb{T}^d}\langle h_{\epsilon}\rangle\mathop{}\!\mathrm{d} x.
	\end{equation}
 Proceeding with the macro-micro (fluid-kinetic) decomposition, we define the orthogonal complement of the projection~$\langle \cdot\rangle$ of $h_{\epsilon}$ as 
	$$h^\perp_{\epsilon}(t,x,v):=h_{\epsilon}(t,x,v)-\langle h_{\epsilon}\rangle(t,x).$$ 
	In this framework, the local mass $\langle h_{\epsilon}\rangle$ is the macroscopic (fluid) part and the complement $h^\perp_{\epsilon}$ is the microscopic (kinetic) part. 
Besides, taking the bracket $\langle\cdot\rangle$ after multiplying the equation in \eqref{epnk} with $1$ and $v$ leads to the following macroscopic equations, 
	\begin{align}
	&{\epsilon}\partial_t\langle h_{\epsilon}\rangle+\nabla_x\cdot\langle vh_{\epsilon}\rangle=0,  \label{hydroFP1} \\
	&{\epsilon}\partial_t\langle vh_{\epsilon}\rangle + \nabla_x\cdot\langle v^{\otimes2}h_{\epsilon}\rangle
	=-\frac{1}{{\epsilon}}\langle h_{\epsilon}\rangle^\beta \langle vh_{\epsilon}\rangle,   \label{hydroFP2}
	\end{align}
	where $\langle vh_{\epsilon}\rangle$ and $\langle v^{\otimes2}h_{\epsilon}\rangle$ represent the local momentum and the stress tensor, respectively. 
	
	
\subsection{Long time behavior}\label{longtimesection}
Our aim is to establish the (uniform-in-${\epsilon}$) exponential decay towards the equilibrium $M_0$ for \eqref{epnk}. 
In particular, when ${\epsilon}=1$, it sets up the exponential convergence in each order derivative based on the smoothness a priori estimates given in Subsection~\ref{smoothsection}. 
	
We remark that the classical coercive method is not applicable in our case to obtain the convergence to equilibrium due to the degeneracy of the ellipticity of the spatially inhomogeneous equation. Indeed, the Poincaré inequality only produces a spectral gap on the orthogonal complement of the projection $\langle\cdot\rangle$; see \eqref{hypoestimate0} below. 
As we mentioned in Subsection~\ref{background}, there are several ways to achieve the long time asymptotics. 
We will mainly follow the argument presented in \cite{EGKM} (see also \cite{KGH}) in a simpler scenario. 
It would also allow us to see some similarity among \cite{EGKM,DMS,He}. 


\begin{proposition}\label{longtime}
Let the function $\lambda_t:\mathbb{R}_+\rightarrow[0,\Lambda]$ with the derivative $\lambda_t'\le0$ on $\mathbb{R}_+$. 
If $h_{\epsilon}$ is a solution to \eqref{epnk} in $\mathbb{R}_+\times\mathbb{T}^d\times\mathbb{R}^d$, associated with the initial data $0\leh_{{\epsilon},\rm in}\le\Lambda$, satisfying 
\begin{equation}\label{hp-prop}
\langle h_{\epsilon}\rangle^\beta(t,x)\ge\lambda_t {\ \ \rm in\ }\mathbb{R}_+\times\mathbb{T}^d {\quad\rm and\quad}
\int_{\mathbb{R}_+}\left(\lambda_t+\lambda_t'\right)\mathop{}\!\mathrm{d} t=\infty,  
\end{equation}
then the solution $h_{\epsilon}$ converges to the state $M_0$ in $L^2(\mathop{}\!\mathrm{d} m)$ as $t\rightarrow\infty$; more precisely, there exists some universal constant $c>0$ such that for any $t>0$, we have 
\begin{equation}\label{longtimeestimate}
\|h_{\epsilon}(t,\cdot,\cdot)-M_0\|^2_{L^2(\mathop{}\!\mathrm{d} m)}
\lesssim \| h_{{\epsilon},\rm in}-M_0\|^2_{L^2(\mathop{}\!\mathrm{d} m)}\exp\left(-c\int_0^t\left(\lambda_s+\lambda_s'\right)\mathop{}\!\mathrm{d} s\right).  
\end{equation}
\end{proposition}
	
	
\begin{proof}
Since the velocity mean of the microscopic part vanishes, $\langle h^\perp_{\epsilon}\rangle=0$, using the equation~\eqref{epnk} and the Poincaré inequality yields that
\begin{align}\label{hypoestimate0}
\frac{1}{2}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\| h_{\epsilon}-M_0\|^2_{L^2(\mathop{}\!\mathrm{d} m)}
&=\frac{1}{{\epsilon}^2} \left(\langle h_{\epsilon}\rangle^{\beta}\mathscr{L}_{\scriptscriptstyle{\rm OU}} h_{\epsilon},\, h_{\epsilon}-M_0\right) =-\frac{1}{{\epsilon}^2}\left(\langle h_{\epsilon}\rangle^{\beta}\nabla_vh^\perp_{\epsilon},\, \nabla_vh^\perp_{\epsilon}\right)\nonumber\\
&\le-\frac{\lambda_t}{{\epsilon}^2}\| \nabla_v h^\perp_{\epsilon}\|^2_{L^2(\mathop{}\!\mathrm{d} m)}
\lesssim -\frac{\lambda_t}{{\epsilon}^2}\| h^\perp_{\epsilon}\|^2_{L^2(\mathop{}\!\mathrm{d} m)}. 
\end{align}
Now we have to recover a new entropy that would give some bound on the projection~$\langle h_{\epsilon}\rangle-M_0$. 
		
For every test function $v\cdot\Psi(t,x)\mu$, with a vector-valued function $\Psi\in H^1_{t,x}(\mathbb{R}^{+}\times\mathbb{T}^d,\mathbb{R}^d)$, we write the weak formulation of \eqref{epnk} as follows	
\begin{equation*}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\Psi,\, h_{\epsilon}\right)=\frac{1}{{\epsilon}} \left( v^{\otimes2}:\nabla_x\Psi,\, h_{\epsilon}\right) +\left(v\cdot\partial_t\Psi,\, h_{\epsilon}\right) 
+\frac{1}{{\epsilon}^2}\left(\langle h_{\epsilon}\rangle^\beta \mathscr{L}_{\scriptscriptstyle{\rm OU}} v\cdot\Psi,\, h_{\epsilon}\right). 
\end{equation*}
Taking the macro-micro decomposition into account, we obtain from the above expression 
\begin{equation}\label{FPpert2}
\begin{split}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\Psi,\, h^\perp_{\epsilon}\right)
=&\frac{1}{{\epsilon}}\left(|v_1|^2{\text{\rm tr}}(\nabla_x\Psi),\, \langle h_{\epsilon}\rangle-M_0\right)+\frac{1}{{\epsilon}}\left(v^{\otimes2}:\nabla_x\Psi,\, h^\perp_{\epsilon}\right)\\
&+\left(v\cdot\partial_t\Psi,\, h^\perp_{\epsilon}\right)
-\frac{1}{{\epsilon}^2}\left(\langle h_{\epsilon}\rangle^\beta v\cdot\Psi,\, h^\perp_{\epsilon}\right). 
\end{split}
\end{equation}
Let us now introduce an auxiliary function $u(t,x)$: for any fixed $t\in\mathbb{R}_+$, $u(t,x)$ is defined as the solution of the following elliptic equation under the compatibility condition \eqref{FP2cond}, 
\begin{equation}\label{hypoPoisson0} 
	-\Delta_x u=\langle h_{\epsilon}\rangle-M_0 {\quad\rm in\ }\mathbb{T}^d, 
\end{equation}
whose elliptic estimate states 
\begin{equation}\label{hypoPoisson1}
\|\nabla_x u\|_{L^2_x}+\|\nabla^2_x u\|_{L^2_x}
\lesssim \|\langle h_{\epsilon}\rangle-M_0\|_{L^2_x}. 
\end{equation}
Besides, observing that $\langle vh_{\epsilon}\rangle$=$\langle vh^\perp_{\epsilon}\rangle$, from \eqref{hydroFP1}, we get 
\begin{equation*}
{\epsilon}\partial_t\langle h_{\epsilon}\rangle+\nabla_x\cdot\langle vh^\perp_{\epsilon}\rangle=0. 
\end{equation*}
Combining this macroscopic relation with \eqref{hypoPoisson0}, we have
\begin{equation*}
\int_{\mathbb{T}^d}|\nabla_x\left(\partial_tu\right)|^2
		=\int_{\mathbb{T}^d}\partial_tu\,\partial_t\langle h_{\epsilon}\rangle
		=-\frac{1}{{\epsilon}}\int_{\mathbb{T}^d}\partial_tu\,\nabla_x\cdot\langle vh^\perp_{\epsilon}\rangle
		=\frac{1}{{\epsilon}}\int_{\mathbb{T}^d}\nabla_x\left(\partial_tu\right)\cdot\langle vh^\perp_{\epsilon}\rangle. 
\end{equation*}
It then follows from H\"older's inequality that 
\begin{equation}\label{hypoPoisson2}
\|\nabla_x(\partial_tu)\|_{L^2_x}
\le\frac{1}{{\epsilon}}\|\langle vh^\perp_{\epsilon}\rangle\|_{L^2_x}
\lesssim \frac{1}{{\epsilon}}\| h^\perp_{\epsilon}\|_{L^2(\mathop{}\!\mathrm{d} m)}. 
\end{equation}
Choosing $\Psi=\nabla_xu$ in \eqref{FPpert2} yields
\begin{equation*}\label{recover}
\begin{split}
-\frac{1}{{\epsilon}}\left(|v_1|^2\Delta_xu,\langle h_{\epsilon}\rangle-M_0\right)
&+\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\nabla_xu,\, h^\perp_{\epsilon}\right)\\
&\lesssim \left(\frac{1}{{\epsilon}}\| \nabla^2_xu\|_{L^2_x}
+\|\nabla_x\left(\partial_tu\right)\|_{L^2_x}
+\frac{1}{{\epsilon}^2}\|\nabla_xu\|_{L^2_x}\right)\| h^\perp_{\epsilon}\|_{L^2(\mathop{}\!\mathrm{d} m)}.
\end{split}
\end{equation*}
Applying \eqref{hypoPoisson0}, \eqref{hypoPoisson1} and \eqref{hypoPoisson2}, we have
\begin{equation*}
\begin{split}
\frac{1}{{\epsilon}}\|\langle h_{\epsilon}\rangle-M_0\|^2_{L^2_x}
+\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\nabla_xu,\, h^\perp_{\epsilon}\right)
&\lesssim \frac{1}{{\epsilon}^2}\|\langle h_{\epsilon}\rangle-M_0\|_{L^2_x}\| h^\perp_{\epsilon}\|_{L^2(\mathop{}\!\mathrm{d} m)}
+\frac{1}{{\epsilon}}\| h^\perp_{\epsilon}\|_{L^2(\mathop{}\!\mathrm{d} m)}^2.
\end{split}
\end{equation*}
By the Cauchy–Schwarz inequality, we arrive at 
\begin{equation}\label{necessario}
\|\langle h_{\epsilon}\rangle-M_0\|^2_{L^2_x}
	+{\epsilon}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\nabla_xu,\, h^\perp_{\epsilon}\right)
\lesssim \frac{1}{{\epsilon}^2}\| h^\perp_{\epsilon}\|_{L^2(\mathop{}\!\mathrm{d} m)}^2.
\end{equation}

Then, \eqref{necessario} combined with \eqref{hypoestimate0} implies that 
\begin{align*}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathcal{E}_{\epsilon}(t)
&\lesssim -\frac{1-{\delta}}{{\epsilon}^2}\| h^\perp_{\epsilon}\|_{L^2(\mathop{}\!\mathrm{d} m)}^2
  -{{\delta}\lambda_t}\|\langle h_{\epsilon}\rangle-M_0\|^2_{L^2_x}
  +{\delta}{\epsilon}\lambda_t'\left(v\cdot\nabla_xu,\, h^\perp_{\epsilon}\right)\\
&\le -{\delta}\lambda_t\| h_{\epsilon}-M_0\|^2_{L^2(\mathop{}\!\mathrm{d} m)} 
  -{\delta}\lambda_t'\left|\left(v\cdot\nabla_xu,\, h^\perp_{\epsilon}\right)\right|,
\end{align*}
where the constant $\delta\in\big(0,\frac{1}{2}\big)$ will be determined and the modified entropy $\mathcal{E}_{\epsilon}$ is defined by 
$$\mathcal{E}_{\epsilon}(t):= \| h_{\epsilon}-M_0\|^2_{L^2(\mathop{}\!\mathrm{d} m)} +{{\delta}{\epsilon}\lambda_t}\left(v\cdot\nabla_xu,\,  h^\perp_{\epsilon}\right).$$ 	
Noticing \eqref{hypoPoisson1} also implies that  
\begin{equation}\label{equw}
\left|\left(v\cdot\nabla_xu,\,  h^\perp_{\epsilon}\right)\right|
\lesssim \|\langle h_{\epsilon}\rangle-M_0\|_{L^2_x}\| h^\perp_{\epsilon}\|_{L^2(\mathop{}\!\mathrm{d} m)}
\le \| h_{\epsilon}-M_0\|^2_{L^2(\mathop{}\!\mathrm{d} m)}. 
\end{equation}
It means that the modified entropy $\mathcal{E}_{\epsilon}$ is equivalent (independent of ${\epsilon}$) to the square of the $L^2(\mathop{}\!\mathrm{d} m)$-distance between $h_{\epsilon}$ and $M_0$, when the constant $\delta>0$ is sufficiently small. 
		
Hence, we have 
\begin{equation*}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t} \mathcal{E}_{\epsilon}(t) 
\lesssim -\left(\lambda_t+\lambda_t'\right) \mathcal{E}_{\epsilon}(t). 
\end{equation*}
The conclusion \eqref{longtimeestimate} then follows from Grönwall's inequality and the equivalence between $\mathcal{E}_{\epsilon}(t)$ and $\| h_{\epsilon}(t,\cdot,\cdot)-M_0\|^2_{L^2(\mathop{}\!\mathrm{d} m)}$. 
\end{proof}
		





We pointed out that the elliptic estimate \eqref{hypoPoisson1} for the Poisson equation \eqref{hypoPoisson0} used in the above proof resulting from a Poincaré type inequality essentially relies on the compactness of the spatial domain. 
It was shown in \cite{BDM} that the related elliptic estimate can be recovered by applying the Nash inequality \cite{Nash} when the spatial domain is the whole space $\mathbb{R}^d$, whose argument is under an abstract setting. 
Inspired by the proof of Proposition~\ref{longtime} above, we are also able to make the construction of \cite{BDM} precise to see that the argument still works for the nonlinear equation \eqref{epnk}. 
We remark that the following algebraic decay rate is optimal in the sense it is the same as in the linear case; see Appendix A of \cite{BDM}. 

\begin{proposition}\label{longwhole}
Assume that the initial data $h_{{\epsilon},\rm in}$ is valued in $[\lambda,\Lambda]$, and satisfies $h_{{\epsilon},\rm in}-M_1\in  L^1(\mathbb{R}^{2d},\mathop{}\!\mathrm{d} m)$ for some universal constant $M_1>0$. 
Let the function $h_{\epsilon}$ valued in $[\lambda,\Lambda]$ be a solution to \eqref{epnk} in $\mathbb{R}_+\times\mathbb{R}^{2d}$ associated with $h_{\epsilon}|_{t=0}=h_{{\epsilon},\rm in}$. Then, for any $t>0$, 
\begin{equation*}
\|h_{\epsilon}-M_1\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}
\lesssim \left(1+\|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}\right) t^{-\frac{d}{4}}. 
\end{equation*}
\end{proposition}


\begin{proof}
By the same derivation of \eqref{hypoestimate0} and \eqref{FPpert2} as in the proof of Proposition~\ref{longtime}, we have the microscopic coercivity
\begin{align}\label{hypo0w}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\| h_{\epsilon}-M_1\|^2_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}
\lesssim -\frac{1}{{\epsilon}^2}\| h^\perp_{\epsilon}\|^2_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}, 
\end{align}
and the identity from the macro-micro decomposition that 
\begin{equation}\label{hypo1w}
\begin{split}
-\left(|v_1|^2\Delta_x w,\, \langle h_{\epsilon}\rangle-M_0\right)_W &+{\epsilon}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\nabla_xw,\, h^\perp_{\epsilon}\right)_W
=\left(v^{\otimes2}:\nabla_x^2w,\, h^\perp_{\epsilon}\right)_W\\
&+{\epsilon}\left(v\cdot\partial_t\nabla_xw,\, h^\perp_{\epsilon}\right)_W-\frac{1}{{\epsilon}}\left(\langle h_{\epsilon}\rangle^\beta v\cdot\nabla_xw,\, h^\perp_{\epsilon}\right)_W, 
\end{split}
\end{equation}
where $(\cdot,\cdot)_W$ denotes the $L^2(\mathbb{R}^{2d},\mathop{}\!\mathrm{d} m)$ inner product, and the function $w(t,x)\in L_t^\infty([0,T];L_x^1\cap L_x^2(\mathbb{R}^d))$ is chosen to be the solution of the following ellipitc equation associated with the constant $\Theta:=\langle |v_1|^2\rangle$ and the macroscopic source $\langle h_{\epsilon}\rangle-M_1$, 
\begin{equation}\label{poissonw}
w-\Theta\Delta_x w=\langle h_{\epsilon}\rangle-M_1 {\quad\rm in\ }\mathbb{R}^d. 
\end{equation}
The elliptic estimate is derived by integrating \eqref{poissonw} against $-\Theta\Delta_x w$ so that 
\begin{equation}\label{ellipticw}
\begin{split}
\Theta\|\nabla_xw\|_{L_x^2(\mathbb{R}^d)}^2 + \Theta^2\|\nabla_x^2w\|_{L_x^2(\mathbb{R}^d)}^2
&= \left(-\Theta\Delta_x w,\, \langle h_{\epsilon}\rangle-M_1\right)_W\\
&= \left(\langle h_{\epsilon}\rangle-M_1 - w,\, \langle h_{\epsilon}\rangle-M_1\right)_W=:\mathcal{A}. 
\end{split}
\end{equation}
It also follows from the same derivation as \eqref{hypoPoisson2} that 
\begin{equation*}
\|\nabla_x(\partial_tw)\|_{L_x^2(\mathbb{R}^d)} \lesssim \frac{1}{{\epsilon}}\|h^\perp_{\epsilon}\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}. 
\end{equation*}
Gathering the above two estimates with \eqref{hypo1w}, we obtain 
\begin{equation*}
\begin{split}
\mathcal{A} + {\epsilon} \frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\nabla_xw,\, h^\perp_{\epsilon}\right)_W
\lesssim \frac{1}{{\epsilon}} \mathcal{A}^\frac{1}{2}\| h^\perp_{\epsilon}\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}
+\|h^\perp_{\epsilon}\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2, 
\end{split}
\end{equation*}
which implies from the Cauchy–Schwarz inequality that 
\begin{equation*}
\mathcal{A} + {\epsilon} \frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\left(v\cdot\nabla_xw,\, h^\perp_{\epsilon}\right)_W
\lesssim \frac{1}{{\epsilon}^2} \| h^\perp_{\epsilon}\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2 . 
\end{equation*}
Denoting the modified entropy $\mathscr{E}_{\epsilon}(t):=\|h_{\epsilon}-M_1\|_{L^2({\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m})}^2+{\delta}{\epsilon}\left(v\cdot\nabla_xw,\, h^\perp_{\epsilon}\right)_W$ with a sufficient small constant ${\delta}>0$, and using \eqref{hypo0w}, we conclude that  
\begin{equation}\label{eew}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathscr{E}_{\epsilon}(t)
\lesssim  -\mathcal{A} -\| h^\perp_{\epsilon}\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2 . 
\end{equation}


Recalling the similar estimate \eqref{equw}, we see that $\mathscr{E}_{\epsilon}$ is equivalent to the square of the $L^2(\mathbb{R}^{2d},\mathop{}\!\mathrm{d} m)$-distance between $h_{\epsilon}$ and $M_1$. 
It thus suffices to recover $\langle h_{\epsilon}\rangle-M_1$ by means of $\mathcal{A}$. 
By \eqref{poissonw} and the convexity of $|\cdot|$, we know that $|w|$ is a subsolution in the sense that 
\begin{equation*}
|w|-\Theta\Delta_x |w| \le |\langle h_{\epsilon}\rangle-M_1| {\quad\rm in\ }\mathbb{R}^d, 
\end{equation*}
and hence $\|w\|_{L_x^1(\mathbb{R}^d)} \le \|\langle h_{\epsilon}\rangle-M_1\|_{L_x^1(\mathbb{R}^d)}$.  
With the aid of Corollary~\ref{existwhole}, 
\begin{equation*}
\|w\|_{L_x^1(\mathbb{R}^d)} \le \|h-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}
\le \|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}. 
\end{equation*}
Applying \eqref{ellipticw} and the Nash inequality $\|w\|_{L_x^2(\mathbb{R}^d)}^{d+2}\lesssim \|w\|_{L_x^1(\mathbb{R}^d)}^2\|\nabla_x w\|_{L_x^2(\mathbb{R}^d)}^d$, we then acquire 
\begin{align*}
\|\langle h_{\epsilon}\rangle-M_1\|_{L_x^2(\mathbb{R}^d)}^2
&\lesssim \mathcal{A} +\|w\|_{L_x^2(\mathbb{R}^d)}^2
\lesssim \mathcal{A} + \|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^{\frac{4}{d+2}} \|\nabla_xw\|_{L_x^2(\mathbb{R}^d)}^{\frac{2d}{d+2}} \\
&\lesssim  \Big(\mathcal{A}^\frac{2}{d+2} + \|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^{\frac{4}{d+2}}\Big)  \mathcal{A}^{\frac{d}{d+2}}. 
\end{align*}
Now that $\|h_{\epsilon}-M_1\|_{L_x^2(\mathbb{R}^d)}^2\le (\Lambda+M_1)\|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}$, 
in the both cases of $\mathcal{A} + \|h^\perp_{\epsilon}\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2 \lessgtr \|h_{\epsilon}-M_1\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2$, we conclude that 
\begin{align*}
\|h_{\epsilon}-M_1\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2
\lesssim \Big(1+\|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^{\frac{4}{d+2}}\Big) 
(\mathcal{A}+\|h^\perp_{\epsilon}\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2 )^{\frac{d}{d+2}}. 
\end{align*}
Combining this with \eqref{eew} and the equivalence between $\mathscr{E}_{\epsilon}$ and $\|h_{\epsilon}-M_1\|_{L^2(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2$, we have 
\begin{equation*}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathscr{E}_{\epsilon}(t) 
\lesssim -\Big(1+\|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^{\frac{4}{d}}\Big)^{\!-1}
\mathscr{E}_{\epsilon}(t) ^{1+\frac{2}{d}}, 
\end{equation*}
Since $\mathscr{E}_{\epsilon}(0)\lesssim \|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}$, we arrive at 
\begin{align*}
\mathscr{E}_{\epsilon}(t)
&\lesssim \Big[\mathscr{E}_{\epsilon}(0)^{-\frac{2}{d}} + \Big(1+\|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^{\frac{4}{d}}\Big)^{\!-1} t\Big]^{-\frac{d}{2}}\\
&\lesssim \Big(1+\|h_{{\epsilon},\rm in}-M_1\|_{L^1(\mathbb{R}^{2d},\,\mathop{}\!\mathrm{d} m)}^2\Big) t^{-\frac{d}{2}}. 
\end{align*}
The proof is thus complete. 
\end{proof}






As far as the case ${\epsilon}=1$ is concerned, we conclude the result of convergence to equilibrium. 



\begin{proof}[Proof of Theorem \ref{wellpose} (ii)] 
Consider $g:=\mu^\frac{1}{2}h$. In view of Proposition~\ref{existence} and Corollary~\ref{existwhole} with the assumption on initial data, we know that $\lambda\le \mu^{-\frac{1}{2}}g\le\Lambda$ in $\mathbb{R}_+\times{\Omega}_x\times\mathbb{R}^d$, for ${\Omega}_x=\mathbb{T}^d$ or $\mathbb{R}^d$. 
By using Proposition~\ref{longtime} to $h=\mu^{-\frac{1}{2}}g$ with $\lambda_t=\lambda$ and ${\Omega}_x=\mathbb{T}^d$, we have an universal constant $c>0$ such that 
\begin{equation*}\label{gglong}
\big\| g(t)-M_0{\mu}^\frac{1}{2}\big\|_{L^2(\mathbb{T}^d\times\mathbb{R}^d)}\lesssim e^{-2ct}.
\end{equation*}
Combining this with the Sobolev embedding and the interpolation, we derive that, for any $k\in\mathbb{N}$ with $k\ge d$, 
\begin{align*}
\big\|g(t)-M_0\mu^\frac{1}{2}\big\|_{{\mathcal C}^{k}(\mathbb{T}^d\times\mathbb{R}^d)}
&\lesssim_k \big\|g(t)-M_0\mu^\frac{1}{2}\big\|_{H^{2k}(\mathbb{T}^d\times\mathbb{R}^d)}\\
&\lesssim_k \big\|g(t)-M_0\mu^\frac{1}{2}\big\|_{H^{4k}(\mathbb{T}^d\times\mathbb{R}^d)}^\frac{1}{2}
\big\| g(t)-M_0\mu^\frac{1}{2}\big\|_{L^2(\mathbb{T}^d\times\mathbb{R}^d)}^\frac{1}{2}. 
\end{align*}
Since the $H^{4k}$-norm on the right hand side is bounded due to the global regularity estimate given by Proposition~\ref{smoothness}, we obtain the exponential convergence to equilibrium in each order derivative. 

The asserted result in the case ${\Omega}_x=\mathbb{R}^d$ is a direct consequence of Proposition~\ref{longwhole}. 
As a side remark, one is also able to upgrade the long time convergence to higher order derivatives of solutions by means of the global regularity estimate and interpolation as above; it yet gives the algebraic decay rate that is not optimal. 
\end{proof}




	
\subsection{Finite time asymptotics}\label{finitetimesection}
The study of macroscopic dynamics for the nonlinear kinetic model~\eqref{epnk} in this subsection relies on the regularity of the target equation~\eqref{limiteq}. 
On account of this, let us begin with mentioning some standard results for the equation~\eqref{limiteq} without proof. 
If the initial data satisfies $\lambda\le{\rho}_{\rm in}\le\Lambda$, then such bounds are preserved along times, $\lambda\le{\rho}\le\Lambda$, in the same spirit as Lemma~\ref{Gaussian}. Combining the parabolic De~Giorgi-Nash-Moser theory with the Schauder theory, we know that the solution ${\rho}$ is smooth for any positive time. 
We state the a priori estimate precisely as follows, where its behavior near the initial time is taken into account in view of the standard scaling technique. 

\begin{lemma}\label{porouslemma}
Let $\rho_{\rm in}\in {\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)$ valued in $[\lambda,\Lambda]$ with ${\alpha}_0\in(0,1)$, and ${\rho}$ be the solution to \eqref{limiteq} in $\mathbb{R}_+\times\mathbb{T}^d$. 
Then, there is some universal constant ${\alpha}\in(0,1)$ such that 
\begin{equation}\label{porous0}
\left\|{\rho}\right\|_{\mathcal{C}^{\alpha}(\mathbb{R}_+\times\mathbb{T}^d)}
\lesssim 1+\left\|\rho_{\rm in}\right\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}.
\end{equation}
Moreover, there exists some constant $C_{\rho}>0$ depending only on universal constants and $\left\|\rho_{\rm in}\right\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}$, such that for any $t\in(0,1]$ and $x\in\mathbb{T}^d$, we have 
	\begin{equation}\label{porous}
	{t}^\frac{1-{\alpha}}{2}\left|\nabla_x{\rho}(t,x)\right|
	+{t}^\frac{2-{\alpha}}{2}\left|\partial_t{\rho}(t,x)\right|
	+{t}^\frac{2-{\alpha}}{2}\left|\nabla^2_x{\rho}(t,x)\right|
	+{t}^\frac{3-{\alpha}}{2}\left|\partial_t\nabla_x{\rho}(t,x)\right|\le C_{\rho}; 
	\end{equation}
and for any $t\ge1$, we have 
	\begin{equation}\label{porous'}
	\left\|\nabla_x{\rho}(t,\cdot)\right\|_{L^\infty(\mathbb{T}^d)}
	+\left\|\partial_t{\rho}(t,\cdot)\right\|_{L^\infty(\mathbb{T}^d)}
	+\left\|\nabla^2_x{\rho}(t,\cdot)\right\|_{L^\infty(\mathbb{T}^d)}
	+\left\|\partial_t\nabla_x{\rho}(t,\cdot)\right\|_{L^\infty(\mathbb{T}^d)}\lesssim 1. 
	\end{equation}
\end{lemma}

We measure the distance between solutions to the scaled nonlinear kinetic model~\eqref{epfnk} and solutions to the fast diffusion equation~\eqref{limiteq} by the relative phi-entropy functional $\mathcal{H}_{\beta}$ (see Definition~\ref{phi-entropy}). 
The following lemma shows the effectiveness of the relative phi-entropy for measuring $L^2$-distance, by virtue of the uniform convexity of $\varphi_{\beta}$. 
It can be seen as a simple version of the Csiszár-Kullback inequality on the relative entropy. 
We give its statement below with a proof taken from \cite[Proposition 2.1]{DL} for the sake of completeness. 
\begin{lemma}\label{entropydist}
Let $h_1$ and $h_2$ be two functions valued in $[0,\Lambda]$. Then, we have
\begin{equation}\label{entropydist1}
		\mathcal{H}_{\beta}(h_1|h_2)\ge \left(1-\frac{{\beta}}{2}\right)\Lambda^{-{\beta}} \|h_1-h_2\|^2_{L^2\left(\mathop{}\!\mathrm{d} m\right)}. 
\end{equation}
If we additionally assume the lower bound that $ h_1,h_2\ge\lambda$, then 
\begin{equation*}\label{entropydist2}
		\mathcal{H}_{\beta}(h_1|h_2)\le \left(1-\frac{{\beta}}{2}\right)\lambda^{-{\beta}} \|h_1-h_2\|^2_{L^2\left(\mathop{}\!\mathrm{d} m\right)}. 
\end{equation*}
\end{lemma}
	
\begin{proof}
Since $\varphi_{\beta}(1)=\varphi_{\beta}'(1)=0$ and ${\beta}\in[0,1]$, for any $z\in\mathbb{R}_+$, there exists $\xi_z\in\mathbb{R}_+$ lying between $1$ and $z$ so that 
\begin{equation*}
\varphi_{\beta}(z)=\frac{1}{2}\varphi''_{\beta}(\xi_z)(z-1)^2=\frac{2-{\beta}}{2}\xi_z^{-{\beta}}(z-1)^2. 
\end{equation*}
Since $\min\{z,1\}\le\xi_z\le\max\{z,1\}$, we have 
\begin{align*}
\int_{\mathbb{T}^d\times\mathbb{R}^d} \max\left\{h_1^{-{\beta}},h_2^{-{\beta}}\right\}|h_1-h_2|^2\mathop{}\!\mathrm{d} m \le&\; \frac{2}{2-{\beta}}\mathcal{H}_{\beta}(h_1|h_2)\\
\le&\; \int_{\mathbb{T}^d\times\mathbb{R}^d} \min\left\{h_1^{-{\beta}},h_2^{-{\beta}}\right\}|h_1-h_2|^2\mathop{}\!\mathrm{d} m,
\end{align*}
which implies the desired results by using the boundedness of $h_1,h_2$. 
\end{proof}



Let us consider the finite time diffusion asymptotics. 
\begin{proposition}\label{finitetime}
Let $\rho_{\rm in}\in {\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)$ valued in $[\lambda,\Lambda]$ with ${\alpha}_0\in(0,1)$, and the sequence of functions $\left\{h_{{\epsilon},\rm in}\right\}_{{\epsilon}\in(0,1)}\subset {\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)$ satisfying 
		\begin{equation*}
		\langleh_{{\epsilon},\rm in}\rangle\ge\lambda{\ \ \rm in\ }\mathbb{T}^d {\quad\rm and\quad} 
		0\leh_{{\epsilon},\rm in} \le \Lambda {\ \ \rm in\ }\mathbb{T}^d\times\mathbb{R}^d. 
		\end{equation*}
Let $h_{\epsilon}$ be the solutions of \eqref{epnk} associated with these initial data. 
Then, there exist some universal constants ${\alpha}\in(0,1)$, $C>0$, and some constant $C_\rho>0$ depending only on universal constants and $\left\|\rho_{\rm in}\right\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}$, such that for any ${\epsilon}\in(0,1)$ and for any $t\in[\underline{T},1]$ with $\underline{T}\in(0,1)$, the following estimate holds, 
\begin{align}\label{finitelimit}
\mathcal{H}_{\beta}(h_{\epsilon}|\rho)(t)
\le C_{\rho} \mathcal{H}_{\beta}(h_{\epsilon}|\rho)(\underline{T})+C_{\rho} {\epsilon} \left({t}^{\frac{{\alpha}-1}{2}} +{\epsilon}{t}^{\frac{{\alpha}-2}{2}}\right), 
\end{align}
where $\rho(t,x)$ is the solution to \eqref{limiteq} associated with the initial data~$\rho_{\rm in}$; 
and for any $t\ge1$, we have 
\begin{align}\label{finitelimit'}
\mathcal{H}_{\beta}(h_{\epsilon}|\rho)(t)
\le \left[\mathcal{H}_{\beta}(h_{\epsilon}|\rho)(t)(1)+ C{\epsilon}\left(1+t^\frac{1}{2}\right)\right]e^{C t}. 
\end{align}
\end{proposition}
	
	
\begin{proof}
For ${\beta}\in[0,1)$, the phi-entropy of $h_{\epsilon}$ relative to $\rho$ reads 
\begin{align*}
\mathcal{H}_{\beta}(h_{\epsilon}|\rho)=\mathcal{H}_{\beta}(h_{\epsilon}|1)-\mathcal{H}_{\beta}({\rho}|1)
-\frac{2-{\beta}}{1-{\beta}}\left(\langle h_{\epsilon}\rangle-{\rho},\;{\rho}^{1-{\beta}}-1\right).
\end{align*}
As far as the entropy $\mathcal{H}_{\beta}(h_{\epsilon}|1)$ is concerned, the entropy dissipation is derived by the equation~\eqref{epnk}, integration by parts and using Hölder's inequality $\langle\nabla_vh_{\epsilon}\rangle^2 \le \langle h_{\epsilon}\rangle^{{\beta}}\langle h_{\epsilon}^{-{\beta}}|\nabla_vh_{\epsilon}|^2\rangle$, 
\begin{align}\label{entropydissipation}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathcal{H}_{\beta}(h_{\epsilon}|1)
&=\frac{2-{\beta}}{1-{\beta}}\left(h^{1-{\beta}},\;h_t\right)
=-\frac{2-{\beta}}{{\epsilon}^2} \left(h_{\epsilon}^{-{\beta}}\nabla_vh_{\epsilon},\; \langle h_{\epsilon}\rangle^{\beta}\nabla_vh_{\epsilon}\right)\nonumber\\
&\le-\frac{2-{\beta}}{{\epsilon}^2}\|\langle\nabla_vh_{\epsilon}\rangle\|^2_{L^2_x}=-\frac{2-{\beta}}{{\epsilon}^2}\|\langle vh_{\epsilon}\rangle\|^2_{L^2_x}.
\end{align}
In view of the limiting equation~\eqref{limiteq}, we have 
\begin{align}\label{entropy2}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathcal{H}_{\beta}({\rho}|1)
=\frac{2-{\beta}}{1-{\beta}}\left({\rho}^{1-{\beta}},\;\partial_t{\rho}\right)
=-(2-{\beta})\left({\rho}^{-{\beta}}\nabla_x{\rho},\;{\rho}^{-{\beta}}\nabla_x{\rho}\right).
\end{align}
A direct computation with the macroscopic equation~\eqref{hydroFP1} and the equation~\eqref{limiteq} leads to
\begin{align}\label{entropy3}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}&\left(\langle h_{\epsilon}\rangle-{\rho},\;{\rho}^{1-{\beta}}-1\right)
=\left({\rho}^{1-{\beta}},\;\partial_t\langle h_{\epsilon}\rangle\right)
 +\left((1-{\beta}){\rho}^{-{\beta}}\langle h_{\epsilon}\rangle-(2-{\beta}){\rho}^{1-{\beta}},\;\partial_t{\rho}\right)\nonumber\\
&=\frac{1-{\beta}}{{\epsilon}}\left({\rho}^{-{\beta}}\nabla_x{\rho},\;\langle vh_{\epsilon}\rangle\right)
 -\left((1-{\beta})\nabla_x\big({\rho}^{-{\beta}}\langle h_{\epsilon}\rangle\big)-(2-{\beta}){\rho}^{-{\beta}}\nabla_x{\rho},\;{\rho}^{-{\beta}}\nabla_x{\rho}\right)
\end{align}
The evolution of $\mathcal{H}_{\beta}(h_{\epsilon}|\rho)$ is then estimated by combining \eqref{entropydissipation}, \eqref{entropy2} and \eqref{entropy3}, 
\begin{align*}
\frac{1}{2-{\beta}}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathcal{H}_{\beta}(h_{\epsilon}|\rho)
\le& -\frac{1}{{\epsilon}^2}\|\langle vh_{\epsilon}\rangle\|^2_{L^2_x}
  -\frac{1}{{\epsilon}} \left(\langle vh_{\epsilon}\rangle,\; {\rho}^{-{\beta}}\nabla_x\rho\right)
+\left(\nabla_x\big({\rho}^{-{\beta}}\langle h_{\epsilon}\rangle-{\rho}^{1-{\beta}}\big),\;{\rho}^{-{\beta}}\nabla_x\rho\right)\\
=&-\left\|{\epsilon}^{-1}\langle vh_{\epsilon}\rangle+\rho^{-{\beta}}\nabla_x{\rho}\right\|^2_{L^2_x}
   +\left({\epsilon}^{-1}\langle vh_{\epsilon}\rangle,\;\rho^{-{\beta}}\nabla_x{\rho}\right)\\
&+\left(\nabla_x\langle h_{\epsilon}\rangle+ 
{\beta}\big(1-{\rho}^{-1}\langle h_{\epsilon}\rangle\big)\nabla_x{\rho},\;\rho^{-2{\beta}}\nabla_x\rho\right)
\end{align*}
We remark that the above inequality also holds for ${\beta}=1$ by a similar computation. 
Abbreviate $Q_{\epsilon}:={\epsilon}^{-1}\langle vh_{\epsilon}\rangle+\rho^{-{\beta}}\nabla_x\rho$, $R_{\epsilon}:=-\nabla_x\cdot\langle v\otimes\nabla_vh_{\epsilon}\rangle-{\epsilon}\partial_t\langle vh_{\epsilon}\rangle$ and write the macroscopic equation \eqref{hydroFP2} in the form of
$\nabla_x\langle h_{\epsilon}\rangle = -{\epsilon}^{-1}\langle h_{\epsilon}\rangle^{{\beta}}\langle vh_{\epsilon}\rangle + R_{\epsilon}$. It then turns out that 
\begin{align}\label{retiveen0}
\frac{1}{2-{\beta}}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathcal{H}_{\beta}(h_{\epsilon}|\rho)
\le& -\|Q_{\epsilon}\|^2_{L^2_{x}} 
+\left(\big(1-{\rho}^{-{\beta}}\langle h_{\epsilon}\rangle^{{\beta}}\big)Q_{\epsilon},\;\rho^{-{\beta}}\nabla_x\rho\right)
+\left(R_{\epsilon},\;\rho^{-2{\beta}}\nabla_x\rho\right)  \nonumber\\
&+\left({\rho}^{-{\beta}}\langle h_{\epsilon}\rangle^{{\beta}}-{\beta}{\rho}^{-1}\langle h_{\epsilon}\rangle -1+{\beta} 
,\;\rho^{-2{\beta}}\left|\nabla_x{\rho}\right|^2\right)  \nonumber\\
\le&\;  2\left\|\rho^{-1-{\beta}}|\nabla_x\rho|\right\|_{L^\infty_{t,x}}^2 \|\langle h_{\epsilon}\rangle-\rho\|^2_{L^2_x}
+\left(R_{\epsilon},\;\rho^{-2{\beta}}\nabla_x\rho\right),  
\end{align}
where for the second inequality, we used the Cauchy–Schwarz inequality 
\begin{equation*}
2\left(\big(1-{\rho}^{-{\beta}}\langle h_{\epsilon}\rangle^{{\beta}}\big)Q_{\epsilon},\;\rho^{-{\beta}}\nabla_x\rho\right)
\le \|Q_{\epsilon}\|^2_{L^2_{x}} + \left(\big|1-{\rho}^{-{\beta}}\langle h_{\epsilon}\rangle^{{\beta}}\big|^2,\;\big|\rho^{-{\beta}}\nabla_x\rho\big|^2\right),
\end{equation*}
and the following two elementary inequalities, with ${\beta}\in[0,1]$, 
\begin{equation*}
		\big|z^{{\beta}}-1\big|\le|z-1| {\quad\rm and\quad} 
		\big|z^{{\beta}}-{\beta} z-1+{\beta}\big|\le|z-1|^2, {\quad\rm for\ any\ } z\in\mathbb{R}_+. 
\end{equation*}
In view of Hölder's inequality and the inequality~\eqref{entropydist1} given in Lemma~\ref{entropydist}, we know that
\begin{align*}
\|\langle h_{\epsilon}\rangle-\rho\|^2_{L^2_x}\le \|h_{\epsilon}-\rho\|^2_{L^2(\mathop{}\!\mathrm{d} m)}
\le 2\Lambda^{{\beta}}\mathcal{H}_{\beta}(h_{\epsilon}|\rho). 
\end{align*}
Combining this with \eqref{retiveen0} and \eqref{porous}, we derive that, for any $t\in(0,1]$, 
\begin{align}\label{retiveen1}
\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}\mathcal{H}_{\beta}(h_{\epsilon}|\rho)
\le C_{\rho}{t}^{{\alpha}-1} \mathcal{H}_{\beta}(h_{\epsilon}|\rho)
	+2\left(R_{\epsilon},\;\rho^{-2{\beta}}\nabla_x\rho\right), 
\end{align}
where the constants ${\alpha}\in(0,1)$ and $C_{\rho}>0$ are provided by Lemma~\ref{porouslemma}. 

We point out that, after integrating in time, the remainder term in \eqref{retiveen1} involving $R_{\epsilon}$ is of order $O({\epsilon})$ due to the control of the entropy production and the regularity of the limiting equation. Indeed, 
\begin{align*}
\int_0^t\left(R_{\epsilon},\;\rho^{-2{\beta}}\nabla_x\rho\right)
=& \int_0^t\left(\langle v\otimes\nabla_vh_{\epsilon}\rangle,\;\nabla_x\big(\rho^{-2{\beta}}\nabla_x\rho\big)\right)
  +{\epsilon}\int_0^t\left(\langle vh_{\epsilon}\rangle,\;\partial_t\big(\rho^{-2{\beta}}\nabla_x\rho\big)\right) \\
&-{\epsilon}\left(\langle vh_{\epsilon}\rangle,\;\rho^{-2{\beta}}\nabla_x\rho\right)(t) 
  +{\epsilon}\left(\langle vh_{\epsilon}\rangle,\;\rho^{-2{\beta}}\nabla_x\rho\right)(0)\\
\lesssim&\; \left(\left\|\nabla_x\big(\rho^{-2{\beta}}\nabla_x\rho\big)\right\|_{L^\infty_{t,x}} +{\epsilon}\left\|\partial_t\big(\rho^{-2{\beta}}\nabla_x\rho\big)\right\|_{L^\infty_{t,x}}\right) 
\int_0^t\left\|\langle vh_{\epsilon}\rangle\right\|_{L^2_x}\\
&+{\epsilon}\left\|\rho^{-2{\beta}}\nabla_x\rho\right\|_{L^\infty_{t,x}}\|\langle vh_{\epsilon}\rangle\|_{L^\infty\left([0,T];L^2_x\right)}. 
\end{align*}
It then follows from \eqref{porous}, \eqref{entropydissipation} and the global upper bound of $h_{\epsilon}$ (Lemma~\ref{Gaussian}) that, for any $t\in(0,1]$, 
\begin{align*}
\int_0^t\left(R_{\epsilon},\rho^{-2{\beta}}\nabla_x\rho\right)
\le&\; C_{\rho}\left({t}^{\frac{{\alpha}-1}{2}}+{\epsilon}{t}^{\frac{{\alpha}-2}{2}}\right)
\left(\int_0^t\left\|\langle vh_{\epsilon}\rangle\right\|_{L^2_{x}}^2\right)^\frac{1}{2} 
+C_{\rho}{\epsilon}{t}^{\frac{{\alpha}-1}{2}} \nonumber\\
\le&\; C_{\rho}\left({\epsilon}{t}^{\frac{{\alpha}-1}{2}} +{\epsilon}^2{t}^{\frac{{\alpha}-2}{2}}\right) \sup\nolimits_{s\in[0,t]}\sqrt{\mathcal{H}_{\beta}(h_{\epsilon}|1)(s)} 
+C_{\rho}{\epsilon}{t}^{\frac{{\alpha}-1}{2}} \nonumber\\
\le&\; C_{\rho}\left({\epsilon}{t}^{\frac{{\alpha}-1}{2}} +{\epsilon}^2{t}^{\frac{{\alpha}-2}{2}}\right)
+C_{\rho}{\epsilon}{t}^{\frac{{\alpha}-1}{2}}. 
\end{align*}	
Combining this estimate with \eqref{retiveen1}, as well as Grönwall's inequality, we conclude \eqref{finitelimit}. 
Besides, we arrive at \eqref{finitelimit'}, if we apply Lemma~\ref{porouslemma} with \eqref{porous'} instead of \eqref{porous} in above argument. This completes the proof. 
\end{proof}
	
	
	
	
We are now in a position to conclude the global in time diffusion asymptotics. 
\begin{proof}[Proof of Theorem \ref{limit}]
We are going to combine Proposition~\ref{longtime}, \ref{finitetime} with a delicate analysis on the relative entropy around the initial time to get Theorem~\ref{limit}. 
The analysis is based on the barrier function method. 
Let us assume the constant ${\alpha}\in(0,1)$ provided by Proposition~\ref{finitetime}. 

\medskip\noindent{\textbf{Step 1.}}	Pointwise estimate. \\
Let us fix ${\delta}\in(0,1)$, $(x_1,v_1)\in\mathbb{T}^d\times B_{R}$ with $R>0$, and consider the function 
\begin{equation*}
\overline{h}(t,x,v):=C_1t+C_2\left(|x-x_1-{\epsilon}^{-1}tv|^2+|v-v_1|^2\right)
\end{equation*}
where the constants $C_1,C_2>0$ are to be determined.  
For any $t\le\frac{{\epsilon}{\delta}}{4(1+R)}$, we have
\begin{equation*}
\overline{h} \ge C_2\left(|x-x_1|^2+|v-v_1|^2-2{\epsilon}^{-1}t|x-x_1||v|\right)\ge\frac{C_2{\delta}^2}{2}=\Lambda
{\quad\rm on\ } \partial B_{\delta}(x_1,v_1), 
\end{equation*} 
where we chose $C_2:=2{\delta}^{-2}\Lambda$. 
For any $(x,v)\in B_{\delta}(x_1,v_1)$, 
\begin{equation*}
		|\langle h\rangle^{\beta} \mathscr{L}_{\scriptscriptstyle{\rm OU}}\overline{h}|
		\lesssim |\Delta_v\overline{h}|+ |v\cdot\nabla_v\overline{h}|
		\lesssim {\delta}^{-2}(1+R^2)\left(1+{\epsilon}^{-2}t^2\right). 
\end{equation*}
Therefore, for any $t\le\frac{{\epsilon}{\delta}}{4(1+R)}$ and $(x,v)\in B_{\delta}(x_1,v_1)$, 
\begin{equation*}
\left(\partial_t+{\epsilon}^{-1}v\cdot\nabla_x-{\epsilon}^{-2}\langle h\rangle^{\beta}\mathscr{L}_{\scriptscriptstyle{\rm OU}}\right)\overline{h}
\ge C_1-C_0{\epsilon}^{-2}{\delta}^{-2}(1+R^2)\left(1+{\epsilon}^{-2}t^2\right) \ge 0, 
\end{equation*}
where the constant $C_0>0$ is universal and we chose $C_1:=2C_0{\epsilon}^{-2}{\delta}^{-2}(1+R^2)$. 
Then, the maximum principle implies that, for any $t\le\frac{{\epsilon}{\delta}}{4(1+R)}$ and $(x,v)\in B_{\delta}(x_1,v_1)$, 
\begin{equation}\label{coninitialbarrier}
|h_{\epsilon}(t,x,v)-h_{{\epsilon},\rm in}(x_1,v_1)|
\le \overline{h}(t,x,v)+\sup\nolimits_{B_{\delta}(x_1,v_1)}|h_{{\epsilon},\rm in}(x,v)-h_{{\epsilon},\rm in}(x_1,v_1)|. 
\end{equation}
In particular, for any $t\le\frac{{\epsilon}{\delta}}{4(1+R)}$ and $(x_1,v_1)\in\mathbb{T}^d\times B_{R}$, 
\begin{equation}\label{limitbarrier}
|h_{\epsilon}(t,x_1,v_1)-h_{{\epsilon},\rm in}(x_1,v_1)|
\lesssim {\epsilon}^{-2}{\delta}^{-2}(1+R^2)t +\|h_{{\epsilon},\rm in}\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)}{\delta}^{\alpha}. 
\end{equation}

As far as the solution ${\rho}$ to the limiting equation~\eqref{limiteq} is concerned, using the H\"older estimate~\eqref{porous0} in Lemma~\ref{porouslemma}, we  derive that, for any $t\in\mathbb{R}_+$, 
\begin{align}\label{limitholder}
\|{\rho}(t)-{\rho}_{\rm in}\|_{L^\infty(\mathbb{T}^d)}
\lesssim \left(1+\|{\rho}_{\rm in}\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}\right)t^{\alpha}. 
\end{align}


\medskip\noindent{\textbf{Step 2.}}	Estimate of the relative entropy around the initial time. \\
Let us restrict our attention on the time interval $[0,\underline{T}]$ with $\underline{T}\in(0,1)$ to be determined. 
Compute the relative entropy $\mathcal{H}_{\beta}(h_{\epsilon}|\rho)$ in terms of initial data as follows, 
\begin{equation}\label{limitinitial0}
\begin{split}
\mathcal{H}_{\beta}(h_{\epsilon}|\rho)
  =\, \mathcal{H}_{\beta}(h_{{\epsilon},\rm in}|\rho_{\rm in})
   +\int_{\mathbb{T}^d\times\mathbb{R}^d}\left[\left(\varphi_{\beta}(h_{\epsilon})-\varphi_{\beta}(h_{{\epsilon},\rm in})\right)
   -\left(\varphi_{\beta}({\rho})-\varphi_{\beta}({\rho}_{\rm in})\right)\right]\mathop{}\!\mathrm{d} m&\\
-\int_{\mathbb{T}^d\times\mathbb{R}^d}\left[\left(\varphi'_{\beta}({\rho})-\varphi'_{\beta}({\rho}_{\rm in})\right)(h_{{\epsilon},\rm in}-{\rho}_{\rm in}) 
 +\varphi'_{\beta}({\rho})(h_{\epsilon}-h_{{\epsilon},\rm in}) -\varphi'_{\beta}({\rho})({\rho}-{\rho}_{\rm in})
 \right]\mathop{}\!\mathrm{d} m&.
\end{split}
\end{equation}
Consider a truncation in $v$ for the integrals on the right hand side. 
Since $h_{\epsilon},h_{{\epsilon},\rm in},{\rho},{\rho}_{\rm in}$ are all bounded from above (Lemma~\ref{Gaussian}), we have 
\begin{equation}\label{limitbarrierc}
\int_{\mathbb{T}^d\times B^c_R}\left[|\varphi_{\beta}(h_{\epsilon})-\varphi_{\beta}(h_{{\epsilon},\rm in})|
+|\varphi_{\beta}({\rho})-\varphi_{\beta}({\rho}_{\rm in})|\right]\mathop{}\!\mathrm{d} m
\lesssim \int_{B^c_{R}}\mathop{}\!\mathrm{d}{\mu}\lesssim R^{-4}. 
\end{equation}
Observe that for any $a,b\in(0,\Lambda]$, there exists $\xi\in\mathbb{R}_+$ lying between $a$ and $b$ such that 
$\varphi_{\beta}(a^2)-\varphi_{\beta}(b^2)=2\xi\varphi'_{\beta}(\xi^2)(a-b)$. 
Meanwhile, for any $\xi\in(0,\Lambda]$, $|\xi\varphi'_{\beta}(\xi^2)|\lesssim1$. 
Thus, 
\begin{equation*}
|\varphi_{\beta}(h_{\epsilon})-\varphi_{\beta}(h_{{\epsilon},\rm in})|\lesssim |h_{\epsilon}-h_{{\epsilon},\rm in}|^\frac{1}{2}
{\quad\rm and\quad}
|\varphi_{\beta}({\rho})-\varphi_{\beta}({\rho}_{\rm in})|\lesssim |{\rho}-{\rho}_{\rm in}|^\frac{1}{2}. 
\end{equation*}
Set $R:={\epsilon}^{-\frac{\eta}{4}}$, ${\delta}:={\epsilon}^{\frac{\eta}{4}}$ and $\underline{T}:=\frac{1}{8}{\epsilon}^{2+2\eta}$ for some constant $\eta\in(0,1)$ to be determined. 
In this setting, $\underline{T}\le\frac{{\epsilon}{\delta}}{4\langle R\rangle}$. 
It then follows from \eqref{limitbarrier}, \eqref{limitholder} and \eqref{limitbarrierc} that, for any $t\le\underline{T}$, 
\begin{align}\label{limitinitial1}
\int_{\mathbb{T}^d\times \mathbb{R}^d}&\left[\left|\varphi_{\beta}(h_{\epsilon})-\varphi_{\beta}(h_{{\epsilon},\rm in})\right|
+\left|\varphi_{\beta}({\rho})-\varphi_{\beta}({\rho}_{\rm in})\right|\right]\mathop{}\!\mathrm{d} m\nonumber\\
&\ \lesssim R^{-4} +\int_{\mathbb{T}^d\times B_R}\left[|h_{\epsilon}-h_{{\epsilon},\rm in}|^\frac{1}{2} +|{\rho}-{\rho}_{\rm in}|^\frac{1}{2}\right]\mathop{}\!\mathrm{d} m\nonumber\\
&\ \lesssim {\epsilon}^{\eta} + {\epsilon}^{\frac{\eta}{2}}+\|h_{{\epsilon},\rm in}\|^\frac{1}{2}_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)}{\epsilon}^{\frac{{\alpha}\eta}{8}}
+\Big(1+\|{\rho}_{\rm in}\|^\frac{1}{2}_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}\Big){\epsilon}^{{\alpha}+{\alpha}\eta}.  
\end{align}
Besides, $\|h_{\epsilon}-h_{{\epsilon},\rm in}\|_{L^1(\mathbb{T}^d\times B^c_R,\,\mathop{}\!\mathrm{d} m)}\lesssim R^{-4}$ and $|\varphi'_{\beta}|\lesssim1$ on $[\lambda,\Lambda]$. 
Therefore, combining \eqref{limitbarrier}, \eqref{limitholder} with the inequality~\eqref{entropydist1} given in Lemma~\ref{entropydist} yields that, for any $t\le\underline{T}$, 
\begin{align}\label{limitinitial2}
\int_{\mathbb{T}^d\times\mathbb{R}^d}&\left[\left|\left(\varphi'_{\beta}({\rho})-\varphi'_{\beta}({\rho}_{\rm in})\right)(h_{{\epsilon},\rm in}-{\rho}_{\rm in})\right| 
 +\left|\varphi'_{\beta}({\rho})(h_{\epsilon}-h_{{\epsilon},\rm in})\right| +\left|\varphi'_{\beta}({\rho})({\rho}-{\rho}_{\rm in})\right|
 \right]\mathop{}\!\mathrm{d} m\nonumber\\
&\ \lesssim  \|h_{{\epsilon},\rm in}-{\rho}_{\rm in}\|_{L^2\left(\mathbb{T}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m\right)} 
+\|h_{\epsilon}-h_{{\epsilon},\rm in}\|_{L^1\left(\mathbb{T}^d\times\mathbb{R}^d,\,\mathop{}\!\mathrm{d} m\right)} +\|{\rho}-{\rho}_{\rm in}\|_{L^1(\mathbb{T}^d)}\nonumber\\
&\ \lesssim  \mathcal{H}_{\beta}^\frac{1}{2}(h_{{\epsilon},\rm in}|{\rho}_{\rm in})+ R^{-4} +\|h_{\epsilon}-h_{{\epsilon},\rm in}\|_{L^1\left(\mathbb{T}^d\times B_R,\,\mathop{}\!\mathrm{d} m\right)} +\|{\rho}-{\rho}_{\rm in}\|_{L^1(\mathbb{T}^d)}\nonumber\\
&\ \lesssim {{\epsilon}'}^\frac{1}{2}
+{\epsilon}^{\eta}+\|h_{{\epsilon},\rm in}\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)}{\epsilon}^{\frac{{\alpha}\eta}{4}}
+\left(1+\|{\rho}_{\rm in}\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}\right){\epsilon}^{2{\alpha}+2{\alpha}\eta}. 
\end{align}
Plugging \eqref{limitinitial1} and \eqref{limitinitial2} into the expression \eqref{limitinitial0}, we derive that, for any $t\le\underline{T}$, 
\begin{equation}\label{limitinitial}
\mathcal{H}_{\beta}(h_{\epsilon}|\rho)(t)
\lesssim  
\left(1+\|h_{{\epsilon},\rm in}\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)}+\|{\rho}_{\rm in}\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}\right) ({\epsilon}+{\epsilon}')^{\frac{{\alpha}\eta}{8}}  .  
\end{equation}


\medskip\noindent{\textbf{Step 3.}}	Conclusion. \\
Recall that we have chosen $\underline{T}=\frac{1}{8}{\epsilon}^{2+2\eta}$. 
In view of \eqref{limitinitial} and the estimate~\eqref{finitelimit} given in Lemma~\ref{finitetime}, one may optimize in $\eta$ to get the result. 
For simplicity, we pick $\eta:=\frac{{\alpha}}{4-2{\alpha}}$ so that $\underline{T}=\frac{1}{8}{\epsilon}^{\frac{4-{\alpha}}{2-{\alpha}}}$ 
and ${\epsilon}\Big(\underline{T}^{\frac{{\alpha}-1}{2}}+{\epsilon}\underline{T}^{\frac{{\alpha}-2}{2}}\Big)\lesssim{\epsilon}^\frac{{\alpha}}{2}$. 
It turns out that for any $t\in[0,1]$, 
\begin{align}\label{conclu}
\mathcal{H}_{\beta}(h_{\epsilon}|\rho)(t)
\le C_{\rho}\mathcal{H}_{\beta}(h_{\epsilon}|\rho)(\underline{T})
+C_{\rho}{\epsilon}\left(\underline{T}^{\frac{{\alpha}-1}{2}}+{\epsilon}\underline{T}^{\frac{{\alpha}-2}{2}}\right)
\le C_* ({\epsilon}+{\epsilon}')^{\frac{{\alpha}\eta}{8}},  
\end{align}
where the constant $C_{\rho}>0$ is provided in Lemma~\ref{finitetime} and the constant $C_*>0$ depends only on universal constants, $\left\|\rho_{\rm in}\right\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d)}$ and $\left\|h_{{\epsilon},\rm in}\right\|_{{\mathcal C}^{{\alpha}_0}(\mathbb{T}^d\times\mathbb{R}^d)}$.
Then, using the estimate \eqref{finitelimit'} given in Lemma~\ref{finitetime} with \eqref{entropydist1}, we arrive at point (i) of Theorem~\ref{limit}. 


As for point (ii) of Theorem~\ref{limit}, applying \eqref{finitelimit'}, together with \eqref{entropydist1} and \eqref{conclu}, for any $t\in[1,\overline{T}]$, we have 
\begin{align*}
\|h_{\epsilon}(t)-\rho(t)\|_{L^2(\mathop{}\!\mathrm{d} m)}^2
\lesssim& \left[\mathcal{H}_{\beta}(h_{\epsilon}|\rho)(1)+ {\epsilon}\Big(1+\overline{T}^\frac{1}{2}\Big)\right] e^{C\overline{T}}\\
\lesssim& \left[C_*({\epsilon}+{\epsilon}')^{\frac{{\alpha}\eta}{8}} +{\epsilon}\Big(1+(-\iota\log({\epsilon}+{\epsilon}'))^\frac{1}{2}\Big)\right]
({\epsilon}+{\epsilon}')^{-\frac{{\alpha}\eta}{16}}
\lesssim C_*({\epsilon}+{\epsilon}')^{\frac{{\alpha}\eta}{16}}. 
\end{align*}
where we picked $\overline{T}:=-\iota\log({\epsilon}+{\epsilon}')$ with $\iota:=\frac{{\alpha}\eta}{16C}$. 
Finally, using Proposition~\ref{longtime} under the additional assumption that $h_{{\epsilon},\rm in}\ge\lambda$, we know the long time behavior that there is some universal constant $c>0$ such that for any $t\ge\overline{T}$, 
\begin{equation*}
\|h_{\epsilon}(t)-\rho(t)\|_{L^2(\mathop{}\!\mathrm{d} m)}
\le \|h_{\epsilon}(t)-M_0\|_{L^2(\mathop{}\!\mathrm{d} m)}+\|{\rho}(t)-M_0\|_{L_x^2}
\lesssim e^{-c\overline{T}}= ({\epsilon}+{\epsilon}')^{c\iota}.  
\end{equation*}
The proof now is complete.  
\end{proof}



