# Extreme sparsity gives rise to functional specialization

Gabriel Béna
Department of Electrical and Electronic Engineering
Imperial College London
London, UK
<EMAIL_ADDRESS>
&Dan F. M. Goodman
Department of Electrical and Electronic Engineering
Imperial College London
London, UK
<EMAIL_ADDRESS>

###### Abstract

Modularity of neural networks – both biological and artificial – can be
thought of either structurally or functionally, and the relationship between
these is an open question. We show that enforcing structural modularity via
sparse connectivity between two dense sub-networks which need to communicate
to solve the task leads to functional specialization of the sub-networks, but
only at extreme levels of sparsity. With even a moderate number of
interconnections, the sub-networks become functionally entangled. Defining
functional specialization is in itself a challenging problem without a
universally agreed solution. To address this, we designed three different
measures of specialization (based on weight masks, retraining and correlation)
and found them to qualitatively agree. Our results have implications in both
neuroscience and machine learning. For neuroscience, it shows that we cannot
conclude that there is functional modularity simply by observing moderate
levels of structural modularity: knowing the brain’s connectome is not
sufficient for understanding how it breaks down into functional modules. For
machine learning, using structure to promote functional modularity – which may
be important for robustness and generalization – may require extremely narrow
bottlenecks between modules.

##  1 Introduction

Modularity of neural networks is a bit like the notion of beauty in art:
everyone agrees that it’s important, but nobody can say exactly what it means.
Modularity has been widely observed in the brain (15; 20; 4; 14), suggested to
be important in evolution (18; 5; 10) and learning (7), and found in trained
artificial neural networks (8). One definition is that a network is modular if
it can be broken down into composable subnetworks that can operate
independently (1). This definition has a structural and functional component.
Structural, because we partition neurons in the network into different
modules. Functional, because these modules are then required to carry out
independent functions.

Structural notions of modularity can be easily measured by constructing a
graph of neurons as nodes and synapses as edges (3), for example the $Q$
modularity measure of how many edges are within modules compared to between
modules (16). Functional measures of modularity are much less straightforward.
How can you quantitatively measure whether a sub-network is independent or
composable? Do different proposed measures of functional modularity agree? Do
they agree with a structural measure?

We investigate these questions with a simple network architecture and set of
tasks designed to allow us to control the degree of structural modularity. We
then compare different measures of functional specialization of these modules
to certain sub-tasks. Reassuringly, we find that the three functional measures
qualitatively agree. However, we also find that only extreme levels of
structural modularity give rise to functional modularity. These results can
serve as a guideline for the design of artificial neural networks that want to
introduce functional modularity. However, they may also pose difficulties for
the interpretation of biological, connectomic data, since our results imply
that even relatively strong structural modularity may have no functional
bearing (20).

We describe the architecture and tasks in section 2, the measures of
functional specialization in section 3, and the results of our experiments in
section 4. We discuss the implications and limitations of these results in
section 5.

##  2 Architecture and tasks

Our approach to investigating the relationship between structural and
functional modularity is to explicitly design a network with a controllable
degree of structural modularity. We then train it on a task designed to be
composed of clearly defined sub-tasks, and measure the degree of functional
modularity in the trained network. The simplest possible network that could be
described as modular would have only two modules. Therefore, our architecture
consists of two sub-networks, each one being densely and recurrently
connected. We add sparse interconnections between neurons in these two sub-
networks, with a degree of sparsity ranging from a single connection in each
direction between the sub-networks, right up to dense inter-connectivity. With
dense inter-connectivity, the network ceases to be structurally modular as
every neuron is connected to every other neuron. Each sub-network receives an
MNIST digit (13) as input, and the task is to return one digit or the other
depending on whether the parity of the two digits is the same or different.
With this task, the two networks are required to communicate, but are able to
specialize in recognizing their own digit because it is sufficient to
communicate only the parity of their own digit to the other sub-network.

In this section we describe the architecture, tasks and training in more
detail. In section 3 we describe the measures of functional modularity.

###  2.1 Tasks

We designed a nested task to allow us to measure functional specialization of
sub-networks.

Sub-network $j$ receives as input an instance of digit
$D_{j}\in\left\\{0,...,9\right\\}$. We define the parity difference as
$P=(D_{0}+D_{1})\bmod 2$. The value of $P$ is 0 if the parity of the two
digits is the same, or 1 if it is different. The network needs to return digit
0 if the parity is the same, or digit 1 if different. Solving this task
requires only communicating a single parity bit between the sub-networks.

$\mathrm{Target}=PD_{0}+(1-P)D_{1}$ (1)

The sub-tasks are to predict only one of the input digits.

$\mathrm{Target}_{k}=D_{k},\quad k\in\left\\{0,1\right\\}$ (2)

Note that we exclude the case $D_{0}=D_{1}$ because the target is then
ambiguous. This case has parity difference $P=0$. To resolve the imbalance in
the distribution of parity differences, we also exclude one case where the
parity difference is $P=1$ by requiring $D_{1}-D_{0}\bmod 10\neq 1$. We
combine these two exclusions by requiring that $(D_{1}-D_{0})\bmod
10\not\in\\{0,1\\}$.

###  2.2 Architectures

The network consists of two dense recurrent sub-networks each receiving their
own separate input, with sparse connections between the two sub-networks (fig.
1(a)). Each sub-network $n\in\left\\{0,1\right\\}$ is composed of a single
hidden layer $\textbf{h}^{(n)}$ of $N=100$ neurons, receiving input through
weights $W^{(n)}_{ih}$ and bias $\beta^{(n)}_{ih}$ and recurrently connected
through weights $W^{(n)}_{hh}$ and bias $\beta^{(n)}_{hh}$. The two sub-
networks are then connected through sparse weights
$W_{sparse}^{0\xrightarrow{}1}$, $W_{sparse}^{1\xrightarrow{}0}$.

(a) Standard Architecture

(b) Bottleneck Architecture

Figure 1: Two network architectures used in this study. In both cases, the
network consists of two densely, recurrently connected sub-networks with
sparse interconnections. Each sub-network receives a separate MNIST digit as
input. In the bottleneck architecture (b), the output of each sub-network is
passed through a narrow bottleneck.

The digits are presented for $T=5$ time steps. Each sub-network is presented
with an input $x^{(n)}$. We write $m=1-n$ to be the index of the other sub-
network. The hidden state $h_{t}^{(n)}$ of sub-network
$n\in\left\\{0,1\right\\}$ at time-step $t$ is then:

$h_{t}^{(n)}=\tanh\left(W^{(n)}_{ih}x^{(n)}+\beta^{(n)}_{ih}+W^{(n)}_{hh}h_{t-1}^{(n)}+\beta^{(n)}_{hh}+W_{sparse}^{m\xrightarrow{}n}h_{t-1}^{(m)}\right).$
(3)

Note that this formulation could be easily adapted for other types of sub-
networks like LSTMs (19) or leaky RNN cells.

Each sub-network is densely connected to a readout layer $\mathbf{r}^{(n)}$
through weights $W^{(n)}_{hr}$ and bias $\beta^{(n)}_{hr}$:

$\mathbf{r}^{(n)}=W^{(n)}_{hr}\mathbf{h}^{(n)}_{T}+\beta^{(n)}_{hr}.$

Given the two readouts layer, we compute the actual decision of the global
model using a competitive _max_ decision. For each input sample, the sub-
network with the highest “certainty” (largest value) takes the decision.

The output layer of the whole network $\mathbf{r}^{out}$ is defined to be the
output layer of the sub-network with the largest maximum value. Let the index
of this sub-network be

$\mu=\operatorname{argmax}_{n}\max\mathbf{r}^{(n)}.$

Then the output layer is defined to be

$\mathbf{r}^{out}=\mathbf{r}^{(\mu)}.$

The prediction of the network is then

$d^{out}=\operatorname{argmax}\mathbf{r}^{out}.$

We use this architecture for the majority of the results shown. For one of the
functional modularity measures, we introduce an additional bottleneck layer
$\mathbf{b}^{(n)}$ in each sub-network . This bottleneck is composed of an
extra, small hidden layer of 5 neurons between the sub-network and the readout
(fig. 1(b)). In that case, the connection from state to output decomposes in
weight matrices $W^{(n)}_{hb}$, $W^{(n)}_{br}$ and biases $\beta^{(n)}_{hb}$,
$\beta^{(n)}_{br}$, and the equations become:

$\displaystyle\mathbf{b}^{(n)}$
$\displaystyle=W^{(n)}_{hb}\mathbf{h}^{(n)}_{T}+\beta^{(n)}_{hb}$
$\displaystyle\mathbf{r}^{(n)}$
$\displaystyle=W^{(n)}_{br}\mathbf{b}^{(n)}+\beta^{(n)}_{br}.$

In all cases, weights are initialized following a uniform distribution
$\mathcal{U}(-0.1,0.1)$

###  2.3 Structural modularity

We define the fraction of connections between two sub-networks as
$p\in[1/N^{2},1]$. The same fraction of connections is used in each direction.
The smallest value is $p=1/N^{2}$ corresponding to a single connection in each
direction, and the largest fraction $p=1$ corresponds to $N^{2}$ connections
in each direction (all-to-all connectivity).

We adapt the $Q$ metric for structural modularity (16) to a directed graph:

$Q=\frac{1}{M}\sum_{ij}(A_{ij}-P_{ij})\delta_{g_{i}g_{j}},$ (4)

where $M$ is the total number of edges in the network, $A$ is the adjacency
matrix, $\delta$ is the Kronecker delta, $g_{i}$ is the group index of node
$i$, and $P_{ij}$ is the probability of a connection between nodes $i$ and $j$
if we randomized the edges respecting node degrees. For our network, we can
analytically compute (appendix A):

$Q=\frac{1}{2}\cdot\frac{1-p}{1+p}.$ (5)

This varies from $Q=0$ when $p=1$ (all nodes connected, so no modularity) to
$Q=1/2$ for $p=0$ (no connections between the sub-networks, so perfect
modularity). Note that for $g$ equally sized groups in a network, the maximum
value that $Q$ can attain is $1-1/g$.

###  2.4 Training

The sub-networks are trained with gradient descent using the ADADELTA (22)
optimizer. The sparse connections in-between sub-networks are trained using
the Deep-R algorithm (2). Each weight is assigned a constant sign $s_{ij}$ and
a associated parameter $\phi_{ij}$. Any weight $w_{ij}$ is considered active
only if the corresponding parameter $\phi_{ij}>0$. If not the weight is said
to be dormant. For non-negative $\phi_{ij}$, the corresponding weight is thus
given by $w_{ij}=s_{ij}\phi_{ij}$. For others, weights are set to 0. Ensuring
sparsity of the weight matrix comes from ensuring that only the desired number
of values of $\mathbf{\phi}$ are positive at all times. We start with the
desired number of positive parameters based on the sparsity $p$, picked at
random. Only those parameters then undergo a gradient descent step, using a
separate ADAM optimizer (11), followed by a random walk step. Whenever one of
these becomes negative, effectively setting its weight to 0, we pick at random
a new one to become positive. This ensures that only the desired number of
weights is active at all times. The gradient descent reinforces important
weights by increasing their corresponding parameters. The random walk in
parameter space continuously explores possible connection patterns and allows
the network to ’rewire’ in case of changing objectives. Effectively, the
parameter responsible for the random walk is annealed from its initial value
to 0 when $p$ goes from 0 to 1. Note that our results are not substantially
changed if we remove the random walk exploration but we found that Deep-R
greatly improves performances at very high sparsities.

##  3 Functional modularity

Defining functional modularity is challenging. We want to capture the notion
that a sub-network is capable of doing a sub-task independently, which we call
functional specialization. However, if that network has never been trained
directly on the sub-task, and has always had input from neurons outside the
sub-network, it is unlikely to be able to do so without some re-training,
which changes the function of the sub-network. We considered three ways of
addressing this issue, leading to three quantitative measures of
specialization, described below. None of the measures are, individually,
entirely satisfactory. However, we show in section 4 that all three measures
qualitatively agree, suggesting that despite potential individual problems
with these measures, they are measuring something meaningful.

We first define three separate metrics, and then give a standard normalization
to turn them into comparable specialization measures in section 3.4.

###  3.1 Bottleneck metric

One straightforward approach to measuring specialization of a sub-network on a
particular sub-task would be to directly evaluate how well it performs on that
sub-task. In our case, that would be monitoring individual performance of sub-
networks when asked to predict either their own or the other sub-network’s
corresponding digits. However, without modification of the network, the sub-
networks would continue to predict digits based on the global parity as they
were initially trained to do. Re-training only the linear readout of sub-
networks while leaving the recurrent weights fixed allows for good performance
on this sub-task. However, due to the simplicity of the task, re-training the
whole hidden state to output readout is sufficiently powerful to allow for
excellent performance even if the recurrent weights are random and untrained,
so this cannot measure anything meaningful about the trained weights. The
issue is that the full state of the sub-network contains a large amount of
information about the original stimulus, and not just the result of any
computations carried out by the sub-network, which is the part we are
interested in.

To effectively constrain the amount of information readily available to the
readout layer, we introduced a narrow bottleneck of only 5 neurons in the
architecture, between the sub-network and the readout (fig. 1(b)). The state
of the bottleneck neurons is sufficient to carry all the information about the
digit class but insufficient to represent much information about the stimulus
itself. Re-training the bottleneck-output readout then provides us with a
measure of performance of the sub-network on the sub-task of classifying the
digit. This approach was inspired by the concept of the information bottleneck
(21) as it minimizes the mutual information between the bottleneck state and
the stimulus, while maximizing the mutual information between the bottleneck
state and the digit class.

Concretely, for sub-network $n$ on sub-task $k$, the bottleneck metric
$\mathcal{M}(n,k)$ is the mean accuracy after:

  1. 1.

Training the network on the main task.

  2. 2.

Freezing all parameters except the weights $W_{br}^{(n)}$ and biases
$\beta_{br}^{(n)}$ between the bottleneck and readout of sub-network $n$.

  3. 3.

Re-initializing and retraining only those weights and biases on the task of
recognizing digit $k$.

###  3.2 Weight mask metric

Another promising approach to measure specialization comes from the recently
discovered idea of weight masking. In order to determine which weights
contribute to which part of the overall computation of a network, a binary
mask is applied onto the network’s weights, effectively setting most of them
to 0. This mask is trained to maximize performance of the masked network on a
target objective, while the actual values of weights remain frozen. This idea
was applied to find high performing sub-networks of randomly initialized
networks (9; 17) and more recently to quantify functional modularity (6). Our
method is similar to this latter approach but uses a different algorithm (17)
to train and discover the masks, and does not exploit the resulting masks in
the same manner. The masks are trained as follows: each weight is fixed but is
attributed a differentiable score, with only the top $q\%$ scoring weights
being used in the forward pass. The training of the mask consists of adjusting
those scores to maximize the performance of the masked network on the desired
objective. Once a mask is discovered for a particular sub-task, the analysis
of the performance of the masks, what weights are retained and the
attributions of those weights among the sub-networks serve as metrics for
specialization.

It should be noted that binary weight masks can be used to discover high
performing sub-networks even in randomly initialized networks (supermasks). We
verified that the performance of the discovered masks were substantially
better when applied to already trained networks, confirming that these could
indeed be used to measure specialization.

The weight mask metric of sub-network $n$ on sub-task $k$ is the proportion of
weights, after training of the weight mask on predicting digit $k$, that are
present in sub-network $n$. Let $\theta=\theta^{0}\cup\theta^{1}$ to be the
set of parameters of the whole network trained on the full task, where
$\theta^{n}$ are the parameters in sub-network $n$. We define
$\mathcal{S}^{q}_{k}\subseteq\theta$ to be the subset of parameters selected
by a mask trained on predicting digit $k$ using only $q\%$ of parameters from
$\theta$. The weight mask metric is then defined as follows, where $|\cdot|$
indicates cardinality:

$\mathcal{M}(n,k)=\frac{|\mathcal{S}^{q}_{k}\cap\theta^{n}|}{|\mathcal{S}^{q}_{k}|}=\frac{|\mathcal{S}^{q}_{k}\cap\theta^{n}|}{q|\theta|}$
(6)

In this work, the weight mask is not applied to input weights $W_{ih}$ and
biases $\beta_{ih}$. We chose $q$ to be as high as possible while ensuring
$\mathcal{M}(n,n)=1$ in the most extreme case of sparsity. In our case, this
led to selecting 5% of the total model’s weights, $q=0.05$. We also ensured
that this value still led to good performance in the masks discovered.

###  3.3 Correlation metric

Our final approach is to directly examine hidden layers’ activities, and how
each sub-network’s activity depends on each digit. One way of doing this is to
examine the correlation of hidden states of sub-networks when only varying a
single digit out of the input pair. A sub-network highly specialized on its
corresponding digit should exhibit high correlation between states
corresponding to inputs where only the second, irrelevant digit is varying. To
measure the correlation we simply used the mean Pearson correlation
coefficient.

For sub-network $n$ on sub-task $k$, we compute the mean Pearson correlation
coefficient between the hidden states of sub-network $n$ for pairs of examples
$x$ and $x^{\prime}$ where digit $k$ of the two examples is the same. Writing
$\operatorname{corr}(x,y)$ for the Pearson correlation coefficient of vectors
$x$ and $y$:

$\mathcal{M}(n,k)=\mathbb{E}_{x,x^{\prime}}[\;\operatorname{corr}(h_{n}(x),h_{n}(x^{\prime})\;|\;\mbox{digit
$k$ of $x$ $=$ digit $k$ of $x^{\prime}$}\;].$ (7)

###  3.4 Specialization Measure

We have now defined our three metrics. Each of those is applied to each sub-
network for each sub-task. We write $\mathcal{M}(n,k)$ for a metric of sub-
network $n$ on task $k$. Our specialization measure is then defined as

$\mathcal{M}_{spec}(n)=\frac{\mathcal{M}(n,n)-\mathcal{M}(n,1-n)}{\mathcal{M}(n,n)+\mathcal{M}(n,1-n)}.$
(8)

So $\mathcal{M}(n,n)$ is sub-network $n$ on the task of recognizing its own
digit, and $\mathcal{M}(n,1-n)$ is sub-network $n$ recognizing the digit of
the other sub-network.

This specialization measure can range from 1, which corresponds to a fully
specialized sub-network on its corresponding sub-task, to -1, which indicates
full specialization on the other sub-network’s corresponding sub-task.

##  4 Results

Figure 2 shows how each of the three measures of functional specialization
discussed in section 3 varies with the structural modularity of the network,
measured either with the sparsity $p$ (left) or with the $Q$ modularity
measure (right). In each case, the qualitative pattern is the same: a sparser,
more structurally modular network leads to a greater degree of functional
specialization. This is reassuring in two senses.

Firstly, it shows that our three very different measures of functional
specialization qualitatively agree, suggesting they are measuring something
meaningful. Secondly, it shows that as we would hope, imposing structural
modularity is sufficient to induce functional modularity.

However, we also note that only extreme levels of sparsity or structural
modularity lead to a high degree of functional specialization across all three
measures. For the bottleneck and correlation metrics, we require $Q>0.49$ (or
sparsity $p<0.01$) before we start to see a rise in functional specialization.
Note that the highest value possible for $Q$ is 0.5 when there are only two
equally-sized sub-networks. This trend is less extreme for the weight mask
metric, where we start to see a rise from the start, but functional
specialization according to this measure still depends strongly non-linearly
on structural modularity, without substantially increasing until around
$Q>0.4$ ($p<0.11$).

In other words, across all measures, a high degree of structural modularity
(e.g. $Q=0.35$) is consistent with little to no functional specialization.

Figure 2: Functional specialization of sub-networks for different degrees of
structural modularity. Each row shows one of the three measures of functional
specialization defined in section 3. The data shown is the same on the left
and right columns, only the horizontal axis changes. On the left, we show the
degree of sparsity between the two sub-networks on a logarithmic scale. On the
right, we show the widely used $Q$ modularity measure. Lines indicate means
and shaded regions indicate one standard deviation around the mean over 10
repeats.

##  5 Discussion

We have demonstrated two substantial results. The first is that, by using a
carefully hand crafted example it is possible to systematically relate
structural and functional modularity. We showed that three entirely different
measures of functional specialization qualitatively agree in this case. This
allows us some confidence that these measures do reflect an important
underlying property of these networks, even if individually they have
limitations (see below). Furthermore, each measure of functional modularity
increases monotonically with increasing levels of structural modularity. This
reassuringly matches our intuition.

The second result, however, is that we require a very high degree of
structural modularity to ensure functional modularity. Or, put another way, a
moderately high degree of structural modularity can be seen in a network that
has almost no functional modularity whatsoever. This has implications both for
the design of artificial networks and for the analysis of biological networks.

For designing networks, it means that imposing a moderate degree of structural
modularity is unlikely to be sufficient, on its own, to ensure any degree of
functional modularity. If we want functional modularity, we either need to use
extreme levels of structural modularity, imposed via very sparse connections
between sub-networks, or – perhaps better – we need to use other methods to
induce this functional modularity.

In terms of analyzing biological data, our results demonstrate that we should
not conclude any degree of functional modularity simply by observing a
moderate degree of structural modularity. This potentially poses problems for
approaches such as connectomics (20), which assume that knowing the structural
properties of networks will give us good constraints on their functional
properties.

Despite good qualitative agreement between the measures of functional
modularity, each one individually has significant limitations individually. In
both bottleneck and weight-mask measures, for example, we effectively retrain
part of the network to perform a sub-task, and this training requires making
choices of parameters such as learning rate, number of epochs, etc. The
resulting metric depends on these choices.

Another limitation of our findings comes from the task design. In order to
easily witness and measure specialization, we designed a simple global task
that could decompose into sub-tasks. The information that needs to be
transmitted between the sub-networks, however, is minimal (a single bit). This
could explain why we only see specialization at extreme levels of sparsity. An
interesting future work would be to craft a global task where the amount of
information transmission needed among sub-networks in order to solve the task
could be independently controlled. We could then examine how the relationship
between sparsity (structural modularity) and specialization (functional
modularity) depends on the required communication bandwidth. Would it shift,
but preserve the overall pattern?

In the end modularity is really like beauty in art. If you give it time, and
different perspectives, you can start to discern its ever-changing shape.

## References

  * [1] Mohammed Amer and Tomás Maul.  A review of modularization techniques in artificial neural networks.  Artificial Intelligence Review, 52(1):527–561, Apr 2019.
  * [2] Guillaume Bellec, David Kappel, Wolfgang Maass, and Robert Legenstein.  Deep rewiring: Training very sparse deep networks, 2018.
  * [3] Edward T. Bullmore and Danielle S. Bassett.  Brain Graphs: Graphical Models of the Human Brain Connectome.  Annual Review of Clinical Psychology, 7(1):113–140, March 2011\.  Publisher: Annual Reviews.
  * [4] Zhang J. Chen, Yong He, Pedro Rosa-Neto, Jurgen Germann, and Alan C. Evans.  Revealing modular architecture of human brain structural networks by using cortical thickness from MRI.  Cerebral Cortex (New York, N.Y.: 1991), 18(10):2374–2381, October 2008.
  * [5] Jeff Clune, Jean-Baptiste Mouret, and Hod Lipson.  The evolutionary origins of modularity.  Proceedings of the Royal Society B: Biological Sciences, 280(1755):20122863, Mar 2013.
  * [6] Róbert Csordás, Sjoerd van Steenkiste, and Jürgen Schmidhuber.  Are neural nets modular? inspecting functional modularity through differentiable weight masks, 2021.
  * [7] Kai Olav Ellefsen, Jean-Baptiste Mouret, and Jeff Clune.  Neural Modularity Helps Organisms Evolve to Learn New Skills without Forgetting Old Skills.  PLOS Computational Biology, 11(4):e1004128, April 2015.  Publisher: Public Library of Science.
  * [8] Daniel Filan, Stephen Casper, Shlomi Hod, Cody Wild, Andrew Critch, and Stuart Russell.  Clusterability in Neural Networks.  arXiv:2103.03386 [cs], March 2021.  arXiv: 2103.03386.
  * [9] Jonathan Frankle and Michael Carbin.  The lottery ticket hypothesis: Finding sparse, trainable neural networks, 2019.
  * [10] Nadav Kashtan and Uri Alon.  Spontaneous evolution of modularity and network motifs.  Proceedings of the National Academy of Sciences, 102(39):13773–13778, September 2005.  Publisher: National Academy of Sciences Section: Physical Sciences.
  * [11] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * [12] Simon Kornblith, Mohammad Norouzi, Honglak Lee, and Geoffrey Hinton.  Similarity of neural network representations revisited, 2019.
  * [13] Yann LeCun and Corinna Cortes.  MNIST handwritten digit database.  2010\.
  * [14] David Meunier, Sophie Achard, Alexa Morcom, and Ed Bullmore.  Age-related changes in modular organization of human brain functional networks.  NeuroImage, 44(3):715–723, February 2009.
  * [15] V. B. Mountcastle.  The columnar organization of the neocortex.  Brain: A Journal of Neurology, 120 ( Pt 4):701–722, April 1997\.
  * [16] M. E. J. Newman.  Modularity and community structure in networks.  Proceedings of the National Academy of Sciences, 103(23):8577–8582, May 2006.
  * [17] Vivek Ramanujan, Mitchell Wortsman, Aniruddha Kembhavi, Ali Farhadi, and Mohammad Rastegari.  What’s hidden in a randomly weighted neural network?, 2020.
  * [18] C. Redies and L. Puelles.  Modularity in vertebrate brain development and evolution.  BioEssays, 23(12):1100–1111, 2001.
  * [19] Haşim Sak, Andrew Senior, and Françoise Beaufays.  Long short-term memory based recurrent neural network architectures for large vocabulary speech recognition, 2014.
  * [20] Olaf Sporns, Giulio Tononi, and Rolf Kötter.  The Human Connectome: A Structural Description of the Human Brain.  PLOS Computational Biology, 1(4):e42, September 2005.  Publisher: Public Library of Science.
  * [21] Naftali Tishby and Noga Zaslavsky.  Deep learning and the information bottleneck principle.  2015 IEEE Information Theory Workshop, ITW 2015, 03 2015.
  * [22] Matthew D. Zeiler.  Adadelta: An adaptive learning rate method, 2012.

##  Appendix A Appendix

###  A.1 Q measure derivation

[16] defines the modularity measure $Q$ for an undirected graph. We modify the
definition for a directed graph as:

$Q=\frac{1}{M}\sum_{ij}(A_{ij}-P_{ij})\delta_{g_{i}g_{j}}=\frac{1}{M}\sum_{ij}\left(A_{ij}-\frac{k_{i}^{out}k_{j}^{in}}{M}\right)\delta_{g_{i}g_{j}}$
(9)

where $M$ is the total number of edges in the network, $A$ is the adjacency
matrix, $\delta$ is the Kronecker delta, $g_{i}$ is the group index of node
$i$, $k_{i}^{out}/k_{j}^{in}$ are the out-degree/in-degree of node $i/j$
respectively, used to compute $P_{ij}$: the probability of a connection
between nodes $i$ and $j$ if we randomized the edges respecting node degrees.
For our network, we have 2 sub-networks (with group index 0 and 1
respectively) of N neurons, densely connected, with p% active inter-
connections between the two. From this we get:

$\displaystyle M=2N^{2}(1+p)$ $\displaystyle\forall i\in[1,2N],\quad
g_{i}=\begin{cases}0,&\text{if }i\in[0,N-1]\\\ 1,&\text{if
}i\in[N,2N]\end{cases}$ $\displaystyle\forall(i,j)\in[0,2N]^{2},\quad
A_{i,j}=\begin{cases}1,&\text{if }\delta_{g_{i}g_{j}}=1\\\ 0\text{ or
}1,&\text{if }\delta_{g_{i}g_{j}}=0\end{cases}$
$\displaystyle\forall(i,j)\in[0,2N]^{2},\quad\begin{cases}k_{i}^{out}=\sum_{j^{\prime}}A_{ij^{\prime}}\approx
N(1+p)\\\ k_{j}^{in}=\sum_{i^{\prime}}A_{i^{\prime}j}\approx
N(1+p)\end{cases}$ And
$\displaystyle\begin{cases}\sum_{i}k_{i}^{out}=\sum_{i}\sum_{j^{\prime}}A_{ij^{\prime}}=N^{2}(1+p)\\\
\sum_{j}k_{j}^{in}=\sum_{j}\sum_{i^{\prime}}A_{i^{\prime}j}=N^{2}(1+p)\end{cases}$
Thus : $\displaystyle
Q=\frac{1}{2N^{2}(1+p)}\sum_{ij}\left(A_{ij}-\frac{(1+p)}{2}\right)\delta_{g_{i}g_{j}}$
$\displaystyle Q=\frac{1}{2N^{2}(1+p)}(2N^{2})\left(1-\frac{(1+p)}{2}\right)$
$\displaystyle\boxed{Q=\frac{1-p}{2(1+p)}}$

###  A.2 Training parameters

Here we present all training parameters used in the global model training
(parity task), the discovery of the weight masks and the bottleneck-readout
retraining. When two optimizers or learning rates are specified, they
correspond to the parameters for the sub-networks/sparse-connections
respectively. When a parameters is indicated as a range, it was geometrically
decreased in between the two bounds when inter-connections sparsity spanned
from 0 to 1. In the case of the global training, we had to train both
architectures (figs. 1(a) and 1(b)), we thus indicates parameters (bottleneck
size and number of epochs) for the absence/presence of the bottleneck.

| Global training | Mask training | Bottleneck re-training
---|---|---|---
Number of sub-networks | $2$ | $2$ | $2$
Number of input neurons | $2*784$ | $2*784$ | $2*784$
Number of hidden neurons | $2*100$ | $2*100$ | $2*100$
Bottleneck size | -/5 | - | $5$
Number of classes | $10$ | $10$ | $10$
Number epochs | $5/7$ | $1$ | $5$
Optimizer | Adadelta/Adam | SGD | Adadelta
Learning rate | $0.5/[0.1-0.01]$ | $0.1$ | $0.5$
Betas | $0.9/0.95$ | - | $0.9$
Table 1: Network and training parameters

###  A.3 Raw metrics

In order to measure specialization, we used a normalized difference between
sub-tasks over three different metrics. This difference however isn’t a fool-
proof way to arrive at out measure. No matter the raw performance of a sub-
network for a given metric, if it achieves the same for both sub-tasks then
our specialization measure is equal to 0. The question that a sub-network
performing perfectly for both sub-task should get the same measure as one
performing very poorly on both is thus open. Fortunately this was not the case
as shown in the raw metrics (fig. 3), and the quantity we measured is
meaningful. We also show the results of the correlation metric when using the
CKA distance [12] (note that this metric is defined as a distance, i.e. 0 when
highly correlated).

Figure 3: Raw metrics.

###  A.4 Global model performance

In this study, we measure the specialization of sub-networks for varying
sparsity of inter-connections. It is to be noted that this specialization
depends on the overall performance of the global model, which in turns should
depend on said sparsity of connections. We made sure that every trained model
showed good performance so that measures on functional modularity were
meaningful (fig. 4).

Figure 4: Trained models performance on global task.


#  ${Q}\bar{Q}$ $({Q}\in{\\{b,c}\\})$ spectroscopy using the modified
Rovibrational model

Zheng-Yuan Fang1,2<EMAIL_ADDRESS>Ya-Rong Wang1,2
<EMAIL_ADDRESS>Cheng-Qun Pang1,2,3111Corresponding author
<EMAIL_ADDRESS>1College of Physics and Electronic Information Engineering,
Qinghai Normal University, Xining 810000, China
2Joint Research Center for Physics, Lanzhou University and Qinghai Normal
University, Xining 810000, China
3Lanzhou Center for Theoretical Physics, Key Laboratory of Theoretical Physics
of Gansu Province, Lanzhou University, Lanzhou, Gansu 730000, China

###### Abstract

Mass spectra of quarkonium systems can be described by different
phenomenological potentials. In the present work, the resonance states of
heavy quarkonium like ($c\bar{c}$ and $b\bar{b}$) are considered as the
rovibrational states. We study a parameterized rovibrational model derived
from the empirical solution of the nonrelativistic Schrödinger equation with
Morse potential, the corrections are composed of colour hyperfine interaction
and spin-orbit interaction of mesons. We obtain the high excited state mass
spectra of charmonium and bottomonium comparing the results in reasonable
agreement with the present experimental data.

###### pacs:

14.40.Be, 12.38.Lg, 13.25.Jx

##  I Introduction

Since the discovery of charm quarks and bottom quarks in 1974 and 1977,
respectively, heavy quarkonium has become an influential and attractive
research field because its physical processes cover the whole energy range of
Quantum Chromodynamics (QCD). This energy range provides us an excellent place
to study the properties of perturbative and non-perturbative QCD Aubert  _et
al._ (1974); Herb  _et al._ (1977); Innes  _et al._ (1977). With the
improvement of experimental accuracy, Particle Data Group (PDG) has collected
dozens of heavy quarkonium particles. The study of hadronic energy spectrum is
one of the main tasks at present. We classify it and then study higher radial
and orbital excited states of heavy quarkonium. In previous work, a series of
phenomenological potential models can use specific techniques to describe the
mass spectra of quarkonium and are in good agreement with the experimental
values Fischer  _et al._ (2015); Lakhina and Swanson (2006); Radford and Repko
(2007); Soni  _et al._ (2018); Ebert  _et al._ (2011); Deng  _et al._ (2017);
Li and Chao (2009); Shah  _et al._ (2012a); Bai-Qing and Kuang-Ta (2009).
Abundant experimental data can accurately fit different parametric potential
models, and these models are consistent in the scale range of observable
states by solving the Schrödinger equation. For example, Barnes et al, used
the standard coulomb potential and linear potential to calculate the energy
spectrum and electromagnetic radiation decay of charmonium Barnes  _et al._
(2005). Considering the color screening effect, Jun-Zhang Wang et al,
calculated the energy spectrum and two-body strong decays of $c\bar{c}$ and
$b\bar{b}$ states, respectively Wang  _et al._ (2019, 2018). The meson family
high excited state spectroscopy describes the most familiar Reggie trajectory
Ebert  _et al._ (2013); Chen (2018, 2021). In particular, the latest
discussion on the overall framework of hadronic spectrum also provides a basis
for our calculation Albaladejo  _et al._ (2021).

Due to the abundant experimental data of heavy mesons, the rovibrational model
derived from the empirical solution of the nonrelativistic Schrödinger
equation with morse potential is only used to describe the excited spectrum of
hadrons rather than to calculate the effective potential between quarks.
Bernardo and Bastos calculated the mass of partial $\pi$, $K$, $N$, $\Sigma$
orbital and radial resonance states, which were in good agreement with the
experimental results Bernardo  _et al._ (2021). Rovibrational model takes into
account the effect of colour hyperfine interaction and spin-orbit interaction,
and has more systematic description of quark model. In the model, the hadronic
family has the resonance states with the same isospin type corresponding to
the fitting parameters of the same set of experimental particles, and we
further predict the mass spectra of other high excited hadrons.

In 1966, Paulie noticed that more than half of the baryon resonance can be
described as the vibration and rotation of the semi-rigid rotor Pauling
(1966). Bernardo and Bastos consider the non-harmonic and semi-rigid rotation
correction of the vibrational rotation. However, when there are many families
of hadron resonances, the predicted values are somewhat different from those
predicted by the traditional quark model. In this work, we consider the
rovibrational model with spin-orbit coupling, spin-spin coupling and tensor
items to calculate the mass spectra of charmonium and bottomonium, where
$Y(4220)$ is calculated as $\psi(4^{3}S_{1})$ Wang  _et al._ (2019); Deng  _et
al._ (2017); Shah  _et al._ (2012a); Li and Chao (2009); Barnes  _et al._
(2005). Our work is better in comparison with available experimental than
other theoretical values.

This paper is aimed to give a systemic study of excited states of charmonium
and bottomonium mass spectroscopy. Verify the feasibility of the model. Morse
potential method of diatomic molecular mass spectrum is mature, after
correction to calculate the double heavy quark particle mass spectrum also
achieved coincident with the experimental.

This paper is organized as follows. In Sec. II, the models employed in this
work are briefly reviewed. The mass spectrum of charmonium and bottomonium
will be performed in Sec. III. The paper ends with a conclusion in Sec. IV.

##  II Models employed in the work

In this work, the modified rovibrational model is utilized to calculate the
mass spectrum of the charmonium and bottomonium meson family. A comparison of
calculated rovibrational states with the experimental and other
phenomenological potential models results is predicted. In the following,
these models will be illustrated briefly.

###  II.1 The modified rovibrational model

In 1929, the exact solutions of the Schrödinger equation for the motion of
atomic nuclei in diatomic molecules are given Morse (1929). It is noteworthy
that approximate solutions of Schrödinger equation with some diatomic
molecular interactions using Nikiforov-Uvarov(NU) method can obtained the
energy eigenvalues and the total normalized wave function Okon  _et al._
(2017). Hesham Mansour et al. choose a phenomenological potential in the
framework of the nonrelativistic Schrödinger equation with relativistic
corrections using the Nikiforov-Uvarov (NU) method obtained the meson mass
spectroscopy Mansour and Gamal (2022). The energy level behavior of diatomic
molecules can be described by the stationary Schrödinger equation with the
internuclear potential in the form of morse potential, which considers the
rotation-vibration coupling including anharmonic and semi-rigid rotor
corrections Bernardo  _et al._ (2021). In case of charmonium and bottomonium
meson, the quark and antiquark separated by a distance $r$ . The radial
equation is given as follows:

$\begin{split}\frac{d^{2}R(r)}{dr^{2}}+\frac{L(L+1)R(r)}{r^{2}}=\frac{-8\pi^{2}m}{h^{2}}[E-D_{e}-D_{e}e^{-2\alpha(r-r_{e})}\\\
+2D_{e}e^{-\alpha(r-r_{e})}]R(r),\end{split}$ (2.1)

where $R(r)$ is the radial wave function, $E$ is the non-relativistic total
energy of the system, $h$ is the Planck’s constant, $m$ is the mass of the
particle, $L$ is the orbital quantum number,
$D_{e}-D_{e}e^{-2\alpha(r-r_{e})}+2D_{e}e^{-\alpha(r-r_{e})}$ is the popular
morse potential, $D_{e}$ and $\alpha$ are parameters characteristic to each
particle. $r_{e}$ is equated with the equilibrium internuclear separation for
diatomic molecules. The morse potential will be used here to describe only the
excited spectrum of hadrons and not to calculate the effective potential of
the quarks. Comparison of Cornell potential model and Godfrey–Isgur (GI) model
for calculating charmonium and bottomonium mass spectrum with considered
colour hyperfine interaction and spin-orbit in teraction. Using the Bernardo
and Bastos solution Bernardo  _et al._ (2021), which includes the
rovibrational coupling and corrections composed of spin-orbit coupling term,
spin-spin coupling term and tensor term for states with $L>0$, the mass of the
resonances $M(n,L)$ can be calculated as

$\begin{split}M(n,L)=D_{e}+\hbar\omega(n+\frac{1}{2})-\hbar\omega\chi_{e}(n+\frac{1}{2})^{2}+B_{rot}L(L+1)\\\
-D_{rot}L^{2}(L+1)^{2}-\alpha_{e}(n+\frac{1}{2})L(L+1)+\frac{1}{2}a_{e}(S(S+1)-\frac{3}{2})\\\
+\frac{1}{2}b_{e}(J(J+1)-S(S+1)-L(L+1))+c_{e}\bm{T},\end{split}$ (2.2)

where $n$ is the radial quantum number, $\hbar\omega(n+\frac{1}{2})$ is the
harmonic contribution, $\hbar\omega\chi_{e}(n+\frac{1}{2})^{2}$ is the
harmonic correction, $B_{rot}L(L+1)$ is the rigid rotor contribution,
$D_{rot}L^{2}(L+1)^{2}$ represents the semi-rigid rotor correction, and the
rovibrational coupling is given by $\alpha_{e}(n+\frac{1}{2})L(L+1)$, which is
the Coriolis effect. Consideration of orbital and radial high excited states
in the charmonium and bottomonium system.
$\frac{1}{2}a_{e}(S(S+1)-\frac{3}{2})$ and $c_{e}\bm{T}$ represents the colour
hyperfine interactions correction, which the GI nonrelativistic quark model,
there are typical forms of contact interaction and tensor interaction Godfrey
and Isgur (1985). $\frac{1}{2}b_{e}(J(J+1)-S(S+1)-L(L+1))$ denotes the spin-
orbit interaction due to one-gluon exchange correction. $\bm{S_{Q}}$ and
$\bm{S_{\bar{Q}}}$ are the spins of quarks in meson component, and the total
spin is $\bm{S}=\bm{S_{Q}}+\bm{S_{\bar{Q}}}$. The spin-orbit operator is
diagonal in a $|\bm{J,S,L}\rangle$ basis, spin-dependent terms with the matrix
elements includes

$\displaystyle\langle\bm{S_{Q}}\cdot\bm{S_{\bar{Q}}\rangle}=\frac{1}{2}(S(S+1)-\frac{3}{2}).$
(2.3)

We can get the relationship between spin and orbital coupling by $6j$ symbol,
the spin-orbit interaction correction terms with the matrix elements includes

$\displaystyle\langle\bm{L\cdot{S}}\rangle=\frac{1}{2}(J(J+1)-S(S+1)-L(L+1)).$
(2.4)

The tensor operator $\bm{T}$ has nonvanishing diagonal matrix elements only
between $L>0$ spin-triplet states, $\bm{r}$ is any unit vector. which are
Barnes  _et al._ (2005)

$\displaystyle\bm{T}=(\frac{\bm{S_{Q}\cdot{r}S_{\bar{Q}}\cdot{r}}}{r^{2}}-\frac{1}{3}\bm{S_{Q}\cdot
S_{\bar{Q}}}),$ (2.5)
$\displaystyle\langle^{3}L_{J}|\bm{T}|^{3}L_{J}\rangle=\left\\{\begin{array}[]{l}-\frac{1}{6(2L+3)},J=L+1;\\\
+\frac{1}{6},J=L;\\\ -\frac{L+1}{6(2L-1)},J=L-1.\end{array}\right.$ (2.9)

The parameter $B_{rot}=\frac{\hbar^{2}}{2mr_{e}^{2}}$ determined using the
experimental mass correlation particles, $D_{e}$ and coefficients
$\hbar\omega$, $\hbar\omega\chi_{e}$, $D_{rot}$, $B_{rot}$, $\alpha_{e}$,
$a_{e}$, $b_{e}$, $c_{e}$ as a function are determined by fitting the
corresponding hadrons system experimental results of the radial and orbital
quantum numbers, respectively. Due to the calculation of highly excited hadron
spectrum, and the lowest mass ground state has significant influence on the
system. It is noteworthy that the ground state $\eta_{c}(1^{1}S_{0})$ of
charmonium and the ground state $\eta_{b}(1S)$ of bottomonium do not
participate in the fitting of the high excited state results more in line with
the experimental values. In this case, in order to make the whole fitting
systematic, we choose to fit both with the orbital and radial excited states
that have been observed experimentally. The behavior of the rovibrational
parameters of Eq. (2.2) is similar to that of diatomic molecule states, and
for the meson system with double heavy quarks, the theoretical results
obtained after considered the influence of the colour hyperfine interaction
and spin-orbit interaction of mesons on the high excited state was better
agreement with the experimental results. This is also the original intention
to plus these perturbations. For molecules and hadrons, the behavior
parameters of the corresponding meson system depend on the experimental
results of the system. For $L=0$, there is no rotation effect, the
contribution of rovibrational model are mainly concentrated in $D_{e}$ and the
harmonic contribution, the anharmonic correction and the perturbation effects
are relatively weak. For $L>0$, which overall calculation results of the
hadrons system are in agreement with the experimental data.

##  III Numerical results and phenomenological analysis

We have determined the free parameters of the model by minimizing the quantity

$\displaystyle\chi^{2}=\sum\limits_{i}(\frac{E^{th}_{i}-M^{exp}_{i}}{M^{exp}_{i}})^{2}.$
(3.1)

The parameters of charmonium fitted from fifteen experimental values.
Experimentally observed mass values of bottomonium which are used to calculate
the minimu $\chi^{2}$ and to obtain the model parameters listed in the third
column of Table 1. Applying the modify rovibrational model and the parameters
in Table 1, the mass spectra of the charmonium and bottomonium can be
obtained, as shown in Table 3 and Table 2.

Table 1: Parameters of the modified rovibrational model, which are determined by fitting the charmonium and bottomonium system experimental data listed in PDG averages.  Parameter | Value $(c\bar{c})$ | Value $(b\bar{b})$
---|---|---
$D_{e}$ (MeV) | 2712.52 | 9171.35
$\hbar\omega$ (MeV) | 767.05 | 597.70
$\hbar\omega\chi_{e}$ (MeV) | 95.43 | 48.89
$B_{rot}$ (MeV) | 273.02 | 306.82
$D_{rot}$ (MeV) | 21.20 | 25.54
$\alpha_{e}\ $(MeV) | 39.12 | 69.77
$a_{e}$ (MeV) | 40.90 | 8.55
$b_{e}$ (MeV) | 12.27 | 13.07
$c_{e}$ (MeV) | 134.42 | 26.81

####  III.0.1 The spectra of charmonium meson excitations

The spectra of charmonium meson orbital and radial excitations are calculated,
and the values are listed in Table 2. All experimental masses are taken from
PDG averages Zyla  _et al._ (2020). The third radial excited state of
$\eta_{c}(3S)$ has a mass of 4.003 GeV, which is smaller than the result of GI
model and close to that reported in Ref Barnes  _et al._ (2005). The
difference between theoretical results for radial excited states
$\eta_{c}(3S)$ and $\eta_{c}(4S)$ are about 0.2 GeV. For the ground state of
$\psi(1^{3}D_{1})$ close to the experimental values, its second and third
radial excited states ($\psi(2^{3}D_{1})$ and $\psi(3^{3}D_{1})$) have the
mass of 4.122 GeV and 4.273 GeV, respectively. When evaluating the theoretical
results of all states, the qualitative evaluation errors selection of all
states are $\mathcal{V}^{Er}_{i}$=5.0 MeV, which are larger than their
respective experimental uncertainty. The reason is that the experimental error
of these particles is relatively small and the distribution is unevenly.
Finally, we obtain that error of the whole system is $\chi^{2}/{N}$ is 55.8,
which is consistent with the similar theory model Shah  _et al._ (2012a). Our
work is better in comparison with available experimental than other values in
the literatures Fischer  _et al._ (2015); Lakhina and Swanson (2006); Soni
_et al._ (2018); Ebert  _et al._ (2011); Deng  _et al._ (2017); Li and Chao
(2009); Shah  _et al._ (2012a); Barnes  _et al._ (2005); Wang  _et al._
(2019); Mansour and Gamal (2022); Shah  _et al._ (2012b); Patel  _et al._
(2016).

Table 2: The mass spectra of the experimental missing orbital and radial excited states have been predicted. Here number of experimental values ($N=15$), we list the error analysis is $\chi^{2}=$$\sum\limits_{i}$$(\frac{E^{th}_{i}-M^{exp}_{i}}{\mathcal{V}^{Er}_{i}})^{2}$, where $Y(4220)$ is calculated as $\psi(4^{3}S_{1})$ state. The values of underline in $\psi(4^{3}S_{1})$ state denote the predicted results of different phenomenological models after fitting the adopt experimental data of underline in parentheses, respectively. The parentheses in the reciprocal second column show the relative error between the theoretical and experimental values after fitting each experimental data of charmonium in the model. All results are in units of MeV.  $\begin{array}[]{cccccccccccccccc}\hline\cr\hline\cr\hline\cr\hline\cr\text{State}&\text{This work}&\text{NU \cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{mansour2022meson}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{ebert2011spectroscopy}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{PhysRevD.86.034015}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{PhysRevD.95.034026}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{soni2018q}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{PhysRevD.79.094004}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{PhysRevD.75.074031}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{PhysRevD.86.034015}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{PhysRevD.74.014012}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Patel_2016}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{JZ \cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{wang2019constructing}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{GI \cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{barnes2005higher}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{PDG \cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Zyla:2020zbs}{\@@citephrase{(}}{\@@citephrase{)}}}}&\text{Error}\\\ \hline\cr\eta_{c}(1^{1}S_{0})&{3041}&3041&2930&2981&2984&2989&2979&2980&2980&3088&2979&2981&2982&2984\ (1.9$\%$)&5.0\\\ \psi(1^{3}S_{1})&{3075}&3140&3110&3096&3097&3094&3097&3097&3097&3168&3096&3096&3098&3097\ (0.71$\%$)&5.0\\\ \eta_{c}(2^{1}S_{0})&{3618}&3661&3680&3635&3637&3602&3623&3597&3633&3669&3600&3642&3623&3638\ (0.55$\%$)&5.0\\\ \psi(2^{3}S_{1})&{3651}&3702&3680&3685&3679&3681&3673&3685&3690&3707&3680&3683&3676&3686\ (0.95$\%$)&5.0\\\ \eta_{c}(3^{1}S_{0})&{4003}&4135&-&3989&4004&4058&3991&4014&3992&4067&4011&4013&4064&-&\\\ \psi(3^{3}S_{1})&{4037}&4050&3800&4039&4030&4129&4022&4095&4030&4094&4077&4035&4100&4039\ (0.05$\%$)&5.0\\\ \eta_{c}(4^{1}S_{0})&{4197}&4414&-&4401&4264&4448&4250&4433&4244&4398&4397&4260&4225&-&\\\ \psi(4^{3}S_{1})&{4231}&$\text@underline{4419}$&-&$\text@underline{4427}$&4281&$\text@underline{4514}$&4273&$\text@underline{4477}$&4273&$\text@underline{4420}$&$\text@underline{4454}$&4274&4225&4230\ ($\text@underline{4421}$)&5.0\\\ h_{c}(1^{1}P_{1})&{3467}&3518&3430&3525&3526&3470&3519&3527&3524&3536&3536&3538&3517&3525\ (1.65$\%$)&5.0\\\ \chi_{c0}(1^{3}P_{0})&{3435}&3414&3320&3413&3415&3428&3433&3416&3392&3448&3488&3464&3445&3415\ (0.58$\%$)&5.0\\\ \chi_{c1}(1^{3}P_{1})&{3512}&3504&3490&3511&3521&3468&3510&3508&3491&3520&3514&3530&3510&3511\ (0.03$\%$)&5.0\\\ \chi_{c2}(1^{3}P_{2})&{3515}&3489&3550&3555&3553&3480&3556&3558&3570&3564&3565&3571&3550&3556\ (1.15$\%$)&5.0\\\ h_{c}(2^{1}P_{1})&{3962}&3824&3750&3926&3916&3943&3908&3960&3922&3950&3996&3933&3956&-&\\\ \chi_{c0}(2^{3}P_{0})&{3933}&3765&3830&3870&3848&3897&3842&3844&3845&3870&3947&3896&3916&3918\ (0.38$\%$)&5.0\\\ \chi_{c1}(2^{3}P_{1})&{4013}&3808&3670&3906&3914&3938&3901&3940&3902&3934&3972&3929&3953&-&\\\ \chi_{c2}(2^{3}P_{2})&{4010}&3915&-&3949&3937&3955&3937&3994&3949&3976&4021&3952&3979&3927\ (2.11$\%$)&5.0\\\ h_{c}(3^{1}P_{1})&{4269}&4137&-&4337&4139&4344&4184&-&4137&4291&-&4200&4318&-&\\\ \chi_{c0}(3^{3}P_{0})&{4240}&4080&-&4301&4146&4296&4131&-&4192&4214&-&4177&4292&-&\\\ \chi_{c1}(3^{3}P_{1})&{4320}&4121&3910&4319&4192&4338&4178&-&4178&4275&-&4197&4317&-&\\\ \chi_{c2}(3^{3}P_{2})&{4317}&4151&-&4354&4211&4358&4208&4212&4316&-&4213&4292&-&\\\ h_{c}(4^{1}P_{1})&{4385}&4416&-&4744&-&4704&-&-&-&-&-&4389&-&-&\\\ \chi_{c0}(4^{3}P_{0})&{4356}&4362&-&4698&-&4653&-&-&-&-&-&4374&-&-&\\\ \chi_{c1}(4^{3}P_{1})&{4436}&4401&-&4728&-&4696&-&-&-&-&-&4387&-&-&\\\ \chi_{c2}(4^{3}P_{2})&{4434}&4230&-&4763&-&4718&-&-&-&-&-&4398&-&-&\\\ \eta_{c2}(1^{1}D_{2})&{3799}&3478&3807&3807&3805&3765&3796&3824&3802&3803&3796&3848&3837&-&\\\ \psi({1^{3}D_{1}})&{3781}&3402&3739&3783&3792&3775&3787&3804&3729&3789&3792&3830&3819&3778\ (0.08$\%$)&5.0\\\ \psi_{2}({1^{3}D_{2}})&{3850}&3461&3550&3795&3807&3772&3798&3824&3788&3804&3794&3848&3838&3822\ (0.73$\%$)&5.0\\\ \psi_{3}({1^{3}D_{3}})&{3861}&3514&3869&3813&3808&3755&3799&3831&3844&3809&3798&3859&3849&-&\\\ \eta_{c2}(2^{1}D_{2})&{4141}&3828&-&4196&4108&4182&4099&4191&4105&4158&4224&4137&4208&-&\\\ \psi({2^{3}D_{1}})&{4122}&3756&-&4105&4095&4188&4089&4164&4057&4143&4222&4125&4194&4159\ (0.88$\%$)&5.0\\\ \psi_{2}({2^{3}D_{2}})&{4192}&3812&-&4190&4109&4188&4100&4189&4095&4159&4223&4137&4208&-&\\\ \psi_{3}({2^{3}D_{3}})&{4203}&3863&-&4220&4112&4176&4103&4202&4132&4167&4425&4144&4217&-&\\\ \eta_{c2}(3^{1}D_{2})&{4291}&4141&-&3549&4336&4553&4326&-&4330&-&-&4343&-&-&\\\ \psi({3^{3}D_{1}})&{4273}&4072&-&4507&4324&4555&4317&-&4293&-&-&4334&-&-&\\\ \psi_{2}({3^{3}D_{2}})&{4342}&4125&-&4544&4337&4557&4327&-&4322&-&-&4343&-&-&\\\ \psi_{3}({3^{3}D_{3}})&{4354}&4174&-&4574&4340&4549&4331&-&4351&-&-&4348&-&-&\\\ \hline\cr\chi^{2}/{N}&55.8&1250&&17.5&33.1&89.4&38.5&49.7&61.1&72.4&64.4&28.8&29.1&\\\ \hline\cr\hline\cr\hline\cr\hline\cr\end{array}$ Table 3: The mass spectra of bottomonium mesons, here number of experimental values ($N$=18), we list the error analysis is $\chi^{2}=$$\sum\limits_{i}$$(\frac{E^{th}_{i}-M^{exp}_{i}}{\mathcal{V}^{Er}_{i}})^{2}$. The parentheses in the reciprocal second column show the relative error between the theoretical and experimental values after fitting each experimental data of bottomonium in the model. All results are in units of MeV. State | $n^{2S+1}L_{J}$ | This work | NU Mansour and Gamal (2022) | Godfrey and Moats (2015) | Segovia  _et al._ (2016) | Bai-Qing and Kuang-Ta (2009) | JZ Wang  _et al._ (2018) | GI Barnes  _et al._ (2005) | PDG Zyla  _et al._ (2020) | Error
---|---|---|---|---|---|---|---|---|---|---
$\eta_{b}(1S)$ | $1^{1}S_{0}$ | 9451 | 9436 | 9402 | 9455 | 9389 | 9398 | 9394 | 9399 (0.55$\%$) | 5.0
$\eta_{b}(2S)$ | $2^{1}S_{0}$ | 9951 | 9991 | 9976 | 9990 | 9987 | 9989 | 9975 | 9999 (0.48$\%$) | 5.0
$\eta_{b}(3S)$ | $3^{1}S_{0}$ | 10354 | 10139 | 10336 | 10330 | 10330 | 10336 | 10333 | – | 5.0
$\eta_{b}(4S)$ | $4^{1}S_{0}$ | 10658 | 10324 | 10523 | – | 10595 | 10597 | 10616 | – | 5.0
$\eta_{b}(5S)$ | $5^{1}S_{0}$ | 10865 | 10498 | 10869 | – | 10817 | 10597 | 10860 | – | 5.0
$\eta_{b}(6S)$ | $6^{1}S_{0}$ | 10973 | 10662 | 11097 | – | 11011 | 10991 | 11079 | – | 5.0
$\Upsilon(1S)$ | $1^{3}S_{1}$ | 9458 | 9491 | 9465 | 9502 | 9460 | 9463 | 9459 | 9460 (0.02$\%$) | 5.0
$\Upsilon(2S)$ | $2^{3}S_{1}$ | 9958 | 10013 | 10003 | 10015 | 10016 | 10017 | 10004 | 10023 (0.65$\%$) | 5.0
$\Upsilon(3S)$ | $3^{3}S_{1}$ | 10361 | 10328 | 10354 | 10349 | 10351 | 10356 | 10354 | 10355 (0.06$\%$) | 5.0
$\Upsilon(4S)$ | $4^{3}S_{1}$ | 10665 | 10546 | 10635 | 10607 | 10611 | 10612 | 10633 | 10579 (0.81$\%$) | 5.0
$\Upsilon(10860)$ | $5^{3}S_{1}$ | 10872 | 10826 | 10878 | 10818 | 10831 | 10822 | 10875 | 10882 (0.09$\%$) | 5.0
$\Upsilon(11020)$ | $6^{3}S_{1}$ | 10980 | 10971 | 11102 | 10995 | 10988 | 11001 | 11092 | 11003 (0.2$\%$) | 5.0
$\Upsilon(7S)$ | $7^{3}S_{1}$ | 10991 | – | – | – | – | 11157 | 11294 | – | 5.0
$h_{b}(1P)$ | $1^{1}P_{1}$ | 9893 | 9879 | 9882 | 9879 | 9903 | 9894 | 9881 | 9899 (0.06$\%$) | 5.0
$h_{b}(2P)$ | $2^{1}P_{1}$ | 10254 | 10222 | 10250 | 10240 | 10256 | 10259 | 10250 | 10260 (0.06$\%$) | 5.0
$h_{b}(3P)$ | $3^{1}P_{1}$ | 10516 | 10417 | 10541 | 10516 | 10529 | 10530 | 10540 | – |
$h_{b}(4P)$ | $4^{1}P_{1}$ | 10681 | 10324 | 10790 | – | 10757 | 10751 | 10790 | – |
$h_{b}(5P)$ | $5^{1}P_{1}$ | 10748 | 10498 | 11016 | – | 10955 | 10938 | 11013 | – |
$\chi_{b0}(1P)$ | $1^{3}P_{0}$ | 9867 | 9843 | 9847 | 9855 | 9865 | 9858 | 9845 | 9859 (0.08$\%$) | 5.0
$\chi_{b0}(2P)$ | $2^{3}P_{0}$ | 10227 | 10196 | 10226 | 10221 | 10226 | 10235 | 10225 | 10231 (0.04$\%$) | 5.0
$\chi_{b0}(3P)$ | $3^{3}P_{0}$ | 10490 | 10134 | 10552 | 10500 | 10502 | 10513 | 10521 | – |
$\chi_{b0}(4P)$ | $4^{3}P_{0}$ | 10655 | 10319 | 10775 | – | 10732 | 10736 | 10773 | – |
$\chi_{b0}(5P)$ | $5^{3}P_{0}$ | 10722 | 10494 | 11004 | – | 10933 | 10926 | 10998 | – |
$\chi_{b1}(1P)$ | $1^{3}P_{1}$ | 9893 | 9874 | 9876 | 9874 | 9897 | 9889 | 9875 | 9893 (0$\%$) | 5.0
$\chi_{b1}(2P)$ | $2^{3}P_{1}$ | 10254 | 10217 | 10246 | 10236 | 10251 | 10255 | 10246 | 10256 (0.01$\%$) | 5.0
$\chi_{b1}(3P)$ | $3^{3}P_{1}$ | 10516 | 10138 | 10538 | 10513 | 10524 | 10527 | 10537 | 10512 (0.04$\%$) | 5.0
$\chi_{b2}(1P)$ | $1^{3}P_{2}$ | 9914 | 9891 | 9897 | 9886 | 9918 | 9910 | 9896 | 9912 (0.01$\%$) | 5.0
$\chi_{b2}(2P)$ | $2^{3}P_{2}$ | 10274 | 10230 | 10261 | 10246 | 10269 | 10269 | 10261 | 10269 (0.05$\%$) | 5.0
$\chi_{b2}(3P)$ | $3^{3}P_{2}$ | 10537 | 10141 | 10550 | 10521 | 10540 | 10539 | 10549 | – |
$\chi_{b2}(4P)$ | $4^{3}P_{2}$ | 10702 | 10325 | 10798 | – | 10767 | 10758 | 10797 | – |
$\chi_{b2}(5P)$ | $5^{3}P_{2}$ | 10767 | – | – | – | – | 10944 | 11020 | – |
$\Upsilon(1D)$ | $1^{3}D_{2}$ | 10164 | 10113 | 10147 | 10122 | 10151 | 10162 | 10147 | 10164 (0$\%$) | 5.0
$\Upsilon_{2}(2D)$ | $2^{3}D_{2}$ | 10245 | 9943 | 10449 | 10418 | 10438 | 10450 | 10449 | – |
$\Upsilon_{1}(1D)$ | $1^{3}D_{1}$ | 10198 | 9729 | 10138 | 10117 | 10145 | 10153 | 10137 | – |
$\Upsilon_{1}(2D)$ | $2^{3}D_{1}$ | 10210 | 9938 | 10441 | 10414 | 10432 | 10442 | 10441 | – |
$\Upsilon_{3}(1D)$ | $1^{3}D_{3}$ | 10198 | 9739 | 10115 | 10127 | 10156 | 10170 | 10155 | – |
$\Upsilon_{3}(2D)$ | $2^{3}D_{3}$ | 10279 | 9947 | 10455 | 10422 | 10442 | 10456 | 10455 | – |
$\chi^{2}/{N}$ |  | $38.8$ | 351.6 | 35.7 | $32.4$ | 10 | 11.3 | 31.4 |  |
Figure 1: Fitted curves of the $\psi$, $\chi_{c1}$, $h_{b}$, $\chi_{b2}$
radial $(M,n)$ resonances calculated by the modified rovibrational model, here
the radial quantum number is $n=n_{r}+1$, solid-line curve show clearly the
mass variation trend when different radial quantum number $n$ is taken. The
open circle and the filled geometry are the theoretical and experimental
values, respectively.

####  III.0.2 The spectra of bottomonium meson excitations

Table 3 show that have been established theoretical and experimental values.
The results of the thirteenth column are calculated using the modified GI
relativistic quark model with screening effects in bottomonium system Wang
_et al._ (2018), all experimental masses are taken from PDG averages Zyla  _et
al._ (2020). Although GI model successfully studied the bottomonium spectrum,
the qualitative evaluation errors selection of all states are
$\mathcal{V}^{Er}_{i}$=5.0 MeV, comparing $\chi^{2}/{N}$ values, one can
easily see that the fitted modified rovibrational model improves the whole
description of the bottomonium spectrum. Since the screening potential can
lower the energy of the high excited states. In the words, the modified GI
model well improves the accuracy of the whole bottomonium system. The
calculation results of other models are listed for comparison in the
literatures Bai-Qing and Kuang-Ta (2009); Barnes  _et al._ (2005); Wang  _et
al._ (2018); Mansour and Gamal (2022); Godfrey and Moats (2015); Segovia  _et
al._ (2016).

##  IV conclusion

This paper presents analysis of mass spectra of the radial and orbital
resonances of charmonium and bottomonium system, using the modified
rovibrational model. Rovibrational model study on diatomic molecular spectra
show satisfactory agreement with experimental values for many molecules. Until
now, the rovibrational model is useful for classifying the hadron resonances,
especially for charmonium and bottomonium system with double heavy quarks. The
modified Godfrey-Isgur quark model and Cornell potential model with the
screening effect traditionally include spin-spin interaction and orbital-spin
interaction. The rovibrational model of diatomic molecular mass spectra
includes the vibration and rotation effects of two atoms, so it is spin-
independent. For small scale of hadrons, spin and orbital dependence must be
considered. According to the comparison of the fitted curves of the $\psi$,
$\chi_{c1}$, $h_{b}$, $\chi_{b2}$ radial $(M,n)$ resonances calculated by the
modified rovibrational, the calculation of too high excited state energy
spectrum has some limitations. The reason is that the interaction of morse
potential decreases at high excited states. It is noteworthy that many
hadronic molecules are mainly ground states, and the modified rovibrational
model has certain reference significance.

##  V ACKNOWLEDGMENTS

This work is supported by the National Natural Science Foundation of China
under Grants No. 11965016, the projects funded by Science and Technology
Department of Qinghai Province (No. 2020-ZJ-728).

## References

  * Aubert  _et al._ (1974) J.-J. Aubert, U. Becker, P. Biggs, J. Burger, M. Chen, G. Everhart, P. Goldhagen, J. Leong, T. McCorriston, T. Rhoades, _et al._ , Physical Review Letters 33, 1404 (1974).
  * Herb  _et al._ (1977) S. Herb, D. Hom, L. Lederman, J. Sens, H. Snyder, J. Yoh, J. Appel, B. Brown, C. Brown, W. R. Innes, _et al._ , Physical Review Letters 39, 252 (1977).
  * Innes  _et al._ (1977) W. R. Innes, J. Appel, B. Brown, C. Brown, K. Ueno, T. Yamanouchi, S. Herb, D. Hom, L. Lederman, J. Sens, _et al._ , Physical Review Letters 39, 1240 (1977).
  * Fischer  _et al._ (2015) C. S. Fischer, S. Kubrak, and R. Williams, The European Physical Journal A 51, 1 (2015).
  * Lakhina and Swanson (2006) O. Lakhina and E. S. Swanson, Phys. Rev. D 74, 014012 (2006).
  * Radford and Repko (2007) S. F. Radford and W. W. Repko, Phys. Rev. D 75, 074031 (2007).
  * Soni  _et al._ (2018) N. Soni, B. Joshi, R. Shah, H. Chauhan, and J. Pandya, The European Physical Journal C 78, 1 (2018).
  * Ebert  _et al._ (2011) D. Ebert, R. N. Faustov, and V. O. Galkin, The European Physical Journal C 71, 1 (2011).
  * Deng  _et al._ (2017) W.-J. Deng, H. Liu, L.-C. Gui, and X.-H. Zhong, Phys. Rev. D 95, 034026 (2017).
  * Li and Chao (2009) B.-Q. Li and K.-T. Chao, Phys. Rev. D 79, 094004 (2009).
  * Shah  _et al._ (2012a) M. Shah, A. Parmar, and P. C. Vinodkumar, Phys. Rev. D 86, 034015 (2012a).
  * Bai-Qing and Kuang-Ta (2009) L. Bai-Qing and C. Kuang-Ta, Communications in Theoretical Physics 52, 653 (2009).
  * Barnes  _et al._ (2005) T. Barnes, S. Godfrey, and E. Swanson, Physical Review D 72, 054026 (2005).
  * Wang  _et al._ (2019) J.-Z. Wang, D.-Y. Chen, X. Liu, and T. Matsuki, Physical Review D 99, 114003 (2019).
  * Wang  _et al._ (2018) J.-Z. Wang, Z.-F. Sun, X. Liu, and T. Matsuki, The European Physical Journal C 78, 1 (2018).
  * Ebert  _et al._ (2013) D. Ebert, R. Faustov, and V. Galkin, Physics of Atomic Nuclei 76, 1554 (2013).
  * Chen (2018) J.-K. Chen, The European Physical Journal C 78, 1 (2018).
  * Chen (2021) J.-K. Chen, The European Physical Journal A 57, 1 (2021).
  * Albaladejo  _et al._ (2021) M. Albaladejo, L. Bibrzycki, S. M. Dawid, C. Fernandez-Ramirez, S. Gonzalez-Solis, A. N. H. Blin, A. W. Jackura, V. Mathieu, M. Mikhasenko, V. I. Mokeev, _et al._ , arXiv preprint arXiv:2112.13436 (2021).
  * Bernardo  _et al._ (2021) D. L. Bernardo, C. C. Bastos, and A. C. Pavão, Chinese Physics C 45, 084104 (2021).
  * Pauling (1966) L. Pauling, Proceedings of the National Academy of Sciences of the United States of America 56, 1676 (1966).
  * Morse (1929) P. M. Morse, Physical review 34, 57 (1929).
  * Okon  _et al._ (2017) I. B. Okon, O. Popoola, and C. N. Isonguyo, Advances in High Energy Physics 2017 (2017).
  * Mansour and Gamal (2022) H. Mansour and A. Gamal, Results in Physics , 105203 (2022).
  * Godfrey and Isgur (1985) S. Godfrey and N. Isgur, Physical Review D 32, 189 (1985).
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Shah  _et al._ (2012b) M. Shah, A. Parmar, and P. Vinodkumar, Physical Review D 86, 034015 (2012b).
  * Patel  _et al._ (2016) S. Patel, P. C. Vinodkumar, and S. Bhatnagar, Chinese Physics C 40, 053102 (2016).
  * Godfrey and Moats (2015) S. Godfrey and K. Moats, Phys. Rev. D 92, 054034 (2015).
  * Segovia  _et al._ (2016) J. Segovia, P. G. Ortega, D. R. Entem, and F. Fernández, Phys. Rev. D 93, 074027 (2016).