# Perfect synchronization in complex networks with higher order interactions

Sangita Dutta1 Prosenjit Kundu2<EMAIL_ADDRESS>Pitambar Khanra3
Chittaranjan Hens4 Pinaki Pal1<EMAIL_ADDRESS>1Department of
Mathematics, National Institute of Technology, Durgapur 713209, India
2Dhirubhai Ambani Institute of Information and Communication Technology,
Gandhinagar, Gujarat, 382007, India  3Department of Mathematics, University at
Buffalo, State University of New York, Buffalo, USA  4Center for Computational
Natural Science and Bioinformatics, International Institute of Informational
Technology, Gachibowli, Hyderabad 500032, India

###### Abstract

We propose a framework for achieving perfect synchronization in complex
networks of Sakaguchi-Kuramoto oscillators in presence of higher order
interactions (simplicial complexes) at a targeted point in the parameter
space. It is achieved by using an analytically derived frequency set from the
governing equations. The frequency set not only provides stable perfect
synchronization in the network at a desired point, but also proves to be very
effective in achieving high level of synchronization around it compared to the
choice of any other frequency sets (Uniform, Normal etc.). The proposed
framework has been verified using scale-free, random and small world networks.
In all the cases, stable perfect synchronization is achieved at a targeted
point for wide ranges of the coupling parameters and phase-frustration. Both
first and second order transitions to synchronizations are observed in the
system depending on the type of the network and phase frustration. The
stability of perfect synchronization state is checked using the low
dimensional reduction approach. The robustness of the perfect synchronization
state obtained in the system using the derived frequency set is checked by
introducing a Gaussian noise around it.

##  I Introduction

Higher order structures, such as three- and four-way interactions in addition
to pairwise interactions, are widespread in neurological, biological,
ecological, and sociological systemsMajhi  _et al._ (2022); Battiston  _et
al._ (2020); Petri  _et al._ (2021); Alvarez-Rodriguez  _et al._ (2021);
Iacopini  _et al._ (2019); Martignon  _et al._ (1995); Yu  _et al._ (2011);
Giusti  _et al._ (2016); Reimann  _et al._ (2017); Sizemore  _et al._ (2018).
In ecological systems such higher order interactions (HOI) where three or more
species interact with each other can stabilize large ecological communities
Chatterjee  _et al._ (2022); Grilli  _et al._ (2017); Bairey  _et al._ (2016).
On the other hand, recent studies suggest the mesoscopic organization of brain
through higher order interactions allows efficient information processing and
offers useful guidelines for performing complex tasks Petri  _et al._ (2014);
Giusti  _et al._ (2015); Sizemore  _et al._ (2018); Ganmor  _et al._ (2011).
Motivated by these practical implications, the network science community is
also concentrating on understanding the various types of collective behavior
ranging from synchronization to epidemic spreading that can go beyond
conventional paired interactions Battiston  _et al._ (2020); Petri  _et al._
(2021); Gambuzza  _et al._ (2021); Skardal and Arenas (2019).

Synchronization observed in the flocking pattern of birds, or rhythmic
flashing of fireflies, can be modelled with the interacting nonlinear
dynamical units Strogatz (2004); Pikovsky  _et al._ (2003); Kuramoto (2003);
Newman (2010); Cohen and Havlin (2010); Barrat  _et al._ (2008). The classic
Kuramoto dynamics, one of the celebrated dynamical model used for studying the
synchronization Kuramoto (1984); Acebrón  _et al._ (2005), encodes the phase
evolution of each node. In a complex network, nodes with non-identical natural
frequencies are connected to one another via paired linkages. Depending on the
coupling configuration or particular frequency design, such complicated
pairwise connections could cause the entire system to reveal “continuous”,
“discontinuous”, or “optimal” synchronization Kuramoto (1984); Gómez-Gardenes
_et al._ (2011); Leyva  _et al._ (2012); Gómez-Gardenes  _et al._ (2007);
Zhang  _et al._ (2015); Coutinho  _et al._ (2013); Ichinomiya (2004); Skardal
and Arenas (2019). However, a slight phase-lag between the oscillators may
destabilize or erode the stable synchronization states, and destroy the
switch-like discontinuous synchronization transition for certain cases
Chamlagai  _et al._ (2022); Kundu  _et al._ (2018, 2017, 2020); Kundu and Pal
(2019); Sakaguchi and Kuramoto (1986); Skardal  _et al._ (2015). Against this
backdrop, we explore the impact of HOI in a network of coupled phase-lag
oscillators. Particularly, in presence of HOI, and fixed coupling strength, we
seek a suitable frequency set which can lead the system to a perfect Kundu
_et al._ (2018) global synchronization state where all phases are in unison.

The HOI can induce abrupt synchronization transitions between incoherent state
and coherent state without any correlations between network structure and
dynamical functions. These synchronized states are stable even in the presence
of repulsive pairwise coupling due to such higher order interaction Skardal
and Arenas (2020). In a system, three-way interactions in addition to pairwise
interactions can cause abrupt desynchronization transitions without any abrupt
synchronisation transitions and extensive multistable partially synchronized
statesSkardal and Arenas (2019) may appear.

Recent advances suggest that, higher order interaction may be vital and play
an important role in general oscillator systems like brain dynamics Petri  _et
al._ (2014); Giusti  _et al._ (2016); Sizemore  _et al._ (2018), collaboration
networks Patania  _et al._ (2017), social contagion network Iacopini  _et al._
(2019) etc. In particular, interaction in 2- simplexes are important to
describe correlation in neuronal activity in brain Cui (2017) providing a
missing link between the structure and dynamics. In fact, considering higher
order interaction is useful where different types of correlation exists
between the nodes in a coupled oscillator systems. Despite these findings,
enhancement of synchronization in phase frustrated dynamics with higher order
interaction has not been explored so far.

Recently, it has been demonstrated that a trade-off between pairwise and
higher order interactions can result in a higher level of synchronization
(optimal synchronization) in relatively weaker coupling strength Skardal  _et
al._ (2021) by using an optimal frequency set obtained from the dominant
eigenvector of the composite Laplacian, which encodes both pairwise and
triangular interactions in the absence of phase frustration. However, in spite
of the reported nontrivial effects of phase frustration on the dynamics of
pairwise networked system Sakaguchi and Kuramoto (1986); Omel’Chenko and
Wolfrum (2012); Kundu  _et al._ (2017); Kundu and Pal (2019); Kundu  _et al._
(2018); Brede and Kalloniatis (2016), the role of it is least explored in
networks in presence of higher order interactions.

Thus, in this paper, we focus on the topic of targeted global synchronization
in phase frustrated complex networks in presence of pairwise as well as higher
order interactions and ask the question ‘Can we determine a set of frequencies
to achieve perfect synchronization at a targeted point of the parameter space
in such networks?’ In the process, we employ the analytical approaches
reported in Skardal  _et al._ (2014); Kundu  _et al._ (2018) and develop a
framework based on linear theory to achieve the goal of perfect
synchronization at a targeted point. Analytical approach eventually leads to
the determination of a set of frequencies involving the structural and
dynamical properties of the network for the achievement of perfect
synchronization. Detailed numerical demonstration with different networks
following the analytic finally confirms the achievement of targeted perfect
synchronization using the derived frequency.

The organization of the paper is as follows. Sections II and III describe the
model and the analytical frame work for deriving a frequency set for
achievement of perfect synchronization at a targeted point. Subsequent section
IV presents the detailed numerical results on perfect synchronization in
different networks. Then the stability and robustness of the synchronization
states are presented in the sections V and VI. A general discussion and
conclusions are made in the section VII. We proceed with the model description
in the next section.

##  II Model Description

We consider different networks of coupled Sakaguchi-Kuramoto Sakaguchi and
Kuramoto (1986); Kuramoto (2003) phase oscillators of size $N$ with higher
order interactions along with pairwise interactions Skardal and Arenas (2020,
2019); Adhikari  _et al._ (2022). Dynamics of each oscillators in the network
is governed by the equation

$\displaystyle\dot{\theta_{i}}$ $\displaystyle=$
$\displaystyle\omega_{i}+K_{1}\sum_{j=1}^{N}A_{ij}\sin(\theta_{j}-\theta_{i}-\alpha)$
$\displaystyle+\frac{K_{2}}{2}\sum_{j=1}^{N}\sum_{l=1}^{N}B_{ijl}\sin(2\theta_{j}-\theta_{l}-\theta_{i}-\beta),$
$\displaystyle\hskip 113.81102pti=1,2,\dots,N,$

where $\theta_{i}$ is the phase and $\omega_{i}$ is the intrinsic natural
frequency of the $i$-th oscillator. $\alpha$ and $\beta$ act as the
frustration terms in the system corresponding to the pairwise and the triadic
interactions respectively. $K_{1}$ and $K_{2}$ are the coupling strengths
associated with $1$-simplex (pairwise interaction) and $2$-simplex (triadic
interaction) respectively. $A_{ij}$ is the $ij$-th element of the adjacency
matrix $A=(A_{ij})_{N\times N}$ associated to the 1-simplex, where $A_{ij}=1$,
if $i$-th and $j$-th node are connected and $A_{ij}=0$, otherwise. Similarly,
$B_{ijl}$ is the adjacency tensor associated to $2$-simplex where $B_{ijl}=1$
if there is a triadic connection between the $i$th, $j$th and $l$th nodes and
$B_{ijl}=0$, otherwise. Note that the networks under consideration is
undirected. As a result, we have $A_{ij}=A_{ji}$ and
$B_{ijl}=B_{ilj}=B_{jil}=B_{lji}=B_{lij}=B_{jli}$ for all admissible $i$ and
$j$. Here we choose $2$ simplicial complex, which leads to the consideration
of the triadic and pairwise interactions only. By the definition of simplicial
complex, every link in that triangle will be in the simplicial complex. So,
the adjacency tensor can be written in the form $B_{ijl}=A_{ij}A_{jl}A_{li}$.

To quantify the level of synchronization in the network we use the Kuramoto
order parameter $r$ given by

$re^{i\psi(t)}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}},$ (2)

where $\psi(t)$ is the average phase of the oscillators at time $t$. The value
of the order parameter lies between $0$ and $1$. $r=0$ indicates the system is
in the incoherent state, while $r=1$ indicates the fully synchronized state of
the system. The coupled equations (II) is used to develop an analytical
framework for the achievement of perfect synchronization in the network. The
details are described in the next section.

##  III Analytical Framework

Here we describe the analytical framework used to derive a frequency set for
achieving perfect synchronization at a targeted point in the parameter space
following the approach reported in Skardal  _et al._ (2014); Kundu  _et al._
(2018). Linearization of the system (II) about the synchronized state
($|\theta_{j}-\theta_{i}|\rightarrow 0$) leads to the equations

$\displaystyle\dot{\theta_{i}}$ $\displaystyle=$
$\displaystyle\omega_{i}+[-K_{1}k_{i}^{(1)}\sin\alpha-
K_{2}k_{i}^{(2)}\sin\beta]-K_{1}\cos\alpha[k_{i}^{(1)}\theta_{i}-\sum_{j=1}^{N}A_{ij}\theta_{j}]-K_{2}\cos\beta[k_{i}^{(2)}\theta_{i}-\sum_{j=1}^{N}A_{ij}(\sum_{l=1}^{N}A_{jl}A_{li})\theta_{j}$
(3) $\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{j=1}^{N}A_{ji}(\sum_{l=1}^{N}A_{il}A_{lj})\theta_{j}]~{}(i=1,2,\dots,N).$

The above equations can be written in vector form as,

$\dot{{\bf{\theta}}}=\omega+d-L{\bf{\theta}},$ (4)

where $d=-K_{1}k^{(1)}\sin\alpha-K_{2}k^{(2)}\sin\beta$ and $L=K_{1}\cos\alpha
L^{(1)}+K_{2}\cos\beta L^{(2)}$ is the composite Laplacian consisting of
Laplacian of pairwise interaction and Laplacian of triadic interaction
$L^{(1)}$ and $L^{(2)}$ respectively. Here, $L^{(1)}$ and $L^{(2)}$ are
defined as

$\displaystyle L^{(1)}=D^{(1)}-A^{(1)}~{}\mathrm{and}~{}\hskip
11.38092ptL^{(2)}=D^{(2)}-(A^{(2)}-\frac{A^{(2)T}}{2}),$

where $A^{(1)}$ denotes the adjacency matrix $A$,
$D^{(1)}=$diag$(k_{1}^{(1)},k_{2}^{(1)},\dots,k_{N}^{(1)})$,
$A^{(2)}=A*(A^{2})^{T}$,
$D^{(2)}=$diag$(k_{1}^{(2)},k_{2}^{(2)},\dots,k_{N}^{(2)})$ with $*$ denoting
the Hadamard product. Also, $k_{i}^{(1)}=\sum_{j=1}^{N}A_{ij}$ is the degree
of $i$-th node considering the pairwise interaction only and
$k_{i}^{(2)}=\frac{1}{2}\sum_{j=1}^{N}\sum_{l=1}^{N}B_{ijl}$ is the degree of
$i$-th node while considering the triadic interactions.

Now, we derive the steady state solution by suitable varying the reference
frame as $\omega\rightarrow\omega-\Omega$, where $\Omega$ is the group angular
velocity and putting $\psi=0$. Therefore $\dot{\theta}=0$ in Eq. (4) gives the
stationary points as,

$\theta^{*}=L^{\dagger}(\omega+d),$ (5)

where $L^{\dagger}$ is the pseudo inverse of the Laplacian $L$, which is of
the form

$L^{\dagger}=\sum_{j=2}^{N}\lambda_{j}^{-1}v_{j}v_{j}^{T},$ (6)

where $0=\lambda_{1}<\lambda_{2}\leq\dots\leq\lambda_{N}$ are the eigenvalues
of $L$ and $v_{j}$ $(j=1,2,\dots,N)$ are the eigenvectors of $L$ with respect
to the eigenvalues $\lambda_{j}$.

The order parameter in Eq. (2) then approximated as

$r\approx 1-\frac{\parallel\theta^{*}\parallel^{2}}{2N},$ (7)

where

$\displaystyle\parallel\theta^{*}\parallel^{2}$ $\displaystyle=$
$\displaystyle\langle\theta^{*},\theta^{*}\rangle$ $\displaystyle=$
$\displaystyle\langle L^{\dagger}(\omega+d),L^{\dagger}(\omega+d)\rangle$
$\displaystyle=$ $\displaystyle\sum_{j=2}^{N}\lambda_{j}^{-2}\langle
v_{j},\omega+d\rangle^{2}.$

Now, by substituting the above expression in the equation (7) we get the order
parameter as

$\displaystyle r$ $\displaystyle=$ $\displaystyle
1-\frac{1}{2N}\sum_{j=2}^{N}\lambda_{j}^{-2}\langle v_{j},\omega+d\rangle^{2}$
(8) $\displaystyle=$ $\displaystyle 1-\frac{1}{2}J(\omega,L),$

where $J(\omega,L)=\frac{1}{N}\sum_{j=2}^{N}\lambda_{j}^{-2}\langle
v_{j},\omega+d\rangle^{2}$ is known as synchrony alignment function (SAF).
Here we note that in Eq. (8), $J(\omega,L)\rightarrow 0$ limit leads to $r=1$,
which denotes the perfect synchronization state. Therefore, for the
achievement of perfect synchronization at a given point, we choose

$\displaystyle\omega+d=0.$ (9)

Thus the frequency set for the achievement of perfect synchronization at a
targeted point in the parameter space is given by

$\omega=K_{1}^{(p)}k^{(1)}\sin\alpha+K_{2}^{(p)}k^{(2)}\sin\beta,$ (10)

where $K_{1}^{(p)}$ and $K_{2}^{(p)}$ are the targeted coupling strength for
pairwise and triadic interaction respectively. We now proceed for the
numerical verification of the proposed framework for achieving perfect
synchronization at targeted coupling strengths for fixed phase-lags $\alpha$
and $\beta$, considering different networks.

##  IV Numerical Verification

The achievement of perfect synchronization in complex networks at a targeted
point in the parameter space using the frequency set (analytically derived by
the method described in the previous section) is numerically verified in this
section. For that purpose, we numerically simulate the Sakaguchi-Kuramoto
model with higher order interactions on top of four different complex networks
of same size ($N=10^{3}$). The considered networks are one Erdős-Rényi (ER)
Erdős and Rényi (1959, 1960, 1961) network; two scale-free(SF) Barabási and
Bonabeau (2003) networks with exponents $\gamma=2.5$ and $3.2$; and one small
world (SW) Watts and Strogatz (1998) network. The SF network with exponent
$2.5$ is of mean degree $\langle k\rangle=6$ and each of the other three
networks have mean degree $\langle k\rangle=8$. Numerical simulations are
performed using $4$-th order Runge-Kutta (RK4) method with step size $\delta
t=0.01$ for sufficiently long time after removing the transients. For next few
numerical illustrations we take $\alpha=\beta$ and set the target point at
$K_{1}^{(p)}=0.1$ and $K_{2}^{(p)}=0.5$. In each case, we numerically continue
the solutions both in forward and backward directions.

Figure 1: Perfect synchronization for different networks. Numerically computed
order parameter $r$ as a function of pairwise coupling strength $K_{1}$ shown
for phase lags $\alpha$=0.1(blue), 0.3(magenta), 0.5(cyan), 0.7(red). Triadic
coupling $K_{2}$ is set to 0.5. Green dot indicates the targeted perfect
synchronization at $K_{1}=0.1$ with $r=1$. (a) ER network for different values
of $\alpha$. $r=1$ is achieved at $K_{1}=0.1$ both in forward and backward
simulations. (b) SF network with $\gamma=2.5$ attain $r=1$ except for forward
simulations at $\alpha=0.7$. (c) SF with $\gamma=3.2$ shows perfect
synchronization for all four lags in forward and backward simulations. (d) SW
network achieves perfect synchronization in forward and backward simulations
only for $\alpha=0.1$. For the other lags backward paths only attain $r=1$.
Figure 2: Synchronization diagram for different frequency set. Order parameter
$r$ as a function of pairwise coupling strength $K_{1}$. $K_{2}$ is fixed to
0.5 and phase lag is taken 0.5. Green, red and blue curve is associated with
the natural frequency chosen from normal distribution, uniform distribution
and derived for perfect synchronization respectively.

First we fix $K_{2}=0.5$ and vary $K_{1}$ around the target point
$K_{1}^{(p)}=0.1$ for the achievement of perfect synchronization in each of
the considered networks for four different values of
$\alpha=\beta=0.1,0.3,0.5,\rm{and}~{}0.7$. For each $\alpha$, we use the
frequency set derived by the equation (10), and obtain the statistically
steady solutions after integrating the system for long time. Although here we
have chosen $\alpha=\beta$, using the derived frequency set, perfect
synchronization can be achieved at a targeted point for $\alpha\neq\beta.$
This fact will be illustrated towards the end of this section.

Starting from $K_{1}=0$, the solution is then continued in the forward
direction by changing the value of $K_{1}$ in small steps and using the last
values of the phase variables of the previous simulation as the initial
condition for the present simulation. In the process, we continue the solution
till $K_{1}=0.3$ and then continue the solution similarly in the backward
direction by gradually reducing the value of $K_{1}$ in steps of $\delta
K_{1}=0.001$ write up to the starting value of $K_{1}$. The variation of the
synchronization order parameter ($r$) with $K_{1}$ for fixed $K_{2}=0.5$
obtained from the simulation data are shown in the FIG. 1. From FIG. 1(a) and
(c) it is observed that in case of ER and SF network with $\gamma=3.2$,
perfect synchronization is achieved at the targeted point (marked by green
dot) during forward as well backward continuation for all values of $\alpha$.
However, for other three cases (FIG. 1(b) and (d)) perfect synchronization is
achieved at the targeted point during backward continuation for all considered
$\alpha$, while, for forward continuation, it is achieved at the targeted
point for lower values of $\alpha$. Interesting to note here that, for the
choice of analytically derived frequency set hysteresis is observed both for
SF and SW networks and large hysteresis loop is formed for higher $\alpha$.
This is counter intuitive, and special result in presence of higher order
interactions, because previous studies with only pairwise interactions suggest
the negative impact of frustration parameter on the hysteresis loop Kundu  _et
al._ (2017). Therefore, from the above discussion it is apparent that using
the analytically derived frequency set, perfect synchronization can be
achieved in a wide variety of networks for large values of the frustration
parameter at a targeted point in the parameter space.

Figure 3: Synchronization of SK model with HOI for different networks.
Synchronization diagrams describe the order parameter $r$ as a function of
triadic coupling $K_{2}$ for fixed $K_{1}=0.1$. The blue, magenta, cyan and
red curves presents the transitions for $\alpha=$0.1, 0.3, 0.5, 0.7
respectively. The green dot denote the perfect synchronization. (a) ER network
achieve $r=1$ for all four lags. (b) SF network with $\gamma=2.5$ attain $r=1$
except for $\alpha=0.7$(forward path). (c) SF network with $\gamma=3.2$ attain
perfect synchronization for all $\alpha$ in both paths. (d) SW network attain
$r=1$ only during forward transition(except $\alpha=0.1$).

Further, the analytically derived frequency which has functional dependence
both on the pairwise as well as triadic interactions not only provide perfect
synchronization at a targeted point in the parameter space but also helps in
achieving very high level of synchronization around the targeted point
compared to other standard frequencies. FIG. 2 shows the comparison of the
level of synchronization in the considered networks for different frequency
sets. Here, simulations are performed with uniform and normal frequencies
along with the derived frequency set. In each case, the derived frequency set
clearly provide substantially higher level of synchronization around the
targeted point.

Next we fix the pairwise coupling strength $K_{1}=0.1$ and vary the triadic
coupling strength $K_{2}$ around the targeted point ($K_{2}^{(p)}=0.5$) and
perform numerical simulations using the SK model on the same complex networks
considered above. As before, in this case also we numerically continue the
solutions both in forward and backward directions by changing $K_{2}$ in the
steps of $\delta{K_{2}}=0.005$ for different values of $\alpha.$ The variation
of $r$ with $K_{2}$ as obtained from the simulation data in different cases
are shown in the FIG. 3 which clearly shows the impact of triadic interaction
in the system. It is observed that for all four networks perfect
synchronization is achieved at the targeted coupling with different transition
paths. ER and SF networks show continuous transition to synchronization for
small $\alpha$. Whereas, at $\alpha=0.7$ both of the SF networks exhibit small
hysteric window around the perfect synchronization point. For higher value of
$\alpha$, the system (phase-lag oscillators in SW network) gradually reaches
to the highest point ($r=1$), and then it sharply drops to zero. We have also
observed that during the backward transition, the perfect as well as higher
degree of synchronization is achieved only for small range of coupling. The
reason is as follows: during the backward transition, the system is already in
fully incoherent states, thus it tries to remain in the incoherent state if we
decrease the coupling strength. However, when it comes closer to
$K_{2}^{(p)}=0.5$, the coupling strength dependent frequency drives the entire
system to reach the perfect (or near perfect) synchronization regime. Further
decrease of coupling strength reduces the level or degree of synchronization
and drives the entire system to be fully desycnhronized. These results confirm
that using the analytically derived frequency set one can achieve perfectly
synchronized state either in the forward transition path or in the backward
transition path. Subsequently we analyze the analytically derived frequency
set.

At this point let us recall that the derived frequency set promotes very high
level of synchronization around the targeted point as the pairwise coupling
strength $K_{1}$ varies for fixed targeted triadic coupling strength (see FIG.
2). Next to check for the similar behavior of the level of synchronization we
perform simulation by vary $K_{2}$ around the targeted point for fixed
$K_{1}=0.1$. FIG. 4 shows that in this case also the derived frequency
provides very high level of synchronization compared to other frequency sets.

Figure 4: Synchronization diagram for different frequency set for parameter
$K_{2}$. Order parameter $r$ is plotted varying the higher order coupling
strength $K_{2}$ for fixed $K_{1}=0.1$ and phase lag 0.5. The system
synchronize for the derived perfect frequency set (blue) in all taken networks
and avoids synchronization in normal (green) and uniform (red) frequency
distributions. Figure 5: Distribution of the frequency for perfect
synchronization using the four different networks namely ER(blue),
SF($\gamma=2.5$)(red), SF($\gamma=3.2$)(purple), SW(brown). In the first row
the frequency is calculated for phase lag $\alpha=0.1$ and the second row the
phase lag $\alpha=0.7$. We take the pairwise coupling strength
$K_{1}^{(p)}=0.1$ and the higher order coupling strength $K_{2}^{(p)}=0.5$.

Now to understand the contributions of pairwise and triadic interactions in
the derived frequency for achieving perfect synchronization at a targeted
point, we analyze the frequency distributions corresponding to the pairwise
($\omega_{1}=K_{1}^{(p)}k^{(1)}\sin\alpha$) and triadic interactions
($\omega_{2}=K_{2}^{(p)}k^{(2)}\sin\alpha$) in the expression (10) along with
the total derived frequency distribution. FIG. 5 shows the distribution of the
derived frequency for perfect synchronization from the equation (10) for four
different networks and two phase lag values ($0.1$ and $0.7$). The figure also
displays the distributions of $\omega_{1}$ and $\omega_{2}$ separately in the
insets of each panel. It is observed that the range of the frequencies
increases substantially with lag. As the sets $\omega_{1}$ and $\omega_{2}$
are linearly correlated with the pairwise degree $k^{(1)}$ and triadic degree
$k^{(2)}$, the distributions also looks like the distribution of the pairwise
and triadic degrees respectively. Since ER network have small number of
triangles, from the distribution of $\omega_{2}$ it is clear that the effect
of triangular interactions is less than the effect of pairwise interactions
(FIG. 5(a)-(b)). For the SF networks also there are more pairwise interactions
than triadic (FIG. 5(c)-(f)). For these networks perfect synchronization is
attained at the targeted coupling up to lag $0.7$. Only in case of SW network
there are higher number of triadic interactions (FIG. 5(g)-(h)), as a reason
the effect of triadic interaction in the constructed frequency is much higher
compared to ER or SF networks. This effect of triadic interaction increases
the heterogeneity in the system and prevent the system to synchronize Chutani
_et al._ (2021), although the predicted frequency is able to overcome this
constraint and reach the global synchrony for a small vicinity of the
parameter space.

As mentioned earlier, although we have presented all the results by setting
$\alpha=\beta$, the proposed scheme works in the general case where
$\alpha\neq\beta.$ To check that we take $\alpha=0.2$ and $\beta=0.1$ and
perform the numerical simulation to achieve perfect synchronization for
$K_{1}^{(p)}=0.1$ and $K_{2}^{(p)}=0.5$ using the derived frequency set given
by the Eq. (10). FIG. 6 shows the variation of the order parameter computed
from the simulation data around the targeted point $K_{1}=0.1$ for fixed
$K_{2}=0.5$ for four considered networks. Also the variation of $r$ with
$K_{2}$ around the targeted point for fixed $K_{1}=0.1$ for these networks are
shown in FIG. 7. From these two figures also the achievement of perfect
synchronization at the targeted point for different phase lags associated with
the pairwise and triadic interactions is quite apparent. Next we perform low
dimensional reduction of the system to understand the stability of the
achieved synchronization state using the derived frequency set.

Figure 6: Synchronization profile for the derived frequency set at
$\alpha=0.2$ and $\beta=0.1$. Order parameter $r$ as a function of pairwise
coupling strength $K_{1}$. $K_{2}$ is fixed to 0.5. Perfect synchronization
(green dot) is achieved at the targeted coupling $K_{1}=0.1$ for all
considered networks. Figure 7: Synchronization profile for the derived
frequency set at $\alpha=0.2$ and $\beta=0.1$. Order parameter $r$ vs $K_{2}$
for the fixed pairwise coupling $K_{1}=0.1$. All considered networks attain
$r=1$ (green dot) at targeted $K_{2}=0.5$ both in forward(red) and
backward(blue) path.

##  V Low dimensional reduction

To analyze the stability of the synchronized state around the targeted
parameter values we reduce the networked system in Eq. (II) to it’s low
dimensional form using the collective coordinate approach Gottwald (2015);
Pinto and Saa (2015); Brede and Kalloniatis (2016). In this approach the phase
of each oscillators is approximated by their intrinsic frequency using

$\theta_{i}(t)=\chi(t)\omega_{i},$ (11)

where $\chi(t)$ is the collective coordinate, which is time dependent. To
ensure the validity of this approach we aim to minimize the error given by

$\displaystyle\epsilon_{i}(\chi)$ $\displaystyle=$
$\displaystyle\dot{\chi}\omega_{i}-\omega_{i}-K_{1}\sum_{j=1}^{N}A_{ij}\sin(\chi(\omega_{j}-\omega_{i})-\alpha)-\frac{K_{2}}{2}$
(12)
$\displaystyle\sum_{j=1}^{N}\sum_{l=1}^{N}B_{ijl}\sin(\chi(2\omega_{j}-\omega_{l}-\omega_{i})-\beta).$

This error will be minimum if it is orthogonal to the tangent space of the
solution space given by Eq. (11), which is spanned by
$\frac{\partial\theta_{i}}{\partial\chi}=\omega_{i}$Gottwald (2015); Pinto and
Saa (2015). Now, projecting this error to the specified subspace and using the
orthogonality property we obtain one-dimensional evolution equation for
$\chi(t)$

$\frac{d\chi}{dt}=g(\chi),$ (13)

where $g(\chi)$ is given by

$\displaystyle g(\chi)$ $\displaystyle=$ $\displaystyle
1+\frac{K_{1}}{\sigma^{2}}\sum_{i=1}^{N}\omega_{i}\sum_{j=1}^{N}A_{ij}\sin(\chi(\omega_{j}-\omega_{i})-\alpha)+\frac{K_{2}}{2\sigma^{2}}$
(14)
$\displaystyle\sum_{i=1}^{N}\omega_{i}\sum_{j=1}^{N}\sum_{l=1}^{N}B_{ijl}\sin(\chi(2\omega_{j}-\omega_{l}-\omega_{i})-\beta)$

and $\sigma^{2}=\sum_{i=1}^{N}\omega_{i}^{2}$. The one-dimensional
differential Eq. (13) will have a stable equilibrium if $g(\chi)=0$ and
$g^{\prime}(\chi)<0$. At the equilibrium $\chi$ is independent of time and
hence Eq.(11) indicates that all $\theta_{i}$’s will be time independent,
which implies the system will exhibit a phase-locked solution. If such fixed
point occurs at $\chi=0$, all phases will be zero, resulting perfect
synchronization.

Figure 8: Variation of $g$ obtained in equation(14) as a function of $\chi$
under perfect frequency for ER network. (a) Coupling $K_{2}$ is set to 0.5 and
the value of $K_{1}$ are taken as 0.06(cyan), 0.08(green), 0.092(magenta),
0.1(blue), 0.12(red). (b) Coupling $K_{1}$ is fixed to 0.1 and the value of
$K_{2}$ varies as 0.01(cyan), 0.15(green), 0.28(magenta), 0.5(blue),
0.8(red).(c) Derived frequency (blue), normally (red) and uniformly (green)
distributed frequency is plotted, showing that only derived frequency reaches
the stable synchronization state.

In FIG. 8(a), and 8(b) we show the behaviour of $g(\chi)$ under derived
frequency for perfect synchronization $\omega$ (targeted at
$K_{1}=0.1,K_{2}=0.5$ respectively), for a set of values of $K_{1}$ and
$K_{2}$ in ER network. We now fix $K_{1}$ to $0.1$ and vary $K_{2}$ in
FIG.8(b). For $K_{2}=0.01$ and $0.15$, $g(\chi)$ does not intersect
$g(\chi)=0$. Here, $g(\chi)$ touches the line $g(\chi)=0$ for the first time
at $K_{2}=0.28$ which provides a stable phase locked solution. From this value
of $K_{2}$ onward $g(\chi)$ crossed the line $g(\chi)=0$ with a negative
slope, which indicates that system has a stable phase-locked solution (as seen
in FIG.3(a)). As per expectation the system exhibits perfect synchronization
at $K_{2}=0.5$ ($g(\chi)=0$ and $g^{\prime}(\chi)<0$ at $\chi=0$). For larger
values of $K_{2}$ the synchronization remains stable. We do the analysis for
other two networks (SF and SW) which we found very similar.
In FIG. 8(a) we fixed $K_{2}$ to $0.5$ and choose $K_{1}$ as
$0.06,~{}0.08,~{}0.092,~{}0.1,~{}0.12$. For $K_{1}=0.06$ and $0.08$ $g(\chi)$
does not intersect $g(\chi)=0$, i.e stable synchronization is not attended at
this point. At $0.092$, $g(\chi)$ touches $g(\chi)=0$, that means stable
synchronization occurs for the first time. As the $K_{1}$ value increases from
$0.092$ the crossing point tends to zero and finally at $K_{1}=0.1$(targeted
value), $g(\chi)$ is zero and $g^{\prime}(\chi)<0$ at $\chi=0$, i.e the
perfect synchronization is obtained. For the values of $K_{1}>0.1~{}(0.12)$
the solution remains stable, as we see it in FIG.1(a).

FIG.8(c) shows the relation between $\chi$ and $g(\chi)$ for three types of
natural frequency $\omega$, namely, frequency for perfect synchronization,
normally distributed frequency and uniformly distributed frequency. Here the
underline network is ER network. We find that the choice of $\omega$ which
induce perfect synchronization, derived in Eq.(10) gives stable equilibrium at
$\chi=0$. For the other two frequency set (normal and uniform distributed)
$g(\chi)$ even never intersect $g(\chi)=0$. This shows that the
synchronization is stable only for the choice of frequency set derived in
Eq.(10) among the considered set of frequencies.

##  VI Robustness of the frequency set

The robustness of the derived frequency set ($\omega$) for perfect
synchronization is now checked by adding a perturbation to $\omega$ in the
form of Gaussian noise such that $\omega_{i}=\omega_{i}+\delta\omega_{i}$,
where $\delta\omega_{i}$ is drawn from a normal distribution
$N(0,\sigma\omega_{i})$, whose mean is $0$ and standard deviation is
$\sigma\omega_{i}$, which represents the multiplicative noise proportional to
$\omega_{i}$. Due to the effect of the perturbation in the optimal frequency,
the system will deviate from the perfect synchronized state. We define this
deviation from synchronization state as synchronization loss and define as
$\rho=1-r$, where $r$ is the synchronization order parameter defined in Eq. 2.
Now, from Eq. (7) we get

$\rho=\frac{1}{2N}\parallel\theta\parallel^{2}\sim\frac{1}{2}{\mathrm{Var}}(\theta),$
(15)

where ${\mathrm{Var}}(\theta)$ denotes the variance of $\theta$. Substituting
the perturbed $\omega$ in Eq. (5), results

$\theta=L^{\dagger}\delta\omega.$

Since $L^{\dagger}$ and $\delta\omega$ are approximately independent to each
other, we can write

$\mathrm{Var}(\theta)=(L^{\dagger})^{2}\mathrm{Var}(\delta\omega)=(L^{\dagger})^{2}\omega^{2}\sigma^{2}.$
(16)

Above equation confirms that $\mathrm{Var}(\theta)\propto\sigma^{2}$, where
the proportion constant is $(L^{\dagger})^{2}\omega^{2}$, which is independent
of $\sigma$ . Using this fact, we reach to a relation between
$\mathrm{Var}(\theta)$ and $\sigma$ as $\mathrm{Var}(\theta)\sim\sigma^{2}$.
Therefore from (15) we get

$\rho\sim\sigma^{2}.$ (17)

Figure 9: Synchronization error$(\rho)$ under deviated frequency. We compute
synchronization loss(blue) after adding noise $N(0,\sigma\omega_{i})$ to the
perfect frequency set. For $\sigma<1$, this error $\rho$ follows
$\rho\sim\sigma^{2}$(blue). For $\sigma\geq 1$, synchronization fully lost,
i.e $\rho\rightarrow 1$.

A numerical verification is done in FIG.9 using 4 networks, such as ER, SF
with two exponents($\gamma=2.5$ and $3.2$) and SW networks. It presents $\rho$
as a function of $\sigma$. In all four cases, this figures show that for small
$\sigma$, $\rho$ is small, means the deviation from synchronization is small
when the deviation of the frequency from the designed frequency set is small.
As the value of $\sigma$ increases (the deviation from the perfect frequency
set increases), $\rho$ also increases following the rule derived in Eq. 17.
After a certain value of $\sigma$ the value of $\rho$ saturates at $1$ because
at this stage the high deviation in the frequency set lead to a complete
incoherent (desynchronized state) with $r\rightarrow 0$ which lead to
$\rho\rightarrow 1$.

##  VII Discussion

In this article, we developed a general mathematical framework to derive the
natural frequencies of the nodes of a complex network which can ensure perfect
synchronization at considerably lower coupling strength. We presented a
synchrony alignment function that measures the interplay between network
structure and oscillator heterogeneity and allows to get a set of frequency
for perfect synchronization. Focusing on the Sakaguchi-Kuramoto model with
higher order interaction, we described the effect of pairwise coupling as well
as higher order coupling to reach the perfectly synchronized state at a
targeted coupling strength for four different types of networks. We found that
the analytically derived natural frequencies involve both structural and
dynamical information of the phase frustrated Kuramoto model with higher order
interaction. The synchronization is promoted by a strong alignment of the
frequency vector with the most dominant Laplacian eigenvectors and the pseudo
inverse operator. In all the cases, we found that, our derived frequency can
help the oscillators to reach to the perfectly synchronized state where as any
other frequency set does not reach perfect synchronization state even with
much higher coupling strength. Later, we derive a low dimensional model to
analyze the stability of the perfectly synchronized state using the collective
coordinate approach. This provides a clear understanding of the stability of
the synchronization state in different coupling strength and different
frequency sets. We also checked the robustness of the derived frequency
$\omega$ by adding a small perturbation in the form of Gaussian noise. This
shows that even we add some noise in the derived frequency set, the
synchronization state is not lost unless the frequency deviation is too high.

ACKNOWLEDGMENTS

S.D. is supported by the INSPIRE program of the DST, India (Code No.
IF190605).

## References

  * Majhi  _et al._ (2022) S. Majhi, M. Perc, and D. Ghosh, Journal of the Royal Society Interface 19, 20220043 (2022).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Physics Reports 874, 1 (2020).
  * Petri  _et al._ (2021) G. Petri, S. Musslick, B. Dey, K. Özcimder, D. Turner, N. K. Ahmed, T. L. Willke, and J. D. Cohen, Nature Physics 17, 646 (2021).
  * Alvarez-Rodriguez  _et al._ (2021) U. Alvarez-Rodriguez, F. Battiston, G. F. de Arruda, Y. Moreno, M. Perc, and V. Latora, Nature Human Behaviour 5, 586 (2021).
  * Iacopini  _et al._ (2019) I. Iacopini, G. Petri, A. Barrat, and V. Latora, Nature Communications 10, 1 (2019).
  * Martignon  _et al._ (1995) L. Martignon, H. Von Hassein, S. Grün, A. Aertsen, and G. Palm, Biological Cybernetics 73, 69 (1995).
  * Yu  _et al._ (2011) S. Yu, H. Yang, H. Nakahara, G. S. Santos, D. Nikolić, and D. Plenz, Journal of Neuroscience 31, 17514 (2011).
  * Giusti  _et al._ (2016) C. Giusti, R. Ghrist, and D. S. Bassett, Journal of Computational Neuroscience 41, 1 (2016).
  * Reimann  _et al._ (2017) M. W. Reimann, M. Nolte, M. Scolamiero, K. Turner, R. Perin, G. Chindemi, P. Dłotko, R. Levi, K. Hess, and H. Markram, Frontiers in Computational Neuroscience , 48 (2017).
  * Sizemore  _et al._ (2018) A. E. Sizemore, C. Giusti, A. Kahn, J. M. Vettel, R. F. Betzel, and D. S. Bassett, Journal of Computational Neuroscience 44, 115 (2018).
  * Chatterjee  _et al._ (2022) S. Chatterjee, S. Nag Chowdhury, D. Ghosh, and C. Hens, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 103122 (2022).
  * Grilli  _et al._ (2017) J. Grilli, G. Barabás, M. J. Michalska-Smith, and S. Allesina, Nature 548, 210 (2017).
  * Bairey  _et al._ (2016) E. Bairey, E. D. Kelsic, and R. Kishony, Nature Communications 7, 1 (2016).
  * Petri  _et al._ (2014) G. Petri, P. Expert, F. Turkheimer, R. Carhart-Harris, D. Nutt, P. J. Hellyer, and F. Vaccarino, Journal of The Royal Society Interface 11, 20140873 (2014).
  * Giusti  _et al._ (2015) C. Giusti, E. Pastalkova, C. Curto, and V. Itskov, Proceedings of the National Academy of Sciences 112, 13455 (2015).
  * Ganmor  _et al._ (2011) E. Ganmor, R. Segev, and E. Schneidman, Proceedings of the National Academy of sciences 108, 9679 (2011).
  * Gambuzza  _et al._ (2021) L. V. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti, Nature Communications 12, 1 (2021).
  * Skardal and Arenas (2019) P. S. Skardal and A. Arenas, Physical Review Letters 122, 248301 (2019).
  * Strogatz (2004) S. Strogatz,  _Sync: The emerging science of spontaneous order_ (Penguin UK, 2004).
  * Pikovsky  _et al._ (2003) A. Pikovsky, M. Rosenblum, J. Kurths, and J. Kurths,  _Synchronization: a universal concept in nonlinear sciences_ , Vol. 12 (Cambridge University Press, 2003).
  * Kuramoto (2003) Y. Kuramoto,  _Chemical oscillations, waves, and turbulence_ (Courier Corporation, 2003).
  * Newman (2010) M. Newman,  _Networks: an introduction_ (Oxford University Press, 2010).
  * Cohen and Havlin (2010) R. Cohen and S. Havlin,  _Complex networks: structure, robustness and function_ (Cambridge University Press, 2010).
  * Barrat  _et al._ (2008) A. Barrat, M. Barthelemy, and A. Vespignani,  _Dynamical processes on complex networks_ (Cambridge University Press, 2008).
  * Kuramoto (1984) Y. Kuramoto, in  _Chemical oscillations, waves, and turbulence_ (Springer, 1984) pp. 111–140.
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, Reviews of Modern Physics 77, 137 (2005).
  * Gómez-Gardenes  _et al._ (2011) J. Gómez-Gardenes, S. Gómez, A. Arenas, and Y. Moreno, Physical Review Letters 106, 128701 (2011).
  * Leyva  _et al._ (2012) I. Leyva, R. Sevilla-Escoboza, J. Buldú, I. Sendiña-Nadal, J. Gómez-Gardenes, A. Arenas, Y. Moreno, S. Gómez, R. Jaimes-Reategui, and S. Boccaletti, Physical Review Letters 108, 168702 (2012).
  * Gómez-Gardenes  _et al._ (2007) J. Gómez-Gardenes, Y. Moreno, and A. Arenas, Physical Review Letters 98, 034101 (2007).
  * Zhang  _et al._ (2015) X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Physical Review Letters 114, 038701 (2015).
  * Coutinho  _et al._ (2013) B. Coutinho, A. Goltsev, S. Dorogovtsev, and J. Mendes, Physical Review E 87, 032106 (2013).
  * Ichinomiya (2004) T. Ichinomiya, Physical Review E 70, 026116 (2004).
  * Chamlagai  _et al._ (2022) P. R. Chamlagai, D. Taylor, and P. S. Skardal, Physical Review E 106, 034202 (2022).
  * Kundu  _et al._ (2018) P. Kundu, C. Hens, B. Barzel, and P. Pal, EPL (Europhysics Letters) 120, 40002 (2018).
  * Kundu  _et al._ (2017) P. Kundu, P. Khanra, C. Hens, and P. Pal, Physical Review E 96, 052216 (2017).
  * Kundu  _et al._ (2020) P. Kundu, P. Khanra, C. Hens, and P. Pal, EPL (Europhysics Letters) 129, 30004 (2020).
  * Kundu and Pal (2019) P. Kundu and P. Pal, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 013123 (2019).
  * Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, Progress of Theoretical Physics 76, 576 (1986).
  * Skardal  _et al._ (2015) P. S. Skardal, D. Taylor, J. Sun, and A. Arenas, Physical Review E 91, 010802 (2015).
  * Skardal and Arenas (2020) P. S. Skardal and A. Arenas, Communications Physics 3, 1 (2020).
  * Patania  _et al._ (2017) A. Patania, G. Petri, and F. Vaccarino, EPJ Data Science 6, 1 (2017).
  * Cui (2017) Y. Cui, Inc., Redwood City, CA, United States (2017).
  * Skardal  _et al._ (2021) P. S. Skardal, L. Arola-Fernández, D. Taylor, and A. Arenas, Physical Review Research 3, 043193 (2021).
  * Omel’Chenko and Wolfrum (2012) E. Omel’Chenko and M. Wolfrum, Physical Review Letters 109, 164101 (2012).
  * Brede and Kalloniatis (2016) M. Brede and A. C. Kalloniatis, Physical Review E 93, 062315 (2016).
  * Skardal  _et al._ (2014) P. S. Skardal, D. Taylor, and J. Sun, Physical Review Letters 113, 144101 (2014).
  * Adhikari  _et al._ (2022) S. Adhikari, J. G. Restrepo, and P. S. Skardal, arXiv preprint arXiv:2208.00909 (2022).
  * Erdős and Rényi (1959) P. Erdős and A. Rényi, Publicationes Mathematicae (Debrecen) 6, 290 (1959).
  * Erdős and Rényi (1960) P. Erdős and A. Rényi, Publ. Math. Inst. Hung. Acad. Sci 5, 17 (1960).
  * Erdős and Rényi (1961) P. Erdős and A. Rényi, Acta Mathematica Hungarica 12, 261 (1961).
  * Barabási and Bonabeau (2003) A.-L. Barabási and E. Bonabeau, Scientific American 288, 60 (2003).
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998).
  * Chutani  _et al._ (2021) M. Chutani, B. Tadić, and N. Gupte, Physical Review E 104, 034206 (2021).
  * Gottwald (2015) G. A. Gottwald, Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 053111 (2015).
  * Pinto and Saa (2015) R. S. Pinto and A. Saa, Physical Review E 92, 062801 (2015).


Composition of q-entropies and hyperbolic orthogonality

NIKOLAOS KALOGEROPOULOS †

Center for Research and Applications
of Nonlinear Systems (CRANS),
University of Patras, Patras 26500, Greece.

Abstract

We point out that the q-entropy composition for independent events has exactly
the same form as the Pythagorean theorem in hyperbolic geometry. We justify
the formal relation of hyperbolic geometry with the q-entropy through the use
of the $\kappa$-entropy, which is directly related to the hyperboloid model of
hyperbolic space. We comment on the relation between orthogonality in this
form of the Pythagorean theorem and the independence of the probability
distributions appearing in the q-entropy composition through the use of the
Dvoretzky-Rogers lemma.

Keywords: q-entropy, $\kappa$-entropy, entropy composition, Pythagorean
theorem, hyperbolic geometry.

† Electronic mail<EMAIL_ADDRESS>

##  Introduction

For more than three decades, the q-entropy (also called “Tsallis entropy”)
[1], which is a one-parameter family of entropies, has been a subject of
investigation in Statistical Mechanics. Despite some progress, it is probably
fair to say that its dynamical foundations and domain of applicability in
Statistical Physics, if any, are still shrouded in mystery. As a result,
speculation abounds, and people have turned mostly to data fittings in an
attempt to justify the use of q-entropy, primarily through its extremal
distributions under some constraints, the q-exponentials [1]. One would,
ideally, like to have a microscopic basis for the use of the q-entropy which
would crucially differentiate it from the well-established
Boltzmann/Gibbs/Shannon (BGS) entropy.

A way to determine differences between the BGS and the q-entropy is to compare
their axiomatic formulations [2], [3], [4], [5]. From this viewpoint, the
major difference, at the axiomatic level, between the q- and the BGS entropy
is in their addition/composition property for probabilistically independent
occurrencies [1].

Several authors have used, in elucidating aspects of q-entropy, this axiomatic
approach as starting point, or the more general but related approach based on
the algebraic properties of composition of entropies, as for instance [6],
[7], [8], [9]. Our investigation is somewhat related, but quite different in
that it relies heavily on geometric, rather than on algebraic, structures
related to the q-entropy in order to elucidate its composition law. We follow
a different path from information geometry [10], as we do not rely on
information theory or statistics in any way, even though one can clearly
discern some rare similarites between our treatments.

We point out that there is a formal analogy between the Pythagorean theorem in
hyperbolic space and the q-entropy composition for independent occurrencies.
We trace back the origin of such a “hyperbolic behavior” in the composition of
the q-entropy to the relation between the q-entropy and the $\kappa$-entropy
[11]. This latter entropy has been motivated by and relies on properties of
the Lorentz transformations [12], [13], [14], [15] which are in turn related
to hyperbolic spaces through the hyperboloid model of hyperbolic geometry.

In Section 2, we derive the expression of the Pythagorean theorem for
hyperbolic spaces which we need in this work. In Section 3, we comment on the
relation between the q- and the $\kappa$\- entropies which justify the
hyperbolic nature of the q-entropy composition law. In Section 4, we present
formal statements mainly from Convex Geometry, which attempt to justify the
relation between orthogonality and probabilistic independence. Section 5
contains further comments related to this work and a general outlook toward
independence and additivity.

##  The Pythagorean theorem for hyperbolic spaces

There is little doubt that the Pythagorean theorem is probably the most
important metric relation in all of Euclidean Geometry. It forms the basis for
its extensions in numerous directions. The most relevant for Physics have been
the quadratic forms determining infinitesimal distances in Riemannian and
semi-Riemannian geometry, the former of which is used in systems requiring
positive-definite metrics where the latter are usually employed in relation to
Special or General Relativity or theories which are based on or extend their
formalism.

What is probably less known is that the Pythagorean theorem has an analogue in
hyperbolic geometry. This is somewhat surprising as most people tend to think
about the Pythagorean theorem as the staple of Euclidean metrics, something
however which does not preclude a similar relation to hold for other
geometries. From the viewpoint of Riemannian geometry, the hyperbolic space is
the unique simply connected $n$-dimensional manifold with constant sectional
curvature $-1$. It should be mentiorned that hyperbolic geometry arose in the
early 19th century in a synthetic/axiomatic fashion as an attempt to determine
whether Euclid’s fifth axiom was independent or could be derived from the
first four, a quest which lasted almost two millenia [16]. One can appreciate
the strong visualization and excellent geometric intuition provided about
hyperbolic geometry when it is seen in its axiomatic/synthetic form. The
synthetic approach is realized by employing one of the several well-known
models of hyperbolic space [17].

We will not attempt to provide a proof of the hyperbolic form of the
Pythagorean theorem here, since it involves a relatively long process which
even though is very instructive from a geometric viepoint, it offers very
little for our purposes. Moreover such a proof can easily be found in the
print and electronic literature. Instead we are content to start from the
conventionally accepted form of the hyperbolic Pythagorean theorem which
states that for a hyperbolic triangle ABC, with the right angle being at C,

$\cosh(c)=\cosh(a)\cosh(b)$ (1)

where $a$ is the length of the side opposite to angle $A$, and similarly for
$b$ and $c$. Squaring both sides and using the obvious identity

$\cosh^{2}(x)=\sinh^{2}(x)+1$ (2)

after a simplification, we get

$\sinh^{2}(c)=\sinh^{2}(a)+\sinh^{2}(b)+\sinh^{2}(a)\sinh^{2}(b)$ (3)

Setting $d_{a}=\sinh(a),\ d_{b}=\sinh(b),\ d_{c}=\sinh(c)$, we finally get

$d_{c}^{2}=d_{a}^{2}+d_{b}^{2}+d_{a}^{2}d_{b}^{2}$ (4)

which is the sought-after form of the Pythagorean theorem in hyperbolic space
which we will be using in the rest of this work. Upon closer inspection it
will turn out that (4) has the same form as the q-entropy composition property
for independent occurrences as will be pointed out in the next Section.

##  Hyperbolicity of q-entropies through $\kappa$-entropies

The q-entropy (also known as “Tsallis entropy”) [1] for a discrete sample
space $\mathfrak{X}$ whose elements are indexed by a set $I$ with
probabilities $p_{i},\ i\in I$ is defined as

$\mathcal{S}_{q}[\\{p_{i}\\}]\ =\ k_{B}\cdot\frac{1}{q-1}\left(1-\sum_{i\in
I}\ p_{i}^{q}\right)$ (5)

where $k_{B}$ stands for the Boltzmann constant to be set to one for
convenience, and where $q\in\mathbb{R}$ is a constant called the “entropic
parameter”. In practice $q$ is usually restricted to a proper subset of
$\mathbb{R}$ in most cases of physical interest. One immediately sees that for
two independent sets of occurrences $A,B\subset\mathfrak{X}$, namely if

$p_{A\cap B}=p_{A}\cdot p_{B}$ (6)

then

$\mathcal{S}_{q}[p_{A\cap
B}]=\mathcal{S}_{q}[p_{A}]+\mathcal{S}_{q}[p_{B}]+(1-q)\mathcal{S}_{q}[p_{A}]\mathcal{S}_{q}[p_{B}]$
(7)

and defining

$\mathcal{S}^{\prime}_{q}=\frac{1}{1-q}\ \mathcal{S}_{q}$ (8)

the composition property (7) becomes

$\mathcal{S}^{\prime}_{q}[p_{A\cap
B}]=\mathcal{S}^{\prime}_{q}[p_{A}]+\mathcal{S}^{\prime}_{q}[p_{B}]+\mathcal{S}^{\prime}_{q}[p_{A}]\mathcal{S}^{\prime}_{q}[p_{B}]$
(9)

which has exactly the same form as the Pythagorean theorem in hyperbolic space
(4).

One question which naturally arises is whether the q-entropy composition
property (9) has anything of substance to do with the Pythagorean theorem in
hyperbolic space, except for their formal similarity. This question can be
further broken down to two, largely indepedent of each other, questions

  1. 1.

What does the q-entropy have to do with hyperbolic space?

  2. 2.

What does orthogonality have to do with the composition property (9)?

It seems that the fact that the q-entropy composition has a formal similarity
to Pythagorean theorem in hyperbolic space has been pointed out before in the
context of information geometry [10]. But so far as we know there has not been
any attempt to understand _why_ the composition property (9) has the same form
as the hyperbolic version of the Pythagorean theorem (4), and instead it has
been taken at face value. The information geometry practitioners focused
instead on creating and further developing a rather elaborate formalism of
dual Riemannian geometric structures and escort distributions [10]. Therefore
the present work may help address this gap in our understanding of q-entropy,
in our opinion.

We begin by stating that the q-entropy is not a distance function as one might
be tempted to assume at first glance. This is indeed true for most entropic
forms which have been proposed over the years [18], including the BGS entropy,
despite occasional comments in the literature about the opposite. Usually such
entropic forms fail to obey the triangle inequality. Moreover one sees such
entropies as ways of measuring the difference between two probability
distributions, in order to establish coordinate reparametrization invariance,
akin to the Kullback-Leibler, the Bregman etc divergences. Then such entropic
forms fail to obey even the symmetry property in their arguments, which is
required of distance functions. Hence, one should be careful not to take the
analogy too far or draw any invalid conclusions by the formal similarity
between (4) and (9).

The answer to question 1 comes through the relation of the q- and the
$\kappa$\- entropies. The $\kappa$-entropy $\mathcal{R}_{\kappa}[p_{i}]$ is a
one-parameter set of entropies defined [12], [13], [14], [15] for discrete
sets of occurrences $p_{i},\ i\in I$ by

$\mathcal{R}_{\kappa}[\\{p_{i}\\}]\ =\ k_{B}\sum_{i\in I}\
\frac{p_{i}^{1-\kappa}-p_{i}^{1+\kappa}}{2\kappa}$ (10)

The $\kappa$-entropy purports to generalize the BGS entropy to a special
relativistic context and is partially motivated by the $\kappa$ distributions
which have been observed in the heliosphere, and which have been a subject of
intense investigation over several decades [19] especially in plasma
Astrophysics.

At this point, and for comparison purposes to the q-entropy composition (9)
only, we digress to state that the composition property of the
$\kappa$-entropy for independent probabilities is not as easy to obtain as the
composition of q-entropies. It was eventually found to be [20]

$\mathcal{R}_{\kappa}[p_{A\cap B}]\ =\
\mathcal{R}_{\kappa}[p_{A}]\left(\frac{1}{\gamma_{\kappa}}\mathcal{R}_{\kappa}\left[p_{B}/\eta_{\kappa}\right]-\mathcal{R}_{\kappa}[p_{B}]\right)+\mathcal{R}_{\kappa}[p_{B}]\left(\frac{1}{\gamma_{\kappa}}\mathcal{R}_{\kappa}\left[p_{A}/\eta_{\kappa}\right]-\mathcal{R}_{\kappa}[p_{A}]\right)$
(11)

where

$\gamma_{\kappa}=\frac{1}{\sqrt{1-\kappa^{2}}}$ (12)

is a rapidity-like factor, and

$\eta_{\kappa}=\left(\frac{1+\kappa}{1-\kappa}\right)^{\frac{1}{2\kappa}}$
(13)

The $\kappa$-entropy has a functional form which is directly descended from
the Lorentz transformations in Special Relativity [12], [13]. Since the
Lorentz group is a normal subgroup of the Poincaré group, which is the group
of isometries of Minkowski space-time, following the viewpoint of F. Klein’s
Erlangen Programme, one can see that properties of the $\kappa$-entropy are
formally related to properties of Minkowski space-time.

The relation between the q- and $\kappa$-entropies is readily seen to be

$\mathcal{R}_{\kappa}[p_{i}]=\frac{1}{2}\left(\mathcal{S}_{1+\kappa}[p_{i}]+\mathcal{S}_{1-\kappa}[p_{i}]\right)$
(14)

This relation can be inverted to give the q-entropies in terms of the
$\kappa$-entropies, but the explicit form is not necessary for our argument.
All that matters is that there is a simple relation between these two entropic
forms. As was mentioned before, the $\kappa$-entropy reflects properties of
Minkowski space-time. To be more concrete, let
$(x_{1},x_{2},\ldots,x_{n+1})\in\mathbb{R}^{n+1}$ endowed with the quadratic
form

$ds^{2}=dx_{1}^{2}+dx^{2}_{2}+\ldots+dx_{n}^{2}-dx_{n+1}^{2}$ (15)

Consider the $n$-dimensional hypersurface defined as the “sphere of unit
imaginary radius” namely

$\mathbb{H}^{n}=\\{(x_{1},x_{2},\ldots,x_{n},x_{n+1})\in\mathbb{R}^{n+1}:\ \
x_{1}^{2}+x_{2}^{2}+\ldots x_{n}^{2}-x_{n+1}^{2}=-1,\ \ x_{n+1}>0\\}$ (16)

This single branch hyperboloid has a space-like, namely a Riemannian, metric
and a careful examination of its properties proves that it is actually
isometric to the $n$-dimensional hyperbolic space $\mathbb{H}^{n}$, hence the
symbol we used for it in (16). This is actually the well-known hyperboloid
model [16], [17] of hyperbolic space.

To finish the argument, the hyperbolic space appears as a subset of Minkowski
space-time, the latter of which has some common properties, formaly of course,
with the $\kappa$-entropy. On the other hand, the q-entropy is a combination
of $\kappa$-entropies. Therefore, it is not totally surprising that some of
the formal properties of the hyperbolic space are inherited by the q-entropy.
The composition property of the q-entropy is one of them. This justifies the
formal resemblance of the hyperbolic Pythagorean theorem (4) and the q-entropy
composition property (9), so far as the word “hyperbolic” is concerned. The
justification of the word “orthogonal” will be dealt with in the next Section.

Before finishing this Section, it may worth making a couple of comments. In
some of our prior works [21], [22] we noticed that the formal resemblance
between the q-entropy composition property and elements of sub-Riemannian
geometry related to the Heisenberg, or more generally, to nilpotent groups.
The problem with Carnot groups, or more generally Carnot-Carathéodory spaces
[23] is that they have a lot of structure, so in a sense they are very
special. On the other hand, manifolds or more general spaces of negative
curvature are plentiful. An example of this is dramatically expressed by W.
Thurston’s work on geometrization [16] which showed that among all
3-dimensional manifolds most of them are hyperbolic. Let us not forget that
something similar occurs with surfaces (manifolds of dimension 2) since
surfaces of any genus greater than one admit a hyperbolic metric. But these
are low-dimensional results which are not of great interest if one works with
configuration or phase spaces of systems with many degrees of freedom, which
are of main interest in Statistical Mechanics.

A more pertinent result in the use of hyperbolic structures valid for
manifolds of high dimension though can be found in [24] where it is proved
that a closed smooth manifold of any dimension can be given a metric of
negatively pinched curvature close to a hyperbolic one, even without
singularities. This when combined with the Mostow rigidity, which states that
two hyperbolic manifolds of dimension greater than two which have isomorphic
fundamental groups are actually isometric [25], makes us suspect that
hyperbolicity in high dimensional manifolds may be more frequent than one
might have naively thought.

The result of all this is that from the viewpoint of geometry, the composition
of q-entropies which has the same form as the hyperbolic Pythagorean theorem
may be quite a bit more frequent than one might believe at a first glance.
Since this composition is the distinguishing property between the BGS and
q-entropies, the latter may be applicable to more physical systems than might
have initially been suspected.

##  Orthogonality and probabilistic independence

The question we wish to address in this Section, is question 2 in the list of
the previous Section, namely, what is the reason for the formal orthogonality
between the joint probabilty distribution and its marginals appearing in the
composition law (9). Since the standard hyperbolic metric and the Euclidean
one are conformal [16], [17] the angles between lines in the hyperbolic plane
$\mathbb{H}^{2}$ and in the Euclidean plane $\mathbb{R}^{2}$ are equal. So it
is sufficient to make an argument for orthogonality for probability
distributions in $\mathbb{R}^{2}$.

Let us consider as sample space a finite set
$\mathfrak{X}=\\{a_{1},a_{2},a_{3},\ldots,a_{n}\\}$ in order to avoid details
related to infinities and to continuous distributions that would obscure the
argument by making it far more technical but would add nothing of essence.
Assign probabilities $p_{1},\ p_{2},\ p_{3},\ldots$ to each element $a_{1},\
a_{2},\ a_{3},\ldots,a_{n}$ of $\mathfrak{X}$ respectively. Obviously

$\sum_{i=1}^{n}\ p_{i}=1,\ \ \ p_{i}\geq 0$ (17)

These probabilities form the vertices of a probability simplex
$\Delta^{n-1}\subset\mathbb{R}^{n}$ when the value of each one of them is
placed on the axes of an orthonormal, Cartesian, coordinate system with
respect to the Euclidean metric of $\mathbb{R}^{n}$. Consider
$A,B\subset\mathfrak{X}$. Then the outcomes of $A$ and $B$ are
probabilistically independent if (6) holds. Does the probabilistic
independence imply orthogonality? The answer, in general, is negative.
However, these two concepts start having some common points if we consider the
space $\mathfrak{X}$ having a large cardinality as is the case of interest in
Statistical Mechanics.

To see how the relation between probabilistc independence and orthogonality
occurs, one can start by moving the barycenter of $\Delta^{n-1}$ to the origin
of $\mathbb{R}^{n}$ by a translation, or more generally an affine
transformation. It may be instructive, and more familiar from a physical
viewpoint, to consider the probabilities $p_{i},\ i=1,\ldots,n$ as unit point
masses placed on the respective axes at distances $p_{i}$ from the origin of
the Cartesian coordinate system of $\mathbb{R}^{n}$ which we use. After the
translation, the resulting simplex, let us call it $\Delta_{\ast}^{n-1}$, is
the boundary of a convex body $\mathcal{K}$ whose center is now at the origin,
and whose vertices have coordinates $p_{i}^{\prime},\ i=1,2,\ldots,n$, which
means that

$\sum_{i=1}^{n}\ p_{i}^{\prime}=0$ (18)

Let us indicate the Euclidean inner product by $\langle\ ,\ \rangle$. Let
$e_{1},\ldots,e_{n}$ indicate the unit vectors along the respective coordinate
axes of the probabilities $p_{1}^{\prime},\ldots,p_{n}^{\prime}$. Then, for
vectors $u,v,w\in\mathbb{R}^{n}$ let us indicate the tensor product by
$u\otimes v$, namely

$(u\otimes v)(w)=\langle v,w\rangle u$ (19)

In case the vector $u$ has unit length, then (19) implies that $u\otimes u$ is
the orthogonal projection to the linear subspace spanned by $u$. Going back to
the physical interpretation of the $p_{i}$, we can see that the “inertia
tensor” of the simplex $\Delta_{\ast}^{n-1}$ is not only diagonal, but even
more so it is the identity, namely

$\sum_{i=1}^{n}\ p_{i}^{\prime}\ e_{i}\otimes e_{i}=1_{n}$ (20)

where $1_{n}$ indicates the identity map on $\mathbb{R}^{n}$. Then a theorem
of F. John [26], states that there is an ellipsoid of maximum volume contained
within $\Delta_{\ast}^{n-1}$ which is actually a Euclidean ball. Then, the
inertia condition (20) implies that such $e_{i}$ would behave like an
orthonormal set since the inner product of any two vectors can be written as

$\langle v,w\rangle=\sum_{i=1}^{n}\ p_{i}^{\prime}\ \langle
v,e_{i}\rangle\langle e_{i},w\rangle$ (21)

for $v,\ w\in\mathbb{R}^{n}$.

It turns out that John’s theorem is valid for a general convex body
$\mathcal{K}$ obeying conditions (18) and (20) and not just for a simplex like
$\Delta_{\ast}^{n-1}.$ In that case, one takes as unit vectors $f_{i},\
i=1,2,\ldots,n$ the normalizations of the vectors which are intersections
$\Delta_{\ast}^{n-1}\cap B_{n}$ where $B_{n}$ indicates the maximum volume
Euclidean ball. To be consistent with the formulation of the previous
paragraph, the coefficients $c_{i}^{\prime},\ i=1,2,\ldots,n$ corresponding to
$p_{i}^{\prime}$ in (18), and (20) are all taken to be positive. But then it
is no longer true that the unit vectors $f_{i},\ i=1,2,\ldots,n$ are
orthogonal. The fact that it is possible to find an orthonormal basis
$\\{z_{i},\ i=1,\ldots,n\\}$ which is not “too far” from the basis $\\{f_{i},\
i=1,2,\ldots,n\\}$ of $\mathbb{R}^{n}$ is the content of the Dvoretzky-Rogers
lemma [26] which states that there is $\\{u_{i},\
i=1,2,\ldots,n\\}\subset\\{f_{i},\ i=1,2,\ldots,n\\}$, where $u_{i}$ belongs
to the linear span of $z_{1},z_{2},\ldots,z_{i}$, such that

$\sqrt{\frac{n-i-1}{n}}\leq\langle u_{i},z_{i}\rangle\leq 1,\ \ \ \ \ \
i=1,2,\ldots n$ (22)

Consider now two subsets $A,B\subset\mathfrak{X}$ and $A\cap B$ be their
intersection as in (6). Assume that their intersection is a connected set, for
simplicity. If it is not, then one can apply this reasoning to each connected
component separately. The intersection of $A$ and $B$ will be a convex subset
of the simplex $\Delta^{n-1}$. Then, this convex subset upon a possible affine
transformation, satisfies the conditions of John’s theorem, therefore it has a
maximum volume inscribed ball for which the Dvoretzky-Rogers lemma applies.
Hence, there are nearly orthonormal vectors forming a basis of $A\cap B$. We
can extend this basis to an orthonormal basis of $A$ and $B$ respectively by a
judicious choice of basis vectors for each one set, regardless of what choice
we make in the other, since $A$ and $B$ do not intersect outside of $A\cap B$.
In essence we can extend the near orthonormal basis provided by the Dvoretzky-
Rogers lemma to an almost orthonormal basis of $A$ and $B$. Hence the sets $A$
and $B$ appearing in (6) are almost orthogonal to $A\cap B$. The word “almost”
appearing everywhere above disappears when one takes $n\rightarrow\infty$ as
is assumed in many cases in Statistical Physics when one take the
thermodynamic limit. Then (22) implies that the majority of vectors $u_{i}$
and $z_{i}$ become collinear hence $A\cap B$ is the orthogonal intersection of
sets $A$ and $B$ for the vast majority of their elements. Hence one can state,
in this particular context, that probabilistic independence implies
orthogonality, as we wanted to ascertain, answering question 2 of the previous
Section.

##  Conclusions and discussion

In this work we pointed out that the q-entropy composition property has the
exact same form as the Pythagorean theorem in hyperbolic space. We justified
the appearance of hyperbolicity in this composition by ascribing it to the
simple relation between the q- and the $\kappa$\- entropies where the latter
are based on Lorentz trasformations, hence on Minkowski space-time, and by
subsequently using the hyperboloid model of hyperbolic space to reach this
conclusion. We justified the appearance of orthogonality through the use of
the Dvoretzky-Rogers lemma applied to convex sets in high dimensional linear
spaces.

In a recent work [27], motivated by the appearence of $\kappa$ distributions
in astrophysical plasmas, the departure of the q-entropy composition from a
simple addition seems to be considered as a drawback, as this is how one may
interpret the word “defect” in the title of that work. If one accepts the
ordinary/usual addition as fundamental and all other forms of addition are
compared to it, then the authors may be right. From the same perspective, if
one accepts the BGS entropy and its composition as the “reference definition”
and “reference composition” relation of entropies to which all other entropic
forms should be compared, and be considered deviations from, then the word
“defect” may be applicable.

What the present work shows, in our opinion, is that the q-entropy, and we
suspect many other entropies, are fundamentally different from the BGS one, as
much as hyperbolic geometry is different from the Euclidean one. The usual
addition is not “better” or “more correct” than a generalized addition such as
the one expressed in (4), (9); they are just different. We should not consider
one as a variation of the other. They apply to different systems and serve
different purposes. This is clearly along the lines of [1] who repeatedly
stresses that the q-entropy is not just a variation of the BGS entropy, as it
might be inferred from its functional form, but a fundamentally different
entropy which is applicable to different systems, and which happens to have as
an appropriate limit the BGS entropy. In our work we show that the BGS and
q-entropies have the same relation as Euclidean to hyperbolic geometries; they
are just different.

Given all this, the major, in our opinion, question remains; to which physical
systems are the q- and $\kappa$\- entropies applicable and what is the
underlying dynamics that such entropies describe at the macroscopic level? We
believe that despite its difficulty, this may be a worthy subject for a future
line of investigation.

## References

  * [1] C. Tsallis, _Introduction to Nonextensive Statistical Mechanics: Approaching a Complex World, Second Edition_ , Springer Nature, Switzerland (2023).
  * [2] R.J.V. Santos,  _Generalization of Shannon’s theorem for Tsallis entropy_ , J. Math. Phys. 38(8), 4104-4107 (1997).
  * [3] S. Abe,  _Axioms and uniqueness theorem for Tsallis entropy_ , Phys. Lett. A 271(1-2) 74-79 (2000).
  * [4] S. Abe,  _General pseudoadditivity of composable entropy prescribed by the existence of equilibrium_ , Phys. Rev. E 63 061105 (2001).
  * [5] A. Enciso, P. Tempesta,  _Uniqueness and characterization theorems for generalized entropies_ , J. Stat. Mech. Theory Exp. 2017 123101 (2017).
  * [6] P. Tempesta,  _Beyond the Shannon-Khinchin formulation: The composability axiom and the universal group property_ , Ann. Phys. 365, 180-197 (2016).
  * [7] P. Tempesta,  _Formal Groups and Z-Entropies_ , Proc. Roy. Soc. A 472, 20160143 (2016).
  * [8] H.J. Jensen, P. Tempesta,  _Group Entropies: From Phase Space Geometry to Entropy Functionals via Group Theory_ , Entropy 20, 804 (2018).
  * [9] H.J. Jensen, P. Tempesta,  _Group Structure as a Foundation for Entropies_ , Entropy 26, 266 (2024).
  * [10] S. Amari,  _Information Geometry and Its Applications_ , Applied Mathematical Sciences 194, Springer, Japan (2020).
  * [11] T. Wada, A. M. Scarfone,  _On the Kaniadakis distributions applied in statistical physics and natural sciences_ , Entropy 25(2), 292 (2023).
  * [12] G. Kaniadakis,  _Statistical mechanics in the context of special relativity_ , Phys. Rev. E 66, 056125 (2002).
  * [13] G. Kaniadakis,  _Statistical mechanics in the context of special relativity II_ , Phys. Rev. E 72, 036108 (2005).
  * [14] G. Kaniadakis,  _Power-law tailed statistical distributions and Lorentz transformations_ , Phys. Lett. A 375, 356-359 (2011).
  * [15] G. Kaniadakis,  _Relativistic Roots of $\kappa$-entropy_, Entropy 26(5), 406 (2024).
  * [16] W.P. Thurston,  _Three-Dimensional Geometry and Topology, Volume 1_ , S. Levy (Ed.), Princeton Math. Series (35, Princeton Univ. Press, Princeton NJ, USA (1997).
  * [17] J.W. Cannon, W.J. Floyd, R. Kenyon, W.R. Parry,  _Hyperbolic Geometry_ , pp. 59-115 in _Flavors of Geometry_ , S. Levy (Ed.), Math. Sci. Res. Inst. Publ. 31, Berkeley, USA (1997).
  * [18] V.M. Ilić, J. Korbel, S. Gupta, A.M. Scarfone,  _An overview of generalized entropic forms_ , Europhys. Lett. 133, 50005 (2021).
  * [19] K. Scherer, E. Hudisic, M. Lazard, H. Fichtner,  _The $\kappa$-cookbook: a novel generalizing approach to unify $\kappa$-like distributions for plasma partice modeling_ Month. Notices Roy. Astron. Soc. 497, 1738-1756 (2020).
  * [20] G. Kaniadakis, A.M. Scafone, A.Sparavigna, T. Wada,  _Composition of $\kappa$-entropy for statistically independent systems_, Phys. Rev. E 95, 052112 (2017).
  * [21] N. Kalogeropoulos,  _Tsallis entropy composition and the Heisenberg group_ , Int. J. Geom. Methods Mod. Phys. 10, 1350032 (2013).
  * [22] A.J. Creaco, N. Kalogeropoulos, _Nilpotence in Physics: the case of Tsallis entropy_ , J. Phys. Conf. Ser. 410, 012148 (2013)
  * [23] M. Gromov,  _Carnot-Carathéodory spaces seen from within_ , pp. 79-233 in  _Sub-Riemannian Geometry_ , A. Bellaiche, J. -J. Risler (Eds.), Progress in Mathematics 144, Birkhäuser, Basel, Switzerland (1996).
  * [24] P. Ontaneda,  _Riemannian hyperbolization_ , Publ. Math. IHES 131, 1-72 (2020).
  * [25] G.D. Mostow  _Strong Rigidity of Locally Symmetric Spaces_ , Princeton Univ. Press, Princeton NJ, USA (1973).
  * [26] S. Artstein-Avidan, A. Giannopoulos, V.D. Milman,  _Asymptotic Geometric Analysis, Part I_ , Math. Surveys Monog. 202, Amer. Math. Soc., Providence, RI, USA (2015).
  * [27] G. Livadiotis, D.J. McComas,  _Entropy defect: Algebra and thermodynamics_ , Europhys. Lett. 144, 21001 (2023).