# Back to the Roots:
Tree-Based Algorithms for Weakly Supervised Anomaly Detection

Thorben Finke<EMAIL_ADDRESS>Institut für Theoretische
Teilchenphysik und Kosmologie, RWTH Aachen, 52074 Aachen, Germany  Marie Hein
<EMAIL_ADDRESS>Institut für Theoretische Teilchenphysik und
Kosmologie, RWTH Aachen, 52074 Aachen, Germany  Gregor Kasieczka
<EMAIL_ADDRESS>Institut für Experimentalphysik, Universität
Hamburg, 22761 Hamburg, Germany  Center for Data and Computing in Natural
Sciences (CDCS), 22607 Hamburg, Germany  Michael Krämer<EMAIL_ADDRESS>aachen.de Institut für Theoretische Teilchenphysik und Kosmologie, RWTH
Aachen, 52074 Aachen, Germany  Alexander Mück<EMAIL_ADDRESS>Institut für Theoretische Teilchenphysik und Kosmologie, RWTH Aachen, 52074
Aachen, Germany  Parada Prangchaikul  Institut für Experimentalphysik,
Universität Hamburg, 22761 Hamburg, Germany  Tobias Quadfasel
<EMAIL_ADDRESS>Institut für Experimentalphysik, Universität
Hamburg, 22761 Hamburg, Germany  David Shih<EMAIL_ADDRESS>NHETC,
Dept. of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854, USA
Manuel Sommerhalder<EMAIL_ADDRESS>Institut für
Experimentalphysik, Universität Hamburg, 22761 Hamburg, Germany

###### Abstract

Weakly supervised methods have emerged as a powerful tool for model-agnostic
anomaly detection at the Large Hadron Collider (LHC). While these methods have
shown remarkable performance on specific signatures such as di-jet resonances,
their application in a more model-agnostic manner requires dealing with a
larger number of potentially noisy input features. In this paper, we show that
using boosted decision trees as classifiers in weakly supervised anomaly
detection gives superior performance compared to deep neural networks. Boosted
decision trees are well known for their effectiveness in tabular data
analysis. Our results show that they not only offer significantly faster
training and evaluation times, but they are also robust to a large number of
noisy input features. By using advanced gradient boosted decision trees in
combination with ensembling techniques and an extended set of features, we
significantly improve the performance of weakly supervised methods for anomaly
detection at the LHC. This advance is a crucial step towards a more model-
agnostic search for new physics.

##  I Introduction

The search for new physics at the Large Hadron Collider (LHC) requires
powerful and model-agnostic anomaly detection methods [1, 2]. Weakly
supervised approaches [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] have proven
effective in identifying specific signatures, such as di-jet resonances, where
signal and background can be distinguished using a limited set of hand-crafted
features. In these approaches, a classifier is generally trained to
distinguish a signal region from a background-dominated control region or
background template. In interesting applications, the signal fraction in the
signal region is small, and the classifier needs to identify a small number of
signal events in an overwhelming background using only the labels for the
signal region and background template. Thus, even a relatively simple
supervised classification problem becomes increasingly challenging in the
weakly supervised setting as the signal fraction decreases. Extending these
methods to a wider range of new physics models without prior knowledge of the
essential discriminative features requires handling a larger number of
potentially noisy input features. This is a significant challenge as the
performance of traditional deep neural networks can degrade under such
conditions.

In this paper we explore the use of boosted decision trees (BDTs) as
classifiers in weakly supervised anomaly detection at the LHC. Boosted
decision trees are known for their strength in analyzing tabular data [14]. In
particular, on small and medium training sets, BDTs outperform deep learning
methods, see for example the reviews [15, 14] and references therein. In the
case of weakly supervised anomaly detection, only the signal events separate
the two classes. Therefore, while the total training set may be large, the
small number of signal data may favor BDTs. In addition, BDTs are generally
less affected by uninformative features [14] (see also Appendix B). In a high-
dimensional input feature space for truly model-agnostic anomaly detection,
many of the input features will inevitably be uninformative for a given signal
model. As such, robustness to uninformative features is an important property
of a model-agnostic method. This robustness and the efficient training and
evaluation times make BDTs an interesting alternative to deep neural networks.

An important aspect of improving performance in weakly supervised learning is
ensembling, which harnesses the collective power of multiple classifiers to
achieve improved classification accuracy. Although ensembling can be applied
to various machine learning methods, the fast training and evaluation time of
BDTs allows ensembling to be used most efficiently.

We study the potential superiority of boosted decision trees over deep neural
networks for weakly supervised anomaly detection using modern gradient
boosting techniques and ensembling methods. Having obtained classifiers that
scale well to a larger number of input features, we investigate whether weak
supervision techniques can take advantage of additional quantities describing
the substructure of jets.

The remainder of this paper is structured as follows: Section II describes
weak anomaly detection and modern boosted decision trees; Section III
introduces the LHC Olympics (LHCO) benchmark dataset used; Section IV
discusses the empirical results of our studies; and Section V presents a
summary and outlook for future work. In Appendix A we compare the performance
of different BDT architectures, in Appendix B we analyze the robustness of
BDTs with respect to noisy features, and in Appendix C we discuss the impact
of ensembling in more detail.

##  II Methods

###  II.1 Optimal anomaly score and its approximations

According to the Neyman-Pearson lemma [16], the most powerful model-agnostic
anomaly score is the likelihood ratio between data and background events

$R_{\rm
optimal}(\boldsymbol{x})=\frac{p_{\mathrm{data}}(\boldsymbol{x})}{p_{\mathrm{bg}}(\boldsymbol{x})}\,,$
(1)

where $\boldsymbol{x}$ is a set of features that describe the events (e.g.
kinematics, high-level variables etc.). This is monotonic with the signal-to-
background likelihood ratio for any, not necessarily specified, signal model.
We will refer to (1) as the optimal anomaly score.

In practice, in any realistic high energy physics context, it is not possible
to achieve the optimal anomaly score, for several reasons. First, the true
background density is generally unknown—our simulations of the Standard Model
and the detector response are imperfect, so at best we have an approximation
to $p_{\rm bg}({\boldsymbol{x}})$. Furthermore, the densities themselves are
intractable and at best we have samples from them. Finally, the data density
is also unknown and must also be approximated from samples.

Existing methods [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] attempt to approximate
the optimal anomaly score through binary classifiers that are trained to
distinguish between samples drawn from the events within a signal region (SR)
and samples drawn from a data-driven background template. To date, this
classification task has been exclusively performed using deep neural networks.
This work will instead focus on the benefits of using boosted decision trees.

So far, methods have primarily focused on the case of resonant anomaly
detection, where the SR is a window in a resonant variable such as dijet
invariant mass, and the background template is derived somehow from
neighboring control regions. However, the idea of the optimal anomaly score
and approximating it with a classifier is more general and not necessarily
limited to resonant anomaly detection, see e.g. [17].

In order to separate out the different approaches to deriving background
templates from the performance of the binary classifier itself—the issue of
interest in this work—we will focus on a version of the optimal anomaly score,
previously termed the idealized anomaly detector (IAD) in [9], which is a
binary classifier trained on SR data and a perfectly modeled background
template. In other words, we presume that the background events in the SR and
the background template events are drawn from the exact same distribution.
Focusing on the IAD in our work ensures that any relative improvement in
significance between approaches is solely due to an improvement in the
classification task itself.

###  II.2 Boosted decision trees

Gradient Boosted Decision Trees (GBDTs) are established as the best performing
methods on tabular data. In recent years, LightGBM [18] has become the state
of the art, mainly due to its fast training and evaluation times achieved by
histogramming the inputs [19]. We use the HistGradientBoostingClassifier
implementation in scikit-learn [20], which is based on LightGBM, because it is
easy to use and performs well. We use default hyperparameters unless stated
otherwise, i.e. a learning rate of $0.1$, a maximum number of leaf nodes per
tree of 31, a maximum of 255 bins per feature, and early stopping with a
patience of 10 iterations. The maximum number of iterations is increased with
respect to the default to 200 to ensure that all training is stopped at the
minimum validation loss and not at the maximum number of trees. The GBDT
implementation uses subsampling for individual trees such that the predictions
can vary for difficult classification tasks like weakly supervised anomaly
detection with small signal fractions.

To further stabilize and improve performance, our classifier uses an ensemble
of $N$ independent runs of the BDT by averaging the $N$ predictions. For this
ensemble, the split between training and validation is randomized. This
results in a further increase in performance compared to a fixed split. Unless
otherwise specified, our classifier uses an ensemble of $N=50$ independent
runs.

We compare the performance of our default GBDT implementation with other BDT
architectures in Appendix A.

###  II.3 Neural Network

We use the neural network architecture of Ref. [9], however we have
implemented it in Tensorflow [21] using Keras [22]. The network is a fully
connected neural network with a binary cross-entropy loss. It consists of 3
hidden layers, each containing 64 nodes and using ReLU activation. The output
layer employs the softmax activation. Training is conducted over 100 epochs
using the Adam [23] optimizer, a learning rate of $10^{-3}$, and a batch size
of 128. Early stopping with a patience of 10 epochs is used in analogy to the
BDT. The NN classifier also applies an ensemble of $N=50$ independent runs
with a randomized validation split.

##  III data

###  III.1 The dataset

Our studies are performed using the R&D dataset [24] from the LHC Olympics
2020 (LHCO) [1]. For the most part, we focus on the original R&D signal model,
$Z^{\prime}\rightarrow X(\rightarrow qq)Y(\rightarrow qq)$, with masses
$m_{Z^{\prime}}=$3.5\text{\,}\mathrm{TeV}$$,
$m_{X}=$500\text{\,}\mathrm{GeV}$$ and $m_{Y}=$100\text{\,}\mathrm{GeV}$$. In
this signal model, the boosted dijets have 2-prong substructure. In Section
IV.4, we further investigate an alternative signal model also included with
the LHCO R&D dataset, where the $X$ and $Y$ particles decay to 3-prong
substructure, i.e. $Z^{\prime}\rightarrow X(\rightarrow qqq)Y(\rightarrow
qqq)$ [24].

The dataset contains $1\,000\,000$ QCD dijet events for the background and
$100\,000$ signal events, both simulated with Pythia 8 [25, 26] and Delphes
3.4.1 [27]. Reconstructed particles are clustered into jets with FastJet [28]
using the anti-$k_{T}$ algorithm [29] with a distance parameter of $R=1$.
Additionally, a trigger of $p_{\mathrm{T}}>$1.2\text{\,}\mathrm{TeV}$$ must be
passed for all events.

Of the background events, approximately $120\,000$ fall in the SR between
$3.3\text{\,}\mathrm{TeV}$ and $3.7\text{\,}\mathrm{TeV}$. We also use
$612\,858$ additional QCD background events in the SR, which can be found here
[30], and which have been employed previously in Ref. [9]. These events are
used for the background template for the IAD, for the supervised classifier,
and for the construction of the test dataset.

Unless otherwise specified, $1000$ signal events are injected into the
$1\,000\,000$ background events of the R&D dataset. As in Ref. [9], the same
signal events are used in different classifiers unless $S/B$ is varied. This
allows for better comparability between BDT and NN results, as the high
computational cost of the NN makes it prohibitive to train multiple
classifiers for each setup (see also Section III.3). This results in 772
signal events in the SR and corresponds to $S/B=6\times 10^{-3}$ and
$S/\sqrt{B}=2.2$.

For training and validation, a total of approximately 272k SR background
events in the background template are trained against 120k SR data events for
the IAD and 54k SR signal events for the supervised classifier. In both cases
a 50 % validation split is used, which is randomized between different runs as
mentioned above. Testing is performed on 340k SR background and 20k SR signal
events.

###  III.2 Features

The features used for training are similar to those used in Ref. [9]. We
select the two jets with the highest $p_{T}$. As features we use the invariant
mass of the lighter of the two jets, $m_{J_{1}}$, the difference in jet mass
between the two jets, $\Delta m_{J}=m_{J_{2}}-m_{J_{1}}$, and several features
based on their n-subjettiness [31, 32]. The original R&D dataset consisted of
$\tau_{1}$, $\tau_{2}$ and $\tau_{3}$ subjettiness features computed with
angular weighting parameter $\beta=1$. In this work, we also computed
additional subjettiness features ($\tau_{n}$ up to $n=9$) for both jets, and
using three different values for $\beta$ ($\beta=0.5$, $\beta=1$ and
$\beta=2$). Thus, a total of $54$ subjettiness features are available
considering both jets. To investigate the impact of the subjettiness features
on the classification performance, we consider various feature sets as listed
in Tab. 1. The subjettiness ratios are defined as
$\tau_{ij}\equiv\tau_{i}/\tau_{j}$.

Name | # features | Features
---|---|---
Baseline | 4 | $\\{m_{J_{1}},\,\Delta m_{J},\,\tau_{21}^{\beta=1,J_{1}},\,\tau_{21}^{\beta=1,J_{2}}\\}$
Extended 1 | 10 | $\\{m_{J_{1}},\,\Delta m_{J},\,\tau_{N,N-1}^{\beta=1,J_{1}},\,\tau_{N,N-1}^{\beta=1,J_{2}}\\}$
for $2\leq N\leq 5$
Extended 2 | 12 | $\\{m_{J_{1}},\,\Delta m_{J},\,\tau_{N}^{\beta=1,J_{1}},\,\tau_{N}^{\beta=1,J_{2}}\\}$
for $N\leq 5$
Extended 3 | 56 | $\\{m_{J_{1}},\,\Delta m_{J},\,\tau_{N}^{\beta,J_{1}},\,\tau_{N}^{\beta,J_{2}}\\}$
for $N\leq 9$ and $\beta\in\\{0.5,1,2\\}$
Table 1: Subjettiness feature sets considered for training. Full training
feature sets always include $m_{J_{1}}$ and $\Delta m_{J}$ as well. Details of
the observables are given in the text.

###  III.3 Metric

Throughout this work we use the significance improvement characteristic (SIC)
as our main metric. We define

$\textrm{SIC}=\frac{\epsilon_{S}}{\sqrt{\epsilon_{B}}}\,$ (2)

as the ratio of the fraction of correctly identified signal events
$\epsilon_{S}$ to the square root of the fraction of background events
misidentified as signal $\sqrt{\epsilon_{B}}$. This is an approximation of the
expected sensitivity gain over an inclusive search, assuming a dominant
Poisson uncertainty in the background event count. For very low background and
signal efficiencies, the statistical error of the SIC value increases, causing
large fluctuations. As the uncertainty of the background efficiency dominates,
a cut-off on the relative statistical error of the background efficiency at
20% is introduced.

For the BDT, all results show the median and 68 % confidence intervals of 10
independent classifiers each employing our ensembling procedure. As training
the NNs is computationally intensive, the results of only one NN classifier
are shown. However, for the large ensemble of $N=50$ independent trainings per
classifier, the variance is expected to be small; this has been verified by
examining 10 different NN classifiers with $N=10$ for all results presented
here.

##  IV Results

###  IV.1 Stability under noisy features

We first investigate our so-called baseline setup which considers the four
features $m_{J_{1}}$, $\Delta m_{J}$, $\tau_{21}^{1}$ and $\tau_{21}^{2}$ and
is close to the one originally considered in [9]. We confirm that BDTs perform
at least on par with NNs in this setup as shown in Fig. 1. The IAD BDT
outperforms the IAD NN for all values of signal efficiency and achieves a
maximum SIC value of 15. As expected, supervised training using signal and
background labels outperforms the weakly supervised IAD, with a difference in
maximum achievable SIC of about 2. For the relatively simple supervised
classification task, BDT and NN achieve the same performance.

Figure 1: Significance improvement characteristic (SIC) curve, eq. (2), for
the baseline scenario with four features introduced in [9]. We compare the
idealized anomaly detector (IAD) and a fully supervised setting for the BDT
and the NN classifier.

Figure 2: The impact of uninformative features on the NN and BDT classifiers.
We show the SIC curves of the IAD NN/BDT classifiers employing the four
baseline features and additional 1, 2, 5, 10, 30 and 50 input features drawn
from Gaussian noise (left/right panel). For 30 and 50 Gaussian noise features,
BDT classifiers with $N=100$ are used instead of the usual $N=50$ ensembling.

Next, we test the stability of NN- and BDT-based idealized anomaly detection
in the presence of noisy features. Since in anomaly detection the potential
signal is not known, it is expected that some input features do not carry
relevant information for its discrimination. We mimic this by adding
additional input features that consist of Gaussian noise with mean zero and a
standard deviation of one. The rest of the architecture is unchanged and the
noise features are drawn from the same normal distribution for signal and
background. In Fig. 2 (left panel) we observe a drastic drop in SIC to a
maximum value below 4 when only five Gaussian noise (5G) features are included
in the idealized neural network classifier. The performance of the BDT IAD, on
the other hand, is much more stable against uninformative features. This is
illustrated in Fig. 2 (right panel), where we show that the BDT IAD still
performs significantly better than random for up to 50 Gaussian noise
features. The robustness of the BDT with respect to uninformative features is
analyzed in more detail in Appendix B.

Ensembling, as introduced in Section II.2, is important for the performance
and stability of the BDT classifier in the presence of noisy features. For the
NN, on the other hand, the improvement due to ensembling is much less
significant (and much more expensive to obtain). As the size of the error band
suggests, the BDT classifier would even benefit from a larger ensemble with
$N>100$ when many noisy features are added. The effect of ensembling is
discussed in more detail in Appendix C.

###  IV.2 Expanding the pool of features

Having identified the BDTs as a method that is stable under the inclusion of
potentially noisy features, we can explore how additional physics features can
improve anomaly detection performance. In Fig. 3, we compare the SIC curves
for the IAD in the baseline setup as in [9] with the IAD employing various
extended feature sets using additional jet substructure information. The
extended sets are given in Table 1, and include additional subjettiness ratios
(extended set 1), individual subjettiness features up to $\tau_{5}$ (extended
set 2), and a set of 54 subjettiness features computed with different angular
weighting parameters $\beta$ (extended set 3). We show results for the NN IAD
classifier (Fig. 3, left panel) and the BDT IAD classifier (Fig. 3, right
panel).

Figure 3: The impact of physics features on the NN and BDT classifiers. We
show the SIC curves of the IAD NN with the four baseline features and the
extended feature sets based on subjettiness variables as introduced in Tab. 1
(left panel), and the corresponding results for the IAD BDT (right panel).

As expected from the study of Gaussian noise, Section IV.1, the NN IAD
classifier is very sensitive to the selection of features. The inclusion of
additional subjettiness ratios leads to a dramatic reduction in performance,
while the inclusion of individual subjettiness features improves the
classification properties.

The BDT IAD classifier, on the other hand, is characterized not only by its
higher performance, but also by its stability under the inclusion of distinct
features. Additional subjettiness ratios and, in particular, the inclusion of
individual subjettiness features lead to a dramatic improvement in
performance. Thus, unlike the NN, adding more and more features does not
affect the performance of the BDT.

From our studies in Section IV.1, it is clear that adding more and more noisy
features will eventually degrade performance. We observed a slight drop in
performance for the BDT when we added even more features (beyond extended set
3) based on the subjettinesses. However, again the NN classifier suffers much
more in this case.

###  IV.3 Dependence on amounts of signal and background

So far we have only considered a fixed number of signal and background events
in the SR. To detect a potential signal, it is of course not essential whether
our classifier reaches ${\rm SIC}=20$ or ${\rm SIC}=40$. Much more important
is the minimum signal fraction that is still detected by the classifier for a
given signal region size. Therefore, it is important to investigate the
potential sensitivity of the classifiers for a range of signal fractions and
signal region sizes.

Figure 4 shows the maximum SIC value as a function of the number of signal
events in the SR for the default fixed number of 120k background events.

Figure 4: The impact of physics features on the sensitivity improvement. We
show the maximum SIC values as a function of the number of signal events in
the SR, with the number of background events held fixed. The BDT IAD
performance is shown for the baseline feature set and the extended feature set
3. The NN performance is shown for the baseline set for comparison.

Achieving significance improvement for rarer signals will translate into
analysis sensitivity for lower cross-sections. For the baseline set of
features, the BDT plateaus at a SIC of 15, which is reached with a minimum of
approximately 500 signal events, corresponding to $S/B=4\times 10^{-3}$. The
same SIC value can be reached with the extended feature set 3 with less than
half the number of events, reaching a plateau SIC value of about 40.

Figure 5: The maximum SIC values as a function of the number of signal events
$N_{sig}$ in the SR and the number of background events $N_{bkg}$ in the SR
and the background template. The BDT IAD performance is shown for the baseline
feature set (upper panel) and the more inclusive extended set 3 (lower panel).

Figure 5 shows the maximum improvement in significance for the IAD BDT as a
function of the number of signal and background events in the signal region.
In contrast to the previous results, for which the background template is
larger than the signal region (see Section III.1) in accordance with the
studies in Ref. [9], in Fig. 5 the size of the background template is equal to
the number of background events $N_{bkg}$ in the SR such that a larger range
of $N_{bkg}$ may be scanned with the existing data. The top panel shows the
results for the baseline setup, and the bottom panel shows the results for the
inclusive set of all subjettiness features (extended set 3). In both panels it
can be seen that the maximum significance improvement for a given signal
significance $S/\sqrt{B}$ is largely independent of the absolute number of
background events. Hence, for larger datasets, the same significance
improvement is achieved with a smaller absolute signal fraction. Comparing the
top and bottom panels, the dramatic increase in performance of the BDT
classifier is evident, with SIC values increasing from 15 for the small
feature set to about 40 for the extended feature set. Even more importantly,
the minimal significance in the original dataset for a discovery with high
significance is decreased quite a bit.

###  IV.4 Three-prong signal

Figure 6: Same as Fig. 3 but for a signal model with three-prong jets.

So far, we have focused on the LHCO R&D signal consisting of 2-pronged jets.
Here we will compare the performance of NN vs. BDT based classifiers across
our different feature sets for a different signal (also from the LHCO R&D
dataset) consisting of 3-pronged jets. Clearly, an anomaly detector that
claims to be both high performing and model agnostic should be able to find
both signals with high sensitivity.

Figure 6 summarizes our results on the 3-prong signal. We see that with the
baseline feature set, the NN is essentially no better than random on the
3-prong signal. The $\tau_{21}$ features of the baseline set offer very little
discrimination power for the 3-prong signal (i.e. are essentially
uninformative noise features), so they are evidently degrading the performance
of the NN. Meanwhile, the BDT is able to achieve a modest significance
improvement, presumably by leveraging the mass features and not being confused
by the uninformative $\tau_{21}$ features.

The comparison between NN and BDT is even more dramatic with extended set 1.
The NN classifier again shows its noise sensitivity: while extended set 1
includes the separating $\tau_{32}$ subjettiness ratios for each jet, it also
includes higher subjettiness ratios. These are essentially uninformative noise
features for the 3-prong signal, and thereby prevent better classification
with the NN. However, once again, the BDT does not suffer from this problem,
and achieves a significant performance gain relative to the baseline features.

Using the subjettinesses in extended sets 2 and 3 instead of the ratios yields
significantly better NN and BDT classifiers, with the BDT again achieving the
best overall performance. Comparing with Fig. 3, we see that NN and BDT
classifiers trained on extended sets 2 and 3 perform well as anomaly detectors
on both 2-prong and 3-prong signals and can reasonably be called “model
agnostic”.

##  V Conclusion

In this paper we have investigated the role of different machine learning
methods (deep learning vs. boosted decision trees) on classification
performance in weakly supervised anomaly detection. Normally, with the large
amount of data and high-dimensional, structured feature sets we have in high
energy physics, deep learning wins over BDTs hands down, and this has been
widely presumed in the literature since the deep learning revolution came to
our field. However, weakly supervised anomaly detection presents special
challenges, for which, surprisingly, deep learning may not be ideally suited.
Although the amount of data is nominally large, in this anomaly detection
context one is trying to detect minute differences between data and
background—the amount of signal is small, and so in a sense the effective
dataset size is small. Furthermore, in the anomaly detection problem, one
wants to be model agnostic, so one wants a broad selection of features, and
for a given signal model maybe only a small subset of them are actually
informative (able to discriminate between signal and background). The rest of
the features are uninformative noise. It is well known and well studied in the
ML literature that BDTs can outperform deep learning in these settings—tabular
data, small to medium dataset sizes, and in the presence of noisy,
uninformative features, see e.g. [14].

In this work, we have confirmed this in the context of weakly supervised
anomaly detection, applied to the setting of the LHC. We have seen that adding
by hand a number of uninformative Gaussian noise features to both signal and
background completely destroys the performance of an NN classifier, but the
anomaly detection ability of a BDT classifier remains robust for at least up
to $\sim 50$ Gaussian noise features. We then considered the performance of
the idealized anomaly detector with the addition of more physics-motivated
features (additional n-subjettinesses $\tau_{n=3,\dots,9}$, although other
families of features such as e.g. energy flow polynomials [33] would be
possible as well). The NN is unstable with respect to additional features and
its performance can deteriorate or improve depending on which features are
added. Meanwhile the BDT’s performance is more stable and consistently
improves as more features are added (up to a certain point—see Fig. 2).
Finally, we have shown that with more n-subjettiness features, both NNs and
BDTs can in principle detect both 2-prong and 3-prong signals, better
realizing the original promise of model-agnostic idealized anomaly detection.

One issue we encountered is the lack of reliable schemes for choosing anomaly
detection methods. We are interested in such model-agnostic metrics for model
selection; this is work in progress. Also, the issue of $m_{JJ}$ sculpting was
left for future work—this important issues should be studied in more detail
and something like LaCathode [11] may be needed to mitigate the sculpting
problems.

Furthermore, this work was focused on the idealized situation where the
background templates are assumed to be perfectly modeled. This allowed us to
fairly and evenly compare the performance of BDTs and NNs across different
feature sets. However, in any realistic scenario, the background template must
be derived from data (e.g. fully data-driven side-band interpolation as in [6,
9, 10, 11] or simulation-assisted interpolation as in [7, 8, 12]). An
essential next step in this study of weakly supervised, model-agnostic anomaly
detection is to consider the performance of BDT and NN classifiers using
various feature sets that are realistically derived from data.

The prospects for this study are very bright. These techniques (BDT
classifiers) can be easily applied to real LHC data, searches based on CWoLa
hunting by ATLAS already exist [34], and their physical performance should be
significantly improved using BDTs and larger feature sets. Although it was not
the primary objective of this work, we also observe that BDTs are much more
efficient in their use of computational resources than NNs—they can be trained
and evaluated in a fraction of the time. This could be enormously beneficial
to real searches (such as [34]) that use weak supervision, where a very large
number of classifiers typically need to be trained and evaluated, especially
for the estimation of systematic uncertainties.

Note added: While this work was being completed, we learned of the work of
[35] who also studied the applications of BDTs to weakly supervised anomaly
detection at the LHC.

## Acknowledgements

TF, MH, MK, and AM would like to thank Thea Klaeboe Aarrestad for sparking our
interest in boosted decision trees. TF is supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant 400140256
- GRK 2497: The physics of the heaviest particles at the Large Hadron
Collider. The research of MH, MK and AM is supported by the DFG under grant
396021762 - TRR 257: Particle Physics Phenomenology after the Higgs Discovery.
GK, MS, and TQ acknowledge support by the DFG under Germany’s Excellence
Strategy 390833306 – EXC 2121: Quantum Universe. DS is supported by DOE grant
DOE-SC0010008. Computations were performed with computing resources granted by
RWTH Aachen University under project rwth0934.

## Code

The code for this paper can be found at https://github.com/uhh-pd-
ml/treebased_anomaly_ detection.

##  Appendix A Comparison of different BDT architectures

Here, we investigate the performance of other tree-based algorithms in
addition to the HGBDT used in the main text. The considered models are:

  * •

HGBDT (Histogrammed Gradient BDT): A histogrammed version of gradient boosted
decision trees. Gradient boosting trains subsequent models on the residuals of
the previous ensemble state. The HistGradientBoostingClassifier implementation
in scikit-learn [20] is used.

  * •

Random forest: Multiple independent decision trees are built based on
subsampling and use a majority vote for classification. The
RandomForestClassifier implementation in scikit-learn [20] is used.

  * •

Adaboost: A boosting algorithm where decision trees are iteratively added to
improve classification results. Weights are introduced for the trees as well
as for the samples to focus on misclassified samples. The AdaBoostClassifier
implementation in scikit-learn [20] is used.

  * •

ROOT TMVA BDT: A BDT implementation that is frequently used in Particle
Physics applications. Here, we us version 6.28.4 [36]. Using the default
settings, this implementation is also based on Adaboost.

For the weakly supervised classification both with and without additional
noise features, the results are summarized in Fig. 7. As in the main text ten
classifiers are used for creating the SIC curve and the corresponding error
bands. However, instead of using ensembling with $N=50$, we here use $N=10$
because not all models are fast to train. For the ROOT TMVA model, we used the
default settings. The optimal parameters of the other models were tuned based
on the validation loss of the weak classification task without noise using the
optuna software package [37]. For the GBDT, we find that the default
parameters described in Section II.2 are close to optimal such that we do not
use any optimization in the main text.

Figure 7: Comparison of different tree-based algorithms for the baseline setup
(left panel) as well as for the baseline setup with additional 10 Gaussian
noise features (right panel). We show results for the HGBDT used in the main
text, a random forest (RF), an Adaboost classifier as well as the BDT
classifier implementation in ROOT TMVA [36].

Without noise, all tree-based algorithms show similar performance, in
particular for low signal efficiencies. When adding 10 Gaussian noise
features, however, a clear performance gap can be seen: While the Adaboost,
random forest and ROOT TMVA models drop significantly, the histogrammed GBDT
retains much of its original performance. Due to the superiority of the
histogrammed GBDT in these noise studies, we use it throughout this work.
These results are compatible with what is generally considered the state of
the art in DT based models [19].

##  Appendix B Uninformative features, rotational invariance, and tabular data

Figure 8: The impact of a random rotation on the BDT IAD and NN IAD using the
baseline feature set and the baseline features with additional 3 Gaussian
noise features. SIC curves are compared with and without a random rotation
applied for NN models (left panel) and BDT models (right panel).

In [14], a key argument for decision tree (DT) based models outperforming deep
learning on tabular data is their robustness to uninformative features. This
robustness stems from the fact that tabular data is not invariant under
rotation, and therefore an algorithm with the same property should be used to
learn from it. Since NNs are rotationally invariant [38], they must learn the
ideal feature orientation in an increasingly high-dimensional input space and
then identify the most informative features. The information in the original
orientation of the data, which in our case is physics based, is lost. On the
other hand, DT-based models are not rotationally invariant and therefore only
operate on the correct orientation with respect to the physics of the event.

Since the data used in this work is tabular, we expect the same behavior in
the weakly supervised setting: While the performance of BDT-based models
should break down significantly under a random rotation of the data, NN-based
models should still show a similar performance. However, if no rotation is
applied, BDT-based models should outperform NNs, and the performance gap is
expected to increase significantly as more uninformative features are used.

Figure 8 shows the performance of both BDT and NN with and without a random
rotation applied to the input features. As in Appendix A, the medians and
error bands are based on ten classifiers using ensembling with $N=10$. We
observe the expected drop in performance for a BDT trained on rotated
features, where the maximum SIC value on the baseline features is halved. When
three noise features are added, the BDT is unable to maintain better than
random performance. A slight drop in SIC of about 2 is also observed for the
NN on the rotated features. With the addition of noise features, the drop is
slightly larger, but the NN still remains better than random.

##  Appendix C Ensembling

Figure 9: The impact of ensembling on the baseline setup with four features
and additional 10 Gaussian noise features for both NN (left panel) and BDT
(right panel). The gray lines show individual SICs of single BDTs/NNs.
Additionally, the median performance of the individual BDTs/NNs and the
performance of the ensemble classifiers are shown. Note that for the BDT the
ensemble classifier is one of the classifiers used to calculate the median
shown in Fig. 2.

In general, the use of ensembles of different classifiers is a well-
established technique in many machine learning applications. A single BDT is
already a form of ensemble learning, as boosting combines different decision
trees. Often, different ML models are combined in an ensemble to exploit their
different strengths. Typically, the ensemble outperforms each of the
individual models. This can be seen for example in the comparison of the top
taggers in [39].

In weakly supervised anomaly detection, the goal is to lower the threshold for
detecting a small fraction of signal in an overwhelming background. Around
this threshold, the performance of classifiers becomes unstable. The
statistical fluctuation in the assignment of signal events to the training or
validation set can be large enough to have a significant impact on the quality
of the classification. In addition, subsampling is used in our BDT
implementation, resulting in variability in performance. For NNs, the random
initialization of the network has a similar effect. The large variance in
individual BDT or NN predictions close to the signal fraction threshold
requires the use of an ensemble classifier as described in Sections II.2 and
II.3.

For the BDT classifier with an ensembling of $N=50$, the performance is
significantly increased even with respect to the best individual training, as
shown in Fig. 9 for the default baseline setup with 10 additional Gaussian
noise features. The variance of the ensemble classifier is quite small as
shown in Fig. 2. For the NN, an increase in the performance of the ensemble
classifier is also observed. However, since no single run achieves a max SIC
value significantly above 2, the NN ensemble cannot match the performance of
the BDT ensemble.

In addition, considering the very different computational costs of training
BDTs and NNs, the BDT ensemble not only performs better, but is also much more
cost effective.

## References

  * [1] G. Kasieczka _et al._ , arXiv:2101.08320 [hep-ph] .
  * [2] T. Aarrestad _et al._ , arXiv:2105.14027 [hep-ph] .
  * Metodiev  _et al._ [2017] E. M. Metodiev, B. Nachman, and J. Thaler, JHEP 10, 174 (2017), arXiv:1708.02949 [hep-ph] .
  * Collins  _et al._ [2018] J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. Lett. 121, 241803 (2018), arXiv:1805.02664 [hep-ph] .
  * Collins  _et al._ [2019] J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. D99, 014038 (2019), arXiv:1902.02634 [hep-ph] .
  * Nachman and Shih [2020] B. Nachman and D. Shih, Phys. Rev. D 101, 075042 (2020), arXiv:2001.04990 [hep-ph] .
  * Andreassen  _et al._ [2020] A. Andreassen, B. Nachman, and D. Shih, Phys. Rev. D 101, 095004 (2020), arXiv:2001.05001 [hep-ph] .
  * Benkendorfer  _et al._ [2021] K. Benkendorfer, L. L. Pottier, and B. Nachman, Phys. Rev. D 104, 035003 (2021), arXiv:2009.02205 [hep-ph] .
  * Hallin  _et al._ [2022] A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel, M. Schlaffer, D. Shih, and M. Sommerhalder, Phys. Rev. D 106, 055006 (2022), arXiv:2109.00546 [hep-ph] .
  * Raine  _et al._ [2023] J. A. Raine, S. Klein, D. Sengupta, and T. Golling, Front. Big Data 6, 899345 (2023), arXiv:2203.09470 [hep-ph] .
  * Hallin  _et al._ [2023] A. Hallin, G. Kasieczka, T. Quadfasel, D. Shih, and M. Sommerhalder, Phys. Rev. D 107, 114012 (2023), arXiv:2210.14924 [hep-ph] .
  * Golling  _et al._ [2023] T. Golling, S. Klein, R. Mastandrea, and B. Nachman, Phys. Rev. D 107, 096025 (2023), arXiv:2212.11285 [hep-ph] .
  * [13] T. Golling, G. Kasieczka, C. Krause, R. Mastandrea, B. Nachman, J. A. Raine, D. Sengupta, D. Shih, and M. Sommerhalder, arXiv:2307.11157 [hep-ph] .
  * [14] L. Grinsztajn, E. Oyallon, and G. Varoquaux, arXiv:2207.08815 [cs.LG] .
  * Borisov  _et al._ [2022] V. Borisov, T. Leemann, K. Seßler, J. Haug, M. Pawelczyk, and G. Kasneci, IEEE Transactions on Neural Networks and Learning Systems , 1 (2022).
  * Neyman and Pearson [1933] J. Neyman and E. S. Pearson, Phil. Trans. Roy. Soc. Lond. A 231, 289 (1933).
  * Finke  _et al._ [2022] T. Finke, M. Krämer, M. Lipp, and A. Mück, JHEP 08, 015 (2022), arXiv:2204.11889 [hep-ph] .
  * Ke  _et al._ [2017] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu, Advances in neural information processing systems 30, 3146 (2017).
  * [19] H. Carlens, “State of competitive machine learning in 2022,” https://mlcontests.com/state-of-competitive-data-science-2022.
  * Pedregosa  _et al._ [2011] F. Pedregosa _et al._ , Journal of Machine Learning Research 12, 2825 (2011).
  * Abadi  _et al._ [2015] M. Abadi _et al._ , “TensorFlow: Large-scale machine learning on heterogeneous systems,” (2015), software available from https://www.tensorflow.org/.
  * Chollet  _et al._ [2015] F. Chollet _et al._ , “Keras,” https://keras.io (2015).
  * [23] D. P. Kingma and J. Ba, arXiv:1412.6980 .
  * Kasieczka  _et al._ [2019] G. Kasieczka, B. Nachman, and D. Shih, “R&d dataset for lhc olympics 2020 anomaly detection challenge,” https://zenodo.org/record/6466204 (2019).
  * Sjöstrand  _et al._ [2006] T. Sjöstrand, S. Mrenna, and P. Z. Skands, JHEP 05, 026 (2006), arXiv:hep-ph/0603175 [hep-ph] .
  * Sjöstrand  _et al._ [2008] T. Sjöstrand, S. Mrenna, and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008), arXiv:0710.3820 [hep-ph] .
  * de Favereau  _et al._ [2014] J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaitre, A. Mertens, and M. Selvaggi (DELPHES 3), JHEP 02, 057 (2014), arXiv:1307.6346 [hep-ex] .
  * Cacciari  _et al._ [2012] M. Cacciari, G. P. Salam, and G. Soyez, Eur. Phys. J. C72, 1896 (2012), arXiv:1111.6097 [hep-ph] .
  * Cacciari and Salam [2006] M. Cacciari and G. P. Salam, Phys. Lett. B641, 57 (2006), arXiv:hep-ph/0512210 [hep-ph] .
  * Shih [2021] D. Shih, “Additional qcd background events for lhco2020 r&d (signal region only),” https://zenodo.org/record/5759086 (2021).
  * Thaler and Van Tilburg [2011] J. Thaler and K. Van Tilburg, JHEP 03, 015 (2011), arXiv:1011.2268 [hep-ph] .
  * Thaler and Van Tilburg [2012] J. Thaler and K. Van Tilburg, JHEP 02, 093 (2012), arXiv:1108.2701 [hep-ph] .
  * Komiske  _et al._ [2019] P. T. Komiske, E. M. Metodiev, and J. Thaler, JHEP 01, 121 (2019), arXiv:1810.05165 [hep-ph] .
  * Aad  _et al._ [2020] G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 125, 131801 (2020), arXiv:2005.02983 [hep-ex] .
  * Freytsis  _et al._ [2023] M. Freytsis, M. Perelstein, and Y. C. San, (2023), arXiv:2310.xxxxx [hep-ph] .
  * Brun and Rademakers [1997] R. Brun and F. Rademakers, Nucl. Inst. & Meth. in Phys. Res. A 389, 81 (1997), proceedings AIHENP’96 Workshop, Lausanne, Sep. 1996.
  * [37] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, “Optuna: A next-generation hyperparameter optimization framework,” arXiv:1907.10902 .
  * Ng [2004] A. Y. Ng, in  _Proceedings of the Twenty-First International Conference on Machine Learning_, ICML ’04 (Association for Computing Machinery, New York, NY, USA, 2004) p. 78.
  * Butter  _et al._ [2019] A. Butter _et al._ , SciPost Phys. 7, 014 (2019), arXiv:1902.09914 [hep-ph] .


# Synchronizability in randomized weighted simplicial complexes

S. Nirmala Jenifer  Department of Physics, Bharathidasan University,
Tiruchirappalli 620024, Tamil Nadu, India  Dibakar Ghosh<EMAIL_ADDRESS>Physics and Applied Mathematics Unit, Indian Statistical Institute, 203 B. T.
Road, Kolkata - 700108, India  Paulsamy Muruganandam<EMAIL_ADDRESS>Department of Physics, Bharathidasan University, Tiruchirappalli 620024, Tamil
Nadu, India

###### Abstract

We present a formula for determining synchronizability in large, randomized
and weighted simplicial complexes. This formula leverages eigenratios and
costs to assess complete synchronizability under diverse network topologies
and intensity distributions. We systematically vary coupling strengths
(pairwise and three-body), degree and intensity distributions to identify the
synchronizability of these simplicial complexes of the identical oscillators
with natural coupling. We focus on randomized weighted connections with
diffusive couplings and check synchronizability for different cases. For all
these scenarios, eigenratios and costs reliably gauge synchronizability,
eliminating the need for explicit connectivity matrices and eigenvalue
calculations. This efficient approach offers a general formula for
manipulating synchronizability in diffusively coupled identical systems with
higher-order interactions simply by manipulating degrees, weights, and
coupling strengths. We validate our findings with simplicial complexes of
Rössler oscillators and confirm that the results are independent of the number
of oscillators, connectivity components and distributions of degrees and
intensities. Finally, we validate the theory by considering a real-world
connection topology using chaotic Rössler oscillators.

##  I Introduction

Complex systems arise when a large number of dynamical units interact with
each other, giving rise to emergent properties distinct from those of the
individual subsystems [1]. Such systems are ubiquitous, found naturally in
entities like the brain and man-made ones like the internet or financial
markets [2]. The study and modelling of these systems have been a longstanding
focus of researchers [3, 4, 5]. Complex networks represent the dynamical
units, and their interactions as nodes and links are a common way to model
such systems. However, these networks typically only account for pairwise
interactions, whereas many complex systems involve interactions among three or
more units. To gain a more comprehensive understanding of these systems, one
must incorporate higher-order networks, such as hypergraphs and simplicial
complexes, to capture these higher-order interactions [6].

Simplicial complexes are used to model complex systems with higher-order
interactions [4, 7, 8, 9]. A simplicial complex consists of $d$-simplices,
where $d$ is the dimension of the simplex. A $0$-dimensional simplex is a
node, a $1$-dimensional simplex is a link, a $2$-dimensional simplex is a
triangle, and so on. A $d$-dimensional simplicial complex consists of
simplices up to dimension $d$. They provide a good representation of complex
systems with higher-order interactions [10, 11, 7, 12]. In most studies, these
simplicial complexes are unweighted, meaning all the links and triangles have
the same weight [13, 14]. However, it has been observed that in most real-
world systems, this is not feasible. For instance, in collaboration networks,
if there is a three-author paper, not every pair of authors has to produce
research papers by themselves. But in unweighted simplicial complexes, if
there is a three-author paper, then there is also the presence of all the
possible two-author research papers, and this leads to information loss, which
one could avoid by introducing proper weights for links [15]. Courtney and
Bianconi [16] have developed a model to distribute weights to links in terms
of bare affinity weights and topological weights. Every simplex has bare-
affinity weights, while the topological weights represent whether a particular
link has a contribution other than being part of a triangle. Therefore,
weighted simplicial complexes provide a more realistic representation of
higher-order networks.

After accurately modelling complex systems, we can better understand their
emergent properties, such as synchronization. As a result, we can control the
synchronization that occurs in natural or artificial systems. Synchronization
occurs when individual dynamical systems adjust their properties to have
common dynamics [3]. Neuronal synchronization can cause epilepsy, while
synchronization from a healthy to an infected state can cause an epidemic
outbreak, and so on. Synchronization can be both constructive and destructive.
However, it is possible to promote synchronization in desired systems and
inhibit it where it leads to destruction. This can be achieved by determining
under what circumstances the system goes into synchronization and
desynchronization and whether synchronization is possible. Two main parameters
are used to calculate the synchronizability of the system: eigenratio and the
cost of the connections of nodes in the system. The lower the eigenratio and
cost, the system has more synchronizability [17].

To calculate the eigenratio and costs, we need to construct the connectivity
matrices, such as adjacency matrices and Laplacian matrices, and then find the
eigenvalues of these matrices. Writing connectivity matrices for large,
complex systems, such as connection networks, collaboration networks, and
neuronal networks, can be challenging, and collecting the necessary data is
not an easy task. It becomes even more difficult when we include higher-order
interactions, as we need to construct the adjacency tensors for these
interactions. Even for three-body interactions, the adjacency tensor contains
$N\times N\times N$ elements. In 2006, Zhou et al. [17] developed an approach
to calculate eigenratios and costs based on the heterogeneity of the system
intensities. This approach significantly reduced the computational cost of
constructing connectivity matrices and evaluating their eigenvalues. We extend
this approach to randomized weighted simplicial complexes. As a result, we can
determine the synchronizability parameters solely from the heterogeneity of
the intensities, specifically the maximum and minimum intensities of a node as
part of links ($d=1$), triangles ($d=2$), and so on and from their coupling
strengths. In this paper, we discuss the mathematical formulation, derive
general formulas for determining the synchronization without explicitly
computing the eigenvalues of huge matrices and corroborate with numerical
results. We validate the results by applying the general formula to Rössler
oscillators and real-world connection networks. We check the emergence of
complete synchronization in diffusively coupled identical Rössler oscillators,
and we show that the results are independent of the number of oscillators in
the networks, the components through which they are connected and also the
distributions of the degrees and intensities.

##  II Analytical Results

In order to derive the general formulas for the eigenratio and the cost, we
first write the dynamical equations for the weighted simplicial complex. Then,
we deduce the corresponding variational equations and modify these equations
in terms of the effective matrix $M$. The eigenvalues of the effective matrix
determine the stability of the synchronized state. Since the effective matrix
is a zero row sum matrix, the first eigenvalue will be zero, corresponding to
the mode along the synchronization manifold. We can find the synchronizability
of the system from the ratio between $\lambda_{N}$ and $\lambda_{2}$. For
simplicity, we consider the simplicial complex of dimension $2$, which we can
extend to any dimension.

We consider a simplicial complex of identically coupled oscillators, described
by the following equations,

$\displaystyle\mathbf{\dot{x}}_{i}=$
$\displaystyle\mathbf{f}(\mathbf{x}_{i})+\sigma_{1}\sum_{j=1}^{N}a_{ij}^{(1)}\omega_{ij}^{(1)}\mathbf{g}^{(1)}(\mathbf{x}_{i},\mathbf{x}_{j})$
$\displaystyle+\sigma_{2}{\sum_{j=1}^{N}}{\sum_{k=1}^{N}}a_{ijk}^{(2)}\omega_{ijk}^{(2)}\mathbf{g}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k}),$
(1)

where $\mathbf{x}_{i}$ is the state vector of $i$-th oscillator of dimension
$m$, and $\mathbf{f}(\mathbf{x}_{i})$ represents the dynamics of the uncoupled
oscillators. ${\omega_{ij}}^{(1)}$ and $\omega_{ijk}^{(2)}$ are the
topological weights of links and triangles, respectively. $a_{ij}^{(1)}$ and
$a_{ijk}^{(2)}$ are the elements of the adjacency matrix $A^{(1)}$ and
adjacency tensor $A^{(2)}$. $a_{ij}^{(1)}=1$ if the nodes $i$ and $j$ form a
link, or $a_{ij}^{(1)}=0$ otherwise. Also $a_{ijk}^{(2)}=1$ if the nodes $i$,
$j$, and $k$ form a triangle, or $a_{ijk}^{(2)}=0$ otherwise. Here
$\sigma_{1}$ and $\sigma_{2}$ are the coupling strengths of the links
(pairwise) and triangles (non-pairwise), respectively. Here,
$\mathbf{g}^{(1)}(\mathbf{x}_{i},\mathbf{x}_{j})$ and
$\mathbf{g}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k})$ are the
synchronization noninvasive functions [7] and are chosen in the forms
$\mathbf{g}^{(1)}(\mathbf{x}_{i},\mathbf{x}_{j})=\mathbf{h}^{(1)}(\mathbf{x}_{j})-\mathbf{h}^{(1)}(\mathbf{x}_{i})$
and
$\mathbf{g}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k})=\mathbf{h}^{(2)}(\mathbf{x}_{j},\mathbf{x}_{k})-\mathbf{h}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{i})$.
At the synchronization state, they tend to zero, making the dynamics of the
system resemble that of the uncoupled oscillators.

For the above diffusive coupling functions, we can rewrite the above equation
as

$\displaystyle\mathbf{\dot{x}}_{i}=$
$\displaystyle\,\mathbf{f}(\mathbf{x}_{i})+\sigma_{1}\sum_{j=1}^{N}a_{ij}^{(1)}\omega_{ij}^{(1)}\left[\mathbf{h}^{(1)}(\mathbf{x}_{j})-\mathbf{h}^{(1)}(\mathbf{x}_{i})\right]$
$\displaystyle\,+\sigma_{2}{\sum_{j=1}^{N}}{\sum_{k=1}^{N}}a_{ijk}^{(2)}\omega_{ijk}^{(2)}\left[\mathbf{h}^{(2)}(\mathbf{x}_{j},\mathbf{x}_{k})-\mathbf{h}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{i})\right],$
(2)

where $\mathbf{h}^{(1)}(\mathbf{x}_{j})$ and
$\mathbf{h}^{(2)}(\mathbf{x}_{j},\mathbf{x}_{k})$ are coupling functions that
couples nodes in links and triangles, respectively. The choice of coupling
functions affects the synchronizability of the system [7]. So we have to
choose the coupling functions in a way that the oscillators tend to
synchronize and we can write,

$\displaystyle
S_{i}^{(1)}=\sum_{j=1}^{N}{a_{ij}}^{(1)}\omega_{ij}^{(1)}\,\;\;\mbox{and}\;\;S_{i}^{(2)}=\frac{1}{2}\sum_{j=1}^{N}\sum_{k=1}^{N}a_{ijk}^{(2)}\omega_{ijk}^{(2)}.$

Here $S_{i}^{(1)}$ and $S_{i}^{(2)}$ are the intensities of node $i$ to form
links and triangles, i.e., the number of weighted links and weighted triangles
incident on node $i$. For randomized simplicial complexes with
$K^{(1)}_{\mbox{min}}\gg 1$, we can use mean-field approximation and this
equation can be rewritten as,

$\displaystyle\dot{\mathbf{x}}_{i}=$
$\displaystyle\,\mathbf{f}(\mathbf{x}_{i})+\sigma_{1}\frac{S_{i}^{(1)}}{K_{i}^{(1)}}\sum_{j=1}^{N}a_{ij}^{(1)}\left(\mathbf{h}^{(1)}(\mathbf{x}_{j})-\mathbf{h}^{(1)}(\mathbf{x}_{i})\right)$
$\displaystyle+\sigma_{2}\frac{S_{i}^{(2)}}{K_{i}^{(2)}}{\sum_{j=1}^{N}}{\sum_{k=1}^{N}}a_{ijk}^{(2)}\left(\mathbf{h}^{(2)}(\mathbf{x}_{j},\mathbf{x}_{k})-\mathbf{h}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{i})\right).$
(3)

Here, $K_{i}^{(1)}$ is the degree of pairwise interactions, i.e., the total
number of links that a node participates in. $K_{i}^{(2)}$ is the total number
of triangles that a node participates in, also called the degree for three-
body interactions. Now, we consider pairwise and non-pairwise local mean
fields as

$\displaystyle\mathbf{\overline{H}}^{(1)}{(\mathbf{x}_{i})}=\frac{1}{K_{i}^{(1)}}\sum_{j=1}^{N}a_{ij}^{(1)}\mathbf{h}^{(1)}(\mathbf{x}_{j})$

due to the interaction between the nodes through links and

$\displaystyle\mathbf{\overline{H}}^{(2)}{(\mathbf{x}_{i},\mathbf{x}_{i}})=\frac{1}{K_{i}^{(2)}}{\sum_{j=1}^{N}}{\sum_{k=1}^{N}}a_{ijk}^{(2)}\mathbf{h}^{(2)}(\mathbf{x}_{j},\mathbf{x}_{k})$

due to the interaction between nodes through triangles. Then, Eq. (3) can be
written in terms of mean field as,

$\displaystyle\mathbf{\dot{x}}_{i}=$
$\displaystyle\,\mathbf{f}(\mathbf{x}_{i})+\sigma_{1}S_{i}^{(1)}\left(\mathbf{\overline{H}}^{(1)}(\mathbf{x_{i}})-\mathbf{h^{(1)}}(\mathbf{x}_{i})\right)$
$\displaystyle+\sigma_{2}{S_{i}^{(2)}}\left(\mathbf{\overline{H}}^{(2)}{(\mathbf{x}_{i},\mathbf{x}_{i}})-\mathbf{h}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{i})\right).$
(4)

So, the interaction between nodes in a large randomized weighted simplicial
complex can be approximated as the interaction of a node with a mean field.
This mean field is not only generated by the interaction of links, as in the
case of complex networks, but also by the interaction of nodes as triangles.

According to the condition for natural coupling, close to the synchronized
state, the interaction between nodes in links and triangles will be similar
since, upon synchronization, the two nodes that are in the same state can be
considered as a single node, and the three-body interactions reduce to the
two-body interactions [7]. The above assumptions enable us to write
$\mathbf{h}^{(2)}(\mathbf{x},\mathbf{x})=\mathbf{h}^{(1)}(\mathbf{x})$ and
$\overline{\mathbf{H}}^{(2)}{(\mathbf{x},\mathbf{x})}=\overline{\mathbf{H}}^{(1)}{(\mathbf{x})}=\mathbf{H(x)}$.
Equation (4) can be rewritten as,

$\displaystyle\mathbf{\dot{x}}_{i}=\mathbf{f}(\mathbf{x}_{i})+\left(\sigma_{1}S_{i}^{(1)}+\sigma_{2}S_{i}^{(2)}\right)\left(\mathbf{H}(\mathbf{x})-\mathbf{h}^{(1)}(\mathbf{x}_{i})\right).$
(5)

The variational equation of the above Eq. (5) is of the form

$\displaystyle\delta\mathbf{\dot{x}}_{i}=$
$\displaystyle\left[J\mathbf{f}(\mathbf{x}^{s})-\left(\sigma_{1}S_{i}^{(1)}+\sigma_{2}S_{i}^{(2)}\right)J\left(\mathbf{h}^{(1)}(\mathbf{x}^{s})\right)\right]\delta\mathbf{x}_{i},$
(6)

where $\mathbf{x}^{s}$ is the synchronization state. The eigenvalues are
approximately equal to $\sigma_{1}S_{i}^{(1)}+\sigma_{2}S_{i}^{(2)}$,
$i=1,2,\ldots,N$. From the above equation, we can write the eigenratio [17]
as,

$\displaystyle
R\approx\frac{\sigma_{1}S_{\mbox{max}}^{(1)}+\sigma_{2}S_{\mbox{max}}^{(2)}}{\sigma_{1}S_{\mbox{min}}^{(1)}+\sigma_{2}S_{\mbox{min}}^{(2)}},$
(7)

where
$\sigma_{1}S_{\mbox{max}}^{(1)}+\sigma_{2}S_{\mbox{max}}^{(2)}\geq\mbox{max}\\{\sigma_{1}S_{1}^{(1)}+\sigma_{2}S_{1}^{(2)},\sigma_{1}S_{2}^{(1)}+\sigma_{2}S_{2}^{(2)},\cdots,\sigma_{1}S_{N}^{(1)}+\sigma_{2}S_{N}^{(2)}\\}$
will give the maximum bound for eigenvalues and similar for minimum
$\sigma_{1}S_{\mbox{min}}^{(1)}+\sigma_{2}S_{\mbox{min}}^{(2)}\leq\mbox{min}\\{\sigma_{1}S_{1}^{(1)}+\sigma_{2}S_{1}^{(2)},\sigma_{1}S_{2}^{(1)}+\sigma_{2}S_{2}^{(2)},\cdots,\sigma_{1}S_{N}^{(1)}+\sigma_{2}S_{N}^{(2)}\\}$.
These maximum and minimum bounds may not necessarily belong to any particular
node index $i$. It further simplifies the problem as it is sufficient to
consider only the maximum and minimum values of the intensities to calculate
the eigenratio and it eliminates the necessity to calculate all the
eigenvalues. It would be of great advantage when $N$ is very large.

For the case of pairwise interaction only ($\sigma_{2}=0.0$), the eigenratio
$R$ is very similar to the previous result [17]. For non-zero pairwise
coupling strength ($\sigma_{1}\neq 0.0,$), we can rewrite the above equation
as,

$\displaystyle
R\approx\frac{S_{\mbox{max}}^{(1)}+\frac{\sigma_{2}}{\sigma_{1}}S_{\mbox{max}}^{(2)}}{S_{\mbox{min}}^{(1)}+\frac{\sigma_{2}}{\sigma_{1}}S_{\mbox{min}}^{(2)}}.$
(8)

In the similar way, we can derive the eigenratio for the $d$-dimensional
simplical complex as,

$\displaystyle
R\approx\frac{S_{\mbox{max}}^{(1)}+\frac{\sigma_{2}}{\sigma_{1}}S_{\mbox{max}}^{(2)}+\cdots+\frac{\sigma_{d}}{\sigma_{1}}S_{\mbox{max}}^{(d)}}{S_{\mbox{min}}^{(1)}+\frac{\sigma_{2}}{\sigma_{1}}S_{\mbox{min}}^{(2)}+\cdots+\frac{\sigma_{d}}{\sigma_{1}}S_{\mbox{min}}^{(d)}}.$
(9)

Thus, the eigenratio of a randomized simplicial complex depends on the
coupling strengths (pairwise and non-pairwise), and maximum and minimum
intensities of the nodes. If we assume all the non-pairwise coupling strengths
are identical to the pairwise coupling strength (i.e.,
$\sigma_{2}=\sigma_{3}=\cdots=\sigma_{d}=\sigma_{1}$), then the eigenratio $R$
is independent on the coupling strengths and depends on the maximum and
minimum intensities as observed for pairwise weighted random networks [17].

To find the tight bounds of the above mean-field approximation, we can write
Eq. (5) as,

$\displaystyle\mathbf{\dot{x}_{i}}=$
$\displaystyle\,\mathbf{F}(\mathbf{x}_{i})-\sigma_{1}\sum_{j=1}^{N}G_{ij}^{(1)}\mathbf{H}^{(1)}(\mathbf{x}_{j})-\sigma_{2}\sum_{j=1}^{N}G_{ij}^{(2)}\mathbf{H}^{(2)}(\mathbf{x}_{j}),$
(10)

where $G_{ij}^{(1)}$ and $G_{ij}^{(2)}$ are the elements of matrices
$G^{(1)}=S^{(1)}D^{(1)-1}L^{(1)}$ and $G^{(2)}=S^{(2)}D^{(2)-1}L^{(2)}$,
respectively. Here $L^{(1)}$ and $L^{(2)}$ are generalized Laplacian matrices,
while $S$ and $D$ are diagonal matrices of strengths and generalized degrees,
respectively. It is a good point that we have separated weights from topology.
We can write normalized Laplacian matrices as
$\overline{L}^{(1)}=D^{-1}L^{(1)}$ and $\overline{L}^{(2)}=D^{-1}L^{(2)}$,
then $G^{(1)}=S^{(1)}\overline{L}^{(1)}$ and
$G^{(2)}=S^{(2)}\overline{L}^{(2)}$. Thus, the largest and smallest
eigenvalues of the matrices $G^{(1)}$ and $G^{(2)}$ are bounded by the
eigenvalues of $\overline{L}^{(1)}$ and $\overline{L}^{(2)}$. We can write the
linearized variational Eq. (10) by following the procedure outlined in [18]
and is given by

$\displaystyle\delta\dot{\mathbf{x}}_{i}=\left(J\mathbf{f}(\mathbf{x}^{s})-\sum_{j=1}^{N}\left[\sigma_{1}G_{ij}^{(1)}+\sigma_{2}G_{ij}^{(2)}\right]J\left(\mathbf{h}^{(1)}(\mathbf{x}^{s})\right)\right)\delta{\mathbf{x}}_{i}.$
(11)

We can define an effective matrix $M$ as,

$\displaystyle M_{ij}=G^{(1)}_{ij}+\frac{\sigma_{2}}{\sigma_{1}}G^{(2)}_{ij}.$
(12)

Equation (11) can be rewritten, in terms of effective matrix $M$, as

$\displaystyle\delta\mathbf{\dot{x}}_{i}=\left[J\mathbf{f}\left(\mathbf{x}^{s}\right)-\sigma_{1}\sum_{j=1}^{N}M_{ij}J\left(\mathbf{h}^{(1)}\mathbf{x}^{s}\right)\right]\delta\mathbf{x}_{i}.$
(13)

The eigenvalues of the matrix $M$ depend on the ratio of coupling strengths,
degrees and intensities, and the eigenratio of this matrix can be expressed as
$R=\lambda_{N}/\lambda_{2},$ where $\lambda_{2}$ and $\lambda_{N}$ are the 2nd
and $N$-th eigenvalues of $M$ such that
$0=\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq\cdots\leq\lambda_{N}$.

Another measure of synchronizability is the cost $C$ involved in the coupling
of nodes in a simplicial complex, and it is the total strength of connections
of all nodes. It can be written as,

$\displaystyle
C=\sigma_{1}\sum_{i=1}^{N}{S_{i}}^{(1)}+\sigma_{2}\sum_{i=1}^{N}S_{i}^{(2)},$
(14)

By using the eigenvalues of $M$ in Eq. (13), the normalized cost can be
written as,

$\displaystyle C_{0}=\frac{C}{N\alpha_{1}}=\frac{\Omega}{\lambda_{2}},$ (15)

where $\alpha_{1}=\sigma_{1}\lambda_{2}(M)$ and mean intensity

$\displaystyle\Omega=\frac{1}{N}\left(\sum_{i=1}^{N}S_{i}^{(1)}+\frac{\sigma_{2}}{\sigma_{1}}\sum_{i=1}^{N}S_{i}^{(2)}\right).$
(16)

From Eq. (6), the normalized cost will be,

$\displaystyle
C_{0}\approx\frac{\Omega}{{S_{\mbox{min}}^{(1)}+rS_{\mbox{min}}^{(2)}}}$ (17)

where $r={\sigma_{2}}/{\sigma_{1}}$. We can verify Eq. (8) and Eq. (17)
numerically by computing the eigenvalues of the effective matrix $M$, and we
can obtain an explicit bound from the eigenvalues of $\overline{L}^{(1)}$ and
$\overline{L}^{(2)}$. Using this bound, we can derive formulas for $R$ and
$C_{0}$ as a function of $S_{\mbox{max}}^{(1)}$, $S_{\mbox{min}}^{(1)}$ and
the generalized degrees.

We can demonstrate that the upper and lower bounds of the nonzero eigenvalues
of the effective matrix $M$ are given by the eigenvalues $\mu^{(1)}_{l}$ and
$\mu^{(2)}_{l}$ of the matrices $G^{(1)}$ and $G^{(2)}$, respectively, as
follows,

$\displaystyle{S_{\mbox{min}}^{(1)}\mu_{2}^{(1)}+rS_{\mbox{min}}^{(2)}\mu_{2}^{(2)}}\leq\lambda_{2}\leq{S_{\mbox{min}}^{(1)}+rS_{\mbox{min}}^{(2)}},$
(18)
$\displaystyle{S_{\mbox{max}}^{(1)}+rS_{\mbox{max}}^{(2)}}\leq\lambda_{N}\leq{S_{\mbox{max}}^{(1)}\mu_{N}^{(1)}+rS_{\mbox{max}}^{(2)}\mu_{N}^{(2)}}.$
(19)

For sufficiently random simplicial complexes, the spectra of the matrices
$G^{(1)}$ and $G^{(2)}$ tend to follow the semicircle law [17] and we can
write $\mu_{2}^{(1)}\approx 1-2/\sqrt{K^{(1)}}$, $\mu_{N}^{(1)}\approx
1+2/\sqrt{K^{(1)}}$, $\mu_{2}^{(2)}\approx 1-2/\sqrt{K^{(2)}}$ and
$\mu_{N}^{(2)}\approx 1+2/\sqrt{K^{(2)}}$, provided that
$K_{\mbox{min}}^{(d)}\gg\sqrt{K^{(d)}}$. Here $K^{(d)}$ is the mean degree of
the $d$-simplex. By introducing these in Eq. (18) and Eq. (19), we get the
following approximations for the bounds of $R$ and $C_{0}$ as

$\displaystyle\frac{S_{\mbox{max}}^{(1)}+rS_{\mbox{max}}^{(2)}}{S_{\mbox{min}}^{(1)}+rS_{\mbox{min}}^{(2)}}$
$\displaystyle\leq R$
$\displaystyle\leq\frac{1+\frac{2}{\sqrt{K^{(1)}}}}{1-\frac{2}{\sqrt{K^{(1)}}}}\left(\frac{S_{\mbox{max}}^{(1)}+r\frac{1+\frac{2}{\sqrt{K^{(2)}}}}{1+\frac{2}{\sqrt{K^{(1)}}}}S_{\mbox{max}}^{(2)}}{S_{\mbox{min}}^{(1)}+r\frac{1-\frac{2}{\sqrt{K^{(2)}}}}{1-\frac{2}{\sqrt{K^{(1)}}}}S_{\mbox{min}}^{(2)}}\right).$
(20) $\displaystyle\frac{\Omega}{S_{\mbox{min}}^{(1)}+rS_{\mbox{min}}^{(2)}}$
$\displaystyle\leq C_{0}$
$\displaystyle\leq\frac{1}{1-\frac{2}{\sqrt{K^{(1)}}}}\left(\frac{\Omega}{S_{\mbox{min}}^{(1)}+r\frac{1-\frac{2}{\sqrt{K^{(2)}}}}{1-\frac{2}{\sqrt{K^{(1)}}}}S_{\mbox{min}}^{(2)}}\right).$
(21)

For a given value of $K$, the synchronizability of randomly weighted
simplicial complexes with a large value of $K_{\mbox{min}}^{(1)}$ is expected
to follow a general formula as follows,

$\displaystyle
R=A_{R}\left(\frac{S_{\mbox{max}}^{(1)}+B_{R_{1}}S_{\mbox{max}}^{(2)}}{S_{\mbox{min}}^{(1)}+B_{R_{2}}S_{\mbox{min}}^{(2)}}\right),$
(22)

and

$\displaystyle
C_{0}=A_{C}\left(\frac{\Omega}{S_{\mbox{min}}^{(1)}+B_{R_{2}}S_{\mbox{min}}^{(2)}}\right),$
(23)

where $A_{R}=\frac{1+2/\sqrt{K^{(1)}}}{1-2/\sqrt{K^{(1)}}}$ and
$A_{C}=\frac{1}{1-2/\sqrt{K^{(1)}}}$. Here, $A_{R}\rightarrow 1$ and
$A_{C}\rightarrow 1$ in the limit $K^{(1)}\rightarrow\infty$. Also,
$B_{R_{1}}=r\frac{1+2/\sqrt{K^{(2)}}}{1+2/\sqrt{K^{(1)}}}$ and
$B_{R_{2}}=r\frac{1-2/\sqrt{K^{(2)}}}{1-2/\sqrt{K^{(1)}}}$.
$B_{R_{1}}\rightarrow r$ and $B_{R_{2}}\rightarrow r$ in the limit
$K^{(d)}\rightarrow\infty(d=1,2)$. They are well agreed with Eqs. (8) and
(17).

For $N$-dimensional simplicial complex,

$\displaystyle
R=A_{R}\left(\frac{S_{\mbox{max}}^{(1)}+\sum_{d=2}^{N}B_{R_{1}}^{(d)}S_{\mbox{max}}^{(N)}}{S_{\mbox{min}}^{(1)}+\sum_{d=2}^{N}B_{R_{2}}^{(d)}S_{\mbox{min}}^{(N)}}\right)$
(24)

and

$\displaystyle
C_{0}=A_{C}\left(\frac{\Omega}{S_{\mbox{min}}^{(1)}+\sum_{d=2}^{N}B_{C}^{(d)}S_{\mbox{min}}^{(N)}}\right),$
(25)

where

$\displaystyle
B_{R_{1}}^{(d)}=r\frac{1+2/\sqrt{K^{(d)}}}{1+2/\sqrt{K^{(1)}}},\;\;B_{R_{2}}^{(d)}=r\frac{1-2/\sqrt{K^{(d)}}}{1-2/\sqrt{K^{(1)}}}.$
(26)

##  III Numerical Results

###  III.1 Verification of the general formula

We compare the expression of $R$ as a function of eigenvalue ratio
${\lambda_{N}}/{\lambda_{2}}$ and $C_{0}$ as a function of the ratio
$\Omega/\lambda_{2}$ and found that they are almost identical when
$K^{(1)}_{\mbox{max}}\gg 1$.

To conduct our analysis, the parameter values are set as follows: $N=10^{3}$,
$\sigma_{1}=0.001$, $\sigma_{2}=0.01$. We arbitrarily fix the values of the
minimum and maximum intensities of the node for pairwise
$(S^{(1)}_{\mbox{min}},S^{(1)}_{\mbox{max}})$ and non-pairwise
$(S^{(2)}_{\mbox{min}},S^{(2)}_{\mbox{max}})$ interactions. In each iteration
of the simulation, we add $100$ to $S^{(1)}_{\mbox{max}}$ and keep all other
values constant. The simulation is run $10^{3}$ times. We first consider a
simplicial complex with a node $i$ connected to all other nodes, resulting in
$K^{(1)}_{\mbox{max}}=N-1$. The maximum number of triangles incident on a node
is $(N-1)\times(N-2)$, and the maximum number of triangles incident on a link
is $(N-2)$. Using this structure, we generate the Laplacian matrices $L^{(1)}$
and $L^{(2)}$. Finally, we uniformly distribute $S$ among the nodes. All the
triangles are considered as nodes having three-body interactions.

We calculate the eigenratio $R$ from the effective matrix $M$ using Eq. (12)
and compare it with Eq. (22). We then plot the values of $R$ against the
corresponding values of the eigenratio in terms of $S$ calculated from the
general formula, which reveals a linear relationship between $R$ and the
eigenratio as a function of $S$.

Figure 1: (a) $R$ (calculated using the eigenvalues of matrix $M$) vs $S$
(eigenratio calculated using the general formula), and (b)
$\Omega/\lambda_{2}$ vs $C_{0}$, expressed in terms of intensities, as
computed using the general formula for different coupling schemes and
distributions: (AA) all-to-all coupling, (PL) power law distribution, (PLICD)
power law distribution with intensities correlated to degrees, (Uni) uniform
distribution, and ($R$) $R$ computed using the eigenvalues of $M$.

Similarly, we calculate the normalized cost $C_{0}$ using Eq. (12) and compare
it with Eq. (23). We then plot the values of $\Omega/\lambda_{2}$ against the
corresponding values of $C_{0}$, which reveals a linear relationship between
them, which is shown dashed-purple line with diamond symbol (AA) in Fig. 1.

One may note that the most realistic networks are scale-free networks, where
the degrees are distributed according to a power law [19, 20]. This power-law
distribution means that only a small portion of nodes has very high degrees,
while a majority of nodes have relatively small degrees. Therefore, we
construct a simplicial complex with power-law distributed degrees and
intensities; however, they are not related. We achieve this by distributing
degrees and intensities to nodes using a power-law distribution. In Fig. 1, we
also show these results by the dotted-red line with stars (PL) and the green-
dashed line with empty circles (PLICD).

Next, we consider the scale-free distribution with intensities correlated to
degrees. We construct this by distributing the intensities, which are
proportional to degrees. Hence, the node with the maximum degree will receive
the maximum intensity, and the node with the minimum degree will receive the
minimum intensity. We apply this procedure to both three-body interactions and
two-body interactions. Subsequently, we plot the graphs for $R$ and $C_{0}$
corresponding to this distribution. The results demonstrate a linear
relationship between the values of the eigenratio and cost calculated from the
general formula. This linear relationship also holds for the uniform
distribution, as depicted by the dash-dotted pink line with square symbols
(Uni) in Fig. 1. From Fig. 1, we can conclude that the general formula holds
true regardless of the distribution of degrees and intensities. Further, we
observe that the values obtained from the general formula are nearly identical
to that obtained by finding the eigenvalues of the effective matrix $M$ (blue
dashed-line in Fig. 1).

As we can see from Fig. 1(b), the cost of connection is dependent on the
distribution. The cost is very low for the power law distribution of degrees
and power law distribution of intensities correlated to degrees. It is the
network topology of scale free networks that arise spontaneously in natural
and man made systems [19]. So, the cost of connection is low for the more
realistic networks. The cost of connection varies as the topology changes.

###  III.2 Effects of coupling strengths

Next, we study effect of the coupling strengths of pairwise and three-body
interactions. We fix the values of the intensities as
$S^{(1)}_{\mbox{min}}=10$, $S^{(1)}_{\mbox{max}}=1000$,
$S^{(2)}_{\mbox{min}}=10$ and $S^{(2)}_{\mbox{max}}=20$, and and change the
values of $r={\sigma_{2}}/{\sigma_{1}}$. Here, to change the value of $r$, we
fix the pairwise interaction coupling strength $\sigma_{1}=0.001$ and vary the
non-pairwise coupling strength $\sigma_{2}$. It is evident from Fig. 2(a) that
as the coupling strength for three-body interactions increases, the eigenratio
decreases and approaches smaller values and increases the synchronizability of
the system. When $r\approx 10$, the general formula gives the exact values for
the eigenratio. Then, for larger values of $r$, $R$ from both the approximated
and general formulas converge.

Figure 2: Variation of $R$ as a function of coupling strength ratio $r$
calculated from the general formula, the approximated formula and the
eigenvalues of $M$ ($R$), for the simplicial complex when degrees and
intensities are uniformly distributed: (a) and (c):
$S^{(1)}_{\mbox{max}}>S^{(2)}_{\mbox{max}}$, and (b) and (d):
$S^{(2)}_{\mbox{max}}>S^{(1)}_{\mbox{max}}$.

One can easily study the effect of other higher-order interactions using Eq.
(9). From Fig. 2(a), at the small values of $r$, (i.e., $\sigma_{2}$), the
value of the eigenratio approaches to
${S^{(1)}_{\mbox{max}}}/{S^{(1)}_{\mbox{min}}}$. As we increase the $r$, by
increasing $\sigma_{2}$, we can drastically reduce the values of $R$. In Fig.
2(a), ${S^{(1)}_{\mbox{max}}}/{S^{(1)}_{\mbox{min}}}=100$ and it falls down to
$3$ as we increase the three-body interactions. Therefore, higher-order
interactions promote synchronization when the first order intensities (degrees
in the case of unweighted networks) are very heterogeneous. The same effect is
shown in the values of cost as we can see from Fig. 2(c). Then we analyze the
synchronizability when the intensity of three-body interaction becomes more
heterogeneous than the pairwise interaction, i.e., when
$S^{(2)}_{\mbox{max}}>S^{(1)}_{\mbox{max}}$. The values are chosen as follows
as $S^{(2)}_{\mbox{max}}=1000$ and $S^{(1)}_{\mbox{max}}=20$. In this case,
the three-body interaction coupling strength has opposite effect. For a given
value of $S^{(2)}_{\mbox{max}}$, the eigenratio increases as we increase the
three-body interaction coupling strength $\sigma_{2}$. Hence, decrease the
synchronizability. The values of eigenratio from the approximated formula and
the general formula converges when the eigenratio is minimum.

The effect of coupling strength is the same for both cost and the eigenratio
[cf. Fig. 2(b) and 2(d)]. The more heterogeneity there is in the second-order
intensities, the faster the cost increases. Hence, the synchronizability
decreases if we increase the value of $\sigma_{2}$. The effect of coupling
strengths on synchronization depends on the intensity of the network topology.
Based on the above results, if $S^{(1)}_{\mbox{max}}>S^{(2)}_{\mbox{max}}$, an
increase in $\sigma_{2}$ enhances synchronizability. However, the opposite
scenario occurs if $S^{(2)}_{\mbox{max}}>S^{(1)}_{\mbox{max}}$, where a larger
$\sigma_{2}$ reduces synchronizability.

In general, higher-order interactions do not always promote synchronization
[21]. They inhibit the synchronization of the system when the weights of the
non-pairwise interactions are higher than those of the pairwise interactions.
With this knowledge, we can promote or inhibit synchronization by manipulating
the weights and coupling strengths of the pairwise and non-pairwise
interactions.

###  III.3 Behaviour of the constants $A_{R}$, $A_{C}$, $B_{R_{1}}$ and
$B_{R_{2}}$

Next, we analyze the behaviour of the constants $A_{R}$, $A_{C}$, $B_{R_{1}}$,
and $B_{R_{2}}$ as the mean degree changes for a simplicial complex with a
mean degree of $K$. For this purpose, we choose the values of the coupling
strengths as $\sigma_{1}=0.001$ and $\sigma_{2}=0.01$ and vary
${S^{(1)}_{\mbox{max}}/{S^{(1)}_{\mbox{min}}}}=1,2,10,100$ by fixing all nodes
to the same degree as in the case of $K$-regular networks. Using Eqs. (22) and
(23), we calculate the values of $A_{R}$, $A_{C}$, $B_{R_{1}}$, and
$B_{R_{2}}$ for large values of mean degree $K$.

As we can see from Figs. 3(a) and 3(b), the heterogeneity of the intensities
has a small effect on $A_{R}$, and the behaviour of $A_{C}$ appears
independent of the intensities’ heterogeneity as $K$ increases.

Figure 3: Variations of (a) $A_{R}$ (b) $A_{C}$, (c) $B_{R_{1}}$ and (d)
$B_{R_{2}}$ against $K$. When $K$ is very large, $A_{R}$ and $A_{C}$ approach
1 and $B_{R_{1}}$ and $B_{R_{2}}$ approach $r$, so we can use the approximate
formula rather than the general formula.

Both $A_{R}$ and $A_{C}$ quickly approach unity as $K$ increases.

Further, from Figs. 3(c) and 3(d), we observe that the values of $B_{R_{1}}$
and $B_{R_{2}}$ are highly dependent on the heterogeneity of the intensities.
$B_{R_{1}}$ remains close to $r=10$ when the intensities are not very
heterogeneous. As the heterogeneity increases, the values decrease for low
$K$s and slowly approach $r$. The values of $B_{R_{2}}$ remain close to the
upper bound when the intensities are heterogeneous and slowly approach $r$.
When the intensities are homogeneous, i.e.,
${S^{(1)}\mbox{max}}/{S^{(1)}\mbox{min}}=1$, the values of $B_{R_{2}}$ remain
close to $r$ for all $K$.

From the behaviours of $A_{R}$, $A_{C}$, $B_{R_{1}}$ and $B_{R_{2}}$, we can
observe that for $K\gg 1$ and not very heterogeneous networks, the general
formula approaches the approximated formula. With this, we can easily analyze
the synchronizability of the weighted simplicial complexes in the presence of
other higher-order interactions, such as four-body interactions, five-body
interactions etc., from the knowledge of mean degrees, coupling strengths, and
maximum and minimum intensities for the interaction networks.

###  III.4 Verification with Rössler oscillators

To validate our findings, we verify the applicability of the proposed model
(1) or (3) by analyzing complete synchronization in a system of identical
Rössler oscillators. In this case, we consider Eq. (3) to represent the
dynamical equations of a weighted simplicial complex composed of Rössler
oscillators as,

$\displaystyle\dot{x}_{i}=$
$\displaystyle\,-y_{i}-z_{i}+\sigma_{1}\frac{S^{(1)}_{i}}{K^{(1)}_{i}}{\sum_{j=1}^{N}a_{ij}^{(1)}(x_{j}-x_{i})}$
$\displaystyle+\sigma_{2}\frac{S^{(2)}_{i}}{K^{(2)}_{i}}\sum_{j=1}^{N}\sum_{k=1}^{N}a_{ijk}^{(2)}(x_{j}^{2}x_{k}-x_{i}^{3}),$
(27a) $\displaystyle\dot{y}_{i}=$ $\displaystyle\,x_{i}+ay_{i},$ (27b)
$\displaystyle\dot{z}_{i}=$
$\displaystyle\,b+z_{i}(x_{i}-c),\;\;i=1,2,\ldots,N,$ (27c)

where $a=0.2$, $b=0.2$ and $c=9.0$, and $N$ is the total number of
oscillators. We construct a weighted simplicial complex comprising $50$
Rössler oscillators. For a given node $i$, $K^{(1)}_{i}$ and $K^{(2)}_{i}$
represent the degrees of pairwise and three-body interactions, respectively,
indicating the total number of links and triangles to which the node belongs.
$S^{(1)}_{i}$ and $S^{(2)}_{i}$ denote the intensities of pairwise and three-
body interactions, respectively, for node $i$. The oscillators are coupled
through their $x$ components for both pairwise and three-body interactions
(later we will show that the results are independent on the choice of
components). We fix the values of $\sigma_{1}$, $\sigma_{2}$,
$S^{(1)}_{\mbox{max}}$, $S^{(2)}_{\mbox{max}}$, $S^{(2)}_{\mbox{min}}$ and
vary the values of $S^{(1)}_{\mbox{min}}$.

Figure 4: Synchronizability of coupled Rössler systems (27) calculated using
cost, eigenratio, and synchronization error ($E$) for different values of the
coupling strengths,
$(S^{(1)}_{\mbox{max}},S^{(2)}_{\mbox{max}},\sigma_{1},\sigma_{2})$: (a)
$(1,10^{-3},1,0.1)$ (first column), (b) $(10,10^{-2},0.1,10^{2})$ (second
column), and (c) $(100,1,10^{-2},10^{-4})$ (third column) chosen as per
equation (29) and $S^{(2)}_{\mbox{min}}=0.1\times S^{(2)}_{\mbox{max}}$. The
synchronizability of the network is independent of the individual values of
the coupling strengths.

The mean degrees are $K^{(1)}=28.4$ and $K^{(2)}=242.1$. To perform the
numerical analysis, we randomly distribute the degrees and the intensities. We
construct the adjacency tensors as we did previously. We calculate the
synchronisation error using the formula [18],

$\displaystyle
E=\left\langle\left(\frac{1}{N(N-1)}\sum_{i,j=1}^{N}\|\boldsymbol{x}_{j}-\boldsymbol{x}_{i}\|^{2}\right)^{\frac{1}{2}}\right\rangle_{T},$
(28)

where $\boldsymbol{x}$ is the state of the node, $\langle\rangle$ represents
the time average over time $T$.

Then, we calculate the eigenratio and cost using the formulas (22) and (23).
We consider three combinations of coupling strengths and maximum intensities,
and Fig. 4 shows the results. We choose the values for the
$S^{(1)}_{\mbox{max}}$, $S^{(2)}_{\mbox{max}}$, $\sigma_{1}$ and $\sigma_{2}$
as per the relations,

$\displaystyle S^{(1)}_{\mbox{max}}\times\sigma_{1}$ $\displaystyle\in[1,2],$
(29a) $\displaystyle S^{(2)}_{\mbox{max}}\times\sigma_{2}$
$\displaystyle\in[10^{-4},10^{0}].$ (29b)

We choose the parameters from the range of coupling strengths for
synchronization in unweighted simplicial complexes of Rössler oscillators [18]
and the results are shown in Fig. 4. We observe that as the eigenratio and
cost decrease, the synchronization error also decreases, leading to a more
synchronized system. Complete synchronization (i.e., zero synchronization
error) occurs when $S^{(1)}_{\mbox{min}}=0.2\times S^{(1)}_{\mbox{max}}$, or
20% of $S^{(1)}_{\mbox{max}}$. This indicates that we can assess the
synchronizability of complex systems with higher-order interactions using the
eigenratio and cost, employing the general formulas (22) and (23) when the
system’s mean degree is large. Additionally, we observe that a large
$S^{(1)}_{\mbox{max}}/S^{(1)}_{\mbox{min}}$ ratio hinders synchronization.
Therefore, to inhibit synchronization in systems where its effects are
undesirable, we can increase the maximum intensity or decrease the minimum
intensity; conversely, to promote synchronization, we can do the opposite.

We validate our formula by conducting additional tests, wherein we select
degrees using power-law and uniform distributions. The results are shown in
Fig. 5. We distribute the degrees and intensities as powerlaw and random
(P-R), uniform and random (U-R), and powerlaw and uniform (P-U) while
maintaining all other parameters as in Fig. 4(a). The results consistently
demonstrate that synchronization is independent of the connectivity patterns
in simplicial complexes. Therefore, the derived formula simplifies the problem
by reducing the necessity to compute adjacency tensors for each interaction.
Moreover, this formula facilitates a swift and straightforward determination
of when a system synchronizes.

Figure 5: Synchronizability calculated using (a) cost, (b) eigenratio, and (c)
synchronization error ($E$) for different distributions of the degrees and the
intensities: P-R, U-R, and P-U correspond to the degrees and intensities
distributed using power law and random, uniform and random, and power law and
uniform, respectively. The other parameters are $\sigma_{1}=1$,
$\sigma_{2}=0.1$, $S^{(1)}_{\mbox{max}}=1$, $S^{(2)}_{\mbox{max}}=10^{-3}$ and
$S^{(2)}_{\mbox{min}}=0.1\times S^{(2)}_{\mbox{max}}$. System synchronization
is independent of the distributions of degrees and intensities. Figure 6:
Synchronizability calculated using (a) cost, (b) eigenratio, and (c)
synchronization error ($E$) for different system sizes $N=50,70,150,200$ and
$\sigma_{1}=10^{-2}$, $\sigma_{2}=10^{-4}$, $S^{(1)}_{\mbox{max}}=100$,
$S^{(2)}_{\mbox{max}}=1$, and $S^{(2)}_{\mbox{min}}=0.1\times
S^{(2)}_{\mbox{max}}$. The complete synchronization of the system is
independent of the total number of nodes present in the network.

We then increase the number of oscillators to $N=70,100,150$ and $200$ while
maintaining all other parameters as in Fig. 4(b). We observe that
synchronization occurs at the same value of $S^{(1)}_{\mbox{min}}$ as in the
system with $N=50$ oscillators. This suggests that complete synchronization in
weighted simplicial complexes of identical oscillators is independent of the
total number of oscillators in the system. We plot the results in Fig. 6.
Notably, the synchronization errors for $N=150$ and $N=200$ are precisely the
same.

Figure 7: Synchronizability calculated using (a) cost, (b) eigenratio, and (c)
synchronization error ($E$) for the interactions happening between different
components. Here $\diamond$ represents the result for $x-y$ scheme (i.e.,
pairwise interaction through $x$ components and three body interactions
through $y$ components), $\times$ for $y-x$ (pairwise interaction through $y$
components and three-body interactions through $x$ variable) and $\circ$ for
$y-y$ (both interactions through $y$ components). The other parameters are
$\sigma_{1}=1$, $\sigma_{2}=0.1$, $S^{(1)}_{\mbox{max}}=1$,
$S^{(2)}_{\mbox{max}}=10^{-3}$ and $S^{(2)}_{\mbox{min}}=0.1\times
S^{(2)}_{\mbox{max}}$. The synchronizability of the system is independent of
the way the nodes are interacting.

Analyzing large weighted complex networks with higher-order interactions poses
a major challenge, especially due to the presence of a vast number of nodes in
natural and man-made systems. Our general formulas (22) and (23) can
significantly simplify the study of complete synchronization in such systems.

Next, we check the synchronizability by taking three different coupling
functions through different variables. We consider the pairwise interactions
through the $x$ components and three-body interactions through the $y$
components ($x-y$ scheme). We then conduct additional experiments with
pairwise interactions through $y$ components, three-body interactions through
$x$ components ($y-x$ scheme), and both interactions through $y$ components
($y-y$ scheme). All other values remain identical to those used in Fig. 4(a).
Figure 7 demonstrates that the system synchronizes at the same value of
$S^{(1)}_{\mbox{min}}$, indicating that synchronization in weighted simplicial
complexes is independent of the specific interaction modes among oscillators.

###  III.5 Verification with real-world connectivity structure

Up to now, we verified our derived analytical theories on synthetic network
structures. Finally, we verify our formula on real-world connectivity networks
[22]. The system is a collection of ants, where the nodes represent the ants,
and the links represent interactions between them within a day. The weight of
links represents the frequency of interaction. We consider all triangles as
three-body interactions and assign weights to each proportional to the weights
of the links in the triangles. At each time instant, there may be two agents
(nodes) that are not connected, i.e., the whole network is disconnected at
that time. However, the union of all the connections over the entire period of
time gives a connected network, which is a necessary condition for complete
synchronization in a graph. The real-world connectivity network has $N=152$,
$K^{(1)}_{\mbox{max}}=143$, $K^{(1)}_{\mbox{min}}=32$,
$K^{(2)}_{\mbox{max}}=7.1K$ with mean degree $K^{(1)}=102$, average number of
triangles $K^{(2)}=4.1K$, $S^{(1)}_{\mbox{max}}=2064$,
$S^{(2)}_{\mbox{max}}=5101.32$, $S^{(1)}_{\mbox{min}}=187$ and
$S^{(2)}_{\mbox{min}}=387.62$. Previously, the oscillatory dynamics was used
to represent the dynamics of each node in the modelling of opinion formation
in social systems [23] since the opinions are not necessarily stationary
always. In our case, we model the dynamics of the ant interaction network and
place chaotic Rössler oscillators on that network. We chose the values for
$\sigma_{1}$ and $\sigma_{2}$ according to the relations in Eq. (29). We fix
$\sigma_{1}=0.001$ and vary $\sigma_{2}$ from $10^{-5}$. Figure 8 displays the
results.

Figure 8: Synchronizability of real-world connectivity network [22] of $N=152$
oscillators modelled by Rössler systems calculated using (a) cost, (b)
eigenratio, and (c) synchronization error ($E$) by varying the nonpairwise
coupling strength $\sigma_{2}$.

One may note that as the eigenratio and the cost decrease, the synchronization
error also decreases. The system is synchronized when the values of the
eigenratio and the cost are at a minimum for the chosen values of the coupling
strengths. From these results, we can conclude that it is possible to predict
the complete synchronizability of complex systems with higher-order
interactions using our formula.

As a whole, it is intriguing to note that complete synchronization in large
complex networks is independent of the total number of nodes present in the
system and the way they are connected. Synchronization depends on the coupling
strengths and the degree of connection between nodes. If all nodes are
connected with nearly equal strengths, the system is easily synchronized.
However, if some nodes have significantly stronger connections compared to
others, achieving synchronization becomes challenging. By knowing the weights
of links and triangles, we can calculate intensities. Subsequently, coupling
strengths must be chosen based on the relationships mentioned earlier. This
allows us to predict when synchronization will occur, which is applicable to a
broad range of complex systems.

##  IV Conclusion

In this paper, we have analyzed the nature of synchronization in randomized
weighted simplicial complexes of identical oscillators with natural coupling.
We derived general formulas for the eigenratio and the cost and verified these
theoretical results for various network topologies with diffusive couplings
and intensity distributions. Our observations indicated that synchronization
is independent of the structure of simplicial complexes, the way the
oscillators interact with one another and the total number of oscillators
present in the system. We observed complete synchronization for fixed values
of $\sigma_{1}$, $\sigma_{2}$, $S^{(1)}_{\mbox{max}}$, $S^{(2)}_{\mbox{max}}$,
and $S^{(2)}_{\mbox{min}}$, when $S^{(1)}_{\mbox{min}}=0.2\times
S^{(1)}_{\mbox{max}}$ (i.e., $20$ percentage of $S^{(1)}_{\mbox{max}}$).
Notably, the synchronization error values for $N=150$ and $N=200$ are nearly
identical, indicating that synchronisation is independent of the number of
oscillators when $N$ is large. Our general formula provides a good
approximation for complete synchronization in these cases.

Furthermore, our observations highlighted that the effect of coupling
strengths varies depending on the largest intensity. When
$S^{(1)}_{\mbox{max}}>S^{(2)}_{\mbox{max}}$, an increase in $\sigma_{2}$
enhances synchronizability. Conversely, when
$S^{(2)}_{\mbox{max}}>S^{(1)}_{\mbox{max}}$, synchronization decreases with an
increase in $\sigma_{2}$. The behaviours of constants $A_{R}$, $A_{C}$,
$B_{R_{1}}$, and $B_{R_{2}}$ suggest that the approximated formula, containing
only coupling strengths and intensities, can be used.

We have validated our results by applying the general formula to determine the
synchronizability of weighted simplicial complexes of diffusively coupled
identical Rössler oscillators in Sec. III.4. We further verified our results
using real-world connectivity structure where the Rössler oscillators are
placed on each node in Sec. III.5. From numerical simulations in Sec. III.4
and III.5, we see that the results are robust on the number of nodes,
components through which they are coupled and the distributions of degree and
intensities.

In conclusion, we demonstrated that the synchronizability of large randomized
weighted simplicial complexes with natural coupling can be determined using
mean degrees, coupling strengths, and intensities without the need to
construct connectivity matrices or explicitly calculate their eigenvalues.
These formulas substantially reduce computation time and cost when assessing a
system’s tendency to synchronize. In essence, these formulas enable the
prediction and control of synchronizability in complex networks with higher-
order interactions, achieved by manipulating the degrees, weights, and
coupling strengths.

###### Acknowledgements.

The work of S.N.J. and P.M. is supported by MoE RUSA 2.0 (Bharathidasan
University - Physical Sciences).

## References

  * Boccaletti et al. [2006] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Complex networks: Structure and dynamics, Phys. Rep. 424, 175 (2006).
  * Kwapień and Drożdż [2012] J. Kwapień and S. Drożdż, Physical approach to complex systems, Phys. Rep. 515, 115 (2012).
  * Kuramoto [1975] Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, in International Symposium on Mathematical Problems in Theoretical Physics, edited by H. Araki (Springer Berlin Heidelberg, Berlin, Heidelberg, 1975) pp. 420–422.
  * Iacopini et al. [2019] I. Iacopini, G. Petri, A. Barrat, and V. Latora, Simplicial models of social contagion, Nat. Commun. 10, 1 (2019).
  * Strogatz and Stewart [1993] S. H. Strogatz and I. Stewart, Coupled oscillators and biological synchronization, Sci. Am. 269, 102 (1993).
  * Majhi et al. [2022] S. Majhi, M. Perc, and D. Ghosh, Dynamics on higher-order networks: A review, J. R. Soc. Interface 19, 20220043 (2022).
  * Battiston et al. [2021] F. Battiston, E. Amico, A. Barrat, G. Bianconi, G. Ferraz de Arruda, B. Franceschiello, I. Iacopini, S. Kéfi, V. Latora, Y. Moreno, et al., The physics of higher-order interactions in complex systems, Nat. Phys. 17, 1093 (2021).
  * Battiston et al. [2020] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Networks beyond pairwise interactions: structure and dynamics, Phys. Rep. 874, 1 (2020).
  * Boccaletti et al. [2023] S. Boccaletti, P. De Lellis, C. del Genio, K. Alfaro-Bittner, R. Criado, S. Jalan, and M. Romance, The structure and dynamics of networks with higher order interactions, Phys. Rep. 1018, 1 (2023).
  * Grilli et al. [2017] J. Grilli, G. Barabás, M. J. Michalska-Smith, and S. Allesina, Higher-order interactions stabilize dynamics in competitive network models, Nature 548, 210 (2017).
  * Mayfield and Stouffer [2017] M. M. Mayfield and D. B. Stouffer, Higher-order interactions capture unexplained complexity in diverse communities, Nat. Ecol. Evol. 1, 0062 (2017).
  * Alvarez-Rodriguez et al. [2021] U. Alvarez-Rodriguez, F. Battiston, G. F. de Arruda, Y. Moreno, M. Perc, and V. Latora, Evolutionary dynamics of higher-order interactions in social networks, Nat. Hum. Behav. 5, 586 (2021).
  * Anwar and Ghosh [2022] M. S. Anwar and D. Ghosh, Stability of synchronization in simplicial complexes with multiple interaction layers, Phys. Rev. E 106, 034314 (2022).
  * Anwar and Ghosh [2023] M. S. Anwar and D. Ghosh, Neuronal synchronization in time-varying higher-order networks, Chaos 33, 073111 (2023).
  * Baccini et al. [2022] F. Baccini, F. Geraci, and G. Bianconi, Weighted simplicial complexes and their representation power of higher-order network data and topology, Phys. Rev. E 106, 034319 (2022).
  * Courtney and Bianconi [2017] O. T. Courtney and G. Bianconi, Weighted growing simplicial complexes, Phys. Rev. E 95, 062301 (2017).
  * Zhou et al. [2006] C. Zhou, A. E. Motter, and J. Kurths, Universality in the synchronization of weighted random networks, Phys. Rev. Lett. 96, 034101 (2006).
  * 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, Stability of synchronization in simplicial complexes, Nat. Commun. 12, 1255 (2021).
  * Barabási and Bonabeau [2003] A.-L. Barabási and E. Bonabeau, Scale-free networks, Sci. Am. 288, 60 (2003).
  * Adamic and Huberman [2000] L. A. Adamic and B. A. Huberman, Power-law distribution of the world wide web, Science 287, 2115 (2000).
  * Zhang et al. [2023] Y. Zhang, M. Lucas, and F. Battiston, Higher-order interactions shape collective dynamics differently in hypergraphs and simplicial complexes, Nat. Commun. 14, 1605 (2023).
  * Rossi and Ahmed [2015] R. Rossi and N. Ahmed, The network data repository with interactive graph analytics and visualization, AAAI 29, 1 (2015).
  * Pluchino et al. [2006] A. Pluchino, S. Boccaletti, V. Latora, and A. Rapisarda, Opinion dynamics and synchronization in a network of scientific collaborations, Physica A 372, 316 (2006).