# Predicting AI Agent Behavior through Approximation of the Perron-Frobenius
Operator

Shiqi Zhang
Department of Mechanical Engineering
University of California, Riverside
Riverside, CA 92521
<EMAIL_ADDRESS>&Darshan Gadginmath
Department of Mechanical Engineering
University of California, Riverside
Riverside, CA 92521
<EMAIL_ADDRESS>Fabio Pasqualetti
Department of Mechanical Engineering
University of California, Riverside
Riverside, CA 92521
<EMAIL_ADDRESS>
###### Abstract

Predicting the behavior of AI-driven agents is particularly challenging
without a preexisting model. In our paper, we address this by treating AI
agents as nonlinear dynamical systems and adopting a probabilistic perspective
to predict their statistical behavior using the Perron-Frobenius (PF)
operator. We formulate the approximation of the PF operator as an entropy
minimization problem, which can be solved by leveraging the Markovian property
of the operator and decomposing its spectrum. Our data-driven methodology
simultaneously approximates the PF operator to perform prediction of the
evolution of the agents and also predicts the terminal probability density of
AI agents, such as robotic systems and generative models. We demonstrate the
effectiveness of our prediction model through extensive experiments on
practical systems driven by AI algorithms.

##  1 Introduction

The integration of artificial intelligence (AI) models within autonomous
agents has transformed many fields, such as autonomous vehicles and
personalized recommendation systems. A wide range of complex models such as
LLMs, diffusion models, and different neural networks architectures have been
used for the above-mentioned applications. They are trained with data that has
inherent biases and can cause misaligned performance. Autonomous agents
operate in dynamic environments, making decisions based on continuous feedback
that allows them to learn and adapt over time. Therefore, studying the
behavior and alignment of these AI-driven agents is critical for several
reasons. Analyzing their actions can help prevent behaviors that conflict with
human values and ethical standards. Further, understanding their behavior is
essential for enhancing their efficiency and reliability, which is
particularly important in safety-critical applications. These intelligent
models are often complex, high-dimensional, and only partially observable over
short time intervals. This complexity raises the question of what abstract
properties can be efficiently quantified to delineate the boundaries of their
intellectual capabilities. The design of models and understanding of the
properties of AI components embedded within these agents depends crucially on
the ability to analyze the interplay between AI-driven decision-making and the
physical behavior of the closed-loop system. This capability is foundational
for users to perceive, predict, and interact effectively with intelligent
systems. It also provides the theoretical and technical basis for practical
tasks such as decision-making and reinforcement learning.

Considering these challenges, it is important to develop methods that are
capable of harnessing critical information to identify the behavior of AI
components in closed-loop. Such an algorithm needs to provide useful laws and
principles governing the behaviors of these AI-driven agents. Among the
emerging methodologies, there has been a notable increase in modeling these
behaviors as nonlinear dynamical systems. Originating from studies in partial
differential equations (PDEs) and fluid mechanics, techniques such as Dynamic
Mode Decomposition (DMD) and its generalizations have demonstrated significant
capability in revealing the underlying evolutionary laws of AI agents. An
alternative to the dynamic-mode-based method is the statistical perspective.
This approach, drawing from statistical mechanics, models the behavior of AI
agents as stochastic processes. The Perron-Frobenius (PF) operators can offer
a useful tool in analyzing the evolution of stochastic processes with
underlying nonlinear dynamics. Although the application of probabilistic
models to learn and predict the statistical behavior of complex AI agents has
increasingly attracted interest in areas such as autonomous driving, motion
planning, and human-robot interaction, algorithms based on this probabilistic
perspective are not yet to be fully explored.

###  1.1 Literature Review

Here we review some empirical and theoretical models from a statistical
modeling perspective that have been developed with the aim of improving the
prediction of the behavior of agents with AI models in the loop. In [1], the
constrained Ulam Dynamic Mode Decomposition method is presented to approximate
the Perron-Frobenius operators for both the deterministic and the stochastic
systems. [2] provides a numerical approximation of the PF operators using the
finite volume method. These early numerical methods lay very important
foundations for the subsequent development of deep-learning-based methods,
which turn out to be more scalable approaches. For example, [3] provides a
direct learning method by training a neural network-based state transition
function (operator). The works [4, 5, 6] estimate the reachability sets of
neural network-controlled systems. The works [7] and [8] use deep operator
networks (DeepONets) and Fourier neural operators (FNOs) to approximate the
solution trajectories of PDEs. More recently, [9] proposes an approximation
method based on kernel density estimation (KDE) and [10] present a deep
reproducing kernel Hilbert module (deep-RKHM), serving as a deep learning
framework for kernel methods. Beside the above-mentioned empirical methods,
there are also some works that try to learn to statistical behavior from an
optimal-transport-theory perspective. For example, in [11], the authors seek
to recover the parameters in dynamical systems with a single smoothly varying
attractor by approximating the physical invariant measure.

A formal analysis of agents with AI models in feedback from the perspective of
the evolution of the density is currently lacking in the literature.
Reachability analysis needs tracking of every possible trajectory which can be
computationally expensive. The density evolution approach uses a single
quantity that measures the probability of evolution of trajectories and offers
a significant computational advantage. Further, the density perspective allows
a convenient method to verify the alignment of machine learning models. We
seek to address these challenges in this work.

###  1.2 Our Contributions

Our main contributions are as follows:

  * •

AI-driven agents behave in unpredictable ways due to machine learning black
boxes. To the best of our knowledge, we are the first ones to look at this
through the lens of propagation of probability densities and the PF operator.
AI agents are trained with data that has inherent biases. This, coupled with
the structure of machine learning models, can potentially alter the alignment
of the model. To verify the alignment of the model, we predict the asymptotic
behavior of the model by analyzing the terminal stationary density of the AI
agents.

  * •

We propose PISA, a novel and scalable algorithm that can simultaneously
predict the evolution of the densities of AI agents and estimate their
terminal density. Our algorithm is motivated by the spectral decomposition
theorem [12, 13] and provides a theoretical backing for its performance. PISA
simultaneously approximates the action of the Perron-Frobenius operator from
the trajectory data of agents and predicts their asymptotic behavior.

  * •

In our proposed algorithm PISA, the model complexity is indexed by the number
of basis functions. The number of basis functions is a tunable parameter that
can be altered according to the user’s needs. We provide a theoretical
guarantee of the existence of the optimal solution to our operator estimation
problem.

  * •

We numerically verify the effectiveness of PISA in a variety of practical
cases and compare it with existing literature. We first predict the behavior
of unicycle robots driven by a controller based on diffusion models. Then we
analyze the behavior of generative models from the lens of density evolution.
Lastly, we apply PISA in the case of predicting the movement of pedestrians.
We observe that PISA performs significantly better than the existing
literature.

##  2 Problem Formulation

Consider an agent with state $x$ whose dynamics are defined by

$\displaystyle\dot{x}=h(x,u),$ (1)

where $u$ is an external input to the system. With a parameterized machine
learning model as feedback $u=\texttt{ML}_{\theta}(x)$, the system’s dynamics
including the feedback input is given by given by

$\dot{x}=h(x,\texttt{ML}_{\theta}(x))=f(x),$ (2)

where $x(t)\in X\subseteq\mathbb{R}^{M}$ and
$f(\cdot):\mathbb{R}^{M}\mapsto\mathbb{R}^{M}$ is a nonlinear continuous
function. We assume that the nonlinear system (2) is bounded almost
everywhere, that is, $x(t)$ is in the $\mathcal{L}^{\infty}$ space. The
dynamics of the probability density of the state of the system
$\rho(x)\in\mathcal{L}_{1}$ is given by,

$\displaystyle\frac{\partial\rho(x,t)}{\partial
t}=\sum_{m=1}^{M}\frac{\partial(f_{m}^{T}(x)\rho(x,t))}{\partial
x_{m}}=A_{P}\rho(x,t).$ (3)

Here, $A_{P}$ is called the infinitesimal generator of the Perron-Frobenius
operator $P$ corresponding to the nonlinear system (2), and (3) is called the
Liouville equation [12]. It is sufficient that $\rho$ is any non-negative
function $\mathcal{L}_{1}$ as it can serve as a probability density function
after normalization [14]. Practically, the evolution of the system (2) can
only be measured by sampling the trajectory. Given a sampling period of
$\tau>0$, the discrete sampling of the system gives rise to a discretized
dynamics of the system defined (2), at times $t=0,\tau,2\tau,\cdots$. For
brevity, we use the notation $\rho_{k}(x)=\rho(x,k\tau)$. The discretized
evolution of $\rho_{k}(x)$ is given by

$\rho_{k+1}(x)=P_{\tau}\circ\rho_{k}(x),$ (4)

where $P_{\tau}$ is called the Perron-Frobenius operator parameterized by the
sampling period corresponding to the nonlinear dynamical system (2). We
illustrate in Figure 1 how the state trajectory $x(t)$ is coupled with the
probability density $\rho(x,t)$ for the Van der Pol oscillator,

$\displaystyle\dot{x}_{1}=x_{2},\quad\dot{x}_{2}=\mu(1-x_{1}^{2})x_{2}-x_{1}.$

(a) Initial states $x(0)$

(b) Initial density $\rho_{0}(x)$

(c) States $x(1500\tau)$

(d) Density $\rho_{1500}(x)$

Figure 1: Illustration of the relationship between the states and probability
density of the Van der Pol oscillator in a bounded domain. The x-axis in the
figures corresponds to $x_{1}$ and the y-axis in the figures corresponds to
$x_{2}$. (a) Initial states of several agents driven by Van der Pol dynamics.
(b) Initial density of state of agents. (c) States of the agent after time
$t=1500\tau$. (d) Density of the states at time $t=1500\tau$. Brighter colors
in (b) and (d) represent higher probability. The states are sampled with
$t=0.01$.

It is important to note that the Perron-Frobenius operator completely defines
the evolution of the density of the system. Hence, our goal is to analyze the
behavior of the AI-driven agents (2), through the estimation of the PF
operator.

Our second goal is to estimate the asymptotic behavior of AI-driven agents.
Several systems (2) exhibit stationary states asymptotically. For example,
robotics systems are designed to stabilize certain points in the domain.
Another example is a diffusion model which is trained to sample from unknown
target distributions. For systems that exhibit stationary states, there exists
an invariant density $\rho^{*}$ [12] for the PF operator such that

$\displaystyle\rho^{*}=P_{\tau}\circ\rho^{*}.$ (5)

Here agents following (2) reach $\rho^{*}$ asymptotically. We seek to estimate
the terminal density $\rho^{*}$ as it provides a convenient method to assess
the alignment of the AI-driven agents.

Our data consists of the state trajectory of $N$ identical agents governed by
the dynamics (2). The trajectory of these agents are collected from $t=0$ to
$t=T$ with a fixed sampling period $\tau=\frac{T}{K}$. The sampled dataset is
given by $\\{\mathcal{X}_{n}\\}_{n=1}^{N}$ of the state $x$, where
$\mathcal{X}_{n}=[\chi_{0}^{n},\chi_{1}^{n},\chi_{2}^{n},\cdots,\chi_{K}^{n}]\in\mathbb{R}^{M\times(K+1)}.$
Note that the collected data set can also come from a single agent starting
from $N$ different initial states. Given the data
$\\{\mathcal{X}_{n}\\}_{n=1}^{N}$, we seek to estimate the Perron-Frobenius
operator and the terminal density $\rho^{*}$ of the system.

##  3  Prediction Informed by Spectral-decomposition Algorithm (PISA) for
Learning Perron-Frobenius Operators

We present our algorithm to estimate the PF operator in this section. Further,
we predict the asymptotic behavior of the system by estimating the terminal
density of the dynamical system.

###  3.1 Density Estimation using Kernel Density Estimation

We employ Kernel Density Estimation to numerically construct the probability
density $\rho_{k}(x)$ using the data $\\{\mathcal{X}_{n}\\}$. We can view
$\\{\mathcal{X}_{n}\\}$ by iterating with respect to time as
$\\{\mathcal{Y}_{k}\\}_{k=0}^{K}$, where
$\mathcal{Y}_{k}=[\chi_{k}^{1},\chi_{k}^{2},\cdots,\chi_{k}^{N}]$ denotes the
state vectors of $N$ particles at time $t=k\tau$. Using kernel density
estimation [15], we then get an empirical probability distribution estimation
$\rho_{k}(x)$. In this paper, we choose to use the Gaussian kernel for the
estimation of $\rho_{k}$ which is given by

$\displaystyle\rho_{k}(x)=\frac{1}{N\sqrt{\det(2\pi\sigma_{k}^{2}I_{M})}}\sum_{n=1}^{N}e^{-\frac{\left\|x-\chi_{k}^{n}\right\|^{2}}{2\sigma_{k}^{2}}}.$
(6)

It is important to note that any choice of density estimation algorithm can be
used with the data $\\{\mathcal{X}_{n}\\}$ to obtain
$\\{\rho_{k}(x)\\}_{k=0}^{K}$. KDE provides a convenient choice for measuring
the probability $\rho(x)$ at fixed reference points. The choice of reference
points and the parameter $\sigma_{k}$ can be chosen by the user to better
approximate $\rho_{k}$. In this paper, we choose to uniformly sample the
reference points in the domain $X$ to estimate every $\rho_{k}(x)$. We fix a
constant $\sigma_{K}=\sigma$ for simplicity. More specific KDE-related tools
can be used based on the system in consideration and the domain, as enlisted
in [16].

###  3.2 Proposed Algorithm

We approximate the PF operator using the following model.

$\displaystyle\begin{aligned}
P_{\tau}\circ\rho(x)&=\rho(x)-\sum_{i=1}^{l}\left(\frac{1}{l}-\texttt{A}^{i}_{\theta}(\rho)\right)\texttt{G}^{i}_{\gamma}(x).\end{aligned}$
(7)

Here, we are decomposing the action of the PF operator on the density
$\rho(x)$ into $2l$ components as given by $l$ functions
$\texttt{A}^{i}_{\theta}(\rho)$ and $l$ functions
$\texttt{G}^{i}_{\gamma}(x)$. The functions $\texttt{A}^{i}_{\theta}(\rho)$
and $\texttt{G}^{i}_{\gamma}(x)$ are parameterized by $\theta$ and $\gamma$
respectively. This method of decomposing the PF operator is guided by the
spectral decomposition theorem which we elaborate on in Section 5.

Given the decomposition of the PF operator, we propose the following loss
function, guided by the spectral decomposition theorem, to learn the parameter
$\theta$ and $\gamma$.

$\displaystyle\begin{aligned}
L({\theta},\gamma)=&\sum_{k=0}^{K-1}D\left(\rho_{k}(x)-\sum_{i=1}^{l}\left(\dfrac{1}{l}-\texttt{A}^{i}_{\theta}(\rho_{k})\right)\texttt{G}^{i}_{\gamma}(x)\Bigg{\|}\rho_{k+1}(x)\right)+\lambda\sum_{i\neq
j}^{l}\braket{\texttt{G}^{i}_{\gamma}(x),\texttt{G}^{j}_{\gamma}(x)}\\\
&+\mu\sum_{r=1}^{l}D\left(\texttt{G}^{r}_{\gamma}(x)-\sum_{i=1}^{l}\left(\dfrac{1}{l}-\texttt{A}^{i}_{\theta}(\texttt{G}^{r}_{\gamma})\right)\texttt{G}^{i}_{\gamma}(x)\Bigg{\|}\texttt{G}^{r+1}_{\gamma}(x)\right).\end{aligned}$
(8)

We then construct PISA as the following alternating optimization algorithm to
compute $\theta$ and $\gamma$, in which we choose
$\texttt{A}^{i}_{\theta}(\rho)$ and $\texttt{G}^{i}_{\gamma}(x)$ to be outputs
of two distinct neural networks parameterized by $\theta$ and $\gamma$,
respectively.

Data: $l>0$, $\lambda>0$, $\mu>0$; $\rho_{k}(x)$, for $k=0,1,\cdots,K$;
initial values of ${\gamma}$ and ${\theta}$; two small positive thresholds
$\epsilon_{1}$ and $\epsilon_{2}$;

Result: $\gamma$ and $\theta$;

1 $N_{\text{epochs}}\leftarrow 1000$

2 while  _$N_{\text{epochs}}\neq 0$_ do

3

4 Solve the following optimization problem to get $\gamma^{*}$
$\displaystyle\begin{split}&\min_{\gamma}~{}~{}L(\theta,\gamma)\\\
&\text{s.t.}~{}~{}\texttt{G}^{i}_{\gamma}(x)\geq 0\text{ and
}\int\texttt{G}^{i}_{\gamma}(x)dx=1,\text{ for }i=1,\cdots,l;\end{split}$ (9)

5 if  _$\left\|\gamma^{*}-\gamma\right\|\geq\epsilon_{1}$_ then

6 $\gamma\leftarrow\gamma^{*}$

7 end if

8 Solve the following optimization problem to get $\theta^{*}$
$\displaystyle\begin{split}&\min_{\theta}~{}~{}L(\theta,\gamma)\\\
&\text{s.t.}~{}~{}\texttt{A}^{i}_{\theta}(\rho)\geq 0,\text{ for
}i=1,\cdots,l;\end{split}$ (10)

9 if  _$\left\|\theta^{*}-\theta\right\|\geq\epsilon_{2}$_ then

10 $\theta\leftarrow\theta^{*}$

11 end if

12 $N_{\text{epochs}}=N_{\text{epochs}}-1$

13 end while

Algorithm 1 : Prediction Informed by Spectral-decomposition Algorithm (PISA)

An important aspect of PISA is that it can also predict the terminal density
of the PF operator. The estimate of terminal density of $P_{\tau}$ can be
expressed as

$\rho^{*}(x)=\frac{1}{l}\sum_{i=1}^{l}\texttt{G}^{i}_{\gamma}(x).$ (11)

##  4 Numerical Experiments

We present the effectiveness of PISA on different numerical testbeds. We
performed the numerical experiments on a machine with Intel i9-9900K CPU with
128GB RAM and the Nvidia Quadro RTX 4000 GPU. In our numerical experiments, we
compare the performance of PISA with the that of [3]. Particularly, [3]
approximates the PF operator as

$\displaystyle\rho_{k+1}=e^{t\cdot\texttt{NN}_{\delta}(x,t)}\rho_{k}.$

Here, note that $\texttt{NN}_{\delta}$ approximates the Liouville operator
$A_{P}$ given in (3). Then $e^{t\cdot\texttt{NN}_{\delta}}$ is an
approximately linear solution to (4).

###  4.1  Unicycle Model with an NN-Based Controller

We first consider agents following the unicycle dynamics,

$\displaystyle\dot{x}_{1}=u_{1}\cos(x_{3}),\quad\dot{x}_{2}=u_{1}\sin(x_{3}),\quad\dot{x}_{3}=u_{2}.$
(12)

We consider the task of controlling $N$ unicycle agents using a diffusion-
model-based controller developed in [17]. The task is to start from a uniform
distribution and finally reach two Gaussian distributions centered at
$-4\mathbf{1}$ and $4\mathbf{1}$, respectively. With the number of agents as
$N=1000$, we generated trajectories of length $T=8$ seconds with a sampling
period of $\tau=0.01$. The initial and final states of the data are
represented in Figure 2(a). We predicted the terminal density $\rho^{*}$ as
given by (11). The estimated terminal density corresponding to the targets the
agents are trained to reach is depicted in Figure 2(b). In Figure 2(c), we
compare PISA against [3], using the same amount of training samples and the
same size of neural networks for both algorithms. We set $l=5$ for PISA and
used fully connected feed-forward neural networks with 3 hidden layers to
learn $\texttt{A}_{\theta}$, $\texttt{G}_{\gamma}$ and
$\texttt{NN}_{\delta}$.The two algorithms are trained for 1000 iterations and
tested on the testing dataset of length $T=3$ seconds. The y-axis of Figure
2(c) depicts the KL divergence between the predicted density and the true
density in the testing dataset. Initially, it is evident that [3] performs
better than PISA due to the nature of the model. Particularly, [3] chooses a
linear solution to the PDE (4) which explains the better performance initially
but results in a rapid decrease in performance. We see that PISA performs
better by one order of magnitude than [3] over a long time horizon.

$x_{2}$$x_{1}$(a) Dataset

$x_{2}$$x_{1}$(b) Predicted $\rho^{*}$

$D(\hat{\rho}_{t}\|\rho_{t})$time(c) Comparison with [3]

Figure 2: Experiments on the unicycle model data generated using controller
from [17]. (a) Initial states are uniformly distributed, depicted in blue,
final states are in red. (b) The estimated $\rho^{*}$ by PISA. (c) Performance
comparison between PISA and [3] on testing data of length $3$ seconds.
Performance metric is the KL divergence between predicted and true densities.
Our dataset contains $K=800$ sampled instants, wherein we use the first $500$
samples of each trajectory as the training data set and the remaining $300$ as
the testing set.

###  4.2 Predicting Behavior of Score-Based Generative Model

We consider the problem of analysis of the behavior of generative models.
Generative models behave as evolving agents by taking data samples from an
initial noise density and reaching an unknown target density. Generative
models are complex and their behavior can be analyzed through the trajectory
of converting noise samples to data samples. Our task is to analyze the
behavior of generative models from the lens of evolving probability densities
using their sampling trajectories. We particularly consider the case of
diffusion models based on estimating the score [18]. In the diffusion model we
use the data lies in five dimensions. In score-based generative models, new
data samples from unknown target distributions are generated using a bi-
directional scheme. Given data from an unknown target distribution, noise is
sequentially added to data samples using a Stochastic Differential Equation in
the forward process until a desired noise distribution is reached. To sample
new data from the target distribution, this forward process is reversed so
that a sample starts from the noise distribution and ends as a sample from the
target distribution. Particularly, [18] reverses the forward process by
learning the _score_ of the data distribution. We seek to study the behavior
of these score-based generative models by their action on samples in the
reverse process. We show that we could not only predict the behavior of
diffusion models but also potentially identify the target distribution by
estimating $\rho^{*}$.

We consider the following ODE for the score-based generative model,

$\displaystyle\dot{x}=u$ (13)

Here, $x,u\in\mathbb{R}^{5}$. In the forward process
$u~{}\sim\mathcal{N}(\bm{0},I)$, whereas in the reverse process, $u$ denotes
the output of a neural network $\texttt{S}_{\psi}(x,t)$ that estimates the
score. The neural network $\texttt{S}_{\psi}(x,t)$ with five outputs
approximates the score of the forward process. In this experiment, the task is
to sample from Gaussians centered at $-4\bm{1}$ and $3\bm{1}$. The noise
distribution is the uniform distribution over the domain $[-8,8]^{5}$. In the
forward process, we start with $N=3000$ samples from
$\mathcal{N}_{1}(-4\bm{1},0.2I)$ and $\mathcal{N}_{1}(3\bm{1},0.2I)$. The data
samples are diffused in the forward process for a time period of 6 seconds. In
the reverse process, to sample from the desired distributions, we learn the
score as proposed in [18]. Once the score is sufficiently learned using a
neural network, we record $N=3000$ trajectories in the reverse process for a
time period of 6 seconds. The first four seconds of the dataset constitute our
training dataset and the last 2 seconds constitute the testing dataset.

It is evident from Figure 3(a) that PISA makes an accurate estimation of the
target distribution. We can see that the estimated $\rho^{*}$ corresponds to
that of the aforementioned target distribution. This implies that the
diffusion model is well aligned with its intended task. Further, in Figure
3(b), we compare PISA with [3] in predicting the evolution of $\rho_{k}$. We
choose $l=10$ for PISA and feedforward neural networks with 3 hidden layers
for $\texttt{A}^{i}_{\theta}$, $\texttt{G}^{i}_{\gamma}$ and
$\texttt{NN}^{i}_{\delta}$. We use the KL divergence between the predicted
$\hat{\rho}_{k}$ and the true $\rho_{k}$ from the testing dataset. Once again,
we see that [3] performs better initially due to its linear solution to the
PDE but its performance deteriorates rapidly. PISA retains relatively much
better performance over a longer time horizon.

$x_{2}$$x_{1}$(a) Predicted $\rho^{*}$

$D(\hat{\rho}_{t}\|\rho_{t})$time(b) Comparison with [3]

Figure 3: Experiments on the five dimensional score-based generative model
[18]. (a) Predicted terminal density $\rho^{*}$ projected on the first two
dimension. Brighter colours indicate a higher probability. (b) Comparison of
performance of PISA and [3] on testing dataset.

###  4.3 UCY Pedestrian Dataset

We show the effectiveness of PISA on a physical data by applying it to the UCY
pedestrian dataset [19]. Here, the task is to predict the movement of
pedestrians by estimating the evolution of the density of pedestrians. The
dataset consists of videos of pedestrians walking in several regions as
depicted in Figure 4(a). We use the Zara01 subsection of the dataset in our
experiments. We obtained pre-processed data from the code repository of [20],
where the video was processed to obtain the $x$ and $y$ coordinates of the
position of the pedestrians. Here we consider each pederstrian as an agent in
our analysis. We assume that every pedestrian is identical and their movement
is governed by the dynamics given in (2). Given the positions of pedestrians
as depicted in Figure 4(a), we approximate the probability density of the
pedestrians as depicted in Figure 4(b).

In Figure 4(c), we once again compare PISA with the exponential model [3] on
the test data for the first 200 time samples. We choose $l=5$ and feedforward
neural networks with 3 hidden layers for $\texttt{A}^{i}_{\theta}$,
$\texttt{G}^{i}_{\gamma}$ and $\texttt{NN}^{i}_{\delta}$. Here, we see that
the initial time period in which the exponential model works better than PISA
is significantly shorter due to the model inaccuracy. However, PISA continues
to perform well over a longer time horizon. It is also important to note that
both models have significantly higher estimation errors in the testing
performance for this experiment compared to performance in experiments on the
unicycle model and the score-based generative model. This is due to the
stochasticity of the data and the assumption we make about the nature of the
pedestrians. Further, with new pedestrians entering the scene and existing
pedestrians leaving the scene, there are jumps in the probability density
which further reduces the performance of KDE-based methods. Better density
approximation algorithms that are suited for stochastic data and for
incorporating jumps in the probability density can be employed to obtain an
improvement in the performance.

(a) Snapshot from dataset

$x_{2}$$x_{1}$(b) Estimated density $\rho_{k}$

$D(\hat{\rho}_{t}\|\rho_{t})$time(c) Comparison with [3]

Figure 4: Experiments on the UCY pedestrian dataset. (a) A snapshot from the
dataset. (b) Corresponding estimated probability density. (c) Comparison
between PISA and [3] in the estimation of future probability densities.

##  5 Theoretical Foundations of Spectral Decomposition

The Perron-Frobenius operator $P_{\tau}$ is a constrictive Markov operator
[12] that pushes forward the probability distribution $\rho_{k}(x)$ to a
stationary distribution $\rho^{*}(x)$ corresponding to the attractors of
dynamical systems. This evolution of the probability distribution is in fact a
Markov Process. For Markov operators with the constrictive property, we have
the following spectral decomposition theorem.

###### Lemma 1

[12, 13] Let $P$ be a constrictive Markov operator. Then there exists an
integer $l$, two sequences of non-negative functions
$g_{i}(x)\in\mathcal{L}_{1}$ and $h_{i}(x)\in\mathcal{L}_{\infty}$,
$i=1,2,\cdots,l,$ and an operator $Q:\mathcal{L}_{1}\mapsto\mathcal{L}_{1}$
such that for all $\rho(x)\in\mathcal{L}_{1}$, $P\circ\rho(x)$ can be written
in the form

$P\circ\rho(x)=\sum_{i=1}^{l}a_{i}(\rho)g_{i}(x)+Q\circ\rho(x),$ (14)

where

$a_{i}(\rho)=\int\rho(x)h_{i}(x)dx.$

The functions $g_{i}(x)$ and the operator $Q$ have the following properites:

1) Each $g_{i}(x)$ is normalized to one and

$g_{i}(x)g_{j}(x)=0,\quad\text{ for all }i\neq j,$ (15)

i.e., the density functions $g_{i}(x)$ have disjoint supports;

2) For each integer $i$ there exists a unique integer $\alpha(i)$ such that

$P\circ g_{i}(x)=g_{\alpha(i)}(x).$ (16)

where $\alpha(i)\neq\alpha(j)$ for $i\neq j$. Thus, $P$ just permutes the
functions $g_{i}(x)$;

3) Moreover,

$\|P^{n}Q\circ\rho(x)\|\rightarrow 0$ (17)

as $n\rightarrow\infty$ for every $\rho(x)\in\mathcal{L}_{1}$.

Lemma 1 states that the action of the PF operator can be decomposed into $l$
components through the functionals $a_{i}(\rho)$ and the functions $g_{i}(x)$.
Here $l$ is a finite integer that serves as a measure of the model complexity
of PISA. Further, the operator $Q$ captures the effect of the terminal density
on $\rho(x)$. As $t\rightarrow\infty$, the action of $Q$ on $\rho(x)$ decays
to 0. This drives our motivation to use $Q\circ\rho(x)$ as

$\displaystyle Q\circ\rho(x)=\rho(x)-\rho^{*}(x).$ (18)

Further, as $t\rightarrow\infty$, we can see that

$\displaystyle\rho^{*}(x)=\frac{1}{l}\sum_{i}^{l}g_{i}(x).$ (19)

This implies that the density functions $g_{i}$ serve as a basis for the
stationary terminal density $\rho^{*}$. It is easy to verify for (19) that
$P\circ\rho^{*}(x)=\rho^{*}(x)$ through the permutation property. Given the
Lemma 1, (18), and (19), we provide a sufficient condition on the output of
our algorithm PISA.

###### Theorem 1

For systems (2) that have a stationary terminal density, there exists a finite
$l$, an operator $Q$, $l$ non-negative functionals
$\texttt{A}^{i}_{\theta}(\rho)$ and $l$ densities $\texttt{G}_{\gamma}(x)$
such that the loss $L(\theta,\gamma)=0$.

###### Proof 1

We provide a brief overview of the proof. Lemma 1 guarantees the existence of
$l$ functions $a_{i}(\rho)$ and $g_{i}(x)$ that exactly decompose the action
of the PF operator. These are approximated using neural networks
$\texttt{A}^{i}_{\theta}(\rho)$ and $\texttt{G}^{i}_{\gamma}(x)$,
respectively. The cost function $L$ is designed to satisfy the properties of
$a_{i}$ and $g_{i}$. The first term in the cost function $L$ addresses the
propagation of the PF operator. The second term addresses the orthogonality
property of every $g_{i}$ and the last term captures the permutative property
$g_{i}$. $\blacksquare$

###### Remark 1

It should be noticed that the spectral decomposition theorem guarantees the
existence of a decomposition in the form of (14), it yet does not specify the
form of $Q$ and $g_{i}(x)$, as well as the number $l$. Also, there is no
guarantee that the decomposition is unique. Therefore, our choice that
$Q\circ\rho(x)=\rho(x)-\rho^{*}(x)$ is just one of many possibilities. An
alternative for the action of $Q$ is $Q\circ\rho(x)=D(\rho\|\rho^{*})$, which
is also a decreasing quantity. This can be verified by the laws of
thermodynamics which guarantee decreasing of relative entropy. Empirically, we
have seen that $Q\circ\rho(x)=\rho(x)-\rho^{*}(x)$ serves as a better choice
for $Q$. Moreover, Lemma 1 guarantees a finite model complexity through the
number of basis functions $l$. This is a tunable parameter in our algorithm
that can adjusted to achieve a desired accuracy of prediction. As described in
our numerical experiments, it turns out that PISA makes good approximations of
both (7) and (11). Our relatively small neural networks to approximate
$a_{i}(\rho)$ and $g_{i}(x)$ through $\texttt{A}^{i}_{\theta}$ and
$\texttt{G}^{i}_{\gamma}$, respectively, are empirically efficient choices.

##  6 Conclusion

In this study, we explore the task of predicting the behavior of agents
controlled by AI models, using a probabilistic framework to analyze the
evolution of probability densities. Our proposed algorithm, PISA, effectively
estimates the Perron-Frobenius (PF) operator that characterizes the evolution
of these densities, thus enabling predictions of both the short and long-term
behavior of AI-driven agents. This ability to forecast asymptotic behaviors is
critical for assessing whether such agents align with human values and
requirements. Currently, our approach utilizes kernel density estimation to
approximate the probability densities from individual trajectories; however,
this method introduces potential inaccuracies. Optimizing the choice of kernel
functions and adopting more sophisticated density estimation techniques could
significantly enhance our model’s performance. Additionally, a comprehensive
evaluation of the model’s complexity and computational efficiency remains to
be conducted, which will be crucial for practical applications and further
scalability.

## References

  * [1] D. Goswami, E. Thackray, and D. A. Paley.  Constrained ulam dynamic mode decomposition: Approximation of the perron-frobenius operator for deterministic and stochastic systems.  IEEE control systems letters, 2(4):809–814, 2018.
  * [2] R. A. Norton, C. Fox, and M. E. Morrison.  Numerical approximation of the frobenius–perron operator using the finite volume method.  SIAM Journal on Numerical Analysis, 56(1):570–589, 2018.
  * [3] Y. Meng, D. Sun, Z. Qiu, M. T. B. Waez, and C. Fan.  Learning density distribution of reachable states for autonomous systems.  In Proceedings of the 5th Conference on Robot Learning, volume 164 of Proceedings of Machine Learning Research, pages 124–136. PMLR, 08–11 Nov 2022.
  * [4] C. Huang, J. Fan, W. Li, X. Chen, and Q. Zhu.  Reachnn: Reachability analysis of neural-network controlled systems.  ACM Transactions on Embedded Computing Systems (TECS), 18(5s):1–22, 2019.
  * [5] M. Everett, G. Habibi, C. Sun, and J. P. How.  Reachability analysis of neural feedback loops.  IEEE Access, 9:163938–163953, 2021.
  * [6] C. Zhang, W. Ruan, and P. Xu.  Reachability analysis of neural network control systems.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, pages 15287–15295, 2023.
  * [7] Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, and A. Anandkumar.  Physics-informed neural operator for learning partial differential equations.  ACM/JMS Journal of Data Science, 2021.
  * [8] T. De Ryck and S. Mishra.  Generic bounds on the approximation error for physics-informed (and) operator learning.  Advances in Neural Information Processing Systems, 35:10945–10958, 2022.
  * [9] S. Surasinghe, J. Fish, and E. M. Bollt.  Learning transfer operators by kernel density estimation.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 34(2), 2024.
  * [10] Y. Hashimoto, M. Ikeda, and H. Kadri.  Deep learning with kernels through rkhm and the perron-frobenius operator.  Advances in Neural Information Processing Systems, 36, 2024.
  * [11] Y. Yang, L. Nurbekyan, E. Negrini, R. Martin, and M. Pasha.  Optimal transport for parameter identification of chaotic dynamics via invariant measures, 2022.
  * [12] A. Lasota and M. C. Mackey.  Chaos, fractals, and noise: stochastic aspects of dynamics, volume 97.  Springer Science & Business Media, 2013.
  * [13] A. Boyarsky and R. Levesque.  Spectral decomposition for combinations of markov operators.  Journal of mathematical analysis and applications, 132(1):251–263, 1988.
  * [14] P. Bevanda, S. Sosnowski, and S. Hirche.  Koopman operator dynamical models: Learning, analysis and control.  Annual Reviews in Control, 52:197–212, 2021.
  * [15] T. Hastie, R. Tibshirani, and J. H. Friedman.  The elements of statistical learning: data mining, inference, and prediction, volume 2.  Springer, 2009.
  * [16] Yen-Chi C.  A tutorial on kernel density estimation and recent advances, 2017.
  * [17] K. Elamvazhuthi, D. Gadginmath, and F. Pasqualetti.  Denoising diffusion-based control of nonlinear systems, 2024.
  * [18] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole.  Score-based generative modeling through stochastic differential equations.  arXiv preprint arXiv:2011.13456, 2020.
  * [19] A. Lerner, Y. Chrysanthou, and D. Lischinski.  Crowds by example.  In Computer graphics forum, volume 26, pages 655–664. Wiley Online Library, 2007.
  * [20] T. Salzmann, B. Ivanovic, P. Chakravarty, and M. Pavone.  Trajectron++: Dynamically-feasible trajectory forecasting with heterogeneous data.  In Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XVIII 16, pages 683–700. Springer, 2020.

##  7 Supplemental Material

###  7.1 Preliminaries and Mathematical Notations

KL divergence between two probability distributions $f(x)$ and $g(x)$ is
defined as

$D(f\|g)=\int f(x)\log\frac{f(x)}{g(x)}dx,$

which serves as a measure of the distance between the two distributions.
Notice also that the $D(f\|g)=0$ if and only if $f(x)-g(x)=0$ almost
everywhere. Moreover, we denote $\mathbb{R}^{M}$ the vector space of all $M$
dimensional real vectors. The inner product of two vectors (functions) $f(x)$
and $g(x)$ in Hilbert space is defined as

$\braket{f(x),g(x)}=\int f^{\prime}(x)g(x)dx,$

where $f^{\prime}(x)$ denotes the conjugate transpose of $f(x)$. Also, we use
$\|\cdot\|$ to denote the norm of a function (vector) in Banach space. The
$\mathcal{L}_{2}$ norm of a vector $f(x)$ in Hilbert space is defined as

$\|f\|_{2}=\sqrt{\int f^{\prime}(x)f(x)dx}.$

The $\mathcal{L}_{1}$ norm of a vector $f(x)$ in Banach space is defined as

$\|f\|_{1}={\int|f(x)|dx}.$

The $\mathcal{L}_{\infty}$ norm of a vector $f(x)$ in Banach space is defined
as

$\|f\|_{\infty}=\sup_{x}{|f(x)|dx}.$

We also denote $\mathcal{L}_{1}$ the set of all absolutely integrable
functions and $\mathcal{L}_{\infty}$ the set of all functions which are almost
bounded everywhere.

###  7.2 Proposition of Perron-Frobenius Operators

###### Proposition 1

[12] Suppose we have a nonlinear bounded dynamical system (2) and the
corresponding Perron-Frobenius operator $P_{\tau}$. Then it follows that

  * •

$P_{\tau}$ is a linear operator;

  * •

$P_{\tau}\circ\rho(x)$ is non-negative if $\rho(x)$ is non-negative;

  * •

Integral invariance:

$\int P_{\tau}\circ\rho(x)dx=\int\rho(x)dx;$

  * •

$P_{\tau}$ has a fixed point (probability distribution function) such that

$P_{\tau}\circ\rho^{*}(x)=\rho^{*}(x)$

and

$\lim_{n\rightarrow\infty}P_{\tau}^{n}\circ\rho(x)=\rho^{*}(x)$

which is also denoted in some literature that $P_{\tau}$ has a preserved
measure $\mu_{\rho^{*}}$.


# Negative temperature is cool for cooling

Gabriella G. Damas  Instituto de Física, Universidade Federal de Goiás,
74.001-970, Goiânia - GO, Brazil  Rogério J. de Assis  Instituto de Física,
Universidade Federal de Goiás, 74.001-970, Goiânia - GO, Brazil  Departamento
de Física, Universidade Federal de São Carlos, 13.565-905, São Carlos - SP,
Brazil  Norton G. de Almeida  Instituto de Física, Universidade Federal de
Goiás, 74.001-970, Goiânia - GO, Brazil

###### Abstract

In this work, we study an autonomous refrigerator composed of three qubits
[Phys. Rev. Lett. 105, 130401 (2010)] operating with one of the reservoirs at
negative temperatures, which has the purpose of cooling one of the qubits. We
find the values of the lowest possible temperature that the qubit of interest
reaches when fixing the relevant parameters, and we also study the limit for
cooling the qubit arbitrarily close to absolute zero. We thus proceed to a
comparative study showing that reservoirs at effective negative temperatures
are more powerful than those at positive temperatures for cooling the qubit of
interest.

###### pacs:

05.30.-d, 05.20.-y, 05.70.Ln

##  I Introduction

In a 2010 work, Linden et al. (Ref. (Linden  _et al._ , 2010)) studied three
models for cooling qubits or arbitrary quantum systems. The first model uses
only qubits, the second one a qubit and a qutrit with nearest-neighbor
interactions, and the third one a single qutrit. The model using only qubits
is a self-contained or autonomous refrigerator (AR) with three interacting
qubits, one of them being the object to be cooled. As shown in Ref. (Linden
_et al._ , 2010), under the condition of perfect insulation of the qubit to be
refrigerated, this model presents no fundamental cooling limit, and thus,
depending on the parameters used, the qubit of interest can be arbitrarily
cooled towards absolute zero.

In this work, we study the AR consisting only of qubits as proposed by Linden
et al. in Ref. (Linden  _et al._ , 2010), putting it to work with one of the
reservoirs at negative temperature. Negative absolute temperatures were
initially considered in 1951 when Purcell first produced spin states with
population inversion (Purcell and Pound, 1951). In 1956, Ramsey theoretically
studied these states, treating them as thermodynamic equilibrium states
(Ramsey, 1956). More than 60 years after the experiment performed by Purcell,
other experiments involving negative temperatures followed (Carr, 2013; Braun
_et al._ , 2013), which drew attention to this topic. While some authors
defend the existence of negative thermodynamic temperatures, some works
associate them with non-equilibrium states and thus refer to them as effective
(or apparent) negative temperatures, as opposed to equilibrium states at
positive thermodynamic temperatures. (Dunkel and Hilbert, 2013; Sokolov, 2013;
Hilbert  _et al._ , 2014; Campisi, 2015; Cerino  _et al._ , 2015; Frenkel and
Warren, 2015; Swendsen and Wang, 2015; Poulter, 2016; Swendsen and Wang, 2016;
Abraham and Penrose, 2017; Hama  _et al._ , 2018; Struchtrup, 2018). As a
general property, experimental realizations of reservoirs with effective
negative temperatures, such as the ones we will consider in this work,
requires population inversion (Carr, 2013; Braun  _et al._ , 2013; de Assis
_et al._ , 2019; Mendonça  _et al._ , 2020). As we will show, the qubit-based
model AR has advantages when operating with one of the reservoirs at negative
temperature compared to operating only with conventional reservoirs at
positive temperatures.

This work is organized as follows: In Section II, we present the model for the
study of AR considering bosonic reservoirs, for which the temperature is
always positive, and fermionic reservoirs, which can reach negative
temperatures. In Section III, we present our results, showing that the use of
a fermionic reservoir with negative temperatures has an advantage over the use
of only bosonic reservoirs. Finally, in Section IV, we present our conclusions
about the comparative study of the AR involving positive and negative
temperatures.

##  II Model

The AR considered here consists of three qubits, each in contact with a
different reservoir (Linden  _et al._ , 2010). The first qubit, which is the
qubit to be cooled, is put in contact with a cold reservoir at temperature
$T_{c}$. The second qubit, analogously to the classical refrigerator, plays
the role of the spiral and is put in contact with a reservoir at "room"
temperature $T_{r}$, while the third qubit works as the engine staying in
contact with a hot reservoir at temperature $T_{h}$. For the AR to work
properly, considering only reservoirs at positive temperatures, the relations
$T_{c}<T_{r}<T_{h}$ and $E_{3}=E_{2}-E_{1}$ \- see Fig. 1(a), must be
satisfied (Linden  _et al._ , 2010).

Figure 1:  Schematic representation of possible level configuration for the
qubit-based AR cooler. (a) Level configurations used in Ref. (Linden  _et al._
, 2010) for positive temperatures, where the condition $E_{3}=E_{2}-E_{1}$ is
required. (b) A possible level configuration with $E_{3}=E_{2}=E_{1}$ for
negative temperatures, where the condition $E_{3}=E_{2}-E_{1}$ is no longer
needed. In (a) and (b), $\left|g\right\rangle_{k}$ and
$\left|e\right\rangle_{k}$ are the ground and excited states of the $k$th
qubit, respectively.

In the weak coupling limit and the Markovian regime, the dynamics governing
the AR is dictated by the master equation

$\frac{\partial\rho}{\partial t}=-i\left[H_{0}+H_{int},\rho\right]+\\\
\frac{1}{2}\sum_{k=1}^{3}\Gamma_{k}^{\downarrow}\left(2\sigma_{-,k}\rho\sigma_{+,k}-\left\\{\sigma_{+,k}\sigma_{-,k},\rho\right\\}\right)+\\\
\frac{1}{2}\sum_{k=1}^{3}\Gamma_{k}^{\uparrow}\left(2\sigma_{+,k}\rho\sigma_{-,k}-\left\\{\sigma_{-,k}\sigma_{+,k},\rho\right\\}\right),$
(1)

where $H_{0}$ and $H_{int}$ are respectively the free qubits Hamiltonian and
the three-body interaction Hamiltonian, given by

$H_{0}=\frac{1}{2}E_{1}\sigma_{z,1}+\frac{1}{2}E_{2}\sigma_{z,2}+\frac{1}{2}E_{3}\mathbb{\sigma}_{z,3}$
(2)

and

$H_{int}=g\left(\sigma_{-,1}\sigma_{+,2}\sigma_{-,3}+\sigma_{+,1}\sigma_{-,2}\sigma_{+,3}\right).$
(3)

In the above equations, $\sigma_{-,k}$, $\sigma_{+,k}$, and $\sigma_{z,k}$
($k=1,2,3$) are respectively the lowering, raising, and z Pauli operators for
the $k\text{th}$ qubit; $E_{1}$, $E_{2}$, and $E_{3}$ are the energy gaps for
each qubit, and $g$ is the coupling constant. Also, we let
$\Gamma_{k}^{\downarrow}=\gamma_{k}\left(1\pm n_{\alpha_{k}}\right)$, where
the signal "+" ("-") is for a bosonic (fermionic) reservoir, and
$\Gamma_{k}^{\uparrow}=\gamma_{k}n_{\alpha_{k}}$ ($\alpha_{1}=c$,
$\alpha_{2}=r$, and $\alpha_{3}=h$), with $\gamma_{k}$ being the dissipation
rate for the $k$th qubit and $n_{\alpha_{k}}$ being the average excitation
number of the reservoir at temperature $T_{{\alpha_{k}}}$. Note that Eq. (1)
describes the dynamics dictated by either bosonic (Carmichael, 1999; Breuer
and Petruccione, 2007) or fermionic (Artacho and Falicov, 1993; Ghosh  _et
al._ , 2012) thermal reservoirs. For bosonic thermal reservoirs
$n_{\alpha_{k}}=1/\left(\text{e}^{E_{k}/T_{\alpha_{k}}}-1\right)$, while for
fermionic reservoirs
$n_{\alpha_{k}}=1/\left(e^{E_{k}/T_{\alpha_{k}}}+1\right)$. As said before,
negative temperatures are characterized by population inversion, which occurs
when $n_{\alpha_{k}}>1/2$, and it is only possible for fermionic reservoirs.
To find the final temperature of qubit 1, we numerically calculate the
asymptotic equilibrium state of Eq. (1). We perform these numerical
calculations using the quantum optics toolbox (Johansson  _et al._ , 2012,
2013).

##  III Results

To compare the AR operation at positive temperature with its operation at
negative temperature, we must consider the same parameters that characterize
the AR, such as coupling constant, energy gaps and dissipation rates. For
clarity, let us start detailing the positive temperature case by fixing the
parameters used in Ref. (Linden  _et al._ , 2010). Thus, we fixed the energies
gaps $E_{1}=1$, $E_{2}=5$, and $E_{3}=4$, which is a special case of the more
general relation $E_{3}=E_{2}-E_{1}$ to AR work properly, and the temperatures
$T_{c}=1,1.5,2$, and $T_{r}=2$. Finally, we fixed the dissipation rates
$\gamma_{1}=\gamma_{2}=\gamma_{3}=1$ and let $T_{h}$ vary from $1$ to $10$.
Cooling is achieved when $T_{1}-T_{c}<0$, which means that the temperature
$T_{1}$ of the qubit $1$ is lower than the temperature $T_{c}$ of the cold
reservoir - See Fig. 2(a). In Fig. 2(b) we show how the temperature $T_{1}$
drops by increasing $T_{h}$. Note from Fig. 2(b) that the temperature of qubit
1 stabilizes at a constant value no matter how high $T_{c}$ is. We numerically
found that the corresponding lowest temperatures reached by the qubit 1 are
$T_{1}=0.9486$ (when $T_{c}=1$), $T_{1}=1.4054$ (when $T_{c}=1.5$), and
$T_{1}=1.867$ (when $T_{c}=2$). These lowest possible temperatures, when
fixing the set of parameters, are dictated by the population of the ground
state $p_{g,1}$, which is given by
$p_{g,1}=1/\left(1+e^{-E_{1}/T_{1}}\right)$. Solving for temperature, we get

$T_{1}=\frac{E_{1}}{\ln\left(\frac{p_{g,1}}{1-p_{g,1}}\right)}.$ (4)

From this equation, we see that $T_{1}$ can reach any value since
$0.5<p_{g,1}<1.0$. However, due to imperfect insulation, the heat flow to
qubit 1 prevents it from reaching temperatures arbitrarily close to absolute
zero. In fact, in our numerical calculations, we found that if we start at low
temperatures as $T_{c}<0.48$, then $T_{1}-T_{c}>0$, meaning that the
refrigerator no longer works. As we shall see, using the same parameters but
considering one of the reservoirs at negative temperatures, we can cool the
qubit 1 below $T_{c}=0.48$.

Figure 2:  (a) $T_{1}-T_{c}$ versus $T_{h}>0$ for three different values of
$T_{c}$: $T_{c}=1$ (dot green line), $T_{c}=1.5$ (dash red line) and $T_{c}=2$
(solid blue line). Refrigeration occurs for $T_{1}-T_{c}<0$. (b) $T_{1}$
behavior for those three initial conditions given in (a). Note that the lowest
temperature of the qubit $T_{1}$ does not decrease any further, no matter how
high $T_{h}$ becomes. The lowest values reached by $T_{1}$ are $T_{1}=0.9486$
(when $T_{c}=1$), $T_{1}=1.4054$ (when $T_{c}=1.5$) and $T_{1}=1.867$ (when
$T_{c}=2$).

Next, we use the perfect insulation condition for the temperature of qubit 1
to get as close to absolute zero as possible (Linden  _et al._ , 2010). This
condition is obtained by isolating qubit 1 from its reservoir, thus letting
$\gamma_{1}\rightarrow 0$. This condition allows us to reach lower
temperatures for qubit 1 than that shown in Fig. 2(b) and obtain the following
analytical expression for the cooling temperature (Linden  _et al._ , 2010)

$T_{1}=\frac{T_{c}}{1+\frac{E_{3}}{E_{1}}\left(1-\frac{T_{c}}{T_{h}}\right)}.$
(5)

The above formula tell us that, under the condition of perfect insulation, it
is possible to cool down to absolute zero, provided that
$E_{3}/E_{1}\rightarrow\infty$.

Now we turn our attention to negative temperatures. In this case, we must
indicate which reservoir will be at negative temperatures. Knowing that the
engine is modeled by qubit 3, which is in contact with a reservoir at hot
temperature $T_{h}$, it is natural to replace the reservoir of qubit 3 by a
reservoir at $T_{h}<0$ since bodies at negative temperatures are known to be
hotter than any body at positive temperature as it scales, from cold to hot,
according to $+0\text{ K},...,+300\text{ K},...,+\infty\text{ K}$,
$-\infty\text{ K},...,-300\text{ K},...,-0\text{ K}$. In view of that, one
might think that instead of using negative temperatures it would suffice to
use arbitrarily hot reservoirs at positive temperatures. However, we saw from
Fig. 2(b) that this is not the case: no matter how high $T_{h}$ is, there is a
limit achieved by $T_{1}$. Remarkably, we will show that under the same
conditions, i.e., the same set of parameters characterizing the AR, the lowest
possible temperature of qubit 1 cooled using positive temperatures can be
pushed even lower using a hot reservoir at negative temperatures.

As said before, to compare the cooling process at either negative or positive
temperatures, we put the AR to work with the same parameters. The steady state
is found using the Eq. (1), but now considering a fermionic reservoir with
$n_{3}>0.5$ ($T_{h}<0$). In Fig. 3(a), we show the difference $T_{1}-T_{c}$ as
a function of $T_{h}<0$. Interestingly, if we compare with Fig. 2(a) for
$T_{h}>0$, we see that at negative temperatures the steady-states are always
cooled, no matter the initial temperature we choose for $T_{h}<0$. Also,
remembering that for $T_{h}<0$ the higher temperatures ranges according to
$\left(-\infty,0^{-}\right)$ in Fig. 3(b) we see that the higher the
temperature $T_{h}$ the lower is the temperature of qubit 1.

Figure 3:  (a) $T_{1}-T_{c}$ versus $T_{h}<0$ using the same three values for
$T_{c}$ as before: $T_{c}=1$ (dot green line), $T_{c}=1.5$ (dash red line) and
$T_{c}=2$ (solid blue line). Refrigeration occurs for $T_{1}-T_{c}<0$. Note
that the higher the temperature $T_{h}$ the greater the cooling of qubit 1, as
before, but now the steady state of qubit 1 is cooled for all negative
temperatures. (b) $T_{1}$ versus $T_{h}$ for those three initial conditions
given in (a). Again, the lowest temperature of the qubit $T_{1}$ does not
decrease any further by increasing $T_{h}$. However, this time the lower
temperatures are closer to absolute zero than when using positive
temperatures. The corresponding lowest temperatures are $T_{1}=0.7805$ (when
$T_{c}=1$ ), $T_{1}=1.1615$ (when $T_{c}=1.5$) and $T_{1}=1.5568$ (when
$T_{c}=2$).

Note that the temperature reached by qubit 1 is yet given by equation Eq. (4).
Again, imperfect insulation of qubit 1 prevents it from reaching temperatures
arbitrarily close to absolute zero, as also shown in Figs. 3(a) and 3(b),
similarly to what happens in the case of $T_{h}>0$. To the set of parameters
we are using, these values are $T_{1}=0.7805$ (when $T_{c}=1$), $T_{1}=1.1615$
(when $T_{c}=1.5$), and $T_{1}=1.5568$ (when $T_{c}=2$). Interestingly, note
that the lowest temperatures $T_{1}$ in the case $T_{h}<0$ are lower than the
temperatures $T_{1}$ corresponding to the case $T_{h}>0$. That is true for all
reference temperatures $T_{c}$ according to our numerical calculation and as
can be seen in few examples shown in the table of Fig. 4(a). Taking $T_{c}$ as
the reference temperature, we can calculate the percentage decrease of $T_{1}$
in relation to $T_{c}$ for both negative and positive temperatures. In Fig.
4(b) we show the corresponding percentages for both $T_{h}>0$ (blue) and
$T_{h}<0$ (red).

Fig. 4 makes it clear that negative temperatures are more effective for
cooling a qubit, as it drives qubit 1 to lower temperatures by a large
percentage. According to our numerical calculations, for this set of
parameters, which is the same as used for positive temperatures, when the
reference temperature $T_{c}$ is as low as $T_{c}=0.48$ the AR stops cooling
if $T_{h}>0$. However, at $T_{h}<0$ the AR continues to cool down until qubit
1 reaches the limit of $T_{1}=0.0275$, that is, an order of magnitude lower.

We can also ask what happens when we impose the condition for perfect
insulation, $\gamma_{1}\rightarrow 0$. In this case, we can demonstrate that
Eq. (5) remains valid. Alternatively, we can write Eq. (5) as

$T_{1}=\frac{T_{C}}{1+\frac{E_{3}}{E_{1}}(1+\frac{T_{C}}{\left|T_{h}\right|})}.$
(6)

Eq. (6) tells us that, as in the case of positive temperatures, it is possible
to cool a qubit to temperatures arbitrarily close to absolute zero. Note,
however, that for cooling a qubit toward zero, since $T_{c}/T_{h}<1$, while
for positive temperatures, Eq. (5), we need to resort only to the ratio
$E_{3}/E_{1}\rightarrow\infty$, for negative temperatures we can have
$T_{c}>\left|T_{h}\right|$. Thus, it is possible to take advantage of negative
temperatures by choosing $T_{c}/\left|T_{h}\right|\rightarrow\infty$.

Regarding the ratio $E_{3}/E_{1}$, there is an important remark: positive
temperatures require the condition $E_{3}=E_{2}-E_{1}$, which limits the level
structures of the qubits needed for AR to work. On the other hand, negative
temperatures allow arbitrary configurations, such as $E_{3}=E_{2}=E_{1}$, see
Fig. 1(b). Besides, in the perfect insulation condition, if we take for
example $E_{3}/E_{1}=1$, Eq.(5) tells us that it is impossible, starting from
finite and positive temperatures, to cool qubit 1 arbitrarily close to
absolute zero. Indeed, according to Eq.(5), the smallest possible temperature
for $E_{3}/E_{1}=1$ is obtained when $T_{c}/T_{h}$ is very small, such that
the final temperature of qubit 1 is $T_{1}=T_{c}/2$. Nevertheless, for
negative temperatures and in the condition of perfect insulation where Eq.(6)
applies, even if we let $E_{3}/E_{1}=1$ we can still make the ratio
$T_{c}/T_{h}$ very large, thus cooling qubit 1 as close to absolute zero as we
want.

Figure 4: (a) Lowest values for $T_{1}$ obtained when either $T_{h}>0$ or
$T_{h}<0$ for some values of $T_{c}$ (b) Percentage of cooling using both
positive and negative temperatures for the pumping heat at $T_{h}$ and taking
$T_{c}$ as reference. According to our numerical calculations, AR stops
working if $T_{c}<0.48$ ($T_{c}<0.0275$) when $T_{h}>0$ ($T_{h}<0$).

##  IV Conclusion

In this work we consider the autonomous quantum refrigerator making use of
three qubits, one qubit to be cooled and the other two playing the role
analogous to the spiral and the engine in a classical refrigerator, as
proposed in Ref. (Linden  _et al._ , 2010). Making a comparative study of the
functioning of this refrigerator in environments with either positive
temperature or negative temperatures, we show that under the same conditions,
that is, with the same parameters that characterize the autonomous quantum
refrigerator, the use of negative temperature brings advantages in relation to
the use restricted to positive temperatures. According to our numerical
simulations, negative temperatures allow a greater cooling range, surpassing
the positive temperature operation by about 20% for some parameters, see Fig.
4(a). Furthermore, while restricting the environments at positive
temperatures, the steady-state of qubit 1, which ends cooled at temperature
$T_{1}$, depends on the temperature $T_{h}$ of the hot reservoir, see Fig.
2(a), by using the hot environment at negative temperature the cooling of
qubit 1 occurs for all $T_{h}$, see Fig. 3(a). Also, by imposing the condition
of perfect insulation, that is, isolating the qubit 1 to be cooled from its
environment and letting it interact directly with the other qubits, we show
that, as in the case of positive temperatures, the qubit 1 can be cooled to a
temperature $T_{1}$ arbitrarily close to absolute zero. However, even for
perfect insulation, it is still possible to demonstrate an advantage of using
negative temperatures. In fact, for $T_{h}>0$, it is necessary to greatly
increase the ratio between the energy $E_{3}$ of the third qubit and the
energy $E_{1}$ of the qubit 1 to be cooled to bring its temperature $T_{1}$ to
values close to absolute zero. On the other hand, using $T_{h}<0$ is possible
to set the ratio $E_{3}/E_{1}$ to reasonably small values but let the ratio
$T_{c}/\left|T_{h}\right|$ very large, such that the temperature of qubit 1 to
be cooled approaches absolute zero. To summarize, using $T_{h}<0$ to cool a
qubit is more powerful than $T_{h}>0$.

###### Acknowledgements.

We acknowledge financial support from the Brazilian agencies: CAPES, financial
code 001, CNPq, FAPEG, and FAPESP. This work was performed as part of the
Brazilian National Institute of Science and Technology (INCT) for Quantum
Information Grant No. 465469/2014-0.

## References

  * Linden  _et al._ (2010) N. Linden, S. Popescu, and P. Skrzypczyk, Phys. Rev. Lett. 105, 130401 (2010).
  * Purcell and Pound (1951) E. M. Purcell and R. V. Pound, Phys. Rev. 81, 279 (1951).
  * Ramsey (1956) N. F. Ramsey, Phys. Rev. 103, 20 (1956).
  * Carr (2013) L. D. Carr, Science 339, 42 (2013).
  * Braun  _et al._ (2013) S. Braun, J. P. Ronzheimer, M. Schreiber, S. S. Hodgman, T. Rom, I. Bloch, and U. Schneider, Science 339, 52 (2013).
  * Dunkel and Hilbert (2013) J. Dunkel and S. Hilbert, Nature Physics 10, 67 (2013), article.
  * Sokolov (2013) I. M. Sokolov, Nature Physics 10, 7 (2013).
  * Hilbert  _et al._ (2014) S. Hilbert, P. Hänggi, and J. Dunkel, Phys. Rev. E 90, 062116 (2014).
  * Campisi (2015) M. Campisi, Phys. Rev. E 91, 052147 (2015).
  * Cerino  _et al._ (2015) L. Cerino, A. Puglisi, and A. Vulpiani, Journal of Statistical Mechanics: Theory and Experiment 2015, P12002 (2015).
  * Frenkel and Warren (2015) D. Frenkel and P. B. Warren, American Journal of Physics 83, 163 (2015), https://doi.org/10.1119/1.4895828 .
  * Swendsen and Wang (2015) R. H. Swendsen and J.-S. Wang, Phys. Rev. E 92, 020103 (2015).
  * Poulter (2016) J. Poulter, Phys. Rev. E 93, 032149 (2016).
  * Swendsen and Wang (2016) R. H. Swendsen and J.-S. Wang, Physica A: Statistical Mechanics and its Applications 453, 24 (2016).
  * Abraham and Penrose (2017) E. Abraham and O. Penrose, Phys. Rev. E 95, 012125 (2017).
  * Hama  _et al._ (2018) Y. Hama, W. J. Munro, and K. Nemoto, Phys. Rev. Lett. 120, 060403 (2018).
  * Struchtrup (2018) H. Struchtrup, Phys. Rev. Lett. 120, 250602 (2018).
  * de Assis  _et al._ (2019) R. J. de Assis, T. M. de Mendonça, C. J. Villas-Boas, A. M. de Souza, R. S. Sarthour, I. S. Oliveira, and N. G. de Almeida, Phys. Rev. Lett. 122, 240602 (2019).
  * Mendonça  _et al._ (2020) T. M. Mendonça, A. M. Souza, R. J. de Assis, N. G. de Almeida, R. S. Sarthour, I. S. Oliveira, and C. J. Villas-Boas, Phys. Rev. Research 2, 043419 (2020).
  * Carmichael (1999) H. J. Carmichael,  _Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations_ (Springer Science $\&$ Business Media, 1999).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Artacho and Falicov (1993) E. Artacho and L. M. Falicov, Phys. Rev. B 47, 1190 (1993).
  * Ghosh  _et al._ (2012) A. Ghosh, S. S. Sinha, and D. S. Ray, Phys. Rev. E 86, 011138 (2012).
  * Johansson  _et al._ (2012) J. Johansson, P. Nation, and F. Nori, Computer Physics Communications 183, 1760 (2012).
  * Johansson  _et al._ (2013) J. Johansson, P. Nation, and F. Nori, Computer Physics Communications 184, 1234 (2013).


# Bijective Enumeration and Sign-Imbalance for
Permutation Depth and Excedances

Sen-Peng Eu Department of Mathematics
National Taiwan Normal University
Taiwan, ROCChinese Air Force Academy
Taiwan, ROC<EMAIL_ADDRESS>Department of Applied Mathematics
National Pingtung University
Taiwan, ROC Tung-Shan Fu Yuan-Hsun Lo Department of Applied Mathematics
National Pingtung University
Taiwan, ROC<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a simplified variant of Biane’s bijection between permutations and
3-colored Motzkin paths with weight that keeps track of the inversion number,
excedance number and a statistic so-called depth of a permutation. This
generalizes a result by Guay-Paquet and Petersen about a continued fraction of
the generating function for depth on the symmetric group $\mathfrak{S}_{n}$ of
permutations. In terms of weighted Motzkin path, we establish an involution on
$\mathfrak{S}_{n}$ that reverses the parities of depth and excedance numbers
simultaneously, which proves that the numbers of permutations with even and
odd depth (excedance numbers, respectively) are equal if $n$ is even and
differ by the tangent number if $n$ is odd. Moreover, we present some
interesting sign-imbalance results on permutations and derangements, refined
with respect to depth and excedance numbers.

##  1 Introduction

###  1.1 Preliminaries and background

Let $\mathfrak{S}_{n}$ be the symmetric group of permutations on
$[n]:=\\{1,2,\dots,n\\}$. Given a permutation
$\sigma=\sigma(1)\cdots\sigma(n)\in\mathfrak{S}_{n}$, an _inversion_ of
$\sigma$ is a pair $(i,j)$ of indices such that $1\leq i<j\leq n$ and
$\sigma(i)>\sigma(j)$. The _inversion number_ of $\sigma$, denoted by
${\textsf{inv}}(\sigma)$, is defined to be the number of inversions of
$\sigma$. Petersen and Tenner [12] presented a statistic called _depth_ for a
Coxeter group, which is defined in terms of factorizations of the elements
into reflections. In the case of symmetric group, the reflections of
$\mathfrak{S}_{n}$ are the transpositions $(i\,j)$, $1\leq i<j\leq n$, and the
depth of a permutation $\sigma\in\mathfrak{S}_{n}$, denoted by $0pt(\sigma)$,
is defined as

$0pt(\sigma)=\min\left\\{\sum_{r=1}^{k}(j_{r}-i_{r}):\sigma=(i_{1}\,j_{1})(i_{2}\,j_{2})\cdots(i_{k}\,j_{k})\right\\}.$
(1.1)

They obtained a simple formula for the calculation of depth

$0pt(\sigma)=\sum_{\sigma(i)>i}\big{(}\sigma(i)-i\big{)},$ (1.2)

which turns out to be one half of the _total displacement_ of $\sigma$ [10,
Problem 5.1.1.28] (also called _Spearman’s disarray_ for $\mathfrak{S}_{n}$ in
[4]).

For any two sequences $\\{\gamma_{h}\\}_{h\geq 0}$ and
$\\{\lambda_{h}\\}_{h\geq 1}$, let $\mathfrak{F}(\gamma_{h},\lambda_{h})$
denote the series in $z$ defined by

$\mathfrak{F}(\gamma_{h},\lambda_{h})=\dfrac{1}{1-\gamma_{0}z-\dfrac{\lambda_{1}z^{2}}{1-\gamma_{1}z-\dfrac{\lambda_{2}z^{2}}{1-\gamma_{2}z-\dfrac{\lambda_{3}z^{2}}{\ddots}}}}$
(1.3)

These are Jacobi-type continued fractions (also known as $J$-_fractions_).
Flajolet [7] gave an combinatorial interpretation of (1.3) as the generating
function for weighted Motzkin paths.

A _Motzkin path_ of length $n$ is a lattice path from the origin to the point
$(n,0)$ staying weakly above the $x$-axis, using the _up step_ $(1,1)$, _down
step_ $(1,-1)$ and _horizontal step_ $(1,0)$. Let $U$, $D$ and $H$ denote an
up step, a down step and a horizontal step, respectively. Let
${\mathcal{M}}_{n}$ denote the set of Motzkin paths of length $n$.

For a Motzkin path $\mu=x_{1}\cdots x_{n}\in{\mathcal{M}}_{n}$, the _height_
of each step $x_{j}$ is the maximum $y$-coordinate achieved by $x_{j}$. For
$x\in\\{U,D,H\\}$ and a nonnegative integer $h$, let $x^{(h)}$ denote a step
$x$ at height $h$. The _weight_ $\omega(\mu)$ of $\mu$ is defined to be the
product of the weight $\omega(x_{j})$ of each step $x_{j}$ for all $j\in[n]$.
By the _area_ under $\mu$, denoted by ${\textsf{area}}(\mu)$, we mean the area
of the region enclosed by $\mu$ and the $x$-axis. Guay-Paquet and Petersen [9]
defined a weight of $\mu$ by setting

$\omega(x_{i})=\begin{cases}ht^{(2h-1)/2}&\mbox{if $x_{i}=U^{(h)}$ or
$D^{(h)}$;}\\\\[1.72218pt] (2h+1)t^{h}&\mbox{if $x_{i}=H^{(h)}$}\end{cases}$

and established a surjective map $\Phi:\sigma\mapsto\mu$ of $\mathfrak{S}_{n}$
onto ${\mathcal{M}}_{n}$ with $0pt(\sigma)={\textsf{area}}(\mu)$, which yields
the following generating function for permutations with respect to depth:

$\sum_{n\geq
0}\sum_{\sigma\in\mathfrak{S}_{n}}t^{0pt(\sigma)}z^{n}=\mathfrak{F}(\gamma_{h},\lambda_{h})$
(1.4)

with $\lambda_{h}=h^{2}t^{2h-1}$ and $\gamma_{h}=(2h+1)t^{h}$, i.e.,

$\sum_{n\geq
0}\sum_{\sigma\in\mathfrak{S}_{n}}t^{0pt(\sigma)}z^{n}=\dfrac{1}{1-z-\dfrac{tz^{2}}{1-3tz-\dfrac{4t^{3}z^{2}}{1-5t^{2}z-\dfrac{9t^{5}z^{2}}{1-7t^{3}z-\dfrac{16t^{7}z^{2}}{1-\cdots}}}}}.$

Guay-Paquet and Petersen [9] commented that map
$\Phi:\mathfrak{S}_{n}\rightarrow{\mathcal{M}}_{n}$ is due to Foata and
Zeilberger [8] and that the statistic $0pt$ coincides with a statistic,
$\mathsf{Edif}$, studied by Clarke, Steingrímsson and Zeng [3].

Following [14], a _3-colored Motzkin path_ is a Motzkin path with three kinds
of horizontal steps, denoted by $H_{1}$, $H_{2}$ and $H_{3}$. Using a double
labeling scheme, Biane [2] established a bijection between permutations and
3-colored Motzkin paths with weight, which keeps track of the number of
inversions. We refer to an alternative version of this bijection given by
Sokal and Zeng in [14, Section 6.2]. Elizalde [6] used this bijection to study
the joint distribution of multiple statistics on $\mathfrak{S}_{n}$. Using
Biane’s method, we propose a multivariate weight function on the steps of a
3-colored Motzkin path and turn the surjective map
$\Phi:\mathfrak{S}_{n}\rightarrow{\mathcal{M}}_{n}$ of Guay-Paquet and
Petersen into a bijection that gives a generalization of (1.4). (This resolves
a conjecture raised by Petersen in an unpublished note [11].)

###  1.2 Main results

Let ${\mathcal{W}}_{n}$ denote the set of 3-colored Motzkin paths of length
$n$ with a weight function $\omega$ defined by setting $\omega(H_{3}^{(0)})=p$
and

$\left\\{\begin{aligned}
\omega(U^{(h)})&\in\\{st^{2h-1},st^{2h-1}q,\dots,st^{2h-1}q^{h-1}\\};\\\
\omega(D^{(h)})&\in\\{q^{2h-1},q^{2h},\dots,q^{3h-2}\\};\\\
\omega(H^{(h)}_{1})&\in\\{st^{h}q^{h},st^{h}q^{h+1},\dots,st^{h}q^{2h-1}\\};\\\
\omega(H^{(h)}_{2})&\in\\{t^{h}q^{h},t^{h}q^{h+1},\dots,t^{h}q^{2h-1}\\};\\\
\omega(H^{(h)}_{3})&=pt^{h}q^{2h}\end{aligned}\right.$ (1.5)

for all $h\geq 1$.

Let $\sigma=\sigma(1)\cdots\sigma(n)\in\mathfrak{S}_{n}$. We say that an index
$i\in[n]$ is an _excedance_ if $\sigma(i)>i$ and a _fixed point_ if
$\sigma(i)=i$. Let ${\textsf{exc}}(\sigma)$ and ${\textsf{fix}}(\sigma)$
denote the number of excedances and fixed points of $\sigma$, respectively.
For all positive integers $k$, we use the notations for $q$-integer
$[k]_{q}=1+q+\cdots+q^{k-1}$ and $[0]_{q}=0$. One of our main results is the
following bijection.

######  Theorem 1.1.

There is a bijection $\Psi:\mathfrak{S}_{n}\rightarrow{\mathcal{W}}_{n}$ such
that a permutation $\sigma\in\mathfrak{S}_{n}$ with $i$ inversions, $j$ fixed
points, $k$ excedances and a depth of $\ell$ is carried to a 3-colored Motzkin
path $\mu\in{\mathcal{W}}_{n}$ with weight
$\omega(\mu)=q^{i}p^{j}s^{k}t^{\ell}$. Therefore, we have

$\sum_{n\geq
0}\left(\sum_{\sigma\in\mathfrak{S}_{n}}q^{{\textsf{inv}}(\sigma)}p^{{\textsf{fix}}(\sigma)}s^{{\textsf{exc}}(\sigma)}t^{0pt(\sigma)}\right)z^{n}=\mathfrak{F}(\gamma_{h},\lambda_{h})$
(1.6)

with $\lambda_{h}=s[h]_{q}^{2}(qt)^{2h-1}$ and
$\gamma_{h}=\big{(}(1+s)[h]_{q}+pq^{h}\big{)}(qt)^{h}$.

Remarks. By the weight given in (1.5), notice that a horizontal step
$x\in\\{H^{(h)}_{2},H^{(h)}_{3}\\}$ can be distinguished by the $q$-factor of
$\omega(x)$ for each $h\geq 1$. That is, if $\omega(x)$ contains $q^{h+d}$
then $x=H^{(h)}_{2}$ ($H^{(h)}_{3}$, respectively) if $0\leq d\leq h-1$
($d=h$, respectively). Thus, the bijection
$\Psi:\mathfrak{S}_{n}\rightarrow{\mathcal{W}}_{n}$ can still be established
with $H^{(h)}_{2},H^{(h)}_{3}$ combined and the parameter $p$ omitted.

The classical _Euler numbers_ $E_{n}$, defined by

$\displaystyle\sum_{n\geq 0}E_{n}\frac{z^{n}}{n!}$ $\displaystyle=\tan z+\sec
z$
$\displaystyle=1+z+\frac{z^{2}}{2!}+2\frac{z^{3}}{3!}+5\frac{z^{4}}{4!}+16\frac{z^{5}}{5!}+61\frac{z^{6}}{6!}+272\frac{z^{7}}{7!}+1385\frac{z^{8}}{8!}+\cdots,$

count the the number of _alternating permutations_ in $\mathfrak{S}_{n}$,
i.e., $\sigma\in\mathfrak{S}_{n}$ such that
$\sigma_{1}>\sigma_{2}<\sigma_{3}>\cdots\sigma_{n}$. The numbers $E_{2n}$ are
called the _secant numbers_ and the numbers $E_{2n+1}$ are called the _tangent
numbers_.

Based on the bijection $\Psi:\mathfrak{S}_{n}\rightarrow{\mathcal{W}}_{n}$, we
shall establish an involution on $\mathfrak{S}_{n}$ in terms of weighted
Motzkin paths, which reverses the parities of depth and the number of
excedances simultaneously. This proves that the numbers of permutations with
even and odd depth (excedance numbers, respectively) are equal if $n$ is even
and differ by $E_{n}$ ($E_{n}$ up to sign, respectively) if $n$ is odd.

######  Theorem 1.2.

There is an involution $\sigma\mapsto\sigma^{\prime}$ on $\mathfrak{S}_{n}$
satisfying

$0pt(\sigma)-0pt(\sigma^{\prime})={\textsf{exc}}(\sigma)-{\textsf{exc}}(\sigma^{\prime})={\textsf{inv}}(\sigma)-{\textsf{inv}}(\sigma^{\prime})\in\\{1,0,-1\\}$

and resulting in the following identities

  1. (i)

${\displaystyle\sum_{\sigma\in\mathfrak{S}_{n}}(-1)^{0pt(\sigma)}=\begin{cases}E_{n}&\mbox{for
$n$ odd}\\\ 0&\mbox{for $n$ even.}\end{cases}}$

  2. (ii)

${\displaystyle\sum_{\sigma\in\mathfrak{S}_{n}}(-1)^{{\textsf{exc}}(\sigma)}=\begin{cases}(-1)^{\frac{n-1}{2}}E_{n}&\mbox{for
$n$ odd}\\\ 0&\mbox{for $n$ even.}\end{cases}}$

Note that the identity in Theorem 1.2(ii) is a classical result of Euler. When
the fixed points of $\sigma$ are ignored, we present some interesting results
on the sign imbalances of permutations and derangements, refined with respect
to depth and excedance numbers.

######  Theorem 1.3.

For $n\geq 1$, we have

$\sum_{\sigma\in\mathfrak{S}_{n}}(-1)^{{\textsf{inv}}(\sigma)}s^{{\textsf{exc}}(\sigma)}t^{0pt(\sigma)}=(1-st)^{n-1}.$
(1.7)

A _derangement_ of size $n$ is a permutation in $\mathfrak{S}_{n}$ that
contains no fixed point. Let ${\mathcal{D}}_{n}\subset\mathfrak{S}_{n}$ be the
set of derangements of size $n$. Define

$F_{n}=F_{n}(s,t)=\sum_{\sigma\in{\mathcal{D}}_{n}}(-1)^{{\textsf{inv}}(\sigma)}s^{{\textsf{exc}}(\sigma)}t^{0pt(\sigma)}.$
(1.8)

Several of the initial polynomials $F_{n}(s,t)$ are listed below:

$\displaystyle F_{1}(s,t)$ $\displaystyle=0,$ $\displaystyle F_{2}(s,t)$
$\displaystyle=-st,$ $\displaystyle F_{3}(s,t)$ $\displaystyle=s(1+s)t^{2},$
$\displaystyle F_{4}(s,t)$ $\displaystyle=s^{2}t^{2}-s(1+s)^{2}t^{3},$
$\displaystyle F_{5}(s,t)$ $\displaystyle=-2s^{2}(1+s)t^{3}+s(1+s)^{3}t^{4}.$

Collected in powers of $t$, the coefficient $s$-polynomials of $F_{n}(s,t)$
for $2\leq n\leq 9$ are listed in Table 1. We obtain a neat expression for the
sign imbalance of the joint distribution of depth and excedance numbers over
${\mathcal{D}}_{n}$.

Table 1: The $s$-polynomial coefficients of $F_{n}(s,t)$ in powers of $t$ for $2\leq n\leq 9$. | $t$ | $t^{2}$ | $t^{3}$ | $t^{4}$ | $t^{5}$ | $t^{6}$ | $t^{7}$ | $t^{8}$
---|---|---|---|---|---|---|---|---
$F_{2}$ | $-s$ |  |  |  |  |  |  |
$F_{3}$ |  | $s(1+s)$ |  |  |  |  |  |
$F_{4}$ |  | $s^{2}$ | $-s(1+s)^{2}$ |  |  |  |  |
$F_{5}$ |  |  | $-2s^{2}(1+s)$ | $s(1+s)^{3}$ |  |  |  |
$F_{6}$ |  |  | $-s^{3}$ | $3s^{2}(1+s)^{2}$ | $-s(1+s)^{4}$ |  |  |
$F_{7}$ |  |  |  | $3s^{3}(1+s)$ | $-4s^{2}(1+s)^{3}$ | $s(1+s)^{5}$ |  |
$F_{8}$ |  |  |  | $s^{4}$ | $-6s^{3}(1+s)^{2}$ | $5s^{2}(1+s)^{4}$ | $-s(1+s)^{6}$ |
$F_{9}$ |  |  |  |  | $-4s^{4}(1+s)$ | $10s^{3}(1+s)^{3}$ | $-6s^{2}(1+s)^{5}$ | $s(1+s)^{7}$

######  Theorem 1.4.

We have

$\sum_{n\geq
1}\left(\sum_{\sigma\in{\mathcal{D}}_{n}}(-1)^{{\textsf{inv}}(\sigma)}s^{{\textsf{exc}}(\sigma)}t^{0pt(\sigma)}\right)z^{n}=\sum_{k\geq
1}(-1)^{k}\left(\sum_{i=0}^{k-1}\binom{k-1}{i}s^{1+i}(1+s)^{k-1-i}z^{k+1+i}\right)t^{k}.$

## References

  * [1]
  * [2] P. Biane (1993): _Permutations suivant le type d’excédance et le nombre d’inversions et interprétation combinatoire d’une fraction continue de Heine_.  European J. Combin. 14(4), pp. 277–284, 10.1006/eujc.1993.1031.
  * [3] R.J. Clarke, E. Steingrímsson & J. Zeng (1997): _New Euler–Mahonian statistics on permutations and words_.  Adv. Appl. Math. 18(3), pp. 237–270, 10.1006/aama.1996.0506.
  * [4] P. Diaconis & R.L. Graham (1977): _Spearman’s footrule as a measure of disarray_.  J. Roy. Statist. Soc. Ser. B 39(2), pp. 262–268, 10.1111/j.2517-6161.1977.tb01624.x.
  * [5] D. Dumont (1995): _Further triangles of Seidel–Arnold type and continued fractions related to Euler and Springer numbers_.  Adv. Appl. Math. 16(3), pp. 275–296, 10.1006/aama.1995.1014.
  * [6] S. Elizalde (2018): _Continued fractions for permutation statistics_.  Discrete Math. Theor. Comput. Sci. 19(2), 10.23638/DMTCS-19-2-11.  arXiv:https://arxiv.org/abs/1703.08742.
  * [7] P. Flajolet (1980): _Combinatorial aspects of continued fractions_.  Discrete Math. 32(2), pp. 125–161, 10.1016/0012-365X(80)90050-3.
  * [8] D. Foata & D. Zeilberger (1990): _Denert’s permutation statistic is indeed Euler–Mahonian_.  Stud. Appl. Math. 83(1), pp. 31–59, 10.1002/sapm199083131.
  * [9] M. Guay-Paquet & T.K. Petersen (2014): _The generating function for total displacement_.  Electron. J. Combin. 21(3):P3.37, 10.37236/4329.
  * [10] D.E. Knuth (1998): _The Art of Computer Programming, vol. 3_ , second edition.  Addison Wesley Longman Publishing Co., Inc., USA.
  * [11] T.K. Petersen: _A note for Doron: the generation function for total displacement (Spearman’s footrule) and inversion_.  Available at https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/noga12yleFeedback.pdf.
  * [12] T.K. Petersen & B.E. Tenner (2015): _The depth of a permutation_.  J. Combin. 6(1–2), pp. 145–178, 10.4310/JOC.2015.v6.n1.a9.
  * [13] N.J.A. Sloane: _The On-Line Encyclopedia of Integer Sequences_.  Available at http://oeis.org.
  * [14] A.D. Sokal & J. Zeng (2022): _Some multivariate master polynomials for permutations, set partitions, and perfect matchings, and their continued fractions_.  Adv. Appl. Math. 138:102341, 10.1016/j.aam.2022.102341.