\section{Introduction}\label{intro}
\IEEEPARstart{W}{avefield} migration is the art of computing the medium reflectivity from the recorded wavefield passing through the medium. Many migration algorithms are available, such as Kirchhoff migration \cite{Claerbout}, one-way wave equation migration \cite{Mulder}, and Reverse-Time Migration (RTM) \cite{BaysalRTM,Zhou,Zhang}. RTM is one of the common migration algorithms among others. RTM is commonly implemented by applying the adjoint of the Born operator to the recorded data \cite{SeisInv}. However, the inverse of the operator is needed for the true image. Consequently, the migration result suffers from amplitude and resolution issues \cite{Dutta,Tang,Liu}. To address this problem, researchers solve the migration problem with a least-squares solution called Least-Squares Reverse-Time Migration (LSRTM). The least-squares solution is usually applied as an iterative optimization problem.

However, LSRTM is computationally expensive \cite{Dai,Tang,Herrmann,Farshad}. To reduce the computation cost of the LSRTM, one can reduce the computation domain by confining the model to a Region of Interest (ROI) by recording or computing the wavefields at the boundary of this region. The process of migrating for a ROI is called target-oriented migration \cite{Valenciano,Haffinger,Willemsen,Yuan,Zhao,Guo2,Ravasi2}. However, deploying receivers on the boundaries of ROI usually is not possible due to physical obstacles.

Due to the aforementioned accessibility issue, the typical approach for target-oriented migration is to opt for redatuming algorithms and only consider the upper horizontal boundary of the ROI or target \cite{Guo1,Guo2,Wapenaar6,Vargas,Schuster,Ravasi1,Luiken,Shoja3}. The shortcoming of only considering the upper boundary is that any wavefield entering the ROI from the medium below the ROI is unaccounted for, hindering the convergence of the inversion process. Moreover, including the lower boundary in the algorithm can add transmission information to the inversion. However, including the lower boundary in the inversion process has rarely been studied directly. For instance, Cui et al. \cite{Cui} derive a representation with a reciprocity theorem and Marchenko redatuming to include surrounding boundaries in target-oriented full waveform inversion (FWI), Diekmann et al. \cite{Diekman1} use a Marchenko retrieved Green's function of the ROI and insert it in the Lippmann-Schwinger integral to create a linear inversion process, and van der Neut et al. \cite{Joost1} design a target-enclosed imaging algorithm with the help of a reciprocity theorem. Of the above-mentioned papers, only \cite{Joost1} directly studies the consequences of including the lower boundary in the imaging process, and the others only implicitly imply the effects of it.

This paper, which is an extension of \cite{Joost1}, studies the contribution of the lower boundary by introducing a target-enclosed LSRTM algorithm. To derive this algorithm, we start by explaining the LSRTM briefly. Next, we derive a target-enclosed representation for Green's functions on the upper and lower boundaries of the ROI by using the reciprocity theorem. Then, we connect this target-enclosed representation to LSRTM to complete our algorithm. After deriving the required equations, we test our algorithm with numerical examples. First, we use physical receivers at the boundaries of the ROI to check the results in the ideal situation. Then, we briefly introduce Marchenko redatuming to explore the possibility of using virtual receivers in our algorithm. For both physical and virtual receivers cases we use both a homogeneous background and a smooth background velocity as the velocity model for migration. Finally, we finish the paper by discussing the results and providing a conclusion.

\section{\label{sec:2} Theory}
To develop the theory, we start with a brief discussion of LSRTM. Next, a representation of the target-enclosed Green's function is given with the help of the reciprocity theorem. Finally, to derive our target-enclosed LSRTM formulation, we combine LSRTM with the target-enclosed representations. In the entire theory section, we are in the frequency-space domain, and for simplicity, we drop the dependency on angular frequency ($\omega$).

\subsection{Least-Squares Reverse-Time Migration}
We start the explanation of LSRTM by investigating the Born integral for the scattered wavefield by a scattering potential \cite{born_wolf_1999,SeisInv,Peter}:
  \begin{equation}
  P_{sct}(\textbf{x}')=\int_{\mathcal{V}}\gamma_0^2(\textbf{x})G_0(\textbf{x}',\textbf{x})\chi^c(\textbf{x})P_{inc}(\textbf{x})d\textbf{x}.
  \label{scattering}
  \end{equation}
In this equation, $\textbf{x}'$ is the observation location, $\textbf{x}$ is a location inside the computation volume ($\mathcal{V}$), $P_{sct}(\textbf{x}')$ is the scattered pressure field at the observation point, $P_{inc}(\textbf{x})$ is the incident pressure field at the computation point and $G_0(\textbf{x}',\textbf{x})$ is the background Green's function between $\textbf{x}$ and $\textbf{x}'$. Moreover, $\gamma_0(\textbf{x})=\frac{\omega}{c_0(\textbf{x})}$ and $\chi^c(\textbf{x})=1-\frac{c^2_0(\textbf{x})}{c^2_{sct}(\textbf{x})}$ is the propagation velocity perturbation, where $c_0(\textbf{x})$ and $c_{sct}(\textbf{x})$ are the background and scatterer's velocity, respectively.

It is possible to rewrite this equation in matrix form:
  \begin{equation}
    \textbf{P}^{pred}_{sct}(\delta \textbf{m})=\textbf{L} \delta \textbf{m}.
    \label{scatt_mat}
  \end{equation}
Here, $\textbf{L}$ is matrix form of the integral operator of equation \ref{scattering} and $\delta\textbf{m}$ is a vector, containing the perturbation $\chi^c(\textbf{x})$.

To obtain an estimation of $\delta\textbf{m}$, we can apply the adjoint of $\textbf{L}$ to the observed data:
  \begin{equation}
    \delta \textbf{m}^{img}=\textbf{L}^\dagger \textbf{P}^{obs}_{sct},
    \label{imaging}
  \end{equation}
  where $\dagger$ denotes complex conjugate transposition.

 We can go one step further to obtain a high-resolution estimation of the image by minimizing the following objective function \cite{Marquardt}:
  \begin{equation}
    J(\delta \textbf{m})=\frac{1}{2}\|\textbf{P}^{pred}_{sct}(\delta \textbf{m})-\textbf{P}^{obs}_{sct}\|^2_2.
    \label{costfunc}
  \end{equation}
  Different optimization algorithms, such as conjugate gradient, can minimize this function. This optimization problem is known as Least-Squares Reverse-Time Migration (LSRTM).

  \subsection{Target-enclosed representations}\label{TE}
  LSRTM is a computationally expensive algorithm. In order to reduce its computational burden, seismologists usually opt for a target-oriented algorithm, limiting the medium to a smaller region. As we mentioned in section \ref{intro}, target-oriented algorithms usually redatum the data to the upper boundary of the target and ignore any information coming from the lower boundary.

  In this section, the idea is to find a representation that can account for a heterogeneous medium above the upper and below the lower boundary of the target area. To derive this representation, we follow \cite{Joost1}. The starting point of the derivation is the acoustic reciprocity theorem of the convolution type \cite{Wapenaar6}, which connects the wavefields of two different states via
  \begin{equation}
  \int_{\mathfrak{d}\mathcal{V}_u}\rho^{-1}(p^+_A(\partial_3p^-_B)+p^-_A(\partial_3p^+_B))d\textbf{x}= \\
  \int_{\mathfrak{d}\mathcal{V}_l}\rho^{-1}(p^+_A(\partial_3p^-_B)+p^-_A(\partial_3p^+_B))d\textbf{x}.
  \label{reciprocity}
  \end{equation}
  Here, we consider a volume $\mathcal{V}$, which is limited by two infinite horizontal surfaces. States A and B (Fig. \ref{States}) are defined in two different media which are identical inside the volume $\mathcal{V}$ with boundaries denoted by $\mathfrak{d}\mathcal{V}_u$ (upper) and $\mathfrak{d}\mathcal{V}_l$ (lower), and arbitrary outside of this volume. In addition, $p^+$ and $p^-$ are decomposed wavefields on the boundaries where $+$ means downgoing and $-$ means upgoing, and $\partial_3$ is the partial derivative in direction $x_3$ (downward).

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja1a.pdf}
\caption{}
\label{stateA}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja1b.pdf}
\caption{}
\label{stateB}
\end{subfigure}
\caption{ (a) State A and (b) state B. Black stars depict source locations, and black reversed triangles depict receiver location.}
\label{States}
\end{figure}

  To continue the derivation we define Green's function $G(\textbf{x},\textbf{x}_s)$ as the solution of the following Helmholtz eqution:
  \begin{equation}
    \nabla^2 G + k^2 G = -\rho \delta(\textbf{x}-\textbf{x}_s),
     \label{waveeq}
  \end{equation}
  where $k(\textbf{x})=\frac{\omega}{c(\textbf{x})}$ is the wavenumber, and $c$ is the propagation velocity.

    A representation for Green's function at the upper boundary of the target ($\mathfrak{d}\mathcal{V}_u$) can be derived by defining state A in the actual medium and state B in a medium identical to medium A inside the volume $\mathcal{V}$ (target) and homogeneous outside of it. We denote the Green's functions of state A with $G(\textbf{x},\textbf{x}_s)$ and Green's function of state B with $G_{tar}(\textbf{x},\textbf{x}_u)$. In state A we define $p_A^\pm=G^\pm(\textbf{x},\textbf{x}_s)$, where $\textbf{x}_s$ is a location at Earth's surface $\mathfrak{d}\mathcal{V}_0$. For state B we have $p_B^\pm=G_{tar}^\pm(\textbf{x},\textbf{x}_u)$. Here, $\textbf{x}_u$ is located at the upper boundary $\mathfrak{d}\mathcal{V}_u$ and $G_{tar}^\pm(\textbf{x},\textbf{x}_u)$ is stimulated by an impulsive point source at $\textbf{x}_u$. Since the half-space above the target is homogeneous in state B, the vertical derivative of the $G_{tar}^+(\textbf{x},\textbf{x}_u)$ at $\mathfrak{d}\mathcal{V}_u$ is
  \begin{equation}
    \lim_{x_{u,3}\to x_{u,3^+}}\partial_3 G_{tar}^+(\textbf{x},\textbf{x}_u)=-\frac{\rho}{2}\delta(\textbf{x}_{\rm H}-\textbf{x}_{u,{\rm H}}).
  \label{dG+}
  \end{equation}
  Here, ${x_{u,3}\to x_{u,3^+}}$ means the limit from below the boundary, where $\textbf{x}_H$ and $\textbf{x}_{u,H}$ denote the horizontal coordinates of $\textbf{x}$ and $\textbf{x}_u$, respectively. In addition, since in state B the medium below the target is reflection-free, $G^-_{tar}(\textbf{x},\textbf{x}_u)$ and its derivative disappear at the lower boundary $\mathfrak{d}\mathcal{V}_l$.
  By substituting all of the ingredients into equation \ref{reciprocity}, the following can be reached:
  \begin{equation}
  G^-(\textbf{x}_u,\textbf{x}_s)=\int_{\mathfrak{d}\mathcal{V}_u}G^+(\textbf{x},\textbf{x}_s)\frac{2\partial_3}{\rho(\textbf{x})}G^-_{tar}(\textbf{x},\textbf{x}_u)d\textbf{x} \\
  +\int_{\mathfrak{d}\mathcal{V}_l}G^-(\textbf{x},\textbf{x}_s)\frac{2\partial_3}{-\rho(\textbf{x})}G^+_{tar}(\textbf{x},\textbf{x}_u)d\textbf{x}.
  \label{representation}
  \end{equation}
This equation is the base for our target-enclosed LSRTM derivation. The first integral on the right-hand side of Equation \ref{representation} accounts for anything entering the medium from the upper boundary, and the second integral accounts for anything that comes from the lower boundary.

  \subsection{Target-enclosed LSRTM}
To merge this representation with LSRTM, we use
  \begin{equation}
  \begin{split}
  \int_{\mathfrak{d}\mathcal{V}}p^\pm_A(\textbf{x})\partial_3p^\mp_B(\textbf{x}) d\textbf{x} = -\int_{\mathfrak{d}\mathcal{V}}p^\mp_B(\textbf{x})\partial_3p^\pm_A(\textbf{x}) d\textbf{x}
  \end{split}
  \label{transpose}
  \end{equation}
 \cite{Keesbook}, and we use Green's function reciprocity:
  \begin{equation}
    G^-_{tar}(\textbf{x}'_u,\textbf{x}_u) = G^-_{tar}(\textbf{x}_u,\textbf{x}'_u)
  \label{Gupreciprocity}
  \end{equation}
  \begin{equation}
    G^+_{tar}(\textbf{x}_l,\textbf{x}_u) = G^-_{tar}(\textbf{x}_u,\textbf{x}_l),
  \label{Glowreciprocity}
  \end{equation}
 where $\textbf{x}'_u$ is an element of $\mathfrak{d}\mathcal{V}_u$ and $\textbf{x}_l$ is an element of $\mathfrak{d}\mathcal{V}_l$. We rewrite Equation \ref{representation} as follows:
  \begin{equation}
    G^-(\textbf{x}_u,\textbf{x}_s)W(\omega)=\int_{\mathfrak{d}\mathcal{V}_u}G^-_{tar}(\textbf{x}_u,\textbf{x}'_u)S_u(\textbf{x}'_u,\textbf{x}_s)d\textbf{x}'_{u} \\
    +\int_{\mathfrak{d}\mathcal{V}_l}G^-_{tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l},
  \label{representation_alt}
  \end{equation}
where $W(\omega)$ is the source signature,
  \begin{equation}
  S_u(\textbf{x}'_u,\textbf{x}_s)=\frac{2\partial'_{3,u}}{-\rho(\textbf{x}'_u)}G^+(\textbf{x}'_u,\textbf{x}_s)W(\omega)
  \label{Source_udf}
  \end{equation}
is the dipole source from the upper boundary of the target, which accounts for reflections from above the upper boundary of the target. Further
  \begin{equation}
    S_l(\textbf{x}_l,\textbf{x}_s)=\frac{2\partial_{3,l}}{\rho(\textbf{x}_l)}G^-(\textbf{x}_l,\textbf{x}_s)W(\omega)
   \label{Source_ldf}
   \end{equation}
is the dipole source term from the lower boundary of the target, which accounts for the reflections generated below the target. Fig. \ref{integrals} represents the right-hand side of Equation \ref{representation_alt}.

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja2a.pdf}
\caption{}
\label{integralup}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja2b.pdf}
\caption{}
\label{integrallow}
\end{subfigure}
\caption{(a) The first integral on the right-hand side of Equation \ref{representation_alt}, and (b) the second integral on the right-hand side of Equation \ref{representation_alt}.}
\label{integrals}
\end{figure}

As said above, the second integral on the right-hand side of equation \ref{representation_alt} is the contribution of the medium below the target to the data. This integral can be split into two terms: 1) the arrival from the lower boundary to the upper one in the background, and 2) the forward scatterings inside the target region. In mathematical terms:
  \begin{equation}
  \int_{\mathfrak{d}\mathcal{V}_l}G^-_{tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l}= \\
  \int_{\mathfrak{d}\mathcal{V}_l}G^-_{0,tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l} \\
  +\int_{\mathfrak{d}\mathcal{V}_l}G^-_{sct,tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l}.
    \label{lower}
  \end{equation}
Here, $G^-_{0,tar}(\textbf{x}_u,\textbf{x}_l)$ is the Green's function in the background model of the target, and $G^-_{sct,tar}(\textbf{x}_u,\textbf{x}_l)$ is the Green's function that contains the scattered events. By substituting equation \ref{lower} in equation \ref{representation_alt} and taking the background contribution to the left-hand side we end up with:
  \begin{equation}
  G^-(\textbf{x}_u,\textbf{x}_s)W(\omega)-\int_{\mathfrak{d}\mathcal{V}_l}G^-_{0,tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l}= \\
  \int_{\mathfrak{d}\mathcal{V}_u}G^-_{sct,tar}(\textbf{x}_u,\textbf{x}'_u)S_u(\textbf{x}'_u,\textbf{x}_s)d\textbf{x}'_{u} \\
  +\int_{\mathfrak{d}\mathcal{V}_l}G^-_{sct,tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l}.
  \label{modifired_representation_df}
  \end{equation}
Importantly, the Green's function $G^-_{tar}(\textbf{x}_u,\textbf{x}'_u)$ in the first integral on the right-hand side of equation \ref{representation_alt} is the scattered Green's function inside the target, so we rename it to $G^-_{sct,tar}(\textbf{x}_u,\textbf{x}'_u)$.

To obtain our target-enclosed LSRTM algorithm, we take three steps: First, we assign
  \begin{equation}
  P^{obs}_{sct,TE}(\textbf{x}_u,\textbf{x}_s)= \\
  G^-(\textbf{x}_u,\textbf{x}_s)W(\omega)-\int_{\mathfrak{d}\mathcal{V}_l}G^-_{0,tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l},
    \label{PobsTE}
  \end{equation}
where "TE" stands for "Target-Enclosed". Second we compute the scattered Green's functions ($G^-_{sct,tar}(\textbf{x}_u,\textbf{x}'_u)$ and $G^-_{sct,tar}(\textbf{x}_u,\textbf{x}_l)$) with equation \ref{scattering}
  \begin{equation}
    G^-_{sct,tar}(\textbf{x}_u,\textbf{x}'_u) = \\
  \int_{\mathcal{V}}\gamma_0^2(\textbf{x})G_{0,tar}(\textbf{x}_u,\textbf{x})\chi^c(\textbf{x})G_{0,tar}(\textbf{x},\textbf{x}'_u)d\textbf{x},
  \end{equation}
and
  \begin{equation}
    G^-_{sct,tar}(\textbf{x}_u,\textbf{x}_l) = \\
    \int_{\mathcal{V}}\gamma_0^2(\textbf{x})G_{0,tar}(\textbf{x}_u,\textbf{x})\chi^c(\textbf{x})G_{0,tar}(\textbf{x},\textbf{x}_l)d\textbf{x},
  \end{equation}
  where $\textbf{x}$ is a location inside the target volume. Using this in the right-hand side of Equation \ref{modifired_representation_df}, we obtain
  \begin{equation}
  \int_{\mathfrak{d}\mathcal{V}_u}G^-_{sct,tar}(\textbf{x}_u,\textbf{x}'_u)S_u(\textbf{x}'_u,\textbf{x}_s)d\textbf{x}'_{u} \\
  +\int_{\mathfrak{d}\mathcal{V}_l}G^-_{sct,tar}(\textbf{x}_u,\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l}= \\
  \int_{\mathcal{V}}\gamma_0^2(\textbf{x})G_{0,tar}(\textbf{x}_u,\textbf{x})\chi^c(\textbf{x})P_{inc,TE}(\textbf{x},\textbf{x}_s)d\textbf{x},
  \end{equation}
 where
  \begin{equation}
  P_{inc,TE}(\textbf{x},\textbf{x}_s)= \\
  \int_{\mathfrak{d}\mathcal{V}_u}G_{0,tar}(\textbf{x},\textbf{x}'_u)S_u(\textbf{x}'_u,\textbf{x}_s)d\textbf{x}'_{u} \\
  +\int_{\mathfrak{d}\mathcal{V}_l}G_{0,tar}(\textbf{x},\textbf{x}_l)S_l(\textbf{x}_l,\textbf{x}_s)d\textbf{x}_{l}.
  \end{equation}
Finally, we complete our derivation by assigning
  \begin{equation}
  P^{pred}_{sct,TE}(\textbf{x}_u,\textbf{x}_s)= \\
  \int_{\mathcal{V}}\gamma_0^2(\textbf{x})G_{0,tar}(\textbf{x}_u,\textbf{x})\chi^c(\textbf{x})P_{inc,TE}(\textbf{x},\textbf{x}_s)d\textbf{x}.
  \label{PpredTE}
  \end{equation}

Our approach computes an incident wavefield, which contains all of the information from the surrounding medium of the target of interest. This means our approach is not limited to one kind of parametrization, and it can be implemented for any parametrization choice. Finally, we can solve the following least-squares problem:
  \begin{equation}
    J(\delta \textbf{m})=\frac{1}{2}\|\textbf{P}^{pred}_{sct,TE}(\delta \textbf{m})-\textbf{P}^{obs}_{sct,TE}\|^2_2.
    \label{costfun_modified}
  \end{equation}

\subsection{Marchenko Green's function retrieval}

  We obtained a target-enclosed LSRTM algorithm in the previous section. Nevertheless, in most real-world situations, one does not have physical access to the boundaries of the region of interest. The alternative to physical receivers inside the medium is to create virtual receivers with redatuming. Marchenko redatuming is a state-of-the-art data-driven approach that can compute Green's functions at any depth level with all orders of multiple reflections from the reflection response at the surface and a smooth background model of the medium.

  To summarise, these redatumed Green's functions are retrieved by iteratively solving the Marchenko-type representations. These representations are \cite{Wapenaar6}

  \begin{equation}
    G^-_{Mar}(\textbf{x}_{v},\textbf{x}_{s}) = \int_{{\mathfrak{d}\mathcal{V}_0}} R(\textbf{x}_s,\textbf{x}'_s)f_1^+(\textbf{x}'_s,\textbf{x}_{v}) \,d\textbf{x}'_s  \\
     -f_1^-(\textbf{x}_{s},\textbf{x}_{v}),
    \label{Mar_r-}
\end{equation}

and

\begin{equation}
    G^+_{Mar}(\textbf{x}_{v},\textbf{x}_{s}) = -\int_{{\mathfrak{d}\mathcal{V}_0}} R(\textbf{x}_s,\textbf{x}'_s)f_1^-(\textbf{x}'_s,\textbf{x}_{v})^* \,d\textbf{x}'_s \\
     + f_1^+(\textbf{x}_{s},\textbf{x}_{v})^* .
    \label{Mar_r+}
\end{equation}

Here, $\mathfrak{d}\mathcal{V}_0$ is the surface, $\textbf{x}_s$ and $\textbf{x}'_s$ are locations at the surface, and $\textbf{x}_v$ is a virtual location on an arbitrary depth. Moreover, $f^\pm_1$ are upgoing (-) and downgoing (+) parts of the focusing function. In addition, $R(\textbf{x}_s,\textbf{x}'_s)$ is the reflection response at the surface which is related to the upgoing Green's function of the medium via

\begin{equation}
    R(\textbf{x}_s,\textbf{x}'_s) = \frac{2\partial'_{3,s}}{-\rho(\textbf{x}'_s)} G^-(\textbf{x}_s,\textbf{x}'_s).
\label{R_to_r}
\end{equation}

Here, $\partial'_{3,s}$ is the vertical partial derivative at $\textbf{x}'_s$. We refer to \cite{Thorbecke2} for a comprehensive explanation of the derivation and numerical algorithms for solving these equations.

Consequently, we can substitute Marchenko Green's functions with the target boundaries Green's functions as follows:

\begin{equation}
    G^+(\textbf{x}_{u},\textbf{x}_{s}) \approx G^+_{Mar}(\textbf{x}_{u},\textbf{x}_{s}),
\label{Mar_Gp}
\end{equation}

\begin{equation}
    G^-(\textbf{x}_{u},\textbf{x}_{s}) \approx G^-_{Mar}(\textbf{x}_{u},\textbf{x}_{s}),
\label{Mar_Gm_up}
\end{equation}

and

\begin{equation}
    G^-(\textbf{x}_{l},\textbf{x}_{s}) \approx G^-_{Mar}(\textbf{x}_{l},\textbf{x}_{s}).
\label{Mar_Gm_low}
\end{equation}

\section{\label{sec:3} Numerical results}

\subsection{\label{sec:3.1} single-sided algorithm vs. double-sided algorithm}

In this section, we aim to visualize the performance of target-enclosed LSRTM. A model with dimensions of 1000 $m$ by 650 $m$ is designed as shown in Fig. \ref{model}. The spatial grid sampling is 5 $m$ in both directions. The target region consists of a rectangular velocity anomaly embedded in a homogeneous background and a constant density. A total of 201 sources are placed at the surface, and 402 receivers are on the target's upper (250 $m$) and lower (550 $m$) boundaries. The required Green's functions and wavefields are computed by a finite-difference algorithm \cite{Thorbecke1} and a Ricker wavelet with a dominant frequency of 30 $Hz$, where the recording time sampling of the receivers is set to 4 $ms$. According to the theory section, the wavefields at the receiver positions are decomposed into upgoing and downgoing components. Fig. \ref{obs} shows the upgoing component of the data at the upper boundary with a source located at $\textbf{x}_s=(0,0)$, corresponding to the first term on the right-hand side of Equation \ref{PobsTE}. We study the effects of including the lower boundary with two different cases: 1) a homogeneous background velocity and 2) a smooth background velocity for migration.

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\includegraphics[width=\textwidth]{Figures/Shoja3a.pdf}
\caption{}
\label{velocity_model}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\includegraphics[width=\textwidth]{Figures/Shoja3b.pdf}
\caption{}
\label{density_model}
\end{subfigure}
\caption{(a) Velocity model and (b) density model. The blue stars at the surface are source locations, and the red dots are the boundaries of the target.}
\label{model}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.4\textwidth]{Figures/Shoja4.pdf}
\caption{Upgoing component of data at the upper boundary with a source located at $\textbf{x}_s=(0,0)$. This is the data corresponding to the first term on the right-hand side of Equation \ref{PobsTE}.}
\label{obs}
\end{figure}

\subsubsection{\label{sec:3.1.1} Homogeneous background velocity}

In this section, we show the results using a homogeneous background velocity in the target area. The true perturbation model of this case is shown in Fig. \ref{pert_homogen}. We design two scenarios to demonstrate the performance and consequences of the target-enclosed algorithm. For the first scenario, we only consider the upper boundary of the target in the inversion process and completely ignore the lower boundary contribution. We call this scenario "single-sided algorithm". For this scenario, the observed data is the same as Fig. \ref{obs}. For the second one, we include the lower boundary contribution, which is the "double-sided algorithm" explained in the previous section. The observed data for this scenario is shown in Fig. \ref{obs_enc_homogen}, which corresponds to the left-hand side of Equation \ref{PobsTE}, whereas Fig. \ref{obs} shows $G^-(\textbf{x}_u,\textbf{x}_s)W(\omega)$. The black arrows in Figures \ref{obs} and \ref{obs_enc_homogen} indicate the total reflection from the reflector below the lower boundary and the forward scattered part of it, respectively. Figures \ref{pred_or_homogen} and \ref{pred_enc_homogen} show the predicted data of single-sided and double-sided algorithms after 30 iterations of LSRTM, respectively. A detailed investigation of Fig. \ref{data_homogen} proves that the double-sided algorithm can predict the forward scattered event that is passing through the perturbation which is indicated by a black arrow in Fig. \ref{pred_enc_homogen}.

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja5.pdf}
\caption{True perturbation model in a homogeneous background.}
\label{pert_homogen}
\end{figure}

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja6a.pdf}
\caption{}
\label{obs_enc_homogen}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja6b.pdf}
\caption{}
\label{pred_or_homogen}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja6c.pdf}
\caption{}
\label{pred_enc_homogen}
\end{subfigure}
\caption{Homogeneous background velocity case, data domain: (a) Upgoing component of double-sided observed data at the upper boundary (Equation \ref{PobsTE}), (b) predicted data of single-sided algorithm after 30 iterations of LSRTM, and (c) predicted data of double-sided algorithm after 30 iterations of LSRTM. All wavefields are recorded at the upper boundary of the target (250 $m$).}
\label{data_homogen}
\end{figure}

To move our investigation further, we show the imaging results of both algorithms in Fig. \ref{img_homogen}. Figures \ref{RTM_or_homogen} and \ref{RTM_Enc_homogen} show the RTM images of both approaches. As we can see, the RTM result of the double-sided algorithm faintly reveals the long wavelength part of the model. Moving to the LSRTM results in Figures \ref{LSRTM_or_homogen} and \ref{LSRTM_Enc_homogen}, we observe an interesting outcome. The LSRTM result of the single-sided algorithm shows that it cannot recover the long wavelength part of the model (Fig. \ref{LSRTM_or_homogen}). In contrast, the double-sided algorithm can incorporate the information embedded inside the forward scattered field (Fig. \ref{LSRTM_Enc_homogen}), and it recovers the long wavelength parts of the volume perturbation. Moreover, Fig. \ref{horiz_homogen} shows the horizontal cross-section of the retrieved perturbation. In this figure, we can see the double-sided approach can recover the vertical boundaries of the perturbation. Nevertheless, since the background velocity for the migration is not updated during LSRTM, the fit of the reflected event from below the target is not accurate.

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja7a.pdf}
\caption{}
\label{RTM_or_homogen}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja7b.pdf}
\caption{}
\label{RTM_Enc_homogen}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja7c.pdf}
\caption{}
\label{LSRTM_or_homogen}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja7d.pdf}
\caption{}
\label{LSRTM_Enc_homogen}
\end{subfigure}
\caption{Homogeneous background velocity case, image domain: (a) RTM image of the single-sided algorithm, (b) RTM image of the double-sided algorithm, (c) LSRTM image of the single-sided algorithm after 30 iterations, and (d) LSRTM image of the double-sided algorithm after 30 iterations.}
\label{img_homogen}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja8.pdf}
\caption{Horizontal cross-section at the depth of 400 m of the retrieved perturbation with a homogeneous migration velocity.}
\label{horiz_homogen}
\end{figure}

To conclude this section, we compare the cost functions of both approaches in Fig. \ref{cost_homogen}. The cost function of target-oriented LSRTM shows a slow convergence rate. In comparison, the target-enclosed approach includes the extra information coming from the lower boundary, so its cost function converges faster and to a lower minimum.

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja9.pdf}
\caption{Cost function comparison of homogeneous background velocity case.}
\label{cost_homogen}
\end{figure}

\subsubsection{\label{sec:3.1.2} Smooth background velocity}

This section uses a smooth background velocity for migration. The setup is exactly the same as before, except for the background velocity. The perturbation model for this case is shown in Fig. \ref{pert_smooth}. Again, we do the same two scenarios as before, i.e. "single-sided" and "double-sided." For the single-sided scenario, the observed data is the same as before (Fig. \ref{obs}). However, since the right-hand side of the equation \ref{PobsTE} is computed in a different background velocity, the last primary in the observed data for the double-sided approach (Fig. \ref{obs_enc_smooth}) is slightly different.

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja10.pdf}
\caption{True perturbation model in a smooth background.}
\label{pert_smooth}
\end{figure}

Comparing the results of both approaches shows (Fig. \ref{data_smooth}) that the double-sided approach (Fig. \ref{pred_enc_smooth}) can predict the reflected event coming from the lower boundary. Since we use a smooth background velocity here, this prediction is more accurate than the previous section's results (Fig. \ref{pred_enc_homogen}). In the image domain (Fig. \ref{img_smooth}), the double-sided approach recovers a faint image of the vertical sides of the rectangle (Fig. \ref{LSRTM_Enc_smooth}) whereas the single-sided results in a more standard image (Fig. \ref{LSRTM_or_smooth}). Moreover, Fig. \ref{horiz_smooth} shows the horizontal cross-section of the retrieved perturbation. In this figure, we can see the double-sided approach can recover the vertical boundaries of the perturbation. The double-sided image is more comparable to the true perturbation in Fig. \ref{pert_smooth}. Finally, investigating the cost functions (Fig. \ref{cost_smooth}) of these approaches shows that the double-sided approach converges better since it can predict the event coming from below the target.

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja11a.pdf}
\caption{}
\label{obs_enc_smooth}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja11b.pdf}
\caption{}
\label{pred_or_smooth}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja11c.pdf}
\caption{}
\label{pred_enc_smooth}
\end{subfigure}
\caption{Smooth background velocity case, data domain: (a) Upgoing component of double-sided observed data at the upper boundary (Equation \ref{PobsTE}), (b) predicted data of single-sided algorithm after 30 iterations of LSRTM, and (c) predicted data of double-sided algorithm after 30 iterations of LSRTM. All wavefields are recorded at the upper boundary of the target (250 $m$).}
\label{data_smooth}
\end{figure}

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja12a.pdf}
\caption{}
\label{RTM_or_smooth}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja12b.pdf}
\caption{}
\label{RTM_Enc_smooth}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja12c.pdf}
\caption{}
\label{LSRTM_or_smooth}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Shoja12d.pdf}
\caption{}
\label{LSRTM_Enc_smooth}
\end{subfigure}
\caption{Smooth background velocity case, image domain: (a) RTM image of the single-sided algorithm, (b) RTM image of the double-sided algorithm, (c) LSRTM image of the single-sided algorithm after 30 iterations, and (d) LSRTM image of the double algorithm after 30 iterations.}
\label{img_smooth}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja13.pdf}
\caption{Horizontal cross-section at the depth of 400 m of the retrieved perturbation with a smooth migration velocity.}
\label{horiz_smooth}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{Figures/Shoja14.pdf}
\caption{Cost function comparison of smooth background velocity case.}
\label{cost_smooth}
\end{figure}

\subsection{\label{sec:3.2} Virtual receivers}

Here, we use the same setup as before but replace the Green's functions at the boundaries with their Marchenko counterparts. In other words, we create virtual receivers with the help of Marchenko redatuming. For a study about the benefits of using Marchenko redatuming instead of a more conventional redatuming algorithm for target-oriented LSRTM we refer to \cite{Shoja3}. In this section, we only focus on including the lower boundary by utilizing virtual receivers created by Marchenko redatuming, and instead of using "double-sided", we use the "target-enclosed" term. Moreover, we show the results for both homogeneous and smooth background velocities in this section.

\subsubsection{\label{sec:3.2.1} Homogeneous background velocity}
Fig. \ref{virobs} shows the data obtained by Marchenko redatuming. In Fig. \ref{virdata_homogen} the observed data calculated by Equation \ref{PobsTE} (Fig. \ref{Marchenko_obs_homogen}) and the predicted data after 30 iterations (Fig. \ref{Marchenko_pred_homogen}) are shown. A comparison between Fig. \ref{Marchenko_pred_homogen} and Fig. \ref{pred_enc_homogen} shows that even in the presence of redatuming error, such as limited aperture and lack of certain parts of the wavelength spectrum, our algorithm can predict acceptable data.
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{Figures/Shoja15.pdf}
    \caption{Marchenko redatumed data with virtual receivers at the upper boundary (250 m) and a source located at $\textbf{x}_s=(0,0)$.}
    \label{virobs}
\end{figure}

\begin{figure}
\centering
\begin{subfigure}{0.4\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja16a.pdf}
\caption{}
\label{Marchenko_obs_homogen}
\end{subfigure}
\hfill
\begin{subfigure}{0.4\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja16b.pdf}
\caption{}
\label{Marchenko_pred_homogen}
\end{subfigure}
\caption{Homogeneous background velocity case with virtual receivers, data domain: (a) Upgoing component of Marchenko redatumed data at the upper boundary (Equation \ref{PobsTE}), (b) predicted data of target-enclosed LSRTM after 30 iterations at the upper boundary. All wavefields are redatumed to the upper boundary of the target (250 $m$).}
\label{virdata_homogen}
\end{figure}

In the image domain, Fig. \ref{virimg_homogen} shows the RTM (Fig. \ref{Marchenko_RTM_homogen}) and LSRTM (Fig. \ref{Marchenko_LSRTM_homogen}) image resulting from Marchenko redatumed data. Comparing Fig. \ref{Marchenko_LSRTM_homogen} with Fig. \ref{LSRTM_Enc_homogen} reveals our target-enclosed algorithm with redatumed data as input can not recover the long wavelength part of the model. This is due to the fact that the direct arrival of the Marchenko-based Green's function is incorrect since it is computed in the background model. Consequently, the forward-scattered waveforms responsible for the long wavelength updates in Figure 7d can no longer be utilized.
\begin{figure}
\centering
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja17a.pdf}
\caption{}
\label{Marchenko_RTM_homogen}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja17b.pdf}
\caption{}
\label{Marchenko_LSRTM_homogen}
\end{subfigure}
\caption{Homogeneous background velocity case with virtual receivers, image domain: (a) RTM image of the target-enclosed algorithm with Marchenko wavefields, (b) LSRTM image of the target-enclosed algorithm with Marchenko wavefields after 30 iterations.}
\label{virimg_homogen}
\end{figure}

\subsubsection{\label{sec:3.2.2} Smooth background velocity}

In Fig. \ref{virdata_smooth} the observed data calculated by Equation \ref{PobsTE} (Fig. \ref{Marchenko_obs_smooth}) and the predicted data after 30 iterations (Fig. \ref{Marchenko_pred_smooth}) are shown. Similar to the homogeneous case, a comparison between Fig. \ref{virdata_smooth} and Fig. \ref{data_smooth} shows that our algorithm successfully predicts the data.
\begin{figure}
\centering
\begin{subfigure}{0.4\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja18a.pdf}
\caption{}
\label{Marchenko_obs_smooth}
\end{subfigure}
\hfill
\begin{subfigure}{0.4\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja18b.pdf}
\caption{}
\label{Marchenko_pred_smooth}
\end{subfigure}
\caption{Smooth background velocity case with virtual receivers, data domain: (a) Upgoing component of Marchenko redatumed data at the upper boundary (Equation \ref{PobsTE}), (b) predicted data of target-enclosed LSRTM after 30 iterations at the upper boundary. All wavefields are redatumed to the upper boundary of the target (250 $m$).}
\label{virdata_smooth}
\end{figure}

In the image domain, Fig. \ref{virimg_smooth} shows the RTM (Fig. \ref{Marchenko_RTM_smooth}) and LSRTM (Fig. \ref{Marchenko_LSRTM_smooth}) image resulting from Marchenko redatumed data. Comparing Fig. \ref{Marchenko_LSRTM_smooth} with Fig. \ref{LSRTM_Enc_smooth} shows that redatumed data reveals an acceptable perturbation model. However, the faint recovered vertical interfaces are not presented in Fig. \ref{Marchenko_LSRTM_smooth} since forward-scattered waveforms are not processed correctly in retrieved Marchenko Green's function at the lower boundary.

\begin{figure}
\centering
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja19a.pdf}
\caption{}
\label{Marchenko_RTM_smooth}
\end{subfigure}
\hfill
\begin{subfigure}{0.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{Figures/Shoja19b.pdf}
\caption{}
\label{Marchenko_LSRTM_smooth}
\end{subfigure}
\caption{Smooth background velocity case with virtual receivers, image domain: (a) RTM image of the target-enclosed algorithm with Marchenko wavefields, (b) LSRTM image of the target-enclosed algorithm with Marchenko wavefields after 30 iterations.}
\label{virimg_smooth}
\end{figure}

\section{Discussion}
In Section \ref{sec:2}, we develop a theory for target-enclosed LSRTM that can limit the computation domain by confining the target between two boundaries. Equations \ref{PobsTE} and \ref{PpredTE} enable us to account for any wavefield entering the target region by including the upper and lower boundaries in the inversion process.

Further, in section \ref{sec:3.1}, several numerical tests are designed to demonstrate the advantage of incorporating the lower boundary in the conventional target-oriented LSRTM. From the data point of view, our double-sided target-enclosed LSRTM, compared to a conventional single-sided target-oriented LSRTM, not only removes the background arrival from the lower boundary to the upper boundary but also can predict the forward scattered field inside the target (Fig. \ref{data_homogen}, Fig. \ref{data_smooth}, and Fig. \ref{virdata_homogen}). Additionally, a comparison between the resulting images of both algorithms with a homogeneous migration velocity (Fig. \ref{img_homogen}), shows that single-sided target-oriented LSRTM can not update the long wavelength part of the model. In contrast, the double-sided target-enclosed LSRTM updates the image according to the Born inversion criteria by integrating the forward scattered wavefield information. In the case of a smooth migration velocity, our algorithm recovers a higher-resolution image and a faint estimation of the vertical sides of the rectangle (Fig. \ref{img_smooth}).

Moreover, in section \ref{sec:3.2}, we investigate the possibility of using virtual receivers created by Marchenko redatuming in our target-enclosed algorithm. Fig. \ref{virimg_homogen} shows that in the case of a homogeneous background migration velocity, it is hardly possible to update the long wavelengths with virtual receiver data, and only short wavelength parts of the perturbation are recovered. Moreover, for the smooth background case, our algorithm increases the image's resolution by updating the short wavelengths in this case.

\section{Conclusion}

This paper proposes a target-enclosed algorithm that can account for the wavefields entering the target region from the upper and lower boundaries of the region. The three main advantages of this paper's algorithm are 1) It significantly reduces the computational domain by limiting it to a smaller domain, 2) it removes interactions with the part of the medium above the upper boundary, and 3) it can incorporate the transmission information from the lower boundary to the upper one.

Nevertheless, our algorithm has some limitations. First, we need physical access to the lower boundary of the target to deploy receivers at the boundaries of the target. Second, we need a background model of the target that can predict the arrival time from the lower boundary to the upper. It is possible to overcome the first limitation by using virtual seismology methods such as Marchenko redatuming to create virtual receivers around the target region \cite{Broggini,Barrera,Diekman2,Dukalski,Ravasi1,Staring,Wapenaar6,Shoja3} as we showed in a numerical example in section \ref{sec:3.2}. For the second one, a reformulation of the target-enclosed LSRTM is possible to make it compatible with full waveform inversion to update the background velocity model \cite{Shoja}.

The need for high-resolution images is increasing daily, which demands more computational power. This paper's proposed target-enclosed LSRTM can produce less computationally demanding high-resolution images by focusing on a relatively small target of interest and including all the interactions between this region and the outside environment.

\section*{conflict of interest}
All authors declare that they have no conflicts of interest.

\bibliographystyle{IEEEtran}


\subsection{\@startsection                        %
  {subsection}{2}{0mm
  {1.5\baselineskip
  {1\baselineskip
  {\indent\normalfont\normalsize\itshape}
\renewcommand*{\@seccntformat}[1]
  \normalfont\bfseries\csname the#1\endcsname\;-               %
}                                                              %
\renewcommand\subsubsection{\@startsection                     %
  {subsubsection}{2}{0mm
  {1.5\baselineskip
  {1\baselineskip
  {\indent\normalfont\normalsize\texttt}
\makeatother                                                   %
                                                               %

\begin{document}
\thispagestyle{empty}

$ $ \vspace {-2.2cm}
\begin{center}
\rule{8.8cm}{0.5pt}
\\[-5pt]
\rule{8.8cm}{0.5pt}
\end{center}
\vspace {1.7cm}

\begin{center}
{\sc\large Giuseppe Toscani},  \
{\sc\large Mattia Zanella}   \
\end{center}
\vspace {1.1cm}

\centerline{\large{\textbf{On a kinetic description of Lotka-Volterra dynamics}}}

\renewcommand{\thefootnote}{\fnsymbol{footnote}}

\renewcommand{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}

\vspace{0.6cm}
\begin{center}
\begin{minipage}[t]{11cm}
\small{
\noindent \textbf{Abstract.}
Owing to the analogies between the problem of wealth redistribution with taxation in a multi-agent society, we introduce and discuss a kinetic model describing  the statistical distributions in time of the sizes of groups of biological systems with prey-predator dynamic. While the evolution of the mean values is shown to be driven by a classical Lotka-Volterra system of differential equations, it is shown that the time evolution of the probability distributions of the size of groups of the two interacting species is heavily dependent both on a kinetic redistribution operator and the degree of randomness present in the system. Numerical experiments are given to clarify the time-behavior of the distributions of groups of the species.
\medskip

\noindent \textbf{Keywords.}
Kinetic models, Boltzmann equation, Wealth distribution, Taxation and redistribution, Lotka-Volterra equations.
\medskip

\noindent \textbf{Mathematics~Subject~Classification:}
91B60, 82C40, 35B40.

}
\end{minipage}
\end{center}

\bigskip

\section{Introduction}

In one of the scientific collaborations between one of the present authors with Professor Giampiero Spiga in the area of kinetic theory of multi-agent systems, the problem of taxation and consequent redistribution of taxed resources in a western society had been studied  \cite{BST}. The purpose of the work was the possibility of understanding how, in the case of a fixed mean wealth, redistribution operators can affect the tails of the equilibrium density of a kinetic model for wealth exchanges.

In this area, there have been several attempts to model simple economies in which a wealth exchange process produces in time a steady profile of wealth distribution similar to that observed in real economies. While the most popular approach  relies on methods borrowed from statistical mechanics for interacting particle systems \cite{Cha,ChCh,CCM,CCS,CPT,DY,DMT,DMT2,Hay,IKR,MT,TTZ,Sla} (cf. also the books \cite{CYC,PT}), other modelling settings that might be considered to constitute a mesoscopic approach is based on generalized Lotka-Volterra models \cite{MBR,SR}. These dynamical systems, also known as the predator-prey equations,  are frequently used to describe the evolution of biological systems in which two species interact, one as a predator and the other as prey \cite{Murray}. In particular, denoting by $x(t)\ge 0$ and $y(t)\ge0$  the size of the populations at time $t\ge 0$, their time evolution is supposed to be influenced by predation, growth and death processes. Therefore, this dynamics is classically described by the following pair of ordinary differential equations
\begin{equation} \label{LV}
\begin{cases}
\frac{d}{dt}x(t) = \alpha x(t) -\beta \,x(t)y(t) ,\\
 \frac{d}{dt} y(t)= - \delta\, y(t) + \gamma \,x(t)y(t) ,
 \end{cases}
\end{equation}
where the preys' population is characterized by the rate of Malthusian growth $\alpha\ge 0$ and by a predation rate $\beta\ge 0$. Furthermore, the predators' population is subject to losses due to death  with rate  $\delta\ge 0$ and on the growth due to predation defined by the parameter $\gamma\ge0$.

While less popular, the approach proposed in \cite{MBR,SR} is capable to establish an interesting parallel between prey-predator dynamics and economic modelling. This analogy, as documented by \cite{Gan}, goes back to the Italian economist Giuseppe Palomba, who first had used  Lotka-Volterra equations \fer{LV} in the 1939 book \cite{Pal}, and to Goodwin, who considered applications of Lotka-Volterra system to economy in the 1967 book \cite{Good}.

In this paper, we provide a kinetic interpretation of the introduced dynamics by means of a system of Boltzmann-type equations. The interacting species are characterized from a statistical perspective according to their size. Hence, the temporal evolution of the statistical distributions of two species, classically denoted as prey and predators,  is subject to variation due to birth and death as well as to interactions between the two families (predation). The building blocks of the model are given by the redistribution operator introduced  in \cite{BST}, which allows for a general representation of the birth and death process of the involved species, and a kinetic bilinear operator of Boltzmann type which quantifies the interaction between species.
Interestingly, the structure of the interaction is similar to that recently introduced by the authors in connection with the study of kinetic models suitable for describing several patterns in compartmental epidemiology connected with the transition between susceptible and infected individuals \cite{DPaTZ, DPTZ}.

The linear operator introduced in \cite{BST} is sufficiently flexible to redistribute to agents a constant amount of wealth independently of the wealth itself, or to redistribute proportionally (or inversely proportionally) to their wealth. Hence, in the new application to a biological system, it allows the number of births  to be quantified in various ways and in terms of the size of the groups of the species.

In more details, the paper is organized as follows: in Section \ref{sect:2} we introduce a kinetic model describing the evolution of the sizes of groups of biological systems obeying a prey-predator dynamic. The transitions depend on elementary interactions determining the evolution in time of probability densities described by a system of non-Maxwellian Boltzmann-like equations. A redistribution operator depending only on the first two moments of the distribution is then considered. In Section \ref{sec:3} we recall basic results on the grazing limit in kinetic theory and we provide a Fokker-Planck description of the dynamics.  Finally, in Section \ref{sec:4} we provide several numerical tests to verify the consistency of the approach with the classical Lotka-Volterra model and we observe the emergence of fat-tailed distributions, in the group of preys, for suitable choices of the redistribution operator.

\section{The model}\label{sect:2}

Following the general approach of \cite{PT} we consider a multiagent system composed by two species, identified as preys and predators,  that are characterized by their number of individuals $x\ge 0$ (respectively $y\ge0$). We denote by $f(x,t)$ (respectively $g(y,t)$) be the statistical distributions at a certain time $t \ge 0$ of the group of preys (respectively predators), classified in terms of their sizes $x,y \in \mathbb R_+$.
Therefore, $f(x,t)dx$ indicates the fraction of the population preys with size in $[x,x + dx]$ at time $t\ge0$. Similarly, $g(y,t)dy$ is the fraction of the populations of predators characterized by a size in $[y,y+dy]$ at time $t\ge0$. Moreover, without loss of generality we fix the initial masses of the distributions $f$ and $g$ equal to one. Therefore, the mean sizes of the populations are given by
\[
m_{f}(t) = \int_{\mathbb R_+}xf(x,t)dx,\qquad m_{g}(t) = \int_{\mathbb R_+}yg(y,t)dy.
\]
The goal is to characterize the time evolution of the pair of probability densities $f(x,t)$, $g(y,t)$ subject to changes including  births and deaths events, and obeying to objective mechanical and behavioral principles when interacting each other.

Let us first characterize the interaction rules between the two species. To achieve this, let us fix a certain time $t= t_0$, considered as the initial time of the observations,  and let $x_0,y_0 $ be the numbers of preys and predators in two groups at time $ t_0 $. Then, we evaluate the statistical distribution of the size of the groups at regular time intervals $ t_{n+1} $, with $n$ nonnegative integer, $ n\in \mathbb{N}_+$, consequent to  the elementary variations $x_n \to x_{n+1}$ and  $y_n \to y_{n+1}$ of the groups in each time interval. Since a prey-predator  interaction leads to a decrease in the number of preys, and, in the presence of a sufficiently high number of preys, to an increase of the number of predators, this change can be modeled as follows
\begin{equations}\label{coll1}
 x_{n+1} &= x_n  -  \Phi_\epsilon(y_n) x_n +   \eta_\epsilon(y_n) x_n. \\
 y_{n+1} &= y_n  + \Psi_\epsilon(x_n) y_n +  \tilde\eta_\epsilon(x_n)  y_n.
\end{equations}
Thus, at each time step, the number  of preys and predators is modified by two different mechanisms, expressed in mathematical terms by two multiplicative terms, both parameterized by a small positive parameter $\epsilon \ll 1$, measuring the intensities of the single variations.

To clarify, the  function $\Phi_\epsilon(\cdot)$ in the first equation of \fer{coll1} quantifies the predation rate, and is consistently assumed to be increasing with respect to the size of the predator population. Likewise, the function $\Psi_\epsilon(\cdot)$ in the second equation of \fer{coll1} quantifies the growth of the predators, increasing with respect to the prey population at least when the number of preys is sufficiently consistent.  In reason of the previous remarks, in what follows we  fix
\begin{equation}\label{Phi}
\Phi_\epsilon(y) =\epsilon\beta \frac y{1+y}; \quad
\Psi_\epsilon(x) = \epsilon \gamma \frac {x-\mu}{1+x}.
\end{equation}
In \fer{Phi}, $\beta$, $\gamma$ and $\mu$ are positive constants. The constants $\beta$ and $\gamma$ measure the maximal predation rate, and, respectively, the maximal growth rate of predators. Last, the constant $\mu \ge 1$ measures the lowest size of the group of preys which is compatible with a growth of the group of predators. These constants are required to satisfy the bounds
\begin{equation}\label{bb}
\epsilon\beta < 1; \quad \epsilon\gamma\mu < 1.
\end{equation}
These conditions ensure that the predation rate is always smaller than one, and that even in absence of preys, the number of predators does not collapse to zero in a single interaction.

In \fer{coll1}, random fluctuations of the groups due to unknown factors  in the environment, like immigration, natural events,  epidemics, etc.  are  expressed by  the independent random variables $\eta_\epsilon(x), \tilde\eta_\epsilon(x)$,  of zero mean and bounded variance, given by
\begin{equation}\label{var1}
\langle \eta_\epsilon^2(y) \rangle  = \epsilon\sigma\frac y{1+y}; \quad  \langle \tilde\eta_\epsilon^2 (x)\rangle =    \epsilon\tilde\sigma \frac {x}{1+x}.
\end{equation}
Note that, consistently with the deterministic variations, the variance of the random fluctuations of one class  is an increasing function of the size of the  other class.

Starting from the elementary interactions \fer{coll1} we can easily write the kinetic equations describing the evolution in time of the pair of population densities.
In a suitable scaling, the evolution of the densities  is
quantitatively described by a system of non-Maxwellian Boltzmann-like equations \cite{MTZ}
\begin{equation}
  \label{kine}
  \begin{split}
  \frac{\partial f}{\partial t} = Q_1(f,g), \\  \frac{\partial g}{\partial t} = Q_2(g,f),
  \end{split}
\end{equation}
where the bilinear operators on the right-hand side of both equations quantify the variations of the population densities due to the elementary interactions of type \fer{coll1} where the instantaneous microscopic  transitions $x\to x^\prime$ and $y\to y^\prime$ are defined as follows
\[
\begin{split}
x^\prime &= x - \Phi_\epsilon(y)x + x \eta_\epsilon(y) \\
y^\prime &= y + \Psi_\epsilon(x)y+y\eta_\epsilon(x).
\end{split}
\]
The homogeneous Boltzmann equations \eqref{kine} can be fruitfully written in weak form. It corresponds to say that the
solution to \fer{kine} satisfies, for all smooth functions $\varphi$  the integro-differential equations
 \begin{equations}
  \label{kine-w}
\frac{d}{dt}\int_{\mathbb{R}_+}\varphi(x)f(x,t)\,dx &=
 \int_{\mathbb{R}_+^2} \kappa_1^\epsilon(x,y)  \left\langle \varphi(x^\prime)-\varphi(x)\right\rangle
f(x,t) g(y,t) \,dx \, dy, \\
\frac{d}{dt}\int_{\mathbb{R}_+} \varphi(y)g(y,t)\,dy &=
 \int_{\mathbb{R}_+^2} \kappa_2^\epsilon(x,y)  \left\langle \varphi(y^\prime)-\varphi(y)\right\rangle
f(x,t) g(y,t) \,dx \, dy.
 \end{equations}
 In \fer{kine-w} the expression $\langle\cdot\rangle$ denotes mathematical expectation with respect to the random variables $\eta_\epsilon$, $\tilde \eta_\epsilon$, and the functions $\kappa_i^\epsilon(x,y)$, $i =1,2$ are the interaction frequencies. In what follows, to maintain the connection between the kinetic model and the Lotka--Volterra system, we fix these frequencies in the form
 \begin{equation}\label{fre}
 \kappa_1^\epsilon(x,y) =\kappa^\epsilon(y)= \frac{1+y}\epsilon; \quad \kappa_2^\epsilon(x,y) = \kappa^\epsilon(x) = \frac{1+x}\epsilon.
 \end{equation}
 Within this choice, the frequency of interactions in the first equation grows with the number of predators, while in the second equation grows with the number of preys.

 By choosing $\varphi(x) =1$ in \fer{kine-w} we can easily observe that  $f(x,t)$ and $g(y,t)$ {remain probability densities} for any  time $t \ge0$, if they are so initially
 \begin{equation}\label{m0}
\int_{\mathbb{R}_+} f(x, t)\,dx = \int_{\mathbb{R}_+} g(y,t)\,dy = 1 .
 \end{equation}

Then, by choosing $ \varphi(x) = x$ in both equations \fer{kine-w} one easily obtains that the evolution of the mean values is defined by the coupled system
\begin{equation}\label{mean2}
\begin{split}
\frac{dm_f(t)}{dt} &= -\beta \,m_f(t)m_g(t); \\
 \frac{dm_g(t)}{dt} &= - \gamma\mu\, m_g(t) + \gamma \,m_f(t)m_g(t).
\end{split}
\end{equation}
Thus, if the functions $\Phi_\epsilon$ and $\Psi_\epsilon$ in \fer{coll1} are given by \fer{Phi}, and the collision frequencies in \fer{kine-w} are equal to \fer{fre}, the mean values of the densities satisfy the Lotka--Volterra system \fer{LV} in which the rate of birth $\alpha$ of the preys is set equal to zero.

We may obtain a kinetic interpretation of the birth rate by complementing the kinetic system \fer{kine} with an additional operator  which accounts for the variation of densities due to births of individuals in the groups. As already discussed in the introduction, we will resort to the redistribution operator introduced in \cite{BST}. In this paper, wealth redistribution mimics the action of a government that use the proceeds of taxation to improve the financial situation of agents based on their income.

According to \cite{BST}, for a given probability density function $h(x,t)$, with mean $m(t)$, we  assume the birth operator of the form
 \begin{equation}\label{redis}
 R_\chi^\alpha(h)(x,t) = \alpha \frac\partial{\partial x} \left[ \left(\chi x - (\chi + 1)m(t)\right)h(x,t)\right].
  \end{equation}
 The weight factor multiplying the distribution function inside the square brackets in \fer{redis}  involves in the mechanism the moments of order zero and one. Such a weight function adds one real parameter $\chi$ to the dynamics. This new parameter characterizes the type of redistribution.  The parameter $\alpha>0$ is determined by the constraint that, while the redistribution operator preserves the total mass, the mean number of individuals increases by a factor $\alpha$. In fact, for each value of the constant $\chi\in \mathbb R$, we have
 \begin{equation}\label{res1}
  \int_{\mathbb{R}_+}x R_\chi^\alpha(h)(x,t) \, dx = \alpha\, m(t),
 \end{equation}
and, provided $h(x,t)$ satisfies in addition the condition $h(0,t) = 0$, we also get
 \begin{equation}\label{res0}
  \int_{\mathbb{R}_+}R_\chi^\alpha(h)(x,t) \, dx = 0.
 \end{equation}
We can observe that the condition $h(0,t) = 0$ is not necessary in the special case $\chi = -1$, which reproduces an
anti-drift operator introduced by Slanina in  \cite{Sla} to mimic wealth redistribution effects.

As shown in \cite{BST}, the effects of the operator $R_\chi^\alpha$ on the birth process on a group described by the density function $f$ are closely related to the values of the parameter $\chi$. The case in which the percentage of births is equally distributed independently of the size of the group  is achieved with the special option $\chi=0$.
In all other cases the births are selective, and may correspond to some partition strategy.   {In the present case, as far as the preys are concerned,  it is suitable to take into account positive values of the parameter $\chi$, which correspond to impose that the percentage  of births is inversely proportional to the size of the group.  On the contrary, we assume $\chi=0$ for the predators group.}

Taking into account the previous discussion, the full system is obtained by adding to equations \fer{kine-w} the birth operators defined in \fer{redis}. Therefore, we have
 \begin{equations}
  \label{kine-wd}
\frac{d}{dt}\int_{\mathbb{R}_+}\varphi(x)f(x,t)\,dx = & \int_{\mathbb{R}_+}\varphi(x)R_\chi^\alpha(f)(x,t)\,dx \, +\\
+ &  \int_{\mathbb{R}_+^2} \kappa_1^\epsilon(x,y) \left\langle \varphi(x^\prime)-\varphi(x)\right\rangle
f(x,t) g(y,t) \,dx \, dy, \\
\frac{d}{dt}\int_{\mathbb{R}_+}\varphi(y)g(y,t)\,dy =& \int_{\mathbb{R}_+}\varphi(y)R_{0}^\nu(g)(y,t)\,dy\, + \\
+ &  \int_{\mathbb{R}_+^2} \kappa_2^\epsilon(x,y) \left\langle \varphi(y^\prime)-\varphi(y)\right\rangle
f(x,t) g(y,t) \,dx \, dy.
 \end{equations}
In \fer{kine-wd}, the constant $\nu$ which defines the percentage of the mean number of births of the group of predators, is assumed to satisfy the bound $ \nu < \gamma\mu$. In this case, resorting to property  \fer{res1}, it is immediate to conclude that the evolution of the mean numbers $m_f(t)$ and $m_g(t)$ satisfies the Lotka-Volterra system \fer{LV}, where now $\delta = \gamma\mu -\nu >0$.

\section{The limit of grazing predation rate}\label{sec:3}

We investigate now the situation in which most of the interactions between species correspond to a very small variation of the sizes of the groups, but at the same time it is possible to keep trace at the macroscopic level of all phenomena affecting predation rules. This kind of asymptotic analysis is referred to as \emph{quasi-invariant limit}, and it allows to prove that under proper assumptions on the interaction parameters,  the kinetic Boltzmann system may be approximated by a system of Fokker-Planck equations, allowing in various cases an analytical investigation of the shape of the solution. This limit corresponds to consider in \fer{kine-w} a Taylor expansion up to the second order of the differences
\[
\varphi(x^\prime)-\varphi(x), \quad \varphi(y^\prime)-\varphi(y),
\]
and subsequently to consider $\epsilon \to 0$, see \cite{PT,T} for further details. It can be proven that, in this asymptotic regime, the kinetic system \fer{kine-w} is well-approximated by the following system of  Fokker-Planck equations
 \begin{equations}
  \label{FP}
\frac{\partial f(x,t)}{\partial t} &= \frac{\partial}{\partial x}\Bigg\{ \frac{\sigma m_g(t)}2 \frac{\partial(x^2 f(x,t))}{\partial x}+ \\
 & +  \left[ (\beta m_g(t) +\alpha\chi) x - \alpha(\chi +1)m_f(t)\right]f(x,t)\Bigg\} \\
\frac{\partial g(y,t)}{\partial t} &= \frac{\partial}{\partial y}\Bigg\{\frac{\tilde\sigma m_f(t)}2 \frac{\partial(y^2 g(y,t))}{\partial y} + \\ & +  \left[ \gamma(\mu -m_f(t))y -\nu m_g(t)\right]g(y,t)\Bigg\}.
 \end{equations}
complemented by no-flux boundary conditions \cite{FPTT}.
We can observe that the evolution of the mean values of \fer{FP} is consistent with the one obtained at the Boltzmann level  in \fer{kine-w}, so that the mean values $m_f(t)$, $m_g(t)$ satisfy the Lotka--Volterra system \fer{LV}.

\section{Numerical tests}\label{sec:4}

In this section we provide several numerical tests to show the consistency of the proposed approach with the classical Lotka-Volterra dynamical system. Furthermore, we will show numerically how the kinetic approach allows additional information to be obtained on the large time behaviour of the densities of preys and predators in dependence of the parameters in the birth operators $R_\chi^\alpha(\cdot)$, $R_0^\nu(\cdot)$. For the approximation of the system of Fokker-Planck equations in \fer{FP} we will adopt the structure preserving numerical strategy developed in \cite{PZ}. These methods are capable to reproduce large time statistical properties of the solution density with arbitrary accuracy together with the preservation of its main physical properties, like its positivity.  We complement  the model \fer{FP}  with no-flux boundary conditions.

In all the following tests we will fix the parameter as reported in Table \ref{tab:1}.

\begin{table}
\begin{center}
\begin{tabular}{ |c| c| c| }
\hline
Parameter & Value & Meaning \\
\hline\hline
 $\alpha$ & 1.00 & Growth rate preys \\
 $\beta$ & 0.50 & Predation rate preys \\
 $\mu$ & 10.00 & Lowest size preys   \\
 $\gamma$ &0.15 & Growth rate predators \\
  $\nu$ & 1.00 & Birth rate predators\\
 $\delta$ & $\gamma\mu-\nu>0$ & Death rate predators\\
 $\sigma^2$ & $10^{-3}$ & Diffusion coefficient preys \\
 $\tilde \sigma^2$ & $10^{-3}$ & Diffusion coefficient predators \\
 \hline
\end{tabular}
 \caption{Parameters chosen in the numerical tests. }
 \label{tab:1}
 \end{center}
\end{table}

\subsection{Test 1. Evolution of the mean values}

We consider the Fokker-Planck model \fer{FP} with initial distribution
\begin{equation}\label{eq:fg0}
\begin{split}
f(x,0) =
\begin{cases}
C_f \exp\left (-\dfrac{(x-m_f^0)^2}{2\kappa_f}\right) & x \ge 0 \\
0 & x<0
\end{cases} \\
g(y,0) =
\begin{cases}
C_g \exp\left (-\dfrac{(y-m_g^0)^2}{2\kappa_g}\right) & y \ge 0 \\
0 & y<0,
\end{cases}
\end{split}
\end{equation}
with $C_f,C_g>0$ normalization constants and with $m_f = 4$, $m_g = 3$ and $\kappa_f = \kappa_g = \frac{1}{10}$.

In Figure \ref{fig:evolution} we show the evolution of the mean values $m_f$, $m_g$ of \fer{FP} over the time horizon $[0,15]$ for an increasing number of gridpoints $N = 201,401,801$ for the discretization of the interval $[0,L]$, $L= 10$. The choice of parameters has been defined in Table \ref{tab:1} and we further considered $\chi = 0$. We highlight in red the numerical solution of the Lotka-Volterra model obtained in the same regime of parameters. We can observe how, consistently with the analytical considerations in Section \ref{sec:3}, we can obtain numerically a good approximation of the original dynamical system in terms of the mean values of the distributions of preys $f(x,t)$ and predators $g(y,t)$.

\begin{figure}
\centering
\includegraphics[scale = 0.325]{figure/evo_mf}
\includegraphics[scale = 0.325]{figure/evo_mg}
\caption{\textbf{Test 1}. Comparison of the evolution of the mean values for the populations of preys $f(x,t)$ and predators $g(y,t)$ with the Lotka-Volterra dynamical system \fer{LV}. The Fokker-Planck system has been approximated by means of a 4th order SP methods with RK4 time integration and an increasing number of gridpoints $N = 201,401,801$ over the domain $[0,10]$. We considered $\Delta t = \Delta x^2= \Delta y^2$. }
\label{fig:evolution}
\end{figure}

\subsection{Test 2. Emergence of fat-tailed distributions}

In this test we show the evolution of the approximated distributions $f(x,t)$ and $g(y,t)$ solution to \fer{FP} for varying values of the parameters $\chi \in \mathbb R$. All the other parameters have been chosen in agreement with Table \ref{tab:1}. We consider the initial distributions defined in \fer{eq:fg0} and we solve the system of Fokker-Planck equations \fer{FP} with $N = 801$ gridpoints in the interval $[0,10]$, and time frame $[0,15]$ with $\Delta t = \Delta x^2$.

In Figure \ref{fig:surf} we provide the total evolution of the distribution functions in the cases $\chi = -1$ (top row) and $\chi = 1$ (bottom row). As specified in Section \ref{sect:2} the two choices are coherent with two different redistribution policies. Indeed, for all $\chi >0$ the number of new births is distributed in small values of $x\ge 0$, whereas for $\chi<0$ the new births are concentrated in high values of $x \ge 0$. Coherently with \cite{BST}, we may observe how the preys population may exhibit a rather different behaviour for large times. Indeed, fat-tail-type distributions appears in the case $\chi = -1$.

To clarify this behaviour, we present in Figure \ref{fig:3times} a plot of the distributions of preys and predators at three fixed times $t = 0,1,8$. We can observe how for $\chi= -1$ the distribution of preys develops fat tails in finite time. This behaviour is further observable in loglog scale, see Figure \ref{fig:loglog}.

\begin{figure}
\centering
\includegraphics[scale = 0.325]{figure/f_chim1}
\includegraphics[scale = 0.325]{figure/g_chim1} \\
\includegraphics[scale = 0.325]{figure/f_chip1}
\includegraphics[scale = 0.325]{figure/g_chip1}
\caption{\textbf{Test 2}. Evolution of the distributions $f(x,t)$ and $g(y,t)$ over the time interval $[0,15]$ and for $x,y \in [0,10]$. Approximation of the system of Fokker-Planck equations obtained with 4th order SP scheme with RK4 time integration, $N = 801$ gridpoints and $\Delta t = \Delta x^2 = \Delta y^2$. We considered $\chi = -1$ (top row) or $\chi = 1$ (bottom row) in \fer{FP}.}
\label{fig:surf}
\end{figure}

\begin{figure}
\includegraphics[scale = 0.325]{figure/f_time_chim1}
\includegraphics[scale = 0.325]{figure/g_time_chim1}\\
\includegraphics[scale = 0.325]{figure/f_time_chi1}
\includegraphics[scale = 0.325]{figure/g_time_chi1}
\caption{\textbf{Test 2}. Distributions of preys $f(x,t)$ (left) and predators $g(y,t)$ (right) at three different time steps $t = 0,1,8$, $x,y \in [0,10]$. We considered $\chi = -1$ (top row) or $\chi = 1$ (bottom row) in \fer{FP}.}
\label{fig:3times}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale = 0.325]{figure/semilog_fg_chim1}
\includegraphics[scale = 0.325]{figure/semilog_fg_chip1}
\caption{\textbf{Test 2}. Distributions of preys $f(x,t)$ and predators $g(y,t)$ at time $T = 15$ in the cases $\chi = -1 $ (left) or $\chi = 1$ (right). }
\label{fig:loglog}
\end{figure}

\section{Conclusions}

Motivated by the similarities between the problem of wealth redistribution with taxation in a multi-agent society, we introduced in this paper a kinetic description of  the statistical distributions in time of the sizes of groups of biological systems obeying  a prey-predator dynamic. The kinetic system is built in such a way that the evolution of the mean values is driven by a classical Lotka-Volterra system of differential equations. Among other characteristics, it is shown that the time evolution of the probability distributions of the size of groups of the two interacting species is heavily dependent on the  kinetic redistribution operator quantifying the births. Numerical experiments clarify the time-behavior of the distributions of groups of the species, and highlight the possibility that the distribution of the sizes of preys could exhibit fat tails.

 \vskip 2cm
\vspace{0.5cm} \indent {\it
A\,c\,k\,n\,o\,w\,l\,e\,d\,g\,m\,e\,n\,t\,s.\;} This work has been written within the activities of GNFM group of INdAM (National Institute of
High Mathematics). GT wish to acknowledge partial support by IMATI, Institute for Applied Mathematics and Information Technologies ``Enrico Magenes'', Via Ferrata 5 Pavia, Italy. MZ acknowledges the support of MUR-PRIN2020 Project No.2020JLWP23 (Integrated Mathematical Approaches to Socio–Epidemiological Dynamics).

\bigskip
\begin{center}


\section{Introduction}
Memristors, also known as resistive random access memories (RRAM) are a new type of memory technology fully embeddable in CMOS, providing a compact nonvolatile, and fast memory \cite{zidan2018future}.
These devices provide fantastic opportunities to integrate logic and memory tightly and allow low-power computing, in particular for Artificial Intelligence models and neuromorphic computing \cite{yu2018neuro}.
Unfortunately, the behavior of memristors is highly complex and partly stochastic \cite{majumdar2021model}: device models do not provide an accurate prediction of their behavior.
It is therefore essential to prototype computing concepts involving memristors experimentally.
However, appropriate platforms are extremely complex to fabricate due to the need to co-integrate commercial CMOS and  memristor devices on the same die. In this work, we designed, fabricated, and tested a prototyping platform, associating an array of 8,192 hafnium-oxide-based memristors and a collection of CMOS periphery circuits.
Our platform is multi-paradigm, which permits prototyping a wide range of both digital and analog projects.

\vspace{-0.25cm}
\section{Description of the Die}
A photograph and a layout view of our integrated circuit are presented in Figs.~1a-b. An electron microscopy image of a memristor in the backend of line of our hybrid memristor/CMOS process is shown in Fig.~1d. A commercial foundry fabricated the CMOS part (including the backend up to metal layer 4), using a 130-nanometer process. Afterward, we deposited the memristors on top of metal 4 using atomic layer deposition, and a fifth layer of metal.

Our integrated circuit embeds periphery circuitry enabling the use of memristors within two modes. Fig.~1c shows a simplified schematic of the circuit. It uses consistent color codes: blue-colored blocks are digital-mode circuits, designed using thin oxide low-power transistors and supplied by digital nominal voltage (except for level shifters), and orange-colored blocks are analog-mode circuits, designed using thick oxide transistors to be compatible with high voltages.

The digital mode circuits (Figs.~1h-k) consist of: row and column decoders to select devices based on input addresses, level shifters on each row and column, which shift digital nominal voltage to higher voltages required to form and program memristors, and  precharge sense amplifiers, with a logic-in-memory feature, at each column  \cite{zhao2014synchronous}. These sense amplifiers allow reading the binary states of memory cells in a highly energy-efficient fashion while optionally performing logic operations at the same time. The complementary approach of \cite{bocquet2018} is used in our array for reducing the bit error rate.

When activating the analog mode (Figs.~1g), digital circuits are deactivated and the memristors array connections are switched to the analog circuitry. In this mode, shift registers configure input multiplexers permitting direct access to the analog state  of memristors, using low-resistance transmission gates. Each word line, bit line, and source line is then connected to the ground or to one of two analog InOut Pads, which can be connected to external equipment, e.g., Keysight B1530, a pulse source and measurement unit widely used to characterize memory devices.

The memristor array and all analog and mixed-signal circuits were designed in a full custom fashion, based on an extensive work on memristor characterization, modeling, and simulation. All digital circuits were placed and routed automatically using an HDL description and a Cadence Encounter flow provided by the foundry. Then, all circuits of the system were assembled manually and routed automatically using a Cadence Encounter flow developed in-house using a homemade abstract view of the memory array.

\vspace{-0.25cm}
\section{Uses of the Platform}

To make the system re-configurable for different projects, we developed the  experimental setup of Fig.~1l: a PCB routes a microcontroller unit and  measurement equipment with our packaged die. Python scripts control the measurement.

Optimizing read and programming strategies using the digital mode can allow the successful implementation of digital applications. Memristors feature a complex interplay between programming energy, reading speed, read disturb effects, and device endurance,  which our platform allows understanding. Fig.~1m shows an endurance study example. A memristor is programmed repeatedly using the digital mode circuits, and the memristor resistance is checked regularly using the analog mode and reported in Fig.~1m. It shows that the memristor resistance starts to degrade after $10^9$ cycles, concurrently with the emergence of bit errors seen by the reads after each programming using sense amplifiers (not shown in Figures). We observed that memristor endurance can vary between $10^3$ and $10^9$ cycles depending on programming conditions.

The analog mode of the platform can be used to prototype computing concepts where memristors are used in an analog fashion, e.g., as artificial synapses in machine learning or neuromorphic circuits \cite{yu2018neuro}. Fig.~1e shows measurements on a memristor in our platform when applying a succession of 15,000 1$V$ 1.5-$\mu$s programming pulses: the memristor resistance progressively increases, a feature that permits the memristor to implement a synaptic learning rule. This use is particularly appealing due to its compactness, but the imperfections of memristors (thermal and random telegraph noise,  cycle-to-cycle, and device-to-device variability) pose challenges that make it necessary to test ideas experimentally. Our platform supports prototyping various neuromorphic experiments, targeting inference, deterministic or probabilistic learning \cite{dalgaty2021situ}.

\vspace{-0.25cm}
\section{Conclusion}
We have designed, fabricated, and tested a flexible multi-paradigm platform to prototype and optimize digital and/or analog computing concepts, based on a hybrid CMOS/memristor integrated circuit.
We are currently using it to validate multiple digital logic-in-memory and analog neuromorphic concepts, and plan to make the platform available to other research groups.

\vspace{-0.25cm}

\begin{acks}
\vspace{-0.1cm}
This work was supported by  ERC Grant NANOINFER (715872) and ANR grant NEURONIC (ANR-18-CE24-0009).
\end{acks}
\vspace{-0.25cm}
\bibliographystyle{ACM-Reference-Format}