# CovarianceNet: Conditional Generative Model for Correct Covariance
Prediction in Human Motion Prediction

Aleksey Postnikov1,2 Aleksander Gamayunov1,2 Gonzalo Ferrer2 1 The authors are
with the Sberbank Robotics Laboratory, Moscow, Russia.
<EMAIL_ADDRESS>2Skolkovo Institute of Science and
Technology, Moscow, Russia<EMAIL_ADDRESS>
###### Abstract

The correct characterization of uncertainty when predicting human motion is
equally important as the accuracy of this prediction. We present a new method
to correctly predict the uncertainty associated with the predicted
distribution of future trajectories. Our approach, CovariaceNet, is based on a
Conditional Generative Model with Gaussian latent variables in order to
predict the parameters of a bi-variate Gaussian distribution. The combination
of CovarianceNet with a motion prediction model results in a hybrid approach
that outputs a uni-modal distribution. We will show how some state of the art
methods in motion prediction become overconfident when predicting uncertainty,
according to our proposed metric and validated in the ETH data-set [1].
CovarianceNet correctly predicts uncertainty, which makes our method suitable
for applications that use predicted distributions, e.g., planning or decision
making.

##  I Introduction

Human Motion Prediction, during the last years, has received the attention of
the research community from different fields: intelligent vehicles, pattern
recognition, graphics, robotics, etc. The motivation to understand and predict
human motion is immense and it has a deep impact in related topics, such as,
decision making, path planning, autonomous navigation, surveillance, tracking,
scene understanding, anomaly detection, etc.

The problem of forecasting where pedestrians will be in the near future is,
however, ill-posed by nature: human beings tend to be unpredictable on their
decisions and motion is neither exempt of it. These random nature in motion
brings an open challenge to prediction algorithms, where algorithms are
desired to be accurate and correctly grasp the uncertainty associated with
their predictions.

To this end, multiple benchmarks have been created and released [1, 2, 3],
providing common grounds to test and evaluate. Most modern motion prediction
algorithms focus on accurate prediction of agent position errors on these
benchmarks. Nonetheless, the precision due to this inherent uncertainty is
equally important, and this paper is an effort to research on this direction.
Prediction algorithms should address this issue as well: it provides a high
degree of interpretability by estimating the associated uncertainty and it
might be of some use for consequent algorithms making use of prediction
information, e.g., planning. For example, Fig. 1 shows two predictions, of
similar predicted error, but different uncertainty estimation.

Figure 1: Comparison of two predicted trajectories in the ETH dataset[1].
CovarianceNet in cyan stars, blue and green ellipses, - mean, 1 sigma, 3
sigma, respectively. Trajectron++ [4] in yellow stars, yellow and red
ellipses, red ellipses - mean, 1 sigma, 3 sigma, respectively. The environment
is challenging, since there is a pedestrian in the way.

In this paper, we focus on predicting uncertainty and human motion prediction.
We propose a hybrid approach shown on Fig. 2, consisting of deep model for
Goals Prediction and model-based trajectory prediction, complemented by a
modern neural-net approach to predict motion as a uni-modal Gaussian
distribution. Several works provide a multi-modal distribution for motion
prediction such as [4], achieving excellent results in benchmarks. However,
combining a mixture distribution into robot planning approaches (for instance)
requires careful considerations. Sampling based techniques require extra
attention when considering multi-modality since low-probability prediction
outcomes might result in dangerous outcomes, when evaluated under a risk
perspective [5].

The contributions of this work are:

  * •

A Conditional Generative Model to correctly predict covariances;

  * •

Hybrid approach combining a model-based motion prediction from the Social
Force Model(SFM) [6] and our implementation of a learning-based Goal
Prediction;

  * •

Metric to measure the correctness of the uncertainty prediction by counting
and averaging the number of times that the prediction lies in the iso-contours
of a bi-variate Gaussian.

Figure 2:  Our Hybrid approach for Human Motion Prediction and Uncertainty
Estimation. We combine customized deep goal prediction model[7] with SFM[6],
for mean pedestrian pose prediction and Covariance Net uses Conditional
Variational Autoencoder[8] deep generative model for uncertainty prediction.
Figure 3: CovarianceNET architecture. CovarianceNET is a part of our proposed
hybrid approach for Human Motion Prediction and Uncertainty Estimation. The
input consists of the spatial coordinates over the past $T_{obs}$ seconds of
all agents in the scene. We use LSTM-based encoder for the agents’ history.
Neighbours’ impact on the predicted motion is encoded by adding an attention
module. We use the CVAE latent variable framework for diverse but realistic
uncertainty prediction[8].

##  II Related Work

Motion prediction has been studied by the robotics community mainly motivated
by its direct relation to robot navigation in social environments, embedded at
its deepest core. Examples of robot navigation include Human-aware approaches
[9, 10], considering prediction into planning [11, 12, 13, 14] or the effect
of planning to prediction [15, 16, 17].

Other fields have studied human motion prediction, where multiple previous
works study the problem of how to forecast pedestrian trajectories and predict
their future behaviors. Broadly, human motion prediction can be divided in two
classes: model-based that creates an empirical model of these transition
functions such as the SFM[6] and its variants[18, 19] or models proposed from
the graphics community [20, 21]. And Learning-based approaches [4, 22, 23, 7,
24, 25] that are becoming the dominant paradigm in human motion prediction, as
well as in other topics due to its unrivaled results. Both of these approaches
have the same input and output data and predict the next state of the
pedestrians in $dt$ time.

Helbing and Molnár [6] have proposed a model based approach for modelling
pedestrians’ behaviour. Authors have shown that pedestrian motion can be
described as a sum of social forces that depend on agents destination point,
other pedestrians, static object (borders of buildings, walls, streets,
obstacles).

Alahi et al. [22] propose a learning-based method with social pooling, where
the Long Short-term Memory(LSTM)[26] states of neighboring agents were pooled
based on their locations in a 2-D grid to form social tensor.

Messaoud et al. [23] apply a multihead attention mechanism to the social
tensor to directly relate distant vehicles and extract a context
representation.

Recently, generative approaches have emerged as state-of-the-art trajectory
forecasting methods due to recent advancements in deep generative models. They
have caused a shift from focusing on predicting the single best trajectory to
producing a distribution of potential future trajectories. Most works in this
category use a deep recurrent backbone architecture with a latent variable
model, such as a Conditional Variational Autoencoder (CVAE) [8] or a
Generative Adversarial Network (GAN) [27].

Trajectron++ [4, 28] \- a multi-agent behavior prediction model that accounts
for the dynamics of the agents, produces predictions possibly conditioned on
potential future robot trajectories which can effectively use heterogeneous
data about the surrounding environment.

Jein et al. [29] proposes a discrete residual flow to recursively updates the
predicted distribution over a pedestrian’s future position in the form of
occupancy maps.

Mahgalam et al. [30] propose to address human trajectory prediction by
modeling intermediate stochastic goals proposing a socially compliant,
endpoint conditioned variational auto-encoder with a novel self-attention
based social pooling layer.

In our previous work [31], we propose a new method for motion prediction -
HSFM-$Σ$NN that combines two different approaches: a feed-forward network
whose layers are model-based transition functions using the Headed Social
Force Model(HSFM) and a Neural Network for covariance prediction.
CovarianceNet is the next step in this work, where we have improved the
approach substantially.

Gal et al., [32] proposed a simple yet effective method for probabilistic
interpretation of dropout which allows to obtain model uncertainty out of
existing deep learning models.

Guo et al. [33] introduced a temperature scaling, confidence calibration
method, that can effectively correct the miscalibration in modern deep neural
networks.

##  III Method

###  III-A Problem formulation

The position of a generic agent $i$ at time $t$ is represented by
$\bm{x}_{i}^{t}=(x_{i},y_{i})_{t}$, where $(x_{i},y_{i})_{t}$ are the
coordinates of agents in the global reference system at the instant of time
$t$. The agent’s trajectory is defined as
$X_{i}^{1:T}=\\{\bm{x}_{i}^{1},...,\bm{x}_{i}^{T}\\}$ from timestamp $1$ to
$T$. We aim to generate plausible trajectory distributions for a time-varying
number of interacting agents. Every trajectory is split into observed and
future: given certain number $T_{obs}$ of observed positions, we seek a
distribution over all agents’ future states for the next $T_{pred}$ time steps
which is denoted as $p(X_{i}^{T_{obs}+1:T_{pred}}|X_{i}^{1:T_{obs}})$.

Our approach is visualized in Fig. 2. The first block is the model-based SFM
[6], a method for motion prediction the mean positions of the agents’ future
states.

$\textit{sfm}^{T_{obs}+1:T_{pred}}_{1:N}=\textit{SFM}(X^{1:T_{obs}}_{1:N},\hat{X}^{T}_{goal_{1:N}})$
(1)

At a higher level the SFM needs to infer future possible destinations of
pedestrians [34]. In particular, we have customized the learning-based
approach by [7] and goals are predicted as:

$\hat{X}^{T}_{goal_{1:N}}=\phi(X^{1:T_{obs}}_{1:N}),$ (2)

where $\phi$ \- deep multihead-attention based model that trained to predict
the position of pedestrian at timestamp $T_{pred}$, $X^{1:T_{obs}}_{1:N}$ \-
the trajectory (scene history) which contains the states of all $N$ agents at
timestamp $t$, $X^{t}_{1:N}=\\{X^{t}_{1},..,X^{t}_{N}\\}$, $P_{goal_{1:N}}$ \-
the position of the pedestrian at timestamp $T_{pred}$, $\textit{sfm}^{t}_{n}$
\- predicted position of $n$ agent at timestamp $t$. Here we emphasize that in
this work our main goal is to show correct covariance prediction, thus we are
not showing implementation details of Goal Prediction model and SFM.

The CovarianceNet network is one of the sub-blocks of our approach, depicted
in Fig. 2 and it is shown in detail in Fig. 3. CovarianceNet is designed to
complement the predictions by the SFM with accurate $x,y$ variances and its
correlation coefficient, such that the bi-variate Gaussian constructed from
the mean and covariance matrix accurately assesses the potential trajectory
distribution of pedestrians, and statistically ensure the properties of a bi-
variate Gaussian distribution.

###  III-B Agent encoder

For each pedestrian in the scene, the current position and all observed
previous positions are known. In our work (see Fig. 3), we are using the Agent
encoder to encode information about agent position and the Neighbours encoder
for encoding information about the influence of the agent and all his
neighbors to each other. Knowing the history of the agent’s movement, one can
calculate velocities and accelerations as a additional source of data.

The observed trajectory of each agent $X^{1:T_{obs}}_{1:N}$ is encoded using
an LSTM encoder

$h_{traj_{i}}^{t}=\textit{LSTM}(X^{1:T_{obs}}_{1:N},h_{traj_{i}}^{t-1},W_{t_{enc}})$
(3)

Here, $h_{traj_{i}}^{t}$ is the hidden state vector of the $i^{th}$ agent at
time $t$. All LSTM encoders share the same weights $W_{t_{enc}}$.

Additionally, we use velocities and accelerations encoded in the same fashion
as in (3) and sum them into a single encoded vector:

$h_{acc_{i}}^{t}=\textit{LSTM}(\frac{\partial}{\partial
t}X^{1:T_{obs}}_{1:N},h_{i_{vel}}^{t-1},W_{v_{enc}})$ (4)
$h_{vel_{i}}^{t}=\textit{LSTM}(\frac{\partial^{2}}{\partial
t^{2}}X^{1:T_{obs}}_{1:N},h_{i_{vel}}^{t-1},h_{i_{acc}}^{t-1},W_{a_{enc}})$
(5) $E_{hist_{i}}^{t}=h_{traj_{i}}^{t}+h_{vel_{i}}^{t}+h_{acc_{i}}^{t};$ (6)

where $h^{t}_{i_{traj}}$ \- encoded trajectory, $h^{t}_{i_{vel}}$ \- encoded
velocities, $h^{t}_{i_{acc}}$ \- encoded accelerations and $E^{t}_{hist}$ \-
full encoded agent’s history.

During the training phase, as it is shown in Fig.3, to produce a latent
distribution additionally leveraged future trajectory encoding
$E^{t}_{fut_{i}}$ for time indexes ${T_{obs}+1...T_{pred}}$, in the same
manner as in Eq.(3).

###  III-C Neighbours encoding

To model neighboring agents’ influence on the modeled agent, all agents’
neighbours in the scene are processed with the additive attention module[35]
and are encoded by the LSTM cell to produce a single neighboring agents’
influence vector $E_{neigh_{i}}^{t}$.

$C_{i}^{t}=\textit{Attention}(X^{1:T_{obs}}_{neigh_{i}},X^{1:T_{obs}}_{i})$
(7) $E_{neigh_{i}}^{t}=\textit{LSTM}(C_{i}^{t},E_{neigh_{i}}^{t-1},W_{neigh})$
(8) $E_{scene_{i}}^{t}=E_{neigh_{i}}^{t}+E_{hist_{i}}^{t}$ (9)

where $C_{i}^{t}$ \- context vectors, $E_{neigh_{i}}^{t}$ \- encoded
neighbours influence on the modeled agent, $E_{scene_{i}}^{t}$ is a full
encoded scene history.

###  III-D Conditional Variational Autoencoder

The network backbone of our approach (Fig. 3) is realized as a version of the
Conditional Variational Autoencoder[8] (CVAE).

The CVAE architecture can be divided in two parts: the prior and the posterior
networks. Both prior and posterior distributions are assumed to be Normal
distributions. The parameters of the prior are computed by the prior network
which only takes the encoded history as input. The posterior parameters are
determined from both the encoded history and the target trajectory. The prior
distribution is a Normal distribution, denoted as
$p_{\phi_{i}}(z|X^{1:T_{obs}}_{i})=\mathcal{N}(\mu_{prior_{i}}(x),\sigma^{2}_{prior_{i}}(x))$.

$p_{\phi_{i}}(z|X^{1:T_{obs}}_{i})=\mathcal{N}(\mu_{prior_{i}},\sigma^{2}_{prior_{i}})$
(10) $\mu_{prior_{i}}=\textit{MLP}(E_{scene_{i}})$ (11)
$\sigma_{prior_{i}}=\textit{MLP}(E_{scene_{i}})$ (12)

During training, the latent variable $z$ will be sampled from the posterior
distribution. Specifically, it takes both the encoded past and the encoded
future trajectories information $X_{i}$ passed through an MLP to obtain a
latent mean $\mu_{latent}$ and a latent sigma $\sigma_{latent}$ to output a
latent distribution
$q_{\phi_{i}}(z|X^{T_{obs}+1:T_{pred}}_{i},X^{1:T_{obs}}_{i})$.

$q_{\phi_{i}}(z|X^{T_{obs}+1:T_{pred}}_{i},X^{1:T_{obs}}_{i})=\mathcal{N}(\mu_{latent_{i}},\sigma^{2}_{latent_{i}})$
(13) $\mu_{latent_{i}}=\textit{MLP}(E_{scene_{i}},E_{fut_{i}})$ (14)
$\sigma_{latent{i}}=\textit{MLP}(E_{scene_{i}},E_{fut_{i}})$ (15)

In our work, we use two fully connected layers with dropout and sigmoid
activation function between them as base multi layer perceptron (MLP).

###  III-E Decoder

The decoder models the probability of a target trajectory
$X^{T_{obs}+1..T_{pred}}$ given the latent variable $z$ sampled from the
latent distribution, encoded node history $E_{scene_{i}}^{t}$ and SFM
trajectory predictions $\textit{sfm}^{t}_{i}$.

We make an assumption that the target distribution
$p(\bm{x}^{t}|X^{1:T_{obs}})$ is a bi-variate Gaussian:

$p_{i}(\bm{x}^{t}_{i}|X^{1:T_{obs}}_{i})=N(\mu^{t}_{i},\Sigma^{t}_{i})$ (16)
$\mu^{t}_{i}=\begin{bmatrix}\mu_{x_{i}}^{t}\\\ \mu_{y_{i}}^{t}\\\
\end{bmatrix};\Sigma^{t}_{i}=\begin{bmatrix}\sigma_{x_{i}}^{t^{2}}&\rho\sigma_{x_{i}}^{t}\sigma_{y_{i}}^{t}\\\
\rho\sigma_{x_{i}}^{t}\sigma_{y_{i}}^{t}&\sigma_{y_{i}}^{t^{2}}\\\
\end{bmatrix}.$ (17)

For each timestamp of prediction, we use encoded information about current
agent state, passed through the MLP as an input to the GRU layer with
$E_{scene}$ as the initial state of the GRU cell, witch recurrently outputs
the new hidden states for each agent.

$h^{t}_{gru_{i}}=\textit{GRU}([\textit{sfm}^{t-1}_{i},z],h^{t-1}_{gru_{i}})$
(18)

GRU hidden states are then used to predict the parameters of a bi-variate
Gaussian distribution $\mathcal{N}(\mu,\sigma,\rho)$ standard deviations
$\sigma$ and correlation coefficient $\rho$, while the mean values $\mu$ are
taken from the SFM prediction $\mu_{i}=\textit{sfm}^{t}_{i}$.

$\rho^{t}_{i}=\textit{MLP}(h^{t}_{gru_{i}})$ (19)
$\sigma^{t}_{i}=\textit{MLP}(h^{t}_{gru_{i}})$ (20)

The entire model is trained by maximising the sequential evidence lower-bound
(ELBO):

$L_{nll}(W)=-\sum^{N}_{i=1}\sum^{T_{pred}}_{t=T_{obs}+1}log(p^{t}_{i}(\bm{x}^{t}_{i}|\mu^{t}_{i},\sigma^{t}_{i},\rho^{t}_{i}))$
(21)

$L^{n}_{kl}=\sum^{N}_{i=1}D_{kl}(q_{\phi_{i}}(z|X^{T_{obs}+1:T_{pred}}_{i},X^{1:T_{obs}}_{i})||N(0,1))$
(22) $\textit{Loss}=\alpha L^{n}_{nll}(W)+L^{n}_{kl}$ (23)

It is important to mention that in our case, when we utilize Goal Predictor
and SFM[6] as model for mean poses prediction we do not optimize $\mu^{t}_{i}$
parameter from (21).

Figure 4: Evaluation of PPEI. Left: results for $1\sigma$. Right: results for
$3\sigma$.

##  IV  Evaluation

In our work, we use the combination of goal predictor and SFM[6] modules as a
base predictor for mean future poses of agents and combine it with the
CovarianceNet in order to predict accurate uncertainties. We compare results
of our method with other three popular approaches for evaluating uncertainty
on publicly-available ETH[1] pedestrian dataset. It consists of real world
human trajectories with rich multi-human interaction scenarios. This dataset
is a standard benchmark in the field as it contain challenging behaviors such
as couples walking together, groups crossing each other, and groups forming
and dispersing. Leave-one-out strategy was used for evaluation, where the
model is trained on four datasets and evaluated on the held-out fifth. An
observation length of 8 timesteps (3.2s) and a prediction horizon of 12
timesteps (4.8s) is used for evaluation. In this section we will compare
results of our method with the following methods:

###  IV-A Evaluated Methods

####  IV-A1 Covariance Forward-Propagation (FP)

In our work we have used SFM for prediction forces acting on all agents in
scene and Human locomotion model which integrates predicted forces to the
state variables of a pedestrians 2D poses at timestamp $t+1$, in the following
generalized way:

$x_{t+1}=T(x_{t})$ (24)

The transition function $T()$ defined by the SFM (24) is a non-linear
differentiable function (by construction). The simplest method for covariance
estimation is using the first-order Taylor expansion:

$x_{t+1}=T(\mu_{t})+G_{t}(x_{t}-\mu_{t})$ (25)

where $\mu_{t}$ is the current state estimate and $G_{t}$ is the Jacobian
matrix $\partial T/\partial x$ evaluated at $\mu_{t}$. From here, we apply
Covariance Propagation of a Gaussian random variable
($x_{t}\sim\mathcal{N}(\mu_{t},\Sigma_{t})$) over a linear function:

$x_{t+1}\sim\mathcal{N}\big{(}T(\mu_{t}),G_{t}\cdot\Sigma_{t}\cdot
G_{t}^{\top}\big{)}$ (26)

####  IV-A2 Social-LSTM [22]

One of the most influential works in direction of Human Motion Distribution
Prediction based on Neural Networks, where each agent is modeled with an LSTM
and nearby agents’ hidden states are pooled at each timestep using a proposed
social pooling operation.

####  IV-A3 Trajectron++ [4]

One of the latest works in this direction is named Trajectron++[4] has four
configurations of predictions, of which we used Most Likely, which gives the
best results of ADE and FDE. Trajectron++ calculates the result of predictions
as a Gaussian Mixture Model (GMM) witch contains 25 Gaussian distributions.

###  IV-B Motion Prediction Evaluation

We use Euclidean distance errors, Average Displacement Error (ADE) and Final
Displacement Error (FDE) to evaluate the accuracy of then mentioned approaches
in a motion prediction task. Metrics are formulated as:

$ADE^{t}=\frac{\sum_{j=1}^{N}||\bm{x}^{t}_{j}-\mu^{t}_{j}||_{2}}{N}$ (27)
$FDE=\frac{\sum_{j=1}^{N}{||\bm{x}^{T_{pred}}_{j}-\mu^{T_{pred}}_{j}||_{2}}}{N}$
(28)

where N - number of processed pedestrians, $\bm{x}_{j}^{t}$ \- ground truth
position of $\textit{j}^{th}$ pedestrian at timestamp t, $\textit{T}_{pred}$
\- prediction horizon, $\mu$ \- predicted mean position.

The results in Table I shows that modern deep approaches are superior to
model-based approaches in terms of ADE and FDE. The Trajectron’s++[4] approach
takes into consideration a multimodal distribution, which results in the best
Euclidean errors. Still, the contributions of this paper are on the covariance
prediction, and combination of SFM and Goal Predictor could be substituted by
any modern prediction method.

###  IV-C Covariance Estimation Evaluations

We propose to evaluate the accuracy of the covariance prediction methods by
the following metric, the Part of Prediction Errors Inside:

$\textit{PPEI}_{\alpha}=E\\{\bm{1}(||\bm{x}-\mu||_{\Sigma}<\alpha)\\},\quad\alpha=\\{1,3\\}.$
(29)

where $\bm{1}()$ is the indicator function and we simply average the number of
times that our prediction lies inside the $1\sigma$ and $3\sigma$ ellipsoids
from the ground truth position by using the Mahalanobis distance.

Figure 5: Mahalanobis error (in red) distances and ADE(in blue) for our proposed method. Mahalanobis solid line is median values, and colored intervals are .25 and .75 percentiles. TABLE I: Comparison of our method against previously published methods on the ETH dataset [1]. Both ADE and FDE are reported in meters and presented here for reference(*is not valid since ground truth last pedestrian pose is used as input). Our main contribution is accurately predicted uncertainties, which can be measured as part of errors inside $1\sigma$ $3\sigma$ interval and its deviation from the ideal 2D Gaussian. Method | ADE | FDE |  % errors inside $1\sigma$ ($\Delta$ from expected) |  % errors inside $3\sigma$ ($\Delta$ from expected) | median Mahalanobis Distance
---|---|---|---|---|---
Social-LSTM[22] | 1.09 | 2.35 | 4.7$\pm$11(35) | 22 $\pm$11(76.7) | 3.4$\pm$5.66
Trajectron++ (Distribution) [4] | 0.71 | 1.66 |  50$\pm$32(10.3) | 56$\pm$30(42.4) | –
SFM + FP | 0.98* | 0.2* | 8.2$\pm$15(31.5) | 11.6$\pm$23(87.2) | –
CovarianceNet | 1.39 | 2.18 | 51.2$\pm$0.3(11.5) | 93.7$\pm$0.01(5) | 1.21$\pm$2.16
Figure 6: PPEI Evaluation of CovarianceNet method and
Trajectron++(Distribution configuration) qualities of predicted uncertainties
peaking Trajectron++ Gaussian with the highest probability (out of GMM 25
Gaussians component).

Figure 4 shows the $\textit{PPEI}_{1}$ and $\textit{PPEI}_{3}$ for each of the
methods described above. Our proposed method achieves a superior performance
(Fig 4) compared to the previously published methods in terms of performance
of predicted uncertainties.

It can be seen in Fig. 4 and in Table I that modern deep approaches output
accurate ADE, FDE but become overconfident when predicting covariances of
distributions, while proposed method with sequentially trained Goal Predictor
and CovarianceNet models are capable to predict statistically correct
covariances. Despite the fact that Trajectron’s++ average $PPEI_{1}$ over all
horizon time is close to ideal, it clearly can be observed a high deviation,
ranging from a high $PPEI_{1}$ (under confident) at initial horizon times to
low values at the final horizon time (over confident), which is shown with
high variance, $PPEI_{1}$=$50$ $\pm 32$ at Table I, while our proposed method
produces more consistent predictions $PPEI_{1}$=$51.2$ $\pm 0.3$.

The forward propagation (FP) method collapses and provides poor results due to
vanishing gradients over multiple iterations. Trajectron shows a decreasing
value of PPEI, becoming clearly overconfident. Social-LSTM also provides a
overconfident covariance estimations, with mean $1\sigma$ PPEI=4.7% and
$3\sigma$ PPEI=22%.

In Fig. 5 we show the median, 25 and 75 percentiles of the Mahalanobis
distances for the predicted distributions of CovarianceNet and the Average
Displacement Error(ADE) over prediction horizon time for predicted mean
trajectories, with dashed red line we show median Mahalanobis Distance (MD)
for ideal bi-variate Gaussian. Our projected median MD for the entire
prediction horizon is 1.21, which is close to the bi-variate Normal
distribution, with a median MD of 1.17.

We use Trajectron++[4] in Distribution configuration to estimate probability
distribution. We used only one Gaussian distribution from GMM for the metric
at each prediction step, which had the highest probability. An example of such
a prediction and that interpretations can be seen in Fig. 1. This evaluation
interpretations have similar to previous interpretations results, with
overconfident distribution for large prediction horizons and under confident
for smaller time horizons.

Fig 6 shows that individual Gaussians, composing GMM, is far from ideal, which
potentially can lead to unexpected results when a person is outside the
probabilistic estimate of the pedestrian position.

##  V Conclusions

We have presented a new hybrid method, CovarianceNet, that combines model-
based Human Motion Prediction with a neural network approach for covariance
prediction. Our approach brings an efficient calculation of motion prediction,
by using the SFM recursively and calculates the uncertainty associated with
this prediction by using a conditional deep generative model CVAE with
Gaussian latent variables.

We have evaluated our results in the ETH dataset and compared with state-of-
the art approaches. While Social-LSTM and Tajectron methods show better
accuracy when evaluating the prediction error, they are overconfident on their
predicted distributions, according to our proposed metric to measure
uncertainty prediction. Our method, CovarianceNet, is able to predict
correctly its uncertainty, thereby, our approach is better suited for
applications requiring accurate prediction of distributions.

## References

  * [1] S. Pellegrini, A. Ess, K. Schindler, and L. Van Gool, “You’ll never walk alone: Modeling social behavior for multi-target tracking,” in 2009 IEEE 12th International Conference on Computer Vision, pp. 261–268, IEEE, 2009\.
  * [2] L. Leal-Taixé, M. Fenzi, A. Kuznetsova, B. Rosenhahn, and S. Savarese, “Learning an image-based motion context for multiple people tracking,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 3542–3549, 2014.
  * [3] H. Caesar, V. Bankiti, A. H. Lang, S. Vora, V. E. Liong, Q. Xu, A. Krishnan, Y. Pan, G. Baldan, and O. Beijbom, “nuscenes: A multimodal dataset for autonomous driving,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 11621–11631, 2020.
  * [4] T. Salzmann, B. Ivanovic, P. Chakravarty, and M. Pavone, “Trajectron++: Multi-agent generative trajectory forecasting with heterogeneous data for control,” arXiv preprint arXiv:2001.03093, 2020.
  * [5] D. Mehta, G. Ferrer, and E. Olson, “Backprop-mpdm: Faster risk-aware policy evaluation through efficient gradient optimization,” in 2018 IEEE International Conference on Robotics and Automation (ICRA), pp. 1740–1746, IEEE, 2018.
  * [6] D. Helbing and P. Molnár, “Social force model for pedestrian dynamics,” Physical review E, vol. 51, no. 5, p. 4282, 1995.
  * [7] F. Giuliari, I. Hasan, M. Cristani, and F. Galasso, “Transformer networks for trajectory forecasting,” arXiv preprint arXiv:2003.08111, 2020.
  * [8] K. Sohn, H. Lee, and X. Yan, “Learning structured output representation using deep conditional generative models,” in Advances in neural information processing systems, pp. 3483–3491, 2015.
  * [9] E. A. Sisbot, L. F. Marin-Urias, R. Alami, and T. Simeon, “A human aware mobile robot motion planner,” IEEE Transactions on Robotics, vol. 23, no. 5, pp. 874–883, 2007.
  * [10] G. Ferrer, A. G. Zulueta, F. H. Cotarelo, and A. Sanfeliu, “Robot social-aware navigation framework to accompany people walking side-by-side,” Autonomous robots, vol. 41, no. 4, pp. 775–793, 2017.
  * [11] B. D. Ziebart, N. Ratliff, G. Gallagher, C. Mertz, K. Peterson, J. A. Bagnell, M. Hebert, A. K. Dey, and S. Srinivasa, “Planning-based prediction for pedestrians,” in IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 3931–3936, IEEE, 2009.
  * [12] C. Fulgenzi, A. Spalanzani, and C. Laugier, “Probabilistic motion planning among moving obstacles following typical motion patterns,” in IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 4027–4033, IEEE, 2009.
  * [13] M. Kuderer, H. Kretzschmar, C. Sprunk, and W. Burgard, “Feature-based prediction of trajectories for socially compliant navigation.,” in Robotics: science and systems, 2012.
  * [14] P. Trautman, J. Ma, R. M. Murray, and A. Krause, “Robot navigation in dense human crowds: Statistical models and experimental studies of human–robot cooperation,” The International Journal of Robotics Research, vol. 34, no. 3, pp. 335–356, 2015.
  * [15] G. Ferrer and A. Sanfeliu, “Proactive kinodynamic planning using the extended social force model and human motion prediction in urban environments,” in 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 1730–1735, IEEE, 2014.
  * [16] G. Ferrer and A. Sanfeliu, “Anticipative kinodynamic planning: multi-objective robot navigation in urban and dynamic environments,” Autonomous Robots, vol. 43, no. 6, pp. 1473–1488, 2019.
  * [17] D. Sadigh, S. Sastry, S. A. Seshia, and A. D. Dragan, “Planning for autonomous cars that leverage effects on human actions.,” in Robotics: Science and Systems, vol. 2, Ann Arbor, MI, USA, 2016.
  * [18] F. Zanlungo, T. Ikeda, and T. Kanda, “Social force model with explicit collision prediction,” EPL (Europhysics Letters), vol. 93, no. 6, p. 68005, 2011.
  * [19] F. Farina, D. Fontanelli, A. Garulli, A. Giannitrapani, and D. Prattichizzo, “Walking ahead: The headed social force model,” PloS one, vol. 12, no. 1, p. e0169734, 2017.
  * [20] J. Van Den Berg, S. J. Guy, M. Lin, and D. Manocha, “Reciprocal n-body collision avoidance,” in Robotics research, pp. 3–19, Springer, 2011.
  * [21] S. J. Guy, J. Chhugani, C. Kim, N. Satish, M. Lin, D. Manocha, and P. Dubey, “Clearpath: highly parallel collision avoidance for multi-agent simulation,” in Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp. 177–187, 2009.
  * [22] A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, L. Fei-Fei, and S. Savarese, “Social LSTM: Human trajectory prediction in crowded spaces,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 961–971, 2016.
  * [23] K. Messaoud, I. Yahiaoui, A. Verroust-Blondet, and F. Nashashibi, “Non-local social pooling for vehicle trajectory prediction,” in 2019 IEEE Intelligent Vehicles Symposium (IV), pp. 975–980, IEEE, 2019.
  * [24] A. Mohamed, K. Qian, M. Elhoseiny, and C. Claudel, “Social-stgcnn: A social spatio-temporal graph convolutional neural network for human trajectory prediction,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 14424–14432, 2020.
  * [25] S. Yi, H. Li, and X. Wang, “Pedestrian behavior understanding and prediction with deep neural networks,” in European Conference on Computer Vision, pp. 263–279, Springer, 2016.
  * [26] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [27] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in Advances in neural information processing systems, pp. 2672–2680, 2014.
  * [28] B. Ivanovic and M. Pavone, “The trajectron: Probabilistic multi-agent trajectory modeling with dynamic spatiotemporal graphs,” in Proceedings of the IEEE International Conference on Computer Vision, pp. 2375–2384, 2019\.
  * [29] A. Jain, S. Casas, R. Liao, Y. Xiong, S. Feng, S. Segal, and R. Urtasun, “Discrete residual flow for probabilistic pedestrian behavior prediction,” in Conference on Robot Learning, pp. 407–419, PMLR, 2020.
  * [30] K. Mangalam, H. Girase, S. Agarwal, K.-H. Lee, E. Adeli, J. Malik, and A. Gaidon, “It is not the journey but the destination: Endpoint conditioned trajectory prediction,” arXiv preprint arXiv:2004.02025, 2020.
  * [31] A. Postnikov, A. Gamayunov, and G. Ferrer, “HSFM-SigmaNN: Combining a feedforward motion prediction network and covariance prediction,” arXiv preprint arXiv:2009.04299, 2020.
  * [32] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in international conference on machine learning, pp. 1050–1059, 2016.
  * [33] C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger, “On calibration of modern neural networks,” in International Conference on Machine Learning, pp. 1321–1330, PMLR, 2017.
  * [34] G. Ferrer and A. Sanfeliu, “Bayesian human motion intentionality prediction in urban environments,” Pattern Recognition Letters, vol. 44, pp. 134–140, 2014.
  * [35] D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by jointly learning to align and translate,” arXiv preprint arXiv:1409.0473, 2014\.


# Transition rates of non-interacting quantum particles from the Widom
insertion formula

M. Hoyuelos<EMAIL_ADDRESS>Instituto de Investigaciones Físicas de Mar
del Plata (IFIMAR-CONICET), Departamento de Física, Facultad de Ciencias
Exactas y Naturales, Universidad Nacional de Mar del Plata, Funes 3350, 7600
Mar del Plata, Argentina

###### Abstract

Transition rates among different states in a system of quantum particles in
contact with a heat reservoir are generally obtained using heuristic
arguments. If $\bar{n}_{i}$ is the average occupation number of the final
state, the factor $1\mp\bar{n}_{i}$ appears, with minus sign for fermions and
plus sign for bosons. It is shown that this factor can be related to the
difference of the chemical potential from that of an ideal classical mixture;
this difference is formally equivalent to the excess chemical potential in a
classical system of interacting particles. Using this analogy, Widom’s
insertion formula is used in the calculation of transition rates. The result
allows an alternative derivation of quantum statistics from the condition that
transition rates depend only on the number of particles in the target energy
level. Instead, if transition rates depend on the particle number only in the
origin level, the statistics of ewkons is obtained; this is an exotic
statistics that can be applied to the description of dark energy and that is
associated to sub-distinguishable particles.

##  I Introduction

Rate equations, that govern the time evolution of average occupation numbers
$\bar{n}_{i}$ of quantum states with energy $\epsilon_{i}$ of a system of
particles in contact with a reservoir at temperature $T$, have important
applications in the description of, for example, Fermi gas, phonons or Bose-
Einstein condensation banyai2 , see also brenig ; banyai . They have the form
of the master equation,

$\frac{\partial\bar{n}_{i}}{\partial
t}=\sum_{j}(\bar{n}_{j}W_{\bar{n}_{j},\bar{n}_{i}}-\bar{n}_{i}W_{\bar{n}_{i},\bar{n}_{j}})$
(1)

with transition rates given by

$W_{\bar{n}_{j},\bar{n}_{i}}=\nu
e^{-\beta(\epsilon_{i}-\epsilon_{j})}(1\mp\bar{n}_{i})$ (2)

where minus sign is for fermions and plus sign for bosons, and $\nu$ is a jump
frequency; $W_{\bar{n}_{j},\bar{n}_{i}}$ is the transition rate from state $j$
to state $i$; in both cases $W_{\bar{n}_{j},\bar{n}_{i}}$ depends on the
occupation number in the destination state, $\bar{n}_{i}$, but does not depend
on $\bar{n}_{j}$. The factor $1\mp\bar{n}_{i}$ is generally justified using
heuristic arguments, see (brenig, , Ch. 32) or (banyai, , Sec. 8.2 and 21.2).
These rates reproduce the correct average occupation numbers in equilibrium,
that is, the Fermi-Dirac and Bose-Einstein distributions.

Transition rates proportional to $1\mp\bar{n}_{i}$ are used in kania  to
obtain a Fokker-Planck equation for the occupation number in a continuous
energy space, see also kania0 ; kania2  and (kadanoff, , Sec. 6.8). Fokker-
Planck equations for fermions and bosons can also be derived from the free
energy functional, see frank0  and (frank, , Sec. 6.5.4); in particular, the
so-called drift form of the Fokker-Planck equation involves a mobility that
behaves as $1\mp\bar{n}_{i}$. A Langevin approach has been developed in frank1
. Fundamental applications are the metal electron and black body radiation
models (frank, , Sec. 6.5.6). It can be considered that $\bar{n}_{i}$
represents the number of particles in a degenerate energy level $i$ that
includes $g_{i}$ different states; in this case, the mentioned factors for
fermions and bosons are $1\mp\bar{n}_{i}/g_{i}$.

In this paper, a method is proposed for the derivation of transition rates for
quantum systems of non-interacting particles. It is based on a previous work
for classical systems with interacting particles dimuro . The main ingredients
are the detailed balance relationship and the Widom insertion formula. From
detailed balance, transition rates can be related with the insertion energy,
that is, the energy needed to add one particle. And from the Widom insertion
formula widom , the insertion energy is related to the excess chemical
potential. A separation of terms with different orders of an extensive
quantity is used to finally obtain the transition rate in terms of the excess
chemical potential, the energy difference and an undetermined jump frequency.
In the case of non-interacting quantum particles that is analyzed here, the
excess chemical potential, $\mu_{\text{q}}$, represents quantum effects
instead of interactions. In both cases (classical with interactions or quantum
without interactions), it gives the departure of the chemical potential from
the classical ideal case.

The generalization to the quantum case of the derivation of transition rates
$W_{\bar{n}_{1},\bar{n}_{2}}$ between generic levels 1 and 2, developed in
Ref. dimuro , is presented in Sec. II. When the excess chemical potential for
fermions or bosons is used in the result for the transition rate, Eq. (2) is
obtained. In Sec. III it is shown that, before assuming a specific form for
$\mu_{\text{q}}$, if the transition rate depends only on the particle number
in the destination level, that is, if $W_{\bar{n}_{1},\bar{n}_{2}}$ depends
only on $\bar{n}_{2}$, then the known statistics of Bose-Einstein, Fermi-Dirac
and Maxwell-Boltzmann are obtained. On the other hand, Sec. IV shows that if
$W_{\bar{n}_{1},\bar{n}_{2}}$ depends only on the occupation number in the
origin level, $\bar{n}_{1}$, then the exotic statistics of ewkons is obtained;
ewkons have been introduced in hoyuelos-sisterna  to describe dark energy.
Summary and conclusions are presented in Sec. V

##  II Theory

Let us consider the quantum formulation of non-interacting particles in a
volume $V$, in equilibrium with a reservoir at temperature $T$ and chemical
potential $\mu$. There is an undetermined number of energy levels identified
with index $i$. Each level, with energy $\epsilon_{i}$, has $n_{i}$ particles.
Some basic concepts are reviewed in the next paragraphs to make the notation
clear.

###  II.1 Basic formulae

The canonical partition function for $N$ particles is

$\mathcal{Z}_{N}=\text{tr
}e^{-\beta\hat{H}}={\sum_{\\{n_{i}\\}}}^{\prime}\prod_{i}\mathcal{Z}_{n_{i}}$
(3)

where $\beta=(k_{B}T)^{-1}$ and $\hat{H}$ is the free Hamiltonian; see, for
example, (greiner, , Chap. 12) or (pathria, , Chap. 6). The sum
$\sum_{\\{n_{i}\\}}^{\prime}$ is over all $n_{i}$ that satisfy
$\sum_{i}n_{i}=N$, and $\mathcal{Z}_{n_{i}}$ is the canonical partition
function for level $i$. For fermions or bosons, we know that
$\mathcal{Z}_{n_{i}}$ is equal to $e^{-\beta\epsilon_{i}n_{i}}$ and the set of
possible values of $\\{n_{i}\\}$ has to be determined for each case. For
classical particles,

$\mathcal{Z}_{n_{i}}^{\text{cl}}=e^{-\beta\epsilon_{i}n_{i}}/n_{i}!,$ (4)

and the canonical partition function of the whole system is
$\mathcal{Z}_{N}^{\text{cl}}=(V/\lambda^{3})^{N}/N!$, where $\lambda$ is the
thermal de Broglie wavelength (pathria, , p. 147).

The grand partition function is

$\mathcal{Q}(\mu,T,V)=\text{tr
}e^{\beta(\mu\hat{N}-\hat{H})}=\sum_{N=0}^{\infty}\mathcal{Z}_{N}e^{\beta\mu
N}=\prod_{i}\mathcal{Q}_{i}$ (5)

where $\hat{N}$ is the total number of particles operator and

$\mathcal{Q}_{i}=\sum_{n_{i}}e^{\beta\mu n_{i}}\mathcal{Z}_{n_{i}}$ (6)

is the grand partition function for level $i$, where the sum is on all allowed
values of $n_{i}$ (0 to $\infty$ for bosons or classical particles; 0 or 1 for
fermions). Each element of the sum in Eq. (6) is proportional to the
probability $P_{n_{i}}$ of having $n_{i}$ particles:

$P_{n_{i}}=\frac{e^{\beta\mu n_{i}}\mathcal{Z}_{n_{i}}}{\mathcal{Q}_{i}}.$ (7)

This probability is used below when introducing detailed balance and
transition rates.

###  II.2 Average number of particles and the Widom insertion formula

If volume $V$ is large, single particle energy levels are very close to each
other. It is considered that one level with energy $\epsilon_{i}$ actually
encompasses a group of $g_{i}$ levels that are very close; number $g_{i}$ also
includes a possible degeneracy. Then, the average number of particles with
energy $\epsilon_{i}$ is given by

$\bar{n}_{i}=g_{i}\frac{1}{\beta}\frac{\partial\ln\mathcal{Q}_{i}}{\partial\mu}.$
(8)

Number $g_{i}$, and $\bar{n}_{i}$, are extensive quantities, proportional to
$V$; $g_{i}$ is used later as an expansion parameter.

Different statistics are treated here in a unified way. In order to do that,
the effective energy, $\phi_{n_{i}}$, defined as

$e^{-\beta\phi_{n_{i}}}=\frac{\mathcal{Z}_{n_{i}}}{\mathcal{Z}_{n_{i}}^{\text{cl}}},$
(9)

is introduced as a parameter to measure non-classical effects:
$\phi_{n_{i}}=0$ for classical particles, and different from zero if
$\mathcal{Z}_{n_{i}}\neq\mathcal{Z}_{n_{i}}^{\text{cl}}$. For fermions,
$\phi_{n_{i}}$ diverges if $n_{i}>1$.

Using Eqs. (6), (9) and (4) to calculate the average number of particles (8),
we get (subscript $i$ in $n_{i}$ can be removed from the sum index to simplify
the notation at this stage)

$\displaystyle\frac{\bar{n}_{i}}{g_{i}}$
$\displaystyle=\frac{1}{\mathcal{Q}_{i}}\sum_{n=0}^{\infty}\frac{n}{n!}e^{-\beta(\epsilon_{i}-\mu)n}e^{-\beta\phi_{n}}$
$\displaystyle=\frac{e^{-\beta(\epsilon_{i}-\mu)}}{\mathcal{Q}_{i}}\sum_{n=1}^{\infty}\frac{1}{(n-1)!}e^{-\beta[\phi_{n}+(\epsilon_{i}-\mu)(n-1)]}$
$\displaystyle=\frac{e^{-\beta(\epsilon_{i}-\mu)}}{\mathcal{Q}_{i}}\sum_{m=0}^{\infty}\frac{1}{m!}e^{-\beta(\phi_{m+1}-\phi_{m})}e^{-\beta[\phi_{m}+(\epsilon_{i}-\mu)m]}$
$\displaystyle=e^{-\beta(\epsilon_{i}-\mu)}\langle
e^{-\beta\Delta\phi_{n}}\rangle,$ (10)

with $\Delta\phi_{n}=\phi_{n+1}-\phi_{n}$. As expected, if $\phi_{n}=0$ then
classical statistics is recovered. Now, we define $\mu_{\text{q}}$ as

$e^{-\beta\mu_{\text{q}}}=\langle e^{-\beta\Delta\phi_{n}}\rangle.$ (11)

It can be interpreted as a correction to the chemical potential representing
non-classical effects; it is an intensive quantity of order
$\mathcal{O}(g_{i}^{0})$. We use the following notation: when $\mu_{\text{q}}$
is written without subscript $n$, it is evaluated at the average number of
particles $\bar{n}_{i}$; instead, $\mu_{\text{q},n}$ is evaluated at a
specific number of particles $n$. Its derivative respect to the number of
particles, $\mu_{\text{q},n}^{\prime}=\frac{\partial\mu_{\text{q},n}}{\partial
n}$, is $\mathcal{O}(g_{i}^{-1})$.

From the average particle number in each case, it is known that

$e^{-\beta\mu_{\text{q}}}=\left\\{\begin{array}[]{ll}1&\text{for classical
particles,}\\\
\frac{1}{1+e^{-\beta(\epsilon_{i}-\mu)}}=1-\bar{n}_{i}/g_{i}&\text{for
fermions,}\\\
\frac{1}{1-e^{-\beta(\epsilon_{i}-\mu)}}=1+\bar{n}_{i}/g_{i}&\text{for
bosons.}\end{array}\right.$ (12)

We will use $\mu_{\text{q}}$ in order to treat all possible statistics in a
unified way, as mentioned before. Eq. (11) is formally equivalent to the Widom
insertion formula widom  (see also (hansen2, , p. 30)) since $\mu_{\text{q}}$
is interpreted as the excess chemical potential and $\Delta\phi_{n}$ as the
insertion energy (the interaction energy needed to insert one particle) in a
classical system of interacting particles. In the next subsection, the
differences $\phi_{n_{2}+1}-\phi_{n_{2}}$ and $\phi_{n_{1}}-\phi_{n_{1}-1}$,
for energy levels 1 and 2, are needed. They can be written in terms of
$\mu_{\text{q},n_{i}}$ using (11) and following the same formal steps that are
described in Ref. dimuro :

$\displaystyle\phi_{n_{2}+1}-\phi_{n_{2}}$
$\displaystyle=\mu_{\text{q},n_{2}}-\frac{1}{2\beta}\frac{\Gamma^{\prime}_{n_{2}}}{\Gamma_{n_{2}}}+\frac{\mu^{\prime}_{\text{q},n_{2}}}{2}+\text{h.t.}$
(13) $\displaystyle\phi_{n_{1}}-\phi_{n_{1}-1}$
$\displaystyle=\underbrace{\vphantom{\frac{1}{\Gamma_{n_{2}}}}\mu_{\text{q},n_{1}}}_{\mathcal{O}(1)}-\underbrace{\frac{1}{2\beta}\frac{\Gamma^{\prime}_{n_{1}}}{\Gamma_{n_{1}}}-\frac{\mu^{\prime}_{\text{q},n_{1}}}{2}}_{\mathcal{O}(g_{i}^{-1})}+\;\text{h.t.}$
(14)

where $\Gamma_{n_{i}}=1+\beta n_{i}\mu_{\text{q},n_{i}}^{\prime}$, and higher
order terms of $1/g_{i}$ are represented by “h.t.”.

###  II.3 Transition rates

Let us use labels 1 and 2 for two generic energy levels in an initial state
with $n_{1}$ and $n_{2}$ particles. After a jump of one particle form level 1
to level 2, the final state has $n_{1}-1$ and $n_{2}+1$ particles. Let us call
$W_{AB}$ the transition rate from state $A=\\{n_{1},n_{2}\\}$ to state
$B=\\{n_{1}-1,n_{2}+1\\}$, and $W_{BA}$ the rate of the inverse process. The
system is in equilibrium and detailed balance is satisfied:

$P_{A}\,W_{A,B}=P_{B}\,W_{B,A},$ (15)

with $P_{A}$ and $P_{B}$ the probabilities of states $A$ and $B$. These
probabilities are $P_{A}=P_{n_{1}}P_{n_{2}}$ and
$P_{B}=P_{n_{1}-1}P_{n_{2}+1}$ with $P_{n_{i}}$ given by (7). Then,

$\mathcal{Z}_{n_{1}}\mathcal{Z}_{n_{2}}\,W_{A,B}=\mathcal{Z}_{n_{1}-1}\mathcal{Z}_{n_{2}+1}\,W_{B,A},$
(16)

From Eqs. (16) and (9) we have

$\frac{W_{A,B}}{W_{B,A}}=\frac{e^{-\beta(\phi_{n_{2}+1}-\phi_{n_{2}})}}{e^{-\beta(\phi_{n_{1}}-\phi_{n_{1}-1})}}\frac{\mathcal{Z}_{n_{1}-1}^{\text{cl}}\mathcal{Z}_{n_{2}+1}^{\text{cl}}}{\mathcal{Z}_{n_{1}}^{\text{cl}}\mathcal{Z}_{n_{2}}^{\text{cl}}}=\frac{e^{-\beta(\phi_{n_{2}+1}-\phi_{n_{2}})}}{e^{-\beta(\phi_{n_{1}}-\phi_{n_{1}-1})}}\frac{n_{1}}{n_{2}+1}e^{-\beta(\epsilon_{2}-\epsilon_{1})},$
(17)

where Eq. (4) for $\mathcal{Z}_{n_{i}}^{\text{cl}}$ was used.

The rate $W_{A,B}$ corresponds to the jump of any of the $n_{1}$ particles
from level 1 to level 2. The transition rates per particle in the origin level
are $W_{n_{1},n_{2}}=W_{A,B}/n_{1}$ and
$W_{n_{2}+1,n_{1}-1}=W_{B,A}/(n_{2}+1)$, where the order of subindices
indicates the jump direction. Then,

$W_{n_{1},n_{2}}\,e^{-\beta(\epsilon_{1}+\phi_{n_{1}}-\phi_{n_{1}-1})}=W_{n_{2}+1,n_{1}-1}\,e^{-\beta(\epsilon_{2}+\phi_{n_{2}+1}-\phi_{n_{2}})}.$
(18)

The next step is to use Eqs. (13) and (14) for $\phi_{n_{2}+1}-\phi_{n_{2}}$
and $\phi_{n_{1}}-\phi_{n_{1}-1}$ and to separate different orders. Orders
$g_{1}^{-1}$ and $g_{2}^{-1}$ are equivalent since $g_{1}$ and $g_{2}$ are
proportional to $V$. Details of these calculations are described in Ref.
dimuro  in the context of interacting classical particles, where the excess
chemical potential, $\mu_{\text{ex},n_{i}}$, has to be replaced by
$\mu_{\text{q},n_{i}}+\epsilon_{i}$. The result is

$W_{n_{1},n_{2}}=\nu
e^{-\beta(\epsilon_{2}-\epsilon_{1})/2}\frac{e^{\beta\mu_{\text{q},n_{1}}/2}}{(1+\beta
n_{1}\mu_{\text{q},n_{1}}^{\prime})^{1/2}}\frac{e^{-\beta\mu_{\text{q},n_{2}}/2}}{(1+\beta
n_{2}\mu_{\text{q},n_{2}}^{\prime})^{1/2}}.$ (19)

where $\nu$ is, in general, a function of the sum $n_{1}+n_{2}$; it represents
the strength of the coupling with the heat reservoir and determines the time
scale of the transition rate. This general approach is useful to determine the
form of transition rates associated to different statistics, that are
represented by $\mu_{\text{q},n_{i}}$, but leaves factor $\nu$ undetermined.

Using (12) for $e^{-\beta\mu_{\text{q},n_{i}}}$ in (19), we get the transition
rates for the known statistics,

$W_{n_{1},n_{2}}=\left\\{\begin{array}[]{ll}\nu
e^{-\beta(\epsilon_{2}-\epsilon_{1})/2}&\text{classical particles}\\\ \nu
e^{-\beta(\epsilon_{2}-\epsilon_{1})/2}(1-n_{2}/g_{2})&\text{fermions}\\\ \nu
e^{-\beta(\epsilon_{2}-\epsilon_{1})/2}(1+n_{2}/g_{2})&\text{bosons.}\end{array}\right.$
(20)

The result for fermions has an immediate interpretation: the transition rate
is proportional to the number of states available in the destination energy
level, $g_{2}-n_{2}$, due to the Pauli exclusion principle. As mentioned in
the introduction, this transition rates for non-interacting quantum particles
can be found in the literature, usually derived through heuristic arguments;
see for example (brenig, , Ch. 32) and (banyai, , Sec. 8.2 and 21.2).
Reproduction of the known results confirms the validity of the procedure
leading to Eq. (19).

##  III Dependence on destination level

From (20), it can be seen that transition rates depend only on the number of
particles in the destination level, $n_{2}$. If the condition that
$W_{\bar{n}_{1},\bar{n}_{2}}$ depends only on $\bar{n}_{2}$ is assumed, then
statistics of fermions, bosons, and classical particles should be obtained.
This assertion is verified as follows. Considering that Eq. (19) is evaluated
at the average particle numbers, the condition implies that the factor that
depends on $\bar{n}_{1}$ must be constant, that is

$\frac{e^{\beta\mu_{\text{q}}/2}}{(1+\beta\bar{n}\mu_{\text{q}}^{\prime})^{1/2}}=1,$
(21)

where subindex 1 is removed in $\bar{n}_{1}$ to lighten the notation. The
constant on the right side is set equal to 1 taking into account the case of
classical particles, for which $\mu_{\text{q}}=0$. Solutions of this equation
are given by

$e^{-\beta\mu_{\text{q}}}=1+\kappa\bar{n},$ (22)

where $\kappa$ is a constant. Comparing with (12), $\kappa=0$ corresponds to
classical particles, $\kappa=-1/g$ to fermions and $\kappa=1/g$ to bosons
(subindex 1 is removed in $g_{1}$ for simplicity).

##  IV Exotic statistics

The previous analysis, in which it has been shown that the assumption of
transition rates depending only on $\bar{n}_{2}$ is enough to derive known
statistics, suggests the possibility of exotic statistics that may be deduced
from the theory when time is reversed. The reasoning is as follows. According
to the Feynman-Stueckelberg interpretation, antiparticles are viewed as
negative energy modes of the quantum field that propagate backward in time. An
hypothetical particle (not necessarily an anti-fermion or an anti-boson) that
propagates backward in time would have transition rates that, instead of
depending on the number of particles in the destination level, depend on the
number in the origin level, $\bar{n}_{1}$. The term that depends on $n_{2}$ in
(19) must be constant. The possibility of time reversed particles is mentioned
here as a motivation to consider transition rates in which the roles of origin
and destination levels are interchanged. The conjecture that there are
particles that meet this condition is proposed without assuming that they are
particles or antiparticles. If that is the case, statistics should satisfy the
condition

$\frac{e^{-\beta\mu_{\text{q}}/2}}{(1+\beta\bar{n}\mu_{\text{q}}^{\prime})^{1/2}}=1,$
(23)

where subindex 2 in $\bar{n}_{2}$ was removed for simplicity. The only
difference with Eq. (21) is that the exponential has the opposite sign.
Solutions are

$e^{-\beta\mu_{\text{q}}}=\frac{1}{1-\kappa/\bar{n}},$ (24)

with $\kappa$ a constant. Using this solution for $e^{-\beta\mu_{\text{q}}}$
in the expression for the mean number of particles, Eq. (10), we obtain

$\bar{n}/g=e^{-\beta(\epsilon-\mu)}+\kappa/g.$ (25)

As before, $\kappa=0$ corresponds to classical particles. The other case that
is considered here is $\kappa=g$, that corresponds to the so-called statistics
of ewkons hoyuelos-sisterna ; hoyuelos1 ; hoyuelos2 . Calling
$\bar{n}_{\text{ewk}}=\bar{n}/g$ the average number of ewkons in one of the
$g$ states, then

$\bar{n}_{\text{ewk}}=e^{-\beta(\epsilon-\mu)}+1.$ (26)

What are the features of particles that produce ewkon statistics? This
question is addressed in the next subsection.

###  IV.1 Features of ewkons

What happens when particles are identical but not fully distinguishable or
indistinguishable? In this subsection, it is shown that an intermediate
category of particles, that can be called sub-distinguishable, have ewkon
statistics. In order to do that, the basic framework is presented first.

There are $N$ particles in $L$ states; each state has energy $\epsilon_{i}$
with $i=1\cdots L$. The Hamiltonian of non-interacting particles is
$\hat{H}=\sum_{i=1}^{L}\hat{n}_{i}\epsilon_{i}$ where $\hat{n}_{i}$ is the
particle number operator of state $i$. Degeneracy is taken into account
implicitly since the values of $\epsilon_{i}$ may be the same for different
$i$ (index $i$ represents a set of quantum numbers that determines the state).
Note that this definition of number $n_{i}$ is different from the one in the
previous sections; it does not encompass the number of particles in degenerate
states or in close energy levels when the volume is large, it is an intensive
quantity whose mean value is the average number of particles of the previous
sections divided by the corresponding $g_{i}$.

Using the eigenstates of the number operators, the system’s state is given by

$|n_{1}\cdots n_{L}\rangle=|n_{1}\rangle\cdots|n_{L}\rangle.$ (27)

As mentioned in Sec. II.1, the canonical partition function is

$\mathcal{Z}_{N}={\sum_{\\{n_{i}\\}}}^{\prime}\langle n_{1}\cdots
n_{L}|e^{-\beta\hat{H}}|n_{1}\cdots
n_{L}\rangle={\sum_{\\{n_{i}\\}}}^{\prime}\prod_{i}\mathcal{Z}_{n_{i}}.$ (28)

Including a super-index $\mp$ to indicate symmetric or anti-symmetric states,
the partition function of state $i$ with $n_{i}$ particles is

$\mathcal{Z}_{n}^{\mp}={}^{\mp}\langle
n|e^{-\beta\hat{n}\epsilon}|n\rangle^{\mp}=e^{-\beta n\epsilon}w^{\mp}$ (29)

where $w^{\mp}$ is the statistical weight factor; $w^{+}=1$ for bosons and
$w^{-}=\delta_{n,0}+\delta_{n,1}$ for fermions; sub-index $i$ was omitted for
simplicity. Let us call $|i_{j}\rangle$ the normalized state of particle $j$.
The $n$ non-interacting particles are described by the number sate $|n\rangle$
that, in turn, is given by the symmetrization or anti-symmetrization of one-
particle product states:

$|n\rangle^{\mp}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{P}(\mp 1)^{p}P|i_{1}\cdots
i_{n}\rangle,$ (30)

where the sum is over possible permutations $P$, $\mathcal{N}$ is the
normalization factor and $p$ is the parity of permutation $P$, with
$P|i_{1}\cdots i_{n}\rangle=|i_{P1}\cdots i_{Pn}\rangle$. Since the $n$
particles are in the same state, all permutations are equal. The normalization
is $\mathcal{N}=n!^{2}$ for bosons or fermions (in the last case the
normalization is 1 since $n$ takes values 0 or 1).

Instead, for classical distinguishable particles, the statistical wight factor
is

$w_{\text{dist}}=1/n!,$ (31)

an expression that leads to the Maxwell-Boltzmann distribution. The usual
justification that can be found in textbooks is that factor $1/n!$ is obtained
in the dilute limit of quantum statistical mechanics. There are other possible
justifications that include the assumption that identical classical particles
should be treated as permutable; more extensive discussions on this
fundamental subject can be found in, for example, saunders ; dieks ; darrigol
; saunders2 . In any case, the statistical weight factor is 1 over the number
of distinguishable configurations, $n!$.

There are two extreme situations: all particles indistinguishable with one
possible configuration, or all particles distinguishable with $n!$ possible
configurations. In the first case, $|i_{P1}\cdots i_{Pn}\rangle=|i_{1}\cdots
i_{n}\rangle$ holds for any $P$ because $i_{1}=\cdots=i_{n}$. In the second
case, the equality $\\{i_{P1}\cdots i_{Pn}\\}=\\{i_{1}\cdots i_{n}\\}$ holds
only if $P$ is the identity. (Ket notation, $|i_{1}\cdots i_{n}\rangle$, is
used here only for indistinguishable particles; in other cases, curly brackets
are used.)

Figure 1: System of four particles arranged on a ring. A cyclic rotation
generates an indistinguishable configuration (right) due to rotation symmetry.
In cycle notation, operator $C$ is represented by $(1,4,3,2)$, meaning that
element 1 goes to 4, 4 to 3, 3 to 2 and 2 to 1.

The purpose is to analyze the possibility of an intermediate situation.
Consider a system of $n$ particles that are arranged along a ring with
rotation symmetry. The arrangement order can be distinguished, but cyclic
rotations produce indistinguishable configurations, see Fig. 1. The number of
distinguishable configurations is $(n-1)!$; since this number is reduced from
$n!$ to $(n-1)!$, the term “sub-distinguishable particles” is proposed. The
resulting statistical weight factor,

$w_{\text{ewk}}=1/(n-1)!,$ (32)

corresponds to statistics of ewkons. The same factor is obtained not only for
cyclic rotations but for any cyclic permutation $C$ that satisfies $C^{n}=I$,
where $I$ is the identity (gross, , p. 29). More generally, sub-
distinguishable particles satisfy, by definition, the following condition: if
$n>0$ and if $P$ is equal to a cyclic permutation $C$ applied $j$ times, that
is $P=C^{j}$ (with $j=1,\cdots,n$), then $\\{i_{P1}\cdots
i_{Pn}\\}=\\{i_{1}\cdots i_{n}\\}$. In other words, configurations that differ
by cyclic permutations are indistinguishable. The partition function,

$\mathcal{Z}_{n,\text{ewk}}=\frac{e^{-\beta\epsilon n}}{(n-1)!}$ (33)

generates the statistics of ewkons for identical sub-distinguishable
particles. Since an ideal gas of ewkons has a negative relation between
pressure and energy density, its main application is to dark energy hoyuelos-
sisterna ; hoyuelos1 ; hoyuelos2 .

The (one-level) grand partition function for ewkons is

$\displaystyle\mathcal{Q}_{\text{ewk}}$ $\displaystyle=\sum_{n}e^{\beta\mu
n}\mathcal{Z}_{n}=\sum_{n=1}^{\infty}\frac{e^{-\beta(\epsilon-\mu)n}}{(n-1)!}=\sum_{m=0}^{\infty}\frac{e^{-\beta(\epsilon-\mu)(m+1)}}{m!}$
$\displaystyle=\exp[-\beta(\epsilon-\mu)+e^{-\beta(\epsilon-\mu)}]$ (34)

and the average number of particles is given by Eq. (26),
$\bar{n}_{\text{ewk}}=e^{-\beta(\epsilon-\mu)}+1$.

From Eqs. (10) (with $\bar{n}_{\text{ewk}}=\bar{n}/g$) and (11), the excess
chemical potential is given by
$e^{-\beta\mu_{\text{q}}}=\bar{n}_{\text{ewk}}e^{\beta(\epsilon-\mu)}$, so

$e^{-\beta\mu_{\text{q}}}=\frac{1}{1-1/\bar{n}_{\text{ewk}}}.$ (35)

with the condition $\bar{n}_{\text{ewk}}>1$.

Using (35) in Eq. (19), the transition rate from state 1 to state 2, with
$\bar{n}_{1}$ and $\bar{n}_{2}$ average particle numbers (subscript ‘ewk’ is
removed) is

$W_{\bar{n}_{1},\bar{n}_{2}}=\nu
e^{-\beta(\epsilon_{2}-\epsilon_{1})/2}(1-1/\bar{n}_{1}).$ (36)

It depends only on the particle number in the origin level. If there is only
one particle, it can not leave its state since the transition rate is zero (it
is an absorbing state). Starting from an initial configuration in which
$n_{i}\geq 1\;\forall i$, this condition remains valid for any subsequent
time. If the initial configuration has empty states, it irreversibly evolves
to the equilibrium configuration in which all states are occupied with at
least one particle, with additional particles provided by the reservoir, that
keeps the chemical potential $\mu$. The problem of a divergent total particle
number (or total energy) is solved by including an upper limit for the energy,
$\epsilon_{\text{max}}$ (hoyuelos2, , Sec. 5). The limit can be assumed as a
feature of the density of states or can be a consequence of a limitation of
the reservoir to provide additional particles (without the limit, the number
of particles that the reservoir has to transfer to the system diverges).

The procedure includes the possibility of super-distinguishable particles,
called genkons in hoyuelos-sisterna ; $n$ particles of this kind have $(n+1)!$
distinguishable configurations. This case corresponds to $\kappa=-g$ in (25),
resulting an average number of particles
$\bar{n}_{\text{gen}}=e^{-\beta(\epsilon-\mu)}-1$. An interesting symmetry of
ewkons and genkons is that they have average particle numbers equal to the
inverse of that of fermions and bosons, respectively, if the sign of
$\epsilon-\mu$ is inverted. Nevertheless, super-distinguishable particles are
not further discussed because, since $n$ can be equal to $-1$, they lack, so
far, a satisfactory physical interpretation.

###  IV.2 Toy model

Simplified models are useful for a better understanding of different particle
features. For example, particles with hard-core interaction in a one-
dimensional lattice, with jump rates between neighboring sites that include an
external force, reproduce the Fermi-Dirac statistics kania ; suarez . Hard-
core interaction plays the role of the Pauli exclusion principle. An
equivalent model for ewkons is proposed in the next paragraphs.

A one-dimensional lattice of size $L$ represents a discrete energy space; the
energy of site $i$ is $\epsilon_{i}=i\Delta\epsilon$, with constant energy
gaps, $\Delta\epsilon=\epsilon_{i+1}-\epsilon_{i}\geq 0$, between any two
neighboring sites. At $t=0$, $N$ particles are randomly distributed in the $L$
sites. A particle in site $i$ jumps to the right or to the left with the rates
corresponding to ewkons:

$\displaystyle W_{\text{r}}$ $\displaystyle=\nu
e^{-\beta\Delta\epsilon/2}(1-1/n_{i})$ $\displaystyle W_{\text{l}}$
$\displaystyle=\nu e^{\beta\Delta\epsilon/2}(1-1/n_{i})$ (37)

Only jumps to nearest neighbors are allowed. In a small time interval $\Delta
t=1/(2N\nu e^{\beta\Delta\epsilon/2})$, the jump probabilities are

$\displaystyle P_{\text{r}}$
$\displaystyle=e^{-\beta\Delta\epsilon}(1-1/n_{i})/(2N)$ $\displaystyle
P_{\text{l}}$ $\displaystyle=(1-1/n_{i})/(2N).$ (38)

In a Monte Carlo (MC) simulation, a particle is randomly chosen with
probability $1/N$ and then jumps to the right or to the left with
probabilities $e^{-\beta\Delta\epsilon}(1-1/n_{i})/2$ and $(1-1/n_{i})/2$
respectively; a MC step is a time interval, equal to $N\Delta t$, in which
each particle has, on average, one chance to jump. Like fermions, ewkons have
an effective repulsive interaction. A particle can jump to another state only
if there is another particle with it. The more particles are in the same state
(or site), more probable it is to jump. Figure 2 shows numerical results of
$\bar{n}_{i}$ against position $i$, with zero current condition at the borders
(see the figure caption for the simulation parameters). After a transient, the
average number of particles converges to Eq. (26):
$\bar{n}_{\text{ewk}}=e^{-\beta(\epsilon-\mu)}+1$, represented by the blue
curve. The value of $e^{-\beta\mu}$ is obtained from the condition
$N=\sum_{i=1}^{L}\bar{n}_{i}$.

If $N<L$, the equilibrium distribution is not reached, and the final state is
a frozen configuration, because there is 0 or 1 particle in each site; the
final configuration depends on the initial condition.

Figure 2: Average number of particles, $\bar{n}_{i}$, against lattice position
$i$ with transition rates for ewkons, Eq. (IV.2). The initial condition is
random. After a time of $15000$ MC steps the equilibrium distribution is
reached. Crosses are averages taken during $5000$ MC steps. Parameters:
$N=200$, $L=100$, $\beta\Delta\epsilon=0.1$. The blue curve is the equilibrium
distribution for ewkons, Eq. (26).

##  V Conclusions

The Widom insertion formula relates the excess chemical potential with the
energy needed to insert one particle. The formula has been previously applied
to derive a general expression of transition rates in a classical system of
interacting particles dimuro . These calculations are adapted for a quantum
system of non-interacting particles at temperature $T$, for which the excess
chemical potential, $\mu_{\text{q}}$, represents quantum effects (it is zero
for the classical ideal mixture). From detailed balance, transition rates can
be related to insertion energies that, in turn, are related to the excess
chemical potential trough the Widom insertion formula. The result for the
transition rate $W_{n_{1},n_{2}}$ is given in terms of the energy difference,
$\epsilon_{2}-\epsilon_{1}$, and the excess chemical potential in both levels,
see Eq. (19). The expression evaluated at average particle numbers reproduces
known results for fermions and bosons that include a factor
$1\mp\bar{n}_{2}/g_{2}$ respectively. With this factors in the transition
rates, the master equation (or Fokker-Planck equation if the energy space is
continuous) gives an irreversible evolution towards the Fermi-Dirac or Bose-
Einstein distributions.

The expression obtained for transition rates provides an alternative procedure
to derive known statistics. If it is assumed that the transition rate depends
only on the number of particles in the destination level, then Fermi-Dirac,
Bose-Einstein and Maxwell-Boltzmann statistics are obtained. This result
prompts an exploration in a time-reversed situation, where origin and
destination levels are interchanged. If the transition rate depends only on
the number of particles in the origin level (instead of destination level),
then the statistics of ewkons is derived.

The statistical weight factor for distinguishable particles is $1/n!$, while
for ewkons it is $1/(n-1)!$. The number of distinguishable configurations for
ewkons is reduced respect to classical particles; for this reason, the term
‘sub-distinguishable particles’ is also used. A sufficient condition to obtain
ewkon statistics is to consider that particle configurations that differ by a
cyclic permutation are indistinguishable; this condition reduces the number of
distinguishable configurations to $(n-1)!$. Ewkon statistics, like Maxwell-
Boltzmann statistics, is outside the scope of the spin-statistics theorem,
since the theorem applies to indistinguishable particles.

In summary, Widom’s insertion formula provides useful information for
determining transition rates in a quantum system of non-interacting particles.
The form of the result suggest the possibility of sub-distinguishable
particles with ewkon statistics when transition rates depend only on the
number of particles in the origin level. It has been shown in previous works
hoyuelos-sisterna ; hoyuelos2  that an ideal gas of ewkons has a negative
relation between pressure and energy density, a feature that makes them
appropriate for the description of dark energy.

## Acknowledgments

This work was partially supported by Consejo Nacional de Investigaciones
Científicas y Técnicas (CONICET, Argentina, PIP 112 201501 00021 CO).

## References

  * (1) L. Bányai, P. Gartner, O. M. Schmitt, H. Haug, Condensation kinetics for bosonic excitons interacting with a thermal phonon bath, Phys. Rev. B 61 (2000) 8823–8834.
  * (2) W. Brenig, Statistical Theory of Heat, Springer, Berlin, 1989.
  * (3) L. Bányai, Lectures on Non-Equilibrium Theory of Condensed Matter, World Scientific, Singapore, 2006.
  * (4) G. Kaniadakis, P. Quarati, Classical model of bosons and fermions, Phys. Rev. E 49 (1994) 5103.
  * (5) G. Kaniadakis, P. Quarati, Kinetic equation for classical particles obeying an exclusion principle, Phys. Rev. E 48 (1993) 4263.
  * (6) G. Kaniadakis, Nonlinear kinetics underlying generalized statistics, Physica A 296 (2001) 405.
  * (7) L. P. Kadanoff, Statistical physics: statics, dynamics and renormalization, World Scientific, Singapore, 2000.
  * (8) T. D. Frank, A. Daffertshofer, Nonlinear Fokker-Planck equations whose stationary solutions make entropy-like functionals stationary, Physica A 272 (1999) 497.
  * (9) T. D. Frank, Nonlinear Fokker-Planck Equations, Springer, Berlin, 2005.
  * (10) T. D. Frank, A langevin approach for the microscopic dynamics of nonlinear fokker-planck equations, Physica A 301 (2001) 52.
  * (11) M. D. Muro, M. Hoyuelos, Application of the Widom insertion formula to transition rates in a lattice, Phys. Rev. E in press (2021).
  * (12) B. Widom, Some topics in the theory of fluids, J. Chem. Phys. 39 (1963) 2808.
  * (13) M. Hoyuelos, P. Sisterna, Quantum statistics of classical particles derived from the condition of a free diffusion coefficient, Phys. Rev. E 94 (2016) 062115\.
  * (14) W. Greiner, L. Neise, H. Stöcker, Thermodynamics and Statistical Mechanics, Springer, New York, 1997.
  * (15) R. K. Pathria, P. D. Beale, Statistical Mechanics, 3rd Edition, Elsevier, Amsterdam, 2011.
  * (16) J.-P. Hansen, I. R. McDonald, Theory of Simple Liquids: with Applications to Soft Matter, Academic Press, Oxford, 2013.
  * (17) M. Hoyuelos, From creation and annihilation operators to statistics, Physica A 490 (2018) 944.
  * (18) M. Hoyuelos, Exotic statistics for particles without spin, J. Stat. Mech.: Theory Exp. 2018 (2018) 073103.
  * (19) S. Saunders, Indistinguishability, in: R. Batterman (Ed.), The Oxford Handbook of Philosophy of Physics, Oxford University Press, 2013.
  * (20) D. Dieks, The Gibbs paradox and particle individuality, Entropy 20 (2018) 466\.
  * (21) O. Darrigol, The Gibbs paradox: Early history and solutions, Entropy 20 (2018) 443.
  * (22) S. Saunders, The concept ‘indistinguishable’, Studies in History and Philosophy of Science Part B: Studies in History and Philosophy of Modern Physics 71 (2020) 37.
  * (23) J. L. Gross, Combinatorial Methods with Computer Applications, Chapman & Hall/CRC, Boca Raton, 2008.
  * (24) G. Suárez, M. Hoyuelos, H. Mártin, Mean-field approach for diffusion of interacting particles, Phys. Rev. E 92 (2015) 062118.