# Triaxial-shape dynamics in the low-lying excited $0^{+}$ state: Role of the
collective mass

Kouhei Washiyama<EMAIL_ADDRESS>Center for Computational
Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan  Research
Center for Superheavy Elements, Kyushu University, Fukuoka 819-0395, Japan
Kenichi Yoshida<EMAIL_ADDRESS>Research Center for Nuclear
Physics, Osaka University, Ibaraki, Osaka 567-0047 Japan  RIKEN Nishina Center
for Accelerator-Based Science, Wako, Saitama 351-0198, Japan

###### Abstract

Background

Non-yrast states in neutron-rich nuclei are being investigated experimentally.
These states reveal various aspects and details of the nuclear structure, such
as the fluctuation around the axially symmetric shape.

Purpose

The beyond-mean-field effects in neutron-rich nuclei with $N\simeq 28$ are
investigated. We focus on the role of collective mass in triaxial-shape
dynamics.

Method

We employ the five-dimensional quadrupole collective Hamiltonian method with
the potential obtained in a constrained Hartree–Fock–Bogoliubov approach with
a Skyrme energy-density functional and the collective-mass functions obtained
by the cranking approximation. The method includes triaxial deformations.

Results

We find that 42Mg, 40Si, 44S, and 46S show $\gamma$-soft: A flat behavior in
the potential energy surface along the triaxial deformation. Their low-lying
spectra show a strong nucleus dependence, while those obtained with a
collective mass assumed as constant are similar to each other. The energy
ratio $E(0_{2}^{+})/E(2_{1}^{+})$ and the $B(E2)$ ratio $B(E2;0_{2}^{+}\to
2_{1}^{+})/B(E2;2_{1}^{+}\to 0_{1}^{+})$ show a unique property of the
$0_{2}^{+}$ state, while the energy and $B(E2)$ ratios in neutron-deficient
$\gamma$-soft nuclei with $N=78$ do not depend on nucleus so much.

Conclusions

Low-lying spectra are determined by not only the potential energy but also the
collective mass. We clarify the important role of the collective mass in low-
energy dynamics in the neutron-rich $N\sim 28$ nuclei.

##  I Introduction

Atomic nuclei exhibit various shapes according to the neutron number, the
proton number, and the excitation energy. Since the ingredients are finite and
their orbital motion is described by quantum mechanics, it is essential to
consider the fluctuation in shape. The five-dimensional quadrupole collective
Hamiltonian (5DCH) [1, 2, 3, 4, 5] as a function of the quadrupole deformation
parameters $\beta$ and $\gamma$ has often been employed in describing low-
energy states. The parameters in the model are the potential energy and the
collective masses.

The shape dynamics is governed by the potential energy at first glance. In
nuclei near the magic numbers, the potential energy surface (PES) shows the
existence of the local minimum at the spherical configuration. As the system
moves away from the magic numbers, it becomes deformed, where one mostly has
the local minimum at the prolately-deformed configuration. Some nuclei show a
soft PES against the triaxial deformation. Consequently, the low-lying
$2^{+}_{2}$ state or the so-called $\gamma$ vibration shows up.

The vibration in the triaxial deformation is not always harmonic. An ideal
situation where the PES against the $\gamma$ direction is flat is investigated
by the Wilets–Jean model [6]. In this model, the mass parameters are assumed
to be constant. A characteristic feature of the low-lying states is that the
$2_{2}^{+}$ state degenerates with the $4^{+}_{1}$ state and is lower than the
$0_{2}^{+}$ state among two-phonon states in a spherical harmonic-oscillator
potential and that the $0_{2}^{+}$ is degenerate with the $3_{1}^{+}$,
$4_{2}^{+}$, and $6_{1}^{+}$ states. In addition to the flatness in the
$\gamma$ direction, another ideal situation where the PES against the $\beta$
direction is flat is investigated in terms of the E(5) critical point symmetry
[7, 8, 9, 10]. The PES in the $\beta$ direction is described by an infinite
square-well potential. The degeneracy of the $2_{2}^{+}$ and $4_{1}^{+}$ is
the same as in the Wilets–Jean model, while the $0_{2}^{+}$ state is not
necessarily degenerate with the $3_{1}^{+}$, $4_{2}^{+}$, and $6_{1}^{+}$
states because of the fluctuation in the $\beta$ direction.

Neutron-rich nuclei around $N=28$ have attracted interest both experimentally
[11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] and theoretically [22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33]. The authors in Refs. [23, 32, 28] found
that the breaking of $N=28$ magicity in 44S is due to a flat potential, which
brings about a wide configuration mixing in the $\beta$–$\gamma$ deformation
space. In a neighboring nucleus 43S, the coexistence of prolately-,
triaxially-, and oblately-deformed states was predicted due to the breaking of
$N=28$ magicity in Ref. [26]. In $N=26$ and $30$ isotones, it was found that
including the triaxial degree of freedom lowers the energy of the $2_{2}^{+}$
state [27]. The shape coexistence of prolate and oblate configurations and
that of oblate and spherical configurations were predicted in 40Mg and 42Si,
respectively [28]. The $\gamma$ vibration was predicted to appear for the
prolate configuration in 40Mg [33]. Those studies have shown that the triaxial
deformation plays an important role in low-energy dynamics in $N\sim 28$
nuclei.

Not only the PES but also the mass parameters may play a role in describing
the low-energy dynamics in a collective Hamiltonian approach. The mass
parameters, as well as the potential energy, depend on the deformations. The
terrain of the PES and the deformation-dependence of mass parameters are
sensitively determined by the shell effect. Thus by changing the
neutron/proton number involving neutron-rich and neutron-deficient nuclei, one
can investigate the role of the mass parameters in low-energy dynamics.
Therefore, we study in the present work the low-energy dynamics governed by
triaxial deformation in neutron-rich nuclei around $N=28$, putting a focus on
the role of the mass parameters in the collective Hamiltonian.

The paper is organized as follows. In Sec. II, we briefly explain the 5DCH
method. In Sec. III, we show the results and discuss the roles of the mass
parameters. Section IV summarizes the paper.

##  II Method

We briefly explain the five-dimensional quadrupole collective Hamiltonian
method. For details, we refer to Refs. [1, 2, 3, 4, 5].

The collective Hamiltonian reads

$H=T_{\text{vib}}+T_{\text{rot}}+V(\beta,\gamma),$ (1)

with the vibrational and rotational kinetic energies,

$\displaystyle T_{\text{vib}}$
$\displaystyle=\frac{1}{2}D_{\beta\beta}(\beta,\gamma)\dot{\beta}^{2}+D_{\beta\gamma}(\beta,\gamma)\dot{\beta}\dot{\gamma}+\frac{1}{2}D_{\gamma\gamma}(\beta,\gamma)\dot{\gamma}^{2},$
(2) $\displaystyle T_{\text{rot}}$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{3}\mathcal{J}_{k}(\beta,\gamma)\omega_{k}^{2}.$
(3)

The functions $D_{\beta\beta}(\beta,\gamma)$, $D_{\beta\gamma}(\beta,\gamma)$,
and $D_{\gamma\gamma}(\beta,\gamma)$ denote the vibrational masses and
$\mathcal{J}_{k}(\beta,\gamma)=4\beta^{2}D_{k}(\beta,\gamma)\sin^{2}(\gamma-2\pi
k/3)$ and $\omega_{k}$ denote the rotational moments of inertia and rotational
angular velocities in the body-fixed frame of a nucleus. After quantizing the
Hamiltonian, we obtain the collective Schrödinger equation as

$\displaystyle[\hat{T}_{\text{vib}}+\hat{T}_{\text{rot}}+V(\beta,\gamma)]\Psi_{\alpha
IM}(\beta,\gamma,\Omega)=E_{\alpha I}\Psi_{\alpha IM}(\beta,\gamma,\Omega),$
(4)

where $E_{\alpha I}$ and $\Psi_{\alpha IM}(\beta,\gamma,\Omega)$ are the
excitation energies and the collective wave functions with the total angular
momentum $I$, its $z$ component $M$ in the laboratory frame, and $\alpha$
distinguishing the states with the same $I$ and $M$. The collective wave
functions are functions of $\beta$, $\gamma$, and the three Euler angles
$\Omega$ and are written as

$\displaystyle\Psi_{\alpha
IM}(\beta,\gamma,\Omega)=\sum_{K=\text{even}}^{I}\Phi_{\alpha
IK}(\beta,\gamma)\braket{\Omega}{IMK},$ (5)

with $\braket{\Omega}{IMK}$ being linear combinations of the Wigner rotational
wave functions and $K$ being the $z$ component of $I$ in the body-fixed frame.
The vibrational wave functions $\Phi_{\alpha IK}(\beta,\gamma)$ are normalized
as

$\int_{0}^{\infty}d\beta\int_{0}^{\pi/3}d\gamma|G(\beta,\gamma)|^{1/2}\\\
\times\sum_{K=\text{even}}^{I}\Phi^{*}_{\alpha
IK}(\beta,\gamma)\Phi_{\alpha^{\prime}IK}(\beta,\gamma)=\delta_{\alpha\alpha^{\prime}},$
(6)

where the volume element $|G(\beta,\gamma)|^{1/2}$ is given by

$|G(\beta,\gamma)|^{1/2}=2\beta^{4}\sqrt{W(\beta,\gamma)R(\beta,\gamma)}\sin
3\gamma,$ (7)

with
$W(\beta,\gamma)=\\{D_{\beta\beta}(\beta,\gamma)D_{\gamma\gamma}(\beta,\gamma)-[D_{\beta\gamma}(\beta,\gamma)]^{2}\\}\beta^{-2}$
and
$R(\beta,\gamma)=D_{1}(\beta,\gamma)D_{2}(\beta,\gamma)D_{3}(\beta,\gamma)$.
Using the collective wave functions (5), the reduced quadrupole transition
probability is given by

$B(E2;\alpha I\to\alpha^{\prime}I^{\prime})=(2I+1)^{-1}|\braket{\alpha
I}{\hat{\mathcal{M}}(E2)}{\alpha^{\prime}I^{\prime}}|^{2},$ (8)

where $\hat{\mathcal{M}}(E2)$ is the electric quadrupole operator.

The collective potential $V(\beta,\gamma)$ is obtained by solving the
constrained Hartree–Fock–Bogoliubov (CHFB) equation constructed from a Skyrme
energy density functional (EDF). The vibrational masses and the rotational
moments of inertia are calculated by the so-called cranking approximation [34,
35, 36] with the quasiparticle states obtained by the CHFB equation, hereafter
denoted as the cranking mass.

We solve the CHFB equation with the two-basis method [37, 38] in the three-
dimensional Cartesian mesh of $R=12.4$ fm with the mesh size of $0.8$ fm. The
single-particle states at positive energies are cut off so as to become the
equivalent quasiparticle energy of $E_{\text{QP}}\approx 30$ MeV, which gives
a good convergence in both the collective potential and the cranking mass in
this paper. We used the SkM∗ EDF [39] and the pairing EDF proposed in Ref.
[40], which depends on the isoscalar and isovector densities. Note that, since
we used a different pairing cutoff scheme [41, 42] from the one in Ref. [40],
we refitted the pairing strength to reproduce the neutron pairing gap of
156Dy. For the collective Hamiltonian in the $\beta$–$\gamma$ plane, we employ
a triangular mesh of $\Delta\beta\approx 0.035$ in the region $0<\beta<0.6$
and $0^{\circ}<\gamma<60^{\circ}$, consisting of about 200 mesh points.

##  III Results and discussion

Figure 1:  Potential energy surface in the $\beta$–$\gamma$ plane of Mg, Si,
and S isotopes with $N=26$, 28, and 30 obtained from constrained DFT
calculations with the SkM∗ EDF.  Figure 2:  Potential energy as a function of
$\gamma$ at a fixed $\beta$ for the selected nuclei. The potential is shifted
to $V(\gamma=0^{\circ})=0$.

Figure 1 shows the PESs in the $\beta$–$\gamma$ plane for the Mg, Si, and S
isotopes with $N=26$, 28, and 30 calculated by the constrained HFB. In the
obtained PESs, one sees the minimum at the prolate configuration in 38Mg and
42S and at the oblate configuration in 42,44Si. The PES of 40Mg has two local
minima at the prolate and oblate sides. One can see that 42Mg, 40Si, 44S, and
46S possess a unique property that the potential shows a flat behavior from
the spherical point to a certain value of $\beta$, and along the $\gamma$
direction, namely $\gamma$ soft.

To see clearly a $\gamma$-soft behavior of PESs, Fig. 2 shows the potential as
a function of $\gamma$ at a fixed $\beta$ for the nuclei mentioned above. In
this figure, the potential is shifted to become $V(\gamma=0^{\circ})=0$. For
each nucleus, the value of $\beta$ is chosen as its mean value obtained with
the ground-state collective wave function ($\alpha=1$, $I=K=0$),

$\braket{\beta}=\int d\beta
d\gamma|G(\beta,\gamma)|^{1/2}|\Phi_{100}(\beta,\gamma)|^{2}\beta,$ (9)

and indicated in the figure. The change in potential along the $\gamma$
direction is 0.2–0.4 MeV.

Figure 3:  Low-lying excitation spectra and $B(E2)$ values in units of $e^{2}$
fm4 of 44S (top) and 46S (middle) with the cranking mass (left) and the
constant mass (right). Those of the E(5)–$\beta^{4}$ and Wilets–Jean models
are shown at the bottom. In the E(5)–$\beta^{4}$ and Wilets–Jean models, the
value of $D$ in the constant mass is determined by fitting $E(0_{2}^{+})=4$
MeV and the $B(E2)$ values are normalized to $B(E2;2_{1}^{+}\to 0_{1}^{+})=50$
$e^{2}$ fm4.

In what follows, we investigate the low-lying states unique in $\gamma$-soft
nuclei. The left panels of Fig 3 shows low-lying spectra of 44S (top) and 46S
(middle) obtained with the cranking mass. In the spectra, $I\leq 6$ in the
$0_{1}^{+}$ band and $I\leq 4$ in the $0_{2}^{+}$ and $2_{2}^{+}$ bands are
plotted and the $B(E2)$ values larger than 1 W.u. are shown. We find that the
spectra between 44S and 46S are very different, especially the order of
$4_{1}^{+}$, $2_{2}^{+}$, and $0_{2}^{+}$, even though their PESs are very
similar to each other as shown in Fig. 1. The energy ratios defined as
$R_{0/2}=E(0_{2}^{+})/E(2_{1}^{+})$, $R_{2/2}=E(2_{2}^{+})/E(2_{1}^{+})$, and
$R_{4/2}=E(4_{1}^{+})/E(2_{1}^{+})$ are $R_{0/2}=1.65$ and 2.45,
$R_{2/2}=1.86$ and 1.89, and $R_{4/2}=2.34$ and 2.49 for 44S and 46S,
respectively. To understand the reason for this difference, we focus on the
role of the vibrational and rotational masses on excitation spectra. To this
end, we employ the Bohr Hamiltonian with a constant mass [43], where the
vibrational and rotational masses obtained by the cranking approximation is
replaced by a constant value, neglecting the $\beta$–$\gamma$ dependence.
Namely [43],

$D_{\beta\beta}=D_{\gamma\gamma}/\beta^{2}=D_{1}=D_{2}=D_{3}\equiv D,\quad
D_{\beta\gamma}=0.$ (10)

The right panels in Fig. 3 show the spectra of 44S and 46S obtained with the
constant mass. The value of $D$ is determined to reproduce the $0_{2}^{+}$
energy obtained with the cranking mass. Clearly, the patterns of the spectra
obtained with the constant mass are similar to each other. The energy ratios
are $R_{0/2}=2.24$ and 2.24, $R_{2/2}=2.07$ and 2.08, and $R_{4/2}=2.06$ and
2.06 for 44S and 46S, respectively. Note that how to determine the value of
$D$ does not change these energy ratios, though the absolute value of energy
changes. If $D$ is determined to give the $2_{2}^{+}$ energy with the cranking
mass, we obtain $R_{0/2}=2.23$, $R_{2/2}=2.07$, and $R_{4/2}=2.06$ for 44S.

Figure 4:  Same as Fig. 1, but for the neutron-deficient $N=78$ 134Ba, 136Ce,
138Nd, 140Sm, and 142Gd nuclides.  Figure 5:  Same as Fig. 2, but for the
$N=78$ nuclei.

To see whether the low-lying spectra are sensitive to the collective mass in
other mass regions, we study low-energy dynamics of neutron-deficient $N=78$
nuclei 134Ba, 136Ce, 138Nd, 140Sm, and 142Gd, whose PESs show flat in the
$\beta$ direction and $\gamma$ soft similarly to those of 44S and 46S. Figures
4 and 5 show how PESs of those nuclei are $\gamma$ soft. The change in PESs is
0.2–0.6 MeV. As the proton number decreases, the PES becomes flatter.

Figure 6:  $R_{0/2}$ (filled square), $R_{2/2}$ (filled circle) and $R_{4/2}$
(filled triangle) for the selected nuclei with the cranking mass (a) and
constant mass (b). The ratios obtained with the E(5)–$\beta^{4}$ and
E(5)–$\beta^{6}$ models and the Wilets–Jean model (WJ) are included in (b).
Note that the $R_{0/2}$ value for the Wilets–Jean model is 3.9. The panel (c)
shows the ratios of the available experimental data from [14] for 44S, from
[20] for 40Si, and from [44] for the $N=78$ nuclei.

Figure 6 shows the $R_{4/2}$, $R_{2/2}$, and $R_{0/2}$ values of 42Mg, 40Si,
44S, and 46S and those of the selected $N=78$ nuclei. The top panel shows the
ratios obtained with the cranking mass, while the bottom depicts those with
the constant mass. In the light nuclei with the cranking mass, $R_{0/2}$
strongly depends on the nucleus, while the variation of $R_{4/2}$ and
$R_{2/2}$ is small. The ratios in $N=78$ nuclei are almost constant,
$R_{2/2}\sim 2.0$ and $R_{0/2}\sim R_{4/2}\sim 2.2$–2.3 by changing nucleus.
Compared with the cranking mass case, the variation of $R_{0/2}$ in the light
nuclei becomes significantly small in the constant mass case. In the $N=78$
nuclei, $R_{4/2}$, $R_{2/2}$, and $R_{0/2}$ values are almost constant around
2.1, 2.1, and 2.5, respectively.

Properties of the low-energy states in $\gamma$-soft nuclei have often been
discussed in view of the E(5) critical point symmetry [7, 8, 9]. The E(5)
symmetry is realized in the Bohr Hamiltonian with an infinite square-well
potential in $\beta$ and a constant in $\gamma$. The $\beta^{2n}$ form of the
potential instead of the infinite-well potential is introduced to describe
realistic systems [10, 9]. In Fig. 6, the energy ratios in the
E(5)–$\beta^{4}$ and E(5)–$\beta^{6}$ models are included and denoted as
$\beta^{4}$ and $\beta^{6}$, respectively. The energy ratios in the light
$N\sim 28$ nuclei and $N=78$ nuclei with the constant mass are close to those
in the E(5)–$\beta^{4}$ and E(5)–$\beta^{6}$ models. The Wilets–Jean model,
which describes $\gamma$-flat PES, gives higher energy ratios than the ones
with the constant mass and with the E(5)–$\beta^{4}$ and E(5)–$\beta^{6}$
models. In the bottom panels of Fig. 3, the low-lying spectra of the
E(5)–$\beta^{4}$ and the Wilets–Jean models are shown up to seniority three.
These spectra are obtained by $V(\beta,\gamma)=C\beta^{4}$ with $C=80$ MeV for
E(5)–$\beta^{4}$ and by $V(\beta,\gamma)=C(\beta-0.3)^{2}$ with $C=200$ MeV
for the Wilets–Jean model, and with the constant mass $D$ giving
$E(0_{2}^{+})=4$ MeV. The pattern of low-lying spectra in the E(5)–$\beta^{4}$
is close to those of 44S and 46S with the constant mass: The degeneracy of the
$2_{2}^{+}$ and $4_{1}^{+}$ states and that of the $0_{3}^{+}$, $3_{1}^{+}$,
$4_{2}^{+}$, and $6_{1}^{+}$ states.

The energy ratios obtained from the available experimental data in neutron-
rich $N\sim 28$ nuclei are $R_{4/2}=1.86$, $R_{2/2}=1.63$ and $R_{0/2}=1.03$
in 44S [14], and $R_{4/2}=2.56$ in 40Si [20], as shown in Fig. 6(c). The
energy ratio $R_{0/2}$ with the cranking mass shows a strong nucleus
dependence in neutron-rich $N\sim 28$ nuclei, and one sees a sudden drop in
44S although the measured value is even lower than the calculated one. In the
$N=78$ nuclei, the energy ratios of the available experimental data in Fig.
6(c) are almost constant, though some of the $R_{0/2}$ values are not
available. These energy ratios are similar to those with both the cranking and
constant mass cases in $N=78$. More experimental data on the low-lying $0^{+}$
state will give insight into roles of the collective mass.

Figure 7:  Decomposition of the $0_{2}^{+}$ energy to the vibrational kinetic
energy $T_{\text{vib}}$ and the potential energy $V$ divided by the total
energy with the cranking mass for the selected nuclei.

To further investigate the origin of the different behaviors of $R_{0/2}$ in
the $N\sim 28$ and $N=78$ nuclei, we decompose the energy of the $0_{2}^{+}$
state to the vibrational kinetic energy and potential energy. These energies
are calculated with the vibrational wave function of the $0_{2}^{+}$ state as

$E=\int d\beta
d\gamma|G(\beta,\gamma)|^{1/2}\Phi^{*}_{200}(\beta,\gamma)\hat{E}\Phi_{200}(\beta,\gamma),$
(11)

where $\hat{E}=\hat{T}_{\text{vib}}$ or $V(\beta,\gamma)$. Figure 7 shows the
vibrational kinetic energy $T_{\text{vib}}$ and the potential energy $V$
divided by their sum $E_{\text{total}}=\braket{T_{\text{vib}}}+\braket{V}$,
the absolute value of the $0_{2}^{+}$ energy, for the $N\sim 28$ and $N=78$
nuclei. Note that there is no contribution of rotation in the $I=0$ states. In
the $N\sim 28$ nuclei, $T_{\text{vib}}$ is larger than $V$. A large
$T_{\text{vib}}$ with the deformation-dependent cranking mass enhances the
nucleus dependence of the $0_{2}^{+}$ energy. Furthermore, we found that the
ratios of $T_{\text{vib}}$, $T_{\text{rot}}$, and $V$ to the total energy in
the $0_{1}^{+}$ and $2_{1}^{+}$ states are almost the same and do not depend
on the nucleus so much. Therefore, we observe the strong nucleus-dependent
$R_{0/2}$ values. In $N=78$, the two energies are close to each other and do
not much depend on the nucleus. Thus, we have obtained that the $R_{0/2}$
values do not depend on the nucleus even with the deformation-dependent
cranking mass.

Figure 8:  $B(E2)$ ratio for the selected nuclei with the cranking mass (a)
and constant mass (b) together with E(5)–$\beta^{4}$ and E(5)–$\beta^{6}$
models and the Wilets–Jean model (WJ).

Figure 9:  Vibrational wavefunctions $|\Psi_{\alpha 00}(\beta,\gamma)|^{2}$
multiplied by $\beta^{4}|W(\beta,\gamma)R(\beta,\gamma)|^{1/2}$ of the
$0_{1}^{+}$ (top) and $0_{2}^{+}$ (bottom) states with the cranking and
constant masses for 44S and 46S. The right two columns show those with the
E(5)–$\beta^{4}$ and Wilets–Jean models.

The similar feature to the energy ratios is indeed obtained in the $B(E2)$
ratios. Figure 8 shows the $B(E2)$ ratios $RB_{0/2}=B(E2;0_{2}\to
2_{1})/B(E2;2_{1}\to 0_{1})$, $RB_{2/2}=B(E2;2_{2}\to 2_{1})/B(E2;2_{1}\to
0_{1})$, and $RB_{4/2}=B(E2;4_{1}\to 2_{1})/B(E2;2_{1}\to 0_{1})$ for the
selected nuclei with the cranking mass (top) and the constant mass (bottom).
The $RB_{0/2}$ value in the light $N\sim 28$ nuclei with the cranking mass
shows a strong dependence on the nucleus. The change in the $RB_{0/2}$ value
in the light nuclei with the constant mass is relatively moderate. In the
$N=78$ nuclei, however, all the $B(E2)$ ratios considered do not depend on the
nucleus with both the cranking mass and the constant mass. The $B(E2)$ ratios
in both the light $N\sim 28$ nuclei and $N=78$ nuclei with the constant mass
are close to those in the E(5)–$\beta^{4}$ and E(5)–$\beta^{6}$ models. These
results imply that the $0_{2}^{+}$ state in the light nuclei is sensitive to
the collective mass. We are going to investigate the $0_{2}^{+}$ state in
terms of the collective wave functions below. Furthermore, we conclude that
the sensitivity of the collective mass to the low-lying spectra is strong in
the light neutron-rich nuclei around $N=28$ uniquely.

We discuss the structure of the collective wave functions (WFs) of the
$0_{1}^{+}$ and $0_{2}^{+}$ states. First, we look into the WFs with the
constant mass shown in Fig. 9. Here the WFs are multiplied by
$\beta^{4}\sqrt{W(\beta,\gamma)R(\beta,\gamma)}$, the volume element without
$\sin 3\gamma$. The WF of the $0_{1}^{+}$ state is spread over along the
$\gamma$ direction. This is indeed expected from the $\gamma$-soft property in
the PESs. The WF of $0_{2}^{+}$ has a node along the $\beta$ direction. This
is common in 44S and 46S. As expected from the PESs in 44S and 46S, the WFs
with the constant mass look similar to those in the E(5)–$\beta^{4}$ model
shown in the fifth column in Fig. 9. Next, we discuss the WFs with the
cranking mass. The WF of the $0_{1}^{+}$ state is spread along the $\gamma$
direction as obtained with the constant mass and peaks at $\gamma\sim
20^{\circ}$ in 44S and at the oblate side in 46S. The WF of the $0^{+}_{2}$
state in 44S is more or less similar to that obtained with the constant mass,
although the localization around the prolate side is strong. In 46S, the
structure of the WF of $0_{2}^{+}$ is different from that obtained with the
constant mass: The WF has two peaks at the prolate and oblate sides. The two-
peak structure in the collective WF of the $0^{+}_{1}$ and $0^{+}_{2}$ states
is a typical feature of shape coexistence, such as in 72Kr [45]. The
collective WF of the only $0^{+}_{2}$ state in 46S looks similar to the
collective WFs of the $0^{+}_{1}$ and $0^{+}_{2}$ states in 72Kr. Thus, the
$0^{+}_{1}$ and $0^{+}_{2}$ states in 46S are not interpreted as a usual shape
coexistence. This result is rather similar to the one obtained in the
Wilets–Jean model shown in the sixth column in Fig. 9.

##  IV Summary

We have investigated the role of the mass parameters in the collective
Hamiltonian for the triaxial-shape dynamics in neutron-rich nuclei with
$N\simeq 28$. The PESs are obtained by the constrained HFB method with a
Skyrme-type EDF. We found that the PES in 42Mg, 40Si, 44S, and 46S possesses a
topography similar to each other and is soft against triaxial deformation. The
low-lying spectra obtained by assuming the mass parameters as constant are
similar. However, the spectra obtained considering the deformation dependence
of the mass parameters with the cranking approximation show characteristic
features. The second $0^{+}$ state is sensitive to the treatment of the mass
parameters. The energy ratio $R_{0/2}$ and the $B(E2)$ ratio $RB_{0/2}$ show a
strong nucleus dependence. The dependence of $R_{0/2}$ and $RB_{0/2}$ on the
nucleus in neutron-deficient $N=78$ nuclei, which also exhibit the
$\gamma$-soft nature, is less pronounced. We clarified the unique role of the
collective mass in $N\sim 28$ nuclei.

## Acknowledgments

The authors thank M. Kimura and Y. Suzuki for the valuable discussions. This
work was supported by JSPS KAKENHI (Grants No. JP19K03824 and No. JP19K03872)
and the JSPS/NRF/NSFC A3 Foresight Program “Nuclear Physics in the 21st
Century.”

## References

  * Libert  _et al._ [1999] J. Libert, M. Girod, and J.-P. Delaroche, Microscopic descriptions of superdeformed bands with the Gogny force: Configuration mixing calculations in the ${A}\sim 190$ mass region, Phys. Rev. C 60, 054301 (1999).
  * Próchniak  _et al._ [2004] L. Próchniak, P. Quentin, D. Samsoen, and J. Libert, A self-consistent approach to the quadrupole dynamics of medium heavy nuclei, Nucl. Phys. A 730, 59 (2004).
  * Nikšić  _et al._ [2009] T. Nikšić, Z. P. Li, D. Vretenar, L. Próchniak, J. Meng, and P. Ring, Beyond the relativistic mean-field approximation. III. Collective Hamiltonian in five dimensions, Phys. Rev. C 79, 034303 (2009).
  * Delaroche  _et al._ [2010] J.-P. Delaroche, M. Girod, J. Libert, H. Goutte, S. Hilaire, S. Péru, N. Pillet, and G. F. Bertsch, Structure of even-even nuclei using a mapped collective Hamiltonian and the D1S Gogny interaction, Phys. Rev. C 81, 014303 (2010).
  * Matsuyanagi  _et al._ [2016] K. Matsuyanagi, M. Matsuo, T. Nakatsukasa, K. Yoshida, N. Hinohara, and K. Sato, Microscopic derivation of the Bohr-Mottelson collective Hamiltonian and its application to quadrupole shape dynamics, Phys. Scripta 91, 063014 (2016), arXiv:1606.08547 [nucl-th] .
  * Wilets and Jean [1956] L. Wilets and M. Jean, Surface oscillations in even-even nuclei, Phys. Rev. 102, 788 (1956).
  * Iachello [2000] F. Iachello, Dynamic symmetries at the critical point, Phys. Rev. Lett. 85, 3580 (2000).
  * Casten and Zamfir [2000] R. F. Casten and N. V. Zamfir, Evidence for a possible E(5) symmetry in ${}^{134}\mathrm{Ba}$, Phys. Rev. Lett. 85, 3584 (2000).
  * Casten [2006] R. F. Casten, Shape phase transitions and critical-point phenomena in atomic nuclei, Nature Phys. 2, 811 (2006).
  * Bonatsos  _et al._ [2004] D. Bonatsos, D. Lenis, N. Minkov, P. P. Raychev, and P. A. Terziev, Sequence of potentials interpolating between the U(5) and E(5) symmetries, Phys. Rev. C 69, 044316 (2004).
  * Sohler  _et al._ [2002] D. Sohler, Z. Dombrádi, J. Timár, O. Sorlin, F. Azaiez, F. Amorini, M. Belleguic, C. Bourgeois, C. Donzaud, J. Duprat, D. Guillemaud-Mueller, F. Ibrahim, J. A. Scarpaci, M. Stanoiu, M. J. Lopez, M. G. Saint-Laurent, F. Becker, F. Sarazin, C. Stodel, G. Voltolini, S. M. Lukyanov, V. Maslov, Y.-E. Penionzhkevich, M. Girod, S. Péru, F. Nowacki, G. Sletten, R. Lucas, C. Theisen, D. Baiborodin, Z. Dlouhy, J. Mrazek, C. Borcea, A. Bauchet, C. J. Moore, and M. J. Taylor, Shape evolution in heavy sulfur isotopes and erosion of the $N=28$ shell closure, Phys. Rev. C 66, 054302 (2002).
  * Longfellow  _et al._ [2021] B. Longfellow, D. Weisshaar, A. Gade, B. A. Brown, D. Bazin, K. W. Brown, B. Elman, J. Pereira, D. Rhodes, and M. Spieker, Quadrupole collectivity in the neutron-rich sulfur isotopes ${}^{38,40,42,44}\mathrm{S}$, Phys. Rev. C 103, 054309 (2021).
  * Force  _et al._ [2010] C. Force, S. Grévy, L. Gaudefroy, O. Sorlin, L. Cáceres, F. Rotaru, J. Mrazek, N. L. Achouri, J. C. Angélique, F. Azaiez, B. Bastin, R. Borcea, A. Buta, J. M. Daugas, Z. Dlouhy, Z. Dombrádi, F. De Oliveira, F. Negoita, Y. Penionzhkevich, M. G. Saint-Laurent, D. Sohler, M. Stanoiu, I. Stefan, C. Stodel, and F. Nowacki, Prolate-spherical shape coexistence at $N=28$ in ${}^{44}\mathrm{S}$, Phys. Rev. Lett. 105, 102501 (2010).
  * Santiago-Gonzalez  _et al._ [2011] D. Santiago-Gonzalez, I. Wiedenhöver, V. Abramkina, M. L. Avila, T. Baugher, D. Bazin, B. A. Brown, P. D. Cottle, A. Gade, T. Glasmacher, K. W. Kemper, S. McDaniel, A. Rojas, A. Ratkiewicz, R. Meharchand, E. C. Simpson, J. A. Tostevin, A. Volya, and D. Weisshaar, Triple configuration coexistence in ${}^{44}\mathrm{S}$, Phys. Rev. C 83, 061305(R) (2011).
  * Gaudefroy  _et al._ [2009] L. Gaudefroy, J. M. Daugas, M. Hass, S. Grévy, C. Stodel, J. C. Thomas, L. Perrot, M. Girod, B. Rossé, J. C. Angélique, D. L. Balabanski, E. Fiori, C. Force, G. Georgiev, D. Kameda, V. Kumar, R. L. Lozeva, I. Matea, V. Méot, P. Morel, B. S. Nara Singh, F. Nowacki, and G. Simpson, Shell erosion and shape coexistence in ${}_{16}^{43}\mathrm{S}_{27}$, Phys. Rev. Lett. 102, 092501 (2009).
  * Mijatović  _et al._ [2018] T. Mijatović, N. Kobayashi, H. Iwasaki, D. Bazin, J. Belarge, P. C. Bender, B. A. Brown, A. Dewald, R. Elder, B. Elman, A. Gade, M. Grinder, T. Haylett, S. Heil, C. Loelius, B. Longfellow, E. Lunderberg, M. Mathy, K. Whitmore, and D. Weisshaar, Lifetime measurements and triple coexisting band structure in ${}^{43}\mathrm{S}$, Phys. Rev. Lett. 121, 012501 (2018).
  * Longfellow  _et al._ [2020] B. Longfellow, D. Weisshaar, A. Gade, B. A. Brown, D. Bazin, K. W. Brown, B. Elman, J. Pereira, D. Rhodes, and M. Spieker, Shape changes in the $N=28$ island of inversion: Collective structures built on configuration-coexisting states in ${}^{43}\mathrm{S}$, Phys. Rev. Lett. 125, 232501 (2020).
  * Campbell  _et al._ [2006] C. M. Campbell, N. Aoi, D. Bazin, M. D. Bowen, B. A. Brown, J. M. Cook, D.-C. Dinca, A. Gade, T. Glasmacher, M. Horoi, S. Kanno, T. Motobayashi, W. F. Mueller, H. Sakurai, K. Starosta, H. Suzuki, S. Takeuchi, J. R. Terry, K. Yoneda, and H. Zwahlen, Measurement of excited states in ${}^{40}\mathrm{Si}$ and evidence for weakening of the $N=28$ shell gap, Phys. Rev. Lett. 97, 112501 (2006).
  * Bastin  _et al._ [2007] B. Bastin, S. Grévy, D. Sohler, O. Sorlin, Z. Dombrádi, N. L. Achouri, J. C. Angélique, F. Azaiez, D. Baiborodin, R. Borcea, C. Bourgeois, A. Buta, A. Bürger, R. Chapman, J. C. Dalouzy, Z. Dlouhy, A. Drouard, Z. Elekes, S. Franchoo, S. Iacob, B. Laurent, M. Lazar, X. Liang, E. Liénard, J. Mrazek, L. Nalpas, F. Negoita, N. A. Orr, Y. Penionzhkevich, Z. Podolyák, F. Pougheon, P. Roussel-Chomaz, M. G. Saint-Laurent, M. Stanoiu, I. Stefan, F. Nowacki, and A. Poves, Collapse of the $N=28$ shell closure in ${}^{42}\mathrm{S}\mathrm{i}$, Phys. Rev. Lett. 99, 022503 (2007).
  * Takeuchi  _et al._ [2012] S. Takeuchi, M. Matsushita, N. Aoi, P. Doornenbal, K. Li, T. Motobayashi, H. Scheit, D. Steppenbeck, H. Wang, H. Baba, D. Bazin, L. Càceres, H. Crawford, P. Fallon, R. Gernhäuser, J. Gibelin, S. Go, S. Grévy, C. Hinke, C. R. Hoffman, R. Hughes, E. Ideguchi, D. Jenkins, N. Kobayashi, Y. Kondo, R. Krücken, T. Le Bleis, J. Lee, G. Lee, A. Matta, S. Michimasa, T. Nakamura, S. Ota, M. Petri, T. Sako, H. Sakurai, S. Shimoura, K. Steiger, K. Takahashi, M. Takechi, Y. Togano, R. Winkler, and K. Yoneda, Well developed deformation in ${}^{42}\mathrm{Si}$, Phys. Rev. Lett. 109, 182501 (2012).
  * Crawford  _et al._ [2019] H. L. Crawford, P. Fallon, A. O. Macchiavelli, P. Doornenbal, N. Aoi, F. Browne, C. M. Campbell, S. Chen, R. M. Clark, M. L. Cortés, M. Cromaz, E. Ideguchi, M. D. Jones, R. Kanungo, M. MacCormick, S. Momiyama, I. Murray, M. Niikura, S. Paschalis, M. Petri, H. Sakurai, M. Salathe, P. Schrock, D. Steppenbeck, S. Takeuchi, Y. K. Tanaka, R. Taniuchi, H. Wang, and K. Wimmer, First spectroscopy of the near drip-line nucleus ${}^{40}\mathrm{Mg}$, Phys. Rev. Lett. 122, 052501 (2019).
  * Rodríguez-Guzmán  _et al._ [2002] R. Rodríguez-Guzmán, J. L. Egido, and L. M. Robledo, Quadrupole collectivity in $N\approx 28$ nuclei with the angular momentum projected generator coordinate method, Phys. Rev. C 65, 024304 (2002).
  * Rodríguez and Egido [2011] T. R. Rodríguez and J. L. Egido, Configuration mixing description of the nucleus ${}^{44}\mathrm{S}$, Phys. Rev. C 84, 051307(R) (2011).
  * Egido  _et al._ [2016] J. L. Egido, M. Borrajo, and T. R. Rodríguez, Collective and single-particle motion in beyond mean field approaches, Phys. Rev. Lett. 116, 052502 (2016).
  * Yao  _et al._ [2011] J. M. Yao, H. Mei, H. Chen, J. Meng, P. Ring, and D. Vretenar, Configuration mixing of angular-momentum-projected triaxial relativistic mean-field wave functions. II. microscopic analysis of low-lying states in magnesium isotopes, Phys. Rev. C 83, 014308 (2011).
  * Kimura  _et al._ [2013] M. Kimura, Y. Taniguchi, Y. Kanada-En’yo, H. Horiuchi, and K. Ikeda, Prolate, oblate, and triaxial shape coexistence, and the lost magicity of $N=28$ in ${}^{43}\mathrm{S}$, Phys. Rev. C 87, 011301(R) (2013).
  * Suzuki and Kimura [2021] Y. Suzuki and M. Kimura, Triaxial deformation and the disappearance of the $N=28$ shell gap, Phys. Rev. C 104, 024327 (2021).
  * Suzuki  _et al._ [2022] Y. Suzuki, W. Horiuchi, and M. Kimura, Erosion of $N=28$ shell closure: Shape coexistence and monopole transition, Prog. Theor. Exp. Phys. 2022, 063D02 (2022).
  * Utsuno  _et al._ [2012] Y. Utsuno, T. Otsuka, B. A. Brown, M. Honma, T. Mizusaki, and N. Shimizu, Shape transitions in exotic Si and S isotopes and tensor-force-driven Jahn-Teller effect, Phys. Rev. C 86, 051301(R) (2012).
  * Utsuno  _et al._ [2015] Y. Utsuno, N. Shimizu, T. Otsuka, T. Yoshida, and Y. Tsunoda, Nature of isomerism in exotic sulfur isotopes, Phys. Rev. Lett. 114, 032501 (2015).
  * Yoshida [2022] K. Yoshida, Cranked Skyrme-Hartree-Fock-Bogoliubov approach for a mean-field description of nuclear rotations near the drip line, Phys. Rev. C 105, 024313 (2022).
  * Li  _et al._ [2011] Z. P. Li, J. M. Yao, D. Vretenar, T. Nikšić, H. Chen, and J. Meng, Energy density functional analysis of shape evolution in $N=28$ isotones, Phys. Rev. C 84, 054304 (2011).
  * Yoshida [2009] K. Yoshida, Skyrme-QRPA calculations for low-lying excitation modes in deformed neutron-rich nuclei, Eur. Phys. J. A 42, 583 (2009), arXiv:0902.3053 [nucl-th] .
  * Inglis [1956] D. R. Inglis, Nuclear moments of inertia due to nucleon motion in a rotating well, Phys. Rev. 103, 1786 (1956).
  * Beliaev [1961] S. T. Beliaev, Concerning the calculation of the nuclear moment of inertia, Nucl. Phys. 24, 322 (1961).
  * Girod and Grammaticos [1979] M. Girod and B. Grammaticos, The zero-point energy correction and its effect on nuclear dynamics, Nucl. Phys. A 330, 40 (1979).
  * Gall  _et al._ [1994] B. Gall, P. Bonche, J. Dobaczewski, H. Flocard, and P.-H. Heenen, Superdeformed rotational bands in the mercury region. A cranked Skyrme-Hartree-Fock-Bogoliubov study, Z. Phys. A 348, 183 (1994).
  * Terasaki  _et al._ [1995] J. Terasaki, P.-H. Heenen, P. Bonche, J. Dobaczewski, and H. Flocard, Superdeformed rotational bands with density dependent pairing interactions, Nucl. Phys. A 593, 1 (1995).
  * Bartel  _et al._ [1982] J. Bartel, P. Quentin, M. Brack, C. Guet, and H.-B. Håkansson, Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force, Nucl. Phys. A 386, 79 (1982).
  * Yamagami  _et al._ [2009] M. Yamagami, Y. R. Shimizu, and T. Nakatsukasa, Optimal pair density functional for the description of nuclei with large neutron excess, Phys. Rev. C 80, 064301 (2009).
  * Bonche  _et al._ [1985] P. Bonche, H. Flocard, P.-H. Heenen, S. J. Krieger, and M. S. Weiss, Self-consistent study of triaxial deformations: Application to the isotopes of Kr, Sr, Zr and Mo, Nucl. Phys. A 443, 39 (1985).
  * Bonche  _et al._ [2005] P. Bonche, H. Flocard, and P.-H. Heenen, Solution of the Skyrme HF+BCS equation on a 3D mesh, Comput. Phys. Commun. 171, 49 (2005).
  * Bohr [1952] A. Bohr, Mat. Fys. Medd. Dan. Vid. Selsk. 26, 14 (1952).
  * [44] National Nuclear Data Center, Brookhaven National Laboratory, https://www.nndc.bnl.gov/nudat3/.
  * Sato and Hinohara [2011] K. Sato and N. Hinohara, Shape mixing dynamics in the low-lying states of proton-rich Kr isotopes, Nucl. Phys. A 849, 53 (2011).


11institutetext: University of Oklahoma, Norman OK 73019, USA 11email:
<EMAIL_ADDRESS>
# Observe Locally, Classify Globally: Using GNNs to Identify Sparse Matrix
Structure

Khaled Abdelaal 11 0000-0001-8963-7523 Richard Veras 11 0000-0003-2633-3391

###### Abstract

The performance of sparse matrix computation highly depends on the matching of
the matrix format with the underlying structure of the data being computed on.
Different sparse matrix formats are suitable for different structures of data.
Therefore, the first challenge is identifying the matrix structure before the
computation to match it with an appropriate data format. The second challenge
is to avoid reading the entire dataset before classifying it. This can be done
by identifying the matrix structure through samples and their features. Yet,
it is possible that global features cannot be determined from a sampling set
and must instead be inferred from local features. To address these challenges,
we develop a framework that generates sparse matrix structure classifiers
using graph convolutional networks. The framework can also be extended to
other matrix structures using user-provided generators. The approach achieves
97% classification accuracy on a set of representative sparse matrix shapes.

###### Keywords:

Sparse Matrix Graph Neural Networks Classification

##  1 Introduction

(a)

(b)

| Diagonal | Random | Rand+Diag | Kronecker
---|---|---|---|---
Diagonal | 500 | 0 | 0 | 0
Random | 0 | 490 | 10 | 0
Rand+Diag | 0 | 40 | 460 | 0
Kronecker | 0 | 20 | 0 | 480

(c)

Figure 1: Efficacy of our classifier framework at determining structure when
the data is permuted. (a) is an off-diagonal matrix, (b) is a re-labelled
variant of (a), and (c) is the confusion matrix for the classifier framework
on re-labelled matrices similar to (b). By using GCNs our approach is
invariant to node labelling and achieves around 97% accuracy.

Sparse matrices represent a fundamental building block used throughout the
field of scientific computing in applications, such as graph analytics,
machine learning, fluid mechanics, and finite element analysis [13, 6]. Such
matrices appear as operands in numerous fundamental computational kernels such
as sparse matrix-vector multiplication (SpMV), Cholesky factorization, LU
factorization, sparse matrix -dense matrix multiplication, and matricized
tensor times Khatri-Rao product (MTTKRP) among others. Building efficient
algorithms for this class of kernels mainly depends on the storage format used
for the sparse matrix as observed in different studies [4, 1]. A variety of
such formats are proposed in literature [9, 10]. Hence, it is crucial to
identify the structure of the matrix to choose the ideal sparse format, and
eventually tailor the algorithm to that format to optimize the workload
performance. However, identifying the structure of the matrix is not always
trivial. Figure 1(a) shows a spy plot of an off-diagonal sparse matrix, and
Figure 1(b) shows the same matrix, with some of the original row indices and
column indices re-labelled. It is less obvious for the latter figure to
provide an insight of the original structure of the non-zeros within the
matrix. Additionally, in case of huge sparse matrices, we might only have
access to samples of the matrix. This could be because of computational or
storage restrictions, or missing data. In these two cases (re-labelling and
sub-sampling), we need efficient techniques to recognize the shape of the
input matrix.

To tackle mentioned issues, we propose a framework to identify sparse matrices
structures, using graph neural networks. Figure 1(c) shows the confusion
matrix for the proposed framework using four sample classes on re-labelled
variants. The framework design is modular, allowing users to easily augment it
with new structures generators or feature sets. The main contributions of this
paper are as follows:

  * •

Proposing a novel, modular Graph Neural Network framework to accurately
predict the shapes of sparse matrices, including partial samples, and re-
labelled variants of original matrices.

  * •

Presenting a new balanced synthetic dataset for structured sparse matrices.

  * •

Providing a performance analysis of graph-level classification on sparse
matrices, using different feature sets.

  * •

Introducing two new compact and efficient feature sets for matrices as graphs,
namely: Linear and Exponential Binned One-Hot Degree Encoding.

The rest of this paper is organized as follows: Section 2 introduces the
necessary background, Section 3 details the design of the proposed framework,
Section 4 discusses the evaluation and results of the framework, while Section
5 describes related work. Finally, Section 6 summarizes the findings of the
paper.

##  2 Background

###  2.1 Graph Neural Networks

Graph neural networks (GNNs) [17] are a class of deep learning models that
operate on graphs or networks. Unlike traditional neural networks that operate
on structured data such as images or sequences, GNNs can handle arbitrary
graph structures with varying node and edge attributes, enabling them to learn
powerful representations of graph-structured data. The key idea behind GNNs is
to iteratively update node embeddings by aggregating information from the
embeddings of their neighbors through the "graph convolution" operation. By
stacking multiple layers of graph convolution and non-linear activation
functions, GNNs can learn hierarchical representations of the graph that
capture both local and global information.

###  2.2 Structured Matrices

Several common structures are observed in sparse matrices, such as:

Diagonal all non-zeros are located on the main or a secondary diagonal. This
structure represents a 1D mesh and commonly appears in various scientific and
engineering applications.

Random the non-zero elements are randomly distributed across the matrix, with
variable density. Such matrices have no specific identifiable structure.

Kronecker Graphs [11] are a class of synthetic graphs that have been widely
used to model real-world networks, and are generated by recursively applying
the Kronecker product of a small base graph with itself. Let $A$ and $B$ be
two matrices. Then, their Kronecker product $A\otimes B$ is given by

$A\otimes B=\begin{pmatrix}a_{11}B&\cdots&a_{1n}B\\\ \vdots&\ddots&\vdots\\\
a_{m1}B&\cdots&a_{mn}B\end{pmatrix}$ (1)

where $a_{ij}$ are the entries of $A$. We use these three classes of
structures, and a combination of them, as a representative set that can be
combined to form more complex relationships [16, 20]. Our framework is not
limited to only these structures, and they serve as an example to evaluate its
performance.

(a) Diagonal

(b) Random

(c) Random+Diagonal

(d) Kronecker Graph

Figure 2: Global Degree Distribution for samples in each matrix (graph) class
studied in this paper. In our approach we classify the shape based on local
views from sampled data.

Degree as a representative node feature Figure 2 illustrates that one can
accurately distinguish between the different classes based on the degree
distribution of the representative graph. For example, for Diagonal matrices
(Figure 2(a)), the degree for all nodes is low, and is either constant or
linear across all nodes. Kronecker graphs follow a power-law degree
distribution, with only a few nodes having many connections (high degree) and
most of the nodes having relatively few connections (low degree). However,
only the local per-node degree view may be immediately available, and not the
global graph view. An example of such a case is only having a sample of the
graph and not the entire graph due to storage or computational limitations.
The power of GNNs can be leveraged to carry out the required task: the
prediction of the sample matrix structure.

##  3 Framework Description

Figure 3: High-Level overview of the framework. It consists of three main
phases: dataset generation, where the synthetic sparse matrices are generated,
prepared as graphs, and have feature set attached. Then, the GNN model
training using 5-fold cross validation to capture the model performance, and
then generate a trained model instances, that is used later in the inference
phase.

The goal of the proposed framework is to predict the structure of the input
sparse matrix through its classification into one of the configured target
classes. We use diagonal, random, diagonal+random, and Kronecker graph as
examples of these classes to evaluate the performance of the framework. New
structures can be seamlessly integrated. Figure 3 shows a high-level overview
of the proposed framework. It consists of three main stages: Dataset
generation, Model Training, and Inference. A synthetic dataset is generated
using different generators for different shapes of matrices, which are then
represented as graphs. In the training phase, we use GNN with 5-fold cross
validation to evaluate the model performance. Finally, the trained model
instance is used for later inference.

###  3.1 Dataset Generation

We generate a balanced dataset of 40K synthetic sparse matrices, covering the
four sample classes through individual generators. Each of the generators
returns a Coordinate (COO) representation for the matrix, excluding the actual
non-zero values. The COO representation is then used as the adjacency list to
build the graph representation.

###  3.2 Feature Set Selection

A per-node feature vector is necessary for the graph neural network to
classify matrices. Node degree can be calculated for rows/columns in input
matrices from their graph representation.

One-Hot Degree Encoding uses a number of features equal to the maximum degree
+ 1. A limitation of this encoding is that it requires the knowledge of the
maximum degree in the entire dataset before training. Also, the required
storage is proportional to the maximum degree recorded in the dataset, which
increases memory requirements and reduces maximum possible batch size during
training. Moreover, it poses complications during inference if the input
matrix has a degree greater than the maximum degree in the training set. In
our dataset, the maximum training set degree is 7710, so the length of one-hot
encoding feature vector per node is 7711, limiting the maximum batch size on
GPU to only 1 graph.

Local Degree Profile (LDP) [3] captures the local structural information of
nodes in their immediate neighborhood. LDP is calculated for each node as a
five-feature vector: the node degree, the minimum degree of its neighbors, the
maximum degree of its neighbors, the mean degree for its neighbors, and the
standard deviation of the degrees of each neighbors. LDP features are easy to
compute for any given graph. Additionally, the number of features per node is
fixed, regardless of the used training data. LDP incurs low storage overhead.

Linear Binned One-Hot Degree Encoding (LBOH) We implement a modified version
of one-hot encoding, to address its limitations. LBOH works by having a fixed
number of buckets for representing one-hot degrees. Buckets ranges are
designed as follows: a set individual sequential buckets from 0 t $\alpha$
(inclusive) where $\alpha$ is a small integer (less than 10). Then, we add a
set of buckets with fixed step $\beta$ from $\alpha$: $(\alpha+\beta)$,
$(\alpha+2\beta)$, …, $(\alpha+k\beta)$ where $(\alpha+k\beta)$ is the maximum
degree threshold. Any degree greater than $(\alpha+k\beta)$ is mapped to the
final bucket.

Figure 4: An example of finding the linear binned one-hot degree encoding for
a node with degree = 5, where the parameters for the encoding scheme are
$\alpha=5$, $\beta=3$, and $k=2$. Degree 5 is mapped to its associated bucket
(5 to 7), then the bucket index (5) is represented using one-hot encoding (1
at the position where the value 5 exists, 0 otherwise).

Figure 4 shows an example of LBOH encoding. As opposed to One-Hot Encoding,
LBOH provides a fixed number of features regardless of the maximum degree in
the training dataset. At the inference stage, only the values of $\alpha$,
$\beta$, and $k$ are needed.

Exponential Binned One-Hot Degree Encoding (EBOH) The main difference between
EBOH and LBOH is the kind of step between buckets ranges. Instead of a linear
step in LBOH, EBOH uses an exponential step to cover more degree values with a
small number of features. First, the value of $\alpha$ is chosen such that
$1\leq\alpha\leq 3$. Then for the buckets, a sequential one-to-one mapping is
performed for values 0 through $2^{\alpha}$. For the following buckets, the
upper bound (inclusive) is $2^{\alpha+i}$ where $i\in[1,k]$ and
$k\in\mathbb{N}^{*}$.

Figure 5: An example of finding the exponential binned one-hot degree encoding
for a node with degree = 7880, where the parameters for the encoding scheme
are $\alpha=2$ and $k=3$. Degree 7880 is mapped to its associated bucket (33
to $\infty$), then the bucket index (8) is represented using one-hot encoding
(1 at the position where the value 8 exists, 0 otherwise).

Figure 5 shows an example of EBOH encoding. EBOH encoding still provides the
benefit of having the number of features independent of the maximum degree in
the training set.

###  3.3 The Graph Neural Network Architecture

InputGraph ConvolutionReLUGraph ConvolutionReLUGraph ConvolutionGlobal Mean
PoolingDropoutLinearOutput

Figure 6: Graph Neural Network architecture.

To identify the structure of the input matrix, the matrix is viewed as the
adjacency list of a graph, enabling the use of machine learning methods
designed for graph data. GNNs provide additional benefits such as allowing the
use of matrices (graphs) of arbitrary sizes as input, Also, GNNs are agnostic
to node ordering. This powerful property enables re-labelling or permuting
nodes in a graph representing a sparse matrix, while maintaining accurate
predictions. The machine learning task of interest is graph-level prediction
since a single label (class) is needed for an entire graph (matrix). The GNN
architecture is shown in Figure 6. The hidden layers are three graph
convolution layers and one linear (output) layer. Graph convolution is an
operation where node embeddings are iteratively generated as the aggregations
from the node neighborhoods. This operation is used to capture complex
features of the graph. The first convolution layer aggregates information from
the local neighborhood of each node. This operation is repeated in subsequent
convolution layers in order to propagate information to increasingly larger
neighborhoods. By the end of three convolution layers, the model has learned a
hierarchical representation of the graph, where the features at each layer
capture increasingly complex structural patterns. The learned representation
so far is "node embeddings". Then, learned node embeddings are reduced into a
single graph embedding using a global mean pooling operation (called readout
layer). Samples are randomly dropped out to reduce overfitting. Finally, a
linear classifier is applied to the graph embedding.

⬇

def generateDiagRandom(size, threshold=2):

""" A function to generate a Diag+Random square matrix """

tuples = [(x,y) for x in range(size) for y in range(size) if
(random.randint(0,10) <= threshold or x == y)]

# seperate tuples into two lists: the row array and the column array

coo_rep = list(map(list, zip(*tuples)))

return coo_rep[0], coo_rep[1], [size, size]

⬇

def process(self):

catMap = [...., {

# Number of instances to generate for this class

’num_iter’:10000,

# Name of the generator function

’generator’:generateDiagRandom,

# A string list of required generator function param

’gen_params’:[’random.randint(MIN_DIM_SIZE,MAX_DIM_SIZE)’] }]

Figure 7:  The two steps needed to add a new class to the classifier
framework. First (top), create a new generator function in the generators
file, and second (bottom), add a dictionary entry to catMap list in the
process method of the dataset class.

Modularity New shapes of matrices can be easily integrated in our framework.
To achieve this, two steps are needed as shown in Figure 7: (1) write a
generator for that new shape, and (2) add an entry to the categories (shapes)
map in the dataset class for this shape, containing the number of dataset
instances to generate, the name of the generator function, and the different
required parameters. The generator is required to return the COO
representation excluding values, and the matrix dimensions. After generating
the new dataset instances for this class, one does not need to re-train the
entire model again. Transfer learning [21] can be used to replace the last
layer of the trained model with a new layer that has the appropriate number of
outputs, after introducing the new shape(s). Then, the weights of all previous
layers are frozen and only the new layer is trained. Another aspect of
modularity in our framework is the ability to seamlessly attach different
feature sets. Feature sets are only computed when the graph is queried. To
implement a new feature set, a modification to the get method of the dataset
is needed. This method first reads in the graph file from disk, calculates the
new feature set, and attaches it to the graph.

##  4 Analysis

We run a set of experiments to evaluate the accuracy of our approach in
detecting structures, using the different feature representation discussed in
Section 3.2.

###  4.1 Evaluation

Experimental Setup Table 1 shows the experimental setup and learning
parameters used in the experiments. We use PyTorch Geometric [8] for the GNN.

Table 1: Experimental setup and Training parameters used in the experiments. * Batch size used for traditional one-hot encoding is 1. Component | Specification | Parameter | Value
---|---|---|---
GPU | NVIDIA RTX A6000 | Optimizer | Adam
GPU Memory | 48 GB GDDR6 | Learning Rate | 0.01
CUDA Version | 11.8 | Error Criterion | Cross Entropy
Main Memory | 64 GB DDR4 | Batch Size | 256*
Operating System | Ubuntu 22.04 | Cross Validation Folds | 5

Evaluation Metrics To evaluate the prediction accuracy of the framework, four
derived metrics are used: accuracy, precision, recall, and F1-score. We report
per-class and overall accuracy and F1 score numbers, since the latter is the
harmonic mean of precision and recall. Using both accuracy and F1-score helps
provide a more comprehensive evaluation of the framework’s performance.
Accuracy gives an overall view of how well the classifier is performing, while
the F1-score provides insights into its ability to correctly classify positive
instances.

###  4.2 Results

Table 2: Performance of the classifier for different degree representations

| Degree Representation
---|---
| One-Hot Encoding | LDP | LBOH | EBOH
| Accuracy | F1 Score | Accuracy | F1 Score | Accuracy | F1 Score | Accuracy | F1 Score
Diagonal | 1.0 | 0.90 | 1.0 | 0.97 | 1.0 | 1.0 | 1.0 | 1.0
Random | 0.90 | 0.91 | 0.64 | 0.76 | 0.92 | 0.95 | 0.95 | 0.96
Random+Diagonal | 0.86 | 0.99 | 0.98 | 0.83 | 0.96 | 0.94 | 0.97 | 0.96
Kronecker | 0.90 | 0.94 | 0.90 | 0.93 | 0.98 | 0.97 | 0.96 | 0.98
Overall | 0.90 | 0.90 | 0.88 | 0.88 | 0.97 | 0.97 | 0.97 | 0.98

Classification Performance Table 2 shows the accuracy and F1 score for the
classifier using the different degree representations discussed in Section
3.2. Performance results show that both LBOH, and EBOH provide high prediction
accuracy of around 97% and a F1 score of around 98%. On the other hand,
traditional one-hot encoding exhibits a lower accuracy of around 90%. One-Hot
Encoding requires a significantly large number of features per node (7711),
limiting the training batch size on the A6000 GPU to only one graph. This
forces the optimizer to adjust the neural network weights very frequently,
hence hurting the overall accuracy. Using LDP as a feature set exhibits
variant model performance across folds depending on the validation set being
used. In some folds, LDP provides high accuracy of around 97% to 98% similar
to EBOH. In other folds, LDP fails to converge to an acceptable loss value,
and ends up with an accuracy of around 74% on the last few epochs. This
performance variance across folds deems LDP unfit for the purposes of our
application. It significantly fails in two classes: Random and Kronecker. It
predicts Random matrices as Random+Diagonal for more than 32.5% of the
instances. This is likely due to the prevalence of the local degree neighbor
summary features (the last four LDP features) instead of focusing on the node
degree. This eventually results in failing to discover the global hierarchical
structures in the matrix. LDP still shows perfect accuracy in case of diagonal
matrices since almost all nodes in the matrix’s graph have the same degree.
LDP prediction quality for Kronecker graphs is also lower than other evaluated
feature sets (around 81% in some folds) for the same reasons.

(a)

(b)

Figure 8: Cross Entropy Loss across different folds in 5-fold cross validation
training using (a) EBOH, and (b) LDP feature set.

Figure 8(a) demonstrates the validation loss across the 5 different folds for
EBOH. It shows almost no variance in the loss across the different folds,
indicating the stability of the model’s performance across folds. On the other
hand, Figure 8(b) shows the validation loss for LDP and illustrates that the
loss does not converge in 2 out of 5 folds.

Classifying Sub-samples and Re-labelled Subgraphs To test the efficacy of GNN
on both aspects, we generate 200 new matrices: 50 for each of the four
classes. For each of them, we generate 10 subgraphs and 10 re-labelled
variants. To generate the subgraphs, we use uniform random node sampling
(URNS) [12]: nodes are randomly selected with uniform probability, as well as
the edges connecting the selected nodes. Re-labelling of a graph $G$ simply
renames the nodes $V$ of the graph, and produces a new graph $G^{\prime}$ with
the same size and degree distribution of the original graph $G$. Figure 9
shows an example of both URNS and random re-labelling.

Table 3: Accuracy comparison for node sampling, node re-labelling, and original graphs using EBOH feature set. Class |  Node Sampling |  Node Re-labelling |  Original Graphs
---|---|---|---
Diagonal |  1 |  1 |  1
Random |  0.83 |  0.98 |  0.98
Random+Diagonal |  0.92 |  0.92 |  0.92
Kronecker |  0.94 |  0.96 |  0.96
Overall |  0.92 |  0.97 |  0.97

Table 3 shows the model’s performance on subgraphs and re-labelled variants as
compared to original full graphs. The table shows that re-labelling node has
no impact on the classification accuracy; it shows the same overall accuracy
of around 97% which is observed for the original graphs. This is expected
because the arrangement of nodes in a graph is irrelevant, since the graph has
the same degree distribution. For subgraphs generated using URNS of larger
graphs, the overall accuracy drops to around 92%. The reason being that random
node sampling can alter the degree distribution of the graph. The random
choice of nodes can result in either isolated nodes (no edges) or much lower
degree nodes as compared to the original graph. This affects the accuracy
specially for complex shapes such as random and Kronecker graphs. One way to
reduce the accuracy loss for samples is to use a more sophisticated graph
sampling technique rather than randomly selecting nodes or edges.

$a$$b$$c$$d$$e$$f$Original graph$b$$c$$d$Subgraph$c$$a$$b$$f$$e$$d$Re-labelled
graph

Figure 9: Example of generating a random sub-sample and a re-labelled variant
of a graph. The original graph (left) contains six nodes. Using URNS, a random
subgraph (middle) of three nodes is generated. A re-labelled variant (right)
is generated using a random 1:1 mapping between the original and new node
labels.

##  5 Related Work

Prediction on Sparse Matrices Several studies investigated the use of machine
learning to predict the optimal sparse format for SpMV on CPU and GPU [23, 19,
18, 15, 14, 2]. Our framework does not directly predict the best sparse
format, instead, we only predict the structure of the input matrix. This
allows de-coupling the sparsity pattern from the sparse format, following the
argument adopted by AlphaSparse [7] since our framework also allows the
seamless integration of new classes. Existing techniques collect a set of
features from each matrix such as: the number of diagonals, the ratio of true
diagonals to total diagonals, the (maximum) number of non-zeros per row, the
variation of the number of nonzeros per row, the ratio of nonzeros in DIA and
ELL data structures, and a factor or power-law distribution. We only need to
calculate one feature per node: its degree. Also, [23] uses a CNN approach to
treat matrices as images, and in order to fix the size of the matrix, they
normalize input matrices into fixed size blocks, losing partial matrix
information in the process. In contrast, our approach handles arbitrary sizes
of matrices, without losing precision, leveraging the power of Graph Neural
Networks. We can optionally sample large matrices and maintain high prediction
accuracy. An additional benefit to our framework is that it is order in-
variant, since matrices are represented as graphs.

Graph Representation for Learning Representing non-attribute graphs is an open
problem [5]. Common approaches employ graph properties such as node degree,
more specifically a one-hot encoding of the degree [22]. One-hot encoding
suffers from numerous limitations (Section 3.2). LDP [3] provides a compact
representation for graph using five features per node. Although the
computation of such feature vector is efficient, using LDP results in
unreliable model performance for our task (Section 4.2). Both our
representations (LBOH and EBOH) outperform one-hot encoding and LDP while
addressing their shortcomings.

##  6 Summary

In this paper, we proposed a GNN based framework to classify structured sparse
matrices. We introduced two novel non-attribute graph representations based on
node degrees: LBOH, and EBOH. We evaluated the efficacy of our framework on a
synthetic, balanced dataset of matrices that we generated containing random
matrices from four sample classes: diagonal, random, random+diagonal, and
Kronecker graphs. Performance results demonstrate a high classification
accuracy of 97% for the framework when using our feature sets: LBOH and EBOH.
They also show high accuracy of 92% and 97% on random node subsamples and re-
labelled variants respectively. Our framework is modular, allowing the
inclusion of additional classes with minimal user effort. Future endeavors
target the automatic generation of the optimal sparse data format and
algorithm for sparse matrix kernels, using the obtained prediction results
from the current framework.

## References

  * [1] Bell, N., Garland, M.: Implementing sparse matrix-vector multiplication on throughput-oriented processors. In: Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis. pp. 1–11 (2009)
  * [2] Benatia, A., Ji, W., Wang, Y., Shi, F.: Sparse matrix format selection with multiclass svm for spmv on gpu. In: 2016 45th International Conference on Parallel Processing (ICPP). pp. 496–505 (2016)
  * [3] Cai, C., Wang, Y.: A simple yet effective baseline for non-attribute graph classification. arXiv preprint arXiv:1811.03508 (2018)
  * [4] Choi, J.W., Singh, A., Vuduc, R.W.: Model-driven autotuning of sparse matrix-vector multiply on gpus. Association for Computing Machinery, New York, NY, USA (2010)
  * [5] Cui, H., Lu, Z., Li, P., Yang, C.: On positional and structural node features for graph neural networks on non-attributed graphs. Association for Computing Machinery, New York, NY, USA (2022)
  * [6] Davis, T.A., Hu, Y.: The university of florida sparse matrix collection 38(1) (2011)
  * [7] Du, Z., Li, J., Wang, Y., Li, X., Tan, G., Sun, N.: Alphasparse: Generating high performance spmv codes directly from sparse matrices. In: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis. SC ’22, IEEE Press (2022)
  * [8] Fey, M., Lenssen, J.E.: Fast graph representation learning with PyTorch Geometric. In: ICLR Workshop on Representation Learning on Graphs and Manifolds (2019)
  * [9] Filippone, S., Cardellini, V., Barbieri, D., Fanfarillo, A.: Sparse matrix-vector multiplication on gpgpus 43(4) (2017)
  * [10] Langr, D., Tvrdík, P.: Evaluation criteria for sparse matrix storage formats. IEEE Transactions on Parallel and Distributed Systems 27(2), 428–440 (2016)
  * [11] Leskovec, J., Chakrabarti, D., Kleinberg, J., Faloutsos, C., Ghahramani, Z.: Kronecker graphs: An approach to modeling networks. Journal of Machine Learning Research 11(33), 985–1042 (2010)
  * [12] Leskovec, J., Faloutsos, C.: Sampling from large graphs. Association for Computing Machinery, New York, NY, USA (2006)
  * [13] Leskovec, J., Krevl, A.: SNAP Datasets: Stanford large network dataset collection. http://snap.stanford.edu/data (Jun 2014)
  * [14] Li, J., Tan, G., Chen, M., Sun, N.: Smat: An input adaptive auto-tuner for sparse matrix-vector multiplication 48(6) (2013)
  * [15] Li, K., Yang, W., Li, K.: Performance analysis and optimization for spmv on gpu using probabilistic modeling. IEEE Transactions on Parallel and Distributed Systems 26(1), 196–205 (2015)
  * [16] Puschel, M., Moura, J., Johnson, J., Padua, D., Veloso, M., Singer, B., Xiong, J., Franchetti, F., Gacic, A., Voronenko, Y., Chen, K., Johnson, R., Rizzolo, N.: Spiral: Code generation for dsp transforms. Proceedings of the IEEE 93(2), 232–275 (2005)
  * [17] Scarselli, F., Gori, M., Tsoi, A.C., Hagenbuchner, M., Monfardini, G.: The graph neural network model. IEEE Transactions on Neural Networks 20(1), 61–80 (2009)
  * [18] Su, B.Y., Keutzer, K.: Clspmv: A cross-platform opencl spmv framework on gpus. Association for Computing Machinery, New York, NY, USA (2012)
  * [19] Tan, G., Liu, J., Li, J.: Design and implementation of adaptive spmv library for multicore and many-core architecture 44(4) (2018)
  * [20] Van Loan, C.: Computational frameworks for the fast Fourier transform. SIAM (1992)
  * [21] Weiss, K., Khoshgoftaar, T.M., Wang, D.: A survey of transfer learning. Journal of Big data 3(1), 1–40 (2016)
  * [22] Xu, K., Hu, W., Leskovec, J., Jegelka, S.: How powerful are graph neural networks? In: International Conference on Learning Representations (2019)
  * [23] Zhao, Y., Li, J., Liao, C., Shen, X.: Bridging the gap between deep learning and sparse matrix format selection. Association for Computing Machinery, New York, NY, USA (2018)