In this paper, we present a novel spectral renormalization exponential integrator method for solving gradient flow problems. Our method is specifically designed to simultaneously satisfy  discrete analogues of the energy dissipation laws and achieve high-order accuracy in time. To accomplish this, our method first incorporates the energy dissipation law into the target gradient flow equation by introducing a time-dependent spectral renormalization (TDSR) factor. Then, the coupled equations are discretized using the spectral approximation in space and the exponential time differencing (ETD) in time. Finally, the resulting fully discrete nonlinear system is decoupled and solved using the Picard iteration at each time step. Furthermore, we introduce an extra enforcing term into the system for updating the TDSR factor, which greatly relaxes the time step size restriction of the proposed method and enhances its computational efficiency. Extensive numerical tests  with various gradient flows  are also presented  to demonstrate the accuracy and effectiveness of our method as well as its high efficiency when combined with an adaptive time-stepping strategy for long-term simulations.

Gradient flows, Energy dissipation, Time-stepping, Spectral renormalization,  Exponential integrator

35K55, 65M12, 65M15, 65F30

§ INTRODUCTION

Partial differential equations (PDEs) with gradient flow structure are commonly used to model physical phenomena in many scientific and engineering fields, including materials science [6, 11, 17]  and fluid dynamics [7, 25]. These PDE systems are derived from the energy variational principle of total free energy in different Sobolev spaces. As a result, gradient flow models typically take the following general form:

∂ϕ/∂t=-_H E[ϕ],  ∈Ω, t> 0,
where $\phi(\x,t)$ denotes the scalar-valued phase function defined over a domain $\Omega\subseteq \mathbb{R}^{d}$ ($d =1,2,3$) at time $t$, $E[\phi]$  is the free energy functional associated with the corresponding physical problem, and
$\mbox{grad}_{H} E[\phi]$ is the functional derivative of $E$ with respect to $\phi$ in the Sobolev space $H$.
This gradient flow model satisfies the energy dissipation law:

=(_H E[ϕ], ∂ϕ/∂t)_H
=-_H E[ϕ]_H^2,
where $(\cdot,\cdot)_{H}$ and $\|\cdot\|_H$ are the inner product and associated norm of the Sobolev space $H$, respectively.
This implies that the phase solution $\phi$ evolves in a way that decreases the energy functional $E[\phi]$ over time.
For evolution gradient-flows-structure equations arising in physical applications, it is important to design numerical schemes that can enforce the corresponding physics into simulations, particularly for the energy dissipation law (<ref>).
Many efforts have been devoted to developing energy stable numerical methods in the long-standing and active research field of gradient flows. These methods include, but are not limited to, the convex splitting method [9, 10, 18, 19], the linear stabilization method [32, 34, 43, 44], the Lagrangian multiplier method [8, 24], the Invariant Energy Quadratization (IEQ) method [45, 48], and the Scalar Auxiliary Variable (SAV) method [41, 42].

The convex splitting method was first introduced by Elliott and Stuart in [18] to numerically study a class of semilinear parabolic  equations. The main idea of the convex splitting method for gradient flow models is to split the free energy functionals into two parts: the convex part and the concave one. The convex and concave parts are then treated implicitly and explicitly, respectively, to derive unconditionally energy stable or energy dissipative numerical schemes, as seen in [10, 19, 23] and the references cited therein. Although first- and second-order unconditionally stable convex splitting schemes have been obtained for some specific energy functionals of gradient flows, there is no general convex-concave splitting form of the energy functional that allows for high-order time-stepping schemes that are unconditionally energy stable.  Another widely-used approach is the linear stabilization method, which adds one or two linear stabilization terms into the corresponding time-stepping schemes based on backward differentiation formulation (BDF), Crank-Nicolson, or exponential time differencing (ETD) methods to greatly relax the time step size constraints. In this type of method, the nonlinear term is usually treated explicitly, leading to a linear system solved at each time step. Although unconditional or conditional stability of the linear stabilization method for some gradient flows  have been successfully established (see e.g., [14, 34, 35, 36]), there does not exist a general framework for stability analysis of these existing schemes, particularly for high-order ones.
Recently, the SAV approach [42] and its variants [5, 27, 31] have been developed to design unconditionally energy stable linear schemes for gradient flow problems. In particular, SAV schemes with high-order accuracy and energy stability have been of great interest, due to the long-time nature of the gradient flow coarsening process. In [5], an arbitrarily high-order extrapolated and linearized Runge–Kutta SAV (RK–SAV) method was constructed for the Allen–Cahn and Cahn–Hilliard  equations. The proposed RK-SAV method is unconditionally stable with respect to a modified discrete energy, and the corresponding error estimate was also rigorously derived. Huang et al. [30] introduced a novel SAV approach to construct an implicit-explicit linear and unconditionally energy stable BDF$k$ $(1 \leq k \leq 5)$ method for gradient flows. However, its unconditional stability only indicates dissipation in time of the auxiliary variable without information on the phase variables.

The steady-state spectral renormalization method was introduced by Ablowitz and Musslimani for the first time in [4] to compute self-localized states of nonlinear waveguides. Since then, it has been utilized in a variety of contexts, such as nonlinear optics [1, 2], Bose-Einstein condensation [26], and water waves [3]. Built on this idea, a type of time-dependent spectral renormalization (TDSR) approach was developed in [12, 13] to incorporate intrinsic physics in the form of conservation laws or dissipation rate equations in the development of numerical schemes. This approach has achieved the impressive capacity for accurately and efficiently capturing certain relevant physical properties regardless of the accuracy of the solution in [13], particularly preserving certain conservation laws of the model. However, it also encounters convergence issues when solving the nonlinear systems of the  physical evolution equation and the  ordinary differential equation involving the renormalization factors with some iterative methods. This leads to the efficiency problem of time step sizes not being large in order to ensure the convergence of the iteration. Such convergence issue often becomes particularly severe for cases where the dissipative properties (such as the energy) of the model are enforced in the TDSR approach. Furthermore, for the dissipative model, the $p$-th order TDSR scheme is numerically stable only in the sense that the error of the computed energy is of order $p$ in time, which doesn't indicate the dissipation property of the model at the discrete level.

A common feature of the gradient flow problem (<ref>) is that its evolution process often takes a considerably long time to reach the steady state. Moreover,  it usually undergoes both fast and slow energy-changing stages throughout the evolution process, which imposes the need for adaptive time stepping during the simulation. Therefore, there is a high demand for the development of high-order, structure-preserving, efficient numerical schemes with variable time steps for the gradient flow (<ref>).  To address these issues, we propose in
this paper a high-order accurate and energy dissipative method, called TDSR-ETD, for solving general gradient flow problems, which combines the spectral renormalization method for handling the energy dissipation law and the exponential integrator methods for accurate and stable time integration. As motivated by the idea of [12, 13], we further apply the proposed  TDSR-ETD  method to a broader class of phase-field models with energy dissipation properties.

The rest of the paper is organized as follows. Section 2 provides a detailed illustration of the proposed TDSR-ETD method for the $L^2$ (Allen-Cahn type) and $H^{-1}$  (Cahn-Hilliard type) gradient flows with respect to a classic free energy functional, under the periodic or homogeneous Neumann boundary conditions. This section includes discussions of the corresponding numerical schemes, solution algorithms, and physical properties, such as energy dissipation and mass conservation. In Section 3, we further explore the application of our TDSR-ETD method to other two types of gradient flow problems, the molecular beam epitaxial model and the phase-field crystal model. Section 4 presents extensive numerical experiments and comparison tests to demonstrate the accuracy and efficiency of the proposed method. Finally, some concluding remarks are given in Section 5.

§ THE SPECTRAL RENORMALIZATION EXPONENTIAL INTEGRATOR METHOD

To illustrate the proposed spectral renormalization exponential integrator method, we  take the gradient flow model (<ref>) with respect to the following classic free energy functional:

where  the parameter $\varepsilon>0$ is related to the interfacial width and  $F: \mathbb{R}\rightarrow  \mathbb{R}$ denotes a nonlinear potential function. Two types of boundary conditions usually will be considered for the above gradient flows, the periodic boundary condition and the homogenous Neumann boundary condition, respectively.
Taking the Sobolev space $H$ to be
$L^2(\Omega)$ or $H^{-1}(\Omega)$ in (<ref>)  gives us the following time-dependent PDEs:

=𝒢_H μ,              ∈Ω, t> 0,

μ=-ε^2Δϕ+f(ϕ),    ∈Ω, t> 0,

with the initial value  $\phi(\x,0)=\phi_{0}(\x)$ for any $\x\in\overline\Omega$, where $f=F'$ and

-I,     ,

We call the above equations (<ref>)  the $L^2(\Omega)$ or  $H^{-1}(\Omega)$  gradient flows with respect to the energy functional (<ref>), respectively (equivalently, the Allen-Cahn type  equation [6]  or the Cahn-Hilliard type equation [11], respectively).

Inspired by the work of [12, 13],
the key idea of the proposed spectral renormalization exponential integrator method is to introduce an extra scalar variable, that is the so-called TDSR factor $R(t)$, to incorporate the energy dissipation law (<ref>) into  the gradient flow problem (<ref>). In particular, we rewrite the solution of the gradient flow problem (<ref>)

\begin{equation}\label{ttt}
\phi(\x,t)=R(t)\psi(\x,t),
\end{equation}

and thus $\mu\big(R\psi\big)=-\varepsilon^{2}\Delta(R\psi)+f(R\psi)$.
Without loss of generality, we focus our following discussion  on
the two-dimensional case $(d=2,\;\x=(x,y))$ and assume $\Omega:=(0,L)^{2}$, but all results derived below can be straightforwardly  extended to the three-dimensional case. Let $T>0$ be a given terminal time and $\{\Dt_{n}=t_{n}-t_{n-1}>0\}_{n=1}^K$ be a general partition of the time interval $[0,T]$ such that $t_0=0$ and $\sum_{n=1}^{K} \Dt_n=T$.

§.§ The $L^2$ gradient flow with the periodic boundary condition

Based on the $L^2$ gradient flow  equation (<ref>) (i.e., $\mathcal{G}_H = -I$) and the energy dissipation law (<ref>), we have the following coupled system for

\begin{align}
&\dps\frac{\partial (R\psi)}{\partial t}=\varepsilon^{2}\Delta(R\psi)-f(R\psi),\qquad (\x,t)\in \Omega\times(0,T],\label{111}\\
&\dps\frac{dE\big[R\psi\big]}{dt}=-\|\mu\big(R\psi\big)\|^{2},\qquad\quad \quad t\in(0,T],\label{prob1-b}
\end{align}

with the initial value  $\dps R(0)=R_{0}$ and $\psi(\x,0)=\phi_{0}(\x)/R_{0}$ for any $\x\in\overline\Omega$.
It is easy to see that if $R_0 = 1$, then the  coupled system (<ref>) has a unique solution of $R(t) \equiv 1$ and $\psi(\x,t) = \phi(\x,t)$, where $\phi(\x,t)$ is the solution of the $L^2$ gradient flow equation (<ref>). Thus, we always set the initial condition $R_{0}=1$ in what flows. For numerical stabilization, a linear splitting  [15] is often applied to the gradient flow equation (<ref>) so that we deal with a transformed equation as

\begin{equation}\label{prob1a}
\dps\frac{\partial (R\psi)}{\partial t}=(\varepsilon^{2}\Delta-s)(R\psi)-(f(R\psi)-sR\psi),\
\end{equation}

where $s\geq0$ is a constant stabilizing parameter.

Fourier spectral discretization in space
Let us apply the Fourier spectral method for the spatial discretization of the system (<ref>).
Note that other space discretization methods can also be used, such as finite difference, finite element or finite volume methods.
The trial function space of the Fourier spectral method is defined as:
X_N(Ω):=span { e^i2(kx+ly)π/L, -N≤k,l≤N},
where $N$ is a positive  integer and $i=\sqrt{-1}$. Let $\Pi_{FS}$ denote the projection operator from $L^{2}(\Omega)$ to $X_{N}(\Omega)$ by

$$\Pi_{FS}v:=\dps\sum_{k,l=-N}^{N}(\widehat{\Pi}_{FS}v)_{k,l}e^{i2(kx+ly)\pi/L},\quad \forall\, v\in L^{2}(\Omega)$$

$$(\widehat{\Pi}_{FS}v)_{k,l}=\frac{(v,e^{i2(kx+ly)\pi/L})}{\|e^{i2(kx+ly)\pi/L}\|^{2}},\quad -N\leq k,l\leq N.$$

Then, we obtain the semi-discretized (in space) system of  (<ref>) as: given $\dps R(0)=1$ and $\psi_{N}(0)=\Pi_{FS}\phi_{0},$
find $R(t)$ and $\psi_{N}(t)\in X_{N}(\Omega)$ for any  $t\in(0,T]$ such that

\begin{align}
&\dps\frac{\partial(R\psi_{N})}{d t}=(\varepsilon^{2}\Delta-s)(R\psi_{N})-\Pi_{FS}\mathcal{N}(R\psi_N),\\
\end{align}

where $\mathcal{N}(R\psi_N)=f(R\psi_N)-sR\psi_N.$
Correspondingly,  $\phi_{N}(\x,t):=R(t)\psi_{N}(\x,t)$ is a semi-discrete approximation to the solution of the $L^{2}$ gradient flow equation (<ref>).
For any $\psi_{N}(t)\in X_{N}(\Omega),$ we can express it as

where $\widehat{\Psi}(t) = (\widehat{\Psi}_{k,l}(t))$ is time-dependent coefficient matrix of dimension $(2N+1)\times(2N+1)$.
Substituting the above expression into the semi-discrete system (<ref>) and using the Duhamel's principle for  $t\in(t_{n},t_{n+1})$ ($n\geq 0$ and $t_0=0$), we obtain

R(t_n+1)Ψ_k,l(t_n+1)=e^_n+1ℒ_k,lR(t_n)Ψ_k,l(t_n)-(I^Π_FS𝒩_1)_k,l,  -N≤k,l≤N,

ℒ_k,l    =-ε^2[(2kπ/L)^2+(2lπ/L)^2]-s,

(I^Π_FS𝒩_1)_k,l    =∫_t_n^t_n+1e^(t_n+1-τ)ℒ_k,l(Π_FS𝒩(R(τ)ψ_N(τ)))_k,ldτfor $-N\leq k,l\leq N$, and

ETD multistep approximation in time
For any function $u(t)$ defined on $[0,T]$,  denote $P_{r,n}u(t)$ as its Lagrange interpolation polynomial of degree $r$ using the values of $u(t)$ at $t_{n+1},t_{n},\cdots,t_{n+1-r},$
then we have
P_r,nu(t)=∑_j=-1^r-1ω_r,j(t)u(t_n-j),  t∈[t_n,t_n+1]
with $\omega_{r,j}(t)= \prod_{l=-1,l\neq j}^{r-1}\frac{t-t_{n-l}}{t_{n-j}-t_{n-l}}.$ As below, we list some of the polynomial $P_{r,n}u(t)$ up to $r=2$:

P_1,nu(t)=ηu(t_n+1)+(1-η) u(t_n),

P_2,nu(t)=(γ_n+1η+1)η/1+γ_n+1 u(t_n+1)+(1-η)(1+γ_n+1η) u(t_n))+γ^2_n+1(η-1)η/1+γ_n+1 u(t_n-1),
where $\eta=(t-t_{n})/\Dt_{n+1}$ and $\gamma_{n+1}=\Dt_{n+1}/\Dt_{n}$. Then, the ETD Multistep approximations (in the spirit of Adam-Moulton) of  the integrations $(I^{u}_{1})_{k,l}$ and $I^{u}_{2}$ are given respectively by
(I^u_1)_k,l   ≈∑_j=-1^r-1u_k,l(t_n-j)∫_t_n^t_n+1e^(t_n+1-τ)ℒ_k,lω_r,j(τ)dτ

:=∑_j=-1^r-1α_k,l^(r,j)u_k,l(t_n-j),  -N≤k,l≤N,

I^u_2   ≈∑_j=-1^r-1u(t_n-j)∫_t_n^t_n+1ω_r,j(τ)dτ:=∑_j=-1^r-1β^(r,j)u(t_n-j).
Below  we give some values of $\{\alpha^{(r,j)}_{k,l}\}_{k,l=-N}^N$ and $\{\beta^{(r,j)}_{k,l}\}_{k,l=-N}^N$ for $r=0,1,2$:

α^1,-1_k,l=α_1, α^1,0_k,l=α_0-α_1;

α^2,-1_k,l=α_1+γ_n+1α_2/1+γ_n+1,  α^2,0_k,l=α_0-α_1+γ_n+1(α_1-α_2),  α^2,1_k,l=-γ_n+1^2/1+γ_n+1(α_1-α_2)

α_0=-1/ℒ_k,l(1-e^ℒ_k,l_n+1),   α_1=-_n+1-α_0/ℒ_k,l_n+1,  α_2=-_n+1-2α_1/ℒ_k,l_n+1,     ℒ_k,l≠0,

α_0=_n+1,   α_1=_n+1/2,  α_2=_n+1/3,

β^2,-1=_n+1(2γ_n+1+3)/6(1+γ_n+1),  β^2,0=3+γ_n+1/6_n+1,

Then, the fully discrete TDSR-ETD scheme with $(r+1)$-th order in time for the continuous  system (<ref>) (and the semi-discrete problem (<ref>))  reads: for $n=r,r+1,\cdots, K-1$, find $R^{n+1}$ and $\widehat{\Psi}^{n+1}$ such that

\begin{align}
&\dps R^{n+1}\widehat{\Psi}^{n+1}_{k,l}=e^{\Dt_{n+1}\mathcal{L}_{k,l}}R^{n}\widehat{\Psi}^{n}_{k,l}-\Big(\mathcal{I}^{\widehat{\Pi}_{FS}\mathcal{N},r}_{1,n+1}\Big)_{k,l},\quad -N\leq k,l\leq N,\label{full-ds1}\\
&\dps E\big[R^{n+1}\psi^{n+1}_{N}\big]-E\big[R^{n}\psi^{n}_{N}\big]=-\max\{\mathcal{I}^{\mu,r}_{2,n+1},0\},\label{full-ds2}
\end{align}

\begin{align}
&\Big(\mathcal{I}^{\widehat{\Pi}_{FS}\mathcal{N},r}_{1,n+1}\Big)_{k,l}\dps \;=\sum_{j=-1}^{r-1}\alpha_{k,l}^{(r,j)}\Big(\widehat{\Pi}_{FS}\mathcal{N}\big(R^{n-j}\psi^{n-j}_{N}\big)\Big)_{k,l},\label{int_num1}\\
\end{align}

Note that $\mathcal{I}^{\mu,r}_{2,n+1}$ approximates $I^{\mu}_{2}$ with $(r+1)$-th order in time and $I^{\mu}_{2}\ge 0$, therefore

$$|\max\{\mathcal{I}^{\mu,r}_{2,n+1},0\}-I^{\mu}_{2}|\leq  |\mathcal{I}^{\mu,r}_{2,n+1}-I^{\mu}_{2}|,$$

and consequently $\max\{\mathcal{I}^{\mu,r}_{2,n+1},0\}$ is also an  $(r+1)$-th order approximation of $I^{\mu}_{2}$. Furthermore, when $r=1$ and $2$,
it is easy to verify that $\beta^{r,j}>0$ for $j=-1,\cdots,r-1$ and thus $\max\{\mathcal{I}^{\mu,r}_{2,n+1},0\}=\mathcal{I}^{\mu,r}_{2,n+1}$.

To start  the $(r+1)$-th order TDSR-ETD scheme (<ref>),
the needed values of $\{(R^{k},\widehat{\Psi}^{k})\}_{k=1}^r$  can be computed
using  the TDSR-ETD schemes with $r$-th order or some other methods.

The proposed TDSR-ETD scheme (<ref>) is a nonlinear scheme and has to be solved with some iterative algorithms at each time step. When a physically valid solution of  (<ref>) exists, due the nonnegativity of $\max\{\mathcal{I}^{\mu,r}_{2,n+1},0\}$, it follows directly from (<ref>) and $\phi^{n}_{N}:=R^{n}\psi^{n}_{N}$ that
E[ϕ^n+1_N]≤E[ϕ^n_N],  n≥0
i.e., the TDSR-ETD scheme (<ref>) is  dissipative with respect to the original energy (<ref>) in the discrete sense. However,
it usually requires the time step size not be too large to ensure  the existence of the solution and the convergence of the iterative solution process for the scheme (<ref>).

Extra forcing term for the TDSR factor and Picard iteration
We know $R^{n}$ is an approximation to $1$ since $R(t)\equiv1$ and $\psi(\x,t)=\phi(\x,t)$ is the solution of the continuous  problem (<ref>).
In order to enforce that the numerical approximation for the TDSR factor $R(t)$ remains close to 1 and relax the solvability of the coupled system, we additionally incorporate  the following forcing term for $R(t)$

dR^2/dt=0,   t> 0,
into  (<ref>) to obtain a slightly modified system as follows:

\begin{align}
&\dps\frac{\partial (R\psi)}{\partial t}-\varepsilon^{2}\Delta(R\psi)+f(R\psi)=0,\qquad (\x,t)\in \Omega\times(0,T],\\
&\dps\frac{d\big(E\big[R\psi\big]+\theta R^2\big)}{dt}=-\|\mu\big(R\psi\big)\|^{2},\qquad\; t\in(0,T],\label{prob12-2}
\end{align}

where $\theta\geq0$ is a constant enforcing parameter.
Note that the above system  is again equivalent to the original problem (<ref>) and the coupled system (<ref>).
The motivation is to simultaneously enforce $R(t)$ to stay close to $1$ and (slightly) relax the energy dissipation law so that the time step size restriction for updating  $R(t)$ can be improved.
By following the similar derivations for  (<ref>), the fully discrete $(r+1)$-th order in time TDSR-ETD scheme for the revised system (<ref>) reads:  for $n=r,r+1,\cdots, K-1$, find $R^{n+1}$ and $\widehat{\Psi}^{n+1}$ such that

\begin{align}
&\dps R^{n+1}\widehat{\Psi}^{n+1}_{k,l}=e^{\Dt_{n+1}\mathcal{L}_{k,l}}R^{n}\widehat{\Psi}^{n}_{k,l}-\Big(\mathcal{I}^{\widehat{\Pi}_{FS}\mathcal{N},r}_{1,n+1}\Big)_{k,l},\label{full-ds-enf1}\\
&\dps E\big[R^{n+1}\psi^{n+1}_{N}\big]-E\big[R^{n}\psi^{n}_{N}\big]+\theta((R^{n+1})^2-(R^{n})^2)=-\max\{\mathcal{I}^{\mu,r}_{2,n+1},0\}.\label{full-ds-enf2}
\end{align}

Then, we have the following modified discrete energy dissipation law  for the TDSR-ETD scheme (<ref>)

E[ϕ^n+1_N]+θ[(R^n+1)^2-1]≤E[ϕ^n_N]+θ[(R^n)^2-1],   n≥0.
We note that the modified discrete energy $E\big[\phi^{n}_{N}\big]+\theta \big[(R^{n})^2-1\big]$ is of $(r+1)$th-order in time approximation to the original discrete energy $E\big[\phi^{n}_{N}\big]$ when $R^{n}$ approximates 1 with $(r+1)$th-order accuracy. Moreover,  $\{E\big[\phi^{n}_{N}\big]\}_{n\geq 0}$ is uniformly bounded from above  if and only if  the modified discrete energy is also bounded from above.

Finally, from $t_n$ to $t_{n+1}$, the solution $(R^{n+1},\widehat{\Psi}^{n+1})$ of (<ref>)  can be efficiently solved with the Picard iteration and the decoupling of the system is as follows:

Algorithm 1.
Set $m=0$, $R^{n+1,(0)}=R^{n}$ and $\widehat{\Psi}^{n+1,(0)} = \widehat{\Psi}^{n}$. Then

(i) compute $R^{n+1,(m+1)}$ by using Newton's method with initial guess $R^{n+1,(m)}$ to solve the scalar nonlinear system (derived from (<ref>))

\begin{align*}%\label{sol1}
&E\big[R^{n+1,(m+1)}\psi^{n+1,(m)}_{N}\big]+\theta (R^{n+1,(m+1)})^2=E\big[R^{n}\psi^{n}_{N}\big]+\theta (R^{n})^2\nonumber\\
\end{align*}

(ii) compute   $\widehat{\Psi}^{n+1,(m+1)}$ by solving the linear system (derived from (<ref>))

\begin{align*}%\label{sol2}
\end{align*}

(iii) set $m = m+1$. If not convergent, repeat Steps (i) and (ii); otherwise  return $R^{n+1}=R^{n+1,(m+1)}$ and $\widehat{\Psi}^{n+1} = \widehat{\Psi}^{n+1,(m+1)}$.

The convergence criterion is set to be $\|\psi^{n+1,(m)}_{N} - \psi^{n+1,(m+1)}_{N}\|_{\infty}\leq \hat\epsilon$ for some pre-given tolerance $\hat\epsilon>0$ in our numerical experiments.

§.§ The $H^{-1}$ gradient flow with homogeneous Neumann boundary condition

Next we turn to present the TDSR-ETD method for the $H^{-1}$ gradient flow problem (<ref>) (i.e., $\mathcal{G}_H = \Delta$).
We similarly obtain the following coupled system for

\begin{align}
&\dps\frac{\partial (R\psi)}{\partial t}=\Delta(\varepsilon^{2}\Delta(R\psi)-f(R\psi)),\qquad (\x,t)\in \Omega\times(0,T],\\
&\dps\frac{d\big(E\big[R\psi\big]+\theta R^2\big)}{dt}=-\|\nabla\mu\big(R\psi\big)\|^{2},\quad t\in(0,T],\label{prob12-b}
\end{align}

with the initial value  $\dps R(0)=1$ and $\psi(\x,0)=\phi_{0}(\x)$ for any $\x\in\overline\Omega$.
Instead of the periodic boundary condition discussed in the previous subsection, we now consider the homogeneous Neumann boundary condition

where $\n$ denotes the outward unit normal vector on the boundary $\partial\Omega$.

We will use Legendre Galerkin (spectral) discretization in space for the coupled system (<ref>) due to the homogeneous Neumann boundary condition.
Denote by $Q_{N}(\Omega)$ the space of polynomials of degree  less than or equal to $N$ with respect to each variable, and define
$Q^0_{N}(\Omega):=\{ v\in Q_{N}(\Omega): \frac{\partial v}{\partial \n}|_{\partial\Omega}=0\}$.
Let $\Pi_{LP}$ denote the usual interpolation operator from $L^{2}(\Omega)$ into $Q_{N}(\Omega)$.
Then we obtain the semi-discretized (in space) problem of  (<ref>):
given $\dps R(0)=1$ and $\psi_{N}(0)=\Pi_{LP}\phi_{0}$,
find $R(t)$, $\psi_{N}(t)\in Q^0_{N}(\Omega)$, and $\mu_{N}(t) \in Q^0_{N}(\Omega)$ for any $q_{N}\in Q^0_{N}(\Omega)$ such that

\begin{align}
&\dps\Big(\frac{\partial(R\psi_{N})}{\partial t},q_{N}\Big)=-(\nabla\mu_{N},\nabla q_{N}),\label{sem-lp-a} \\
&\dps(\mu_{N},q_{N})=\varepsilon^{2}(\nabla(R\psi_{N}),\nabla q_{N})+s(R\psi_{N},q_{N})+\big(\Pi_{LP}\mathcal{N}(R\psi_{N}),q_{N}\big),\label{sem-lp-b}\\
&\dps\frac{d\big(E\big(R\psi_N\big)+\theta R^2\big)}{dt}=-\|\nabla\mu\big(R\psi_{N}\big)\|^{2}.\label{sem-lp-c}
\end{align}

Let us express $\psi_{N}(\x,t)$ and $\mu_{N}(\x,t)$ as

ψ_N(,t)=∑_k,l=0^N-2Ψ_k,l(t)h_k(x)h_l(y),    μ_N(,t)=∑_k,l=0^N-2μ_k,l(t)h_k(x)h_l(y),

and $L_{k}(\cdot)$ denotes  the Legendre polynomial of degree $k$.
Then we deduce from (<ref>) and (<ref>) the following ordinary differential equation system:

where $M=(M_{k,l})$, $S=(S_{k,l})$ and $\widehat{\mathcal{N}}$ are  $(N-1)\times(N-1)$ matrices with the elements given by
M_k,l=(h_l,h_k),  S_k,l=(h'_l,h'_k), (𝒩)_k,l=(Π_LP𝒩(Rψ_N),h_l(x)h_k(y))
for $0\leq k,l\leq N-2.$
It is readily seen  [39] that the mass matrix $M$ is a symmetric positive definite matrix
and the stiffness matrix $S$ is diagonal with positive elements.
Therefore, the ODE system (<ref>) can be rewritten as

+ M^-1𝒩 M^-1.

Let $\Lambda$ be the diagonal matrix whose diagonal elements are  $\{\lambda_{k}>0\}_{k=0}^{N-2}$, the eigenvalues of $M^{-1}S$, and $P$ be the corresponding eigenvector matrix which is orthonormal, i.e.,

Define $\overline{\Psi}(t)=P^T\widehat{\Psi}(t)P$ and $\overline{\mu}(t)=P^T\widehat{\mu}(t)P$, then we derive from (<ref>) and  (<ref>) that

+ H,

$H(t):=(MP)^{-1}\widehat{\mathcal{N}}(t) (MP)^{-T}$.
Moreover, (<ref>) can be equivalently rewritten as

d(RΨ_k,l)/dt=-(λ_k+λ_l)μ_k,l,  0≤k,l≤N-2.

+ (H)_k,l,  0≤k,l≤N-2.

Therefore, we obtain an equivalent system for (<ref>) as follows:

\begin{align}
&\dps \frac{d\big(R\overline{\Psi}_{k,l}\big)}{dt}=-(\lambda_{k}+\lambda_{l})(\varepsilon^{2}\big(\lambda_{k}+\lambda_{l})+s)R\overline{\Psi}_{k,l}-\mathcal{H}_{k,l}, \quad 0\leq k,l\leq N-2,\\
&\dps\frac{d\big(E\big[R\psi_N\big]+\theta R^2\big)}{dt}=-\|\nabla\mu\big(R(t)\psi_{N}\big)\|^{2},
\end{align}

with $\mathcal{H}_{k,l}(t):=(\lambda_{k}+\lambda_{l})(H(t))_{k,l}$ and $\dps \psi_{N}(t)\dps=\sum_{k,l=0}^{N-2}(P\overline{\Psi}(t)P^{T})_{k,l}h_{k}(x)h_{l}(y).$

Similar to the process of deriving the fully discrete scheme (<ref>), by applying  Duhamel's principle to (<ref>), we can obtain the
fully discrete TDSR-ETD scheme with $r+1$ order in time for the $H^{-1}$ gradient flow (<ref>) as follows:

for $n = r, r + 1,\cdots, K-1$, find $R^{n+1}$ and $\overline{\Psi}^{n+1}$ such that

\begin{align}
&\dps R^{n+1}\overline{\Psi}^{n+1}_{k,l}=e^{\Dt_{n+1}\mathcal{L}_{k,l}}R^{n}\overline{\Psi}^{n}_{k,l}-\Big(\mathcal{I}^{\mathcal{H},r}_{1,n+1}\Big)_{k,l},\label{full-Ne1}\\
% &\dps \psi_{N}^{n+1}\dps=\sum_{k,l=0}^{N-2}(P\overline{\Psi}^{n+1}P^{T})_{k,l}h_{k}(x)h_{l}(y),\label{full-Ne2}\\
&\dps E\big[R^{n+1}\psi^{n+1}_{N}\big]-E\big[R^{n}\psi^{n}_{N}\big]+\theta((R^{n+1})^2-(R^{n})^2)=-\max\{\mathcal{I}^{\nabla\mu,r}_{2,n+1},0\}, \label{full-Ne3}
\end{align}

with $\psi_{N}^{n+1}=\sum_{k,l=0}^{N-2}(P\overline{\Psi}^{n+1}P^{T})_{k,l}h_{k}(x)h_{l}(y)$, where

\begin{equation*}%\label{int_num2}
\mathcal{L}_{k,l}=-(\lambda_{k}+\lambda_{l})[\varepsilon^{2}\big(\lambda_{k}+\lambda_{l})+s],\quad
\Big(\mathcal{I}^{\mathcal{H},r}_{1,n+1}\Big)_{k,l}\dps=\sum_{j=-1}^{r-1}\alpha_{k,l}^{(r,j)}\big(\mathcal{H}^{n-j}\big)_{k,l},\\
\end{equation*}

for $0\leq k,l\leq N-2$, and
From $t_n$ to $t_{n+1}$, the solution $(R^{n+1},\overline{\Psi}^{n+1})$ of (<ref>)  can  be efficiently solved using Algorithm 1 in the same way.

The following modified discrete energy dissipation law is again deduced  for the TDSR-ETD scheme (<ref>)

\begin{equation}
E\big[\phi^{n+1}_{N}\big]+\theta \big[(R^{n+1})^2-1\big]\leq E\big[\phi^{n}_{N}\big]+\theta\big[ (R^{n})^2-1\big], \quad n\geq 0,
\end{equation}

with $\phi^{n}_{N}=R^{n}\psi^{n}_{N}$.
In addition, the $H^{-1}$ gradient flow (<ref>) is mass-conserved along the time (i.e., $\frac{d}{dt}\int_{\Omega}\phi d\x=0$),
it is easy to check the  TDSR-ETD scheme (<ref>) conserves the mass in the discrete sense, that is

\begin{equation}
\int_{\Omega}\phi^{n}_{N}(\x)d\x = \int_{\Omega}\phi^{0}_{N}(\x)d\x, \quad\forall\,n\geq 0.
\end{equation}

§ APPLICATION TO OTHER TYPES OF GRADIENT FLOWS

In addition to the above two classic gradient flows (Allen-Cahn and Cahn-Hilliard equations) with respect to (<ref>), the proposed TDSR-ETD method can be naturally applied to many other types of gradient flow problems. In this section, we  will discuss its application to  two challenging phase field models, the molecular beam epitaxial (MBE) model [20, 33] and the phase-field crystal (PFC) model [16, 17].

§.§ The MBE model

In this subsection, we aim to construct efficient energy dissipative numerical schemes  for the MBE model equipped with periodic boundary conditions via the TDSR-ETD method.
The MBE model can be viewed as the $L^{2}$-gradient flow with respect to the free energy functional defined by

where the two commonly used nonlinear potential functionals $F_M: \mathbb{R}^{d}\rightarrow \mathbb{R}$ are given by

Then, the MBE model reads

∂ϕ/∂t=-(ε^2Δ^2ϕ+f_M(∇ϕ)),   (,t)∈Ω×(0,T],
with initial condition $\phi(\x,0)=\phi_{0}(x)$ for any $\x\in\overline{\Omega}$, where $f_M(\nabla\phi)={\delta F_M(\nabla\phi)}/{\delta\phi}$.
It satisfies the following energy dissipation law:
Setting $\phi(\x,t)=R(t)\psi(\x,t)$ as before, and then we can rewrite the MBE model (<ref>) into the following
equivalent coupled system for $(R(t),\psi(\x,t))$:

\begin{align}
&\dps\frac{\partial (R\psi)}{\partial t}=-\varepsilon^{2}\Delta^{2}(R\psi)+s\Delta(R\psi)-\mathcal{N}_{M}(R\nabla\psi), \quad(\x,t)\in \Omega\times(0,T],\\
&\dps\frac{d\big(E_M\big[R\psi\big]+\theta R^2\big)}{dt}=-\|\mu_{M}\big(R\psi\big)\|^{2},\qquad\quad\qquad\quad t\in(0,T],
\dps
\end{align}

with $R(0)=1$ and $\psi(\x,0)=\phi_{0}(\x)$ for any $\x\in\overline{\Omega}$, where $\mathcal{N}_{M}(\cdot)$ and $\mu_{M}(\cdot)$ are given respectively by

\begin{equation*}
\mathcal{N}_{M}(\nabla(R\psi))=f_M(\nabla(R\psi))+s\nabla\cdot\nabla(R\psi),\quad
\mu_{M}\big(R\psi\big)=\varepsilon^{2}\Delta^{2}(R\psi)+f_M(\nabla(R\psi)).
% \end{eqnarray*}
\end{equation*}

Let us apply the Fourier spectral method for the spatial discretization of the system (<ref>) as done in subsection (<ref>)
and set $\psi_{N}(0)=\Pi_{FS}\phi_{0}$. For  any $\psi_N\in X_N$ there is a representation (<ref>).
By following a similar process discussed in the previous section, we obtain the fully-discrete $(r+1)$-th order in time TDSR-ETD scheme for the MBE model (<ref>) as follows:

for $n = r, r + 1, \cdots, K-1$, find $R^{n+1}$ and ${\widehat\Psi}^{n+1}$ such that

\begin{align}
&\dps R^{n+1}\widehat{\Psi}^{n+1}_{k,l}=e^{\Dt_{n+1}(\mathcal{L}_{M})_{k,l}}R^{n}\widehat{\Psi}^{n}_{k,l}-\Big(\mathcal{I}^{\widehat\Pi_{FS}\mathcal{N}_{M},r}_{1,n+1}\Big)_{k,l},\label{full-ds-mbe1}\\
&\dps E_M\big[R^{n+1}\psi^{n+1}_{N}\big]-E_M\big[\phi^{n}_{N}\big]+\theta((R^{n+1})^2-(R^{n})^2)=-\max\{\mathcal{I}^{\mu_{M},r}_{2,n+1},0\},\label{full-ds-mbe2}
\end{align}

(ℒ_M)_k,l    =-ε^2[(2kπ/L)^2+(2lπ/L)^2]^2-s[(2kπ/L)^2+(2lπ/L)^2],

(ℐ^Π_FS𝒩_M,r_1,n+1)_k,l    =∑_j=-1^r-1α_k,l^(r,j)(Π_FS𝒩_M(∇(R^n-jψ^n-j_N)))_k,l,
for $-N\leq k,l\leq N$, and
\mathcal{I}^{\mu_M,r}_{2,n+1}\dps\;=\sum_{j=-1}^{r-1}\beta^{(r,j)}\|\mu_M\big(R^{n-j}\psi_{N}^{n-j}\big)\|^{2}.
Obviously, the TDSR-ETD scheme (<ref>) is energy dissipative in the following sense:

\begin{equation}
E_M\big[\phi^{n+1}_{N}\big]+\theta \big[(R^{n+1})^2-1\big]\leq E_M\big[\phi^{n}_{N}\big]+\theta \big[(R^{n})^2-1\big], \quad n\geq 0,
\end{equation}

where $\phi^{n}_{N} = R^n\psi^n_N$.

§.§ The PFC model

The  PFC model was first proposed in [16, 17] and has been frequently used for modelling crystal growth at the atomic scale in both space and diffusive time scales.
It allows for the nucleation of crystallites at arbitrary locations and orientations while containing elastic and plastic deformations.
In this subsection, we consider the following phase field crystal equation with the periodic boundary condition:

=Δμ_p(ϕ),                   (,t)∈Ω×(0,T],

μ_p(ϕ)=(Δ+σ)^2ϕ+ϕ^3-δϕ,  (,t)∈Ω×(0,T],

with initial condition $\phi(\x,0)=\phi_{0}(\x),$ where $\sigma$ and $\delta$ are two positive constants such that $0<\delta<\sigma^{2}$.
The above PFC model is mass conservative. Moreover, it satisfies the following energy dissipation law:
d/dt E_p[ϕ]=-∇μ_p(ϕ)^2≤0,
Setting $\phi(\x,t)=R(t)\psi(\x,t)$ as before, and then we can rewrite the PFC model (<ref>) into the following
equivalent coupled system for $(R(t),\psi(\x,t))$:

\begin{align}
&\dps\frac{\partial (R\psi)}{\partial t}
=\Delta(\Delta+\sigma)^{2}(R\psi)+s\Delta(R\psi)-\mathcal{N}_{p}(R\psi),\qquad (\x,t)\in \Omega\times(0,T],\\\
\frac{d}{dt} \big(E_p[R\psi]+\theta R^2\big)=-\|\nabla\mu_{p}(R\psi)\|^{2},\qquad\qquad\quad\qquad\quad t\in(0,T],
\dps
\end{align}

with $R(0)=1, \psi(\x,0)=\phi_{0}(\x)$ for any $\x\in\overline{\Omega}$,

Let us again apply the Fourier spectral method for the spatial discretization of (<ref>).
Then we obtain the $(r+1)$-th order in time TDSR-ETD scheme for the PFC model (<ref>) as follows: for $n = r, r + 1, \cdots, K-1$, find $R^{n+1}$ and ${\widehat\Psi}^{n+1}$ such that

\begin{align}
&\dps R^{n+1}\widehat{\Psi}^{n+1}_{k,l}=e^{\Dt_{n+1}\big(\mathcal{L}_{p}\big)_{k,l}}R^{n}\widehat{\Psi}^{n}_{k,l}-\Big(\mathcal{I}^{\widehat\Pi_{FS}\mathcal{N}_{p},r}_{1,n+1}\Big)_{k,l},\label{full-ds-pfc1}\\
&\dps E_p\big[R^{n+1}\psi^{n+1}_{N}\big]-E_p\big[\phi^{n}_{N}\big]+\theta((R^{n+1})^2-(R^{n})^2)=-\max\{\mathcal{I}^{\mu_{p},r}_{2,n+1},0\},\label{full-ds-pfc2}
\end{align}

(ℒ_p)_k,l    =-((2kπ/L)^2+(2lπ/L)^2)[(-(2kπ/L)^2-(2lπ/L)^2+σ)^2+s],

(ℐ^Π_FS𝒩_p,r_1,n+1)_k,l    =∑_j=-1^r-1α_k,l^(r,j)(Π_FS𝒩_p(R^n-jψ^n-j_N))_k,l,
for $-N\leq k,l\leq N$, and
\mathcal{I}^{\mu_M,r}_{2,n+1}\dps:=\sum_{j=-1}^{r-1}\beta^{(r,j)}\|\mu_p\big(R^{n-j}\psi_{N}^{n-j}\big)\|^{2}.
It is easy to check that the TDSR-ETD scheme (<ref>) is mass conservative and energy dissipative in the following sense:
E_p[ϕ^n+1_N]+θ[(R^n+1)^2-1]≤E_p[ϕ^n_N]+θ[(R^n)^2-1],  n≥0,
where $\phi^{n}_{N} = R^n\psi^n_N$.

§ NUMERICAL EXPERIMENTS

In this section,  various numerical experiments and comparisons are presented to validate the theoretical results derived for the proposed TDSR-ETD schemes. In the following experiments, the tolerance error for the Picard iteration in Algorithm 1 at each time step is always set to be $\hat\epsilon = 10^{-7}$ unless specified otherwise.
All the errors are computed under the $L^{\infty}$-norm measure throughout the numerical tests. Specifically, we test the TDSR-ETD schemes (<ref>)  with $r=1$ and $r=2$, which are second-order and third-order accurate in time and thus referred to as TDSR-ETD2 and  TDSR-ETD3 respectively.

§.§ Convergence tests

Consider the gradient flow equation (<ref>) with the double-well potential $F(\phi) = \frac14(\phi^2-1)^2$:

∈Ω, t>0,

ϕ(,0)=ϕ_0(),                       ∈Ω,

where $\mathcal{G}_{H}=-I$ (the $L^2$ gradient flow, i.e., Allen-Cahn equation) or $\mathcal{G}_{H}=\Delta$ (the $H^{-1}$ gradient flow,  i.e., Cahn-Hilliard equation) as defined in (<ref>), and

The periodic  and the homogeneous Neumann boundary conditions are imposed for the above $L^2$ and $H^{-1}$ gradient flows, respectively. We set the interfacial width parameter $\varepsilon^2=0.01$ with the computational domain $\Omega=(0,2\pi)^2$ for the $L^2$ gradient flow, and $\varepsilon^2=2.5\times 10^{-3}$ with $\Omega=(-1,1)^2$ for the $H^{-1}$ gradient flow.

(a) $\phi$

(b)  $R$

(Example <ref>) Plots of the numerical errors of $\phi$  and $R$  with respect  to  the time step size $\tau$ produced by the TDSR-ETD2 and TDSR-ETD3 schemes (<ref>)  for the $L^2$ gradient flow. For the purpose of comparison, the corresponding result of $\phi$ produced by the SAV-BDF2 scheme is added in the left panel.

We employ the proposed TDSR-ETD schemes (<ref>) and (<ref>) respectively, to compute the numerical solutions of (<ref>) for these two types of gradient flows. The Fourier spectral method with $128\times128$ basis modes is applied in the spatial discretization of the $L^2$ gradient flow, and the spatial discretization used  for the $H^{-1}$ gradient flow is the Legendre spectral method with $256\times256$ basis modes. The uniform time step is used and we set the stabilizing parameter $s=2$ as done in [15] for this test.
It holds $R(t)\equiv1$ for the continuous problem but there are no exact solutions of $\phi$ or $\psi$ available for the tested two gradient flows. Thus we take the approximate solutions produced by the TDSR-ETD3 scheme with the time step size $\Dt=10^{-5}$ as the reference solution for computing  numerical solution errors on $\phi$ and $R$ at $T=1$.

For the $L^{2}$ gradient flow ($\mathcal{G}_{H}=-I$, Allen-Cahn equation), we choose a set of uniformly refined time step sizes $\tau=0.1\times 2^{-k}, k=0,1,\cdots,4$.  We set $\hat\epsilon = 10^{-10}$ for the Picard iteration in this test to catch $R$ very accurately to produce the expected order of convergence for $R$, especially for the high-order scheme TDSR-ETD3. Figure <ref> presents the numerical errors of $\phi$  and $R$ computed by the TDSR-ETD2 and  TDSR-ETD3 schemes as functions of the time step size in log-log scale, in which two choices of  the forcing parameter are tested, $\theta=0$ and $\theta=10$ respectively. Moreover, for comparison, we also display in Figure <ref>-(a) the error behavior of $\phi$  by the SAV-BDF2 scheme.
When  $\theta=0$, we find that  TDSR-ETD2 always fails due to non-convergence of Algorithm 1 and TDSR-ETD3 works only for the time step size $\tau= 1/160$. When $\theta=10$, both TDSR-ETD2 and TDSR-ETD3 work well for all tested time step sizes.
This implies that the introduction of the enforcing parameter $\theta$ does greatly improve the time step size restriction for the TDSR-ETD schemes.
It is observed as expected  that the desired second-order accuracy in time of TDSR-ETD2  and third-order accuracy  of  TDSR-ETD3 are achieved for  numerical solutions of both $\phi$ and $R$, and the numerical errors of TDSR-ETD3 are among 10 to 100 times smaller than those of TDSR-ETD2, which demonstrate the advantage of higher-order accuracy. In addition, we see  from Figure <ref> that the numerical errors of TDSR-ETD2 always remain about 10 times smaller than those of SAV-BDF2 under the same time step size, which indicates TDSR-ETD2 performs much better than SAV-BDF2 even though they are both second-order accurate in time.

(a) TDSR-ETD2 with $\theta=10, 10^{2},10^{3}$ and $10^{4}$ and SAV-BDF2. Left: $\phi$, right: R.

(b) TDSR-ETD3 with $\theta=10, 10^{2},10^{3}$ and $10^{4}$. Left: $\phi$, right: R.
(Example <ref>) Plots of the numerical errors of  $\phi$  and $R$ with respect to the time step size $\tau$ produced by the TDSR-ETD2 and TDSR-ETD3 schemes (<ref>)   for the $H^{-1}$ gradient flow. For the purpose of comparison, the corresponding result of $\phi$ produced by the SAV-BDF2 scheme is added in the left panel of (a).

(a) $\tau=5\times 10^{-3}$

(b) $\tau=10^{-3}$

(Example <ref>) Evolution  of the number of Picard iterations of Algorithm 1 at each time step (up to the time $T=5$) for the TDSR-ETD3 scheme (<ref>)  with uniform time stepping  for the $H^{-1}$ gradient flow.

Next, we investigate the effect of the enforcing parameter $\theta$ on the convergence behaviors of $\phi$ and $R$, and on the needed numbers of  Picard iterations for convergence  in Algorithm 1 for the TDSR-ETD2 and TDSR-ETD3 schemes by testing the $H^{-1}$ gradient flow ($\mathcal{G}_{H}=\Delta$, Cahn-Hilliard equation). We choose several different enforcing parameters as $\theta=10, 10^{2}, 10^{3},$ and $10^{4}$ and a set of uniformly refined time step sizes $\tau=0.01\times 2^{-k}, k=0,1,\cdots,6$. Figure <ref> presents the error behaviors of the two TDSR-ETD schemes with respect to the time step sizes. Moreover, we also plot in Figure <ref>-(a)  the numerical errors  by the SAV-BDF2 scheme to make a comparison. It is shown as expected that TDSR-ETD2 and TDSR-ETD3 are of second-order and third-order accuracy in time for both $\phi$ and $R$, respectively, for all tested cases. Moreover, TDSR-ETD2 again achieves about 10 times better accuracy than SAV-BDF2 under the same time step size. In addition, we observe  from Figure <ref> that the accuracy of the numerical solutions for $\phi$ isn't  affected by the choice of these four values for $\theta$ but larger $\theta$ forces the numerical solution of $R$ to be closer to the exact solution $1$.
The evolution in time of the number of Picard iterations of Algorithm 1 at each time step (up to a longer  time $T=5$ for better illustration) is plotted in Figure <ref> for the TDSR-ETD3 scheme with two time step sizes $\tau=5\times 10^{-3}$ and $\tau=10^{-3}$. It shows that a larger forcing parameter $\theta$ can  reduce the number of Picard iterations and is thus more efficient, particularly for the case of larger time step size such as $\tau=5\times 10^{-3}$ used in the simulation.
Therefore, we  will by default set $\theta=10^{4}$ in all remaining numerical experiments.

Consider the following MBE model with the periodic boundary condition:

∂ϕ/∂t+ε^2Δ^2 ϕ+f_M(∇ϕ)=0,                  ∈(0,2π)^2, t>0,

ϕ(,0)=0.1(sin3xsin2y+sin5xsin5y),       ∈(0,2π)^2,

where $\varepsilon^{2}=0.01$ and $f_M(\nabla\phi)=\delta F_M(\nabla\phi)/\delta\phi$ with  $F_M$ defined in (<ref>).

We adopt the TDSR-ETD scheme (<ref>) with $128\times128$ Fourier modes and uniform time steps, and set the stabilizing parameter $s=2$ and $s=1/8$ for the model with slope selection and without slope selection [32],  respectively.
Moreover, the numerical solution computed by the TDSR-ETD3 scheme with the  small  time step $\tau=10^{-5}$ is used as a reference solution at $T=1$ for calculating the numerical solution errors.
Figure <ref> plots the numerical errors of $\phi$ with respect to the time step size $\tau$ in the log-log scale for the MBE model with and without slope selection, respectively. The results are again in good agreement with the expected convergence rate in time, i.e., second order  for  TDSR-ETD2 and third order for TDSR-ETD3.

(a) With slope selection

(b)  Without slope selection

(Example <ref>) Plots of the numerical errors of $\phi$  with respect  to  the time step size $\tau$ produced by the TDSR-ETD2 and TDSR-ETD3 schemes  (<ref>) for the MBE model.

Consider following PFC model (<ref>) with the periodic boundary condition and $\sigma=1$ and $\delta=0.025$:

=Δ((Δ+σ)^2ϕ+ϕ^3-δϕ),        ∈(0,32)^2, t>0,

ϕ(,0)=sin(πx/16)cos(πy/16),         ∈(0,32)^2.

We adopt  the TDSR-ETD scheme (<ref>) with  $256\times256$ Fourier modes and uniform time stepping. We set the stabilizing parameter  $s=\delta$.
To validate the temporal accuracy, a reference solution of the PFC model (<ref>) at $T=1$ is computed by the TDSR-ETD3 scheme with a small enough time step size $\Dt=10^{-5}$.
In Figure <ref>, the numerical error as a function of the time step size is plotted in the log-log scale.
It is observed that TDSR-ETD2 and TDSR-ETD3  again achieve the expected second order and third order accuracy in time, respectively.

(Example <ref>) Plots of the numerical errors of  $\phi$  and $R$ with respect to the time step size $\tau$ produced by the TDSR-ETD2 and TDSR-ETD3 schemes (<ref>) for the PFC model.

§.§ Combined with a time-adaptive strategy

One of the main advantages of the energy stable or energy dissipative numerical schemes is that they can be easily and naturally applied in a time-adaptive strategy.
In this subsection, we investigate the efficiency of the proposed TDSR-ETD schemes combined with a time-adaptive strategy.
There already exist several efficient time-adaptive strategies to the energy stable schemes,
including the ones [21, 40] according to local relative errors, and those [28, 29, 38] based on the changing rates of free energy.
It is known that the coarsening dynamics usually experience large changing rates of free energy in some time intervals and small energy changing rates during some others.
Therefore, we here adopt the following robust time-adaptive strategy based on the energy variation [38]:

where $\Dt_{min},\Dt_{max}$ are predetermined minimum and maximum time step sizes,
$\gamma$ is a positive constant to be determined.
Obviously, the adaptive time-stepping scheme will automatically select small time step sizes when the energy changing rate is large, and large time step sizes when small.  We will make use of the TDSR-ETD3 schemes with $\theta=10^{4}$ in combination with the above time-adaptive strategy (<ref>).

(a) The fixed time step size $\tau=0.1$

(b)  The time-adaptive strategy  (<ref>)   with $\tau_{min}=10^{-4}, \tau_{max}=0.1$, and $\gamma=10^{5}$

(c) The fixed time step size $\tau=10^{-4}$
Snapshots of the simulated coarsening process produced by the TDSR-ETD3 scheme (<ref>)  for  the $H^{-1}$ gradient flow at $t=0.1, 10, 20, 100, 200$, respectively.

§.§.§ Coarsening dynamic by the Cahn-Hilliard equation

We first consider  a long-time coarsening dynamic process governed by  $H^{-1}$ gradient flow (<ref>) (i.e., the Cahn-Hilliard equation) with the double-well potential, subject to the periodic boundary condition with $\varepsilon^{2}=0.01$ and the domain $\Omega=(0,2\pi)^{2}$.
The initial condition   is given by random data ranging from $-0.05$ to $0.05$ in the domain, i.e., $\phi(x,y,0) = 0.05\,(2\,\text{rand}(\cdot)-1)$ where
$\text{rand}(\cdot)$ randomly generates a number between 0 and 1.
The simulation of the coarsening dynamics  is performed by the TDSR-ETD3 scheme with $128\times128$ Fourier modes.
Figure <ref> presents a comparison on the numerical solutions (up to $T=200$) produced by
the uniform time steps with $\tau=0.1$, the time-adaptive strategy (<ref>) with parameters $\Dt_{min}=10^{-4},\Dt_{max}=0.1,$ and $\gamma=10^{5}$, and the uniform time stepping  with $\tau=10^{-4}$.
It is observed that TDSR-ETD3 with  the  large time step size  $\tau=0.1$  leads to an inaccurate phase transition process, while
the time-adaptive strategy gives the correct coarsening process which matches very well with the results produced by  the small time step size  $\Dt=10^{-4}$.
The evolution of the time step size used by TDSR-ETD3 with the time-adaptive strategy (<ref>) is shown in Figure <ref>-(a).
We observe that  the large time step size $\tau=0.1$ is chosen during the time interval $[40,200]$, which clearly demonstrates the efficiency of TDSR-ETD3 with the time-adaptive strategy.
The evolutions of
the number of Picard iterations at each time step are  displayed in Figure <ref>-(b) for all the above three different types of time stepping schemes.
It shows that TDSR-ETD3 with the time-adaptive strategy also significantly  reduces the number of Picard iterations at each time step compared to the case of fixed large time step size $\Dt=0.1$ while still maintaining  the same accuracy as the case of the fixed small time step size $\Dt=10^{-4}$. Except for a very short initial time period and the time interval $[16,40]$, the number of Picard iterations is only either 2 or 3. In addition, the evolutions of the simulated modified and original free energies are plotted in Figure <ref>-(c) and (d), respectively, which show they are almost identical, and the time-adaptive  strategy and  the uniform time stepping  with $\tau=10^{-4}$ both give accurate results.

(a) The adaptive time step size

(b) The number of Picard iterations

(c) The modified energy

(d) The original energy

Evolutions of   (a) the  time step size by the time-adaptive strategy (<ref>), (b) the number of Picard iterations, (c) the modified energy $E\big[\phi^{n}_{N}\big]+\theta \big[(R^{n})^2-1]$, and (d) the original energy  $E\big[\phi^{n}_{N}\big]$ produced by the TDSR-ETD3 scheme (<ref>)  with  the time-adaptive strategy (<ref>)  for the $H^{-1}$ gradient flow.

§.§.§ Thin film epitaxial growth

We next numerically investigate the power laws of the height growth function $\phi$ governed by the MBE model (<ref>) with slope selection or without slope selection.
The domain  and the width parameter are set to be $\Omega=(0,2\pi)^{2}$ and $\varepsilon =0.03$, respectively. The initial condition is a random state with values ranging from $-0.001$ to $0.001$, i.e., $\phi(x,y,0) = 0.001(2\,\text{rand}(\cdot)-1)$.
The thin film epitaxial growth process usually requires a long time simulation before reaching a steady state, and the free energy (<ref>) of the MBE model usually undergoes some very fast and slow rates of change during the process. Therefore, robust, high-order, and energy stable numerical schemes are crucial for such simulation, in order to guarantee the accuracy  and easily be  employed with time-adaptive strategies to speed up the computations.
We take the TDSR-ETD3 scheme (<ref>) and  the time-adaptive strategy (<ref>) with $\Dt_{min}=10^{-5},\Dt_{max}=10^{-2}$ and $\gamma=100$. The Fourier spectral method with $128\times128$ basis modes is used for spatial discretization.
Figure <ref> presents the evolutions of the energy and the time step size for the models with slope section and without slope selection, respectively.
It is observed that the energy of the MBE model with slope selection  dissipates closely at a level of $O(t^{-1/3})$ and the one of the case without slope selection is approximately $O(-\log_{10}(t)),$ which are consistent well with those reported in [38, 44].
We also see that a large time step size $\Dt=10^{-2}$ is automatically chosen by the time-adaptive strategy in the later period of the simulations.
In Figure <ref>, we display some  snapshots of the height function $\phi$
and its Laplacian $\Delta\phi$ for both models obtained around the times $t=1, 50, 100,$ and $500$.

(a) With slope selection. Left: energy, right: time step size.

(b) Without slope selection. Left: energy, right: time step size

Evolutions of the simulated energy  and the time step size  by the TDSR-ETD3 scheme (<ref>)  with  the time-adaptive strategy (<ref>)  for  the MBE model.

Isolines of the simulated height function $\phi$  and its
Laplacian $\Delta \phi$  for the MBE model with slope selection (left two columns) and without slope selection (right two columns) around the times $t=1$, $50$, $100$, and $500$ (from top to bottom).

§.§.§ Crystal growth in a supercooled liquid

Now we simulate the polycrystal growth governed by the PFC model (<ref>) in a supercooled liquid.

Let us first consider the evolution of three crystallites with different orientations in two dimensions [47], whose initial condition is given by
ϕ(x_l,y_l,0)=0.285+0.446(cos(0.66/√(3)y_l)cos(0.66x_l)-0.5cos(1.32/√(3)y_l)),   l=1,2,3,
where $x_l$ and $y_{l}$ define a local system of cartesian coordinates that is oriented with the crystallite lattice. The corresponding parameters are set to be $\sigma=1$ and $\delta=0.25$.
As shown in the first snapshot of Figure <ref>, the three initial crystallites in different orientations are located in the blocks with a length of $40$ in the domain $\Omega=(0,800)^{2},$  which are obtained via different definitions of the local coordinates $(x_{l},y_{l})$ using the following affine transformation of the global coordinates $(x,y)$:

where $\alpha=-\pi/4, 0, \pi/4$ for $l=1,2,3$, respectively. Some similar numerical examples  can  be found in [17, 22, 46].
We adopt the TDSR-ETD3 scheme (<ref>) with $1024\times1024$ Fourier modes for spatial discretization
and the time-adaptive strategy (<ref>) with $\Dt_{min}=0.02, \Dt_{max}=1$ and $\gamma=10$. The simulation of the polycrystal growth lasts up to $T=2000$ in a supercooled liquid.
In Figure <ref>, we plot the evolutions of the simulated energy and the  time step sizes, which again  validate the stability and the efficiency of the proposed scheme (<ref>) with the time-adaptive strategy (<ref>).
Figure <ref> displays  snapshots  of the crystal growth  around the times $t=0, 100, 200,300,400,500,600,700,800,900,1000$, and $2000$.
The numerical  results match very well with those reported in  [47],  which was produced by the relaxed SAV-BDF2 scheme with the fixed small time step size $\tau=0.02$.
The growth of the crystalline phase and the movement of the well-defined crystal-liquid interface are clearly observed, and
furthermore, it shows that different arrangements of crystallites may lead to defects and dislocations. Such crystal growth phenomenon has also been observed in [17, 22, 46].

Evolutions of the simulated energy (left) and the time step size (right) by the TDSR-ETD3 scheme (<ref>) with the time-adaptive strategy (<ref>) for the growth of three crystallites governed by the 2D PFC model.

Snapshots of the simulated growth of three crystallites  governed by the 2D PFC model around the times $t=0, 100, 200, 300, 400, 500, 600, 700, 800,$ $900, 1000$, and $2000$ respectively.

Next we simulate the crystal growth in three dimensions with $\Omega=(0,50)^{3}$ and initial data $\phi(x,y,z,0)=\overline{\phi}+0.01 \text{rand}\,(\cdot)$ where $\overline{\phi}$ is a constant.
The parameters  are set to be $\sigma=1$ and $\delta=0.56$ and $64\times64\times64$ Fourier modes are applied to the spatial discretization in the  TDSR-ETD3 scheme (<ref>). In Figure <ref>,  snapshots of the iso-surfaces of the $\phi=0$ and the density field $\phi$ at the time $t=3000$ are displayed for the PFC model with different values of
$\overline{\phi}=-0.2,-0.35$, and $-0.43$. These results are consistent with those reported in [22, 37, 47].

Snapshots of the simulated  iso-surfaces of $\phi=0$  (top row) and   density field $\phi$ (bottom row) at $T=3000$ for the 3D PFC model with $\overline{\phi}=-0.2, -0.35$, and $-0.43$ (from left to right), respectively.

§ CONCLUSIONS

In this paper we have developed and analyzed a family of  high-order energy dissipative TDSR-ETD schemes for solving general gradient flows, utilizing a combination of the time-dependent spectral renormalization approach and the exponential integrator method. By incorporating the TDSR factor into the system, the discrete energy dissipation laws for the proposed schemes can be easily derived. In addition, the introduction of a new  enforcing term for updating  the TDSR factor significantly improves the time-step restriction to the numerical solution of  the resulting nonlinear systems. Extensive numerical experiments in two and three dimensions have been conducted to demonstrate the accuracy and efficiency of the proposed TDSR-ETD method, especially when adopting a time-adaptive strategy based on energy variation. It is also worth noting that the complexity and nonlinearity of the proposed schemes require nontrivial error analysis, which will require further deeper study.
Moreover, the nonlinear term was treated implicitly in our proposed TDSR-ETD schemes, leading to the need of solving a nonlinear system at each time step. Therefore, it would also be  very interesting and important to investigate whether the proposed TDSR approach can be further extended to  construct  linear high-order schemes with similar energy dissipative properties.

§ ACKNOWLEDGMENTS

We would like to express our sincere appreciation to Prof. Ziad Musslimani of Florida State University for his inspiring work and presentations on  the spectral renormalization method in physics, which stimulated the main idea behind our research.

[1]
M. J. Ablowitz, N. Antar, İ. Bakırtaş, and B. Ilan, Band-gap boundaries and fundamental solitons in complex two-dimensional
nonlinear lattices, Phys. Rev. A, 81 (2010), p. 033834.

[2]
M. J. Ablowitz, N. Antar, İ. Bakırtaş, and B. Ilan, Vortex and dipole solitons in complex two-dimensional nonlinear lattices,
Phys. Rev. A, 86 (2012), p. 033804.

[3]
M. J. Ablowitz, A. S. Fokas, and Z. H. Musslimani, On a new
non-local formulation of water waves, J. Fluid Mech., 562 (2006),
pp. 313–343.

[4]
M. J. Ablowitz and Z. H. Musslimani, Spectral renormalization method
for computing self-localized solutions to nonlinear systems, Opt. lett., 30
(2005), pp. 2140–2142.

[5]
G. Akrivis, B. Li, and D. Li, Energy-decaying extrapolated RK–SAV
methods for the Allen–Cahn and Cahn–Hilliard equations, SIAM J. Sci.
Comput., 41 (2019), pp. A3703–A3727.

[6]
S. M. Allen and J. W. Cahn, A microscopic theory for antiphase
boundary motion and its application to antiphase domain coarsening, Acta
Metall. Mater., 27 (1979), pp. 1085–1095.

[7]
D. M. Anderson, G. B. Mcfadden, and A. A. Wheeler, Diffuse-interface
methods in fluid mechanics, Annu. Rev. Fluid Mech., 30 (1997), pp. 139–165.

[8]
S. Badia, F. Guill$\acute{e}$n-Gonz$\acute{a}$lez, and J. V.
Guti$\acute{e}$rrez-Santacreu, Finite element approximation of nematic
liquid crystal flows using a saddle-point structure, J. Comput. Phys., 230
(2011), pp. 1686–1706.

[9]
J. W. Barrett and J. F. Blowey, Finite element approximation of a
model for phase separation of a multi–component alloy with non-smooth free
energy, Numer. Math., 77 (1997), pp. 1–34.

[10]
A. Baskaran, J. S. Lowengrub, C. Wang, and S. M. Wise, Convergence
analysis of a second order convex splitting scheme for the modified phase
field crystal equation, SIAM J. Numer. Anal., 51 (2013), pp. 2851–2873.

[11]
J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System.
I. Interfacial Free Energy, J. Chem. Phys., 28 (1958), pp. 258–267.

[12]
S. Chandramouli, A. Farhat, and Z. H. Musslimani, Time-dependent
duhamel renormalization method with multiple conservation and dissipation
laws, Nonlinearity, 35 (2022), p. 1286.

[13]
J. T. Cole and Z. H. Musslimani, Time-dependent spectral
renormalization method, Physica D: Nonlinear Phenomena, 358 (2017),
pp. 15–24.

[14]
Q. Du, L. Ju, X. Li, and Z. Qiao, Stabilized linear semi-implicit
schemes for the nonlocal Cahn–Hilliard equation, J. Comput. Phys., 363
(2018), pp. 39–54.

[15]
Q. Du, L. Ju, X. Li, and Z. Qiao, Maximum bound principles for a
class of semilinear parabolic equations and exponential time differencing
schemes, SIAM Rev., 63 (2021), pp. 317–359.

[16]
K. R. Elder and M. Grant, Modeling elastic and plastic deformations
in nonequilibrium processing using phase field crystals, Phys. Rev. E, 70
(2004), p. 051605.

[17]
K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, Modeling
elasticity in crystal growth, Phys. Rev. Lett., 88 (2002), p. 245701.

[18]
C. M. Elliott and A. M. Stuart, The global dynamics of discrete
semilinear parabolic equations, SIAM J. Numer. Anal., 30 (1993),
pp. 1622–1663.

[19]
D. J. Eyre, Unconditionally Gradient Stable Time Marching the
Cahn-Hilliard Equation, Mater. Res. Soc. Sympos. Proc., 529 (1998),
pp. 39–46.

[20]
L. Golubović, Interfacial coarsening in epitaxial growth models
without slope selection, Phys. Rev. Lett., 78 (1997), p. 90.

[21]
H. Gomez and T. J. Hughes, Provably unconditionally stable,
second-order time-accurate, mixed variational methods for phase-field
models, J. Comput. Phys., 230 (2011), pp. 5310–5327.

[22]
H. Gomez and X. Nogueira, An unconditionally energy-stable method
for the phase field crystal equation, Comput. Methods Appl. Mech. Eng.,
249-252 (2012), pp. 52–61.

[23]
Z. Guan, J. S. Lowengrub, C. Wang, and S. M. Wise, Second order
convex splitting schemes for periodic nonlocal Cahn-Hilliard and Allen-Cahn
equations, J. Comput. Phys., 277 (2014), pp. 48–71.

[24]
F. Guillén-González and G. Tierra, On linear schemes for a
Cahn–Hilliard diffuse interface model, J. Comput. Phys., 234 (2013),
pp. 140–171.

[25]
M. E. Gurtin and D. Polignone, Two-phase binary fluids and
immiscible fluids described by an order parameter, Math. Models Methods
Appl. Sci., 6 (1996), pp. 815–831.

[26]
M. A. Hoefer and B. Ilan, Dark solitons, dispersive shock waves, and
transverse instabilities, Multiscale Model. Simul., 10 (2012), pp. 306–341.

[27]
D. Hou, M. Azaiez, and C. Xu, A variant of scalar auxiliary variable
approaches for gradient flows, J. Comput. Phys., 395 (2019), pp. 307–332.

[28]
D. Hou, L. Ju, and Z. Qiao, A linear second-order maximum bound
principle-preserving BDF scheme for the Allen-Cahn equation with a general
mobility, to be appeared in Math. Comput., arXiv:2211.00852v2,  (2023),
pp. 1–28.

[29]
D. Hou and Z. Qiao, An implicit–explicit second order BDF scheme
with variable steps for gradient flows, J. Sci. Comput., 94 (2023), p. 39.

[30]
F. Huang, J. Shen, and Z. Yang, A highly efficient and accurate new
scalar auxiliary variable approach for gradient flows, SIAM J. Sci.
Comput., 42 (2020), pp. A2514–A2536.

[31]
M. Jiang, Z. Zhang, and J. Zhao, Improving the accuracy and
consistency of the scalar auxiliary variable (SAV) method with relaxation,
J. Comput. Phys., 456 (2022), p. 110954.

[32]
L. Ju, X. Li, Z. Qiao, and H. Zhang, Energy stability and error
estimates of exponential time differencing schemes for the epitaxial growth
model without slope selection, Math. Comput., 87 (2021), pp. 1859–1885.

[33]
B. Li and J. Liu, Thin film epitaxy with or without slope
selection, European J. Appl. Math., 14 (2003), pp. 713–743.

[34]
D. Li, Z. Qiao, and T. Tang, Characterizing the stabilization size
for semi–implicit Fourier–spectral method to phase field equations, SIAM
J. Numer. Anal., 54 (2016), pp. 1653–1681.

[35]
X. Li, Z. Qiao, and C. Wang, Convergence analysis for a stabilized
linear semi-implicit numerical scheme for the nonlocal Cahn–Hilliard
equation, Math. Comput., 90 (2021), pp. 171–188.

[36]
X. Li, Z. Qiao, and C. Wang, Stabilization parameter analysis of a
second–order linear numerical scheme for the nonlocal Cahn–Hilliard
equation, IMA J. Numer. Anal., 00 (2022), pp. 1–26.

[37]
N. Provatas, J. A. Dantzig, B. Athreya, P. Chan, P. Stefanovic,
N. Goldenfeld, and K. R. Elder, Using the phase-field crystal method in
the multi-scale modeling of microstructure evolution, Jom, 59 (2007),
pp. 83–90.

[38]
Z. Qiao, Z. Zhang, and T. Tang, An adaptive time-stepping strategy
for the molecular beam epitaxy models, SIAM J. Sci. Comput., 33 (2011),
pp. 1395–1414.

[39]
J. Shen, T. Tang, and L. L. Wang, Spectral methods, Algorithms,
Analysis and Applications, Springer, Berlin, 2010.

[40]
J. Shen, T. Tang, and J. Yang, On the maximum principle preserving
schemes for the generalized Allen–Cahn equation, Commun. Math. Sci., 14
(2016), pp. 1517–1534.

[41]
J. Shen and J. Xu, Convergence and error analysis for the scalar
auxiliary variable(SAV) schemes to gradient flows, SIAM J. Numer. Anal., 56
(2018), pp. 2895–2912.

[42]
J. Shen, J. Xu, and J. Yang, The scalar auxiliary variable (SAV)
approach for gradient flows, J. Comput. Phys., 353 (2018), pp. 407–416.

[43]
J. Shen and X. Yang, Numerical approximations of Allen-Cahn and
Cahn-Hilliard equations, Discrete Conti. Dyn. Syst. Ser. A, 28 (2010),
pp. 1669–1691.

[44]
C. Xu and T. Tang, Stability analysis of large time-stepping methods
for epitaxial growth models, SIAM J. Numer. Anal., 44 (2006),
pp. 1759–1779.

[45]
X. Yang, Linear, first and second-order, unconditionally energy
stable numerical schemes for the phase field model of homopolymer blends, J.
Comput. Phys., 327 (2016), pp. 294–316.

[46]
X. Yang and D. Han, Linearly first- and second-order,
unconditionally energy stable schemes for the phase field crystal model, J.
Comput. Phys., 88 (2017), pp. 1116–1134.

[47]
Y. Zhang and J. Shen, A generalized SAV approach with relaxation
for dissipative systems, J. Comput. Phys., 464 (2022), p. 111311.

[48]
J. Zhao, Q. Wang, and X. Yang, Numerical approximations for a phase
field dendritic crystal growth model based on the invariant energy
quadratization approach, Int. J. Numer. Meth. Eng., 110 (2017),
pp. 279–300.


* [25] G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, Phys. Lett. B 735 (2014), 90-91 [arXiv:1403.7512 [hep-ph]].
  * [26] P. Athron, C. Balázs, D. H. J. Jacob, W. Kotlarski, D. Stöckinger and H. Stöckinger-Kim, JHEP 09 (2021), 080 [arXiv:2104.03691 [hep-ph]].
  * [27] R. H. Parker, C. Yu, W. Zhong, B. Estey and H. Müller, Science 360 (2018), 191 [arXiv:1812.04130 [physics.atom-ph]].
  * [28] C. Yu, W. Zhong, B. Estey, J. Kwan, R. H. Parker and H. Müller, Annalen Phys. 531 (2019) no.5, 1800346.
  * [29] L. Morel, Z. Yao, P. Cladé and S. Guellati-Khélifa, Nature 588 (2020) no.7836, 61-65.
  * [30] G. F. Giudice, P. Paradisi and M. Passera, JHEP 11 (2012), 113 [arXiv:1208.6583 [hep-ph]].
  * [31] H. Davoudiasl and W. J. Marciano, Phys. Rev. D 98 (2018) no.7, 075011 [arXiv:1806.10252 [hep-ph]].
  * [32] Y. Kahn, G. Krnjaic, S. Mishra-Sharma and T. M. P. Tait, JHEP 05 (2017), 002 [arXiv:1609.09072 [hep-ph]].
  * [33] A. Crivellin, M. Hoferichter and P. Schmidt-Wellenburg, Phys. Rev. D 98 (2018) no.11, 113002 [arXiv:1807.11484 [hep-ph]].
  * [34] J. Liu, C. E. M. Wagner and X. P. Wang, JHEP 03 (2019), 008 [arXiv:1810.11028 [hep-ph]].
  * [35] B. Dutta and Y. Mimura, Phys. Lett. B 790 (2019), 563-567 [arXiv:1811.10209 [hep-ph]].
  * [36] X. F. Han, T. Li, L. Wang and Y. Zhang, Phys. Rev. D 99 (2019) no.9, 095034 [arXiv:1812.02449 [hep-ph]].
  * [37] M. Endo and W. Yin, JHEP 08 (2019), 122 [arXiv:1906.08768 [hep-ph]].
  * [38] J. Kawamura, S. Raby and A. Trautner, Phys. Rev. D 100 (2019) no.5, 055030 [arXiv:1906.11297 [hep-ph]].
  * [39] M. Abdullah, B. Dutta, S. Ghosh and T. Li, Phys. Rev. D 100 (2019) no.11, 115006 [arXiv:1907.08109 [hep-ph]].
  * [40] M. Badziak and K. Sakurai, JHEP 10 (2019), 024 [arXiv:1908.03607 [hep-ph]].
  * [41] A. E. Cárcamo Hernández, S. F. King, H. Lee and S. J. Rowley, Phys. Rev. D 101 (2020) no.11, 115016 [arXiv:1910.10734 [hep-ph]].
  * [42] G. Hiller, C. Hormigos-Feliu, D. F. Litim and T. Steudtner, Phys. Rev. D 102 (2020) no.7, 071901 [arXiv:1910.14062 [hep-ph]].
  * [43] C. Cornella, P. Paradisi and O. Sumensari, JHEP 01 (2020), 158 [arXiv:1911.06279 [hep-ph]].
  * [44] A. E. Cárcamo Hernández, Y. Hidalgo Velásquez, S. Kovalenko, H. N. Long, N. A. Pérez-Julve and V. V. Vien, Eur. Phys. J. C 81 (2021) no.2, 191 [arXiv:2002.07347 [hep-ph]].
  * [45] N. Haba, Y. Shimizu and T. Yamada, PTEP 2020 (2020) no.9, 093B05 [arXiv:2002.10230 [hep-ph]].
  * [46] I. Bigaran and R. R. Volkas, Phys. Rev. D 102 (2020) no.7, 075037 [arXiv:2002.12544 [hep-ph]].
  * [47] S. Jana, V. P. K. and S. Saad, Phys. Rev. D 101 (2020) no.11, 115037 [arXiv:2003.03386 [hep-ph]].
  * [48] L. Calibbi, M. L. López-Ibáñez, A. Melis and O. Vives, JHEP 06 (2020), 087 [arXiv:2003.06633 [hep-ph]].
  * [49] C. H. Chen and T. Nomura, Nucl. Phys. B 964 (2021), 115314 [arXiv:2003.07638 [hep-ph]].
  * [50] J. L. Yang, T. F. Feng and H. B. Zhang, J. Phys. G 47 (2020) no.5, 055004 [arXiv:2003.09781 [hep-ph]].
  * [51] I. Doršner, S. Fajfer and S. Saad, Phys. Rev. D 102 (2020) no.7, 075007 [arXiv:2006.11624 [hep-ph]].
  * [52] C. Hati, J. Kriewald, J. Orloff and A. M. Teixeira, JHEP 07 (2020), 235 [arXiv:2005.00028 [hep-ph]].
  * [53] L. Delle Rose, S. Khalil and S. Moretti, Phys. Lett. B 816 (2021), 136216 [arXiv:2012.06911 [hep-ph]].
  * [54] A. E. C. Hernández, S. F. King and H. Lee, Phys. Rev. D 103 (2021) no.11, 115024 [arXiv:2101.05819 [hep-ph]].
  * [55] A. Bodas, R. Coy and S. J. D. King, Eur. Phys. J. C 81 (2021) no.12, 1065 [arXiv:2102.07781 [hep-ph]].
  * [56] X. F. Han, T. Li, H. X. Wang, L. Wang and Y. Zhang, Phys. Rev. D 104 (2021) no.11, 115001 [arXiv:2104.03227 [hep-ph]].
  * [57] A. E. C. Hernández, S. Kovalenko, M. Maniatis and I. Schmidt, JHEP 10 (2021), 036 [arXiv:2104.07047 [hep-ph]].
  * [58] D. Borah, M. Dutta, S. Mahapatra and N. Sahu, Phys. Rev. D 105 (2022) no.1, 015029 [arXiv:2109.02699 [hep-ph]].
  * [59] L. T. Hue, A. E. Cárcamo Hernández, H. N. Long and T. T. Hong, “Heavy singly charged Higgs bosons and inverse seesaw neutrinos as origins of large $(g-2)_{e,\mu}$ in two higgs doublet models,” arXiv:2110.01356 [hep-ph].
  * [60] P. Langacker, Rev. Mod. Phys. 81 (2009), 1199-1228 [arXiv:0801.1345 [hep-ph]].
  * [61] W. Altmannshofer, S. Gori, S. Profumo and F. S. Queiroz, JHEP 12 (2016), 106 [arXiv:1609.04026 [hep-ph]].
  * [62] F. C. Correia and S. Fajfer, Phys. Rev. D 94 (2016) no.11, 115023 [arXiv:1609.00860 [hep-ph]].
  * [63] F. C. Correia and S. Fajfer, JHEP 10 (2019), 278 [arXiv:1905.03867 [hep-ph]].
  * [64] F. C. Correia and S. Fajfer, JHEP 10 (2019), 279 [arXiv:1905.03872 [hep-ph]].
  * [65] W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, Phys. Rev. D 89 (2014), 095033 [arXiv:1403.1269 [hep-ph]].
  * [66] A. Crivellin, G. D’Ambrosio and J. Heeck, Phys. Rev. Lett. 114 (2015), 151801 [arXiv:1501.00993 [hep-ph]].
  * [67] A. Crivellin, G. D’Ambrosio and J. Heeck, Phys. Rev. D 91 (2015) no.7, 075006 [arXiv:1503.03477 [hep-ph]].
  * [68] D. Aristizabal Sierra, F. Staub and A. Vicente, Phys. Rev. D 92 (2015) no.1, 015001 [arXiv:1503.06077 [hep-ph]].
  * [69] A. Crivellin, L. Hofer, J. Matias, U. Nierste, S. Pokorski and J. Rosiek, Phys. Rev. D 92 (2015) no.5, 054013 [arXiv:1504.07928 [hep-ph]].
  * [70] A. Celis, J. Fuentes-Martin, M. Jung and H. Serodio, Phys. Rev. D 92 (2015) no.1, 015007 [arXiv:1505.03079 [hep-ph]].
  * [71] D. Bhatia, S. Chakraborty and A. Dighe, JHEP 03 (2017), 117 [arXiv:1701.05825 [hep-ph]].
  * [72] J. F. Kamenik, Y. Soreq and J. Zupan, Phys. Rev. D 97 (2018) no.3, 035002 [arXiv:1704.06005 [hep-ph]].
  * [73] C. H. Chen and T. Nomura, Phys. Lett. B 777 (2018), 420-427 [arXiv:1707.03249 [hep-ph]].
  * [74] J. E. Camargo-Molina, A. Celis and D. A. Faroughy, Phys. Lett. B 784 (2018), 284-293 [arXiv:1805.04917 [hep-ph]].
  * [75] L. Darmé, K. Kowalska, L. Roszkowski and E. M. Sessolo, JHEP 10 (2018), 052 [arXiv:1806.06036 [hep-ph]].
  * [76] S. Baek and C. Yu, JHEP 11 (2018), 054 [arXiv:1806.05967 [hep-ph]].
  * [77] A. Biswas and A. Shaw, JHEP 05 (2019), 165 [arXiv:1903.08745 [hep-ph]].
  * [78] B. C. Allanach and J. Davighi, Eur. Phys. J. C 79 (2019) no.11, 908 [arXiv:1905.10327 [hep-ph]].
  * [79] A. Crivellin, C. A. Manzari, M. Alguero and J. Matias, Phys. Rev. Lett. 127 (2021) no.1, 011801 [arXiv:2010.14504 [hep-ph]].
  * [80] A. Crivellin, C. A. Manzari, W. Altmannshofer, G. Inguglia, P. Feichtinger and J. M. Camalich, “Towards excluding a light $Z^{\prime}$ explanation of $b\to s\ell^{+}\ell^{-}$,” arXiv:2202.12900 [hep-ph].
  * [81] M. Raggi [NA48/2], Nuovo Cim. C 38 (2016) no.4, 132 [arXiv:1508.01307 [hep-ex]].
  * [82] A. Anastasi, D. Babusci, G. Bencivenni, M. Berlowski, C. Bloise, F. Bossi, P. Branchini, A. Budano, L. Caldeira Balkeståhl and B. Cao, et al. Phys. Lett. B 750 (2015), 633-637 [arXiv:1509.00740 [hep-ex]].
  * [83] D. Androić et al. [Qweak], Nature 557 (2018) no.7704, 207-211 [arXiv:1905.08283 [nucl-ex]].
  * [84] S. Baek, N. G. Deshpande, X. G. He and P. Ko, Phys. Rev. D 64 (2001), 055006 [arXiv:hep-ph/0104141 [hep-ph]].
  * [85] E. Ma, D. P. Roy and S. Roy, Phys. Lett. B 525 (2002), 101-106 [arXiv:hep-ph/0110146 [hep-ph]].
  * [86] D. W. P. Amaral, D. G. Cerdeno, A. Cheek and P. Foldenauer, Eur. Phys. J. C 81 (2021) no.10, 861 [arXiv:2104.03297 [hep-ph]].
  * [87] J. Heeck and W. Rodejohann, Phys. Rev. D 84 (2011), 075007 [arXiv:1107.5238 [hep-ph]].
  * [88] P. Foldenauer, Phys. Rev. D 99 (2019) no.3, 035007 [arXiv:1808.03647 [hep-ph]].
  * [89] S. N. Gninenko and N. V. Krasnikov, Phys. Lett. B 783 (2018), 24-28 [arXiv:1801.10448 [hep-ph]].
  * [90] I. Holst, D. Hooper and G. Krnjaic, Phys. Rev. Lett. 128 (2022) no.14, 141802 [arXiv:2107.09067 [hep-ph]].
  * [91] G. y. Huang, F. S. Queiroz and W. Rodejohann, Phys. Rev. D 103 (2021) no.9, 095005 [arXiv:2101.04956 [hep-ph]].
  * [92] H. Banerjee, B. Dutta and S. Roy, JHEP 03 (2021), 211 [arXiv:2011.05083 [hep-ph]].
  * [93] A. Kamada, K. Kaneta, K. Yanagi and H. B. Yu, JHEP 06 (2018), 117 [arXiv:1805.00651 [hep-ph]].
  * [94] T. Araki, S. Hoshino, T. Ota, J. Sato and T. Shimomura, Phys. Rev. D 95 (2017) no.5, 055006 [arXiv:1702.01497 [hep-ph]].
  * [95] T. Araki, F. Kaneko, Y. Konishi, T. Ota, J. Sato and T. Shimomura, Phys. Rev. D 91 (2015) no.3, 037301 [arXiv:1409.4180 [hep-ph]].
  * [96] S. Patra, S. Rao, N. Sahoo and N. Sahu, Nucl. Phys. B 917 (2017), 317-336 [arXiv:1607.04046 [hep-ph]].
  * [97] K. Harigaya, T. Igari, M. M. Nojiri, M. Takeuchi and K. Tobe, JHEP 03 (2014), 105 [arXiv:1311.0870 [hep-ph]].
  * [98] M. Frank, Y. Hiçyılmaz, S. Mondal, Ö. Özdal and C. S. Ün, JHEP 10 (2021), 063 [arXiv:2107.04116 [hep-ph]].
  * [99] G. Bhattacharyya, M. Dutra, Y. Mambrini and M. Pierre, Phys. Rev. D 98 (2018) no.3, 035038 [arXiv:1806.00016 [hep-ph]].
  * [100] G. Arcadi, P. Ghosh, Y. Mambrini, M. Pierre and F. S. Queiroz, JCAP 11 (2017), 020 [arXiv:1706.04198 [hep-ph]].
  * [101] G. Arcadi, Y. Mambrini, M. H. G. Tytgat and B. Zaldivar, JHEP 03 (2014), 134 [arXiv:1401.0221 [hep-ph]].
  * [102] E. Dudas, Y. Mambrini, S. Pokorski and A. Romagnoni, JHEP 08 (2009), 014 [arXiv:0904.1745 [hep-ph]].
  * [103] E. J. Chun, J. C. Park and S. Scopel, JHEP 02 (2011), 100 [arXiv:1011.3300 [hep-ph]].
  * [104] M. T. Frandsen, F. Kahlhoefer, S. Sarkar and K. Schmidt-Hoberg, JHEP 09 (2011), 128 [arXiv:1107.2118 [hep-ph]].
  * [105] A. Alves, S. Profumo and F. S. Queiroz, JHEP 04 (2014), 063 [arXiv:1312.5281 [hep-ph]].
  * [106] R. Foot, X. G. He, H. Lew and R. R. Volkas, Phys. Rev. D 50 (1994), 4571-4580 [arXiv:hep-ph/9401250 [hep-ph]].
  * [107] J. Heeck, Phys. Lett. B 758 (2016), 101-105 [arXiv:1602.03810 [hep-ph]].
  * [108] W. Altmannshofer, C. Y. Chen, P. S. Bhupal Dev and A. Soni, Phys. Lett. B 762 (2016), 389-398 [arXiv:1607.06832 [hep-ph]].
  * [109] A. J. Buras, A. Crivellin, F. Kirk, C. A. Manzari and M. Montull, JHEP 06 (2021), 068 [arXiv:2104.07680 [hep-ph]].
  * [110] A. Crivellin, S. Najjari and J. Rosiek, JHEP 04 (2014), 167 [arXiv:1312.0634 [hep-ph]].
  * [111] R. K. Ellis, B. Heinemann, J. de Blas, M. Cepeda, C. Grojean, F. Maltoni, A. Nisati, E. Petit, R. Rattazzi and W. Verkerke, et al. “Physics Briefing Book: Input for the European Strategy for Particle Physics Update 2020,” arXiv:1910.11775 [hep-ex].
  * [112] J. P. Delahaye, M. Diemoz, K. Long, B. Mansoulié, N. Pastrone, L. Rivkin, D. Schulte, A. Skrinsky and A. Wulzer, “Muon Colliders,” arXiv:1901.06150 [physics.acc-ph].
  * [113] K. Long, D. Lucchesi, M. Palmer, N. Pastrone, D. Schulte and V. Shiltsev, Nature Phys. 17 (2021) no.3, 289-292 [arXiv:2007.15684 [physics.acc-ph]].
  * [114] J. De Blas et al. [Muon Collider], “The physics case of a 3 TeV muon collider stage,” arXiv:2203.07261 [hep-ph].
  * [115] D. Banerjee et al. [NA64], Phys. Rev. Lett. 120 (2018) no.23, 231802 [arXiv:1803.07748 [hep-ex]].
  * [116] D. Banerjee et al. [NA64], Phys. Rev. D 101 (2020) no.7, 071101 [arXiv:1912.11389 [hep-ex]].
  * [117] J. P. Lees et al. [BaBar], Phys. Rev. Lett. 113 (2014) no.20, 201801 [arXiv:1406.2980 [hep-ex]].
  * [118] P. L. Anthony et al. [SLAC E158], Phys. Rev. Lett. 95 (2005), 081601 [arXiv:hep-ex/0504049 [hep-ex]].
  * [119] M. Deniz et al. [TEXONO], Phys. Rev. D 81 (2010), 072001 [arXiv:0911.1597 [hep-ex]].
  * [120] P. Vilain et al. [CHARM-II], Phys. Lett. B 302 (1993), 351-355.
  * [121] S. R. Mishra et al. [CCFR], Phys. Rev. Lett. 66 (1991), 3117-3120.
  * [122] K. I. Aoki, Z. Hioki, M. Konuma, R. Kawabe and T. Muta, Prog. Theor. Phys. Suppl. 73 (1982), 1-225.
  * [123] D. Espriu, J. Manzano and P. Talavera, Phys. Rev. D 66 (2002), 076002 [arXiv:hep-ph/0204085 [hep-ph]].
  * [124] A. Denner, Fortsch. Phys. 41 (1993), 307-420 [arXiv:0709.1075 [hep-ph]].
  * [125] A. Freitas, JHEP 04 (2014), 070 [arXiv:1401.2447 [hep-ph]].
  * [126] P. A. Zyla et al. [Particle Data Group], PTEP 2020 (2020) no.8, 083C01.
  * [127] H. H. Patel, Comput. Phys. Commun. 197 (2015), 276-290 [arXiv:1503.01469 [hep-ph]].
  * [128] G. Aad et al. [ATLAS], Phys. Rev. D 90 (2014) no.7, 072010 [arXiv:1408.5774 [hep-ex]].
  * [129] R. Akers et al. [OPAL], Z. Phys. C 67 (1995), 555-564.
  * [130] A. M. Baldini et al. [MEG], Eur. Phys. J. C 76 (2016) no.8, 434 [arXiv:1605.05081 [hep-ex]].
  * [131] A. M. Baldini et al. [MEG II], Eur. Phys. J. C 78 (2018) no.5, 380 [arXiv:1801.04688 [physics.ins-det]].
  * [132] B. Aubert et al. [BaBar], Phys. Rev. Lett. 104 (2010), 021802 [arXiv:0908.2381 [hep-ex]].
  * [133] E. Kou et al. [Belle-II], PTEP 2019 (2019) no.12, 123C01 [erratum: PTEP 2020 (2020) no.2, 029201] [arXiv:1808.10567 [hep-ex]].
  * [134] U. Bellgardt et al. [SINDRUM], Nucl. Phys. B 299 (1988), 1-6.
  * [135] A. Blondel, A. Bravar, M. Pohl, S. Bachmann, N. Berger, M. Kiehn, A. Schoning, D. Wiedner, B. Windelband and P. Eckert, et al. “Research Proposal for an Experiment to Search for the Decay $\mu\to eee$,” arXiv:1301.6113 [physics.ins-det].
  * [136] K. Hayasaka, K. Inami, Y. Miyazaki, K. Arinstein, V. Aulchenko, T. Aushev, A. M. Bakich, A. Bay, K. Belous and V. Bhardwaj, et al. Phys. Lett. B 687 (2010), 139-143 [arXiv:1001.3221 [hep-ex]].
  * [137] A. Abada et al. [FCC], Eur. Phys. J. C 79 (2019) no.6, 474.
  * [138] W. H. Bertl et al. [SINDRUM II], Eur. Phys. J. C 47 (2006), 337-346.
  * [139] T. M. Nguyen [DeeMe], PoS FPCP2015 (2015), 060.
  * [140] B. E. Krikler [COMET], “An Overview of the COMET Experiment and its Recent Progress,” arXiv:1512.08564 [physics.ins-det].
  * [141] Y. Kuno, Presentation at the Flavour Session of the CERN Council Open Symposium on the Update of the European Strategy for Particle Physics. Granada, 13-16 May 2019.
  * [142] R. Abramishvili et al. [COMET], PTEP 2020 (2020) no.3, 033C01 [arXiv:1812.09018 [physics.ins-det]].
  * [143] L. Bartoszek et al. [Mu2e], “Mu2e Technical Design Report,” arXiv:1501.05241 [physics.ins-det].
  * [144] G. Aad et al. [ATLAS], Phys. Rev. D 90 (2014) no.7, 072010 [arXiv:1408.5774 [hep-ex]].
  * [145] R. Akers et al. [OPAL], Z. Phys. C 67 (1995), 555-564.
  * [146] A. Abada, M. E. Krauss, W. Porod, F. Staub, A. Vicente and C. Weiland, JHEP 11 (2014), 048 [arXiv:1408.0138 [hep-ph]].
  * [147] T. Hahn and M. Perez-Victoria, Comput. Phys. Commun. 118 (1999), 153-165 [arXiv:hep-ph/9807565 [hep-ph]].
  * [148] G. Passarino and M. J. G. Veltman, Nucl. Phys. B 160 (1979), 151-207.
  * [149] J. Hisano, T. Moroi, K. Tobe and M. Yamaguchi, Phys. Rev. D 53 (1996), 2442-2459 [arXiv:hep-ph/9510309 [hep-ph]].
  * [150] R. Kitano, M. Koike and Y. Okada, Phys. Rev. D 66 (2002), 096002 [erratum: Phys. Rev. D 76 (2007), 059902] [arXiv:hep-ph/0203110 [hep-ph]].
  * [151] T. Fukuyama, Y. Mimura and Y. Uesaka, Phys. Rev. D 105 (2022) no.1, 015026 [arXiv:2108.10736 [hep-ph]].
  * [152] L. Willmann, P. V. Schmidt, H. P. Wirtz, R. Abela, V. Baranov, J. Bagaturia, W. H. Bertl, R. Engfer, A. Grossmann and V. W. Hughes, et al. Phys. Rev. Lett. 82 (1999), 49-52 [arXiv:hep-ex/9807011 [hep-ex]].
  * [153] H. Albrecht et al. [ARGUS], Z. Phys. C 53 (1992), 367-374.
  * [154] A. Anastassov et al. [CLEO], Phys. Rev. D 55 (1997), 2559-2576 [erratum: Phys. Rev. D 58 (1998), 119904].
  * [155] B. Aubert et al. [BaBar], Phys. Rev. Lett. 105 (2010), 051602 [arXiv:0912.0242 [hep-ex]].
  * [156] A. Pich, I. Boyko, D. Dedovich and I. I. Bigi, Int. J. Mod. Phys. A 24S1 (2009), 715-737.
  * [157] Y. S. Amhis et al. [HFLAV], Eur. Phys. J. C 81 (2021) no.3, 226 [arXiv:1909.12524 [hep-ex]].
  * [158] A. Sirlin, Phys. Rev. D 22 (1980), 971-981.
  * [159] M. Lindner, M. Platscher and F. S. Queiroz, Phys. Rept. 731 (2018), 1-82 [arXiv:1610.06587 [hep-ph]].
  * [160] F. Jegerlehner and A. Nyffeler, Phys. Rept. 477 (2009), 1-110 [arXiv:0902.3360 [hep-ph]].
  * [161] J. P. Leveille, Nucl. Phys. B 137 (1978), 63-76.
  * [162] S. Iguro, Y. Omura and M. Takeuchi, JHEP 09 (2020), 144 [arXiv:2002.12728 [hep-ph]].
  * [163] A. Abada, J. Kriewald and A. M. Teixeira, Eur. Phys. J. C 81 (2021) no.11, 1016 [arXiv:2107.06313 [hep-ph]].
  * [164] S. Schael et al. [ALEPH], Eur. Phys. J. C 49 (2007), 411-437 [arXiv:hep-ex/0609051 [hep-ex]].
  * [165] G. Aad et al. [ATLAS], Phys. Rev. Lett. 115 (2015) no.3, 031801 [arXiv:1503.04430 [hep-ex]].
  * [166] G. Aad et al. [ATLAS], JHEP 03 (2020), 145 [arXiv:1910.08447 [hep-ex]].
  * [167] A. M. Sirunyan et al. [CMS], JHEP 05 (2020), 033 [arXiv:1911.03947 [hep-ex]].
  * [168] A. M. Sirunyan et al. [CMS], JHEP 07 (2021), 208 [arXiv:2103.02708 [hep-ex]].
  * [169] G. Aad et al. [ATLAS], Phys. Lett. B 796 (2019), 68-87 [arXiv:1903.06248 [hep-ex]].
  * [170] S. Iguro, K. A. Mohan and C. P. Yuan, Phys. Rev. D 101 (2020) no.7, 075011 [arXiv:2001.09079 [hep-ph]].
  * [171] C. Aimè, A. Apyan, M. Attia Mahmoud, N. Bartosik, F. Batsch, A. Bertolin, M. Bonesini, D. Buttazzo, M. Casarsa and M. G. Catanesi, et al. “Muon Collider Physics Summary,” arXiv:2203.07256 [hep-ph].
  * [172] D. Stratakis et al. [International Muon Collider], “A Muon Collider Facility for Physics Discovery,” arXiv:2203.08033 [physics.acc-ph].
  * [173] J. Blümlein and K. Schönwald, Mod. Phys. Lett. A 37 (2022) no.07, 2230004 [arXiv:2202.08476 [hep-ph]].
  * [174] E. Accomando, A. Belyaev, J. Fiaschi, K. Mimasu, S. Moretti and C. Shepherd-Themistocleous, JHEP 01 (2016), 127 [arXiv:1503.02672 [hep-ph]].
  * [175] M. Endo, S. Iguro and T. Kitahara, JHEP 06 (2020), 040 [arXiv:2002.05948 [hep-ph]].
  * [176] S. Schael et al. [ALEPH, DELPHI, L3, OPAL, SLD, LEP Electroweak Working Group, SLD Electroweak Group and SLD Heavy Flavour Group], Phys. Rept. 427 (2006), 257-454 [arXiv:hep-ex/0509008 [hep-ex]].
  * [177] A. Tewsley-Booth [Muon g-2], “The First Measurement of the Muon Anomalous Magnetic Moment from the Fermilab Muon $g-2$ Collaboration,” arXiv:2202.11172 [hep-ex].
  * [178] A. Abada, A. M. Teixeira, A. Vicente and C. Weiland, JHEP 02 (2014), 091 [arXiv:1311.2830 [hep-ph]].