# A New Reduced Basis Method for Parabolic Equations Based on Single-
Eigenvalue Acceleration

Qijia Zhai, Qingguo Hong, Xiaoping Xie

###### Abstract

In this paper, we develop a new reduced basis (RB) method, named as Single
Eigenvalue Acceleration Method (SEAM), for second order parabolic equations
with homogeneous Dirichlet boundary conditions. The high-fidelity numerical
method adopts the backward Euler scheme and conforming simplicial finite
elements for the temporal and spatial discretizations, respectively. Under the
assumption that the time step size is sufficiently small and time steps are
not very large, we show that the singular value distribution of the high-
fidelity solution matrix $U$ is close to that of a rank one matrix. We select
the eigenfunction associated to the principal eigenvalue of the matrix
$U^{\top}U$ as the basis of Proper Orthogonal Decomposition (POD) method so as
to obtain SEAM and a parallel SEAM. Numerical experiments confirm the
efficiency of the new method.

Keywords: Reduced basis method; Proper orthogonal decomposition; Singular
value; Second order parabolic equation

##  1 Introduction

The Reduced Basis (RB) method is a type of model order reduction approach for
numerical approximation of problems involving repeated solution of
differential equations generated in engineering and applied sciences. It was
first proposed in [2] for the analysis of nonlinear structures in the 1970s,
and later extended to many other problems such as partial differential
equations (PDEs) with multiple parameters or different initial conditions [9,
4, 41, 10, 11, 37, 12, 39], PDE-constrained parametric optimization and
control problems [16, 34, 17, 18], and inverse problems [24, 31]. It should be
mentioned that the work in [38, 48] has led to a decisive improvement in the
computational aspects of RB methods, owing to an efficient criterion for the
selection of the basis functions, a systematic splitting of the computational
procedure into an offline (parameter-independent) phase and an online
(parameter-dependent) phase, and the use of posteriori error estimates that
guarantee certified numerical solutions for the reduced problems. The RB
methods that include the offline and online phases as their essential
constituents have become the most widely used ones.

The Proper Orthogonal Decomposition (POD) method, combined with the Galerkin
projection method, is a typical RB method. It uses a set of orthonormal bases,
which can represent the known data in the sense of least squares, to linearly
approximate the target variables so as to obtain a low-dimensional approximate
model. Since it is optimal in the least square sense, the POD method has the
property of completely relying, without making any a priori assumption, on the
data.

The POD method, whose predecessor was initially presented by K.Pearson[36] in
1901 as an eigenvector analysis method and was originally conceived in the
framework of continuous second-order processes by Berkooz[5], has been widely
used in many fields with different appellations. In the singular value
analysis and sample identification, the method is called the Karhumen-Loeve
expansion [21]. In statistics it is named as the principal component analysis
(PCA) [19]. In geophysical fluid dynamics and meteorological sciences, it is
called the empirical orthogonal function method (EOF)[33, 42]. We can also see
the widespread applications of POD method in fluid dynamics and viscous
structures[3, 8, 15, 25, 35, 40, 43, 44, 45], optimal fluid control problems
[20, 30], numerical analysis of PDEs [26, 22, 23, 1, 6, 27, 29, 28] and
machine learning [32, 47, 13, 7, 46].

In this paper, we develop a new RB/POD method, named as Single Eigenvalue
Acceleration Method (SEAM), for a full discretization, using backward Euler
scheme for temporal discretization and continuous simplicial finite elements
for spatial discretization, of second order parabolic equations with
homogeneous Dirichlet boundary conditions. The idea of SEAM is inspired by a
POD numerical experiment for a one-dimensional heat conduction problem:

$\left\\{\begin{array}[]{ll}\frac{\partial u}{\partial
t}=\frac{\partial^{2}u}{\partial x^{2}},&x\in(0,1),t\in(0,0.1],\\\
u(x,0)=\sin(4\pi x),&x\in(0,1),\\\
u(0,t)=u(1,t)=0,&t\in(0,0.1].\end{array}\right.$ (1)

When we focus on the numerical solution matrix $U$ of (1), where each column
consists of the value of the numerical solution at a node, and the number of
columns is the same as the number of time steps (We divide $(0,T]$ into 1000
equal parts and $D$ into 99 equal parts to get the high-fidelity numerical
solution; see Figure 1), we observe something interesting: the principal
singular value is much larger than other singular values of the numerical
solution matrix! Here we recall that the singular values of $U$ are the
arithmetic square roots of the eigenvalues of $U^{\top}U$. Notice that when we
use the POD method the choice of basis functions is often empirical. The
observation tells us that we only need to choose one basis function when
solving this equation with POD.

We show in Figure 1 , Figure 2 and Table 1 the SEAM (POD with one basis
function) solution and computational details, respectively. We can see that
when SEAM is used, the computational time can be saved greatly, and the
relative error between the SEAM solution and the high-fidelity solution is
very small.

(a) High-fidelity numerical solution

(b) SEAM solution

Figure 1: High-Fidelity numerical solution and SEAM solution for (1)

(a) The first 20 eigenvalues of $U^{\top}U$ (from big to small): the 1st one
is 1007.63, and the 2nd one is 8.7e-9

(b) Relative error between high-fidelity solution and SEAM solution

Figure 2:  The first 20 eigenvalues of $U^{\top}U$ and relative error between SEAM solution and high-fidelity solution of (1)  High-fidelity model | SEAM
---|---
Number of FEM d.o.fs (each time step) | 98 |  Number of SEAM d.o.fs (each time step) | 1
|  | d.o.fs reduction | 98:1
FE solution time | 1.73s | SEAM solution time | 762ms
SEAM-relative error in $L^{2}$-norm | 1.2e-8
Table 1: Computational details for the high-fidelity / SEAM solutions of (1).

Based on the above observation from the numerical experiment of 1d case, we
investigate SEAM for d-dimensional (d=2,3) second order parabolic problems
with homogeneous Dirichlet boundary conditions. Under the assumption that the
time step size is sufficiently small and the time steps are not very large, we
show that the singular value distribution of the high-fidelity solution matrix
is very close to a reference matrix of rank one. Then we select the
eigenfunction associated to the principal eigenvalue as the POD basis so as to
obtain SEAM.

The rest of the paper is arranged as follows. Section 2 introduces the model
problem and its full discretization. Section 3 gives the SEAM algorithm and
its parallel version. Section 4 discusses the singular value distribution of
numerical solution matrix, and Section 5 reports some numerical experiments.
Finally, Section 6 gives concluding remarks.

##  2 Model problem and fully discrete scheme

Let $D\subset\mathcal{R}^{d}$ ($d=1,2,3$) be an open, bounded domain, and set
$D_{T}:=D\times(0,T]$ for constant $T>0$. We consider the following second-
order parabolic problem: find $u=u(x,t)$ such that

$\begin{cases}u_{t}-\nabla\cdot(\alpha(x)\nabla u)+c(x)u=f,&\text{ in
}D_{T},\\\ u=u_{0},&\text{ on }D\times\\{0\\},\\\ u=0,&\text{ on }\partial
D\times[0,T],\end{cases}$ (2)

where $\alpha(x)=\left(\alpha_{k,j}(x)\right)_{d\times d}$ with
$\alpha_{k,j}(x)=\alpha_{j,k}(x)\in W^{1,\infty}(D)$, $c(x)\in L^{\infty}(D)$,
and $f=f(x,t):D_{T}\rightarrow\mathcal{R}$ and
$u_{0}(x):D\rightarrow\mathcal{R}$ are given functions.

Let $\mathcal{T}_{h}=\bigcup\\{T\\}$ be a conforming shape-regular simplicial
decomposition of the domain $D$, and let $\mathbb{V}_{h}\subset H_{0}^{1}(D)$
be a space of standard finite elements of arbitrary order with respect to
$\mathcal{T}_{h}$. Denote by $\left\\{x_{k}:\ k=1,2,\cdots,M\right\\}$ the set
of all interior nodes of $\mathcal{T}_{h}$, and by $\\{\varphi_{k}(x):\
k=1,2,\cdots,M\\}$ the nodal basis of $\mathbb{V}_{h}$ with

$\varphi_{k}(x_{j})=\delta_{kj}\quad\text{for }k,j=1,2,\cdots,M.$

Then we have

$\mathbb{V}_{h}=\text{span}\left\\{\varphi_{k}:\ k=1,\cdots,M\right\\}.$

We divide the time interval $(0,T)$ into a uniform grid with nodes
$t_{n}=n\tau$ and step size $\tau=T/N$ for $n=0,1,\ldots,N$. For convenience
we set $g^{n}:=g(t_{n})$ for any function $g(t)$.

Let $\pi_{h}u_{0}\in\mathbb{V}_{h}$ be the piecewise interpolation of $u_{0}$.
We consider the following fully discrete scheme for (2):

For $n=1,2,...,N$, find $u_{n,h}\in\mathbb{V}_{h}$ such that

$\frac{1}{\tau}\left(u_{n,h}-u_{n-1,h},v\right)_{D}+a\left(u_{n,h},v\right)=\left(f_{n},v\right)_{D},\quad\forall
v\in\mathbb{V}_{h},$ (3)

with initial data $u_{0,h}=\pi_{h}u^{0}$. Here $(\cdot,\cdot)_{D}$ denotes the
$L^{2}$-inner product on $D$, and $a\left(u,v\right):=\left(\alpha\nabla
u,\nabla v\right)_{D}+\left(cu_{n,h},v\right)_{D}$.

Introduce the numerical solution vector $U_{n}$, the mass matrix
$\mathcal{M}$, the total stiffness matrix $\mathcal{S}$, and the load vector
$F$, defined respectively by

$U_{n}:=\left(u_{n,h}(x_{1}),\cdots,u_{n,h}(x_{M})\right)^{\top},\quad
n=0,1,\cdots,N,$
$\mathcal{M}:=\left((\varphi_{k},\varphi_{j})_{D}\right)_{M\times
M},\quad\mathcal{S}:=\left(a(\varphi_{k},\varphi_{j})\right)_{M\times M},\quad
F:=\left(\left(f,\varphi_{1}\right),\cdots,\left(f,\varphi_{M}\right)\right)^{\top}.$

Notice that $\mathcal{M}$ and $\mathcal{S}$ are both symmetric and positive
definite. Since
$u_{n,h}(x)=\sum\limits_{k=1}^{M}u_{n,h}(x_{k})\varphi_{k}(x),$ for a given
$U_{0}$ the system (3) is of the matrix form

$(\mathcal{M}+\tau\mathcal{S})U_{n}=\mathcal{M}U_{n-1}+\tau F.$ (4)

##  3 Single-Eigenvalue Acceleration Method

In this section, we give the SEAM algorithm, based on POD, and a parallel
version. Let $U_{n}$ be the high-fidelity numerical solution vector defined by
(4) for each $n$, and let

$U:=\left(U_{0},U_{1},\cdots,U_{n}\right),\quad n=0,1,\cdots,N,$ (5)

be the numerical solution matrix of (2). Denote

$X:=U^{\top}U,$ (6)

and let

$\lambda_{0}(X)\geq\lambda_{1}(X)\geq\cdots\geq\lambda_{n}(X)\geq 0$

be the eigenvalues of $X$, then the singular values of $U$ are

$\sqrt{\lambda_{0}(X)},\sqrt{\lambda_{1}(X)},\cdots,\sqrt{\lambda_{n}(X)}.$

Set

$\widehat{\mathcal{U}}:=\operatorname{span}\left\\{U_{0},U_{1},\cdots,U_{n}\right\\},\quad\hat{d}:=\operatorname{dim}\widehat{\mathcal{U}}.$

Without loss of generality, we assume that $1\leq\hat{d}<<M$. The POD method
is to find the standard orthogonal basis functions $\beta_{j}$
$(j=1,2,\cdots,p)$ of $\widehat{\mathcal{U}}$ such that for any $p$, $1\leq$
$p\leq\hat{d}$, the mean square error between $U_{n}$ $(n=0,1,\cdots)$ and its
orthogonal projection onto the space
$\text{span}\\{\beta_{1},\cdots,\beta_{p}\\}$ is minimized on average, i.e.

$\min_{\left\\{\beta_{j}\right\\}_{j=1}^{p}\subset\widehat{\mathcal{U}}}\sum_{k=0}^{n}\left\|{U}_{k}-\sum_{j=1}^{p}\langle{U}_{k},\beta_{j}\rangle\beta_{j}\right\|_{2}^{2}.$
(7)

A solution $\left\\{\beta_{j}\right\\}_{j=1}^{p}$ to (7) is called a POD-basis
of rank $p$.

We will show later (cf. Theorem 4.1 and Corollary 4.1) that the non-zero
eigenvalues of $X$ except for the principal eigenvalue are close to zero if
$n$ is not very large and the time step size $\tau$ is small enough. Based
upon the theoretical result, we choose $p=1$ in the POD method (7). This means
that we only need to use the biggest eigenvalue $\lambda_{0}(X)$ and the
associated eigenvector $\mathbf{b}_{0}=(b_{0},b_{1},\cdots,b_{n})^{T}$ to
obtain the rank 1 POD basis (cf. [43])

$\beta_{1}=\frac{1}{\sqrt{\lambda_{0}}}\sum_{k=0}^{n}b_{k}{U}_{k}.$

For convenience, we name the POD method with $p=1$ as Single Eigenvalue
Acceleration Method (SEAM), and the corresponding SEAM scheme for the full
discretization (4) is as follows:

Given $\alpha_{0}=\beta_{1}^{\top}U_{0}\in\Re,$ for $k=1,2,\cdots,n,$ the SEAM
approximation solution $U^{seam}_{k}$ at $k$-th time step is obtained by

$\left\\{\begin{array}[]{l}\beta_{1}^{\top}(\mathcal{M}+\tau\mathcal{S})\beta_{1}\alpha_{k}=\beta_{1}^{\top}\mathcal{M}\beta_{1}\alpha_{k-1}+\tau\beta_{1}^{\top}F,\\\
U^{seam}_{k}=\alpha_{k}\beta_{1}.\end{array}\right.$ (8)

Note that the coefficients
$\beta_{1}^{\top}(\mathcal{M}+\tau\mathcal{S})\beta_{1}$ and
$\beta_{1}^{\top}\mathcal{M}\beta_{1}$ in (8) are positive constants.

By Theorem 4.1 the model-order-reduction error of SEAM (8) is (cf. [39])

$\sum_{k=0}^{n}\left\|{U}_{k}-U^{seam}_{k}\right\|_{2}^{2}=\sum_{k=1}^{n}\lambda_{k}(X)\leq
Cn^{2}\tau,$

where $C$ is a positive constant independent of $n$ and $\tau$. This means
that the SEAM algorithm may be inaccurate when $n^{2}\tau$ is not small
enough. In other words, we can not expect the desired accuracy if $n$ is
large. To overcome such a weakness, we will give a parallel SEAM below (cf.
Remark 4.3).

For simplicity of notation, we assume $N+1=(\tilde{n}+1)(n+1)$. We divide the
total numerical solution matrix as

$\left(U_{0},U_{1},\cdots,U_{N}\right)=\left(\mathfrak{U}_{0},\mathfrak{U}_{1},\cdots,\mathfrak{U}_{\tilde{n}}\right),$

with the submatrices

$\mathfrak{U}_{k}:=(U_{k(n+1)},U_{k(n+1)+1},\cdots,U_{k(n+1)+n}),\quad
k=0,1,\cdots,\tilde{n}.$ (9)

Then the parallel SEAM can be described as following:

Algorithm 1 Parallel SEAM

1:for $k=0:\tilde{n}$ do

2: $\beta_{k}\leftarrow POD(\mathfrak{U}_{k})$

3:end for

4:$\alpha_{k,0}=\beta_{k}^{\top}U_{k(n+1)},\quad k=0,1,\cdots,\tilde{n}$

5:for $j=1:n$ do

6:
$\beta_{k}^{\top}(\mathcal{M}+\tau\mathcal{S})\beta_{k}\alpha_{k,j}=\beta_{k}^{\top}\mathcal{M}\beta_{k}\alpha_{k,j-1}+\tau\beta_{k}^{\top}F,\quad
k=0,1,\cdots,\tilde{n}$

7:end for

8:output
$\mathfrak{U}^{seam}_{k}=\beta_{k}\left(\alpha_{k,0},\alpha_{k,1},\cdots,\alpha_{k,n}\right),\quad
k=0,1,\cdots,\tilde{n}.$

##  4 Singular Value Distribution of Numerical Solution Matrix

In this section, we apply a perturbation technique to investigate the singular
value distribution of the numerical solution matrix $U$ in (5). We need the
following result due to Hoffman and Wielandt [14]:

###### Lemma 4.1

If $A$ and $E$ are two $(n+1)\times(n+1)$ symmetric matrices, then there hold

$\sum_{k=0}^{n}\left(\lambda_{k}(A+E)-\lambda_{k}(A)\right)^{2}\leq\|E\|_{\mathfrak{F}}^{2}$

and

$\lambda_{n}(E)\leq\lambda_{k}(A+E)-\lambda_{k}(A)\leq\lambda_{0}(E),\quad
0\leq k\leq n,$

where $\lambda_{n}(E)\leq\lambda_{n-1}(E)\leq\cdots\leq\lambda_{0}(E)$ are the
eigenvalues of $E$, and
$\|E\|_{\mathfrak{F}}=\sqrt{\sum\limits_{k=0}^{n}\sum\limits_{j=0}^{n}E_{kj}^{2}}$
is the Frobenius norm of $E$.

Let us turn back to the matrix equation (4). Since

$\mathcal{M}+\tau\mathcal{S}=\mathcal{M}^{1/2}(I+\tau\mathcal{M}^{-1/2}\mathcal{S}\mathcal{M}^{-1/2})\mathcal{M}^{1/2},$

we rewrite (4) as

$U_{n}=\mathcal{M}^{-1/2}\left(I+\tau\bm{\mathcal{A}}\right)^{-1}\mathcal{M}^{1/2}U_{n-1}+\tau\mathcal{M}^{-1/2}\left(I+\tau\bm{\mathcal{A}}\right)^{-1}\mathcal{M}^{-1/2}F$
(10)

with $\bm{\mathcal{A}}=\mathcal{M}^{-1/2}\mathcal{S}\mathcal{M}^{-1/2}$, for a
given $U_{0}$ and $n=1,2,...,N$.

In what follows we assume that the time step $\tau\in(0,1)$ is small enough
such that

$\max\\{\tau\|\bm{\mathcal{A}}\|_{2},\|I-\tau\bm{\mathcal{A}}\|_{2}\\}<1.$
(11)

By noticing that $\bm{\mathcal{A}}$ is SPD, this assumption is reasonable.

Let $\rho(\cdot)$ represent the spectral radius of a matrix. As (11) implies
$\rho(\tau\bm{\mathcal{A}})<1$, we have

$\left(I+\tau\bm{\mathcal{A}}\right)^{-1}=\sum_{k=0}^{\infty}(-1)^{k}(\tau\bm{\mathcal{A}})^{k}=I-\tau\bm{\mathcal{A}}+O\left(\tau^{2}\right).$
(12)

In the case of $F=0$, (10) is reduced to

$\left\\{\begin{array}[]{l}\tilde{U}_{0}=U_{0},\\\
\widetilde{U}_{n}=\mathcal{M}^{-1/2}\left(I+\tau\bm{\mathcal{A}}\right)^{-1}\mathcal{M}^{1/2}\widetilde{U}_{n-1},\quad
1\leq n\leq N.\end{array}\right.$ (13)

Inspired by (12), for $n=0,1,2,\cdots$ we define $U_{n}^{*}$ by

$\left\\{\begin{array}[]{l}U_{0}^{*}=U_{0},\\\
U_{n}^{*}=\mathcal{M}^{-1/2}\left(I-\tau\bm{\mathcal{A}}\right)\mathcal{M}^{1/2}U_{*}^{n-1}=\mathcal{M}^{-1/2}\left(I-\tau\bm{\mathcal{A}}\right)^{n}\mathcal{M}^{1/2}U_{0}^{*},\end{array}\right.$
(14)

and define $\bar{U}_{n}^{*}$ by

$\left\\{\begin{array}[]{l}\bar{U}_{0}^{*}=U_{0},\\\
\bar{U}_{n}^{*}=\mathcal{M}^{-1/2}\left(I-\tau\|\bm{\mathcal{A}}\|_{2}I\right)\mathcal{M}^{1/2}\bar{U}_{*}^{n-1}=\left(1-\tau\|\bm{\mathcal{A}}\|_{2}\right)^{n}\bar{U}_{0}^{*}.\end{array}\right.$
(15)

Denote

$\widetilde{U}:=(\widetilde{U}_{0},\widetilde{U}_{1},\cdots,\widetilde{U}_{n}),\quad
U_{*}:=\left(U_{0}^{*},{U}_{1}^{*},\cdots,{U}_{n}^{*}\right),\quad\bar{U}_{*}:=\left(\bar{U}_{0}^{*},\bar{U}_{1}^{*},,\cdots,\bar{U}_{n}^{*}\right),$
$\widetilde{X}:=\widetilde{U}^{\top}\widetilde{U},\quad\bar{X}_{*}:=\bar{U}_{*}^{\top}\bar{U}_{*}.$

###### Remark 4.1

We easily see that the rank of $\bar{U}_{*}$ is at most 1, which implies that
$\bar{U}_{*}$ has at most one non-zero singular value, and thus $\bar{X}_{*}$
has at most one nonzero eigenvalue.

In the sequel we use $a\lesssim b$ to denote $a\leq Cb$ for simplicity, where
$C$ is a generic positive constant independent of $n$ and $\tau$ and may be
different at its each occurrence.

The following lemma shows that $\bar{X}_{*}$ can be viewed as a perturbation
of $\widetilde{X}$ in some sense.

###### Lemma 4.2

Assume that $\tau$ is small enough such that (11) holds, then we have

$\displaystyle\|\widetilde{X}-\bar{X}_{*}\|_{\mathfrak{F}}^{2}\lesssim
n^{4}\tau^{2}.$ (16)

#### Proof:

From (13)-(15) we get

$\displaystyle\widetilde{T}_{*}:=\widetilde{U}-U_{*}=\left(\bm{0},\tilde{\theta}_{1}U_{0},\tilde{\theta}_{2}U_{0}\cdots,\tilde{\theta}_{n}U_{0}\right),$
$\displaystyle\bar{T}_{*}:=U_{*}-\bar{U}_{*}=\left(\bm{0},\bar{\theta}_{1}U_{0},\bar{\theta}_{2}U_{0},\cdots,\bar{\theta}_{n}U_{0}\right),$

where

$\displaystyle\tilde{\theta}_{n}$
$\displaystyle:=\mathcal{M}^{-1/2}\left[\left(I-\tau\bm{\mathcal{A}}+O\left(\tau^{2}\right)\right)^{n}-\left(I-\tau\bm{\mathcal{A}}\right)^{n}\right]\mathcal{M}^{1/2}$
$\displaystyle=\mathcal{M}^{-1/2}\left[\sum\limits_{k=1}^{n}\binom{n}{k}\left(I-\tau\bm{\mathcal{A}}\right)^{n-k}O\left(\tau^{2k}\right)\right]\mathcal{M}^{1/2},$
$\displaystyle\bar{\theta}_{n}$
$\displaystyle:=\mathcal{M}^{-1/2}\left[\left(I-\tau\bm{\mathcal{A}}\right)^{n}-\left(I-\tau\|\bm{\mathcal{A}}\|_{2}I\right)^{n}\right]\mathcal{M}^{1/2}$
$\displaystyle=\tau\mathcal{M}^{-1/2}\left[\left(\|\bm{\mathcal{A}}\|_{2}I-\bm{\mathcal{A}}\right)\sum\limits_{k=1}^{n}\left(1-\tau\|\bm{\mathcal{A}}\|_{2}\right)^{k-1}\left(I-\tau\bm{\mathcal{A}}\right)^{n-k}\right]\mathcal{M}^{1/2}.$

Thus,

$\displaystyle\widetilde{X}-\bar{X}_{*}$ $\displaystyle=$
$\displaystyle(\widetilde{U}^{\top}\widetilde{U}-U_{*}^{\top}U_{*})+(U_{*}^{\top}U_{*}-\bar{U}_{*}^{\top}\bar{U}_{*})$
(17) $\displaystyle=$
$\displaystyle\left((U_{*}+\widetilde{T}_{*})^{\top}(U_{*}+\widetilde{T}_{*})-U_{*}^{\top}U_{*}\right)+\left((\bar{U}_{*}+\bar{T}_{*})^{\top}(\bar{U}_{*}+\bar{T}_{*})-\bar{U}_{*}^{\top}\bar{U}_{*}\right)$
$\displaystyle=$
$\displaystyle(U_{*}^{\top}\widetilde{T}_{*}+\widetilde{T}_{*}^{\top}U_{*}+\widetilde{T}_{*}^{\top}\widetilde{T}_{*})+(\bar{U}_{*}^{\top}\bar{T}_{*}+\bar{T}_{*}^{\top}\bar{U}_{*}+\bar{T}_{*}^{\top}\bar{T}_{*})$
$\displaystyle=:$ $\displaystyle E_{*}+\bar{E}_{*}.$

From the assumption (11) and the definitions of $U_{*}$, $\bar{U}_{*}$,
$\widetilde{T}_{*}$ and $\bar{T}_{*}$, we have the following estimates:

$\displaystyle\|U_{*}\|_{\mathfrak{F}}^{2}$ $\displaystyle\leq$
$\displaystyle\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\sum_{k=0}^{n}\|I-\tau\bm{\mathcal{A}}\|_{2}^{2k}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle(n+1)\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2},$
$\displaystyle\|\bar{U}_{*}\|_{\mathfrak{F}}^{2}$ $\displaystyle\leq$
$\displaystyle\sum_{k=0}^{n}(1-\tau\|\bm{\mathcal{A}}\|_{2})^{2k}\|U_{0}\|_{2}^{2}\leq(n+1)\|U_{0}\|_{2}^{2},$
$\displaystyle\|\widetilde{T}_{*}\|_{\mathfrak{F}}^{2}$ $\displaystyle\leq$
$\displaystyle\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\sum_{k=1}^{n}\left[\sum\limits_{j=1}^{k}\binom{k}{j}\|I-\tau\bm{\mathcal{A}}\|_{2}^{k-j}O\left(\tau^{2j}\right)\right]^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\sum_{k=1}^{n}\left[\sum\limits_{j=1}^{k}\binom{k}{j}O\left(\tau^{2j}\right)\right]^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\sum_{k=1}^{n}\left[\sum\limits_{j=1}^{k}O\left(\tau^{j}\right)\right]^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}$
$\displaystyle\lesssim$ $\displaystyle
n\tau^{2}\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2},$

and

$\displaystyle\|\bar{T}_{*}\|_{\mathfrak{F}}^{2}$
$\displaystyle\leq\tau^{2}\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\|\|\bm{\mathcal{A}}\|_{2}I-\bm{\mathcal{A}}\|_{2}^{2}\sum_{k=1}^{n}\left[\sum\limits_{j=1}^{k}\left(1-\tau\|\bm{\mathcal{A}}\|_{2}\right)^{j-1}\|I-\tau\bm{\mathcal{A}}\|^{k-j}\right]^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}$
$\displaystyle\leq\tau^{2}\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\|\|\bm{\mathcal{A}}\|_{2}I-\bm{\mathcal{A}}\|_{2}^{2}(\sum_{k=1}^{n}k^{2})\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}$
$\displaystyle\lesssim
n^{3}\tau^{2}\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\|\|\bm{\mathcal{A}}\|_{2}I-\bm{\mathcal{A}}\|_{2}^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2},$

Consequently, we obtain

$\displaystyle\|E_{*}\|_{\mathfrak{F}}^{2}$
$\displaystyle=\|U_{*}^{\top}\widetilde{T}_{*}+\widetilde{T}_{*}^{\top}U_{*}+\widetilde{T}_{*}^{\top}\widetilde{T}_{*}\|_{\mathfrak{F}}^{2}$
$\displaystyle\leq
8\|U_{*}\|_{\mathfrak{F}}^{2}\|\widetilde{T}_{*}\|_{\mathfrak{F}}^{2}+2\|\widetilde{T}_{*}\|_{\mathfrak{F}}^{4}$
$\displaystyle\lesssim
n^{2}\tau^{2}\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{4}\|\mathcal{M}^{1/2}\|_{2}^{4}\|U_{0}\|_{2}^{4}$
$\displaystyle\lesssim n^{2}\tau^{2},$
$\displaystyle\|\bar{E}_{*}\|_{\mathfrak{F}}^{2}$
$\displaystyle=\|\bar{U}_{*}^{\top}\bar{T}_{*}+\bar{T}_{*}^{\top}\bar{U}_{*}+\bar{T}_{*}^{\top}\bar{T}_{*}\|_{\mathfrak{F}}^{2}$
$\displaystyle\leq
8\|\bar{U}_{*}\|_{\mathfrak{F}}^{2}\|\bar{T}_{*}\|_{\mathfrak{F}}^{2}+2\|\bar{T}_{*}\|_{\mathfrak{F}}^{4}$
$\displaystyle\lesssim\left(n^{4}\tau^{2}\|U_{0}\|_{2}^{2}+n^{6}\tau^{4}\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\|\|\bm{\mathcal{A}}\|_{2}I-\bm{\mathcal{A}}\|_{2}^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}\right)$
$\displaystyle\quad\cdot\|\mathcal{M}^{-1/2}\|_{\mathfrak{F}}^{2}\|\|\bm{\mathcal{A}}\|_{2}I-\bm{\mathcal{A}}\|_{2}^{2}\|\mathcal{M}^{1/2}\|_{2}^{2}\|U_{0}\|_{2}^{2}$
$\displaystyle\lesssim n^{4}\tau^{2},$

which, together with (17), yields

$\|\widetilde{X}-\bar{X}_{*}\|_{\mathfrak{F}}^{2}\leq\|E_{*}\|_{\mathfrak{F}}^{2}+\|\bar{E}_{*}\|_{\mathfrak{F}}^{2}\lesssim
n^{2}\tau^{2}+n^{4}\tau^{2}\lesssim n^{4}\tau^{2}.$

$\square$

Based on Lemma 4.2, we can get an error estimate between ${X}=U^{\top}U$
$\bar{X}=\bar{U}^{\top}\bar{U}$ and :

###### Lemma 4.3

Assume that $\tau$ is small enough such that (11) holds, then we have

$\displaystyle\|X-\bar{X}_{*}\|_{\mathfrak{F}}^{2}\lesssim n^{4}\tau^{2}.$
(18)

#### Proof:

From (10) and (13) it follows

$U_{n}-\widetilde{U}_{n}=\tau\mathcal{M}^{-1/2}\sum\limits_{k=1}^{n}\left(I+\tau\bm{\mathcal{A}}\right)^{-k}\mathcal{M}^{-1/2}F,\quad
n=1,2,\cdots.$ (19)

Recalling that
$\widetilde{U}=(\widetilde{U}_{0},\widetilde{U}_{1},\cdots,\widetilde{U}_{n})$
and $U=\left(U_{0},{U}_{1},\cdots,{U}_{n}\right),$ we have

$T_{*}:=U-\widetilde{U}=\tau\left(\bm{0},\theta_{1}F,\cdots,\theta_{n}F\right),$

where

$\displaystyle{\theta}_{n}$
$\displaystyle:=\mathcal{M}^{-1/2}\sum_{k=1}^{n}\left(I+\tau\bm{\mathcal{A}}\right)^{-k}\mathcal{M}^{-1/2}.$

Then we obtain

$\displaystyle{X}-\widetilde{X}$
$\displaystyle=(\widetilde{U}+T_{*})^{\top}(\widetilde{U}+T_{*})-\widetilde{U}^{\top}\widetilde{U}$
$\displaystyle=\widetilde{U}^{\top}T_{*}+T_{*}^{\top}\widetilde{U}+T_{*}^{\top}T_{*}.$

By (11) it is easy to show that

$\|\widetilde{U}\|_{\mathfrak{F}}^{2}\lesssim
n,\quad\|T_{*}\|_{\mathfrak{F}}^{2}\lesssim
n\tau^{2}\|F\|_{\mathfrak{F}}^{2}.$

As a result, we get

$\displaystyle\|{X}-\widetilde{X}\|_{\mathfrak{F}}^{2}$ $\displaystyle\lesssim
n^{2}\tau^{2}\|F\|_{\mathfrak{F}}^{2}+n^{2}\tau^{4}\|F\|_{\mathfrak{F}}^{4}\lesssim
n^{2}\tau^{2},$

which, together with the triangle inequality and Lemma 4.2, yields the desired
estimate. $\square$

Using Remark 4.1 and Lemmas 4.1 and 4.2, we immediately obtain the following
main conclusion for the full discretization (10):

###### Theorem 4.1

Assume that $\tau$ is small enough such that (11) holds, then we have

$\displaystyle\left(\lambda_{0}({X})-\lambda_{0}(\bar{X}_{*})\right)^{2}+\sum\limits_{k=1}^{n}\lambda_{k}({X})^{2}\lesssim
n^{4}\tau^{2}.$ (20)

This theorem directly yields the following result:

###### Corollary 4.1

Assume that $\tau$ is small enough such that (11) holds, and that
$n=O(\tau^{-\frac{1}{2}+\epsilon})$ for some $\epsilon\in(0,\frac{1}{2})$,
then we have

$\displaystyle\left(\lambda_{0}(X)-\lambda_{0}(\bar{X}_{*})\right)^{2}+\sum_{k=1}^{n}\lambda_{k}^{2}(X)\lesssim\tau^{4\epsilon}.$
(21)

###### Remark 4.2

From Corollary 4.1, we can see that each singular value of the high-fidelity
solution matrix ${U}$ is a small perturbation of the corresponding singular
value of the reference matrix $\bar{U}_{*}$, provided that
$n=O(\tau^{-\frac{1}{2}+\epsilon})$. Since $\bar{U}_{*}$ has at most one non-
zero singular value, all the singular values of $U$, except for the largest
one, are then close to zero.

###### Remark 4.3

From the analysis we know that Corollary 4.1 holds not only for
$U=(U_{0},U_{1},\cdots,U_{n})$, but also for
$U=(U_{k},U_{k+1},\cdots,U_{k+n})$ with $k=0,1,\cdots$. This motivates the
parallel SEAM in Section 3.

##  5 Numerical Results

This section provides some numerical experiments to verify the performance of
SEAM for the fully discrete scheme (10) using continuous linear finite
elements. All the tests are performed on a 12th-Gen Intel 3.20 GHz Core i9
computer.

###  5.1 Two-dimensional tests

We consider the following two-dimensional problem:

$\left\\{\begin{array}[]{ll}\frac{\partial u}{\partial
t}=\alpha_{1}(x,y)\frac{\partial^{2}u}{\partial
x^{2}}+\alpha_{2}(x,y)\frac{\partial^{2}u}{\partial y^{2}}-c(x,y)u+f,&\text{
in }D\times(0,T],\\\ u=u_{0},&\text{ on }D\times\\{0\\},\\\ u=0,&\text{ on
}\partial D\times[0,T],\end{array}\right.$ (22)

where $D=(0,1)\times(0,1)$ and $T=1$. To obtain the high-fidelity numerical
solutions, we first divide the spatial domain $D$ uniformly into $32\times
32=1024$ squares, then divide each square into $2$ isosceles right triangles.
For the temporal subdivision, the coefficients $\alpha_{i}$, the initial data
$u_{0}(x,y)$ and the source term $f$, we consider three situations:

  * S1:

$\tau=1\times 10^{-4}$, $\alpha_{1}(x,y)=\alpha_{2}(x,y)=c(x,y)=1$,
$u_{0}(x,y)=\sin(\pi xy)$, and $f=0,\ xy$;

  * S2:

$\tau=1\times 10^{-4}$, $\alpha_{1}(x,y)=x^{2}$, $\alpha_{2}(x,y)=y^{2}$,
$c(x,y)=\pi^{2}(1-2x^{2}y^{2})$, $u_{0}(x,y)=\sin(\pi x)\sin(\pi y)$, and
$f=0,\ 10$;

  * S3:

$\tau=2.5\times 10^{-3}$, $\alpha_{1}(x,y)=x^{2}$, $\alpha_{2}(x,y)=y^{2}$,
$c(x,y)=\pi^{2}(1-2x^{2}y^{2})$, $u_{0}(x,y)=\sin(\pi x)\sin(\pi y)$, and
$f=0$.

In S1 and S2, we take $n=\tilde{n}=100$. In S3, we take $n=\tilde{n}=20$.
Figures 3 to 8 show the distribution of eigenvalues of
$\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$ for different $f$ in situations
S1,S2 and S3, where the numerical solution matrix $\mathfrak{U}_{i}$ is
defined by (9) for $i=1,2,\cdots,\tilde{n}$. In each figure, the label ’No.s
of time interval’ means the $\tilde{n}$ time intervals, and the label ’No.s of
eigenvalues’ means the first 5 eigenvalues in each time interval.

Figures 9 to 14 demonstrate the high-fidelity numerical solutions and the
corresponding SEAM solutions of (22).

Tables 2 gives some computational details as well as the relative errors
between the the SEAM solution $u_{SEAM}$ and the high-fidelity numerical
solution $u_{fem}$, with

$\text{SEAM-
Error}_{L^{2}}:=\left(\int_{D\times(0,T]}(u_{fem}-u_{SEAM})^{2}d{\bf
x}dt\right)^{\frac{1}{2}}/\left(\int_{D\times(0,T]}u_{fem}^{2}d{\bf
x}dt\right)^{\frac{1}{2}}.$

From Figures 3 to 8, we have the following observations on the distribution of
eigenvalues:

  * •

In each case the principal eigenvalue of
$\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$ is much larger than other
eigenvalues. For example, in the case S2 with $f=xy$, for the first 100 steps,
the principal eigenvalue is 2373.15 and the second largest one is 0.000256;

  * •

The eigenvalues of $\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$ except the
principal one are close to zero when $f=0$ and, though being not so close to
zero, are still very small. These are consistent with Corollary 4.1;

  * •

As $i$ increases from $1$ to $\tilde{n}$, the principal eigenvalues of
$\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$ gradually decrease. For example, in
the case S1 with $f=0$, the principal eigenvalues decrease from 3376.2 to
38.0. In the case S2 with $f=10$, the principal eigenvalues decrease from
2389.5 to 1258.0. And in the case S3, the principal eigenvalues decrease from
4495.9 to 328.3.

From Figures 9 to 14 and Tables 2, we have the following observations on the
SEAM solutions:

  * •

SEAM gives accurate numerical solutions for all cases;

  * •

SEAM is much faster than the original high-fidelity numerical method, due to
the remarkable reduction of the size of the discrete model from $M=961$ to
$p=1$.

(a) The 1th-20th time intervals

(b) The 21th-40th time intervals

(c) The 41th-60th time intervals

(d) The 61th-80th time intervals

(e) The 81th-100th time intervals

Figure 3: The first 5 eigenvalues of $\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$
(from big to small) for S1, $f=0$.

(a) The 1th-20th time intervals

(b) The 21th-40th time intervals

(c) The 41th-60th time intervals

(d) The 61th-80th time intervals

(e) The 81th-100th time intervals

Figure 4: The first 5 eigenvalues of $\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$
(from big to small) for S1 with $f=xy$.

(a) The 1th-20th time intervals

(b) The 21th-40th time intervals

(c) The 41th-60th time intervals

(d) The 61th-80th time intervals

(e) The 81th-100th time intervals

Figure 5: The first 5 eigenvalues of $\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$
(from big to small) for S2 with $f=0$.

(a) The 1th-20th time intervals

(b) The 21th-40th time intervals

(c) The 41th-60th time intervals

(d) The 61th-80th time intervals

(e) The 81th-100th time intervals

Figure 6: The first 5 eigenvalues of $\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$
(from big to small) for S2 with $f=10$.

(a) The 1th-20th time intervals

(b) The 21th-40th time intervals

(c) The 41th-60th time intervals

(d) The 61th-80th time intervals

(e) The 81th-100th time intervals

Figure 7: The first 5 eigenvalues of $\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$
(from big to small) for S2 with $f=xy$. Figure 8: The first 5 eigenvalues of
$\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$ (from big to small) for S3 with
$f=0$.

(a) $t=0.25$

(b) $t=0.50$

(c) $t=0.75$

(d) $t=1.00$

High-fidelity numerical solution

(e) $t=0.25$

(f) $t=0.50$

(g) $t=0.75$

(h) $t=1.00$

SEAM solution

Figure 9: High-fidelity / SEAM solutions of (22) at different time for S1 with
$f=0$.

(a) $t=0.25$

(b) $t=0.50$

(c) $t=0.75$

(d) $t=1.00$

High-fidelity numerical solution

(e) $t=0.25$

(f) $t=0.50$

(g) $t=0.75$

(h) $t=1.00$

SEAM solution

Figure 10: High-fidelity / SEAM solutions of (22) at different time for S1
with $f=xy$.

(a) $t=0.25$

(b) $t=0.50$

(c) $t=0.75$

(d) $t=1.00$

High-fidelity numerical solution

(e) $t=0.25$

(f) $t=0.50$

(g) $t=0.75$

(h) $t=1.00$

SEAM solution

Figure 11: High-fidelity / SEAM solutions of (22) at different time for S2
with $f=0$.

(a) $t=0.25$

(b) $t=0.50$

(c) $t=0.75$

(d) $t=1.00$

High-fidelity numerical solution

(e) $t=0.25$

(f) $t=0.50$

(g) $t=0.75$

(h) $t=1.00$

SEAM solution

Figure 12: High-fidelity / SEAM solutions of (22) at different time for S2
with $f=10$.

(a) $t=0.25$

(b) $t=0.50$

(c) $t=0.75$

(d) $t=1.00$

High-fidelity numerical solution

(e) $t=0.25$

(f) $t=0.50$

(g) $t=0.75$

(h) $t=1.00$

SEAM solution

Figure 13: High-fidelity / SEAM solutions of (22) at different time for S2
with $f=xy$.

(a) $t=0.25$

(b) $t=0.50$

(c) $t=0.75$

(d) $t=1.00$

High-fidelity numerical solution

(e) $t=0.25$

(f) $t=0.50$

(g) $t=0.75$

(h) $t=1.00$

SEAM solution

Figure 14: High-fidelity / SEAM solutions of (22) at different time for S3 with $f=0$. High-fidelity model | SEAM
---|---
Number of FEM d.o.fs (each $t_{n}$) | $M=961$ |  Number of SEAM d.o.fs (each $t_{n}$) | $p=1$
|  | d.o.fs reduction | $961:1$
FE solution time ($S1,f=0$) | 8419s | SEAM solution time | 12s
FE solution time ($S1,f=xy$) | 8432s | SEAM solution time | 14s
FE solution time ($S2,f=0$) | 8510s | SEAM solution time | 12s
FE solution time ($S2,f=10$) | 8529s | SEAM solution time | 11s
FE solution time ($S2,f=xy$) | 8540s | SEAM solution time | 13s
FE solution time ($S3,f=0$) | 301s | SEAM solution time | 3s
SEAM-Error${}_{L^{2}}$($S1,f=0$) |  | $1.2e{-8}$ |
SEAM-Error${}_{L^{2}}$($S1,f=xy$) |  | $1.4e{-8}$ |
SEAM-Error${}_{L^{2}}$($S2,f=0$) |  | $2.1e{-8}$ |
SEAM-Error${}_{L^{2}}$($S2,f=10$) |  | $1.6e{-8}$ |
SEAM-Error${}_{L^{2}}$($S2,f=xy$) |  | $1.2e{-8}$ |
SEAM-Error${}_{L^{2}}$($S3,f=0$) |  | $1.4e{-8}$ |
Table 2: Computational details for the high-fidelity / SEAM solutions.

###  5.2 A three-dimensional test

We consider the following three-dimensional problem:

$\left\\{\begin{array}[]{ll}\frac{\partial u}{\partial
t}=\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}+\frac{\partial^{2}u}{\partial z^{2}},&\text{ in }D\times(0,T],\\\
u=\sin(2\pi x)\sin(2\pi y)\sin(2\pi z),&\text{ on }D\times\\{0\\},\\\
u=0,&\text{ on }\partial D\times[0,T],\end{array}\right.$ (23)

where $D=(0,1)\times(0,1)\times(0,1)$ and $T=1$. To obtain the high-fidelity
numerical solution, we first divide the spatial domain $D$ uniformly into
$32\times 32\times 32=32768$ cubes, then divide each cube into $6$
tetrahedrons, and for the temporal subdivision we still take $\tau=2.5\times
10^{-3}$ and $n=\tilde{n}=20$ in the full discretization (3).

Figure 15 shows the distribution of eigenvalues of
$\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$, which is consistent with Corollary
4.1. Figure 16 shows the high-fidelity / SEAM solutions of (23), and Table 3
gives some computational details as well as the relative errors between the
the SEAM solution and the high-fidelity numerical solution. We can see that
the SEAM algorithm yields an accurate approximation solution and is much
faster than the high-fidelity numerical method, due to the remarkable
reduction of the size of the discrete model from $M=29397$ to $p=1$.

Figure 15: The first 5 eigenvalues of
$\mathfrak{U}_{i}^{\top}\mathfrak{U}_{i}$ (from big to small) for (23).

(a) $t=0.25$

(b) $t=0.50$

(c) $t=0.75$

(d) $t=1.00$

High-fidelity numerical solution

(e) $t=0.25$

(f) $t=0.50$

(g) $t=0.75$

(h) $t=1.00$

SEAM solution

Figure 16: High-fidelity / SEAM solutions of (23) at different time. High-fidelity model | SEAM
---|---
Number of FEM d.o.fs (each $t_{n}$) | $M=29397$ |  Number of SEAM d.o.fs (each $t_{n}$) | $p=1$
|  | d.o.fs reduction | $29397:1$
FE solution time | 40320s | SEAM solution time | 8s
$\text{SEAM-Error}_{L^{2}}$ | $1.2e{-8}$
Table 3: Computational details for the high-fidelity / SEAM solutions of (23).

##  6 Conclusion

We have studied the singular value distribution of the numerical solution
matrix corresponding to a full discretization of second order parabolic
equations that uses the backward Euler scheme in time and continuous
simplicial finite elements in space. Based on the property that the singular
value distribution of the matrix formed by the high-fidelity solution is
similar to a rank one matrix under the assumption that the time step size is
sufficiently small, we have proposed the SEAM algorithm by choosing the
eigenvector associated to the largest singular value as the POD basis. The
resulting RB scheme at each time step is a linear equation with only one
unknown. Numerical experiments have demonstrated the remarkable efficiency of
SEAM.

##  7 Acknowledgement

This work was supported in part by the National Natural Science Foundation of
China (12171340) and National Key R&D Program of China (2020YFA0714000).

## References

  * [1] D. Ahlman, F. Söderlund, J. Jackson, A. Kurdila, and W. Shyy, Proper orthogonal decomposition for time-dependent lid-driven cavity flows, Numerical Heat Transfer, Part B: Fundamentals, 42 (2002), pp. 285–306.
  * [2] B. O. Almroth, P. Stern, and F. A. Brogan, Automatic choice of global shape functions in structural analysis, AIAA Journal, 16 (1978), pp. 525–528.
  * [3] N. Aubry, P. Holmes, J. L. Lumley, and E. Stone, The dynamics of coherent structures in the wall region of a turbulent boundary layer, Journal of Fluid Mechanics, 192 (1988), pp. 115–173.
  * [4] A. Barrett and G. Reddien, On the reduced basis method, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 75 (1995), pp. 543–549.
  * [5] G. Berkooz, Observations on the Proper Orthogonal Decomposition, Springer New York, New York, NY, 1992, pp. 229–247.
  * [6] Y. Cao, J. Zhu, Z. Luo, and I. Navon, Reduced-order modeling of the upper tropical pacific ocean model using proper orthogonal decomposition, Computers & Mathematics with Applications, 52 (2006), pp. 1373–1386.
  * [7] W. Chen, J. S. Hesthaven, B. Junqiang, Y. Qiu, Z. Yang, and Y. Tihao, Greedy nonintrusive reduced order model for fluid dynamics, AIAA Journal, 56 (2018), pp. 4927–4943.
  * [8] D. T. Crommelin and A. J. Majda, Strategies for model reduction: Comparing different optimal bases, Journal of the Atmospheric Sciences, 61 (2004), pp. 2206–2217.
  * [9] J. L. Eftang, D. J. Knezevic, and A. T. Patera, An hp certified reduced basis method for parametrized parabolic partial differential equations, Mathematical and Computer Modelling of Dynamical Systems, 17 (2011), pp. 395–422.
  * [10] J. P. Fink and W. C. Rheinboldt, On the error behavior of the reduced basis technique for nonlinear finite element approximations, Zamm-zeitschrift Fur Angewandte Mathematik Und Mechanik, 63 (1983), pp. 21–28.
  * [11] J. P. Fink and W. C. Rheinboldt, Solution manifolds and submanifolds of parametrized equations and their discretization errors, Numer. Math., 45 (1984), p. 323–343.
  * [12] M. D. Gunzburger, Finite Element Methods for Viscous Incompressible Flows, Computer Science and Scientific Computing, Academic Press, San Diego, 1989\.
  * [13] J. Hesthaven and S. Ubbiali, Non-intrusive reduced order modeling of nonlinear problems using neural networks, Journal of Computational Physics, 363 (2018), pp. 55–78.
  * [14] A. Hoffman and H. Wielandt, The variation of the spectrum of a normal matrix, Duke Math. J., 20 (1953), pp. 37–39.
  * [15] P. Holmes, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge University Press, 2012.
  * [16] L. Iapichino, S. Ulbrich, and S. Volkwein, Multiobjective pde-constrained optimization using the reduced-basis method, Advances in Computational Mathematics, 43 (2017), pp. 945–972.
  * [17] K. Ito and S. Ravindran, A reduced basis method for control problems governed by pdes, in Control and estimation of distributed parameter systems, Springer, 1998, pp. 153–168.
  * [18] K. Ito and S. S. Ravindran, A reduced-order method for simulation and control of fluid flows, Journal of computational physics, 143 (1998), pp. 403–425.
  * [19] I. T. Jolliffe, Principal component analysis, Springer, 2011.
  * [20] R. D. Joslin, M. D. Gunzburger, R. A. Nicolaides, G. Erlebacher, and M. Y. Hussaini, Self-contained automated methodology for optimal flow control, AIAA Journal, 35 (1997), pp. 816–824.
  * [21] F. Keinosuke, Introduction to Statistical Pattern Recognition, Academic Press, 2 ed., 1990.
  * [22] Kunisch and Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numerische Mathematik, 90 (2001), pp. 117–148.
  * [23] Kunisch and Volkwein, Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics, SIAM Journal on Numerical Analysis, 40 (2002), pp. 492–515.
  * [24] C. Lieberman, K. Willcox, and O. Ghattas, Parameter and state model reduction for large-scale statistical inverse problems, SIAM Journal on Scientific Computing, 32 (2010), pp. 2523–2542.
  * [25] J. Lumley, Transition and Turbulence, Academic Press, 1981.
  * [26] Z. Luo and G. Chen, Proper Orthogonal Decomposition Methods for Partial Differential Equations, Mathematics in Science and Engineering, Academic Press, 2019.
  * [27] Z. Luo, J. Chen, J. Zhu, R. Wang, and I. M. Navon, An optimizing reduced order fds for the tropical pacific ocean reduced gravity model, International Journal for Numerical Methods in Fluids, 55 (2007), pp. 143–161.
  * [28] Z. Luo, R. Wang, and J. Zhu, Finite difference scheme based on proper orthogonal decomposition for the nonstationary navier-stokes equations, Science in China Series A: Mathematics, 50 (2007), pp. 1186–1196.
  * [29] Z. Luo, J. Zhu, R. Wang, and I. Navon, Proper orthogonal decomposition approach and error estimation of mixed finite element methods for the tropical pacific ocean reduced gravity model, Computer Methods in Applied Mechanics and Engineering, 196 (2007), pp. 4184–4195.
  * [30] H. V. Ly and H. T. Tran, Proper orthogonal decomposition for flow calculations and optimal control in a horizontal cvd reactor, Quarterly of Applied Mathematics, 60 (2002), pp. 631–656.
  * [31] Y. Maday, A. T. Patera, J. D. Penn, and M. Yano, A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics, International Journal for Numerical Methods in Engineering, 102 (2015), pp. 933–965.
  * [32] L. Mainini and K. Willcox, Surrogate modeling approach to support real-time structural assessment and decision making, AIAA Journal, 53 (2015), p. 1612 – 1626.  Cited by: 39.
  * [33] A. J. Majda, I. Timofeyev, and E. Vanden-Eijnden, Systematic strategies for stochastic mode reduction in climate, Journal of the Atmospheric Sciences, 60 (2003), pp. 1705–1722.
  * [34] A. Manzoni and S. Pagani, A certified rb method for pde-constrained parametric optimization problems, Communications in Applied and Industrial Mathematics, 10 (2019), pp. 123–152.
  * [35] P. Moin and R. D. Moser, Characteristic-eddy decomposition of turbulence in a channel, Journal of Fluid Mechanics, 200 (1989), pp. 471–509.
  * [36] K. Pearson, Liii. on lines and planes of closest fit to systems of points in space, Philosophical Magazine Series 1, 2, pp. 559–572.
  * [37] J. S. Peterson, The reduced basis method for incompressible viscous flow calculations, SIAM Journal on Scientific and Statistical Computing, 10 (1989), pp. 777–786.
  * [38] C. Prudhomme, D. V. Rovas, K. Veroy, L. Machiels, Y. Maday, A. T. Patera, and G. Turinici, Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods , Journal of Fluids Engineering, 124 (2001), pp. 70–80.
  * [39] A. Quarteroni, A. Manzoni, and F. Negri, Reduced Basis Methods for Partial Differential Equations An Introduction, Cham Springer International Publishing, 2016.
  * [40] M. Rajaee, S. K. F. Karlsson, and L. Sirovich, Low-dimensional description of free-shear-flow coherent structures and their dynamical behaviour, Journal of Fluid Mechanics, 258 (1994), pp. 1–29.
  * [41] W. C. Rheinboldt, On the theory and error estimation of the reduced basis method for multi-parameter problems, Nonlinear Analysis-theory Methods & Applications, 21 (1993), pp. 849–858.
  * [42] F. M. Selten, Baroclinic empirical orthogonal functions as basis functions in an atmospheric model, Journal of the Atmospheric Sciences, 54 (1997), pp. 2099–2114.
  * [43] Sirovich and Lawrence, Turbulence and the dynamics of coherent structures. Parts I-III., Quarterly of Applied Mathematics, 45 (1987), pp. 561–590.
  * [44] Sirovich and Lawrence, Turbulence and the dynamics of coherent structures. ii. symmetries and transformations, Quarterly of Applied Mathematics, 45 (1987), pp. 573–582.
  * [45] Sirovich and Lawrence, Turbulence and the dynamics of coherent structures. iii. dynamics and scaling, Quarterly of Applied Mathematics, 45 (1987), pp. 583–590.
  * [46] R. Swischuk, L. Mainini, B. Peherstorfer, and K. Willcox, Projection-based model reduction: Formulations for physics-based machine learning, Computers & Fluids, 179 (2019), pp. 704–717.
  * [47] E. Ulu, R. Zhang, and L. B. Kara, A data-driven investigation and estimation of optimal topologies under variable loading configurations, Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization, 4 (2016), pp. 61–72.
  * [48] K. Veroy, C. Prud’homme, D. Rovas, and A. Patera, A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations, 16th AIAA Computational Fluid Dynamics Conference, (2003).


#  $Z$ boson mixing and the mass of the $W$ boson

Yu-Pan Zeng<EMAIL_ADDRESS>School of Electronics and Information
Engineering, Guangdong Ocean University, Zhanjiang 524088, China  School of
Physics, Sun Yat-sen University, Guangzhou 510275, China  Chengfeng Cai
<EMAIL_ADDRESS>School of Physics, Sun Yat-sen University, Guangzhou
510275, China  Yu-Hang Su<EMAIL_ADDRESS>School of Physics, Sun Yat-
sen University, Guangzhou 510275, China  Hong-Hao Zhang
<EMAIL_ADDRESS>School of Physics, Sun Yat-sen University, Guangzhou
510275, China

###### Abstract

We explore the possibility of explaining the $W$ boson mass with an extra
gauge boson mixing with the $Z$ boson at tree level. Extra boson mixing with
the $Z$ boson will change the expression of the $Z$ boson mass, thus altering
the $W$ boson mass. We explore two models in this work. We find that in the
Derivative Portal Dark Matter model, there are parameters space which can give
the observed $W$ boson mass, as well as the observed Dark Matter relic
density. These parameters space can also fulfill the constraints from the
electroweak oblique parameters and Dark Matter indirect detection. In the U(1)
extension model, the kinetic mixing between extra boson and $B$ boson can also
give the observed $W$ boson mass. However, to fulfill electroweak oblique
parameters fit the kinetic mixing in the U(1) model can only contribute about
$27~{}\mathrm{MeV}$ extra mass to the Standard Model $W$ boson mass. Both
models indicate the extra vector boson with the best fit mass around
$120~{}\mathrm{GeV}$.

##  I Introduction

Recently the Collider Detector at Fermilab (CDF) Collaboration has measured
the mass of the $W$ boson to be $80.4335\pm 0.0094~{}\mathrm{GeV}$ Aaltonen
_et al._ (2022), which is deviated from Standard Model (SM) prediction of
$80.357\pm 0.006~{}\mathrm{GeV}$ Zyla  _et al._ (2020) and which seems to
indicate new physics beyond SM. There are lots of works appeared to discuss
this topic Cirigliano  _et al._ (2022); Borah  _et al._ (2022); Chowdhury  _et
al._ (2022); Arcadi and Djouadi (2022); Zhang and Feng (2022); Mondal (2022);
Nagao  _et al._ (2022); Kanemura and Yagyu (2022); Kawamura  _et al._ (2022);
Péli and Trócsányi (2022); Ghoshal  _et al._ (2022); Perez  _et al._ (2022);
Zheng  _et al._ (2022); Ahn  _et al._ (2022); Heo  _et al._ (2022); Crivellin
_et al._ (2022); Endo and Mishima (2022); Du  _et al._ (2022a); Cheung  _et
al._ (2022); Di Luzio  _et al._ (2022a); Balkin  _et al._ (2022); Biekötter
_et al._ (2022); Krasnikov (2022); Paul and Valli (2022); Babu  _et al._
(2022); Di Luzio  _et al._ (2022b); Bagnaschi  _et al._ (2022); Heckman
(2022); Lee and Yamashita (2022); Cheng  _et al._ (2022); Bahl  _et al._
(2022); Song  _et al._ (2022); Asadi  _et al._ (2022); Athron  _et al._
(2022a); Sakurai  _et al._ (2022); Fan  _et al._ (2022a); Zhu  _et al._
(2022a); Arias-Aragón  _et al._ (2022); Cacciapaglia and Sannino (2022);
Blennow  _et al._ (2022); Lu  _et al._ (2022); Strumia (2022); Athron  _et
al._ (2022b); Yang and Zhang (2022); de Blas  _et al._ (2022); Tang  _et al._
(2022); Du  _et al._ (2022b); Campagnari and Mulders (2022); Zhu  _et al._
(2022b); Fan  _et al._ (2022b); Han  _et al._ (2022); Algueró  _et al._
(2022); Ghorbani and Ghorbani (2022); Yuan  _et al._ (2022); Bhaskar  _et al._
(2022). In this work we will explore physics beyond SM which can give the
observed mass of the $W$ boson at tree level.

In the SM, the mass of the $W$ boson and the $Z$ boson are given by the Higgs
mechanism. Since the $Z$ boson is combination of the $B$ boson and the $W^{3}$
boson, which is a component of the gauge triplet $W^{i}$, the mass of the $W$
boson and the $Z$ boson are connected. Therefore it is difficult to solely
change the mass of the $W$ boson. One way to alter the mass of the $W$ boson
is to mix the $Z$ boson with an extra vector boson. Mix the $Z$ boson with
another boson will inevitably alter the mass expression of the $Z$ boson which
may alter the value of the $\mathrm{SU}(2)_{L}$ gauge coupling and thus the
mass of the $W$ boson. There are usually two kinds of mixing: direct mixing in
mass matrix and kinetic mixing. Though the normalization of the kinetic mixing
terms will result in mass mixing, we will consider two models in this work:
the Derivative Portal Dark Matter (DPDM) model Zeng (2022) and the U(1) model
Holdom (1991); Lao  _et al._ (2020). In these two models, the extra gauge
boson are connected to the SM through kinetic mixing to the $Z$ boson and the
$B$ boson respectively. The kinetic mixing will alter the mass expression of
the $Z$ boson and thus the mass of the $W$ boson at tree level. Since
electroweak oblique parameters have a strong constraint on electroweak
physics, we will consider the electroweak oblique parameters constraint on
these models. For the DPDM model, we also consider constraints from the
observed Dark Matter (DM) relic density and DM indirect detection.

This work is structured as follows: In Sec. II we generally discuss the
mechanism whereby the mixing between an extra boson and the $Z$ boson changes
the mass of the $W$ boson. In Sec. III we explore two models and discuss their
capability of altering the $W$ mass. We also explore constraints from
electroweak oblique parameters, DM relic density, and DM indirect detection.
We conclude in Sec. IV

##  II general discussion of prediction of the mass of the $W$ boson

In this section we will discuss in general how an extra boson mixing with the
$Z$ boson changes the mass of the $W$ boson. To see this we first write down
the mass of the $W$ boson $m_{W}$ and the mass of the $Z$ boson $m_{Z}$ given
by SM:

$\displaystyle m_{W}^{2}=\frac{1}{4}g^{2}v^{2},\
m_{Z}^{2}=\frac{1}{4}(g^{2}+g^{\prime 2})v^{2},$ (1)

where $g$ and $g^{\prime}$ are the gauge couplings of
$\mathrm{SU}(2)_{\mathrm{L}}$ and $\mathrm{U}(1)_{\mathrm{Y}}$ gauge symmetry.
And $v$ is the vacuum expectation value (vev) of the Higgs boson. When
choosing the Fermi coupling constant $G_{F}$, the $Z$ boson mass $m_{Z}$ and
the fine-structure constant $\alpha$ as input parameters, the $W$ boson mass
will then be determined. Because the parameters involved in the $W$ boson mass
can be determined by adding the following equations:

$\displaystyle G_{F}=\frac{1}{\sqrt{2}v^{2}},\
e=\sqrt{4\pi\alpha}=\frac{gg^{\prime}}{\sqrt{g^{2}+g^{\prime 2}}}.$ (2)

Going beyond the SM, we will mix the $Z$ boson with another vector boson.
After that the real mass of the $Z$ will be the square root of one of the
eigenvalues of the following mass matrix:

$\displaystyle\begin{pmatrix}\frac{1}{4}(g^{2}+g^{\prime 2})v^{2}&b\\\
b&a\end{pmatrix},$ (3)

where we have used $a$ and $b$ to denote some general mass terms. The
eigenvalues of the mass matrix can be written as:

$\displaystyle
m_{Z,Z^{\prime}}^{2}=\frac{1}{2}\left(\frac{1}{4}(g^{2}+g^{\prime
2})v^{2}+a\pm\sqrt{\left(\frac{1}{4}(g^{2}+g^{\prime
2})v^{2}+a\right)^{2}-a(g^{2}+g^{\prime 2})v^{2}+4b^{2}}\right).$ (4)

Define $c=\frac{1}{4}(g^{2}+g^{\prime 2})v^{2}$, then we will get compact form
of
$m_{Z,Z^{\prime}}^{2}=\frac{1}{2}\left(a+c\pm\sqrt{(a-c)^{2}+4b^{2}}\right)$.
We can see the heavier mass of $m_{Z,Z^{\prime}}$ will be bigger than both $a$
and $c$, and the lighter mass of $m_{Z,Z^{\prime}}$ will be smaller than both
$a$ and $c$. Therefore in order to have a bigger $c$, since observation of the
$W$ mass indicates larger $g$, the value of $a$ must be lager than $c$.
Therefore the mass of the $Z$ boson should correspond to the minus sign in Eq.
(4). Adopting the input parameters as $G_{F}=1.1663787\times
10^{-5}~{}\mathrm{GeV}^{-2},\ m_{Z}=91.1876~{}\mathrm{GeV},\ \alpha\approx
1/128$ Zyla  _et al._ (2020), we can draw a blue band which saturates the
observed mass of the $W$ boson in $3\sigma$ confidence level in Fig. 1.

Figure 1: Band which gives the $W$ mass between $80.4053$ and $80.4617$.

Actually we can calculate the analytic relation between $a$ and $b$ by taking
the mass of $W$ boson $m_{W}$ as an input parameter. From Eq. (4) we can
write:

$\displaystyle b^{2}$ $\displaystyle=c(a-m_{Z}^{2})+m_{Z}^{4}-m_{Z}^{2}a$ (5)
$\displaystyle=\frac{4m_{W}^{4}}{4m_{W}^{2}-e^{2}v^{2}}(a-m_{Z}^{2})+m_{Z}^{4}-m_{Z}^{2}a.$

Then we can give constraint on models beyond SM according to Eq. (5). Actually
the above discussion does not take loop corrections from SM into
consideration. Considering the loop corrections from SM we should replace
$m_{W}$ in Eq. (5) with $m_{W}-\delta m_{W}$, where $\delta m_{W}$ represents
the loop corrections to $m_{W}$ from SM.

##  III models beyond SM

In this section we will explore two models beyond SM which mix the $Z$ boson
with an extra vector boson and might give the observed $W$ boson mass. We also
consider other constraints like electroweak oblique parameters constraint, DM
relic density constraint, and DM indirect detection constraint.

###  III.1 Derivative Portal Dark Matter

The DPDM model extends the SM with an extra vector boson which links the dark
sector and the SM through its kinetic mixing with the $Z$ boson. The relevant
Lagrangian of the DPDM model can be written as Zeng (2022):

$\displaystyle\mathcal{L}=$
$\displaystyle-\frac{1}{4}Z^{\mu\nu}Z_{\mu\nu}-\frac{1}{4}Z^{\prime\mu\nu}Z^{\prime}_{\mu\nu}-\frac{\epsilon}{2}Z^{\mu\nu}Z_{\mu\nu}^{\prime}$
$\displaystyle+\sum\limits_{f}Z_{\mu}\bar{f}\gamma^{\mu}(g_{V}-g_{A}\gamma^{5})f+g_{\chi}Z_{\mu}^{\prime}\bar{\chi}\gamma^{\mu}\chi$
$\displaystyle+\frac{1}{2}m_{Z}^{2}Z_{\mu}Z^{\mu}+\frac{1}{2}m_{Z^{\prime}}^{2}Z_{\mu}^{\prime}Z^{\prime\mu}-m_{\chi}\bar{\chi}\chi.$

After normalization of the kinetic terms, the kinetic mixing between $Z$ and
$Z^{\prime}$ actually result in mass mixing between them. The kinetic terms of
the Lagrangian can be normalized by:

$\displaystyle K=\begin{pmatrix}-k_{1}&k_{2}\\\ k_{1}&k_{2}\end{pmatrix},$ (7)

where $k_{1}=1/\sqrt{2-2\epsilon}$ and $k_{2}=1/\sqrt{2+2\epsilon}$. The
normalization will result in the following mass matrix between the two vector
bosons:

$\displaystyle\begin{pmatrix}k_{1}^{2}M_{1}&k_{1}k_{2}M_{2}\\\
k_{1}k_{2}M_{2}&k_{2}^{2}M_{1}\end{pmatrix},$ (8)

where $M_{1}=m_{Z}^{2}+m_{Z^{\prime}}^{2}$ and
$M_{2}=m_{Z^{\prime}}^{2}-m_{Z}^{2}$. One can use an orthogonal matrix $O$ to
diagonalize the mass matrix, and $O$ can be defined as

$\displaystyle O=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix},\ \text{with}\tan
2\theta=\frac{2k_{1}k_{2}M_{2}}{(k_{2}^{2}-k_{1}^{2})M_{1}}.$ (9)

Therefore according to Eq. (5) we can give constraint on $m_{Z^{\prime}}$ and
$\epsilon$ as:

$\displaystyle
k_{1}k_{2}M_{2}=\sqrt{(k_{1}^{2}M_{1}-m_{\hat{Z}}^{2})(k_{2}^{2}M_{1}-m_{\hat{Z}}^{2})},$
(10)

where we have used $m_{\hat{Z}}$ to represent the experiment observed mass of
$Z$ boson, which is meant to distinguish from $m_{Z}$. Also we can use the
measured mass of $Z^{\prime}$ boson $m_{\hat{Z}^{\prime}}$ to reformulate Eq.
(10):

$\displaystyle
m_{Z}^{2}=\frac{1}{8k_{1}^{2}k_{2}^{2}}(m_{\hat{Z}}^{2}+m_{\hat{Z}^{\prime}}^{2})-\sqrt{\frac{1}{64k_{1}^{4}k_{2}^{4}}(m_{\hat{Z}}^{2}+m_{\hat{Z}^{\prime}}^{2})^{2}-\frac{1}{4k_{1}^{2}k_{2}^{2}}m_{\hat{Z}}^{2}m_{\hat{Z}^{\prime}}^{2}}.$
(11)

Apart from giving mass to the $W$ boson, we will also calculate the tree level
$S,T,U$ constraints on this model. The neutral-current coupling between the
$Z$ boson and SM fermions in the DPDM model can be written as:

$\displaystyle L_{NC,\hat{Z}ff}$
$\displaystyle=\sum\limits_{f}(-k_{2}\sin\theta-
k_{1}\cos\theta)\hat{Z}_{\mu}\bar{f}\gamma^{\mu}(g_{V}-g_{A}\gamma^{5})f$ (13)
$\displaystyle=\sum\limits_{f}(-k_{2}\sin\theta-
k_{1}\cos\theta)\hat{Z}_{\mu}\bar{f}\gamma^{\mu}\frac{e}{s_{w}c_{w}}(T^{3}_{f}\frac{1-\gamma^{5}}{2}-Q_{f}s_{w}^{2})f,$

where $\hat{Z}_{\mu}$ is the mass eigenstate of the $Z$ boson. We see the
above coupling is the same as that in the SM except for an extra factor
$(-k_{2}\sin\theta-k_{1}\cos\theta)$, and the form of the charged-current in
the DPDM is the same as that in the SM. Using the effective-lagrangian
techniques given by Burgess  _et al._ (1994):

$\displaystyle\mathcal{L}_{CC,Wff}=-\frac{e}{\sqrt{2}\hat{s}_{w}}(1-\frac{\alpha
S}{4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}+\frac{\hat{c}_{w}^{2}\alpha
T}{2(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}+\frac{\alpha
U}{8\hat{s}_{w}^{2}})\sum\limits_{ij}V_{ij}\bar{f}_{i}\gamma^{\mu}\gamma_{L}f_{j}W_{\mu}^{\dagger}+\mathrm{c.c.}$
(14)
$\displaystyle\mathcal{L}_{NC,\hat{Z}ff}=\frac{e}{\hat{s}_{w}\hat{c}_{w}}(1+\frac{\alpha
T}{2})\sum\limits_{f}\bar{f}\gamma^{\mu}[T^{3}_{f}\frac{1-\gamma^{5}}{2}-Q_{f}(\hat{s}_{w}^{2}+\frac{\alpha
S}{4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}-\frac{\hat{c}_{w}^{2}\hat{s}_{w}^{2}\alpha
T}{\hat{c}_{w}^{2}-\hat{s}_{w}^{2}})]f\hat{Z}_{\mu},$ (15)

where $\hat{s}_{w}=\sin\hat{\theta}_{w}$ and
$\hat{c}_{w}=\cos\hat{\theta}_{w}$ and they are defined by:

$\displaystyle\hat{s}_{w}\hat{c}_{w}m_{\hat{Z}}=s_{w}c_{w}\frac{1}{2}\sqrt{g^{2}+g^{\prime
2}}v=\frac{1}{2}ev=s_{w}c_{w}m_{Z}.$ (16)

Now we can write $S,\ T$ and $U$ in the DPDM model as

$\displaystyle\alpha
T=2(\frac{\hat{s}_{w}\hat{c}_{w}}{s_{w}c_{w}}(-k_{2}\sin\theta-
k_{1}\cos\theta)-1)$ (17) $\displaystyle\alpha
S=4\hat{c}_{w}^{2}\hat{s}_{w}^{2}\alpha
T+4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})(s_{w}^{2}-\hat{s}_{w}^{2})$ (18)
$\displaystyle\alpha
U=8\hat{s}_{w}^{2}(\frac{\hat{s}_{w}}{s_{w}}-1+\frac{\alpha
S}{4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}-\frac{\hat{c}_{w}^{2}\alpha
T}{2(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}).$ (19)

Then we constrain the DPDM model with global fit results given by table five
of Ref. de Blas  _et al._ (2022):

$\displaystyle S=0.005\pm 0.097,\ T=0.04\pm 0.12,\ U=0.134\pm 0.087,$ (20)

with the correlation coefficient $\rho_{ST}=0.91,\ \rho_{SU}=-0.65,\
\rho_{TU}=-0.88$.

The DPDM model can naturally escape stringent constraint from DM direct
detection due to a cancellation mechanism Cai  _et al._ (2022); Zeng (2022):
the scattering between DM and SM fermions mediated by the $Z$ and $Z^{\prime}$
bosons will cancel out in the zero momentum transfer limit (i.e., the
scattering amplitude is proportional to the transferred momentum). Since in
$t$ channel the derivative of mediators in momentum space is proportional to
the transferred momentum, models where dark sector linked to SM by the
derivative of mediators will possess cancellation mechanism. In the DPDM model
both spin dependent and spin independent direct detection interaction are
mediated by derivative of mediators, therefore both process are suppressed in
the DPDM model. In Fig. 2 we have drawn the constraints from observed DM relic
density, the observed $W$ mass and the electroweak oblique parameters.

Figure 2: The lightblue area is excluded by Planck experiment Aghanim  _et
al._ (2020). The blue line gives the observed DM relic density. The red line
gives the observed $W$ at tree level. The green line has taken the SM model
loop corrections into consideration and gives the observed $W$ mass, with the
dashed green lines correspond to the $3\sigma$ upper and lower deviation. The
red star gives the best fit of electroweak oblique parameters $STU$, and the
purple line corresponds to $\Delta\chi^{2}=6.18$ with respect to the best fit.

In Fig. 2 the red line gives the observed $W$ boson mass solely. The green
line gives the observed $W$ boson mass with SM loop corrections taken into
consideration. The dashed green lines correspond to the $3\sigma$ mass
deviated from the $W$ boson mass (i.e. $80.4335\pm 3\times
0.0094~{}\mathrm{GeV}$). The blue line saturates the observed DM relic
density, while the light blue area is excluded by Planck experiment Aghanim
_et al._ (2020). The valley of the blue line arises from the enhancement of
the DM annihilation by $Z^{\prime}$ resonance. The DM relic density is
calculated in settings $m_{\chi}=60~{}\mathrm{GeV},\ g_{\chi}=0.1$ by
numerical tools: FeynRules 2 Alloul  _et al._ (2014), MadGraph Alwall  _et
al._ (2014), and MadDM Ambrogi  _et al._ (2019). We use the following
definition of $\chi^{2}$ to fit $m_{\hat{Z}^{\prime}}$ and $\epsilon$ through
$STU$:

$\displaystyle\chi^{2}=XCov^{-1}X^{T},$ (21)

where

$\displaystyle X=\begin{pmatrix}S-0.005&T-0.04&U-0.134\end{pmatrix}$ (22)
$\displaystyle Cov=\begin{pmatrix}0.097^{2}&\rho_{ST}\times 0.097\times
0.12&\rho_{SU}\times 0.097\times 0.087\\\ \rho_{ST}\times 0.097\times
0.12&0.12^{2}&\rho_{TU}\times 0.12\times 0.087\\\ \rho_{SU}\times 0.097\times
0.087&\rho_{TU}\times 0.12\times 0.087&0.087^{2}\end{pmatrix}.$ (23)

The degree of freedom (d.o.f) of the $\chi^{2}$ is 3. The red star represents
the best fit of $STU$: $m_{\hat{Z}^{\prime}}=116.63~{}\mathrm{GeV},\
\epsilon=0.025,\ \chi^{2}=3.21$. $\Delta\chi^{2}=6.18$ with respect to the
best fit value is denoted by the purple line. From Fig. 2 we see that the red
star lies in the area circled by green lines, which means that the global fit
of $STU$ has encoded the information of the $W$ boson mass. Also it is clear
that the best fit point meets the observed DM relic density. The purple line
indicates there is large area which can give explanation to the $W$ boson
mass. To make the parameters fall into the purple circle
$m_{\hat{Z}^{\prime}}$ should satisfy $102~{}\mathrm{GeV}\lesssim
m_{\hat{Z}^{\prime}}\lesssim 155~{}\mathrm{GeV}$. To make area where
$m_{\hat{Z}^{\prime}}\lesssim 114~{}\mathrm{GeV}$ and
$m_{\hat{Z}^{\prime}}\gtrsim 132~{}\mathrm{GeV}$ not excluded by Planck
experiment, one can change the DM mass $m_{\chi}$ and thus the annihilation
resonance area will move accordingly. On the other hand, one can increase the
extra gauge coupling $g_{\chi}$ or simply not introduce dark matter in this
model. In the parameters setting adopted by Fig. 2, it is hard to find area
which is not excluded by DM indirect detection. However, by switching
$m_{\chi}$ and $g_{\chi}$ into free parameters and fixing
$m_{\hat{Z}^{\prime}}$ and $\epsilon$, one can find large areas that escape
constraints from $W$ boson mass, DM relic density and DM indirect detection at
the same time. For example, fixing $m_{\hat{Z}^{\prime}}$ and $\epsilon$ to
the best fit values obtained above
($m_{\hat{Z}^{\prime}}=116.63~{}\mathrm{GeV},\ \epsilon=0.025$), we draw
constraints on the DPDM model from DM relic density and DM indirect detection
in Fig. 3.

Figure 3: Constraints from DM relic density and DM indirect detection, where
light blue areas are excluded by the Planck experiment Aghanim  _et al._
(2020) and yellow area is excluded by the Fermi-LAT experiment Albert  _et
al._ (2017). Blue lines correspond to the observed DM relic density.

When the DM mass is larger than the $W$ boson mass, DM annihilate largely into
$W$ bosons. When the DM mass is smaller than the $W$ boson mass, DM annihilate
largely into $s$ quarks. Therefore we use these two annihilation channels to
show the DM indirect detection constraints on the DPDM model. The DM
annihilation cross section and DM indirect detection constraints are obtained
from MadDM and Fermi-LAT experiment Ambrogi  _et al._ (2019); Albert  _et al._
(2017). In Fig. 3 the light blue areas and the yellow area are excluded by the
Planck experiment and Fermi-LAT experiment respectively, and the blue lines
correspond to the observed DM relic density.

Note that collider experiments have set strong constraints on extra vector
boson Aaboud  _et al._ (2017); Aad  _et al._ (2019); Zyla  _et al._ (2020).
For $m_{\hat{Z}^{\prime}}<209~{}\mathrm{GeV}$, the LEP-II experiment implies
the couplings between $Z^{\prime}$ boson and the SM fermions are smaller than
or of order $10^{-2}$ Zyla  _et al._ (2020). The couplings between the
$Z^{\prime}$ boson and the SM fermions in the DPDM model can be written as:

$\displaystyle L_{NC,\hat{Z}^{\prime}ff}$
$\displaystyle=\sum\limits_{f}\hat{Z}_{\mu}^{\prime}\bar{f}\gamma^{\mu}(g_{V}^{\prime}-g_{A}^{\prime}\gamma^{5})f$
(26)
$\displaystyle=\sum\limits_{f}(-k_{1}\sin\theta+k_{2}\cos\theta)\hat{Z}_{\mu}^{\prime}\bar{f}\gamma^{\mu}(g_{V}-g_{A}\gamma^{5})f$
$\displaystyle=\sum\limits_{f}(-k_{1}\sin\theta+k_{2}\cos\theta)\hat{Z}_{\mu}^{\prime}\bar{f}\gamma^{\mu}\frac{e}{s_{w}c_{w}}(T^{3}_{f}\frac{1-\gamma^{5}}{2}-Q_{f}s_{w}^{2})f.$

Taking $f$ to be electron we can calculate $g_{V}^{\prime}=0.0005$ and
$g_{A}^{\prime}=0.0073$ for $m_{\hat{Z}^{\prime}}=116.63~{}\mathrm{GeV},\
\epsilon=0.025$. Therefore in the parameters setting we considered the DPDM
model is safe from the LEP constraint. For hadron colliders the upper limits
of couplings between the $Z^{\prime}$ boson and quarks are of order
$10^{-1}$(see summary of bounds in Fig. 88.2 of Ref. Zyla  _et al._ (2020)),
which is relatively weaker constraints since the $Z^{\prime}$ coupling to the
leptons and quarks are similar and of the same order.

For clarity, we take $m_{\hat{Z}^{\prime}}=116.63~{}\mathrm{GeV},\
\epsilon=0.025,\ m_{\chi}=1000~{}\mathrm{GeV},\ g_{\chi}=0.21$ as a benchmark
point and summarize phenomenological constraints for the benchmark point in
Table 1, where the direct detection constraint on DM-xenon scattering events
is calculated with the same procedure as Ref. Zeng  _et al._ (2022). Here for
illustration we only considered the spin-independent direct detection
scattering events, which is extremely small as expected by the cancellation
mechanism of the DPDM model.

Table 1: Summary of phenomenological constraints for the benchmark point $m_{\hat{Z}^{\prime}}=116.63~{}\mathrm{GeV},\ \epsilon=0.025,\ m_{\chi}=1000~{}\mathrm{GeV},\ g_{\chi}=0.21$. | $m_{W}$ |  $\chi^{2}$ of $STU$ | $\Omega_{DM}h^{2}$
---|---|---|---
model value | 80.4136 | 3.21 | 0.1235
constraint |  $80.4335\pm 0.0094~{}\mathrm{GeV}$ Aaltonen  _et al._ (2022) | d.o.f=3 | 0.1200$\pm$ 0.0012 Aghanim  _et al._ (2020)
| $<\sigma_{ann}v>$ |  $Z^{\prime}$-electron couplings | DM-Xe scattering events
model value | $4.5\times 10^{-26}~{}\mathrm{cm}^{3}/\mathrm{s}$ |  $g_{V}^{\prime}=0.0005$ and $g_{A}^{\prime}=0.0073$ | $1.2\times 10^{-9}$
constraint |  $2.3\times 10^{-25}~{}\mathrm{cm}^{3}/\mathrm{s}$ Ambrogi  _et al._ (2019); Albert  _et al._ (2017) |  $\sim O(10^{-2})$ Zyla  _et al._ (2020) | 7.9 Meng  _et al._ (2021)

###  III.2 U(1) model

In the U(1) model there is a gauge boson of an extra
$\mathrm{U}(1)_{\mathrm{X}}$ gauge symmetry which connects to the gauge boson
of SM $\mathrm{U}(1)_{\mathrm{Y}}$ symmetry through kinetic mixing. In this
section we will adopt the same model setting as Ref. Lao  _et al._ (2020).
Then the kinetic mixing terms can be written as:

$\displaystyle\mathcal{L}_{\mathrm{K}}=-\frac{1}{4}B^{\mu\nu}B_{\mu\nu}-\frac{1}{4}X^{\mu\nu}X_{\mu\nu}-\frac{\epsilon}{2}B^{\mu\nu}X_{\mu\nu},$
(27)

where $B_{\mu}$ and $X_{\mu}$ are the gauge fields of
$\mathrm{U}(1)_{\mathrm{Y}}$ and $\mathrm{U}(1)_{\mathrm{X}}$ gauge symmetry.
After Higgs filed getting its vev there will be mass mixing term between
$B_{\mu}$ and $W^{3}_{\mu}$, and then the matrix of $(W_{\mu}^{3},\ B_{\mu},\
X_{\mu})$ can be denoted as:

$\displaystyle\frac{1}{2}\begin{pmatrix}W^{3\mu}&B^{\mu}&X^{\mu}\end{pmatrix}\begin{pmatrix}g^{2}v^{2}/4&-gg^{\prime}v^{2}/4&0\\\
-gg^{\prime}v^{2}/4&g^{\prime 2}v^{2}/4&0\\\
0&0&g_{x}^{2}v_{s}^{2}\end{pmatrix}\begin{pmatrix}W_{\mu}^{3}\\\ B_{\mu}\\\
X_{\mu}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{pmatrix}W^{3\mu}&B^{\mu}&X^{\mu}\end{pmatrix}K^{-1T}OO^{T}K^{T}\begin{pmatrix}g^{2}v^{2}/4&-gg^{\prime}v^{2}/4&0\\\
-gg^{\prime}v^{2}/4&g^{\prime 2}v^{2}/4&0\\\
0&0&g_{x}^{2}v_{s}^{2}\end{pmatrix}KOO^{T}K^{-1}\begin{pmatrix}W_{\mu}^{3}\\\
B_{\mu}\\\ X_{\mu}\end{pmatrix}$ (28) $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{pmatrix}A^{\mu}&\hat{Z}^{\mu}&\hat{Z}^{\prime\mu}\end{pmatrix}\begin{pmatrix}0&0&0\\\
0&m_{\hat{Z}}^{2}&0\\\
0&0&m_{\hat{Z}^{\prime}}^{2}\end{pmatrix}\begin{pmatrix}A_{\mu}\\\
\hat{Z}_{\mu}\\\ \hat{Z}^{\prime}_{\mu}\end{pmatrix},$

where $g_{x}$ is gauge coupling of the $\mathrm{U}(1)_{\mathrm{X}}$ gauge
symmetry and $v_{s}$ is the vev of an dark scalar which gives mass to
$X_{\mu}$. In Eq. (28) we have used $K$ to normalize the kinetic terms of
$B_{\mu}$ and $X_{\mu}$ and used $O$ to diagonalize the mass matrix and
transform the fields to their mass eigenstates. The masses of the two massive
vector boson $\hat{Z}$ and $\hat{Z}^{\prime}$ are:

$\displaystyle m_{\hat{Z},\hat{Z}^{\prime}}^{2}$
$\displaystyle=\frac{1}{8}(g^{2}v^{2}+g^{\prime 2}k_{1}^{2}v^{2}+g^{\prime
2}k_{2}^{2}v^{2}+4g_{x}^{2}k_{1}^{2}v_{s}^{2}+4g_{x}^{2}k_{2}^{2}v_{s}^{2}$
$\displaystyle\pm\sqrt{\left(g^{2}v^{2}+\left(k_{1}^{2}+k_{2}^{2}\right)\left(g^{\prime
2}v^{2}+4g_{x}^{2}v_{s}^{2}\right)\right)^{2}-16g_{x}^{2}v^{2}v_{s}^{2}\left(g^{2}\left(k_{1}^{2}+k_{2}^{2}\right)+4g^{\prime
2}k_{1}^{2}k_{2}^{2}\right)}),$

where we denote the masses and mass eigenstates with hat in order to keep
consistent with the DPDM model. Note that the kinetic mixing between $B_{\mu}$
and $X_{\mu}$ will not change the form of the electric charge $e$. The
definition of electric can be extracted from couplings between the photon and
the Higgs doublet. Which in this model will be:

$\displaystyle
e=g[KO]_{11}=g^{\prime}[KO]_{21}=\frac{2g^{\prime}k_{2}}{\sqrt{1+\frac{4g^{\prime
2}k_{2}^{2}}{g^{2}}+\frac{k_{2}^{2}}{k_{1}^{2}}}}=\frac{gg^{\prime}}{\sqrt{g^{2}+g^{\prime
2}}},$ (30)

where $[KO]_{ij}$ represents the element which lies in the $i^{th}$ row and
the $j^{th}$ column of matrix $KO$. The neutral-current coupling between $Z$
boson and SM fermions in the U(1) model can be written as:

$\displaystyle L_{NC,\hat{Z}ff}$
$\displaystyle=\sum\limits_{f}\hat{Z}_{\mu}\bar{f}\gamma^{\mu}(g_{V}-g_{A}\gamma^{5})f,$
(31) $\displaystyle\mathrm{with}\ g_{V}$
$\displaystyle=g_{A}+g^{\prime}[KO]_{22}Q_{f},\
g_{A}=\frac{T_{f}^{3}}{2}(-g^{\prime}[KO]_{22}+g[KO]_{12}).$ (32)

From the above expression we can read $S,T,U$ as:

$\displaystyle\alpha
T=\frac{2\hat{s}_{w}\hat{c}_{w}(-g^{\prime}[KO]_{22}+g[KO]_{12})}{e}-2$ (33)
$\displaystyle\alpha
S=\frac{-4g^{\prime}[KO]_{22}(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}{-g^{\prime}[KO]_{22}+g[KO]_{12}}-4\hat{s}_{w}^{2}(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})+4\hat{c}_{w}^{2}\hat{s}_{w}^{2}\alpha
T$ (34) $\displaystyle\alpha
U=8\hat{s}_{w}^{2}(\frac{\hat{s}_{w}}{s_{w}}-1+\frac{\alpha
S}{4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}-\frac{\hat{c}_{w}^{2}\alpha
T}{2(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}).$ (35)

Now we can give a line which predicts the observed $W$ boson mass in this
model in Fig. 4.

Figure 4: The red line gives the observed $W$ boson mass at tree level. The
green line has taken the SM model loop corrections into consideration and
gives the observed $W$ boson mass, with the dashed green lines correspond to
the $3\sigma$ upper and lower deviation. The red star gives the best fit of
the electroweak oblique parameters $STU$, and the purple lines correspond to
$\Delta\chi^{2}=6.18$.

In Fig. 4 we also use red dashed line to show that the U(1) model can solely
give the observed mass of the $W$ boson. The green line takes the SM loop
corrections into consideration and gives the observed $W$ boson mass, with the
dashed green lines correspond to the $3\sigma$ upper and lower deviation. We
use the same definition of $\chi^{2}$ as Eq. (21). The red star being the best
fit of electroweak oblique parameters $STU$ :
$m_{\hat{Z}^{\prime}}=133.65~{}\mathrm{GeV},\ \epsilon=0.048,\
\chi^{2}=24.94$, with the purple lines corresponding to $\Delta\chi^{2}=6.18$.
From Fig. 4 we see that the electroweak oblique parameters results do not fall
into the area circled by the green lines which represents the direct
calculation of $m_{W}$. This means that the parameters space which gives the
observed $W$ boson mass can not give corresponding electroweak couplings that
fit the electroweak oblique parameters nicely. Our results shows that the best
fit of the U(1) model can only give about $27~{}\mathrm{MeV}$ extra mass to
the SM $W$ boson mass. See Appendix A for comparison of electroweak oblique
parameters between the DPDM model and the U(1) model.

##  IV Conclusion

In this work we have explored the possibility of altering the $W$ boson mass
at tree level through mixing between an extra gauge boson and the $Z$ boson.
We first gave general discussion of the effects from mixing extra vector boson
with the $Z$ boson, then explored two realistic models: the DPDM model and the
U(1) model. In the DPDM model the extra gauge boson mixes with the $Z$ boson
through the kinetic mixing between the extra boson and the $Z$ boson, while in
the U(1) model the extra gauge boson mixes with the $Z$ boson through the
kinetic mixing between the extra boson and the $B$ boson. Apart from giving
the $W$ boson mass, we also discussed the electroweak oblique parameters
constraints for both models, and explored DM relic density and DM indirect
detection constraints for the DPDM model. We find that in both model the best
fit value for the extra vector boson mass is around $120~{}\mathrm{GeV}$.
While the best fit of the U(1) model can only contribute $27~{}\mathrm{MeV}$
extra mass to the SM $W$ boson mass, the best fit of the DPDM model can give
the observed $W$ boson mass as well as the observed DM relic density. The DPDM
model can also escape stringent DM direct detection and the best fit of the
DPDM can saturate the constraints from DM indirect detection and rough
estimation of collider bounds. Detailed collider search for the DPDM model
seemed interesting and is leaved for future works.

###### Acknowledgements.

This work is supported in part by the National Natural Science Foundation of
China (NSFC) under Grant Nos. 11875327 and 11905300, the China Postdoctoral
Science Foundation under Grant No. 2018M643282, the Fundamental Research Funds
for the Central Universities, and the Sun Yat-Sen University Science
Foundation.

## Note added

During the finalizing of this manuscript, we noticed that Ref. Zhang and Feng
(2022) appeared on arxiv. Ref. Zhang and Feng (2022) discusses explanation of
the $W$ boson mass with U(1) dark matter model as well as several
phenomenology constraints on DM. Our work discusses models with an extra gauge
boson which can explain the $W$ boson mass. Apart from the DPDM model, we also
discussed the U(1) model, but in different scenarios.

##  Appendix A Comparison of $STU$ between the DPDM model and the U(1) model

To see why the $STU$ behave differently in the DPDM model and the U(1) model
analytically, we compare the expression and structure of the $STU$ between
these two models in this section.

In Sec. III.2 we gave the general expression of $STU$ of the U(1) model, here
we will show the exact expression of these variables. To do this we take $K$
in Eq. (28) as:

$\displaystyle K=\begin{pmatrix}1&0&0\\\ 0&1&-t_{\epsilon}\\\
0&0&1/c_{\epsilon}\end{pmatrix},$ (36)

where $c_{\epsilon}=\sqrt{1-\epsilon^{2}}$ and
$t_{\epsilon}=\epsilon/c_{\epsilon}$. Then $O$ in Eq. (28) can be denoted as:

$\displaystyle O=O_{w}O_{\xi}=\begin{pmatrix}s_{w}&-c_{w}&0\\\
c_{w}&s_{w}&0\\\ 0&0&1\end{pmatrix}\begin{pmatrix}1&&\\\ &c_{\xi}&s_{\xi}\\\
&-s_{\xi}&c_{\xi}\end{pmatrix},$ (37)

where $s_{w},\ c_{w},\ s_{\xi}$ and $c_{\xi}$ are shorthand notations for
$\sin\theta_{w},\ \cos\theta_{w},\ \sin\xi$ and $\cos\xi$ respectively.
Substituting $K$ and $O$ into the expression of the $STU$ of the U(1) model
and after simplification we arrive:

$\displaystyle\alpha T$
$\displaystyle=-2c_{\xi}\frac{s_{w}c_{w}}{\hat{s}_{w}\hat{c}_{w}}-2$ (38)
$\displaystyle\alpha S$
$\displaystyle=4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})(s_{w}^{2}-\hat{s}_{w}^{2}+\frac{s_{w}^{2}c_{w}^{2}-\hat{s}_{w}^{2}\hat{c}_{w}^{2}}{s_{w}^{2}})+4\hat{c}_{w}^{2}\hat{s}_{w}^{2}\alpha
T$ (39) $\displaystyle\alpha U$
$\displaystyle=8\hat{s}_{w}^{2}(\frac{\hat{s}_{w}}{s_{w}}-1+\frac{\alpha
S}{4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}-\frac{\hat{c}_{w}^{2}\alpha
T}{2(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})})$ (41)
$\displaystyle=8\hat{s}_{w}^{2}(\frac{\hat{s}_{w}}{s_{w}}-1+(s_{w}^{2}-\hat{s}_{w}^{2}+\frac{s_{w}^{2}c_{w}^{2}-\hat{s}_{w}^{2}\hat{c}_{w}^{2}}{s_{w}^{2}}))-4\hat{s}_{w}^{2}\hat{c}_{w}^{2}\alpha
T.$

For comparison we write down the $STU$ of the DPDM model:

$\displaystyle\alpha T$
$\displaystyle=2(\frac{\hat{s}_{w}\hat{c}_{w}}{s_{w}c_{w}}(-k_{2}\sin\theta-
k_{1}\cos\theta)-1)$ (42) $\displaystyle\alpha S$
$\displaystyle=4\hat{c}_{w}^{2}\hat{s}_{w}^{2}\alpha
T+4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})(s_{w}^{2}-\hat{s}_{w}^{2})$ (43)
$\displaystyle\alpha U$
$\displaystyle=8\hat{s}_{w}^{2}(\frac{\hat{s}_{w}}{s_{w}}-1+\frac{\alpha
S}{4(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})}-\frac{\hat{c}_{w}^{2}\alpha
T}{2(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})})$ (45)
$\displaystyle=8\hat{s}_{w}^{2}(\frac{\hat{s}_{w}}{s_{w}}-1+(s_{w}^{2}-\hat{s}_{w}^{2}))-4\hat{s}_{w}^{2}\hat{c}_{w}^{2}\alpha
T.$

We see that the $STU$ expression of these two model are similar. Knowing that
$\Delta s=\hat{s}_{w}-s_{w}$ is small, we can expand $S$ and $U$ to $O(\Delta
s)$ and neglect the higher order terms. Then for the U(1) model

$\displaystyle\alpha T$
$\displaystyle=-2c_{\xi}\frac{s_{w}c_{w}}{\hat{s}_{w}\hat{c}_{w}}-2$ (46)
$\displaystyle\alpha S$
$\displaystyle=-8\hat{c}_{w}^{2}(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})\frac{\Delta
s}{\hat{s}_{w}}+4\hat{c}_{w}^{2}\hat{s}_{w}^{2}\alpha T$ (47)
$\displaystyle\alpha U$
$\displaystyle=-8\hat{s}_{w}^{2}(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})\frac{\Delta
s}{\hat{s}_{w}}-4\hat{s}_{w}^{2}\hat{c}_{w}^{2}\alpha T,$ (48)

while for the DPDM model

$\displaystyle\alpha T$
$\displaystyle=2(\frac{\hat{s}_{w}\hat{c}_{w}}{s_{w}c_{w}}(-k_{2}\sin\theta-
k_{1}\cos\theta)-1)$ (49) $\displaystyle\alpha S$
$\displaystyle=4\hat{c}_{w}^{2}\hat{s}_{w}^{2}\alpha
T-8\hat{s}_{w}(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})\Delta s$ (50)
$\displaystyle\alpha U$
$\displaystyle=8\hat{s}_{w}(\hat{c}_{w}^{2}-\hat{s}_{w}^{2})\Delta
s-4\hat{s}_{w}^{2}\hat{c}_{w}^{2}\alpha T.$ (51)

It is interesting to see that the first order expression of $S$ and $U$ in the
DPDM model are the opposite of each other. Actually $\hat{s}_{w}$ and
$\hat{c}_{w}$ can be calculated from Eq. (16), therefore we can further
simplify the expression of the $STU$ as

$\displaystyle T$
$\displaystyle=-\frac{2c_{\xi}s_{w}c_{w}}{\alpha\hat{s}_{w}\hat{c}_{w}}-\frac{2}{\alpha},\
S=-866.17\Delta s+0.716T,\ U=-263.76\Delta s-0.716T$ (52)

and

$\displaystyle T$
$\displaystyle=(\frac{2\hat{s}_{w}\hat{c}_{w}(-k_{2}\sin\theta-
k_{1}\cos\theta)}{s_{w}c_{w}\alpha}-\frac{2}{\alpha})$ (53) $\displaystyle S$
$\displaystyle=0.716T-263.76\Delta s,\ U=263.76\Delta s-0.716T.$ (54)

$\Delta s$ represents the deviation of $s_{w}$ from its SM value
$\hat{s}_{w}$. Taking $s_{w}=0.48269$, which gives the desired $g$ (thus the
desired $m_{W}$), we can obtain $\Delta s=0.00046$. Therefore the $S$ and $U$
of the U(1) model and the DPDM model can be written as $S=0.716T-0.398,\
U=-0.716T-0.121$ and $S=0.716T-0.121,\ U=-0.716T+0.121$ respectively. From Eq.
(38) we see the $T$ in the U(1) model will always be negative, which means in
the U(1) model the $S$ will be smaller than $-0.398$ when the $W$ boson mass
is satisfied. In this case the $S$ deviates largely from its central value
$0.005$, therefore to fulfill the electroweak oblique parameters fit the
region of large $W$ boson mass enhancement will be disfavored (as we see in
Fig. 4 that the best fit of the U(1) model can only provide about
$27~{}\mathrm{MeV}$ extra mass to the SM $W$ boson mass). As a comparison, the
$S$ of the DPDM model does not deviate too much from its central value. Also
the $U$ of the DPDM model has positive value, and thus it is easier to
approach it large central value $0.134$. Therefore it is reasonable to see the
DPDM model to fulfill the oblique electroweak parameters fit and the $W$ boson
mass simultaneously.

To further understand why the contours of these two model behave differently,
we further write down the $T$ of the U(1) model as:

$\displaystyle
T=\frac{2}{\alpha}\sqrt{\frac{t_{\epsilon}^{2}s_{w}^{2}}{(x^{2}-1)^{2}+t_{\epsilon}^{2}s_{w}^{2}}}x-\frac{2}{\alpha},$
(55)

where $x=\frac{s_{w}c_{w}}{\hat{s}_{w}\hat{c}_{w}}$. Noticing $(x^{2}-1)^{2}$
is much smaller than $t_{\epsilon}^{2}s_{w}^{2}$ in the region we studied, we
can expand $T$ to $O(\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}})$ as:

$\displaystyle
T=\frac{2}{\alpha}(1-\frac{1}{2}\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}})x-\frac{2}{\alpha}.$
(56)

From Fig. 4 we see that the upper purple line is similar and close to the
lower dashed green line. All points in the lower dashed green line give the
same $W$ boson mass, which means the same $g$ and the same $s_{w}$. Therefore
we expect points in the purple line have similar $s_{w}$. If along the purple
line the $T$ are almost the same, then $S$ and $U$ will also have little
changes along the purple line since they are determined by $s_{w}$ and $T$.
Then it is reasonable to see the $\chi^{2}$ contour (the purple line) and the
$m_{W}$ contour (the dashed green line) behave similarly. Taking
$s_{w}=0.48286$ (the lower dashed green line value) as an approximation and
$\epsilon=0.2$ as a benchmark, we obtain
$\frac{1}{2}\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}}=3.64\times
10^{-5}$ and $x=1-4.21\times 10^{-4}$. We see $T$ can be rewritten as:

$\displaystyle T$
$\displaystyle=\frac{2}{\alpha}(-\frac{1}{2}\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}}-(1-x)+\frac{1}{2}\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}}(1-x))$
(58)
$\displaystyle\approx\frac{2}{\alpha}(-\frac{1}{2}\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}}-(1-x)).$

$\frac{1}{2}\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}}$ decreases as
$\epsilon$ increases, and when $\epsilon=0.5$ it decreases to $4.55\times
10^{-6}$. We see the $T$ do change little when $s_{w}$ is fixed and $\epsilon$
increases. When $\epsilon$ decreases,
$\frac{1}{2}\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}s_{w}^{2}}$ can be comparable
or larger than $4.21\times 10^{-4}$, therefore to keep $\chi^{2}$ unchanged
$s_{w}$ should change accordingly. From Fig. 4 it is obvious to see $s_{w}$ in
the lower left part of the upper purple line do increases. For the DPDM model,
$T$ can be written as:

$\displaystyle
T=\frac{2}{\alpha}\frac{1}{\sqrt{1+\frac{(x^{2}-1)^{2}}{t_{\epsilon}^{2}}}}x-\frac{2}{\alpha}.$
(59)

The expression of $T$ in the DPDM model is similar to that of the U(1) model,
however the $\Delta\chi^{2}$ of the DPDM model is different from that of the
U(1) model, which means in the case of the DPDM model the minor change of the
$T$ will affect the $\chi^{2}$ value distinctively. To see this clearly we
will write down the expression of $\chi^{2}$ explicitly:

$\displaystyle\chi^{2}$
$\displaystyle=1494.83S^{2}+S(-3620.67T-2228.11U+428.445)+2500.26T^{2}$ (61)
$\displaystyle+T(3445.65U-643.634)+1415.9U^{2}-506.147U+45.7134.$

We see $\chi^{2}=45.7134$ corresponds to the SM case. For the U(1) model the
$\chi^{2}$ can be approximated as:

$\displaystyle\chi^{2}$ $\displaystyle=7.10963\times 10^{8}\Delta s^{2}+\Delta
s(-53071T-237605)$ (63) $\displaystyle+75.2286T^{2}+25.5335T+45.7134.$

While for the DPDM model the $\chi^{2}$ can be approximated as:

$\displaystyle\chi^{2}$ $\displaystyle=3.57505\times 10^{8}\Delta s^{2}+\Delta
s(-77147.6T-246508)$ (65) $\displaystyle+75.2286T^{2}+25.5335T+45.7134.$

We see the $\chi^{2}$ in the DPDM model depends more on $T$ than that of the
U(1) model. Also the $\chi^{2}$ value of the contour in the DPDM model is one
third of that in the U(1) model. Therefore the $\chi^{2}$ in the DPDM model
will vary more drastically as $T$ varies than that in the U(1) model,
resulting in a different contour shape.

## References

  * Aaltonen  _et al._ (2022) T. Aaltonen _et al._ (CDF), Science 376, 170 (2022).
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Cirigliano  _et al._ (2022) V. Cirigliano, W. Dekens, J. de Vries, E. Mereghetti, and T. Tong, (2022), arXiv:2204.08440 [hep-ph] .
  * Borah  _et al._ (2022) D. Borah, S. Mahapatra, D. Nanda, and N. Sahu, (2022), arXiv:2204.08266 [hep-ph] .
  * Chowdhury  _et al._ (2022) T. A. Chowdhury, J. Heeck, S. Saad, and A. Thapa, (2022), arXiv:2204.08390 [hep-ph] .
  * Arcadi and Djouadi (2022) G. Arcadi and A. Djouadi, (2022), arXiv:2204.08406 [hep-ph] .
  * Zhang and Feng (2022) K.-Y. Zhang and W.-Z. Feng, (2022), arXiv:2204.08067 [hep-ph] .
  * Mondal (2022) P. Mondal, (2022), arXiv:2204.07844 [hep-ph] .
  * Nagao  _et al._ (2022) K. I. Nagao, T. Nomura, and H. Okada, (2022), arXiv:2204.07411 [hep-ph] .
  * Kanemura and Yagyu (2022) S. Kanemura and K. Yagyu, (2022), arXiv:2204.07511 [hep-ph] .
  * Kawamura  _et al._ (2022) J. Kawamura, S. Okawa, and Y. Omura, (2022), arXiv:2204.07022 [hep-ph] .
  * Péli and Trócsányi (2022) Z. Péli and Z. Trócsányi, (2022), arXiv:2204.07100 [hep-ph] .
  * Ghoshal  _et al._ (2022) A. Ghoshal, N. Okada, S. Okada, D. Raut, Q. Shafi, and A. Thapa, (2022), arXiv:2204.07138 [hep-ph] .
  * Perez  _et al._ (2022) P. F. Perez, H. H. Patel, and A. D. Plascencia, (2022), arXiv:2204.07144 [hep-ph] .
  * Zheng  _et al._ (2022) M.-D. Zheng, F.-Z. Chen, and H.-H. Zhang, (2022), arXiv:2204.06541 [hep-ph] .
  * Ahn  _et al._ (2022) Y. H. Ahn, S. K. Kang, and R. Ramos, (2022), arXiv:2204.06485 [hep-ph] .
  * Heo  _et al._ (2022) Y. Heo, D.-W. Jung, and J. S. Lee, (2022), arXiv:2204.05728 [hep-ph] .
  * Crivellin  _et al._ (2022) A. Crivellin, M. Kirk, T. Kitahara, and F. Mescia, (2022), arXiv:2204.05962 [hep-ph] .
  * Endo and Mishima (2022) M. Endo and S. Mishima, (2022), arXiv:2204.05965 [hep-ph] .
  * Du  _et al._ (2022a) X. K. Du, Z. Li, F. Wang, and Y. K. Zhang, (2022a), arXiv:2204.05760 [hep-ph] .
  * Cheung  _et al._ (2022) K. Cheung, W.-Y. Keung, and P.-Y. Tseng, (2022), arXiv:2204.05942 [hep-ph] .
  * Di Luzio  _et al._ (2022a) L. Di Luzio, M. Nardecchia, and C. Toni, (2022a), arXiv:2204.05945 [hep-ph] .
  * Balkin  _et al._ (2022) R. Balkin, E. Madge, T. Menzo, G. Perez, Y. Soreq, and J. Zupan, (2022), arXiv:2204.05992 [hep-ph] .
  * Biekötter  _et al._ (2022) T. Biekötter, S. Heinemeyer, and G. Weiglein, (2022), arXiv:2204.05975 [hep-ph] .
  * Krasnikov (2022) N. V. Krasnikov, (2022), arXiv:2204.06327 [hep-ph] .
  * Paul and Valli (2022) A. Paul and M. Valli, (2022), arXiv:2204.05267 [hep-ph] .
  * Babu  _et al._ (2022) K. S. Babu, S. Jana, and V. P. K., (2022), arXiv:2204.05303 [hep-ph] .
  * Di Luzio  _et al._ (2022b) L. Di Luzio, R. Gröber, and P. Paradisi, (2022b), arXiv:2204.05284 [hep-ph] .
  * Bagnaschi  _et al._ (2022) E. Bagnaschi, J. Ellis, M. Madigan, K. Mimasu, V. Sanz, and T. You, (2022), arXiv:2204.05260 [hep-ph] .
  * Heckman (2022) J. J. Heckman, (2022), arXiv:2204.05302 [hep-ph] .
  * Lee and Yamashita (2022) H. M. Lee and K. Yamashita, (2022), arXiv:2204.05024 [hep-ph] .
  * Cheng  _et al._ (2022) Y. Cheng, X.-G. He, Z.-L. Huang, and M.-W. Li, (2022), arXiv:2204.05031 [hep-ph] .
  * Bahl  _et al._ (2022) H. Bahl, J. Braathen, and G. Weiglein, (2022), arXiv:2204.05269 [hep-ph] .
  * Song  _et al._ (2022) H. Song, W. Su, and M. Zhang, (2022), arXiv:2204.05085 [hep-ph] .
  * Asadi  _et al._ (2022) P. Asadi, C. Cesarotti, K. Fraser, S. Homiller, and A. Parikh, (2022), arXiv:2204.05283 [hep-ph] .
  * Athron  _et al._ (2022a) P. Athron, M. Bach, D. H. J. Jacob, W. Kotlarski, D. Stöckinger, and A. Voigt, (2022a), arXiv:2204.05285 [hep-ph] .
  * Sakurai  _et al._ (2022) K. Sakurai, F. Takahashi, and W. Yin, (2022), arXiv:2204.04770 [hep-ph] .
  * Fan  _et al._ (2022a) J. Fan, L. Li, T. Liu, and K.-F. Lyu, (2022a), arXiv:2204.04805 [hep-ph] .
  * Zhu  _et al._ (2022a) B.-Y. Zhu, S. Li, J.-G. Cheng, R.-L. Li, and Y.-F. Liang, (2022a), arXiv:2204.04688 [astro-ph.HE] .
  * Arias-Aragón  _et al._ (2022) F. Arias-Aragón, E. Fernández-Martínez, M. González-López, and L. Merlo, (2022), arXiv:2204.04672 [hep-ph] .
  * Cacciapaglia and Sannino (2022) G. Cacciapaglia and F. Sannino, (2022), arXiv:2204.04514 [hep-ph] .
  * Blennow  _et al._ (2022) M. Blennow, P. Coloma, E. Fernández-Martínez, and M. González-López, (2022), arXiv:2204.04559 [hep-ph] .
  * Lu  _et al._ (2022) C.-T. Lu, L. Wu, Y. Wu, and B. Zhu, (2022), arXiv:2204.03796 [hep-ph] .
  * Strumia (2022) A. Strumia, (2022), arXiv:2204.04191 [hep-ph] .
  * Athron  _et al._ (2022b) P. Athron, A. Fowlie, C.-T. Lu, L. Wu, Y. Wu, and B. Zhu, (2022b), arXiv:2204.03996 [hep-ph] .
  * Yang and Zhang (2022) J. M. Yang and Y. Zhang, (2022), arXiv:2204.04202 [hep-ph] .
  * de Blas  _et al._ (2022) J. de Blas, M. Pierini, L. Reina, and L. Silvestrini, (2022), arXiv:2204.04204 [hep-ph] .
  * Tang  _et al._ (2022) T.-P. Tang, M. Abdughani, L. Feng, Y.-L. S. Tsai, and Y.-Z. Fan, (2022), arXiv:2204.04356 [hep-ph] .
  * Du  _et al._ (2022b) X. K. Du, Z. Li, F. Wang, and Y. K. Zhang, (2022b), arXiv:2204.04286 [hep-ph] .
  * Campagnari and Mulders (2022) C. Campagnari and M. Mulders, Science 376, abm0101 (2022).
  * Zhu  _et al._ (2022b) C.-R. Zhu, M.-Y. Cui, Z.-Q. Xia, Z.-H. Yu, X. Huang, Q. Yuan, and Y. Z. Fan, (2022b), arXiv:2204.03767 [astro-ph.HE] .
  * Fan  _et al._ (2022b) Y.-Z. Fan, T.-P. Tang, Y.-L. S. Tsai, and L. Wu, (2022b), arXiv:2204.03693 [hep-ph] .
  * Han  _et al._ (2022) X.-F. Han, F. Wang, L. Wang, J. M. Yang, and Y. Zhang,  (2022), 10.1088/1674-1137/ac7c63, arXiv:2204.06505 [hep-ph] .
  * Algueró  _et al._ (2022) M. Algueró, A. Crivellin, C. A. Manzari, and J. Matias, (2022), arXiv:2201.08170 [hep-ph] .
  * Ghorbani and Ghorbani (2022) K. Ghorbani and P. Ghorbani, (2022), arXiv:2204.09001 [hep-ph] .
  * Yuan  _et al._ (2022) G.-W. Yuan, L. Zu, L. Feng, Y.-F. Cai, and Y.-Z. Fan, (2022), arXiv:2204.04183 [hep-ph] .
  * Bhaskar  _et al._ (2022) A. Bhaskar, A. A. Madathil, T. Mandal, and S. Mitra, (2022), arXiv:2204.09031 [hep-ph] .
  * Zeng (2022) Y.-P. Zeng, (2022), arXiv:2203.09462 [hep-ph] .
  * Holdom (1991) B. Holdom, Phys. Lett. B 259, 329 (1991).
  * Lao  _et al._ (2020) J. Lao, C. Cai, Z.-H. Yu, Y.-P. Zeng, and H.-H. Zhang, Phys. Rev. D 101, 095031 (2020), arXiv:2003.02516 [hep-ph] .
  * Burgess  _et al._ (1994) C. Burgess, S. Godfrey, H. König, D. London, and I. Maksymyk, Physical Review D 49, 6115 (1994).
  * Cai  _et al._ (2022) C. Cai, Y.-P. Zeng, and H.-H. Zhang, JHEP 01, 117 (2022), arXiv:2109.11499 [hep-ph] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Alloul  _et al._ (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. 185, 2250 (2014), arXiv:1310.1921 [hep-ph] .
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Ambrogi  _et al._ (2019) F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, L. Mantani, O. Mattelaer, and G. Mohlabeng, Phys. Dark Univ. 24, 100249 (2019), arXiv:1804.00044 [hep-ph] .
  * Albert  _et al._ (2017) A. Albert _et al._ (Fermi-LAT, DES), Astrophys. J. 834, 110 (2017), arXiv:1611.03184 [astro-ph.HE] .
  * Aaboud  _et al._ (2017) M. Aaboud _et al._ (ATLAS), JHEP 10, 182 (2017), arXiv:1707.02424 [hep-ex] .
  * Aad  _et al._ (2019) G. Aad _et al._ (ATLAS), Phys. Lett. B 796, 68 (2019), arXiv:1903.06248 [hep-ex] .
  * Zeng  _et al._ (2022) Y.-P. Zeng, X. Xiao, and W. Wang, Phys. Lett. B 824, 136822 (2022), arXiv:2108.11381 [hep-ph] .
  * Meng  _et al._ (2021) Y. Meng _et al._ (PandaX-4T), Phys. Rev. Lett. 127, 261802 (2021), arXiv:2107.13438 [hep-ex] .