11institutetext: University of Helsinki, Finland, 11email:
<EMAIL_ADDRESS>
# Approximation algorithms for confidence bands for time series

Nikolaj Tatti  0000-0002-2087-5360

###### Abstract

Confidence intervals are a standard technique for analyzing data. When applied
to time series, confidence intervals are computed for each time point
separately. Alternatively, we can compute confidence bands, where we are
required to find the smallest area enveloping $k$ time series, where $k$ is a
user parameter. Confidence bands can be then used to detect abnormal time
series, not just individual observations within the time series. We will show
that despite being an NP-hard problem it is possible to find optimal
confidence band for some $k$. We do this by considering a different problem:
discovering regularized bands, where we minimize the envelope area minus the
number of included time series weighted by a parameter $\alpha$. Unlike normal
confidence bands we can solve the problem exactly by using a minimum cut. By
varying $\alpha$ we can obtain solutions for various $k$. If we have a
constraint $k$ for which we cannot find appropriate $\alpha$, we demonstrate a
simple algorithm that yields
$\mathit{\mathcal{O}}\mathopen{}\left(\sqrt{n}\right)$ approximation guarantee
by connecting the problem to a minimum $k$-union problem. This connection also
implies that we cannot approximate the problem better than
$\mathit{\mathcal{O}}\mathopen{}\left(n^{1/4}\right)$ under some (mild)
assumptions. Finally, we consider a variant where instead of minimizing the
area we minimize the maximum width. Here, we demonstrate a simple
2-approximation algorithm and show that we cannot achieve better approximation
guarantee.

##  1 Introduction

Confidence intervals are a common tool to summarize the underlying
distribution, and to indicate outlier behaviour. In this paper we will study
the problem of computing confidence intervals for time series.

Korpela et al. [11] proposed a notion for computing confidence intervals:
instead of computing point-wise confidence intervals, the authors propose
computing confidence bands. More formally, given $n$ time series $T$, we are
asked to find $k$ time series $U\subseteq T$ that minimize the envelope area,
that is, the sum $\sum_{i}\left(\max_{t\in U}t(i)\right)-\left(\min_{t\in
U}t(i)\right)$. The benefit, as argued by Korpela et al. [11], of using
confidence bands instead of point-wise confidence intervals is better family-
wise error control: if we were to use point-wise intervals we can only say
that a time series _at some fixed point_ is an outlier and require a
correction for multiple testing (such as Bonferroni correction) if we want to
state with a certain probability that the _whole_ time series is normal.

In this paper we investigate the approximation algorithms for finding
confidence bands. While Korpela et al. [11] proved that finding the optimal
confidence band is an NP-hard problem, they did not provide any approximation
algorithms nor any inapproximability results.

We will first show that despite being an NP-hard problem, we can solve the
problem for _some_ $k$. We do this by considering a different problem, where
instead of having a hard constraint we have an objective function that prefers
selecting time series as long as they do not increase the envelope area too
much. The objective depends on the parameter $\alpha$, larger values of
$\alpha$ allow more increase in the envelope area. We will show that this
problem can be solved exactly in polynomial time and that each $\alpha$
correspond to a certain value of $k$. We will show that there are at most
$n+1$ of such bands, and that we can discover all of them in polynomial time
by varying $\alpha$.

Next, we provide a simple algorithm for approximating confidence bands by
connecting the problem to the weighted $k$-MinUnion problem. We will provide a
variant of an algorithm by Chlamtáč et al. [2] that yields $\sqrt{n}+1$
guarantee. We also argue that—under certain conjecture—we cannot approximate
the problem better than $\mathit{\mathcal{O}}\mathopen{}\left(n^{1/4}\right)$.

Finally, we consider a variant of the problem where instead of minimizing the
envelope area, we minimize the width of the envelope, that is, we minimize the
maximum difference between the envelope boundaries. We show that a simple
algorithm can achieve 2-approximation. This approximation provides interesting
contrast to the inapproximability results when minimizing the envelope area.
Surprisingly this guarantee is tight: we will also show that the there is no
polynomial-time algorithm with smaller guarantee unless
$\textbf{P}=\textbf{NP}$.

The remainder of the paper is organized as follows. We define the optimization
problems formally in Section 2. We solve the regularized band problem in
Section 3, approximate minimization of envelope area in Section 4, and
approximate minimization of envelope width in Section 5. Section 6 is devoted
to the related work. We present our experiments in Section 7 and conclude with
discussion in Section 8.

##  2 Preliminaries and problem definitions

Assume that we are given time series $T$ with each time series ${f}:{D}\to{}$
mapping from domain $D$ to a real number. We will often write
$n={\left|T\right|}$ to be the number of given time series, and
$m={\left|D\right|}$ to mean the size of the domain.

Given a set of time series $T$, we define the upper and lower _envelopes_ as

$\mathit{ub}\mathopen{}\left(T,i\right)=\max_{t\in
T}t(i)\quad\text{and}\quad\mathit{\ell
b}\mathopen{}\left(T,i\right)=\min_{t\in T}t(i)\quad.$

Our main goal is to find $k$ time series that minimize the envelope area.

######  Problem 1 (SumBand)

Given a set of $n$ time series $T=\left(t_{1},\ldots,t_{n}\right)$, an integer
$k\leq n$, and a time series $x\in T$ find $k$ time series $U\subseteq T$
containing $x$ minimizing

$\mathit{s_{1}}\mathopen{}\left(U\right)=\sum_{i}\mathit{ub}\mathopen{}\left(U,i\right)-\mathit{\ell
b}\mathopen{}\left(U,i\right)\quad.$

We will refer to $U$ as _confidence bands_.

Note that the we also require that we must specify at least one sequence $x\in
T$ that must be included in the input whereas the original definition of the
problem given by Korpela et al. [11] did not require specifying $x$. As we
will see later, this requirement simplifies the computational problem. On the
other hand, if we do not have $x$ at hand, then we can either test every $t\in
T$ as $x$, or we can use the mean or the median of $T$. We will use the latter
option as it does not increase the computational complexity and at the same
time is a reasonable assumption. Note that in this case most likely $x\notin
T$, so we define $T^{\prime}=T\cup\left\\{x\right\\}$, increase
$k^{\prime}=k+1$, and solve SumBand for $T^{\prime}$ and $k^{\prime}$ instead.

We can easily show that the area function
$\mathit{s_{1}}\mathopen{}\left(\cdot\right)$ is a submodular function for all
non-empty subsets, that is,

$\mathit{s_{1}}\mathopen{}\left(U\cup\left\\{t\right\\}\right)-\mathit{s_{1}}\mathopen{}\left(U\right)\leq\mathit{s_{1}}\mathopen{}\left(W\cup\left\\{t\right\\}\right)-\mathit{s_{1}}\mathopen{}\left(W\right),$

where $U\supseteq W\neq\emptyset$. In other words, adding $t$ to a larger set
$U$ increases the cost less than adding $t$ to $W$.

We also consider a variant of SumBand where instead of minimizing the area of
the envelope, we will minimize the maximum width.

######  Problem 2 (InfBand)

Given a set of $n$ time series $T=\left(t_{1},\ldots,t_{n}\right)$, an integer
$k\leq n$, and a time series $x\in T$, find $k$ time series $U\subseteq T$
containing $x$ minimizing

$\mathit{s_{\infty}}\mathopen{}\left(U\right)=\max_{i}\mathit{ub}\mathopen{}\left(U,i\right)-\mathit{\ell
b}\mathopen{}\left(U,i\right)\quad.$

We will show that we can 2-approximate InfBand and that the ratio is tight.

Finally, we consider a regularized version of SumBand, where instead of
requiring that the set has a minimum size $k$, we add a term
$-\alpha{\left|U\right|}$ into the objective function. In other words, we will
favor larger sets as long as the area
$\mathit{s_{1}}\mathopen{}\left(U\right)$ does not increase too much.

######  Problem 3 (RegBand)

Given a set of $n$ time series $T=\left(t_{1},\ldots,t_{n}\right)$, a number
$\alpha>0$, and a time series $x\in T$, find a subset $U\subseteq T$
containing $x$ minimizing

$\mathit{s_{\mathit{reg}}}\mathopen{}\left(U;\alpha\right)=\mathit{s_{1}}\mathopen{}\left(U\right)-\alpha{\left|U\right|}\quad.$

In case of ties, use ${\left|U\right|}$ as a tie-breaker, preferring larger
values.

We refer to the solutions of RegBand as regularized bands. It turns out that
RegBand can be solved in polynomial time. Moreover, the solutions we obtain
from RegBand will be useful for approximating SumBand.

##  3 Regularized bands

In this section we will list useful properties of of RegBand, show how can we
solve RegBand in polynomial time for a single $\alpha$, and finally
demonstrate how we can discover _all_ regularized bands by varying $\alpha$.

###  3.1 Properties of regularized bands

Our first observation is that the output of RegBand also solves SumBand for
certain size constraints.

###### Proposition 1

Assume time series $T$ and $\alpha>0$. Let $U$ be a solution to
$\textsc{RegBand}(\alpha)$. Then $U$ is also a solution for SumBand with
$k={\left|U\right|}$.

The proof of this proposition is trivial and is omitted.

Our next observation is that the solutions to RegBand form a chain.

###### Proposition 2

Assume time series $T$ and $0<\alpha<\beta$. Let $V$ be a solution to
$\textsc{RegBand}(T,\alpha)$ and let $U$ be a solution to
$\textsc{RegBand}(T,\beta)$. Then $V\subseteq U$.

###### Proof

Assume otherwise. Let $W=V\setminus U$. Due to the optimality of $V$,

$0\geq\mathit{s_{\mathit{reg}}}\mathopen{}\left(V;\alpha\right)-\mathit{s_{\mathit{reg}}}\mathopen{}\left(V\cap
U;\alpha\right)=\mathit{s_{1}}\mathopen{}\left(V\right)-\mathit{s_{1}}\mathopen{}\left(V\cap
U\right)-\alpha{\left|W\right|}\quad.$

Since $\mathit{s_{1}}$ is a submodular function, we have

$\mathit{s_{1}}\mathopen{}\left(V\right)-\mathit{s_{1}}\mathopen{}\left(V\cap
U\right)=\mathit{s_{1}}\mathopen{}\left(W\cup(V\cap
U)\right)-\mathit{s_{1}}\mathopen{}\left(V\cap
U\right)\geq\mathit{s_{1}}\mathopen{}\left(W\cup
U\right)-\mathit{s_{1}}\mathopen{}\left(U\right)\quad.$

Combining these inequalities leads to

$\begin{split}0&\geq\mathit{s_{1}}\mathopen{}\left(W\cup
U\right)-\mathit{s_{1}}\mathopen{}\left(U\right)-\alpha{\left|W\right|}\\\
&\geq\mathit{s_{1}}\mathopen{}\left(W\cup
U\right)-\mathit{s_{1}}\mathopen{}\left(U\right)-\beta{\left|W\right|}\\\
&=\mathit{s_{\mathit{reg}}}\mathopen{}\left(W\cup
U;\beta\right)-\mathit{s_{\mathit{reg}}}\mathopen{}\left(U;\beta\right),\end{split}$

which contradicts the optimality of $U$.∎

This property is particularly useful as it allows clean visualization: the
envelopes resulting from different values of $\alpha$ will not intersect.
Moreover, it allows us to stored all regularized bands by simply storing, per
each time series, the index of the largest confidence band containing the time
series.

Interestingly, this result does not hold for SumBand.

###### Example 1

Consider 4 constant time series $t_{1}=0$, $t_{2}=-1$ and $t_{3}=t_{4}=2$. Set
the seed time series $x=t_{1}$. Then the solution for SumBand with $k=2$ is
$\left\\{t_{1},t_{2}\right\\}$ and the solution SumBand with $k=3$ is
$\left\\{t_{1},t_{3},t_{4}\right\\}$.

###  3.2 Computing regularized band for a single $\alpha$

Our next step is to solve RegBand in polynomial time. Note that since
$\mathit{s_{1}}\mathopen{}\left(\cdot\right)$ is submodular, then so is
$\mathit{s_{\mathit{reg}}}\mathopen{}\left(\cdot\right)$. Minimizing
submodular function is solvable in polynomial-time [15]. Solving RegBand using
a generic solver for minimizing submodular functions is slow, so instead we
will solve the problem by reducing it to a minimum cut problem. In such a
problem, we are given a weighted directed graph $G=(V,E,W)$, two nodes, say
$\theta,\eta\in V$, and ask to partition $V$ into $X\cup Y$ such that
$\theta\in X$ and $\eta\in Y$ minimizing the total weight of edges from $X$ to
$Y$.

In order to define $G$ we need several definitions. Assume we are given $n$
time series $T$, a real number $\alpha$ and a seed time series $x\in T$. Let
$m$ be the size of the domain. For $i\in[m]$, we define
$p_{i}=\left\\{t_{j}(i)\mid j\in[n]\right\\}$ to be the set (with no
duplicates) sorted, smallest values first. In other words, $p_{ij}$ is the
$j$th smallest distinct observed value in $T$ at $i$. Let $P$ be the
collection of all $p_{i}$.

We also define $c_{ij}$ to be the number of time series at $i$ smaller than or
equal to $p_{ij}$, that is, $c_{ij}={\left|\left\\{\ell\in[n]\mid
t_{\ell}(i)\leq p_{ij}\right\\}\right|}$. We also write $c_{i0}=0$.

We are now ready to define our graph. We define a weighted directed graph
$G=(V,E,W)$ as follows. The nodes $V$ have three sets $A$, $B$, and $C$. The
set $A$ has ${\left|P\right|}$ nodes, a node $a_{ij}\in A$ corresponding to
each entry $p_{ij}\in P$. The set $B=\left\\{b_{j}\right\\}$ has $n$ nodes,
and the set $C$ has two nodes, $\theta$ and $\eta$. Here, $\theta$ acts as a
source node and $\eta$ acts as a terminal node.

The edges and the weights are as follows: For each $a_{ij}\in A$ such that
$p_{ij}>x(i)$, we add an edge $(a_{i(j-1)},a_{ij})$ with the weight

$w(a_{i(j-1)},a_{ij})=n-c_{i(j-1)}+\frac{m}{\alpha}(p_{i(j-1)}-x(i))\quad.$

For each $a_{ij}\in A$ such that $p_{ij}<x(i)$, we add an edge
$(a_{i(j+1)},a_{ij})$ with the weight

$w(a_{i(j+1)},a_{ij})=c_{ij}+\frac{m}{\alpha}(x(i)-p_{i(j+1)})\quad.$

For each $a_{ij}\in A$ such that $p_{ij}=x(i)$, we add an edge
$(\theta,a_{ij})$ with the weight $\infty$. For each $i\in[m]$ and
$\ell={\left|p_{i}\right|}$, we add two edges $(a_{i\ell},\eta)$ and
$(a_{i1},\eta)$ with the weights

$w(a_{i\ell},\eta)=\frac{m}{\alpha}(p_{i\ell}-x(i))\quad\text{and}\quad
w(a_{i1},\eta)=\frac{m}{\alpha}(x(i)-p_{i1})\quad.$

In addition, for each $i\in[m]$, $\ell\in[n]$, let $j$ be such that
$p_{ij}=t_{\ell}(i)$ and define two edges $(a_{ij},b_{\ell})$ and
$(b_{\ell},a_{ij})$ with the weights,

$\displaystyle w(a_{ij},b_{\ell})$ $\displaystyle=1$ $\displaystyle
w(b_{\ell},a_{ij})$ $\displaystyle=\infty\quad.$

Our next proposition states the minimum cut of $G$ also minimizes RegBand.

###### Proposition 3

Let $X,Y$ be a $(\theta,\eta)$-cut of $G$ with the optimal cost. Define
$f(i)=\min_{j}\left\\{p_{ij}\mid a_{ij}\in X\right\\}$ and
$g(i)=\max_{j}\left\\{p_{ij}\mid a_{ij}\in X\right\\}$.

Then the cost of the cut is equal to

$nm-m{\left|\left\\{j\mid b_{j}\in
X\right\\}\right|}+\frac{m}{\alpha}\sum_{i}g(i)-f(i)\quad.$

Moreover, if $b_{\ell}\in X$, then $g(i)\leq b_{\ell}(i)\leq f(i)$, for all
$i$.

###### Proof

The last claim follows immediately as otherwise there is a cross-edge with
infinite cost making the cut suboptimal.

Define $u(i)=\operatorname*{arg\,min}_{j}\left\\{p_{ij}\mid a_{ij}\in
X\right\\}$ and $v(i)=\operatorname*{arg\,max}_{j}\left\\{p_{ij}\mid a_{ij}\in
X\right\\}$ to be the indices yielding $f$ and $g$. Define also

$d_{i}={\left|\left\\{j\mid u(i)\leq t_{j}(i)\leq
v(i)\right\\}\right|}=c_{iv(i)}-c_{i(u(i)-1)}$

to be the number of time series between $u(i)$ and $v(i)$ at $i$.

Note that $a_{ij}\in X$ whenever $u(i)\leq j\leq v(j)$ as otherwise we can
move $a_{ij}$ to $X$ and decrease the cost.

The cut consists of the cross-edges originating from $a_{iv(i)}$ and
$a_{iu(i)}$, and cross-edges between $A$ and $B$. The cost of the former is
equal to

$\begin{split}&\sum_{i}n-c_{iv(i)}+m\frac{p_{iv(i)}-x(i)}{\alpha}+c_{i(u(i)-1)}+m\frac{x(i)-p_{iu(i)}}{\alpha}\\\
&\qquad=nm+\sum\frac{m}{\alpha}(g(i)-f(i))-\sum_{i}d_{i}\end{split}$

while the cost of the latter is

$\begin{split}\sum_{i}{\left|\left\\{j\mid a_{ij}\in X,b_{j}\notin
X\right\\}\right|}&=\sum_{i}{\left|\left\\{j\mid u(i)\leq t_{j}(i)\leq v(i)\in
X,b_{j}\notin X\right\\}\right|}\\\ &=\sum_{i}d_{i}-{\left|\left\\{j\mid
u(i)\leq t_{j}(i)\leq v(i)\in X,b_{j}\in X\right\\}\right|}\\\
&=\sum_{i}d_{i}-m{\left|\left\\{j\mid b_{j}\in X\right\\}\right|}\quad.\\\
\end{split}$

Combining the two equations proves the claim.∎

###### Corollary 1

Let $U^{\prime}$ be the solution to $\textsc{RegBand}(\alpha)$. Let $(X,Y)$ be
a minimum $(\theta,\eta)$-cut of $G$. Set $U=\left\\{t_{\ell}\mid b_{\ell}\in
X\right\\}$. Then
$\mathit{s_{\mathit{reg}}}\mathopen{}\left(U;\alpha\right)=\mathit{s_{\mathit{reg}}}\mathopen{}\left(U^{\prime};\alpha\right)$.

###### Proof

Proposition 3 states that the cost of the minimum cut is
$nm+\frac{m}{\alpha}\mathit{s_{\mathit{reg}}}\mathopen{}\left(U;\alpha\right)$.

Construct a cut $(X^{\prime},Y^{\prime})$ from $U^{\prime}$ by setting
$X^{\prime}$ to be the nodes from $A$ and $B$ that correspond to the time
series $U^{\prime}$. The proof of Proposition 3 now states that the cut is
equal to
$nm+\frac{m}{\alpha}\mathit{s_{\mathit{reg}}}\mathopen{}\left(U^{\prime};\alpha\right)$.

The optimality of $(X,Y)$ proves the claim.∎

We may encounter a pathological case, where we have multiple cuts with the
same optimal cost. RegBand requires that in such case we use largest solution.
This can be enforced by modifying the weights: first scale the weights so that
they are all multiples of $nm+1$, then add 1 to the weight of each
$(\theta,\alpha_{ij})$. The cut with the modified graph yields the largest
band with the optimal cost.

The constructed graph $G$ has $\mathit{\mathcal{O}}\mathopen{}\left(nm\right)$
nodes and $\mathit{\mathcal{O}}\mathopen{}\left(nm\right)$ edges.
Consequently, we can compute the minimum cut in
$\mathit{\mathcal{O}}\mathopen{}\left((nm)^{2}\right)$ time [13]. In practice,
solving minimum cut is much faster.

###  3.3 Computing all regularized bands

Now that we have a method for solving $\textsc{RegBand}(\alpha)$ for a fixed
$\alpha$, we would like to find solutions for all $\alpha$. Note that
Proposition 2 states that we can have at most $n+1$ different bands.

We can enumerate the bands with the divide-and-conquer approach given in
Algorithm 1. Here, we are given two, already discovered, regularized bands
$U\subsetneq V$, and we try to find a middle band $W$ with $U\subsetneq
W\subsetneq V$. If $W$ exists, we recurse on both sides. To enumerate all
bands, we start with $\textsc{EnumReg}(\left\\{x\right\\},V)$.

1
$\gamma\leftarrow\frac{\mathit{s_{1}}\mathopen{}\left(V\right)-\mathit{s_{1}}\mathopen{}\left(U\right)}{{\left|V\right|}-{\left|U\right|}}-\frac{\Delta}{n^{2}}$;

2 $W\leftarrow$ solution to $\textsc{RegBand}(\gamma)$;

3 if  _$U\neq W$_ then

4 report $W$;

5 $\textsc{EnumReg}(U,W)$; $\textsc{EnumReg}(W,V)$;

6

Algorithm 1 $\textsc{EnumReg}(U,V)$ finds all regularized bands between $U$
and $V$

The following proposition proves the correctness of the algorithm: during each
split we will always find a new band if such exist.

###### Proposition 4

Assume time series $T$ with $n$ time series. Let $\left\\{U_{i}\right\\}$ be
all the possible regularized confidence bands ordered using inclusion. Define

$\Delta=\min\left\\{{\left|t(i)-u(i)\right|}\mid t,u\in T,i,t(i)\neq
u(i)\right\\}\quad.$

Let $i<j$ be two integers and define

$\gamma=\frac{\mathit{s_{1}}\mathopen{}\left(U_{j}\right)-\mathit{s_{1}}\mathopen{}\left(U_{i}\right)}{{\left|U_{j}\right|}-{\left|U_{i}\right|}}-\frac{\Delta}{n^{2}}\quad.$

Let $U_{\ell}$ be the solution for $\textsc{RegBand}(\gamma)$. Then
$i\leq\ell<j$. If $j>i+1$, then $i<\ell$, otherwise $\ell=i$.

For simplicity, let us define
$f(x,y)=\frac{\mathit{s_{1}}\mathopen{}\left(U_{y}\right)-\mathit{s_{1}}\mathopen{}\left(U_{x}\right)}{{\left|U_{y}\right|}-{\left|U_{x}\right|}}$.

In order to prove the result we need the following technical lemma.

###### Lemma 1

Assume time series $T$ with $n$ time series. Let $\left\\{U_{i}\right\\}$ be
all the possible regularized confidence bands ordered using inclusion. Let
$\alpha>0$. Let $U_{i}$ be the solution for $\textsc{RegBand}(\alpha)$. Then
$f(i-1,i)\leq\alpha<f(i,i+1)$.

###### Proof

Due to the optimality of $U_{i}$,

$\mathit{s_{1}}\mathopen{}\left(U_{i}\right)-\alpha{\left|U_{i}\right|}=\mathit{s_{\mathit{reg}}}\mathopen{}\left(U_{i};\alpha\right)<\mathit{s_{1}}\mathopen{}\left(U_{i+1}\right)-\alpha{\left|U_{i+1}\right|}\quad.$

Solving for $\alpha$ gives us the right-hand side of the claim. Similarly,

$\mathit{s_{1}}\mathopen{}\left(U_{i}\right)-\alpha{\left|U_{i}\right|}=\mathit{s_{\mathit{reg}}}\mathopen{}\left(U_{i};\alpha\right)\leq\mathit{s_{1}}\mathopen{}\left(U_{i-1}\right)-\alpha{\left|U_{i-1}\right|}\quad.$

Solving for $\alpha$ gives us the left-hand side of the claim.∎

######  Proof (of Proposition 4)

It is straightforward to see that Lemma 1 implies that $f(a,b)\leq f(x,y)$ for
$a\leq x$ and $b\leq y$. Moreover, the equality holds only if $x=a$ and $y=b$,
in other cases $f(a,b)+\frac{\Delta}{n^{2}}\leq f(x,y)$.

If $\ell\geq j$, then Lemma 1 states that $f(i,j)\leq
f(\ell-1,\ell)\leq\gamma$, which contradicts the definition of $\gamma$. Thus
$\ell<j$.

Since $f(i,j)-f(i-1,i)\geq\frac{\Delta}{n^{2}}$, we have $f(i-1,i)\leq\gamma$.
If $\ell<i$, then Lemma 1 states that $\gamma<f(i-1,i)$, which is a
contradiction. Thus, $\ell\geq i$.

If $j=i+1$, then immediately $\ell=i$.

Assume that $j>i+1$. Since $f(i,j)-f(i,i+1)\geq\frac{\Delta}{n^{2}}$, we have
$f(i,i+1)\leq\gamma$. If $\ell=i$, then according to Lemma 1
$\gamma<f(i,i+1)$, which is a contradiction. Thus, $\ell>i$. ∎

Lemma 1 reveals an illuminating property of regularized bands, namely each
band minimizes the ratio of additional envelope area and the number of new
time series.

###### Proposition 5

Let $U$ be a regularized band. Define
$g(X)=\frac{\mathit{s_{1}}\mathopen{}\left(X\right)-\mathit{s_{1}}\mathopen{}\left(U\right)}{{\left|X\right|}-{\left|U\right|}}$.
Let $V\supsetneq U$ be the adjacent regularized band. Then
$g(V)=\min_{X\supsetneq U}g(X)$.

###### Proof

Let $O=\operatorname*{arg\,min}_{X\supsetneq U}g(X)$, and set $\beta=g(O)$. We
will prove that $g(V)\leq\beta$. Let $W=\textsc{RegBand}(\beta)$. Let $\alpha$
be the parameter for which $U=\textsc{RegBand}(\alpha)$. Assume that
$\alpha\geq\beta$. We can rewrite the equality $\beta=g(O)$ as

$0=\mathit{s_{\mathit{reg}}}\mathopen{}\left(O;\beta\right)-\mathit{s_{\mathit{reg}}}\mathopen{}\left(U;\beta\right)\geq\mathit{s_{\mathit{reg}}}\mathopen{}\left(O;\alpha\right)-\mathit{s_{\mathit{reg}}}\mathopen{}\left(U;\alpha\right),$

which violates the optimality of $U$. Thus $\alpha<\beta$. Proposition 2
states that $U\subseteq W$. Moreover, due to submodularity,

$\mathit{s_{\mathit{reg}}}\mathopen{}\left(O\cup
W;\beta\right)-\mathit{s_{\mathit{reg}}}\mathopen{}\left(W;\beta\right)\leq\mathit{s_{\mathit{reg}}}\mathopen{}\left(O\cup
U;\beta\right)-\mathit{s_{\mathit{reg}}}\mathopen{}\left(U;\beta\right)=0,$

which due to the optimality of $W$ implies that $O\subseteq W$. Thus $W\neq U$
and $V\subseteq W$. Lemma 1, possibly applied multiple times, shows that
$g(V)\leq g(W)\leq\beta$. ∎

Proposition 2 states that there are at most $n+1$ bands. Queries done by
EnumReg yield the same band at most twice. Thus, EnumReg performs at most
$\mathit{\mathcal{O}}\mathopen{}\left(n\right)$ queries, yielding
computational complexity of
$\mathit{\mathcal{O}}\mathopen{}\left(n^{3}m^{2}\right)$. In practice, EnumReg
is faster: the number of bands is significantly smaller than $n$ and the
minimum cut solver scales significantly better than
$\mathit{\mathcal{O}}\mathopen{}\left(n^{2}m^{2}\right)$. Moreover, we can
further improve the performance with the following observation: Proposition 2
states that when processing $\textsc{EnumReg}(U,V)$, the bands will be between
$U$ and $V$. Hence, we can ignore the time series that are outside $V$, and we
can safely replace $U$ with its envelope $\mathit{\ell
b}\mathopen{}\left(U\right)$ and $\mathit{ub}\mathopen{}\left(U\right)$.111We
need to make sure that the envelope is always selected. This can be done by
connecting $\theta$ to the envelope with edges of infinite weight.

##  4 Discovering confidence bands minimizing $\mathit{s_{1}}$

In this section, we will study SumBand. Korpela et al. [11] showed that the
problem is NP-hard. We will argue that we can approximate the problem and
establish a (likely) lower bound for the approximation guarantee.

1 $\left\\{B_{i}\right\\}\leftarrow\textsc{EnumReg}(\left\\{x\right\\},T)$;

2

3$j\leftarrow$ largest index for which ${\left|B_{j}\right|}\leq k$;

4 if ${\left|B_{j}\right|}\leq k-\sqrt{n}$ then $W\leftarrow B_{j+1}\setminus
B_{j}$ else $W\leftarrow T\setminus B_{j}$;

5

6$U\leftarrow B_{j}$;

7 greedily add $k-{\left|U\right|}$ entries from $W$ to $U$, minimizing
$\mathit{s_{1}}$ at each step;

8 return $U$;

Algorithm 2 $\textsc{FindSum}(T,k,x)$, approximates SumBand

As a starting point, note that SumBand is an instance of $k$-MinUnion,
weighted minimum $k$-union problem. In $k$-MinUnion we are given $n$ sets over
a universe with weighted points, and ask to select $k$ sets minimizing the
weighted union. In our case, the universe is the set $P$ described in Section
3, the weights are the distances between adjacent points, and a set consists
of all the points between a time series and $x$.

The _unweighted_ $k$-MinUnion problem has several approximation algorithms: a
simple algorithm achieving
$\mathit{\mathcal{O}}\mathopen{}\left(\sqrt{n}\right)$ guarantee by Chlamtáč
et al. [2] and an algorithm achieving lower approximation guarantee of
$\mathit{\mathcal{O}}\mathopen{}\left(n^{1/4}\right)$ by Chlamtáč et al. [3].
We will use the former algorithm due to its simplicity and the fact that it
can be easily adopted to handle weights.

The pseudo-code for the algorithm is given in Algorithm 2. The algorithm first
looks for the largest possible regularized band, say $B_{j}$, whose size at
most $k$. The remaining time series are then selected greedily from a set of
candidates $W$. The set $W$ depends on how many additional time series is
needed: if we need at most $\sqrt{n}$ additional time series, we set $W$ to be
the remaining time series $T\setminus B_{j}$, otherwise we select the time
series from the next regularized band, that is, we set $W=B_{j+1}\setminus
B_{j}$.

###### Proposition 6

FindSum yields $\sqrt{n}+1$ approximation guarantee.

###### Proof

Let $O$ be the optimal solution for $\textsc{SumBand}(k)$, and let
$r=\mathit{s_{1}}\mathopen{}\left(O\right)$. Let $U$ be the output of FindSum.
Assume that $B_{j}\neq O$, as otherwise we are done. We split the proof in two
cases.

First, assume that ${\left|B_{j}\right|}\leq k-\sqrt{n}$. Since
$\mathit{s_{1}}$ is submodular we have $\mathit{s_{1}}\mathopen{}\left(O\cup
B_{j}\right)-\mathit{s_{1}}\mathopen{}\left(B_{j}\right)\leq\mathit{s_{1}}\mathopen{}\left(O\right)-\mathit{s_{1}}\mathopen{}\left(\left\\{x\right\\}\right)=r$,
leading to

$\frac{\mathit{s_{1}}\mathopen{}\left(B_{j+1}\right)-\mathit{s_{1}}\mathopen{}\left(B_{j}\right)}{{\left|B_{j+1}\right|}-{\left|B_{j}\right|}}\leq\frac{\mathit{s_{1}}\mathopen{}\left(O\cup
B_{j}\right)-\mathit{s_{1}}\mathopen{}\left(B_{j}\right)}{{\left|O\cup
B_{j}\right|}-{\left|B_{j}\right|}}\leq\frac{r}{k-{\left|B_{j}\right|}}\leq\frac{r}{\sqrt{n}},$

where the first inequality is due to Proposition 5. Rearranging the terms
gives us

$\mathit{s_{1}}\mathopen{}\left(B_{j+1}\right)-\mathit{s_{1}}\mathopen{}\left(B_{j}\right)\leq\frac{r({\left|B_{j+1}\right|}-{\left|B_{j}\right|})}{\sqrt{n}}\leq
r\frac{n}{\sqrt{n}}=r\sqrt{n}\quad.$

Finally,

$\begin{split}\mathit{s_{1}}\mathopen{}\left(U\right)&=\mathit{s_{1}}\mathopen{}\left(B_{j}\right)+(\mathit{s_{1}}\mathopen{}\left(U\right)-\mathit{s_{1}}\mathopen{}\left(B_{j}\right))\\\
&\leq\mathit{s_{1}}\mathopen{}\left(B_{j}\right)+(\mathit{s_{1}}\mathopen{}\left(B_{j+1}\right)-\mathit{s_{1}}\mathopen{}\left(B_{j}\right))\leq\mathit{s_{1}}\mathopen{}\left(B_{j}\right)+r\sqrt{n}\leq
r(1+\sqrt{n}),\\\ \end{split}$

where the last inequality is implied by Proposition 1 and the fact that
${\left|B_{j}\right|}\leq k$.

Assume that ${\left|B_{j}\right|}>k-\sqrt{n}$, and let
$q=k-{\left|B_{j}\right|}$. Note that $q<\sqrt{n}$. Let $c_{1},\ldots,c_{q}$
be the additional time series added to $U$. Write
$U_{i}=B_{j}\cup\left\\{c_{1},\ldots,c_{i}\right\\}$.

Let $c_{i}^{\prime}$ be the closest time series to $x$ outside $U_{i-1}$. Note
that
$\mathit{s_{1}}\mathopen{}\left(\left\\{x,c_{i}^{\prime}\right\\}\right)-\mathit{s_{1}}\mathopen{}\left(\left\\{x\right\\}\right)=\mathit{s_{1}}\mathopen{}\left(\left\\{x,c_{i}^{\prime}\right\\}\right)\leq
r$ for $i=1,\ldots,q$ as otherwise $r$ has to be larger. In addition,
Proposition 1 and ${\left|B_{j}\right|}\leq k$ imply that
$\mathit{s_{1}}\mathopen{}\left(B_{j}\right)\leq r$. Consequently,

$\begin{split}\mathit{s_{1}}\mathopen{}\left(U_{q}\right)&=\mathit{s_{1}}\mathopen{}\left(B_{j}\right)+\sum_{i=1}^{q}\mathit{s_{1}}\mathopen{}\left(U_{i}\right)-\mathit{s_{1}}\mathopen{}\left(U_{i-1}\right)\\\
&\leq\mathit{s_{1}}\mathopen{}\left(B_{j}\right)+\sum_{i=1}^{q}\mathit{s_{1}}\mathopen{}\left(\left\\{x,c_{i}^{\prime}\right\\}\right)-\mathit{s_{1}}\mathopen{}\left(\left\\{x\right\\}\right)\leq(1+\sqrt{n})r,\end{split}$

where the first inequality is due to the submodularity of $\mathit{s_{1}}$. ∎

FindSum resembles greatly the algorithm given by Chlamtáč et al. [2] but has
few technical differences: we select $B_{j}$ as our starting point whereas the
algorithm by Chlamtáč et al. [2] constructs the starting set by iteratively
finding and adding sets with the smallest average area,
$\mathit{s_{1}}\mathopen{}\left(X\right)/{\left|X\right|}$, that is, solving
the problem given in Proposition 5.222The original algorithm is described
using set/graph terminology but we use our terminology to describe the
differences. Such sets can be found with a linear program. Proposition 5
implies that both approaches result in the same set $B_{j}$ but our approach
is faster.333The computational complexity of the state-of-the-art linear
program solver is
$\mathit{\mathcal{O}}\mathopen{}\left((nm)^{2.37}\log(nm/\delta)\right)$,
where $\delta$ is the relative accuracy [4]. We may need to solve
$\mathit{\mathcal{O}}\mathopen{}\left(n\right)$ such problems, leading to a
total time of
$\mathit{\mathcal{O}}\mathopen{}\left(n(nm)^{2.37}\log(nm/\delta)\right)$.
Moreover, this modification allows us to prove a tighter approximation
guarantee: the authors prove that their algorithm yields $2\sqrt{n}$ guarantee
whereas we show that we can achieve $\sqrt{n}+1$ guarantee. Additionally, we
select additional time series iteratively by selecting those time series that
result in the smallest increase of the current area, whereas the original
algorithm would simply select time series that are closest to
$\left\\{x\right\\}$.

Chlamtáč et al. [3] argued that under some mild but technical conjecture there
is no polynomial-time algorithm that can approximate $k$-MinUnion better than
$\mathit{\mathcal{O}}\mathopen{}\left(n^{1/4}\right)$. Next we will show that
we can reduce $k$-MinUnion to SumBand while preserving approximation.

###### Proposition 7

If there is an $f(n)$-approximation polynomial-time algorithm for SumBand,
then there is an $f(n+1)$-approximation polynomial-time algorithm for
$k$-MinUnion.

###### Proof

Assume that we are given an instance of $k$-MinUnion with $n$ sets
${\mathcal{S}}=\left(S_{1},\ldots,S_{n}\right)$. Let $D=\bigcup_{i}S_{i}$ be
the union of all $S_{i}$.

Define $T$ containing $n+1$ time series over the domain $D$. The first $n$
time series correspond to the sets $S_{i}$, that is, given $i\in D$, we set
$t_{j}(i)=1$ if $i\in S_{j}$, and $0$ otherwise. The remaining single time
series, named $x$, is set to be $0$.

Assume that we have an algorithm estimating $\textsc{SumBand}(T,x,k+1)$, and
let $U$ be the output of this algorithm. Note that since $x\in U$, we have
$\mathit{\ell b}\mathopen{}\left(U,i\right)=0$.

Let ${\mathcal{V}}$ be the subset of ${\mathcal{S}}$ corresponding to the non-
zero time series in $U$. Let $C=\bigcup_{S\in{\mathcal{V}}}S$ be the union of
sets in ${\mathcal{V}}$. Since $\mathit{\ell b}\mathopen{}\left(U,i\right)=0$,
and $\mathit{ub}\mathopen{}\left(U,i\right)=1$ if and only if $i\in C$, we
have $\mathit{s_{1}}\mathopen{}\left(U\right)={\left|C\right|}$. ∎

The above result implies that unless the conjecture suggested by Chlamtáč et
al. [3] is false, we cannot approximate SumBand better than
$\mathit{\mathcal{O}}\mathopen{}\left(n^{1/4}\right)$. This proposition holds
even if we replace $\mathit{s_{1}}\mathopen{}\left(\cdot\right)$ with an
$\ell_{p}^{p}$ norm, $\sum_{i}{\left|t(i)-u(i)\right|}^{p}$, where $1\leq
p<\infty$, or any norm that reduces to hamming distance if $t$ is a binary
sequence and $u$ is 0. Interestingly, we will show in the next section that we
can achieve a tighter approximation if we use $\mathit{s_{\infty}}$.

##  5 Discovering confidence bands minimizing $\mathit{s_{\infty}}$

In this section we consider the problem InfBand. Namely, we will show that a
straightforward algorithm 2-approximates the problem, and more interestingly
we show that the guarantee is tight.

The algorithm for $\textsc{InfBand}(T,x,k)$ is simple: we select $k$ time
series that are closest to $x$ according to the norm
$\left\|t(i)-x(i)\right\|_{\infty}=\max_{i}{\left|t(i)-x(i)\right|}$. We will
refer to this algorithm as FindInf.

It turns out that this simple algorithm yields 2-approximation guarantee.

###### Proposition 8

FindInf yields 2-approximation for InfBand.

###### Proof

Let $U$ be the optimal solution for InfBand. Let $V$ be the result produced by
FindInf. Define $c=\max_{t\in V}\left\|t-x\right\|_{\infty}$. Then

$c=\max_{t\in V}\left\|t-x\right\|_{\infty}\leq\max_{t\in
U}\left\|u-x\right\|_{\infty}\leq\mathit{s_{\infty}}\mathopen{}\left(U\right),$

where the first inequality holds since $V$ contains the closest time series
and the second inequality holds since $x\in U$.

Let $i$ be the index such that
$\mathit{s_{\infty}}\mathopen{}\left(V\right)=\mathit{ub}\mathopen{}\left(V,i\right)-\mathit{\ell
b}\mathopen{}\left(V,i\right)$. Then

$\mathit{s_{\infty}}\mathopen{}\left(V\right)=\mathit{ub}\mathopen{}\left(V,i\right)-\mathit{\ell
b}\mathopen{}\left(V,i\right)=(\mathit{ub}\mathopen{}\left(V,i\right)-t(i))+(t(i)-\mathit{\ell
b}\mathopen{}\left(V,i\right))\leq 2c\quad.$

Thus, $\mathit{s_{\infty}}\mathopen{}\left(V\right)\leq 2c\leq
2\mathit{s_{\infty}}\mathopen{}\left(U\right)$, proving the claim. ∎

While FindInf is trivial, surprisingly it achieves the best possible
approximation guarantee for a polynomial-time algorithm.

###### Proposition 9

There is no polynomial-time algorithm for InfBand that yields $\alpha<2$
approximation guarantee unless $\textbf{P}=\textbf{NP}$.

###### Proof

To prove the claim we will show that we can solve $k$-Clique in polynomial
time if we can $\alpha$-approximate InfBand with $\alpha<2$. Since $k$-Clique
is an NP-complete problem, this is a contradiction unless
$\textbf{P}=\textbf{NP}$.

The goal of $k$-Clique is given a graph $G=(V,E)$ with $n$ nodes and $m$ edges
to detect whether there is a $k$-clique, a fully connected subgraph with $k$
nodes, in $G$. We can safely assume that $G$ has no nodes that are fully-
connected.

Fix an order for nodes $V=\left(v_{1},\ldots,v_{n}\right)$ and let $F$ be all
the edges that are not in $E$, that is,
$F=\left\\{(v_{x},v_{y})\mid(v_{x},v_{y})\notin E,x<y\right\\}$.

Next, we will define an instance of InfBand. The set of time series
$T=\left(t_{1},\ldots,t_{n}\right)\cup\left\\{x\right\\}$ consists of $n$ time
series $t_{i}$ corresponding to the node $v_{i}$, and a single time series $x$
which we will use a seed. We set the domain to be $F$. Each time series
$t_{i}$ maps an element of $e=(v_{x},v_{y})\in F$ to an integer,

$t_{i}(e)=1,\ \text{if}\ i=x,\quad t_{i}(e)=-1,\ \text{if}\ i=y,\quad
t_{i}(e)=0,\ \text{otherwise}\quad.$

We also set $x=0$. First note that since $t_{i}(e)$ is an integer between $-1$
and $1$, the score $\mathit{s_{\infty}}\mathopen{}\left(U\right)$ is either
$0$, $1$, or $2$ for any $U\subseteq T$.

Since we do not have any fully-connected nodes in $G$, there is no non-zero
$t_{i}$ in $T$. Since $x\in U$ for any solution of InfBand, then
$\mathit{s_{\infty}}\mathopen{}\left(U\right)=0$ implies ${\left|U\right|}=1$.

Let $W\subseteq V$ be a subset of nodes, and let $U$ be the corresponding time
series. We claim that $\mathit{s_{\infty}}\mathopen{}\left(U\right)=1$ if and
only if $W$ is a clique. To prove the claim, first observe that if
$v_{i},v_{j}\in W$ such that $e=(v_{i},v_{j})\in F$, then $t_{i}(e)=1$ and
$t_{j}(e)=-1$, thus $\mathit{s_{\infty}}\mathopen{}\left(U\right)=2$. On the
other hand, if $W$ is a clique, then for every $t_{i},t_{j}\in U$ and $e\in F$
such that $t_{i}(e)\neq 0$, we have $t_{j}(e)=0$ since otherwise
$(v_{i},v_{j})\notin E$. Thus,
$\mathit{s_{\infty}}\mathopen{}\left(U\right)=1$ if and only if $W$ is a
clique.

Let $O$ be the solution for $\textsc{InfBand}(T,k+1,x)$. Note that
$\mathit{s_{\infty}}\mathopen{}\left(O\right)=1$ if and only if $G$ has a
$k$-clique, and $\mathit{s_{\infty}}\mathopen{}\left(O\right)=2$ otherwise.

Let $S$ be the output of $\alpha$-approximation algorithm. Since $k>1$, we
know that $\mathit{s_{\infty}}\mathopen{}\left(O\right)$ is either 1 or 2. If
$\mathit{s_{\infty}}\mathopen{}\left(O\right)=2$, then
$\mathit{s_{\infty}}\mathopen{}\left(S\right)=2$. If
$\mathit{s_{\infty}}\mathopen{}\left(O\right)=1$, then
$\mathit{s_{\infty}}\mathopen{}\left(S\right)\leq\alpha\mathit{s_{\infty}}\mathopen{}\left(O\right)<2\times
1$. Thus, $\mathit{s_{\infty}}\mathopen{}\left(S\right)=1$. In summary,
$\mathit{s_{\infty}}\mathopen{}\left(O\right)=\mathit{s_{\infty}}\mathopen{}\left(S\right)$.

We have shown that $\mathit{s_{\infty}}\mathopen{}\left(S\right)=1$ if and
only if $G$ has a $k$-clique. This allows us to detect $k$-clique in $G$ in
polynomial time proving our claim. ∎

##  6 Related work

Confidence bands are envelopes for which confidence intervals of individual
points hold simultaneously. Davison and Hinkley [5], Mandel and Betensky [12]
proposed a non-parametric approach for finding simultaneous confidence
intervals. Here, time series are ordered based on its _maximum_ value, and
$\alpha$-confidence intervals are obtained by removing $\alpha/2$ portions
from each tail. Note that unlike SumBand and InfBand this definition is not
symmetric: if we flip the sign of time series we may get a different interval.

There is a strong parallel between finding regularized bands and finding dense
subgraphs. Proposition 5 states that the inner-most regularized band has the
smallest average envelope area, or alternatively it has the highest ratio of
time series per envelope area. A related graph-theoretical concept is a dense
subgraph, a subgraph $H$ of a given subgraph $G$ with the largest ratio
${\left|E(H)\right|}/{\left|V(H)\right|}$. The method proposed by Goldberg [7]
for finding dense subgraphs in polynomial time is based on maximizing
${\left|E(H)\right|}-\alpha{\left|V(H)\right|}$ and selecting $\alpha$ to be
as small as possible without having an empty solution. Moreover, Tatti [16]
extended the notion of dense subgraphs to density-friendly core decomposition,
which essentially consists of the subgraphs minimizing
${\left|E(H)\right|}-\alpha{\left|V(H)\right|}$ for various values of
$\alpha$, the algorithm for finding the decomposition is similar to the
algorithm for enumerating all regularized bands. In addition, Tsourakakis [17]
extended the notion of dense subgraphs to triangle-density and hypergraphs,
and also used minimum cut to find the solutions. As pointed out in Section 4
is that we can view time series as sets of points in $P$. In fact, the minimum
cut used in Section 3 share some similarities with the minimum cut proposed by
Tsourakakis [17]. Finally, the algorithm proposed by Korpela et al. [11] to
find confidence bands resembles the algorithm by Charikar [1] for
approximating the densest subgraph: in the former we delete the time series
that reduce the envelope area the most while in the latter we delete vertices
that have the smallest degree.

We assume that we are given a seed time series $x$. If such series is not
given then we need to test every $t\in T$ as a seed. If we consider a special
case of $k=2$, then the problem of finding regularized band reduces to the
closest pair problem: find two time series with the smallest distance: a well-
studied problem in computational geometry. A classic approach by Rabin [14],
Dietzfelbinger et al. [6], Khuller and Matias [10] allows to solve the closest
pair problem in $\mathit{\mathcal{O}}\mathopen{}\left(n\right)$ time but the
analysis treats the size of the domain, $m$, as a constant; otherwise, the
computational complexity has an exponential factor in $m$ and can be only used
for very small values of $m$. For large values of $m$, Indyk et al. [9]
proposed an algorithm for solving the closest pair problem minimizing
$\mathit{s_{1}}\mathopen{}\left(\cdot\right)$ in
$\mathit{\mathcal{O}}\mathopen{}\left(n^{2.687}\right)$ time and minimizing
$\mathit{s_{\infty}}\mathopen{}\left(\cdot\right)$ in
$\mathit{\mathcal{O}}\mathopen{}\left(n^{2.687}\log\Delta\right)$ time, where
$\Delta$ is the width of the envelope of the whole data.

##  7 Experimental evaluation

In this section we describe our experimental evaluation.

We implemented EnumReg and FindSum using C++ and used a laptop with Intel Core
i5 (2.3GHz) to conduct our experiments.444The code is available at
https://version.helsinki.fi/DACS As a baseline we used the algorithm by
Korpela et al. [11], which we will call Peel. We implemented Peel also with
C++, and modified it to make sure that the seed time series $x$ is always
included. Finally, we implemented FindInf with Python. In all algorithms we
used the median as the seed time series.

Table 1: Basic characteristics of the datasets and performance measures of the algorithms. Here, $n$ stands for the number of time series, $m$ stands for the domain size, ${\left|B_{1}\right|}$ is the size of the smallest non-trivial regularized band, ${\left|\mathcal{B}\right|}$ is the number of regularized bands, and time is the required time to execute EnumReg in seconds. The scores $\mathit{s_{1}}$ for the algorithms FindSum, FindInf, and Peel are normalized with the envelope area of the whole data and multiplied by 100.  |  |  |  |  |  |  $\mathit{s_{1}}$ for $k=\left\lfloor 0.9n\right\rfloor$ |  $\mathit{s_{1}}$ for $k=\left\lfloor 0.95n\right\rfloor$
---|---|---|---|---|---|---|---
_Dataset_ | $n$ | $m$ | ${\left|B_{1}\right|}$ | ${\left|\mathcal{B}\right|}$ | Time | Sum | Peel | Inf | Sum | Peel | Inf
Milan | $245$ | $12$ | $209$ | $17$ | $0.03$ | $70.34$ | $72.49$ | $74.1$ | $75.31$ | $76.99$ | $78.45$
Power | $1\,417$ | $24$ | $1\,102$ | $56$ | $3.68$ | $70.94$ | $72.83$ | $77.06$ | $78.89$ | $81.17$ | $82.31$
ECG-normal | $1\,507$ | $253$ | $1\,289$ | $72$ | $39.44$ | $51.72$ | $52$ | $72.97$ | $57.22$ | $57.51$ | $73.15$
ECG-pvc | $520$ | $253$ | $484$ | $19$ | $6.67$ | $80.28$ | $80.02$ | $91.97$ | $83.84$ | $83.92$ | $95.69$
Table 2: Scores $\mathit{s_{\infty}}$ of discovered confidence bands. The scores are normalized with the envelope width of the whole data and multiplied by 100.  |  $\mathit{s_{\infty}}$ for $k=\left\lfloor 0.9n\right\rfloor$ |  $\mathit{s_{\infty}}$ for $k=\left\lfloor 0.95n\right\rfloor$
---|---|---
_Dataset_ | Sum | Peel | Inf | Sum | Peel | Inf
Milan | $72.54$ | $79.78$ | $67.35$ | $78.04$ | $79.78$ | $73.95$
Power | $79.08$ | $82.16$ | $73.13$ | $82.16$ | $98.71$ | $79.08$
ECG-normal | $64.78$ | $64.78$ | $54.81$ | $65.64$ | $64.78$ | $57.39$
ECG-pvc | $93.24$ | $93.24$ | $66.41$ | $93.24$ | $93.24$ | $81.9$

$0$$50$$100$$150$$200$$250$$-1$$0$$1$$2$med.SumInf

$0$$50$$100$$150$$200$$250$$-2$$0$$2$

Figure 1: Envelopes for ECG-normal (left) and ECG-pvc (right) and
$k=\left\lfloor 0.9n\right\rfloor$.

Datasets: We used 4 real-world datasets as benchmark datasets. The first
dataset, Milan, consists of monthly averages of maximum daily temperatures in
Milan between the years 1763–2007.555https://www.ncdc.noaa.gov/ The second
dataset, Power, consists of hourly power consumption (variable
global_active_power) of a single household over almost 4 years, a single time
series representing a
day.666http://archive.ics.uci.edu/ml/datasets/Individual+household+electric+power+consumption
Our last 2 datasets ECG-normal and ECG-pvc are heart beat data [8]. We used
MLII data of a single patient (id 106) from the MIT-BIH arrhythmia
database,777https://physionet.org/content/mitdb/1.0.0/ and split the
measurements into normal beats (ECG-normal) and abnormal beats with premature
ventricular contraction (ECG-pvc). Each time series represent measurements
between $-300$ms and $400$ms around each beat. The sizes of the datasets are
given in Table 1.

Results: First let us consider EnumReg. From the results given Table 1 we see
that the number of distinct regularized bands ${\left|\mathcal{B}\right|}$ is
low: about 4%–7% of $n$, the number of time series. Having so few bands in
practice reduces the computational cost of EnumReg since the algorithm tests
at most $2{\left|\mathcal{B}\right|}$ values of $\alpha$ Interestingly, the
smallest non-trivial band $B_{1}$ is typically large, containing about 70%–90%
of the time series. Note that Proposition 5 states that $B_{1}$ has the
smallest ratio of
$\mathit{s_{1}}\mathopen{}\left(B_{1}\right)/{\left|B_{1}\right|}$. For our
benchmark datasets, $B_{1}$ is large suggesting that most time series are
equally far away from the median while the remaining the time series exhibit
outlier behaviour.

The algorithms are fast for our datasets: Table 1 show that EnumReg requires
at most 40 seconds. Additional steps required by FindSum are negligible,
completing in less than a second. The baseline algorithm is also fast,
requiring less than a second to complete.

Let us now compare FindSum against Peel. We compared the obtained areas by
both algorithms with $k=\left\lfloor 0.9n\right\rfloor$ and $k=\left\lfloor
0.95n\right\rfloor$. We see from the results in Table 1, that FindSum performs
slightly better than Peel. The improvement in score is modest, 1%–2%. We
conjecture that in practice Peel is close to the optimal, so any improvements
are subtle. Interestingly, enough Peel performs better than FindSum for ECG-
pvc and $\gamma=0.1$. The reason for this is that the inner band $B_{1}$
contains more than 90% of the time series. In such a case FindSum will reduce
to a simple greedy method, starting from $\left\\{x\right\\}$. Additional
testing revealed that Peel outperforms FindSum when
$k\leq{\left|B_{1}\right|}$ about 50%–90%, depending on the dataset,
suggesting that whenever $k\leq{\left|B_{1}\right|}$ it is probably better to
run both algorithms and select the better envelope.

Next let us compare FindInf against the other methods. The results in Tables
1–2 show that FindInf yields inferior $\mathit{s_{1}}$ scores but superior
$\mathit{s_{\infty}}$ scores. This is expected as FindInf optimizes
$\mathit{s_{\infty}}$ while FindSum and Peel optimize $\mathit{s_{1}}$. The
differences are further highlighted in the envelopes for ECG datasets shown in
Figure 1: FindInf yields larger envelopes but provides a tighter bound under
the peak (R wave).

##  8 Concluding remarks

In this paper we consider the approximation algorithms for discovering
confidence bands. Namely, we proposed a practical algorithm that approximates
SumBand with a guarantee of
$\mathit{\mathcal{O}}\mathopen{}\left(n^{1/2}\right)$. We also argued that the
lower bound for the guarantee is most likely
$\mathit{\mathcal{O}}\mathopen{}\left(n^{1/4}\right)$. In addition, we showed
that we can 2-approximate InfBand, a variant of SumBand problem, with a simple
algorithm and that the guarantee is tight.

Our experiments showed that FindSum outperforms the original baseline method
for large values of $k$, that is, as long as $k$ is larger than the smallest
regularized band. Our experiments suggest that this condition usually holds,
if we are interested, say in, 90%–95% confidence.

Interesting future line of work is to study the case for time series with
multiple modes, that is, a case where instead of a single seed time series, we
are given a set of time series, and we are asked to find confidence bands
around each seed.

## References

  * Charikar [2000] Charikar, M.: Greedy approximation algorithms for finding dense components in a graph.  APPROX (2000)
  * Chlamtáč et al. [2018] Chlamtáč, E., Dinitz, M., Konrad, C., Kortsarz, G., Rabanca, G.: The densest k-subhypergraph problem.  SIAM Journal on Discrete Mathematics 32(2), 1458–1477 (2018)
  * Chlamtáč et al. [2017] Chlamtáč, E., Dinitz, M., Makarychev, Y.: Minimizing the union: Tight approximations for small set bipartite vertex expansion.  In: Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 881–899. SIAM (2017)
  * Cohen et al. [2021] Cohen, M.B., Lee, Y.T., Song, Z.: Solving linear programs in the current matrix multiplication time.  Journal of the ACM (JACM) 68(1), 1–39 (2021)
  * Davison and Hinkley [1997] Davison, A.C., Hinkley, D.V.: Bootstrap methods and their application.  Cambridge university press (1997)
  * Dietzfelbinger et al. [1997] Dietzfelbinger, M., Hagerup, T., Katajainen, J., Penttonen, M.: A reliable randomized algorithm for the closest-pair problem.  Journal of Algorithms 25(1), 19–51 (1997)
  * Goldberg [1984] Goldberg, A.V.: Finding a maximum density subgraph.  University of California Berkeley Technical report (1984)
  * Goldberger et al. [2000] Goldberger, A.L., Amaral, L.A., Glass, L., Hausdorff, J.M., Ivanov, P.C., Mark, R.G., Mietus, J.E., Moody, G.B., Peng, C.K., Stanley, H.E.: Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals.  circulation 101(23), e215–e220 (2000)
  * Indyk et al. [2004] Indyk, P., Lewenstein, M., Lipsky, O., Porat, E.: Closest pair problems in very high dimensions.  In: International Colloquium on Automata, Languages, and Programming. pp. 782–792. Springer (2004)
  * Khuller and Matias [1995] Khuller, S., Matias, Y.: A simple randomized sieve algorithm for the closest-pair problem.  Information and Computation 118(1), 34–37 (1995)
  * Korpela et al. [2014] Korpela, J., Puolamäki, K., Gionis, A.: Confidence bands for time series data.  Data mining and knowledge discovery 28(5), 1530–1553 (2014)
  * Mandel and Betensky [2008] Mandel, M., Betensky, R.A.: Simultaneous confidence intervals based on the percentile bootstrap approach.  Computational statistics & data analysis 52(4), 2158–2165 (2008)
  * Orlin [2013] Orlin, J.B.: Max flows in $O(nm)$ time, or better.  In: Proceedings of the forty-fifth annual ACM symposium on Theory of computing. pp. 765–774 (2013)
  * Rabin [1976] Rabin, M.O.: Probabilistic algorithms.  In: Traub, J.F. (ed.) Algorithms and Complexity: New Directions and Recent Results. Academic Press New York (1976)
  * Schrijver [2000] Schrijver, A.: A combinatorial algorithm minimizing submodular functions in strongly polynomial time.  Journal of Combinatorial Theory, Series B 80(2), 346–355 (2000)
  * Tatti [2019] Tatti, N.: Density-friendly graph decomposition.  ACM Transactions on Knowledge Discovery from Data (TKDD) 13(5), 1–29 (2019)
  * Tsourakakis [2015] Tsourakakis, C.: The k-clique densest subgraph problem.  In: Proceedings of the 24th international conference on world wide web. pp. 1122–1132 (2015)


# Possible $P$\- wave $D_{s}\bar{D}_{s0}(2317)$ molecular state
$Y^{{}^{\prime}}(4274)$

Hong Qiang Zhu  College of Physics and Electronic Engineering, Chongqing
Normal University, Chongqing 401331,China  Yin Huang111corresponding author
<EMAIL_ADDRESS>Asia Pacific Center for Theoretical Physics, Pohang
University of Science and Technology, Pohang 37673, Gyeongsangbuk-do, South
Korea  School of Physical Science and Technology, Southwest Jiaotong
University, Chengdu 610031,China

###### Abstract

Stimulated by the measurement of the $J\psi\phi$ decay model of $Y(4274)$ by
the LHCb Collaboration, we consider a possible interpretation of this state as
a hadron molecular-a bound state of $D_{s}$ and $\bar{D}_{s0}(2317)$ mesons.
Using effective Lagrangian approach, we calculate the two-body strong decay
channels
$Y(4274)\to{}J/\psi\phi,\chi_{c0}\eta,\chi_{c0}\eta,D^{*}_{s}\bar{D}_{s}$,
$D\bar{D}^{*}$, $K\bar{K}^{*}$, and $\phi\phi$ through hadronic loops and
three-body decays into $\pi^{0}{}D_{s}\bar{D}_{s}$. In comparison with the
LHCb data, our results show that $Y(4274)$ cannot be assigned to be a
$D_{s}\bar{D}_{s0}(2317)$ molecular state. The calculated partial decay widths
with $J^{P}=1^{+}$ $D_{s}\bar{D}_{s0}$ molecular state picture indicate that
allowed decay modes, $\chi_{c0}\eta$ and $\chi_{c1}\eta$, may have the
smallest branching ratio and are of the order of 0.0 MeV. Future experimental
measurements of such two processes can be quite useful to test the different
interpretations of the $Y(4274)$. If $P-$wave $D_{s}\bar{D}_{s0}$ molecular
exist [we marked as $Y^{{}^{\prime}}(4274)$], the total decay is at the order
of 1.06-1.84 MeV, which seems to be within the reach of the current
experiments such as Belle II. In addition, the calculated partial decay widths
indicate that allowed decay mode, $D\bar{D}^{*}$, may have the biggest
branching ratio. The experimental measurements for this strong decay process
could be a crucial to observe such a new state $Y^{{}^{\prime}}(4274)$.

##  I Introduction

Thanks to the great progress of the experiment in the past several decades,
many hadrons that cannot be ascribed into the simple $\bar{q}q$ configuration
for mesons or $qqq$ configuration for baryons have been reported Zyla:2020zbs
. For example, various hidden-charm pentaquark were observed in the
$J/\psi{}p(\Lambda)$ invariant mass from the heavy baryon decay
$\Lambda_{b}^{+}\to{}K^{-}J/\psi{}p$ ($\Xi_{b}^{-}\to{}K^{-}J/\psi{}\Lambda$)
at LHCb Collaboration, $P_{c}(4312,4440,4450)$ and $P_{cs}(4459)$ Aaij:2019vzc
; LHCb:2020jpq . Their confirmation, and determination of their quantum
numbers, would allow new insights into the binding mechanisms present in
multiquark systems, and help improve understanding of QCD in the
nonperturbative regime.

In 2017, a charmonium-like meson named $Y(4274)$ was observed again by the
LHCb collaboration in the analysis of the $B^{+}\to{}J/\psi\phi{}K^{+}$
reaction LHCb:2016axx ; LHCb:2016nsl . The observed resonance masses, widths,
and favorable quantum numbers are

$\displaystyle M$ $\displaystyle=4273.3\pm{}8.3^{+17.2}_{-3.6}~{}~{}~{}{\rm
MeV},$ $\displaystyle\Gamma$ $\displaystyle=56\pm
11^{+8.0}_{-11}~{}~{}~{}~{}~{}{\rm MeV},~{}~{}J^{PC}=1^{++},$ (1)

respectively, which are consistent with the early CDF Collaboration Yi:2010aa
report

$\displaystyle M$ $\displaystyle=4274.4^{+8.4}_{-6.7}(stat)~{}~{}~{}{\rm
MeV},$ $\displaystyle\Gamma$
$\displaystyle=32.3^{+21.9}_{-15.3}(stat)~{}~{}~{}~{}~{}{\rm MeV}.$ (2)

Since the statistic is not enough, its spin-parity quantum number was not
confirmed by the CDF Collaboration. However, the isospin of this state is zero
and it contains at least four valence quarks from the observed $J/\psi{}\phi$
decay mode.

Following the discovery of the $Y(4274)$, several theoretical studies have
been performed. In the QCD sum rules approach, based on the analysis of the
mass spectrum, the $Y(4274)$ can be interpreted as the $S$-wave
$cs\bar{c}\bar{s}$ state with spin-parity $J^{P}=1^{+}$ Chen:2016oma ;
Wang:2016dcb . The compact tetraquark model, implemented by Stancu, can also
describe the $Y(4274)$ Stancu:2009ka , while only one $J^{P}=1^{+}$ state
exists. While in Ref. Lu:2016cwr  the mass of the $Y(4274)$ was studied in the
relativized quark model and it was shown that the $Y(4274)$ cannot be
explained as a tetraquark state, however, it can be a good candidate of the
conventional $\chi_{c1}(3^{3}P_{1})$ state. And the $\chi_{c1}(3^{3}P_{1})$
explain for $Y(4274)$ was also proposed in Ref. Gui:2018rvv . In the context
of the QCD two-point sum rule method by taking into account the quark, gluon,
and mixed vacuum condensates, Ref. Agaev:2017foq  assigned the $Y(4274)$ as a
sextet [$cs$][$\bar{c}\bar{s}$] diquark-antidiquark state with spin-parity
$J^{P}=1^{+}$. Based on the spin-spin interaction, Maiani et al. suggest that
the $Y(4274)$ may have quantum number $0^{++}$ or $2^{++}$ Maiani:2016wlq ,
which contradicts the experimental observation LHCb:2016axx ; LHCb:2016nsl .
Moreover, a detailed calculation is performed by Zhu Zhu:2016arf  where the
$Y(4274)$ may be described simultaneously by adding the up and down quark
components.

Although the studies of Refs. Chen:2016oma ; Wang:2016dcb ; Stancu:2009ka ;
Agaev:2017foq ; Lu:2016cwr ; Gui:2018rvv ; Maiani:2016wlq ; Zhu:2016arf  seem
to indicate that this state is a compact tetraquark state or diquark-
antidiquark state, $Y(4274)$ might still be a $D_{s}\bar{D}_{s0}(2317)$
hadronic molecule state. Since the mass of $Y(4274)$ is about 12 MeV below the
threshold of $D_{s}\bar{D}_{s0}$ ($m_{D_{s0}}=2317.8\pm{}0.6$ MeV and
$m_{D_{s}}=1968.34\pm{}0.07$ MeV Zyla:2020zbs ), it is reasonable to regard it
is a bound state of $D_{s}\bar{D}_{s0}$. The idea comes from molecular state
interpretation of deuteron due to deuteron mass is a little below the
corresponding threshold and exhibit sizable spatial extension. And because the
quantum numbers of $\bar{D}_{s0}$ and $D_{s}$ are $J^{P}=0^{+}$ and
$J^{P}=0^{-}$, respectively, to form a bound state with quantum number
$J^{P}=1^{+}$, the coupling between $Y(4274)$ and its constituents should be a
$P$ wave. Indeed, it is shown in Ref. He:2016pfa  that the interaction between
a $D_{s}$ meson and a $\bar{D}_{s0}$ meson is strong enough to form a bound
state with a mass about 4274 MeV.

From Ref He:2016pfa , $Y(4274)$ may be a molecular state. However, currently,
we cannot fully exclude other possible explanations such as a compact
pentaquark state Chen:2016oma ; Wang:2016dcb ; Stancu:2009ka ; Agaev:2017foq ;
Lu:2016cwr ; Gui:2018rvv ; Maiani:2016wlq ; Zhu:2016arf . Further research is
required to distinguish whether it is a molecular or compact multi-quark
state. One way to distinguish the two scenarios is to study the allowing
strong decay widths of the $Y(4274)$ baryon due to the strong decay almost
saturates the total strong decay width. In the present paper we consider
possible strong decay modes using an effective Lagrangian approach by assuming
that $Y(4274)$ is a hadronic molecule state of $D_{s}$ and $\bar{D}_{s0}$.

This work is organized as follows. The theoretical formalism is explained in
Sec. II. The predicted partial decay widths are presented in Sec. III,
followed by a short summary in the last section.

##  II FORMALISM AND INGREDIENTS

Besides the $J/\psi{}\phi$ decay model, which other decay is allowed. We first
find that the transition from $Y(4274)$ to final states composed of purely
neutral state $A\bar{A}$ are strictly forbidden by the conservation of the $c$
parity. Thus, the decay of $Y(4274)$ into $\eta_{c}\eta$, $\chi_{c0}\eta$,
$\chi_{c1}\eta$, $\phi\phi$, $\eta\eta$, $D^{*}_{s}\bar{D}_{s}$,
$D\bar{D}^{*}$, $K\bar{K}^{*}$ and $\pi^{0}\bar{D}_{s}D_{s}$ are allowed by
considering appropriate phase space Zyla:2020zbs . However, the transitions
$Y(4274)\to{}\eta_{c}\eta$ and $\eta\eta$ are also strictly forbidden by the
conservation of angular momentum. In this work, we will calculate
$J/\psi{}\phi$, $\chi_{c0}\eta$, $\chi_{c1}\eta$, $D^{*}_{s}\bar{D}_{s}$,
$D\bar{D}^{*}$, $K\bar{K}^{*}$ and $\phi\phi$ strong decay patterns of
$p$-wave $D_{s}\bar{D}_{s0}$ molecular state within the effective Lagrangians
approach, and find the relation between $D_{s}\bar{D}_{s0}$ molecular state
and $Y(4274)$ by comparing with the experiment observation.

Before introducing the theoretical framework, we need construct the flavor
functions for $D_{s}\bar{D}_{s0}$ system with definite $I(J^{P{\cal{C}}})$.
Since the $Y(4274)$ carry quantum numbers $I(J^{P{\cal{C}}})=0(1^{++})$, the
flavor function for a definite charge parity ${\cal{C}}=1$ can be easy
obtained Liu:2013rxa

$\displaystyle|D_{s}\bar{D}_{s0}\rangle=\frac{1}{\sqrt{2}}[D_{s}^{+}D_{s0}^{-}-D_{s}^{-}D_{s0}^{+}].$
(3)

Considering the quantum number $J^{P}=1^{+}$ and the flavor function,
$Y(4274)$ should couple to its components dominantly via $P-$ wave, and the
corresponding effective Lagrangian is in the form Ma:2010xx

$\displaystyle{\cal{L}}_{Y(4274)}(x)$
$\displaystyle=g_{YD_{s}\bar{D}_{s0}}Y^{\mu}(x)\int{}d^{4}y\Phi(y^{2})$
$\displaystyle\times{}\frac{1}{\sqrt{2}}[D^{+}_{s}(x+\omega_{D^{-}_{s0}}y)\overleftrightarrow{\partial_{\mu}}D^{-}_{s0}(x-\omega_{D^{+}_{s}}y)$
$\displaystyle-D^{-}_{s}(x+\omega_{D^{+}_{s0}}y)\overleftrightarrow{\partial_{\mu}}D^{+}_{s0}(x-\omega_{D^{-}_{s}}y)],$
(4)

where $\omega_{D_{s}}=m_{D_{s}}/(m_{D_{s}}+m_{D_{s0}})$ and
$\omega_{D_{s0}}=m_{D_{s0}}/(m_{D_{s}}+m_{D_{s0}})$. In the Lagrangian, an
effective correlation function $\Phi(y^{2})$ is introduced to describe the
distribution of the two constituents, the $D_{s}$ and the $\bar{D}_{s0}$, in
the hadronic molecular $Y(4274)$ state. The introduced correlation function
also makes the Feynman diagrams ultraviolet finite. Here we choose the Fourier
transformation of the correlation to be a Gaussian form in the Euclidean space
Faessler:2007gv ; Faessler:2007us ; Dong:2008gb ; Dong:2009uf ; Dong:2009yp ;
Dong:2017rmg ; Dong:2014ksa ; Dong:2014zka ; Dong:2013kta ; Dong:2013iqa ;
Dong:2013rsa ; Dong:2012hc ; Dong:2011ys ; Dong:2010xv ; Dong:2010gu ;
Dong:2009tg ; Dong:2017gaw ; Yang:2021pio ; Zhu:2020jke  ,

$\displaystyle\Phi(p^{2})\doteq\exp(-p_{E}^{2}/\Lambda^{2})$ (5)

with $\Lambda$ being the size parameter which characterizes the distribution
of the components inside the molecule. The value of $\Lambda$ could not be
determined from first principles, therefore it should better be determined by
experimental data. It is usually chosen to be about $1$ GeV, which depends on
experimental total widths Faessler:2007gv ; Faessler:2007us ; Dong:2008gb ;
Dong:2009uf ; Dong:2009yp ; Dong:2017rmg ; Dong:2014ksa ; Dong:2014zka ;
Dong:2013kta ; Dong:2013iqa ; Dong:2013rsa ; Dong:2012hc ; Dong:2011ys ;
Dong:2010xv ; Dong:2010gu ; Dong:2009tg ; Dong:2017gaw ; Yang:2021pio ;
Zhu:2020jke . In this work, we vary $\Lambda$ in a range of 0.9 GeV
$\leq\Lambda\leq{}1.10$ GeV to study whether the $Y(4274)$ can be interpreted
as $P-$ wave molecule composed of $D_{s}\bar{D}_{s0}$.

Figure 1: Self-energy of the $Y(4274)$ state.

The coupling constant $g_{YD_{s}\bar{D}_{s0}}$ in Eq. (4) can be computed by
the compositeness condition Weinberg:1962hj ; Salam:1962ap , which indicates
that the renormalization constants of a composite particle wave function
should be zero, i.e.,

$\displaystyle Z_{Y}=1-\frac{d\Sigma^{T}_{Y}}{dk_{0}}|_{k_{0}=m_{Y}}=0,$ (6)

where the $\Sigma^{T}_{Y}$ is the transverse part of the mass operator and
relates to its mass operator via the relation.

$\displaystyle\Sigma^{\mu\nu}_{Y}(p)=(g_{\mu\nu}-p^{\mu}p^{\nu}/p^{2})\Sigma^{T}_{Y}+....$
(7)

The concrete forms of the mass operators of the $Y(4274)$ corresponding to the
diagrams in Fig. 1 is

$\displaystyle\Sigma^{\mu\nu}_{Y}(k_{0})$
$\displaystyle=\frac{g^{2}_{YD_{s}\bar{D}_{s0}}}{2}\int\frac{d^{4}k_{1}}{(2\pi)^{4}}\\{\Phi^{2}[(k_{1}-k_{0}\omega_{D^{+}_{s}})^{2}]\frac{i}{k_{1}^{2}-m^{2}_{D^{+}_{s}}}$
$\displaystyle\times\frac{i}{(k_{1}-k_{0})^{2}-m^{2}_{D^{-}_{s0}}}+\Phi^{2}[(k_{1}-k_{0}\omega_{D^{-}_{s}})^{2}]\frac{i}{k_{1}^{2}-m^{2}_{D^{-}_{s}}}$
$\displaystyle\times\frac{i}{(k_{1}-k_{0})^{2}-m^{2}_{D^{+}_{s0}}}\\}(k_{0}-2k_{1})_{\mu}(k_{0}-2k_{1})_{\nu},$
(8)

where $k_{0}^{2}=m^{2}_{Y}$ with $k_{0},m_{Y}$ denoting the four momenta and
mass of the $Y(4274)$, respectively, $k_{1}$, $m_{D_{s}}$, and $m_{D_{s0}}$
are the four-momenta, the mass of the $D_{s}$ meson, and the mass of the
$\bar{D}_{s0}$ meson, respectively. With above preparations, we can obtain the
coupling constant of the $D_{s}\bar{D}_{s0}$ molecule to its components

$\displaystyle\frac{1}{g^{2}_{YD_{s}\bar{D}_{s0}}}$
$\displaystyle=\frac{m_{Y}}{8\pi^{2}{\cal{S}}}\int_{0}^{\infty}d\alpha\int_{0}^{\infty}d\beta\sum_{i=1}^{2}$
$\displaystyle\times\exp[-\frac{1}{\Lambda^{2}}(-{\cal{F}}_{i}m^{2}_{Y}+{\cal{H}}_{i}+\frac{{\cal{C}}_{i}^{2}m_{Y}^{2}}{z})](\frac{{\cal{F}}_{i}}{z^{3}}-\frac{{\cal{C}}^{2}_{i}}{z^{4}})$
(9)

where
${\cal{F}}_{1}=(2\omega^{2}_{D_{s}^{+}}+\beta),{\cal{F}}_{2}=(2\omega^{2}_{D_{s}^{-}}+\beta),{\cal{H}}_{1}=\alpha{}m^{2}_{D_{s}^{+}}+\beta{}m^{2}_{D^{-}_{s0}}$,
${\cal{H}}_{2}=\alpha{}m^{2}_{D_{s}^{-}}+\beta{}m^{2}_{D^{+}_{s0}}$,
${\cal{C}}_{1}=(2\omega_{D_{s}^{+}}+\beta)$,${\cal{C}}_{2}=(2\omega_{D_{s}^{-}}+\beta)$,
$z=2+\alpha+\beta$, and ${\cal{S}}=1.0$ GeV.

###  II.1 The decay $Y(4274)\to{}J/\psi{}\phi$

Since the $Y(4274)$ was observed in the $J/\psi{}\phi$ invariant mass, we
first calculate the $J/\psi{}\phi$ two-body decay width of the $Y(4274)$ via
the triangle diagrams show in Fig. 2.

Figure 2: Feynman diagrams for the $Y(4274)\to{}J/\psi\phi$ decay processes.
We also show the definitions of the
kinematics($k_{0},k_{1},k_{2},p_{1},p_{2},q$,and $q^{{}^{\prime}}$) used in
the calculation.

The hadronic decay of the $D_{s}\bar{D}_{s0}$ molecular state into
$J/\psi\phi$ mediated by the exchange of the $D_{s}^{*}$ meson. To evaluate
the diagrams, in addition to the Lagrangian in Eq. 4, the following effective
Lagrangians, responsible for vector meson $V(=J/\psi,\phi)$ coupling to
$D_{s}^{*}\bar{D}_{s0}$ are needed as well Ma:2010xx

$\displaystyle{\cal{L}}_{D_{s}^{*}\bar{D}_{s0}V}=g_{V\bar{D}_{s0}D_{s}^{*}}[D^{-}_{s0}D^{*+}_{s;\mu\nu}-D^{+}_{s0}D^{*-}_{s;\mu\nu}]V^{\mu\nu},$
(10)

where $V^{\mu\nu}=\partial^{\mu}V^{\nu}-\partial^{\nu}V^{\mu}$. The coupling
constant $g_{J/\psi\bar{D}_{s0}D_{s}^{*}}=0.225$ GeV-1 and
$g_{\phi\bar{D}_{s0}D_{s}^{*}}=0.135$ GeV-1 are got from Ref. Ma:2010xx .

To compute the $D_{s}D_{s}^{*}V$ vertices, we also need the following
effective Lagrangian Oset:2002sh ; Bramon:1992kr

$\displaystyle{\cal{L}}_{VVP}=\frac{G}{\sqrt{2}}\epsilon^{\mu\nu\alpha\beta}\langle\partial_{\mu}V_{\nu}\partial_{\alpha}V_{\beta}P\rangle,$
(11)

where $G=3h^{2}/(4\pi^{2}f)$ with $h=-G_{V}m_{\rho}/(\sqrt{2}f^{2})$,
$f=0.093$ GeV, $G_{V}=0.069$ GeV Oset:2002sh ; Bramon:1992kr , and
$m_{\rho}=0.775$ GeV. $V_{\mu}$ and $P$ are standard $SU(4)$ matrices
constructed with 16-plet of the vector meson contain $\rho$ and the 16-plet of
pseudoscalar mesons containing the $\pi$, respectively.

$V_{\mu}=\left(\begin{array}[]{cccc}\frac{1}{\sqrt{2}}(\rho^{0}+\omega)&\rho^{+}&K^{*+}&\bar{D}^{*0}\\\
\rho^{-}&\frac{1}{\sqrt{2}}(-\rho^{0}+\omega)&K^{*0}&-D^{*-}\\\
K^{*-}&\bar{K}^{*0}&\phi&D^{*-}_{s}\\\ D^{*0}&-D^{*+}&D^{*+}_{s}&J/\psi\\\
\end{array}\right)_{\mu},$ (12)
$P=\sqrt{2}\left(\begin{array}[]{cccc}\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}+\frac{\eta^{{}^{\prime}}}{\sqrt{3}}&\pi^{+}&K^{+}&\bar{D}^{0}\\\
\pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}+\frac{\eta^{{}^{\prime}}}{\sqrt{3}}&K^{0}&-D^{-}\\\
K^{-}&\bar{K}^{0}&-\sqrt{\frac{2}{3}}\eta+\frac{1}{\sqrt{3}}\eta^{{}^{\prime}}&D^{-}_{s}\\\
D^{0}&-D^{+}&D^{+}_{s}&\eta_{c}\\\ \end{array}\right).$ (13)

With above effective Lagrangians, we can get the decay amplitudes
corresponding to the diagrams in Fig. 2

$\displaystyle{\cal{M}}_{a}$
$\displaystyle=\frac{G}{\sqrt{2}}g_{J/\psi{}\bar{D}_{s0}D^{*}_{s}}g_{YD_{s}\bar{D}_{s0}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{-}_{s}}-k_{2}\omega_{D_{s0}^{+}})^{2}]$
$\displaystyle\times{}(p_{1\mu}g_{\nu\sigma}-p_{1\nu}g_{\mu\sigma})(q_{\mu}g_{\nu\eta}-q_{\nu}g_{\mu\eta})(g_{\eta\lambda}-\frac{q_{\eta}q_{\lambda}}{m^{2}_{D^{*+}_{s}}})$
$\displaystyle\times{}\epsilon^{\tau\lambda\alpha\beta}q_{\tau}p_{2\alpha}(k_{1}^{\varphi}-k_{2}^{\varphi})\epsilon^{Y}_{\varphi}(k_{0})\epsilon^{*\sigma}_{J/\psi}(p_{1})\epsilon^{*\phi}_{\beta}(p_{2})$
$\displaystyle\times{}\frac{1}{q^{2}-m^{2}_{D^{*+}_{s}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{-}_{s0}}},$
(14) $\displaystyle{\cal{M}}_{b}$
$\displaystyle=\frac{G}{\sqrt{2}}g_{\phi\bar{D}_{s0}D^{*}_{s}}g_{YD_{s}\bar{D}_{s0}}\int\frac{d^{4}q^{{}^{\prime}}}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{-}_{s}}-k_{2}\omega_{D_{s0}^{+}})^{2}]$
$\displaystyle\times{}(p_{2\lambda}g_{\sigma\tau}-p_{2\sigma}g_{\lambda\tau})(q^{{}^{\prime}}_{\lambda}g_{\sigma\eta}-q^{{}^{\prime}}_{\sigma}g_{\lambda\eta})(g_{\nu\eta}-\frac{q^{{}^{\prime}}_{\nu}q^{{}^{\prime}}_{\eta}}{m^{2}_{D^{*+}_{s}}})$
$\displaystyle\times{}\epsilon^{\mu\nu\alpha\beta}q^{{}^{\prime}}_{\mu}p_{1\alpha}(k_{1}^{\varphi}-k_{2}^{\varphi})\epsilon^{Y}_{\varphi}(k_{0})\epsilon^{*\phi}_{\tau}(p_{2})\epsilon^{*J/\psi}_{\beta}(p_{1})$
$\displaystyle\times{}\frac{1}{q^{{}^{\prime}2}-m^{2}_{D^{*+}_{s}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{-}_{s0}}},$
(15) $\displaystyle{\cal{M}}_{c}$
$\displaystyle={\cal{M}}_{a}(D^{+}_{s0}\to{}D^{-}_{s0},D^{-}_{s}\to{}D^{+}_{s},D^{*+}_{s}\to{}D^{*-}_{s}),$
(16) $\displaystyle{\cal{M}}_{d}$
$\displaystyle={\cal{M}}_{b}(D^{+}_{s0}\to{}D^{-}_{s0},D^{-}_{s}\to{}D^{+}_{s},D^{*+}_{s}\to{}D^{*-}_{s}).$
(17)

###  II.2 The decay $Y(4274)\to{}\chi_{c0}\eta$ and $\chi_{c1}\eta$

In this section, we compute the other possible decays of $Y(4274)$ with
$D_{s}\bar{D}_{s0}$ molecular. Fig. 4 shows the hadronic decay of the
$D_{s}\bar{D}_{s0}$ molecular state into $\chi_{c0}\eta$ and $\chi_{c1}\eta$
mediated by the exchange of $D_{s}$ meson.

Figure 3: Diagrams for $Y(4274)\to{}\eta\chi_{c0/c1}$ decay processes.

The ingredients need are the $\chi_{c0}D^{+}_{s}D^{-}_{s}$ and
$\chi_{c1}D^{+}_{s}D^{-}_{s}$ Lagrangians Colangelo:2003sa ; Casalbuoni:1996pg

$\displaystyle{\cal{L}}_{\chi_{c1}D^{+}_{s}D^{-}_{s}}=-ig_{\chi_{c1}D^{+}_{s}D^{-}_{s}}(D_{s}^{+}\partial_{\mu}D_{s}^{-}-\partial_{\mu}D_{s}^{+}D_{s}^{-})\chi_{c1}^{\mu},$
(18)
$\displaystyle{\cal{L}}_{\chi_{c0}D_{s}^{+}D_{s}^{-}}=-g_{\chi_{c0}D_{s}^{+}D_{s}^{-}}\chi_{c0}D_{s}^{+}D_{s}^{-}.$
(19)

The coupling constants $g_{\chi_{c1}D^{+}_{s}D^{-}_{s}}$ and
$g_{\chi_{c0}D_{s}^{+}D_{s}^{-}}$ can be fixed from the heavy quark field
theory Casalbuoni:1996pg .

Moreover, the effective Lagrangian, responsible for the coupling of $D_{s0}$
to $\eta{}D_{s}$, is needed as well Huang:2019qmw

$\displaystyle{\cal{L}}_{D_{s0}D_{s}\eta}=g_{D_{s0}D_{s}\eta}D_{s0}D_{s}\eta,$
(20)

where the coupling constant is found to be $g_{D_{s0}D_{s}\eta}=6.40$
Gamermann:2006nm .

Thus, we can obtain the following amplitudes for the decays
$Y(4274)\to{}\chi_{c0}\eta$

$\displaystyle{\cal{M}}_{a}$
$\displaystyle=\frac{g_{\eta\bar{D}_{s0}D_{s}}g_{YD_{s}\bar{D}_{s0}}g_{D_{s}^{-}D_{s}^{+}\chi_{c0}}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{-}_{s}}-k_{2}\omega_{D_{s0}^{+}})^{2}]$
$\displaystyle\times{}(k_{1}^{\varphi}-k_{2}^{\varphi})\epsilon^{Y}_{\varphi}(k_{0})\frac{1}{q^{2}-m^{2}_{D_{s}^{+}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{+}_{s0}}},$
(21) $\displaystyle{\cal{M}}_{b}$
$\displaystyle=-\frac{g_{\eta\bar{D}_{s0}D_{s}}g_{YD_{s}\bar{D}_{s0}}g_{D_{s}^{-}D_{s}^{+}\chi_{c0}}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{+}_{s}}-k_{2}\omega_{D_{s0}^{-}})^{2}]$
$\displaystyle\times{}(k_{1}^{\varphi}-k_{2}^{\varphi})\epsilon^{Y}_{\varphi}(k_{0})\frac{1}{q^{2}-m^{2}_{D_{s}^{-}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{+}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{-}_{s0}}},$
(22)

and $Y(4274)\to{}\chi_{c1}\eta$

$\displaystyle{\cal{M}}_{a}$
$\displaystyle=\frac{g_{\eta\bar{D}_{s0}D_{s}}g_{YD_{s}\bar{D}_{s0}}g_{D_{s}^{-}D_{s}^{+}\chi_{c1}}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{-}_{s}}-k_{2}\omega_{D_{s0}^{+}})^{2}]$
$\displaystyle\times{}(k_{1}^{\varphi}-k_{2}^{\varphi})\epsilon^{Y}_{\varphi}(k_{0})(k_{2}^{\mu}-q^{\mu})\epsilon^{\chi_{c1}}_{\mu}(p_{2})$
$\displaystyle\times\frac{1}{q^{2}-m^{2}_{D_{s}^{+}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{+}_{s0}}},$
(23) $\displaystyle{\cal{M}}_{b}$
$\displaystyle=-\frac{g_{\eta\bar{D}_{s0}D_{s}}g_{YD_{s}\bar{D}_{s0}}g_{D_{s}^{-}D_{s}^{+}\chi_{c1}}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{+}_{s}}-k_{2}\omega_{D_{s0}^{-}})^{2}]$
$\displaystyle\times{}(k_{1}^{\varphi}-k_{2}^{\varphi})\epsilon^{Y}_{\varphi}(k_{0})(k_{2}^{\mu}-q^{\mu})\epsilon^{\chi_{c1}}_{\mu}(p_{2})$
$\displaystyle\times\frac{1}{q^{2}-m^{2}_{D_{s}^{-}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{+}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{-}_{s0}}}.$
(24)

The minus sign in Eqs. (27) and (24) come from the flavor function that shown
in Eq. (3). Then, we find
${\cal{M}}[Y(4274)\to{}\chi_{c0/c1}\eta]_{Total}={\cal{M}}_{a}+{\cal{M}}_{b}=0$.

###  II.3 The decay $Y(4274)\to{}D^{*}_{s}\bar{D}_{s}$

Fig. 4 shows the hadronic decay of the $D_{s}\bar{D}_{s0}$ molecular state
into $D^{*}_{s}\bar{D}_{s}$ mediated by the exchange of $\eta$ meson.

Figure 4: Diagrams for $Y(4274)$ decay processes:(a,d)
$Y(4274)\to{}D^{*+}_{s}D^{-}_{s}$, and (b,c)
$Y(4274)\to{}D^{*-}_{s}D^{+}_{s}$.

The ingredients need are the vector($V$)-pseudoscalar($P$)-pseudoscalar($P$)
Lagrangian

$\displaystyle{\cal{L}}_{VPP}=-ig\langle[P,\partial_{\mu}P]V^{\mu}\rangle,$
(25)

The coupling $g$ is fixed from the strong decay width of $K^{*}\to{}K\pi$.
With the help of Eq. (25), the two-body decay width
$\Gamma(K^{*+}\to{}K^{0}\pi^{+})$ is related to $g$ as

$\displaystyle\Gamma(K^{*+}\to{}K^{0}\pi^{+})=\frac{g^{2}}{6\pi{}m^{2}_{K^{*+}}}{\cal{P}}^{3}_{\pi{}K^{*}}=\frac{2}{3}\Gamma_{K^{*+}},$
(26)

where ${\cal{P}}_{\pi{}K^{*}}$ is the three-momentum of the $\pi$ in the rest
frame of the $K^{*}$. Using the experimental strong decay width
($\Gamma_{K^{*+}}=50.3\pm{}0.8$ MeV) and the masses of the particles needed in
the present work Zyla:2020zbs  we obtain $g=4.61$.

Thus, we can obtain the following amplitudes for the decay
$Y(4274)\to{}D^{*}_{s}\bar{D}_{s}$

$\displaystyle{\cal{M}}_{a}$
$\displaystyle=\frac{-gg_{V\bar{D}_{s0}D_{s}^{*}}g_{YD_{s}\bar{D}_{s0}}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{-}_{s}}-k_{2}\omega_{D_{s0}^{+}})^{2}]$
$\displaystyle\times{}(p_{1\mu}g_{\nu\eta}-p_{1\nu}g_{\mu\eta})(q_{\mu}g_{\nu\lambda}-q_{\nu}g_{\mu\lambda})(g_{\lambda\alpha}-\frac{q_{\lambda}q_{\alpha}}{m^{2}_{V}})$
$\displaystyle\times{}(k_{2\alpha}+p_{2\alpha})(k_{1\omega}-k_{2\omega})\epsilon^{*\eta}_{D^{*+}_{s}}(p_{1})\epsilon^{\omega}_{Y}(k_{0})$
$\displaystyle\times{}\frac{1}{q^{2}-m^{2}_{V}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{+}_{s0}}},$
(27) $\displaystyle{\cal{M}}_{b}$
$\displaystyle={\cal{M}}_{b}(D^{+}_{s0}\to{}D^{-}_{s0},D^{-}_{s}\to{}D^{+}_{s},D^{*+}_{s}\to{}D^{*-}_{s}),$
(28) $\displaystyle{\cal{M}}_{c}$
$\displaystyle={\cal{M}}_{f}(D^{-}_{s0}\to{}D^{+}_{s0},D^{+}_{s}(k_{2})\to{}D^{-}_{s}(k_{2}),D^{*+}_{s}\to{}D^{*-}_{s})$
(29) $\displaystyle{\cal{M}}_{d}$
$\displaystyle=-\frac{gg_{YD_{s}\bar{D}_{s0}}g_{D_{s0}D_{s}\eta}}{\sqrt{6}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{+}_{s}}-k_{2}\omega_{D_{s0}^{-}})^{2}]$
$\displaystyle\times{}(q_{\mu}-k_{2\mu})(k_{1\nu}-k_{2\nu})\epsilon^{Y}_{\nu}(k_{0})\epsilon^{*\mu}_{D^{*+}_{s}}(p_{1})$
$\displaystyle\times\frac{1}{q^{2}-m^{2}_{\eta}}\frac{1}{k_{2}^{2}-m^{2}_{D^{+}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{-}_{s0}}},$
(30)

###  II.4 The decay $Y(4274)\to{}D_{s}\bar{D}_{s}\pi^{0}$

Now we turn to $D_{s}\bar{D}_{s}\pi^{0}$ three-body decay channel of
$Y(4274)$. Under the $D_{s}\bar{D}_{s0}$ molecular state assignment, $Y(4274)$
first dissociates into $D_{s}^{+}D_{s0}^{-}$ or $D_{s}^{-}D_{s0}^{+}$. Then,
the decay $Y(4274)\to{}D_{s}^{+}D_{s}^{-}\pi^{0}$ occur via the transitions
$D_{s0}^{\pm}\to{}D^{\pm}_{s}\pi^{0}$, where $D_{s0}^{\pm}$ decay into
$D^{\pm}_{s}\pi^{0}$ by considering the $\eta-\pi^{0}$ mixing mechanism
BaBar:2003oey ; Liu:2006jx . The relevant Feynman diagrams ars shown in Fig.
5.

Figure 5: Feynman diagrams for $Y(4274)\to{}D_{s}\bar{D}_{s}\pi^{0}$ decay
processes. We also show the definitions of the
kinematics($k_{0},p_{1},p_{2},p_{3}$) used in the calculation.

The Lagrangian including $\eta-\pi^{0}$ mixing for $\pi^{0}D_{s0}D_{s}$ have
been constructed in Ref. He:2011ed  and in the form

$\displaystyle{\cal{L}}_{\pi^{0}D_{s0}D_{s}}$
$\displaystyle=g_{\pi^{0}D_{s0}D_{s}}\pi^{0}D_{s0}D_{s},$ (31)

where the coupling constant $g_{\pi^{0}D_{s0}D_{s}}$ can be extracted by the
relation

$\displaystyle\Gamma(D_{s0}^{\pm}\to{}D_{s}^{\pm}\pi^{0})=g^{2}_{\pi^{0}D_{s0}D_{s}}\frac{|\vec{{\cal{P}}}_{\pi^{0}}|}{8\pi{}m^{2}_{D^{\pm}_{s0}}}.$
(32)

In above, ${\cal{P}}_{\pi^{0}}$ is the three-momentum of $\pi$ in the rest
frame of $D_{s0}$.

At present, the experimental partial decay width on
$D_{s0}^{\pm}\to{}D_{s}^{\pm}\pi^{0}$ is absent. We only find the absolute
branching fraction $B(D_{s0}^{\pm}\to{}D_{s}^{\pm}\pi^{0})$ is measured as
$1.00^{+0.00}_{-0.14}$(stat)${}^{+0.00}_{-0.14}$(syst) BESIII:2017vdm  and the
upper limit on $D_{s0}$ width is 3.8 MeV at the 95$\%$ confidence level (C.L.)
BaBar:2006eep . In this work, we focus on whether $Y(4274)$ can be a $P$\-
wave $D_{s}\bar{D}_{s0}$ molecular state and hope that the theoretical maximum
decay width cannot be compared with the experimental data. Thus, we use upper
limit $\Gamma{(D_{s0})}=3.8$ MeV to determine the coupling constant
$g_{\pi^{0}D_{s0}D_{s}}=1.3124$ GeV and in further calculation.

The general expression of the decay amplitude
$Y(4274)\to{}D_{s}\bar{D}_{s}\pi^{0}$ is

$\displaystyle{\cal{M}}$
${}_{a/b}=\frac{g_{\pi^{0}D_{s0}D_{s}}g_{YD_{s}\bar{D}_{s0}}}{\sqrt{2}}(q_{\mu}-p_{1\mu})\epsilon^{\mu}_{Y}(k_{0})$
$\displaystyle\times\Phi[(p_{1}\omega_{\bar{D}_{s0}}-q\omega_{D_{s}})^{2}]\frac{1}{q^{2}-m^{2}_{D_{s0}}+im_{D_{s0}}\Gamma_{D_{s0}}},$
(33)

where $q=p_{2}+p_{3}$.

###  II.5 The decay $Y(4274)\to{}D\bar{D}^{*}$, $K\bar{K}^{*}$, and $\phi\phi$

In this section, we calculate the two-body Okubo-Zweig-Iizuka (OZI) allowed
strong decays $Y(4274)\to{}D\bar{D}^{*}$, $K\bar{K}^{*}$, and $\phi\phi$. The
processes is described as a quark-antiquark pair $c\bar{c}$ or $s\bar{s}$
annihilation in the initial state. Meanwhile, a light quark-antiquark pair is
created and then regroups into two outgoing hadrons by a quark rearrangement
process. The decay of $Y(4274)$ into other channels, such as
$Y(4274)\to{}\pi\bar{\pi}$, are ignored because in these processes, which
involve the creation or annihilation of two $\bar{q}q$ $(q=u,d,s,c)$ quark
pairs, are usually strongly suppressed. The quark-level diagrams are depicted
in Fig. 6.

Figure 6: The $Y(4274)$ decay process via the OZI mechanism. $q=u,d$ quarks
for diagram (a) and $q=u,d,s$ quarks for diagram (b).

And corresponding hadron-level diagrams are in Fig. 7.

Figure 7: Diagrams for $Y(4274)$ decay processes:(a,b)
$Y(4274)\to{}D\bar{D}^{*}$,(c,d) $Y(4274)\to{}K\bar{K}^{*}$, and (e,f)
$Y(4274)\to{}\phi\phi$.

Besides the Lagrangians above, the effective Lagrangian of $KDD_{s0}$ vertex
is also needed Huang:2019qmw

$\displaystyle{\cal{L}}_{KDD_{s0}}$ $\displaystyle=g_{KDD_{s0}}KDD_{s0},$ (34)

where the coupling constant
$g_{KDD_{s0}}=\sqrt{2}g_{K^{-}D^{0}D^{-}_{s0}}=10.21$ is obtained from the
coupling of the $D_{s0}$ to the $DK$ channel in isospin $I=0$ Gamermann:2006nm
. Putting all pieces together, we obtain the amplitudes for
$Y(4274)\to{}D\bar{D}^{*}$, $K\bar{K}^{*}$, and $\phi\phi$ which correspond to
the diagrams in Fig. 7,

$\displaystyle{\cal{M}}_{a}$
$\displaystyle=\frac{g_{KDD_{s0}}gg_{Y}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{-}_{s}}-k_{2}\omega_{D_{s0}^{+}})^{2}]$
$\displaystyle\times{}(q_{\mu}-k_{2\mu})(k_{1\nu}-k_{2\nu})\epsilon_{Y}^{\nu}(k_{0})\epsilon^{*\mu}_{D^{*-(0)}}(p_{2})$
$\displaystyle\times\frac{1}{q^{2}-m^{2}_{K^{0(+)}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{+}_{s0}}},$
(35) $\displaystyle{\cal{M}}_{b}$
$\displaystyle=-\frac{g_{KDD_{s0}}gg_{Y}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{+}_{s}}-k_{2}\omega_{D_{s0}^{-}})^{2}]$
$\displaystyle\times{}(q_{\mu}-k_{2\mu})(k_{1\nu}-k_{2\nu})\epsilon_{Y}^{\nu}(k_{0})\epsilon^{*\mu}_{D^{*+(0)}}(p_{2})$
$\displaystyle\times\frac{1}{q^{2}-m^{2}_{\bar{K}^{0(-)}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{+}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{-}_{s0}}},$
(36) $\displaystyle{\cal{M}}_{c}$
$\displaystyle={\cal{J}}\frac{g_{KDD_{s0}}gg_{Y}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{-}_{s}}-k_{2}\omega_{D_{s0}^{+}})^{2}]$
$\displaystyle\times{}(q_{\mu}-k_{2\mu})(k_{1\nu}-k_{2\nu})\epsilon_{Y}^{\nu}(k_{0})\epsilon^{*\mu}_{K^{*-(0)}}(p_{2})$
$\displaystyle\times\frac{1}{q^{2}-m^{2}_{D^{0(+)}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{+}_{s0}}},$
(37) $\displaystyle{\cal{M}}_{d}$
$\displaystyle={\cal{J}}\frac{g_{KDD_{s0}}gg_{Y}}{\sqrt{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{+}_{s}}-k_{2}\omega_{D_{s0}^{-}})^{2}]$
$\displaystyle\times{}(q_{\mu}-k_{2\mu})(k_{1\nu}-k_{2\nu})\epsilon_{Y}^{\nu}(k_{0})\epsilon^{*\mu}_{K^{*+(0)}}(p_{2})$
$\displaystyle\times\frac{1}{q^{2}-m^{2}_{D^{0(-)}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{+}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{-}_{s0}}},$
(38) $\displaystyle{\cal{M}}_{e}$
$\displaystyle=-\frac{Gg_{Y}g_{\phi\bar{D}_{s0}D_{s}^{*}}}{2}\int\frac{d^{4}q}{(2\pi)^{4}}\Phi[(k_{1}\omega_{D^{+}_{s}}-k_{2}\omega_{D_{s0}^{-}})^{2}]$
$\displaystyle\times{}(q^{\mu}g^{\nu\eta}-q^{\nu}g^{\mu\eta})(p_{1\mu}g_{\nu\rho}-p_{1\nu}g_{\mu\rho})(-g^{\eta\lambda}+q^{\eta}q^{\lambda}/m^{2}_{D_{s}^{*+}})$
$\displaystyle\times{}\epsilon^{\tau\lambda\alpha\beta}q_{\tau}p_{2\alpha}(k_{1\sigma}-k_{2\sigma})\epsilon_{\phi}^{*\rho}(p_{1})\epsilon_{\phi}^{*\beta}(p_{2})\epsilon_{Y}^{\sigma}(k_{0})$
$\displaystyle\times{}\frac{1}{q^{2}-m^{2}_{D_{s}^{*+}}}\frac{1}{k_{2}^{2}-m^{2}_{D^{-}_{s}}}\frac{1}{k_{1}^{2}-m^{2}_{D^{+}_{s0}}},$
(39) $\displaystyle{\cal{M}}_{f}$
$\displaystyle={\cal{M}}_{e}(D_{s0}^{+}\to{}D_{s0}^{-},D_{s}^{-}\to{}D_{s}^{+},D_{s}^{*+}\to{}D_{s}^{*-}),$
(40)

where ${\cal{J}}$=1 and $-1$ are for $D^{0}(D^{0})$ and $D^{+}(D^{-})$
exchange, respectively.

Once the amplitudes are determined, the corresponding partial decay widths can
be obtained, which reads,

$\displaystyle\Gamma(Y(4274)\to MB)$
$\displaystyle=\frac{1}{24\pi}\frac{|\vec{p}_{1}|}{m^{2}_{Y}}\overline{|{\cal{M}}|^{2}},$
(41) $\displaystyle\Gamma(Y(4274)\to D_{s}\bar{D}_{s}\pi^{0})$
$\displaystyle=\frac{1}{3(2\pi)^{5}}\frac{1}{16m^{2}_{Y}}\overline{|{\cal{M}}|^{2}}|\vec{p}^{*}_{3}||\vec{p}_{1}|$
$\displaystyle\times{}dm_{\pi^{0}D_{s}^{\pm}}d\Omega^{*}_{p_{3}}d\Omega_{p_{1}},$
(42)

where $J$ is the total angular momentum of the Y(4274) state, the
$|\vec{p}_{1}|$ is the three-momenta of the decay products in the center of
mass frame, the overline indicates the sum over the polarization vectors of
the final hadrons, and $MB$ denotes the decay channel of $MB$, i.e.,
$J/\psi{}\phi$, $\chi_{c0}\eta$, $\chi_{c1}\eta$, $D^{*}_{s}\bar{D}_{s}$,
$D\bar{D}^{*}$, $K\bar{K}^{*}$ and $\phi\phi$. In Eq. 42, the
$\vec{p}^{*}_{3}$ and $\Omega^{*}_{p_{3}}$ are the momentum and angle of the
particle $D^{\pm}_{s}$ in the rest frame of $D^{\pm}_{s}$ and $\pi^{0}$,
respectively, and $\Omega_{p_{1}}$ is the angle of $D_{s}^{\pm}$ in the rest
frame of the decaying particle $Y(4274)$. $m_{\pi^{0}D_{s}^{\pm}}$ is the
invariant mass for $D^{\pm}_{s}$ and $\pi^{0}$ and must meet
$m_{D^{\pm}_{s0}}+m_{\pi^{0}}\leq{}m_{\pi^{0}D_{s}^{\pm}}\leq{}m_{Y}-m_{D^{\pm}_{s}}$.

##  III RESULTS and Discussions

In this work, we study the strong decays of the $Y(4274)$ to the two-body
final states $J/\psi\phi,\chi_{c0}\eta,\chi_{c1}\eta,D^{*}_{s}\bar{D}_{s}$,
$D\bar{D}^{*}$, $K\bar{K}^{*}$, $\phi\phi$ and three-body decay into
$\pi^{0}D_{s}\bar{D}_{s}$ assuming that $Y(4274)$ is a
$D_{s}\bar{D}_{s0}(2317)$ molecular state. In order to obtain the decay width
shown in Figs. 2, 4, and 5, the coupling constant $g_{YD_{s}\bar{D}_{s0}}$
should be computed first.

Figure 8: The $\Lambda$ dependence of the coupling $g_{YD_{s}\bar{D}_{s0}}$
estimated from compositeness condition.

According to the compositeness condition that we introduced in Eq. 6,
$\Lambda$ dependence of the coupling constant $g_{YD_{s}\bar{D}_{s0}}$ is
computed. With a value of cutoff $\Lambda=0.9-1.1$ GeV, the corresponding
coupling constants shown in Fig 8. We note that they decrease slowly with the
increase of cut-off, and the coupling constant is almost independent of
$\Lambda$, where $Y(4274)$ is $P-$wave $D_{s}\bar{D}_{s0}$ molecular state.
According to the studies in Refs Faessler:2007gv ; Faessler:2007us ;
Dong:2008gb ; Dong:2009uf ; Dong:2009yp ; Dong:2017rmg ; Dong:2014ksa ;
Dong:2014zka ; Dong:2013kta ; Dong:2013iqa ; Dong:2013rsa ; Dong:2012hc ;
Dong:2011ys ; Dong:2010xv ; Dong:2010gu ; Dong:2009tg ; Dong:2017gaw ;
Yang:2021pio ; Zhu:2020jke , a typical value of $\Lambda=1.0$ GeV is often
employed. Thus, in this work we take $\Lambda$ = 1.0 GeV and the corresponding
coupling constants are $g_{YD_{s}\bar{D}_{s0}}=13.34^{+1.11}_{-0.89}$ GeV,
which the error reflects variation in $\Lambda$ from $0.9$ to $1.1$ GeV.

Figure 9: (Color line) Partial decay widths of the $Y(4274)\to{}J/\psi{}\phi$
(black solid line), $Y(4274)\to{}\pi^{0}D_{s}\bar{D}_{s}$ (blue dot line),
$Y(4274)\to{}D_{s}^{*}\bar{D}_{s}$ (red dash line),
$Y(4274)\to{}\chi_{c0}\eta$ (violet short dot line),
$Y(4274)\to{}D\bar{D}^{*}$ (magenta dash dot line), $Y(4274)\to{}K\bar{K}^{*}$
(olive dash dot dot line), $Y(4274)\to{}\chi_{c1}\eta$ (purple short dash dot
line), and $Y(4274)\to{}\phi\phi$ (navy short dash line).

Once the coupling constant $g_{YD_{s}\bar{D}_{s0}}=13.34^{+1.11}_{-0.89}$ GeV
are determined, the decay widths of $Y(4274)$ can be calculated
straightforwardly. In Fig. 9, we show the partial decay widths of
$Y(4274)\to{}J/\psi\phi,\chi_{c0}\eta,\chi_{c1}\eta,D^{*}_{s}\bar{D}_{s}$,
$\pi^{0}D_{s}\bar{D}_{s}$, $D\bar{D}^{*}$, $K\bar{K}^{*}$, and $\phi\phi$ as a
function of cutoff parameter $\Lambda$. We find that the estimated two-body
decay widths increases with increase of cut-off and are all insensitive to
cut-off parameter $\Lambda$. While the $Y(4274)\to{}\pi^{0}D_{s}\bar{D}_{s}$
three-body decay decreases, but very slowly. We also find that the partial
decay width is the largest for transition $Y(4274)\to{}D\bar{D}^{*}$. Thanks
the flavor symmetry of the wave function that shown in Eq. (3), the
$Y(4274)\to{}\chi_{c0}\eta$ and $\chi_{c1}\eta$ two-body decay widths are of
the order of about 0.0 MeV. We also note that the three-body transition
strength is quite small, and the decay width is of the order of 0.158-0.190
MeV. The small three-body decay width can be easy understand due to the decay
$Y(4274)\to D_{s}\bar{D}_{s}\pi^{0}$ is an isospin violation process.

Two reasons can help us to understand why transition
$Y(4274)\to{}D\bar{D}^{*}$ provides the dominant contribution. First is that
transition $Y(4274)\to{}D\bar{D}^{*}$ is $s$-wave decay, and the lowest
angular momentum gives the dominant contribution. Larger
$Y(4274)\to{}D\bar{D}^{*}$ decay can also be understood due to the main
component of $D_{s0}$ is $DK$, and the coupling constant related to this
vertex is larger than the others. The same $D_{s0}DK$ coupling also exist in
the $Y(4274)\to{}K\bar{K}^{*}$ reaction. However, its partial decay width is
small. A possible explanation is that a light quark-antiquark pair creation or
annihilation in Fig. 6 is easier than that of a heavy quark-antiquark pair.

Figure 10: (Color line) The total decay width of the $Y(4274)$. The cyan bands
denote the experimental total width LHCb:2016axx ; LHCb:2016nsl .

We also show the dependence of the total decay width on cutoff $\Lambda$ in
Fig. 10. In present calculation, we vary $\Lambda$ from 0.9 to 1.1 GeV. In
this $\Lambda$ range, the total decay width increases, and predicted decay
width $\Gamma_{Y}$=1.25-2.0 MeV is much smaller than the experimental width,
which disfavors $Y(4274)$ in a $D_{s}\bar{D}_{s0}$ molecular picture. If we
increase $\Lambda$ to higher values, the total widths of $Y(4274)$ cannot be
reproduced until a much larger $\Lambda$ value of about 9.6 is adopted.
Unfortunately, there are no such studies on taking $\Lambda=9.6$ or higher,
which is reasonable. Hence, the assignment as a $P-$wave $D_{s}\bar{D}_{s0}$
molecular state is impossible for $Y(4274)$ based on the total decay width
experimentally measured. This is quite different from the conclusion in Ref.
He:2016pfa  that the interaction between a $D_{s}$ meson and a $\bar{D}_{s0}$
meson is strong enough to form a bound state with a mass about 4274 MeV, which
can be associated to $Y(4274)$. Comparing our results with those in Ref.
He:2016pfa , it seems that a study of the spectroscopy alone does not give a
complete picture of its nature.

Combination our results in Fig. 10 with conclusion that in Ref. He:2016pfa ,
we predict a $P-$wave $D_{s}\bar{D}_{s0}$ molecular that we marked as
$Y^{{}^{\prime}}(4274)$ may exist. Taking a typical value
$\Gamma(D^{\pm}_{s0}\to\pi{}D^{\pm}_{s})=79.3$ KeV Faessler:2007gv , where
$D_{s0}$ is assumed to be a $DK$ bound state, the corresponding partial decay
widths are $\Gamma[Y^{{}^{\prime}}(4274)\to{}D_{s}^{*}\bar{D}_{s}]=0.20-0.33$
MeV, $\Gamma[Y^{{}^{\prime}}(4274)\to{}J/\psi\phi]=0.048-0.090$,
$\Gamma[Y^{{}^{\prime}}(4274)\to{}\pi^{0}{}D_{s}\bar{D}_{s}]=0.0066-0.0080$
MeV, $\Gamma[Y^{{}^{\prime}}(4274)\to{}D\bar{D}^{*}]=0.76-1.30$ MeV,
$\Gamma[Y^{{}^{\prime}}(4274)\to{}K\bar{K}^{*}]=0.065-0.14$ MeV,
$\Gamma[Y^{{}^{\prime}}(4274)\to{}\chi_{c0}\eta]=0.0$ MeV,
$\Gamma[Y^{{}^{\prime}}(4274)\to{}\chi_{c1}\eta]=0.0$ MeV, and
$\Gamma[Y^{{}^{\prime}}(4274)\to{}\phi\phi]=0.0089-0.018$ MeV, which yields a
total decay width of 1.06-1.84 MeV. And find that the transition
$Y^{{}^{\prime}}(4274)\to{}D\bar{D}^{*}$ is the main decay channel, almost
saturating the total width. The experimental measurements for this strong
decay process could be a crucial to observe such a state.

##  IV summary

In this work, inspired by the studies in Ref. He:2016pfa  that showed the
likely existence of a $D_{s}\bar{D}_{s0}$ bound state, we have studied its
partial decay widths into $Y(4274)\to{}J/\psi\phi$, $\chi_{c0}\eta$,
$\chi_{c1}\eta$, $\phi\phi$, $\eta\eta$, $D^{*}_{s}\bar{D}_{s}$,
$D\bar{D}^{*}$, $K\bar{K}^{*}$ and $\pi^{0}\bar{D}_{s}D_{s}$. These decays
involve the treatment of $Y(4274)$ state as a quasi-bound state of
$D_{s}\bar{D}_{s0}$ and utilizing the Weinberg compositeness condition to
determine the corresponding coupling. Our studies find the $P-$wave
$D_{s}\bar{D}_{s0}$ assignments for $Y(4274)$ is not allowed, it may be a
compact tetraquark state or diquark-antidiquark state Chen:2016oma ;
Wang:2016dcb ; Stancu:2009ka ; Agaev:2017foq ; Lu:2016cwr ; Gui:2018rvv ;
Maiani:2016wlq ; Zhu:2016arf . If $Y(4274)$ is a compact tetraquark state or
diquark-antidiquark state Chen:2016oma ; Wang:2016dcb ; Stancu:2009ka ;
Agaev:2017foq ; Lu:2016cwr ; Gui:2018rvv ; Maiani:2016wlq ; Zhu:2016arf , we
suggest experimental to search for $Y(4274)$ in the
$Y(4274)\to{}\chi_{c0}\eta$ and $Y(4274)\to{}\chi_{c1}\eta$ reactions that the
partial decay widths are of the order of 0.0 MeV by assuming $Y(4274)$ is a
$D_{s}\bar{D}_{s0}$ molecular state. Theoretical investigations on decay modes
and further experimental information on partial decay widths will be helpful
to distinguish which the inner structure of the $Y(4274)$ state is possible.

However, a $P-$wave $D_{s}\bar{D}_{s0}$ molecular with the total decay width
is at the order of 1.06 to 1.84 MeV [we marked as $Y^{{}^{\prime}}(4274)$] is
found. The predicted decay width seems to suggest that it is possible to
observe such a state at Belle or Belle II, e.g., via the inclusive invariant
mass distribution $D\bar{D}^{*}$, which is the largest transition. On the
other hand, its production yields at these experimental setups remain to be
studied.

If the $Y(4274)$ could be a $P-$wave $D_{s}\bar{D}_{s0}$ molecular. It
naturally leads us to think about whether there are contributions come from
the light pseudoscalar meson exchanges $(\pi,K,\eta,\eta^{\prime})$ and the
vector meson exchanges $(\rho,\omega,K^{*})$…, or there exist other decay
mechanisms for the few-body systems. However, only $J/\psi$, $\phi$, $\eta$,
$K$, $D$, $D_{s}$, and $D_{s}^{*}$ exchange are included because there is no
information on studies about other meson exchange in the vertices of charm-
strange meson. Moreover, we always think the two-body decay modes of the
multi-quark states are usually the dominant ones. Hence, the currently
calculation is enough to explain the $Y(4274)$ cannot be a $P-$wave
$D_{s}\bar{D}_{s0}$ molecular.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China
under Grant No.12104076, the Science and Technology Research Program of
Chongqing Municipal Education Commission (Grant No. KJQN201800510), and the
Opened Fund of the State Key Laboratory on Integrated Optoelectronics
(GrantNo. IOSKL2017KF19). Yin Huang want to thanks the support from the
Development and Exchange Platform for the Theoretic Physics of Southwest
Jiaotong University under Grants No.11947404 and No.12047576, the Fundamental
Research Funds for the Central Universities(Grant No. 2682020CX70), and the
National Natural Science Foundation of China under Grant No.12005177.

## References

  * (1) P. A. Zyla et al. [Particle Data Group], PTEP 2020, 083C01 (2020).
  * (2) R. Aaij et al. [LHCb], Phys. Rev. Lett. 122, 222001 (2019).
  * (3) R. Aaij et al. [LHCb], Sci. Bull. 66, 1278-1287 (2021).
  * (4) R. Aaij et al. [LHCb], Phys. Rev. Lett. 118, 022003 (2017).
  * (5) R. Aaij et al. [LHCb], Phys. Rev. D 95,012002 (2017).
  * (6) K. Yi [CDF], PoS ICHEP2010, 182 (2010) arXiv:1010.3470.
  * (7) H. X. Chen, E. L. Cui, W. Chen, X. Liu and S. L. Zhu, Eur. Phys. J. C 77, 160 (2017).
  * (8) Z. G. Wang, Eur. Phys. J. C 77, 174 (2017).
  * (9) F. Stancu, J. Phys. G 37, 075017 (2010) [erratum: J. Phys. G 46, 019501 (2019)].
  * (10) Q. F. Lü and Y. B. Dong, Phys. Rev. D 94,074007 (2016).
  * (11) L. C. Gui, L. S. Lu, Q. F. Lü, X. H. Zhong and Q. Zhao, Phys. Rev. D 98, 016010 (2018).
  * (12) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 95, 114003 (2017).
  * (13) L. Maiani, A. D. Polosa and V. Riquer, Phys. Rev. D 94, 054026 (2016).
  * (14) R. Zhu, Phys. Rev. D 94, 054009 (2016).
  * (15) J. He, Phys. Rev. D 95, 074004 (2017).
  * (16) Y. R. Liu, Phys. Rev. D 88, 074008 (2013) doi:10.1103/PhysRevD.88.074008 [arXiv:1304.7467 [hep-ph]].
  * (17) Y. L. Ma, Estimates for $X(4350)$ Decays from the Effective Lagrangian Approach, Phys. Rev. D 82, 015013 (2010).
  * (18) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 81, 074011 (2010).
  * (19) A. Faessler, T. Gutsche, V. E. Lyubovitskij and Y. L. Ma, Phys. Rev. D 76, 014005 (2007).
  * (20) A. Faessler, T. Gutsche, V. E. Lyubovitskij and Y. L. Ma, Phys. Rev. D 76, 114008 (2007).
  * (21) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 77, 094013 (2008).
  * (22) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, J. Phys. G 38, 015001 (2011).
  * (23) Y. Dong, A. Faessler, T. Gutsche, S. Kovalenko and V. E. Lyubovitskij, Phys. Rev. D 79, 094013 (2009).
  * (24) Y. Dong, A. Faessler, T. Gutsche, Q.F. Lü and V. E. Lyubovitskij, Phys. Rev. D 96, 074027 (2017).
  * (25) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 90, 094001 (2014).
  * (26) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 90, 074032 (2014).
  * (27) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 89, 034018 (2014).
  * (28) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 88, 014030 (2013).
  * (29) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Few Body Syst. 54, 1011 (2013).
  * (30) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, J. Phys. G 40, 015002 (2013).
  * (31) Y. Dong, A. Faessler, T. Gutsche, S. Kumano and V. E. Lyubovitskij, Phys. Rev. D 83, 094005 (2011).
  * (32) Y. Dong, A. Faessler, T. Gutsche, S. Kumano and V. E. Lyubovitskij, Phys. Rev. D 82, 034035 (2010).
  * (33) Y. Dong, A. Faessler, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 81, 014006 (2010).
  * (34) Y. Dong, A. Faessler and V. E. Lyubovitskij, Prog. Part. Nucl. Phys. 94, 282 (2017).
  * (35) F. Yang, Y. Huang and H. Q. Zhu, [arXiv:2107.13267 [hep-ph]].
  * (36) H. Zhu, N. Ma and Y. Huang, Eur. Phys. J. C 80, 1184 (2020).
  * (37) S. Weinberg, Phys. Rev. 130, 776 (1963).
  * (38) A. Salam, Nuovo Cim. 25, 224 (1962).
  * (39) E. Oset, J. R. Pelaez and L. Roca, Phys. Rev. D 67, 073013 (2003).
  * (40) A. Bramon, A. Grau and G. Pancheri, Phys. Lett. B 283, 416 (1992).
  * (41) P. Colangelo, F. De Fazio and T. N. Pham, Phys. Rev. D 69, 054023 (2004).
  * (42) R. Casalbuoni, A. Deandrea, N. Di Bartolomeo, R. Gatto, F. Feruglio and G. Nardulli, Phys. Rept. 281, 145-238 (1997).
  * (43) Y. Huang, M. Z. Liu, Y. W. Pan, L. S. Geng, A. Martínez Torres and K. P. Khemchandani, Phys. Rev. D 101,014022 (2020).
  * (44) D. Gamermann, E. Oset, D. Strottman and M. J. Vicente Vacas, Phys. Rev. D 76, 074016 (2007).
  * (45) B. Aubert et al. [BaBar], Phys. Rev. Lett. 90, 242001 (2003).
  * (46) X. Liu, Y. M. Yu, S. M. Zhao and X. Q. Li, Eur. Phys. J. C 47, 445-452 (2006).
  * (47) J. He and X. Liu, Eur. Phys. J. C 72, 1986 (2012).
  * (48) M. Ablikim et al. [BESIII], Phys. Rev. D 97, 051103 (2018).
  * (49) B. Aubert et al. [BaBar], Phys. Rev. D 74, 032007 (2006).