# Representing individual electronic states in crystals for machine learning
quasiparticle band structures

Nikolaj Rørbæk Knøsgaard  Computational Atomic-scale Materials Design (CAMD),
Department of Physics, Technical University of Denmark, 2800 Kgs. Lyngby
Denmark  Kristian Sommer Thygesen  Computational Atomic-scale Materials Design
(CAMD), Department of Physics, Technical University of Denmark, 2800 Kgs.
Lyngby Denmark

## Abstract

We address the problem of representing quantum states of electrons in a solid
for the purpose of machine leaning state-specific electronic properties.
Specifically, we construct a fingerprint based on energy decomposed operator
matrix elements (ENDOME) and radially decomposed projected density of states
(RAD-PDOS), which are both obtainable from a standard density functional
theory (DFT) calculation. Using such fingerprints we train a gradient boosting
model on a set of 46k G0W0 quasiparticle energies. The resulting model
predicts the self-energy correction of states in materials not seen by the
model with a mean absolute error of 0.14 eV. By including the material’s
calculated dielectric constant in the fingerprint the error can be further
reduced by 30%, which we find is due to an enhanced ability to learn the
correlation/screening part of the self-energy. Our work paves the way for
accurate estimates of quasiparticle band structures at the cost of a standard
DFT calculation.

##  1 Introduction

The electronic band structure is one of the most fundamental and important
characteristics of a crystalline solid. It relates the quantum mechanical
energy levels of an electron in the solid to its (crystal) momentum and
provides the basis for describing and understanding a range of materials
properties. As a consequence, the accurate prediction of electronic band
structures represents a corner-stone problem of computational condensed matter
physics.

Density functional theory (DFT)[1] with semi-local exchange-correlation
functionals[2] is the standard method for solving the electronic structure
problem of materials from first principles. However, the DFT single-particle
energies do not in general provide an accurate model for the electronic band
structure.[3] Instead, the gold standard for band structure calculations is
represented by the GW self-energy method[4], which provides the true
quasiparticle (QP) band structure, i.e. it goes beyond a mean field
description by explicitly accounting for exchange and many-body screening
effects.[5, 6] In Ref. [7] the mean absolute error on the calculated band gap
relative to experimental references for a set of ten simple semiconductors and
insulators was found to be 2.05 eV for DFT-LDA and 0.31 eV for non-
selfconsistent G0W0@LDA[huser2013quasiparticle]. Very similar results have
been found in other studies.[8, 9] The improved accuracy of the GW method
comes at the price of a significantly more involved methodology and a much
higher computational cost. In practice, this means that GW calculations are
limited to small-scale studies of relatively simple materials.

Recently, machine learning (ML) has attracted widespread interest as a means
to predict materials properties without performing expensive quantum
mechanical calculations.[10, 11, 12, 13, 14, 15] In the context of band gap
predictions, Zhou _et al._ trained a support vector machine on 3896
experimental band gaps using a representation based only on elemental
properties of the constituent atoms.[16] Rajan _et al._ used different
regressions methods to predict band gaps of MXene crystals using a training
set of 76 G0W0 band gaps and a representation encoding atomic and structural
properties.[17] Liang _et al._ used a representation based on atomic ionicity
descriptors to predict GW band gaps of a set of 2D semiconductors.[18] In all
these previous studies, the ML model was trained to predict the size of the
band gap rather than the full $k$-resolved band structure. Thereby, important
information is missed including the type of the band gap (direct or indirect),
the curvature of the valence and conduction bands at the extrema points
(effective masses), and the position and dispersion of other bands away from
the band gap. Predicting the full band structure directly from the atomic
structure of the material is a daunting challenge that, although possible in
principle, would require highly sophisticated ML models and immense amounts of
training data.

Here we take a different approach, in which the output from a DFT calculation
is taken as input to a ML model to predict the full GW band structure. The
philosophy behind our approach is that standard DFT calculations are
computationally very cheap, in particular compared to GW, and although they do
not directly produce the desired precision, they hold the gist of the
material’s genome and thus should provide an excellent starting point for
accurate property predictions. In our scheme, the rich, but unmanageable,
information contained in the DFT wave functions is encoded into low
dimensional fingerprints via energy resolved orbital projections and operator
matrix elements. These state-specific electronic fingerprints provide a
description of the ”local environment” of a given electronic eigenstate in the
infinite dimensional Hilbert space, and are thus analogue to the well known
fingerprints used to describe atoms in chemical environments [19].

Using a data set of 286 G0W0 band structures of non-magnetic 2D semiconductors
comprising a total of 46.000 $(\varepsilon_{nk}^{\mathrm{QP}},k)$ pairs, we
train a gradient boosting algorithm to predict the G0W0 correction of an
eigenstate from its DFT fingerprint. The method achieves a mean absolute error
(MAE) of 0.14 eV for individual band energies and 0.18 eV for the band gap.
These deviations are significantly smaller than the typical size of the G0W0
corrections and also lower than the accuracy of the G0W0 method itself. The
model can be further and significantly improved by adding the static
electronic polarisability to the fingerprint. A SHAP feature analysis reveals
that the inclusion of the polarisability allows the ML model to distinguish
between materials with similar PBE band structures but different dielectric
screening properties, which is directly related to the size of the GW
correction.

##  2 Results

Figure 1a shows an example of a PBE (orange) and G0W0 (green) band structure
for monolayer MoS2 (note that spin-orbit interactions are not included
throughout this work). It is clear that there are significant differences
between the two descriptions. First of all, G0W0 yields a QP band gap of 2.53
eV in good agreement with the experimental value of 2.5 eV[20] while PBE
yields a significantly smaller band gap of 1.58 eV. It can also be noted that
unoccupied bands are shifted up in energy while occupied bands are shifted
down. This is in fact a general trend across all the materials in the data set
and it leads to a double peak in the histogram of G0W0 corrections with the
peak of negative (positive) corrections corresponding to occupied (empty)
bands, see Figure 1b. The absolute values of the G0W0 corrections range from 0
to 3 eV with an average value of 1.17 eV, see the histogram in Figure 1c.
Returning to the band diagram in panel (a) we further note that not all the
bands are shifted by the same amount - even when disregarding the different
sign for occupied/empty bands. Although for most materials, all the occupied
bands experience similar, though material specific, shifts and the same holds
for the empty bands, there are several examples, like MoS2, where this is not
the case. Therefore, an accurate prediction of G0W0 corrections for general
bands requires a representation that not only encodes the occupation of the
state, but also information about the energy and shape of the wave function
and its relation to other relevant states of the crystal.

The ENDOME and RAD-PDOS representations, defined in the Methods Section, are
attempts to generalise the notion of the ”local environment” of an atom, which
has been successfully employed to represent solids and molecules in machine
learning studies, to the case of an electronic state. The ENDOME fingerprint
represents the local environment of an energy eigenstate $|nk\rangle$ in terms
of operator matrix elements between the state itself and other eigenstates of
the crystal, $|\langle nk|\hat{A}|n^{\prime}k^{\prime}\rangle|^{2}$. These
matrix elements are arranged on an grid as a function of the energy difference
$\varepsilon_{nk}-\varepsilon_{n^{\prime}k^{\prime}}$, and their sign is used
to encode the occupation of the ”final state” $|n^{\prime}k^{\prime}\rangle$.
With the ENDOME fingerprint two states are thus considered similar if they
have similar matrix elements with other states of similar occupations and
energy displacements. In this work we include matrix elements for the position
operator, momentum operator, and Laplacian operator. Since we exclusively
consider 2D materials in the present work, the fingerprints are split into in-
plane and out-of-plane components for the position operator (labeled $xy$ and
$z$, respectively) and the momentum operator (labeled $p_{\text{xy}}$ and
$p_{\text{z}}$). The RAD-PDOS fingerprint is a correlation function in energy
and radial distance between the atomic orbital projections (onto angular
momentum channels $s,p,$ and $d$) of the reference eigenstate and all other
eigenstates of the crystal. Figure 2 visualises the two types of fingerprints
for three different electronic states of MoS2.

Any reasonable fingerprint should comply with certain general requirements[13]
of which _invariance_ and _simplicity_ are the most fundamental. In the
present context, this means that the fingerprint should be invariant with
respect to the choice of unit cell (number of primitive cells, rotations and
translations), the gauge used for the Bloch wave functions, and that it should
be computationally cheap to generate compared to a full G0W0 calculation. Both
the ENDOME and RAD-PDOS fingerprints clearly fulfill these requirements.
Besides the invariance and simplicity conditions, the fingerprints should also
be _unique_ such that two different systems (here electronic states) are not
mapped to the same fingerprint, and they should be _descriptive_ such that
systems with similar properties are close in fingerprint space. The
interpretation and quantitative assessment of notions such as ”different
systems” and ”similar properties” are obviously problem dependent. This fact
can make it difficult for problem independent fingerprints like the ENDOME and
RAD-PDOS to meet these requirements in general. This is, however, not a
principal problem, and can usually be solved by increasing the size of the
training data set, at least as long as the fingerprints are complex and
flexible enough to capture the variations in the considered systems that are
relevant to the specific learning problem.

An impression of the descriptiveness of the fingerprints can be obtained from
Figure 3, which shows 2-dimensional projections of the ENDOME-$p_{xy}$ and
RAD-PDOS-$dd$ fingerprints using $t$-distributed stochastic neighbor embedding
(tSNE) color coded by the GW corrections. It is clear that data points, which
are close in pxy-space have similar GW corrections. The $dd$ fingerprint is
also descriptive for some data points, but there is also a large blob of data
points that are indistinguishable in fingerprint space but have very different
GW corrections. Not unexpectedly, these points correspond to the subset of
materials without valence $d$-electrons, which results in all-zero $dd$
fingerprint vectors. The tSNE plots for the other components of the ENDOME and
RAD-PDOS fingerprints look similar.

To predict the state-specific G0W0 corrections to the PBE eigenvalues of 2D
semiconductors, we use the XGBoost package [21] to build a machine learning
model based on a gradient boosting algorithm for decision tree ensembles. The
G0W0 data set was described and analysed in detail in Ref. [22]. We split the
data set into a training set of 228 randomly selected materials (37851
electronic states) and a test set consisting of the remaining 58 materials
(8766 electronic states). As objective function we use the mean absolute error
(MAE) between the predicted and actual G0W0 corrections. The electronic states
are represented by the ENDOME and RAD-PDOS fingerprints supplemented by a set
of extra features consisting of the occupation of the state ($f_{nk}=0,1$),
its distance to the Fermi energy ($E_{F}-\varepsilon_{nk}$), the PBE band gap
of the material ($E_{\mathrm{gap}}$), and the static averaged in-plane and
out-of-plane polarisabilities of the material ($\alpha_{x}+\alpha_{y}$ and
$\alpha_{z}$). The averaged in-plane polarisability is used to ensure
invariance of the feature with respect to rotations of the 2D material in the
plane, which is important for materials with in-plane anisotropy. The effect
of including the polarisabilities in the fingerprint has been analysed
separately (see later discussion).

The results of the model together with relevant baselines for assessing its
performance, are summarised in Table 1. The first row shows the estimated
accuracy of our target G0W0 data relative to experiments based on previous
reports in the literature[7, 8, 9]. Experimental data for individual
band/state QP energies are scarce and subject to significant uncertainties,
and thus do not represent a meaningful reference. The remaining rows of the
Table show the mean absolute error (MAE) on the band gap and individual state
energies for different approximate methods versus G0W0. The MAE on state
energies is evaluated over all the bands for which G0W0 data is available,
namely the 8 highest valence bands (VB) and 4 lowest conduction bands (CB).
The second and third rows are straightforward comparisons of band energies
from PBE and HSE06 with G0W0, respectively. The fourth row shows the MAE
between G0W0 and PBE after the occupied and unoccuped PBE energies have been
rigidly shifted (by applying a ”scissors operator”) to match the valence band
maximum (VBM) and conduction band minimum (CBM) of the G0W0 band structure.
From this it follows that the lowest possible MAE on individual band energies
obtainable with a model trained to predict only the VBM and CBM energies, is
0.17 eV. The last two rows of the table shows the MAE on the test set obtained
with the XGBoost model (see below for more details). Improved performance for
the band gap can be obtained by training the model only on the highest valence
and lowest conduction band (last row); however, such a restriction on the
training data reduces the prediction accuracy for bands further away from the
band gap. The numbers marked by (*) refer to the MAE obtained when the static
polarisability of the materials is included in the fingerprint (see later
discussion).

In the following, unless stated otherwise, results refer the case where the
model has been trained on all bands (8VB + 4CB) and with the static
polarizabilities included in the fingerprint.

Figure 4a shows a parity plot of the predicted vs. true values for the train
and test set. The evaluation yields MAEs of 0.05 eV and 0.11 eV for train and
test set, respectively. To test for potential bias of the model, the residual
distributions are plotted in Figure 4b, showing that both the train and test
set have residuals distributed evenly around 0 eV. To estimate the effect of
adding more data to the train set, a learning curve is shown in Figure 4c. The
learning curve is calculated by continuously adding more materials to the
training set while evaluating the performance on a constant test set. The test
set MAE decreases significantly up to $\approx 50$ materials after which the
learning curve flattens considerably, although still presenting a slightly
decreasing MAE. This suggests that a generalizable model can be trained using
a rather limited number of materials, though it should be noted that
overfitting issues decrease with the amount of materials in the training set.
In general, it is difficult to assess whether the learning ability of the
model is limited by the flexibility of the model/fingerprint or by the noise
level in the data set. We do stress, however, that the numerical precision of
the G0W0 corrections is not expected to be much better than 0.05 eV due to
errors introduced by e.g. plane wave extrapolation and linearisation of the
self-energy, see [22]. This could explain (part of) the finite prediction
error of the model.

The data used to train and evaluate the ML model represent states/energies
evaluated at discrete uniformly distributed $k$-points of the Brillouin zone.
However, the resulting ML model can of course be used to predicting the G0W0
energy corrections of states at arbitrary $k$-points and thereby generate
full, densely sampled band structures. Figure 5 shows examples of ML generated
band structures for PtO2, SbClTe, GeS2, and CaCl2, which are all test set
materials. For comparison, the PBE and the true discrete G0W0 energies are
also shown. Overall, the ML bands closely interpolates the true G0W0 energies.
In cases where the ML bands deviate, e.g. the cinduction bands of CaCl2, they
still present a better description than PBE. Interestingly, the ML model is
able to deviate from a scissors operator that would ascribe the same
corrections to all occupied and all unoccupied bands, respectively. This is
for example clear in the PtO2 band structure where the four conduction bands
are shifted by different amounts.

The ML state energies can be translated into ML band gaps by simply
calculating the vertical difference between conduction band minimum and
valence band maximum. Figure 6 shows parity plots of the predicted band gaps
vs. G0W0 band gaps for a ML model trained on all bands and a ML model trained
only on valence and conduction bands. Due to the discreteness of the original
G0W0 data, the ML band gap has been evaluated on the same states (discrete
$k$-points) that define the G0W0 gap. The PBE and HSE06 data are also shown as
baselines. Only data from the test set has been used for the comparison. The
PBE and HSE06 functionals systematically underestimate the band gaps leading
to MAEs of 1.70 eV and 0.85 eV, respectively. The ML model trained on all
bands achieves a MAE on the band gap of 0.18 eV, while training the ML model
only on valence and conduction bands reduces the band gap MAE to 0.15 eV, but
at the cost of increasing the MAE on the individual state energies across all
bands from 0.11 to 0.22 eV.

Methods | Target property MAE
---|---
| Band gap (eV) | State energies (eV)
G0W0 vs. experiment |  $\approx$ 0.3 | N/A
PBE vs. G0W0 | 1.70 | 1.17
HSE06 vs. G0W0 | 0.85 | 0.47
PBE with ideal scissor-operator vs. G0W0 | 0 | 0.17
ML (8VB+4CB) vs. G0W0 | 0.23 ,$0.18^{(*)}$ | 0.14 , $0.11^{(*)}$
ML (VB+CB) vs. G0W0 | 0.18 , $0.15^{(*)}$ | 0.31, $0.22^{(*)}$
Table 1: Summary of results. The table shows the mean absolute error (MAE) on
the band gap and individual state energies for G0W0 versus experiments and
different approximate methods versus G0W0, respectively. The MAE on state
energies is always evaluated for the 8 highest valence bands (VB) and 4 lowest
conduction bands (CB). ML(X) refers to the test set MAE of the gradient
boosting model after training on all bands (8VB+4CB) or only the highest
valence and lowest conduction band (VB+CB), respectively. The values marked by
(*) are obtained after training the model with the static polarisability of
the materials included as extra features in the fingerprint.

##  3 Discussion

Often the evaluation of a machine learning model stops after considering the
overall performance in terms of an objective function like the MAE. However,
important insight may be gained by analysing how the model responds to
different features in the input data. This is particularly important when
devising new types of fingerprints. To extract information about the role of
the different features composing the fingerprint vectors used in the present
work, a feature importance analysis is performed using a feature subset hold-
out method. The features are grouped at two different levels: The first level
has four groups, namely the RAD-PDOS components, the ENDOME components, the
extra features covering the PBE gap, occupation number, distance to the Fermi
level, and finally the in-plane and out-of-plane polarizabilities. The second
level breaks the RAD-PDOS and ENDOME components further down into their
individual $ll^{\prime}$ angular momentum blocks and operator matrix elements,
respectively. The analysis is carried out in two complementary ways where a
group of features is either used exclusively or dropped from the full
fingerprint when training the ML model.

Figure 7 shows the test set MAE on individual state energies for the various
feature groups with the all-feature baseline indicated by the vertical black
line. Focusing first on panel (a), the analysis shows that both the RAD-PDOS
and ENDOME perform well by themselves, though not as well as the full
fingerprint. The extra features, in particular the polarisabilities, are
unable to produce an accurate ML model. The poor performance of the
polarisability-only feature is unsurprising as this feature is fully material
specific and not even able to distinguishing between occupied and unoccupied
states. Panel (b) shows the same analysis when the feature groups are broken
further down. When used alone, the $pp$, $ss$ and $sp$ components of the RAD-
PDOS perform best followed by the various operator matrix elements of the
ENDOME. An interesting observation is that at this level of feature grouping,
almost any group of features can be dropped without increasing the MAE, except
for the in-plane polarisability, $\alpha_{xy}$, which results in a significant
27% increase of the MAE from 0.11 eV to 0.14 eV. This reveals a clear feature
synergy since $\alpha_{xy}$ in itself does not have any predictive ability
unless it is combined with other features (see below). In general, there seems
to be some redundant information in the various fingerprint components since
dropping any of the feature sets, at least at the second level of grouping,
does not affect the test score by much. In some cases, the model might even
gain performance when dropping some features (not visible on the scale of the
plot). This suggests that a feature selection might improve the performance of
the model further.

The role of the $\alpha_{xy}$ feature and its synergy with other features is
further investigated using the general feature importance method SHAP, which
is a game theoretic approach to explain the output of any machine learning
model [23]. SHAP builds an explanation model on top of a ML model which
relates the output from the ML model to the importance of individual features
for each predicted output. The SHAP values for a given feature can thus be
interpreted as the direct effect of that feature on the model output, i.e. the
difference between the model’s prediction when used with and without that
particular feature in the input. Figure 8a shows the SHAP values for
$\alpha_{xy}$ as a function of $\alpha_{xy}$. Only states from the test set
are shown in Figure 8, and the color code in panel (a) reflects the occupancy
of the state. The plot shows a surprisingly clear trend: The SHAP values for
occupied states increase consistently and monotonously for increasing
$\alpha_{xy}$ while the opposite trend is seen for the empty states. In the
following we present a physical explanation for this observation.

The G0W0 correction can be split into two terms with distinctly different
physical origin: $\Delta
E^{\mathrm{QP}}_{nk}=(v_{nk}^{\mathrm{x}}-v^{\mathrm{xc}}_{nk})+\Delta_{nk}^{\mathrm{scr}}$.
The first term (in parenthesis) represents the difference between the local
xc-potential (in this case the PBE potential) and the nonlocal exact exchange
potential while the last term accounts for the interaction of the
electron/hole with its own polarisation cloud. The first term is typically
negative for occupied states and positive for unoccupied states (Hartree-Fock
typically opens the PBE gap), but its magnitude depends on the detailed shape
of the wave functions of the system. In particular, this term can be quite
different for different states of the same material. Moreover, one does not
expect the size of this term to correlate with the material’s static
polarisability and thus it should not be captured by the $\alpha_{xy}$-SHAP
values. The second term is always positive for occupied states (hole
quasiparticles) and negative for unoccupied states (electron quasiparticles)
because the Coulomb interaction of the bare particle with its oppositely
charged polarisation cloud will always stabilise the quasiparticle, thus
shifting occupied states up and empty states down in energy[24, 25, 26]. Now,
the shape and size of the polarisation cloud does not depend on the detailed
shape of the wave function, but is largely governed by the (microscopic)
polarisability of the material. Therefore, on purely physical grounds, the
static macroscopic polarisability, $\alpha_{xy}$, is expected to provide a
good descriptor for $\Delta_{nk}^{\mathrm{scr}}$: A large value of
$\alpha_{xy}$ signals high screening ability of the material and therefore
large QP polarisation clouds, which in turn will yield a large
$\Delta_{nk}^{\mathrm{scr}}$ (with opposite signs for occupied/empty states).
This is exactly what is seen in Figure 8a. By subtracting the
$\alpha_{xy}$-SHAP values for the states at the CBM and VBM, we obtain the
$\alpha_{xy}$-SHAP values for the band gap correction, see Figure 8b. These
show that the $\alpha_{xy}$ feature increases the band gap in materials with
low screening and decreases the band gap in materials with high screening.
Again, this is perfectly in line with the physical understanding of screening-
induced renormalisation of the band gaps[24, 25, 26].

It can be noted that the $\alpha_{xy}$-SHAP values for the state energies and
band gaps are significantly larger than the change in the MAE upon
including/dropping $\alpha_{xy}$ from the feature set, see Figure 7b. For
example, the $\alpha_{xy}$-SHAP values for the band gap range from -0.50 to
1.25 eV while the MAE decreases by 0.03 eV when $\alpha_{xy}$ is included.
This is due to the redundant information carried by the feature set. When the
model is trained without $\alpha_{xy}$ as feature, other features can, to a
large extent, provide the same information. For example, the PBE band gap
alone correlates fairly well with $\alpha_{xy}$, see inset of Figure 8b. To
test this hypothesis, we have carried out the same SHAP analysis for
$E_{g}^{\mathrm{PBE}}$ on a model trained with and without $\alpha_{xy}$ in
the feature set. The analysis shows that when $\alpha_{xy}$ is used to train
the model, the $E_{g}^{\mathrm{PBE}}$-SHAP values are fairly low (below $\pm
0.1$ eV) and do not show any clear trends. In contrast, when $\alpha_{xy}$ is
not included in the fingerprint, the $E_{g}^{\mathrm{PBE}}$-SHAP values are
very similar to the $\alpha_{xy}$-SHAP values shown in Figure 7, although the
values are slightly smaller and the trend less pronounced. This shows that in
the absence of $\alpha_{xy}$ the model uses $E_{g}^{\mathrm{PBE}}$ to encode
similar information. However, the model also finds that $\alpha_{xy}$ provides
a better description of $\Delta_{nk}^{\mathrm{scr}}$ than does
$E_{g}^{\mathrm{PBE}}$, which is why the SHAP values of $E_{g}^{\mathrm{PBE}}$
are dwarfed by those of $\alpha_{xy}$ when both features are available for
learning.

In summary, we have introduced two different methods to generate fingerprints
of individual electronic states based on information available from a standard
DFT ground state calculation (eigenvalues and wave functions). Although the
current work has focused on states in periodic 2D crystals, the methods can be
straightforwardly used to fingerprint states in 3D crystals as well as non-
periodic structures like molecules or surfaces. The fingerprints were used to
train a decision-tree based ML model to predict the G0W0 corrections to the
PBE band structure of a 2D semiconductor. The model achieves a MAE of 0.14 eV
for individual state energies, which is reduced to 0.11 eV when the static
polarisability is included in the fingerprint. For the band gap, the MAE is
0.15-0.23 eV depending on whether the model is trained on all bands or only
the valence/conduction bands and whether or not the static polarisability is
included in the fingerprint. This level of precision is highly encouraging
considering that the noise on the employed G0W0 data for individual state
energies could be on the order of 0.05 eV and that the accuracy of the G0W0
method itself, when evaluated against experimental band gaps, is about 0.3 eV.
Since the bottleneck of the computations is the self-consistent DFT
calculation (in particular the structural relaxation if performed), the method
enables GW-quality band structures at the cost of a DFT calculation.

##  4 Methods

This section describes the definition and generation of the Energy Decomposed
Operator Matrix Elements (ENDOME) and Radially Decomposed Projected Density Of
States (RAD-PDOS) fingerprints. In addition, the G0W0 band structure data set
is presented along with a description of the employed machine learning model.

###  4.1 Electronic state fingerprints

The ENDOME fingerprint is based on operator matrix elements between electronic
states (here assumed to be Bloch states of a periodic crystal)

$\displaystyle A_{nk,n^{\prime}k^{\prime}}=|\langle
nk|\hat{A}|n^{\prime}k\rangle|^{2}$ (1)

where $\hat{A}$ is some operator. For a reference state $|nk\rangle$ with
energy $\varepsilon_{nk}$, the ENDOME fingerprint is defined as

$\displaystyle
m^{A}_{nk}(E)=\sum_{n^{\prime}k^{\prime}}A_{nk,n^{\prime}k^{\prime}}G\left(E-(\varepsilon_{nk}-\varepsilon_{n^{\prime}k});\delta_{E}\right)\exp\left(-\alpha_{E}E\right)\operatorname{sign}(E_{F}-\varepsilon_{n^{\prime}k^{\prime}}),$
(2)

where $G(x;\delta)$ is a Gaussian of width $\delta$ centered at $x=0$. This
function encodes the matrix element between the reference state and all other
states at an energy distance of $E$ from the reference state. In principle,
any operator can be used to create fingerprints, but in this study we include
the position operators ($x,y,z$), the momentum operators
($\nabla_{x},\nabla_{y},\nabla_{z}$), and the Laplace operator ($\nabla^{2}$).
These operators are all diagonal in the $k$ index. In addition, we include the
all-one matrix, $A_{nk,n^{\prime}k^{\prime}}=1$, which essentially yields the
density of states (DOS) translated to the energy of the reference state,
$\varepsilon_{nk}$.

In practice, the function $m^{A}_{nk}(E)$ is represented on a uniformly spaced
energy grid with 50 energy points from -10 to 10 eV around the reference
state. Since we consider 2D materials, the in-plane ($x$ and $y$) components
of both the position and momentum operators are collected into a single
fingerprint vector (i.e. $m^{xy}_{nk}=m^{x}_{nk}+m^{y}_{nk}$ and similarly for
the momentum operator) while the out-of-plane $z$ component is treated
separately. For a given reference state, the ENDOME fingerprint thus consists
of six 50-dimensional vectors resulting in a total of 300 features.

The RAD-PDOS encodes the electronic structure in terms of the density of
states projected onto atomic orbitals. Specifically, a correlation function in
energy and radial distance is defined as

$\displaystyle\rho_{nk}^{\nu\nu^{\prime}}(E,R)=\dfrac{1}{N_{e}}\sum_{n^{\prime}k^{\prime}aa^{\prime}}$
$\displaystyle\rho_{nk}^{a\nu}\rho_{n^{\prime}k^{\prime}}^{a^{\prime}\nu^{\prime}}G\left(R-|R_{a}-R_{a^{\prime}}|;\delta_{R}\right)\exp\left(-\alpha_{R}R\right)G\left(E-(\varepsilon_{nk}-\varepsilon_{n^{\prime}k^{\prime}});\delta_{E}\right)$
(3)
$\displaystyle\times\exp\left(-\alpha_{E}E\right)\operatorname{sign}(E_{F}-\varepsilon_{n^{\prime}k^{\prime}})$

where $N_{e}$ is the number of electrons in the system, $a$ and $a^{\prime}$
denote atoms in the primitive unit cell and the entire crystal, respectively,
and $\nu$ and $\nu^{\prime}$ denote atomic orbitals. The atomic projections
are given by

$\displaystyle\rho_{nk}^{a\nu}=|\langle\psi_{nk}|a\nu\rangle|^{2}$ (4)

The functions $\rho_{nk}^{\nu\nu^{\prime}}(E,R)$ are represented on a uniform
$(E,R)$-grid of size $25\times 20$ spanning the intervals from -10 to 10 eV
(centered around the reference energy $\varepsilon_{nk}$) and 0 to 5 Å,
respectively. For the Gaussian smearing functions we use $\delta_{E}=0.3$ eV
and $\delta_{R}=0.25$ Å, respectively. For a given state, the RAD-PDOS
fingerprint consists of six 2D grids of 500 points each resulting in a total
of 3000 features.

Figure 2 shows examples of ENDOME and RAD-PDOS fingerprints for three
different states at the $K$-point of MoS2. Note that some of the RAD-PDOS
fingerprints are qualitatively similar (e.g. sp and pp) but the scales differ
by about an order of magnitude. This is due to the fact that the density of
states projected onto $s$ and $p$ orbitals have similar dependence on energy.

###  4.2 The G0W0 data set

The data set comprises quasiparticle (QP) energies from 286 G0W0 band
structures of non-magnetic 2D semiconductors covering 14 different crystal
structures and 52 chemical elements. The QP energies have been obtained from
plane-wave-based one-shot G0W0@PBE calculations with full frequency
integration and were produced as a part of the Computational 2D Materials
Database (C2DB)[27, 28]. The data set has been described and analysed in
detail in [22].

The QP energies of the data set have been calculated under the standard
assumption that the G0W0 self-energy can be treated within first-order
perturbation theory and linearised around the non-interacting reference
energy, $\omega=\varepsilon_{nk}$, leading to the expression

$\displaystyle E^{\mathrm{QP}}_{nk}\approx\epsilon_{nk}+\
Z\mathrm{Re}\left[\langle\psi_{nk}|\Sigma(\epsilon_{nk})|\psi_{nk}\rangle\right]$
(5)

where

$\displaystyle
Z=\left(1-\left.\dfrac{\partial\Sigma}{\partial\omega}\right|_{\omega=\epsilon_{nk}}\right)^{-1}$
(6)

is the QP weight and $\psi_{nk}$ is the PBE wave function with eigenvalues
$\epsilon_{nk}$. In practice, the G0W0 correction to the PBE energies, $\Delta
E_{nk}^{\mathrm{QP}}=E^{\mathrm{QP}}_{nk}-\epsilon_{nk}$, were used as targets
for the machine learning model.

To ensure the highest data quality, the original data set was filtered such
that only states with QP weight between 0.7 and 1.0 were kept. As shown in
Ref. [22] the MAE on the QP correction of such states due to the linearisation
of the QP equation is 0.04 eV.

###  4.3 Machine learning model

The machine learning model is built using a gradient boosting method from the
XGBoost distribution based on decision trees in an ensemble [21]. The choice
of XGBoost as learning algorithm is based on its generality and good
performance across multiple machine learning applications and the ability of
training on large amounts of data. For this specific purpose of predicting
G0W0 QP energies, a gaussian process regression method has also been tested
resulting in similar prediction accuracy.

A train and test set is created using a random $80/20\%$ split on the material
level which results in a train set of 228 materials (37851 QP energies) and a
test set of 58 materials (8766 QP energies). Hyperparameters of the learning
algorithm ($\text{max depth}=5$, $\text{learning rate}=0.15$ and $\text{number
of estimators}=60$) (XXX: should we mention what they are? Does this include
hyper parameters defining the fingerprints?)are tuned using a grid search
method with a 5-fold cross-validation of the $80\%$ train set. The performance
of the machine learning is based on the mean absolute error (MAE) of the
$20\%$ test set.

Figure 1: (a) Example of PBE and G0W0 band structures of monolayer MoS2. The
prediction target data is the difference in energy between the PBE and G0W0
energies. (b) Histogram of the G0W0 corrections for all states in all
materials. (c) Histogram of the absolute values of the G0W0 corrections with a
mean of 1.17 eV. Figure 2: Visualization of electronic state fingerprints for
MoS2. Top left panel shows the PBE band structure. Top right panels show
ENDOME fingerprints of the conduction band minimum and valence band maximum
states for the K-point. The panels below show six RAD-PDOS fingerprints for
combinations of s, p and d orbitals.  Figure 3: a) tSNE components of ENDOME
pxy. b) tSNE components of RAD-PDOS dd fingerprints color-coded with the GW
corrections. For pxy, states with similar GW corrections are also close in
fingerprint space. In b) a large amount of the states with both positive and
negative GW corrections have similar distances in fingerprint space,
corresponding to the materials without d-electrons where the RAD-PDOS dd
fingerprint will be all zeros. Figure 4: Left: Parity plot showing the ML
predicted vs. true values of the GW correction for individual states for the
train and test set. The MAE’s of the train and test set are 0.05 and 0.11 eV,
respectively. Mid: Histograms of the prediction residuals of the train and
test set. Right: Learning curve for the ML model showing validation MAE as
function of number of materials/states in training set. Figure 5: Examples of
band structures for four materials from the test set. Both PBE and GW band
structures are shown along with the ML predictions.  Figure 6: Parity plots
for predicted bandgaps vs. GW bandgaps for PBE and HSE06 and two different ML
models predicting GW corrections for either all bands (MAE = 0.18 eV) or only
valence and conduction bands (MAE = 0.15 eV) which significantly outperform
PBE and HSE06 with MAE’s of 1.70 and 0.85 eV, respectively.  Figure 7: Feature
analysis of ML model. Solid bars refers to a ML model using only the specific
features while the shaded bars are for a ML model without these features.
Left: High-level feature groups. Right: Low-level feature groups.  Figure 8:
SHAP values for $\alpha_{xy}$ for the prediction of GW correction energies
colorcoded by occupancy. For materials with a low polarisability the ML model
predicts a more negative GW correction for the occupied states and a more
positive correction for the unoccpied states. For materials with a high
polarisability the occupied states are predicted with a more positive
correction when using the polarisability as a feature while the unoccupied
states are only weakly affected.

##  5 Data and code availability

The structures of the materials used for this study can be found in the
C2DB[29]. The dataset generated for this study is available at
https://gitlab.com/knosgaard/electronic-structure-fingerprints along with the
Python code used to compute the fingerprints.

## References

  * [1] Walter Kohn and Lu Jeu Sham.  Self-consistent equations including exchange and correlation effects.  Physical review, 140(4A):A1133, 1965.
  * [2] John P Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Phys. Rev. Lett., 77(18):3865, 1996.
  * [3] RW Godby, M Schlüter, and LJ Sham.  Accurate exchange-correlation potential for silicon and its discontinuity on addition of an electron.  Phys. Rev. Lett., 56(22):2415, 1986.
  * [4] Lars Hedin.  New method for calculating the one-particle green’s function with application to the electron-gas problem.  Phys. Rev., 139(3A):A796, 1965.
  * [5] Dorothea Golze, Marc Dvorak, and Patrick Rinke.  The GW compendium: A practical guide to theoretical photoemission spectroscopy.  Front. Chem., 7, 2019.
  * [6] Ferdi Aryasetiawan and Olle Gunnarsson.  The GW method.  Rep. Prog. Phys., 61(3):237, 1998.
  * [7] Falco Hüser, Thomas Olsen, and Kristian S Thygesen.  Quasiparticle gw calculations for solids, molecules, and two-dimensional materials.  Phys. Rev. B Condens. Matter, 87(23):235132, 2013.
  * [8] M Shishkin and G Kresse.  Self-consistent GW calculations for semiconductors and insulators.  Phys. Rev. B Condens. Matter, 75(23):235102, 2007.
  * [9] Dmitrii Nabok, Andris Gulans, and Claudia Draxl.  Accurate all-electron G0W0 quasiparticle energies employing the full-potential augmented plane-wave method.  Phys. Rev. B Condens. Matter, 94(3):035118, 2016.
  * [10] Logan Ward, Ankit Agrawal, Alok Choudhary, and Christopher Wolverton.  A general-purpose machine learning framework for predicting properties of inorganic materials.  npj Computational Materials, 2(1):1–7, 2016.
  * [11] Luca M Ghiringhelli, Jan Vybiral, Sergey V Levchenko, Claudia Draxl, and Matthias Scheffler.  Big data of materials science: critical role of the descriptor.  Physical review letters, 114(10):105503, 2015.
  * [12] Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller, and O Anatole Von Lilienfeld.  Fast and accurate modeling of molecular atomization energies with machine learning.  Physical review letters, 108(5):058301, 2012.
  * [13] Felix Faber, Alexander Lindmaa, O. Anatole von Lilienfeld, and Rickard Armiento.  Crystal structure representations for machine learning models of formation energies.  International Journal of Quantum Chemistry, 115(16):1094–1101, 2015\.
  * [14] Haoyan Huo and Matthias Rupp.  Unified representation of molecules and crystals for machine learning, 2018.
  * [15] Peter Bjørn Jørgensen, Murat Mesta, Suranjan Shil, Juan Maria García Lastra, Karsten Wedel Jacobsen, Kristian Sommer Thygesen, and Mikkel N Schmidt.  Machine learning-based screening of complex molecules for polymer solar cells.  The Journal of chemical physics, 148(24):241735, 2018.
  * [16] Ya Zhuo, Aria Mansouri Tehrani, and Jakoah Brgoch.  Predicting the band gaps of inorganic solids by machine learning.  The Journal of Physical Chemistry Letters, 9(7):1668–1673, 2018\.  PMID: 29532658.
  * [17] Arunkumar Chitteth Rajan, Avanish Mishra, Swanti Satsangi, Rishabh Vaish, Hiroshi Mizuseki, Kwang-Ryeol Lee, and Abhishek K Singh.  Machine-learning-assisted accurate band gap predictions of functionalized mxene.  Chemistry of Materials, 30(12):4031–4038, 2018.
  * [18] Jiechun Liang and Xi Zhu.  Phillips-inspired machine learning for band gap and exciton binding energy prediction.  The journal of physical chemistry letters, 10(18):5640–5646, 2019\.
  * [19] Sandip De, Albert P Bartók, Gábor Csányi, and Michele Ceriotti.  Comparing molecules and solids across structural and alchemical space.  Physical Chemistry Chemical Physics, 18(20):13754–13769, 2016.
  * [20] AR Klots, AKM Newaz, Bin Wang, D Prasai, H Krzyzanowska, Junhao Lin, D Caudel, NJ Ghimire, J Yan, BL Ivanov, et al.  Probing excitonic states in suspended two-dimensional semiconductors by photocurrent spectroscopy.  Scientific Reports, 4(1):1–7, 2014.
  * [21] Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.
  * [22] Asbjørn Rasmussen, Thorsten Deilmann, and Kristian S. Thygesen.  Towards fully automated gw band structure calculations: What we can learn from 60.000 self-energy evaluations.  npj Computational Materials, 7(1):22, Jan 2021.
  * [23] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 4765–4774. Curran Associates, Inc., 2017.
  * [24] Jeffrey B Neaton, Mark S Hybertsen, and Steven G Louie.  Renormalization of molecular electronic levels at metal-molecule interfaces.  Physical Review Letters, 97(21):216405, 2006.
  * [25] Juan Maria Garcia-Lastra, Carsten Rostgaard, Angel Rubio, and Kristian Sommer Thygesen.  Polarization-induced renormalization of molecular levels at metallic and semiconducting surfaces.  Physical Review B, 80(24):245427, 2009.
  * [26] Per S Schmidt, Christopher E Patrick, and Kristian S Thygesen.  Simple vertex correction improves g w band energies of bulk and two-dimensional crystals.  Physical Review B, 96(20):205206, 2017.
  * [27] Sten Haastrup, Mikkel Strange, Mohnish Pandey, Thorsten Deilmann, Per S Schmidt, Nicki F Hinsche, Morten N Gjerding, Daniele Torelli, Peter M Larsen, Anders C Riis-Jensen, et al.  The computational 2d materials database: high-throughput modeling and discovery of atomically thin crystals.  2D Materials, 5(4):042002, 2018.
  * [28] MN Gjerding, A Taghizadeh, A Rasmussen, S Ali, F Bertoldo, T Deilmann, UP Holguin, NR Knøsgaard, M Kruse, S Manti, et al.  Recent progress of the computational 2d materials database (c2db).  arXiv:2102.03029, 2021.
  * [29] Computational 2D Materials Database (C2DB).  https://cmr.fysik.dtu.dk/c2db/c2db.html.  Accessed: 2021-07-01.


# Role of the scalar $f_{0}(980)$ in the process
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$

Han Zhang  School of Physics and Microelectronics, Zhengzhou University,
Zhengzhou, Henan 450001, China  Yun-He Lyu  School of Physics and
Microelectronics, Zhengzhou University, Zhengzhou, Henan 450001, China  Li-
Juan Liu<EMAIL_ADDRESS>School of Physics and Microelectronics,
Zhengzhou University, Zhengzhou, Henan 450001, China  En Wang
<EMAIL_ADDRESS>School of Physics and Microelectronics, Zhengzhou
University, Zhengzhou, Henan 450001, China

###### Abstract

Based on the BESIII measurements on the reaction of
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$, we investigate this process by
considering the $S$-wave pseudoscalar-pseudoscalar interaction within the
unitary chiral approach, and the contributions from the intermediate
resonances $f_{0}(1370)$ and $f_{2}(1270)$. Our calculation could reasonably
reproduce the experimental data, and our results imply that the $f_{0}(980)$,
dynamically generated from the $S$-wave pseudoscalar-pseudoscalar interaction,
plays an important role in this process, and the contributions from the
intermediate resonances $f_{0}(1370)$ and $f_{2}(1270)$ are also necessary.
The more precise measurements of this process in future could shed light on
the nature of the $f_{0}(1370)$ and $f_{2}(1270)$.

##  I Introduction

In last decades, the constituent quark model successfully explains the
composition of the most mesons [1]. However, the properties of the light
scalar mesons are difficult to be described within the constituent quark
model. For the scalar mesons with masses less than 1 GeV, $f_{0}(500)$,
$a_{0}(980)$, $f_{0}(980)$, and $K(700)$, there are different explanations,
such as the $q\bar{q}$ states, tetraquark states, or molecular states [2, 3,
4, 5]. In addition, the nature of the isospin-zero scalar mesons with larger
masses, $f_{0}(1370)$, $f_{0}(1500)$, and $f_{0}(1710)$, are also difficult to
be fully understood [1]. Recently, the BESIII and Belle/Belle II have
accumulated lots of experimental information about the charmed hadrons decays,
which provides an important platform to explore the internal structures of
those scalar mesons [6, 9, 7, 8].

For instance, the BESIII Collaboration have observed the scalar $S(1710)$111It
should be stressed that BESIII does not distinguish between the $a_{0}(1710)$
and $f_{0}(1710)$ in the $D^{+}_{s}\to K^{0}_{S}K^{0}_{S}\pi^{+}$, and denotes
the combined state as $S(1710)$ [10]. with the mass of $1723\pm 11\pm 2$ MeV
and width of $140\pm 14\pm 4$ MeV in the process of $D_{s}^{+}\to
K^{0}_{S}K^{0}_{S}\pi^{+}$ [10]. Later, the similar structure was also
observed in the process of $D_{s}^{+}\to K^{+}K^{0}_{S}\pi^{0}$ by the BESIII
Collaboration, associated with a new scalar state $a_{0}(1817)$ with the mass
of $1817\pm 8\pm 20$ MeV and width of $97\pm 22\pm 15$ MeV [11]. The $\sim
100$ MeV mass discrepancy between $S(1710)$ and $a_{0}(1817)$ causes the doubt
whether they are isospin partners. However, the line-shapes of the
$a_{0}(1817)$ and $S(1710)$ in the $K^{+}K^{0}_{S}$ and $K^{0}_{S}K^{0}_{S}$
invariant mass distributions are so similar that the measured mass discrepancy
of the $a_{0}(1817)$ and $S(1710)$ may be due to their peaks’ being very close
to the threshold of the $K\bar{K}$ mass spectrum. We have investigated both
reactions, and found that the BESIII measurements of both reactions could be
well reproduced by regarding the intermediate state $a_{0}(1817)$ as the
$K^{*}\bar{K}^{*}$ molecular state [12, 13], which is also supported by the
study of Ref. [14].

The BESIII Collaboration recently performed an amplitude analysis of the
process $D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$, and determined the branching
fraction $\mathcal{B}$($D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$) = (2.8 $\pm$ 0.4
$\pm$ 0.4)$\%$ [15], which improves the precision with a factor of two
compared to the results of the CLEO Collaboration [16]. The significant
$f_{0}(980)$ signal and an enhancement structure around 1300 MeV were observed
in the $\pi^{0}\pi^{0}$ invariant mass spectrum, and the later structure could
be associated with the intermediate resonances $f_{0}(1370)$ and
$f_{2}(1270)$. The scalar $f_{0}(980)$ could be explained as the $K\bar{K}$
molecular state dynamically generated from the $S$-wave pseudoscalar-
pseudoscalar interaction [17, 18], which was supported by many works [6, 9, 7,
8]. In addition, the resonance $f_{0}(1370)$ is a broad state, and its mass
and width are not well established. There are various explanations for the
nature of $f_{0}(1370)$, such as the $n\bar{n}$ state, molecular state [19,
20, 21], $s\bar{s}$ [22], and glueball [24, 23]. In Refs. [19, 20, 21],
$f_{2}(1270)$ is also explained as the bound system of two vector mesons,
which is challenged in Ref. [25]. Therefore, the study of the process
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$ should be helpful to deepen our
understanding about the nature of the $f_{0}(980)$, $f_{0}(1370)$, and
$f_{2}(1270)$.

In this work, we will investigate the process
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$ by considering the scalar $f_{0}(980)$
generated from the $S$-wave pseudoscalar-pseudoscalar interaction within the
chiral unitary approach and the contributions from the intermediate resonances
$f_{0}(1370)$ and $f_{2}(1270)$.

This paper is organized as follows. In Sec. II, we will present the mechanism
for the process $D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$, and results and
discussions will be given in Sec. III, followed by a short summary in the last
section.

##  II Formalism

In this section, we will present the mechanism of the process
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$. This process will happen via three steps,
the weakly decay, the hadronization, and the final state interactions [9, 7,
12, 13]. Firstly, the $c$ quark of the initial $D_{s}^{+}$ weakly decays into
an $s$ quark and a $W^{+}$ boson, and subsequently the $W^{+}$ boson decays
into a $\bar{d}u$ quark pair. Then all the quarks, together with the quark
pair $q\bar{q}$ (=$u\bar{u}$ \+ $d\bar{d}$ \+ $s\bar{s}$ ) created from the
vacuum with the quantum numbers $J^{PC}=0^{++}$, hadronize into the hadrons,
which could be classified as the $W^{+}$ internal emission of Figs. 1(a-b) and
the $W^{+}$ external emission of Fig. 1(c-d).

For the $W^{+}$ internal emissions of Figs. 1(a) and 1(b), the $s\bar{d}$ or
$u\bar{s}$ hadronize into the $\bar{K}^{0}$ or $K^{+}$, and the hadronization
of the other quarks could be expressed as

$\displaystyle\sum_{i}u(\bar{q}_{i}q_{i})\bar{s}=\sum_{i}M_{1i}M_{i3}=(M^{2})_{13},$
(1)
$\displaystyle\sum_{i}s(\bar{q}_{i}q_{i})\bar{d}=\sum_{i}M_{3i}M_{i2}=(M^{2})_{32},$
(2)

where $i=1,~{}2,~{}3$ correspond to the $u$, $d$, and $s$ quarks,
respectively, and the $M$ is the $q\bar{q}$ matrix,

$\displaystyle M=\left(\begin{array}[]{ccc}u\bar{u}&u\bar{d}&u\bar{s}\\\
d\bar{u}&d\bar{d}&d\bar{s}\\\ s\bar{u}&s\bar{d}&s\bar{s}\end{array}\right)\,.$
(6)

Within the $SU(3)$ flavor symmetry, the matrix $M$ can be written in terms of
pseudoscalar mesons as [7]

$\displaystyle
M=\left(\begin{array}[]{ccc}\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}+\frac{\eta^{\prime}}{\sqrt{6}}&\pi^{+}&K^{+}\\\
\pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}+\frac{\eta^{\prime}}{\sqrt{6}}&K^{0}\\\
K^{-}&\bar{K}^{0}&-\frac{\eta}{\sqrt{3}}+\frac{2\eta^{\prime}}{\sqrt{6}}\end{array}\right)\,.$
(10)

Figure 1: The quark level diagrams for the $D_{s}^{+}$ weak decays. (a) The
internal emission of $D_{s}^{+}\to\bar{K}^{0}u\bar{s}$ and hadronization of
the $u\bar{s}$ through $q\bar{q}$ with vacuum quantum numbers. (b) The
internal emission of $D_{s}^{+}\to K^{+}s\bar{d}$ and hadronization of the
$s\bar{d}$ through $q\bar{q}$ with vacuum quantum numbers. (c) The external
emission of $D_{s}^{+}\to\pi^{+}s\bar{s}$ and hadronization of the $s\bar{s}$
through $q\bar{q}$ with vacuum quantum numbers. (d) The external emission of
$D_{s}^{+}\to u\bar{d}\eta$ and hadronization of the $u\bar{d}$ through
$q\bar{q}$ with vacuum quantum numbers.

Since the $\eta^{\prime}$ has a large mass and does not play a role in the
generation of $f_{0}(980)$, we ignore the $\eta^{\prime}$ component in this
work. Then the Eqs. (1-2) can be rewritten as,

$\displaystyle(M^{2})_{13}=\frac{1}{\sqrt{2}}K^{+}\pi^{0}+\pi^{+}K^{0},$ (11)
$\displaystyle(M^{2})_{32}=K^{-}\pi^{+}-\frac{1}{\sqrt{2}}\bar{K^{0}}\pi^{0}.$
(12)

For the $W^{+}$ external emission of Fig. 1(c), the quark $u\bar{d}$ of the
$W^{+}$ decay hadroinze into $\pi^{+}$, and the $s\bar{s}$ pair, together with
the created $q\bar{q}$ pair, hadronize into the states as,

$\displaystyle\sum_{i}s(\bar{q}_{i}q_{i})\bar{s}$ $\displaystyle=$
$\displaystyle\sum_{i}M_{3i}M_{i3}=(M^{2})_{33}$ (13) $\displaystyle=$
$\displaystyle K^{+}K^{-}+K^{0}\bar{K}^{0}+\frac{1}{3}\eta\eta.$

For the $W^{+}$ external emission of Fig. 1(d), the $s\bar{s}$ pair could
hadronize into the $\eta$ meson, and the quark $u\bar{d}$ of the $W^{+}$
decay, together with the created $q\bar{q}$ pair, hadroinze into
$\pi^{+}\eta$, which contributes to the process $D_{s}^{+}\to\pi^{+}\eta\eta$.
Thus, we have,,

$\displaystyle\sum_{i}u(\bar{q}_{i}q_{i})\bar{d}(s\bar{s})$ $\displaystyle=$
$\displaystyle\sum_{i}M_{1i}M_{i2}M_{33}$ (14) $\displaystyle=$
$\displaystyle(M^{2})_{12}\left(-\frac{1}{\sqrt{3}}\eta\right)=-\frac{2}{3}\pi^{+}\eta\eta.$

Then the processes of the $D^{+}_{s}$ decaying into all possible states could
be expressed as

$\displaystyle H^{(a)}$ $\displaystyle=$ $\displaystyle
V_{cs}V_{ud}\left(\frac{1}{\sqrt{2}}\pi^{0}K^{+}+\pi^{+}K^{0}\right)\bar{K}^{0},$
(15) $\displaystyle H^{(b)}$ $\displaystyle=$ $\displaystyle
V_{cs}V_{ud}\left(\pi^{+}K^{-}-\frac{1}{\sqrt{2}}\pi^{0}\bar{K}^{0}\right)K^{+},$
(16) $\displaystyle H^{(c)}$ $\displaystyle=$ $\displaystyle
CV_{cs}V_{ud}\left(K^{+}K^{-}+K^{0}\bar{K^{0}}+\frac{1}{3}\eta\eta\right)\pi^{+},$
(17) $\displaystyle H^{(d)}$ $\displaystyle=$ $\displaystyle
CV_{cs}V_{ud}\left(-\frac{2}{3}\eta\eta\right)\pi^{+},$ (18)

where $V_{ud}$ and $V_{cs}$ are the CKM matrix elements. Since the external
emission of $W^{+}$ bosons is color-favored relative to the $W^{+}$ internal
emission, an extra color factor $C$ can be introduced to account for the
relative weight of the $W^{+}$ external emission with respect to the $W^{+}$
internal emission. For the $W^{+}$ external emission, the $u\bar{d}$ quark
pair from the $W^{+}$ decay can form the color singlet $\pi^{+}$, and the $u$
and $\bar{d}$ have three choices of colors, while for the $W^{+}$ internal
emission, the $u$, $\bar{d}$, $s$ and $\bar{s}$ quarks from the $W^{+}$ decay
have fixed colors. Thus, the factor $C$ is taken to be 3 in this work [26, 27,
28, 29]. Now we have all the possible components after the
hadronization,222Here we neglect the components $\pi^{0}\bar{K}^{0}K^{+}$.
Indeed, the interaction of the $\bar{K}^{0}K^{+}$ to $\pi^{+}\pi^{0}$ can be
given by the intermediate resonances, which is shown very small by BESIII
[15].

$\displaystyle H=$ $\displaystyle H^{(a)}+H^{(b)}+H^{(c)}+H^{(d)}$
$\displaystyle=$ $\displaystyle
V_{p}V_{cs}V_{ud}\left\\{C(K^{+}K^{-}+K^{0}\bar{K}^{0}+\frac{1}{3}\eta\eta)\pi^{+}\right.$
$\displaystyle\left.+\bar{K^{0}}K^{0}\pi^{+}+K^{+}K^{-}\pi^{+}\right\\}$
$\displaystyle=$ $\displaystyle
V_{p}V_{cs}V_{ud}\left\\{(C+1)\left(K^{+}K^{-}\pi^{+}+K^{0}\bar{K}^{0}\pi^{+}\right)\right.$
$\displaystyle\left.-\frac{C}{3}\eta\eta\pi^{+}\right\\},$ (19)

where the $V_{p}$ is the factors of the production vertices containing all the
dynamics. After the preliminary weakly decay, the meson pairs of $K^{+}K^{-}$,
$K^{0}\bar{K}^{0}$, and $\eta\eta$ could undergo the $S$-wave final state
interaction to give rise to the $\pi^{0}\pi^{0}$ final state, where the scalar
meson $f_{0}(980)$ could be dynamically generated as shown in Fig. 2.

Figure 2: The S-wave final state interaction of the meson pairs.

The amplitude of the $S$-wave pseudoscalar-pseudoscalar interaction,
generating the scalar $f_{0}(980)$, now can be written as,

$\displaystyle{\mathcal{M}}_{f_{0}(980)}$ $\displaystyle=$ $\displaystyle
V_{p}\left[(C+1)G_{K^{+}K^{-}}t_{K^{+}K^{-}\to\pi^{0}\pi^{0}}\right.$ (20)
$\displaystyle+(C+1)G_{K^{0}\bar{K}^{0}}t_{K^{0}\bar{K}^{0}\to\pi^{0}\pi^{0}}$
$\displaystyle\left.-\frac{C}{3}G_{\eta\eta}t_{\eta\eta\to\pi^{0}\pi^{0}}\right],$

where $G_{i}$ is the loop function of two-meson propagator, and $t_{i\to j}$
is the transition amplitude of the $i$-channel to $j$-channel, both of which
are the functions of the $\pi^{0}\pi^{0}$ invariant mass $M_{\pi^{0}\pi^{0}}$.
The loop function is given by,

$G_{i}=i\int\frac{d^{4}q}{(2\pi)^{4}}\,\frac{1}{(P-q)^{2}-m_{1}^{2}+i\epsilon}\,\frac{1}{q^{2}-m^{2}_{2}+i\epsilon},$
(21)

where $m_{1}$ and $m_{2}$ are the masses of the two mesons in the loop of the
$i$-channel, and $P$ and $q$ are the four-momenta of the two-meson system and
the second meson, respectively. The Mandelstam invariant
$s=P^{2}=M_{\pi^{0}\pi^{0}}^{2}$. The loop function of Eq. (21) is
logarithmically divergent, and there are two methods to solve this singular
integral, either using the three-momentum cut-off method, or the dimensional
regularization method. The choice of a particular regularization scheme does
not, of course, affect our argumentation. In this work, we performed the
integral for $q$ in Eq. (21) with a cut-off $|\vec{q}_{\rm max}|=600$ MeV [31,
30]. The transition amplitude $t_{i\to j}$ can be obtained by solving the
Bethe-Salpeter equation in coupled channels,

$T=[1-VG]^{-1}V,$ (22)

where $V$ is a $5\times 5$ matrix of the interaction kernel, we take five
channels $\pi^{+}\pi^{-}$, $\pi^{0}\pi^{0}$, $K^{+}K^{-}$, $K^{0}\bar{K}^{0}$,
and $\eta\eta$. The explicit expressions of the $5\times 5$ matrix elements in
$S$-wave are given by [32, 30, 31, 33]

$\displaystyle
V_{11}=-\frac{1}{2f^{2}}s,~{}~{}~{}V_{12}=-\frac{1}{\sqrt{2}f^{2}}(s-m^{2}_{\pi}),$
$\displaystyle V_{13}=-\frac{1}{4f^{2}}s,~{}~{}V_{14}=-\frac{1}{4f^{2}}s,$
$\displaystyle
V_{15}=-\frac{1}{3\sqrt{2}f^{2}}m^{2}_{\pi},~{}~{}~{}V_{22}=-\frac{1}{2f^{2}}m^{2}_{\pi},$
$\displaystyle
V_{23}=-\frac{1}{4\sqrt{2}f^{2}}s,~{}~{}~{}V_{24}=-\frac{1}{4\sqrt{2}f^{2}}s,$
$\displaystyle
V_{25}=-\frac{1}{6f^{2}}m^{2}_{\pi},~{}~{}V_{33}=-\frac{1}{2f^{2}}s,~{}~{}~{}V_{34}=-\frac{1}{4f^{2}}s,$
$\displaystyle
V_{35}=-\frac{1}{12\sqrt{2}f^{2}}(9s-6m^{2}_{\eta}-2m^{2}_{\pi}),$
$\displaystyle
V_{44}=-\frac{1}{2f^{2}}s,~{}~{}V_{45}=-\frac{1}{12\sqrt{2}f^{2}}(9s-6m^{2}_{\eta}-2m^{2}_{\pi}),$
$\displaystyle V_{55}=-\frac{1}{18f^{2}}(16m^{2}_{K}-7m^{2}_{\pi}),$ (23)

where $f=93$ MeV is the pion decay constant, and $s$ is invariant mass square
of the meson-meson system, $m_{\pi}$, $m_{K}$, and $m_{\eta}$ are the masses
of the pion, kaon, and $\eta$ mesons, respectively [1]. The unitary
normalization $|\eta\eta>\to\frac{1}{\sqrt{2}}|\eta\eta>$ and
$|\pi^{0}\pi^{0}>\to\frac{1}{\sqrt{2}}|\pi^{0}\pi^{0}>$ has been taken easily
account for the identify of the particle when using the loop function $G$
without an extra factor [30].

In addition to the scalar $f_{0}(980)$, the BESIII has observed the
enhancement structure around 1300 MeV in the $\pi^{0}\pi^{0}$ invariant mass
distribution, which could be associated with the resonances $f_{0}(1370)$ and
$f_{2}(1270)$. Hence, we will also take into account the contributions from
the intermediate resonances $f_{0}(1370)$ and $f_{2}(1270)$.

For the contribution of $f_{0}(1370)$ in the decay of
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$, we describe it using the Breit-Winger
form,

$\displaystyle\mathcal{M}_{f_{0}(1370)}=\frac{\alpha\times
m^{2}_{f_{0}(1370)}}{m^{2}_{\pi^{0}\pi^{0}}-M^{2}_{f_{0}(1370)}+i\Gamma_{f_{0}(1370)}M_{f_{0}(1370)}},$
(24)

where $\alpha$ is the strength of the $f_{0}(1370)$. Considering that the mass
and width of $f_{0}(1370)$ have large uncertainties [1], we fix its mass to be
1300 MeV, the center position of the enhancement structure in the
$\pi^{0}\pi^{0}$ invariant mass distribution of BESIII measurements [15], and
take the width as a free parameter.

Taking into account that the $f_{2}(1270)$ couples to the $\pi^{0}\pi^{0}$ in
$D$-wave, we can write the contribution of this resonance as follows [34],

$\displaystyle\mathcal{M}_{f_{2}(1270)}=\frac{\beta\times(3{\rm
cos}^{2}\theta-1)\times{\tilde{p}_{\pi^{0}}^{2}}}{M^{2}_{\pi^{0}\pi^{0}}-M^{2}_{f_{2}(1270)}+i{M_{f_{2}(1270)}}\Gamma_{f_{2}(1270)}},$
(25)

where $\beta$ stands for the strength of the $D$-wave amplitudes. The mass and
width of the $f_{2}(1270)$ are $M_{f_{2}(1270)}=1275.5$ MeV and
$\Gamma_{f_{2}(1270)}=186.7$ GeV, taken from the Review of Particle Physics
[1]. $\tilde{p}_{\pi^{0}}$ is the momentum of $\pi^{0}$ in the
$\pi^{0}\pi^{0}$ rest frame,

$\displaystyle\tilde{p}_{\pi^{0}}=\frac{\lambda^{1/2}\left(M^{2}_{\pi^{0}\pi^{0}},m^{2}_{\pi^{0}},M^{2}_{\pi^{0}}\right)}{2M_{\pi^{0}\pi^{0}}}.$
(26)

The parameter $\theta$ is the angle between the momentum of $\pi^{0}$ and
$\pi^{+}$ in the rest frame of the $\pi^{0}\pi^{0}$ system [34],

$\displaystyle{\rm
cos}\theta=\frac{M^{2}_{\pi^{+}\pi^{0}}-M^{2}_{D^{+}_{s}}-M^{2}_{\pi^{0}}+2\tilde{p}^{0}_{D^{+}_{s}}\tilde{p}^{0}_{\pi^{0}}}{2\tilde{p}_{\pi^{+}}\tilde{p}_{\pi^{0}}},$
(27)

where $\tilde{p}^{0}_{D_{s}^{+}}$ ($\tilde{p}^{0}_{\pi^{0}}$) is the energy of
$D_{s}^{+}$ ($\pi^{0}$) in the $\pi^{0}\pi^{0}$ rest frame, and
$\tilde{p}_{\pi^{+}}=\tilde{p}_{D_{s}^{+}}$ is the $\pi^{+}$ ($D_{s}^{+}$)
momentum in this same frame. We give the explicit forms for those variables
below,

$\displaystyle\tilde{p}_{\pi^{+}}$ $\displaystyle=$
$\displaystyle\tilde{p}_{D^{+}_{s}}=\frac{\lambda^{1/2}\left(M^{2}_{D^{+}_{s}},M^{2}_{\pi^{0}\pi^{0}},m^{2}_{\pi^{+}}\right)}{2M_{\pi^{0}\pi^{0}}},$
$\displaystyle\tilde{p}^{0}_{D^{+}_{s}}$ $\displaystyle=$
$\displaystyle\sqrt{M^{2}_{D^{+}_{s}}+\tilde{p}^{2}_{D^{+}_{s}}},$
$\displaystyle\tilde{p}^{0}_{\pi^{0}}$ $\displaystyle=$
$\displaystyle\frac{M_{\pi^{0}\pi^{0}}}{2},$ (28)

where $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2yz-2xz$.

Consequently, the total amplitude of $D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$ can
be described by

$\displaystyle{\mathcal{M}}$ $\displaystyle=$ $\displaystyle a_{\rm
bg}+{\mathcal{M}}_{f_{0}(980)}+{\mathcal{M}}_{f_{0}(1370)}+{\mathcal{M}}_{f_{2}(1270)},$
(29)

where the constant $a_{\rm bg}$ is the background contribution333In Ref. [15],
BESIII has considered the background contributions from the
$D_{s}^{+}\to\pi^{+}\pi^{0}\eta$ where $\pi^{0}\eta$ is misreconstructed as
$\pi^{0}\pi^{0}$. In addition, there are also the background from the non-
resonant contributions.. As a result, the amplitude of Eq. (29) depends on the
two independent invariant masses $M_{\pi^{0}\pi^{0}}$ and
$M_{\pi^{+}\pi^{0}}$, the double differential width for the process
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$ is given by,

$\frac{d^{2}\Gamma}{dM_{\pi^{0}\pi^{0}}dM_{\pi^{+}\pi^{0}}}=\frac{M_{\pi^{0}\pi^{0}}M_{\pi^{+}\pi^{0}}}{(2\pi)^{3}8m^{3}_{D^{+}_{s}}}|\mathcal{M}|^{2}.$
(30)

One can obtain the ${d\Gamma}/{dM_{\pi^{0}\pi^{0}}}$ and
${d\Gamma}/{dM_{\pi^{+}\pi^{0}}}$ by integrating the Eq. (30) over the other
invariant mass variable with relations as follows,

$\frac{d\Gamma}{dM_{\pi^{0}\pi^{0}}}=\int\frac{M_{\pi^{0}\pi^{0}}M_{\pi^{+}\pi^{0}}}{(2\pi)^{3}8m^{3}_{D^{+}_{s}}}|\mathcal{M}|^{2}dM_{\pi^{+}\pi^{0}}.$
(31)

With a given $M_{\pi^{0}\pi^{0}}$, the upper and lower bounds of the
$M_{\pi^{+}\pi^{0}}$ are,

$\displaystyle(M^{2}_{\pi^{+}\pi^{0}})_{\rm max}$ $\displaystyle=$
$\displaystyle\left(E^{*}_{\pi^{+}}+E^{*}_{\pi^{0}}\right)^{2}-$ (32)
$\displaystyle\left(\sqrt{E^{*2}_{\pi^{+}}-m^{*2}_{\pi^{+}}}-\sqrt{E^{*2}_{\pi^{0}}-m^{*2}_{\pi^{0}}}\right)^{2},$
$\displaystyle(M^{2}_{\pi^{+}\pi^{0}})_{\rm min}$ $\displaystyle=$
$\displaystyle\left(E^{*}_{\pi^{+}}+E^{*}_{\pi^{0}}\right)^{2}-$ (33)
$\displaystyle\left(\sqrt{E^{*2}_{\pi^{+}}-m^{*2}_{\pi^{+}}}+\sqrt{E^{*2}_{\pi^{0}}-m^{*2}_{\pi^{0}}}\right)^{2},$

here $E_{\pi^{+}}^{*}$ and $E^{*}_{\pi^{0}}$ are the energies of $\pi^{+}$ and
$\pi^{0}$ in the $\pi^{0}\pi^{0}$ rest frame, respectively,

$\displaystyle E^{*}_{\pi^{+}}$ $\displaystyle=$
$\displaystyle\frac{m^{2}_{D^{+}_{s}}-M^{2}_{\pi^{0}\pi^{0}}-m^{2}_{\pi^{+}}}{2M_{\pi^{0}\pi^{0}}},$
$\displaystyle E^{*}_{\pi^{0}}$ $\displaystyle=$
$\displaystyle\frac{M^{2}_{\pi^{0}\pi^{0}}-m^{2}_{\pi^{0}}+m^{2}_{\pi^{0}}}{2M_{\pi^{0}\pi^{0}}}.$
(34)

Similarly, we can obtain the $\pi^{+}\pi^{0}$ invariant mass distribution.

##  III Numerical results and discussion

Figure 3: The $\pi^{0}\pi^{0}$ (a) and $\pi^{+}\pi^{0}$ (b) mass distributions
for the process $D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$. The experimental datas of
BESIII are represented by points with error bars, and the fitting results of
theoretical calculations by the solid black lines. The colored dashed lines
indicate the components of the theoretical model.

Figure 4: The Dalitz plots of ‘$M^{2}_{\pi^{+}\pi^{0}}$’ vs.
‘$M^{2}_{\pi^{0}\pi^{0}}$’ (a) and ‘$M^{2}_{\pi^{+}\pi^{0}}$’ vs.
‘$M^{2}_{\pi^{+}\pi^{0}}$’ (b) for the process
$D^{+}_{s}\to\pi^{+}\pi^{0}\pi^{0}$.

In our model, there are five free parameters, (1) $V_{p}$, the global
normalization of $f_{0}(980)$ in Eq. (20), (2) the constant $a_{\rm bg}$ as
the background contribution, (3) $\alpha$ corresponding to the strength of the
$f_{0}(1370)$ amplitude in Eq. (24), (4) the width $\Gamma_{f_{0}(1370)}$, (5)
$\beta$ as the strength of the $f_{2}(1270)$ amplitude in Eq. (25). In order
to present our numerical results, we fit the free parameters to the BESIII
measurements of the $\pi^{0}\pi^{0}$ and $\pi^{+}\pi^{0}$ invariant mass
distributions. Here we have $V_{p}=200$, $a_{\rm bg}=300$,
$\Gamma_{f_{0}(1370)}=236$ MeV, $\alpha=85$, and $\beta=250$.444The free
parameters $V_{p}$, $a_{\rm bg}$, $\alpha$, and $\beta$ could be complex,
which will affect the interferences of the different contributions.
Considering our results could reasonably describe the experimental results, we
take them to real in order to reduce the numbers of the free parameters.

Then we have calculated the $\pi^{0}\pi^{0}$ and $\pi^{+}\pi^{0}$ invariant
mass distributions as shown in Fig. 3. The black points with error bars
labeled as ‘BESIII data’ are the BESIII data taken from Ref. [15], and the
solid black curves labeled as ‘Total’ are our theoretical results for the
total contributions of Eq. (29). In addition, we present the contributions
from the $f_{0}(980)$, $f_{0}(1370)$, $f_{2}(1270)$, and the background, which
are labeled as ‘$f_{0}(980)$’, ‘$f_{0}(1370)$’, ‘$f_{2}(1270)$’, and
‘Background’, respectively. One can find a significant cusp signal around 980
MeV in the $\pi^{0}\pi^{0}$ invariant mass distribution, which can be
associated with the scalar $f_{0}(980)$ dynamically generated from the
$S$-wave pseudoscalar-pseudoscalar interaction. The enhancement structure
around 1300 MeV could be well described, and mainly comes from the
contributions of the $f_{0}(1370)$ and $f_{2}(1270)$. Indeed, the mass and
width of the $f_{0}(1370)$ have large uncertainties, and it is difficult to
extract the exact properties of $f_{0}(1370)$ due to the overlap with the
signal of $f_{2}(1270)$. Thus, the more precise of the $\pi^{0}\pi^{0}$ mass
distribution of $D^{+}_{s}\to\pi^{+}\pi^{0}\pi^{0}$ should be more helpful to
shed light on the nature of the resonance $f_{0}(1370)$ and $f_{2}(1270)$.

For the $\pi^{+}\pi^{0}$ invariant mass distribution, our results are in
agreement with the BESIII results considering the large experimental
uncertainties. Indeed, one could find a structures around 1.6 GeV, which could
be due to resonance $\rho_{(}1700)$. Since the BESIII measurements have large
uncertainties, we do not consider the possible contribution from the resonance
$\rho(1700)$ in this work, and one could perform the complicate calculations
when more precise data is available.

We also present the Dalitz plots of ‘$M^{2}_{\pi^{+}\pi^{0}}$’ vs.
‘$M^{2}_{\pi^{0}\pi^{0}}$’ and ‘$M^{2}_{\pi^{+}\pi^{0}}$’ vs.
‘$M^{2}_{\pi^{+}\pi^{0}}$’ for the process $D^{+}_{s}\to\pi^{+}\pi^{0}\pi^{0}$
in Figs. 4(a) and 4(b), and one can easily find the contributions of the
$f_{0}(980)$, $f_{0}(1370)$, and $f_{2}(1270)$. BESIII has reported the Dalitz
plot of ‘$M^{2}_{\pi^{+}\pi^{0}}$’ vs. ‘$M^{2}_{\pi^{+}\pi^{0}}$’, and one can
find our results are consistent with the BESIII measurements.

##  IV Conclusions

Motivated by the recent BESIII measurments about the decay
$D_{s}^{+}\to\pi^{+}\pi^{0}\pi^{0}$ [15], we have investigated this Cabibbo-
favored process by considering the $S$-wave pseudoscalar-pseudoscalar
interactions within the chiral unitary approach, which dynamically generates
the scalar $f_{0}(980)$. In addition, the contributions from the intermediate
resonances $f_{0}(1370)$ and $f_{2}(1270)$ are taken into account.

we have calculated the $\pi^{0}\pi^{0}$ and $\pi^{+}\pi^{0}$ mass
distributions, and find a peak around 980 MeV, which could be associated with
the scalar $f_{0}(980)$, and an enhancement structure around 1300 MeV, mainly
due to the intermediate resonances $f_{0}(1370)$ and $f_{2}(1270)$. Our
results of both the $\pi^{0}\pi^{0}$ and $\pi^{+}\pi^{0}$ invariant mass
distributions are in agreement with the BESIII measurements. In addition, we
have predicted the Dalitz plots of ‘$M^{2}_{\pi^{+}\pi^{0}}$’ vs.
‘$M^{2}_{\pi^{0}\pi^{0}}$’ and ‘$M^{2}_{\pi^{+}\pi^{0}}$’ vs.
‘$M^{2}_{\pi^{+}\pi^{0}}$’.

###### Acknowledgements.

We warmly thank Prof. Ju-Jun Xie and Prof. Li-Sheng Geng for useful
discussions. This work is supported by the National Natural Science Foundation
of China under Grant No. 12192263, the Natural Science Foundation of Henan
under Grand No. 222300420554, the Project of Youth Backbone Teachers of
Colleges and Universities of Henan Province (2020GGJS017), the Youth Talent
Support Project of Henan (2021HYTP002), and the Open Project of Guangxi Key
Laboratory of Nuclear Physics and Nuclear Technology, No.NLK2021-08.

## References

  * [1] R. L. Workman et al. [Particle Data Group], Review of Particle Physics, PTEP 2022, 083C01 (2022)
  * [2] G. ’t Hooft, G. Isidori, L. Maiani, A. D. Polosa and V. Riquer, A Theory of Scalar Mesons, Phys. Lett. B 662 (2008), 424-430
  * [3] F. E. Close and N. A. Tornqvist, Scalar mesons above and below 1-GeV, J. Phys. G 28 (2002), R249-R267
  * [4] G. Y. Wang, S. C. Xue, G. N. Li, E. Wang and D. M. Li, Strong decays of the higher isovector scalar mesons, Phys. Rev. D 97 (2018) no.3, 034030
  * [5] C. R. Deng, J. L. Ping and F. Wang, Dynamical Study of the Light Scalar Mesons below 1 GeV in a Flux-tube Model, Chin. Phys. C 37, 033101 (2013)
  * [6] G. Y. Wang, N. C. Wei, H. M. Yang, E. Wang, L. S. Geng and J. J. Xie, Roles of $a_{0}(980)$, $\Lambda(1670)$, and $\Sigma(1385)$ in the ${\Lambda}_{c}^{+}\to{\eta}{\Lambda}{\pi}^{+}$ decay, Phys. Rev. D 106 (2022) no.5, 056001
  * [7] M. Y. Duan, J. Y. Wang, G. Y. Wang, E. Wang and D. M. Li, Role of scalar $a_{0}(980)$ in the single Cabibbo suppressed process $D^{+}\rightarrow\pi^{+}\pi^{0}\eta$, Eur. Phys. J. C 80 (2020) no.11, 1041
  * [8] Z. Wang, Y. Y. Wang, E. Wang, D. M. Li and J. J. Xie, The scalar $f_{0}(500)$ and $f_{0}(980)$ resonances and vector mesons in the single Cabibbo-suppressed decays $\Lambda_{c}\to pK^{+}K^{-}$ and $p\pi^{+}\pi^{-}$, Eur. Phys. J. C 80 (2020) no.9, 842
  * [9] J. Y. Wang, M. Y. Duan, G. Y. Wang, D. M. Li, L. J. Liu and E. Wang, The $a_{0}(980)$ and $f_{0}(980)$ in the process $D_{s}^{+}\to K^{+}K^{-}{\pi}^{+}$, Phys. Lett. B 821, 136617 (2021)
  * [10] M. Ablikim et al. [BESIII], Study of the decay $D_{s}^{+}\to K_{S}^{0}K_{S}^{0}\pi^{+}$ and observation an isovector partner to $f_{0}(1710)$, Phys. Rev. D 105 (2022) no.5, L051103
  * [11] M. Ablikim et al. [BESIII], Observation of an $a_{0}$-like State with Mass of 1.817 GeV in the Study of $D_{s}^{+}\to K_{S}^{0}K^{+}{\pi}^{0}$ Decays, Phys. Rev. Lett. 129 (2022) no.18, 182001
  * [12] X. Zhu, H. N. Wang, D. M. Li, E. Wang, L. S. Geng and J. J. Xie, Further understanding the nature of $a_{0}(1710)$ in the $D^{+}_{s}\to\pi^{0}K^{+}K^{0}_{S}$ decay, Phys. Rev. D 107 (2023) no.3, 034001
  * [13] X. Zhu, D. M. Li, E. Wang, L. S. Geng and J. J. Xie, Theoretical study of the process $D_{s}^{+}\to{\pi}^{+}K_{S}^{0}K_{S}^{0}$ and the isovector partner of $f_{0}(1710)$, Phys. Rev. D 105 (2022) no.11, 116010
  * [14] L. R. Dai, E. Oset and L. S. Geng, The $D_{s}^{+}\rightarrow\pi^{+}K_{S}^{0}K_{S}^{0}$ reaction and the $I=1$ partner of the $f_{0}(1710)$ state, Eur. Phys. J. C 82 (2022) no.3, 225
  * [15] M. Ablikim et al. [BESIII], Amplitude analysis and branching fraction measurement of the decay $D_{s}^{+}\to{\pi}^{+}{\pi}^{0}{\pi}^{0}$, JHEP 01, 052 (2022)
  * [16] P. Naik et al. [CLEO], Measurement of the Pseudoscalar Decay Constant $f(D_{s})$ Using $D_{s}^{+}\to\tau^{+}\nu,\tau^{+}\to\rho^{+}\bar{\nu}$ Decays, Phys. Rev. D 80, 112004 (2009)
  * [17] J. A. Oller and E. Oset, Chiral symmetry amplitudes in the S wave isoscalar and isovector channels and the $\sigma$, f0(980), a0(980) scalar mesons, Nucl. Phys. A 620 (1997), 438-456 [erratum: Nucl. Phys. A 652 (1999), 407-409]
  * [18] J. A. Oller, E. Oset and J. R. Pelaez, Meson D meson interaction in a nonperturbative chiral approach, Phys. Rev. 59 (1999), 074001 [erratum: Phys. Rev. D 60 (1999), 099906; erratum: Phys. Rev. D 75 (2007), 099903]
  * [19] R. Molina, D. Nicmorus and E. Oset, The $\rho\rho$ interaction in the hidden gauge formalism and the $f_{0}(1370)$ and $f_{2}(1270)$ resonances, Phys. Rev. D 78 (2008), 114018
  * [20] L. S. Geng and E. Oset, Vector meson-vector meson interaction in a hidden gauge unitary approach, Phys. Rev. D 79 (2009), 074009
  * [21] C. García-Recio, L. S. Geng, J. Nieves, L. L. Salcedo, E. Wang and J. J. Xie, Low-lying even parity meson resonances and spin-flavor symmetry revisited, Phys. Rev. D 87 (2013) no.9, 096006
  * [22] L. S. Celenza, S. f. Gao, B. Huang, H. Wang and C. M. Shakin, Covariant confinement model for the calculation of the properties of scalar mesons, Phys. Rev. C 61 (2000), 035201
  * [23] P. Minkowski and W. Ochs, Identification of the glueballs and the scalar meson nonet of lowest mass, Eur. Phys. J. C 9 (1999), 283-312
  * [24] S. Janowski, F. Giacosa and D. H. Rischke, Is $f_{0}(1710)$ a glueball?, Phys. Rev. D 90 (2014) no.11, 114005
  * [25] D. Gülmez, U. G. Meißner and J. A. Oller, A chiral covariant approach to $\rho\rho$ scattering, Eur. Phys. J. C 77 (2017) no.7, 460
  * [26] L. L. Wei, H. S. Li, E. Wang, J. J. Xie, D. M. Li and Y. X. Li, Search for a $D\bar{D}$ bound state in the $\Lambda_{b}\rightarrow\Lambda D\bar{D}$ process, Phys. Rev. D 103, 114013 (2021)
  * [27] L. R. Dai, J. M. Dias and E. Oset, Disclosing $D^{*}\bar{D}^{*}$ molecular states in the $B_{c}^{-}\rightarrow\pi^{-}J/\psi\omega$ decay, Eur. Phys. J. C 78, no.3, 210 (2018)
  * [28] Y. Zhang, E. Wang, D. M. Li and Y. X. Li, Search for the $D^{*}\bar{D}^{*}$ molecular state $Z_{c}(4000)$ in the reaction $B^{-}\rightarrow J/\psi\rho^{0}K^{-}$, Chin. Phys. C 44, no.9, 093107 (2020)
  * [29] L. R. Dai, G. Y. Wang, X. Chen, E. Wang, E. Oset and D. M. Li, The $B^{+}\rightarrow J/\psi\omega K^{+}$ reaction and $D^{\ast}\bar{D}^{\ast}$ molecular states, Eur. Phys. J. A 55, no.3, 36 (2019)
  * [30] W. H. Liang and E. Oset, $B^{0}$ and $B^{0}_{s}$ decays into $J/\psi$ $f_{0}(980)$ and $J/\psi$ $f_{0}(500)$ and the nature of the scalar resonances, Phys. Lett. B 737, 70-74 (2014)
  * [31] J. M. Dias, F. S. Navarra, M. Nielsen and E. Oset, $f_{0}(980)$ production in $D_{s}^{+}\rightarrow\pi^{+}\,\pi^{+}\,\pi^{-}$ and $D_{s}^{+}\rightarrow\pi^{+}\,K^{+}\,K^{-}$ decays, Phys. Rev. D 94, no.9, 096002 (2016)
  * [32] D. Gamermann, E. Oset, D. Strottman and M. J. Vicente Vacas, Dynamically generated open and hidden charm meson systems, Phys. Rev. D 76 (2007), 074016
  * [33] H. A. Ahmed, Z. Y. Wang, Z. F. Sun and C. W. Xiao, Study $B^{0}_{(s)}$ decays into $\phi$ and a scalar or vector meson, Eur. Phys. J. C 81 (2021), 695.
  * [34] E. Wang, H. X. Chen, L. S. Geng, D. M. Li and E. Oset, Hidden-charm pentaquark state in $\Lambda^{0}_{b}\to J/\psi p\pi^{-}$ decay, Phys. Rev. D 93, no.9, 094001 (2016)