# Kinetic models for systems of interacting agents with multiple microscopic
states

Marzia Bisi Department of Mathematical, Physical and Computer Sciences,
University of Parma, Italy<EMAIL_ADDRESS>Nadia Loy Department of
Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca
degli Abruzzi 24, 10129 Torino, Italy<EMAIL_ADDRESS>
###### Abstract

We propose and investigate general kinetic models with transition
probabilities that can describe the simultaneous change of multiple
microscopic states of the interacting agents. These models can be applied to
many problems in socio-economic sciences, where individuals may change both
their compartment and their characteristic kinetic variable, as for instance
kinetic models for epidemics or for international trade with possible
transfers of agents. Mathematical properties of our kinetic model are proved,
as existence and uniqueness of a solution for the Cauchy problem in suitable
Wasserstein spaces. The quasi-invariant asymptotic regime, leading to simpler
kinetic Fokker-Planck-type equations, is investigated and commented on in
comparison with other existing models. Some numerical tests are performed in
order to show time evolution of distribution functions and of meaningful
macroscopic fields, even in case of non-constant interaction probabilities.

Keywords: Boltzmann equation, Markov process, multi-agent system, socio-
economic modelling

##  1 Introduction

In the literature of kinetic models for multi-agent systems, there is an
increasing interest in phenomena where the agents are characterized by a
multiple microscopic state and they are, in particular, divided into
subpopulations.

Kinetic theory for multi-agent systems has its roots in the classical kinetic
theory related to the Boltzmann equation for the description of a rarefied gas
[11], in which individuals are molecules identified by a microscopic state $v$
that is the velocity and that changes because of binary interactions. The
classical kinetic theory for gas dynamics has been generalized to various
kinds of interacting systems, where the microscopic state is not necessarily
the velocity, providing reasonable mathematical models for many socio-economic
problems, as the evolution of wealth distribution [13, 40], the opinion
formation [38, 39], the pedestrian or vehicular traffic dynamics [22, 23],
birth and death processes [25, 30] and many others. Also in this field, models
describing the interaction of different populations through a system of
Boltzmann equations have been proposed for instance in [21, 19] for wealth
exchanges, in [20] for opinion formation in presence of leaders, in [8] for
multilane traffic models.

The classical Boltzmann equation has been extended some decades ago to
mixtures of different gaseous species [12, 27], even in presence of chemical
reactions [24], and also consistent BGK approximations have been proposed and
investigated [2, 5, 34]. In these models, when describing bimolecular chemical
reactions, a given binary interaction between molecules may simultaneously
cause both the change of the velocity and transfers of the involved molecules
to different species. Also models in which the molecule is characterized by
the belonging to a species, a molecular velocity and another inner variable
modelling the internal energy have been proposed [9].

In the field of kinetic equations for multi-agent systems applied to socio-
economic phenomena, a very interesting class of models in which individuals
have a multiple microscopic state corresponds to models in which the total
population is divided into subgroups and each agent is also described by a
physical quantity (wealth, opinion, viral load, etc.). Each group is then
characterized by a distribution function depending on the given microscopic
physical quantity that can be exchanged both with individuals of the same
subgroup and with individuals of a different subgroup, according to suitable
interaction rules. For simplicity, interactions causing exchanges of goods and
the ones giving rise to a change of subgroup of one agent are often modelled
separately, by means of different kinetic operators [19, 21, 20, 18, 31, 17].
Formally, models belonging to this class have been derived in [30], where the
authors describe multi-agent systems in which the agents are characterized by
a double microscopic state: a physical quantity $v$, changing because of
binary interactions, and a label $x$ denoting the subgroup of the agent,
changing as a consequence of a Markovian process. The two stochastic processes
for the evolution of $v$ and $x$ are independent and occur with different
frequencies, giving thus rise to different operators, where the one relevant
to the variation of $v$ may be the sum of inter- or intra- species
interactions.

However, a class of models worth to be investigated is the one in which the
agent changes the microscopic quantity $v$ and the subgroup simultaneously as
a consequence of the same binary interaction. An example is given by the
aforementioned bimolecular chemical reactions in gas mixtures. The classical
kinetic description for chemically reacting gases was introduced in the
pioneering work [36] where each gas species has a different distribution
function of the microscopic velocity of its molecules. A pair of molecules
belonging to potentially different species may interact exchanging both the
velocity and the species. Such bimolecular reaction is described by means of a
collisional operator that involves both distributions of reactants and
products, and the change of the velocity is included by means of a
transformation with unit Jacobian like in the classical strong Boltzmann
equation. In the literature of kinetic theory for socio-economic sciences, the
recent paper [4] describes the trade among different subpopulations, living in
different countries, taking into account also possible transfers of
individuals from one country to another, by means of suitable Boltzmann-type
operators similar to the ones modelling bimolecular chemical reactions in gas
mixtures. In this model, then, binary interactions between individuals may
lead to both an exchange of wealth and to a transfer to another subpopulation
as a consequence of the same binary interaction, i.e. the microscopic state
identifying the wealth and the one related to the label, that denotes the
subpopulation, change simultaneously.

Another topic that has gained much interest in recent years also in kinetic
theory, is the modelling of the spread of an epidemic: in this respect,
compartmental Boltzmann models allowing the passage of individuals from an
epidemiological compartment to another have been proposed, essentially of SIR
type, where a susceptible individual could become infected and then removed
because of healing or death [18, 15, 16]. A different kinetic description of
infectious diseases consists in modelling interactions among different types
of human cells, including the immune cells [31, 17]. For example, in some
epidemic models, susceptible (carrying a vanishing viral load $v$) and
infected individuals (with $v>0$) interact exchanging the quantity $v$, and as
soon as the susceptible individual’s viral load becomes positive because of
the binary interaction, then he/she becomes infected [15, 16]. In these works,
the authors, similarly to [30], start from a microscopic description in which
each agent is characterized by the microscopic quantity $v$ and by the label
$x$ denoting, in this case, the compartment. The microscopic dynamics is then
described through discrete in time stochastic processes in which the new
microscopic physical quantity and label are modelled through Markov-type jump
processes governed by suitable transition probabilities [15, 16]. As a
consequence, the kinetic model implementing the prescribed microscopic
dynamics is a kinetic equation with an operator in which the kernel is related
to the transition probability.

This formulation of the Boltzmann equation is well known also in kinetic
theory for a single gas. Indeed, besides the classical Boltzmann operator in
which the kernel has a proper cross section taking into account the
intermolecular potential (depending on the relative speed and on the impact
angles) [11], other different forms have been used in the literature. The most
common one is the so–called Waldmann representation [42], showing in the
kernel the probability distribution of the collision process transforming pre-
collision velocities $(v,w)$ into the post-collision ones
$(v^{\prime},w^{\prime})$, and the Boltzmann integral over the unit sphere is
replaced by integrals over post-collision velocity variables. An analogous
scattering kernel formulation replaces Waldmann kernel by its integral over
the velocity of the partner molecule [37]. The equivalence between these
kinetic equations has been proved in [7] for microscopic interactions that
conserve the average and the energy. In the literature of kinetic equations
for socio-economic sciences, the authors in [30] show the equivalence between
the collision-like Boltzmann equation and Markovian jump-processes described
by transition probabilities that can be related to the Waldmann
(probabilistic) representation of the Boltzmann equation.

In this paper we will present kinetic models for socio-economic problems in
which agents have a multiple microscopic state, starting from a microscopic
stochastic process ruled by transition probabilities that allow to describe
the simultaneous change of all the microscopic variables and by a microscopic
state-dependent interaction frequency. In the case in which the agent state is
given by a microscopic quantity $v$ and by a label $x$ denoting the subgroup,
we will show that this approach has some advantages with respect to the
classical collision model [4]. Indeed, as explained also in [4], the
construction of the gain term of the Boltzmann operators requires the
invertibility of the collision process, that obviously holds in gas-dynamics
(because of conservation of total momentum and energy in each collision), but
not in human interactions, that are also influenced by non-deterministic
(random) effects. This invertibility property is not needed in the operator
with a transition probability in the kernel, because each single interaction
has its own probability, not related to the reverse process. Moreover,
realistic situations with non constant interaction probabilities are easier to
manage in the stochastic Boltzmann formulation, therefore this approach could
have many applications in kinetic modelling of social sciences. For these
reasons in this paper we present a formal and organic treatment of kinetic
equations involving a microscopic stochastic process that simultaneously
changes several internal states of the interacting agents (typically, their
compartment and the value of their kinetic variable). In more detail, the
paper is organized as follows.

In Section 2 we formally derive the general form of a kinetic model
implementing a microscopic dynamics in which each agent is characterized by a
set $\boldsymbol{z}\in\Omega\subset\mathbb{R}^{d}$ of microscopic states which
may change simultaneously in each binary interaction, that is described by a
transition probability and ruled by a microscopic state-dependent frequency.
Even though the procedure is quite classical, stating the discrete in time
stochastic process will be useful for writing the Nanbu-Babovski Algorithm for
simulating the microscopic dynamics. Then we revise and establish the relation
with some well-known models such as the collision-like Boltzmann equation and
the kinetic equations describing transfers among different groups due to
binary interactions. Finally, we explicitly derive the kinetic equation for a
multi-agent system in which a binary interaction causes simultaneously a
transfer and an exchange of a microscopic quantity. In Section 3, mathematical
properties of the Cauchy problem associated to our general Boltzmann equation
are discussed, proving existence and uniqueness of a solution in suitable
Wasserstein spaces. Then, the quasi-invariant limit commonly used to
investigate socio-economic kinetic models is adapted to our general frame,
allowing to derive suitable Fokker-Planck equations with additional terms
taking into account transfers of agents. Section 4 is devoted to the
investigation of a specific kinetic model fitting into our general framework,
describing international trade with possible transfers of individuals:
evolution of number density and mean wealth of each country are computed from
the kinetic model, the quasi-invariant limit is performed, and analogies and
differences with respect to analogous models for a single population [13] are
discussed, with particular reference to the Pareto index of steady
distributions. In Section 5 we show some numerical tests, simulating our
kinetic equations by means of a Nanbu-Babovski Monte Carlo algorithm
implementing the discrete in time stochastic process presented in Section 2:
the evolution of distribution functions and of macroscopic quantities are
commented on for varying parameters. Section 6 contains some concluding
remarks and perspectives.

##  2 Kinetic models for binary interactions processes

In this section, we provide a formal derivation of kinetic equations
implementing binary interactions among agents whose microscopic state is a
vector, described by means of Markovian processes, where the interaction
frequency depends on the microscopic state of the interacting agents. Then, we
illustrate the relation to well-known kinetic models for binary interactions
leading to exchange of physical quantities (the collision-like Boltzmann
equation) and to label-switch processes, also named transfers. We eventually
present a general framework for describing, through kinetic equations,
microscopic binary interaction processes leading to both exchanges of a
physical quantity and label-switches.

###  2.1 Kinetic models with transition probabilities

Let us consider a large system of agents described by a microscopic state
$\boldsymbol{z}\in\Omega\subset\mathbb{R}^{d}$. We shall suppose that the
change of the microscopic state (of all its components simultaneously) is due
to stochastic binary interactions. A probabilistic description of such
interactions may be given by means of transition probability functions

$T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y})>0,\qquad\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})>0\qquad\forall\boldsymbol{z},\boldsymbol{y}\in\Omega,\
t>0,$ (1)

namely the conditional probabilities that, given a binary interaction between
an agent $\boldsymbol{z}$ and an agent $\boldsymbol{y}$, the first changes
into $\boldsymbol{z}^{\prime}$ while the second into
$\boldsymbol{y}^{\prime}$, respectively. Such a microscopic description may be
assimilated to a Markov-type jump process. In order for
$T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y}),\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})$
to be conditional probability densities, they have to satisfy the following
further property:

$\int_{\Omega}T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y})\,d\boldsymbol{z}^{\prime}=1,\quad\int_{\Omega}\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})\,d\boldsymbol{y}^{\prime}=1\qquad\forall\boldsymbol{z},\boldsymbol{y}\in\Omega,\
t>0.$ (2)

The binary interactions may happen with a frequency
$\lambda_{\boldsymbol{z}\boldsymbol{y}}$, namely the frequency of the binary
interactions between two agents having microscopic states
$\boldsymbol{z},\boldsymbol{y}$ depends on the microscopic states themselves.
We remark that the two transition probabilities $T$ and $\tilde{T}$ are given
in order to take into account for possible asymmetries in the binary
interactions. We remark that the symmetry of the binary interactions is here
expressed by

$T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y})=\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y}),$
(3)

and
$\lambda_{\boldsymbol{z}\boldsymbol{y}}=\lambda_{\boldsymbol{y}\boldsymbol{z}}$.
As classically done [32], a kinetic description of the multi-agent system can
be derived by introducing discrete in time stochastic processes. Let
$\boldsymbol{Z}_{t},\boldsymbol{Y}_{t}\in\Omega$ be random variables
describing the microscopic state of two agents at time $t>0$. Let
$f=f(\boldsymbol{z},t)$ be the probability density function associated to our
multi-agent system, i.e. the probability density function of the random
variable of a given agent $\boldsymbol{Z}_{t}$, thus satisfying

$\int_{\Omega}f(\boldsymbol{z},t)\,d\boldsymbol{z}=1.$ (4)

During a sufficiently small time $\Delta{t}>0$ the agents may or may not
change their state $\boldsymbol{Z}_{t},\boldsymbol{Y}_{t}$ depending on
whether a binary interaction takes place or not. We express this discrete-in-
time random process as

$\displaystyle\boldsymbol{Z}_{t+\Delta{t}}$
$\displaystyle=(1-\Theta)\boldsymbol{Z}_{t}+\Theta\boldsymbol{Z}_{t}^{\prime},$
(5) $\displaystyle\boldsymbol{Y}_{t+\Delta{t}}$
$\displaystyle=(1-\Theta)\boldsymbol{Y}_{t}+\Theta\boldsymbol{Y}_{t}^{\prime},$

where $\boldsymbol{Z}_{t}^{\prime},\,\boldsymbol{Y}_{t}^{\prime}$ are random
variables describing the new microscopic state of $\boldsymbol{Z}_{t}$ and
$\boldsymbol{Y}_{t}$ respectively after a binary interaction and having joint
probability density functions
$g=g(\boldsymbol{Z}_{t}^{\prime}=\boldsymbol{z}^{\prime};\boldsymbol{Z}_{t}=\boldsymbol{z},\boldsymbol{Y}_{t}=\boldsymbol{y}),\tilde{g}=\tilde{g}(\boldsymbol{Y}_{t}^{\prime}=\boldsymbol{y}^{\prime};\boldsymbol{Z}_{t}=\boldsymbol{z},\boldsymbol{Y}_{t}=\boldsymbol{y})$,
while $\Theta\in\\{0,\,1\\}$ is a Bernoulli random variable, which we assume
to be independent of all the other random variables appearing in (5),
discriminating whether a binary interaction takes place ($\Theta=1$) or not
($\Theta=0$) during the time $\Delta{t}$. In particular, we set the
probability to change the microscopic state

$\operatorname{Prob}(\Theta=1)=\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\Delta{t},$
(6)

where $\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}$ is the interaction
frequency between agents with microscopic states $\boldsymbol{Z}_{t}$ and
$\boldsymbol{Y}_{t}$. Notice that, for consistency, we need
$\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\Delta{t}\leq 1$.

Let now $\phi=\phi(\boldsymbol{z})$ be an observable quantity defined on
$\boldsymbol{z}\in\Omega$. From (5)-(6), together with the assumed
independence of $\Theta$, we see that the mean variation rate of $\phi$ in the
time interval $\Delta{t}$ satisfies

$\displaystyle\frac{\langle\phi(\boldsymbol{Z}_{t+\Delta{t}})\rangle-\langle\phi(\boldsymbol{Z}_{t})\rangle}{\Delta{t}}+\frac{\langle\phi(\boldsymbol{Y}_{t+\Delta{t}})\rangle-\langle\phi(\boldsymbol{Y}_{t})\rangle}{\Delta{t}}=$
$\displaystyle=\frac{\langle(1-\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\Delta
t)\phi(\boldsymbol{Z}_{t})\rangle+\Delta
t\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Z}_{t}^{\prime})\rangle-\langle\phi(\boldsymbol{Z}_{t})\rangle}{\Delta{t}}$
$\displaystyle+\frac{\langle(1-\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\Delta
t)\phi(\boldsymbol{Y}_{t})\rangle+\Delta
t\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Y}_{t}^{\prime})\rangle-\langle\phi(\boldsymbol{Y}_{t})\rangle}{\Delta
t}$

where $\langle C_{t}\rangle$ denotes the average of the random variable
$C_{t}$ with respect to its probability density function. Whence, we deduce
the instantaneous time variation of the average of $\phi$ in the limit
$\Delta{t}\to 0^{+}$ as

$\frac{d}{dt}\langle\phi(\boldsymbol{Z}_{t})\rangle=\dfrac{1}{2}\Big{(}\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Z}_{t}^{\prime})\rangle+\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Y}_{t}^{\prime})\rangle-\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Z}_{t})\rangle-\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Y}_{t})\rangle\Big{)}$
(7)

where we used the fact that
$\langle\phi(\boldsymbol{Z}_{t})\rangle=\langle\phi(\boldsymbol{Y}_{t})\rangle$
that implies
$\langle\phi(\boldsymbol{Z}_{t})\rangle+\langle\phi(\boldsymbol{Y}_{t})\rangle=2\langle\phi(\boldsymbol{Z}_{t})\rangle$.

We now specify the gain terms as

$\displaystyle\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Z}_{t}^{\prime})\rangle=\int_{\Omega}\int_{\Omega^{2}}\phi(\boldsymbol{z}^{\prime})\lambda_{\boldsymbol{z}\boldsymbol{y}}g(\boldsymbol{z}^{\prime};\boldsymbol{z},\boldsymbol{y})\,d\boldsymbol{z}d\boldsymbol{y}\,d\boldsymbol{z}^{\prime},$
$\displaystyle\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Y}_{t}^{\prime})\rangle=\int_{\Omega}\int_{\Omega^{2}}\phi(\boldsymbol{y}^{\prime})\lambda_{\boldsymbol{z}\boldsymbol{y}}\tilde{g}(\boldsymbol{y}^{\prime};\boldsymbol{z},\boldsymbol{y})\,d\boldsymbol{z}d\boldsymbol{y}\,d\boldsymbol{y}^{\prime},$

where $g$ and $\tilde{g}$ are the joint probability density functions of
$\boldsymbol{Z}_{t}^{\prime}$ and $\boldsymbol{Y}_{t}^{\prime}$, respectively,
and of the samples of the random variables
$\boldsymbol{Z}_{t}=\boldsymbol{z},\boldsymbol{Y}_{t}=\boldsymbol{y}$ at time
$t$. The probability density functions $g$ and $\tilde{g}$ are defined as

$g(\boldsymbol{z}^{\prime};\boldsymbol{z},\boldsymbol{y})=T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y})\,f_{2}(\boldsymbol{z},\boldsymbol{y},t)\qquad\tilde{g}(\boldsymbol{y}^{\prime};\boldsymbol{z},\boldsymbol{y})=\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})\,f_{2}(\boldsymbol{z},\boldsymbol{y},t),\,$
(8)

where $f_{2}(\boldsymbol{z},\boldsymbol{y},t)$ is the joint distribution of
the couple $(\boldsymbol{z},\boldsymbol{y})$ at time $t$. As typically done in
kinetic theory, we assume propagation of chaos, i.e. $\boldsymbol{z}$ and
$\boldsymbol{y}$ are independently distributed, which allows us to perform the
factorization
$f_{2}(\boldsymbol{z},\boldsymbol{y},t)=f(\boldsymbol{z},t)f(\boldsymbol{y},t)$,
so that we can write

$g(\boldsymbol{z}^{\prime};\boldsymbol{z},\boldsymbol{y})=T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y})\,f(\boldsymbol{z},t)f(\boldsymbol{y},t)\qquad\tilde{g}(\boldsymbol{y}^{\prime};\boldsymbol{z},\boldsymbol{y})=\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})\,f(\boldsymbol{z},t)f(\boldsymbol{y},t).$

It is immediate to verify that $g$ and $\tilde{g}$ are probability density
functions thanks to (2) and (4). Analogously, the loss terms can be naturally
written as

$\displaystyle\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Z}_{t})\rangle=\int_{\Omega^{2}}\phi(\boldsymbol{z})\lambda_{\boldsymbol{z}\boldsymbol{y}}f(\boldsymbol{z},t)f(\boldsymbol{y},t)\,d\boldsymbol{z}d\boldsymbol{y},$
$\displaystyle\langle\lambda_{\boldsymbol{Z}_{t}\boldsymbol{Y}_{t}}\phi(\boldsymbol{Y}_{t})\rangle=\int_{\Omega^{2}}\phi(\boldsymbol{y})\lambda_{\boldsymbol{z}\boldsymbol{y}}f(\boldsymbol{z},t)f(\boldsymbol{y},t)\,d\boldsymbol{z}d\boldsymbol{y}.$

Therefore Eq. (7) can be stated as

$\displaystyle\frac{d}{dt}\int_{\Omega}f(\boldsymbol{z},t)\phi(\boldsymbol{z})\,d\boldsymbol{z}$
$\displaystyle=\dfrac{1}{2}\int_{\Omega}\int_{\Omega^{2}}\lambda_{\boldsymbol{z}\boldsymbol{y}}T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y})\left(\phi(\boldsymbol{z}^{\prime})-\phi(\boldsymbol{z})\right)f(\boldsymbol{z},t)f(\boldsymbol{y},t)\,d\boldsymbol{z}d\boldsymbol{y}d\boldsymbol{z}^{\prime}$
(9)
$\displaystyle+\dfrac{1}{2}\int_{\Omega}\int_{\Omega^{2}}\lambda_{\boldsymbol{z}\boldsymbol{y}}\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})\left(\phi(\boldsymbol{y}^{\prime})\,-\phi(\boldsymbol{y})\right)f(\boldsymbol{z},t)f(\boldsymbol{y},t)\,d\boldsymbol{z}d\boldsymbol{y}d\boldsymbol{y}^{\prime},$

where we have used (2) in order to write the loss terms.

In the following, we shall illustrate three meaningful examples of kinetic
models describing binary interaction processes by means of transition
probabilities: $i)$ binary interactions causing exchange of physical
quantities; $ii)$ binary interactions leading to transfers of individuals;
$iii)$ binary interactions leading to both exchange of physical quantities and
transfers of individuals.

###  2.2 Boltzmann-type description of classical binary interaction dynamics

Let us consider the case in which the microscopic state of the agent is a non-
negative physical quantity $\boldsymbol{z}=v\in\Omega=\mathbb{R}_{+}$.
Extensions to negative and possibly also bounded microscopic states are mostly
a matter of technicalities. In general, as classically done in kinetic theory
[32], if $v,w\in\mathbb{R}_{+}$ denote the pre-interaction states of any two
interacting agents, their post-interaction states $v^{\prime},w^{\prime}$ will
be given by general interaction rules in the form

$v^{\prime}=I(v,w)+D(v,w)\eta,\qquad
w^{\prime}=\tilde{I}(v,w)+\tilde{D}(v,w)\eta_{*}$ (10)

where $\eta$ and $\eta_{*}$ are independent random variables satisfying
$\langle\eta\rangle=\langle\eta_{*}\rangle=0,\langle\eta^{2}\rangle=\langle\eta_{*}^{2}\rangle=1$,
namely with zero average and unitary variance. It is known that an aggregate
description of the (sole) binary interaction dynamics inspired by the
principles of statistical mechanics can be obtained by introducing a
probability density function $f=f(v,t)\geq 0$ such that $f(v,t)dv$ gives the
proportion of agents having at time $t$ a microscopic state comprised between
$v$ and $v+dv$. Such a probability density function satisfies a Boltzmann-type
kinetic equation, which in weak form reads

$\frac{d}{dt}\int_{\mathbb{R}_{+}}f(v,t)\phi(v)\,dv=\dfrac{\lambda}{2}\langle\int_{\mathbb{R}_{+}^{2}}\big{(}\phi(v^{\prime})+\phi(w^{\prime})-\phi(v)-\phi(w)\big{)}f(v,t)f(w,t)\,dvdw\rangle$
(11)

where $\lambda$ is the interaction frequency that we here assume to be
independent of the microscopic states of the agents.

On the other hand, if we want to describe the binary interactions through
transition probabilities (1) with
$\boldsymbol{z}^{\prime}=v^{\prime},\boldsymbol{y}^{\prime}=w^{\prime}$, we
have that (9) can be rewritten as

$\displaystyle\frac{d}{dt}\int_{\mathbb{R}_{+}}f(v,t)\phi(v)\,dv$
$\displaystyle=\dfrac{\lambda}{2}\int_{\mathbb{R}_{+}^{2}}\left(\int_{\mathbb{R}_{+}}\phi(v^{\prime})T(v^{\prime}|v,w)\,dv^{\prime}-\phi(v)\right)f(v,t)f(w,t)\,dvdw$
(12)
$\displaystyle+\dfrac{\lambda}{2}\int_{\mathbb{R}_{+}^{2}}\left(\int_{\mathbb{R}_{+}}\phi(w^{\prime})\tilde{T}(w^{\prime}|v,w)\,dw^{\prime}-\phi(w)\right)f(v,t)f(w,t)\,dvdw.$

We can define the first two statistical moments of the distribution function
$f$ as:

$M(t):=\int_{\mathbb{R}_{+}}vf(t,\,v)\,dv,\qquad
E(t):=\int_{\mathbb{R}_{+}}v^{2}f(t,\,v)\,dv,$

that represent the average and energy, respectively. As done in [29] in the
symmetric case, in order to establish a relation between (11)-(10) and (12),
we investigate primarily the trend of the statistical moments as prescribed by
the two different models. Setting $\phi(v)=v,v^{2}$ in (12) and (11)-(10)
yields the evolution equations for $M,E$ for a system of agents obeying the
microscopic dynamics expressed in terms of transition probabilities (1) or
interaction rules (10), respectively. By comparing the evolution equations of
$M$ and $E$ prescribed by the two different kinetic models, we see that the
evolution is the same if we choose

$\displaystyle I(v,\,w)=V_{T}(v,\,w),\quad
D(v,\,w)=\sqrt{E_{T}(v,\,w)-V_{T}^{2}(v,\,w)}=:D_{T}(v,\,w),$
$\displaystyle\tilde{I}(v,\,w)=V_{\tilde{T}}(v,\,w),\quad\tilde{D}(v,\,w)=\sqrt{E_{\tilde{T}}(v,\,w)-V_{\tilde{T}}^{2}(v,\,w)}=:D_{\tilde{T}}(v,\,w)$

where

$V_{T}(v,\,w):=\int_{\mathbb{R}_{+}}v^{\prime}T(v^{\prime}\,|\,v,\,w)\,dv^{\prime},\qquad
E_{T}(v,\,w):=\int_{\mathbb{R}_{+}}v^{\prime
2}T(v^{\prime}\,|\,v,\,w)\,dv^{\prime}$

and

$V_{\tilde{T}}(v,\,w):=\int_{\mathbb{R}_{+}}w^{\prime}\tilde{T}(w^{\prime}\,|\,v,\,w)\,dw^{\prime},\qquad
E_{\tilde{T}}(v,\,w):=\int_{\mathbb{R}_{+}}w^{\prime
2}\tilde{T}(w^{\prime}\,|\,v,\,w)\,dw^{\prime}$

denote the mean and the energy, respectively, of $T$ and $\tilde{T}$ for a
given pair $(v,\,w)\in{\mathbb{R}_{+}}\times{\mathbb{R}_{+}}$ of pre-
interaction states, while $D_{T}(v,\,w)$ and $D_{\tilde{T}}(v,\,w)$ are the
standard deviations of $T$ and $\tilde{T}$, respectively. Therefore, if
dealing with (11) we can consider the collisions

$v^{\prime}=V_{T}(v,\,w)+D_{T}(v,\,w)\eta,\qquad
w^{\prime}=V_{\tilde{T}}(v,\,w)+D_{\tilde{T}}(v,\,w)\eta_{\ast},$ (13)

and this choice makes formulations (11) and (12) equivalent at the macroscopic
level (at least for the mass, average and energy). As highlighted in [29], in
general, (11) with (10) and (12) are not the same kinetic equation, although
with the choice (13) they account for the same evolution of the first and
second statistical moments of $f$. Nevertheless, if in (12) we take

$T(v^{\prime}|v,w)=\delta\Big{(}v^{\prime}-(V_{T}(v,w)+D_{T}(v,w)\eta)\Big{)},\qquad\tilde{T}(w^{\prime}|v,w)=\delta\Big{(}w^{\prime}-(V_{\tilde{T}}(v,w)+D_{\tilde{T}}(v,w)\eta_{\ast})\Big{)}$
(14)

where in the right-hand side $\delta$ is the Dirac delta, then we can formally
show that (12) becomes exactly (11)-(10). Of course, in this case the right
hand side (12) is meant to be written in brackets $\langle\cdot\rangle$.

###  2.3 Label switch process caused by binary interactions

Let us now consider the case in which $\Omega=\mathcal{I}_{n}=\\{1,...,n\\}$
and the microscopic discrete variable $x\in\mathcal{I}_{n}$ is regarded as a
label, that may denote the belonging of the agent to a certain group or
subpopulation. We assume that label switches, i.e. migrations across
subpopulations, can be caused by binary interactions between agents, causing a
transfer of (potentially) both of them, but in such a way that the total mass
of the agents in the system is conserved. We say that this process is formally
a Markov-type one because the probability to switch from the current labels
$x,y$ to new labels $x^{\prime},y^{\prime}$ does not depend on how the agents
reached previously the labels $x,y$. In particular we denote

$P_{xy}^{x^{\prime}y^{\prime}}:=P(x^{\prime},y^{\prime}|x,y)$ (15)

the conditional probability density function of switching to the groups
$x^{\prime},y^{\prime}$ given the pre-interaction labels $x,y$.

######  Remark 2.1.

Since the variables $x,y$ are discrete, the mapping
$(x^{\prime},y^{\prime})\mapsto P(x^{\prime},y^{\prime}|x,y)$ is a discrete
probability measure. Consequently, we actually have

$\int_{\mathcal{I}_{n}^{2}}P(x^{\prime},y^{\prime}|x,y)\,dx^{\prime}dy^{\prime}=\sum_{x^{\prime},y^{\prime}\in\mathcal{I}_{n}}P(x^{\prime},y^{\prime}|x,y)=1.$
(16)

If we introduce the probability density function $f=f(x,t)\geq 0$ of the
agents with label $x$ at time $t$, its evolution can be modelled by a kinetic
equation describing a Markov-type jump process:

$\partial_{t}f(x^{\prime},t)=\lambda\left(\int_{\mathcal{I}_{n}^{3}}P_{xy}^{x^{\prime}y^{\prime}}f(x,t)f(y,t)\,dxdydy^{\prime}-f(x^{\prime},t)\right),$
(17)

where $\lambda>0$ is the (constant) switch frequency. In weak form (17) reads

$\frac{d}{dt}\int_{\mathcal{I}_{n}}\psi(x)f(x,t)\,dx=\lambda\int_{\mathcal{I}_{n}^{2}}\int_{\mathcal{I}_{n}^{2}}(\psi(x^{\prime})-\psi(x))P_{xy}^{x^{\prime}y^{\prime}}f(x,t)f(y,t)\,dx^{\prime}dy^{\prime}\,dxdy,$
(18)

where $\psi:\mathcal{I}_{n}\to\mathbb{R}$ is an observable quantity (test
function) defined on $\mathcal{I}_{n}$. Equation (18) can be derived by (9) by
setting $\boldsymbol{z}=x$ and

$P_{xy}^{x^{\prime}y^{\prime}}=\dfrac{T(x^{\prime}|x,y)+\tilde{T}(y^{\prime}|x,y)}{2}.$

Since $x\in\mathcal{I}_{n}$ is discrete, we may conveniently represent the
distribution function $f$ as

$f(x,t)=\sum_{i=1}^{n}f_{i}(t)\delta(x-i),$ (19)

where $\delta(x-i)$ is the Dirac distribution centred in $x=i$ and
$f_{i}=f_{i}(t)\geq 0$ is the probability that an agent is labelled by $x=i$
at time $t$. In this way, we reconcile the weak form (18) with the convention
introduced in Remark 2.1, and (18) actually becomes

$\sum_{i=1}^{n}\psi(i)f_{i}^{\prime}(t)=\lambda\sum_{i,l=1}^{n}\sum_{j,k=1}^{n}(\psi(i)-\psi(j))P_{jk}^{il}f_{j}(t)f_{k}(t).$
(20)

We have conservation of the total mass thanks to (16), as we can verify
setting $\psi=1$ in (18). Using (19) and, then, setting $\psi=1$ in (20), this
corresponds to

$\sum_{i=1}^{n}f_{i}(t)=1.$ (21)

Choosing $\psi$ such that $\psi(s)=1$ for a certain $s\in\mathcal{I}_{n}$ and
$\psi(x)=0$ for all $x\in\mathcal{I}_{n}\setminus\\{s\\}$ we get in particular

$f_{s}^{\prime}=\lambda\left(\sum_{j,k,l=1}^{n}P_{jk}^{sl}f_{j}f_{k}-f_{s}\right),\qquad
s=1,\,\dots,\,n,$ (22)

where we have used (16) and (21). If we allow the interaction frequency
$\lambda_{xy}$ to depend on the labels of the interacting agents, then the
equation becomes

$\sum_{i=1}^{n}\psi(i)f_{i}^{\prime}(t)=\sum_{i,l=1}^{n}\sum_{j,k=1}^{n}(\psi(i)-\psi(j))\lambda_{jk}P_{jk}^{il}f_{j}(t)f_{k}(t)$

so that

$f_{s}^{\prime}=\sum_{j,k,l=1}^{n}\left(\lambda_{jk}P_{jk}^{sl}f_{j}f_{k}-\lambda_{sk}P_{sk}^{jl}f_{s}f_{k}\right),\qquad
s=1,\,\dots,\,n.$ (23)

In general we define

$\beta_{ij}^{kl}:=\lambda_{ij}P_{ij}^{kl}$ (24)

the rate of transfer for a couple from the subgroups $(i,j)$ to $(k,l)$ as a
consequence of a binary interaction between two agents labelled $i$ and $j$.

###  2.4 Interacting particles with label switch and exchange of physical
quantities

Let us now consider the case in which an agent is characterized by a physical
quantity $v\in\mathbb{R}_{+}$ and by a label $x\in\mathcal{I}_{n}$ that,
again, denotes the belonging of the agent to a certain group. Hence, now, the
microscopic state is
$\boldsymbol{z}=(x,v)\in\Omega=\mathcal{I}_{n}\times\mathbb{R}_{+}$. The
present framework allows to describe a situation in which an agent, as a
consequence of a single binary interaction, changes (simultaneously) both the
microscopic quantity $v$ and the label $x$. Agents within the same group, i.e.
with the same label, are assumed to be indistinguishable. We can take into
account the possibility that the interactions among agents with the same label
differ from those among agents with different labels. In general, if
$(x,v),\,(y,w)\in\mathcal{I}_{n}\times\mathbb{R}_{+}$ denote the pre-
interaction states of any two interacting agents, their post-interaction
quantities $v^{\prime},\,w^{\prime}$ will be given by (10), where
$I=I_{xy},\tilde{I}=\tilde{I}_{xy},D=D_{xy},\tilde{D}=\tilde{D}_{xy}$ may
depend on $x,y$. Moreover, agents can also migrate to other subgroups
$x^{\prime},y^{\prime}$ and this microscopic transfer process is described by
the probability function (15).

We now want to derive a kinetic equation for the joint distribution function
$f=f(x,v,t)\geq 0$, such that $f(x,v,t)dv$ gives the proportion of agents
labelled by $x\in\mathcal{I}_{n}$ and having microscopic state comprised
between $v$ and $v+dv$ at time $t$. The discreteness of $x$ allows us to
represent $f$ as [30]

$f(x,v,t)=\sum_{i=1}^{n}f_{i}(v,t)\delta(x-i),$ (25)

where $f_{i}=f_{i}(v,t)\geq 0$ is the distribution function of the microscopic
state $v$ of the agents with label $i$ and, in particular, $f_{i}(v,t)dv$ is
the proportion of agents with label $i$ whose microscopic state is comprised
between $v$ and $v+dv$ at time $t$.

Since both the interactions and the label switching conserve the total mass of
the system, we may assume that $f(x,v,t)$ is a probability distribution,
namely:

$\int_{\mathbb{R}_{+}}\int_{\mathcal{I}_{n}}f(x,v,t)\,dx\,dv=\sum_{i=1}^{n}\int_{\mathbb{R}_{+}}f_{i}(v,t)\,dv=1\qquad\forall\,t>0.$
(26)

Notice, however, that the $f_{i}$’s are in general not probability density
functions because their $v$-integral varies in time due to the label
switching. We denote by

$\rho_{i}(t):=\int_{\mathbb{R}_{+}}f_{i}(v,t)\,dv$ (27)

the mass of the group of agents with label $i$, thus $0\leq\rho_{i}(t)\leq 1$
and

$\sum_{i=1}^{n}\rho_{i}(t)=1\qquad\forall\,t>0.$

Let us also define the first statistical moment of $f_{i}$

$M_{i}(t):=\int_{\mathbb{R}_{+}}v\,f_{i}(v,t)\,dv$

so that the average of the $i-$th group is

$m_{i}(t):=\dfrac{M_{i}(t)}{\rho_{i}(t)}.$

The kinetic evolution equation for $f(x,v,t)$, expressed as in (25), is given
by (9), where now $\boldsymbol{z}=(x,v)$, which has to hold for every
$\phi=\phi(x,v):\mathcal{I}_{n}\times\mathbb{R}_{+}\to\mathbb{R}$. Hence the
evolution equation for $f$ is

$\displaystyle\begin{aligned}
\frac{d}{dt}\sum_{i=1}^{n}&\int_{\mathbb{R}_{+}}\phi(i,v)f_{i}(v,t)\,dv\\\
&=\dfrac{1}{2}\int_{\mathbb{R}_{+}^{2}}\sum_{i,j,k=1}^{n}\lambda_{jk}\int_{\mathbb{R}_{+}}T((i,v^{\prime})|(j,v),(k,w))\big{(}\phi(i,v^{\prime})-\phi(j,v)\big{)}f_{j}(v,t)f_{k}(w,t)\,dvdwdv^{\prime}\\\
&+\dfrac{1}{2}\int_{\mathbb{R}_{+}^{2}}\sum_{l,j,k=1}^{n}\lambda_{jk}\int_{\mathbb{R}_{+}}\tilde{T}((l,w^{\prime})|(j,v),(k,w))\big{(}\phi(l,w^{\prime})-\phi(k,w)\big{)}f_{j}(v,t)f_{k}(w,t)\,dvdwdw^{\prime}.\end{aligned}$
(28)

Choosing $\phi(x,v)=\psi(x)\varphi(v)$ with $\psi$ such that $\psi(s)=1$ for a
certain $s\in\mathcal{I}_{n}$ and $\psi(x)=0$ for all
$x\in\mathcal{I}_{n}\setminus\\{s\\}$, we finally obtain the following system
of equations for the subgroup distributions $f_{s}$

$\displaystyle\begin{aligned}
\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)f_{s}(v,t)\,dv&=\\\
&=\dfrac{1}{2}\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,i=1}^{n}\int_{\mathbb{R}_{+}}\Big{(}\lambda_{jk}\varphi(v^{\prime})T((s,v^{\prime})|(j,v),(k,w))f_{j}(v,t)\\\
&-\lambda_{sk}\varphi(v)T((i,v^{\prime})|(s,v),(k,w))f_{s}(v,t)\Big{)}f_{k}(w,t)\,dvdwdv^{\prime}\\\
&+\dfrac{1}{2}\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\int_{\mathbb{R}_{+}}\Big{(}\lambda_{jk}\varphi(w^{\prime})\tilde{T}((s,w^{\prime})|(j,v),(k,w))f_{k}(w,t)\\\
&-\lambda_{js}\varphi(w)\tilde{T}((l,w^{\prime})|(j,v),(s,w))f_{s}(w,t)\Big{)}f_{j}(v,t)dvdwdw^{\prime}.\end{aligned}$
(29)

In particular, in order to implement the microscopic process (10)-(15), we
choose

$\displaystyle T((x^{\prime},v^{\prime})|(x,v),(y,w))$
$\displaystyle=\langle\int_{\mathcal{I}_{n}}P_{xy}^{x^{\prime}y^{\prime}}(v,w)\delta(v^{\prime}-(I_{xy}(v,w)+D_{xy}(v,w)\eta))dy^{\prime}\rangle,$
(30) $\displaystyle\tilde{T}((y^{\prime},w^{\prime})|(x,v),(y,w))$
$\displaystyle=\langle\int_{\mathcal{I}_{n}}P_{xy}^{x^{\prime}y^{\prime}}(v,w)\delta(w^{\prime}-(\tilde{I}_{xy}(v,w)+\tilde{D}_{xy}(v,w)\eta_{\ast}))dx^{\prime}\rangle,$

where we remark that
$P_{xy}^{x^{\prime}y^{\prime}}=P_{xy}^{x^{\prime}y^{\prime}}(v,w)$ may depend
on the microscopic physical quantities of the interacting agents. Considering
(30), (29) becomes

$\displaystyle\begin{aligned}
\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)f_{s}(v,t)\,dv&=\\\
&=\dfrac{1}{2}\langle\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\Big{(}\beta_{jk}^{sl}\varphi(I_{xy}(v,w)+D_{xy}(v,w)\eta)f_{j}(v,t)\Big{.}\\\
\Big{.}&\phantom{\dfrac{1}{2}\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\beta_{jk}^{sl}\varphi(I_{xy}(v,w)}-\beta_{sk}^{jl}\varphi(v)f_{s}(v,t)\Big{)}f_{k}(w,t)\,dvdw\rangle\\\
&+\dfrac{1}{2}\langle\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,i=1}^{n}\Big{(}\beta_{jk}^{is}\varphi(\tilde{I}_{xy}(v,w)+\tilde{D}_{xy}(v,w)\eta_{\ast})f_{k}(w,t)\Big{.}\\\
\Big{.}&\phantom{\dfrac{1}{2}\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\beta_{jk}^{sl}\varphi(I_{xy}(v,w)}-\beta_{js}^{ik}\varphi(w)f_{s}(w,t)\Big{)}f_{j}(v,t)dvdw\rangle,\end{aligned}$
(31)

where we have used (24). We remark that now $\beta_{ij}^{kl}$ may depend on
the microscopic variables $v$ and $w$ through $P_{ij}^{kl}$. Moreover, as done
in [28] in a linear case, also the interaction frequency may depend on the
microscopic state $v$.

We remark that in this case where we consider both exchanges (10) (with
$I,\tilde{I},D$ and $\tilde{D}$ that may depend on the labels) and transfers
(15), asymmetric binary interactions arise quite commonly, even if the two
processes (10) and (15), separately, are symmetric. This is mainly due to the
fact that the microscopic rule (10) depends on the label of the agents.
Indeed, if we consider a transfer $(i,j)\rightarrow(k,l)$, the reverse
transfer $(k,l)\rightarrow(i,j)$ may occur with a different probability and
with a different interaction law, losing thus the reversibility of the
process, usually assumed in classical Boltzmann descriptions. Specifically, in
gas mixtures the reversibility is guaranteed by conservations of momentum and
total energy, and the post–collision velocities may be uniquely determined in
terms of the pre–collision velocities and of the impact angles [12, 27], even
in presence of chemical reactions [36]. The break of symmetry between the
direct and the reverse collision is known to occur in presence of inelastic
collisions (for instance in granular media [6]) causing a decay in time of the
kinetic energy of the system. For interactions involving human beings the
kinetic approach is much more complicated (even under simplistic assumptions),
and exchanges of goods and transfers among different compartments may be
non–symmetric. Just to give an example, in socio-economic problems the
fraction of the own wealth that each agent is willing to give to the others
may depend also on the proper amount of wealth [3], and moreover a transfer
from a poor country to a rich one might be much more probable than the reverse
transfer [4]. This is why the general approach with generally different
transition probabilities $T$ and $\tilde{T}$ provides a useful tool for a
correct description of this kind of processes. Moreover, it allows to build
more easily exchange and transfer operators for a generic number $n$ of
subpopulations. Indeed, the usual way of extending Boltzmann theory to a set
of $n>1$ constituents consists in building up a set of $n$ Boltzmann
equations, each one for the distribution function of the $i$-th constituent,
with $i=1,\dots,n$ [24, 36, 21]. On the other hand, our transition probability
approach includes the label of the individual compartment into the set of
microscopic states characterizing the individual, dealing thus with only one
kinetic equation, that could be separated into $n$ different equations only
when one needs to compute the pertinent moments of each compartment by
choosing the appropriate test function as done in (29). This turns out to be a
great advantage even from the computational point of view. As we will see in
numerical tests shown in Section 5, in the present approach it is also
straightforward to consider transition rates $\beta_{ij}^{kl}$ explicitly
dependent on the microscopic states, both through the binary interaction
frequency and /or through the transfer probability $P_{ij}^{kl}$,
investigating thus more realistic cases with respect to classical kinetic
descriptions that, for the sake of simplicity, assume constant interaction
probabilities in the kernel of the Boltzmann operators. Furthermore, this
approach allows to include the stochastic contributions $\eta$ and
$\eta_{\ast}$ in the binary interaction rules, as invertibility is not
required as in the construction of the operators.

##  3 Formal study of the kinetic equation with transition probabilities

In this section we intend to revise and illustrate some analytical tools that
are useful for formally studying equation (9). After briefly stating some
results on the existence and uniqueness of the solution, we consider the
quasi-invariant limit in various regimes of equation (31) involving label
switching and exchange of physical quantities.

###  3.1 Basic theory of kinetic models with transition probabilities in
Wasserstein spaces

The strong form of (9) coupled with an initial condition
$f_{0}(\boldsymbol{z})$ defines the following Cauchy problem

$\begin{cases}&\dfrac{\partial}{\partial
t}f(\boldsymbol{z},t)=Q^{+}(f,f)-f(\boldsymbol{z},t)\displaystyle\int_{\Omega}\lambda_{\boldsymbol{z}\boldsymbol{y}}f(\boldsymbol{y},t)\,d\boldsymbol{y},\qquad
t>0,\quad\boldsymbol{z}\in\Omega,\\\
&f(0,\boldsymbol{z})=f_{0}(\boldsymbol{z}),\qquad\boldsymbol{z}\in\Omega,\end{cases}$
(32)

where

$Q^{+}(f,f)=\dfrac{1}{2}\int_{\Omega^{2}}\lambda_{{}^{\prime}\\!\boldsymbol{z}{}^{\prime}\\!\boldsymbol{y}}T(\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y})f({}^{\prime}\\!\boldsymbol{z},t)f({}^{\prime}\\!\boldsymbol{y},t)\,d{}^{\prime}\\!\boldsymbol{z}d{}^{\prime}\\!\boldsymbol{y}+\dfrac{1}{2}\int_{\Omega^{2}}\lambda_{{}^{\prime}\\!\boldsymbol{z}{}^{\prime}\\!\boldsymbol{y}}\tilde{T}(\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{y},{}^{\prime}\\!\boldsymbol{z})f({}^{\prime}\\!\boldsymbol{z},t)f({}^{\prime}\\!\boldsymbol{y},t)\,d{}^{\prime}\\!\boldsymbol{z}d{}^{\prime}\\!\boldsymbol{y},$
(33)

where ${}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y}$ are the pre-
interaction states, with the compatibility condition
$\int_{\Omega}f_{0}(\boldsymbol{z})\,d\boldsymbol{z}=1$ as (4) holds true. We
remark that everything could be written for a generic mass $\rho>0$. Let us
now define

$\bar{\lambda}:=\int_{\Omega}\lambda_{\boldsymbol{z},\boldsymbol{y}}f(\boldsymbol{y},t)\,d\boldsymbol{y}$

that we assume to be constant throughout the whole text (this assumption
includes the case of a constant $\lambda_{\boldsymbol{z}\boldsymbol{y}}$). If
we multiply both sides of the equation by $e^{\bar{\lambda}t}$ and we
integrate in time we get

$\displaystyle
f(\boldsymbol{z},t)=e^{-\bar{\lambda}t}f_{0}(\boldsymbol{z})+\int_{0}^{t}e^{\bar{\lambda}(s-t)}\left[\dfrac{1}{2}\int_{\Omega^{2}}\right.$
$\displaystyle\left.\lambda_{{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y}}T(\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y})f({}^{\prime}\\!\boldsymbol{z},s)f({}^{\prime}\\!\boldsymbol{y},s)\,d{}^{\prime}\\!\boldsymbol{z}d{}^{\prime}\\!\boldsymbol{y}\right.$
(34)
$\displaystyle\left.+\dfrac{1}{2}\int_{\Omega^{2}}\lambda_{{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y}}\tilde{T}(\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{y},{}^{\prime}\\!\boldsymbol{z})f({}^{\prime}\\!\boldsymbol{z},s)f({}^{\prime}\\!\boldsymbol{y},s)\,d{}^{\prime}\\!\boldsymbol{z}d{}^{\prime}\\!\boldsymbol{y}\right]ds,$

where we have used (4). Let us define $(\Omega,d)$ a polish space. Analogously
to what has been done in [23], we see that an appropriate space in which (34)
can be studied is $X:=\mathcal{C}([0,\bar{t}];\mathcal{M}_{+}(\Omega))$, where
$\bar{t}>0$ is a final time and $\mathcal{M}_{+}(\Omega)$ is the space of
positive measures on $\Omega$ having unitary mass. Therefore $f\in X$ is a
continuous mapping as a function of time over $[0,\bar{t}]$ and it is a
positive measure satisfying (4) as a function of the microscopic state
$\boldsymbol{z}\in\Omega$. In particular, $X$ is a complete state with the
distance

$\sup_{t\in[0,\bar{t}]}W_{1}(f(t,\boldsymbol{z}),g(t,\boldsymbol{y}))$

where

$W_{1}(f(t,\boldsymbol{z}),g(t,\boldsymbol{y}))=\inf_{\underline{\mu}\in\Gamma(f,g)}\int_{\Omega^{2}}d(\boldsymbol{z},\boldsymbol{y})\underline{\mu}(t,\boldsymbol{z},\boldsymbol{y})\,d\boldsymbol{z}d\boldsymbol{y}$
(35)

is the 1-Wasserstein distance between $f(t,\cdot)$ and
$g(t,\cdot)\in\mathcal{M}_{+}(\Omega)$, being $\Gamma(f,g)$ the space of the
probability density functions defined on $\Omega^{2}$ having marginals $f$ and
$g$.

As done in [23], we shall always assume that the transition probabilities $T$
and $\tilde{T}$ satisfy the following Lipschitz continuity property.

######  Assumption 3.1.

Let
$T(\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y}),\tilde{T}(\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y})\in\mathscr{P}(\Omega)$
for all ${}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y}\in\Omega$,
where $\mathscr{P}(\Omega)$ is the space of probability measures on $\Omega$.
We assume that there exists $Lip(T)>0$, such that

$W_{1}(T(\cdot|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y}),T(\cdot|{}^{\prime}\\!\boldsymbol{z}_{*},{}^{\prime}\\!\boldsymbol{y}_{*}))\leq\textrm{Lip}(T)(|{}^{\prime}\\!\boldsymbol{z}-{}^{\prime}\\!\boldsymbol{z}_{*}|+|{}^{\prime}\\!\boldsymbol{y}-{}^{\prime}\\!\boldsymbol{y}_{*}|)$

for all
${}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y},{}^{\prime}\\!\boldsymbol{z}_{*},{}^{\prime}\\!\boldsymbol{y}_{*}\in\Omega$
and that the same holds for $\tilde{T}$.

The following result holds.

######  Theorem 3.2.

Let $f_{0}\in\mathcal{M}_{+}(\Omega)$ and let us assume that $T$, $\tilde{T}$
satisfy assumption 3.1, that $W_{1}(f_{0},T),W_{1}(f_{0},\tilde{T})<\infty$,
that $\bar{\lambda}$ is constant and that
$\lambda_{\boldsymbol{z}\boldsymbol{y}}$ is lower and upper bounded, i.e.
$\exists\,\tilde{\lambda}^{m},\tilde{\lambda}^{M}>0$ such that
$0<\tilde{\lambda}^{m}<\lambda_{\boldsymbol{z}\boldsymbol{y}}<\tilde{\lambda}^{M}<\infty$.
Then, there exists a unique $f\in X$ which solves (32), and (32) exhibits
continuous dependence on the initial data.

Moreover, if
$\tilde{\lambda}^{M}\textrm{Lip}(T),\tilde{\lambda}^{M}\textrm{Lip}(\tilde{T})<\dfrac{1}{2}$,
then (32) admits a unique equilibrium distribution $f_{\infty}$, which is a
probability measure on $\Omega$ and which is also globally attractive, i.e.

$\lim_{t\rightarrow\infty}W_{1}(f(t;\cdot);f_{\infty})=0$

for every solution $f$ to (32).

###### Proof.

The proof follows the same steps as done in [23] Appendix A, where the authors
prove the results in the case of a bounded $\Omega$ and, thus, use the dual
form of the 1-Wasserstein distance due to the Rubinstein-Kantorovitch Theorem
[1]. In the present case, as $\Omega$ is arbitrary, we can use in the proof
the definition (35) recalling the hypothesis
$W_{1}(f_{0},T),W_{1}(f_{0},\tilde{T})<\infty$ and
$\tilde{\lambda}^{m}<\lambda_{\boldsymbol{z}\boldsymbol{y}}<\tilde{\lambda}^{M}$.
∎

Moreover, the following Theorem holds in the case of label switching and
exchange of physical quantities.

######  Theorem 3.3.

Let the transition probability distributions have the form

$T(x,v|({}^{\prime}\\!x,{}^{\prime}\\!v),({}^{\prime}\\!y,{}^{\prime}\\!w))=\sum_{i=1}^{n}T_{i}(v|({}^{\prime}\\!x,{}^{\prime}\\!v),({}^{\prime}\\!y,{}^{\prime}\\!w))\delta(x-i)$

where $T_{i}$ satisfies

$|T_{i}(v|({}^{\prime}\\!x,{}^{\prime}\\!v),({}^{\prime}\\!y,{}^{\prime}\\!w))|<Lip(T_{i})(|{}^{\prime}\\!x-{}^{\prime}\\!y|+|{}^{\prime}\\!v-{}^{\prime}\\!w|)$

and let the analogous property hold for $\tilde{T}$ and $\tilde{T}_{i}$. Let
moreover

$f_{0}(x,v)=\sum_{i=1}^{n}f^{(0)}_{i}(v)\,\delta(x-i)$

be a prescribed kinetic distribution function at time $t=0$ over the space of
microscopic states $(x,v)\in\mathcal{I}_{n}\times\mathbb{R}_{+}$ such that
$f_{i}\geq 0$, $\sum_{i=1}^{n}\int_{\mathbb{R}_{+}}f_{i}(v,t)\,dv=1,\,\forall
t>0$. Then the unique solution to (32) is of the form

$f(x,v,t)=\sum_{i=1}^{n}f_{i}(v,t)\,\delta(x-i)$

(analogous to (19)), with coefficients $f_{i}(v,t)$ given by (23) along with
the initial conditions $f_{i}(v,0)=f^{(0)}_{i}(v)$. In addition, it depends
continuously on the initial datum as stated by Theorem 3.2.

###  3.2 Quasi-invariant limit

One of the most interesting issues in the study of kinetic models is the
characterisation of the stationary distributions arising asymptotically for
$t\rightarrow\infty$, which depict the emergent behaviour of the system. The
jump process model (9) hardly allows one to investigate in detail the trend to
equilibrium and the profile of the stationary distributions, and the explicit
expression of the steady state $f_{\infty}$ can be inferred only in particular
cases [14, 35, 23].

It is widely known that the classic collisional Boltzmann equation (11)-(10)
offers several analytical tools, which often permit to explicitly recover
accurate approximations of $f_{\infty}$ by means of suitable asymptotic
procedures. The basic idea of such procedures is to approximate an integro-
differential Boltzmann equation with an appropriate partial differential
equation, more amenable to analytical investigations, at least in some regimes
of the parameters of the microscopic interactions. A prominent framework in
which this type of asymptotic analysis is successfully applied to (11) is that
of the quasi-invariant interactions. This concept was first introduced in the
kinetic literature on multi-agent systems in [13, 38] as a reminiscence of the
grazing collisions studied in the classical kinetic theory, see [41]. This
corresponds to introducing a small parameter $\epsilon$ such that the
microscopic interaction rule can be written as

$v^{\prime}\approx v+\mathcal{O}(\epsilon)$ (36)

and analyzing the dynamics on a longer time scale, setting a new time variable

$\tau:=\epsilon t$ (37)

in order to compensate for the smallness of the interactions.

In this spirit, in [29], that concerns the investigation of a parallelism
between the model (11)-(10) and (12), the authors propose a way to translate
the concept of quasi-invariancy, typically used in the context of collision-
like Boltzmann equations (11)-(10), to the language of transition
probabilities. The idea is the following. Let
${}^{\prime}\\!\boldsymbol{Z},\,\boldsymbol{Z}\in\mathbb{R}_{+}$ be the random
variables representing the pre- and post-interaction states, respectively, of
an agent, and ${}^{\prime}\\!\boldsymbol{Z}_{\ast}\in\mathbb{R}_{+}$ the one
representing the pre-interaction state of the other agent involved in the
interaction. In the probabilistic description via the transition
probabilities, we say that interactions are quasi-invariant if, given
$0<\epsilon\ll 1$,

$\operatorname{Prob}(\left\lvert\boldsymbol{Z}-{}^{\prime}\\!\boldsymbol{Z}\right\rvert>\epsilon\,|\,{}^{\prime}\\!\boldsymbol{Z},\,{}^{\prime}\\!\boldsymbol{Z}_{\ast})\leq\epsilon;$
(38)

in other words, if the post-interaction state is, in probability, close to the
pre-interaction state, so that the interactions produce a small transfer of
microscopic state between the interacting agents.

In the present framework, in order to have a quasi-invariant transition
probability, we can introduce rescaled transition probabilities defined by the
following transform

$T_{\epsilon}(\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y})=\mathcal{F}_{\epsilon}[T](\boldsymbol{z}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y}),\qquad\tilde{T}_{\epsilon}(\boldsymbol{y}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y})=\mathcal{F}_{\epsilon}[\tilde{T}](\boldsymbol{y}|{}^{\prime}\\!\boldsymbol{z},{}^{\prime}\\!\boldsymbol{y})$
(39)

where

$\mathcal{F}_{\epsilon}[T]:\mathscr{P}_{1}(\mathbb{R}_{+})\longmapsto\mathscr{P}_{1}(\mathbb{R}_{+})$

is a family of operators (for $\epsilon>0$) defined on the space of the
probability measures defined on $\mathbb{R}_{+}$. We require that
$\mathcal{F}_{\epsilon}$ satisfies the following three properties:

  * $F1$

$\mathcal{F}_{1}$ is the identity;

  * $F2$

$\lim_{\epsilon\rightarrow
0}W_{1}(\mathcal{F}_{\epsilon}[T],\delta(\boldsymbol{z}^{\prime}-\boldsymbol{z}))=0,\qquad\lim_{\epsilon\rightarrow
0}W_{1}(\mathcal{F}_{\epsilon}[\tilde{T}],\delta(\boldsymbol{z}^{\prime}-\boldsymbol{z}))=0$
(40)

  * $F3$

$\lim_{\epsilon\rightarrow
0}W_{1}(\dfrac{\mathcal{F}_{\epsilon}[T]}{\epsilon},T)=0,\qquad\lim_{\epsilon\rightarrow
0}W_{1}(\dfrac{\mathcal{F}_{\epsilon}[\tilde{T}]}{\epsilon},\tilde{T})=0$ (41)

meaning that $\epsilon=1$ corresponds to the basic regime ($F1$), that for
small values of $\epsilon$ the microscopic state tends not to change ($F2$),
and that on the long time scale (37) the dynamics is ruled by $T$ and
$\tilde{T}$ ($F3$) [29]. An example of properly rescaled transition
probabilities is

$T_{\epsilon}(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y})=(1-\epsilon)\delta(\boldsymbol{z}^{\prime}-\boldsymbol{z})+\epsilon
T(\boldsymbol{z}^{\prime}|\boldsymbol{z},\boldsymbol{y}),\qquad\tilde{T}_{\epsilon}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})=(1-\epsilon)\delta(\boldsymbol{y}^{\prime}-\boldsymbol{y})+\epsilon\tilde{T}(\boldsymbol{y}^{\prime}|\boldsymbol{z},\boldsymbol{y})$
(42)

as introduced in [29], satisfying the properties $F1$, $F2$, $F3$.

In the following we shall investigate the quasi invariant limit in the three
examples illustrated in the previous section.

####  3.2.1 Boltzmann-type description of classical binary interaction
dynamics

The quasi-invariant limit procedure is classically applied to the collisional
Boltzmann equation with microscopic interaction rules (11)-(10) (or
(11)-(13)). Let us introduce a small parameter $0<\epsilon\ll 1$, a time scale
(37) and a corresponding probability density function
$f^{\epsilon}(\tau,v)=f(\tau/\epsilon,v)$. For what we have said in Section 2,
as shown in [29], we have that the quasi-invariant microscopic rule for having
the same evolution of the average and the energy of both $f$ and
$f^{\epsilon}$ on the $t$\- and $\tau$-time scale, respectively, is

$\displaystyle v^{\prime}=V_{T_{\epsilon}}(v,w)+D_{T_{\epsilon}}(v,w)\eta,$

and analogously for the microscopic rule for $w^{\prime}$. If we consider the
quasi-invariant transition probability (42), we have that
$V_{T_{\epsilon}}(v,w)=v+\epsilon(I(v,w)-v)$ and
$D_{T_{\epsilon}}(v,w)=\sqrt{\epsilon(1-\epsilon)(I(v,w)-v)^{2}+\epsilon
D(v,w)^{2}}$ so that the quasi-invariant microscopic rules are

$\displaystyle
v^{\prime}=v+\epsilon(I(v,w)-v)+\sqrt{\epsilon(1-\epsilon)(I(v,w)-v)^{2}+\epsilon
D(v,w)^{2}}\,\eta),$ (43) $\displaystyle
w^{\prime}=w+\epsilon(\tilde{I}(v,w)-w)+\sqrt{\epsilon(1-\epsilon)(\tilde{I}(v,w)-w)^{2}+\epsilon\tilde{D}(v,w)^{2}}\,\eta_{\ast}).$

Therefore, in terms of transition probabilities, in order to recover the same
evolution of the first two moments in the two models (11)-(10) and (12), as
shown in the previous section, we must consider $T_{\epsilon}$ having average
$\epsilon(I(v,w)-v)$ and variance $\epsilon(1-\epsilon)(I(v,w)-v)^{2}+\epsilon
D(v,w)^{2}$ and analogously for $\tilde{T}_{\epsilon}$. In particular, in
order for (11)-(10) and (12) to be the same model an appropriate choice is
(42) with $\boldsymbol{z}=v$ or

$\displaystyle
T_{\epsilon}(v^{\prime}|v,w)=\delta\Big{(}v^{\prime}-(v+\epsilon(I(v,w)-v)+\sqrt{\epsilon(1-\epsilon)(I(v,w)-v)^{2}+\epsilon
D(v,w)^{2}}\,\eta)\Big{)},$ (44)
$\displaystyle\tilde{T}_{\epsilon}(w^{\prime}|v,w)=\delta\Big{(}w^{\prime}-(w+\epsilon(\tilde{I}(v,w)-w)+\sqrt{\epsilon(1-\epsilon)(\tilde{I}(v,w)-w)^{2}+\epsilon\tilde{D}(v,w)^{2}}\,\eta_{\ast})\Big{)}.$

We remark that for $\epsilon=1$ we recover $T$ and $\tilde{T}$, while for
$\epsilon\rightarrow 0$ we have that

$\displaystyle W_{1}(T_{\epsilon}(\cdot|v,w),\delta(\cdot-v))$
$\displaystyle\leq\int_{\mathbb{R_{+}^{2}}}|v^{\prime}-w^{\prime}|T_{\epsilon}(v^{\prime}|v,w)\delta(w^{\prime}-v)\,dv^{\prime}dw^{\prime}$
$\displaystyle=|\epsilon(I(v,w)-v)+\sqrt{\epsilon(1-\epsilon)(I-v)^{2}+\epsilon
D^{2}}\,\eta|\longrightarrow_{\epsilon\rightarrow 0}0$

and it can also be easily verified that $F3$ holds true for $T_{\epsilon}$
defined by (44), and analogously for $\tilde{T}_{\epsilon}$. If we consider
for a moment the symmetric case (3), for simplicity, we have that plugging
(37) into (12) and considering (39) satisfying (41), letting $\epsilon$ to
$0^{+}$ yields

$\partial_{\tau}f=\int_{\mathbb{R}_{+}}\int_{\mathbb{R}_{+}}T(v\,|\,{}^{\prime}\\!v,\,{}^{\prime}\\!w)f(\tau,\,{}^{\prime}\\!v)f(\tau,\,{}^{\prime}\\!w)\,d{}^{\prime}\\!v\,d{}^{\prime}\\!w-f,$
(45)

which is structurally identical to the very general equation (12) and does not
give any further information. Therefore, in spite of the quasi-invariant
structure of the interactions, it is in principle not easier to extract from
(45) any more detailed information about the asymptotic trends.

Instead, if we consider (12) with (44) and (37), then we can perform the
quasi-invariant limit, as in [29], and we obtain letting $\epsilon$ to $0^{+}$

$\displaystyle\begin{aligned}
\partial_{\tau}f&=\frac{1}{2}\partial^{2}_{v}\left\\{\left[\int_{\mathbb{R}_{+}}\left(\left(V_{T}(v,\,w)-v\right)^{2}+D_{T}^{2}(v,\,w)\right)f(\tau,\,w)\,dw\right]f\right\\}\\\
&\phantom{=}-\partial_{v}\left[\left(\int_{\mathbb{R}_{+}}V_{T}(v,\,w)f(\tau,\,w)\,dw-v\right)f\right],\end{aligned}$
(46)

where $V_{T},D_{T}$ are the average and variance of $T$ defined in Section
2.2.

####  3.2.2 Label switch process caused by binary interactions

Let us now consider the transfer process described by the kinetic equation
(23). In order to write a quasi-invariant regime, we can express the fact
that, given a collision, individuals have a small probability of jumping, i.e.

$P_{\epsilon_{xy}}^{x^{\prime}y^{\prime}}=\epsilon
P_{xy}^{x^{\prime}y^{\prime}}\quad\text{if}\quad(x,y)\neq(x^{\prime},y^{\prime}),\quad
P_{\epsilon_{xy}}^{x^{\prime}y^{\prime}}=1-\epsilon\quad\text{if}\quad(x,y)=(x^{\prime},y^{\prime}).$
(47)

Then, considering a long time scale (37) we have that (23) is

$\dfrac{df_{s}^{\epsilon}}{d\tau}=\dfrac{1}{\epsilon}\sum_{j,k,l=1,\\\ j\neq
s,k\neq
l}^{n}\left(P_{\epsilon_{jk}}^{sl}\lambda_{{jk}}f_{j}^{\epsilon}-P_{\epsilon_{sk}}^{jl}\lambda_{sk}f_{s}^{\epsilon}\right)f_{k}^{\epsilon},\qquad
s=1,\,\dots,\,n.$ (48)

Plugging (47) in (48) we obtain

$\dfrac{df_{s}^{\epsilon}}{d\tau}=\sum_{j,k,l=1}^{n}\left(\beta_{jk}^{sl}f_{j}^{\epsilon}f_{k}^{\epsilon}-\beta_{sk}^{jl}f_{s}^{\epsilon}f_{k}^{\epsilon}\right),\qquad
s=1,\,\dots,\,n$

which is structurally identical to the very general equation (23) with $\tau$
instead of $t$, meaning that it is the probability of jumping that rules the
dynamics on the long time scale.

####  3.2.3 Interacting particles with label switch and exchange of physical
quantities

Let us now consider the case in which the binary interactions lead both to a
transfer and to an exchange of the physical quantity $v$. Therefore, for the
two processes we shall consider a quasi-invariant regime given by (36)-(47).
In this case the transition probability (30) may be rescaled as

$T_{\epsilon}((x^{\prime},v^{\prime})|(x,v),(y,w))=P_{\epsilon_{xy}}^{x^{\prime}y^{\prime}}\delta\Big{(}v^{\prime}-(v+\epsilon(I_{xy}-v)+\sqrt{\epsilon(1-\epsilon)(I_{xy}-v)^{2}+\epsilon
D_{xy}^{2}}\,\eta)\Big{)}$ (49)

and analogously for $\tilde{T}_{\epsilon}$. Let us consider the symmetric case
for simplicity of notation, bearing in mind that the asymmetric case can be
treated analogously. Plugging the latter in (31) and considering the re-
scaling (37) and reminding (47), if we let $\epsilon\rightarrow 0^{+}$ we
obtain

$\displaystyle\begin{aligned}
&\frac{d}{d\tau}\int_{\mathbb{R}_{+}}\varphi(v)f_{s}(\tau,v)\,dv=\\\
&=\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1,j\neq s\vee k\neq
l}^{n}\varphi(v)\big{(}\beta_{jk}^{sl}f_{j}(\tau,v)-\beta_{sk}^{jl}f_{s}(\tau,v)\big{)}f_{k}(\tau,w)\,dvdw\end{aligned}$
(50)

that means that the dynamics is ruled by the label switches and gives no
further information. In case of symmetry, it is the symmetric form of (31).

Let us now consider a different regime and, in particular, let us only
consider a quasi-invariant exchange rule, i.e. (36). Let us then rescale (31)
with (37) and let us consider the quasi-invariant transition probability

$T_{\epsilon}((x^{\prime},v^{\prime})|(x,v),(y,w))=P_{{xy}}^{x^{\prime}y^{\prime}}\delta\Big{(}v^{\prime}-(v+\epsilon(I_{xy}-v)+\sqrt{\epsilon(1-\epsilon)(I_{xy}-v)^{2}+\epsilon
D_{xy}^{2}}\,\eta)\Big{)},$ (51)

i.e. the exchange of the physical quantity is actually quasi-invariant, whilst
the label-switch process is not. We obtain

$\displaystyle\begin{aligned}
&\frac{d}{d\tau}\int_{\mathbb{R}_{+}}\varphi(v)f_{s}(\tau,v)\,dv=\\\
&=\langle\dfrac{1}{\epsilon}\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\big{(}\beta_{jk}^{sl}\varphi(v+\epsilon(I_{jk}-v)+\sqrt{\epsilon(1-\epsilon)(I_{jk}-v)^{2}+\epsilon
D_{jk}^{2}}\eta)f_{j}(\tau,v)\\\
&\phantom{=\dfrac{1}{\epsilon}\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\big{(}\beta_{jk}^{sl}\varphi(v+\epsilon(I_{jk}-v)+}-\beta_{sk}^{jl}\varphi(v)f_{s}(\tau,v)\big{)}f_{k}(\tau,w)\,dvdw\rangle.\end{aligned}$
(52)

As $\epsilon$ is small, we can Taylor expand

$\displaystyle\langle\varphi(v+\epsilon(I_{jk}(v,w)-v)+\sqrt{\epsilon(1-\epsilon)(I_{jk}-v)^{2}+\epsilon
D_{jk}^{2}}\eta)\rangle=$
$\displaystyle=\varphi(v)+\epsilon(I_{jk}(v,w)-v)\varphi^{\prime}(v)+\dfrac{1}{2}\varphi^{\prime\prime}(v)\epsilon\left((1-\epsilon)(I_{jk}-v)^{2}+D_{jk}^{2}\right)+\mathcal{O}(\epsilon^{2})$

and plugging the latter in (52) we have that

$\displaystyle\begin{aligned}
&\frac{d}{d\tau}\int_{\mathbb{R}_{+}}\varphi(v)f_{s}(\tau,v)\,dv=\\\
&=\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\beta_{jk}^{sl}\big{(}\varphi^{\prime}(v)(I_{jk}-v)+\dfrac{1}{2}((1-\epsilon)(I_{jk}-v)^{2}+D_{jk}^{2})\varphi^{\prime\prime}(v)\big{)}f_{j}(\tau,v)f_{k}(\tau,w)\,dvdw\\\
&+\dfrac{1}{\epsilon}\int_{\mathbb{R}_{+}}\sum_{j,k,l=1,j\neq
s}^{n}\varphi(v)\big{(}\beta_{jk}^{sl}f_{j}(\tau,v)-\beta_{sk}^{jl}f_{s}(\tau,v)\big{)}\rho_{k}(\tau)\,dv\end{aligned}$

where we remind that $\rho_{k}$ is the mass of the $k$-th population. Let us
now consider an expansion for the probability density function of the whole
population $f$

$f(\tau,x,v)=f^{(0)}(\tau,x,v)+\epsilon
f^{(1)}(\tau,x,v)+\mathcal{O}(\epsilon^{2})$ (53)

where the zero-th and first order moments satisfy

$\displaystyle\int_{\mathcal{I}_{n}\times\mathbb{R}_{+}}fdxdv=\rho^{(0)}:=\int_{\mathcal{I}_{n}\times\mathbb{R}_{+}}f^{(0)}\,dxdv,\qquad\rho^{(1)}:=\int_{\mathcal{I}_{n}\times\mathbb{R}_{+}}f^{(1)}\,dxdv=0,$
(54)
$\displaystyle\int_{\mathcal{I}_{n}\times\mathbb{R}_{+}}fvdxdv=M^{(0)}:=\int_{\mathcal{I}_{n}\times\mathbb{R}_{+}}f^{(0)}v\,dxdv,\qquad
M^{(1)}:=\int_{\mathcal{I}_{n}\times\mathbb{R}_{+}}f^{(1)}v\,dxdv=0.$

For each $f_{i}$ this translates into

$f_{i}(\tau,v)=f_{i}^{(0)}(\tau,v)+\epsilon
f_{i}^{(1)}(\tau,v)+\mathcal{O}(\epsilon^{2})$

and (54) translates to

$\rho^{(0)}=\sum_{i=1}^{n}\rho_{i}^{(0)}=\int_{\mathcal{I}_{n}\times\mathbb{R}_{+}}fdxdv,\qquad\sum_{i=1}^{n}\rho_{i}^{(1)}=0.$
(55)

Comparing equal orders of $\epsilon$ and supposing symmetry, we obtain

$\displaystyle\dfrac{1}{\epsilon}\int_{\mathbb{R}_{+}}\sum_{j,k,l=1,j\neq
s}^{n}\varphi(v)\big{(}\beta_{jk}^{sl}f_{j}^{(0)}(\tau,v)-\beta_{sk}^{jl}f_{s}^{(0)}(\tau,v)\big{)}\rho_{k}^{(0)}(\tau)dv=0,$
(56)

so that

$\displaystyle f_{s}^{(0)}(\tau,v)=\dfrac{\sum_{j,k,l=1,j\neq
s}^{n}\beta_{jk}^{sl}f_{j}^{(0)}(\tau,v)\rho_{k}^{(0)}(\tau)}{\sum_{j,k,l=1,j\neq
s}^{n}\beta_{sk}^{jl}\rho_{k}^{(0)}(\tau)},$ (57)

while, at the first order

$\displaystyle\frac{d}{d\tau}\int_{\mathbb{R}_{+}}\varphi(v)f_{s}^{(0)}(\tau,v)\,dv=$
(58)
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\sum_{j,k,l=1}^{n}\beta_{jk}^{sl}\big{(}\varphi^{\prime}(v)(I_{jk}(v,w)-v)+\frac{1}{2}((I_{jk}(v,w)-v)^{2}+D_{jk}^{2})\varphi^{\prime\prime}(v)\big{)}f_{j}^{(0)}(\tau,v)f_{k}^{(0)}(\tau,w)\,dvdw$
$\displaystyle+\int_{\mathbb{R}_{+}}\sum_{j,k,l=1,j\neq
s}^{n}\varphi(v)\big{(}\beta_{jk}^{sl}f_{j}^{(0)}(\tau,v)-\beta_{sk}^{jl}f_{s}^{(0)}(\tau,v)\big{)}\rho_{k}^{(1)}(\tau)dv$
$\displaystyle+\int_{\mathbb{R}_{+}}\sum_{j,k,l=1,j\neq
s}^{n}\varphi(v)\big{(}\beta_{jk}^{sl}f_{j}^{(1)}(\tau,v)-\beta_{sk}^{jl}f_{s}^{(1)}(\tau,v)\big{)}\rho_{k}^{(0)}(\tau)dv,$

from which we can obtain, resorting to the strong form thanks to integration
by parts, a Fokker-Planck-type equation with a reaction term for each
$f_{s}^{(0)}$ that is

$\displaystyle\partial_{\tau}f_{s}^{(0)}(\tau,v)$ $\displaystyle=$ (59)
$\displaystyle-\partial_{v}\sum_{j,k,l=1}^{n}\beta_{jk}^{sl}\int_{\mathbb{R}_{+}}(I_{jk}(v,w)-v)f_{k}^{(0)}(\tau,w)\,dwf_{j}^{(0)}(\tau,v)$
$\displaystyle+\partial^{2}_{vv}\dfrac{1}{2}\sum_{j,k,l=1}^{n}\beta_{jk}^{sl}\int_{\mathbb{R}_{+}}((I_{jk}(v,w)-v)^{2}+D_{jk}^{2})f_{k}^{(0)}(\tau,w)\,dwf_{j}^{(0)}(\tau,v)$
$\displaystyle+\sum_{j,k,l=1,j\neq
s}^{n}\big{(}\beta_{jk}^{sl}f_{j}^{(0)}(\tau,v)-\beta_{sk}^{jl}f_{s}^{(0)}(\tau,v)\big{)}\rho_{k}^{(1)}(\tau)$
$\displaystyle+\sum_{j,k,l=1,j\neq
s}^{n}\big{(}\beta_{jk}^{sl}f_{j}^{(1)}(\tau,v)-\beta_{sk}^{jl}f_{s}^{(1)}(\tau,v)\big{)}\rho_{k}^{(0)}(\tau).$

The latter reaction terms also involve the first order corrections
$f_{i}^{(1)},\,i=1,...,n$. In order to find univocally the solutions
$f_{s}^{(0)},f_{s}^{(1)}$ to (57) and (59) satisfying (55), we need a number
of conditions (to be looked for example in conserved quantities) that is equal
to the number of degrees of freedom.

##  4 Kinetic model for international trade allowing transfer of individuals

In this section, we are going to rephrase with the current framework the
kinetic model for international trade allowing transfer of individuals
investigated in [4], where the author presents a model of interacting
individuals divided into two subpopulations and allowed, by means of binary
interactions, to exchange wealth and to migrate to the other subgroup. Here,
then, we have that $n=2$, the physical quantity $v\in\mathbb{R}_{+}$ is the
wealth, while the label $x\in\mathcal{I}_{2}$ denotes the subgroup. Note that,
in the present work, we are only going to consider binary interactions giving
rise to both exchanges of the wealth and transfer simultaneously.

###  4.1 From the microscopic to the macroscopic model

For what concerns the exchange of the physical quantity $v$, we are going to
consider simple linear microscopic rules (10) with

$I_{xy}(v,w)=(1-\omega_{x})v+\omega_{y}w,\qquad D_{xy}=\zeta_{xy}v$ (60)

where $\omega_{i}\in[0,1],\,i\in\mathcal{I}_{2}$ and we are dealing with
symmetric interactions. We consider possible transfers given by

$\displaystyle\textrm{(a)}\qquad 1+1\rightarrow 1+2,\qquad\textrm{(b)}\qquad
2+2\rightarrow 1+2,$ (61) $\displaystyle\textrm{(c)}\qquad 1+2\rightarrow
1+1,\qquad\textrm{(d)}\qquad 1+2\rightarrow 2+2,$

therefore only one of the two interacting agents moves to the other subgroup.
The latter implies that the only non-vanishing values of $P_{ij}^{kl}$
correspond to the 4-plets

$(i,j,k,l)\in\\{(1,1,1,2),(1,1,2,1),(2,2,1,2),(2,2,2,1),(1,2,1,1),(1,2,2,2),(2,1,1,1),(2,1,2,2)\\}.$
(62)

The kinetic equation describing this microscopic dynamics is (29) with (30),
where $P_{ij}^{kl}$ is defined according to (62) and the microscopic exchange
dynamics by (60). The evolution of the mass of each subpopulation is given by
setting $\varphi=1$ in (29), with the prescribed dynamics, for $s=1,2$ and it
results in

$\displaystyle\begin{aligned}
\partial_{t}\rho_{1}=(\beta_{22}^{12}\rho_{2}+\beta_{12}^{11}\rho_{1})\rho_{2}-(\beta_{11}^{12}\rho_{1}+\beta_{12}^{22}\rho_{2})\rho_{1}\\\
\partial_{t}\rho_{2}=(\beta_{11}^{12}\rho_{1}+\beta_{12}^{22}\rho_{2})\rho_{1}-(\beta_{22}^{12}\rho_{2}+\beta_{12}^{11}\rho_{1})\rho_{2}\end{aligned}$
(63)

while setting $\varphi=v$ in (29) gives the evolution of the first moments for
$s=1,2$

$\displaystyle\begin{aligned}
\partial_{t}M_{1}=\left(\beta_{12}^{11}\rho_{1}+\beta_{22}^{12}\rho_{2}\right)M_{2}-\left(\beta_{12}^{22}\rho_{2}+\beta_{11}^{12}\rho_{1}\right)M_{1},\\\
\partial_{t}M_{2}=\left(\beta_{12}^{22}\rho_{2}+\beta_{11}^{12}\rho_{1}\right)M_{1}-\left(\beta_{12}^{11}\rho_{1}+\beta_{22}^{12}\rho_{2}\right)M_{2}.\end{aligned}$
(64)

As a consequence, the averages of the wealth of population 1 and 2 evolve as

$\displaystyle\begin{aligned}
\partial_{t}m_{1}=\dfrac{\rho_{2}}{\rho_{1}}\left(\beta_{22}^{12}\rho_{2}+\beta_{12}^{11}\rho_{1}\right)(m_{2}-m_{1}),\\\
\partial_{t}m_{2}=\dfrac{\rho_{1}}{\rho_{2}}\left(\beta_{11}^{12}\rho_{1}+\beta_{12}^{22}\rho_{2}\right)(m_{1}-m_{2}).\end{aligned}$
(65)

We remark that we have assumed symmetry in the interaction rates, i.e.

$\beta_{ii}^{12}=\beta_{ii}^{21},\quad\beta_{12}^{ii}=\beta_{21}^{ii},\quad\forall
i=1,2.$

We observe that the total mass and average

$\bar{\rho}:=\rho_{1}+\rho_{2}=1,\qquad\bar{M}:=M_{1}+M_{2}$

are conserved in time. Regarding the stationary states of the masses, we have
that

$\rho_{2}^{\infty}=\alpha\rho_{1}^{\infty}$

where

$\alpha=\dfrac{-(\beta_{12}^{11}-\beta_{12}^{22})+\sqrt{(\beta_{12}^{11}-\beta_{12}^{22})^{2}+4\beta_{11}^{12}\beta_{22}^{12}}}{2\beta_{22}^{12}},$

and, taking into account that the sum of the two densities is constant, we
have that

$\rho_{1}^{\infty}=\dfrac{\bar{\rho}}{1+\alpha},\qquad\rho_{2}^{\infty}=\dfrac{\alpha\bar{\rho}}{1+\alpha}.$

Therefore

$\rho_{1}^{\infty}=\rho_{2}^{\infty}=\dfrac{\bar{\rho}}{2}\qquad\textrm{if and
only if}\qquad\alpha=1.$

Bearing in mind that $\beta_{ij}^{kl}=P_{ij}^{kl}\lambda_{ij}$, the latter
condition is satisfied if $P_{12}^{22}=P_{12}^{11}=0.5$ and
$\lambda_{11}=\lambda_{22}$, which means that the probability for interacting
agents with different labels and going to the same subgroup is the same, and
the frequency of interaction among agents of the same subgroup is the same for
all subgroups. To this regard, we observe that, as the stationary state only
depends on $\alpha$, then there may be a switch in the population size
($(\rho_{2}^{\infty}-\rho_{1}^{\infty})(\rho_{2}(0)-\rho_{1}(0))<0$) if

$(\rho_{2}(0)-\rho_{1}(0))(\alpha-1)<0.$ (66)

For what concerns the average, the sufficient and necessary condition to be
met at the stationary state is

$m_{1}^{\infty}=m_{2}^{\infty}=:m^{\infty}$ (67)

for every choice of the parameters. The latter implies that
$M_{2}^{\infty}>M_{1}^{\infty}$ if and only if $\alpha>1$. Moreover, because
of conservation of mass and total momentum, we have that

$m^{\infty}=\rho_{1}(0)m_{1}(0)+\rho_{2}(0)m_{2}(0)$ (68)

that implies that the final average wealth is closer to the initial average
wealth of the subgroup that was more populated at $t=0$.

###  4.2 Quasi-invariant limit

If we assume, for simplicity of notation, that the probability of transfer
towards the $i$-th subgroup is independent of the countries of the interacting
agents, i.e.

$\beta_{1}^{2}:=\beta_{11}^{12}=\beta_{12}^{22},\qquad\beta_{2}^{1}:=\beta_{22}^{12}=\beta_{12}^{11}$
(69)

we have that

$\rho_{2}^{\infty}=\alpha\rho_{1}^{\infty},\qquad M_{2}^{\infty}=\alpha
M_{1}^{\infty},$ (70)

where $\alpha=\beta_{1}^{2}/\beta_{2}^{1}$, and, as $\bar{M}$ and
$\bar{\rho}=1$ are conserved quantities, we have that

$\rho_{1}^{\infty}=\dfrac{\bar{\rho}}{1+\alpha},\qquad
M_{1}^{\infty}=\dfrac{\bar{M}}{1+\alpha}.$

Let us now consider the quasi-invariant regime defined by the transition
probability

$T_{\epsilon}((x^{\prime},v^{\prime})|(x,v),(y,w))=P_{{xy}}^{x^{\prime}y^{\prime}}\delta\Big{(}v^{\prime}-(v+\epsilon(I_{xy}(v,w)-v)+\sqrt{\epsilon
D_{xy}(v,w)^{2}}\eta)\Big{)}.$ (71)

The latter, even if it satisfies the requirements F1, F2, F3 prescribed in
Section 3.2, differently from (51), does not guarantee the same evolution of
the energy in the quasi-invariant regime. By expanding the distribution
functions in powers of $\epsilon$, imposing that the globally invariant
quantities (zero-th and first order moments) remain unexpanded, we get that
the constraints (55) are

$\rho_{1}^{(1)}+\rho_{2}^{(1)}=0,\qquad M_{1}^{(1)}+M_{2}^{(1)}=0.$ (72)

Therefore, we have two degrees of freedom and we can determine the values of
$\rho_{1}^{(1)},M_{1}^{(1)}$ as we have two conserved quantities. Following
the same procedure as before, we find that (57) now is

$f_{1}^{(0)}=\dfrac{1}{\alpha}f_{2}^{(0)}\,.$ (73)

Then we have that $\rho_{2}^{(0)}=\alpha\rho_{1}^{(0)}$ and
$M_{2}^{(0)}=\alpha M_{1}^{(0)}$, i.e.
$\rho_{i}^{(0)}=\rho_{i}^{\infty},M_{i}^{(0)}=M_{i}^{\infty},\,i=1,2$, which
means that the masses and averages of order zero correspond to the equilibrium
ones. Therefore

$\rho_{1}^{(0)}=\dfrac{\bar{\rho}}{1+\alpha},\qquad
M_{1}^{(0)}=\dfrac{\bar{M}}{1+\alpha}\,.$ (74)

At the first order (59) for $s=1$ (for $s=2$ an analogous result applies)
specialises into

$\displaystyle\partial_{\tau}f_{1}^{(0)}$
$\displaystyle=-\partial_{v}\Big{(}\beta_{1}^{2}(\omega_{1}M_{1}^{(0)}-\omega_{2}v\rho_{1}^{(0)})f_{1}^{(0)}+\beta_{1}^{2}(\omega_{2}M_{2}^{(0)}-\omega_{2}v\rho_{2}^{(0)})f_{1}^{(0)}\Big{.}$
(75) .
$\displaystyle\phantom{=-\partial_{v}(\beta_{1}^{2}}+\beta_{1}^{2}(\omega_{1}M_{1}^{(0)}-\omega_{1}v\rho_{1}^{(0)})f_{1}^{(0)}+\beta_{2}^{1}(\omega_{2}M_{2}^{(0)}-\omega_{1}v\rho_{2}^{(0)})f_{1}^{(0)}\Big{)}$
$\displaystyle\phantom{=}+\partial^{2}_{v^{2}}\dfrac{1}{2}\Big{(}\beta_{1}^{2}\left[\zeta_{12}^{2}v^{2}\right]f_{1}^{(0)}\rho_{1}+\beta_{1}^{2}\left[\zeta_{22}^{2}v^{2}\right]f_{1}^{(0)}\rho_{2}\Big{.}$
.
$\displaystyle\phantom{=-\partial_{v}(\beta_{1}^{2}}+\beta_{1}^{2}\left[\zeta_{11}^{2}v^{2}\right]f_{1}^{(0)}\rho_{1}^{(0)}+\beta_{2}^{1}\left[\zeta_{12}^{2}v^{2}\right]f_{1}^{(0)}\rho_{2}^{(0)}\Big{)}$
$\displaystyle\phantom{=}+\left(f_{2}^{(1)}\beta_{2}^{1}-f_{1}^{(1)}\beta_{1}^{2}\right)\bar{\rho}$

where use of (73) has been made. Integrating (75) over $\mathbb{R}_{+}$ and
(75) multiplied by $v$ over $\mathbb{R}_{+}$, along with the conditions
$f_{1}(0)=0$ and $\lim_{v\rightarrow+\infty}f_{1}(v)=0$, and remembering (72),
we discover that both $\rho_{1}^{(1)}=\rho_{2}^{(1)}=0$ and
$M_{1}^{(1)}=M_{2}^{(1)}=0$. This implies the fact that both the masses and
the averages of $f_{1}$ and $f_{2}$ are at the equilibrium even at
$\mathcal{O}(\epsilon)$ accuracy. Using relations (73)-(74) in (75), we obtain

$\displaystyle\partial_{\tau}f_{1}^{(0)}$
$\displaystyle=-\dfrac{\beta_{1}^{2}}{1+\alpha}\partial_{v}\left[\Big{(}(\omega_{1}\bar{M}-\omega_{2}v\bar{\rho})+\alpha(\omega_{2}\bar{M}-\omega_{2}v\bar{\rho})+(\omega_{1}\bar{M}-\omega_{1}v\bar{\rho})+(\omega_{2}\bar{M}-\omega_{1}v\bar{\rho})\Big{)}f_{1}^{(0)}\right]$
(76)
$\displaystyle\phantom{=}+\dfrac{\zeta^{2}\beta_{1}^{2}\bar{\rho}}{2(1+\alpha)}\partial^{2}_{v^{2}}\Big{(}v^{2}(3+\alpha)f_{1}^{(0)}\Big{)}$
$\displaystyle\phantom{=}+\left(f_{2}^{(1)}\beta_{2}^{1}-f_{1}^{(1)}\beta_{1}^{2}\right)\bar{\rho}$

where we have also assumed that the stochastic fluctuations are the same in
each kind of interaction, i.e.

$\zeta_{ij}=\zeta,\qquad\forall i,j\in\mathcal{I}_{2}.$

Equation (76) is a Fokker-Planck equation with reaction term, where the
advection-diffusion part (first and second line) only involves $f_{1}^{(0)}$
and its (known) mass and average (74), while the reaction term (third line in
(76)) only depends on $f_{1}^{(1)},f_{2}^{(1)}$. As
$\rho_{1}^{(1)}=\rho_{2}^{(1)}=0$ and $M_{1}^{(1)}=M_{2}^{(1)}=0$, the
reaction term does not influence the mass and average of $f_{1}^{(0)}$. It is
therefore reasonable to look for the stationary solution to the Fokker-Planck
equation without reaction term, i.e.

$\displaystyle\partial_{\tau}\tilde{f}_{1}^{(0)}$
$\displaystyle=-\dfrac{\beta_{1}^{2}}{1+\alpha}\partial_{v}\left[\Big{(}(\omega_{1}\bar{M}-\omega_{2}v\bar{\rho})+\alpha(\omega_{2}\bar{M}-\omega_{2}v\bar{\rho})+(\omega_{1}\bar{M}-\omega_{1}v\bar{\rho})+(\omega_{2}\bar{M}-\omega_{1}v\bar{\rho})\Big{)}\tilde{f}_{1}^{(0)}\right]$
(77)
$\displaystyle\phantom{=}+\dfrac{\zeta^{2}\beta_{1}^{2}\bar{\rho}}{2(1+\alpha)}\partial^{2}_{v^{2}}\Big{(}v^{2}(3+\alpha)\tilde{f}_{1}^{(0)}\Big{)}$

that is (76) where we neglect the third term on the right hand side, as this
does not contribute to a variation of mass and average of $f_{1}^{(0)}$.
Therefore we obtain

$\tilde{f}_{1}^{(0)}=\dfrac{\bar{\rho}}{1+\alpha}v^{-2\left(1+\dfrac{\gamma}{2}\right)}\exp^{-\dfrac{\bar{M}}{\bar{\rho}}\dfrac{\gamma}{v}},\qquad\gamma=\dfrac{B}{D},\quad
B=2\omega_{1}+\omega_{2}(1+\alpha),\quad D=\zeta^{2}\dfrac{3+\alpha}{2},$ (78)

The mass and average can be verified to be $\bar{\rho}/(1+\alpha)$ and
$\dfrac{\bar{M}}{1+\alpha}$ respectively, while the energy is
$\dfrac{\bar{M}}{(1+\alpha)(\gamma-1)}$. Moreover, we can determine the Pareto
index of the first population that is (approximated by)

$PI_{1}=\gamma+1=\dfrac{2\omega_{1}+\omega_{2}(1+\alpha)}{\frac{\zeta^{2}}{2}(3+\alpha)}+1$
(79)

that depends on the trading propensity of both populations
$\omega_{1},\omega_{2}$, on the ratio $\alpha$ that involves the rates
$\beta_{i}^{j}$ and on the stochasticity $\zeta^{2}$. Since, according to
(73), $\tilde{f}_{2}^{(0)}=\alpha\,\tilde{f}_{1}^{(0)}$, both populations have
the same (approximate) Pareto index.

We remark that considering only one population corresponds to setting
$\omega_{1}=\omega_{2}=\omega$ and $\beta_{1}^{2}=\beta_{2}^{1}$. If
$\beta_{2}^{1}=\beta_{1}^{2}$, then there is no reaction term in (76) so that
the stationary state (78) is exact and $\alpha=1$ that implies
$f_{1}^{(0)}=f_{2}^{(0)}$. Moreover, the Pareto index (now exact) is

$\dfrac{2\omega}{\zeta^{2}}+1,$

that coincides with the one commonly obtained from a kinetic model for a
single population [13]. Note that, even keeping $\alpha\not=1$, in the case
$\omega_{1}=\omega_{2}$ the Pareto index of each group coincides with that
relevant to a single population.

##  5 Numerical tests

In this section we present some numerical tests that illustrate the dynamics
of the model that we have introduced in the previous section. We integrate the
kinetic equation (31) numerically using a modified version of the Nanbu-
Babovski Monte Carlo algorithm (see Algorithm 1 in the Appendix). The latter
is based on a direct implementation of the time discrete stochastic
microscopic process (5)-(6)-(8) with (30) as illustrated in Section 2.4 for
$N$ agents, which in the limit $\Delta t\rightarrow 0^{+}$ produces the
kinetic equation (31). In particular, we shall consider the microscopic rules
(60)-(62), with $n=2$. We perform an empirical statistics of the $N$ simulated
agents and define the distribution functions $f_{1}^{MC},f_{2}^{MC}$, their
masses $\rho_{1}^{MC},\rho_{2}^{MC}$ and first moments
$M_{1}^{MC},M_{2}^{MC}$.

###  5.1 Interacting particles with label switch and exchange of physical
quantities

In all numerical tests we consider $\rho_{1}(0)=0.9$ and $\rho_{2}(0)=0.1$,
i.e. the subgroup labelled with $x=1$ is initially more populated. We remark
that the assumed symmetry in the process (62) implies
$P_{11}^{11}=P_{11}^{22}=P_{22}^{22}=P_{22}^{11}=0$ and then
$P_{11}^{12}=P_{11}^{21}=0.5$, being $P_{ij}^{kl}$ a conditional probability.
In the following numerical tests, we shall always consider
$\omega_{1}=\omega_{2}=0.5$, as the values of $\omega_{1}$ and $\omega_{2}$ do
not affect the averages’ evolution and stationary state (65). Moreover, we fix

$\lambda_{11}=1$

and we vary the other parameters.

In the first set of simulations (Test 1 in the following), we choose
$\lambda_{22}=10,\lambda_{12}=1$, $P_{12}^{11}=0.5,P_{12}^{22}=0.5$ and we
consider two different initial conditions for the distributions $f_{1}$ and
$f_{2}$. In Case A we have that the first population is poorer than the second
population at time $t=0$, i.e.

$f_{1}(v,0)=\rho_{1}(0)\mathbb{1}_{[0,1]}(v),\qquad
f_{2}(v,0)=\rho_{2}(0)\dfrac{1}{10}\mathbb{1}_{[5,15]}(v)$ (80)

while in Case B we invert the initial wealths

$f_{1}(v,0)=\rho_{1}(0)\dfrac{1}{10}\mathbb{1}_{[5,15]}(v),\qquad
f_{2}(v,0)=\rho_{2}(0)\mathbb{1}_{[0,1]}(v),$

i.e. the second population is poorer than the first population at time $t=0$.
We report the results in Figure 1. First of all, we observe that this choice
of parameters prescribes $\alpha<1$, that, as showed by the macroscopic
equations (63), implies $\rho_{1}^{\infty}>\rho_{2}^{\infty}$ (see Figure
1(a)). As forecast by theoretical results, the final average wealth
$m^{\infty}=m_{1}=m_{2}$ is closer to the initial wealth of the initially more
populated subgroup: then $m^{\infty}$ is smaller in case A) and larger in case
B) (see Figure 1 (b)). This implies a different behaviour of the first moment
of both $f_{1}$ and $f_{2}$ (see Figure 1 (c)): while in case B) the first
population remains the richer one as it is the one that is initially more
populated, in scenario A), the mean wealth is inverted as the first population
becomes reacher. In each case we compare the evolution of the macroscopic
quantities
$\rho_{1}^{MC},\rho_{2}^{MC},M_{1}^{MC},M_{2}^{MC},m_{1}^{MC},m_{2}^{MC}$ as
prescribed by the microscopic model (5)-(6)-(8)-(30) with (60)-(62) and the
ones whose evolution is given by the derived equations (63)-(64)-(65) for the
macroscopic quantities $\rho_{1},\rho_{2},M_{1},M_{2},m_{1},m_{2}$. Being
$\Delta t=1e-2$ and $N=10^{6}$, we observe a very good agreement between the
solution of the microscopic model and the one of the macroscopic model. The
integration of the kinetic equations also allows to approximate numerically
the distribution functions $f_{1}$ and $f_{2}$ that we report at the
equilibrium in Figure 1(d) in both cases A) and B). In Figure 1 (e)-(f) we
report the time evolution of $f_{1}^{MC}$ and $f_{2}^{MC}$ in case A) with
initial condition (80). In Fig. 1 (e)-(f) we report the time evolution for
$f_{1}^{MC},f_{2}^{MC}$.

Figure 1: Test 1: $\lambda_{11}=1,\lambda_{22}=10,\lambda_{12}=1$,
$P_{12}^{11}=0.5,P_{12}^{22}=0.5$. In all figures we report the time evolution
of the macroscopic quantities as prescribed by the kinetic model and by the
derived macroscopic equations: masses (a), averages (b), first moments (c).
Straight lines correspond to the solution to the macroscopic equations
(63)-(65)-(64), while circles correspond to the solution of the kinetic
equation that we obtain by simulating with the Monte Carlo algorithm 1 the
microscopic model (5)-(6)-(8)-(30) as illustrated in Sec. 2.4 with the
microscopic rules (60)-(62). In all figures we compare the results of the
solutions given the two different initial conditions A and B. In figure (d) we
show the steady states of distributions $f_{1}$, $f_{2}$ in both test cases A)
and B), while in figures (e), (f) we report time evolution of distributions
functions in the test case A).

In the second set of simulations (Test 2) we choose
$\lambda_{11}=1,\lambda_{22}=1,\lambda_{12}=10$, i.e. intra-species
interactions have the same frequency in the two populations, while the inter-
group interactions have a higher frequency. Moreover, we consider three cases
for the inter-group interactions: in case i)
$P_{12}^{11}=0.5,P_{12}^{22}=0.5$, i.e. given an inter-group interaction, the
probability for both agents of transferring is the same, while in case ii)
$P_{12}^{11}=0.2,P_{12}^{22}=0.8$, i.e. the probability of transferring to the
subgroup 2 is higher and iii) $P_{12}^{11}=0.8,P_{12}^{22}=0.2$, i.e. the
probability of transferring to the subgroup 1 is higher. The initial condition
is set as in (80). We observe that in the three cases i), ii) and iii) we have
respectively $\alpha=1,\alpha>1,\alpha<1$, that imply, see (63),
$\rho_{1}^{\infty}=\rho_{2}^{\infty},\rho_{1}^{\infty}<\rho_{2}^{\infty},\rho_{2}^{\infty}<\rho_{1}^{\infty}$,
respectively. In particular, because of (66), in case ii) we have a switch in
the trend of the populations, as population 1 becomes the less populated (see
Fig. 2 (a)). This also implies a different (non monotone) trend of the first
moments as reported in Fig. 2 (b).

Figure 2: Test 2. Solutions to the macroscopic model (63)-(64). In (a) we
report the masses $\rho_{1},\rho_{2}$ and in (b) the first moments
$M_{1},M_{2}$. Population 1 is in red,while population 2 is in blue. The
parameters are $\lambda_{11}=1,\lambda_{22}=1,\lambda_{12}=10$ and we show
three cases: i) (straight lines) $P_{12}^{11}=0.5,P_{12}^{22}=0.5$, ii) (*
marker) $P_{12}^{11}=0.2,P_{12}^{22}=0.8$ and iii) (o marker)
$P_{12}^{11}=0.8,P_{12}^{22}=0.2$.

In Test 3, we consider a switching probability depending on the microscopic
wealth. In this case, it is not immediate, in general, to derive equations for
the macroscopic quantities, unless in special cases, for example $P_{ij}^{kl}$
having a linear dependence on $v$, or by imposing a monokinetic closure [32].
In Figure 3 we have that $\lambda_{11}=.1,\lambda_{22}=1,\lambda_{12}=10$ and
$P_{12}^{11}=0.2\dfrac{1}{4}\left[1-\exp^{-v}+1-\exp^{-w}\right],P_{12}^{22}=0.8\dfrac{1}{4}\left[1-\exp^{-v}+1-\exp^{-w}\right]$.
We also present a comparison with the solution of macroscopic equations
(63)-(65)-(64) (straight lines) where we consider constant switching
probabilities $P_{12}^{11}=0.2,P_{12}^{22}=0.8.$ We observe that both the
microscopic model with $v$-dependent switching probabilities and the
macroscopic model with constant switching probabilities forecast a similar
behavior of the macroscopic quantities in the long run. The microscopic model
with a $v$-dependent switching probability forecasts the same behaviour but
with a delay, and this is due to the fact that the switching probabilities are
smaller than the constant ones. In magenta and green we also present the
results of the simulation of the microscopic model in case
$P_{12}^{11}=0.2\dfrac{1}{2}\left[\exp^{-v}+\exp^{-w}\right],P_{12}^{22}=0.8\dfrac{1}{2}\left[\exp^{-v}+\exp^{-w}\right]$.
In this case we can observe that the convergence is even slower. This is due
to the fact that for large values of the wealth $v$, the switching probability
is very small.

Figure 3: Test 3. Solution to the microscopic model (5)-(6)-(8)-(30) as
illustrated in Sec. 2.4 with the microscopic rules (60)-(62) with the Monte
Carlo algorithm 1 (circles). Here
$\lambda_{11}=1,\lambda_{22}=.1,\lambda_{12}=10$ and
$P_{12}^{11}=0.2\dfrac{1}{2}\left[1-\exp^{-v}+1-\exp^{-w}\right],P_{12}^{22}=0.8\dfrac{1}{2}\left[1-\exp^{-v}+1-\exp^{-w}\right]$.
We also present a comparison with the solution of macroscopic equations
(63)-(65)-(64) (straight lines) where we consider constant switching
probabilities $P_{12}^{11}=0.2,P_{12}^{22}=0.8.$ In magenta and green we also
present the results of the simulation of the microscopic model in case
$P_{12}^{11}=0.2\dfrac{1}{2}\left[\exp^{-v}+\exp^{-w}\right],P_{12}^{22}=0.8\dfrac{1}{2}\left[\exp^{-v}+\exp^{-w}\right]$.

###  5.2 Quasi-invariant regime and Fokker-Planck equation

In this section, we consider the quasi-invariant regime (51), i.e. we analyse
the dynamics on a long time-scale by considering small exchanges of wealth,
while the switching probability is not rescaled. In this framework, we have
seen that it is possible to approximate the leading order of the stationary
solution $f_{1}^{(0)}$ through (78). Here, we compare the stationary state
(78) with the solution $f_{1}^{MC}$ obtained by the numerical integration of
the microscopic process (5)-(6)-(8)-(71) with the microscopic rules (60)-(62),
where in the quasi-invariant regime (71) we have chosen $\epsilon=10^{-3}$. In
Figure 4 we represent the analytical $\tilde{f}_{1}^{(0)}$ as given in (78)
and the approximation of $f_{1}^{MC}$. We can remark that, despite the fact
that $f_{1}^{MC}$ is obtained through a Monte Carlo simulation and
$\tilde{f}_{1}^{(0)}$ is an approximation, the agreement is quite good. In the
right panel, we also represent the numerical approximation of the distribution
function $g_{1}^{MC}$ of the first population, in case we consider a quasi-
invariant transition probability defined by (51). With this choice it is
granted that the evolution of both the average and the energy in the quasi-
invariant regime $\epsilon\ll 1$ is the same as in the standard regime defined
by $\epsilon=1$. In this case it was not possible to determine the stationary
state explicitly like in the quasi-invariant regime leading to (78) that, on
the other hand, only grants that the average is the same for each $\epsilon$.

Figure 4: Solution of the Fokker-Planck equation (76). Left: comparison
between the approximated stationary state $\tilde{f}_{1}^{(0)}$ given by (78)
and the solution $f_{1}^{MC}$ of the microscopic model (5)-(6)-(8)-(71) with
the microscopic rules (60)-(62), with $\epsilon=10^{-3}$ in the quasi-
invariant regime (71) obtained with $N=10^{6},\Delta t=10^{-3}$. Right:
comparison between $f_{1}^{MC}$ and $g_{1}^{MC}$ obtained integrating the
same, with the quasi-invariant regime (51).

##  6 Conclusions

In this paper we have presented a general framework for modeling systems of
interacting particles with multiple microscopic states changing simultaneously
according to a given dynamics. In particular, the microscopic description
relies on Markovian processes described by transition probabilities, as they
depend on the pre-interaction states, and the interaction frequency depends on
the microscopic states. The fact of starting from the microscopic stochastic
process allows to describe in more detail the dynamics, by including
parameters and quantities related to the phenomenon under study that can be
observed. The derivation, through kinetic equations, of macroscopic equations
allows to obtain also at the aggregate level a higher level of detail that is
inherited from the underlying microscopic dynamics.

Under some assumptions, general results concerning well-posedness, existence
and uniqueness of a solution for the Cauchy problem associated to our kinetic
equation have been shown. We have also rephrased the concept of quasi-
invariant limit in the present framework, leading to evolution equations of
Fokker-Planck type.

We have applied the present modeling framework in order to describe systems of
binarily interacting agents characterized by a physical quantity $v$
(representing wealth, or opinion, or viral load, etc.) and by a label $x$
denoting the belonging to a given subgroup. The physical quantity changes
according to binary interaction rules, while the label changes, simultaneously
with the physical quantity, through a switch process caused by the same binary
interaction. In this context, we have seen that the description of the
microscopic process by means of transition probabilities allows us to remove
the reversibility assumption on the interaction rule, modeling thus also
stochasticity in the binary encounters giving rise to transfers (not present
in the paper [4] using classical Boltzmann operators analogous to the reactive
ones). Moreover, in our framework it is easier to consider a non-constant
switching probability, i.e. depending on the microscopic physical quantity. We
have analyzed and discussed various quasi-invariant regimes and performed some
numerical tests showing a very good agreement between the microscopic Monte
Carlo simulations and the derived macroscopic equations.

The modelling framework investigated in this paper is worth to be applied and
generalized to many other problems. As first, the model for international
trade with transfers presented in Section 4 could be extended by adding an
extra independent microscopic process for $v$, describing the exchange of
goods without transfers; this would make the model even more similar to the
kinetic description of gaseous mixtures, where elastic collisions (which do
not change the nature of the particles) coexist with chemical reactions
(changing the species of the reacting particles). Epidemic models based on a
kinetic approach could be improved owing to our stochastic framework with
multiple states as it allows to start from a microscopic description and to
consider independent or simultaneous microscopic stochastic dynamics for the
different variables of the microscopic state. For example, the so-called
“non–conservative” interactions giving rise to the passage from one
compartment to another could be made more realistic taking into account also
the simultaneous change of viral load (of individuals) as done in [15, 16] or
internal activity (of cells). Moreover, the relation with kinetic models with
label switching and gradient descent could be established [10]. Eventually,
applications to situations with many internal states is the final scope of our
framework. It could provide for instance a physically reasonable description
of mixtures of polyatomic gases, with each molecule characterized by its
species label $i$, its velocity ${\bf v}\in\mathbb{R}^{3}$, and its internal
energy that could also be separated into the vibrational part (typically
described by a discrete variable) and the rotational part (typically
approximated by a continuous variable) [26, 9]. Even in econophysics, the
possible influence of the personal knowledge of the market on the strategy
adopted in the trades (as sketched in [33] for a single population) could be
described considering the individual knowledge as an additional microscopic
state, besides the population label and the individual amount of wealth.
Suitable quasi-invariant limits and properties of steady states of such non-
standard kinetic descriptions of various interacting populations are
completely open problems worth to be investigated in future research.

Acknowledgments This research was initiated during the post-doc contract of
N.L. at the Department of Mathematical, Physical and Computer Sciences of
Parma University, funded by the Italian National Research Project “Multiscale
phenomena in Continuum Mechanics: singular limits, off-equilibrium and
transitions” (Prin 2017YBKNCE). The authors also thank the support by
University of Parma, by Politecnico di Torino, by the Italian National Group
of Mathematical Physics (GNFM-INdAM), and by the Italian PRIN Research Project
“Integrated Mathematical Approaches to Socio–Epidemiological Dynamics” (Prin
2020JLWP23, CUP: E15F21005420006).

## Appendix: Nanbu-Babovski algorithm

Data:

  * •

$N\in\mathbb{N}$ total number of agents of the system;

  * •

$N_{1}^{n},\,N_{2}^{n}\in\mathbb{N}$ numbers of agents in $x=1$, $x=2$,
respectively, at time $t^{n}:=n\Delta{t}$ and $v_{1}^{n},\,v_{2}^{n}$ the
microscopic states of agents in $x=1$, $x=2$, respectively, at time $t^{n}$;

1 Fix $\Delta{t}\leq\min\\{\frac{1}{\max\lambda_{ij}}\\}$;

2 for  _$n=0,\,1,\,2,\,\dots$_ do

3 Compute
$\rho_{1}^{MC,n}=\dfrac{N_{1}^{n}}{N},\qquad\rho_{2}^{MC,n}=\dfrac{N_{2}^{n}}{N},\qquad
m_{1}^{MC,n}=\dfrac{1}{N_{1}^{n}}\displaystyle{\sum_{k=1}^{N_{1}^{n}}}v_{k}^{n},\qquad
m_{2}^{MC,n}=\dfrac{1}{N_{2}^{n}}\displaystyle{\sum_{k=1}^{N_{2}^{n}}}v_{k}^{n}$;

4 repeat

5 Pick randomly two agents $(x_{i}^{n},\,v_{i}^{n})$,
$(x_{j}^{n},\,v_{j}^{n})$ with $i\neq j$;

6 for  _$h=i,\,j$_ do

7 Sample
$\Theta\sim\operatorname{Bernoulli}(\lambda_{x_{i}^{n}x_{j}^{n}}\Delta{t})$;

8 if  _$\Theta=1$_ then

9 for  _$\\{x^{\prime}_{i},\,x^{\prime}_{j}\\}\in\mathcal{I}_{n}^{2}$_ do

10 Sample $J\in\\{1,\,0\\}$ with law
$\operatorname{Prob}(J=1)=P_{x_{i}^{n},x_{j}^{n}}^{x_{i}^{\prime}x_{j}^{\prime}},\quad\operatorname{Prob}(J=0)=1-P_{x_{i}^{n},x_{j}^{n}}^{x_{i}^{\prime}x_{j}^{\prime}}$;
if  _$J=1$_ then

11 Set $(x_{i}^{n+1},x_{j}^{n+1})=(x_{i}^{\prime},x_{j}^{\prime})$;
Set $(v_{i}^{n+1},\,v_{j}^{n+1})=(v_{i}^{\prime},v_{j}^{\prime})$ where
$(v_{i}^{\prime},v_{j}^{\prime})$ is given by (10)-(60) and break

12

13

14 else

15 Set $x_{h}^{n+1}=x_{h}^{n},\,v_{h}^{n+1}=v_{h}^{n}$;

16

17

18 until  _no unused pairs of agents are left_ ;

Algorithm 1 Nanbu-Babovski algorithm with mass transfer for model (5)-(6)-(8)-
(60)-(62)

## References

  * [1] L. Ambrosio, N. Gigli, and G. Savarè.  Gradient flows in metric spaces and in the space of probability measures.  Birkhauser Verlag, Basel, 2008.
  * [2] P. Andries, K. Aoki, and B. Perthame.  A consistent bgk-type model for gas mixtures.  J. Stat. Phys., 106:993–1018, 2002.
  * [3] M. Bisi.  Some kinetic models for a market economy.  Boll. Unione Mat. Ital., 10:143–158, 2017.
  * [4] M. Bisi.  Kinetic model for international trade allowing transfer of individuals.  Phil. Trans. A, 380:20210156. (pp. 1–14), 2022.
  * [5] M. Bisi, M. Groppi, and G. Spiga.  Kinetic Bhatnagar-Gross-Krook model for fast reactive mixtures and its hydrodynamic limit.  Phys. Rev. E, 81:036327 (pp. 1–9), 2010.
  * [6] A.V. Bobylev, J.A. Carrillo, and I. Gamba.  On some properties of kinetic and hydrodynamics equations for inelastic interactions.  J. Stat. Phys., 98:743–773, 2000.
  * [7] V.C. Boffi, V. Protopopescu, and G. Spiga.  On the equivalence between the probabilistic, kinetic, and scattering kernel formulations of the Boltzmann equation.  Physica A, 164:400–410, 1990.
  * [8] R. Borsche, A. Klar, and M. Zanella.  Kinetic-controlled hydrodynamics for multilane traffic models.  Physica A, 587:126486 (pp. 1–17), 2022.
  * [9] T. Borsoni, M. Bisi, and M. Groppi.  A general framework for the kinetic modelling of polyatomic gases.  Comm. Math. Phys., 393:215–266, 2022.
  * [10] M. Burger and A. Rossi.  Analysis of kinetic models for label switching and stochastic gradient descent.  Preprint: arXiv.2207.00389, 2022.
  * [11] C. Cercignani.  The Boltzmann Equation and its Applications.  Number 67 in Applied Mathematical Sciences. Springer, New York, 1988.
  * [12] S. Chapman and T.G. Cowling.  The Mathematical Theory of Non-Uniform Gases.  Cambridge University Press, Cambridge, 1970.
  * [13] S. Cordier, L. Pareschi, and G. Toscani.  On a kinetic model for a simple market economy.  J. Stat. Phys., 120(1):253–277, 2005.
  * [14] M. Delitala and T. Lorenzi.  A mathematical model for value estimation with public information and herding.  Kinet. Relat. Models, 7(1):29–44, 2014.
  * [15] R. Della Marca, N. Loy, and A. Tosin.  An SIR-like kinetic model tracking individuals’ viral load.  Networks and Heterogeneous Media, 17:467–494, 2022.
  * [16] R. Della Marca, N. Loy, and A. Tosin.  An SIR model with viral load-dependent transmission.  Preprint: arXiv:2208.12004, 2022.
  * [17] R. Della Marca, M.D.P. Machado Ramos, C. Ribeiro, and A.J. Soares.  Mathematical modelling of oscillating patterns for chronic autoimmune diseases.  Math. Meth. Appl. Sci., 45:7144–7161, 2022.
  * [18] G. Dimarco, L. Pareschi, G. Toscani, and M. Zanella.  Wealth distribution under the spread of infectious diseases.  Phys. Rev. E, 102(2):022303, 2020.
  * [19] B. During.  Multi-species models in econo- and sociophysics.  Econophysics & economics of games, social choices and quantitative techniques (eds. B. Basu, B.K. Chakrabarti, S.R. Chakravarty, K. Gangopadhyay), pages 83–89, 2010.  Dordrecht: Springer.
  * [20] B. During, P. Markowich, J.F. Pietschmann, and M.T. Wolfram.  Boltzmann and Fokker-Planck equations modelling opinion formation in the presence of strong leaders.  Proc. R. Soc. A, 465:3687–3708, 2009.
  * [21] B. During and G. Toscani.  International and domestic trading and wealth distribution.  Commun. Math. Sci., 6:1043–1058, 2008.
  * [22] A. Festa, A. Tosin, and M.T. Wolfram.  Kinetic description of collision avoidance in pedestrian crowds by sidestepping.  Kinet. Relat. Models, 11:491–520, 2018.
  * [23] P. Freguglia and A. Tosin.  Proposal of a risk model for vehicular traffic: A Boltzmann-type kinetic approach.  Commun. Math. Sci., 15(1):213–236, 2017.
  * [24] V. Giovangigli.  Multicomponent Flow Modeling.  Birkhäuser, Boston, 1999.
  * [25] C. D. Greenman and T. Chou.  Kinetic theory of age-structured stochastic birth-death processes.  Phys. Rev. E, 93(1):012112, 2016.
  * [26] G. Herzberg.  Molecular Spectra and Molecular Structure.  Van Nostrand Reinold, New York, 1950.
  * [27] M.N. Kogan.  Rarefied Gas Dynamics.  Plenum Press, New York, 1969.
  * [28] N. Loy, T. Hillen, and K. Painter.  Direction dependent turning leads to anisotropic diffusion and persistence.  European Journal of Applied Mathematics, 33(4):729–765, 2022.
  * [29] N. Loy and A. Tosin.  Markov jump processes and collision-like models in the kinetic description of multi-agent systems.  Commun. Math. Sci., 18(6):1539–1568, 2020.
  * [30] N. Loy and A. Tosin.  Boltzmann-type equations for multi-agent systems with label switching.  Kinet. Relat. Models, 14(5):867–894, 2021.
  * [31] M.D.P. Machado Ramos, C. Ribeiro, and A.J. Soares.  A kinetic model of T-cell autoreactivity in autoimmune diseases.  J. Math. Biol., 79:2005–2031, 2019.
  * [32] L. Pareschi and G. Toscani.  Interacting Multiagent Systems: Kinetic equations and Monte Carlo methods.  Oxford University Press, 2013.
  * [33] L. Pareschi and G. Toscani.  Wealth distribution and collective knowledge: a Boltzmann approach.  Phil. Trans. A, 372:20130396, 2014.
  * [34] M. Pirner.  A review on BGK models for gas mixtures of mono and polyatomic molecules.  Fluids, 6:393, 2021.
  * [35] G. Puppo, M. Semplice, A. Tosin, and G. Visconti.  Kinetic models for traffic flow resulting in a reduced space of microscopic velocities.  Kinet. Relat. Models, 10(3):823–854, 2017.
  * [36] A. Rossani and G. Spiga.  A note on the kinetic theory of chemically reacting gases.  Physica A, 272:563–573, 1999.
  * [37] G. Spiga, T. Nonnenmacher, and V.C. Boffi.  Moment equations for the diffusion of the particles of a mixture via the scattering kernel formulation of the nonlinear Boltzmann equation.  Physica A, 131:431–448, 1985.
  * [38] G. Toscani.  Kinetic models of opinion formation.  Commun. Math. Sci., 4(3):481–496, 2006.
  * [39] G. Toscani, A. Tosin, and M. Zanella.  Opinion modeling on social media and marketing aspects.  Phys. Rev. E, 98:022315, 2018.
  * [40] G. Toscani, A. Tosin, and M. Zanella.  Multiple-interaction kinetic modeling of a virtual-item gambling economy.  Phys. Rev. E, 100:012308, 2019.
  * [41] C. Villani.  On a new class of weak solutions to the spatially homogeneous Boltzmann and Landau equations.  Arch. Ration. Mech. Anal., 143(3):273–307, 1998.
  * [42] L. Waldmann.  Transporterscheinungen in gasen von mittlerem druck.  Handbuch der Physik, 12:295–514, 1958.  S. Flügge ed., Springer Verlag, Berlin.