\section{Introduction}
\label{intro}

The adiabatic piston problem is an old problem in Thermodynamics and Statistical Mechanics. As described by Callen \cite{ca63} many years ago, this is to predict the equilibrium state of a cylinder containing two subsystems (usually two gases) separated by a movable adiabatic piston, i.e. a wall that isolates the two subsystems when kept fixed. The same author, and others \cite{cule79}, stated that Thermodynamics can only provide a partial solution to the problem, namely the equality of pressures.

Lieb \cite{li99} went one step further by pointing out that a Statistical Mechanics solution to the problem is in contradiction with Thermodynamics. Namely, as described by Landau and Lifshitz \cite{lali13} and by Feynmann \cite{fe11}, the relaxation of the system, through the motion of the piston, is towards a global equilibrium, where all subsystems including the piston have the same temperature and the pressures of the gasses are equal. But this implies that the initially hotter subsystem decreases its entropy by doing work to the other subsystem, which may be in contradiction with the Second Principle \cite{glpr71,chlesi02a,grma13}. According to Lieb and Yngvason \cite{li99,liyn98,liyn99} the adiabatic piston problem can be solved by reconsidering the principle of maximizing the entropy in the absence of constrains.

Most works on the adiabatic piston use simple models, such as non-interacting particles. In this context, rigorous results can be found in the literature. If the initial condition is of a fixed piston surrounded by two ideal gases in equilibrium with the same pressures and temperatures, the dynamics of the piston occurs in several stages \cite{chle02}: the initial condition becomes unstable after a short time and the piston start oscillating \cite{si99,nesi04,wr07}, the amplitude of the oscillations attenuates exponentially, and a final thermalization leads to the equilibration of the whole system. For more general initial conditions, where the two gases may be in different states, the evolution is in two stages \cite{chlesi02a,grpale03,cachle04}: a first fast evolution towards mechanical equilibrium, where the two gases have the same pressure, is followed by a much slower relaxation towards thermal equilibrium, dominated by the energy flux through the piston. Moreover, if the effect of recollisions (when a particle collide with the piston more than once) can be removed under some limits, the motion of the piston converges to an Ornstein-Uhlenbeck process \cite{chlesi02a,ho71,dugole81}. See \cite{chlesi02a} for an attempt to include recollisions.

The last stage of the dynamics, under the condition of mechanical equilibrium but with the temperatures of the gases still different, has attained much attention. As a matter of fact, the existence of this stage is closely related to the original controversy of the adiabatic piston problem. An interesting model, that keeps the system always in this last stage, considers two semiinfinite ideal gases in each sides of the piston having the same pressure and temperature difference. Even the fact that the net force to the piston is zero, the system reaches a steady state where the mean velocity of the piston is constant, different from zero, and towards the hotter gas \cite{grpi99,pigr99,grfr99,pi01,ch04,itsa15}. When the mass of the piston goes to infinity, its mean velocity vanishes \cite{pi02,grpa02,grpale02}. See also \cite{grpale04} for a kinetic description of the model in close relation with Thermodynamics.

From a mesoscopic and macroscopic scales perspective, the major challenge when dealing with the adiabatic piston is to derive  closed equations for the piston and the gases. This has been done rigorously \cite{lepisi00,chlesi02} and approximately \cite{crdise96,gr99,cepapivu07,gi10,mavake05,magaba06} for the case of non-interacting and similar models. For more realistic ones, such as models with hard spheres or disks, most of the existing results are numerical \cite{hure06,whgove02}, or use rude approximations \cite{brkh12}. As in the case of non-interacting particles, the numerical simulations show an exponential relaxation of the system towards thermal equilibrium \cite{kevama00,fosh13}.

So far, the numerous theoretical results have not been compared against experiments. This is partly because the last stage of the evolution, actually the interesting one for the adiabatic piston problem, is too much slow to be observed in macroscopic system. As far as we know, there is only one experimental work that can be relevant for the present discussion \cite{leda10}. It provides experimental measurements on the thermalization of a piston surrounded by a few hundreds of grains  in a very dense configuration, close to the jamming transition. What we learn from this experiment, and also from the theory done in the dilute regime \cite{brkh11}, is that the equilibration in the inelastic or granular case is completely different from the classic or elastic one. Specifically, when there is no energy source supplying the dissipation of the inelastic collisions, each of the gases reaches a state close to the homogeneous cooling state \cite{kh18}, and the ``thermal'' equilibration occurs when the cooling rates (or rates of energy loss) of both gases and the piston are the same.  The new ``equilibrium'' criterion, as opposite of the temperatures being the same for the elastic case, explains the emergence of non-equilibrium phase transitions \cite{brrebr05,hure06a,brkh10}. See also \cite{brkh11,brru04,brru08,brru09,brru09,brru10} for theoretical and numerical studies relevant for experimental situations.

The works on the adiabatic piston in contact with granular matter suggest the need to include at least two new ingredients in the study of the (elastic) adiabatic piston. Firstly, if we are interested in making some contact with reality, a systematic derivation of a theory for the piston in contact with a finite number of particles is needed. Secondly, the theory should include interactions between particles. The aim of the present work is to propose a mesoscopic description of a realistic model of the adiabatic piston. The model considers a large but finite number of interacting particles, and reduces to the case of non-interacting particles when the diameters of the particles tend to zero.

The organization of the work is as follow. The model is presented in Sec. \ref{sec:2} where we  also show, by computing some relevant quantities of the system at thermal equilibrium, the importance of keeping the spatial correlations between the gases and the piston. Section \ref{sec:3} contains the main results of the work, namely the derivation of a closed set of kinetic equations for the gases, conditioned to a position of the piston, and the piston. It is shown that the equations support the equilibrium solution of Sec. \ref{sec:2}. As an application to the results of Sec. \ref{sec:3}, in Sec. \ref{sec:4} we demonstrate the $\mathcal H$-theorem, which implies the evolution of the system towards thermal equilibrium from any initial condition. Finally, Sec. \ref{sec:5} includes a discussion together with some conclusions.

\section{Model and equilibrium \label{sec:2}}
\subsection{Model}
The system is a $d$-dimensional cylinder of length $L$ and section $S$ divided into two parts by a moving wall of section $S$, zero width, and mass $m_p$ (the piston). The normal direction of the piston remains parallel to the axis of the cylinder, taken as the $X$-axis. The system contains $N_i$ $d-$dimensional hard spheres with masses $m_i$ and diameter $\sigma_i$, to the left of the piston for $i=1$ and to the right of the piston for $i=2$. See figure \ref{fig:1} for a sketch of the system.

\begin{figure}[!h]
  \centering
  \includegraphics[width=.5\textwidth]{fig1.eps}
  \caption{A representation of the system for $d=2$. }
  \label{fig:1}
\end{figure}

The microstate of the system is given by the set of positions and velocities as $\{\mathbf r_1,\dots,\mathbf r_{N},x_p;\mathbf v_1,\dots,\mathbf v_{N},v_p\}\equiv \{\mathbf R,x_p;\mathbf V,v_p\}$ or the set of positions and momenta $\Gamma\equiv \{\mathbf r_1,\dots,\mathbf r_{N},x_p;\mathbf p_1,\dots,\mathbf p_{N},v_p\}\equiv \{\mathbf R,x_p;\mathbf P,p_p\}$, where $\mathbf p_i=m_i\mathbf v_i$, $\mathbf p_p=m_p\mathbf v_p$,  and $N=N_1+N_2$ is the total number of particles. In the latter expressions, the indexes $i=1,\dots,N_1$ stand for the particles to the left of the piston and $i=N_1+1,\dots,N$ for the particles to the right.

The state of the system changes because of the free motion of the particles and the piston, as well as of instantaneous collisions among particles and between particles and the piston. The collisions conserve energy and linear momentum. More specifically, for two particles colliding with velocities $\mathbf v_k$ and $\mathbf v_l$, their postcollisional velocities, denoted by a prime, are
\begin{eqnarray}
  \label{eq:coll}
    && \mathbf{v}_k'=\mathbf{v}_k-[(\mathbf{v}_k-\mathbf{v}_l)\cdot \hat{\boldsymbol{\sigma}}]\hat{\boldsymbol{\sigma}}, \\
  \label{eq:coll2}
    && \mathbf{v}_l'=\mathbf{v}_l+[(\mathbf{v}_k-\mathbf{v}_l)\cdot \hat{\boldsymbol{\sigma}}]\hat{\boldsymbol{\sigma}},
\end{eqnarray}
where $\hat{\boldsymbol{\sigma}}$ is a unit vector pointing from the center of particle $k$ to the center of particle $l$ when in contact. The collision rule for a particle with mass $m_i$ and velocity $\mathbf v$ and the piston with velocity $v_p$  is
\begin{eqnarray}
  &&v_x'= v_x-\frac{2m_p}{m_i+m_p}(v_x-v_p), \\
  &&\mathbf v_\perp'= \mathbf v_\perp, \\
  &&v_p'= v_p-\frac{2m_i}{m_i+m_p}(v_p-v_x),
\end{eqnarray}
when $v_x>v_p$ if the particle is on the left ($i=1$) and $v_x<v_p$ if the particle is on the right ($i=2$). The symbol $\perp$ represents the component of the velocity normal to the $X$-axis. An alternative characterization of the dynamics, more suitable when using the phase space representation, is provided by the hamiltonioan
\begin{equation}
  \label{eq:ham}
  \mathcal H =\left\{
    \begin{split}
      & \sum_{i=1}^{N_1}\frac{p_i^2}{2m_1}+\sum_{i=N_1+1}^{N}\frac{p_i^2}{2m_2}+\frac{p_p^2}{2m_p}, \quad (\mathbf R,x_p)\in \mathcal V, \\
      & +\infty, \qquad (\mathbf R,x_p)\notin \mathcal V,
    \end{split}
  \right.
\end{equation}
where $\mathcal V$ is the volume of the phase space accessible to the system (no overlaps).

The main focus in this work is on a mesoscopic description of the system. First, we consider the canonical equilibrium and the probability density of the phase space $\rho(\mathbf R,x_p;\mathbf P,p_p)$ and related probability densities. Later, a non-equilibrium description based on Kinetic Theory will be proposed and the distribution functions for the gases and the piston, to be defined later, will be used instead. Many of the results to be obtained at equilibrium, although of trivial derivation, will deserve as a reference and a guide for a more general out-of-equilibrium description.

\subsection{Canonical equilibrium}

If the system is big enough or if it is in contact with a thermal bath with temperature $T$, i.e. if the walls of the container vibrates in equilibrium with temperature $T$, then the probability density of the whole system is given by the canonical distribution \cite{pabe11}
\begin{equation}
  \rho=\frac{e^{-\beta \mathcal H}}{h^{dN+1}\mathcal Z},
\end{equation}
where $h$ is the Planck constant, $\beta\equiv \frac{1}{k_BT}$ with $k_B$ the Boltzmann constant, $\mathcal H$ is the hamiltonian of Eq. \eqref{eq:ham}, and $\mathcal Z$ the partition function. The latter is difficult to compute in general, due to the exclusion effects. However, for $N(\sigma_1+\sigma_2) \ll L$, that is if we neglect the volume of particles, an approximation to be taken along the work, the canonical distribution is well approximated as
\begin{eqnarray}
  \nonumber
  \mathcal Z=&&h^{-(dN+1)} \int d\mathbf R dx_p d\mathbf P dp_p  e^{-\beta \mathcal H} \\ \simeq&& \nonumber \frac{(LS)^{N}L}{\lambda_{1}^{dN_1}\lambda_{2}^{dN_2}\lambda_{p}} \int_{0}^1 ds \left[\int_0^s dx\right]^{N_1}\left[\int_s^1 dy\right]^{N_2} \\
  =&& \frac{(LS)^{N}L}{\lambda_{1}^{dN_1}\lambda_{2}^{dN_2}\lambda_{p}} \frac{N_1!N_2!}{(N+1)!}.
\end{eqnarray}
where $\lambda_{i}\equiv \sqrt{\frac{\beta h^2}{2\pi m_i}}$ is the thermal length associated to a mass $m_i$, $i=1,2,p$.

From the probability density we can obtain other relevant quantities, such as the probability density of a particle to the left $\rho_1$, to the right $\rho_2$, and of the piston $\rho_p$. For the first one, we have
\begin{eqnarray}
  \nonumber    \rho_1(\mathbf r,\mathbf p)\equiv &&\int \prod\limits_{i=2}^{N}d\mathbf r_id\mathbf p_i  dx_p dp_p  \rho(\mathbf R, x_p;\mathbf p,p_p)\\   \nonumber \simeq &&\frac{\lambda_{1}^de^{-\beta\frac{p^2}{2m_1}}}{h^dLS}\frac{(N+1)!}{N_1!N_2!}  \int_{x/L}^1 ds\left[\int_0^s dx\right]^{N_1-1}\left[\int_s^1 dy\right]^{N_2}\\ =&&\frac{(N+1)\lambda_{1}^d}{N_1h^dLS}\left[1-I\left(\frac{x}{L},N_1,N_2+1\right)\right]e^{-\beta\frac{p^2}{2m_1}},
\end{eqnarray}
where $I(x;a,b)$ is the regularized incomplete beta function. Proceeding analogously,
\begin{eqnarray}
  && \rho_2(\mathbf r,\mathbf p)\simeq \frac{(N+1)\lambda_{2}^d}{N_2h^dLS}I\left(\frac{x}{L},N_1+1,N_2\right)e^{-\beta\frac{p^2}{2m_2}}, \\
  \label{eq:rhop}
  && \rho_p(x,p)\simeq \frac{(N+1)!\lambda_{p}}{N_1!N_2!hL}\left(\frac{x}{L}\right)^{N_1}\left(1-\frac{x}{L}\right)^{N_2}e^{-\beta\frac{p^2}{2m_p}}.
\end{eqnarray}
These expressions describe two non-homogeneous gases and a fluctuating piston, as shown in Fig. \ref{fig:1} for two representative cases. Only for $N_1,N_2\to \infty$ (regardless the mass of the piston) are the two gases spatially homogeneous, with the position of the piston fixed (although with a fluctuating velocity).

Despite the previous results, the system is spatially homogeneous at equilibrium, namely the probability of finding any particle (including the piston) in a given position is uniform. This can be seen from the density of particles, defined as
\begin{eqnarray}
  && n_i(\mathbf r,t)\equiv N_i \int d\mathbf p\ \rho_i(\mathbf r,\mathbf p, t), \qquad i=1,2, \\
  \label{eq:den1}
  && n_p(x,t)\equiv \int d p\ \rho_p(x,p, t),
\end{eqnarray}
for the particles on the left ($i=1$), the particles on the right ($i=2$), and the piston. At thermal equilibrium, each of the quantities depends on $x$, however it is easily seen that the global density $n$ is spatially uniform
\begin{equation}
n(x)\equiv  n_1(x)+n_2(x)+n_p(x)=\frac{N+1}{LS}.
\end{equation}

\begin{figure}[!h]
  \centering
  \includegraphics[width=.45\textwidth]{fig2_1.eps}
  \includegraphics[width=.45\textwidth]{fig2_2.eps}
  \caption{Probability densities for particles and the piston at equilibrium for different number of particles.}
  \label{fig:2}
\end{figure}

Similar arguments allow us to show that the pressure is  is $\frac{N+1}{LS\beta}$ along the $X$-axis and $\frac{N+1}{LS\beta}$ along any other normal direction. The calculation can be made using the partition function (then we should change the volume by fixing $S$ in the first case and fixing $L$ in the second one) or by computing the net flux of linear momentum across an imaginary $(d-1)$-surface.

\subsection{Conditional probabilities}

As already shown, the fluctuations of the position of the piston make the probability densities for the gases to be non-homogenous at thermal equilibrium. Interestingly, the conditional probability densities of the gases to a position of the piston, $\rho_i(\mathbf r,\mathbf p|x_p)$ for $i=1,2$, are spatially homogenous at equilibrium. Using the definition of the conditional probabilities, and neglecting volume exclusion effects, we have
\begin{eqnarray}
 \nonumber
  \rho_1(\mathbf r,\mathbf p|x_p)\equiv &&\frac{\rho_1(\mathbf r,x_p,\mathbf p)}{n_p(x_p)} \\ = \nonumber  &&\int_{x_i\le x_p}\prod_{i=1}^{N_1-1} d\mathbf r_id\mathbf p_i\int_{x_j\ge x_p}\prod_{j=N_1+1}^{N} d\mathbf r_jd\mathbf p_jdp_p\rho(\mathbf R,x_p;\mathbf P,p_p) \\ = && \frac{\lambda_1^d}{h^dSx_p}e^{-\beta \frac{p^2}{2m_1}}, \qquad x\le x_p,
\end{eqnarray}
and
\begin{equation}
  \rho_2(\mathbf r,\mathbf p|x_p)=\frac{\lambda_2^d}{h^dS(L-x_p)}e^{-\beta \frac{p^2}{2m_2}}, \qquad x\ge x_p.
\end{equation}

There are two remarkable aspects of the previous results. On the one hand, the conditional distributions characterize the equilibrium state as if the piston where fixed. On the other hand, at thermal equilibrium the piston always ``feels'' homogeneous gases. The latter implies that we can also compute the equilibrium probability density of the piston as if it where in equilibrium under the force $F(x_p)$ exerted by the gases. Namely, at equilibrium,
\begin{equation}
  \begin{split}
    F&=\int_{\frac{p_x}{m_1}>\frac{p_p}{m_p}}d\mathbf r_\perp d\mathbf p dp_p\ \left|\frac{p_x}{m_1}-\frac{p_p}{m_p}\right| \rho_1((x_p,\mathbf r_\perp),\mathbf p|x_p)\frac{\lambda_{p}}{h}e^{-\beta\frac{p_p^2}{2m_p}}\Delta^-p\\ &  -\int_{\frac{p_x}{m_2}<\frac{p_p}{m_p}}d\mathbf r_\perp d\mathbf p dp_p\ \left|\frac{p_p}{m_p}-\frac{p_x}{m_2}\right| \rho_2((x_p,\mathbf r_\perp),\mathbf p|x_p)\frac{\lambda_{p}}{h}e^{-\beta\frac{p^2}{2m_p}}\Delta^+p,
  \end{split}
\end{equation}
where
\begin{equation}
  \label{eq:19}
  \begin{split}
    \Delta^-p=& \frac{2m_1m_p}{m_1+m_p}\left(\frac{p_x}{m_1}-\frac{p_p}{m_p}\right);  \\
    \Delta^+p=& \frac{2m_2m_p}{m_2+m_p}\left(\frac{p_p}{m_p}-\frac{p_x}{m_2}\right),
  \end{split}
\end{equation}
are the change of the linear momentum of the piston due to collisions with the gases. Operating,
\begin{eqnarray}
  \nonumber
  F=&&\frac{2m_1m_p}{m_1+m_p}\frac{N_1\lambda_1^d\lambda_{p}}{h^{d+1}x_p}\int_{\frac{p_x}{m_1}>\frac{p_p}{m_p}} d\mathbf p dp_p\ \left(\frac{p_x}{m_1}-\frac{p_p}{m_p}\right)^2 e^{-\beta\left( \frac{p^2}{2m_1}+\frac{p_p^2}{2m_p}\right)} \\ \nonumber && -\frac{2m_2m_p}{m_2+m_p}\frac{N_2\lambda_2^d\lambda_{p}}{h^{d+1}(L-x_p)}\int_{\frac{p_x}{m_2}<\frac{p_p}{m_p}} d\mathbf p dp_p\ \left(\frac{p_p}{m_p}-\frac{p_x}{m_2}\right)^2 e^{-\beta\left( \frac{p^2}{2m_2}+\frac{p_p^2}{2m_p}\right)} \\ =&&\frac{1}{\beta}\left[\frac{N_1}{x_p}-\frac{N_2}{(L-x_p)}\right].
\end{eqnarray}
The associated ``effective'' potential $U(x)$ verifies $F(x)=-\frac{dU(x)}{dx}$ and is, up to an additive constant,
\begin{equation}
  U(x)=-\frac{1}{\beta}\ln\left[x^{N_1}(L-x)^{N_2}\right].
\end{equation}
Hence, the desired probability density of the piston is
\begin{equation}
  \rho_p(x,p)=\frac{1}{\mathcal Z_p}\exp\left[-\beta\left(\frac{p^2}{2m_p}+U(x)\right)\right]=\frac{x^{N_1}(L-x)^{N_2}}{\mathcal Z_p}e^{-\beta\frac{p^2}{2m_p}},
\end{equation}
with $\mathcal Z_p$ a normalization constant. The latter expression of $\rho_p$ coincides with Eq. \eqref{eq:rhop}.

\section{Kinetic description\label{sec:3}}

As shown in the previous section, even at equilibrium, there are spatial correlations between the gases and the piston. These correlations should be kept in order to have an accurate description of a finite, out-of-equilibrium system. Moreover, a easier characterization of the thermal equilibrium of the gases is possible using the conditional distributions. For these reasons, in this section we derive equations for the conditional distribution functions $f_i(\mathbf r,\mathbf v,t|x_p)$ for $i=1,2$ of the gases, as well as for the distribution function of the piston $f_p(x_p,v_p,t)$.

\subsection{Preliminary definitions}

The distribution function of the gases $f_i(\mathbf r,\mathbf v,t)$ are defined such as $f_i(\mathbf r,\mathbf v,t)d\mathbf rd\mathbf v$ is the mean number of particles on the left ($i=1$) or on the right ($i=2$) with position between $\mathbf r$ and $\mathbf r+d\mathbf r$ and velocity between $\mathbf v$ and $\mathbf v+d\mathbf v$ at time $t$. They are related with the probability densities $\rho_i(\mathbf r,\mathbf p,t)$ as
\begin{equation}
  f_i(\mathbf r,\mathbf v,t)=N_im_i^{d}\rho_i(\mathbf r,\mathbf v/m_i,t), \quad i=1,2.
\end{equation}
The distribution function of the piston $f_p(x_p,v_p,t)$ is the probability density of finding the piston with position $x_p$ and velocity $v_p$ at time $t$, now
\begin{equation}
  f_p(x_p, v_p,t)=m_p\rho_p(x_p,v_p/m_p,t).
\end{equation}

We define the conditional distribution of the gases as
\begin{equation}
  \label{eq:cond}
  f_i(\mathbf r,\mathbf v,t|x_p)\equiv \frac{f_i(\mathbf r,\mathbf v,x_p,t)}{n_p(x_p,t)}, \quad i=1,2,
\end{equation}
where $f_i(\mathbf r,\mathbf v,x_p,t)$ is the joint probability distribution for a particle of gas $i$ with position $\mathbf r$ and velocity $\mathbf v$ and the piston with position $x_p$, regardless its velocity, at time $t$. It is convenient to express the new probability density in terms of the more general joint probability $f_{ip}(\mathbf r,\mathbf v,x_p,v_p,t)$ where the velocity of the piston is also specified,
\begin{equation}
  f_{i}(\mathbf r,\mathbf v,x_p,t)\equiv \int dv_p\ f_{ip}(\mathbf r,\mathbf v,x_p,v_p,t).
\end{equation}
The other quantity in Eq. \eqref{eq:cond} is the density of the piston $n_p(x_p,t)$ already defined in Eq. \eqref{eq:den1}, now given in terms of $f_p$ as
\begin{equation}
  n_p(x_p,t)=\int dv_p \ f_p(x_p,v_p,t).
\end{equation}

Finally, it is convenient for forthcoming discussions to define the conditional probability distribution of a position and a velocity of the gases to a position and velocity of the piston as
\begin{equation}
  f_i(\mathbf r,\mathbf v,t|x_p,v_p)\equiv \frac{f_{ip}(\mathbf r,\mathbf v,x_p,v_p,t)}{f_p(x_p,v_p,t)}, \quad i=1,2.
\end{equation}
It is forth observing that $f_i(\mathbf r,\mathbf v,t|x_p)\ne \int f_i(\mathbf r,\mathbf v,t|x_p,v_p)$ in general, due to velocity correlations. An exception is when the system is at thermal equilibrium.

\subsection{Set of Boltzmann-like kinetic equations}

\subsubsection{Kinetic equation for the piston}
If we were to proceed as usual, we should write down a set of equations for the distribution functions, of the gases and the piston, from the Liouville equation. This would result in a non-closed set of equations, since higher order distribution function would appear. The difference in the present case with respect to other problems is that the conventional molecular chaos hypothesis does not work, and new approximations are needed.

Take first the kinetic equation for $f_p(x_p,v_p,t)$, resulting form integrating the Liouville equation over all variables except the position and velocity of the system,
\begin{equation}
  \label{eq:pist}
  \partial_t f_p+v_p \partial_{x_p}f_p=J_{p1}[v_p|f_1,f_p]+J_{p2}[v_p|f_2,f_p].
\end{equation}
It is nothing but a probabilistic balance equation, where the lhs accounts for the free-of-collision motion of the piston and the rhs accounts for the collisions with the gases. The latter includes the collision operators,
\begin{eqnarray}
  \nonumber J_{p1}=&&\int d\mathbf r_\perp\int_{v_x^*>v_p^*}d\mathbf v\ (v_x^*-v_p^*)f_{1p}((x_p,\mathbf r_\perp),\mathbf v^*,x_p,v_p^*,t)\\ && -\int d\mathbf r_\perp\int_{v_x>v_p}d\mathbf v\ (v_x-v_p)f_{1p}((x_p,\mathbf r_\perp),\mathbf v,x_p,v_p,t), \\ \nonumber J_{p2}=&&\int d\mathbf r_\perp\int_{v_x^*<v_p^*}d\mathbf v\ (v_p^*-v_x^*)f_{2p}((x_p,\mathbf r_\perp),\mathbf v^*,x_p,v_p^*,t)\\ && -\int d\mathbf r_\perp\int_{v_x<v_p}d\mathbf v\ (v_p-v_x)f_{2p}((x_p,\mathbf r_\perp),\mathbf v,x_p,v_p,t),
\end{eqnarray}
which are functions of $v_p$ and functionals of the distribution functions. The asterisk on the velocities indicates precollisional values,
\begin{equation}
  \begin{split}
    v_x^*\equiv & v_x-\frac{2m_p}{m_i+m_p}(v_x-v_p), \\
    \mathbf v_\perp^*\equiv & \mathbf v_\perp, \\
    v_p^*\equiv & v_p-\frac{2m_i}{m_i+m_p}(v_p-v_x),
  \end{split}
\end{equation}
for $i=1$ or $2$. The functional dependence of the operator can be made more clear if we replace the joint probabilities using the following relation
\begin{equation}
  f_{ip}((x_p,\mathbf r_\perp),\mathbf v,x_p,v_p,t)=f_{i}((x_p,\mathbf r_\perp),\mathbf v,t|x_p,v_p) f_p(x_p,v_p,t).
\end{equation}

\paragraph{Assumption 1.} The usual molecular chaos hypothesis for precollisions is translated in this work to
\begin{equation}
  f_{i}((x_p,\mathbf r_\perp),\mathbf v,t|x_p,v_p)\simeq f_{i}((x_p,\mathbf r_\perp),\mathbf v,t|x_p),
\end{equation}
for precollisions. That is, we keep the correlations in the position but remove precollisional velocity correlations.

Observe that the functional dependence of the collision operators are evaluated when a particle is at the position of the piston, dissregarding the diameter of the particle. Since we are interested in a mesoscopic description where the relevant spatial variation of the distribution function occur over a mean free path of the particles, the simplification is expected to be accurate for the gases dilute enough.

\subsubsection{Kinetic equations for the gases}

To complete the kinetic description, we need equations for the conditional probabilities $f_i(\mathbf r,\mathbf v,t|x_p)$. For that purpose, and according to the definition of Eq. \eqref{eq:cond}, we need equations for $n_p(x_p,t)$ and $f_{ip}(\mathbf r,\mathbf v,x_p,v_p,t)$.

Integrating Eq. \eqref{eq:pist}, we get
\begin{eqnarray}
  \label{eq:np}
  \nonumber
    \partial_tn_p(x_p,t)= && \int dv_p \partial_tf_p(x_p,v_p,t)=\int dv_p\left(-v_p\partial_{x_p}f_p+J_{p1}+J_{p2}\right) \\ =&&
    -\partial_{x_p}\left[u_p(x_p,t)n_p(x_p,t)\right],
\end{eqnarray}
where
\begin{equation}
  u_p(x_p,t)=\frac{1}{n_p(x_p,t)}\int dv_p\ v_p f_p(x_p,v_p,t)
\end{equation}
is the mean velocity of the piston. In Eq.~\eqref{eq:np} we have used the conservation of the number of particles in collisions,
\begin{equation}
    \int dv_p J_{p1}=\int dv_p J_{p2}=0.
\end{equation}

The kinetic equation for the joint probability $f_{1p}(\mathbf r,\mathbf v,x_p,v_p,t)$, resulting from the Liouville equation after a proper integration, is
\begin{equation}
  \label{eq:f1p}
  \partial_t f_{1p}=-(\mathbf v\cdot \partial_{\mathbf r}+v_p\partial_{x_p})f_{1p}+J_{1p,11'}+J_{1p,p1'}+J_{1p,p2}
\end{equation}
where the collision operators take into account the different mechanisms the distribution function $f_{1p}$ can change through collisions. For the operator $J_{1p,11'}$ the collision is between two particles on the left of the piston, for $J_{1p,p1'}$ the collision is between the piston and a second particle on the left, and $J_{1p,p2}$ accounts for collisions between the piston and a particle on the right. These collision terms involve three-body distribution functions. There is still one kind of collision to consier, involving a particle with position and velocity given by the argument of $f_{1p}$ and the piston, when they are in contact ($x=x_p)$ and about to collide ($v_x\ge v_p$). In this case, the usual procedure is to impose a boundary condition to $f_{1p}$, which is a consequence of the conservation of probability:
\begin{equation}
  \label{eq:f1pbc}
  (v_p'-v_x')f_{1p}((x_p,\mathbf r_\perp),\mathbf v',x_p,v'_p,t)=(v_x-v_p)f_{1p}((x_p,\mathbf r_\perp),\mathbf v,x_p,v_p,t),
\end{equation}
for  $v_x\ge v_p$. The primed variables are given by the collision rule in Eqs. \eqref{eq:coll}-\eqref{eq:coll2}. For the gas on the right, an analogous equation and boundary condition result.

Before proceeding with the derivation of the equation for $f_i(\mathbf r,\mathbf v,t|x_p)$, it is worth focusing on some consequences of the boundary condition in Eq. \eqref{eq:f1pbc}. If we integrate it over velocities for $v_x\ge v_p$, after a change of variable and taking into account that $dv_xdv_p=dv_x'dv_p'$, we get
\begin{equation}
  \label{eq:bvel}
  n_p(x_p,t)n_i(x_p,t)u_p(x_p,t)=n_p(x_p,t)n_i(x_p,t)u_{i,x}(x_p,t), \quad i=1,2,
\end{equation}
with the mean velocity of the gases $\mathbf u_i(\mathbf r,t)$ ($i=1,2$) defined as
\begin{equation}
  n_i(\mathbf r,t)\mathbf u_i(\mathbf r,t)\equiv \int d\mathbf v\ \mathbf v f_i(\mathbf r,\mathbf v,t), \quad i=1,2.
\end{equation}
For nonzero densities, the velocity of the gases and the piston in contact coincide.

Integrating Eq. \eqref{eq:f1p} over the velocity of the piston, we have
\begin{equation}
  \partial_t \int dv_p f_{1p}= -\mathbf v\cdot \partial_{\mathbf r}\int dv_p f_{1p}-\partial_{x_p}\int dv_p\ v_p f_{1p}+\int dv_p\ J_{1p,11'},
\end{equation}
where we have used the fact that the collision operators conserve the number of particles. Using last equation and Eq. \eqref{eq:np} with the temporal derivative of Eq. \eqref{eq:cond}, after some manipulations, we arrive at
\begin{eqnarray}
  \nonumber
  \partial_t f_1(\mathbf r,\mathbf v,t|x_p)=&&f_1(\mathbf r,\mathbf v,t|x_p)\frac{\partial_{x_p}\left[u_p(x_p,t)n_p(x_p,t)\right]}{n_p(x_p,t)}-\frac{1}{n_p(x_p,t)}\mathbf v\cdot \partial_{\mathbf r}\int dv_p f_{1p}\\ && -\frac{1}{n_p(x_p,t)}\partial_{x_p}\int dv_p\ v_p f_{1p}+\frac{1}{n_p(x_p,t)}\int dv_p\ J_{1p,11'},
\end{eqnarray}
which is not a closed equation yet.

\paragraph{Assumptions 2 and 3.} In order to close the equation we include additional approximations  that go beyond the standard molecular chaos hypothesis, namely
\begin{eqnarray}
  \label{eq:ap1}
  \nonumber
  \int dv_p\ v_pf_{1p}(\mathbf r,\mathbf v,x_p,v_p,t)=&& \int dv_p\ v_pf_{1p}(\mathbf r,\mathbf v,t|x_p,v_p)f_p(x_p,v_p,t)\\ \simeq&& \nonumber f_{1p}(\mathbf r,\mathbf v,t|x_p)\int dv_p\ v_pf_p(x_p,v_p,t)\\= && n_p(x_p,t)u_p(x_p,t)f_{1p}(\mathbf r,\mathbf v,t|x_p),
\end{eqnarray}
and
\begin{equation}
  \label{eq:ap2}
  \int dv_p\ J_{1p,11'}\simeq  J_{11} \int dv_p\ f_{1p}(\mathbf r,\mathbf v,t|x_p)=n_p(x_p,t)J_{11},
\end{equation}
with $J_{11}$ the collision operator of two particles on the left, which is a function of $\mathbf v$ and a functional of $f_1(\mathbf r,\mathbf v,t|x_p)$.  For hard-sphere collisions, the latter approximation is
\begin{eqnarray}
  \nonumber
\int dv_p\ J_{1p,11'}    &&=\int dv_p\ f_p(x_p,v_p,t) \int d\mathbf v_2 d\hat{\boldsymbol{\sigma}}\ \Theta [(\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}}](\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}} \\ && \nonumber \quad \times  (b^{-1}-1) f_{1}(\mathbf r,\mathbf v,t|x_p,v_p)f_{1}(\mathbf r,\mathbf v_2,t|x_p,v_p) \\ \nonumber &&\simeq \int dv_p\ f_p(x_p,v_p,t) \int d\mathbf v_2 d\hat{\boldsymbol{\sigma}}\ \Theta [(\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}}](\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}} \\ && \nonumber \quad \times  (b^{-1}-1) f_{1}(\mathbf r,\mathbf v,t|x_p)f_{1}(\mathbf r,\mathbf v_2,t|x_p)\\ &&= n_p(x_p,t)J_{11},
\end{eqnarray}
where $b^{-1}$ is the restitution operator, which replaces the velocities of the colliding particles by their precollisional values, and
\begin{eqnarray}
  \nonumber
  J_{11}= \int d\mathbf v_2 d\hat{\boldsymbol{\sigma}}\ && \Theta [(\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}}](\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}} \\ && \times (b^{-1}-1) f_{1}(\mathbf r,\mathbf v,t|x_p)f_{1}(\mathbf r,\mathbf v_2,t|x_p).
\end{eqnarray}

Finally, we arrive at a closed equation for $f_1(\mathbf r,\mathbf v,t|x_p)$,
\begin{equation}
  \label{eq:eqf1}
  \partial_t f_1(\mathbf r,\mathbf v,t|x_p) \simeq-\mathbf v\cdot \partial_{\mathbf r}f_1(\mathbf r,\mathbf v,t|x_p)-u_p(x_p,t)\partial_{x_p}f_1(\mathbf r,\mathbf v,t|x_p) +J_{11}.
\end{equation}
With similar approximation, we also obtain the kinetic equation for the gas on the right of the piston,
\begin{equation}
  \label{eq:eqf2}
  \partial_t f_2(\mathbf r,\mathbf v,t|x_p) \simeq-\mathbf v\cdot \partial_{\mathbf r}f_2(\mathbf r,\mathbf v,t|x_p)-u_p(x_p,t)\partial_{x_p}f_2(\mathbf r,\mathbf v,t|x_p) +J_{22},
\end{equation}
with
\begin{eqnarray}
  \nonumber
  J_{22}= \int d\mathbf v_2 d\hat{\boldsymbol{\sigma}}\ && \Theta [(\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}}](\mathbf{v}-\mathbf{v}_2)\cdot \hat{\boldsymbol{\sigma}} \\ && \times (b^{-1}-1) f_{2}(\mathbf r,\mathbf v,t|x_p)f_{2}(\mathbf r,\mathbf v_2,t|x_p).
\end{eqnarray}

For the boundary conditions, we have to differentiate from particles in contact with the walls of the cylinder, denoted by $\partial V$, and the piston. In the first case, the conservation of probability imposes
\begin{equation}
  \label{eq:bcdv}
    f_i(\mathbf r,\mathbf v,t|x_p)=f_i(\mathbf r,\mathbf v',t|x_p), \quad \mathbf r\in \partial V, \quad i=1,2,
  \end{equation}
  for $\mathbf v\cdot \hat{\boldsymbol{n}}=-\mathbf v'\cdot \hat{\boldsymbol{n}}$, where $\hat{\boldsymbol{n}}$ is a unit vector normal to $\partial V$.

  \paragraph{Assumption 4.} For the boundary conditions at the position of the piston, different approximations are possible. We take the one resulting from Eq. \eqref{eq:f1pbc}, after integration over the velocity of the piston and assuming $f_1(\mathbf r,\mathbf v,t|x_p,v_p) \simeq f_1(\mathbf r,\mathbf v,t|x_p) $ inside the integrals,
\begin{eqnarray}
  \label{eq:bc1}
\nonumber
  && \int_{v_x>v_p} dv_p\ (v_p'-v_x')f_{1}((x_p,\mathbf r_\perp),\mathbf v',t|x_p)f_p(x_p,v_p',t) \\ && \qquad \simeq  \int_{v_x>v_p}\ (v_x-v_p)f_{1}((x_p,\mathbf r_\perp),\mathbf v,t|x_p)f_p(x_p,v_p,t).
\end{eqnarray}
and
\begin{eqnarray}
  \label{eq:bc2}
\nonumber
  && \int_{v_x<v_p} dv_p\ (v_p'-v_x')f_{2}((x_p,\mathbf r_\perp),\mathbf v',t|x_p)f_p(x_p,v_p',t) \\ && \qquad \simeq  \int_{v_x<v_p}\ (v_x-v_p)f_{2}((x_p,\mathbf r_\perp),\mathbf v,t|x_p)f_p(x_p,v_p,t).
\end{eqnarray}
These approximate boundary conditions have the advantage, if compared to others, to conserve the exact relation of Eq. \eqref{eq:bvel}, which is of crucial importance for the conservation of particles, as seen later.

\subsection{Canonical equilibrium}

The proposed system of equations \eqref{eq:pist},\eqref{eq:eqf1}-\eqref{eq:eqf2} with boundary conditions \eqref{eq:bc1}-\eqref{eq:bc2} has the canonical equilibrium as a solution, given by
\begin{eqnarray}
  \label{eq:f1eq}
  f_1((x,\mathbf r_\perp),\mathbf v,t|x_p)=&&\frac{N_1}{x_pS}\frac{m_1^d\lambda_1^d}{h^d}e^{-\beta\frac{1}{2}m_1v^2}, \quad x\le x_p, \\
  \label{eq:f2eq}
    f_2((x,\mathbf r_\perp),\mathbf v,t|x_p)=&&\frac{N_1}{(L-x_p)S}\frac{m_2^d\lambda_2^d}{h^d}e^{-\beta\frac{1}{2}m_2v^2}, \quad x\ge x_p, \\
  \label{eq:fpeq}
  f_p(x_p,v_p,t)=&&n_p(x_p)\frac{m_p\lambda_p}{h}e^{-\beta\frac{1}{2}m_pv_p^2},
\end{eqnarray}
with
\begin{equation}
  \label{eq:nxp}
n_p(x_p)=\frac{(N+1)!}{N_1!N_2!L}\left(\frac{x}{L}\right)^{N_1}\left(1-\frac{x}{L}\right)^{N_2}.
\end{equation}

To demonstrate the latter assertion, take $J_{p1}$ at canonical equilibrium. Since $v_x^*-v_p^*=-(v_x-v_p)$ and $m_1v_x^{*2}+m_pv_p^{*2}=m_1v_x^{2}+m_pv_p^{2}$, after integrating over $\mathbf v_\perp$, it is
\begin{eqnarray}
  \nonumber
  J_{p1}[v_p|f_p,f_1]=&&-\frac{N_1m_1\lambda_1m_p\lambda_p}{x_ph^2}\rho_p(x_p)\int_{v_x<v_p}dv_x\ (v_x-v_p)e^{-\beta\frac{1}{2}m_pv_p^2-\beta\frac{1}{2}m_1v_x^2}\\   \nonumber
 &&  -\frac{N_1m_1\lambda_1m_p\lambda_p}{x_ph^2}\rho_p(x_p)\int_{v_x>v_p}dv_x\ (v_x-v_p)e^{-\beta\frac{1}{2}m_pv_p^2-\beta\frac{1}{2}m_1v_x^2} \\ =&&\frac{N_1m_p\lambda_p}{x_ph}n_p(x_p)v_pe^{-\beta\frac{1}{2}m_pv_p^2}=\frac{N_1}{x_p}v_pf_p(x_p,v_p,t).
\end{eqnarray}
Similarly,
\begin{equation}
  \label{eq:31}
  J_{p2}[v_p|f_p,f_2]=-\frac{N_2}{L-x_p}v_pf_p(x_p,v_p,t).
\end{equation}
Hence,
\begin{equation}
  \frac{d}{dx_p}n_p(x_p)=\frac{N_1}{x_p}-\frac{N_2}{L-x_p}v_pf_p(x_p,v_p,t),
\end{equation}
which has expression \eqref{eq:nxp} as the unique normalized solution.

The considerations for the conditional probabilities are straightforward. Since $J_{11}=J_{22}=0$ for gaussian distributions, Eqs. \eqref{eq:eqf1} and \eqref{eq:eqf2} clearly support the equilibrium solutions.

\subsection{Some comments and remarks}

The equations for the distribution function of the piston \eqref{eq:pist} and the conditional distribution functions for the gases \eqref{eq:eqf1}-\eqref{eq:eqf2} constitute one of the main results of the present work. The chose of the conditional distribution as the function to describe the gases has been motivated by the analysis of the canonical equilibrium, since there we see the need of including spatial correlations between gases and the piston, and the equation of the piston where the conditional probabilities emerge in a natural way.

The derivation of the kinetic equations used several approximations, namely the molecular chaos for inter-particle and piston-particle collisions (assumption 1), and a generalization of molecular chaos for particle-piston collisions (assumptions 2-4). The molecular chaos hypothesis neglects pre-collisional velocity correlations, as usual when dealing with the Boltzmann equations. The generalized molecular chaos hypothesis goes one step further by also removing ``some'' post-collisional velocity correlations, as specified by Eqs. \eqref{eq:ap1}-\eqref{eq:ap2} and \eqref{eq:bc1}-\eqref{eq:bc2}.  The approximate expression of Eqs. \eqref{eq:ap1}-\eqref{eq:ap2} are expected to be correct as far as $f_i(\mathbf r,\mathbf v,t|v_p,x_p)\simeq f_i(\mathbf r,\mathbf v,t|x_p)$ inside the integrals. This, in turn, is expected to be true in a general situation with $x$ not very close to $x_p$, since collisions between particles and the piston, as well as the collisions among particles, tend to destroy the velocity correlations. When $x\simeq x_p$, which is the case of Eqs. \eqref{eq:bc1}-\eqref{eq:bc2}, the removal of correlations is not that obvious. In fact, the boundary conditions  \eqref{eq:bc1}-\eqref{eq:bc2} are strictly wrong. However, if the section of the piston $S$ is big enough $S/\sigma_i^{d-1}\gg 1$, we could find a regime of dilute gases where the postcollisional velocity correlations decay to zero on a distance much smaller than the mean free path of the particles, meaning that the boundary conditions become correct in the mesoscale. An exact analysis of this considerations needs a formulation with the functions $f_i(\mathbf r,\mathbf v,t|v_p,x_p)$, which is beyond the scope of the present work.

Observe that the structure of the kinetic equations for the conditional distributions of the gases  is different from the one expected from naive considerations, as terms proportional to the mean velocity of the piston appear. Although the presence of these terms can be justified by means of physical considerations, i.e they are important in order to compute the work done by the piston correctly, they can compromise some key properties of the equation. Namely, for a given physical initial condition, the kinetic equations for the conditional distribution functions should provide positive and normalized solutions for all later times. Both aspects could be indirectly demonstrated from the equation of $f_{ip}$ and $n_p$ or directly. As an example, we show in the sequel that the kinetic equations \eqref{eq:pist}-\eqref{eq:eqf1}-\eqref{eq:eqf2} preserve the normalization of distribution functions.

Consider first the equation of the density of the piston Eq. \eqref{eq:np}. Integrating over $x_p\in(0,L)$, using the divergence theorem, and the fact that the piston can not leave the system, $n_p(0,t)u_p(0,t)=n_p(L,t)u_p(L,t)=0$, it turns out that $\int dx_p\ n_p(x_p,t)$ is a conserved quantity. Hence, the normalization of $f_p$ is ensured, provided it is initially normalized.

Consider now the gas on the left. After integration of Eq. \eqref{eq:eqf1} over $\mathbf v$, we face several terms:
\begin{equation}
  \int d\mathbf v\ \partial_t f_1(\mathbf r,\mathbf v,t|x_p)=\partial_tn_1(\mathbf r,t|x_p),
\end{equation}
with
\begin{equation}
  n_1(\mathbf r,t|x_p)\equiv \int d\mathbf v\ f_1(\mathbf r,\mathbf v,t|x_p)
\end{equation}
being the conditional density;
\begin{equation}
  \int d\mathbf v\ \mathbf v\cdot \partial_{\mathbf r}f_1(\mathbf r,\mathbf v,t|x_p)=\partial_{\mathbf r}\cdot[n_1(\mathbf r,t|x_p)\mathbf u_1(\mathbf r,t|x_p)],
\end{equation}
with the conditional velocity defined as
\begin{equation}
  \mathbf u_1(\mathbf r,t|x_p)\equiv \frac{1}{n_1(\mathbf r,t|x_p)}\int d\mathbf v \ \mathbf v f_1(\mathbf r,\mathbf v,t|x_p);
\end{equation}
\begin{equation}
  \int d\mathbf v\ u_p(x_p,t)\partial_{x_p}f_1(\mathbf r,\mathbf v,t|x_p)=u_p(x_p,t)\partial_{x_p}n_1(\mathbf r,t|x_p);
\end{equation}
and
\begin{equation}
 \int d\mathbf v\ J_{11}=0.
\end{equation}
Hence, the balance equation for the conditional density is
\begin{equation}
  \partial_{t}n_1(\mathbf r,t|x_p)=-\partial_{\mathbf r}\cdot[n_1(\mathbf r,t|x_p)\mathbf u_1(\mathbf r,t|x_p)]-u_p(x_p,t)\partial_{x_p}n_1(\mathbf r,t|x_p),
\end{equation}
and a similar one for the conditional density of the gas on the right of the piston $n_2(\mathbf r,t|x_p)$. Integrating over $x\in(0,x_p)$ and all possible values of $\mathbf r_\perp$:
\begin{eqnarray}
  \nonumber
  \partial_t N_1(t|x_p)=&& -\left. \int d\mathbf r_\perp\ n_1(\mathbf r,t|x_p) u_{1,x}(\mathbf r,t|x_p)\right|_{x=x_p} \\ &&-u_p(x_p,t)\int d\mathbf r_\perp \int_0^{x_p}dx \partial_{x_p}n_1(\mathbf r,t|x_p),
\end{eqnarray}
with
\begin{equation}
  N_1(t|x_p)\equiv \int d\mathbf r\ n_1(\mathbf r,t|x_p),
\end{equation}
We have used the condition
\begin{equation}
  \mathbf u_1(x\in \partial V,t|x_p)\cdot \hat{\boldsymbol{n}}=0,
\end{equation}
with $\hat{\boldsymbol{n}}$ a unit vector normal to boundary of the cylinder $\partial V$, which is a direct consequence of the kinetic boundary condition \eqref{eq:bcdv}. The quantity $u_{1,x}$ is the $X$ component of $\mathbf u_1(\mathbf r,t|x_p)$. Now,
\begin{equation}
  \int_0^{x_p}dx \partial_{x_p}n_1(\mathbf r,t|x_p)=\partial_{x_p}N_1(t|x_p)-n_1(x_p,t|x_p).
\end{equation}
Hence,
\begin{eqnarray}
  \nonumber
  \partial_t N_1(x_p,t)=&&-\left. \int d\mathbf r_\perp\ n_1(\mathbf r,t|x_p) [u_{1,x}(\mathbf r,t|x_p) -u_p(x_p,t)]\right|_{x=x_p}\\ && -u_p(x_p,t)\partial_{x_p}N_1(t|x_p).
\end{eqnarray}
The first term on the rhs is zero, using the kinetic boundary condition at the piston \eqref{eq:bc1},
\begin{eqnarray}
  \nonumber
  && \int_{v_x>v_p} dv_xdv_p(v_x-v_p)f_1(x_p,v_x,t|x_p)f_p(x_p,v_p,t)\\ && +\int_{v_x>v_p} dv_xdv_p (v_x'-v_p')f_1(x_p,v_x',t|x_p)f_p(x_p,v_p',t)\simeq 0
\end{eqnarray}
which, making a change of variable on the second integral, turns
\begin{eqnarray}
  \nonumber
  && \int dv_xdv_p(v_x-v_p)f_1(x_p,v_x,t|x_p)f_p(x_p,v_p,t)\simeq 0 \\ \nonumber &&\Rightarrow n_1(x_p,t|x_p)n_p(x_p,t)(u_{1,x}-u_{p})\simeq 0\\ && \Rightarrow n_1(x_p,t|x_p)(u_{1,x}-u_{p})\simeq 0,
\end{eqnarray}
if $n_p(x_p,t)\ne 0$. Hence, the balance equation for $N_1$ results
\begin{equation}
  \partial_t N_1(t|x_p)=-u_p(x_p,t)\partial_{x_p}N_1(t|x_p).
\end{equation}
From the nature of this equation we infer that if the conditional distribution function for the gas on the left is normalized to $N_1$  for all allowed values of $x_p$, that is $N_1(0|x_p)=N_1$, then $N_1(t|x_p)=N_1$ for $t\ge 0$.

 Similar considerations allow us to conclude that the conditional distribution of the gas on the right of the piston is also normalized to $N_2$ for all times.

\section{$\mathcal{H}$-theorem\label{sec:4}}
In this section we demonstrate the $\mathcal H$-theorem for the system of kinetic equations derived along the previous section. It states that the $\mathcal H$ function, to be defined below, mush approach a limit where the distribution functions are of thermal equilibrium given by Eqs. \eqref{eq:f1eq}-\eqref{eq:fpeq}. We first prove $\mathcal H$ is a decreasing function, then that its time derivative is zero at thermal equilibrium, and finally that it is bounded from below by its value at equilibrium.

\subsection{Definition of $\mathcal H$}
Following a recent work \cite{magabr18}, take the $\mathcal{H}$ function as
\begin{equation}
  \mathcal H\equiv \int d\Gamma \rho_N(\Gamma,t)\ln\left[\frac{\rho_N(\Gamma,t)}{\rho_{01}^{N_1}\rho_{02}^{N_2}\rho_{0p}}\right],
\end{equation}
with $\rho_N$ the probability density of the whole system, $\Gamma\equiv (\mathbf R,x_p;\mathbf P,p_p)$, and $\rho_{0i}$ constants that make the expression dimensionless. Removing velocity correlations, as usual, it is
\begin{eqnarray}
  \nonumber
  \rho_N(\mathbf R,x_p;\mathbf P,p_p,t)&&=\rho_N(\mathbf R;\mathbf P,t|p_p,x_p)\rho_p(x_p,p_p,t)\\ \nonumber && \simeq \rho_N(\mathbf R;\mathbf P,t|x_p)\rho_p(x_p,p_p,t)\\ && \simeq \prod_{i=1}^{N_1}\rho_1(\mathbf r_i,\mathbf p_i,t|x_p)\prod_{j=1}^{N_2}\rho_2(\mathbf r_j,\mathbf p_j,t|x_p)\rho_p(x_p,p_p,t),
\end{eqnarray}
and the $\mathcal H$-function turns
\begin{eqnarray}
  \nonumber
  \mathcal H \simeq && N_1 \int d\mathbf rdx_pd\mathbf pdp_p\ \rho_1(\mathbf r,\mathbf p,t|x_p)\rho_p (x_p,p_p,t)\ln\left[\frac{\rho_1(\mathbf r,\mathbf p,t|x_p)}{\rho_{01}}\right]  \\  \nonumber && + N_2 \int d\mathbf rdx_pd\mathbf pdp_p\ \rho_2(\mathbf r,\mathbf p,t|x_p)\rho_p (x_p,p_p,t)\ln\left[\frac{\rho_2(\mathbf r,\mathbf p,t|x_p)}{\rho_{02}}\right] \\ && + \int dx_p dp_p\ \rho_p(x_p,p_p,t)\ln\left[\frac{\rho_p(x_p,p_p,t)}{\rho_{0p}}\right].
\end{eqnarray}
In terms of the distribution functions,
\begin{eqnarray}
  \label{eq:etaf}
  \nonumber
  \mathcal H \simeq && \int d\mathbf rdx_pd\mathbf vdv_p\ f_1(\mathbf r,\mathbf v,t|x_p)f_p (x_p,v_p,t)\ln\left[\frac{f_1(\mathbf r,\mathbf v,t|x_p)}{f_{01}}\right]  \\ \nonumber && + \int d\mathbf rdx_pd\mathbf vdv_p\ f_2(\mathbf r,\mathbf v,t|x_p)f_p (x_p,v_p,t)\ln\left[\frac{f_2(\mathbf r,\mathbf v,t|x_p)}{f_{02}}\right] \\ && + \int dx_p dv_p\ f_p(x_p,v_p,t)\ln\left[\frac{f_p(x_p,v_p,t)}{f_{0p}}\right],
\end{eqnarray}
with the new constants $f_{0i}$ having the dimensions of their respective distribution functions.

\subsection{Proof of $\frac{d}{dt}\mathcal H \le 0$}
Now we compute the time derivative of $\mathcal H$ as given by its last expression. For the sake of simplicity on the notation, we omit the obvious dependence of the different functions.

The time derivative of the first term  of Eq. \eqref{eq:etaf} is
\begin{eqnarray}
  \nonumber
    \frac{d}{dt}\int d\mathbf rdx_pd\mathbf v\ f_1n_p\ln\left(\frac{f_1}{f_{01}}\right)= && \int d\mathbf rdx_pd\mathbf v\ \left[1+\ln\left(\frac{f_1}{f_{01}}\right)\right]n_p\partial_tf_1 \\ && + \int d\mathbf rdx_pd\mathbf v\ f_1\ln\left(\frac{f_1}{f_{01}}\right)\partial_tn_p,
\end{eqnarray}
which is, after using the equations for $f_1$ and $n_p$ \eqref{eq:eqf1} and \eqref{eq:np},
\begin{eqnarray}
  \nonumber
  && \int d\mathbf rdx_pd\mathbf v\ \left[1+\ln\left(\frac{f_1}{f_{01}}\right)\right]n_p\left(-\mathbf v\cdot \partial_{\mathbf r}f_1-u_p\partial_{x_p}f_1+J_{11}\right) \\ &&- \int d\mathbf rdx_pd\mathbf v\ f_1\ln\left(\frac{f_1}{f_{01}}\right)\partial_{x_p}\left(u_pn_p\right).
\end{eqnarray}
After some algebra,
\begin{eqnarray}
  && \left[1+\ln\left(\frac{f_1}{f_{01}}\right)\right]n_p \left(-\mathbf v\cdot \partial_{\mathbf r}f_1\right) = -\mathbf v\cdot \partial_{\mathbf r}\left[\ln\left(\frac{f_1}{f_{01}}\right)n_p f_1\right], \\
  \nonumber
  && \left[1+\ln\left(\frac{f_1}{f_{01}}\right)\right]n_p \left(- u_p \partial_{x_p}f_1\right) =-\partial_{x_p}\left[\ln\left(\frac{f_1}{f_{01}}\right)n_pu_pf_1\right]\\ && \qquad +f_1\ln\left(\frac{f_1}{f_{01}}\right)\partial_{x_p}(u_pn_p),
\end{eqnarray}
and their respective integrals include
\begin{eqnarray}
  \nonumber
  -&&\int d\mathbf v \int_0^L dx_p\int d\mathbf r_\perp\int_0^{x_p} d x\ \mathbf v\cdot \partial_{\mathbf r}\left[\ln\left(\frac{f_1}{f_{01}}\right)n_p f_1\right]\\ && =-\int d\mathbf v \int_0^L dx_p\int d\mathbf r_\perp\ v_x\left[\ln\left(\frac{f_1}{f_{01}}\right)n_p f_1\right]_{x=x_p}, \\
\nonumber    -&&\int d\mathbf r_\perp\int_0^L dx_p\int_0^{x_p} d x\ \partial_{x_p}\left[\ln\left(\frac{f_1}{f_{01}}\right)n_pu_pf_1\right]\\  && = \int d\mathbf r \left[\ln\left(\frac{f_1}{f_{01}}\right)n_pu_pf_1\right]_{x_p=x}.
\end{eqnarray}
After some manipulations, we finally arrive at
\begin{eqnarray}
  \nonumber
  \frac{d}{dt}\int d\mathbf rdx_pd\mathbf v\ f_1n_p \ln\left(\frac{f_1}{f_{01}}\right)=&& \int d\mathbf r d\mathbf v\ (u_p-v_x)\left.\ln\left(\frac{f_1}{f_{01}}\right)n_pf_1\right|_{x_p=x}\\ &&+ \int d\mathbf rdx_pd\mathbf v\ \ln\left(\frac{f_1}{f_{01}}\right)n_pJ_{11},
\end{eqnarray}
and similarly
\begin{eqnarray}
  \nonumber
  \frac{d}{dt}\int d\mathbf rdx_pd\mathbf v\ f_2n_p \ln\left(\frac{f_2}{f_{02}}\right)=& &\int d\mathbf r d\mathbf v\ (v_x-u_p)\left.\ln\left(\frac{f_2}{f_{02}}\right)n_pf_2\right|_{x_p=x}\\ &&+ \int d\mathbf rdx_pd\mathbf v\ \ln\left(\frac{f_2}{f_{02}}\right)n_pJ_{22}.
\end{eqnarray}
It is still possible to rewrite the terms involving the piston on the last two expression in a form which will be useful soon. After some manipulations and the use of the kinetic boundary conditions \eqref{eq:bc1}-\eqref{eq:bc2}, we have
\begin{eqnarray}
  \nonumber
    && \int d\mathbf r d\mathbf v\ (u_p-v_x)\left.\ln\left(\frac{f_1}{f_{01}}\right)n_pf_1\right|_{x_p=x}\\ && = \int d\mathbf r \int_{v_x>v_p}d\mathbf v dv_p\ (v_x-v_p)f_1f_p\left.\ln\left(\frac{f_1'}{f_{1}}\right)\right|_{x_p=x}
\end{eqnarray}
and
\begin{eqnarray}
   \nonumber
  && \int d\mathbf r d\mathbf v\ (u_x-v_p)\left.\ln\left(\frac{f_1}{f_{01}}\right)n_pf_1\right|_{x_p=x}\\ && = \int d\mathbf r \int_{v_x<v_p}d\mathbf v dv_p\ (v_p-v_x)f_2f_p\left.\ln\left(\frac{f_2'}{f_{2}}\right)\right|_{x_p=x},
\end{eqnarray}
where the primer on the functions indicates that the velocity of the respective arguments has to be replaced by its postcollision value.

For the time derivative of the last term of Eq.~\eqref{eq:etaf}, we have
\begin{equation}
  \label{eq:79}
  \begin{split}
\frac{d}{dt}  \int & dx_p dv_p\ f_p(x_p,v_p,t)\ln\left[\frac{f_p(x_p,v_p,t)}{f_{0p}}\right]\\ &=\int dx_p dv_p\ \left(1+\ln\frac{f_p}{f_{0p}}\right)\left[-v_p\partial_{x_p}f_p+J_{p1}+J_{p2}\right] \\ & =\int dx_p dv_p\ \ln\frac{f_p}{f_{0p}}\left(J_{p1}+J_{p2}\right),
  \end{split}
\end{equation}
where we have used the conservation of particles in collisions and the fact that $f_p(0,v_p,t)=f_p(L,v_p,t)=0$. Using the expressions for the collision operators and after some manipulations,
\begin{equation}
  \int dx_p dv_p\ \ln\frac{f_p}{f_{0p}} J_{p1}=\int d\mathbf r \int_{v_x>v_p}\left. d\mathbf v dv_p\  (v_x-v_p)f_1f_p\ln\frac{f_p'}{f_{p}}\right|_{x_p=x}
\end{equation}
and
\begin{equation}
  \int dx_p dv_p\ \ln\frac{f_p}{f_{0p}} J_{p2} =\int d\mathbf r \int_{v_x<v_p}\left. d\mathbf v dv_p\  (v_p-v_x)f_2f_p\ln\frac{f_p'}{f_{p}}\right|_{x_p=x}.
\end{equation}

Collecting the different terms,
\begin{equation}
  \begin{split}
    \frac{d}{dt}\mathcal H&= \int d\mathbf r \int_{v_x>v_p}d\mathbf v dv_p\ (v_x-v_p)f_1f_p\left.\ln\left(\frac{f_1'f_p'}{f_{1}f_p}\right)\right|_{x_p=x} \\ &+ \int d\mathbf r \int_{v_x<v_p}d\mathbf v dv_p\ (v_p-v_x)f_2f_p\left.\ln\left(\frac{f_2'f_p'}{f_{2}f_p}\right)\right|_{x_p=x}\\ &+ \int d\mathbf rdx_pd\mathbf v\ n_p\left[ \ln\left(\frac{f_1}{f_{01}}\right)J_{11}+ \ln\left(\frac{f_2}{f_{02}}\right)J_{22}\right].
  \end{split}
\end{equation}
As it is well known, the last term on the rhs is negative \cite{mc89}. The first term on the rhs is, using the inequality $x\ln\frac{y}{x}\le y-x$ for positive numbers $x$ and $y$,
\begin{equation}
  \begin{split}
    &\int d\mathbf r \int_{v_x>v_p}d\mathbf v dv_p\ (v_x-v_p)f_1f_p\left.\ln\left(\frac{f_1'f_p'}{f_{1}f_p}\right)\right|_{x_p=x} \\ &\le \int d\mathbf r \int_{v_x>v_p}d\mathbf v dv_p\ (v_x-v_p)\left.\left(f_1'f_p'-f_{1}f_p\right)\right|_{x_p=x},
  \end{split}
\end{equation}
which implies, using the fact that $v_x-v_p=v_p'-v_x'$ and the kinetic boundary condition,
\begin{equation}
  \begin{split}
    &\int d\mathbf r \int_{v_x>v_p}d\mathbf v dv_p\ (v_x-v_p)f_1f_p\left.\ln\left(\frac{f_1'f_p'}{f_{1}f_p}\right)\right|_{x_p=x} \\ &\le \int d\mathbf r \int_{v_x>v_p}d\mathbf v dv_p\ (v_x-v_p)\left. f_1f_p\right|_{x_p=x} \\ & - \int d\mathbf r \int_{v_x>v_p}d\mathbf v dv_p\ (v_x-v_p)\left.f_{1}f_p\right|_{x_p=x}=0.
  \end{split}
\end{equation}
Analogously,
\begin{equation}
  \begin{split}
    &\int d\mathbf r \int_{v_x<v_p}d\mathbf v dv_p\ (v_p-v_x)f_2f_p\left.\ln\left(\frac{f_2'f_p'}{f_{2}f_p}\right)\right|_{x_p=x} \le 0.
  \end{split}
\end{equation}

Hence,
\begin{equation}
  \frac{d}{dt}\mathcal H\le 0.
\end{equation}

Following the argument in \cite{mc89}, for $\frac{d}{dt}\mathcal H$ to be zero, we infer the need of the distribution function to be that of the canonical distribution, Eqs. \eqref{eq:f1eq}-\eqref{eq:fpeq}.

\subsection{Proof of $\mathcal H\ge \mathcal H_0$}
Let $\mathcal H_0$ be the function $\mathcal H$ evaluated at thermal equilibrium, with Eqs. \eqref{eq:f1eq}-\eqref{eq:fpeq} and \eqref{eq:etaf}. It is
\begin{eqnarray}
  \label{eq:difh0}
  \nonumber
  \mathcal H-\mathcal H_0= && \int d\mathbf rdx_pd\mathbf vdv_p\ f_1f_p\ln\left[\frac{f_1}{f_{1,eq}}\right]  + \int d\mathbf rdx_pd\mathbf vdv_p\ f_2f_p\ln\left[\frac{f_2}{f_{2,eq}}\right] \\ \nonumber && + \int dx_p dv_p\ f_p\ln\left[\frac{f_p}{f_{p,eq}}\right]\\ \nonumber && +\int d\mathbf rdx_pd\mathbf vdv_p\ (f_{1,eq}f_{p,eq}-f_1f_p)\ln\left[\frac{f_{1,eq}}{f_{01}}\right] \\ \nonumber && + \int d\mathbf rdx_pd\mathbf vdv_p\ (f_{2,eq}f_{p,eq}-f_2f_p)\ln\left[\frac{f_{2,eq}}{f_{02}}\right] \\ && + \int dx_p dv_p\ (f_{p,eq}-f_p)\ln\left[\frac{f_{p,eq}}{f_{0p}}\right],
\end{eqnarray}
where the suffix $eq$ denotes the equilibrium distributions.

The first term on the rhs is
\begin{eqnarray}
\label{eq:inec1}
\nonumber
\int d\mathbf rdx_pd\mathbf vdv_p\ && f_1f_p\ln\left[\frac{f_1}{f_{1,eq}}\right]= \\ && \int d\mathbf rdx_pd\mathbf vdv_p\ f_{1,eq}f_p\left\{1+\frac{f_1}{f_{1,eq}}\left[\ln\left(\frac{f_1}{f_{1,eq}}\right)-1\right]\right\},
\end{eqnarray}
where we have used the fact that
\begin{equation}
  \int d\mathbf rdx_pd\mathbf vdv_p\ f_{1,eq}f_p=\int d\mathbf rdx_pd\mathbf vdv_p\ f_{1}f_p.
\end{equation}
Since $1+x(\ln x-1)\ge 0$ for $x\ge 0$, the expression in Eq. \eqref{eq:inec1} is positive. A similar reasoning shows that the second and third terms of Eq. \eqref{eq:difh0} are also positive.

The last three terms of Eq. \eqref{eq:difh0} are
\begin{eqnarray}
  \label{eq:3t}
  -\frac{\beta}{2}&&\left[\int d\mathbf rdx_pd\mathbf vdv_p\ (f_{1,eq}f_{p,eq}-f_1f_p)m_1v_1^2\right. \\ \nonumber && + \int d\mathbf rdx_pd\mathbf vdv_p\ (f_{2,eq}f_{p,eq}-f_2f_p)m_2v_2^2\\ && + \left.\int dx_p dv_p\ (f_{p,eq}-f_p)m_pv_p^2\right],
\end{eqnarray}
where use has been made of the form of the equilibrium distribution function and the fact that
\begin{eqnarray}
  && \int d\mathbf rdx_pd\mathbf vdv_p\ (f_{i,eq}f_{p,eq}-f_if_p)=0, \quad i=1,2, \\
  && \int dx_p dv_p\ (f_{p,eq}-f_p)=0.
\end{eqnarray}
Hence, Eq. \eqref{eq:3t} is
\begin{equation}
  -\frac{\beta}{2}\left(E_{eq}-E\right)=0,
\end{equation}
where $E_{eq}$ is the mean energy at thermal equilibrium, with is the same as the mean energy $E$ for any other state, since collisions conserve energy.

This way, $\mathcal H\ge \mathcal H_0$.

\section{Discussion and conclusions \label{sec:5}}
In this work, we have proposed a mesoscopic description of a system made of a cylinder including a finite number of hard spheres in $d$ dimensions divided by the adiabatic piston. We started by computing some equilibrium properties of the system as given by the canonical distribution. The analysis showed that the fluctuations of the position of the piston make the probability distribution of the gases to be spatially inhomogenous. However, a direct calculation also showed that the conditional distributions of the gases to a position of the piston are spatially homogenous, which provide a direct characterization of the global equilibrium. Furthermore, it can be shown that the conditional distributions of the gases to a position and a velocity of the piston depend on the position but not on the velocity, which reflect the presence of spatial correlation as well as the absence of velocity correlations at thermal equilibrium.

The derivation of the kinetic equations depend on few assumptions. As it is usual in Kinetic Theory, we have used the molecular chaos hypothesis by removing velocity correlations when two particles or a particle and the piston are about to collide. In addition, we have also remove some velocity correlations between particles and the piston after a collision. As already discussed, this is expected to be a good approximation if the piston suffers from many collisions in a mean free time of the surrounding particles, when the gases are dilute enough. As a result of the approximations, we ended up with a closed system of equations which solved the adiabatic piston problem, that is to say describes the evolution of the system towards thermal equilibrium where the two gases have the same pressure and all temperatures, including that of the piston, are equal. As a main difference with respect to other theories one may find in the literature, the one here is not restricted to the thermodynamic limit, nor to small particle to piston masses ratio, and include short-ranged collisions between particles.

Many of the results of the present work are easily generalized along different directions. Even though the kinetic description is for a system of hard spheres, more general interactions among particle can be considered, by replacing the collision operators $J_{11}$ and $J_{22}$ with the appropriate ones. Moreover, for most the of the collision operators conserving the number of particles, linear momentum, and energy, the $\mathcal H$-theorem still holds. More in general, we can also include dissipation, thermostats, vibrating walls, more pistons as in \cite{br12,cacesavu17}, and so on.

An interesting aspect not addressed in this work is to precise the conditions under which the assumptions are expected to be correct. This would need a formulation including velocity correlations, which are seems to be present for small system even at thermal equilibrium \cite{cegrsavivu14}. This aspect is left for future works.

\begin{acknowledgements}
I dedicate this work to the memory of Mar\'ia Jos\'e Ruiz Montero.
\end{acknowledgements}

\bibliographystyle{spphys}


\section{Introduction}
An important early result in quantum mechanics is the so-called avoided crossing, also anti-crossing theorem.
It dates back to the work of J. von Neumann, and E. Wigner~\cite{vonNeumann1993} -- more accurately,  two years  earlier, to the work of F. Hund~\cite{Hund1927};
J. von Neumann and E. Wigner   proved
 and generalized
the intuition of F. Hund.
Their finding  is as follows  (see also \cite{Suzuki2013, cmu2}):

~~\\
Suppose the
Hermitian matrix~$H$ depends on a number of parameters $\kappa_1, \kappa_2, \cdots$. Then {\it when one can change only
one or two parameters it is in general not possible to reach an intersection of the eigenvalues of $H$}. This number goes to one if $H$ is a real Hermitian matrix; that is,
{\it in
this case it suffices to be able to change two real parameters to make two eigenvalues degenerate. }

~~\\
The role of anti-crossings in adiabatic quantum computations is well known~\cite{Farhi472, Vazirani}.
This is expressed as the total adiabatic evolution time $T$ being inversely proportional to the square of the energy difference between the two lowest energies of
the given Hamiltonian
\begin{equation}\label{H}
    H(s) = (1-s) H_{initial}  + s H_{final},
\end{equation}
for $0\leq s=t/T\leq 1$.
 This energy gap, for many instances\footnote{By an instance, we mean  a time dependent Hamiltonian; thus, we  include the choice of the initial Hamiltonian in our definition--and not just the problem Hamiltonian.}, is very small,  and the enlargement of this gap is the  focal point in current research in the area of adiabatic quantum computations. The goal of our paper is to present a general principled methodology to study anti-crossings that generalizes previous studies (that are ad-hoc) and also provides guidelines for new constructions.

\subsection{Connection to Morse functions} Let us start by  making the connection with Morse theory \cite{PMIHES_1988__68__99_0}.
Suppose $\lambda_1(s)$ and $\lambda_2(s)$ are two parallel planar curves; i.e., their tangent vectors span  a two dimensional vector space  (e.g., two eigenvalues of the time dependent Hamiltonian $H(s)$). Let us assume that $\lambda_1(s)$ and $\lambda_2(s)$ are twice differentiable on a compact interval $I\subset \mathbb [0,1]$ and satisfying $\lambda_1(s)<\lambda_2(s)$ for all $s$ in $I$. Suppose at some point $s^* \in I$, the two curves  converge to each other (without intersecting) and then diverge again; that is, $\lambda_1(s)$ is convex on $I$ while $\lambda_1(s)$ is concave. Without  loss of generality, we can write
\begin{equation}
\lambda_1(s) = s^2 +\varepsilon \mbox{ and }\lambda_2(s) = -s^2-\varepsilon,
\end{equation}
 with $\varepsilon>0$.
The spectral gap $\Delta(s)=\lambda_2(s)-\lambda_1(s)=2s^2+2\varepsilon$
is minimized at the anti-crossing $s^*=0$; see Figure 1. Now, we embed the two planar curves into the three dimensional Euclidean space
$\mathbb E^3$ and define the function
$$
    f(s, \lambda) = (\lambda - \lambda_1(s))  (\lambda - \lambda_2(s)).
$$
The graph of the function $f$ is the saddle surface (also depicted in Figure 1). With this embedding, the anti-crossing point $s^*$ is mapped into a non degenerate critical point $(s^*, \lambda^*)$ of $f$ with a negative  Gaussian curvature \footnote{Because $f\in \mathbb Q[s, \lambda]$, a point $p$
in $M$ is a non degenerate critical point of the function $f$ if and only if the remainder  ${\sf NF}_\mathcal B (Hessian(f))$  is not zero at $p$, where $\mathcal B$
is a Groebner basis for the ideal generated by the two polynomials $f_s$ and $f_\lambda$.}.
As we shall show, this embedding also constrains {\em what} can be done to the picture on the right in order to eliminate the anti-crossings and provides  guidelines to {\em how} it can be done. The key constructions are gradient flows \cite{Hurtubise}, the  Euler characteristic, and the Gauss-Bonnet theorem \cite{doCarmo} which distributes the curvature consistently with  the Euler characteristic.
\begin{figure}[h]
    \centering
    \subfloat
    {{\includegraphics[width=3cm]{pics/pic1.png} }}%
    \qquad
        \centering
    \subfloat
    {{\includegraphics[width=3cm]{pics/pic2.png} }}%
    \caption{{\small  Anti-crossing happens at saddle points. Note that the picture on the left can be recovered by intersecting the saddle surface (right) with
the horizontal plane $z=0$. The resulting curve is commonly referred to as Dupin indicatrix.}}
    \label{Chimera}
  \end{figure}

\subsection{Outline of the paper}
In section  2, we describe the general methodology for changing the dynamics around anti-crossing. We start by reviewing the connection with Morse theory and add more details than the simple connection we made in the Introduction. We then explain that the central object of our investigation
is one (or more) parameter families of Hamiltonians and review the appropriate mathematical machinery, i.e.,  Cerf theory.   The third section discusses
two concrete illustrations of the general theory; these two situations are taken from published literature. The first situation repels the ground state from the rest of the spectrum while the second situation weakens first order quantum phase transitions (QPT) into a higher order. The role of the Gauss-Bonnet theorem is elucidated.   We conclude with some suggestions for future research.

\section{The general theory}
We briefly review the Morse theoretical depiction of the adiabatic evolution of the Hamiltonian (\ref{H}) that we have introduced in \cite{cmu2}. We then extend this depiction to
parametrized families of Hamiltonians.  For Morse theorists, this extension is referred to as Cerf theory.
\subsection{Morse theory}
To the Hamiltonian (\ref{H}), we assign the Morse function $f: M\rightarrow [0,1],$ given by the characteristics polynomial (aka secular function)
\begin{equation}\label{detH}
    f(s, \lambda) = det(H(s)-\lambda I),
\end{equation}
where $I$ is the identity $n\times n-$matrix. The surface $M$, aka cobordism in the language of Morse theory, is defined by the (non degenerate) critical points
of $f$ using handlebodies decomposition  \cite{cmu2}. This   interplay  between $f$ and the topology (homology) of the cobordism $M$ is the apex of Morse theory, and can be expressed in more precise terms as follows: the trajectories of the gradient flow
 \begin{system}\label{gradientFlowEnergies}
	\frac{d}{d \tau} \lambda(\tau) &=& - \partial_\lambda f (\lambda(\tau) s(\tau)),\\[3mm]
	\frac{d}{d \tau} s(\tau) &=& - \partial_s f (\lambda(\tau), s(\tau)),
\end{system}%
triangulate $M$ and the homology of this complex is isomorphic to the ``usual" singular homology of~$M$.  One can then work on the cobordism $M$ to investigate $f$,
or vice versa. This paper is aligned with the former.

~~\\
If we Taylor series
$f$ around each of its critical points (and rotate axes to coincide with principal curvature directions), we get a simple local expression of $f$:
\begin{lemma}[Morse lemma - Revisited]\label{MorseLemma}
Given $p=(s_p, \lambda_p)\in cr(f),$ i.e., a non degenerate critical point of the smooth function $f: M\rightarrow [0,1]$. We have
\begin{equation}\label{MorseEq}
	f(s, \lambda) = f(p) + K_1(p) (s-s_p)^2 + K_2(p) (\lambda-\lambda_p)^2 + higher\, order\, terms,
\end{equation}
with $K_i(p)$ being the principal curvature of $M$ at $p$.
\end{lemma}
The lemma, which we use repeatedly, gives the two energies around the critical point $p$ as  solutions of the equation (\ref{MorseEq}). In particular, the energy difference between the two
 is given by
\begin{equation}\label{Delta}
	\Delta(s) = 2\,{\frac {\sqrt {- K_2(p) \left( {\it f(p)}+ K_1(p){s}^{2}
 \right) }}{{K_2(p)}}}
\end{equation}
(after translating $p$ to the origin--it is important to notice that translations, as well as rotations- above,  are isometries, i.e., distance preserving transformations). This gives
\begin{equation}\label{minDelta}
	\mathrm{min}_s\Delta(s) = 2\,{\frac {\sqrt {- K_2(p) \left( {\it f(p)}
 \right) }}{{K_2(p)}}}.
\end{equation}
The square root is always real-valued: if $p$ is a saddle, i.e., an anti-crossing, then $K_1(p)<0<K_2(p)$ and $f(p)<0$. If
$p$ is a maximum then $K_1(p)<0>K_2(p)$ and $f(p)>0$. If $p$ is a minimum then $K_1(p)>0<K_2(p)$ and $f(p)<0$.

~~\\ Remarkably, the different local descriptions of $f$ can be glued together. The fact that these local depictions can be glued together is the incarnation of the globality of Morse theory. In the language of {\it functors}, this can be presented  as follows: Let $ \mathcal O_M$ denote the category of standard open sets of $M$
sorted by inclusion and $C^\infty(M)$ denote the category of real-valued smooth functions sorted by restrictions.
\begin{proposition}
The functor $\mathcal F_M:  \mathcal O_M \rightarrow C^\infty(M)$, defined by the lemma for each open set, is a sheaf.
\end{proposition}
For readers familiar with this language, the proof is almost evident, particularly when  we take into account that critical points are isolated.

\subsection{Cerf theory}
As mentioned before, the central object in this paper is not the mere Hamiltonian (\ref{H}) but  the one-parameter families of Hamiltonians
\begin{equation}\label{Hb}
    H^b(s) = H(s) + b H_{enhancement}.
\end{equation}
This gives rise to a one-parameter family of Morse functions $f^b: M^b\rightarrow [0,1]$ and to the homotopy:
\begin{equation}\label{h}
h:M^0\times [0, 1] \rightarrow M^1,
\end{equation}
with $h(s, \lambda, 0) = f^0(s, \lambda)$
and $h(s, \lambda, 1) = f^1(s, \lambda)$.
An important implication of  the homotopy (\ref{h}) is that the  Euler characteristic is given independently of $b$ by the Gauss-Bonnet formula:
\begin{equation}
    \int_{m \in M^b} K^b d m  = 2\pi \chi.
\end{equation}
Because the curvature is essentially dumped at the critical points, the above formula constrains  the quantities $\{K^b(p^b), p^b\in cr(f^b)\}$, and,  by (\ref{Delta}), the energy differences, simultaneously.  In fact, for each critical point $p^b\in cr(f^b)$, we have
\begin{equation}
	\mathrm{min}_s\Delta_{p^b}(s) = 2\,{\frac {\sqrt {- K_2^ b(p^b) \left( {\it f^b(p^b)}
 \right) }}{{K_2^ b}}}
\end{equation}
while at the same time we also have
\begin{equation}
	\sum_{p^b\in cr(f^b)}  K_{1}^b(p^b)K_{2}^ b (p^b)  \int _{V^b(p^b)} dm= 2\pi \chi,
\end{equation}
where  the neighbourhoods $V^b(p^b)$ are the graph of $f^b$ (given by the  lemma).  Equivalently, we are constrained (in our quest of enhancement)  by the commutativity of the following diagrams:

 \begin{center}
\begin{tikzcd}
f^0 \arrow[to=Z, "{h(- ,b)}"] \arrow[to=2-2]
& {} \\
& \qquad  \qquad    H_*(M^b, \mathbb Z)=H_*(M^0, \mathbb Z)\\
|[alias=Z]|  f^b  \arrow[to=2-2]
\end{tikzcd}
\end{center}
where $f^b\in \Gamma \mathcal{F}_{M_b}$ is the unique global section of the functor $\mathcal{F}_{M_b}$.
In particular, as $b$ changes, the set of critical points changes but consistently with the invariant homology groups. This means that critical points might move around and  new critical points
might appear (born), or old ones disappear (death), by pairs of (saddle, minimum) or  (saddle, maximum), which amounts to adding, or deleting, a cylinder (surface with $\chi=0$), thus not affecting the topology.

\section{Two illustrative examples}
We give two illustrations of the general methodology described above, i.e, homotopically altering the dynamics around the anti-crossing by adding critical points to its vicinity   (through the process of birth/death of critical pairs).  In the first example, we push away the first excited state from the ground state by adding degeneracy, whilst in the second, we change a first order quantum phase transition (QPT) into a second order by adding {\it non-stoquasticity}. The role of the Gauss-Bonnet theorem is also exhibited in both cases.

\subsection{Enhancement by degeneracy: Repelling the excited state}
This section revisits the work of Dixon \cite{Dickson1}.

~\\
The main idea is described in Figure \ref{pantsAndCaps}. Initially, the two lowest energies are attracted to each other due to the presence  of an anti-crossing (as in (a)).
The addition of a minimum  (or maximum as depicted, by making the ground state degenerate at $s=1$, will repel the rest of the spectrum from the ground state, simply because level sets around a minimum  (or maximum) are ellipses ((b) and (c)). See also Figure \ref{dickson}.

   \begin{figure}[h]
    \centering
    \subfloat  []
    {{\includegraphics[width=5cm]{pics/SaddleOnly.png} }
        \centering
    \subfloat  [ ]
    {{\includegraphics[width=5cm]{pics/SaddleAndMax.png} } }%
      \subfloat  [ ]
    {{\includegraphics[width=4cm]{pics//SaddleAndMax2D.png} } }%
    \caption{ (a): The initial spectrum lies on  the so-called pair of pants. The two energies are dragged closer by the saddle point i.e., anti-crossing. (b) and (c):
    The maximum repels the ground state from the rest of the spectrum.  }
    \label{pantsAndCaps}
  \end{figure}

   \begin{figure}[h]
    \centering
    \subfloat  []
    {{\includegraphics[width=4cm]{pics/dickson1.png} }
        \centering
    \subfloat  [ ]
    {{\includegraphics[width=4cm]{pics/dickson2.png} } }%
    \caption{ Two examples that use  the explicit deformation given in  \cite{Dickson1}.  For instance, (b) is obtained by replacing the problem Hamiltonian
    $H_{final} = -\sigma^1$  with $H_{final} + b/2 (\sigma^1+1 ) (1-\sigma^2+1) + b/2 (1-\sigma^1 ) (\sigma^3+1)$.
      }
    \label{dickson}
  \end{figure}

~\\
To make this more concrete,  consider the following situation where we initially have (Morse function of the Grover search adiabatic Hamiltonian \cite{cerf, cmu2})
\begin{equation}
	f(s, \lambda) =\left( {2}^{-N}-1 \right)  \left( {s}^{2} -s\right) +{\lambda}^{2}-
\lambda =
 \left( \lambda- \lambda_2(s) \right) \times \left( \lambda-\lambda_1(s) \right).
\end{equation}

~~\\
Similar to \cite{Dickson1}, we make the ground state  degenerate at $s=1$ by defining the one parameter family of functions  ($b>0$):
\begin{equation}
	f^b(s, \lambda) = \left( \lambda- \lambda_2(s)-b) \right) \times \left( \lambda-\lambda_1(s) \right) \times \left( \lambda-  b) \lambda_1(s) -\lambda_1(s)\right).
\end{equation}
Here $f^0(s, \lambda) = f(s, \lambda)  \times  \left( \lambda-\lambda_1(s) \right).$ Let us, for instance, fix $b$ to 1 and $N$ to~5.   The Morse function $f^1(s, \lambda)$ now has two critical points in the rectangular region $[1/2-\epsilon/2, 1/2+\epsilon/2] \times [0, R]$ where $R$ is a large number (see Figure \ref{mexican}):
\begin{equation}
	p^1= (s = 1/2, \lambda = 0.78),\, q^1= (s = 1/2, \lambda = 2.05).
\end{equation}

~~\\
Figures  \ref{mexican} and \ref{twoSpec} also show that $p^1$ is  a maximum while $q^1$ is an anti-crossing point (a saddle).
When $\epsilon$ is infinitesimally small (it makes sense to take it infinitesimally small, because the addition of $a$ is   affecting ony the $\lambda$ direction; hence, we would like to concentrate on this direction), the two principal curvatures are essentially the same, i.e., $K_2$ (recall that the principal curvatures are by definition the maximum and minimum  of the curvatures of all curves passing through the given point, and because the surface shrinks to a curve, the two are equal).  The application of the Gauss Bonnet formula on  the region $M\cap [1/2-\epsilon/2, 1/2+\epsilon/2] \times [0, R]$,
with $M$ being the surface depicted in (b) above, gives
\begin{equation}
	K_2(p^1)(R-\frac{2}{3}R^3 K_2(p^1) ) + K_2(q^1)(R-\frac{2}{3}c^3K_2(q^1) )  +  O(\epsilon) = 0,
\end{equation}
which gives the functional dependence between $K_2(p^1)$ and $K_2(q^1)$.
Returning to our instance, a direct calculation of the two principal curvatures confirms this dependence:
$
	 K_2(p_1) = - K_2(q_1) = -3.81.
$

  \begin{figure}[h]
    \centering

    {{\includegraphics[width=2cm]{pics/mexican.png} }
    \caption{Gauss-Bonnet theorem constrains the energy differences.}
    \label{mexican}
  \end{figure}

 \begin{figure}[h]
    \centering

    {{\includegraphics[width=4cm]{pics/dickson0.png} }
    \caption{The spectrum of $f^0$ is in purple. The spectrum  of $f^1$ is in green; the ground state is degenerate at $s=1$. Both spectra share the same ground state. The green spectrum  has three energies. The first two energies are separated by a maximum $p_1= (s = 1/2, \lambda = 0.78)$; figuratively, because level sets around $p_1$ are ellipses, the two energies repel each other,  and this keeps the ground state away from the rest of the spectrum. }
    \label{twoSpec}
  \end{figure}

\subsection{From first order QPT to second order}
When the system endures a first order QPT, the spectral gap between the two lowest energies is exponentially small, yielding exponentially  large evolution time. The reason is quite simple. In  a first order QPT, the two local minima that the ground state occupies before and after the transition are too far apart (in the Hamming distance sense), so the tunnelling is exponentially  slow.   This singularity is expressed mathematically by the discontinuity of the first derivative of  the free energy.  A higher order QPT comes with less severe shrinking of the spectral gap.

~~\\
The bare bone idea of this section is depicted in the picture below. Initially, we have an anti-crossing in the region $M^0\cap B$. The angle
$\theta $
between the two asymptotes is given by
\begin{equation}
\theta= \pi -\arccos \left( {\frac {K_1+K_2}{K_1-K_2}} \right) = \pi-\arccos \left( {\frac {m}{\sqrt {{m}^{2}-4\,K}}} \right).
 \end{equation}
Here $m$ is the mean curvature.   For first order QPTs,  this angle is large, i.e., approaches
$\pi$. The figure on the right shows the birth of a pair (min, saddle)  in the vicinity of the anti-crossing.  The effect of this homotopical  change is that the new angle
$\theta$ is reduced.    As before, the curvatures of the three critical points are related with the relation
$
	 \int K = -2\pi.
$

 \begin{figure}[h]
    \centering
    \subfloat  []
    {{\includegraphics[width=4cm]{pics/1stOrderQPT.png} }

    \subfloat  []
    {{\includegraphics[width=6cm]{pics/2ndOrderQPT.png} }
    \caption{A first order QPT becomes second order by adding the pair (min, saddle). The angle between the two asymptotes, which was initially large (a), is now reduced (b).}
    \label{pattern}
  \end{figure}

~~\\
Let us make this more concrete by considering the following, somewhat simplified, scenario:
 \begin{figure}[h]
    \centering

    {{\includegraphics[width=9cm]{pics/someWhatSimple.png} }
    \caption{Illustrative scenario.}
    \label{charts}
  \end{figure}
 ~~\\
The Morse function is given by
\begin{equation}
	f^\epsilon(s, \lambda) =  (\lambda-\lambda^*(s)) (\lambda + \lambda^*(s))
\end{equation}
with
\begin{equation}
	\lambda^*(s) = \left( \frac{2}{3}\,{\frac {{s}^{2}}{{\epsilon}^{2}}}-\frac{1}{6}\,{\frac {{s}^{4}}{{
\epsilon}^{4}}} \right)  \left( a_{{1}}-2\,a_{{0}}+a_{{-1}} \right) +
 \left( \frac{1}{6}\,{\frac {{s}^{3}}{{\epsilon}^{3}}}-\frac{2}{3}\,{\frac {s}{
\epsilon}} \right)  \left( a_{{-1}}-a_{{1}} \right) +a_{{0}}
\end{equation}
and $a_{i}=1/{K^\epsilon(p_i)}$ for $i=-1,0,1$. Fix $B= [-2\epsilon,2\epsilon] \times [-2a_0, 2a_0]\times [-R, R]$,
with $R$ a large enough positive number-- precisely,  ${R\geq max(\{|f^\epsilon(p_i)|,\, i=-1,0,1  \})}$. Working out
Gauss-Bonnet formula on $M\cap B$ yields the dependancy between the different curvatures:
\begin{eqnarray}\nonumber
	-2\pi &=& {\frac {32}{2835}}\,a_{{0}} \left( -2929\,{a_{{0
}}}^{2}+896\,{a_{{-1}}}^{2}+896\,{a_{{1}}}^{2}-32\,a_{{1}}a_{{0}}-32\,
a_{{0}}a_{{-1}}+256\,a_{{1}}a_{{-1}} \right) \epsilon\\\nonumber
&+& {\frac {16}{2835
}}\, \left( -10489\,{a_{{0}}}^{2}+896\,{a_{{1}}}^{2}-32\,a_{{1}}a_{{0}
}+256\,a_{{1}}a_{{-1}}-32\,a_{{0}}a_{{-1}}+896\,{a_{{-1}}}^{2}
 \right) \epsilon\,R\\
 & -& \frac{16}{3}\,{a_{{0}}}^{3}R.
\end{eqnarray}
A better example would be  to have two different parameters instead of just $\epsilon$, so the distance between the minimum  and the saddle is independent of their distance from the original saddle. The computation generalizes easily but gives complicated expressions.

\subsubsection{The case of the ferromagnetic $p-$spin model }
Here we pay a particular attention to the very instructive model of the  ferromagnetic $p-$spin. See the work by Nishimori and Takada \cite{10.3389/fict.2017.00002}. The time dependent Hamiltonian
is given simply by
\begin{equation}
	H(s) = -s N\left ( \frac{1}{N} \sum_{i=1}^N \sigma_i^z\right)^p - (1-s) \sum_{i=1}^N \sigma_i^x,
\end{equation}
where the $p-$th power refers to the matrix power; i.e., composition of linear operators on the Hilbert space ${\mathbb C^2}^{\otimes N}$. A comprehensive
investigation into the Hamiltonian $H(s)$ can be found in \cite{Bapst_2012}. This Hamiltonian exhibits a first order QPT for $p\geq 3$;  this can be shown in various ways--for instance,  by taking the thermodynamical limit $N=\infty$, in which case one gets  a classical (continuous) spin of a single particle. Figure \ref{thermo} shows
the continuity and discontinuity of the first derivative of the classical energy $e(p,s)$ when $p=2$ and $p=5$ respectively.

 \begin{figure}[h]
    \centering
    \subfloat  []
    {{\includegraphics[width=4cm]{pics/ep2.png} }
        \centering
    \subfloat  [ ]
    {{\includegraphics[width=4cm]{pics/ep5.png} } }%
    \\
      \subfloat  [ ]
    {{\includegraphics[width=4cm]{pics/ep5nonStoquastic.png} } }%
    \caption{ (a): The minimum of the energy $e(p=2)$ is continuous function of $s$; this minimum is  given by $\theta=\pi -\arccos \left({\frac {-1+s}{2s}} \right)$. The phase transition is second order. (b):
    The minimum of the energy $e(p=5)$ is discontinuous at $s=0.38$. The phase transition is first order. (c) After the addition of the non-stoquastic term with $b =0.1$, the  minimum of the energy is now a continuous function of~$s$. }
    \label{thermo}
  \end{figure}
 ~\\
For finite $N$, the system is also reduced into a much smaller space. This is because of the commutativity of $H(s)$ with the total spin
$S$. Interestingly, this reduction can be generalized using the theory of irreducible representations \cite{cmurep}, but here we recall the usual reduction.  Conveniently, the
 lowest part of the spectrum of the total operator is also the lowest part of the spectrum of $H(s=0,b)$. And because of the anti-crossing theorem, this is also true for $H(s, b)$. This lowest part is given by the states with maximum spin $N/2$ i.e.,  states  $|N/2, m\rangle$, for
 $0\leq m\leq N$,  and defined with
\begin{eqnarray}
	S_z |N/2, m\rangle  &=& m |N/2, m\rangle\\
	S^2 |N/2, m\rangle  &=& N/2 (N/2+1) |N/2, m\rangle.
\end{eqnarray}
In addition we have
\begin{equation}
	S_{\pm} |N/2, m\rangle =  \sqrt{N/2(N/2+1) - m(m\pm1)  } \,  |N/2, m \pm 1\rangle
\end{equation}
where $S_{\pm} =S_x \pm S_y$ are the spin raising the lowering operators. All of this  continues to hold with
the new non-stoquastic Hamiltonian  (\cite{10.3389/fict.2017.00002}):
\begin{equation}
	H(s, b) =  -s b N\left ( \frac{1}{N} \sum_{i=1}^N \sigma_i^z\right)^p
	+ s(1-b) N  \left( \frac{1}{N} \sum_{i=1}^N \sigma_i^x\right)^k
	 - (1-s) \sum_{i=1}^N \sigma_i^x.
\end{equation}
A direct calculation gives the matrix elements:
\begin{eqnarray}\nonumber
	H(s,b)_{m, m} &=&  s \left( -bN \left( 1-2\,{\frac {m-1}{N}} \right) ^{p}+ \left( 1-b
 \right)  \left( 2\,m-1-2\,{\frac { \left( m-1 \right) ^{2}}{N}}
 \right)  \right)
\\[3mm] \nonumber
H(s,b)_{m, m+1} &=& H(s,b)_{m+1, m}  = - \left( 1-s \right) \sqrt { \left( N-m+1 \right) m}\\[3mm] \nonumber
H \left( s,b \right) _{{m,m+2}}&=& H(s,b)_{m+2, m} ={\frac {s \left( 1-b \right) \sqrt {m
 \left( m+1 \right)  \left( N-m+1 \right)  \left( N-m \right) }}{N}}.\nonumber
\end{eqnarray}
Other matrix elements are zero.
Sending the non-stoquasticity parameter $b$ to zero amounts to the homotopical reduction of the angle $\theta$ described in the beginning.  This
pattern   is reproduced in Figure \ref{Nicespec} with the two lowest energies; see also Figure \ref{angles}. Figure \ref{dance} shows how critical points are born and dead
as $b$ changes without affecting the Euler characteristic. As the parameter  $b$ decreases towards zero, the number of critical points increases, forcing the total curvature to redistribute itself on a larger set of critical
points.~\begin{figure}[h]
    \centering

    {{\includegraphics[width=3cm]{pics/angles.png} }
    \caption{ Decrease of the angle between the two asymptotes around the anti-crossing for $N=7$ and $p=5$. In red, the asymptotes at $b=1$. In green, asymptotes at $b=0.1$.  }
    \label{angles}
  \end{figure}

  \newpage

 \begin{figure}[H]
    \centering
    \subfloat  []
    {{\includegraphics[width=11cm]{pics/specN7p5mu1.png} }

    \quad
    \subfloat  []
    {{\includegraphics[width=14cm, height=7cm]{pics/specN7p5mu01bis.png} }
    \caption{Spectrum of $H(s, b)$ for $N=7$ and $p=5$. (a) $b =1$ (b)  $b =0.1$. The pattern of Figure \ref{pattern} is reproduced with the two lowest energies (at the left most saddle).}
    \label{Nicespec}
  \end{figure}

  \newpage

\begin{figure}[H]
    \centering
    \subfloat  []
    {{\includegraphics[width=8cm, height=6cm]{pics/N7p5.png} }
        \centering
    \subfloat  [ ]
    {{\includegraphics[width=7cm,  height=6cm]{pics/N7p3.png} } }%
    ~~\\
      \subfloat  [ ]
    {{\includegraphics[width=5cm]{pics/N7p5nberofcpts.png} } }%
    \qquad    \quad \quad \,
      \subfloat  [ ]
    {{\includegraphics[width=5cm]{pics/N7p3nberofcpts.png} } }%
    \caption{ (a): The evolution of the critical points for $p=5$ and $N=7$. In red are the saddle points (anti-crossings). In green, are the maxima, and in blue, the minima.  Euler characteristic  is constant  $\chi=-7$ for all $b$ in $\mathbb R$. (b): Same for $p=3.$  (c) and(d):
    The number of critical points (with $N=7$) for $p=5$ and $p=3$ respectively. This number is minimized at $b =1$ and maximized at $b =0$.}
    \label{dance}
  \end{figure}

 \newpage
 \newpage

\section{Conclusion}
The paper presents a general methodology for enhancing the computational power of adiabatic quantum computations. The central idea is to homotopically deform the dynamics around the anti-crossing to weaken its effect.  We have also illustrated the key role the Gauss-Bonnet theorem plays in this process. As illustrated by the two examples, the modifications provide speedup either by repelling the excited state or by transitioning from a first order QPT to the second, and so provide guidance on how to design enhancements. One interesting future direction is to perform such a deformation using degenerate critical points (hence, Conley theory \cite{cmu2}),  as such points are no longer constrained by  the Gauss-Bonnet theorem. This gives more freedom on how to control the evolution.

\bibliographystyle{plain}