# Computational multiscale methods for parabolic wave approximations in
heterogeneous media

Eric Chung , Yalchin Efendiev , Sai-Mang Pun , and Zecheng Zhang Department of
Mathematics, The Chinese University of Hong Kong, Shatin, Hong KongDepartment
of Mathematics and Institute of Scientific Computing, Texas A&M University,
College Station, TX 77843, USADepartment of Mathematics, Texas A&M University,
College Station, TX 77843, USADepartment of Mathematics, Purdue University,
West Lafayette, IN 47906, USA

###### Abstract

In this paper, we develop a computational multiscale to solve the parabolic
wave approximation with heterogeneous and variable media. Parabolic wave
approximation is a technique to approximate the full wave equation. One
benefit of the method is that: one wave propagation direction can be taken as
an evolution direction, and we then can discretize it using a classical scheme
like Backward Euler. Consequently, we obtain a set of quasi-gas-dynamic (QGD)
models with different heterogeneous permeability fields. Then, we employ
constraint energy minimization generalized multiscale finite element method
(CEM-GMsFEM) to perform spatial discretization for the problem. The resulting
system can be solved by combining the central difference in time evolution.
Due to the variable media, we apply the technique of proper orthogonal
decomposition (POD) to further the dimension of the problem and solve the
corresponding model problem in the POD space instead of in the whole
multiscale space spanned by all possible multiscale basis functions. We prove
the stability of the full discretization scheme and give the convergence
analysis of the proposed approximation scheme. Numerical results verify the
effectiveness of the proposed method.

##  1 Introduction

Parabolic wave approximation have been used to approximate wave equations with
a preferred direction [2, 3]. To be more specific, we study the approximation
of the full wave equation:

$\displaystyle\rho\partial_{tt}u-\nabla\cdot(\mu\nabla u)=0.$ (1)

Here, $\nabla$ denotes the gradient operator in $\mathbb{R}^{3}$; $\rho$ and
$\mu$ are positive functions in $\mathbb{R}^{3}$. We then consider the
parabolic approximation of (1) in $D\times\Omega\subseteq\mathbb{R}^{3}$: for
$(z,x)\in D\times\Omega$ with the boundary $\Gamma:=\partial(D\times\Omega)$,
find $v=\sigma^{1/2}u$ such that

$\displaystyle\begin{split}c^{-1}\partial_{tt}v+\partial_{t}(\partial_{z}v)-\frac{1}{2}\nabla_{x}\cdot\left(c\nabla_{x}v\right)&=0&\quad\text{in
}D\times\Omega\times(0,T],\\\ v(z,x,0)&=v_{0}(z,x)&\quad\text{in
}D\times\Omega,\\\ \partial_{t}v(z,x,0)&=v_{1}(z,x)&\quad\text{in
}D\times\Omega,\\\ v(z,x,t)&=g(z,x,t)&\quad\text{in
}\Gamma\times(0,T].\end{split}$ (2)

Here, $\nabla_{x}$ denotes the gradient operator defined in the bounded domain
$\Omega\subset\mathbb{R}^{2}$; $D=[d_{1},d_{2}]\subset\mathbb{R}$ is a bounded
domain; $v_{0}$ and $v_{1}$ are initial conditions; $g$ is boundary condition
and $T>0$ is a given terminal time.

Wave equations of this type have a propagation direction $z$ which plays a
role of time [4]. Many real-world problems can be solved by the parabolic
approximation; in particular in the areas of geology [22], under water
acoustics [5, 6, 43, 44] and optics [23, 32, 36]. There are two benefits for
using the parabolic approximation instead of the full wave equation [2, 3]:
(i) it is easier to be realized and computationally more efficient; and (ii)
it makes it possible to use the approximation as an evolution equation in $z$
direction. The second property makes it possible to solve the problem as an
equation evolving in $z$ direction; we hence can discretize ($\partial_{z}v$)
by applying classical difference scheme. Consequently, we obtain a quasi-gas-
dynamic (QGD) equation which have been thoroughly studied in literature [11,
12, 13].

It is relatively common in the geological problems [22] that the medium under
consideration is highly variable and heterogeneous; that is, the function
field $c$ is non-homogeneous given a cross section in $z$ and is fast-changing
in $z$ direction. This brings in two difficulties. The first to mention is the
multiscale property brought by the heterogeneous field. The second difficulty
is the intense work in solving a QGD model given a $z$ cross-section in the
$z$ evolution. Throughout the work, we study (2) with variable and
heterogeneous media. We will derive scheme in discretizing the time and $z$
evolution and also provide our approach to solve the problems mentioned
before.

As we have discussed, the model has heterogeneous in $z$ direction. Directly
solving this problem on fine mesh can capture the multiscale features;
however, this is computationally intense and this issue becomes exacerbate
when people are solving time-dependent problems. Therefore, many methods which
solve the multiscale problems on coarser mesh have been proposed. These
include homogenization-based approaches [8, 9, 10, 29, 40, 41, 47], multiscale
finite element methods [31, 35, 37, 38], generalized multiscale finite element
methods (GMsFEM) [14, 15, 18, 19, 21, 25, 30], constraint energy minimizing
GMsFEM (CEM-GMsFEM) [16, 17], nonlocal multi-continua (NLMC) approaches [20],
metric-based upscaling [45], heterogeneous multiscale method [1, 24],
localized orthogonal decomposition (LOD) [33, 42], equation free approaches
[46, 48, 49], computational continua [26, 27, 28], hierarchical multiscale
method [7, 34, 50], and so on. Some of these approaches, such as
homogenization-based approaches, are designed for problems with scale
separation. In this work, we apply the CEM-GMsFEM [16, 17] and provide the
convergence analysis of our proposed scheme based on the coarse mesh
convergence results of the CEM-GMsFEM.

The second difficulty of the problem is the variable media. If we discretize
$z$ evolution using some classical difference scheme, each $z$ level is a QGD
model with heterogeneous. This model can be solved in the framework of CEM-
GMsFEM [11]; however the coarse scale basis evaluation is time consuming; in
particular, this process will be repeated for each level of $z$. We hence
proposed the proper orthogonal decomposition (POD) technique [39] which is
target to find low dimensional subspace such that the error of the orthogonal
projection is minimized in the sense of the norm induced by the inner product
of the original space. To be more specific, the proposed method can be
summarized as follows:

We proposed a method to solve a parabolic-wave model with variable and
heterogeneous media. We first apply the backward Euler scheme to
discretization the term $\partial_{z}v$. This leads to a set of two
dimensional (space) QGD models (3); we call this level of discretization the
quasi-time scheme and an unconditional stability result is established. To
further discretize the problem, we apply the central scheme to deal with the
time derivative $v_{tt}$, $v_{t}$; and then use CEM-GMsFEM in space on coarse
scale to capture the heterogeneous brought by the media. We then prove that
the full discretization scheme is stable in an energy norm under some CFL
condition by using an inverse inequality in the multiscale space.

The key of the CEM-GMsFEM method is to construct the CEM basis. The standard
procedure is first to build the auxiliary multiscale basis by solving local
spectral problems in coarse mesh; we then can construct the CEM basis by
evaluating a set of energy minimization problems. Due to the variable media,
we need to construct a set of basis for each QGD model and the solution of the
full discretized scheme is in the space of all multiscale basis. This is time
consuming; and we hence can apply the POD technique to find the best set of
orthogonal basis in the sense of $L_{2}$ minimization; that is, the projection
error of the full discretized solution onto the POD basis is optimal in the
norm induced by the original space. In practice, we collect CEM basis for some
QGD models and then construct POD basis of the space spanned by all CEM basis;
the POD models will finally be solved by using POD basis. A convergence
analysis of the POD approximation is established and the numerical results
prove the algorithm is successful.

The remainder of the paper is organized as follows. In Section 2, we present
some preliminaries of the model problem and briefly overview the framework of
proper orthogonal decomposition. Section 3 is devoted to the multiscale
methods and we will briefly overview the construction of multiscale basis
function within the framework of the CEM-GMsFEM. In Section 4, we present a
complete analysis of the proposed computational multiscale method. We then
present some numerical results to demonstrate the efficiency of the proposed
method in Section 5. Concluding remarks are drawn in Section 6.

##  2 Preliminaries

In this section, we present some preliminaries of the model problem. For
simplicity, we assume the homogeneous boundary condition $g=0$ is equipped in
the model problem (2). The extension of inhomogeneous case is straightforward.
The system (2) is the one we shall consider throughout the remainder of this
paper and we shall develop computational multiscale method for efficiently
simulating the problem. We remark that under appropriate regularity
assumptions on initial and boundary conditions, the problem (2) has a unique
solution such that

$t\to v(z,x,t)\in W^{1,\infty}(0,T;L^{2}(D\times\Omega))\cap
L^{\infty}(0,T;L^{2}(D\times\Omega)),$ $z\to v(z,x,t)\in
L^{\infty}(D,H^{1}(0,T;L^{2}(\Omega))).$

The last property enables us to consider $z$ as an evolution direction. The
result is an application of the semigroup theory and the Hille-Yoshida
theorem. See [3, Section 4] for more details.

Instead of the PDE formulation (2), we consider its corresponding variational
formulation. In the following, we treat $z$ as an evolution direction and use
backward Euler method to discretize the term $\partial_{z}v$. We divide the
domain $D=[d_{1},d_{2}]$ along $z$-direction into $K$ pieces. We write, for
$k=0,1,\cdots,K$,

$v_{k}=v(z_{k})\quad\text{with}\quad z_{k}=d_{1}+k\Delta z\quad\text{and}\quad
d_{2}=d_{1}+K\Delta z.$

The quasi-time discretization of (2) reads: find $\theta_{k}\in V$ for
$k=1,\cdots,K-1$ such that

$\displaystyle\left(\ddot{\theta}_{k},w;z_{k}\right)_{c}+\frac{1}{\Delta
z}\left(\dot{\theta}_{k}-\dot{\theta}_{k-1},w\right)+\frac{1}{2}a(\theta_{k},w;z_{k})=0\quad\text{
for all }w\in V:=H_{0}^{1}(\Omega)$ (3)

or equivalently

$\Delta
z\left(\dot{\theta}_{k},w;z_{k}\right)_{c}+\left(\dot{\theta}_{k},w\right)+\frac{\Delta
z}{2}a(\theta_{k},w;z_{k})=\left(\dot{\theta}_{k-1},w\right)\quad\text{ for
all }w\in V,$

where $(v,w;z_{k})_{c}:=\int_{\Omega}c^{-1}(z_{k})vw\leavevmode\nobreak\ dx$,
$(v,w):=\int_{\Omega}vw\leavevmode\nobreak\ dx$, and
$a(v,w;z_{k}):=\int_{\Omega}c(z_{k})\nabla v\cdot\nabla w\leavevmode\nobreak\
dx$. Here, we denote $\dot{v}:=\partial_{t}v$ and $\ddot{v}:=\partial_{tt}v$.
We assume that $c\in L^{\infty}(D\times\Omega)$ and we denote
$c_{\max}:=\left\|c\right\|_{L^{\infty}(D\times\Omega)}$. Employing Galerkin’s
method and the method of energy estimate, one can show the well-posedness of
the variational formulation (3). Denote $\left\|v\right\|:=\sqrt{(v,v)}$,
$\left\|v\right\|_{c(z_{k})}:=\sqrt{(v,v;z_{k})_{c}}$, and
$\left\|v\right\|_{a(z_{k})}:=\sqrt{a(v,v;z_{k})}$. We establish the following
stability estimate for the quasi-time discretization (3).

######  Lemma 2.1.

Let $\\{v_{k}\\}_{k=0}^{K}\subseteq V$ solve the equation (3). Then, the
following stability estimate holds:

$\left\|\dot{v_{K}}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta
z\sum_{k=1}^{K}\mathcal{E}_{k}(v_{k}(T))\lesssim\left\|\dot{v_{0}}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta
z\sum_{k=1}^{K}\mathcal{E}_{k}(v_{k}(0)),$

where we denote
$\mathcal{E}_{k}(v):=\left\|\dot{v}\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v\right\|_{a(z_{k})}^{2}$.

######  Remark.

The term
$\mathcal{E}_{k}(v_{k}(0))=\left\|v_{1}(z_{k})\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{0}(z_{k})\right\|_{a(z_{k})}^{2}$
on the right-hand side of the inequality above is determined by the initial
conditions.

###### Proof of Lemma 2.1.

Taking $w=\dot{v_{k}}$, we have

$\displaystyle\begin{split}&\quad\frac{\Delta
z}{2}\frac{d}{dt}\left(\left\|\dot{v}_{k}\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}\right\|_{a(z_{k})}^{2}\right)+\left\|\dot{v}_{k}\right\|^{2}-(\dot{v}_{k-1},\dot{v}_{k})\\\
&=\Delta
z\left(\ddot{v}_{k},\dot{v}_{k};z_{k}\right)_{c}+(\dot{v}_{k}-\dot{v}_{k-1},\dot{v_{k}})+\frac{\Delta
z}{2}a(v_{k},\dot{v_{k}};z_{k})=0.\end{split}$

Thus, we have

$\quad\frac{\Delta
z}{2}\frac{d}{dt}\left(\left\|\dot{v}_{k}\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}\right\|_{a(z_{k})}^{2}\right)+\left\|\dot{v}_{k}\right\|^{2}=(\dot{v}_{k},\dot{v}_{k-1})\leq\left\|\dot{v}_{k}\right\|\left\|\dot{v}_{k-1}\right\|.$

Multiplying by $2$ and integrating over $(0,T]$, we have

$\displaystyle\begin{split}&{\Delta
z}\left(\left\|\dot{v}_{k}(T)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(T)\right\|_{a(z_{k})}^{2}\right)+2\int_{0}^{T}\left\|\dot{v}_{k}\right\|^{2}dt\\\
&\leq
2\int_{0}^{T}\left\|\dot{v}_{k}\right\|\left\|\dot{v}_{k-1}\right\|dt+\Delta
z\left(\left\|\dot{v}_{k}(0)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(0)\right\|_{a(z_{k})}^{2}\right)\\\
&\leq\int_{0}^{T}\left\|\dot{v}_{k}\right\|^{2}dt+\int_{0}^{T}\left\|\dot{v}_{k-1}\right\|^{2}dt+\Delta
z\left(\left\|\dot{v}_{k}(0)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(0)\right\|_{a(z_{k})}^{2}\right).\\\
\end{split}$

Therefore, we have

$\displaystyle\begin{split}\left\|\dot{v}_{k}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}&+{\Delta
z}\left(\left\|\dot{v}_{k}(T)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(T)\right\|_{a(z_{k})}^{2}\right)\\\
&\leq\left\|\dot{v}_{k-1}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta
z\left(\left\|\dot{v}_{k}(0)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(0)\right\|_{a(z_{k})}^{2}\right).\end{split}$

Summing over $k=1,\cdots,K$, we obtain the desired result and this completes
the proof. ∎

###  2.1 The Proper Orthogonal Decomposition

In this section, we briefly introduce the proper orthogonal decomposition
(POD) method. This method aims to generate optimally ordered orthogonal basis
functions in the least squares sense for a given set of theoretical,
experimental, or computational data. Reduced-order models or surrogate models
are then obtained by truncating this set of optimal basis functions, providing
considerable computational savings over the original high-dimensional
problems.

Let $X$ be a real Hilbert space endowed with inner product $(\cdot,\cdot)_{X}$
and norm $\left\|\cdot\right\|_{X}$. We set
$\mathcal{V}:=\operatorname{span}\\{y_{1},y_{2},\cdots,y_{n}\\}$ with each
$y_{i}\in X$ for $i\in\\{1,\cdots,n\\}$. We refer to $\mathcal{V}$ as ensemble
consisting of the snapshots $\\{y_{i}\\}_{i=1}^{n}$, at least one of which is
assumed to be non-zero. Let $\\{\psi_{k}\\}_{k=1}^{\mathcal{N}}$ denote a set
of orthonormal basis functions of $\mathcal{V}$ with
$\mathcal{N}:=\text{dim}(\mathcal{V})\leq n$. Then, each member of the
ensemble can be expressed as

$y_{j}=\sum_{k=1}^{\mathcal{N}}(y_{j},\psi_{k})_{X}\psi_{k}$

for each $j\in\\{1,\cdots,n\\}$. The POD method consists in choosing the
orthonormal basis functions such that for every
$\ell\in\\{1,\cdots,\mathcal{N}\\}$ the mean square error between the elements
$y_{j}$ (for any $j\in\\{1,\cdots,n\\}$), and the corresponding $\ell$-th
partial sum is minimized on average:

$\displaystyle\begin{split}&\min_{\\{\psi_{k}\\}_{k=1}^{\ell}}\frac{1}{n}\sum_{j=1}^{n}\left\|y_{j}-\sum_{k=1}^{\ell}(y_{j},\psi_{k})_{X}\psi_{k}\right\|_{X}^{2}\\\
&\text{subject to}\quad(\psi_{k},\psi_{t})=\delta_{kt}\quad\text{for any
}k,t\in\\{1,2,\cdots,\ell\\}.\end{split}$ (4)

Here, $\delta_{kt}$ denotes the Kronecker-delta function. A solution
$\\{\psi_{k}\\}_{k=1}^{\ell}$ to (4) is called a POD-basis of rank $\ell$. We
introduce the correlation matrix

$K=\left(\frac{1}{n}(y_{j},y_{i})_{X}\right)\in\mathbb{R}^{n\times n}$

corresponding to the snapshots $\\{y_{j}\\}_{j=1}^{n}$. The matrix $K$ is
positive semi-definite and has rank $\mathcal{N}$. The minimization problem
(4) can be reduced to an eigenvalue problem

$\displaystyle Kv=\lambda v.$ (5)

We sort all the positive eigenvalues in a decreasing order as
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{\mathcal{N}}>0$ and the
associated eigenvectors are denoted by $v_{k}$ with $k=1,\cdots,\mathcal{N}$.
It can be shown that the POD-basis of rank $\ell\in\mathbb{N}^{+}$ with
$\ell\leq\mathcal{N}$ is formed by

$\displaystyle\varphi_{k}=\frac{1}{\sqrt{\lambda_{k}}}\sum_{j=1}^{n}(v_{k})_{j}y_{j}\quad\text{for
}k=1,\cdots,\ell.$ (6)

Here, $(v_{k})_{j}$ is the $j$-th component of the eigenvector $v_{k}$. The
basis functions $\\{\varphi_{k}\\}_{k=1}^{\ell}$ form a POD-basis of rank
$\ell$ and we have the following error formula.

######  Proposition 2.2.

Let $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{\mathcal{N}}>0$ be the
positive eigenvalues of $K$ in (5) and
$v_{1},\cdots,v_{\mathcal{N}}\in\mathbb{R}^{n}$ be the associated
eigenvectors. Then, $\\{\varphi_{k}\\}_{k=1}^{\ell}$ given by (6) forms a set
of POD-basis of rank $\ell$ with $\ell\leq\mathcal{N}$. Moreover, we have the
error formula

$\frac{1}{n}\sum_{j=1}^{n}\left\|y_{j}-\sum_{k=1}^{\ell}(y_{j},\varphi_{k})_{X}\varphi_{k}\right\|_{X}^{2}=\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.$

In practice, we shall make use of the decay property of eigenvalues in
$\lambda_{k}$ and choose the first $\ell$ dominant eigenvalues such that the
ratio
$\zeta:=\frac{\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}}{\sum_{k=1}^{\mathcal{N}}\lambda_{k}}$
is small enough to achieve an expected accuracy, for instance $\zeta=1\%$. One
would prefer the eigenvalues decays as fast as possible so that one can ensure
high accuracy with few POD basis functions.

##  3 Multiscale Method

In this section, we develop the computational multiscale method in order to
solve the parabolic wave approximation. For spatial discretization, we will
apply the CEM-GMsFEM. In particular, for each node $z_{i}\in D$ along the
$z$-direction, we will construct a set of multiscale basis functions in the
spirit of CEM-GMsFEM. To further reduce the dimension of the multiscale space,
we will perform POD procedure related to these CEM basis functions. Once the
multiscale space is constructed, one can use leapfrog scheme to discretize
time derivatives and solve the resulting fully-discretized problem.

###  3.1 Spatial Discretization: CEM-GMsFEM

In this section, we briefly outline the framework of CEM-GMsFEM and present
the construction of the multiscale space. First, we introduce fine and coarse
grids for the computational domain. Let
$\mathcal{T}^{H}=\\{K_{i}\\}_{i=1}^{N}$ be a conforming partition of the
domain $\Omega$ with mesh size $H>0$ defined by

$H:=\max_{K\in\mathcal{T}^{H}}\Big{(}\max_{x,y\in K}\lvert x-y\rvert\Big{)}.$

We refer to this partition as coarse grid. We denote $N\in\mathbb{N}^{+}$ the
total number of coarse elements. Subordinate to the coarse grid, we define the
fine grid partition $\mathcal{T}^{h}$ (with mesh size $h\ll H$) by refining
each coarse element $K\in\mathcal{T}^{H}$ into a connected union of finer
elements. We assume the refinement above is performed such that
$\mathcal{T}^{h}$ is also a conforming partition of the domain $\Omega$.
Denote $N_{c}$ the number of interior coarse grid nodes of $\mathcal{T}^{H}$
and $\\{x_{i}\\}_{i=1}^{N_{c}}$ the collection of interior coarse nodes in the
coarse grid.

####  3.1.1 Spectral Decomposition

We present the construction of the auxiliary multiscale basis functions. Let
$K_{i}\in\mathcal{T}^{H}$ be a coarse block. Define $V(K_{i})$ as the
restriction of the abstract space $V$ on the coarse element $K_{i}$. For each
$z_{k}\in D$, we consider a local spectral problem: Find
$\sigma_{j}^{(i,k)}\in\mathbb{R}$ and $\phi_{j}^{(i,k)}\in V(K_{i})$ such that

$\displaystyle
a_{i}(\phi_{j}^{(i,k)},v;z_{k})=\sigma_{j}^{(i,k)}s_{i}(\phi_{j}^{(i,k)},v;z_{k})\quad\text{for
all}\leavevmode\nobreak\ v\in V(K_{i}).$ (7)

Here, $a_{i}:V(K_{i})\times V(K_{i})$ is a symmetric non-negative definite
bilinear form and $s_{i}:V(K_{i})\times V(K_{i})$ is a symmetric positive
definite bilinear form. We remark that the above problem is solved on the fine
mesh in the actual computations. Based on the analysis, we choose

$a_{i}(v,w;z_{k}):=\int_{K_{i}}c(z_{k})\nabla v\cdot\nabla
w\leavevmode\nobreak\ dx,\quad
s_{i}(v,w;z_{k}):=\int_{K_{i}}\tilde{\kappa}(z_{k})vw\leavevmode\nobreak\ dx,$

where
$\tilde{\kappa}:=\sum_{j=1}^{N_{c}}c(z_{k})\lvert\nabla\chi_{j,k}^{\text{ms}}\rvert^{2}$.
The functions $\\{\chi_{j,k}^{\text{ms}}\\}_{j=1}^{N_{c}}$ are the standard
multiscale finite element basis functions which satisfy the partition of unity
property. More precisely, $\chi_{j,k}^{\text{ms}}$ is the solution of the
following system:

$\displaystyle\nabla\cdot(c(z_{k})\nabla\chi_{j,k}^{\text{ms}})$
$\displaystyle=0\quad$ $\displaystyle\text{in each }K\subset\omega_{j},$
$\displaystyle\chi_{j,k}^{\text{ms}}$ $\displaystyle=g_{j}\quad$
$\displaystyle\text{on }\partial K\setminus\partial\omega_{j},$
$\displaystyle\chi_{j,k}^{\text{ms}}$ $\displaystyle=0\quad$
$\displaystyle\text{on }\partial\omega_{j}.$

The function $g_{j}$ is continuous and linear along the boundary of the coarse
element. We assume that the eigenvalues $\sigma_{j}^{(i,k)}$ are arranged in
ascending order and we pick $\ell_{i,k}\in\mathbb{N}^{+}$ corresponding
eigenfunctions to construct the local auxiliary space
$V_{\text{aux}}^{(i,k)}:=\text{span}\\{\phi_{j}^{(i,k)}:j=1,\cdots,\ell_{i}\\}$.
We assume the normalization
$s_{i}\left(\phi_{j}^{(i,k)},\phi_{j}^{(i,k)};z_{k}\right)=1$. After that, we
define the global auxiliary multiscale space
$V_{\text{aux}}^{k}:=\bigoplus_{i=1}^{N}V_{\text{aux}}^{(i,k)}$. We remark
that the global auxiliary multiscale space is used to construct multiscale
basis functions that are orthogonal to the auxiliary space with respect to the
weighted $L^{2}$ inner product $s(\cdot,\cdot;z_{k})$.

Note that the bilinear form $s_{i}(\cdot,\cdot;z_{k})$ defines an inner
product with norm
$\left\|\cdot\right\|_{s(K_{i};k)}:=\sqrt{s_{i}(\cdot,\cdot;z_{k})}$ in the
local auxiliary space $V_{\text{aux}}^{(i,k)}$. Based on these local inner
products and norms, one can naturally define a new inner product and norm for
the global auxiliary space $V_{\text{aux}}$ as follows: for all $v,w\in
V_{\text{aux}}^{k}$,

$\displaystyle
s(v,w;z_{k}):=\sum_{i=1}^{N}s_{i}(v,w;z_{k})\quad\text{and}\quad\left\|v\right\|_{s(z_{k})}:=\sqrt{s(v,v;z_{k})}.$
(8)

Note that if $\\{\chi_{j,k}^{\text{ms}}\\}_{j=1}^{N_{c}}$ is a set of bilinear
partition of unity, then $\left\|v\right\|_{s(z_{k})}\leq
H^{-1}(\max\\{c\\})^{1/2}\left\|v\right\|$ for any $v\in L^{2}(\Omega)$. In
addition, we define $\pi_{k}:L^{2}(\Omega)\to V_{\text{aux}}^{k}$ as the
projection with respect to the inner product $s(\cdot,\cdot;z_{k})$ such that

$\pi_{k}(u):=\sum_{i=1}^{N}\sum_{j=1}^{\ell_{i}}s_{i}(u,\phi_{j}^{(i,k)};z_{k})\phi_{j}^{(i,k)}\quad\text{for
all }u\in L^{2}(\Omega).$

####  3.1.2 The construction of multiscale basis functions

In this section, we present the construction of the multiscale basis
functions. First, we define an oversampling region for each coarse element.
Specifically, given a non-negative integer $m\in\mathbb{N}$ and a coarse
element $K_{i}$, we define the oversampling region $K_{i,m}\subset\Omega$ such
that

$K_{i,m}:=\left\\{\begin{array}[]{lr}K_{i}&\text{if }m=0,\\\
\displaystyle{\bigcup\\{K:K_{i,m-1}\cap K\neq\emptyset\\}}&\text{if }m\geq
1.\end{array}\right.$

See Figure 1 for an illustration of oversampling region. For simplicity, we
denote $K_{i}^{+}$ the oversampled region $K_{i,m}$ for some nonnegative
integer $m$.

Figure 1: Oversampling region with $m=1$.

Recall that $V(K_{i}^{+})$ is the restriction of $V$ on the coarse patch
$K_{i}^{+}$. Let $V_{0}(K_{i}^{+})$ be the subspace of $V(K_{i}^{+})$ with
zero trace on the boundary $\partial K_{i}^{+}$. For each eigenfunction
$\phi_{j}^{(i,k)}\in V_{\text{aux}}^{k}$, we define the multiscale basis
$\psi_{j,\text{ms}}^{(i,k)}\in V_{0}(K_{i}^{+})$ to be the solution of the
equation:

$\displaystyle
a(\psi_{j,\text{ms}}^{(i,k)},v)+s\left(\pi(\psi_{j,\text{ms}}^{(i,k)}),\pi_{k}(v)\right)=s(\phi_{j}^{(i,k)},v)\quad\text{for
all }v\in V_{0}(K_{i}^{+}).$ (9)

Then, the multiscale space is defined as
$V_{\text{ms}}^{k}:=\text{span}\left\\{\psi_{j,\text{ms}}^{(i,k)}:i=1,\cdots,N,\leavevmode\nobreak\
j=1,\cdots,\ell_{i,k}\right\\}$. By construction, we have
$\text{dim}(V_{\text{ms}}^{k})=\text{dim}(V_{\text{aux}}^{k})$. After that, we
define

$V_{\text{CEM}}:=\text{span}\left\\{\psi_{j,\text{ms}}^{(i,k)}:i=1,\cdots,N,\leavevmode\nobreak\
j=1,\cdots,\ell_{i,k},\leavevmode\nobreak\ k=0,\cdots,K\right\\}.$

We will perform the procedure of POD for this space to construct the
multiscale space that will be used in actual simulation.

###  3.2 Construction of multiscale reduced basis functions using POD

In this section, we present the construction of multiscale reduced basis
functions using the POD technique. First, we define
$\\{v_{\text{CEM},k}\\}_{k=0}^{K}\subset V_{\text{CEM}}$ to be the solution of
the following equation:

$\displaystyle\left(\ddot{v}_{\text{CEM},k},w;z_{k}\right)_{c}+\frac{1}{\Delta
z}\left(\dot{v}_{\text{CEM},k}-\dot{v}_{\text{CEM},k-1},w\right)+\frac{1}{2}a(v_{\text{CEM},k},w;z_{k})=0\quad\text{
for all }w\in V_{\text{CEM}}.$ (10)

In the abstract framework of POD, we set the space
$X=L^{2}(0,T;H^{1}(\Omega))$. We define the inner product associated with $X$
as $\left\langle u,v\right\rangle:=\int_{0}^{T}(\nabla u,\nabla
v)+(u,v)\leavevmode\nobreak\ dt$ for $u,v\in X$. The corresponding norm will
then be defined as
$\left\|\cdot\right\|_{L^{2}(0,T;H^{1}(\Omega))}:=\sqrt{\left\langle\cdot,\cdot\right\rangle}$.
Next, we define the snapshot space $\mathcal{V}$ as the collection of:

$y_{j}:=v_{\text{CEM}}(t,z_{j-1}),\quad
y_{j+K+1}:=\ddot{v}_{\text{CEM}}(t,z_{j-1}),\quad j=1,\cdots,K+1,$

and

$y_{j+2K+2}=\tilde{\partial}\dot{v}_{\text{CEM}}(t,z_{j}):=\frac{\dot{v}_{\text{CEM}}(t,z_{j})-\dot{v}_{\text{CEM}}(t,z_{j-1})}{\Delta
z},\quad j=1,\cdots,K.$

Then, we perform the POD procedure on the snapshot space $\mathcal{V}$ as
described in Section 2.1. It should be noted that the correlation matrix is
defined as $K_{ij}:=\frac{1}{3K+2}(y_{j},y_{i})_{H^{1}}$. We denote
$\\{\psi_{k}\\}_{k=1}^{\mathcal{N}}$ the corresponding POD basis functions
with cardinality $\mathcal{N}\in\mathbb{N}^{+}$. We then define
$V_{\text{POD}}^{\ell}:=\operatorname{span}\\{\psi_{k}\\}_{k=1}^{\ell}$ for a
given positive number $\ell\leq\mathcal{N}$.

###  3.3 Fully Discretization

In this section, we present the fully discretization for the problem (3). We
can further consider the discretization in time. We divide $(0,T]$ into $N$
pieces and write

$v^{n}=v(t^{n})\quad\text{with}\quad t^{n}=n\Delta t\quad\text{and}\quad\Delta
t=\frac{T}{N}.$

Specifically, we use first- and second-order central difference schemes. We
define $\tau:=\Delta t\Delta z$. The fully discretization reads: for
$k=1,\cdots,K-1$ and $n=1,\cdots,N-1$, find $\\{v_{k}^{n}\\}\subset
V_{\text{POD}}^{\ell}$ such that

$\displaystyle\left(\frac{v_{k}^{n+1}-2v_{k}^{n}+v_{k}^{n-1}}{(\Delta
t)^{2}},w;z_{k}\right)_{c}+\left(\frac{v_{k}^{n+1}-v_{k}^{n-1}}{2\tau},w\right)+\frac{1}{2}a(v_{k}^{n},w;z_{k})=\left(\frac{v_{k-1}^{n+1}-v_{k-1}^{n-1}}{2\tau},w\right)$
(11)

for all $w\in V_{\text{POD}}^{\ell}$. In short, we discretize implicitly in
$z$ and explicitly in time. We show that the fully discretization is stable in
time. To this aim, we first recall the inverse inequality for the multiscale
space, which is proved in [11].

######  Proposition 3.1 (Lemma 4.5 in [11]).

Assume that $\\{\chi_{j,k}^{\text{ms}}\\}_{j=1}^{N_{c}}$ is a set of bilinear
partition of unity. For any $v\in V_{\text{CEM}}$, there is a constant
$C_{\text{inv}}>0$ such that

$\displaystyle\left\|\nabla v\right\|_{a}\leq
C_{\text{inv}}H^{-1}\left\|v\right\|.$

Then, we have the following stability result for (11).

######  Lemma 3.2.

Suppose that the following CFL condition

$\frac{1}{c_{\max}}-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta
t)^{2}\geq\delta,$

holds for some $\delta>0$. Then, we have the following stability estimate:

$\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{K}^{n+1}-v_{K}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{N,k}\leq\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{0}^{n+1}-v_{0}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{0,k}.$

###### Proof.

Define the energy as follows:

$\displaystyle\mathcal{E}_{n,k}:=\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta
t}\right\|_{c}^{2}+\frac{1}{2}a(v_{k}^{n+1},v_{k}^{n}).$

Using the inverse inequality for the multiscale functions, one can show that

$\displaystyle\mathcal{E}_{n,k}$
$\displaystyle=\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta
t}\right\|_{c}^{2}+\frac{1}{2}a(v_{k}^{n+1},v_{k}^{n})$
$\displaystyle=\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta
t}\right\|_{c}^{2}+\frac{1}{4}a(v_{k}^{n+1},v_{k}^{n+1})+\frac{1}{4}a(v_{k}^{n},v_{k}^{n})-\frac{1}{4}a(v_{k}^{n+1}-v_{k}^{n},v_{k}^{n+1}-v_{k}^{n})$
$\displaystyle\geq\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta
t}\right\|_{c}^{2}+\frac{1}{4}\left(\left\|v_{k}^{n+1}\right\|_{a}^{2}+\left\|v_{k}^{n}\right\|_{a}^{2}\right)-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta
t)^{2}\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\|^{2}$
$\displaystyle\geq\big{(}\frac{1}{c_{\max}}-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta
t)^{2}\big{)}\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta
t}\right\|^{2}+\frac{1}{4}\left(\left\|v_{k}^{n+1}\right\|_{a}^{2}+\left\|v_{k}^{n}\right\|_{a}^{2}\right)\geq
0.$

Next, taking $w=v_{k}^{n+1}-v_{k}^{n-1}$ in (11) and denoting
$f^{n}_{k-1}=v_{k-1}^{n+1}-v_{k-1}^{n-1}$, we have

$\displaystyle\frac{1}{2\tau}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta
t}\right\|_{c}^{2}-\left\|\frac{v_{k}^{n}-v_{k}^{n-1}}{\Delta
t}\right\|_{c}^{2}+\frac{1}{2}a(v_{k}^{n},v_{k}^{n+1}-v_{k}^{n-1})=\frac{1}{2\tau}\left(f^{n}_{k-1},v_{k}^{n+1}-v_{k}^{n-1}\right).$

Then, it implies that

$\displaystyle\begin{split}\frac{1}{2\tau}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\mathcal{E}_{n,k}-\mathcal{E}_{n-1,k}&=\frac{1}{2\tau}\left(f^{n}_{k},v_{k}^{n+1}-v_{k}^{n-1}\right)\\\
&\leq\frac{1}{2\tau}\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|\\\
&\leq\frac{1}{4\tau}\left(\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|^{2}+\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}\right)\\\
\implies\frac{1}{4\tau}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\mathcal{E}_{n,k}&\leq\frac{1}{4\tau}\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|^{2}+\mathcal{E}_{n-1,k}.\end{split}$

Summing over $n=1,\cdots,N$, we have

$\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\mathcal{E}_{N,k}\leq\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|^{2}+\mathcal{E}_{0,k}.$

Therefore, summing over $k=1,\cdots,K$, we have

$\frac{1}{4\tau}\sum_{n=1}^{N}\left(\sum_{k=1}^{K}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}-\sum_{k=0}^{K-1}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}\right)+\sum_{k=1}^{K}\mathcal{E}_{N,k}\leq\sum_{k=1}^{K}\mathcal{E}_{0,k}$

and

$\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{K}^{n+1}-v_{K}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{N,k}\leq\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{0}^{n+1}-v_{0}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{0,k}.$

This completes the proof. ∎

######  Remark.

It should be noted that the CFL condition depends on the quantity
$c_{\max}:=\left\|c\right\|_{L^{\infty}(D\times\Omega)}$. It also relies on
the constant $C_{\text{inv}}$ which comes from the inverse inequality (see
Lemma 4.5 in [11]). From the proof of that lemma, the constant
$C_{\text{inv}}$ is independent of the ratio of the contrast values in the
permeability. We are using the CEM-GMsFEM, which is a coarse mesh method,
hence the coarse mesh size $H$ will make the CFL much less restrictive. The
idea to improve the CFL condition is to use the implicit scheme in time; more
precisely, one can use $a(v_{k}^{n+1},w;z_{k})$ in (11). This will become one
of our future works.

##  4 Convergence Analysis

In this section, we present the convergence analysis of the proposed POD-based
multiscale method for the parabolic wave equation. Let $v_{\text{CEM}}\in
V_{\text{CEM}}$ be the semi-discretized solution which solves the problem
(10). Throughout this section, we assume that the coercivity and boundedness
of this bilinear form hold. Specifically, there exist two constants $\gamma$
and $\beta$, independent of any $z_{k}\in D$, such that

$\displaystyle\gamma\|u\|_{H_{1}(\Omega)}^{2}\leq
a(u,u;z_{k})\quad\text{and}\quad
a(u,w;z_{k})\leq\beta\|u\|_{H_{1}(\Omega)}\|w\|_{H_{1}(\Omega)}$ (12)

for any $u,w\in V$. We denote $(\cdot,\cdot)_{H^{1}}$ the $H^{1}$-inner
product such that $(u,v)_{H^{1}}:=(\nabla u,\nabla v)+(u,v)$ for any $u,v\in
V$ with its associated norm being
$\left\|\cdot\right\|_{H^{1}}:=\sqrt{(\cdot,\cdot)_{H^{1}}}$.

We have the following lemma for the proper orthogonal decomposition.

######  Lemma 4.1 (cf. Proposition 1 in [39]).

Let $v_{\text{CEM}}\in V_{\text{CEM}}$ be the semi-discretized solution which
solves the problem (10) and $K_{ij}:=\frac{1}{3K+2}(y_{j},y_{i})_{H^{1}}$ be
the correlation matrix and $\lambda_{k}$ be the corresponding eigenvalues
sorted ascending. For any integer $\ell$ with $0<\ell\leq\mathcal{N}$, the
following error formula holds:

$\displaystyle\begin{split}\frac{1}{3K+2}&\left[\int_{0}^{T}\sum_{k=1}^{K+1}\left\|v_{\text{CEM}}(t,z_{k})-\sum_{j=1}^{\ell}\left(v_{\text{CEM}}(t,z_{k}),\psi_{j}\right)_{H^{1}}\psi_{j}\right\|_{H^{1}}^{2}\right.\\\
&+\left.\sum_{k=1}^{K}\left\|\tilde{\partial}\dot{v}_{\text{CEM}}(t,z_{k})-\sum_{j=1}^{\ell}\left(\tilde{\partial}\dot{v}_{\text{CEM}}(t,z_{k}),\psi_{j}\right)_{H^{1}}\psi_{j}\right\|_{H^{1}}^{2}\right.\\\
&\left.+\sum_{k=1}^{K+1}\left\|\ddot{v}_{\text{CEM}}(t,z_{k})-\sum_{j=1}^{\ell}\left(\ddot{v}_{\text{CEM}}(t,z_{k}),\psi_{j}\right)_{H^{1}}\psi_{j}\right\|_{H^{1}}^{2}dt\right]=\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.\end{split}$
(13)

Next, we define the Ritz-projection $P^{\ell}:V\rightarrow
V_{\text{POD}}^{\ell}$ by:

$\displaystyle(P^{\ell}v,\psi)_{H^{1}}=(v,\psi)_{H^{1}}\quad\text{ for all
}\psi\in V_{\text{POD}}^{\ell},$ (14)

for any $v\in V$. We have the following estimate for the projection operator.

######  Lemma 4.2.

Let $\\{z_{k}\\}_{k=1}^{K+1}\subset D$. For any
$\ell\in\\{1,2,...,\mathcal{N}\\}$ and $v\in V$, the projection operator
$P^{\ell}$ satisfies the following estimate:

$\displaystyle\frac{1}{K+1}\sum_{k=1}^{K+1}\int_{0}^{T}\left\|v_{\text{CEM}}(z_{k})-P^{\ell}v_{\text{CEM}}(z_{k})\right\|_{H^{1}}^{2}\leavevmode\nobreak\
dt\lesssim\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.$ (15)

###### Proof.

Let $v\in V$ be arbitrary. By the definition of the projection $P^{\ell}$, we
have

$\displaystyle\begin{split}\left\|v-P^{\ell}v\right\|_{H^{1}}^{2}=(v-P^{\ell}v,v-P^{\ell}v)_{H^{1}}=(v-P^{l}v,v-\psi)_{H^{1}}\leq\left\|v-P^{\ell}v\right\|_{H^{1}}\left\|v-\psi\right\|_{H^{1}}\end{split}$

for all $\psi\in V_{\text{POD}}^{\ell}$. This implies that,

$\displaystyle\left\|v-P^{\ell}v\right\|_{H^{1}}\lesssim\left\|v-\psi\right\|_{H^{1}}\quad\text{
for all }\psi\in V_{\text{POD}}^{\ell}.$ (16)

The lemma follows immediately by Lemma 4.1 and (16). ∎

Clearly, from the proof of Lemma 4.2, the corollary below follows immediately.

######  Corollary 4.3.

Let $v_{\text{CEM}}\in V_{\text{CEM}}$ be the semi-discretized solution which
solves (10).

Then, the following estimates hold:

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\int_{0}^{T}\left\|\tilde{\partial}\dot{v}_{\text{CEM}}(z_{k})-P^{\ell}\tilde{\partial}\dot{v}_{\text{CEM}}(z_{k})\right\|_{H^{1}}^{2}\leavevmode\nobreak\
dt\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}$ (17)

and

$\displaystyle\frac{1}{K+1}\sum_{k=0}^{K}\int_{0}^{T}\left\|\ddot{v}_{\text{CEM}}(z_{k})-P^{\ell}\ddot{v}_{\text{CEM}}(z_{k})\right\|_{H^{1}}^{2}dt\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.$
(18)

Let $\\{U_{k}\\}_{k=0}^{K}\subset V_{\text{POD}}^{\ell}$ be the semi-
discretized solution satisfying the following equation

$\displaystyle(\ddot{U}_{k},\phi;z_{k})_{c}+(\tilde{\partial}\dot{U}_{k},\phi)+\frac{1}{2}a(U_{k},\phi;z_{k})=0$
(19)

for all $\phi\in V_{\text{POD}}^{\ell}$ with appropriate initial condition
$U_{0}$. Here, we denote
$\tilde{\partial}\dot{U}_{k}:=\frac{\dot{U}_{k}-\dot{U}_{k-1}}{\Delta z}$. We
are now able to analyze the error. Assume that $\\{v_{k}^{n}\\}\subset
V_{\text{POD}}^{\ell}$ solves (11), $\\{U_{k}\\}_{k=0}^{K}\subset
V_{\text{POD}}^{\ell}$ solves (19), and the function $v_{\text{CEM}}\in
V_{\text{CEM}}$ solves (10). We decompose the error into three parts:

$v_{k}^{n}-v_{\text{CEM}}(t_{n})=\underbrace{v_{k}^{n}-U_{k}(t_{n})}_{=:\mu_{k}^{n}}+\underbrace{U_{k}(t_{n})-P^{\ell}v_{\text{CEM}}(t_{n})}_{=:\nu_{k}^{n}}+\underbrace{P^{\ell}v_{\text{CEM}}(t_{n})-v_{\text{CEM}}(t_{n})}_{=:\rho_{k}^{n}}.$

We denote $\mu_{k}(t)$, $\nu_{k}(t)$, and $\rho_{k}(t)$ the piecewise linear
functions that interpolates $\\{\mu_{k}^{n}\\}$, $\\{\nu_{k}^{n}\\}$, and
$\\{\rho_{k}^{n}\\}$ in time, respectively. Due to Lemma 4.2, we have

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|\rho_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.$

It should be noted that the term $\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}$
comes from the method of POD and is a typical error term in the POD analysis.
The error decay to $0$ very fast; and by the theory, this error is optimal
[39] since the method of POD solves a minimization problem (4).

Next, using the notation
$\tilde{\partial}\nu_{k}=\frac{\nu_{k}-\nu_{k-1}}{\Delta z}$ for all
$k=1,...,K$ and the equation (10), we obtain

$\displaystyle(\ddot{\nu}_{k},\psi;z_{k})_{c}+(\tilde{\partial}\dot{\nu}_{k},\psi)+\frac{1}{2}a(\nu_{k},\psi;z_{k})$
$\displaystyle=-(P^{\ell}\ddot{v}_{\text{CEM}},\psi;z_{k})_{c}-(\tilde{\partial}P^{\ell}\dot{v}_{\text{CEM}},\psi)-\frac{1}{2}a(v_{\text{CEM}},\psi;z_{k})$
$\displaystyle=(\ddot{v}_{\text{CEM}}-P^{\ell}\ddot{v}_{\text{CEM}},\psi;z_{k})_{c}+(\tilde{\partial}\dot{v}_{\text{CEM}}-\tilde{\partial}P^{\ell}\dot{v}_{\text{CEM}},\psi)$

for any $\psi\in V_{\text{POD}}^{\ell}$. Let us denote

$h_{k}:=\tilde{\partial}\dot{v}_{\text{CEM}}-\tilde{\partial}P^{\ell}\dot{v}_{\text{CEM}}\quad\text{and}\quad
w_{k}:=\ddot{v}_{\text{CEM}}-P^{\ell}\ddot{v}_{\text{CEM}}.$

By Corollary 4.3, it follows that

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|h_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}+\frac{1}{K+1}\sum_{k=0}^{K}\left\|w_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.$

Using the same technique of showing Lemma 2.1, one can obtain an estimate for
$\nu_{k}$:

$\displaystyle\|\dot{\nu}_{k}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta
z\sum_{k=1}^{K}\mathcal{E}_{k}(\nu_{k}(T))\lesssim\Delta
z\sum_{k=1}^{K}\left\|h_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}+\left\|w_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.$

It remains to estimate the term $\mu_{k}$. This error comes from the temporal
discretization and it follows that

$\left\|\mu_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}=\int_{0}^{T}\left\|\mu_{k}\right\|_{H^{1}}^{2}dt=\int_{0}^{T}O(\tau^{2})dt=O(\tau).$

Denote $v_{k}(t)$ the piecewise linear function that interpolates
$\\{v_{k}^{n}\\}$ in time. As a result, we have the following error estimate

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|v_{k}-v_{\text{CEM}}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim
O(\tau)+\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.$ (20)

We remark that the error between the solution $\theta$ (that solves (3)) and
the solution $v_{\text{CEM}}$ has the relation: $\left\|\theta-
v_{\text{CEM}}\right\|_{L^{2}(0,T;H^{1}(\Omega))}=O(H)$. Therefore, the whole
error estimate reads as follows:

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|\theta-
v_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim
O(H^{2})+O(\tau)+\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.$ (21)

##  5 Numerical Experiments

In this section, we present some numerical results to demonstrate the
efficiency of the proposed computational multiscale method. In the experiments
below, we set the spatial domain to be $\Omega=(0,1)^{2}$. We partition the
spatial domain into $10\times 10$ uniform square elements and we refer this
partition to the coarse grid with mesh size $H=\sqrt{2}/10$. Further, we
divide each coarse element into $10\times 10$ uniform square elements and the
corresponding fine grid has resolution of the size $100\times 100$. We refer
this to the fine grid of the spatial domain with mesh size $h=\sqrt{2}/100$.
We compute the numerical solutions $v_{k}^{N}$ at the terminal time for all
$k=1,\cdots,K$ and compare it with the reference solutions $v_{f}(z_{k},x,T)$
for all $k$, which are computed by using the underlying fine grid. We measure
the relative error in $L^{2}$ norm at each layer $z_{k}$, which is defined as
follows:

$e_{2}^{k}:=\frac{\left\|v_{k}^{N}-v_{f}(z_{k},x,T)\right\|}{\left\|v_{f}(z_{k},x,T)\right\|}.$

###  5.1 The First Experiment

In the first experiment, we set the initial condition to be $v(0,x,t)=\sin(\pi
x_{1})\sin(\pi x_{2})\sin(t)$ where $x=(x_{1},x_{2})\in\Omega$; also we set
homogeneous initial condition for all $z$, i.e., $v(z,x,0)=0$ for all
$(z,x)\in D\times\Omega$. The temporal direction and $z$-direction are
discretized as in the scheme (11). We choose the time step to be $\Delta
t=10^{-5}$ and the step size in $z$-direction is $\Delta z=10^{-4}$. We set
$K=30$, i.e., the $z$ direction is partitioned into $30$ levels and
$D=[0,30\Delta z]$.

In this experiment, the pattern of the permeability field $c(x)$ is given in
Figure 2 (Left). In $x_{1}$\- and $x_{2}$-directions, the permeability will
follow the same pattern. However, along the $z$-direction, the permeability
field $c(x)$ has different values of contrast. See Figure 2 (Right) for an
illustration for this layer-structured heterogeneous media. In particular, we
partition the permeability field along the $z$-direction into $3$ different
layers; the value of contrast in the yellow region within each layer are equal
to $10$, $15$, and $20$, respectively. We remark that the patterns keep the
same as demonstrated in Figure 2 (Left).

To obtain the set of POD basis functions, we first construct the CEM basis
functions for each layer and then obtain the reduced basis functions by
performing the POD procedure on the combined basis space. We choose 100 POD
basis functions from the total $300$ CEM basis functions in this experiment.
We remark that the number of the POD reduced basis functions depends on the
decay of the eigenvalues $\\{\lambda_{k}\\}$ corresponding to the POD
construction. In practice, one has to include sufficiently many reduced basis
functions so that the tail of the eigenvalues is smaller than a given
tolerance of accuracy.

Figure 2: Left: The two-dimensional pattern of permeability field in
Experiment $1$. Right: Permeability field at different layers $z_{k}$ has the
same pattern but different values of contrast.

The record of relative error $e_{2}^{k}$ for all $k\in\\{1,\cdots,K\\}$ at the
terminal time is plotted in Figure 3 (Left). The relative error at each layer
is around the level of magnitude of $10^{-3}$. One can observe from the graph
that the proposed numerical scheme is capable for accurately approximating the
fine-scale solution.

Figure 3: The relative error $e_{2}^{k}$. Left: Experiment 1; Right:
Experiment 2.

###  5.2 Experiment 2

In the second experiment, we consider the so-called Marmousi heterogeneous
field. More precisely, given $z_{k}\in D$, $c(z_{k})$ is a Marmousi
permeability field as shown in Figure 4.

Figure 4: Marmousi permeability field in Experiment $2$.

The initial and boundary data are the same as in the first experiment. We set
$T=5\times 10^{-3}$ and $\Delta t=5\times 10^{-6}$. We set $\Delta z=10^{-4}$
and $D=[0,10\Delta z]$ with $K=10$. For each $z_{k}\in D$, the medium
$c(z_{k})$ is obtained by taking the data from a very fine Marmousi field of
size $600\times 600$. To be more specific, we partition the fine field (which
is of size $600\times 600$) into a coarse field of size $100\times 100$ whose
coarse element size is $6\times 6$. Then, in each coarse element, we randomly
pick a value in one of the fine elements in this coarse element; and we
formulate our computational fine grid ($100\times 100$) as combining all
selected elements accordingly.

To get the set of POD basis functions, we first calculate three sets of CEM
basis functions with $k\in\\{0,4,8\\}$ and obtain the POD basis by performing
the POD procedure on the combined set of CEM basis functions. We choose $50$
POD basis functions from the total $300$ CEM basis functions. The record of
relative error for all $z_{k}$ at the terminal time is plotted in Figure 3
(Right). The relative error at each layer is around the level of magnitude of
$10^{-4}$ in this case. In particular, the relative error $e_{2}^{k}$ is about
$3.9\times 10^{-4}$ for $k=10$. This demonstrates the efficiency of the
proposed multiscale method for the simulation of parabolic wave formulation.

##  6 Conclusion

In this paper, we developed a computational multiscale method for simulating
the parabolic wave formulation with highly heterogeneous media. For the
spatial discretization, we employed the recently developed CEM-GMsFEM, which
is proved to be efficient to reduce the dimension of the model in space. We
then combined the technique of proper orthogonal decomposition to further
reduce the dimension along the quasi time direction. A complete analysis of
the proposed algorithm has been provided. Numerical results are provided to
demonstrate the effectiveness and efficiency of the proposed multiscale
method.

## Acknowledgement

## References

  * [1] A. Abdulle and B. Engquist.  Finite element heterogeneous multiscale methods with near optimal computational complexity.  Multiscale Modeling & Simulation, 6(4):1059–1084, 2007.
  * [2] A. Bamberger, B. Engquist, L. Halpern, and P. Joly.  Higher order paraxial wave equation approximations in heterogeneous media.  SIAM Journal on Applied Mathematics, 48(1):129–154, 1988.
  * [3] A. Bamberger, B. Engquist, L. Halpern, and P. Joly.  Parabolic wave equation approximations in heterogenous media.  SIAM Journal on Applied Mathematics, 48(1):99–128, 1988.
  * [4] A. Bamberger, B. Euggnist, L. Helpern, and P. Joly.  The paraxial approximation for the wave equation: some new results.  Adv. Comput. Methods Partial Differential Equations, pages 340–344, 1984.
  * [5] P. Blomgren, G. Papanicolaou, and H. Zhao.  Super-resolution in time-reversal acoustics.  The Journal of the Acoustical Society of America, 111(1):230–248, 2002.
  * [6] H. Brock, R. Buchal, and C. Spofford.  Modifying the sound-speed profile to improve the accuracy of the parabolic-equation technique.  The Journal of the Acoustical Society of America, 62(3):543–552, 1977.
  * [7] D. L. Brown, Y. Efendiev, and V. H. Hoang.  An efficient hierarchical multiscale finite element method for stokes equations in slowly varying media.  Multiscale Modeling & Simulation, 11(1):30–58, 2013.
  * [8] E. Cances, V. Ehrlacher, F. Legoll, and B. Stamm.  An embedded corrector problem to approximate the homogenized coefficients of an elliptic equation.  Comptes Rendus Mathematique, 353(9):801–806, 2015.
  * [9] J. Chen, S. Sun, and Z. He.  Homogenize coupled stokes–cahn–hilliard system to darcy’s law for two-phase fluid flow in porous medium by volume averaging.  Journal of Porous Media, 22(1), 2019.
  * [10] J. Chen, S. Sun, and X. Wang.  Homogenization of two-phase fluid flow in porous media via volume averaging.  Journal of Computational and Applied Mathematics, 353:265–282, 2019\.
  * [11] B. Chetverushkin, E. Chung, Y. Efendiev, S.-M. Pun, and Z. Zhang.  Computational multiscale methods for quasi-gas dynamic equations.  arXiv preprint arXiv:2009.00068, 2020.
  * [12] B. Chetverushkin, A. Saveliev, and V. Saveliev.  Compact quasi-gasdynamic system for high-performance computations.  Computational Mathematics and Mathematical Physics, 59(3):493–500, 2019.
  * [13] B. N. Chetverushkin, N. D’Ascenzo, A. V. Saveliev, and V. Saveliev.  Kinetic model and magnetogasdynamics equations.  Computational Mathematics and Mathematical Physics, 58(5):691–699, 2018.
  * [14] E. Chung, Y. Efendiev, and T. Y. Hou.  Adaptive multiscale model reduction with generalized multiscale finite element methods.  Journal of Computational Physics, 320:69–95, 2016.
  * [15] E. Chung, Y. Efendiev, and C. S. Lee.  Mixed generalized multiscale finite element methods and applications.  Multiscale Modeling & Simulation, 13:338–366, 2014.
  * [16] E. Chung, Y. Efendiev, and W. T. Leung.  Constraint energy minimizing generalized multiscale finite element method.  Computer Methods in Applied Mechanics and Engineering, 339:298–319, 2018.
  * [17] E. Chung, Y. Efendiev, and W. T. Leung.  Constraint energy minimizing generalized multiscale finite element method in the mixed formulation.  Computational Geosciences, 22(3):677–693, 2018.
  * [18] E. T. Chung, Y. Efendiev, and W. T. Leung.  Generalized multiscale finite element methods for wave propagation in heterogeneous media.  Multiscale Modeling & Simulation, 12(4):1691–1721, 2014.
  * [19] E. T. Chung, Y. Efendiev, and W. T. Leung.  Fast online generalized multiscale finite element method using constraint energy minimization.  Journal of Computational Physics, 355:450–463, 2018.
  * [20] E. T. Chung, Y. Efendiev, W. T. Leung, M. Vasilyeva, and Y. Wang.  Non-local multi-continua upscaling for flows in heterogeneous fractured media.  Journal of Computational Physics, 372:22–34, 2018.
  * [21] E. T. Chung, W. T. Leung, and S. Pollock.  Goal-oriented adaptivity for GMsFEM.  Journal of Computational and Applied Mathematics, pages 625–637, 2015.
  * [22] J. F. Claerbout.  Fundamentals of geophysical data processing, volume 274.  Citeseer, 1976.
  * [23] J. D. Cole.  Modern developments in transonic flow.  SIAM Journal on Applied Mathematics, 29(4):763–787, 1975.
  * [24] W. E and B. Engquist.  The Heterogeneous Multiscale Methods.  Communications in Mathematical Sciences, 1(1):87–132, 2003.
  * [25] Y. Efendiev, J. Galvis, and T. Y. Hou.  Generalized multiscale finite element methods (GMsFEM).  Journal of Computational Physics, 251:116–135, 2013.
  * [26] D. Fafalis and J. Fish.  Computational continua for linear elastic heterogeneous solids on unstructured finite element meshes.  International Journal for Numerical Methods in Engineering, 115(4):501–530, 2018.
  * [27] J. Fish and S. Kuznetsov.  Computational continua.  International Journal for Numerical Methods in Engineering, 84(7):774–802, 2010.
  * [28] J. Fish and Z. Yuan.  Multiscale enrichment based on partition of unity.  International Journal for Numerical Methods in Engineering, 62(10):1341–1359, 2005.
  * [29] S. Fu, E. Chung, and G. Li.  Edge multiscale methods for elliptic problems with heterogeneous coefficients.  Journal of Computational Physics, 2019.
  * [30] K. Gao, S. Fu, R. L. Gibson Jr, E. T. Chung, and Y. Efendiev.  Generalized multiscale finite-element method (GMsFEM) for elastic wave propagation in heterogeneous, anisotropic media.  Journal of Computational Physics, 295:161–188, 2015.
  * [31] H. Hajibeygi, D. Kavounis, and P. Jenny.  A hierarchical fracture model for the iterative multiscale finite volume method.  Journal of Computational Physics, 230(4):8729–8743, 2011.
  * [32] A. Hasegawa and F. Tappert.  Transmission of stationary nonlinear optical pulses in dispersive dielectric fibers. i. anomalous dispersion.  Applied Physics Letters, 23(3):142–144, 1973.
  * [33] P. Henning, A. Målqvist, and D. Peterseim.  A localized orthogonal decomposition method for semi-linear elliptic problems.  ESAIM: Mathematical Modelling and Numerical Analysis, 48(5):1331–1349, 2014.
  * [34] V. Hoang and C. Schwab.  High dimensional finite elements for elliptic problems with multiple scales.  Multiscale Modeling & Simulation, 3:168–194, 2004.
  * [35] T. Y. Hou and X.-H. Wu.  A multiscale finite element method for elliptic problems in composite materials and porous media.  Journal of Computational Physics, 134:169–189, 1997.
  * [36] J. Hudson.  A parabolic approximation for elastic waves.  Wave Motion, 2(3):207–214, 1980.
  * [37] P. Jenny, S. Lee, and H. Tchelepi.  Multi-scale finite volume method for elliptic problems in subsurface flow simulation.  Journal of Computational Physics, 187:47–67, 2003.
  * [38] P. Jenny, S. Lee, and H. Tchelepi.  Adaptive multi-scale finite volume method for multi-phase flow and transport in porous media.  Multiscale Modeling & Simulation, 3:30–64, 2004.
  * [39] K. Kunisch and S. Volkwein.  Galerkin proper orthogonal decomposition methods for parabolic problems.  Numerische mathematik, 90(1):117–148, 2001.
  * [40] C. Le Bris, F. Legoll, and A. Lozinski.  An MsFEM type approach for perforated domains.  Multiscale Modeling & Simulation, 12(3):1046–1077, 2014.
  * [41] C. Le Bris, F. Legoll, and F. Thomines.  Multiscale finite element approach for weakly random problems and related issues.  ESAIM: Mathematical Modelling and Numerical Analysis, 48(3):815–858, 2014.
  * [42] A. Målqvist and D. Peterseim.  Localization of elliptic multiscale problems.  Mathematics of Computation, 83(290):2583–2603, 2014.
  * [43] S. T. McDaniel.  Parabolic approximations for underwater sound propagation.  The Journal of the Acoustical Society of America, 58(6):1178–1185, 1975.
  * [44] S. T. McDaniel.  Propagation of normal mode in the parabolic approximation.  The Journal of the Acoustical Society of America, 57(2):307–311, 1975.
  * [45] H. Owhadi and L. Zhang.  Metric-based upscaling.  Communications on Pure and Applied Mathematics, 60:675–723, 2007\.
  * [46] A. Roberts and I. Kevrekidis.  General tooth boundary conditions for equation free modeling.  SIAM Journal on Scientific Computing, 29(4):1495–1510, 2007.
  * [47] A. Salama, S. Sun, M. F. El Amin, Y. Wang, and K. Kumar.  Flow and Transport in Porous Media: A Multiscale Focus.  Geofluids, 2017, 2017.
  * [48] G. Samaey, I. Kevrekidis, and D. Roose.  Patch dynamics with buffers for homogenization problems.  Journal of Computational Physics, 213(1):264–287, 2006.
  * [49] G. Samaey, D. Roose, and I. Kevrekidis.  The gap-tooth scheme for homogenization problems.  Multiscale Modeling & Simulation, 4(1):278–306, 2005.
  * [50] W. C. Tan and V. H. Hoang.  High dimensional finite element method for multiscale nonlinear monotone parabolic equations.  Journal of Computational and Applied Mathematics, 345:471–500, 2019\.


# Dark states with electromagnetic form factors at electron colliders

Yu Zhang  School of Physics, Hefei University of Technology, Hefei
230601,China  Mao Song  School of Physics and Optoelectronics Engineering,
Anhui University, Hefei 230601,China  Liangwen Chen  Institute of Modern
Physics, CAS, Lanzhou 730000, China  Advanced Energy Science and Technology
Guangdong Laboratory, Huizhou 516000, China

###### Abstract

Electromagnetically neutral dark sector particles may feebly interact with
photons through higher dimensional effective operators, such as mass-dimension
5 magnetic and electric dipole moment, and a mass- dimension 6 anapole moment
and charge radius operators. In this work, we use hypercharge gauge field form
factors to treat dark states, which will induce not only electromagnetic form
factors but also the corresponding $Z$ boson operators. Taking a Dirac fermion
$\chi$ as an example, we investigate the probes of searching for such dark
states at current and future $e^{+}e^{-}$ collider experiments including
BESIII, STCF, Belle II and CEPC via monophoton searches. Comparing to current
experiments, we find that electron colliders including BESIII, STCF, Belle II,
which operate with the center-of-mass at several GeV, have leading sensitivity
on the corresponding electromagnetic form factors for the mass-dimension 5
operators with dark states lighter than several GeV, while they can not
provide competitive upper limits for the mass-dimension 6 operators. Future
CEPC operated with the center-of-mass on and beyond the mass of $Z$-boson with
competitive luminosity can probe the unexplored parameter space for dark
states with mass-dimension 5 (6) operators in the mass region of $m\lesssim
100$ GeV (10 MeV $\lesssim m\lesssim 100$ GeV).

##  I Introduction

Primary importance for our understanding of elementary interactions is
shedding light on the dark sector states. Searches for kinetically mixed dark
photons and for new particles as mediators connecting to standard model (SM)
constituent of the prime dark-sector physics cases Essig:2013lka . However,
some dark states sharing a coupling to the SM photon have received
comparatively less attention. In this scenario, even if dark states are
perfectly electromagnetic (EM) neutral, higher-dimensional effective couplings
to the SM photons are still possible. For dark states $\chi$ considered as
Dirac fermion, through various moments, such as magnetic dipole moments (MDM)
and electric dipole moments (EDM) at mass-dimension 5, and anapole moment (AM)
and charge radius interaction (CR) at mass-dimension 6, the couplings to the
photon can be present Chu:2018qrm ; Kavanagh:2018xeh .

In general, by interactions associated with photons, dark states $\chi$ with
EM form factors can be produced, which therefore could be studied by
accelerator-based experiments and stars. In Ref. Chu:2018qrm , $\chi$ pair-
production in electron beams on fixed targets at NA64 NA64:2017vtt , LDMX
LDMX:2018cma , BDX BDX:2016akw , and mQ Prinz:1998ua  has been studied. In
addition, constraints from rare meson ($B$ and $K$) decays Bird:2004ts ;
E949:2004uaj ; BNL-E949:2009dza  and SM precision observables such as
$(g-2)_{\mu}$ and the running of the fine-structure constant are also worked
out. The current limits on and detection prospects of such dark sector
particles $\chi$ are illustrated in Ref. Chu:2020ysb , utilizing high-
intensity proton beams, concretely from LSND LSND:1996jxj , MiniBooNE-DM
MiniBooNEDM:2018cxm , CHARM-II CHARM-II:1989nic ; CHARM-II:1994dzw , and E613
Ball:1980ojt  and projected SHiP SHiP:2015vad  and DUNE DUNE:2015lol ;
DUNE:2017pqt . Via their scattering off electrons in the Forward Liquid Argon
Experiment detector at the LHC Forward Physics Facility, the prospects of
searching for electromagnetically interacting $\chi$ particles have been
studied in Ref. Kling:2022ykt . A detailed astrophysical study of stellar
cooling constraints for the mass of dark states $\chi$ dropping below MeV are
further complemented in Ref. Chu:2019rok ; Chang:2019xva . Using missing
energy searches at colliders, the constraints from $e^{+}e^{-}$ colliders
BaBar Chu:2018qrm  and LEP Fortin:2011hv ; Chu:2018qrm , and from proton-
proton collisions at LHC are also investigated Fortin:2011hv ; Arina:2020mxo .

In this work, we extend the analysis to other electron colliders operated at
the GeV scale, including BESIII Asner:2008nq , Belle II Belle-II:2010dht , and
the proposed Super Tau Charm Factory (STCF) Luo:2019xqt ; Charm-
TauFactory:2013cnj ; Shi:2020nrf , which can probe the light dark states with
the mass less than GeV scale. In order to accurately measure the discovered
Higgs, many high energy electron colliders are proposed, including Circular
Electron Positron Collider (CEPC) CEPCStudyGroup:2018ghi , Future Circular
Collider $e^{+}e^{-}$ (FCC-ee) FCC:2018evy , International Linear Collider
(ILC) ILC:2019gyn , and Compact Linear $e^{+}e^{-}$ Collider (CLIC)
Robson:2018enq . We will take CEPC as an example to investigate the
sensitivity on the dark state with EM form factor via higher-dimensional
moments at future high energy electron colliders.

This paper is organized as follows. In Sec. II, we introduce the interactions
between the dark states and photons. In Sec. III we describe the signal and
backgrounds of probing the dark state with electromagnetic form factor at
electron colliders. In Sec. IV we present the constraints on the corresponding
couplings at BESIII, STCF and Belle II that are operated with the center-of-
mass (CM) energy $\sqrt{s}\ll M_{Z}$. In Sec. V, we consider the constraints
at CEPC that is operated with $\sqrt{s}\geq M_{Z}$. Our conclusions are drawn
in Sec. VI.

##  II Dark states with Electromagnetic form factor interactions

In this work, we consider that the dark state $\chi$ is taken as a Dirac
fermion, which may have interactions with the photon gauge filed $A_{\mu}$ or
its field strength tensor
$F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ as Chu:2018qrm ;
Kavanagh:2018xeh ,

$\displaystyle{\cal
L}_{\chi}=\frac{1}{2}\mu_{\chi}\bar{\chi}\sigma^{\mu\nu}\chi
F_{\mu\nu}+\frac{i}{2}d_{\chi}\bar{\chi}\sigma^{\mu\nu}\gamma^{5}\chi
F_{\mu\nu}-a_{\chi}\bar{\chi}\gamma^{\mu}\gamma^{5}\chi\partial^{\nu}F_{\mu\nu}+b_{\chi}\bar{\chi}\gamma^{\mu}\chi\partial^{\nu}F_{\mu\nu}.$
(1)

Here $e$ is the electric charge; $\mu_{\chi}$ and $d_{\chi}$ are the
dimensional coefficients of the mass-dimension 5 MDM and EDM interactions,
expressed in units of the Bohr magneton $\mu_{B}\equiv e/(2m_{e})$ with
$m_{e}$ being the electron mass, and $\sigma_{\mu\nu}\equiv
i[\gamma^{\mu},\gamma^{\nu}]/2$; $a_{\chi}$ and $b_{\chi}$ are the dimensional
coefficients of the mass-dimension 6 AM and CR interactions.

Below the electroweak scale, the effective couplings between $\chi$ and
photons in Eq. (1) can be induced from more general hypercharge form factor
operators, i.e., with $F_{\mu\nu}$ replaced by the hypercharge field strength
$B_{\mu\nu}\equiv\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}$ and $A_{\mu}$
replaced by the associated gauge boson $B_{\mu}$ Shi:2020nrf ,

$\displaystyle{\cal
L}_{\chi}=\frac{1}{2}\mu_{\chi}^{B}\bar{\chi}\sigma^{\mu\nu}\chi
B_{\mu\nu}+\frac{i}{2}d_{\chi}^{B}\bar{\chi}\sigma^{\mu\nu}\gamma^{5}\chi
B_{\mu\nu}-a_{\chi}^{B}\bar{\chi}\gamma^{\mu}\gamma^{5}\chi\partial^{\nu}B_{\mu\nu}+b_{\chi}^{B}\bar{\chi}\gamma^{\mu}\chi\partial^{\nu}B_{\mu\nu}.$
(2)

Hypercharge form factors are linear combinations of electromagnetic form
factors and the corresponding $Z$-boson operators, weighted by appropriate
factors of the cosine and sine of the Weinberg angle, $c_{W}$ and $s_{W}$.
After electroweak symmetry breaking, it then can be observed that
${\mathcal{C}}_{\chi}={\mathcal{C}}_{\chi}^{B}c_{W}$ with
${\mathcal{C}}_{\chi}=\mu_{\chi},\,d_{\chi},\,a_{\chi},\,b_{\chi}$, and the
coupling between $\chi$ and the $Z$-boson will also be induced. All the
couplings are taken as real.

##  III Signal and Background at $e^{+}e^{-}$ colliders

The cross section for the single photon production from $e^{+}e^{-}$
annihilation, $e^{+}e^{-}\to\chi\bar{\chi}\gamma$, can be approximately
factorized into the process without photon emission,
$e^{+}e^{-}\to\gamma/Z\to\chi\bar{\chi}$, times the improved Altarelli-Parisi
radiator function Nicrosini:1989pn ; Montagna:1995wp ,

$\frac{d^{2}\sigma}{dx_{\gamma}dz_{\gamma}}=H\left(x_{\gamma},z_{\gamma};s\right)\sigma_{0}(s_{\gamma})\,,$
(3)

where the radiator function is

$H\left(x_{\gamma},z_{\gamma};s\right)=\frac{\alpha}{\pi}\frac{1}{x_{\gamma}}\left[\frac{1+(1-x_{\gamma})^{2}}{1-z_{\gamma}^{2}}-\frac{x_{\gamma}^{2}}{2}\right]\,.$
(4)

Here, ${s}$ and ${s_{\gamma}}$ are the square of the center-of-mass (CM)
energies of the $e^{+}e^{-}$ and $\chi\bar{\chi}$ system, respectively, with
$s_{\gamma}=(1-x_{\gamma})s$, $E_{\gamma}$ is the energy of the initial state
radiation (ISR) photon, $x_{\gamma}=2E_{\gamma}/\sqrt{s}$ is the energy
fraction emitted away by ISR, $z_{\gamma}=\cos\theta_{\gamma}$ with
$\theta_{\gamma}$ being the polar angle of the photon. The cross section of
the $\chi$ pair production without ISR $\sigma_{0}$ reads

$\displaystyle\sigma_{0}(s)$ $\displaystyle=$
$\displaystyle\frac{\alpha}{4}\frac{f\left(s\right)}{s^{2}}\sqrt{\frac{s-4m_{\chi}^{2}}{s}}\Bigg{[}c_{W}^{2}+(g_{L}+g_{R})\frac{s(s-M_{Z}^{2})}{(s-M_{Z}^{2})^{2}+M_{Z}^{2}\Gamma_{Z}^{2}}$
(5) $\displaystyle+$
$\displaystyle\frac{1}{2}\frac{1}{c_{W}^{2}}(g_{L}^{2}+g_{R}^{2})\frac{s^{2}}{(s-M_{Z}^{2})^{2}+M_{Z}^{2}\Gamma_{Z}^{2}}\Bigg{]},$

with $g_{L}=-\frac{1}{2}+s_{W}^{2}$ and $g_{R}=s_{W}^{2}$, and $M_{Z}$ and
$\Gamma_{Z}$ are the mass and decay width of the $Z$-boson. The factor $f(s)$
is given as

$\displaystyle\mathrm{MDM}:f\left(s\right)=\frac{2}{3}\mu_{\chi}^{2}s^{2}\left(1+\frac{8m_{\chi}^{2}}{s}\right),$
(6)
$\displaystyle\mathrm{EDM}:f\left(s\right)=\frac{2}{3}d_{\chi}^{2}s^{2}\left(1-\frac{4m_{\chi}^{2}}{s}\right),$
(7)
$\displaystyle\mathrm{AM}:f\left(s\right)=\frac{4}{3}a_{\chi}^{2}s^{3}\left(1-\frac{4m_{\chi}^{2}}{s}\right),$
(8)
$\displaystyle\mathrm{CR}:f\left(s\right)=\frac{4}{3}b_{\chi}^{2}s^{3}\left(1+\frac{2m_{\chi}^{2}}{s}\right).$
(9)

It can be found that, in the limit of $\frac{m_{\chi}^{2}}{s}\to 0$, the
production rates of the $\chi$ pair with the EM form factors for mass-
dimension-5 MDM and EDM operators have the same forms. So do for mass-
dimension-6 AM and CR operators. When $\sqrt{s}\ll M_{Z}$, the production rate
for $\chi$ pair in Eq. (5) tends to be same with the one only considering dark
sector-photon interactions Chu:2018qrm .

For monophoton searches at electron colliders, the backgrounds consist of two
categories: the irreducible background and the reducible background. The
irreducible background arises from the neutrino pair production associated
with one visible photon $e^{+}e^{-}\to\nu\bar{\nu}\gamma$. The reducible
background comes from one visible photon in the final state together with
several other SM particles that cannot be detected because of the detector
limitations. The reducible background will be discussed later in details for
each experiment, since it strongly depends on the detector performance.

##  IV $e^{+}e^{-}$ colliders operated with $\sqrt{s}\ll M_{Z}$

###  IV.1 Belle II

The constraints on the light dark states with electromagnetic form factors
from Belle II via monophoton searches have been investigated in Ref.
Chu:2018qrm , in which the authors follow Ref. Essig:2013vha , and scale up
the BaBar background from Ref. BaBar:2008aby  to an integrated luminosity of
50 ab-1 with the CM energy of 10.57 GeV, employ a constant efficiency cut of
50% in both search regions and take identical geometric cuts for Belle II and
BABAR. Ref. Kou:2018nap  has provided the exact background subtraction in the
monophoton search of Belle II to probe an invisibly decaying dark photon. In
this work, we revisit the constraints from Belle II via monophoton searches
following the strategy in Ref. Kou:2018nap .

At Belle II, the Electromagnetic Calorimeter (ECL), which covers a polar angle
region of $(12.4-155.1)^{\circ}$ and has inefficient gaps between the barrel
and the endcaps for polar angles between $(31.4-32.2)^{\circ}$ and
$(128.7-130.7)^{\circ}$ in the lab frame Kou:2018nap . The photons and
electrons can be detected in the ECL. In the monophoton signature, the
reducible background at Belle II consists of two major parts Liang:2019zkb :
one is mainly due to the lack of polar angle coverage of the ECL near the beam
direction, i.e., $\theta>155.1^{\circ}$ or $\theta<12.4^{\circ}$, which is
referred to as the “bBG”; the other one is mainly owing to the gaps between
the three segments in the ECL detector, i.e., $\theta\in(31.4-32.2)^{\circ}$
or $(128.7-130.7)^{\circ}$, which is referred to as the “gBG”.

The bBG comes from the electromagnetic processes
$e^{+}e^{-}\to\gamma+\not{X}$, dominated by
$e^{+}e^{-}\to\gamma\not{\gamma}(\not{\gamma})$ and
$e^{+}e^{-}\to\gamma\not{e}^{+}\not{e}^{-}$. Here, $\not{X}$ denotes the other
particle(s) in the final state are emitted along the beam directions. Thus,
except the single detected photon all the other final state particles are
emitted along the beam directions with $\theta>155.1^{\circ}$ or
$\theta<12.4^{\circ}$ in the lab frame, which are out of the cover polar angle
region of the ECL.

At the asymmetric Belle II detector, for the monophoton events from reducible
bBG, the maximum energy of the final photon in the CM frame $E_{\gamma}^{m}$,
is given by Liang:2019zkb ; Zhang:2020fiu  (if not exceeding $\sqrt{s}/2$)

$E_{\gamma}^{m}(\theta_{\gamma})=\frac{\sqrt{s}(A\cos\theta_{1}-\sin\theta_{1})}{A(\cos\theta_{1}-\cos\theta_{\gamma})-(\sin\theta_{\gamma}+\sin\theta_{1})},$
(10)

where all angles are given in the CM frame, and
$A=(\sin\theta_{1}-\sin\theta_{2})/(\cos\theta_{1}-\cos\theta_{2})$, with
$\theta_{1}$ and $\theta_{2}$ being the polar angles corresponding to the
edges of the ECL detector. To remove the nasty bBG, the detector cut

$E_{\gamma}^{\mathrm{CM}}>E_{\gamma}^{m}$ (11)

is adopted for the final monophoton (hereafter the “bBG cut”), with
$E_{\gamma}^{\mathrm{CM}}$ being the photon energy in the CM frame.

In the gBG, the monophoton energy can be quite large around
$\theta_{\gamma}\sim 0$ region, because the gaps in the ECL are significantly
away from the beam direction. The simulations for gBG have been carried out by
Ref. Kou:2018nap  to search for an invisibly decaying dark photon. Two
different sets of detector cuts are designed to optimize the detection
efficiency for different masses of the dark photon: the “low-mass cut” and
“high-mass cut”. The “low-mass cut” can be described as $\theta_{\rm min}^{\rm
low}<\theta_{\gamma}^{\rm lab}<\theta_{\rm max}^{\rm low}$, where $\theta_{\rm
min}^{\rm low}$ and $\theta_{\rm max}^{\rm low}$ are the minimum and maximum
angles for the photon in the lab frame, respectively fitted as functions of
Duerr:2019dmv

$\displaystyle\theta_{\mathrm{min}}^{\mathrm{low}}$ $\displaystyle=$
$\displaystyle
5.399^{\circ}(E_{\gamma}^{\mathrm{CM}})^{2}/\mathrm{GeV}^{2}-58.82^{\circ}E_{\gamma}^{\mathrm{CM}}/\mathrm{GeV}+195.71^{\circ},$
(12) $\displaystyle\theta_{\mathrm{max}}^{\mathrm{low}}$ $\displaystyle=$
$\displaystyle-7.982^{\circ}(E_{\gamma}^{\mathrm{CM}})^{2}/\mathrm{GeV}^{2}+87.77^{\circ}E_{\gamma}^{\mathrm{CM}}/\mathrm{GeV}-120.6^{\circ}.$
(13)

The “high-mass cut” can be described as $\theta_{\rm min}^{\rm
high}<\theta_{\gamma}^{\rm high}<\theta_{\rm max}^{\rm high}$, where
$\theta_{\rm min}^{\rm high}$ and $\theta_{\rm max}^{\rm high}$ can be
respectively fitted as functions of Duerr:2019dmv

$\displaystyle\theta_{\mathrm{min}}^{\mathrm{high}}$ $\displaystyle=$
$\displaystyle
3.3133^{\circ}(E_{\gamma}^{\mathrm{CM}})^{2}/\mathrm{GeV}^{2}-33.58^{\circ}E_{\gamma}^{\mathrm{CM}}/\mathrm{GeV}+108.79^{\circ},$
(14) $\displaystyle\theta_{\mathrm{max}}^{\mathrm{high}}$ $\displaystyle=$
$\displaystyle-5.9133^{\circ}(E_{\gamma}^{\mathrm{CM}})^{2}/\mathrm{GeV}^{2}+54.119^{\circ}E_{\gamma}^{\mathrm{CM}}/\mathrm{GeV}-13.781^{\circ}.$
(15)

In order to probe the sensitivity for the light dark states with
electromagnetic form factors at Belle II, we define
$\chi^{2}({\mathcal{C}}_{\chi})\equiv S^{2}/(S+B)$, where $S$ ($B$) is the
number of events in the signal (background) process and
${\mathcal{C}}_{\chi}=\mu_{\chi},\,d_{\chi},\,a_{\chi},\,b_{\chi}$ denotes the
dimensional coefficient of MDM, EDM, AM and CR interactions, respectively. For
background, $B=B_{\rm ir}+B_{\rm re}$ consists of the number of events in
irreducible background $B_{\rm ir}$ and reducible background $B_{\rm re}$. The
number of events $B_{\rm ir}$ ($S$) can be obtained from the irreducible
background (signal) by integrating the differential cross section in the phase
space regions under the related detector cuts, and assuming photon detection
efficiency as 95% Kou:2018nap . It is found that about 300 (25000) gBG events
survived the “low-mass cut” (“high-mass cut”) with 20 fb-1 integrated
luminosity Kou:2018nap  in the reducible background, which will be rescaled
according to the considered luminosity.

By solving $\chi^{2}({\mathcal{C}}_{\chi})-\chi^{2}(0)=2.71$, one can obtain
the 95% confidence level (C.L.) upper limits on the corresponding
electromagnetic form factors for dark states. The upper limits under “low-mass
cut” and “high-mass cut” at Belle II with 50 ab-1 integrated luminosity are
shown in Fig. 1. We can see that the constraints under “low-mass cut” are
better than “high-mass cut” for MDM with $m_{\chi}\lesssim$ 2.6 GeV, EDM with
$m_{\chi}\lesssim$ 2.1 GeV, AM with $m_{\chi}\lesssim$ 1.4 GeV, and CR with
$m_{\chi}\lesssim$ 2.0 GeV. The upper limits of the relevant interaction
coefficients at 50 ab-1 Belle II can down to about $8.6\times
10^{-7}\,\mu_{B}$, and $3.9\times 10^{-5}$ GeV-2 for light dark state with
mass-dimension 5 and mass-dimension 6 operators, respectively.

In order to investigate the effects of the gBG, and compare with other
experiments where detailed simulations with gBG are not available, the 95%
C.L. upper limits without taking gBG into account are also presented in Fig.
1. In this scenario, the reducible background can vanish with the “bBG cut”,
and now the background events are all provided by the irreducible backgrounds
that survived the “bBG cut”. We find that the upper bound under “bBG cut” can
be down to about $1.6\times 10^{-7}\,\mu_{B}$ ($6.1\times 10^{-6}$ GeV-2) for
mass-dimension 5 (6) operators, which is about five (six) times stronger than
the one when gBG is considered under “low-mass cut”.

Figure 1:  The expected 95% C.L. upper limits on the electromagnetic form
factors for mass-dimension 5 operators (left) through MDM (solid) and EDM
(dashed), and mass-dimension 6 operators (right) through AM (solid) and CR
interaction (dashed), at Belle II under “low-mass cut” (blue), “high-mass cut”
(red) and “bBG cut” (black), with 50 ab-1 integrated luminosity.

###  IV.2 BESIII and STCF

The proposed STCF Luo:2019xqt  in China is a symmetric double ring electron-
positron collider. It is the next generation tau charm facility and successor
of the BESIII experiment, and designed to have CM energy ranging from 2 to 7
GeV. At BESIII and STCF, the cut for the final photon: $E_{\gamma}>25\ {\rm
MeV}$ in the barrel ($|z_{\gamma}|<0.8$) or $E_{\gamma}>50\ {\rm MeV}$ in the
end-caps ($0.86<|z_{\gamma}|<0.92$) Ablikim:2017ixv  is applied. Besides, we
use the BESIII detector parameters to analyze the projected constraints from
STCF due to the similarity of these two experiments. Since there is no
released analysis on gBG at BESIII so far as we know, gBG in the monophoton
reducible background at BESIII and STCF is not considered. Without taking gBG
into account, the monophoton reducible background at BESIII and STCF, mainly
arises from $e^{+}e^{-}\to\not{\gamma}\not{\gamma}\gamma$ and
$e^{+}e^{-}\to\not{e}^{+}\not{e}^{-}\gamma$ . At symmetric BESIII and STCF, we
also apply the detector cut Liu:2019ogn :

$E_{\gamma}>E_{\gamma}^{m}(\theta_{\gamma})=\frac{\sqrt{s}}{(1+{\sin\theta_{\gamma}}/{\sin\theta_{b}})},$
(16)

on the final state photon to remove the reducible background, where
$\theta_{b}$ denotes the angle at the boundary of the sub-detectors. Taking
into account the coverage of main drift chamber, electromagnetic calorimeter,
and time-of-flight, we have the polar angel $\cos{\theta_{b}}=0.95$ at BESIII
Liu:2018jdi . The photon detection efficiency is assumed as $100\%$ at both
BESIII and STCF in this work, since photon reconstruction efficiencies are all
more than 99% BESIII:2011ysp  at BESIII.

Since 2012, the monophoton trigger has been implemented at BESIII, and the
corresponding events have been collected with the luminosity of about 28 fb-1
at the CM energy ranging from 2.125 GeV to 4.95 GeV until 2021. We compute the
number of events due to signal ($S$) and backgrounds ($B$) under the applied
cuts, and define $\chi_{\rm
tot}^{2}({\mathcal{C}}_{\chi})=\sum_{i}\chi_{i}^{2}({\mathcal{C}}_{\chi})$,
where $\chi_{i}^{2}({\mathcal{C}}_{\chi})\equiv S_{i}^{2}/(S_{i}+B_{i})$ with
index $i$ denoting each BESIII colliding energy. The expected 95% C.L. upper
limits on the electromagnetic form factors of the light dark fermion $\chi$
according to about 28 fb-1 luminosity collected at BESIII are shown in Fig. 2
by demanding $\chi_{\rm tot}^{2}({\mathcal{C}}_{\chi})=\chi^{2}(0)+2.71$. Fig.
2 also shows the expected 95% C.L. upper limits with assumed 30 ab-1
luminosity at three typical colliding energies, $\sqrt{s}=2,4,7$ GeV, in
future STCF, respectively. We find that, BESIII can probe couplings down to
about $1.1\times 10^{-6}\,\mu_{B}$ for light dark states with mass-dimension 5
operators and down to $1.0\times 10^{-4}$ GeV-2 with mass-dimension 6
operators. With the same luminosity, operated at lower energy, STCF has better
sensitivity in probing the light dark fermion $\chi$ with the electromagnetic
form factors though mass-dimension 5 operators. This is because the monophoton
cross section in small mass $\chi$ production is not very dependent on the CM
energy, while in the background decreases with the increment of the CM energy.
On the contrary, for mass-dimension 6 operators, the production rates of light
dark states are even more sensitive to the center-of-mass energy, thus higher
energy STCF has better sensitivity.

Figure 2:  The expected 95% C.L. upper limits on the electromagnetic form
factors at BESIII and STCF, for mass-dimension 5 operators (left) though MDM
(solid) and EDM (dashed) , and for mass-dimension 6 operators (right) through
AM (solid) and CR interaction (dashed) . The limits from BESIII (black) is
obtained with about 28 fb-1 integrated luminosity collected at the various CM
energies ranging from 2.125 GeV to 4.95 GeV during 2012-2021. The expected
limits from STCF are shown at three typical energy points with 30 ab-1
integrated luminosity for $\sqrt{s}$= 2 (green), 4 (red), and 7 GeV (blue),
respectively.

##  V $e^{+}e^{-}$ colliders operated with $\sqrt{s}\geq M_{Z}$

###  V.1 LEP

The monophoton searches have been investigated carefully by all four LEP
experiments Workman:2022ynf , while L3 Collaboration presents the most
stringent limits on the cross section, both on $Z$-pole L3:1998uub  with an
integrated luminosity of 100 pb-1 at CM energies $\sqrt{s}=89.45-91.34$ GeV,
and off $Z$-pole with 619 pb-1 at $\sqrt{s}=188.6-209.2$ GeV L3:2003yon .
Using the L3 off Z-pole data for the monophoton searches L3:2003yon , the
bounds in the presence of the $\chi$ couplings to only photon via mass-
dimension 5 MDM and EDM operators, and mass-dimension 6 AM and CR operators
have been studied in Refs. Fortin:2011hv ; Chu:2018qrm . In this work, we
revisit the sensitivity at LEP on the $\chi$ couplings to not only photon but
also $Z$ boson via Eq. (2) using the L3 measurements both on L3:1998uub  and
off L3:2003yon  $Z$-pole. For comparison of cross sections at $Z$-pole, we
require photon energy within the range 1 GeV $<E_{\gamma}<$ 10 GeV and the
angular acceptance $45^{\circ}<\theta_{\gamma}<135^{\circ}$, following the
same event acceptance criteria as in Ref. L3:1998uub  with six data subsets.
Similarly, for the off $Z$-pole analysis, the high-energy photon should lie in
the kinematic region $14^{\circ}<\theta_{\gamma}<166^{\circ}$, and
$p_{T}^{\gamma}>0.02\sqrt{s}$, following the same event topology as described
in Ref. L3:2003yon  with eight data subsets. We obtain the bounds on the
couplings from the data subset that leads to the best constraints using
$|\sigma^{\rm SM}+\sigma^{\chi}-\sigma^{\rm exp}|\leq\delta\sigma^{\rm exp}$.
Here $\sigma^{\rm SM}$ is the SM cross section, $\sigma^{\chi}$ represents the
contribution from $\chi\bar{\chi}\gamma$ production, and $\sigma^{\rm
exp}\pm\delta\sigma^{\rm exp}$ denotes the experiment result. We find that the
off $Z$-pole measurement imposes more stringent bound than the $Z$-pole
measurement does, which can be seen in Fig. 3.

Due to the couplings with $Z$ boson via Eq. (2), the dark sectors with
$m_{\chi}<M_{Z}/2$ can now be constrained by invisible $Z$ decay. The $Z$
boson partial decay widths into $\chi\bar{\chi}$ mediated by hypercharge form
factors can be expressed as

$\displaystyle\Gamma_{Z\to\chi\bar{\chi}}=\frac{s_{W}^{2}f(M_{Z}^{2})}{16\pi
M_{Z}}\sqrt{\frac{M_{Z}^{2}-4m_{\chi}^{2}}{M_{Z}^{2}}}.$ (17)

The total width of the $Z$ boson has been measured accurately by the LEP
experiments which place a strong bound on beyond-the-SM contributions
$\Gamma_{Z\to\chi\bar{\chi}}<2.0$ MeV at 95% C.L. ALEPH:2005ab . In Fig. 3, we
show the constraints on the electromagnetic form factors for $\chi$ from the
measurement of invisible $Z$ decay at LEP. It can be found that the limits on
dark states with electromagnetic form factors by invisible $Z$ decay are
stricter than those by monophton searches at LEP for MDM with
$m_{\chi}\lesssim$ 45 GeV, EDM with $m_{\chi}\lesssim$ 40 GeV, AM with
$m_{\chi}\lesssim$ 25 GeV, and CR with $m_{\chi}\lesssim$ 40 GeV.

###  V.2 CEPC

In the following, we will focus on the future CEPC CEPCStudyGroup:2018rmc ;
CEPCStudyGroup:2018ghi . The CEPC, proposed by the Chinese high energy physics
community in 2012, is designed to run primarily at a center-of-mass energy of
240 GeV as a Higgs factory ($H$-mode) with a total luminosity of $20\
\mathrm{ab}^{-1}$ for ten years running. In addition, it will also be operated
on the $Z$-pole as a $Z$ factory ($Z$-mode) with a total luminosity of $100\
\mathrm{ab}^{-1}$ for two years running at $\sqrt{s}=91.2\ \mathrm{GeV}$,
perform a precise $WW$ threshold scan ($WW$-mode) with a total luminosity of
$\sim 6\ \mathrm{ab}^{-1}$ for one year running at $\sqrt{s}\sim$ $160\
\mathrm{GeV}$, and be upgraded to a center-of-mass energy of 360 GeV, close to
the $t\bar{t}$ threshold ($t\bar{t}$-mode) with a total luminosity of $\sim 1\
\mathrm{ab}^{-1}$ for five years CEPCPhysicsStudyGroup:2022uwl . For the
monophoton signature at CEPC, we use the cuts for the final detected photon
following the CEPC CDR CEPCStudyGroup:2018ghi : $|z_{\gamma}|<0.99$ and
$E_{\gamma}>0.1$ GeV. Similar as the symmetric STCF and BESIII, we also apply
the cut in Eq. (16) with $\cos\theta_{b}=0.99$ to remove the reducible
background events Liu:2019ogn .

We use the simple criteria $S^{2}/B=2.71$ to study the 95% C.L. upper bounds
on the couplings at CEPC, which are shown in Figure 3. Here we compute the
limits based on $20\ \mathrm{ab}^{-1}$ data in the $H$-mode, $6\
\mathrm{ab}^{-1}$ data in the $WW$-mode, $100\ \mathrm{ab}^{-1}$ data in the
$Z$-mode, and $1\ \mathrm{ab}^{-1}$ data in the $t\bar{t}$-mode. The $Z$-mode
has the best sensitivity for MDM with $m_{\chi}\lesssim$ 36 GeV, EDM with
$m_{\chi}\lesssim$ 32 GeV, AM with $m_{\chi}\lesssim$ 20 GeV, and CR with
$m_{\chi}\lesssim$ 30 GeV, due to the largest luminosity, which can probe the
couplings down to about $9.0\times 10^{-8}\,\mu_{B}$, and $3.2\times 10^{-7}$
GeV-2 for mass-dimension 5 and 6 operators. The $H$-mode has the best
sensitivity for MDM with 36 GeV $\lesssim m_{\chi}\lesssim$ 91 GeV, EDM with
32 GeV $\lesssim m_{\chi}\lesssim$ 78 GeV, AM with 20 GeV $\lesssim
m_{\chi}\lesssim$ 60 GeV, and CR with 30 GeV $\lesssim m_{\chi}\lesssim$ 79
GeV. And the corresponding couplings can be probed down to about $3.0\times
10^{-7}\,\mu_{B}$, $5.2\times 10^{-7}\,\mu_{B}$, $6.2\times 10^{-7}$ GeV-2 and
$4.6\times 10^{-7}$ GeV-2, respectively, for the case where $m_{\chi}\sim 50$
GeV by the $H$-mode running of CEPC. Although the luminosity of
$t\bar{t}$-mode is only 1/6 of that of $WW$-mode, the upper limits from
$t\bar{t}$-mode are still about twice stronger than that of $WW$-mode for
light dark states with mass-dimension operators, due to the fact that the
production cross sections for $\chi$ are larger and the SM irreducible
background is smaller in the $t\bar{t}$-mode than the $WW$-mode. The
$t\bar{t}$-mode has the best sensitivity for heavy dark states $\chi$. With
$m_{\chi}\sim 100$ GeV, $\mu_{\chi}\sim 7.0\times 10^{-7}\,\mu_{B}$,
$d_{\chi}\sim 1.6\times 10^{-6}\,\mu_{B}$, $a_{\chi}\sim 1.2\times 10^{-6}$
GeV-2, and $b_{\chi}\sim 7.0\times 10^{-7}$ GeV-2 can be probed by the
$t\bar{t}$-mode running of CEPC.

Figure 3:  The expected 95% C.L. upper limits on the electromagnetic form
factors for mass-dimension 5 operators (left) through MDM (solid) and EDM
(dashed) , and for mass-dimension 6 operators (right) through AM (solid) and
CR interaction (dashed) at CEPC in the $Z$-mode with 100 ab-1 luminosity
(black), in the $H$-mode with 20 ab-1 luminosity (blue), and in $WW$-mode with
6 ab-1 luminosity (red), and in $t\bar{t}$-mode with 1 ab-1 luminosity
(green), respectively. The constraints at LEP from monophoton searches using
the on L3:1998uub  (orange) and off L3:2003yon  (skyblue) $Z$-pole data and
the measurement of invisible decay width of the $Z$ boson ALEPH:2005ab
(purple) are also shown with shaded regions.

##  VI Discussion and conclusions

The Belle II limit (red solid line) combines the low-mass and high-mass limits
in Fig. 3, where both the bBG and the gBG are considered.

The landscape of current excluded parameter space in the plane of dark states
mass and coupling to the photon with mass-dimension 5 (left panel) operators
through MDM (solid) and EDM (dashed) and with mass-dimension 6 (right panel)
operators through AM (solid) and CR interaction (dashed) are shown in Fig. 4
by shaded regions, obtaining from terrestrial experiments, such as proton-beam
experiments CHARM-II or E613 Chu:2020ysb , monophoton searches and $Z$-boson
invisible decay at LEP and monojet searches at LHC Arina:2020mxo , and
astrophysics supernovae SN 1987A Chu:2019rok . It should be noted that we only
show several more competitive constraints in Fig. 4, the more complete results
can be found in Refs. Chu:2018qrm ; Chu:2019rok ; Chu:2020ysb ; Kling:2022ykt
. The 95% C.L. constraints on the dark states with electromagnetic form
factors derived above from the electron colliders, BESIII, STCF, Belle-II, and
CEPC are also plotted with lines in Fig. 4.

For the mass-dimension 5 operators, $Z$ boson invisible decay is most
sensitive with $m_{\chi}\lesssim 45$ (40) GeV through MDM (EDM) among the
current constraints from the terrestrial experiments mentioned above.
Monophoton search at the LEP is currently the strongest constraint in the
range of 45 GeV $\lesssim m_{\chi}\lesssim 100$ GeV through MDM. BESIII can
probe new parameter space that is previously unconstrained by other
experiments for mass $\lesssim$ 1 GeV, with 28 fb-1 data collected during
2012-2021. BESIII with the omission of the gBG (28 fb-1) only leads to a
slightly weaker limit than Belle II (50 ab-1) with gBG included for
$m_{\chi}\lesssim 1$ GeV. To compare the sensitivity on the dark states-photon
couplings from different electron colliders, we also present a Belle II limit
(green curve) with gBG omitted. Although the STCF luminosity (30 ab-1) is
lower than Belle II (50 ab-1), STCF has better sensitivity in probing the low-
mass region ($m_{\chi}\lesssim 1$ GeV) than Belle II. This is because STCF is
operated at a lower colliding energy ($\sqrt{s}=4$ GeV) where the monophoton
cross section in SM is smaller than Belle II ($\sqrt{s}=10.58$ GeV), and
$\chi\bar{\chi}\gamma$ production rate is not very dependent on the CM energy
for mass-dimension 5 operators. The about five times of magnitude difference
in sensitivity between the two Belle II limits, the cyan curve and the green
curve in Fig. 4, shows that the control on gBG is very important in probing
the electromagnetic form factors via mass-dimension 5 operators. When the
background due to the gaps in the detectors is neglected, the future CEPC can
give leading constraints than other electron colliders, which can probe the
coupling down to $9.0\times 10^{-8}\,\mu_{B}$. While with about 100 ab-1
luminosity running at 91.2 GeV, the bound only improves about dozens of
percent than STCF and Belle-II with gBG omitted. It implies that the low-
energy electron colliders can secure a place in the future to probe low-mass
light dark states with electromagnetic form factors via mass-dimension 5
operators.

With regard to the mass-dimension 6 operators, the bounds from the mono-jet
search at LHC constrain better than other current experimental sensitivity in
the plotted region in Fig.4, expect for the light dark states $\chi$ with
$m_{\chi}\lesssim 10$ MeV, which are constrained dominantly by astrophysical
bound from SN1987A. The upper limits from low-energy electron colliders, such
as BESIII, STCF and Belle II (except Belle II with gBG omitted), are all
excluded by the monophoton search at LEP. This is because that, for mass-
dimension 6 operators, the production rates of light dark states $\chi$ are
even more sensitive to the CM energy, suggesting that it is unlikely for low-
energy experiments to play any role in the foreseeable future. The high-energy
colliders, such as CEPC, can probe a vast region of the parameter space that
is previously unexplored, for the light dark states with electromagnetic form
factors via mass-dimension 6 operators in the mass region from 20 MeV to 100
GeV. Compared to current LHC bounds, the improvement on upper limits of
couplings is about six times of magnitude for the mass less than 10 GeV.

Figure 4:  The expected 95% C.L. exclusion limits on the electromagnetic form
factors at electron colliders, including BESIII, STCF, Belle II and CEPC for
mass-dimension 5 operators (left) through MDM (solid) and EDM (dashed) , and
for mass-dimension 6 operators (right) through AM (solid) and CR interaction
(dashed). The Belle II limits (cyan lines) combine the low-mass and high-mass
limits in Fig. 1, where both the bBG and the gBG are considered. The other
Belle II limits (green lines) are obtained with the “bBG” cut where the gBG is
omitted. The limits from BESIII (black lines) are obtained with about 28 fb-1
integrated luminosity collected at the various CM energies from 2.125 GeV to
4.95 GeV during 2012-2021. The STCF limits (red lines) are obtained with
$\sqrt{s}=$ 4 GeV and 30 ab-1. The CEPC limits (blue lines) combine the best
limits in $Z$-mode, $H$-mode, and $t\bar{t}$-mode in Fig. 3. The gBG is not
considered at BESIII, STCF and CEPC. The landscape of current leading
constraints are also shown with shaded regions, exploiting from terrestrial
experiments, such as proton-beam experiments CHARM-II or E613 Chu:2020ysb ,
monophoton searches and $Z$-boson invisible decay at LEP and monojet searches
at LHC Arina:2020mxo , and astrophysics supernovae SN 1987A Chu:2019rok .

###### Acknowledgements.

This work was supported in part by the National Natural Science Foundation of
China (Grants No.12105327, No. 11805001), the Natural Science Foundation of
Anhui Province (No.2208085MA10), and the Key Research Foundation of Education
Ministry of Anhui Province of China (No.KJ2021A0061).

## References

  * (1) R. Essig, J. A. Jaros, W. Wester, P. Hansson Adrian, S. Andreas, T. Averett, O. Baker, B. Batell, M. Battaglieri and J. Beacham, et al. [arXiv:1311.0029 [hep-ph]].
  * (2) X. Chu, J. Pradler and L. Semmelrock, Phys. Rev. D 99 (2019) no.1, 015040 doi:10.1103/PhysRevD.99.015040 [arXiv:1811.04095 [hep-ph]].
  * (3) B. J. Kavanagh, P. Panci and R. Ziegler, JHEP 04 (2019), 089 doi:10.1007/JHEP04(2019)089 [arXiv:1810.00033 [hep-ph]].
  * (4) D. Banerjee et al. [NA64], Phys. Rev. D 97 (2018) no.7, 072002 doi:10.1103/PhysRevD.97.072002 [arXiv:1710.00971 [hep-ex]].
  * (5) T. Åkesson et al. [LDMX], [arXiv:1808.05219 [hep-ex]].
  * (6) M. Battaglieri et al. [BDX], [arXiv:1607.01390 [hep-ex]].
  * (7) A. A. Prinz, R. Baggs, J. Ballam, S. Ecklund, C. Fertig, J. A. Jaros, K. Kase, A. Kulikov, W. G. J. Langeveld and R. Leonard, et al. Phys. Rev. Lett. 81 (1998), 1175-1178 doi:10.1103/PhysRevLett.81.1175 [arXiv:hep-ex/9804008 [hep-ex]].
  * (8) C. Bird, P. Jackson, R. V. Kowalewski and M. Pospelov, Phys. Rev. Lett. 93 (2004), 201803 doi:10.1103/PhysRevLett.93.201803 [arXiv:hep-ph/0401195 [hep-ph]].
  * (9) V. V. Anisimovsky et al. [E949], Phys. Rev. Lett. 93 (2004), 031801 doi:10.1103/PhysRevLett.93.031801 [arXiv:hep-ex/0403036 [hep-ex]].
  * (10) A. V. Artamonov et al. [BNL-E949], Phys. Rev. D 79 (2009), 092004 doi:10.1103/PhysRevD.79.092004 [arXiv:0903.0030 [hep-ex]].
  * (11) X. Chu, J. L. Kuo and J. Pradler, Phys. Rev. D 101 (2020) no.7, 075035 doi:10.1103/PhysRevD.101.075035 [arXiv:2001.06042 [hep-ph]].
  * (12) C. Athanassopoulos et al. [LSND], Nucl. Instrum. Meth. A 388 (1997), 149-172 doi:10.1016/S0168-9002(96)01155-2 [arXiv:nucl-ex/9605002 [nucl-ex]].
  * (13) A. A. Aguilar-Arevalo et al. [MiniBooNE DM], Phys. Rev. D 98 (2018) no.11, 112004 doi:10.1103/PhysRevD.98.112004 [arXiv:1807.06137 [hep-ex]].
  * (14) K. De Winter et al. [CHARM-II], Nucl. Instrum. Meth. A 278 (1989), 670 doi:10.1016/0168-9002(89)91190-X
  * (15) P. Vilain et al. [CHARM-II], Phys. Lett. B 335 (1994), 246-252 doi:10.1016/0370-2693(94)91421-4
  * (16) R. Ball, C. T. Coffin, H. R. Gustafson, L. W. Jones, M. J. Longo, T. J. Roberts, B. P. Roe, E. Wang, M. Duffy and G. Fanourakis, et al. eConf C801002 (1980), 172-174 PRINT-81-0319 (MICHIGAN).
  * (17) M. Anelli et al. [SHiP], [arXiv:1504.04956 [physics.ins-det]].
  * (18) R. Acciarri et al. [DUNE], [arXiv:1512.06148 [physics.ins-det]].
  * (19) B. Abi et al. [DUNE], [arXiv:1706.07081 [physics.ins-det]].
  * (20) F. Kling, J. L. Kuo, S. Trojanowski and Y. D. Tsai, [arXiv:2205.09137 [hep-ph]].
  * (21) X. Chu, J. L. Kuo, J. Pradler and L. Semmelrock, Phys. Rev. D 100 (2019) no.8, 083002 doi:10.1103/PhysRevD.100.083002 [arXiv:1908.00553 [hep-ph]].
  * (22) J. H. Chang, R. Essig and A. Reinert, JHEP 03 (2021), 141 doi:10.1007/JHEP03(2021)141 [arXiv:1911.03389 [hep-ph]].
  * (23) J. F. Fortin and T. M. P. Tait, Phys. Rev. D 85 (2012), 063506 doi:10.1103/PhysRevD.85.063506 [arXiv:1103.3289 [hep-ph]].
  * (24) C. Arina, A. Cheek, K. Mimasu and L. Pagani, Eur. Phys. J. C 81 (2021) no.3, 223 doi:10.1140/epjc/s10052-021-09010-1 [arXiv:2005.12789 [hep-ph]].
  * (25) D. M. Asner, T. Barnes, J. M. Bian, I. I. Bigi, N. Brambilla, I. R. Boyko, V. Bytev, K. T. Chao, J. Charles and H. X. Chen, et al. Int. J. Mod. Phys. A 24 (2009), S1-794 [arXiv:0809.1869 [hep-ex]].
  * (26) T. Abe et al. [Belle-II], [arXiv:1011.0352 [physics.ins-det]].
  * (27) X. D. Shi, X. R. Zhou, X. S. Qin and H. P. Peng, JINST 16 (2021) no.03, P03029 doi:10.1088/1748-0221/16/03/P03029 [arXiv:2011.01654 [physics.ins-det]].
  * (28) Q. Luo, W. Gao, J. Lan, W. Li and D. Xu, doi:10.18429/JACoW-IPAC2019-MOPRB031
  * (29) A. E. Bondar et al. [Charm-Tau Factory], Phys. Atom. Nucl. 76 (2013), 1072-1085 doi:10.1134/S1063778813090032
  * (30) J. B. Guimarães da Costa et al. [CEPC Study Group], [arXiv:1811.10545 [hep-ex]].
  * (31) A. Abada et al. [FCC], Eur. Phys. J. ST 228 (2019) no.2, 261-623 doi:10.1140/epjst/e2019-900045-4
  * (32) H. Aihara et al. [ILC], [arXiv:1901.09829 [hep-ex]].
  * (33) A. Robson, P. N. Burrows, N. Catalan Lasheras, L. Linssen, M. Petric, D. Schulte, E. Sicking, S. Stapnes and W. Wuensch, [arXiv:1812.07987 [physics.acc-ph]].
  * (34) O. Nicrosini and L. Trentadue, Phys. Lett. B 231 (1989), 487-491 doi:10.1016/0370-2693(89)90699-0
  * (35) G. Montagna, O. Nicrosini, F. Piccinini and L. Trentadue, Nucl. Phys. B 452 (1995), 161-172 doi:10.1016/0550-3213(95)00399-D [arXiv:hep-ph/9506258 [hep-ph]].
  * (36) R. Essig, J. Mardon, M. Papucci, T. Volansky and Y. M. Zhong, JHEP 11 (2013), 167 doi:10.1007/JHEP11(2013)167 [arXiv:1309.5084 [hep-ph]].
  * (37) B. Aubert et al. [BaBar], [arXiv:0808.0017 [hep-ex]].
  * (38) E. Kou et al. [Belle-II], PTEP 2019 (2019) no.12, 123C01 [erratum: PTEP 2020 (2020) no.2, 029201] doi:10.1093/ptep/ptz106 [arXiv:1808.10567 [hep-ex]].
  * (39) J. Liang, Z. Liu, Y. Ma and Y. Zhang, Phys. Rev. D 102 (2020) no.1, 015002 doi:10.1103/PhysRevD.102.015002 [arXiv:1909.06847 [hep-ph]].
  * (40) Y. Zhang, Z. Yu, Q. Yang, M. Song, G. Li and R. Ding, Phys. Rev. D 103 (2021) no.1, 015008 doi:10.1103/PhysRevD.103.015008 [arXiv:2012.10893 [hep-ph]].
  * (41) M. Duerr, T. Ferber, C. Hearty, F. Kahlhoefer, K. Schmidt-Hoberg and P. Tunney, JHEP 02 (2020), 039 doi:10.1007/JHEP02(2020)039 [arXiv:1911.03176 [hep-ph]].
  * (42) M. Ablikim et al. [BESIII], Phys. Rev. D 96 (2017) no.11, 112008 doi:10.1103/PhysRevD.96.112008 [arXiv:1707.05178 [hep-ex]].
  * (43) Z. Liu, Y. H. Xu and Y. Zhang, JHEP 06 (2019), 009 doi:10.1007/JHEP06(2019)009 [arXiv:1903.12114 [hep-ph]].
  * (44) Z. Liu and Y. Zhang, Phys. Rev. D 99 (2019) no.1, 015004 doi:10.1103/PhysRevD.99.015004 [arXiv:1808.00983 [hep-ph]].
  * (45) M. Ablikim et al. [BESIII], Phys. Rev. D 83 (2011), 112005 doi:10.1103/PhysRevD.83.112005 [arXiv:1103.5564 [hep-ex]].
  * (46) R. L. Workman [Particle Data Group], PTEP 2022 (2022), 083C01
  * (47) M. Acciarri et al. [L3], Phys. Lett. B 431 (1998), 199-208 doi:10.1016/S0370-2693(98)00519-X
  * (48) P. Achard et al. [L3], Phys. Lett. B 587 (2004), 16-32 doi:10.1016/j.physletb.2004.01.010 [arXiv:hep-ex/0402002 [hep-ex]].
  * (49) S. Schael et al. [ALEPH, DELPHI, L3, OPAL, SLD, LEP Electroweak Working Group, SLD Electroweak Group and SLD Heavy Flavour Group], Phys. Rept. 427 (2006), 257-454 doi:10.1016/j.physrep.2005.12.006 [arXiv:hep-ex/0509008 [hep-ex]].
  * (50) [CEPC Study Group], [arXiv:1809.00285 [physics.acc-ph]].
  * (51) H. Cheng et al. [CEPC Physics Study Group], [arXiv:2205.08553 [hep-ph]].