# Chaos and Regularity in the Double Pendulum
with Lagrangian Descriptors

Javier Jiménez López<EMAIL_ADDRESS>Departamento de Física y
Matemáticas, Facultad de Ciencias, Universidad de Alcalá, 28805 Alcalá de
Henares, Madrid, Spain.  V. J. García-Garrido<EMAIL_ADDRESS>(Coresponding author). Departamento de Física y Matemáticas, Facultad de
Ciencias, Universidad de Alcalá, 28805 Alcalá de Henares, Madrid, Spain.

###### Abstract

In this paper we apply the method of Lagrangian descriptors as an indicator to
study the chaotic and regular behavior of trajectories in the phase space of
the classical double pendulum system. In order to successfully quantify the
degree of chaos with this tool, we first derive Hamilton’s equations of motion
for the problem in non-dimensional form, showing that they can be written
compactly using matrix algebra. Once the dynamical equations are obtained, we
carry out a parametric study in terms of the system’s total energy and the
other model parameters (lengths and masses of the pendulums, and gravity), to
determine the extent of the chaotic and regular regions in the phase space.
Our numerical results show that for a given mass ratio, the maximum chaotic
fraction of phase space trajectories is attained when the pendulums have equal
lengths. Moreover, we give a characterization of the growth and decay of chaos
in the system in terms of the model parameters, and explore the hypothesis
that the chaotic fraction follows an exponential law over different energy
regimes.

keywords: Hamiltonian systems, Double pendulum, Chaos indicators, Lagrangian
descriptors, Chaos and Regularity.

##  1 Introduction

The study of dynamical systems and their chaotic behavior has captured the
attention of the scientific community for decades given its prevalence in
natural phenomena and its relevance in technological applications [25]. Among
these systems, the double pendulum stands out as a classical and fascinating
model, not only for its conceptual simplicity but also for the complexity of
its dynamical behavior. This system is a perfect example of a nonlinear
system, which, despite its apparent mechanical simplicity, exhibits a great
variety of dynamical behaviors ranging from periodic motions to chaos.

The purpose of this study is to analyze the fraction of chaotic trajectories
in the phase space of the double pendulum, which is a paradigmatic benchmark
example of chaos in the nonlinear sciences literature [11]. In order to carry
out this task, we first write the model equations in dimensionless form, and
provide a simplified expression for them using matrix algebra. We follow a
similar approach as that developed for other mechanical systems such as, for
example, the Furuta pendulum [18]. This formulation allows for a more detailed
exploration of the degree of chaos present in the system depending on its
characteristic parameters. In a recent article concerned with coupled
pendulums [24], the authors highlight the importance of carrying out complete
parametric studies for such systems, by means of applying different
techniques, as a fundamental step towards obtaining relevant insights about
their rich and complex dynamical behavior. In the problem we address in this
paper, we use chaos indicators derived from the method of Lagrangian
descriptors [8, 26] to quantify chaos in the double pendulum. This trajectory-
based diagnostic tool has proved to be particularly effective for the analysis
of regular and chaotic dynamics, and also for revealing the underlying phase
space structures, in many physical systems [6, 19, 4, 20]. This methodology,
which is straightforward to implement since it only requires the integration
of ensembles of trajectories, has the capability of reducing the computational
time required to study the invariant manifolds that divide the phase space
into regions where trajectories display qualitatively distinct dynamical
behaviors.

The relevance of our study lies not only in its contribution to the
understanding of the double pendulum, but also to the analysis of similar
Hamiltonian systems such as those arising in the motion of celestial bodies
[10], which share fundamental characteristics with the double pendulum. These
connections make the findings of our study potentially applicable in a broader
context. In the existing literature, several approaches have been undertaken
to study chaos in the double pendulum, ranging from experimental works as
those presented in [12, 21] to numerical studies such as [9, 23, 17, 5].
However, the application of Lagrangian descriptors in this problem offers a
new perspective that combines both mathematical rigor with physical intuition.
Recently, a preliminary study of the chaotic fraction in the phase space of
the double pendulum has been performed in [2] by means of the Maximum Lyapunov
Exponent. However, in this paper the authors only address the case where the
masses and lengths of the pendulums are equal to one, and gravity is taken as
that of the Earth. Our work goes further and extends the results obtained in
[2]. We do so by first making the double pendulum model dimensionless and use
matrix algebra to simplify the equations of motion. Once Hamilton’s equations
are written in terms of the mass and length ratios, we then conduct a detailed
parametric study to determine how the masses and lengths of the pendulums
influence the degree of chaos that the system exhibits.

This article is outlined as follows. In Section 2 we develop the non-
dimensional Hamiltonian model of the double pendulum and investigate the
critical points of its potential energy surface. Section 3 is devoted to the
methodology of this work. We explain how one can define different chaos
indicators derived from Lagrangian descriptors in order to study the chaotic
fraction of the phase space as a function of the energy and the model
parameters. In Section 4 we discuss and interpret the main results obtained
from our analysis of the double pendulum system. We finish this work with
Section 5, where we summarize the conclusions and provide some ideas for
future research.

##  2 Hamiltonian Formulation of the Mechanical System

We begin this section by deriving the Hamiltonian formulation and Hamilton’s
equations of motion for the double pendulum. Despite this is a well-known
academic exercise, carried out in many courses and textbooks on Analytical
Mechanics, see e.g. [7], we explain here how to obtain a compact and
dimensionless expression of the problem, written in terms of matrix algebra.
This approach is inspired by the work developed in [18] for the Furuta
pendulum, which has similarities with our model system.

Figure 1: Diagram of the double pendulum system, taken from [1]

The first step is to define the configuration space vector
$\boldsymbol{\theta}=(\theta_{1},\theta_{2})^{T}$ that contains the degrees of
freedom of the mechanical system. If we refer to the diagram shown in Fig. 1,
we have that the position vectors of the masses $m_{1}$ and $m_{2}$ for the
double pendulum system are given by:

$\mathbf{r}_{1}=l_{1}\begin{bmatrix}\sin\theta_{1}\\\
-\cos\theta_{1}\end{bmatrix}\quad,\quad\mathbf{r}_{2}=\mathbf{r}_{1}+\mathbf{r}_{12}=l_{1}\begin{bmatrix}\sin\theta_{1}\\\
-\cos\theta_{1}\end{bmatrix}+l_{2}\begin{bmatrix}\sin\theta_{2}\\\
-\cos\theta_{2}\end{bmatrix}\;,$ (1)

and the velocity vectors:

$\dot{\mathbf{r}}_{1}=l_{1}\dot{\theta}_{1}\begin{bmatrix}\cos\theta_{1}\\\
\sin\theta_{1}\end{bmatrix}\quad,\quad\dot{\mathbf{r}}_{2}=\dot{\mathbf{r}}_{1}+\dot{\mathbf{r}}_{12}=l_{1}\dot{\theta}_{1}\begin{bmatrix}\cos\theta_{1}\\\
\sin\theta_{1}\end{bmatrix}+l_{2}\dot{\theta}_{2}\begin{bmatrix}\cos\theta_{2}\\\
\sin\theta_{2}\end{bmatrix}\;.$ (2)

The kinetic energy of the system is:

$T=\dfrac{1}{2}m_{1}\lVert\dot{\mathbf{r}}_{1}\rVert^{2}+\dfrac{1}{2}m_{2}\lVert\dot{\mathbf{r}}_{2}\rVert^{2}\;.$
(3)

Now, we have that:

$\lVert\dot{\mathbf{r}}_{1}\rVert^{2}=\dot{\mathbf{r}}_{1}\cdot\dot{\mathbf{r}}_{1}=l_{1}^{2}\dot{\theta_{1}}^{2}=\dot{\boldsymbol{\theta}}^{T}\begin{bmatrix}l_{1}^{2}&0\\\
0&0\end{bmatrix}\dot{\boldsymbol{\theta}}\;,$ (4)

where $\dot{\boldsymbol{\theta}}=(\dot{\theta}_{1},\dot{\theta}_{2})^{T}$ is
the vector of generalized velocities. The degrees of freedom of the system are
$\theta_{1}$ and $\theta_{2}$. Moreover:

$\lVert\dot{\mathbf{r}}_{2}\rVert^{2}=\dot{\mathbf{r}}_{12}\cdot\dot{\mathbf{r}}_{12}=\lVert\dot{\mathbf{r}}_{1}\rVert^{2}+\lVert\dot{\mathbf{r}}_{12}\rVert^{2}+2\dot{\mathbf{r}}_{1}\cdot\dot{\mathbf{r}}_{12}\;,$
(5)

and using the expressions:

$\lVert\dot{\mathbf{r}}_{12}\rVert^{2}=\dot{\mathbf{r}}_{12}\cdot\dot{\mathbf{r}}_{12}=l_{2}^{2}\dot{\theta_{2}}^{2}=\dot{\boldsymbol{\theta}}^{T}\\!\begin{bmatrix}0&0\\\
0&l_{2}^{2}\end{bmatrix}\dot{\boldsymbol{\theta}}\;,$ (6)
$2\dot{\mathbf{r}}_{1}\cdot\dot{\mathbf{r}}_{12}=2l_{1}l_{2}\dot{\theta}_{1}\dot{\theta}_{2}\cos\Delta\theta=\dot{\boldsymbol{\theta}}^{T}\\!\begin{bmatrix}0&l_{1}l_{2}\cos\Delta\theta\\\
l_{1}l_{2}\cos\Delta\theta&0\end{bmatrix}\dot{\boldsymbol{\theta}}\;,$ (7)

where $\Delta\theta=\theta_{1}-\theta_{2}$, we arrive at:

$T=\dfrac{1}{2}\dot{\boldsymbol{\theta}}^{T}\\!\begin{bmatrix}(m_{1}+m_{2})l_{1}^{2}&m_{2}l_{1}l_{2}\cos\Delta\theta\\\
m_{2}l_{1}l_{2}\cos\Delta\theta&m_{2}l_{2}^{2}\end{bmatrix}\dot{\boldsymbol{\theta}}\;.$
(8)

This can be rewritten as:

$T=\dfrac{1}{2}m_{2}l_{2}^{2}\,\dot{\boldsymbol{\theta}}^{T}\\!\begin{bmatrix}\beta&\alpha\cos\Delta\theta\\\
\alpha\cos\Delta\theta&1\end{bmatrix}\dot{\boldsymbol{\theta}}=\dfrac{1}{2}m_{2}l_{2}^{2}\dot{\boldsymbol{\theta}}^{T}B(\cos\Delta\theta)\dot{\boldsymbol{\theta}}\;,$
(9)

where we have defined the matrix:

$B(x)=\begin{bmatrix}\beta&\alpha x\\\ \alpha x&1\end{bmatrix}\;,$ (10)

and introduced the dimensionless quantities:

$\alpha=\dfrac{l_{1}}{l_{2}}>0\quad,\quad\beta=\dfrac{m_{1}+m_{2}}{m_{2}}\alpha^{2}=\dfrac{\alpha^{2}}{\mu}>0\quad,\quad\mu=\dfrac{m_{2}}{m_{1}+m_{2}}\;.$
(11)

If we set $\sigma=\dfrac{m_{1}}{m_{2}}$ then:

$\mu=\dfrac{1}{1+\sigma}\quad,\quad\beta={\alpha^{2}}(1+\sigma)\;.$

Notice that by construction, the model parameters satisfy the conditions
$\beta>\alpha^{2}$, $0<\mu<1$ and $0<\sigma<+\infty$. Using this simplified
expression for the kinetic energy, the Lagrangian of the system is:

$\mathcal{L}=T-V=\dfrac{1}{2}m_{2}l_{2}^{2}\,\dot{\boldsymbol{\theta}}^{T}\\!B(\cos\Delta\theta)\,\dot{\boldsymbol{\theta}}+(m_{1}+m_{2})\,l_{1}g\cos\theta_{1}+m_{2}gl_{2}\cos\theta_{2}\;.$
(12)

If we non-dimensionalize time by setting $t=\sqrt{\dfrac{l_{2}}{g}}\,\tau$
then:

$\dot{\boldsymbol{\theta}}=\dfrac{d\boldsymbol{\theta}}{dt}=\dfrac{d\boldsymbol{\theta}}{d\tau}\dfrac{d\tau}{dt}=\sqrt{\dfrac{g}{l_{2}}}\dfrac{d\boldsymbol{\theta}}{d\tau}=\sqrt{\dfrac{g}{l_{2}}}\,\boldsymbol{\theta}^{\prime}\;.$
(13)

This yields the non-dimensional kinetic energy:

$\mathcal{T}=\dfrac{T}{m_{2}l_{2}g}=\dfrac{1}{2}\left(\boldsymbol{\theta}^{\prime}\right)^{T}\\!\\!B(\cos\Delta\theta)\,\boldsymbol{\theta}^{\prime}\;.$
(14)

Similarly, the potential energy of the system is:

$V=-(m_{1}+m_{2})l_{1}g\cos\theta_{1}-m_{2}gl_{2}\cos\theta_{2}=-m_{2}l_{2}g\left[\dfrac{\beta}{\alpha}\cos\theta_{1}+\cos\theta_{2}\right]\;,$
(15)

and if we make it non-dimensional:

$\mathcal{V}=\dfrac{V}{m_{2}l_{2}g}=-\dfrac{\beta}{\alpha}\cos\theta_{1}-\cos\theta_{2}\;.$
(16)

Hence, the dimensionless lagrangian is:

$\mathcal{L}=\dfrac{1}{2}\left(\boldsymbol{\theta}^{\prime}\right)^{T}\\!\\!B(\cos\Delta\theta)\,\boldsymbol{\theta}^{\prime}+\dfrac{\beta}{\alpha}\cos\theta_{1}+\cos\theta_{2}\;.$
(17)

In order to write the Hamiltonian, we need first the conjugate momenta:

$\mathbf{p}=\nabla_{\boldsymbol{\theta}^{\prime}}\mathcal{L}=B(\cos\Delta\theta)\,\boldsymbol{\theta}^{\prime}=(p_{1},p_{2})^{T}\;.$
(18)

Therefore:

$\boldsymbol{\theta}^{\prime}=B^{-1}(\cos\Delta\theta)\,\mathbf{p}\;,$ (19)

where:

$B^{-1}(x)=\dfrac{1}{\beta-\alpha^{2}x^{2}}\begin{bmatrix}1&-\alpha x\\\
-\alpha x&\beta\end{bmatrix}\;.$ (20)

The non-dimensional Hamiltonian is:

$\mathcal{H}(\boldsymbol{\theta},\mathbf{p})=\mathcal{T}(\boldsymbol{\theta},\mathbf{p})+\mathcal{V}\left(\boldsymbol{\theta}\right)=\dfrac{1}{2}\mathbf{p}^{T}B^{-1}(\cos\Delta\theta)\,\mathbf{p}-\dfrac{\beta}{\alpha}\cos\theta_{1}-\cos\theta_{2}\;,$
(21)

and Hamilton’s equations are given by the system:

$\begin{cases}\boldsymbol{\theta}^{\prime}=\nabla_{\mathbf{p}}\mathcal{H}=\nabla_{\mathbf{p}}\mathcal{T}=B^{-1}(\cos\Delta\theta)\,\mathbf{p}\\\\[2.84544pt]
\mathbf{p}^{\prime}=-\nabla_{\boldsymbol{\theta}}\mathcal{H}=-\nabla_{\boldsymbol{\theta}}\mathcal{T}-\nabla_{\boldsymbol{\theta}}\mathcal{V}\end{cases}\;.$
(22)

Now:

$\nabla_{\boldsymbol{\theta}}\mathcal{V}=\begin{bmatrix}\dfrac{\beta}{\alpha}\sin\theta_{1}\\\\[8.5359pt]
\sin\theta_{2}\end{bmatrix}\;,$ (23)

and to calculate $\nabla_{\boldsymbol{\theta}}\mathcal{T}$ we need the
derivative of the matrix $B^{-1}$. From the matrix identity:

$BB^{-1}=\mathbb{I}\;,$ (24)

we can differentiate this expression to obtain:

$\dfrac{dB}{dx}B^{-1}+B\dfrac{d}{dx}\left(B^{-1}\right)=\mathbb{O}\;,$ (25)

where $\mathbb{O}$ is the zero matrix. Therefore:

$\dfrac{d}{dx}\left(B^{-1}\right)=-B^{-1}\dfrac{dB}{dx}B^{-1}\;,$ (26)

and taking the derivative of Eq. (10):

$\dfrac{dB}{dx}=\begin{bmatrix}0&\alpha\\\ \alpha&0\end{bmatrix}\;.$ (27)

This procedure yields the result:

$C(x)=-B^{-1}\dfrac{dB}{dx}B^{-1}=-\dfrac{\alpha}{\left(\beta-\alpha^{2}x^{2}\right)^{2}}\begin{bmatrix}-2\alpha
x&\beta+\alpha^{2}x^{2}\\\ \beta+\alpha^{2}x^{2}&-2\alpha\beta
x\end{bmatrix}\;.$ (28)

which is a symmetric matrix. We arrive at:

$\nabla_{\boldsymbol{\theta}}\mathcal{T}=\dfrac{1}{2}\mathbf{p}^{T}C(\cos\Delta\theta)\,\mathbf{p}\,\nabla_{\boldsymbol{\theta}}\left(\cos\Delta\theta\right)\quad,\quad\nabla_{\boldsymbol{\theta}}\left(\cos\Delta\theta\right)=-\sin\Delta\theta\;\nabla_{\boldsymbol{\theta}}\left(\Delta\theta\right)\;,$
(29)

where $\nabla_{\boldsymbol{\theta}}\left(\Delta\theta\right)=(1,-1)^{T}$. The
resulting equations of motion are:

$\begin{cases}\boldsymbol{\theta}^{\prime}=B^{-1}(\cos\Delta\theta)\,\mathbf{p}\\\
\mathbf{p}^{\prime}=\dfrac{\sin\Delta\theta}{2}\mathbf{p}^{T}C(\cos\Delta\theta)\,\mathbf{p}\begin{bmatrix}1\\\
-1\end{bmatrix}-\begin{bmatrix}\dfrac{\beta}{\alpha}\sin\theta_{1}\\\\[2.84544pt]
\sin\theta_{2}\end{bmatrix}\end{cases}\;.$ (30)

If we introduce at this point the state vector
$\boldsymbol{\sigma}=(\boldsymbol{\theta},\mathbf{p})^{T}$, Hamilton’s
equations can be written in compact form as:

$\boldsymbol{\sigma}^{\prime}=\mathbf{F}(\boldsymbol{\sigma})\;,$ (31)

where $\mathbf{F}=(\mathbf{F}_{1},\mathbf{F}_{2})^{T}$ is the Hamiltonian
vector field whose components are defined by:

$\mathbf{F}_{1}(\boldsymbol{\sigma})=B^{-1}(\cos\Delta\theta)\,\mathbf{p}\quad,\quad\mathbf{F}_{2}(\boldsymbol{\sigma})=\dfrac{\sin\Delta\theta}{2}\mathbf{p}^{T}C(\cos\Delta\theta)\,\mathbf{p}\begin{bmatrix}1\\\
-1\end{bmatrix}-\begin{bmatrix}\dfrac{\beta}{\alpha}\sin\theta_{1}\\\\[5.69046pt]
\sin\theta_{2}\end{bmatrix}\;.$ (32)

The Jacobian matrix determining the variational equations of the linearized
dynamics about the phase space point
$\boldsymbol{\sigma}_{0}=(\boldsymbol{\theta}_{0},\mathbf{p}_{0})^{T}$ has the
block form:

$D_{\boldsymbol{\sigma}}\mathbf{F}(\boldsymbol{\sigma}_{0})=\begin{bmatrix}D_{\boldsymbol{\theta}}\mathbf{F}_{1}(\boldsymbol{\sigma}_{0})&D_{\mathbf{p}}\mathbf{F}_{1}(\boldsymbol{\sigma}_{0})\\\\[5.69046pt]
D_{\boldsymbol{\theta}}\mathbf{F}_{2}(\boldsymbol{\sigma}_{0})&D_{\mathbf{p}}\mathbf{F}_{2}(\boldsymbol{\sigma}_{0})\end{bmatrix}\;,$
(33)

where the subscript indicates the variables with respect to which we
differentiate to obtain each of the matrices involved in the jacobian. After
some algebra, the resulting matrix is:

$D_{\boldsymbol{\sigma}}\mathbf{F}=\begin{bmatrix}-\sin\Delta\theta\;C\left(\cos\Delta\theta\right)\left(\mathbf{p}\otimes\nabla_{\boldsymbol{\theta}}\left(\Delta\theta\right)\right)&B^{-1}(\cos\Delta\theta)\\\\[5.69046pt]
G(\boldsymbol{\theta},\mathbf{p})-H_{\mathcal{V}}&\sin\Delta\theta\left(\nabla_{\boldsymbol{\theta}}\left(\Delta\theta\right)\otimes\mathbf{p}\right)C\left(\cos\Delta\theta\right)\end{bmatrix}\;,$
(34)

such that $G(\boldsymbol{\theta},\mathbf{0})=\mathbb{O}$ and the hessian
matrix of the potential energy is:

$H_{\mathcal{V}}(\boldsymbol{\theta})=\begin{bmatrix}\dfrac{\beta}{\alpha}\cos\theta_{1}&0\\\\[2.84544pt]
0&\cos\theta_{2}\\\\[2.84544pt] \end{bmatrix}\;.$ (35)

The dynamical system in Eq. (30) has 4 equilibria, denoted by
$\boldsymbol{\sigma}_{\text{eq}}=(\theta_{1,\text{eq}},\theta_{2,\text{eq}},0,0)$,
located in configuration space ($\mathbf{p}=\mathbf{0}$) at the points:

$\theta_{1,\text{eq}}=0,\pi\quad,\quad\theta_{2,\text{eq}}=0,\pi\;.$ (36)

and we will denote them as follows:

$\boldsymbol{\sigma}^{1}_{\text{eq}}=(0,0,0,0)\;\;,\;\;\boldsymbol{\sigma}^{2}_{\text{eq}}=(0,\pi,0,0)\;\;,\;\;\boldsymbol{\sigma}^{3}_{\text{eq}}=(\pi,0,0,0)\;\;,\;\;\boldsymbol{\sigma}^{4}_{\text{eq}}=(\pi,\pi,0,0)\;,$
(37)

These points are critical points of the potential energy in Eq. (16), and in
order to get a better understanding of the geometrical features of the
potential energy surface, we have represented the function in Fig. 2, together
with the equilibria of Hamilton’s equations. The energies of the equilibrium
points are given by:

$\begin{cases}\mathcal{H}_{1}=\mathcal{V}_{1}(0,0)=-\alpha(1+\sigma)-1\\\\[8.5359pt]
\mathcal{H}_{2}=\mathcal{V}_{2}(0,\pi)=-\alpha(1+\sigma)+1\\\\[8.5359pt]
\mathcal{H}_{3}=\mathcal{V}_{3}(\pi,0)=\alpha(1+\sigma)-1\\\\[8.5359pt]
\mathcal{H}_{4}=\mathcal{V}_{4}(\pi,\pi)=\alpha(1+\sigma)+1\end{cases}\;.$
(38)

Notice that the energy levels of the equilibrium points of the system satisfy
the following relations, depending on the values of the model parameters. If
$\alpha(1+\sigma)>1$ then:

$\mathcal{H}_{1}<\mathcal{H}_{2}<\mathcal{H}_{3}<\mathcal{H}_{4}\;,$ (39)

whereas in the case $\alpha(1+\sigma)<1$ we have:

$\mathcal{H}_{1}<\mathcal{H}_{3}<\mathcal{H}_{2}<\mathcal{H}_{4}\;.$ (40)

In particular, there is a critical value of the parameters,
$\alpha(1+\sigma)=1$ for which $\mathcal{H}_{2}=\mathcal{H}_{3}$. This means
that, for $\beta=\alpha$, the saddles located at the configuration space
points $(0,\pi)$ and $(\pi,0)$ have the same energy.

Figure 2: Contours of the potential energy function given in (16) for the
double pendulum system with $\alpha=1$ and $\beta=2$. We have marked the
equilibrium points of Hamilton’s equations using magenta squares (saddles),
circles (local minimum) and diamonds (local maximum).

Next, we will carry out a stability analysis of the equilibria, with the
purpose of determining the local behavior of trajectories about these points.
First we define the constant $\gamma=\beta-\alpha^{2}$, and calculate the
Jacobian matrix evaluated at each equilibrium point of the system, which
yields:

$D_{\boldsymbol{\sigma}}\mathbf{F}(\boldsymbol{\sigma}^{1}_{\text{eq}})=\dfrac{1}{\gamma}\begin{bmatrix}0&0&1&-\alpha\\\\[2.84544pt]
0&0&-\alpha&\beta\\\\[2.84544pt]
-\dfrac{\beta\gamma}{\alpha}&0&0&0\\\\[7.11317pt]
0&-\gamma&0&0\end{bmatrix}\quad,\quad
D_{\boldsymbol{\sigma}}\mathbf{F}(\boldsymbol{\sigma}^{2}_{\text{eq}})=\dfrac{1}{\gamma}\begin{bmatrix}0&0&1&\alpha\\\\[2.84544pt]
0&0&\alpha&\beta\\\\[2.84544pt]
-\dfrac{\beta\gamma}{\alpha}&0&0&0\\\\[7.11317pt]
0&\gamma&0&0\end{bmatrix}\;,$ (41)
$D_{\boldsymbol{\sigma}}\mathbf{F}(\boldsymbol{\sigma}^{3}_{\text{eq}})=\dfrac{1}{\gamma}\begin{bmatrix}0&0&1&\alpha\\\\[2.84544pt]
0&0&\alpha&\beta\\\\[2.84544pt]
\dfrac{\beta\gamma}{\alpha}&0&0&0\\\\[7.11317pt]
0&-\gamma&0&0\end{bmatrix}\quad,\quad
D_{\boldsymbol{\sigma}}\mathbf{F}(\boldsymbol{\sigma}^{4}_{\text{eq}})=\dfrac{1}{\gamma}\begin{bmatrix}0&0&1&-\alpha\\\\[2.84544pt]
0&0&-\alpha&\beta\\\\[2.84544pt]
\dfrac{\beta\gamma}{\alpha}&0&0&0\\\\[7.11317pt] 0&\gamma&0&0\end{bmatrix}\;.$
(42)

Analyzing the eigenvalues of these matrices, it can be shown that
$\boldsymbol{\sigma}^{1}_{\text{eq}}$ is stable, with center$\times$center
stability, $\boldsymbol{\sigma}^{2}_{\text{eq}}$ and
$\boldsymbol{\sigma}^{3}_{\text{eq}}$ are index-1 saddles, that is, they have
saddle$\times$center stability. Finally, $\boldsymbol{\sigma}^{4}_{\text{eq}}$
is an index-2 saddle, with saddle$\times$saddle stability.

##  3 Methodology

In this section we describe the methodology that we have followed to quantify
the chaotic and regular fraction of the phase space in the double pendulum
system. There exists many chaos indicators in the literature to carry out this
task [22], however they typically require the solution of the variational
equations of motion together with Hamilton’s equations in order to classify a
trajectory as chaotic or regular. Moreover, the integration time for applying
these techniques is typically of the order $10^{6}$, which would require
symplectic integrators to preserve the energy, or the use of high order Runge-
Kutta numerical schemes. Here we explore the use of the method of Lagrangian
descriptors (LDs) to carry out the classification of ensembles of
trajectories. This new approach has been validated recently on continuous and
discrete dynamical systems by developing several chaos indicators based on the
values of the LD function [8, 26]. These diagnostics have been benchmarked
against the Small Alignment Index method and have shown a success rate for the
classification of trajectories of above $90\%$. Moreover, the LD chaos
indicators only require an integration time of about three orders of magnitude
less than the classical methods. Furthermore, the implementation of LDs is
straightforward, since one only needs to add an extra differential equation to
the system of Eqs. 30. This makes it a suitable tool for all practitioners of
Dynamical Systems Theory and offers an important advantage with respect to
computational time.

The method of Lagrangian descriptors is a trajectory-based diagnostic
technique that was originally developed in the field of Geophysics to analyze
lagrangian transport and mixing processes in the ocean and the atmosphere [14,
16]. A lagrangian descriptor is a scalar function constructed in the following
way. Given a dynamical system in the form:

$\dfrac{d\mathbf{x}}{dt}=\mathbf{f}(\mathbf{x},t)\;.$ (43)

we define a non-negative function
$\mathcal{F}(\mathbf{x}(t;\mathbf{x}_{0}),t)$ that depends on the initial
condition $\mathbf{x}_{0}$ at time $t=t_{0}$. To determine the LD scalar
field, denoted by $\mathcal{L}$, we set an integration time $\tau>0$ and
calculate:

$\mathcal{L}(\mathbf{x}_{0},t_{0},\tau)=\mathcal{L}_{f}(\mathbf{x}_{0},t_{0},\tau)+\mathcal{L}_{b}(\mathbf{x}_{0},t_{0},\tau)\;,$
(44)

where the forward ($\mathcal{L}_{f}$) and backward ($\mathcal{L}_{b}$)
components of the LD function are given by:

$\mathcal{L}_{f}(\mathbf{x}_{0},t_{0},\tau)=\int_{t_{0}}^{t_{0}+\tau}\mathcal{F}(\mathbf{x}(t;\mathbf{x}_{0}),t)\,dt\quad,\quad\mathcal{L}_{b}(\mathbf{x}_{0},t_{0},\tau)=\int_{t_{0}-\tau}^{t_{0}}\mathcal{F}(\mathbf{x}(t;\mathbf{x}_{0}),t)\,dt\;.$
(45)

Hence, the calculation of LDs involve accumulating the value of the function
$\mathcal{F}$ along the trajectory starting at $\mathbf{x}_{0}$ as it evolves
forward and backward in time. In the literature it has been shown rigorously
that the scalar field generated by the method has the capability of
identifying the invariant sets (equilibria, stable and unstable manifolds,
tori, periodic orbits, etc,) that characterize the dynamical behavior of
trajectories in the phase space of the system [15, 13].

In this work we will use the LD whose definition resembles the $p$-norm of
Functional Analysis [13]. Thus, we select as the non-negative function the
expression:

$\mathcal{F}(\mathbf{x},t)=\sum_{i=1}^{n}|f_{i}(\mathbf{x},t)|^{p}\;\;,\;\;0<p\leq
1\;,$ (46)

where $f_{i}$ is the $i$-th component of the vector field that determines the
dynamical system in Eq. 43. In particular we will set $p=1/2$, which is
commonly used because it yields a scalar field where the phase space
structures are nicely detected. Notice that the classification of a trajectory
as chaotic or regular is independent of the fact that we integrate it forward
or backward in time, since both directions must yield the same answer. Hence,
for our analysis we have decided to integrate only forward, focusing on the
output provided by the component $\mathcal{L}_{f}$. For the integration time,
we have chosen a value of $\tau=700$ which is close to the one used in
previous works [8, 26] for the computation of the LD chaos indicators. This
value has been adjusted by looking at several trajectories that are known a
priori to be chaotic and regular for the double pendulum, see Fig. 3. The
numerical method that we have used to solve Hamilton’s equations is an
adaptive Runge-Kutta scheme of order 8 (in particular, the Dormand-Prince 8(9)
algorithm) with an error tolerance of $10^{-8}$, which has achieved a
conservation of the energy with at least 4 significant digits for all the
cases considered in this work.

A) B)

Figure 3: A) Poincaré section in the phase space slice given by Eq. (47) for
the system with $\alpha=1$, $\sigma=1$ and $\mathcal{H}_{0}=20$. We have
marked with an orange point a chaotic initial condition, and with blue a
regular one; B) Chaos indicators defined in Eq. (49) calculated for the
initial conditions in panel A).

In order to study the chaotic and regular fraction of the phase space in the
double pendulum system, we will carry out our simulations on the Poincaré
section defined as follows:

$\Sigma\left(\mathcal{H}_{0}\right)=\left\\{(\theta_{1},\theta_{2},p_{1},p_{2})\in\mathbb{R}^{4}\;\Big{|}\;\mathcal{H}=\mathcal{H}_{0}\;,\;\theta_{2}=0\;,\;\beta
p_{2}-\alpha p_{1}\cos(\Delta\theta)>0\right\\}\;.$ (47)

Note that, given the symmetry of the system with respect to the angle
$\theta_{1}$, we will only consider positive values for this configuration
variable, that is, $\theta_{1}\in[0,\pi]$. To use the chaos indicators derived
from Lagrangian descriptors, we need to obtain the neighbors of the initial
condition we want to analyze. These neighbors are given by:

$\mathbf{y}_{i}^{\pm}=\mathbf{x}_{0}\pm\sigma_{i}\,\mathbf{e}_{i}\;\;,\;\;i=1,\ldots,n\;,$
(48)

where $\mathbf{e}_{i}$ is the $i$-th canonical basis vector in
$\mathbb{R}^{n}$, and $\sigma_{i}$ represents the distance between the central
point $\mathbf{x}_{0}$ and its neighbors on the grid. The value of $n$
represents the dimension of the space where the initial conditions are
selected. In particular, since our phase space slice is a plane, then $n=2$,
and hence each initial condition on the grid has 4 neighbors. For our
simulations we have chosen a value of $\sigma_{i}=10^{-4}$. Using these
points, one can construct different chaos indicators based on LDs. In this
paper we use the following four diagnostics:

$\begin{split}D_{L}^{n}(\mathbf{x}_{0})&=\dfrac{1}{2n\mathcal{L}_{f}(\mathbf{x}_{0})}\sum_{i=1}^{n}\big{|}\mathcal{L}_{f}(\mathbf{x}_{0})-\mathcal{L}_{f}\left(\mathbf{y}_{i}^{+}\right)\big{|}+\big{|}\mathcal{L}_{f}(\mathbf{x}_{0})-\mathcal{L}_{f}\left(\mathbf{y}_{i}^{-}\right)\big{|}\;,\\\
R_{L}^{n}(\mathbf{x}_{0})&=\Bigg{|}1-\dfrac{1}{2n\mathcal{L}_{f}(\mathbf{x}_{0})}\sum_{i=1}^{n}\mathcal{L}_{f}\left(\mathbf{y}_{i}^{+}\right)+\mathcal{L}_{f}\left(\mathbf{y}_{i}^{-}\right)\Bigg{|}\;,\\\\[2.84544pt]
C_{L}^{n}(\mathbf{x}_{0})&=\dfrac{1}{2n}\sum_{i=1}^{n}\dfrac{\big{|}\mathcal{L}_{f}\left(\mathbf{y}_{i}^{+}\right)-\mathcal{L}_{f}\left(\mathbf{y}_{i}^{-}\right)\big{|}}{\sigma_{i}}\;,\\\\[2.84544pt]
S_{L}^{n}(\mathbf{x}_{0})&=\dfrac{1}{n}\sum_{i=1}^{n}\dfrac{\big{|}\mathcal{L}_{f}\left(\mathbf{y}_{i}^{+}\right)-2\mathcal{L}_{f}(\mathbf{x}_{0})+\mathcal{L}_{f}\left(\mathbf{y}_{i}^{-}\right)\big{|}}{\sigma_{i}^{2}}\end{split}$
(49)

where $\mathcal{L}_{f}(\cdot)$ is the LD value calculated for an integration
time $\tau$. Notice that the chaotic or regular nature of a trajectory is
equivalently characterized if the we integrate forward or backward, so only
one of the components of the LD function is required to determine the chaos
indicators. If we set a tolerance $\varepsilon>0$ (that depends on the
dynamical system under study and its parameters, and also on the chaos
indicator chosen), we say that the initial condition $\mathbf{x}_{0}$ gives
rise to a chaotic trajectory when the logarithm of the value of the indicator
at that point is above $\varepsilon$. Otherwise, we classify the trajectory as
regular.

To determine the threshold corresponding to each energy and $\alpha$ and
$\sigma$ parameter values, we have developed an algorithm based on the
histograms generated using the chaos indicators in Eq. 49. Upon examining the
shape of these histograms, we have identified three potential scenarios: the
first is where the entire phase space is regular, leading to only one peak in
the histogram; the second scenario occurs when the phase space contains both
regular and chaotic regions, resulting in two peaks in the histogram, with the
left one corresponding to the regular fraction and the right one to the
chaotic fraction; and the third scenario is when the entire phase space is
chaotic, where we will again observe a single peak in the histogram. In the
first and third scenarios, determining the threshold is relatively
straightforward, but in the second scenario, we have decided to define the
threshold as the value of the chaos indicator corresponding to the lowest
column between the two peaks, which is the minimum value between the two
maximum values.

Even though the four chaos indicators we have considered have a very similar
performance, we have noticed that the $R^{n}_{L}$ indicator is more sensitive
to the initial conditions given to our algorithm and disagrees with the other
three indicators in some of the cases we have studied, while the other three
indicators converge to the same classification in the vast majority of cases,
being $D^{n}_{L}$ and $S^{n}_{L}$ the two best indicators for the cases we
have analyzed.

With this algorithm, we can automate the process of classifying initial
conditions in such a way that we can classify any initial condition by
comparing the value of the chaos indicator of the condition itself with the
threshold we have determined. In Fig. 4, the histograms for the $D^{n}_{L}$
and $S^{n}_{L}$ indicators are displayed with a vertical red line
corresponding to the position of the threshold. From this and following the
described process, we have obtained the initial conditions classified by both
indicators and we can compare how they have classified the trajectories with
the Poincaré section given in Fig. 3 in such a way that it is observed how the
chaos indicators have reconstructed the regular and chaotic regions of the
phase space as we expected, and very accurately. We want to highlight that, as
mentioned at the beginning of this section, the indicators based on Lagrangian
descriptors do not classify all the initial conditions correctly. To make this
clear, we have indicated in the corresponding figures with dots and red arrows
misclassified initial conditions which are regular initial conditions that had
been classified as chaotic. Additionally, we have observed that the
misclassification is more frequent in the regular region than in the chaotic
region of the phase space.

A) B)
C) D)

Figure 4: Classification of a random ensemble of $10^{4}$ initial conditions
carried out by means of the chaos indicators based on Lagrangian descriptors.
The orange dots correspond to chaotic motion, whereas blue dots indicate
regular behavior. For this simulation we have selected an energy of
$\mathcal{H}_{0}=20$, and the model parameters are set to $\alpha=1$ and
$\sigma=1$. A) Histogram showing the distribution of the $D^{n}_{L}$
indicator. We have marked with a vertical red line the threshold that
separates chaotic from regular behavior; B) Poincaré section overlaid with the
classified initial conditions using the $D^{n}_{L}$ indicator. The red dot
shows an initial condition that has not been correctly classified by the
method.; C) and D) Same analysis, but for the $S^{n}_{L}$ indicator.

##  4 Results

This section describes all the results we have obtained from the analysis of
the chaotic fraction of phase space trajectories in the double pendulum system
using the chaos indicators based on LD. In our numerical experiments we have
explored the parameter space $(\alpha,\sigma,\mathcal{H}_{0})$, where
$\alpha=l_{1}/l_{2}$ is the ratio of lengths, $\sigma=m_{1}/m_{2}$ the ratio
of masses and $\mathcal{H}_{0}$ represents the total energy of the system. For
the values of the parameters we have selected $\alpha_{i}=2^{i}$ and
$\sigma_{j}=2^{j}$, where $i,j\in\left\\{-4,-3,-2,-1,0,1,2,3,4\right\\}$, and
we have considered 170 different energy levels for each case
$(\alpha_{i},\sigma_{j})$ of the model parameters. The energies are chosen as
follows: 40 energy levels are distributed uniformly from the energy of the
local minimum of the potential energy surface ($\mathcal{H}_{1}$) up to the
energy of the local maximum of the potential ($\mathcal{H}_{4}$), and the
other 130 energy values are taken above $\mathcal{H}_{4}$ with a unit step.
For each simulation, we have generated an ensemble of $10^{4}$ initial
conditions, sampled randomly on the surface of section given by Eq. (47) using
a standard Monte Carlo approach.

One of the main results that can be obtained from our simulations is that, for
a given value of $\sigma$ (the mass fraction of $m_{1}$ with respect to
$m_{2}$), the maximum chaotic fraction of the phase space is attained when the
lengths of the pendulums are approximately equal, and decreases monotonically
as $\alpha$ (the length ratio of $l_{1}$ with respect to $l_{2}$) moves away
from unity. On the other hand, for a given $\alpha>1$, chaos increases in the
system as $\sigma$ gets smaller, whereas for a fixed $\alpha<1$, chaos
increases with $\sigma$. In particular, this means that if the first pendulum
is longer than the second, chaos is enhanced by making the second mass heavier
than the first one. On the contrary, in the case where the second pendulum is
longer, chaos grows as the mass of the first pendulum gets larger than that of
the second pendulum. All this information about the behavior of the system is
summarized in Fig. 5.

A)  B)

Figure 5: Maximum chaotic fraction of phase space as a function of the model
parameters $\alpha=l_{1}/l_{2}$ and $\sigma=m_{1}/m_{2}$ calculated with the
indicator $S^{n}_{L}$. The range of energy levels used for the simulations is
$[\mathcal{H}_{1},\mathcal{H}_{4}+130]$.

Another interesting feature of the simulations we have carried out is that
sudden peaks and valleys appear in the chaotic fraction of phase space
trajectories. For example, in the case where $\alpha=1$ and $\sigma=8$, this
behavior occurs in a short energy interval as displayed in Fig. 6. As the
system goes from an energy level of $\mathcal{H}_{0}=16.025$ to
$\mathcal{H}_{0}=19.025$, the chaotic fraction increases and decreases
abruptly by about $15\%$. This is a clear indication that significant changes
are taking place in the phase space structure of the system due to
bifurcations and the formation and destruction of KAM tori. These features are
nicely captured by the Poincaré sections displayed below the chaotic fraction
curve in Fig. 6. Notice also that this behavior manifests repeatedly along the
curve, but the peaks become less abrupt each time they occur.

Figure 6: Chaotic fraction ($\mathcal{C}$) versus energy ($\mathcal{H}$)
calculated with the chaos indicator $S^{n}_{L}$ for the model parameters
$\alpha=1$ and $\sigma=8$. The black dotted lines are the energies of the
local minimum ($\mathcal{H}_{1}$), of the index-1 saddles ($\mathcal{H}_{2}$
and $\mathcal{H}_{3}$) and that of the local maximum ($\mathcal{H}_{4}$) of
the potential energy surface of the double pendulum. The panels on the bottom
represent Poincaré sections for the energies marked on the blue curve above
with yellow dots. A sudden peak in the chaotic fraction indicates a relevant
change in the underlying phase space structure of the system.

In addition to the study of the peaks that appear in the chaotic fraction
curve, we have also performed different fits to study the trend of chaos
growth and decay in the system with energy. Our simulations indicate that for
values of the energy larger than that corresponding to the maximum of the
potential, the curve behaves like a decreasing exponential of the form
$\mathcal{A}e^{-\mathcal{B}x}$. For comparison, we have also considered a
linear fit of the form $Mx+N$, but the results obtained from this fit are much
worse than those provided by the exponential model. This seems to indicate
that chaos in the double pendulum system tends to follow an exponential
decrease from the energy of the maximum. For example, if we take the case
where $\alpha=1$ and $\sigma=8$, the $R^{2}$ coefficient for the exponential
fit in this region is approximately $0.9586$, while for the linear fit the
value obtained is $0.8048$, which shows that the linear fit is not suitable
for this section of the curve as it is clearly shown in Fig. 7 B). To support
our claim, we have carried out this exponential fit for a fixed length ratio
$\alpha=1$ and different values of the mass ratio $\sigma$. A summary of the
results obtained is included in the table displayed in Fig. 7 D). It is
important to note here that by using a fit with a decreasing exponential we
have assumed that chaos tends to 0 in the system when the energy becomes
arbitrarily large. However, we do not know if this assumption is true or not,
even though the curve that represents the chaotic fraction as a function of
the energy seems to point in that direction. Needless to say that one needs to
be careful with these claims about the asymptotic behavior of the system. It
is evident that one cannot make generalizations of this behavior to other
dynamical systems. Indeed, there are many Hamiltonian systems in the
literature where the chaotic fraction of the phase space does not decrease
with energy, see e.g. [3].

Another aspect that we have analyzed about the chaotic fraction as a function
of the model parameters is the growth of chaos in the energy regime that
covers the values from that of the minimum of the potential to the energy of
the first index-1 saddle. For this section of the curve, we have also studied
an exponential fit, since in this region the chaotic fraction grows very
rapidly with energy. We have observed that an exponential function is
perfectly suitable for this purpose, as can be seen in Fig. 7 where the fit
has the functional form $\mathcal{A}e^{\mathcal{B}x}$. We have performed a
detailed study for the exponential models in the case of $\alpha=1$ so that
the values of the coefficients and the $R^{2}$ values for the fits are
presented in the table included in Fig. 7 C). Moreover, Fig. 7 A) displays the
section of the chaotic fraction curve where this growth occurs, and we have
overlaid the exponential fit of the data in order to show how well the
exponential fit performs in this approximation. These results are in the line
of those discussed in [2], where they address this behavior but only for the
case where the pendulums have equal masses and lengths (that corresponds to
our case $\alpha=1$ and $\sigma=0$).

A) B)
C)  D)

Figure 7: Fits performed for the chaotic fraction as a function of energy
using exponential functions (red) and linear functions (black) for the model
parameters $\alpha=1$ and $\sigma=8$. A) Exponential fit for the section of
the curve in which the chaotic fraction grows; B) Exponential fit for the
section of the graph where the chaotic fraction decreases; C) Coefficients of
the exponential fit and $R^{2}$ indicator determined for different values of
the model parameters in the section of the curve where the chaotic fraction
increases; D) Coefficients of the exponential fit and $R^{2}$ indicator
determined for different values of the model parameters in the energy regime
where the chaotic fraction decreases.

##  5 Conclusions

In this paper we have quantified the extent of chaos and regularity in the
double pendulum by means of applying chaos indicators derived from LDs. Our
findings show how these diagnostic tools can be used to effectively measure
the degree of chaos in a computationally efficient manner. In order to study
how the chaotic fraction of trajectories varies with energy and the other
model parameters, we have shown that it is important to implement an adaptive
strategy to set the threshold value for the chaos indicator that is required
to accurately classify between chaotic and regular motion. The methodology is
based on determining the location of the minimum between the two peaks that
arise in the histogram of the chaos indicator obtained from the ensemble of
initial conditions used in the simulations. This process is relevant since
there is no universal threshold, as this critical value depends on many
factors such as the integration time, the distance between neighbors, the
model parameters and also the energy of the system.

Our parametric study of the degree of chaos in the double pendulum system has
allowed us to establish how the masses and lengths of the pendulums influence
the dynamics of this classical Hamiltonian system. The numerical analysis we
have carried out reveals that the systems attains the maximum fraction of
chaos when the pendulums have equal lengths. Another interesting feature of
the problem that our simulations highlight is that, if the second pendulum is
longer than the first one, chaos grows as the mass of the first pendulum is
increased with respect to that of the second one. Conversely, if the first
pendulum is longer, then chaos grows as the mass of the second pendulum is
larger than the first one.

Regarding the growth and decay of chaos in the system as a function of energy,
we have found that in the low-energy regime, that is, from the energy at the
minimum of the potential to that of the first index-1 saddle, chaos increases
following an exponential law. This result was also reported in [2], but only
for one of the cases we cover in this work. By performing a curve fitting
analysis, we have also observed from our numerical experiments that in the
energy regime starting at approximately the value of the local maximum of the
potential, the chaotic fraction shows an exponentially decreasing trend. This
behavior seems to indicate that the system will recover full regularity in the
high-energy regime. Interestingly, the reduction of chaos displayed by the
system is not monotone, as abrupt peaks and valleys appear in the chaotic
fraction curve. These sudden variations in the extent of chaos point in the
direction that bifurcations and other significant changes are taking place in
the underlying phase space of the system, as KAM islands are being generated
and destroyed.

The work we have carried out in this paper opens new avenues for future
research and raises many questions. For example, it would be interesting to
analyze how chaos grows and decays in different Hamiltonian systems and, in
particular, determine under what conditions full regularity is restored for
the high-energy regime. Regarding the appearance of abrupt variations in the
chaotic fraction, we believe that this phenomenon could be correlated with the
complexity of the transitions that are taking in the phase space structure. It
might be worth examining this relation by means of measuring the jaggedness of
the chaotic fraction curve using for instance the fractal dimension. Another
aspect of the problem that could be pursued is the development of control
strategies on the system so that chaos varies smoothly with energy, reducing
the sudden peaks and valleys that we have observed in the behavior of the
double pendulum. Our research efforts in the future will focus on addressing
all these questions, with the goal of gaining more insights and shedding some
light on the intertwined chaotic and regular nature of Hamiltonian systems. We
hope that this work could serve as a motivation for the further development of
these chaos diagnostic techniques, as their simple implementation and
computational efficiency make them ideal candidates to explore the chaotic
behavior of any dynamical system, especially those that are high-dimensional
and where the visualization of phase space structure is impractical.

## Acknowledgements

The authors would like to acknowledge the support received for this research
through the grant “Introducción a la Investigación” from the Universidad de
Alcalá.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] D. Assencio.  Double pendulum: Hamiltonian formulation.  https://dassencio.org/46, 2014.  Last accessed 9 February 2024.  URL: https://dassencio.org/46.
  * [2] S. Cabrera, E. D. Leonel, and A. C. Marti.  Regular and chaotic phase space fraction in the double pendulum.  arXiv preprint arXiv:2312.13436, 2023.
  * [3] P. Dahlqvist and G. Russberg.  Existence of stable orbits in the ${\mathit{x}}^{2}$${\mathit{y}}^{2}$ potential.  Phys. Rev. Lett., 65:2837–2838, Dec 1990.  doi:10.1103/PhysRevLett.65.2837.
  * [4] A. Darwish, S. Norouzi, G. Di Labbio, and L. Kadem.  Extracting Lagrangian coherent structures in cardiovascular flows using Lagrangian descriptors.  Physics of Fluids, 33(11):111707, 2021.  doi:10.1063/5.0064023.
  * [5] D. Deleanu.  The dynamics of a double pendulum: Classic and modern approach.  Annals of “Dunarea de Jos” University of Galati, Mathematics, Physics, Theoretical Mechanics, Fascicle II, Year, 3:203–212, 2011.
  * [6] V. J. García-Garrido, A. Ramos, A. M. Mancho, J. Coca, and S. Wiggins.  A dynamical systems perspective for a real-time response to a marine oil spill.  Marine Pollution Bulletin, 112(1):201–210, 2016.  doi:10.1016/j.marpolbul.2016.08.018.
  * [7] H. Goldstein.  Classical Mechanics.  Pearson, 2011.
  * [8] M. Hillebrand, S. Zimper, A. Ngapasare, M. Katsanikas, S. Wiggins, and C. Skokos.  Quantifying chaos using lagrangian descriptors.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(12):123122, 2022.  doi:10.1063/5.0120889.
  * [9] I. Intan Indiati, J. Saefan, and P. Marwoto.  Numerical approach of hamilton equations on double pendulum motion with axial forcing constraint.  Journal of Physics: Conference Series, 739(1):012066, 2016.  doi:10.1088/1742-6596/739/1/012066.
  * [10] K. Kaheman, J. J. Bramburger, N. Kutz, and S. L. Brunton.  Saddle transport and chaos in the double pendulum.  Nonlinear Dynamics, 111(8):7199–7233, 2023.  doi:10.1007/s11071-023-08232-0.
  * [11] H. Korsch, H. Jodl, and T. Hartmann.  Chaos: A Program Collection for the PC.  Springer Berlin Heidelberg, 2007.  doi:10.1007/978-3-540-74867-0.
  * [12] R. B. Levien and S. M. Tan.  Double pendulum: An experiment in chaos.  American Journal of Physics, 61(11):1038–1044, 1993.  doi:10.1119/1.17335.
  * [13] C. Lopesino, F. Balibrea-Iniesta, V. J. García-Garrido, S. Wiggins, and A. M. Mancho.  A theoretical framework for lagrangian descriptors.  International Journal of Bifurcation and Chaos, 27(01):1730001, 2017\.  doi:10.1142/S0218127417300014.
  * [14] J. A. J. Madrid and A. M. Mancho.  Distinguished trajectories in time dependent vector fields.  Chaos, 19:013111, 2009.  doi:10.1063/1.3056050.
  * [15] A. M. Mancho, S. Wiggins, J. Curbelo, and C. Mendoza.  Lagrangian descriptors: A method for revealing phase space structures of general time dependent dynamical systems.  Communications in Nonlinear Science and Numerical, 18:3530–3557, 2013.  doi:10.1016/j.cnsns.2013.05.002.
  * [16] C. Mendoza and A. M. Mancho.  Hidden geometry of ocean flows.  Phys. Rev. Lett., 105:038501, 2010.  doi:10.1103/PhysRevLett.105.038501.
  * [17] K. Miyamoto.  Long term simulations of the double pendulum by keeping the value of hamiltonian constant.  In 2013 8th EUROSIM Congress on Modelling and Simulation, pages 130–135. IEEE, 2013.  doi:10.1109/EUROSIM.2013.33.
  * [18] F. J. Muñoz Almaraz, E. Freire, and J. Galán-Vioque.  Bifurcation behavior of the furuta pendulum.  International Journal of Bifurcation and Chaos, 17(08):2571–2578, 2007.  doi:10.1142/S0218127407018634.
  * [19] S. Naik and S. Wiggins.  Detecting reactive islands in a system-bath model of isomerization.  Phys. Chem. Chem. Phys., 22:17890–17912, 2020.  doi:10.1039/D0CP01362E.
  * [20] S. Raffa, G. Merisio, and F. Topputo.  Finding regions of bounded motion in binary asteroid environment using lagrangian descriptors.  Communications in Nonlinear Science and Numerical Simulation, 121:107198, 2023.  doi:10.1016/j.cnsns.2023.107198.
  * [21] T. Shinbrot, C. Grebogi, J. Wisdom, and J. A. Yorke.  Chaos in a double pendulum.  American Journal of Physics, 60(6):491–499, 1992.  doi:10.1119/1.16860.
  * [22] C. Skokos, G. Gottwald, and J. Laskar.  Chaos Detection and Predictability.  Lecture Notes in Physics. Springer Berlin Heidelberg, 2016.  doi:10.1007/978-3-662-48410-4.
  * [23] T. Stachowiak and T. Okada.  A numerical analysis of chaos in the double pendulum.  Chaos, Solitons & Fractals, 29(2):417–422, 2006.  doi:10.1016/j.chaos.2005.08.032.
  * [24] W. Szumiński.  A new model of variable-length coupled pendulums: from hyperchaos to superintegrability.  Nonlinear Dynamics, pages 1–29, 2024.  doi:10.1007/s11071-023-09253-5.
  * [25] J. M. T. Thompson.  Chaos, fractals and their applications.  International Journal of Bifurcation and Chaos, 26(13):1630035, 2016\.  doi:10.1142/S0218127416300354.
  * [26] S. Zimper, A. Ngapasare, M. Hillebrand, M. Katsanikas, S. Wiggins, and S. Ch.  Performance of chaos diagnostics based on lagrangian descriptors. application to the 4d standard map.  Physica D: Nonlinear Phenomena, 453:133833, 2023.  doi:10.1016/j.physd.2023.133833.


# Why Existing Multimodal Crowd Counting Datasets Can Lead to Unfulfilled
Expectations in Real-World Applications

Martin Thißen Hochschule Darmstadt Prof. Dr. Elke Hergenröther Hochschule
Darmstadt

###### ABSTRACT

More information leads to better decisions and predictions, right? Confirming
this hypothesis, several studies concluded that the simultaneous use of
optical and thermal images leads to better predictions in crowd counting.
However, the way multimodal models extract enriched features from both
modalities is not yet fully understood. Since the use of multimodal data
usually increases the complexity, inference time, and memory requirements of
the models, it is relevant to examine the differences and advantages of
multimodal compared to monomodal models. In this work, all available
multimodal datasets for crowd counting are used to investigate the differences
between monomodal and multimodal models. To do so, we designed a monomodal
architecture that considers the current state of research on monomodal crowd
counting. In addition, several multimodal architectures have been developed
using different multimodal learning strategies. The key components of the
monomodal architecture are also used in the multimodal architectures to be
able to answer whether multimodal models perform better in crowd counting in
general. Surprisingly, no general answer to this question can be derived from
the existing datasets. We found that the existing datasets hold a bias toward
thermal images. This was determined by analyzing the relationship between the
brightness of optical images and crowd count as well as examining the
annotations made for each dataset. Since answering this question is important
for future real-world applications of crowd counting, this paper establishes
criteria for a potential dataset suitable for answering whether multimodal
models perform better in crowd counting in general.

### Keywords

Crowd Counting, Multimodal Learning, RGB-T, Transformer

##  1 Introduction

One of the biggest challenges of crowd counting in real-world applications is
dealing with varying lighting conditions. Since crowd counting can be very
important for event security and crowd monitoring, good performance
independent of lighting conditions is essential for real-world applications.
Especially at night, lighting is often poor, resulting in less contrast and
information in optical images and thus reducing the accuracy of prediction
models. In this case, thermal images are more suitable because they do not
rely on visible light. On the other hand, optical images can contain more
information during the daytime compared to monochrome thermal images due to
their color information. In addition, the environment may heat up during the
day, resulting in lower contrast in thermal images, as human body temperature
is almost constant. Overall, the use of both modalities seems to be symbiotic
and to lead to better results compared to the use of a single modality. Using
multiple modalities to train a model has led to state-of-the-art results in
many cases. In particular, with the rise of transformers [Vas17], where inputs
are transformed into homogeneous tokens, using multiple modalities such as
text or images in a model has become easier. In the area of monomodal crowd
counting, the use of transformers has not been fully explored. To our best
knowledge, with the exception of one work [Tia21], previous research has
focused only on convolutional networks. The use of transformers has tremendous
potential, as previous work [Zha16] [Li18] has often achieved better results
when improving the extraction of multi-scale features. Although existing work
[Liu21] [Pen20] concludes that the use of optical and thermal imagery leads to
better crowd counting predictions, it is not yet fully understood how such
models internally work and how they extract enriched features from both
modalities.

Apart from the lack of understanding of how multimodal models work internally,
it is not entirely understood whether the multimodal approach leads to better
crowd counting results in general or only under certain conditions. Further
research with potential influencing factors such as illumination, distance to
the crowd, or number of people per image is needed to gain more certainty
about whether multimodal crowd counting leads to better predictions in
general. For this reason, in this paper we investigate the impact of using
optical and thermal images simultaneously in crowd counting.

To investigate the impact of using optical and thermal images simultaneously
in crowd counting, we designed a monomodal and several multimodal
architectures consisting of the same key components. When we designed the
monomodal model, we took into account the latest developments in the field of
monomodal crowd counting. In addition, we have developed three multimodal
models that incorporate different strategies of multimodal learning. To allow
a comparison between the monomodal and the multimodal architectures, all key
components of the monomodal architecture are also part of the multimodal
architectures. The goal of this comparison is to find out whether multimodal
models lead to better crowd counting results in general or only under certain
conditions. Since this comparison led to interesting findings, we further
analyzed all the datasets used to compare the models. To this end, we examined
the relationship between the brightness of optical images and the number of
individuals in the image. We also randomly selected a subset of each dataset
and examined how individuals were labeled in the images from both modalities.

In examining the differences between the monomodal and the multimodal
architectures, we found that existing datasets have a bias toward thermal
images. This does not allow us to determine whether multimodal crowd counting
leads to better results in general or only under certain conditions. For this
reason, we have described criteria for a dataset suitable for investigating
the research question.

##  2 Related Work

Monomodal Crowd Counting: Crowd counting has been studied for decades. While a
few works have used thermal images for crowd counting, most works have used
optical images to examine crowd counting. As in other areas, the use of deep
learning models [Wan15] [Fu15] has led to more accurate predictions in crowd
counting. In recent years, the use of a density map-based approach for crowd
counting has become prevalent. Many recent works have addressed the question
of how to deal with scale variations in images. In particular, techniques such
as multi-column models [Zha16] or dilated convolutions [Li18] have been used
to extract multi-scale features from the image. Since such techniques aim to
increase the receptive field of a network, it was no surprise that state-of-
the-art results could be achieved by using a transformer encoder [Vas17] to
extract features [Tia21].

Multimodal Crowd Counting: Multimodal learning is becoming increasingly
relevant in the field of crowd counting. So far, the use of optical and
thermal images [Liu21] [Pen20] [Tan22] [Gu22] as well as the use of optical
and depth images [Lia22] [Lia19] has been investigated. However, depth images
provide only a limited depth range (0 $\sim$ 20 meters), making them
unsuitable for many real-world crowd counting applications [Liu21]. Also, when
using depth images, there is still the problem that less information is
available in poorly illuminated scenes. For this reason, we will focus on the
use of optical and thermal images in this paper. While all work concludes that
the additional use of thermal images leads to better predictions in crowd
counting, it is not fully understood under what circumstances it is beneficial
to complement optical images with thermal images to obtain better predictions.
Previous work has focused primarily on constructing a novel model architecture
that outperforms the state-of-the-art in multimodal crowd counting. While this
approach proves the effectiveness of the models created, it does not allow us
to fully understand how complementary information is extracted from both
modalities.

Multimodal Crowd Counting Datasets: Similar to different multimodal models,
two different datasets [Liu21] [Pen20] consisting of optical and thermal image
pairs have been published in recent years. The dataset published by Peng et
al. [Pen20] was acquired with a drone and contains 3,600 image pairs.
Furthermore, this dataset contains information about distance (scale of
individuals), illumination and crowd count per image pair. The other dataset,
which was published by Liu et al. [Liu21], contains 2,030 image pairs. The
image pairs of this dataset were taken from a normal perspective. Information
on the number of individuals and lighting is available for each image pair.

##  3 Effectiveness of Multimodal Crowd Counting

To allow a comparison between monomodal and multimodal architectures, we first
developed a monomodal architecture. This monomodal model takes into account
recent advances in the field of monomodal crowd counting and its main
components are reused in subsequent multimodal architectures to allow a fair
comparison. Since the constructed monomodal architecture is heavily inspired
by recent advances in monomodal crowd counting and does not incorporate any
new strategies, we only used one monomodal model for comparison.

###  3.1 Monomodal Architecure

The monomodal architecture designed in this work is inspired by the work of
Tian et al. [Tia21] as well as the implementation of the work realized in
[Wan21a]. The CCTrans model designed by Tian et al. [Tia21] achieves state-of-
the-art results on multiple monomodal crowd counting benchmarks [Zha16]
[Wan21b] [Idr18]. Our monomodal architecture is shown in Fig. 1. Instead of
Twins [Chu21], which was used by Tian et al. [Tia21], we used PVTv2 [Wan21c]
as the transformer-based backbone in our architecture. By empirical analysis,
we found that the PVTv2 architecture leads to better results for us. More
specifically, for our monomodal architecture, we used the PVTv2 B0 variant,
which allows shorter training time and requires less computational resources.
However, this leads to slightly worse results compared to other PVTv2 variants
with more parameters. This was acceptable to us, as our primary goal was not
to construct a novel architecture with state-of-the-art results. Furthermore,
we adopted the pyramid feature aggregation and regression head of Tian et al.
[Tia21], but used the convolution kernel sizes used in [Wan21a]. Again,
through empirical analysis, we found that these kernel sizes led to slightly
better results for us.

Figure 1: The architecture of our monomodal model, which is inspired by the
work of Tian et al. [Tia21] as well as the implementation of the work realized
in [Wan21a]. The input image is first transformed into tokens. From these
tokens, features are extracted by a hierarchical transformer-based backbone.
The hierarchical feature maps are then aggregated and finally used by the
regression head to predict the crowd count. The parameter $d$ indicates the
dilatation rate used.

###  3.2 Multimodal Architectures

After the monomodal architecture was designed, three different multimodal
architectures were designed that incorporate different strategies of
multimodal learning. As mentioned before, the key characteristics of the
monomodal architecture are also incorporated in the three different multimodal
architectures. The idea behind this is that when using the same weight
initialization (prior) and the same model properties (which constrain the
hypothesis space), better results can only be explained by more information
provided by the additional modality (data). In this work, we chose to use
early and late fusion as two simple multimodal strategies. These have also
been used in previous work on multimodal crowd counting [Liu21] [Pen20]. In
addition, we apply a more advanced deep fusion strategy using the Information
Aggregation and Distribution Module (IADM) of Liu et al. [Liu21] which has
been shown to be effective for multimodal crowd counting.

Early Fusion Model: With the early fusion strategy, modalities are fused at
the beginning of the model. For this purpose, the constructed monomodal model
was adapted to support 6-channel inputs by changing the amount of filters in
the first layer. Thus, the multimodal early fusion model has the same number
of parameters as the monomodal model.

Late Fusion Model: In contrast to the early fusion strategy, the fusion of
modalities takes place at the end of the model with the late fusion strategy.
The idea here is that features of both modalities are first extracted
individually. Thus, except for the final layer ($1\times 1$ convolution), both
modalities are investigated with the constructed monomodal model individually.
Then, the extracted feature maps from both individual columns are
concatenated. Based on the concatenated feature maps, a density map is then
finally predicted by a $1\times 1$ convolution. Hereby, the late fusion model
requires around twice as many parameters as the monomodal model and the early
fusion model.

Deep Fusion Model: In contrast to the early fusion and late fusion
architectures, the multimodal information exchange in the deep fusion
architecture takes place during feature extraction. For this purpose, a third
column is added to the architecture, which extracts the complementary
information of both modalities. In particular, this is done by using the IADM
of Liu et al. [Liu21]. Through the IADM, information is exchanged between the
modality-specific columns and the cross-modality column. However, this only
takes place during feature extraction in the backbone, as shown in Fig. 2. Of
all the models used in this work, this architecture requires the most
parameters.

Figure 2: The architecture of our deep fusion model. To extract complementary
information and enable exchange between modality-specific and modality-shared
columns, we use the IADM of Liu et al. [Liu21]. In their work, it was shown
that the use of the IADM is effective for multimodal data.

###  3.3 Evaluation

To evaluate the performance of the monomodal model and the three multimodal
models, we used the mean absolute error (MAE) and the root mean squared error
(RMSE). Both of these measures are widely used in crowd counting. The use of
these measures allows comparison of our results with the results of other
work. The mean absolute error and root mean squared error are defined as
follows:

$MAE=\frac{1}{N}\sum_{i=1}^{N}|y_{i}-\hat{y}_{i}|\;,$ (1)

$RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}}\;,$ (2)

where $N$ is the number of image pairs, $y_{i}$ is the ground-truth number of
individuals in image pair $i$, and $\hat{y}_{i}$ is the predicted number of
individuals for image pair $i$.

###  3.4 Training

Overall, our training approach is heavily inspired by the training approach
chosen by Tian et al. [Tia21]. The B0 variant of the PVTv2 architecture was
initialized with pre-trained weights in all experiments. We used random
cropping with a cropping size of 256 for both dimensions and horizontal
flipping with a probability of 50% as augmentation strategies. In addition,
AdamW [Los19] was used as the optimizer and a batch size of 8 was chosen for
training. The learning rate was $1\mathrm{e}{-5}$ in all experiments, but was
increasingly regulated by a weight decay of $1\mathrm{e}{-4}$. Bayesian loss
[Ma19] with a sigma value of 8 was used as a loss function. The models were
trained for 60 epochs in all experiments.

Modality | Architecture | MAE | RMSE
---|---|---|---
RGB | Monomodal | 26.48 | 55.28
T | Monomodal | 15.19 | 28.27
RGB-T | Early Fusion | 14.92 | 25.86
RGB-T | Late Fusion | 13.83 | 25.16
RGB-T | Deep Fusion | 14.32 | 24.64
RGB-T | BL + IADM [Liu21] | 15.61 | 28.18
RGB-T | TAFNet [Tan22] | 12.38 | 22.45
Table 1: Performance of the different architectures on the RGBT-CC [Liu21] dataset. The use of thermal images leads to dramatically better results compared to optical images. Moreover, the multimodal approach leads to better results than the monomodal approach. Modality | Architecture | MAE | RMSE
---|---|---|---
RGB | Monomodal | 10.40 | 16.44
T | Monomodal | 6.70 | 10.20
RGB-T | Early Fusion | 7.41 | 11.43
RGB-T | Late Fusion | 7.01 | 11.18
RGB-T | Deep Fusion | 7.20 | 11.45
RGB-T | MMCCN [Pen20] | 7.27 | 11.45
RGB-T | MFCC [Gu22] | 7.96 | 12.50
Table 2: Performance of the different architectures on the Drone-RGBT [Pen20]
dataset. Surprisingly, using thermal images solely with the monomodal
architecture led to the best result for the Drone-RGBT dataset. In contrast,
using optical images solely with the monomodal architecture leads to
considerably worse results.

###  3.5 Results

The results for all constructed models on both datasets are shown in Tab. 1
and Tab. 2. Three aspects in particular caught our attention, which we
describe in more detail below.

Discrepancy between optical and thermal images in both datasets. One of the
first things we noticed is that the monomodal model performs much better on
thermal images than on optical images, as can be seen in Tab. 1 and Tab. 2.
This holds true for both datasets. Nevertheless, the discrepancy is larger for
the RGBT-CC dataset than for the Drone-RGBT dataset. Since we used the exact
same model and training approach, these results raise the question of whether
thermal images are more suitable for crowd counting in general. Before
investigating this question, we first wanted to gain a better understanding of
both data sets. The investigation is described in more detail in Section 4.

The monomodal model performs better than the multimodal models for the Drone-
RGBT [Pen20] dataset. Contrary to our assumption that the multimodal approach
of using optical and thermal images would lead to better crowd counting
predictions, using thermal images solely led to the best result for the Drone-
RGBT dataset. This result further affirmed our motivation to gain a better
understanding of both datasets. To our best knowledge, we have achieved state-
of-the-art results for the Drone-RGBT dataset using the monomodal
architecture.

IADM [Liu21] seems to be less effective with transformer encoders. Comparing
the three multimodal models, the late fusion model achieves the best results
on both datasets. The deep fusion model, although more complex and shown to be
effective by Liu et al. [Liu21], performs worse in our study than the late
fusion model. Since Liu et al. also compared the IADM to a late fusion model,
the most obvious explanation for this is the use of a transformer encoder in
our work. Liu et al. did not use a transformer encoder in their work.
Nevertheless, a more detailed investigation beyond this work is needed to
better understand why the IADM is less effective when used with transformer
encoders.

##  4 Analysis of Existing Multimodal Crowd Counting Datasets

To understand more profoundly whether thermal images are better for crowd
counting in general, or whether the characteristics of the datasets used lead
to better results on thermal images, we used two different approaches.

###  4.1 Relationship Between Brightness and Crowd Count

First, we investigated the relationship between the brightness of optical
images and the number of individuals. We suspected that many optical images in
both datasets were taken in poorly illuminated environments, which could be
the reason for the discrepancy between thermal and optical images. This would
also be in line with our main motivation to use multimodal data. Since the two
metrics we used consider the counting error and are sensitive to outliers, we
thought it is relevant to investigate the relationship between brightness and
crowd count. To measure the brightness of an optical image, we used the
following equation:

$Brightness=\frac{\sum_{i=1}^{W*H}R_{i}+G_{i}+B_{i}}{3*W*H}\;,$ (3)

where $W$ is the width and $H$ is the height of the optical image. $R_{i}$,
$G_{i}$ and $B_{i}$ represent the three color values of pixel $i$. The
relationship between brightness and crowd count for both datasets are shown in
Fig. 3 and Fig. 4.

Figure 3: Scatter plot showing the relationship between the brightness of
optical images and crowd count in the RGBT-CC dataset. It can be seen that the
RGBT-CC dataset is unbalanced in terms of brightness and crowd count. Many
images with very low brightness have a high crowd count.

Figure 4: Scatter plot showing the relationship between the brightness of
optical images and crowd count in the Drone-RGBT dataset. The relationship
between brightness and crowd count appears very uniform compared to the
distribution of the RGBT-CC dataset.

The RGBT-CC [Liu21] dataset is unbalanced regarding brightness and crowd
count. The RGBT-CC dataset contains many images with very low brightness and
high crowd count, as can be seen in Fig. 3. In comparison, the images in the
Drone-RGBT dataset are much brighter on average and the overall distribution
between brightness and number of individuals is much more balanced, as shown
in Fig. 4. Since both metrics are sensitive to outliers and many optical
images with very low brightness (low information) have a high crowd count in
the RGBT-CC dataset, we assume that this explains the bigger discrepancy for
the RGBT-CC dataset between optical and thermal images.

This finding has serious implications on our research question. Since this
imbalance of the dataset likely affects all trained models and results in
higher activations for thermal input, we believe that the research question
cannot be thoroughly investigated with the RGBT-CC dataset. In particular, we
assume that many optical images with low brightness (low optical information)
and a high crowd count (high error) will cause the model to pay more attention
to thermal images as the counting error is propagated back into the network
during training. In this way, it is difficult to verify whether multimodal
crowd counting leads to better results in general, since a certain condition
(low brightness, high crowd count) has a great impact on the training of the
model as well as the metrics. Nevertheless it is important to say that images
with low brightness and high crowd count are not a problem, but are important
and desirable for training a robust crowd counting model. However, we are
concerned about whether our research question can be fairly investigated due
to an inherent correlation between the number of people and brightness in the
RGBT-CC dataset.

###  4.2 Annotation Sample Analysis

Our second approach to better understand both datasets was to perform a sample
analysis of how the annotations were made. Since both datasets contain two
different modalities recorded with two different cameras, we wanted to
understand if images of both modalities were synchronized and how perspective
changes were handled (because the cameras were probably next to each other
during the recording). We decided to perform the sample analysis of how the
annotations were made since both datasets provide shared annotations for both
modalities. To this purpose, we randomly selected 10% of the image pairs per
dataset and visualized the annotations in the images of both modalities to
verify how the individuals were labeled in each image.

Only thermal images were used to label individuals in both datasets. By
randomly selecting 10% of all image pairs per dataset and visualizing the
annotations for both modalities, we found that both datasets used only the
thermal image to label individuals. Examples of both datasets showing that
only thermal images were used to label individuals are provided in the
Appendix in Fig. 6 and Fig. 6.

All image pairs of the Drone-RGBT dataset were taken at night. We have seen
that the optical images in the Drone-RGBT dataset are on average brighter than
the optical images in the RGBT-CC dataset. However, by looking at the
annotations for each image pair in the Drone-RGBT dataset, we perceived that
all images were taken at night. To gain more confidence in this perception, we
looked at all the optical images in the Drone-RGBT dataset. In this way, we
found that all image pairs in the Drone-RGBT dataset were taken at night.
Nevertheless, many images were taken in environments with much artificial
light, therefore the optical images are on average brighter than those of the
RGBT-CC dataset. The fact that all images were taken at night adds a new
perspective to the results obtained with the Drone-RGBT dataset. This leads to
the assumption that optical images do not provide additional information at
night and that a monomodal approach with thermal images leads to better
results. Further research beyond this paper is needed to validate this
assumption.

For image pairs in the RGBT-CC dataset, individuals were sometimes visible in
one modality but not the other. Liu et al. [Liu21] have already stated in
their work that optical and thermal images in the RGBT-CC dataset are not
strictly aligned because they were captured with different sensors. However,
when examining the annotations, we found that not only were the image pairs
not strictly aligned, but sometimes individuals were visible in one modality
but not the other. Examples for this are provided in the Appendix in Fig. 7.

##  5 Criteria for a Multimodal Crowd Counting Dataset

Because both datasets have some weaknesses that make it difficult to draw
general conclusions about the effectiveness of multimodal crowd counting, we
decided to set criteria for a suitable dataset. Overall, the image pairs
should be taken evenly throughout the day. In this way, the variability of the
two modalities gets extensively covered. Ideally, this would even take into
account different seasons and climate zones. Also, the crowd count per image
pair should be independent of when the image was taken. This allows for an
equal influence of both modalities on the multimodal model during training, as
no modality receives more attention due to a higher counting error. When
labeling individuals, both modalities should be considered so that later
models can learn to extract the information from both modalities and
incorporate it into the prediction (even when one modality contains little
information and the other contains much). Furthermore, the images for both
modalities should be taken simultaneously. In this way, the images of both
modalities are aligned as precisely as possible, which allows the use of the
same annotations for both modalities.

##  6 Is Multimodal Crowd Counting Better in General?

The goal of this work was to find out if the simultaneous use of optical and
thermal images leads to better predictions in crowd counting in general. We
found that existing datasets have a bias toward thermal images, making it
difficult to draw general conclusions about the effectiveness of multimodal
crowd counting. The results on the Drone-RGBT dataset indicate that solely
using thermal images at night results in better predictions than a multimodal
approach. Since the RGBT-CC dataset contains both daytime and nighttime
images, the better predictions with multimodal data seem to indicate that the
multimodal approach leads to better results during the daytime. However, these
assumptions are by no means proven, but could serve as hypotheses for future
research. Furthermore, we encourage the creation of a multimodal dataset in
order to be able to investigate such hypotheses. We have provided criteria for
the creation of such a dataset in the previous Section 5. However, it remains
an open question whether multimodal crowd counting (including technical
challenges like perspective distortion and synchronization between modalities)
is the perfect approach. It could also be the case that two monomodal models
produce better results than one multimodal model. For example, one monomodal
model could be used with optical images during the day and another with
thermal images at night.

##  7 Conclusion

In this work, we found that existing multimodal crowd counting datasets have a
bias toward thermal images. For this reason, we outlined criteria for a
balanced dataset. To our best knowledge, we also obtained state-of-the-art
results on the multimodal Drone-RGBT dataset. Interestingly, for this we used
solely thermal images and the monomodal model constructed in this work.
Considering the results of this work, we encourage the creation of a
multimodal dataset that meets the criteria outlined in this paper. In this
way, we can understand more profoundly whether the simultaneous use of optical
images and thermal images leads to better predictions in crowd counting in
general.

## REFERENCES

  * [Vas17] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. Gomez, L. Kaiser and I. Polosukhin, “Attention is All you Need” in Advances in Neural Information Processing Systems, Curran Associates, Inc., vol. 30, pp. 6000–6010, 2017.
  * [Tia21] Y. Tian, X. Chu and H. Wang (2021), “CCTrans: Simplifying and Improving Crowd Counting with Transformer” unpublished.
  * [Zha16] Y. Zhang, D. Zhou, S. Chen, S. Gao and Y. Ma, “Single-Image Crowd Counting via Multi-Column Convolutional Neural Network” in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016, pp. 589-597.
  * [Li18] Y. Li, X. Zhang and D. Chen, “CSRNet: Dilated Convolutional Neural Networks for Understanding the Highly Congested Scenes” in 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2018, pp. 1091-1100.
  * [Liu21] L. Liu, J. Chen, H. Wu, G. Li, C. Li and L. Lin, “Cross-Modal Collaborative Representation Learning and a Large-Scale RGBT Benchmark for Crowd Counting” in 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2021, pp. 4821-4831.
  * [Pen20] T. Peng, Q. Li and P. Zhu, “RGB-T Crowd Counting from Drone: A Benchmark and MMCCN Network” in Computer Vision - ACCV 2020: 15th Asian Conference on Computer Vision, 2020, pp. 497–513.
  * [Tan22] H. Tang, Y. Wang and L. Chau (2022), “TAFNet: A Three-Stream Adaptive Fusion Network for RGB-T Crowd Counting” unpublished.
  * [Lia22] D. Lian, X. Chen, J. Li, W. Luo and S. Gao, “Locating and Counting Heads in Crowds With a Depth Prior” in in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 44, 2022, pp. 9056–9072.
  * [Lia19] D. Lian, J. Li, J. Zheng, W. Luo and S. Gao, “Density Map Regression Guided Detection Network for RGB-D Crowd Counting and Localization” in 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1821–1830.
  * [Wan21a] F. Wang and F. Taioli (2021) CCTrans: Simplifying and Improving Crowd Counting with Transformer (Code reproduction) [Source code]. https://github.com/wfs123456/CCTrans
  * [Wan21b] Q. Wang, J. Gao, W. Lin and X. Li, “NWPU-Crowd: A Large-Scale Benchmark for Crowd Counting and Localization” in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 43, no. 6, 2021, pp. 2141–2149.
  * [Idr18] H. Idrees et al., “Composition Loss for Counting, Density Map Estimation and Localization in Dense Crowds” in European Conference on Computer Vision, 2018, pp. 544–559.
  * [Chu21] X. Chu, Z. Tian, Y. Wang, B. Zhang, H. Ren, X. Wei, H. Xia and C. Shen H., “Twins: Revisiting the Design of Spatial Attention in Vision Transformers” in Advances in Neural Information Processing Systems, vol. 34, 2021, pp. 9355–9366.
  * [Wan21c] W. Wang, E. Xie, X. Li, D. Fan, K. Song, D. Liang, T. Lu, P. Luo, L. Shao (2021), “PVT v2: Improved Baselines with Pyramid Vision Transformer” unpublished.
  * [Los19] I. Loshchilov and F. Hutter, “Decoupled Weight Decay Regularization” in 7th International Conference on Learning Representations, 2019.
  * [Ma19] Z. Ma, X. Wei, X. Hong and Y. Gong, “Bayesian Loss for Crowd Count Estimation With Point Supervision” in IEEE/CVF International Conference on Computer Vision (ICCV), 2019, pp. 6141–6150.
  * [Gu22] S. Gu and Z. Lian (2022), “A Unified Multi-Task Learning Framework of Real-Time Drone Supervision for Crowd Counting” unpublished.
  * [Wan15] C. Wang, H. Zhang, L. Yang, S. Liu, and X. Cao, “Deep People Counting in Extremely Dense Crowds” in Proceedings of the 23rd ACM international conference on Multimedia, 2015, pp. 1299–1302.
  * [Fu15] M. Fu, P. Xu, X. Li, Q. Liu, M. Ye, and C. Zhu, “Fast crowd density estimation with convolutional neural networks” in Engineering Applications of Artificial Intelligence, vol. 43, 2015, pp. 81–88.

## Appendix

Figure 5: Illustration of the annotations (yellow squares) of the RGBT-CC
[Liu21] dataset for both modalities. It can be seen that the annotations were
created based on the thermal images. For better understanding, examples of
inaccuracies in the annotation of the optical images have been encircled. The
thermal images are also encircled in the corresponding places, but the
annotations are correct there.

Figure 6: The annotations of the Drone-RGBT [Pen20] dataset are shown in
yellow squares. It can be seen that the annotations were made on the basis of
the thermal images. For better understanding, individuals have been encircled
who are easily recognizable in the optical image, but have not been annotated.
In the thermal image, on the other hand, it can be seen that precisely these
individuals are difficult to recognize, which is probably why they were not
annotated.

Figure 7: Two examples from the RGBT-CC dataset where a person is visible in
one modality but not the other. It can be seen that this is due to the time-
delayed capture of the images from both modalities.